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