#include <EmcRec.h>
Public Member Functions | |
EmcRec (const std::string &name, ISvcLocator *pSvcLocator) | |
EmcRec (const std::string &name, ISvcLocator *pSvcLocator) | |
StatusCode | execute () |
StatusCode | execute () |
StatusCode | finalize () |
StatusCode | finalize () |
StatusCode | initialize () |
StatusCode | initialize () |
Private Attributes | |
NTuple::Item< float > | a20Moment |
NTuple::Item< float > | a20Moment |
NTuple::Item< float > | a42Moment |
NTuple::Item< float > | a42Moment |
NTuple::Item< float > | cluster2ndMoment |
NTuple::Item< float > | cluster2ndMoment |
NTuple::Item< float > | dphi1 |
NTuple::Item< float > | dphi1 |
NTuple::Item< float > | e3x3 |
NTuple::Item< float > | e3x3 |
NTuple::Item< float > | e5x5 |
NTuple::Item< float > | e5x5 |
NTuple::Item< float > | enseed |
NTuple::Item< float > | enseed |
NTuple::Item< float > | eseed |
NTuple::Item< float > | eseed |
NTuple::Item< float > | etof2x1 |
NTuple::Item< float > | etof2x1 |
NTuple::Item< float > | etof2x3 |
NTuple::Item< float > | etof2x3 |
EmcRecCluster2ShowerAbs * | fCluster2Shower |
EmcRecCluster2ShowerAbs * | fCluster2Shower |
RecEmcClusterMap | fClusterMap |
bool | fDigiCalib |
EmcRecDigit2Hit | fDigit2Hit |
RecEmcDigitMap | fDigitMap |
int | fEventNb |
EmcRecHit2Cluster | fHit2Cluster |
RecEmcHitMap | fHitMap |
bool | fOnlineMode |
int | fOutput |
std::vector< std::string > | fPositionMode |
std::vector< std::string > | fPositionMode |
RecEmcShowerMap | fShowerMap |
double | fTimeMax |
double | fTimeMin |
bool | fTofEnergy |
NTuple::Item< float > | latMoment |
NTuple::Item< float > | latMoment |
int | m_event |
PropertyMgr | m_propMgr |
IRawDataProviderSvc * | m_rawDataProviderSvc |
IRawDataProviderSvc * | m_rawDataProviderSvc |
NTuple::Tuple * | m_tuple |
NTuple::Tuple * | m_tuple |
NTuple::Item< float > | mpi0 |
NTuple::Item< float > | mpi0 |
NTuple::Item< long > | ncluster |
NTuple::Item< long > | ncluster |
NTuple::Item< long > | ndigi |
NTuple::Item< long > | ndigi |
NTuple::Item< long > | nhit |
NTuple::Item< long > | nhit |
NTuple::Item< long > | nneu |
NTuple::Item< long > | nneu |
NTuple::Item< long > | npart |
NTuple::Item< long > | npart |
NTuple::Item< long > | nphi |
NTuple::Item< long > | nphi |
NTuple::Item< long > | nrec |
NTuple::Item< long > | nrec |
NTuple::Item< long > | nrun |
NTuple::Item< long > | nrun |
NTuple::Item< long > | ntheta |
NTuple::Item< long > | ntheta |
NTuple::Item< float > | phi1 |
NTuple::Item< float > | phi1 |
NTuple::Item< long > | phigap1 |
NTuple::Item< long > | phigap1 |
NTuple::Item< long > | pid |
NTuple::Item< long > | pid |
NTuple::Array< float > | pp1 |
NTuple::Array< float > | pp1 |
NTuple::Array< float > | pp2 |
NTuple::Array< float > | pp2 |
NTuple::Item< float > | secondMoment |
NTuple::Item< float > | secondMoment |
NTuple::Item< float > | theta1 |
NTuple::Item< float > | theta1 |
NTuple::Item< long > | thtgap1 |
NTuple::Item< long > | thtgap1 |
NTuple::Item< float > | tp |
NTuple::Item< float > | tp |
NTuple::Item< float > | tphi |
NTuple::Item< float > | tphi |
NTuple::Item< float > | ttheta |
NTuple::Item< float > | ttheta |
|
00035 : 00036 Algorithm(name, pSvcLocator),fCluster2Shower(0) 00037 { 00038 m_event=0; 00039 fPositionMode.push_back("log"); 00040 fPositionMode.push_back("5x5"); 00041 // Declare the properties 00042 m_propMgr.declareProperty("Output",fOutput=0); 00043 m_propMgr.declareProperty("EventNb",fEventNb=0); 00044 m_propMgr.declareProperty("DigiCalib",fDigiCalib=false); 00045 m_propMgr.declareProperty("TofEnergy",fTofEnergy=false); 00046 m_propMgr.declareProperty("OnlineMode",fOnlineMode=false); 00047 m_propMgr.declareProperty("TimeMin",fTimeMin=0); 00048 m_propMgr.declareProperty("TimeMax",fTimeMax=60); 00049 m_propMgr.declareProperty("PositionMode",fPositionMode); 00050 00051 IJobOptionsSvc* jobSvc; 00052 service("JobOptionsSvc", jobSvc); 00053 jobSvc->setMyProperties("EmcRecAlg", &m_propMgr); 00054 }
|
|
|
|
|
|
Reconstruction begins. State 1: Initialize digit map Retrieve Tof Digi Retrieve EMC digi State 2: DigitMap --> HitMap State 3: HitMap --> ClusterMap State 4: ClusterMap --> ShowerMap State 5: Emc Tof Match State 6: Register to TDS 00137 { 00138 00139 MsgStream log(msgSvc(), name()); 00140 log << MSG::DEBUG << "in execute()" << endreq; 00141 00144 int event, run; 00145 00146 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader"); 00147 if (!eventHeader) { 00148 log << MSG::FATAL << "Could not find Event Header" << endreq; 00149 return( StatusCode::FAILURE); 00150 } 00151 run=eventHeader->runNumber(); 00152 event=eventHeader->eventNumber(); 00153 if(fEventNb!=0&&m_event%fEventNb==0) { 00154 cout<<m_event<<"-------: "<<run<<","<<event<<endl; 00155 } 00156 m_event++; 00157 if(fOutput>=2) { 00158 cout<<"===================================="<<endl; 00159 cout<<"run= "<<run<<"; event= "<<event<<endl; 00160 } 00161 00162 Hep3Vector posG; //initial photon position 00163 #ifndef OnlineMode 00164 RecEmcTofHitMap tofHitMap; 00165 if(!(m_rawDataProviderSvc->isOnlineMode())) { 00166 if(fOutput>=1) { 00167 // Retrieve mc truth 00168 SmartDataPtr<McParticleCol> mcParticleCol(eventSvc(),"/Event/MC/McParticleCol"); 00169 if (!mcParticleCol) { 00170 log << MSG::WARNING << "Could not find McParticle" << endreq; 00171 } else { 00172 HepLorentzVector pG; 00173 McParticleCol::iterator iter_mc = mcParticleCol->begin(); 00174 for (;iter_mc != mcParticleCol->end(); iter_mc++) { 00175 log << MSG::INFO 00176 << " particleId = " << (*iter_mc)->particleProperty() 00177 << endreq; 00178 pG = (*iter_mc)->initialFourMomentum(); 00179 posG = (*iter_mc)->initialPosition().v(); 00180 } 00181 ttheta = pG.theta(); 00182 tphi = pG.phi(); 00183 if(tphi<0) { 00184 tphi=twopi+tphi; 00185 } 00186 tp = pG.rho(); 00187 00188 // Retrieve EMC truth 00189 SmartDataPtr<EmcMcHitCol> emcMcHitCol(eventSvc(),"/Event/MC/EmcMcHitCol"); 00190 if (!emcMcHitCol) { 00191 log << MSG::WARNING << "Could not find EMC truth" << endreq; 00192 } 00193 00194 RecEmcID mcId; 00195 unsigned int mcTrackIndex; 00196 double mcPosX=0,mcPosY=0,mcPosZ=0; 00197 double mcPx=0,mcPy=0,mcPz=0; 00198 double mcEnergy=0; 00199 00200 EmcMcHitCol::iterator iterMc; 00201 for(iterMc=emcMcHitCol->begin();iterMc!=emcMcHitCol->end();iterMc++){ 00202 mcId=(*iterMc)->identify(); 00203 mcTrackIndex=(*iterMc)->getTrackIndex(); 00204 mcPosX=(*iterMc)->getPositionX(); 00205 mcPosY=(*iterMc)->getPositionY(); 00206 mcPosZ=(*iterMc)->getPositionZ(); 00207 mcPx=(*iterMc)->getPx(); 00208 mcPy=(*iterMc)->getPy(); 00209 mcPz=(*iterMc)->getPz(); 00210 mcEnergy=(*iterMc)->getDepositEnergy(); 00211 00212 if(fOutput>=2){ 00213 cout<<"mcId="<<mcId<<"\t"<<"mcTrackIndex="<<mcTrackIndex<<endl; 00214 cout<<"mcPosition:\t"<<mcPosX<<"\t"<<mcPosY<<"\t"<<mcPosZ<<endl; 00215 cout<<"mcP:\t"<<mcPx<<"\t"<<mcPy<<"\t"<<mcPz<<endl; 00216 cout<<"mcDepositEnergy=\t"<<mcEnergy<<endl; 00217 } 00218 } 00219 } 00220 } 00221 00223 SmartDataPtr<TofDigiCol> tofDigiCol(eventSvc(),"/Event/Digi/TofDigiCol"); 00224 if (!tofDigiCol) { 00225 log << MSG::FATAL << "EmcRec could not find Tof digi!!" << endreq; 00226 return( StatusCode::FAILURE); 00227 } 00228 00229 // Tof energy calibration 00230 if(fTofEnergy) { 00231 EmcRecTofDigitCalib tofDigitCalib; 00232 tofDigitCalib.Convert(tofDigiCol, tofHitMap); 00233 } 00234 } 00235 #endif 00236 00238 fDigitMap.clear(); 00239 fHitMap.clear(); 00240 fClusterMap.clear(); 00241 fShowerMap.clear(); 00242 00243 // Get EmcCalibConstSvc. 00244 IEmcCalibConstSvc *emcCalibConstSvc; 00245 StatusCode sc = service("EmcCalibConstSvc", emcCalibConstSvc); 00246 if(sc != StatusCode::SUCCESS) { 00247 cout << "EmcRecDigit2Hit Error: Can't get EmcCalibConstSvc." << endl; 00248 } 00249 00250 SmartDataPtr<EmcDigiCol> emcDigiCol(eventSvc(),"/Event/Digi/EmcDigiCol"); 00251 if (!emcDigiCol) { 00252 log << MSG::FATAL << "Could not find EMC digi" << endreq; 00253 return( StatusCode::FAILURE); 00254 } 00255 00256 EmcDigiCol::iterator iter3; 00257 for (iter3=emcDigiCol->begin();iter3!= emcDigiCol->end();iter3++) { 00258 RecEmcID id((*iter3)->identify()); 00259 double adc=RawDataUtil::EmcCharge((*iter3)->getMeasure(), 00260 (*iter3)->getChargeChannel()); 00261 double tdc=RawDataUtil::EmcTime((*iter3)->getTimeChannel()); 00262 00263 //for cable connection correction 00264 unsigned int nnpart=EmcID::barrel_ec(id); 00265 unsigned int nnthe=EmcID::theta_module(id); 00266 unsigned int nnphi=EmcID::phi_module(id); 00267 00268 int index = emcCalibConstSvc->getIndex(nnpart,nnthe,nnphi); 00269 int ixtalNumber=emcCalibConstSvc->getIxtalNumber(index); 00270 00271 if(run>0&&ixtalNumber>0) { 00272 unsigned int npartNew=emcCalibConstSvc->getPartID(ixtalNumber); 00273 unsigned int ntheNew=emcCalibConstSvc->getThetaIndex(ixtalNumber); 00274 unsigned int nphiNew=emcCalibConstSvc->getPhiIndex(ixtalNumber); 00275 id=EmcID::crystal_id(npartNew,ntheNew,nphiNew); 00276 }//------- 00277 00278 RecEmcDigit aDigit; 00279 aDigit.Assign(id,adc,tdc); 00280 fDigitMap[id]=aDigit; 00281 } 00282 00283 if(fOutput>=2) { 00284 RecEmcDigitMap::iterator iDigitMap; 00285 for(iDigitMap=fDigitMap.begin(); 00286 iDigitMap!=fDigitMap.end(); 00287 iDigitMap++){ 00288 cout<<iDigitMap->second; 00289 } 00290 } 00291 // DigitMap ok 00292 00293 // oooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 00295 fDigit2Hit.Convert(fDigitMap,fHitMap); 00296 if(fOutput>=2) { 00297 RecEmcHitMap::iterator iHitMap; 00298 for(iHitMap=fHitMap.begin(); 00299 iHitMap!=fHitMap.end(); 00300 iHitMap++){ 00301 cout<<iHitMap->second; 00302 } 00303 fDigit2Hit.Output(fHitMap); 00304 } 00305 00306 00308 fHit2Cluster.Convert(fHitMap,fClusterMap); 00309 // 00310 if(fOutput>=2) { 00311 RecEmcClusterMap::iterator iClusterMap; 00312 for(iClusterMap=fClusterMap.begin(); 00313 iClusterMap!=fClusterMap.end(); 00314 iClusterMap++){ 00315 cout<<iClusterMap->second; 00316 } 00317 } 00318 00320 fCluster2Shower->Convert(fClusterMap,fShowerMap); 00321 // 00322 if(fOutput>=2) { 00323 RecEmcShowerMap::iterator iShowerMap; 00324 for(iShowerMap=fShowerMap.begin(); 00325 iShowerMap!=fShowerMap.end(); 00326 iShowerMap++) { 00327 cout<<iShowerMap->second; 00328 } 00329 } 00330 00331 #ifndef OnlineMode 00332 if(!(m_rawDataProviderSvc->isOnlineMode())) { 00334 if(fTofEnergy) { 00335 EmcRecTofMatch tofMatch; 00336 tofMatch.Match(fShowerMap, tofHitMap); 00337 } 00338 00339 //nshower=fShowerMap.size(); 00340 //EmcRecFindTofShower tofShower; 00341 //tofShower.FindShower(tofHitMap, fShowerMap); 00342 } 00343 #endif 00344 00346 EmcRecTDS tds; 00347 tds.RegisterToTDS(fHitMap,fClusterMap,fShowerMap); 00348 if(fOutput>=2) { 00349 tds.CheckRegister(); 00350 } 00351 00352 // oooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 00353 #ifndef OnlineMode 00354 if(!(m_rawDataProviderSvc->isOnlineMode())) { 00355 if(fOutput>=1) { 00356 nrun=run; 00357 nrec=event; 00358 00359 //cout<<"cm="<<cm<<"\tmm="<<mm<<"\tGeV="<<GeV<<"\tMeV="<<MeV<<endl; 00360 ndigi=fDigitMap.size(); 00361 nhit=fHitMap.size(); 00362 ncluster=fClusterMap.size(); 00363 RecEmcShowerVec fShowerVec; 00364 RecEmcShowerMap::iterator iShowerMap; 00365 for(iShowerMap=fShowerMap.begin(); 00366 iShowerMap!=fShowerMap.end(); 00367 iShowerMap++) { 00368 fShowerVec.push_back(iShowerMap->second); 00369 } 00370 sort(fShowerVec.begin(), fShowerVec.end(), greater<RecEmcShower>()); 00371 nneu=fShowerMap.size(); 00372 00373 RecEmcShower aShower; 00374 RecEmcShowerVec::iterator iShowerVec; 00375 iShowerVec=fShowerVec.begin(); 00376 npart=-99; 00377 ntheta=-99; 00378 nphi=-99; 00379 if(iShowerVec!=fShowerVec.end()) { 00380 aShower=*iShowerVec; 00381 RecEmcID id=aShower.getShowerId(); 00382 npart=EmcID::barrel_ec(id); 00383 ntheta=EmcID::theta_module(id); 00384 nphi=EmcID::phi_module(id); 00385 } 00386 pp1[0]=aShower.energy()*aShower.x()/aShower.position().mag(); 00387 pp1[1]=aShower.energy()*aShower.y()/aShower.position().mag(); 00388 pp1[2]=aShower.energy()*aShower.z()/aShower.position().mag(); 00389 pp1[3]=aShower.energy(); 00390 theta1=(aShower.position()-posG).theta(); 00391 phi1=(aShower.position()-posG).phi(); 00392 if(phi1<0) { 00393 phi1=twopi+phi1; 00394 } 00395 thtgap1=aShower.ThetaGap(); 00396 phigap1=aShower.PhiGap(); 00397 RecEmcID seed,nseed; 00398 seed=aShower.getShowerId(); 00399 nseed=aShower.NearestSeed(); 00400 eseed=aShower.eSeed(); 00401 e3x3=aShower.e3x3(); 00402 e5x5=aShower.e5x5(); 00403 etof2x1=aShower.getETof2x1(); 00404 etof2x3=aShower.getETof2x3(); 00405 if(aShower.getCluster()) { 00406 cluster2ndMoment=aShower.getCluster()->getSecondMoment(); 00407 } 00408 secondMoment=aShower.secondMoment(); 00409 latMoment=aShower.latMoment(); 00410 a20Moment=aShower.a20Moment(); 00411 a42Moment=aShower.a42Moment(); 00412 if(nseed.is_valid()) { 00413 enseed=fHitMap[nseed].getEnergy(); 00414 } else { 00415 enseed=0; 00416 } 00417 00418 dphi1=phi1-tphi; 00419 if(dphi1<-pi) dphi1=dphi1+twopi; 00420 if(dphi1>pi) dphi1=dphi1-twopi; 00421 00422 if(iShowerVec!=fShowerVec.end()) { 00423 iShowerVec++; 00424 if(iShowerVec!=fShowerVec.end()) { 00425 aShower=*iShowerVec; 00426 pp2[0]=aShower.energy()*aShower.x()/aShower.position().mag(); 00427 pp2[1]=aShower.energy()*aShower.y()/aShower.position().mag(); 00428 pp2[2]=aShower.energy()*aShower.z()/aShower.position().mag(); 00429 pp2[3]=aShower.energy(); 00430 } 00431 } 00432 00433 //if(abs(dphi1)<0.5&&abs(dphi1)>0.1) fDigit2Hit.Output(fHitMap); 00434 00435 //for pi0 00436 if(fShowerVec.size()>=2) { 00437 RecEmcShowerVec::iterator iShowerVec1,iShowerVec2; 00438 iShowerVec1=fShowerVec.begin(); 00439 iShowerVec2=fShowerVec.begin()+1; 00440 double e1=(*iShowerVec1).energy(); 00441 double e2=(*iShowerVec2).energy(); 00442 double angle=(*iShowerVec1).position().angle((*iShowerVec2).position()); 00443 mpi0=sqrt(2*e1*e2*(1-cos(angle))); 00444 } 00445 m_tuple->write(); 00446 } 00447 } 00448 #endif 00449 00450 return StatusCode::SUCCESS; 00451 }
|
|
|
|
00454 { 00455 EmcRecParameter::Kill(); 00456 if(fCluster2Shower) delete fCluster2Shower; 00457 00458 MsgStream log(msgSvc(), name()); 00459 log << MSG::INFO << "in finalize()" << endreq; 00460 return StatusCode::SUCCESS; 00461 }
|
|
|
|
00057 { 00058 00059 MsgStream log(msgSvc(), name()); 00060 log << MSG::INFO << "in initialize()" << endreq; 00061 00062 // Get RawDataProviderSvc 00063 StatusCode sc = service("RawDataProviderSvc", m_rawDataProviderSvc); 00064 if(sc != StatusCode::SUCCESS) { 00065 log << MSG::ERROR << "EmcRec Error: Can't get RawDataProviderSvc." << endmsg; 00066 } 00067 00068 #ifndef OnlineMode 00069 if(!(m_rawDataProviderSvc->isOnlineMode())) { 00070 m_tuple=0; 00071 StatusCode status; 00072 00073 if(fOutput>=1) { 00074 NTuplePtr nt(ntupleSvc(),"FILE301/n1"); 00075 if ( nt ) m_tuple = nt; 00076 else { 00077 m_tuple=ntupleSvc()->book("FILE301/n1",CLID_ColumnWiseTuple,"EmcRec"); 00078 if( m_tuple ) { 00079 status = m_tuple->addItem ("pid",pid); 00080 status = m_tuple->addItem ("tp",tp); 00081 status = m_tuple->addItem ("ttheta",ttheta); 00082 status = m_tuple->addItem ("tphi",tphi); 00083 status = m_tuple->addItem ("nrun",nrun); 00084 status = m_tuple->addItem ("nrec",nrec); 00085 status = m_tuple->addItem ("nneu",nneu); 00086 status = m_tuple->addItem ("ncluster",ncluster); 00087 status = m_tuple->addItem ("npart",npart); 00088 status = m_tuple->addItem ("ntheta",ntheta); 00089 status = m_tuple->addItem ("nphi",nphi); 00090 status = m_tuple->addItem ("ndigi",ndigi); 00091 status = m_tuple->addItem ("nhit",nhit); 00092 status = m_tuple->addItem ("pp1",4,pp1); 00093 status = m_tuple->addItem ("theta1",theta1); 00094 status = m_tuple->addItem ("phi1",phi1); 00095 status = m_tuple->addItem ("dphi1",dphi1); 00096 status = m_tuple->addItem ("eseed",eseed); 00097 status = m_tuple->addItem ("e3x3",e3x3); 00098 status = m_tuple->addItem ("e5x5",e5x5); 00099 status = m_tuple->addItem ("enseed",enseed); 00100 status = m_tuple->addItem ("etof2x1",etof2x1); 00101 status = m_tuple->addItem ("etof2x3",etof2x3); 00102 status = m_tuple->addItem ("cluster2ndMoment",cluster2ndMoment); 00103 status = m_tuple->addItem ("secondMoment",secondMoment); 00104 status = m_tuple->addItem ("latMoment",latMoment); 00105 status = m_tuple->addItem ("a20Moment",a20Moment); 00106 status = m_tuple->addItem ("a42Moment",a42Moment); 00107 status = m_tuple->addItem ("mpi0",mpi0); 00108 status = m_tuple->addItem ("thtgap1",thtgap1); 00109 status = m_tuple->addItem ("phigap1",phigap1); 00110 status = m_tuple->addItem ("pp2",4,pp2); 00111 } 00112 else { // did not manage to book the N tuple.... 00113 log << MSG::ERROR <<"Cannot book N-tuple:" << long(m_tuple) << endmsg; 00114 return StatusCode::FAILURE; 00115 } 00116 } 00117 } 00118 fCluster2Shower = new EmcRecCluster2Shower; 00119 } else { 00120 fCluster2Shower = new EmcRecFastCluster2Shower; 00121 } 00122 #else 00123 fCluster2Shower = new EmcRecFastCluster2Shower; 00124 #endif 00125 00126 // Get EmcRecParameter 00127 EmcRecParameter& Para=EmcRecParameter::GetInstance(); 00128 Para.SetDigiCalib(fDigiCalib); 00129 Para.SetTimeMin(fTimeMin); 00130 Para.SetTimeMax(fTimeMax); 00131 Para.SetPositionMode(fPositionMode); 00132 00133 return StatusCode::SUCCESS; 00134 }
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|