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