#include <EventAssemblyAlg.h>
Public Member Functions | |
EventAssemblyAlg (const std::string &name, ISvcLocator *pSvcLocator) | |
StatusCode | initialize () |
StatusCode | execute () |
StatusCode | finalize () |
Private Member Functions | |
double | thetaCut (double p, double nSigma, int parId) |
double | phiCut (double p, double nSigma, int parId) |
Private Attributes | |
bool | myMsgFlag |
bool | myActiveCutFlag |
double | myEmcThetaCut |
double | myEmcPhiCut |
double | myEmcThetaNSigmaCut |
double | myEmcPhiNSigmaCut |
int | myParticleId |
int | m_Output |
NTuple::Tuple * | myNTuple1 |
NTuple::Item< double > | myExted |
NTuple::Item< double > | myMatched |
NTuple::Item< double > | myEnergyMatched |
Definition at line 19 of file EventAssemblyAlg.h.
EventAssemblyAlg::EventAssemblyAlg | ( | const std::string & | name, | |
ISvcLocator * | pSvcLocator | |||
) |
Definition at line 23 of file EventAssemblyAlg.cxx.
References m_Output, myActiveCutFlag, myEmcPhiCut, myEmcPhiNSigmaCut, myEmcThetaCut, myEmcThetaNSigmaCut, myMsgFlag, and myParticleId.
00023 :Algorithm(name, pSvcLocator) 00024 { 00025 //Declare the properties 00026 declareProperty("MsgFlag",myMsgFlag=false); 00027 declareProperty("ActiveCutFlag",myActiveCutFlag=true); 00028 declareProperty("EmcThetaMatchCut",myEmcThetaCut=-1); 00029 declareProperty("EmcPhiMatchCut",myEmcPhiCut=-1); 00030 declareProperty("EmcThetaNSigmaCut",myEmcThetaNSigmaCut=5); 00031 declareProperty("EmcPhiNSigmaCut",myEmcPhiNSigmaCut=5); 00032 declareProperty("ParticleId",myParticleId); 00033 declareProperty("Output", m_Output = 0); 00034 }
StatusCode EventAssemblyAlg::execute | ( | ) |
Definition at line 68 of file EventAssemblyAlg.cxx.
References abs, TofID::barrel_ec(), EventModel::EvtRec::Event, EventModel::EvtRec::EvtRecEvent, EventModel::EvtRec::EvtRecTrackCol, Bes_Common::FATAL, genRecEmupikp::i, Bes_Common::INFO, TofID::is_mrpc(), ganga-rec::j, TofID::layer(), m_Output, M_PI, TofID::module(), msgSvc(), myActiveCutFlag, myEmcPhiCut, myEmcPhiNSigmaCut, myEmcThetaCut, myEmcThetaNSigmaCut, myEnergyMatched, myExted, myMatched, myMsgFlag, myNTuple1, myParticleId, TofID::phi_module(), phiCut(), EventModel::Recon::RecEmcShowerCol, EventModel::Recon::RecExtTrackCol, EventModel::Recon::RecMdcDedxCol, EventModel::Recon::RecMdcKalTrackCol, EventModel::Recon::RecMdcTrackCol, EventModel::Recon::RecMucTrackCol, EventModel::Recon::RecTofTrackCol, EvtRecEvent::setTotalCharged(), EvtRecEvent::setTotalNeutral(), EvtRecEvent::setTotalTracks(), thetaCut(), and x.
00069 { 00070 MsgStream log(msgSvc(), name()); 00071 log << MSG::INFO << "in execute()" << endreq; 00072 00073 int numOfChargedTrks = 0; 00074 int numOfKalTrks = 0; 00075 int numOfDedx = 0; 00076 int numOfExt = 0; 00077 int numOfTof = 0; 00078 int numOfShower = 0; 00079 int numOfMatchedShower = 0; 00080 int numOfMucTrks = 0; 00081 00082 // Part 1: Get the event header, print out event and run number 00083 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader"); 00084 if (!eventHeader) 00085 { 00086 log << MSG::FATAL << "Could not find Event Header" << endreq; 00087 return( StatusCode::FAILURE); 00088 } 00089 log << MSG::INFO << ": retrieved event: " << eventHeader->eventNumber() << " run: " << eventHeader->runNumber() << endreq; 00090 00091 //Part 2: Retrieve Rec data 00092 SmartDataPtr<RecMdcTrackCol> recMdcTrackCol(eventSvc(),EventModel::Recon::RecMdcTrackCol); 00093 RecMdcTrackCol::iterator beginOfMdcTrkCol = RecMdcTrackCol::iterator(NULL); 00094 if(!recMdcTrackCol) 00095 { 00096 log<<MSG::INFO<<"Could not find RecMdcTrackCol!"<<endreq; 00097 } 00098 else 00099 { 00100 beginOfMdcTrkCol = recMdcTrackCol->begin(); 00101 numOfChargedTrks = recMdcTrackCol->size(); 00102 } 00103 00104 SmartDataPtr<RecMdcKalTrackCol> mdcKalTrackCol(eventSvc(), EventModel::Recon::RecMdcKalTrackCol); 00105 RecMdcKalTrackCol::iterator beginOfMdcKalTrackCol = RecMdcKalTrackCol::iterator(NULL); 00106 if(!mdcKalTrackCol) 00107 { 00108 log<<MSG::INFO<<"Could not find MdcKalTrackCol!"<<endreq; 00109 } 00110 else 00111 { 00112 numOfKalTrks = mdcKalTrackCol->size(); 00113 beginOfMdcKalTrackCol = mdcKalTrackCol->begin(); 00114 } 00115 00116 00117 SmartDataPtr<RecMdcDedxCol> mdcDedxCol(eventSvc(), EventModel::Recon::RecMdcDedxCol); 00118 RecMdcDedxCol::iterator beginOfDedxCol = RecMdcDedxCol::iterator(NULL); 00119 if(!mdcDedxCol) 00120 { 00121 log<<MSG::INFO<<"Could not find MdcDedxCol!"<<endreq; 00122 } 00123 else { 00124 numOfDedx = mdcDedxCol->size(); 00125 beginOfDedxCol = mdcDedxCol->begin(); 00126 } 00127 00128 00129 SmartDataPtr<RecExtTrackCol> extTrackCol(eventSvc(), EventModel::Recon::RecExtTrackCol); 00130 RecExtTrackCol::iterator beginOfExtTrackCol = RecExtTrackCol::iterator(NULL); 00131 if(!extTrackCol) 00132 { 00133 log<<MSG::INFO<<"Could not find RecExtTrackCol!"<<endreq; 00134 } 00135 else 00136 { 00137 beginOfExtTrackCol = extTrackCol->begin(); 00138 numOfExt = extTrackCol->size(); 00139 } 00140 00141 SmartDataPtr<RecTofTrackCol> TofTrackCol(eventSvc(), EventModel::Recon::RecTofTrackCol); 00142 RecTofTrackCol::iterator beginOfTofTrackCol = RecTofTrackCol::iterator(NULL); 00143 if(!TofTrackCol) 00144 { 00145 log<<MSG::INFO<<"Could not find RecTofTrackCol!"<<endreq; 00146 } 00147 else 00148 { 00149 beginOfTofTrackCol = TofTrackCol->begin(); 00150 numOfTof = TofTrackCol->size(); 00151 } 00152 00153 SmartDataPtr<RecEmcShowerCol> emcRecShowerCol(eventSvc(), EventModel::Recon::RecEmcShowerCol); 00154 RecEmcShowerCol::iterator beginOfEmcTrackCol = RecEmcShowerCol::iterator(NULL); 00155 if(!emcRecShowerCol) 00156 { 00157 log<<MSG::INFO<<"Could not find RecEmcShowerCol!"<<endreq; 00158 } 00159 else 00160 { 00161 beginOfEmcTrackCol = emcRecShowerCol->begin(); 00162 numOfShower = emcRecShowerCol->size(); 00163 } 00164 00165 SmartDataPtr<RecMucTrackCol> mucTrackCol(eventSvc(), EventModel::Recon::RecMucTrackCol); 00166 RecMucTrackCol::iterator beginOfMucTrackCol = RecMucTrackCol::iterator(NULL); 00167 if(!mucTrackCol) 00168 { 00169 log<<MSG::INFO<<"Could not find RecMucTrackCol!"<<endreq; 00170 } 00171 else 00172 { 00173 beginOfMucTrackCol = mucTrackCol->begin(); 00174 numOfMucTrks = mucTrackCol->size(); 00175 } 00176 00177 //Some match flag 00178 00179 int MaxHit = numOfChargedTrks+numOfShower+numOfTof; 00180 00181 vector<bool> dEdxMatched; 00182 vector<bool> kalTrkMatched; 00183 vector<bool> extMatched; 00184 vector<bool> TofMatched; 00185 vector<bool> emcMatched; 00186 vector<bool> mucMatched; 00187 00188 for(int i=0;i<MaxHit;i++) 00189 { 00190 dEdxMatched.push_back(false); 00191 kalTrkMatched.push_back(false); 00192 extMatched.push_back(false); 00193 TofMatched.push_back(false); 00194 emcMatched.push_back(false); 00195 mucMatched.push_back(false); 00196 } 00197 00198 00199 //Part 3: Match tracks 00200 // RecTrackListCol* aNewRecTrkListCol = new RecTrackListCol; 00201 EvtRecTrackCol* aNewEvtRecTrackCol = new EvtRecTrackCol; 00202 if(numOfChargedTrks+numOfShower)//if mdc trackCol and emc trackCol neither exists,there is no match job. 00203 { 00204 for(int i=0;i<numOfChargedTrks;i++)//match charged tracks 00205 { 00206 EvtRecTrack* aNewEvtRecTrack = new EvtRecTrack; 00207 int mdcTrkID = (*(beginOfMdcTrkCol+i))->trackId(); 00208 aNewEvtRecTrack->setTrackId(mdcTrkID); 00209 aNewEvtRecTrack->setMdcTrack(*(beginOfMdcTrkCol+i)); 00210 // aNewEvtRecTrack->setMdcTrkIdx(i); 00211 for(int j=0;j<numOfKalTrks;j++)//match MdcKalTrack 00212 if( !kalTrkMatched[j] && mdcTrkID==(*(beginOfMdcKalTrackCol+j))->trackId() ) 00213 { 00214 aNewEvtRecTrack->setMdcKalTrack(*(beginOfMdcKalTrackCol+j)); 00215 // aNewRecTrkList->setMdcKalTrkIdx(j); 00216 kalTrkMatched[j]=true; 00217 } 00218 00219 for(int j=0;j<numOfDedx;j++)//match dEdx 00220 if( !dEdxMatched[j] && mdcTrkID==(*(beginOfDedxCol+j))->trackId() ) 00221 { 00222 aNewEvtRecTrack->setMdcDedx(*(beginOfDedxCol+j)); 00223 // aNewRecTrkList->setDedxTrkIdx(j); 00224 dEdxMatched[j]=true; 00225 } 00226 00227 for(int j=0;j<numOfExt;j++)//match ExtTrack 00228 if( !extMatched[j] && mdcTrkID==(*(beginOfExtTrackCol+j))->trackId() ) 00229 { 00230 aNewEvtRecTrack->setExtTrack(*(beginOfExtTrackCol+j)); 00231 // aNewRecTrkList->setExtTrkIdx(j); 00232 extMatched[j]=true; 00233 } 00234 for(int j=0;j<numOfTof;j++)//match Tof 00235 if( !TofMatched[j] && mdcTrkID==(*(beginOfTofTrackCol+j))->trackID() ) 00236 { 00237 aNewEvtRecTrack->addTofTrack(*(beginOfTofTrackCol+j)); 00238 // aNewRecTrkList->setTofTrkIdx(j); 00239 TofMatched[j]=true; 00240 } 00241 for(int j=0;j<numOfMucTrks;j++)//match MucTrack 00242 if( !mucMatched[j] && mdcTrkID==(*(beginOfMucTrackCol+j))->trackId() ) 00243 { 00244 aNewEvtRecTrack->setMucTrack(*(beginOfMucTrackCol+j)); 00245 // aNewRecTrkList->setMucTrkIdx(j); 00246 mucMatched[j]=true; 00247 } 00248 00249 if(m_Output==1){ 00250 myExted = -1.; 00251 myMatched = -1.; 00252 myEnergyMatched = -1.; 00253 } 00254 if(aNewEvtRecTrack->isExtTrackValid()&&(aNewEvtRecTrack->extTrack())->emcVolumeNumber()!=-1)//start to match Emc showers 00255 { 00256 if(m_Output==1)myExted = 1.; 00257 //resieve ext information @ Emc 00258 Hep3Vector extPos = (aNewEvtRecTrack->extTrack())->emcPosition(); 00259 double extTheta = extPos.theta(); 00260 double extPhi = extPos.phi(); 00261 00262 double dAngleMin = 9999.; 00263 int indexOfEmcNearest = -1; 00264 00265 for(int j=0;j<numOfShower;j++)//start to find nearest emc shower 00266 { 00267 if(!emcMatched[j]) 00268 { 00269 HepPoint3D emcPot = (*(beginOfEmcTrackCol+j))->position(); 00270 Hep3Vector emcPos(emcPot.x(),emcPot.y(),emcPot.z()); 00271 double dAngle = extPos.angle(emcPos); 00272 if(dAngle<dAngleMin) 00273 { 00274 dAngleMin = dAngle; 00275 indexOfEmcNearest = j; 00276 } 00277 } 00278 } 00279 if(indexOfEmcNearest>=0)//if nearest Emc shower is found 00280 { 00281 HepPoint3D emcPot = (*(beginOfEmcTrackCol+indexOfEmcNearest))->position(); 00282 double theta = emcPot.theta(); 00283 double phi = emcPot.phi(); 00284 double deltaTheta = (extTheta-theta); 00285 double deltaPhi = fmod(extPhi-phi+5*M_PI,2*M_PI)-M_PI; 00286 if(myActiveCutFlag) 00287 { 00288 double tanLamda = (*(beginOfMdcTrkCol+i))->helix(4); 00289 double kappa = (*(beginOfMdcTrkCol+i))->helix(2); 00290 double p = sqrt(1+tanLamda*tanLamda)/fabs(kappa); 00291 double aThetaCut = thetaCut(p,myEmcThetaNSigmaCut,myParticleId); 00292 double aPhiCUt = phiCut(p,myEmcPhiNSigmaCut,myParticleId); 00293 if(fabs(deltaTheta)<aThetaCut && fabs(deltaPhi)<aPhiCUt) 00294 { 00295 aNewEvtRecTrack->setEmcShower(*(beginOfEmcTrackCol+indexOfEmcNearest)); 00296 // aNewRecTrkList->setEmcTrkIdx(indexOfEmcNearest); 00297 emcMatched[indexOfEmcNearest]=true; 00298 numOfMatchedShower++; 00299 if(m_Output==1)myMatched = 1.; 00300 } 00301 } 00302 else if(fabs(deltaTheta)<myEmcThetaCut && fabs(deltaPhi)<myEmcPhiCut)//if the nearest shower is in the fixed match window,it will be matched. 00303 { 00304 aNewEvtRecTrack->setEmcShower(*(beginOfEmcTrackCol+indexOfEmcNearest)); 00305 // aNewRecTrkList->setEmcTrkIdx(indexOfEmcNearest); 00306 emcMatched[indexOfEmcNearest]=true; 00307 numOfMatchedShower++; 00308 if(m_Output==1)myMatched = 1.; 00309 } 00310 } 00311 } 00312 aNewEvtRecTrackCol->push_back(aNewEvtRecTrack); 00313 if(m_Output==1){ 00314 myEnergyMatched = 0.; 00315 myNTuple1->write(); 00316 } 00317 }//end charged track match 00318 00319 int id=0; 00320 for(int i=0;i<numOfShower;i++)//If Emc shower isn't macthed, it will be added to the end of RecTrackListCol 00321 { 00322 if(!emcMatched[i]) 00323 { 00324 EvtRecTrack* aNewEvtRecTrack = new EvtRecTrack; 00325 aNewEvtRecTrack->setTrackId(id+numOfChargedTrks); 00326 aNewEvtRecTrack->setEmcShower(*(beginOfEmcTrackCol+i)); 00327 // aNewRecTrkList->setEmcTrkIdx(i); 00328 00329 //EMC and TOF match 00330 HepPoint3D emcPos((*(beginOfEmcTrackCol+i))->position()); 00331 double emcPhi = emcPos.phi(); 00332 emcPhi = emcPhi < 0 ? emcPhi+CLHEP::twopi : emcPhi; 00333 00334 int module=9999, nphi1=9999, nphi2=9999, nphi1_mrpc=9999; 00335 module = (*(beginOfEmcTrackCol+i))->module(); 00336 if(module==1) { //barrel 00337 nphi1 = (int)(emcPhi*88./CLHEP::twopi); //tof layer1 number 00338 nphi2 = (int)(emcPhi*88./CLHEP::twopi+0.5); //tof layer2 number 00339 nphi2 += 88; 00340 } else { //endcap 00341 nphi1 = (int)(emcPhi*48./CLHEP::twopi+0.5); 00342 00343 00344 //new mrpc detector 00345 double fi_endtof_mrpc = atan2(emcPos.y(),emcPos.x())-3.141592654/2.; 00346 if(fi_endtof_mrpc<0) fi_endtof_mrpc += 2*3.141592654; 00347 nphi1_mrpc=((int)(fi_endtof_mrpc/(3.141593/18)))+1; //We have 36 modules in each endcap (divided again into two layers)! The numbering starts with 1. 00348 00349 if(nphi1_mrpc%2==0){nphi1_mrpc=nphi1_mrpc/2;} //Calculate the correct number as in each endcap we have two layers of modules furnished with numbers 1 to 18. 00350 else {nphi1_mrpc=(nphi1_mrpc+1)/2;} 00351 00352 if(emcPos.z()>0) {(nphi1_mrpc -= 19); nphi1_mrpc=nphi1_mrpc*(-1);} //Consider that for the east endcap the module at x=0;y=max has the number 18 (the numbering is != to west endcap) 00353 // nphi1_mrpc is now the number of the nearest ToF module! 00354 00355 } 00356 00357 int dPhiMin = 9999; 00358 int partid_dPhiMin = 7; 00359 bool mrpcmatch = false; 00360 00361 int indexOfTofNearest = -1; 00362 for(int j=0;j<numOfTof;j++) { //start to find nearest Tof shower 00363 if( !TofMatched[j] ) { 00364 00365 Identifier id((*(beginOfTofTrackCol+j))->tofID()); 00366 bool is_mrpc = TofID::is_mrpc(id); 00367 00368 if(!is_mrpc) 00369 { 00370 int barrel_ec = TofID::barrel_ec(id); 00371 int layer = TofID::layer(id); 00372 int im = TofID::phi_module(id); 00373 im += layer * 88; 00374 00375 if(module!=barrel_ec) continue; 00376 00377 int dPhi; 00378 if(layer == 0) { 00379 dPhi = abs(im - nphi1); 00380 } else { 00381 dPhi = abs(im - nphi2); 00382 } 00383 if(module==1) { 00384 dPhi = dPhi>=44 ? 88-dPhi : dPhi; 00385 } else { 00386 dPhi = dPhi>=24 ? 48-dPhi : dPhi; 00387 } 00388 00389 if(dPhi < dPhiMin) { 00390 dPhiMin = dPhi; 00391 indexOfTofNearest = j; 00392 } 00393 }//close if(!is_mrpc) 00394 else//is_mrpc 00395 { 00396 int barrel_ec = TofID::barrel_ec(id); 00397 int module_mrpc = TofID::module(id); 00398 00399 int dphi; 00400 dphi = abs(module_mrpc-nphi1_mrpc); 00401 00402 if(dphi < dPhiMin) 00403 { 00404 dPhiMin=dphi; 00405 partid_dPhiMin=barrel_ec; 00406 indexOfTofNearest = j; 00407 mrpcmatch = true; 00408 } 00409 }//close else is_mrpc 00410 00411 }//close if( !TofMatched[j] ) 00412 }//close for(int j=0;j<numOfTof;j++) 00413 00414 /* 00415 if(0) 00416 { 00417 cout << "Event Assembly Algorithm" <<endl; 00418 cout << "EMC Shower number = " << i << endl; 00419 cout << "mrpcmatch ? " << mrpcmatch << endl; 00420 cout << "indexOfTofNearest " << indexOfTofNearest <<endl; 00421 cout << "dphiMin " << dPhiMin << endl; 00422 cout << "partid_mrpc " << partid_dPhiMin <<endl; 00423 cout << "partid_emc " << module <<endl <<endl; 00424 } 00425 */ 00426 00427 //match in z-direction 00428 if(indexOfTofNearest>=0) { 00429 HepPoint3D tofPos(0,0.867,(*(beginOfTofTrackCol+indexOfTofNearest))->zrhit()); 00430 double matWindow = 0; 00431 double eShower = (*(beginOfEmcTrackCol+i))->e5x5();//energy of a 5x5 cluster 00432 if(eShower>0) { 00433 matWindow = 0.01277+0.004268/sqrt(eShower); 00434 } 00435 } 00436 00437 00438 if(indexOfTofNearest>=0 && dPhiMin<=1) { //if nearest Tof shower is found 00439 // calculate matching window for barrel 00440 double eShower = (*(beginOfEmcTrackCol+i))->e5x5(); 00441 double matWindow = 0; 00442 if(eShower>0) { 00443 matWindow = 0.01277+0.004268/sqrt(eShower); 00444 } 00445 matWindow = 0.2; //temp code 00446 00447 HepPoint3D tofPos(0,0.867,(*(beginOfTofTrackCol+indexOfTofNearest))->zrhit()); 00448 if(fabs(tofPos.theta()-emcPos.theta())<matWindow || module!=1) { 00449 aNewEvtRecTrack->addTofTrack(*(beginOfTofTrackCol+indexOfTofNearest)); 00450 TofMatched[indexOfTofNearest]=true; 00451 //cout << "EventAssembly : New neutral Toftrack added " << endl; 00452 00453 } 00454 }//close if(indexOfTofNearest>=0 && dPhiMin<=1) 00455 00456 aNewEvtRecTrackCol->push_back(aNewEvtRecTrack); 00457 id++; 00458 }//close if(!) 00459 }//close for(int i=0;i<numOfShower;i++) <---This are the EMC showers 00460 }//end match 00461 00462 //match statistics 00463 if(myMsgFlag) cout<<"Charged tracks: "<<numOfChargedTrks<<" Ext track: "<<numOfExt <<" Shower:"<<numOfShower<<" Matched shower:"<<numOfMatchedShower<<endl; 00464 00465 //Part 4: register RecTrackListCol and EventList in TDS 00466 DataObject *aDataObject1; 00467 eventSvc()->findObject(EventModel::EvtRec::Event, aDataObject1); 00468 if ( aDataObject1 == NULL ) { 00469 EvtRecObject* aRecObject = new EvtRecObject; 00470 StatusCode sc = eventSvc()->registerObject(EventModel::EvtRec::Event, aRecObject); 00471 if ( sc.isFailure() ) { 00472 log << MSG::FATAL << "Could not register EvtRecObject in TDS!" << endreq; 00473 return StatusCode::FAILURE; 00474 } 00475 } 00476 00477 DataObject* aDataObject2; 00478 eventSvc()->findObject(EventModel::EvtRec::EvtRecEvent, aDataObject2); 00479 if ( aDataObject2 == NULL ) { 00480 EvtRecEvent* aRecEvent = new EvtRecEvent; 00481 aRecEvent->setTotalTracks(numOfChargedTrks+numOfShower-numOfMatchedShower); 00482 aRecEvent->setTotalCharged(numOfChargedTrks); 00483 aRecEvent->setTotalNeutral(numOfShower-numOfMatchedShower); 00484 StatusCode sc = eventSvc()->registerObject(EventModel::EvtRec::EvtRecEvent, aRecEvent); 00485 if ( sc.isFailure() ) { 00486 log << MSG::FATAL << "Could not register EvtRecEvent in TDS!" << endreq; 00487 return( StatusCode::FAILURE); 00488 } 00489 } 00490 else { 00491 EvtRecEvent* aRecEvent = dynamic_cast<EvtRecEvent*>(aDataObject2); 00492 aRecEvent->setTotalTracks(numOfChargedTrks+numOfShower-numOfMatchedShower); 00493 aRecEvent->setTotalCharged(numOfChargedTrks); 00494 aRecEvent->setTotalNeutral(numOfShower-numOfMatchedShower); 00495 } 00496 00497 DataObject* aDataObject3; 00498 eventSvc()->findObject(EventModel::EvtRec::EvtRecTrackCol, aDataObject3); 00499 if ( aDataObject3 != NULL ) { 00500 //IDataManagerSvc* dataMgrSvc = dynamic_cast<IDataManagerSvc*> (eventSvc()); 00501 SmartIF<IDataManagerSvc> dataMgrSvc(eventSvc()); 00502 dataMgrSvc->clearSubTree(EventModel::EvtRec::EvtRecTrackCol); 00503 eventSvc()->unregisterObject(EventModel::EvtRec::EvtRecTrackCol); 00504 } 00505 00506 StatusCode sc = eventSvc()->registerObject(EventModel::EvtRec::EvtRecTrackCol, aNewEvtRecTrackCol); 00507 if ( sc.isFailure() ) { 00508 log << MSG::FATAL << "Could not register EvtRecTrackCol in TDS!" << endreq; 00509 return( StatusCode::FAILURE); 00510 } 00511 00512 //Part 5: check RecTrackListCol in TDS 00513 SmartDataPtr<EvtRecTrackCol> evtRecTrackCol(eventSvc(), EventModel::EvtRec::EvtRecTrackCol); 00514 if(!evtRecTrackCol) 00515 { 00516 log<<MSG::FATAL<<"Could not find RecTrackListCol in TDS!"<<endreq; 00517 } 00518 00519 EvtRecTrackCol::iterator itOfRecTrkListCol = evtRecTrackCol->begin(); 00520 if(myMsgFlag) 00521 for(int i=0;itOfRecTrkListCol!=evtRecTrackCol->end();itOfRecTrkListCol++,i++) 00522 { 00523 cout<<"******************** Track "<<i<<": *********************"<<endl; 00524 cout<<"Track ID: "<<(*itOfRecTrkListCol)->trackId()<<endl; 00525 if((*itOfRecTrkListCol)->isMdcTrackValid()) 00526 { 00527 RecMdcTrack* mdcTrk = (*itOfRecTrkListCol)->mdcTrack(); 00528 double kappa = mdcTrk->helix(2); 00529 double tanl = mdcTrk->helix(4); 00530 cout<<"Mdc kappa, tanl: "<<kappa<<", "<<tanl<<endl; 00531 } 00532 if((*itOfRecTrkListCol)->isExtTrackValid()) 00533 { 00534 RecExtTrack* extTrack = (*itOfRecTrkListCol)->extTrack(); 00535 Hep3Vector emcPos = extTrack->emcPosition(); 00536 cout<<"Ext Emc pos:"<<emcPos.x()<<","<<emcPos.y()<<","<<emcPos.z()<<endl; 00537 } 00538 if((*itOfRecTrkListCol)->isEmcShowerValid()) 00539 { 00540 RecEmcShower* emcTrack = (*itOfRecTrkListCol)->emcShower(); 00541 HepPoint3D emcPos = emcTrack->position(); 00542 double x = emcPos.x(); 00543 double y = emcPos.y(); 00544 double z = emcPos.z(); 00545 cout<<"Emc rec pos:"<<x<<","<<y<<","<<z<<endl; 00546 } 00547 } 00548 00549 return StatusCode::SUCCESS; 00550 }
StatusCode EventAssemblyAlg::finalize | ( | ) |
Definition at line 553 of file EventAssemblyAlg.cxx.
References Bes_Common::INFO, and msgSvc().
00554 { 00555 MsgStream log(msgSvc(), name()); 00556 log << MSG::INFO << "in finalize()" << endreq; 00557 00558 return StatusCode::SUCCESS; 00559 }
StatusCode EventAssemblyAlg::initialize | ( | ) |
Definition at line 37 of file EventAssemblyAlg.cxx.
References calibUtil::ERROR, Bes_Common::INFO, m_Output, msgSvc(), myEmcPhiCut, myEmcPhiNSigmaCut, myEmcThetaCut, myEmcThetaNSigmaCut, myEnergyMatched, myExted, myMatched, myNTuple1, and ntupleSvc().
00038 { 00039 MsgStream log(msgSvc(), name()); 00040 log << MSG::INFO << "in initialize()" << endreq; 00041 00042 if(myEmcThetaCut<0) myEmcThetaCut = 0.1745; 00043 if(myEmcPhiCut<0) myEmcPhiCut = 0.2618; 00044 00045 log << MSG::INFO <<"EmcThetaCut,EmcPhiCut: "<<myEmcThetaCut<<", "<<myEmcPhiCut<<endreq; 00046 log << MSG::INFO <<"emcThetaNSigmaCut, emcPhiNSigmaCut: "<<myEmcThetaNSigmaCut<<", "<<myEmcPhiNSigmaCut<<endreq; 00047 00048 if(m_Output == 1){ 00049 NTuplePtr nt1(ntupleSvc(), "FILE998/match"); 00050 if ( nt1 ) myNTuple1 = nt1; 00051 else { 00052 myNTuple1 = ntupleSvc()->book("FILE998/match", CLID_ColumnWiseTuple, "match"); 00053 if ( myNTuple1 ) { 00054 myNTuple1->addItem("exted",myExted); 00055 myNTuple1->addItem("matched",myMatched); 00056 myNTuple1->addItem("Ematched",myEnergyMatched); 00057 } 00058 else { 00059 log << MSG::ERROR << " Cannot book N-tuple:" << long(myNTuple1) << endmsg; 00060 //return StatusCode::FAILURE; 00061 } 00062 } 00063 } 00064 return StatusCode::SUCCESS; 00065 }
double EventAssemblyAlg::phiCut | ( | double | p, | |
double | nSigma, | |||
int | parId | |||
) | [private] |
Definition at line 581 of file EventAssemblyAlg.cxx.
References cut, exp(), and myEmcPhiCut.
Referenced by execute().
00582 { 00583 if(p<0.3) p=0.3; 00584 if(parId<1||parId>3) parId=3; 00585 if(nSigma<=0) nSigma = 5.; 00586 double cut = myEmcPhiCut; 00587 switch(parId) 00588 { 00589 case 1:cut = fabs(exp(-2.187*p-2.945)+0.03236)+(exp(-2.977*p-3.558)+0.005566)*nSigma; break; 00590 case 2:cut = fabs(exp(-2.216*p-2.13)+0.03314)+(exp(-2.436*p-3.392)+0.005049)*nSigma; break; 00591 case 3:cut = fabs(exp(-1.63*p-2.931)+0.03223)+(exp(-3.192*p-2.756)+0.01533)*nSigma; break; 00592 } 00593 00594 return cut; 00595 }
double EventAssemblyAlg::thetaCut | ( | double | p, | |
double | nSigma, | |||
int | parId | |||
) | [private] |
Definition at line 564 of file EventAssemblyAlg.cxx.
References cut, exp(), and myEmcThetaCut.
Referenced by execute().
00565 { 00566 if(p<0.3) p=0.3; 00567 if(parId<1||parId>3) parId=3; 00568 if(nSigma<=0) nSigma = 5.; 00569 double cut = myEmcThetaCut; 00570 switch(parId) 00571 { 00572 case 1:cut = fabs(-0.00159)+(exp(-3.566*p-4.14)+0.006804)*nSigma; break; 00573 case 2:cut = fabs(-0.00145)+(exp(-4.268*p-3.305)+0.01129)*nSigma; break; 00574 case 3:cut = fabs(-0.00161)+(exp(-5.955*p-3.052)+0.01751)*nSigma; break; 00575 } 00576 00577 return cut; 00578 }
int EventAssemblyAlg::m_Output [private] |
Definition at line 40 of file EventAssemblyAlg.h.
Referenced by EventAssemblyAlg(), execute(), and initialize().
bool EventAssemblyAlg::myActiveCutFlag [private] |
double EventAssemblyAlg::myEmcPhiCut [private] |
Definition at line 35 of file EventAssemblyAlg.h.
Referenced by EventAssemblyAlg(), execute(), initialize(), and phiCut().
double EventAssemblyAlg::myEmcPhiNSigmaCut [private] |
Definition at line 38 of file EventAssemblyAlg.h.
Referenced by EventAssemblyAlg(), execute(), and initialize().
double EventAssemblyAlg::myEmcThetaCut [private] |
Definition at line 34 of file EventAssemblyAlg.h.
Referenced by EventAssemblyAlg(), execute(), initialize(), and thetaCut().
double EventAssemblyAlg::myEmcThetaNSigmaCut [private] |
Definition at line 37 of file EventAssemblyAlg.h.
Referenced by EventAssemblyAlg(), execute(), and initialize().
NTuple::Item<double> EventAssemblyAlg::myEnergyMatched [private] |
NTuple::Item<double> EventAssemblyAlg::myExted [private] |
NTuple::Item<double> EventAssemblyAlg::myMatched [private] |
bool EventAssemblyAlg::myMsgFlag [private] |
NTuple::Tuple* EventAssemblyAlg::myNTuple1 [private] |
int EventAssemblyAlg::myParticleId [private] |