StRoot: StMuDSTMaker/COMMON/StMuEmcUtil.cxx Source File

00001 #include "StMuEmcUtil.h"
00002 #include "StEvent.h"
00003 #include "StMessMgr.h"
00004 #include "StEventTypes.h"
00005 #include "StMuEmcCollection.h"
00006 #include "SystemOfUnits.h"
00007 #include "StEmcUtil/geometry/StEmcGeom.h"
00008 #include "StMuEmcTowerData.h"
00009 
00010 ClassImp(StMuEmcUtil)
00011 
00012 #define __EMC_HITS_ID_DIM__ 18000
00013 
00014 
00015 StMuEmcUtil::StMuEmcUtil()
00016     : TObject()
00017 {
00018   for(Int_t i =0;i<4;i++) mGeo[i]=StEmcGeom::getEmcGeom(i+1);
00019 }
00020 StMuEmcUtil::~StMuEmcUtil()
00021 {
00022 }
00023 StMuEmcCollection* StMuEmcUtil::getMuEmc(const StEmcCollection *emccol)
00024 {
00025   if(!emccol) return NULL;
00026   StMuEmcCollection* muEmc=new StMuEmcCollection();
00027   fillMuEmc(muEmc,emccol);
00028   return muEmc;
00029 }  
00030 StEmcCollection* StMuEmcUtil::getEmc(const StMuEmcCollection* muEmc)
00031 {
00032   if(!muEmc) return NULL;
00033   StEmcCollection *emc=new StEmcCollection();
00034   fillEmc(emc,muEmc);
00035   return emc;
00036 }
00037 void StMuEmcUtil::fillMuEmc(StMuEmcCollection *muEmc, const StEmcCollection *emccol)
00038 {
00039   if(!emccol) return;
00040   if(!muEmc) return;
00041       
00042   // starting by hits;    
00043   //cout <<"Filling hits and clusters \n";
00044   for(Int_t d=0; d<8; d++)
00045   {  
00046     Int_t EmcDet=d+1;
00047     
00048     StDetectorId id = static_cast<StDetectorId>(d+kBarrelEmcTowerId);
00049     const StEmcDetector* detector=emccol->detector(id);
00050     if(detector)
00051     {                          
00052       Int_t maxMod = 121;
00053       if(d>3) maxMod = 14;
00054       //cout <<"Filling hits for detetor "<<EmcDet<<endl;
00055       for(Int_t j=1;j<maxMod;j++)
00056       {
00057         const StEmcModule* module = detector->module(j);
00058         if(module)
00059         {
00060           const StSPtrVecEmcRawHit& rawHit=module->hits();
00061           Int_t nhits=(Int_t)rawHit.size();
00062           if(nhits>0)
00063             for(Int_t k=0;k<nhits;k++)
00064             {
00065               Int_t m = rawHit[k]->module();
00066               Int_t e = rawHit[k]->eta();
00067               Int_t s = abs(rawHit[k]->sub());
00068               Int_t adc = rawHit[k]->adc();
00069               Float_t energy = rawHit[k]->energy();
00070               Int_t cal = rawHit[k]->calibrationType();              
00071               Int_t rid;
00072               bool save = true;
00073               if(d<4 && cal>127) save = false;
00074               if(save)
00075               {
00076                 if (d<4) // for the barrel
00077                 {
00078                   mGeo[d]->getId(m,e,s,rid);
00079                 }
00080                 else
00081                 {
00082                   if(getEndcapId(EmcDet,m,e,s,rid)) continue;// on error
00083                 }
00084               
00085                 if(EmcDet == 1 || EmcDet == 5  ) // towers save only ADC
00086                 {
00087                   muEmc->setTowerADC(rid,adc,EmcDet);
00088                 }              
00089                 if(EmcDet==2 || EmcDet == 6) //pre shower
00090                 {
00091                   muEmc->addPrsHit(EmcDet);
00092                   StMuEmcHit* muHit = muEmc->getPrsHit(muEmc->getNPrsHits(EmcDet)-1,EmcDet); 
00093                   muHit->setId(rid);
00094                   muHit->setAdc(adc);
00095                   muHit->setEnergy(energy);
00096                   muHit->setCalType(cal);      
00097                 }
00098                 if(EmcDet==3 || EmcDet==4 || EmcDet==7 || EmcDet==8)
00099                 {
00100                   muEmc->addSmdHit(EmcDet);
00101                   StMuEmcHit* muHit = muEmc->getSmdHit(muEmc->getNSmdHits(EmcDet)-1,EmcDet);          
00102                   muHit->setId(rid);
00103                   muHit->setAdc(adc);
00104                   muHit->setEnergy(energy);
00105                   muHit->setCalType(cal);
00106                 }      
00107               }
00108             }
00109           
00110         }
00111       }
00112       Int_t n_crate=0;
00113       switch (EmcDet) {
00114       case 1:
00115         n_crate=StMuEmcTowerData::nBTowCrates;
00116         break;
00117       case 2:
00118         n_crate=StMuEmcTowerData::nBPrsCrates;
00119         break;
00120       case 3:
00121         n_crate=StMuEmcTowerData::nBSmdCrates;
00122         break;
00123       case 5:
00124         n_crate=StMuEmcTowerData::nETowCrates;
00125         break;
00126       case 6:
00127         n_crate=StMuEmcTowerData::nEPrsCrates;
00128         break;
00129       case 7:
00130         n_crate=StMuEmcTowerData::nESmdCrates;
00131         break;
00132       }
00133       for (Int_t i_crate=1; i_crate<=n_crate; i_crate++) {
00134         muEmc->setCrateStatus(detector->crateStatus(i_crate),i_crate,EmcDet);
00135       }
00136     }
00137   } 
00138   
00139   return;
00140 }
00141 
00142 void StMuEmcUtil::fillEmc(StEmcCollection* emc, const StMuEmcCollection* muEmc)
00143 {
00144   if(!muEmc) return;
00145   if(!emc) return;
00146   //cout <<"FILLING EMC COLLECTION\n";
00147   
00148   for(Int_t i=0;i<8;i++)
00149   {
00150     Int_t det=i+1;
00151     
00152     StDetectorId id = static_cast<StDetectorId>(i+kBarrelEmcTowerId);
00153     int nMod = 120;
00154     if(i>=4) nMod = 13;
00155     StEmcDetector* detector = new StEmcDetector(id, nMod);
00156     emc->setDetector(detector);
00157     // hits
00158     Int_t nh=0;
00159     if (det==1) nh = 4800; 
00160     if (det==5) nh = 720; 
00161     if (det==2 || det ==6) nh=muEmc->getNPrsHits(det);
00162     if (det==3 || det==4 || det==7 || det==8) nh=muEmc->getNSmdHits(det);
00163     //cout <<"Number of hits for detector "<<det<<" = "<<nh<<endl;
00164     for(Int_t j=0;j<nh;j++)
00165     {
00166       Bool_t save = kTRUE;
00167       Int_t m,e,s,rid;
00168       Int_t a=0,cal=0;
00169       Float_t energy=0;
00170       if(det==1 || det==5) // towers have only ADC
00171       {
00172         a = muEmc->getTowerADC(j+1,det);
00173         if(det==1) mGeo[det-1]->getBin(j+1,m,e,s);
00174         else  {
00175           if( getEndcapBin(det,j+1,m,e,s)) continue ;// on error
00176         }
00177         energy = 0;
00178         cal    = 0;
00179         if(a==0) save = kFALSE;
00180       }      
00181       if(det==2 || det ==6) //prs
00182       {
00183         const StMuEmcHit* hit=muEmc->getPrsHit(j,det);
00184         if(hit)
00185         {
00186           rid=hit->getId();
00187           if(det==2) mGeo[det-1]->getBin(rid,m,e,s);
00188           else {
00189             if( getEndcapBin(det,rid,m,e,s)) continue ;// on error
00190           }
00191           a=hit->getAdc();
00192           cal=hit->getCalType();
00193           energy=hit->getEnergy();
00194         } else save = kFALSE;
00195       }
00196       if(det==3 || det==4 || det==7 || det==8) //smd
00197       {
00198         const StMuEmcHit* hit=muEmc->getSmdHit(j,det);
00199         if(hit)
00200         {
00201           rid=hit->getId();
00202           if(det<5) mGeo[det-1]->getBin(rid,m,e,s);
00203           else {
00204             if( getEndcapBin(det,rid,m,e,s)) continue ;// on error
00205           }
00206           a=hit->getAdc();
00207           cal=hit->getCalType();
00208           energy=hit->getEnergy();
00209         } else save = kFALSE;
00210       }
00211       if(save)
00212       {
00213         StEmcRawHit* rawHit=new StEmcRawHit(id,(UInt_t)m,(UInt_t)e,(UInt_t)s,(UInt_t)a,energy);
00214         rawHit->setCalibrationType(cal);
00215         //cout <<"det = "<<det<<"  Hit number "<<j<<"  m = "<<m<<"  e = "<<e<<"  s = "<<s<<"  adc = "<<a<<"  en = "<<energy<<"\n";
00216         detector->addHit(rawHit);
00217       }
00218     }
00219   
00220     Int_t n_crate=0;
00221     switch (det) {
00222     case 1:
00223       n_crate=StMuEmcTowerData::nBTowCrates;
00224       break;
00225     case 2:
00226       n_crate=StMuEmcTowerData::nBPrsCrates;
00227       break;
00228     case 3:
00229     case 4:
00230       n_crate=StMuEmcTowerData::nBSmdCrates;
00231       break;
00232     case 5:
00233       n_crate=StMuEmcTowerData::nETowCrates;
00234       break;
00235     case 6:
00236       n_crate=StMuEmcTowerData::nEPrsCrates;
00237       break;
00238     case 7:
00239     case 8:
00240       n_crate=StMuEmcTowerData::nESmdCrates;
00241       break;
00242     }
00243     for (Int_t i_crate=1; i_crate<=n_crate; i_crate++) {
00244       detector->setCrateStatus(i_crate,muEmc->getCrateStatus(i_crate,det));
00245     }
00246   }
00247   return;
00248 }
00249 
00250 //=================================================
00251 //=================================================
00252 int  StMuEmcUtil::getEndcapId(int d,int m, int e, int s,int &rid) const {
00253   rid=1;
00254   /*  first tower or first strip is default,
00255       I do not like it, but it is the only safe value,JB
00256   */
00257  
00258   TString text;
00259 
00260   if(     m<=0 || m >12 || d<5 || d>8 ) {
00261     text="m<=0 || m >12 || d<5 || d>8 ";
00262     goto abort;
00263   }
00264   
00265   switch (d) {
00266 
00267   case 5: 
00268     if(     e<=0 || e>12 || s<=0 || s>5 ) {
00269       text="e<=0 || e>12 || s<=0 || s>5 ,towers";
00270       goto abort;
00271     } 
00272     rid = 60*(m-1) + 12*(s-1) + e-1;
00273     break;
00274 
00275   case 6: 
00276     if(     e<=0 || e>12 || s<=0 || s>15 ) {
00277       text="e<=0 || e>12 || s<=0 || s>15 ,  pre/post";
00278       goto abort;
00279     } 
00280     rid = 180*(m-1) + 12*(s-1) + e-1;
00281     break;
00282 
00283   case 7: 
00284   case 8:
00285     if(     s!=1 || e<=0 || e>288 ) {
00286       text=" s!=1 || e<=0 || e>288,  SMD";
00287       goto abort;
00288     } 
00289     rid = 300*(m-1) + e-1;
00290     break;
00291   default:;
00292   }
00293 
00294   rid++; // to compensate for counting from 1
00295 
00296   if(     rid<=0 ) {
00297       text=" rid<=0";
00298       goto abort;
00299     } 
00300 
00301   return 0;  // all went OK
00302 
00303  abort:
00304   gMessMgr->Error() <<"StMuEmcUtil::getEndcapId(), FATAL internal error: "
00305                     <<text<<"\n d="<<d<<" m="<<m<<" e="<<" s="<<s<<" rid "<<rid
00306                     <<"\n ENDCAP data may be wrong, " << endm;
00307   return 1;
00308 }
00309 
00310 
00311 int StMuEmcUtil::getEndcapBin(int d,int rid0,int &m, int &e, int &s) const
00312 {
00313   m=e=s=1;
00314   /*  first tower or first strip is default,
00315       I do not like it, but it is the only safe value,JB
00316   */
00317  
00318   TString text;
00319   int rid=rid0-1;
00320   int x;
00321 
00322   if(     rid0<=0 || d<5 || d>8 ) {
00323     text="rid0<=0 || d<5 || d>8 ";
00324     goto abort;
00325   }
00326 
00327   switch (d) {
00328 
00329   case 5: 
00330     m=1+ rid/60;
00331     x=rid%60;
00332     s=1 +x/12;
00333     e=1 + x%12;
00334     if     ( m>12 || s<1 || s>5 || e>12 ) {
00335       text=" m>12 || s<1 || s>5 || e>12 , tower";
00336       goto abort;
00337     }
00338     break;
00339 
00340   case 6:
00341     m=1+ rid/180;
00342     x=rid%180;
00343     s=1 +x/12;
00344     e=1 + x%12;
00345     if     ( m>12 || s<1 || s>15 || e>12 ) {
00346       text=" m>12 || s<1 || s>15 || e>12 , pre/post";
00347       goto abort;
00348     }
00349     break;
00350 
00351   case 7:
00352   case 8:
00353     m=1+ rid/300;
00354     s=1;
00355     e=1 + rid%300;
00356     if     ( m>12 || e>288 || e <1 ) {
00357       text=" m>12 || s<=0 || s>5 || e>12 , smd";
00358       goto abort;
00359     }
00360     break;
00361   default: ;
00362   }
00363   
00364   return 0; // all went OK
00365   
00366  abort:
00367   gMessMgr->Error() <<"StMuEmcUtil::getEndcapBin(), FATAL internal error: "
00368                     <<text<<"\n d="<<d<<" m="<<m<<" e="<<" s="<<s<<" rid0 "<<rid0
00369                     <<"\n ENDCAP data may be wrong, "
00370                     <<" assert() should be here, JB"<<endm;
00371   return 1;
00372 }