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