FairRoot/PandaRoot
PndGemTrackFinderOnHitsTB.cxx
Go to the documentation of this file.
1 // -------------------------------------------------------------------------
2 // ----- PndGemTrackFinderOnHitsTB source file -----
3 // ----- Created 02.06.2009 by R. Karabowicz -----
4 // -------------------------------------------------------------------------
5 
6 #include "PndGemTrackFinderOnHitsTB.h"
7 
8 // Pnd includes
9 #include "FairRootManager.h"
10 #include "FairRunAna.h"
11 #include "FairRuntimeDb.h"
12 #include "FairBaseParSet.h"
13 #include "FairTrackParam.h"
14 #include "FairTrackParP.h"
15 #include "PndGemMCPoint.h"
16 #include "PndGemDigiPar.h"
17 #include "FairRootManager.h"
18 #include "PndDetectorList.h"
19 #include "PndTrack.h"
20 #include "PndTrackCand.h"
21 #include "PndTrackCandHit.h"
22 // ROOT includes
23 #include "TClonesArray.h"
24 #include "TGeoManager.h"
25 #include "TMatrixFSym.h"
26 #include "TStopwatch.h"
27 
28 // C++ includes
29 #include <iostream>
30 #include <iomanip>
31 #include <map>
32 #include <cmath>
33 using std::cout;
34 using std::endl;
35 using std::map;
36 
37 // ----- Default constructor -------------------------------------------
38 PndGemTrackFinderOnHitsTB::PndGemTrackFinderOnHitsTB() :
39  PndGemTrackFinder(),
40  fDigiPar(NULL),
41  fMCTrackArray(NULL),
42  fMCPointArray(NULL),
43  fVerbose(0),
44  fPrimary(0),
45  fSigmaMult(60.),
46  fParThetaA(0.),
47  fParThetaB(0.),
48  fParTheta0(0.),
49  fParTheta1(0.),
50  fParTheta2(0.),
51  fParTheta3(0.),
52  fParRadPhi0(0.),
53  fParRadPhi2(0.),
54  fPrepTime(0.),
55  fSegmTime(0.),
56  fMatchTime(0.),
57  fRemoveTime(0.),
58  fWriteTime(0.),
59  fAllTime(0.),
60  fNofEvents(0),
61  fMCAvailable(-666),
62  fNofClHits(0),
63  fNofExpectedTrackSegments(0),
64  fNofFoundTrackSegments(0)
65 {
66  for ( Int_t ipar = 0 ; ipar < 3 ; ipar++ ) {
67  fParMat0[ipar] = 0.;
68  fParMat1[ipar] = 0.;
69  }
70 }
71 
72 // ----- Destructor ----------------------------------------------------
73 PndGemTrackFinderOnHitsTB::~PndGemTrackFinderOnHitsTB() { }
74 
75 // ----- Init -----------------------------------------------------------
76 void PndGemTrackFinderOnHitsTB::Init() {
77 
78  fNofExpectedTrackSegments = 0;
79  fNofFoundTrackSegments = 0;
80 
81  // Get and check FairRootManager
82  FairRootManager* ioman = FairRootManager::Instance();
83  if( !ioman ) {
84  cout << "-E- "<< GetName() <<"::Init: "
85  << "RootManager not instantised!" << endl;
86  return;
87  }
88 
89  // Get the pointer to the singleton FairRunAna object
90  FairRunAna* ana = FairRunAna::Instance();
91  if(NULL == ana) {
92  cout << "-E- "<< GetName() <<"::Init :"
93  <<" no FairRunAna object!" << endl;
94  return;
95  }
96  // Get the pointer to run-time data base
97  FairRuntimeDb* rtdb = ana->GetRuntimeDb();
98  if(NULL == rtdb) {
99  cout << "-E- "<< GetName() <<"::Init :"
100  <<" no runtime database!" << endl;
101  return;
102  }
103 
104  // Get MCTrack array
105  fMCTrackArray = (TClonesArray*) ioman->GetObject("MCTrack");
106  if( !fMCTrackArray ) {
107  cout << "-I- "<< GetName() <<"::Init: No MCTrack array!"
108  << endl;
109  // return;
110  }
111 
112  // Get PndGemPoint (MCPoint) array
113  fMCPointArray = (TClonesArray*) ioman->GetObject("GEMPoint");
114  if( !fMCPointArray ) {
115  cout << "-I- "<< GetName() <<"::Init: No MCPoint array!"
116  << endl;
117  // return;
118  }
119 
120  fDigiPar = (PndGemDigiPar*)(rtdb->getContainer("PndGemDetectors"));
121  cout << "-I- " << "PndGemTrackFinderOnHitsTB" << "::Init(). There are " << fDigiPar->GetNStations() << " GEM stations." << endl;
122  cout << "-I- " << "PndGemTrackFinderOnHitsTB" << "::Init(). Initialization succesfull." << endl;
123  fParThetaA = fDigiPar->GetTrackFinderOnHits_ParThetaA();
124  fParThetaB = fDigiPar->GetTrackFinderOnHits_ParThetaB();
125 
126  fParTheta0 = fDigiPar->GetTrackFinderOnHits_ParTheta0();
127  fParTheta1 = fDigiPar->GetTrackFinderOnHits_ParTheta1();
128  fParTheta2 = fDigiPar->GetTrackFinderOnHits_ParTheta2();
129  fParTheta3 = fDigiPar->GetTrackFinderOnHits_ParTheta3();
130 
131  fParRadPhi0 = fDigiPar->GetTrackFinderOnHits_ParRadPhi0();
132  fParRadPhi2 = fDigiPar->GetTrackFinderOnHits_ParRadPhi2();
133  for ( Int_t in = 0 ; in < 3 ; in++ ) {
134  fParMat0[in] = fDigiPar->GetTrackFinderOnHits_ParMat0(in);
135  fParMat1[in] = fDigiPar->GetTrackFinderOnHits_ParMat1(in);
136  }
137 
138  std::cout << "-I- "<< GetName() <<": Intialization successfull" << std::endl;
139 }
140 
141 // ----- Private method SetParContainers -------------------------------
142 void PndGemTrackFinderOnHitsTB::SetParContainers() {
143 
144  // Get run and runtime database
145  FairRunAna* run = FairRunAna::Instance();
146  if ( ! run ) Fatal("SetParContainers", "No analysis run");
147 
148  FairRuntimeDb* db = run->GetRuntimeDb();
149  if ( ! db ) Fatal("SetParContainers", "No runtime database");
150 
151  // Get GEM digitisation parameter container
152  fDigiPar = (PndGemDigiPar*)(db->getContainer("PndGemDetectors"));
153 
154 }
155 // -------------------------------------------------------------------------
156 
157 // ----- Public method DoFind ------------------------------------------
158 Int_t PndGemTrackFinderOnHitsTB::DoFind(TClonesArray* hitArray,
159  TClonesArray* trackArray,
160  TClonesArray* trackCandArray) {
161  // cout << "======== PndGemTrackFinderOnHitsTB::Exec(Event = " << fNofEvents << " ) ====================" << endl;
162 
163  fTimer.Start();
164 
165  // Get GEM digitisation parameter container
166  fTrackSegments.clear();
167 
168  // Count events
169  fNofEvents++;
170 
171  if ( fVerbose ) {
172  cout << endl << endl<< endl << endl;
173  cout << "=======================================================" << endl;
174  cout << "-I- Event No: " << fNofEvents << endl;
175  cout << "=======================================================" << endl;
176 
177 
178  cout <<"-I- "<< GetName() <<"::DoFind "<< endl;
179  cout << "-------------------------------------------------------" << endl;
180  cout << " ### Start DoFind" << endl;
181  cout << "-------------------------------------------------------" << endl;
182  }
183 
184  // Check pointers
185  fMCAvailable = kTRUE;
186  if( !fMCTrackArray || !fMCPointArray ) {
187  fMCAvailable = kFALSE;
188  }
189 
190  if( !hitArray ) {
191  cout << "-E- "<< GetName() <<"::DoFind: "
192  << "Hit arrays missing! "<< endl;
193  return -1;
194  }
195 
196  // Initialise control counters
197  //Int_t nNoTrack = 0; //[R.K. 01/2017] unused variable?
198  //Int_t nNoGemPoint = 0; //[R.K. 01/2017] unused variable?
199  //Int_t nNoGemHit = 0; //[R.K. 01/2017] unused variable?
200 
201  // Create pointers to GemHit and GemPoint
202  //PndGemHit* gemHit = NULL; //[R.K. 01/2017] unused variable?
203  //PndGemHit* gemHit2 = NULL; //[R.K. 01/2017] unused variable?
204  //PndTrackCand* gemTrackCand = NULL; //[R.K. 01/2017] unused variable?
205 
206  // Declare variables outside the loop
207  //Int_t trackIndex = 0; // Gem track index //[R.K. 01/2017] unused variable?
208  //Int_t relDetID = -1;//3000; // //[R.K. 01/2017] unused variable?
209 
210  if ( fMCAvailable & fVerbose ) {
211  cout <<"# MC Tracks: "<< fMCTrackArray->GetEntriesFast() << endl;
212  cout <<"# MC Points: "<< fMCPointArray->GetEntriesFast() << endl;
213  }
214 
215  // Number of Gem hits
216  Int_t nGemHits = hitArray->GetEntriesFast();
217  if(fVerbose) cout <<"# GemHits: "<< nGemHits << endl;
218 
219  fPrepTime+=fTimer.RealTime();
220  fTimer.Continue();
221 
222  for ( Int_t istat1 = 0 ; istat1 < fDigiPar->GetNStations() ; istat1++ ) {
223  for ( Int_t istat2 = istat1+1 ; istat2 < fDigiPar->GetNStations() ; istat2++ ) {
224  FindTrackSegments(hitArray,istat1,istat2);
225  }
226  }
227 
228  fSegmTime+=fTimer.RealTime();
229  fTimer.Continue();
230 
231  if ( fVerbose ) {
232  cout << "event " << fNofEvents << " >>> " << fTrackSegments.size() << " track segments. " << endl;
233  if ( fMCAvailable )
234  PrintMCTrackSegments(hitArray);
235  else
236  PrintTrackSegments(hitArray);
237  }
238 
239  Int_t nr = MatchTrackSegments();
240 
241  fMatchTime+=fTimer.RealTime();
242  fTimer.Continue();
243 
244  RemoveCloneTracks(nr);
245 
246  fRemoveTime+=fTimer.RealTime();
247  fTimer.Continue();
248 
249  if ( fVerbose ) {
250  cout << "************************************************" << endl;
251  if ( fMCAvailable )
252  PrintMCTracks(hitArray,nr);
253  else
254  PrintTracks(hitArray,nr);
255  cout << "************************************************" << endl;
256  cout << "finished printing tracks" << endl;
257  }
258 
259  nr = CreateTracks(hitArray,trackArray,trackCandArray,nr);
260 
261  fWriteTime+=fTimer.RealTime();
262  fTimer.Continue();
263 
264  if ( fVerbose ) {
265  cout << "------------------------------------------------" << endl;
266  cout << "!!!!!!!!!!!!!!!!!! " << nr << " tracks have been found" << endl;
267  cout << "------------------------------------------------" << endl;
268  }
269 
270  // cout << " PndGemTrackFinderOnHitsTB::DoFind. Found " << nr << " tracks out of " << hitArray->GetEntries() << " hits." << endl;
271 
272  fAllTime+=fTimer.RealTime();
273 
274  fTimer.Stop();
275 
276  return nr;
277 }
278 // ------------------------------------------------------------
279 
280 // --- Private method to create tracks ------------------------
281 Int_t PndGemTrackFinderOnHitsTB::CreateTracks(TClonesArray* hitArray,
282  TClonesArray* trackArray,
283  TClonesArray* trackCandArray,
284  Int_t nofRecoTracks) {
285 
286  Int_t nofCreatedTracks = 0;
287 
288  const Int_t kNofRecoTracks = nofRecoTracks;
289  const Int_t kNofStatDbl = 2*fDigiPar->GetNStations();
290 
291  Int_t hitIndices[kNofRecoTracks][kNofStatDbl*10]; // RKCHANGE
292  Int_t nofHits[kNofRecoTracks];
293 
294  Double_t meanMom[kNofRecoTracks];
295  Double_t meanPhi[kNofRecoTracks];
296  Double_t meanThe[kNofRecoTracks];
297  Int_t nofTS[kNofRecoTracks];
298 
299  PndTrackCandHit tcHit;
300  PndTrackCand* gemTrackCand;
301  PndGemHit* gemHit;
302  //PndGemHit* gemHit2; //[R.K. 01/2017] unused variable?
303 
304  for ( Int_t itr = 0 ; itr < nofRecoTracks ; itr++ ) {
305  nofHits[itr] = 0;
306  meanMom[itr] = 0.;
307  meanPhi[itr] = 0.;
308  meanThe[itr] = 0.;
309  nofTS[itr] = 0;
310 
311  for ( size_t its = 0 ; its < fTrackSegments.size() ; its++ ) {
312  TrackSegmentTB iterTS = fTrackSegments[its];
313  if ( iterTS.recoTrackIndex != itr ) continue;
314 
315  for ( Int_t ihit = 0 ; ihit < 2 ; ihit++ ) {
316  Bool_t hitAdded = kFALSE;
317  for ( Int_t itrh = 0 ; itrh < nofHits[itr] ; itrh++ ) {
318  if ( hitIndices[itr][itrh] == iterTS.hitIndex[ihit] )
319  hitAdded = kTRUE;
320  }
321  if ( hitAdded ) continue;
322  hitIndices[itr][nofHits[itr]] = iterTS.hitIndex[ihit];
323  nofHits[itr]+=1;
324  }
325 
326  if ( meanPhi[itr]/(nofTS[itr]+1) < 5. && iterTS.trackPhi > 355. ) { iterTS.trackPhi -= 360.; }
327  if ( meanPhi[itr]/(nofTS[itr]+1) > 355. && iterTS.trackPhi < 5. ) { iterTS.trackPhi += 360.; }
328  meanMom[itr] += iterTS.trackMom;
329  meanPhi[itr] += iterTS.trackPhi;
330  meanThe[itr] += iterTS.trackTheta;
331  nofTS[itr] ++;
332  }
333 
334  if ( nofTS[itr] == 0 ) continue;
335 
336  meanMom[itr] = meanMom[itr]/nofTS[itr];
337  meanPhi[itr] = meanPhi[itr]/nofTS[itr];
338  meanThe[itr] = meanThe[itr]/nofTS[itr];
339 
340  gemTrackCand = new((*trackCandArray)[nofCreatedTracks]) PndTrackCand();
341 
342  // cout << "TRACK----------------------------------" << endl;
343 
344  Double_t sumZ = 0.;
345  Double_t sumZSq = 0.;
346  Double_t sumZT = 0.;
347  Double_t sumT = 0.;
348  for ( Int_t itrh = 0 ; itrh < nofHits[itr] ; itrh++ ) {
349  gemHit = (PndGemHit*)hitArray->At(hitIndices[itr][itrh]);
350 
351  // cout << "gemHit " << ih << " AT z = " << gemHit->GetZ() << " cm and t = " << gemHit->GetTimeStamp() << " ns" << endl;
352  gemTrackCand->AddHit(FairRootManager::Instance()->GetBranchId("GEMHit"),
353  hitIndices[itr][itrh],gemHit->GetPosition().Mag());
354 
355  if ( fVerbose )
356  cout << "hit " << itrh << " @ "
357  << gemHit->GetX() << " "
358  << gemHit->GetY() << " "
359  << gemHit->GetZ() << " was created at " << gemHit->GetTimeStamp() << endl;
360  sumZ += gemHit->GetZ();
361  sumZSq += gemHit->GetZ()*gemHit->GetZ();
362  sumZT += gemHit->GetZ()*gemHit->GetTimeStamp();
363  sumT += gemHit->GetTimeStamp();
364 
365  }
366  Double_t time0 = (sumT*sumZSq-sumZ*sumZT)/(nofHits[itr]*sumZSq-sumZ*sumZ);
367  if ( fVerbose )
368  cout << "----> time0 = " << time0 << endl;
369 
370  gemTrackCand->Sort();
371 
372  gemTrackCand->SetTimeStamp(time0);
373 
374  TVector3 pos (0.,0.,0.);
375  TVector3 mom;
376  mom.SetMagThetaPhi(TMath::Abs(meanMom[itr]),meanThe[itr]*TMath::DegToRad(),meanPhi[itr]*TMath::DegToRad());
377 
378  // cout << "TRACK "
379  // << "( " << meanMom[itr] << " , " << meanThe[itr] << " , " << meanPhi[itr] << " ) >>> "
380  // << "( " << mom.Mag() << " , " << mom.Theta()*TMath::RadToDeg() << " , " << mom.Phi()*TMath::RadToDeg() << " )" << endl;
381 
382  Int_t charge = -1;
383  // if ( mom.Mag() < 0. )
384  if ( meanThe[itr] < 0. ) {
385  // cout << "THE CHARGE IS " << 1 << endl;
386  charge = 1;
387  mom.SetX(-mom.X());
388  mom.SetY(-mom.Y());
389  }
390 
391  FairTrackParP firstPar(pos,mom,
392  TVector3(0.5, 0.5, 0.5),
393  0.1*mom,
394  charge,
395  pos,
396  TVector3(1.,0.,0.),
397  TVector3(0.,1.,0.));
398  FairTrackParP lastPar (pos,mom,
399  TVector3(0.5,0.5,0.5),//tempMom[0],tempMom[1],tempMom[2]),
400  0.1*mom,
401  charge,
402  pos,
403  TVector3(1.,0.,0.),
404  TVector3(0.,1.,0.));
405  // cout << "q = " << firstPar->GetQ() << endl;
406 // TClonesArray &pndtracks = *trackArray;
407 // Int_t size = pndtracks.GetEntriesFast();
408 // PndTrack* pndTrack = new(pndtracks[size]) PndTrack(*firstPar, *lastPar, *gemTrackCand);
409  new((*trackArray)[nofCreatedTracks]) PndTrack(firstPar, lastPar, *gemTrackCand);
410 
411  /* PndTrack* checkTrack = (PndTrack*) trackArray->At(nofCreatedTracks);
412  for ( Int_t ih = 0 ; ih < kNofStatDbl ; ih++ ) {
413  if ( hitIndices[itr][ih] == -1 ) continue;
414  checkTrack->AddLink(FairLink("GemHit", hitIndices[itr][ih]));
415  }*/
416 
417 // cout << "now q = " << checkTrack->GetParamFirst().GetQ() << endl;
418 
419  nofCreatedTracks++;
420  }
421 
422  return nofCreatedTracks;
423 }
424 // ------------------------------------------------------------
425 
426 // --- Private method to remove clone track -------------------
427 void PndGemTrackFinderOnHitsTB::RemoveCloneTracks(Int_t nofRecoTracks) {
428  Bool_t printInfo = kFALSE;//TRUE;
429 
430  if ( fVerbose > 4 || printInfo )
431  cout << "Trying to remove clone tracks" << endl;
432 
433  const Int_t kNofRecoTracks = nofRecoTracks;
434  const Int_t kNofStatDbl = 2*fDigiPar->GetNStations();
435 
436  Int_t hitIndices[kNofRecoTracks][kNofStatDbl];
437  Int_t nofHits[kNofRecoTracks];
438  Int_t nofMultiHits[kNofRecoTracks];
439 
440  Double_t meanMom[kNofRecoTracks];
441  Double_t meanPhi[kNofRecoTracks];
442  Double_t meanThe[kNofRecoTracks];
443 
444  Int_t nofComb = 0;
445  for ( Int_t i1 = 0 ; i1 < fDigiPar->GetNStations() ; i1++ ) for ( Int_t i2 = i1+1 ; i2 < fDigiPar->GetNStations() ; i2++ ) nofComb++;
446  const Int_t kNofComb = nofComb;
447 
448  Int_t nofTS[kNofRecoTracks];
449  Double_t valMom[kNofRecoTracks][kNofComb];
450  Double_t valPhi[kNofRecoTracks][kNofComb];
451  Double_t valThe[kNofRecoTracks][kNofComb];
452  Double_t trackCons[kNofRecoTracks];
453 
454  for ( Int_t itr = 0 ; itr < nofRecoTracks ; itr++ ) {
455  for ( Int_t ih = 0 ; ih < kNofStatDbl ; ih++ ) hitIndices[itr][ih] = -1;
456  nofHits[itr] = 0;
457  nofMultiHits[itr] = 0;
458  meanMom[itr] = 0.;
459  meanPhi[itr] = 0.;
460  meanThe[itr] = 0.;
461  nofTS[itr] = 0;
462  for ( Int_t ic = 0 ; ic < kNofComb ; ic++ ) {
463  valMom[itr][ic] = -666.;
464  valPhi[itr][ic] = -666.;
465  valThe[itr][ic] = -666.;
466  }
467  trackCons[itr] = 0.;
468 
469  for ( size_t its = 0 ; its < fTrackSegments.size() ; its++ ) {
470  TrackSegmentTB iterTS = fTrackSegments[its];
471  if ( iterTS.recoTrackIndex != itr ) continue;
472  for ( Int_t istat = 0 ; istat < 2 ; istat++ ) {
473  hitIndices[itr][2*iterTS.stationIndex[istat] ] = iterTS.hitIndex[2*istat ];
474  hitIndices[itr][2*iterTS.stationIndex[istat]+1] = iterTS.hitIndex[2*istat+1];
475  }
476  if ( meanPhi[itr]/(nofTS[itr]+1) < 5. && iterTS.trackPhi > 355. ) { iterTS.trackPhi -= 360.; }
477  if ( meanPhi[itr]/(nofTS[itr]+1) > 355. && iterTS.trackPhi < 5. ) { iterTS.trackPhi += 360.; }
478 
479  valMom[itr][nofTS[itr]] = iterTS.trackMom;
480  valPhi[itr][nofTS[itr]] = iterTS.trackPhi;
481  valThe[itr][nofTS[itr]] = iterTS.trackTheta;
482  meanMom[itr] += iterTS.trackMom;
483  meanPhi[itr] += iterTS.trackPhi;
484  meanThe[itr] += iterTS.trackTheta;
485  nofTS[itr]++;
486  }
487  for ( Int_t isens = 0 ; isens < kNofStatDbl ; isens++ )
488  if ( hitIndices[itr][isens] > -0.5 )
489  nofHits[itr] += 1;
490 
491  meanMom[itr] = meanMom[itr]/nofTS[itr];
492  meanPhi[itr] = meanPhi[itr]/nofTS[itr];
493  meanThe[itr] = meanThe[itr]/nofTS[itr];
494 
495  if ( fVerbose > 4 || printInfo )
496  cout << " MEAN = " << meanMom[itr] << " " << meanPhi[itr] << " " << meanThe[itr] << endl;
497  for ( Int_t its = 0 ; its < nofTS[itr] ; its++ ) {
498  if ( fVerbose > 4 || printInfo )
499  cout << " v" << its << " = " << valMom[itr][its] << " " << valPhi[itr][its] << " " << valThe[itr][its] << endl;
500  trackCons[itr] += TMath::Abs((valMom[itr][its]-meanMom[itr])/meanMom[itr]);
501  trackCons[itr] += TMath::Abs( valPhi[itr][its]-meanPhi[itr]);
502  trackCons[itr] += TMath::Abs( valThe[itr][its]-meanThe[itr]);
503  }
504  trackCons[itr] = trackCons[itr]/(3.*nofTS[itr]);
505  }
506 
507  for ( Int_t itr = 0 ; itr < nofRecoTracks ; itr++ ) {
508  for ( Int_t ih = 0 ; ih < kNofStatDbl ; ih++ ) {
509  if ( hitIndices[itr][ih] == -1 ) continue;
510  Bool_t hitFound = kFALSE;
511  for ( Int_t itr2 = 0 ; itr2 < nofRecoTracks ; itr2++ ) {
512  if ( itr == itr2 ) continue;
513  for ( Int_t ih2 = 0 ; ih2 < kNofStatDbl ; ih2++ ) {
514  if ( hitIndices[itr2][ih2] == -1 ) continue;
515  if ( hitIndices[itr2][ih2] == hitIndices[itr][ih] ) {
516  hitFound = kTRUE;
517  nofMultiHits[itr]++;
518  break;
519  }
520  }
521  if ( hitFound ) break;
522  }
523  }
524 
525  Int_t cloneIndicators = 0;
526  if ( trackCons[itr] > 1 ) cloneIndicators++; // TS's too different
527  if ( nofTS[itr] <= kNofComb*2/3 ) cloneIndicators++; // not enough TS
528  if ( nofMultiHits[itr] >= nofHits[itr]/3 ) cloneIndicators++; // too many hits shared with other tracks
529 
530  if ( cloneIndicators >= 2 ) { // remove this track candidate
531  for ( size_t its = 0 ; its < fTrackSegments.size() ; its++ ) {
532  TrackSegmentTB iterTS = fTrackSegments[its];
533  if ( iterTS.recoTrackIndex != itr ) continue;
534  iterTS.recoTrackIndex = -1;
535  fTrackSegments[its] = iterTS;
536  }
537  }
538  if ( fVerbose > 4 || printInfo ) {
539  cout << " consistency = " << trackCons[itr] << ( trackCons[itr] > 1 ? " YES":"") << " ( > 1 )" << endl;
540  cout << " segments = " << nofTS[itr] << ( nofTS[itr] <= kNofComb*2/3 ? " YES":"") << " ( <= " << kNofComb*2/3 << " )" << endl;
541  cout << " multihits = " << nofMultiHits[itr] << ( nofMultiHits[itr] >= nofHits[itr]/3 ? " YES":"") << " ( >= " << nofHits[itr]/3 << " )" << endl;
542  cout << "TRACK " << itr << " HAS " << nofMultiHits[itr] << " MULTI HITS AND CONSISTENCY = " << trackCons[itr] << (cloneIndicators>=2?" >>> REMOVED!!!! ":"") << endl;
543  }
544  }
545 }
546 // ------------------------------------------------------------
547 
548 // --- Private method to match track segments -----------------
549 Int_t PndGemTrackFinderOnHitsTB::MatchTrackSegments() {
550  // const Int_t kMaxNofSegments = fDigiPar->GetNStations()-1;
551 
552  Bool_t printInfo = kFALSE;//TRUE;
553 
554  Int_t nofRecoTracks = 0;
555 
556  for ( size_t itrc = 0 ; itrc < fTrackSegments.size() ; itrc++ ) {
557  TrackSegmentTB origTS = fTrackSegments[itrc];
558  if ( origTS.recoTrackIndex != -1 ) continue; // this track segment has already been used
559  if ( fVerbose > 4 || printInfo )
560  cout << "MATCH TRACK SEGMENTS TO SEGM " << itrc << " with params: " << origTS.trackMom << " " << origTS.trackPhi << " " << origTS.trackTheta << endl;
561  // if ( origTS.stationIndex[1] != origTS.stationIndex[0]+1 ) continue; // only consider segments of hits on neighbouring stations
562 
563  vector<Int_t> trackSegs;
564  trackSegs.push_back(itrc);
565 
566  Bool_t overrepr = kFALSE;
567  for ( size_t itrc2 = 0 ; itrc2 < fTrackSegments.size() ; itrc2++ ) {
568  if ( itrc2 == itrc ) continue; // do not match segment with itself
569  TrackSegmentTB matchTS = fTrackSegments[itrc2];
570  if ( origTS.stationIndex[0] == matchTS.stationIndex[0] && origTS.stationIndex[1] == matchTS.stationIndex[1] &&
571  origTS.sensorNumber[0] == matchTS.sensorNumber[0] && origTS.sensorNumber[1] == matchTS.sensorNumber[1] ) continue;
572  if ( matchTS.recoTrackIndex != -1 ) continue; // this track segment has already been used
573  // if ( matchTS.stationIndex[1] != matchTS.stationIndex[0]+1 ) continue; // only consider segments of hits on neighbouring stations
574 
575  // check if track matches with any of the track segments in trackSegs
576  Double_t meanMom = 0.;
577  Double_t meanPhi = 0.;
578  Double_t meanThe = 0.;
579  Bool_t matching = kFALSE;
580  for ( size_t its = 0 ; its < trackSegs.size() ; its++ ) {
581  TrackSegmentTB iterTS = fTrackSegments[trackSegs[its]];
582  if ( !matching ) {
583  if ( matchTS.hitIndex[0] == iterTS.hitIndex[0] ) { matching = kTRUE; } //break; } // matchTS matches after iterTS
584  if ( matchTS.hitIndex[1] == iterTS.hitIndex[1] ) { matching = kTRUE; } //break; } // matchTS matches before iterTS
585  if ( matchTS.hitIndex[0] == iterTS.hitIndex[1] ) { matching = kTRUE; } //break; } // matchTS shares hit on first station with iterTS
586  if ( matchTS.hitIndex[1] == iterTS.hitIndex[0] ) { matching = kTRUE; } //break; } // matchTS shares hit on second station with iterTS
587  }
588  meanMom += iterTS.trackMom/(trackSegs.size());
589  meanPhi += iterTS.trackPhi/(trackSegs.size());
590  meanThe += iterTS.trackTheta/(trackSegs.size());
591  }
592  if ( !matching ) continue; // matchTS does not match to any track segment in trackSegs
593 
594  if ( meanPhi < 5. && matchTS.trackPhi > 355. ) matchTS.trackPhi-=360.;
595  if ( meanPhi > 355. && matchTS.trackPhi < 5. ) matchTS.trackPhi+=360.;
596 
597  if ( TMath::Abs(matchTS.trackPhi-meanPhi) > 5. ||
598  TMath::Abs(matchTS.trackTheta-meanThe) > 5. ||
599  TMath::Abs(matchTS.trackMom-meanMom) > 0.3*TMath::Abs(meanMom) ) {
600  if ( fVerbose > 4 || printInfo ) {
601  cout << " mean > " << meanMom << " " << meanPhi << " " << meanThe << endl;
602  cout << " ++++ > " << matchTS.trackMom << " " << matchTS.trackPhi << " " << matchTS.trackTheta << endl;
603  }
604  continue;
605  }
606 
607  // check if there exists in trackSegs a segment that uses different hits on some station
608  Bool_t mismatch = kFALSE;
609  Bool_t mismatchSolved = kFALSE;
610  for ( size_t its = 0 ; its < trackSegs.size() ; its++ ) {
611  TrackSegmentTB iterTS = fTrackSegments[trackSegs[its]];
612  if ( matchTS.stationIndex[0] == iterTS.stationIndex[0] &&
613  matchTS.stationIndex[1] == iterTS.stationIndex[1] &&
614  matchTS.sensorNumber[0] == iterTS.sensorNumber[0] &&
615  matchTS.sensorNumber[1] == iterTS.sensorNumber[1] ) {
616 
617  if ( fVerbose > 4 || printInfo )
618  cout << "there are already " << trackSegs.size() << " segments in this track: " << endl;
619  meanMom = 0.;
620  meanPhi = 0.;
621  meanThe = 0.;
622  for ( size_t its2 = 0 ; its2 < trackSegs.size() ; its2++ ) {
623  if ( its == its2 ) continue; // the track segment in question should not go to mean
624  TrackSegmentTB aveTS = fTrackSegments[trackSegs[its2]];
625  if ( fVerbose > 4 || printInfo )
626  cout << its2 << " (" << trackSegs[its2] << ") > " << aveTS.trackMom << " " << aveTS.trackPhi << " " << aveTS.trackTheta << endl;
627  meanMom += aveTS.trackMom/(trackSegs.size()-1);
628  meanPhi += aveTS.trackPhi/(trackSegs.size()-1);
629  meanThe += aveTS.trackTheta/(trackSegs.size()-1);
630  }
631  if ( fVerbose > 4 || printInfo ) {
632  cout << its << " (" << trackSegs[its] << ") > " << iterTS.trackMom << " " << iterTS.trackPhi << " " << iterTS.trackTheta << endl;
633  cout << " this track conflicts with track " << its << " and is: " << endl;
634  cout << "C (" << itrc2 << ") > " << matchTS.trackMom << " " << matchTS.trackPhi << " " << matchTS.trackTheta << endl;
635  cout << " should decide on mean of other track segments, which is: " << endl;
636  cout << "MEAN > " << meanMom << " " << meanPhi << " " << meanThe << endl;
637  }
638  if ( meanPhi < 5. ) {
639  if ( iterTS.trackPhi > 355. ) iterTS.trackPhi-=360.;
640  if ( matchTS.trackPhi > 355. ) matchTS.trackPhi-=360.;
641  }
642  if ( meanPhi > 355 ) {
643  if ( iterTS.trackPhi < 5. ) iterTS.trackPhi+=360.;
644  if ( matchTS.trackPhi < 5. ) matchTS.trackPhi+=360.;
645  }
646  Bool_t changeSegm = kFALSE;
647  // check if new segment is better than the old
648  if ( TMath::Abs(matchTS.trackPhi-meanPhi) < TMath::Abs(iterTS.trackPhi-meanPhi) ) // && TMath::Abs(matchTS.trackPhi-meanPhi) < 5. )
649  if ( TMath::Abs(matchTS.trackTheta-meanThe) < TMath::Abs(iterTS.trackTheta-meanThe) ) // && TMath::Abs(matchTS.trackTheta-meanThe) < 5. )
650  if ( TMath::Abs(matchTS.trackMom-meanMom) < TMath::Abs(iterTS.trackMom-meanMom) ) // && 0.3*TMath::Abs(meanMom) )
651  changeSegm = kTRUE;
652  // new segment is better, replace it
653  if ( changeSegm ) {
654  if ( fVerbose )
655  cout << "Changing segment " << trackSegs[its] << " with segment " << itrc2 << endl;
656  trackSegs[its] = itrc2;
657  mismatchSolved = kTRUE;
658  continue;
659  }
660  /* if ( TMath::Abs(iterTS.trackPhi-meanPhi) < 5. )
661  if ( TMath::Abs(iterTS.trackTheta-meanThe) < 5. )
662  if ( TMath::Abs(iterTS.trackMom-meanMom) < 0.3*TMath::Abs(meanMom) ) {
663  }
664  mismatch = kTRUE;
665  break;
666  */
667  mismatchSolved = kTRUE;
668  continue;
669  }
670  }
671  if ( mismatch ) {
672  if ( fVerbose > 4 || printInfo )
673  cout << " PROBLEM HERE " << endl;
674  overrepr = kTRUE;
675  break;
676  }
677  if ( mismatchSolved ) { continue; }
678  // matchTS matches with some of trackSegs, and it does not conflict with any of trackSegs
679 
680  trackSegs.push_back(itrc2);
681  }
682  if ( overrepr ) continue;
683  if ( trackSegs.size() == 1 ) continue;
684 
685  if ( fVerbose > 4 || printInfo )
686  cout << "segments: " << endl;
687  for ( size_t its = 0 ; its < trackSegs.size() ; its++ ) {
688  TrackSegmentTB iterTS = fTrackSegments[trackSegs[its]];
689  if ( fVerbose > 4 || printInfo ) {
690  cout << iterTS.stationIndex[0] << "." << iterTS.sensorNumber[0] << " "
691  << iterTS.stationIndex[1] << "." << iterTS.sensorNumber[1] << " > " << flush;
692  for ( Int_t ihit = 0 ; ihit < 2 ; ihit++ )
693  cout << setw(4) << iterTS.hitIndex[ihit] << " " << flush;
694  cout << setw(11) << iterTS.trackMom << " " << setw(11) << iterTS.trackPhi << " " << setw(11) << iterTS.trackTheta << endl;
695  }
696  iterTS.recoTrackIndex = nofRecoTracks;
697  fTrackSegments[trackSegs[its]] = iterTS;
698  }
699  if ( fVerbose > 4 || printInfo )
700  cout << " seems to belong to one track" << endl;
701  nofRecoTracks++;
702  }
703 
704  if ( fVerbose || printInfo )
705  cout << "********* " << nofRecoTracks << " track candidates have been found" << endl;
706 
707  return nofRecoTracks;
708 }
709 // ------------------------------------------------------------
710 
711 // ---- Private function to find track segments on 2 stations --------------------
712 Int_t PndGemTrackFinderOnHitsTB::FindTrackSegments(TClonesArray* hitArray, Int_t stat1Id, Int_t stat2Id) {
713  Bool_t printInfo = kFALSE;//TRUE;
714 
715  PndGemHit* gemHit1;
716  PndGemHit* gemHit2;
717  Double_t gemHit1X = 0., gemHit1Y = 0., gemHit1Z = 0., gemHit1T = 0.;
718  Double_t gemHit2X = 0., gemHit2Y = 0., gemHit2Z = 0., gemHit2T = 0.;
719 
720  Int_t nGemHits = hitArray->GetEntriesFast();
721 
722  for(Int_t iHit1 = 0; iHit1 < nGemHits; iHit1++ ) {
723  // Get the pointer to Gem hit
724  gemHit1 = (PndGemHit*) hitArray->At(iHit1);
725  if ( gemHit1->GetCharge() < 1. ) continue; //RKChange
726  if ( gemHit1->GetStationNr()-1 != stat1Id ) continue;
727  if ( fVerbose > 3 || printInfo )
728  cout << "LOOKING FOR TRACK SEGMENTS STARTING AT " << gemHit1->GetX() << " " << gemHit1->GetY() << " " << gemHit1->GetZ() << endl;
729 
730  gemHit1X = gemHit1->GetX();
731  gemHit1Y = gemHit1->GetY();
732  gemHit1Z = gemHit1->GetZ();
733  gemHit1T = gemHit1->GetTimeStamp();
734 
735  if ( fVerbose > 3 || printInfo )
736  cout << "possible segment " << fTrackSegments.size() << " starts here " << gemHit1X << " " << gemHit1Y << " " << gemHit1Z << endl;
737  Double_t radius = TMath::Sqrt(gemHit1X*gemHit1X+gemHit1Y*gemHit1Y);
738  Double_t pangle = TMath::ACos(gemHit1X/radius);
739  if ( gemHit1Y < 0 )
740  pangle = 2.*TMath::Pi() - pangle;
741  Double_t trackTheta = fParThetaA * radius / gemHit1Z + fParThetaB;
742  if ( fVerbose > 3 || printInfo )
743  cout << " -> with theta of " << trackTheta << " (radius = " << radius << " and phi angle = " << pangle*TMath::RadToDeg() << ")" << endl;
744 
745  for(Int_t iHit2 = 0; iHit2 < nGemHits; iHit2++){
746  gemHit2 = (PndGemHit*) hitArray->At(iHit2);
747  if ( gemHit2->GetCharge() < 1. ) continue; //RKChange
748  if ( gemHit2->GetStationNr()-1 != stat2Id ) continue;
749 
750  gemHit2X = gemHit2->GetX();
751  gemHit2Y = gemHit2->GetY();
752  gemHit2Z = gemHit2->GetZ();
753  gemHit2T = gemHit2->GetTimeStamp();
754 
755  if ( gemHit2T < gemHit1T || gemHit2T > gemHit1T+11. ) continue; // the hits are too far away in time or hit2 is before hit1
756 
757  Double_t radius2 = TMath::Sqrt(gemHit2X*gemHit2X+gemHit2Y*gemHit2Y);
758  Double_t pangle2 = TMath::ACos(gemHit2X/radius2);
759  if ( gemHit2Y < 0 )
760  pangle2 = 2.*TMath::Pi() - pangle2;
761 
762  if ( pangle < TMath::Pi()/2. && pangle2 > TMath::Pi()*3./2. ) pangle2 -= TMath::Pi()*2.;
763  if ( pangle > TMath::Pi()*3./2. && pangle2 < TMath::Pi()/2. ) pangle2 += TMath::Pi()*2.;
764 
765  if ( TMath::Abs(pangle-pangle2)*TMath::RadToDeg() > 40 ) continue;
766 
767  if ( fVerbose > 3 || printInfo ) {
768  cout << " trying to match it with " << gemHit2X << " , " << gemHit2Y << " , " << gemHit2Z << endl;
769  cout << " (radius = " << radius2 << " and phi angle = " << pangle2*TMath::RadToDeg() << ")" << endl;
770  }
771 
772  Double_t dphi = pangle2-pangle;
773  Double_t expectedRad2 = radius / ( TMath::Cos(dphi) - TMath::Sin(dphi) / TMath::Tan(gemHit2Z/(gemHit2Z-gemHit1Z)*dphi) );
774 
775  Double_t zDistRatio = gemHit2Z/gemHit1Z;
776  Double_t zDiffRatio = gemHit1Z/(gemHit2Z-gemHit1Z);
777  Double_t zCuDiff = gemHit2Z*gemHit2Z*gemHit2Z-gemHit1Z*gemHit1Z*gemHit1Z;
778  Double_t zSqDiff = gemHit2Z*gemHit2Z-gemHit1Z*gemHit1Z;
779  Double_t zDiff = gemHit2Z-gemHit1Z;
780  Double_t par0_mom = fParMat0[0]*zCuDiff+fParMat0[1]*zSqDiff+fParMat0[2]*zDiff;
781  Double_t par1_mom = fParMat1[0]*zCuDiff+fParMat1[1]*zSqDiff+fParMat1[2]*zDiff;
782 
783  Double_t expRadUncert = 0.08*radius*zDistRatio; // changed 0.05 to 0.08
784 
785  if ( fVerbose > 3 || printInfo )
786  cout << " -> while expected radius was " << expectedRad2 << " with error of " << expRadUncert << endl;
787 
788  if ( radius2>expectedRad2+expRadUncert || radius2<expectedRad2-expRadUncert ) continue;
789 
790  if ( fVerbose > 4 || printInfo )
791  cout << "STRONG CORRELATION FOR THIS HIT!!!" << endl;
792  // calculate phi and momentum basing on the pangle-pangle2;
793  Double_t trackMomentum = 666.;
794  if ( (pangle-pangle2) != 0 )
795  trackMomentum = (par0_mom+par1_mom/radius)/((pangle-pangle2)*TMath::RadToDeg());;
796 
797  if ( fVerbose > 0 && stat1Id == 0 && stat2Id == 1 ) {
798  cout << "CALCULATED MOMENTUM " << trackMomentum << " FROM ( " << par0_mom << " + " << par1_mom << " / " << radius << " ) / ( " << pangle << " - " << pangle2 << " ) * " << TMath::RadToDeg() << " = ( " << par0_mom+par1_mom/radius << " ) / ( " << pangle-pangle2 << " ) * " << TMath::RadToDeg() << endl;
799  cout << pangle-pangle2 << " from " << TMath::ACos(gemHit1X/radius) << " - " << TMath::ACos(gemHit2X/radius2) << " = ACos( " << gemHit1X << " / " << radius << ") - ACos( " << gemHit2X << " / " << radius2 << " ) " << endl;
800  cout << " " << gemHit1X << " " << gemHit2X << endl;
801  cout << " " << gemHit1Y << " " << gemHit2Y << endl;
802  cout << " " << gemHit1Z << " " << gemHit2Z << endl;
803  }
804  Double_t circRad = TMath::Sqrt( (gemHit2X-gemHit1X)*(gemHit2X-gemHit1X) + (gemHit2Y-gemHit1Y)*(gemHit2Y-gemHit1Y) ) / TMath::Sin(pangle-pangle2) / 2.;
805 
806  Double_t mom_trans = circRad / 196.;
807 
808  Double_t newMomPar = 4.585; // value for stations 2-3
809  if ( TMath::Abs(gemHit1Z-117.) < 5. ) { // station 1
810  if ( TMath::Abs(gemHit2Z-153.) < 5. ) // station 2
811  newMomPar = 5.602; // value for stations 1-2
812  if ( TMath::Abs(gemHit2Z-188.) < 5. ) // station 3
813  newMomPar = 10.188; // value for stations 1-3
814  }
815 
816  trackMomentum = newMomPar/(TMath::RadToDeg()*TMath::Abs(pangle2-pangle));
817 
818  trackTheta = TMath::ATan(mom_trans/trackMomentum);
819 
820  // if ( fVerbose > 3 || printInfo )
821  // cout << "calculated track momentum is " << trackMomentum << endl;
822  Double_t trackPhiAngle = pangle+(pangle-pangle2)*zDiffRatio;
823  if ( trackPhiAngle < 0. ) trackPhiAngle += TMath::Pi()*2.;
824  if ( trackPhiAngle > TMath::Pi()*2. ) trackPhiAngle -= TMath::Pi()*2.;
825 
826  // if ( fVerbose > 3 || printInfo )
827  // cout << "calculated phi is " << trackPhiAngle*TMath::RadToDeg() << endl;
828 
829  // if ( fVerbose ) {
830  // cout << "OK, THE MOMENTUM_" << gemHit1->GetStationNr() << gemHit2->GetStationNr()
831  // << " mom = " << trackMomentum
832  // << " the = " << trackTheta*TMath::RadToDeg()
833  // << " phi = " << trackPhiAngle*TMath::RadToDeg()
834  // << " >>> mom_trans " << mom_trans << endl;
835  // }
836 
837  TrackSegmentTB tempTS;
838  tempTS.stationIndex[0] = stat1Id;
839  tempTS.stationIndex[1] = stat2Id;
840  tempTS.sensorNumber[0] = gemHit1->GetSensorNr();
841  tempTS.sensorNumber[1] = gemHit2->GetSensorNr();
842  tempTS.hitIndex[0] = iHit1;
843  tempTS.hitIndex[1] = iHit2;
844  tempTS.trackMom = trackMomentum;
845  tempTS.trackPhi = trackPhiAngle*TMath::RadToDeg();
846  tempTS.trackTheta = trackTheta *TMath::RadToDeg();
847  tempTS.recoTrackIndex = -1;
848 
849  // cout << "SEGMENT CREATED BY HITS @ " << gemHit1->GetTimeStamp() << " (" << gemHit1->GetZ() << ") AND " << gemHit2->GetTimeStamp() << " (" << gemHit2->GetZ() << ")" << endl;
850 
851  if ( fVerbose > 2 || printInfo ) {
852  cout << " FOUND SEGMENT (stat. " << stat1Id << " & " << stat2Id << "), hits "
853  << iHit1 << ", " << iHit2
854  << " >>> " << trackMomentum << " GeV, " << trackPhiAngle*TMath::RadToDeg() << " deg, " << trackTheta << " deg." << endl;
855  }
856 
857  fTrackSegments.push_back(tempTS);
858 
859  }
860  }
861 
862  if ( fVerbose > 2 || printInfo )
863  cout << "===>>> " << fTrackSegments.size() << " track segments" << endl;
864  return fTrackSegments.size();
865 }
866 // ------------------------------------------------------------
867 
868 // --- Private method to print track segments -----------------
869 void PndGemTrackFinderOnHitsTB::PrintTrackSegments(TClonesArray* hitArray) {
870  PndGemHit* gemHit;
871 
872  //const Int_t kNofGemStations = fDigiPar->GetNStations(); //[R.K. 01/2017] unused variable?
873 
874  for ( size_t itrc = 0 ; itrc < fTrackSegments.size() ; itrc++ ) {
875  TrackSegmentTB tempTS = fTrackSegments[itrc];
876  cout << tempTS.stationIndex[0] << " " << tempTS.stationIndex[1] << " >> segment " << itrc << ": " << flush;
877  for ( Int_t ihit = 0 ; ihit < 2 ; ihit++ ) {
878  gemHit = (PndGemHit*) hitArray->At(tempTS.hitIndex[ihit]);
879  cout << " <" << gemHit->GetX() << "/" << gemHit->GetY() << "> " << flush;
880  cout << setw(3) << tempTS.hitIndex[ihit] << flush;
881 
882  }
883  cout << setw(11) << tempTS.trackMom << " " << setw(11) << tempTS.trackPhi << " " << setw(11) << tempTS.trackTheta << endl;
884  }
885 }
886 // ------------------------------------------------------------
887 
888 // --- Private method to print track segments -----------------
889 void PndGemTrackFinderOnHitsTB::PrintMCTrackSegments(TClonesArray* hitArray) {
890  PndGemHit* gemHit;
891  PndGemMCPoint* mcPoint;
892 
893  const Int_t kNofMCTracks = fMCTrackArray->GetEntriesFast();
894 
895  vector<Int_t> nofFiredStations(kNofMCTracks,0);
896  vector<TVector3> mcTrackMomentum(kNofMCTracks);
897 
898  const Int_t kNofGemPoints = fMCPointArray->GetEntriesFast();
899 
900  const Int_t kNofGemStations = fDigiPar->GetNStations();
901  Int_t nofPointsPerStation[kNofMCTracks][kNofGemStations];
902  for ( Int_t imct = 0 ; imct < kNofMCTracks ; imct++ )
903  for ( Int_t is = 0 ; is < kNofGemStations ; is++ )
904  nofPointsPerStation[imct][is] = 0;
905 
906  for ( Int_t imcp = 0 ; imcp < kNofGemPoints ; imcp++ ) {
907  mcPoint = (PndGemMCPoint*)fMCPointArray->At(imcp);
908 
909  Int_t stationNr = fDigiPar->GetStationNr(mcPoint->GetSensorId());
910 
911  nofPointsPerStation[mcPoint->GetTrackID()][stationNr-1] += 1;
912  }
913 
914  Int_t nofCGM = 0;
915  for ( Int_t imct = 0 ; imct < kNofMCTracks ; imct++ ) {
916  // fMCTrackNofCrossedGemStations[imct] = nofCGM;
917  nofCGM = 0;
918  for ( Int_t is = 0 ; is < kNofGemStations ; is++ )
919  if ( nofPointsPerStation[imct][is] > 0 )
920  nofCGM++;
921  nofFiredStations[imct] = nofCGM;
922  PndMCTrack* mcTrack = (PndMCTrack*) fMCTrackArray->At(imct);
923  mcTrackMomentum[imct] = mcTrack->GetMomentum();
924  }
925 
926  vector<Int_t> nofTrSegments(kNofMCTracks,0);
927  vector<Int_t> nofTrMCId(kNofMCTracks,0);
928  vector<Int_t> segmentMCId(fTrackSegments.size(),-1);
929 
930  for ( size_t itrc = 0 ; itrc < fTrackSegments.size() ; itrc++ ) {
931  for ( Int_t itr = 0 ; itr < kNofMCTracks ; itr++ ) nofTrMCId[itr] = 0;
932  TrackSegmentTB tempTS = fTrackSegments[itrc];
933  cout << tempTS.stationIndex[0] << " " << tempTS.stationIndex[1] << " >> segment " << itrc << ": " << flush;
934  for ( Int_t ihit = 0 ; ihit < 2 ; ihit++ ) {
935  gemHit = (PndGemHit*) hitArray->At(tempTS.hitIndex[ihit]);
936  cout << " <" << gemHit->GetX() << "/" << gemHit->GetY() << "> " << flush;
937  cout << setw(3) << tempTS.hitIndex[ihit] << flush;
938  mcPoint = (PndGemMCPoint*) fMCPointArray->At(gemHit->GetRefIndex());
939  cout << setw(2) << gemHit->GetRefIndex() << "/" << mcPoint->GetTrackID() << ") " << flush;
940 
941  if ( ihit < 2 && nofTrMCId[mcPoint->GetTrackID()] < 1 )
942  nofTrMCId[mcPoint->GetTrackID()] += 1;
943  if ( ihit > 1 && nofTrMCId[mcPoint->GetTrackID()] < 2 )
944  nofTrMCId[mcPoint->GetTrackID()] += 2;
945  }
946  cout << setw(11) << tempTS.trackMom << " " << setw(11) << tempTS.trackPhi << " " << setw(11) << tempTS.trackTheta << endl;
947  Int_t bestMCId = -1;
948  for ( Int_t itr = 0 ; itr < kNofMCTracks ; itr++ ) {
949  if ( nofTrMCId[itr] != 3 ) continue;
950  if ( bestMCId > -1 ) { bestMCId = -1; break; }
951  bestMCId = itr;
952  }
953  cout << " " << flush;
954  if ( bestMCId == -1 ) cout << " -- NO MATCHING MC -----------------" << endl;
955  else {
956  cout << "MC " << setw(3) << bestMCId << " TRUTH: "
957  << setw(11) << mcTrackMomentum[bestMCId].Mag() << " "
958  << setw(11) << TMath::RadToDeg()*mcTrackMomentum[bestMCId].Phi()+(mcTrackMomentum[bestMCId].Phi()>=0.?0:360.) << " "
959  << setw(11) << TMath::RadToDeg()*mcTrackMomentum[bestMCId].Theta() << endl;
960  nofTrSegments[bestMCId] += 1;
961  }
962  segmentMCId[itrc] = bestMCId;
963  }
964 
965  Int_t expNofS = 0;
966  Int_t fndNofS = 0;
967  for ( Int_t imct = 0 ; imct < kNofMCTracks ; imct++ ) {
968  if ( nofTrSegments[imct] == 0 || nofFiredStations[imct] == 0 ) continue;
969  cout << " track " << imct << " fired " << nofFiredStations[imct] << " stations, and " << nofTrSegments[imct] << " segments were created:" << endl;
970 
971  for ( size_t itrc = 0 ; itrc < fTrackSegments.size() ; itrc++ ) {
972  if ( segmentMCId[itrc] != imct ) continue;
973  TrackSegmentTB tempTS = fTrackSegments[itrc];
974  cout << " " << itrc << " " << setw(11) << tempTS.trackMom << " " << setw(11) << tempTS.trackPhi << " " << setw(11) << tempTS.trackTheta << endl;
975  }
976 
977  if ( mcTrackMomentum[imct].Mag() < 0.5 ) continue;
978  Int_t expSegms = 0;
979  for ( Int_t is = 0 ; is < nofFiredStations[imct] ; is++ ) for ( Int_t is2 = is+1 ; is2 < nofFiredStations[imct] ; is2++ ) expSegms++;
980  expNofS += expSegms;
981  fndNofS += nofTrSegments[imct];
982  }
983  fNofExpectedTrackSegments += expNofS;
984  fNofFoundTrackSegments += fndNofS;
985 
986  cout << "********* FOUND " << fndNofS << " OUT OF " << expNofS << " SEGMENTS IN THIS EVENT" << endl;
987  cout << "********* FOUND " << fNofFoundTrackSegments << " OUT OF " << fNofExpectedTrackSegments << " SEGMENTS IN ALL EVENTS" << endl;
988 }
989 // ------------------------------------------------------------
990 
991 // --- Private method to print track candidates ---------------
992 void PndGemTrackFinderOnHitsTB::PrintTracks(TClonesArray* hitArray, Int_t nofRecoTracks) {
993 
994  PndGemHit* gemHit;
995 
996  const Int_t kNofStatDbl = 2*fDigiPar->GetNStations();
997 
998  for ( Int_t itr = 0 ; itr < nofRecoTracks ; itr++ ) {
999  Double_t meanMom = 0.;
1000  Double_t meanPhi = 0.;
1001  Double_t meanThe = 0.;
1002 // meanMom = 0.;
1003 // meanPhi = 0.;
1004 // meanThe = 0.;
1005 
1006  vector<Int_t> hitIndices(kNofStatDbl,-1);
1007  cout << "===================== TRACK " << itr << " ======================" << endl;
1008  Int_t nofTS = 0;
1009  for ( size_t its = 0 ; its < fTrackSegments.size() ; its++ ) {
1010  TrackSegmentTB iterTS = fTrackSegments[its];
1011  if ( iterTS.recoTrackIndex != itr ) continue;
1012  // for ( Int_t ih = 0 ; ih < 4 ; ih++ ) cout << its << " > " << ih << " > " << iterTS.hitIndex[ih] << endl;
1013  hitIndices[2*iterTS.stationIndex[0]-1+iterTS.sensorNumber[0]] = iterTS.hitIndex[0];
1014  hitIndices[2*iterTS.stationIndex[1]-1+iterTS.sensorNumber[1]] = iterTS.hitIndex[1];
1015  cout << "segm " << its << " " << iterTS.trackMom << " " << iterTS.trackPhi << " " << iterTS.trackTheta << " "
1016  << iterTS.stationIndex[0] << "." << iterTS.sensorNumber[0] << ": " << iterTS.hitIndex[0] << " / "
1017  << iterTS.stationIndex[1] << "." << iterTS.sensorNumber[1] << ": " << iterTS.hitIndex[1] << endl;
1018  meanMom += iterTS.trackMom;
1019  meanPhi += iterTS.trackPhi;
1020  meanThe += iterTS.trackTheta;
1021  nofTS++;
1022  }
1023  meanMom = meanMom/nofTS;
1024  meanPhi = meanPhi/nofTS;
1025  meanThe = meanThe/nofTS;
1026 
1027  if ( nofTS == 0 ) { cout << "THIS TRACK WAS REMOVED " << endl; continue; }
1028 
1029  for ( Int_t ihit = 0 ; ihit < kNofStatDbl ; ihit++ ) {
1030  cout << ihit << "/" << hitIndices[ihit] << "/" << flush;
1031  if ( hitIndices[ihit] == -1 ) { cout << "- - " << flush; continue; }
1032  gemHit = (PndGemHit*) hitArray->At(hitIndices[ihit]);
1033  if ( gemHit->GetRefIndex() == -1 ) { cout << "- " << flush; continue;}
1034  cout << setw(2) << gemHit->GetRefIndex() << "_ " << flush;
1035  }
1036  cout << endl;
1037  }
1038 }
1039 // ------------------------------------------------------------
1040 
1041 // --- Private method to print track candidates ---------------
1042 void PndGemTrackFinderOnHitsTB::PrintMCTracks(TClonesArray* hitArray, Int_t nofRecoTracks) {
1043 
1044  PndGemHit* gemHit;
1045  FairMCPoint* mcPoint;
1046 
1047  const Int_t kNofStatDbl = 2*fDigiPar->GetNStations();
1048 
1049  for ( Int_t itr = 0 ; itr < nofRecoTracks ; itr++ ) {
1050  vector<Int_t> nofTrMCId(500,0);
1051 
1052  Double_t meanMom = 0.;
1053  Double_t meanPhi = 0.;
1054  Double_t meanThe = 0.;
1055 // meanMom = 0.;
1056 // meanPhi = 0.;
1057 // meanThe = 0.;
1058 
1059  vector<Int_t> hitIndices(kNofStatDbl,-1);
1060  cout << "===================== TRACK " << itr << " ======================" << endl;
1061  Int_t nofTS = 0;
1062  for ( size_t its = 0 ; its < fTrackSegments.size() ; its++ ) {
1063  TrackSegmentTB iterTS = fTrackSegments[its];
1064  if ( iterTS.recoTrackIndex != itr ) continue;
1065  // for ( Int_t ih = 0 ; ih < 4 ; ih++ ) cout << its << " > " << ih << " > " << iterTS.hitIndex[ih] << endl;
1066  hitIndices[2*iterTS.stationIndex[0]-1+iterTS.sensorNumber[0]] = iterTS.hitIndex[0];
1067  hitIndices[2*iterTS.stationIndex[1]-1+iterTS.sensorNumber[1]] = iterTS.hitIndex[1];
1068  cout << "segm " << its << " " << iterTS.trackMom << " " << iterTS.trackPhi << " " << iterTS.trackTheta << " "
1069  << iterTS.stationIndex[0] << "." << iterTS.sensorNumber[0] << ": " << iterTS.hitIndex[0] << " / "
1070  << iterTS.stationIndex[1] << "." << iterTS.sensorNumber[1] << ": " << iterTS.hitIndex[1] << endl;
1071  meanMom += iterTS.trackMom;
1072  meanPhi += iterTS.trackPhi;
1073  meanThe += iterTS.trackTheta;
1074  nofTS++;
1075  }
1076  meanMom = meanMom/nofTS;
1077  meanPhi = meanPhi/nofTS;
1078  meanThe = meanThe/nofTS;
1079 
1080  if ( nofTS == 0 ) { cout << "THIS TRACK WAS REMOVED " << endl; continue; }
1081 
1082  for ( Int_t ihit = 0 ; ihit < kNofStatDbl ; ihit++ ) {
1083  cout << ihit << "/" << hitIndices[ihit] << "/" << flush;
1084  if ( hitIndices[ihit] == -1 ) { cout << "- - " << flush; continue; }
1085  gemHit = (PndGemHit*) hitArray->At(hitIndices[ihit]);
1086  if ( gemHit->GetRefIndex() == -1 ) { cout << "- - " << flush; continue;}
1087  mcPoint = (FairMCPoint*) fMCPointArray->At(gemHit->GetRefIndex());
1088  cout << setw(2) << gemHit->GetRefIndex() << " _" << mcPoint->GetTrackID() << "_ " << flush;
1089  nofTrMCId[mcPoint->GetTrackID()] += 1;
1090  }
1091  cout << endl;
1092  Int_t bestMCId = -1;
1093  Int_t largestNofTracks = 0;
1094  for ( Int_t itm = 0 ; itm < 500 ; itm++ ) {
1095  if ( nofTrMCId[itm] == largestNofTracks ) { bestMCId = -1; }
1096  if ( nofTrMCId[itm] > largestNofTracks ) { bestMCId = itm; largestNofTracks = nofTrMCId[itm]; }
1097  }
1098  cout << " " << flush;
1099  cout << setw(11) << meanMom << " " << setw(11) << meanPhi << " " << setw(11) << meanThe << endl;
1100  cout << " " << flush;
1101  if ( bestMCId == -1 ) cout << " -- NO MATCHING MC -----------------" << endl;
1102  else {
1103  PndMCTrack* mcTrack = (PndMCTrack*) fMCTrackArray->At(bestMCId);
1104  TVector3 mcmom = mcTrack->GetMomentum();
1105  cout << "MC " << bestMCId << " TRUTH: "
1106  << setw(11) << mcmom.Mag() << " "
1107  << setw(11) << TMath::RadToDeg()*mcmom.Phi()+(mcmom.Phi()>=0.?0:360.) << " "
1108  << setw(11) << TMath::RadToDeg()*mcmom.Theta() << endl;
1109  }
1110  }
1111 }
1112 // ------------------------------------------------------------
1113 
1114 // --- Private method to print info at the end ---------------
1115 void PndGemTrackFinderOnHitsTB::Finish() {
1116  cout << "---------------------------------------------------------------------" << endl;
1117  cout << " >>> TF >>> prep time = " << fPrepTime << "s (get data from input)" << endl;
1118  cout << " >>> TF >>> segm time = " << fSegmTime-fPrepTime << "s (create track segments, " << fSegmTime << ")" << endl;
1119  cout << " >>> TF >>> match time = " << fMatchTime-fSegmTime << "s (match track segments, " << fMatchTime << ")" << endl;
1120  cout << " >>> TF >>> remove time = " << fRemoveTime-fMatchTime << "s (remove clone tracks, " << fRemoveTime << ")" << endl;
1121  cout << " >>> TF >>> write time = " << fWriteTime-fRemoveTime << "s (wrie tracks, " << fWriteTime << ")" << endl;
1122  cout << " >>> TF >>> all time = " << fAllTime-fWriteTime << "s (all time spent in DoFind, " << fAllTime << ")" << endl;
1123  cout << "---------------------------------------------------------------------" << endl;
1124 }
1125 // ------------------------------------------------------------
1126 
1127 ClassImp(PndGemTrackFinderOnHitsTB)
PndGemTrackFinderOnHitsTB::PrintTracks
void PrintTracks(TClonesArray *hitArray, Int_t nofRecoTracks)
Definition: PndGemTrackFinderOnHitsTB.cxx:992
pos
TVector3 pos
Definition: GammaSpectraAnalysis.C:43
PndGemTrackFinderOnHitsTB::fTrackSegments
std::vector< TrackSegmentTB > fTrackSegments
Definition: PndGemTrackFinderOnHitsTB.h:124
fVerbose
int fVerbose
Definition: poormantracks.C:24
PndGemDigiPar::GetTrackFinderOnHits_ParTheta1
Double_t GetTrackFinderOnHits_ParTheta1()
Definition: PndGemDigiPar.h:68
trackArray
TClonesArray * trackArray
Definition: NHitsPerEvent.C:14
run
Int_t run
Definition: autocutx.C:47
PndGemDigiPar::GetTrackFinderOnHits_ParThetaB
Double_t GetTrackFinderOnHits_ParThetaB()
Definition: PndGemDigiPar.h:65
PndGemDigiPar::GetTrackFinderOnHits_ParTheta3
Double_t GetTrackFinderOnHits_ParTheta3()
Definition: PndGemDigiPar.h:70
map
PndTransMap * map
Definition: sim_emc_apd.C:99
CAMath::Sqrt
static T Sqrt(const T &x)
Definition: PndCAMath.h:37
PndGemTrackFinder
Abstract base class for concrete Gem track finding algorithm.
Definition: PndGemTrackFinder.h:25
CAMath::Sin
static T Sin(const T &x)
Definition: PndCAMath.h:42
PndMCTrack::GetMomentum
TVector3 GetMomentum() const
Definition: PndMCTrack.h:78
CAMath::Tan
float Tan(float x)
Definition: PndCAMath.h:165
PndGemDigiPar::GetTrackFinderOnHits_ParMat1
Double_t GetTrackFinderOnHits_ParMat1(Int_t in)
Definition: PndGemDigiPar.h:76
PndGemDigiPar
Digitization Parameter Class for GEM part.
Definition: PndGemDigiPar.h:31
mom
Double_t mom
Definition: plot_dirc.C:14
CAMath::Cos
static T Cos(const T &x)
Definition: PndCAMath.h:43
PndGemTrackFinderOnHitsTB::PrintTrackSegments
void PrintTrackSegments(TClonesArray *hitArray)
Definition: PndGemTrackFinderOnHitsTB.cxx:869
PndGemTrackFinderOnHitsTB::RemoveCloneTracks
void RemoveCloneTracks(Int_t nofRecoTracks)
Definition: PndGemTrackFinderOnHitsTB.cxx:427
PndGemHit::GetSensorNr
Int_t GetSensorNr() const
Definition: PndGemHit.h:83
PndGemDigiPar::GetTrackFinderOnHits_ParTheta0
Double_t GetTrackFinderOnHits_ParTheta0()
Definition: PndGemDigiPar.h:67
PndGemDigiPar::GetTrackFinderOnHits_ParMat0
Double_t GetTrackFinderOnHits_ParMat0(Int_t in)
Definition: PndGemDigiPar.h:75
CAMath::Abs
static T Abs(const T &x)
Definition: PndCAMath.h:39
PndGemTrackFinderOnHitsTB
OnHits track finding algorithm.
Definition: PndGemTrackFinderOnHitsTB.h:40
PndTrackCand::AddHit
void AddHit(UInt_t detId, UInt_t hitId, Double_t rho)
Definition: PndTrackCand.cxx:36
PndGemMCPoint::GetSensorId
Int_t GetSensorId() const
Definition: PndGemMCPoint.h:90
Double_t
Double_t
Definition: SimCompleteLinkDef.h:6
rtdb
FairRuntimeDb * rtdb
Definition: hit_dirc.C:66
PndGemHit::GetCharge
Double_t GetCharge() const
Definition: PndGemHit.h:69
PndGemTrackFinderOnHitsTB::FindTrackSegments
Int_t FindTrackSegments(TClonesArray *hitArray, Int_t stat1Id, Int_t stat2Id)
Definition: PndGemTrackFinderOnHitsTB.cxx:712
is
static int is
Definition: ranlxd.cxx:374
PndGemDigiPar::GetNStations
Int_t GetNStations()
Definition: PndGemDigiPar.h:45
PndGemHit::GetStationNr
Int_t GetStationNr() const
Definition: PndGemHit.h:81
PndGemDigiPar::GetStationNr
Int_t GetStationNr(Int_t sensorId)
Definition: PndGemDigiPar.h:54
PndGemTrackFinderOnHitsTB::PrintMCTracks
void PrintMCTracks(TClonesArray *hitArray, Int_t nofRecoTracks)
Definition: PndGemTrackFinderOnHitsTB.cxx:1042
PndGemHit::GetPosition
TVector3 GetPosition() const
Definition: PndGemHit.h:71
PndGemDigiPar::GetTrackFinderOnHits_ParRadPhi2
Double_t GetTrackFinderOnHits_ParRadPhi2()
Definition: PndGemDigiPar.h:73
PndGemTrackFinderOnHitsTB::CreateTracks
Int_t CreateTracks(TClonesArray *hitArray, TClonesArray *trackArray, TClonesArray *trackCandArray, Int_t nofRecoTracks)
Definition: PndGemTrackFinderOnHitsTB.cxx:281
ClassImp
ClassImp(PndAnaContFact)
PndGemTrackFinderOnHitsTB::DoFind
virtual Int_t DoFind(TClonesArray *hitArray, TClonesArray *trackArray, TClonesArray *trackCandArray)
Definition: PndGemTrackFinderOnHitsTB.cxx:158
Pi
Double_t Pi
Definition: PndGeoHypGeCrystal.cxx:10
PndGemDigiPar::GetTrackFinderOnHits_ParThetaA
Double_t GetTrackFinderOnHits_ParThetaA()
Definition: PndGemDigiPar.h:64
PndGemDigiPar::GetTrackFinderOnHits_ParTheta2
Double_t GetTrackFinderOnHits_ParTheta2()
Definition: PndGemDigiPar.h:69
PndGemDigiPar::GetTrackFinderOnHits_ParRadPhi0
Double_t GetTrackFinderOnHits_ParRadPhi0()
Definition: PndGemDigiPar.h:72
PndGemTrackFinderOnHitsTB::PrintMCTrackSegments
void PrintMCTrackSegments(TClonesArray *hitArray)
Definition: PndGemTrackFinderOnHitsTB.cxx:889
Generated on Wed Apr 3 2019 08:43:57 for FairRoot/PandaRoot by   doxygen 1.8.5