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"
15#include "CLHEP/Matrix/SymMatrix.h"
16#include "CLHEP/Vector/ThreeVector.h"
29string parName[5]={
"e",
"mu",
"pi",
"kaon",
"proton"};
34 myParticleName =
"pi";
39 declareProperty(
"ParticleName",myParticleName);
40 declareProperty(
"GeantGeomOptimization",myGeomOptimization=
true);
41 declareProperty(
"MessageFlag",msgFlag);
42 declareProperty(
"ResultMessageFlag",myResultFlag);
43 declareProperty(
"BFieldOn",myBFieldOn=
true);
44 declareProperty(
"InputTrk",myInputTrk);
45 declareProperty(
"Tof",m_tof=2);
46 declareProperty(
"UseMucKalFilter",myUseMucKal=
true);
47 declareProperty(
"MucWindow",myMucWindow=6);
53 if(myExtTrack)
delete myExtTrack;
59 MsgStream log(
msgSvc(), name());
60 log << MSG::INFO <<
"initialize()" << endreq;
62 myExtTrack =
new Ext_track(msgFlag,myBFieldOn,myGeomOptimization,m_tof,myUseMucKal,myMucWindow);
63 myExtTrack->
Initialization(msgFlag,myBFieldOn,myGeomOptimization,myUseMucKal,myMucWindow);
90 return StatusCode::SUCCESS;
98 MsgStream log(
msgSvc(), name());
99 log << MSG::INFO <<
"execute()" << endreq;
100 int eventNumber, runNumber;
102 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),
"/Event/EventHeader");
105 log << MSG::FATAL <<
"Could not find Event Header" << endreq;
106 return( StatusCode::FAILURE);
108 runNumber = eventHeader->runNumber();
109 eventNumber = eventHeader->eventNumber();
113 cout<<
"TrackExt: ******************* Start a event *******************"<<endl;
114 cout<<
"run= "<<runNumber<<
"; event= "<<eventNumber<<endl;
118 SmartDataPtr<MucDigiCol> mucDigiCol(eventSvc(),
"/Event/Digi/MucDigiCol");
120 log << MSG::FATAL <<
"Could not find MUC digi" << endreq;
121 return( StatusCode::SUCCESS);
131 if(myParticleName==
"e") parID=0;
132 else if(myParticleName==
"mu") parID=1;
133 else if(myParticleName==
"pi") parID=2;
134 else if(myParticleName==
"kaon") parID=3;
135 else if(myParticleName==
"proton"||myParticleName==
"anti_proton") parID=4;
137 if(myInputTrk==
"Mdc")
139 SmartDataPtr<RecMdcTrackCol> aMdcTrackCol(eventSvc(),
"/Event/Recon/RecMdcTrackCol");
142 log << MSG::WARNING <<
"Can't find RecMdcTrackCol in TDS!" << endreq;
143 return( StatusCode::SUCCESS);
148 else if(myInputTrk==
"Kal")
150 SmartDataPtr<RecMdcKalTrackCol> aMdcKalTrackCol(eventSvc(),
"/Event/Recon/RecMdcKalTrackCol");
153 log << MSG::WARNING <<
"Can't find RecMdcKalTrackCol in TDS!" << endreq;
154 return( StatusCode::SUCCESS);
160 log << MSG::WARNING <<
"Wrong type of inputTrk:" << myInputTrk << endreq;
161 return( StatusCode::SUCCESS);
172 for(
int i=0; i<5; i++)
184 double tofInMdc = aExtMdcTrack.
GetTrkTof();
188 cout<<
"Start From:"<<position.x()<<
' '<<position.y()<<
' '
189 <<position.z()<<endl;
191 cout<<
"Start Error matrix:"<<error<<endl;
192 cout<<
"Path before start:"<< pathInMDC << endl;
195 G4String aParticleName(
parName[i]);
197 if(!aParticleName.contains(
"proton"))
199 if(charge>0) aParticleName +=
"+";
200 else aParticleName +=
"-";
204 if(charge>0) aParticleName =
"proton";
205 else aParticleName =
"anti_proton";
210 cout<<
"Charge: "<<charge<<endl;
211 cout<<
"Particle: "<<aParticleName<<endl;
218 extSteppingAction->
Reset();
221 bool m_trackstatus =
false;
222 int trk_startpart = 0;
223 while(!m_trackstatus)
228 {cout<<
"-------has modified more than 20 times---------"<<endl;
break;}
229 if(myExtTrack->
Set(position,
momentum,error,aParticleName,pathInMDC,tofInMdc))
233 m_trackstatus = extSteppingAction->
TrackStop();
236 m_trackstatus =
true;
244 if(msgFlag) cout<<
"will add aExtTrack!"<<endl;
247 if(aExtTrack) aExtTrackCol->add(aExtTrack);
248 else if(msgFlag) cout<<
"No aExtTrack!"<<endl;
252 if(msgFlag) cout<<
"No aExtTrackCol!"<<endl;
254 if(msgFlag) cout<<
"add a aExtTrack!"<<endl;
286 SmartIF<IDataManagerSvc> dataManSvc(eventSvc());
289 DataObject *extTrackCol;
290 eventSvc()->findObject(
"/Event/Recon/RecExtTrackCol",extTrackCol);
291 if(extTrackCol != NULL) {
292 dataManSvc->clearSubTree(
"/Event/Recon/RecExtTrackCol");
293 eventSvc()->unregisterObject(
"/Event/Recon/RecExtTrackCol");
298 StatusCode sc = eventSvc()->registerObject(
"/Event/Recon/RecExtTrackCol", aExtTrackCol);
299 if(sc!=StatusCode::SUCCESS) {
300 log << MSG::FATAL <<
"Could not register RecExtTrackCol in TDS!" << endreq;
301 return( StatusCode::FAILURE);
306 SmartDataPtr<RecExtTrackCol> aExtTrkCol(eventSvc(),
"/Event/Recon/RecExtTrackCol");
308 log << MSG::FATAL <<
"Can't find RecExtTrackCol in TDS!" << endreq;
309 return( StatusCode::FAILURE);
312 RecExtTrackCol::iterator iterOfExtTrk;
315 for(iterOfExtTrk = aExtTrkCol->begin();iterOfExtTrk!=aExtTrkCol->end();iterOfExtTrk++)
319 for(
int i=0; i<5; i++)
322 cout<<
"##########track"<<j<<
": "<<
"("<<i<<
")"<<endl;
323 cout<<
"******TOF1:******"<<endl;
324 cout<<
"VolumeName: "<<(*iterOfExtTrk)->tof1VolumeName(i)<<
"\t"
325 <<
"VolumeNumber: "<<(*iterOfExtTrk)->tof1VolumeNumber(i)<<
"\t"<<endl
326 <<
"Position: "<<(*iterOfExtTrk)->tof1Position(i)<<
"\t"
327 <<
"Momentum: "<<(*iterOfExtTrk)->tof1Momentum(i)<<
"\t"<<endl
328 <<
"Error matrix: "<<(*iterOfExtTrk)->tof1ErrorMatrix(i)
329 <<
"Error z: "<<(*iterOfExtTrk)->tof1PosSigmaAlongZ(i)<<
"\t"
330 <<
"Error Tz: "<<(*iterOfExtTrk)->tof1PosSigmaAlongT(i)<<
"\t"
331 <<
"Error x: "<<(*iterOfExtTrk)->tof1PosSigmaAlongX(i)<<
"\t"
332 <<
"Error y: "<<(*iterOfExtTrk)->tof1PosSigmaAlongY(i)<<endl
333 <<
"Tof: "<<(*iterOfExtTrk)->tof1(i)<<
"\t"
334 <<
"PathOF: "<<(*iterOfExtTrk)->tof1Path(i)
336 cout<<
"******TOF2:******"<<endl;
337 cout<<
"VolumeName: "<<(*iterOfExtTrk)->tof2VolumeName(i)<<
"\t"
338 <<
"VolumeNumber: "<<(*iterOfExtTrk)->tof2VolumeNumber(i)<<
"\t"<<endl
339 <<
"Position: "<<(*iterOfExtTrk)->tof2Position(i)<<
"\t"
340 <<
"Momentum: "<<(*iterOfExtTrk)->tof2Momentum(i)<<
"\t"<<endl
341 <<
"Error matrix: "<<(*iterOfExtTrk)->tof2ErrorMatrix(i)
342 <<
"Error z: "<<(*iterOfExtTrk)->tof2PosSigmaAlongZ(i)<<
"\t"
343 <<
"Error Tz: "<<(*iterOfExtTrk)->tof2PosSigmaAlongT(i)<<
"\t"
344 <<
"Error x: "<<(*iterOfExtTrk)->tof2PosSigmaAlongX(i)<<
"\t"
345 <<
"Error y: "<<(*iterOfExtTrk)->tof2PosSigmaAlongY(i)<<endl
346 <<
"Tof: "<<(*iterOfExtTrk)->tof2(i)<<
"\t"
347 <<
"PathOF: "<<(*iterOfExtTrk)->tof2Path(i)
351 cout<<
"******EMC:******"<<endl
352 <<
"VolumeName: "<<(*iterOfExtTrk)->emcVolumeName(i)<<
"\t"
353 <<
"VolumeNumber: "<<(*iterOfExtTrk)->emcVolumeNumber(i)<<
"\t"<<endl
354 <<
"Position: "<<(*iterOfExtTrk)->emcPosition(i)<<
"\t"
355 <<
"Momentum: "<<(*iterOfExtTrk)->emcMomentum(i)<<
"\t"<<endl
356 <<
"Error matrix: "<<(*iterOfExtTrk)->emcErrorMatrix(i)
357 <<
"Error theta: "<<(*iterOfExtTrk)->emcPosSigmaAlongTheta(i)<<
"\t"
358 <<
"Error phi: "<<(*iterOfExtTrk)->emcPosSigmaAlongPhi(i)<<
"\t"
359 <<
"EMC path: "<<(*iterOfExtTrk)->emcPath(i)
363 cout<<
"******MUC:******"<<endl
364 <<
"VolumeName: "<<(*iterOfExtTrk)->mucVolumeName(i)<<
"\t"
365 <<
"VolumeNumber: "<<(*iterOfExtTrk)->mucVolumeNumber(i)<<endl
366 <<
"Position: "<<(*iterOfExtTrk)->mucPosition(i)<<
"\t"
367 <<
"Momentum: "<<(*iterOfExtTrk)->mucMomentum(i)<<
"\t"<<endl
368 <<
"Error matrix: "<<(*iterOfExtTrk)->mucErrorMatrix(i)
369 <<
"Error z: "<<(*iterOfExtTrk)->mucPosSigmaAlongZ(i)<<
"\t"
370 <<
"Error Tz: "<<(*iterOfExtTrk)->mucPosSigmaAlongT(i)<<
"\t"
371 <<
"Error x: "<<(*iterOfExtTrk)->mucPosSigmaAlongX(i)<<
"\t"
372 <<
"Error y: "<<(*iterOfExtTrk)->mucPosSigmaAlongY(i)
375 cout<<
"*******MUC KALMANFILTER***********"<<endl;
376 cout<<
"Chisq is "<<(*iterOfExtTrk)->MucKalchi2(i)<<endl;
377 cout<<
"Nfit is "<<(*iterOfExtTrk)->MucKaldof(i)<<endl;
378 cout<<
"chiL "<<(*iterOfExtTrk)->MucKalchi2()<<endl;
406 if(msgFlag) cout<<
"****************** End a event! ****************"<<endl<<endl;
534return StatusCode::SUCCESS;
540 MsgStream log(
msgSvc(), name());
541 log << MSG::INFO <<
"finalize()" << endreq;
546 return StatusCode::SUCCESS;
**********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
ObjectVector< RecExtTrack > RecExtTrackCol
void SetTrackId(int trackId)
void SetParType(int aParType=2)
double GetParticleCharge() const
double GetTrackLength() const
bool SetMdcKalTrkCol(RecMdcKalTrackCol *aPointer)
bool SetMdcRecTrkCol(RecMdcTrackCol *aPointer)
const Hep3Vector GetPosition() const
const HepSymMatrix GetErrorMatrix() const
const Hep3Vector GetMomentum() const
void SetMsgFlag(bool aFlag)
void SetExtTrackPointer(RecExtTrack *aExtTrack)
void InfmodMuc(Hep3Vector &pos, Hep3Vector &mom, HepSymMatrix &err)
void Set_which_tof_version(int version)
void SetMucDigiColPointer(MucDigiCol *rawdigicol)
void TrackExtrapotation()
void Initialization(const bool aMsgFlag, const bool Bfield, const bool GeomOptimization, const bool aUseMucKal, const int aMucWindow)
ExtSteppingAction * GetStepAction()
bool Set(const Hep3Vector &xv3, const Hep3Vector &pv3, const HepSymMatrix &err, const std::string &particleName, const double pathInMDC, const double tofInMdc)
TrkExtAlg(const std::string &name, ISvcLocator *pSvcLocator)