EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
EicIdealTrackingCode.cxx
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file EicIdealTrackingCode.cxx
1 //
2 // //AYK, 2013/02/28:
3 //
4 // Well, there was no easy way to inherit from any generic enough PandaRoot class,
5 // so considered to copy over PndFtsTrackerIdeal.h(cxx) and customize/generalize them
6 // to my taste; will actually have to look more attentively into these codes later
7 // and clean them up; in particular figure out how to get rid of "PndTrack" prefices;
8 //
9 
10 #include "PndMCTrack.h"
11 #include "PndTrackCand.h"
12 #include "PndTrack.h"
13 
14 #include "FairTrackParP.h"
15 #include "FairMCPoint.h"
16 #include "FairHit.h"
17 #include "FairTask.h"
18 #include "FairRootManager.h"
19 
20 #include "TObjectTable.h"
21 #include "TClonesArray.h"
22 #include "TParticlePDG.h"
23 #include "TRandom3.h"
24 #include <iostream>
25 
26 #include <EicDetName.h>
27 #include "EicIdealTrackingCode.h"
28 
29 // -----------------------------------------------------------------------------------------------
30 
31 EicIdealTrackingCode::EicIdealTrackingCode():
32  FairTask("EIC Ideal Tracking Code"), fMCTracks(new TClonesArray()),
33  fTrackIds(new TClonesArray()), fMomSigma(0,0,0), fDPoP(0.), fRelative (kFALSE),
34  fVtxSigma(0,0,0), fEfficiency(1.),
35  fTracksArrayName("EicIdealTrack"), pdg(0)
36 {
37  fTrackCands = new TClonesArray("PndTrackCand");
38  fTracks = new TClonesArray("PndTrack");
39  fVerbose = 0;
40 } // EicIdealTrackingCode::EicIdealTrackingCode()
41 
42 // -----------------------------------------------------------------------------------------------
43 
44 EicIdealTrackingCode::~EicIdealTrackingCode()
45 {
46  //FairRootManager *fManager =FairRootManager::Instance();
47  //fManager->Write();
48 } // EicIdealTrackingCode::~EicIdealTrackingCode()
49 
50 // -----------------------------------------------------------------------------------------------
51 
52 void EicIdealTrackingCode::Register()
53 {
54  FairRootManager::Instance()->Register(fTracksArrayName,"EicTrk", fTracks, kTRUE);
55  FairRootManager::Instance()->Register(fTracksArrayName+"Cand","EicTrk", fTrackCands, kTRUE);
56  if(fVerbose>3) Info("Register","Done.");
57 } // EicIdealTrackingCode::Register()
58 
59 // -----------------------------------------------------------------------------------------------
60 
61 int EicIdealTrackingCode::AddDetectorGroup(const char *name)
62 {
63  EicDetectorGroup group(name);
64 
65  fGroups.push_back(group);
66 
67  return 0;
68 } // EicIdealTrackingCode::AddDetectorGroup()
69 
70 // -----------------------------------------------------------------------------------------------
71 
72 InitStatus EicIdealTrackingCode::Init()
73 {
74  if(fVerbose>3) Info("Init","Start initialisation.");
75 
76  FairRootManager *fManager = FairRootManager::Instance();
77 
78  // Get MC track array;
79  fMCTracks = (TClonesArray*)fManager->GetObject("MCTrack");
80  if ( ! fMCTracks ) {
81  std::cout << "-W- EicIdealTrackingCode::Init: No MCTrack array! Needed for MC Truth" << std::endl;
82  return kERROR;
83  }
84 
85  // Get detector-specific MC & digi arrays;
86  for (std::vector<EicDetectorGroup>::iterator it=fGroups.begin(); it!=fGroups.end(); it++)
87  {
88  it->_fMCPoints = (TClonesArray*)fManager->GetObject(it->dname->Name() + "MoCaPoint");
89  if ( ! it->_fMCPoints ) {
90  std::cout << "-W- EicIdealTrackingCode::Init: No " << it->dname->Name() << "Point array!" << std::endl;
91  return kERROR;
92  }
93 
94  //it->_fHits = (TClonesArray *)fManager->GetObject(it->dname->Name() + "DigiHit");
95  it->_fHits = (TClonesArray *)fManager->GetObject(it->dname->Name() + "TrackingDigiHit");
96  if ( ! it->_fHits ) {
97  std::cout << "-W- EicIdealTrackingCode::Init: No " << it->dname->Name() << "Hit array!" << std::endl;
98  return kERROR;
99  }
100 
101  //it->_fBranchID = FairRootManager::Instance()->GetBranchId(it->dname->Name() + "DigiHit");
102  it->_fBranchID = fManager->GetBranchId(it->dname->Name() + "TrackingDigiHit");
103 
104  // And create direct access lookup table entries; yes, all this looks like back
105  // doors into the original FairRoot codes;
106  fManager->AddMoCaLookupEntry(it->_fBranchID, it->_fMCPoints);
107  fManager->AddDigiLookupEntry(it->_fBranchID, it->_fHits);
108  } //for it
109 
110  if(fVerbose>3) Info("Init","Fetched all arrays.");
111 
112  Register();
113 
114  pdg = new TDatabasePDG();
115  if(fVerbose>3) Info("Init","End initialisation.");
116 
117  return kSUCCESS;
118 } // EicIdealTrackingCode::Init()
119 
120 // -----------------------------------------------------------------------------------------------
121 
122 void EicIdealTrackingCode::Exec(Option_t * option)
123 {
124  Reset();
125  if(fVerbose>3) Info("Exec","Start eventloop.");
126 
127  if(fVerbose>4) Info("Exec","Print some array properties");
128 
129  FairHit* ghit = NULL;
130  std::map<Int_t, FairHit*> firstHit;
131  std::map<Int_t, FairHit*> lastHit;
132  FairMCPoint* myPoint=NULL;
133  std::map<Int_t, FairMCPoint*> firstPoint;
134  std::map<Int_t, FairMCPoint*> lastPoint;
135  std::map<Int_t, PndTrackCand*> candlist;
136 
137  // Loop through all detector groups (TRS, FGT, etc);
138  for (std::vector<EicDetectorGroup>::iterator it=fGroups.begin(); it!=fGroups.end(); it++)
139  {
140  for (Int_t ih = 0; ih < it->_fHits->GetEntriesFast(); ih++) {
141  ghit = (FairHit*) it->_fHits->At(ih);
142  if(!ghit) {
143  if(fVerbose>3) Error("Exec","Have no ghit %i, array size: %i",ih, it->_fHits->GetEntriesFast());
144  continue;
145  }
146  Int_t mchitid=ghit->GetRefIndex();
147  //printf("%d\n", mchitid);
148  if(mchitid<0) {
149  if(fVerbose>3) Error("Exec","Have a mcHit %i",mchitid);
150  continue;
151  }
152  myPoint = (FairMCPoint*)(it->_fMCPoints->At(mchitid));
153  if(!myPoint) continue;
154  Int_t trackID = myPoint->GetTrackID();
155  if(trackID<0) continue;
156 
157  if(fVerbose>5) Info("Exec","Have a Hit %i at Track index %i",ih,trackID);
158  if(fVerbose>5) Info("Exec"," --> mchitid %i",mchitid);
159 
160  PndTrackCand* cand=candlist[trackID];
161  if(NULL==cand){
162  if(fVerbose>5) Info("Exec","Create new PndTrack object %i",trackID);
163  cand=new PndTrackCand();
164  cand->setMcTrackId(trackID);
165  if(fVerbose>5) Info("Exec","Create new PndTrack object finished %i",trackID);
166  }
167  if(fVerbose>5) Info("Exec","add the hit %i to trackcand %i",ih,trackID);
168 #if _THINK_
169  // Well, it looks like distance from the IP should be a good guess on hit
170  // ordering parameter; not the best idea -> THINK later, need much better algorithm
171  // (order along the trajectory, whatever this means);
172  //cand->AddHit(it->_fBranchID,ih,ghit->GetZ());
173 #endif
174  //printf("%3d -> %7.2f -> %7.2f %7.2f %7.2f\n", ih, myPoint->GetTime(),
175  // ghit->GetX(), ghit->GetY(), ghit->GetZ());
176  cand->AddHit(it->_fBranchID,ih,sqrt(ghit->GetX()*ghit->GetX() +
177  ghit->GetY()*ghit->GetY() + ghit->GetZ()*ghit->GetZ()));
178  //cand->AddHit(it->_fBranchID, ih, myPoint->GetTime());
179  //cand->AddHit(it->_fBranchID,ih,sqrt(ghit->GetX()*ghit->GetX() + ghit->GetY()*ghit->GetY()));
180  if(!firstHit[trackID] || firstHit[trackID]->GetZ() > ghit->GetZ()) {
181  firstHit[trackID]=ghit;
182  firstPoint[trackID]=myPoint;
183  }
184  if(!lastHit[trackID] || lastHit[trackID]->GetZ() < ghit->GetZ()) {
185  lastHit[trackID]=ghit;
186  lastPoint[trackID]=myPoint;
187  }
188 
189  candlist[trackID] = cand; // set
190  }
191  }
192  if(fVerbose>3) Info("Exec","Insert to TCA (depending on efficiency)");
193 
194  // re-iterate and select by efficiency & number of hits
195  std::map<Int_t, PndTrackCand*>::iterator candit;
196  PndMCTrack *mc=NULL;
197  TVector3 svtx, smom;
198  Int_t charge=0, trackID=-1;
199 
200  for(candit=candlist.begin(); candit!=candlist.end(); ++candit) {
201  PndTrackCand* tcand=candit->second;
202  trackID=candit->first;
203  if(!tcand) {
204  if(fVerbose>3) Warning("Exec","Have no candidate at %i",trackID);
205  continue;
206  }
207  if( tcand->GetNHits() < 3 ) continue;
208  if(0 < fEfficiency && fEfficiency < 1){
209  if(gRandom->Rndm() > fEfficiency) continue;
210  }
211  tcand->Sort();
212  mc = (PndMCTrack*)fMCTracks->At(trackID);
213  if (mc->GetPdgCode()<100000000)
214  charge = (Int_t)TMath::Sign(1.0, ((TParticlePDG*)pdg->GetParticle(mc->GetPdgCode()))->Charge());
215  else charge = 1;
216  tcand->setMcTrackId(trackID);
217  // prepare track parameters
218  firstHit[trackID]->Position(svtx); // set position to first hit
219  SmearFWD(svtx, fVtxSigma);
220  firstPoint[trackID]->Momentum(smom);
221  if (fRelative) fMomSigma=fDPoP*smom;
222  SmearFWD(smom, fMomSigma);
223  //FairTrackParP* firstPar=new FairTrackParP(svtx, smom,
224  // fVtxSigma, fMomSigma,
225  // charge, svtx,
226  // TVector3(1.,0.,0.), TVector3(0.,1.,0.));
227  FairTrackParP firstPar(svtx, smom,
228  fVtxSigma, fMomSigma,
229  charge, svtx,
230  TVector3(1.,0.,0.), TVector3(0.,1.,0.));
231 
232  lastHit[trackID]->Position(svtx);
233  SmearFWD(svtx, fVtxSigma);
234  lastPoint[trackID]->Momentum(smom);
235  SmearFWD(smom, fMomSigma);
236  //FairTrackParP* lastPar=new FairTrackParP(svtx, smom,
237  // fVtxSigma, fMomSigma,
238  // charge, svtx,
239  // TVector3(1.,0.,0.), TVector3(0.,1.,0.));
240  FairTrackParP lastPar(svtx, smom,
241  fVtxSigma, fMomSigma,
242  charge, svtx,
243  TVector3(1.,0.,0.), TVector3(0.,1.,0.));
244 
245  if(fVerbose>3) Info("Exec","Store candidate at %i",trackID);
246  // Creates a new hit in the TClonesArray.
247  if(fVerbose>3) Info("AddTrack","Adding a Track.");
248  TClonesArray &pndtracks = *fTracks;
249  TClonesArray &pndtrackcands = *fTrackCands;
250  // TClonesArray &pndtrackids = *fTrackIds;
251  Int_t size = pndtrackcands.GetEntriesFast();
252  if(pndtracks.GetEntriesFast() != size) {
253  Error("Exec","Arrays out of synchronisation: %i tracks, %i cands. Abort event."
254  ,pndtracks.GetEntriesFast(),pndtrackcands.GetEntriesFast());
255  return;
256  }
257 
258 #if 1
259  PndTrackCand* pndTrackCand = new(pndtrackcands[size]) PndTrackCand(*tcand);
260  PndTrack* pndTrack = new(pndtracks[size])
261  PndTrack(firstPar, lastPar, *tcand,0,0,1,mc->GetPdgCode(),trackID,
262  FairRootManager::Instance()->GetBranchId("MCTrack"));
263  //PndTrack(firstPar, lastPar, *tcand,0,0,1,mc->GetPdgCode(),trackID,
264  // FairRootManager::Instance()->GetBranchId("MCTrack"));
265 #endif
266  }
267 
268  // Clean up candidate list;
269  for(candit=candlist.begin(); candit!=candlist.end(); ++candit)
270  delete candit->second;
271  {
272  //static unsigned counter;
273  //printf("%6d\n", counter++);
274  }
275 
276  if(fVerbose>3) Info("Exec","End eventloop.");
277 } // EicIdealTrackingCode::Exec()
278 
279 // -----------------------------------------------------------------------------------------------
280 
281 #include <FairRun.h>
282 
283 void EicIdealTrackingCode::Finish()
284 {
285  std::cout << " Found "<< fTracks->GetEntriesFast() << " tracks\n";
286 
287  FairRun *fRun = FairRun::Instance();
288 
289  // I guess there is no need to save/restore current directory here?;
290  fRun->GetOutputFile()->cd();
291  Write("EicIdealTrackingCode");
292 } // EicIdealTrackingCode::Finish()
293 
294 // -----------------------------------------------------------------------------------------------
295 
296 void EicIdealTrackingCode::Reset()
297 {
298  //for(unsigned iq=0; iq<fTracks->GetEntriesFast(); iq++) {
299  //delete (PndTrack*)(fTracks->At(iq));
300  //fTracks[iq] = 0;
301  //}
302  //for(unsigned iq=0; iq<fTrackCands->GetEntriesFast(); iq++)
303  //delete fTrackCands->At(iq);
304 
305  //if (fTracks->GetEntriesFast() != 0) fTracks->Clear();
306  if (fTracks->GetEntriesFast() != 0) fTracks->Delete();
307  //if (fTrackCands->GetEntriesFast() != 0) fTrackCands->Clear();
308  if (fTrackCands->GetEntriesFast() != 0) fTrackCands->Delete();
309 } // EicIdealTrackingCode::Reset()
310 
311 // -----------------------------------------------------------------------------------------------
312 
313 void EicIdealTrackingCode::SmearFWD(TVector3 &vec, const TVector3 &sigma)
314 {
315  // gaussian smearing
316  Double_t rannn=0.;
317  rannn = gRandom->Gaus(vec.X(),sigma.X());
318  vec.SetX(rannn);
319 
320  rannn = gRandom->Gaus(vec.Y(),sigma.Y());
321  vec.SetY(rannn);
322 
323  rannn = gRandom->Gaus(vec.Z(),sigma.Z());
324  vec.SetZ(rannn);
325 
326  return;
327 } // EicIdealTrackingCode::SmearFWD()
328 
329 // -----------------------------------------------------------------------------------------------
330 
331 ClassImp(EicDetectorGroup)
332 ClassImp(EicIdealTrackingCode)
333