#include <TrkExtAlg.h>
Public Member Functions | |
| StatusCode | execute () |
| StatusCode | execute () |
| StatusCode | finalize () |
| StatusCode | finalize () |
| StatusCode | initialize () |
| StatusCode | initialize () |
| TrkExtAlg (const std::string &name, ISvcLocator *pSvcLocator) | |
| TrkExtAlg (const std::string &name, ISvcLocator *pSvcLocator) | |
| ~TrkExtAlg () | |
| ~TrkExtAlg () | |
Private Attributes | |
| bool | msgFlag |
| bool | myBFieldOn |
| Ext_track * | myExtTrack |
| Ext_track * | myExtTrack |
| bool | myGeomOptimization |
| string | myInputTrk |
| string | myParticleName |
| bool | myResultFlag |
|
||||||||||||
|
00032 :Algorithm(name, pSvcLocator)
00033 {
00034 myParticleName = "pi";
00035 msgFlag = false;
00036 myResultFlag = false;
00037 myInputTrk= "Kal";
00038
00039 declareProperty("ParticleName",myParticleName);
00040 declareProperty("GeantGeomOptimization",myGeomOptimization=true);
00041 declareProperty("MessageFlag",msgFlag);
00042 declareProperty("ResultMessageFlag",myResultFlag);
00043 declareProperty("BFieldOn",myBFieldOn=true);
00044 declareProperty("InputTrk",myInputTrk);
00045 }
|
|
|
00049 {
00050 if(myExtTrack) delete myExtTrack;
00051 }
|
|
||||||||||||
|
|
|
|
|
|
|
|
|
|
00093 {
00094 MsgStream log(msgSvc(), name());
00095 log << MSG::INFO << "execute()" << endreq;
00096 int eventNumber, runNumber;
00097
00098 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
00099 if (!eventHeader)
00100 {
00101 log << MSG::FATAL << "Could not find Event Header" << endreq;
00102 return( StatusCode::FAILURE);
00103 }
00104 runNumber = eventHeader->runNumber();
00105 eventNumber = eventHeader->eventNumber();
00106 if(msgFlag)
00107 {
00108 cout<<"TrackExt: ******************* Start a event *******************"<<endl;
00109 cout<<"run= "<<runNumber<<"; event= "<<eventNumber<<endl;
00110 }
00111
00112 ExtMdcTrack aExtMdcTrack;
00113 aExtMdcTrack.SetMsgFlag(msgFlag);
00114
00115 bool setTrk=false;
00116
00117 int parID;
00118 if(myParticleName=="e") parID=0;
00119 else if(myParticleName=="mu") parID=1;
00120 else if(myParticleName=="pi") parID=2;
00121 else if(myParticleName=="kaon") parID=3;
00122 else if(myParticleName=="proton"||myParticleName=="anti_proton") parID=4;
00123
00124 if(myInputTrk=="Mdc")
00125 {
00126 SmartDataPtr<RecMdcTrackCol> aMdcTrackCol(eventSvc(),"/Event/Recon/RecMdcTrackCol");
00127 if(!aMdcTrackCol)
00128 {
00129 log << MSG::WARNING << "Can't find RecMdcTrackCol in TDS!" << endreq;
00130 return( StatusCode::SUCCESS);
00131 }
00132 setTrk=aExtMdcTrack.SetMdcRecTrkCol(aMdcTrackCol);
00133
00134 }
00135 else if(myInputTrk=="Kal")
00136 {
00137 SmartDataPtr<RecMdcKalTrackCol> aMdcKalTrackCol(eventSvc(),"/Event/Recon/RecMdcKalTrackCol");
00138 if(!aMdcKalTrackCol)
00139 {
00140 log << MSG::WARNING << "Can't find RecMdcKalTrackCol in TDS!" << endreq;
00141 return( StatusCode::SUCCESS);
00142 }
00143 setTrk=aExtMdcTrack.SetMdcKalTrkCol(aMdcKalTrackCol);
00144 }
00145 else
00146 {
00147 log << MSG::WARNING << "Wrong type of inputTrk:" << myInputTrk << endreq;
00148 return( StatusCode::SUCCESS);
00149 }
00150
00151 RecExtTrackCol *aExtTrackCol = new RecExtTrackCol;
00152 if(setTrk)
00153 {
00154 while(aExtMdcTrack.GetOneGoodTrk())
00155 {
00156 RecExtTrack *aExtTrack = new RecExtTrack;
00157
00158 for(int i=0; i<5; i++)// extrapolate one track with all 5 hypotheses
00159 {
00160 if(aExtMdcTrack.ReadTrk(i))
00161 {
00162 aExtTrack->SetParType(i);
00163
00164 int trackID = aExtMdcTrack.GetTrackID();
00165 aExtTrack->SetTrackId(trackID);
00166 Hep3Vector position = aExtMdcTrack.GetPosition();
00167 Hep3Vector momentum = aExtMdcTrack.GetMomentum();
00168 HepSymMatrix error = aExtMdcTrack.GetErrorMatrix();
00169 double pathInMDC = aExtMdcTrack.GetTrackLength();
00170 double tofInMdc = aExtMdcTrack.GetTrkTof();
00171
00172 if(msgFlag)
00173 {
00174 cout<<"Start From:"<<position.x()<<' '<<position.y()<<' '
00175 <<position.z()<<endl;
00176 cout<<"Start Momentum:"<<momentum.x()<<' '<<momentum.y()<<' '
00177 <<momentum.z()<<endl;
00178 cout<<"Start Error matrix:"<<error<<endl;
00179 cout<<"Path before start:"<< pathInMDC << endl;
00180 }
00181
00182 ExtSteppingAction *extSteppingAction;
00183 extSteppingAction = myExtTrack->GetStepAction();
00184
00185 G4String aParticleName(parName[i]);
00186 double charge = aExtMdcTrack.GetParticleCharge();
00187 if(!aParticleName.contains("proton"))
00188 {
00189 if(charge>0) aParticleName += "+";
00190 else aParticleName += "-";
00191 }
00192 else
00193 {
00194 if(charge>0) aParticleName = "proton";
00195 else aParticleName = "anti_proton";
00196 }
00197
00198 if(msgFlag)
00199 {
00200 cout<<"Charge: "<<charge<<endl;
00201 cout<<"Particle: "<<aParticleName<<endl;
00202 }
00203
00204 if(myExtTrack->Set(position,momentum,error,aParticleName,pathInMDC,tofInMdc))
00205 {
00206 extSteppingAction->SetExtTrackPointer(aExtTrack);
00207 myExtTrack->TrackExtrapotation();
00208 }
00209 }
00210 }
00211
00212 aExtTrack->SetParType(parID);
00213
00214 if(msgFlag) cout<<"will add aExtTrack!"<<endl;
00215 if(aExtTrackCol)
00216 {
00217 if(aExtTrack) aExtTrackCol->add(aExtTrack);
00218 else if(msgFlag) cout<<"No aExtTrack!"<<endl;
00219 }
00220 else
00221 {
00222 if(msgFlag) cout<<"No aExtTrackCol!"<<endl;
00223 }
00224 if(msgFlag) cout<<"add a aExtTrack!"<<endl;
00225
00226 /*
00227 //For Test
00228 if(myFile)
00229 {
00230 (*myFile)<<endPoint.x()<<' '<<endPoint.y()<<' '
00231 <<endPoint.z()<<' '<<endMomentum.x()
00232 <<' '<<endMomentum.y()<<' '<<endMomentum.z()
00233 <<' '<<endErrorMatrix(1,1)<<' '<<endErrorMatrix(2,2)
00234 <<' '<<endErrorMatrix(3,3)<<' '<<endErrorMatrix(4,4)
00235 <<' '<<endErrorMatrix(5,5)<<' '<<endErrorMatrix(6,6)
00236 <<endl;
00237 }
00238 else {
00239 log << MSG::ERROR <<"can't open file" << endreq;
00240
00241 }
00242 */
00243 }//while
00244 }//if
00245
00246 //Register ExtTrackCol to TDS.
00247 /* ReconEvent *aReconEvent = new ReconEvent();
00248 StatusCode sc = eventSvc()->registerObject("/Event/Recon",aReconEvent);
00249 if(sc!=StatusCode::SUCCESS) {
00250 log << MSG::FATAL << "Could not register ReconEvent" <<endreq;
00251 return( StatusCode::FAILURE);
00252 }
00253 */
00254
00255 IDataManagerSvc *dataManSvc = dynamic_cast<IDataManagerSvc*> (eventSvc());
00256
00257 DataObject *extTrackCol;
00258 eventSvc()->findObject("/Event/Recon/RecExtTrackCol",extTrackCol);
00259 if(extTrackCol != NULL) {
00260 dataManSvc->clearSubTree("/Event/Recon/RecExtTrackCol");
00261 eventSvc()->unregisterObject("/Event/Recon/RecExtTrackCol");
00262 }
00263
00264 StatusCode sc = eventSvc()->registerObject("/Event/Recon/RecExtTrackCol", aExtTrackCol);
00265 if(sc!=StatusCode::SUCCESS) {
00266 log << MSG::FATAL << "Could not register RecExtTrackCol in TDS!" << endreq;
00267 return( StatusCode::FAILURE);
00268 }
00269
00270 //Check ExtTrackCol in TDS.
00271 SmartDataPtr<RecExtTrackCol> aExtTrkCol(eventSvc(),"/Event/Recon/RecExtTrackCol");
00272 if (!aExtTrkCol) {
00273 log << MSG::FATAL << "Can't find RecExtTrackCol in TDS!" << endreq;
00274 return( StatusCode::FAILURE);
00275 }
00276
00277 RecExtTrackCol::iterator iterOfExtTrk;
00278 int j=1;
00279 for(iterOfExtTrk = aExtTrkCol->begin();iterOfExtTrk!=aExtTrkCol->end();iterOfExtTrk++)
00280 {
00281 if(myResultFlag)
00282 {
00283 for(int i=0; i<5; i++)
00284 {
00285 //TOF information.
00286 cout<<"##########track"<<j<<": "<<"("<<i<<")"<<endl;
00287 cout<<"******TOF1:******"<<endl;
00288 cout<<"VolumeName: "<<(*iterOfExtTrk)->tof1VolumeName(i)<<"\t"
00289 <<"VolumeNumber: "<<(*iterOfExtTrk)->tof1VolumeNumber(i)<<"\t"<<endl
00290 <<"Position: "<<(*iterOfExtTrk)->tof1Position(i)<<"\t"
00291 <<"Momentum: "<<(*iterOfExtTrk)->tof1Momentum(i)<<"\t"<<endl
00292 <<"Error matrix: "<<(*iterOfExtTrk)->tof1ErrorMatrix(i)
00293 <<"Error z: "<<(*iterOfExtTrk)->tof1PosSigmaAlongZ(i)<<"\t"
00294 <<"Error Tz: "<<(*iterOfExtTrk)->tof1PosSigmaAlongT(i)<<"\t"
00295 <<"Error x: "<<(*iterOfExtTrk)->tof1PosSigmaAlongX(i)<<"\t"
00296 <<"Error y: "<<(*iterOfExtTrk)->tof1PosSigmaAlongY(i)<<endl
00297 <<"Tof: "<<(*iterOfExtTrk)->tof1(i)<<"\t"
00298 <<"PathOF: "<<(*iterOfExtTrk)->tof1Path(i)
00299 <<endl;
00300 cout<<"******TOF2:******"<<endl;
00301 cout<<"VolumeName: "<<(*iterOfExtTrk)->tof2VolumeName(i)<<"\t"
00302 <<"VolumeNumber: "<<(*iterOfExtTrk)->tof2VolumeNumber(i)<<"\t"<<endl
00303 <<"Position: "<<(*iterOfExtTrk)->tof2Position(i)<<"\t"
00304 <<"Momentum: "<<(*iterOfExtTrk)->tof2Momentum(i)<<"\t"<<endl
00305 <<"Error matrix: "<<(*iterOfExtTrk)->tof2ErrorMatrix(i)
00306 <<"Error z: "<<(*iterOfExtTrk)->tof2PosSigmaAlongZ(i)<<"\t"
00307 <<"Error Tz: "<<(*iterOfExtTrk)->tof2PosSigmaAlongT(i)<<"\t"
00308 <<"Error x: "<<(*iterOfExtTrk)->tof2PosSigmaAlongX(i)<<"\t"
00309 <<"Error y: "<<(*iterOfExtTrk)->tof2PosSigmaAlongY(i)<<endl
00310 <<"Tof: "<<(*iterOfExtTrk)->tof2(i)<<"\t"
00311 <<"PathOF: "<<(*iterOfExtTrk)->tof2Path(i)
00312 <<endl;
00313
00314 //EMC information.
00315 cout<<"******EMC:******"<<endl
00316 <<"VolumeName: "<<(*iterOfExtTrk)->emcVolumeName(i)<<"\t"
00317 <<"VolumeNumber: "<<(*iterOfExtTrk)->emcVolumeNumber(i)<<"\t"<<endl
00318 <<"Position: "<<(*iterOfExtTrk)->emcPosition(i)<<"\t"
00319 <<"Momentum: "<<(*iterOfExtTrk)->emcMomentum(i)<<"\t"<<endl
00320 <<"Error matrix: "<<(*iterOfExtTrk)->emcErrorMatrix(i)
00321 <<"Error theta: "<<(*iterOfExtTrk)->emcPosSigmaAlongTheta(i)<<"\t"
00322 <<"Error phi: "<<(*iterOfExtTrk)->emcPosSigmaAlongPhi(i)<<"\t"
00323 <<"EMC path: "<<(*iterOfExtTrk)->emcPath(i)
00324 <<endl;
00325
00326 //MUC information
00327 cout<<"******MUC:******"<<endl
00328 <<"VolumeName: "<<(*iterOfExtTrk)->mucVolumeName(i)<<"\t"
00329 <<"VolumeNumber: "<<(*iterOfExtTrk)->mucVolumeNumber(i)<<endl
00330 <<"Position: "<<(*iterOfExtTrk)->mucPosition(i)<<"\t"
00331 <<"Momentum: "<<(*iterOfExtTrk)->mucMomentum(i)<<"\t"<<endl
00332 <<"Error matrix: "<<(*iterOfExtTrk)->mucErrorMatrix(i)
00333 <<"Error z: "<<(*iterOfExtTrk)->mucPosSigmaAlongZ(i)<<"\t"
00334 <<"Error Tz: "<<(*iterOfExtTrk)->mucPosSigmaAlongT(i)<<"\t"
00335 <<"Error x: "<<(*iterOfExtTrk)->mucPosSigmaAlongX(i)<<"\t"
00336 <<"Error y: "<<(*iterOfExtTrk)->mucPosSigmaAlongY(i)
00337 <<endl;
00338
00339 //Muc Ext hits information
00340 ExtMucHitVec aExtMucHitVec = (*iterOfExtTrk)->GetExtMucHitVec(i);
00341 int numOfMucHits = aExtMucHitVec.size();
00342 cout<<"******MUC hits:"<<numOfMucHits<<"\t"<<i<<"******"<<endl;
00343 for(int j=0;j<numOfMucHits;j++)
00344 {
00345 cout<<"###Muc Hit "<<j<<":###"<<endl
00346 <<"VolumeName: "<<aExtMucHitVec[j].GetVolumeName()<<"\t"
00347 <<"VolumeNumber: "<<aExtMucHitVec[j].GetVolumeNumber()<<"\t"<<endl
00348 <<"Position: "<<aExtMucHitVec[j].GetPosition()<<"\t"
00349 <<"Momentum: "<<aExtMucHitVec[j].GetMomentum()<<"\t"<<endl
00350 <<"Error z: "<<aExtMucHitVec[j].GetPosSigmaAlongZ()<<"\t"
00351 <<"Error Tz: "<<aExtMucHitVec[j].GetPosSigmaAlongT()<<"\t"
00352 <<"Error x: "<<aExtMucHitVec[j].GetPosSigmaAlongX()<<"\t"
00353 <<"Error y: "<<aExtMucHitVec[j].GetPosSigmaAlongY()<<"\t"
00354 <<endl;
00355 }
00356
00357 }
00358 }
00359 j++;
00360
00361 } // loop ExtTrkCol
00362
00363 if(msgFlag) cout<<"****************** End a event! ****************"<<endl<<endl;
00364
00365 /*
00366 //--------- For Ext Test ----------------
00367 // Retrieve mc truth
00368 SmartDataPtr<McParticleCol> mcParticleCol(eventSvc(),"/Event/MC/McParticleCol");
00369 if (!mcParticleCol) {
00370 log << MSG::FATAL << "Could not find McParticle" << endreq;
00371 return( StatusCode::FAILURE);
00372 }
00373
00374 int numOfTrack = mcParticleCol->size();
00375 if(msgFlag) cout<< "numOfMcTrack: " << numOfTrack << endl;
00376 if(numOfTrack!=2) return StatusCode::SUCCESS;
00377
00378 // Retrieve Emc Mc Hits
00379 SmartDataPtr<EmcMcHitCol> emcMcHitCol(eventSvc(),"/Event/MC/EmcMcHitCol");
00380 if (!emcMcHitCol) {
00381 log << MSG::FATAL << "Could not find EMC truth" << endreq;
00382 return( StatusCode::FAILURE);
00383 }
00384
00385 McParticleCol::iterator iter_mc = mcParticleCol->begin();
00386 EmcMcHitCol::iterator iterEmcBegin = emcMcHitCol->begin();
00387 int numOfEmcHits = emcMcHitCol->size();
00388
00389 // ExtTrackCol *aExtTrackCol = new ExtTrackCol;
00390
00391 for (int i = 1;iter_mc != mcParticleCol->end(); iter_mc++, i++) {
00392 bool flag = (*iter_mc)->primaryParticle();
00393 if(!flag) continue;
00394 int particleId = (*iter_mc)->particleProperty();
00395 double charge = particleId/abs(particleId);
00396 unsigned int sFlag = (*iter_mc)->statusFlags();
00397 int trackIdx = (*iter_mc)->getTrackIndex();
00398 HepPoint3D iPos = (*iter_mc)->initialPosition();
00399 HepPoint3D fPos = (*iter_mc)->finalPosition();
00400 HepLorentzVector iFMomentum = (*iter_mc)->initialFourMomentum();
00401 HepLorentzVector fFMomentum = (*iter_mc)->finalFourMomentum();
00402
00403
00404 bool emcHitFlag = false;
00405 double thetaOfEmcHit = 0;
00406 double phiOfEmcHit = 0;
00407 Hep3Vector posOfEmcHit(0,0,0);
00408 //Get Emc Truth
00409 for(int j=0;j<numOfEmcHits;j++)
00410 {
00411 if(trackIdx==(*(iterEmcBegin+j))->getTrackIndex())
00412 {
00413 emcHitFlag = true;
00414 double x = (*(iterEmcBegin+j))->getPositionX();
00415 double y = (*(iterEmcBegin+j))->getPositionY();
00416 double z = (*(iterEmcBegin+j))->getPositionZ();
00417 Hep3Vector vec(x,y,z);
00418 posOfEmcHit = vec;
00419 thetaOfEmcHit = posOfEmcHit.theta();
00420 phiOfEmcHit = posOfEmcHit.phi();
00421 break;
00422 }
00423 }
00424
00425 ExtSteppingAction *extSteppingAction;
00426 extSteppingAction = myExtTrack->GetStepAction();
00427
00428 G4String aParticleName(myParticleName);
00429 if(!aParticleName.contains("proton"))
00430 {
00431 if(charge>0) aParticleName += "+";
00432 else aParticleName += "-";
00433 }
00434 else
00435 {
00436 if(charge>0) aParticleName = "proton";
00437 else aParticleName = "anti_proton";
00438 }
00439
00440 ExtTrack *aExtTrack = new ExtTrack;
00441 aExtTrack->SetTrackID(trackIdx);
00442
00443 Hep3Vector position(iPos);
00444 Hep3Vector momentum(iFMomentum.x(),iFMomentum.y(),iFMomentum.z());
00445 HepSymMatrix error(6,0);
00446
00447 if(myExtTrack->Set(position,momentum,error,aParticleName,0))
00448 {
00449 extSteppingAction->SetExtTrackPointer(aExtTrack);
00450 myExtTrack->TrackExtrapotation();
00451 }
00452
00453 Hep3Vector extEmcPos = aExtTrack->GetEmcPosition();
00454 double volumeNum = aExtTrack->GetEmcVolumeNumber();
00455 double thetaExt = extEmcPos.theta();
00456 double phiExt = extEmcPos.phi();
00457
00458 if(myResultFlag)
00459 {
00460 cout<< "*******Mc Track " << i <<" :******"<<endl;
00461 cout<< "Parimary particle :" << flag <<" particleId = " << (*iter_mc)->particleProperty() << endl;
00462 cout<< "Track Index: " << trackIdx << endl;
00463 cout<< "initialPosition: "<< iPos.x() << "," << iPos.y() << "," << iPos.z() << endl;
00464 cout<< "initialFourMomentum: " <<iFMomentum.x()<<","<<iFMomentum.y()<<","<<iFMomentum.z() << endl;
00465 cout<<"Emc Truth:"<<emcHitFlag<<"!"<<posOfEmcHit<<endl;
00466
00467 cout<<"Ext Emc: "<<aExtTrack->GetEmcVolumeName()<<aExtTrack->GetEmcPosition()<<endl;
00468 }
00469
00470 myCharge = charge;
00471 myEmcHitFlag = (emcHitFlag)? 1.:-1.;
00472 myEmcHitTheta = thetaOfEmcHit;
00473 myEmcHitPhi = phiOfEmcHit;
00474 myEmcVolNum = volumeNum;
00475 myEmcExtTheta = thetaExt;
00476 myEmcExtPhi = phiExt;
00477 myDTheta = myEmcHitTheta-myEmcExtTheta;
00478 myDPhi = myEmcHitPhi-myEmcExtPhi;
00479 while(myDTheta<-1*M_PI) myDTheta+=2.0*M_PI;
00480 while(myDPhi<-1*M_PI) myDPhi+=2.0*M_PI;
00481 while(myDPhi>M_PI) myDPhi-=2.0*M_PI;
00482 while(myDTheta>M_PI) myDTheta-=2.0*M_PI;
00483 myNtuple->write();
00484
00485 if(aExtTrack) delete aExtTrack;
00486 }
00487
00488 //--------- end Ext Test ----------------
00489 */
00490
00491 return StatusCode::SUCCESS;
00492 }
|
|
|
|
|
|
00496 {
00497 MsgStream log(msgSvc(), name());
00498 log << MSG::INFO << "finalize()" << endreq;
00499
00500 // delete myExtTrack;
00501 // myFile->close();
00502
00503 return StatusCode::SUCCESS;
00504 }
|
|
|
|
|
|
00055 {
00056 MsgStream log(msgSvc(), name());
00057 log << MSG::INFO << "initialize()" << endreq;
00058
00059 myExtTrack = new Ext_track(msgFlag,myBFieldOn,myGeomOptimization);
00060 myExtTrack->Initialization(msgFlag,myBFieldOn,myGeomOptimization);
00061 // myFile = new ofstream("ExtData.txt");
00062
00063 /*
00064 //--------- For Ext Test ----------------
00065 NTuplePtr nt(ntupleSvc(),"FILE501/ext");
00066 if ( nt ) myNtuple = nt;
00067 else {
00068 myNtuple=ntupleSvc()->book("FILE501/ext",CLID_ColumnWiseTuple,"TrkExt");
00069 if(myNtuple) {
00070 myNtuple->addItem("charge",myCharge);
00071 myNtuple->addItem("emcHitFlag",myEmcHitFlag);
00072 myNtuple->addItem("emcHitTheta",myEmcHitTheta);
00073 myNtuple->addItem("emcHitPhi",myEmcHitPhi);
00074 myNtuple->addItem("emcVolNum",myEmcVolNum);
00075 myNtuple->addItem("emcExtTheta",myEmcExtTheta);
00076 myNtuple->addItem("emcExtPhi",myEmcExtPhi);
00077 myNtuple->addItem("dTheta",myDTheta);
00078 myNtuple->addItem("dPhi",myDPhi);
00079 }
00080 else { // did not manage to book the N tuple....
00081 log << MSG::ERROR <<"Cannot book N-tuple:" << long(myNtuple) << endmsg;
00082 return StatusCode::FAILURE;
00083 }
00084 }
00085 //-------- end Ext Test -----------------
00086 */
00087 return StatusCode::SUCCESS;
00088 }
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
1.3.9.1