BOSS 7.0.6
BESIII Offline Software System
Loading...
Searching...
No Matches
TrkExtAlg.cxx
Go to the documentation of this file.
1//
2// File: TrkExtAlg.cxx
3// Author: Wang Liangliang
4// Date: 2005.4.4
5//
6#include "GaudiKernel/MsgStream.h"
7#include "GaudiKernel/AlgFactory.h"
8#include "GaudiKernel/ISvcLocator.h"
9#include "GaudiKernel/SmartDataPtr.h"
10#include "GaudiKernel/IDataProviderSvc.h"
11#include "GaudiKernel/IDataManagerSvc.h"
12#include "GaudiKernel/PropertyMgr.h"
13#include "EventModel/EventHeader.h"
14
15#include "CLHEP/Matrix/SymMatrix.h"
16#include "CLHEP/Vector/ThreeVector.h"
17
18#include "ExtEvent/RecExtTrack.h"
19#include "TrkExtAlg/TrkExtAlg.h"
20#include "TrkExtAlg/ExtSteppingAction.h"
21/////////////////////////////////change
22//#include "TrkExtAlg/Helix.h"
23#include "TrackUtil/Helix.h"
24#include "DetVerSvc/IDetVerSvc.h"
25//////////////////////////////////
26#include "TrkExtAlg/ExtMdcTrack.h"
27
28using namespace Event;
29
30string parName[5]={"e","mu","pi","kaon","proton"};
31
32//Constructor
33TrkExtAlg::TrkExtAlg(const string& name, ISvcLocator* pSvcLocator):Algorithm(name, pSvcLocator)
34{
35 myParticleName = "pi";
36 msgFlag = false;
37 myResultFlag = false;
38 myInputTrk= "Kal";
39
40 declareProperty("ParticleName",myParticleName);
41 declareProperty("GeantGeomOptimization",myGeomOptimization=true);
42 declareProperty("MessageFlag",msgFlag);
43 declareProperty("ResultMessageFlag",myResultFlag);
44 declareProperty("BFieldOn",myBFieldOn=true);
45 declareProperty("InputTrk",myInputTrk);
46 declareProperty("UseMucKalFilter",myUseMucKal=true);
47 declareProperty("MucWindow",myMucWindow=6);
48 declareProperty("SetSeed",m_setSeed=false);
49}
50
51//Destructor
52TrkExtAlg::~TrkExtAlg()
53{
54 if(myExtTrack) delete myExtTrack;
55}
56
57//////////////////////
58StatusCode TrkExtAlg::initialize()
59{
60 MsgStream log(msgSvc(), name());
61 log << MSG::INFO << "initialize()" << endreq;
62
63 IDetVerSvc* detVerSvc;
64 StatusCode sc_det = service("DetVerSvc", detVerSvc);
65 if( sc_det.isFailure() ) {
66 log << MSG::FATAL << "can't retrieve DetVerSvc instance" << endreq;
67 return sc_det;
68 }
69
70 m_detVer = detVerSvc->phase();
71 log << MSG::INFO << "** ~~~~~~ ** : retrieved DetectorStage = " << m_detVer << endreq;
72
73 myExtTrack = new Ext_track(msgFlag,myBFieldOn,myGeomOptimization,m_detVer,myUseMucKal,myMucWindow);
74 myExtTrack->Initialization(msgFlag,myBFieldOn,myGeomOptimization,myUseMucKal,myMucWindow);
75 // myFile = new ofstream("ExtData.txt");
76
77 /*
78 //--------- For Ext Test ----------------
79 NTuplePtr nt(ntupleSvc(),"FILE501/ext");
80 if ( nt ) myNtuple = nt;
81 else {
82 myNtuple=ntupleSvc()->book("FILE501/ext",CLID_ColumnWiseTuple,"TrkExt");
83 if(myNtuple) {
84 myNtuple->addItem("charge",myCharge);
85 myNtuple->addItem("emcHitFlag",myEmcHitFlag);
86 myNtuple->addItem("emcHitTheta",myEmcHitTheta);
87 myNtuple->addItem("emcHitPhi",myEmcHitPhi);
88 myNtuple->addItem("emcVolNum",myEmcVolNum);
89 myNtuple->addItem("emcExtTheta",myEmcExtTheta);
90 myNtuple->addItem("emcExtPhi",myEmcExtPhi);
91 myNtuple->addItem("dTheta",myDTheta);
92 myNtuple->addItem("dPhi",myDPhi);
93 }
94 else { // did not manage to book the N tuple....
95 log << MSG::ERROR <<"Cannot book N-tuple:" << long(myNtuple) << endmsg;
96 return StatusCode::FAILURE;
97 }
98 }
99 //-------- end Ext Test -----------------
100 */
101 return StatusCode::SUCCESS;
102}
103
104
105//////////////////////
106StatusCode TrkExtAlg::execute()
107{
108 //cout<<"a new event "<<endl;
109 MsgStream log(msgSvc(), name());
110 log << MSG::INFO << "execute()" << endreq;
111 int eventNumber, runNumber;
112
113 if(m_setSeed==true) CLHEP::HepRandom::setTheSeed(9000);
114 //cout<<"seed "<<CLHEP::HepRandom::getTheSeed();
115
116 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
117 if (!eventHeader)
118 {
119 log << MSG::FATAL << "Could not find Event Header" << endreq;
120 return( StatusCode::FAILURE);
121 }
122 runNumber = eventHeader->runNumber();
123 eventNumber = eventHeader->eventNumber();
124 // if(eventNumber!=5&&eventNumber!=8&&eventNumber!=142) return( StatusCode::SUCCESS);//1900736 //839313
125 if(msgFlag)
126 {
127 cout<<"TrackExt: ******************* Start a event *******************"<<endl;
128 cout<<"run= "<<runNumber<<"; event= "<<eventNumber<<endl;
129 }
130 // cout<<"run= "<<runNumber<<"; event= "<<eventNumber<<endl;
131 /////
132 SmartDataPtr<MucDigiCol> mucDigiCol(eventSvc(),"/Event/Digi/MucDigiCol");
133 if(!mucDigiCol) {
134 log << MSG::FATAL << "Could not find MUC digi" << endreq;
135 return( StatusCode::SUCCESS);
136 }
137 //cout<<"digi col is "<<mucDigiCol->size()<<endl;
138 ////
139 ExtMdcTrack aExtMdcTrack;
140 aExtMdcTrack.SetMsgFlag(msgFlag);
141
142 bool setTrk=false;
143
144 int parID;
145 if(myParticleName=="e") parID=0;
146 else if(myParticleName=="mu") parID=1;
147 else if(myParticleName=="pi") parID=2;
148 else if(myParticleName=="kaon") parID=3;
149 else if(myParticleName=="proton"||myParticleName=="anti_proton") parID=4;
150
151 if(myInputTrk=="Mdc")
152 {
153 SmartDataPtr<RecMdcTrackCol> aMdcTrackCol(eventSvc(),"/Event/Recon/RecMdcTrackCol");
154 if(!aMdcTrackCol)
155 {
156 log << MSG::WARNING << "Can't find RecMdcTrackCol in TDS!" << endreq;
157 return( StatusCode::SUCCESS);
158 }
159 setTrk=aExtMdcTrack.SetMdcRecTrkCol(aMdcTrackCol);
160
161 }
162 else if(myInputTrk=="Kal")
163 {
164 SmartDataPtr<RecMdcKalTrackCol> aMdcKalTrackCol(eventSvc(),"/Event/Recon/RecMdcKalTrackCol");
165 if(!aMdcKalTrackCol)
166 {
167 log << MSG::WARNING << "Can't find RecMdcKalTrackCol in TDS!" << endreq;
168 return( StatusCode::SUCCESS);
169 }
170 setTrk=aExtMdcTrack.SetMdcKalTrkCol(aMdcKalTrackCol);
171 }
172 else
173 {
174 log << MSG::WARNING << "Wrong type of inputTrk:" << myInputTrk << endreq;
175 return( StatusCode::SUCCESS);
176 }
177
178 RecExtTrackCol *aExtTrackCol = new RecExtTrackCol;
179 if(setTrk)
180 {
181 while(aExtMdcTrack.GetOneGoodTrk())
182 {
183 //cout<<"start a track........................"<<endl;
184 RecExtTrack *aExtTrack = new RecExtTrack;
185
186 for(int i=0; i<5; i++)// extrapolate one track with all 5 hypotheses
187 {
188 if(aExtMdcTrack.ReadTrk(i))
189 {
190 aExtTrack->SetParType(i);
191
192 int trackID = aExtMdcTrack.GetTrackID();
193 aExtTrack->SetTrackId(trackID);
194 Hep3Vector position = aExtMdcTrack.GetPosition();
195 Hep3Vector momentum = aExtMdcTrack.GetMomentum();
196 HepSymMatrix error = aExtMdcTrack.GetErrorMatrix();
197 double pathInMDC = aExtMdcTrack.GetTrackLength();
198 double tofInMdc = aExtMdcTrack.GetTrkTof();
199
200 if(msgFlag)
201 {
202 cout<<"Start From:"<<position.x()<<' '<<position.y()<<' '
203 <<position.z()<<endl;
204 cout<<"Start Momentum:"<<momentum.x()<<' '<<momentum.y()<<' '<<momentum.z()<<endl;
205 cout<<"Start Error matrix:"<<error<<endl;
206 cout<<"Path before start:"<< pathInMDC << endl;
207 }
208
209 G4String aParticleName(parName[i]);
210 double charge = aExtMdcTrack.GetParticleCharge();
211 if(!aParticleName.contains("proton"))
212 {
213 if(charge>0) aParticleName += "+";
214 else aParticleName += "-";
215 }
216 else
217 {
218 if(charge>0) aParticleName = "proton";
219 else aParticleName = "anti_proton";
220 }
221
222 if(msgFlag)
223 {
224 cout<<"Charge: "<<charge<<endl;
225 cout<<"Particle: "<<aParticleName<<endl;
226 }
227
228 ExtSteppingAction *extSteppingAction;
229 extSteppingAction = myExtTrack->GetStepAction();
230 extSteppingAction->Set_which_tof_version(m_detVer);//Set the ToFversionnumber. Required for the MRPC
231
232 extSteppingAction->Reset();
233 extSteppingAction->SetMucDigiColPointer(mucDigiCol);
234 extSteppingAction->SetExtTrackPointer(aExtTrack);
235 bool m_trackstatus = false;
236 int trk_startpart = 0;
237 while(!m_trackstatus)
238 {
239
240 trk_startpart++;//just for protection
241 if(trk_startpart>20)
242 {cout<<"-------has modified more than 20 times---------"<<endl;break;}
243 if(myExtTrack->Set(position,momentum,error,aParticleName,pathInMDC,tofInMdc))
244 {
245 myExtTrack->TrackExtrapotation();
246 extSteppingAction->InfmodMuc(position,momentum,error);
247 m_trackstatus = extSteppingAction->TrackStop();
248 }
249 else
250 m_trackstatus = true;
251 }
252 }
253
254 }
255
256 aExtTrack->SetParType(parID);
257
258 if(msgFlag) cout<<"will add aExtTrack!"<<endl;
259 if(aExtTrackCol)
260 {
261 if(aExtTrack) aExtTrackCol->add(aExtTrack);
262 else if(msgFlag) cout<<"No aExtTrack!"<<endl;
263 }
264 else
265 {
266 if(msgFlag) cout<<"No aExtTrackCol!"<<endl;
267 }
268 if(msgFlag) cout<<"add a aExtTrack!"<<endl;
269
270 /*
271 //For Test
272 if(myFile)
273 {
274 (*myFile)<<endPoint.x()<<' '<<endPoint.y()<<' '
275 <<endPoint.z()<<' '<<endMomentum.x()
276 <<' '<<endMomentum.y()<<' '<<endMomentum.z()
277 <<' '<<endErrorMatrix(1,1)<<' '<<endErrorMatrix(2,2)
278 <<' '<<endErrorMatrix(3,3)<<' '<<endErrorMatrix(4,4)
279 <<' '<<endErrorMatrix(5,5)<<' '<<endErrorMatrix(6,6)
280 <<endl;
281 }
282 else {
283 log << MSG::ERROR <<"can't open file" << endreq;
284
285 }
286 */
287 }//while
288 }//if
289
290 //Register ExtTrackCol to TDS.
291 /* ReconEvent *aReconEvent = new ReconEvent();
292 StatusCode sc = eventSvc()->registerObject("/Event/Recon",aReconEvent);
293 if(sc!=StatusCode::SUCCESS) {
294 log << MSG::FATAL << "Could not register ReconEvent" <<endreq;
295 return( StatusCode::FAILURE);
296 }
297 */
298
299 //cout<<"will write in to service "<<endl;
300 SmartIF<IDataManagerSvc> dataManSvc(eventSvc());
301 //IDataManagerSvc *dataManSvc = dynamic_cast<IDataManagerSvc*> (eventSvc());
302
303 DataObject *extTrackCol;
304 eventSvc()->findObject("/Event/Recon/RecExtTrackCol",extTrackCol);
305 if(extTrackCol != NULL) {
306 dataManSvc->clearSubTree("/Event/Recon/RecExtTrackCol");
307 eventSvc()->unregisterObject("/Event/Recon/RecExtTrackCol");
308 }
309
310 // cout<<"write register object "<<endl;
311
312 StatusCode sc = eventSvc()->registerObject("/Event/Recon/RecExtTrackCol", aExtTrackCol);
313 if(sc!=StatusCode::SUCCESS) {
314 log << MSG::FATAL << "Could not register RecExtTrackCol in TDS!" << endreq;
315 return( StatusCode::FAILURE);
316 }
317
318 //Check ExtTrackCol in TDS.
319 //cout<<"Check ExtTrackCol in TDS."<<endl;
320 SmartDataPtr<RecExtTrackCol> aExtTrkCol(eventSvc(),"/Event/Recon/RecExtTrackCol");
321 if (!aExtTrkCol) {
322 log << MSG::FATAL << "Can't find RecExtTrackCol in TDS!" << endreq;
323 return( StatusCode::FAILURE);
324 }
325
326 RecExtTrackCol::iterator iterOfExtTrk;
327 int j=1;
328 //cout<<"aExtTrkCol->size is "<<aExtTrkCol->size()<<" --------------"<<endl;
329 for(iterOfExtTrk = aExtTrkCol->begin();iterOfExtTrk!=aExtTrkCol->end();iterOfExtTrk++)
330 {
331 if(myResultFlag)
332 {
333 for(int i=0; i<5; i++)
334 {
335 //TOF information.
336 cout<<"##########track"<<j<<": "<<"("<<i<<")"<<endl;
337 cout<<"******TOF1:******"<<endl;
338 cout<<"VolumeName: "<<(*iterOfExtTrk)->tof1VolumeName(i)<<"\t"
339 <<"VolumeNumber: "<<(*iterOfExtTrk)->tof1VolumeNumber(i)<<"\t"<<endl
340 <<"Position: "<<(*iterOfExtTrk)->tof1Position(i)<<"\t"
341 <<"Momentum: "<<(*iterOfExtTrk)->tof1Momentum(i)<<"\t"<<endl
342 <<"Error matrix: "<<(*iterOfExtTrk)->tof1ErrorMatrix(i)
343 <<"Error z: "<<(*iterOfExtTrk)->tof1PosSigmaAlongZ(i)<<"\t"
344 <<"Error Tz: "<<(*iterOfExtTrk)->tof1PosSigmaAlongT(i)<<"\t"
345 <<"Error x: "<<(*iterOfExtTrk)->tof1PosSigmaAlongX(i)<<"\t"
346 <<"Error y: "<<(*iterOfExtTrk)->tof1PosSigmaAlongY(i)<<endl
347 <<"Tof: "<<(*iterOfExtTrk)->tof1(i)<<"\t"
348 <<"PathOF: "<<(*iterOfExtTrk)->tof1Path(i)
349 <<endl;
350 cout<<"******TOF2:******"<<endl;
351 cout<<"VolumeName: "<<(*iterOfExtTrk)->tof2VolumeName(i)<<"\t"
352 <<"VolumeNumber: "<<(*iterOfExtTrk)->tof2VolumeNumber(i)<<"\t"<<endl
353 <<"Position: "<<(*iterOfExtTrk)->tof2Position(i)<<"\t"
354 <<"Momentum: "<<(*iterOfExtTrk)->tof2Momentum(i)<<"\t"<<endl
355 <<"Error matrix: "<<(*iterOfExtTrk)->tof2ErrorMatrix(i)
356 <<"Error z: "<<(*iterOfExtTrk)->tof2PosSigmaAlongZ(i)<<"\t"
357 <<"Error Tz: "<<(*iterOfExtTrk)->tof2PosSigmaAlongT(i)<<"\t"
358 <<"Error x: "<<(*iterOfExtTrk)->tof2PosSigmaAlongX(i)<<"\t"
359 <<"Error y: "<<(*iterOfExtTrk)->tof2PosSigmaAlongY(i)<<endl
360 <<"Tof: "<<(*iterOfExtTrk)->tof2(i)<<"\t"
361 <<"PathOF: "<<(*iterOfExtTrk)->tof2Path(i)
362 <<endl;
363
364 //EMC information.
365 cout<<"******EMC:******"<<endl
366 <<"VolumeName: "<<(*iterOfExtTrk)->emcVolumeName(i)<<"\t"
367 <<"VolumeNumber: "<<(*iterOfExtTrk)->emcVolumeNumber(i)<<"\t"<<endl
368 <<"Position: "<<(*iterOfExtTrk)->emcPosition(i)<<"\t"
369 <<"Momentum: "<<(*iterOfExtTrk)->emcMomentum(i)<<"\t"<<endl
370 <<"Error matrix: "<<(*iterOfExtTrk)->emcErrorMatrix(i)
371 <<"Error theta: "<<(*iterOfExtTrk)->emcPosSigmaAlongTheta(i)<<"\t"
372 <<"Error phi: "<<(*iterOfExtTrk)->emcPosSigmaAlongPhi(i)<<"\t"
373 <<"EMC path: "<<(*iterOfExtTrk)->emcPath(i)
374 <<endl;
375
376 //MUC information
377 cout<<"******MUC:******"<<endl
378 <<"VolumeName: "<<(*iterOfExtTrk)->mucVolumeName(i)<<"\t"
379 <<"VolumeNumber: "<<(*iterOfExtTrk)->mucVolumeNumber(i)<<endl
380 <<"Position: "<<(*iterOfExtTrk)->mucPosition(i)<<"\t"
381 <<"Momentum: "<<(*iterOfExtTrk)->mucMomentum(i)<<"\t"<<endl
382 <<"Error matrix: "<<(*iterOfExtTrk)->mucErrorMatrix(i)
383 <<"Error z: "<<(*iterOfExtTrk)->mucPosSigmaAlongZ(i)<<"\t"
384 <<"Error Tz: "<<(*iterOfExtTrk)->mucPosSigmaAlongT(i)<<"\t"
385 <<"Error x: "<<(*iterOfExtTrk)->mucPosSigmaAlongX(i)<<"\t"
386 <<"Error y: "<<(*iterOfExtTrk)->mucPosSigmaAlongY(i)
387 <<endl;
388
389 cout<<"*******MUC KALMANFILTER***********"<<endl;
390 cout<<"Chisq is "<<(*iterOfExtTrk)->MucKalchi2(i)<<endl;
391 cout<<"Nfit is "<<(*iterOfExtTrk)->MucKaldof(i)<<endl;
392 cout<<"chiL "<<(*iterOfExtTrk)->MucKalchi2()<<endl;
393 // cout<<"Pull is "<<(*iterOfExtTrk)->mucKalPull()<<endl;
394
395 // cout<<"Residual is "<<(*iterOfExtTrk)->mucKalResidual()<<endl;
396 //Muc Ext hits information
397 /* ExtMucHitVec aExtMucHitVec = (*iterOfExtTrk)->GetExtMucHitVec();
398 int numOfMucHits = aExtMucHitVec.size();
399 cout<<"******MUC hits:"<<numOfMucHits<<"******"<<endl;
400 for(int j=0;j<numOfMucHits;j++)
401 {
402 cout<<"###Muc Hit "<<j<<":###"<<endl
403 <<"VolumeName: "<<aExtMucHitVec[j].GetVolumeName()<<"\t"
404 <<"VolumeNumber: "<<aExtMucHitVec[j].GetVolumeNumber()<<"\t"<<endl
405 <<"Position: "<<aExtMucHitVec[j].GetPosition()<<"\t"
406 <<"Momentum: "<<aExtMucHitVec[j].GetMomentum()<<"\t"<<endl
407 <<"Error z: "<<aExtMucHitVec[j].GetPosSigmaAlongZ()<<"\t"
408 <<"Error Tz: "<<aExtMucHitVec[j].GetPosSigmaAlongT()<<"\t"
409 <<"Error x: "<<aExtMucHitVec[j].GetPosSigmaAlongX()<<"\t"
410 <<"Error y: "<<aExtMucHitVec[j].GetPosSigmaAlongY()<<"\t"
411 <<endl;
412 }
413 */
414 }
415 }
416 j++;
417
418 } // loop ExtTrkCol
419
420 if(msgFlag) cout<<"****************** End a event! ****************"<<endl<<endl;
421
422 /*
423 //--------- For Ext Test ----------------
424 // Retrieve mc truth
425 SmartDataPtr<McParticleCol> mcParticleCol(eventSvc(),"/Event/MC/McParticleCol");
426 if (!mcParticleCol) {
427 log << MSG::FATAL << "Could not find McParticle" << endreq;
428 return( StatusCode::FAILURE);
429 }
430
431 int numOfTrack = mcParticleCol->size();
432 if(msgFlag) cout<< "numOfMcTrack: " << numOfTrack << endl;
433 if(numOfTrack!=2) return StatusCode::SUCCESS;
434
435 // Retrieve Emc Mc Hits
436 SmartDataPtr<EmcMcHitCol> emcMcHitCol(eventSvc(),"/Event/MC/EmcMcHitCol");
437 if (!emcMcHitCol) {
438 log << MSG::FATAL << "Could not find EMC truth" << endreq;
439 return( StatusCode::FAILURE);
440 }
441
442 McParticleCol::iterator iter_mc = mcParticleCol->begin();
443 EmcMcHitCol::iterator iterEmcBegin = emcMcHitCol->begin();
444 int numOfEmcHits = emcMcHitCol->size();
445
446 // ExtTrackCol *aExtTrackCol = new ExtTrackCol;
447
448 for (int i = 1;iter_mc != mcParticleCol->end(); iter_mc++, i++) {
449 bool flag = (*iter_mc)->primaryParticle();
450 if(!flag) continue;
451 int particleId = (*iter_mc)->particleProperty();
452 double charge = particleId/abs(particleId);
453 unsigned int sFlag = (*iter_mc)->statusFlags();
454 int trackIdx = (*iter_mc)->getTrackIndex();
455 HepPoint3D iPos = (*iter_mc)->initialPosition();
456 HepPoint3D fPos = (*iter_mc)->finalPosition();
457 HepLorentzVector iFMomentum = (*iter_mc)->initialFourMomentum();
458 HepLorentzVector fFMomentum = (*iter_mc)->finalFourMomentum();
459
460
461 bool emcHitFlag = false;
462 double thetaOfEmcHit = 0;
463 double phiOfEmcHit = 0;
464 Hep3Vector posOfEmcHit(0,0,0);
465 //Get Emc Truth
466 for(int j=0;j<numOfEmcHits;j++)
467 {
468 if(trackIdx==(*(iterEmcBegin+j))->getTrackIndex())
469 {
470 emcHitFlag = true;
471 double x = (*(iterEmcBegin+j))->getPositionX();
472 double y = (*(iterEmcBegin+j))->getPositionY();
473 double z = (*(iterEmcBegin+j))->getPositionZ();
474 Hep3Vector vec(x,y,z);
475 posOfEmcHit = vec;
476 thetaOfEmcHit = posOfEmcHit.theta();
477 phiOfEmcHit = posOfEmcHit.phi();
478 break;
479 }
480 }
481
482 ExtSteppingAction *extSteppingAction;
483 extSteppingAction = myExtTrack->GetStepAction();
484
485 G4String aParticleName(myParticleName);
486 if(!aParticleName.contains("proton"))
487 {
488 if(charge>0) aParticleName += "+";
489 else aParticleName += "-";
490 }
491 else
492 {
493 if(charge>0) aParticleName = "proton";
494 else aParticleName = "anti_proton";
495}
496
497ExtTrack *aExtTrack = new ExtTrack;
498aExtTrack->SetTrackID(trackIdx);
499
500Hep3Vector position(iPos);
501Hep3Vector momentum(iFMomentum.x(),iFMomentum.y(),iFMomentum.z());
502HepSymMatrix error(6,0);
503
504if(myExtTrack->Set(position,momentum,error,aParticleName,0))
505{
506 extSteppingAction->SetExtTrackPointer(aExtTrack);
507 myExtTrack->TrackExtrapotation();
508}
509
510Hep3Vector extEmcPos = aExtTrack->GetEmcPosition();
511double volumeNum = aExtTrack->GetEmcVolumeNumber();
512double thetaExt = extEmcPos.theta();
513double phiExt = extEmcPos.phi();
514
515if(myResultFlag)
516{
517 cout<< "*******Mc Track " << i <<" :******"<<endl;
518 cout<< "Parimary particle :" << flag <<" particleId = " << (*iter_mc)->particleProperty() << endl;
519 cout<< "Track Index: " << trackIdx << endl;
520 cout<< "initialPosition: "<< iPos.x() << "," << iPos.y() << "," << iPos.z() << endl;
521 cout<< "initialFourMomentum: " <<iFMomentum.x()<<","<<iFMomentum.y()<<","<<iFMomentum.z() << endl;
522 cout<<"Emc Truth:"<<emcHitFlag<<"!"<<posOfEmcHit<<endl;
523
524 cout<<"Ext Emc: "<<aExtTrack->GetEmcVolumeName()<<aExtTrack->GetEmcPosition()<<endl;
525}
526
527myCharge = charge;
528myEmcHitFlag = (emcHitFlag)? 1.:-1.;
529myEmcHitTheta = thetaOfEmcHit;
530myEmcHitPhi = phiOfEmcHit;
531myEmcVolNum = volumeNum;
532myEmcExtTheta = thetaExt;
533myEmcExtPhi = phiExt;
534myDTheta = myEmcHitTheta-myEmcExtTheta;
535myDPhi = myEmcHitPhi-myEmcExtPhi;
536while(myDTheta<-1*M_PI) myDTheta+=2.0*M_PI;
537while(myDPhi<-1*M_PI) myDPhi+=2.0*M_PI;
538while(myDPhi>M_PI) myDPhi-=2.0*M_PI;
539while(myDTheta>M_PI) myDTheta-=2.0*M_PI;
540myNtuple->write();
541
542if(aExtTrack) delete aExtTrack;
543}
544
545//--------- end Ext Test ----------------
546*/
547
548return StatusCode::SUCCESS;
549}
550
551/////////////////////
552StatusCode TrkExtAlg::finalize()
553{
554 MsgStream log(msgSvc(), name());
555 log << MSG::INFO << "finalize()" << endreq;
556
557 // delete myExtTrack;
558 // myFile->close();
559
560 return StatusCode::SUCCESS;
561}
momentum
**********INTEGER nmxhep !maximum number of particles DOUBLE PRECISION vhep INTEGER jdahep COMMON hepevt $ !serial number $ !number of particles $ !status code $ !particle ident KF $ !parent particles $ !childreen particles $ !four momentum
Definition: BesBdkRc/BesBdkRc-00-00-02/src/fortran/HepEvt.h:23
RecExtTrackCol
ObjectVector< RecExtTrack > RecExtTrackCol
Definition: RecExtTrack.h:209
msgSvc
IMessageSvc * msgSvc()
Definition: ServiceAccessor.h:13
detVerSvc
IDetVerSvc * detVerSvc
Definition: TofRawDataProvider.cxx:31
parName
string parName[5]
Definition: TrkExtAlg.cxx:30
DstExtTrack::SetTrackId
void SetTrackId(int trackId)
Definition: DstExtTrack.h:193
DstExtTrack::SetParType
void SetParType(int aParType=2)
Definition: DstExtTrack.h:194
ExtMdcTrack::GetOneGoodTrk
bool GetOneGoodTrk()
Definition: ExtMdcTrack.cxx:51
ExtMdcTrack::GetParticleCharge
double GetParticleCharge() const
Definition: ExtMdcTrack.cxx:356
ExtMdcTrack::GetTrackLength
double GetTrackLength() const
Definition: ExtMdcTrack.cxx:309
ExtMdcTrack::GetTrkTof
double GetTrkTof() const
Definition: ExtMdcTrack.h:56
ExtMdcTrack::GetTrackID
int GetTrackID()
Definition: ExtMdcTrack.h:50
ExtMdcTrack::SetMdcKalTrkCol
bool SetMdcKalTrkCol(RecMdcKalTrackCol *aPointer)
Definition: ExtMdcTrack.cxx:42
ExtMdcTrack::SetMdcRecTrkCol
bool SetMdcRecTrkCol(RecMdcTrackCol *aPointer)
Definition: ExtMdcTrack.cxx:33
ExtMdcTrack::GetPosition
const Hep3Vector GetPosition() const
Definition: ExtMdcTrack.cxx:227
ExtMdcTrack::ReadTrk
bool ReadTrk(int pid)
Definition: ExtMdcTrack.cxx:77
ExtMdcTrack::GetErrorMatrix
const HepSymMatrix GetErrorMatrix() const
Definition: ExtMdcTrack.cxx:251
ExtMdcTrack::GetMomentum
const Hep3Vector GetMomentum() const
Definition: ExtMdcTrack.cxx:242
ExtMdcTrack::SetMsgFlag
void SetMsgFlag(bool aFlag)
Definition: ExtMdcTrack.h:44
ExtSteppingAction::SetExtTrackPointer
void SetExtTrackPointer(RecExtTrack *aExtTrack)
Definition: ExtSteppingAction.h:44
ExtSteppingAction::InfmodMuc
void InfmodMuc(Hep3Vector &pos, Hep3Vector &mom, HepSymMatrix &err)
Definition: ExtSteppingAction.cxx:1018
ExtSteppingAction::Set_which_tof_version
void Set_which_tof_version(int version)
Definition: ExtSteppingAction.h:50
ExtSteppingAction::SetMucDigiColPointer
void SetMucDigiColPointer(MucDigiCol *rawdigicol)
Definition: ExtSteppingAction.h:56
Ext_track
Definition: Ext_track.h:40
Ext_track::TrackExtrapotation
void TrackExtrapotation()
Definition: Ext_track.cxx:343
Ext_track::Initialization
void Initialization(const bool aMsgFlag, const bool Bfield, const bool GeomOptimization, const bool aUseMucKal, const int aMucWindow)
Definition: Ext_track.cxx:90
Ext_track::GetStepAction
ExtSteppingAction * GetStepAction()
Definition: Ext_track.h:90
Ext_track::Set
bool Set(const Hep3Vector &xv3, const Hep3Vector &pv3, const HepSymMatrix &err, const std::string &particleName, const double pathInMDC, const double tofInMdc)
Definition: Ext_track.cxx:235
IDetVerSvc::phase
virtual int phase()=0
TrkExtAlg::finalize
StatusCode finalize()
Definition: TrkExtAlg.cxx:552
TrkExtAlg::execute
StatusCode execute()
Definition: TrkExtAlg.cxx:106
TrkExtAlg::TrkExtAlg
TrkExtAlg(const std::string &name, ISvcLocator *pSvcLocator)
Definition: TrkExtAlg.cxx:33
TrkExtAlg::initialize
StatusCode initialize()
Definition: TrkExtAlg.cxx:58
Event
Definition: Event.h:21
charge
float charge
Definition: prepare_data.cxx:28