FairRoot/PandaRoot
PndMvdGemTrackFinderOnHits.cxx
Go to the documentation of this file.
1 // -------------------------------------------------------------------------
2 // ----- PndMvdGemTrackFinderOnHits source file -----
3 // ----- Created 02.06.2009 by R. Karabowicz -----
4 // -------------------------------------------------------------------------
5 
6 #include "PndMvdGemTrackFinderOnHits.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 "PndSdsMCPoint.h"
16 #include "PndGemMCPoint.h"
17 #include "PndGemDigiPar.h"
18 #include "FairRootManager.h"
19 #include "PndDetectorList.h"
20 #include "PndTrack.h"
21 #include "PndTrackCand.h"
22 #include "PndTrackCandHit.h"
23 // ROOT includes
24 #include "TClonesArray.h"
25 #include "TGeoManager.h"
26 #include "TMatrixFSym.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 PndMvdGemTrackFinderOnHits::PndMvdGemTrackFinderOnHits() : PndPersistencyTask("MvdGem TrackFinder On Hits", 0) {
39  fMvdPixelHitArray = NULL;
40  fMvdStripHitArray = NULL;
41  fGemHitArray = NULL;
42 
43  fTrackArray = NULL;
44  fTrackCandArray = NULL;
45 
46  fMCTrackArray = NULL;
47  fMCGemPointArray = NULL;
48  fMCMvdPointArray = NULL;
49  fNofEvents = 0;
50 }
51 
52 // ----- Default constructor -------------------------------------------
53 PndMvdGemTrackFinderOnHits::PndMvdGemTrackFinderOnHits(Int_t iVerbose)
54  : PndPersistencyTask("MvdGem TrackFinder On Hits", iVerbose) {
55  fMvdPixelHitArray = NULL;
56  fMvdStripHitArray = NULL;
57  fGemHitArray = NULL;
58 
59  fTrackArray = NULL;
60  fTrackCandArray = NULL;
61 
62  fMCTrackArray = NULL;
63  fMCGemPointArray = NULL;
64  fMCMvdPointArray = NULL;
65  fNofEvents = 0;
66 
67  SetPersistency(kTRUE);
68 }
69 
70 // ----- Destructor ----------------------------------------------------
71 PndMvdGemTrackFinderOnHits::~PndMvdGemTrackFinderOnHits() {
72  fTrackArray ->Delete();
73  fTrackCandArray->Delete();
74 }
75 
76 // ----- Init -----------------------------------------------------------
77 InitStatus PndMvdGemTrackFinderOnHits::Init() {
78 
79  fNofExpectedTrackSegments = 0;
80  fNofFoundTrackSegments = 0;
81 
82  // Get and check FairRootManager
83  FairRootManager* ioman = FairRootManager::Instance();
84  if( !ioman ) {
85  cout << "-E- "<< GetName() <<"::Init: "
86  << "RootManager not instantised!" << endl;
87  return kERROR;
88  }
89 
90  // Get the pointer to the singleton FairRunAna object
91  FairRunAna* ana = FairRunAna::Instance();
92  if(NULL == ana) {
93  cout << "-E- "<< GetName() <<"::Init :"
94  <<" no FairRunAna object!" << endl;
95  return kERROR;
96  }
97  // Get the pointer to run-time data base
98  FairRuntimeDb* rtdb = ana->GetRuntimeDb();
99  if(NULL == rtdb) {
100  cout << "-E- "<< GetName() <<"::Init :"
101  <<" no runtime database!" << endl;
102  return kERROR;
103  }
104 
105  // Get MCTrack array
106  fMCTrackArray = (TClonesArray*) ioman->GetObject("MCTrack");
107  if( !fMCTrackArray )
108  cout << "-I- "<< GetName() <<"::Init: No MCTrack array!" << endl;
109 
110  // Get PndGemPoint (MCPoint) array
111  fMCGemPointArray = (TClonesArray*) ioman->GetObject("GEMPoint");
112  if( !fMCGemPointArray )
113  cout << "-I- "<< GetName() <<"::Init: No MCGemPoint array!" << endl;
114 
115  // Get PndMvdPoint (MCPoint) array
116  fMCMvdPointArray = (TClonesArray*) ioman->GetObject("MVDPoint");
117  if( !fMCMvdPointArray )
118  cout << "-I- "<< GetName() <<"::Init: No MCMvdPoint array!" << endl;
119 
120  // Get Mvd pixel hit Array
121  fMvdPixelHitArray = (TClonesArray*) ioman->GetObject("MVDHitsPixel");
122  if ( !fMvdPixelHitArray ) {
123  cout << "-W- "<< GetName() <<"::Init: No PndMvdPixelHit array!" << endl;
124  return kERROR;
125  }
126  // Get Mvd strip hit Array
127  fMvdStripHitArray = (TClonesArray*) ioman->GetObject("MVDHitsStrip");
128  if ( !fMvdStripHitArray ) {
129  cout << "-W- "<< GetName() <<"::Init: No PndMvdStripHit array!" << endl;
130  return kERROR;
131  }
132  // Get Gem hit Array
133  fGemHitArray = (TClonesArray*) ioman->GetObject("GEMHit");
134  if ( !fGemHitArray ) {
135  cout << "-W- "<< GetName() <<"::Init: No PndGemHit array!" << endl;
136  return kERROR;
137  }
138 
139  // Create and register output PndTrack array
140  fTrackArray = new TClonesArray("PndTrack",100);
141  ioman->Register("MVDGEMTrack", "Gem Tracks", fTrackArray, GetPersistency());
142 
143  // Create and register output PndTrackCand array
144  fTrackCandArray = new TClonesArray("PndTrackCand",100);
145  ioman->Register("MVDGEMTrackCand", "Gem Track Cands", fTrackCandArray, GetPersistency());
146 
147 
148  fDigiPar = (PndGemDigiPar*)(rtdb->getContainer("PndGemDetectors"));
149  cout << "-I- " << "PndMvdGemTrackFinderOnHits" << "::Init(). There are " << fDigiPar->GetNStations() << " GEM stations." << endl;
150  cout << "-I- " << "PndMvdGemTrackFinderOnHits" << "::Init(). Initialization succesfull." << endl;
151  fParThetaA = fDigiPar->GetTrackFinderOnHits_ParThetaA();
152  fParThetaB = fDigiPar->GetTrackFinderOnHits_ParThetaB();
153 
154  fParTheta0 = fDigiPar->GetTrackFinderOnHits_ParTheta0();
155  fParTheta1 = fDigiPar->GetTrackFinderOnHits_ParTheta1();
156  fParTheta2 = fDigiPar->GetTrackFinderOnHits_ParTheta2();
157  fParTheta3 = fDigiPar->GetTrackFinderOnHits_ParTheta3();
158 
159  fParRadPhi0 = fDigiPar->GetTrackFinderOnHits_ParRadPhi0();
160  fParRadPhi2 = fDigiPar->GetTrackFinderOnHits_ParRadPhi2();
161  for ( Int_t in = 0 ; in < 3 ; in++ ) {
162  fParMat0[in] = fDigiPar->GetTrackFinderOnHits_ParMat0(in);
163  fParMat1[in] = fDigiPar->GetTrackFinderOnHits_ParMat1(in);
164  }
165 
166  std::cout << "@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@" << std::endl;
167  std::cout << "-I- "<< GetName() <<": Intialization successfull" << std::endl;
168  std::cout << "@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@" << std::endl;
169 
170  return kSUCCESS;
171 }
172 
173 // ----- Private method SetParContainers -------------------------------
174 void PndMvdGemTrackFinderOnHits::SetParContainers() {
175 
176  // Get run and runtime database
177  FairRunAna* run = FairRunAna::Instance();
178  if ( ! run ) Fatal("SetParContainers", "No analysis run");
179 
180  FairRuntimeDb* db = run->GetRuntimeDb();
181  if ( ! db ) Fatal("SetParContainers", "No runtime database");
182 
183  // Get GEM digitisation parameter container
184  fDigiPar = (PndGemDigiPar*)(db->getContainer("PndGemDetectors"));
185 
186 }
187 // -------------------------------------------------------------------------
188 
189 // ----- Private method ReInit -----------------------------------------
190 InitStatus PndMvdGemTrackFinderOnHits::ReInit() {
191 
192  return kSUCCESS;
193 }
194 // -------------------------------------------------------------------------
195 
196 // ----- Public method DoFind ------------------------------------------
197 void PndMvdGemTrackFinderOnHits::Exec(Option_t*) {
198  if ( fVerbose > 0 ) {
199  cout << "===========================================" << endl;
200  cout << "=============== EVENT " << fNofEvents << " =================" << endl;
201  }
202  fTrackArray ->Delete();
203  fTrackCandArray->Delete();
204 
205  fNofEvents++;
206 
207  fTrackSegments.clear();
208 
209  // Check pointers
210  fMCAvailable = kTRUE;
211  if( !fMCTrackArray || !fMCGemPointArray || !fMCMvdPointArray )
212  fMCAvailable = kFALSE;
213 
214  for ( Int_t itr = 0 ; itr < fMCTrackArray->GetEntriesFast() ; itr++ ) {
215  PndMCTrack* mcTrack = (PndMCTrack*) fMCTrackArray->At(itr);
216  TVector3 mcmom = mcTrack->GetMomentum();
217 // cout << "MC " << itr << " p = : "
218 // << setw(11) << mcmom.Mag() << " GeV, "
219 // << setw(11) << TMath::RadToDeg()*mcmom.Phi()+(mcmom.Phi()>=0.?0:360.) << " deg, "
220 // << setw(11) << TMath::RadToDeg()*mcmom.Theta() << " deg." << endl;
221  }
222  /*
223  for ( Int_t itp = 0 ; itp < fMCMvdPointArray->GetEntriesFast() ; itp++ ) {
224  PndSdsMCPoint* mcPoint = (PndSdsMCPoint*) fMCMvdPointArray->At(itp);
225  cout << " MVD #" << itp << " *" << mcPoint->GetTrackID() << "* "
226  << mcPoint->GetX() << " "
227  << mcPoint->GetY() << " "
228  << mcPoint->GetZ() << " " << endl;
229  }
230  for ( Int_t itp = 0 ; itp < fMCGemPointArray->GetEntriesFast() ; itp++ ) {
231  PndGemMCPoint* mcPoint = (PndGemMCPoint*) fMCGemPointArray->At(itp);
232  cout << " GEM #" << itp << " *" << mcPoint->GetTrackID() << "* "
233  << mcPoint->GetX() << " "
234  << mcPoint->GetY() << " "
235  << mcPoint->GetZ() << " " << endl;
236  }
237  */
238 
239  // Initialise control counters
240  //Int_t nNoTrack = 0; //[R.K. 01/2017] unused variable?
241  //Int_t nNoGemPoint = 0; //[R.K. 01/2017] unused variable?
242  //Int_t nNoGemHit = 0; //[R.K. 01/2017] unused variable?
243 
244  FindTrackSegments();
245  for ( Int_t istat1 = 0 ; istat1 < fDigiPar->GetNStations() ; istat1++ ) {
246  FindTrackSegments(istat1);
247  for ( Int_t istat2 = istat1+1 ; istat2 < fDigiPar->GetNStations() ; istat2++ ) {
248  FindTrackSegments(istat1,istat2);
249  }
250  }
251 
252  if ( fVerbose ) {
253  cout << "event " << fNofEvents << " >>> " << fTrackSegments.size() << " track segments. " << endl;
254  if ( fMCAvailable )
255  PrintMCTrackSegments();
256  else
257  PrintTrackSegments();
258  }
259 
260  Int_t nr = MatchTrackSegments();
261 
262  RemoveCloneTracks(nr);
263 
264  if ( fVerbose ) {
265  cout << "************************************************" << endl;
266  if ( fMCAvailable )
267  PrintMCTracks(nr);
268  else
269  PrintTracks(nr);
270  cout << "************************************************" << endl;
271  cout << "finished printing tracks" << endl;
272  }
273 
274  nr = CreateTracks(nr);
275 
276  if ( fVerbose ) {
277  cout << "------------------------------------------------" << endl;
278  cout << "!!!!!!!!!!!!!!!!!! " << nr << " tracks have been found" << endl;
279  cout << "------------------------------------------------" << endl;
280  }
281 }
282 // ------------------------------------------------------------
283 
284 // --- Private method to create tracks ------------------------
285 Int_t PndMvdGemTrackFinderOnHits::CreateTracks(Int_t nofRecoTracks) {
286  Bool_t printInfo = kFALSE;//kTRUE;
287 
288  Int_t nofCreatedTracks = 0;
289 
290  const Int_t kNofRecoTracks = nofRecoTracks;
291 
292  const Int_t maxNofHits = 100;
293  Int_t hitIndices[kNofRecoTracks][maxNofHits][2]; // detId, hitId for each track hit
294  Int_t nofHits [kNofRecoTracks];
295  Double_t meanMom[kNofRecoTracks];
296  Double_t meanPhi[kNofRecoTracks];
297  Double_t meanThe[kNofRecoTracks];
298  Int_t nofTS [kNofRecoTracks];
299 
300  PndTrackCandHit tcHit;
301  PndTrackCand* trackCand;
302  PndGemHit* gemHit;
303  PndSdsHit* mvdHit;
304 
305  Int_t itr = -1;
306  for ( itr = 0 ; itr < nofRecoTracks ; itr++ ) {
307  nofHits[itr] = 0;
308  meanMom[itr] = 0.;
309  meanPhi[itr] = 0.;
310  meanThe[itr] = 0.;
311  nofTS[itr] = 0;
312  }
313 
314  for ( size_t its = 0 ; its < fTrackSegments.size() ; its++ ) {
315  TrackSegment iterTS = fTrackSegments[its];
316  if ( iterTS.recoTrackIndex == -1 ) continue;
317  itr = iterTS.recoTrackIndex;
318 
319  Bool_t hitInTrack[2] = {kFALSE,kFALSE};
320  for ( Int_t ih1 = 0 ; ih1 < 2 ; ih1++ ) {
321  for ( Int_t ihit = 0 ; ihit < nofHits[itr] ; ihit++ )
322  if ( hitIndices[itr][ihit][0] == iterTS.detId [ih1] &&
323  hitIndices[itr][ihit][1] == iterTS.hitIndex[ih1] ) {
324  hitInTrack[ih1] = kTRUE;
325  break;
326  }
327  if ( hitInTrack[ih1] == kFALSE ) {
328  hitIndices[itr][nofHits[itr]][0] = iterTS.detId [ih1];
329  hitIndices[itr][nofHits[itr]][1] = iterTS.hitIndex[ih1];
330  nofHits[itr]++;
331  }
332  }
333 
334  if ( nofTS[itr] > 0 ) {
335  if ( meanPhi[itr]/nofTS[itr] < 5. && iterTS.trackPhi > 355. ) { iterTS.trackPhi -= 360.; }
336  if ( meanPhi[itr]/nofTS[itr] > 355. && iterTS.trackPhi < 5. ) { iterTS.trackPhi += 360.; }
337  }
338  meanMom[itr] += iterTS.trackMom;
339  meanPhi[itr] += iterTS.trackPhi;
340  meanThe[itr] += iterTS.trackTheta;
341  nofTS[itr] ++;
342  }
343  for ( itr = 0 ; itr < nofRecoTracks ; itr++ ) {
344  if ( nofTS[itr] == 0 ) continue;
345 
346  meanMom[itr] = meanMom[itr]/nofTS[itr];
347  meanPhi[itr] = meanPhi[itr]/nofTS[itr];
348  meanThe[itr] = meanThe[itr]/nofTS[itr];
349 
350  trackCand = new((*fTrackCandArray)[nofCreatedTracks]) PndTrackCand();
351 
352  if ( fVerbose > 3 || printInfo ) {
353  cout << "---------------------------------------------------------" << endl;
354  cout << " Track " << itr << " segments:" << endl;
355 
356  for ( size_t its = 0 ; its < fTrackSegments.size() ; its++ ) {
357  TrackSegment iterTS = fTrackSegments[its];
358  if ( iterTS.recoTrackIndex != itr ) continue;
359  cout << iterTS.trackMom << " GeV/c, " << iterTS.trackTheta << " deg theta, " << iterTS.trackPhi << " deg phi)" << endl;
360  }
361 
362  cout << " track has " << nofHits[itr] << " hits:" << endl;
363  }
364  Int_t nGemHits = 0;
365  Int_t nMvdStrH = 0;
366  Int_t nMvdPixH = 0;
367  for ( Int_t ihit = 0 ; ihit < nofHits[itr] ; ihit++ ) {
368  Int_t detId = kGemHit;
369  if ( hitIndices[itr][ihit][0] < 0 )
370  detId = (-hitIndices[itr][ihit][0])>>27;
371 
372  if ( fVerbose > 3 || printInfo ) {
373  cout << " #" << ihit << ". "
374  << hitIndices[itr][ihit][0] << "."
375  << hitIndices[itr][ihit][1] << " -> "
376  << detId << flush;
377  }
378 
379  if ( detId == kGemHit ) {
380  gemHit = (PndGemHit*)fGemHitArray->At(hitIndices[itr][ihit][1]);
381  trackCand->AddHit(FairRootManager::Instance()->GetBranchId("GEMHit"),
382  hitIndices[itr][ihit][1],gemHit->GetPosition().Mag());
383  nGemHits += 1;
384  if ( fVerbose > 3 || printInfo ) {
385  cout << " ( " << gemHit->GetX()
386  << " " << gemHit->GetY()
387  << " " << gemHit->GetZ()
388  << " ) " << flush;
389  }
390  }
391  else if ( detId == kMVDHitsStrip ) { // FIXME: Use FairRootManager::Instance()->GetBranchId()
392  mvdHit = (PndSdsHit*)fMvdStripHitArray->At(hitIndices[itr][ihit][1]);
393  trackCand->AddHit(FairRootManager::Instance()->GetBranchId("MVDHitsStrip"),hitIndices[itr][ihit][1],mvdHit->GetPosition().Mag());
394  nMvdStrH += 1;
395  if ( fVerbose > 3 || printInfo ) {
396  cout << " ( " << mvdHit->GetX()
397  << " " << mvdHit->GetY()
398  << " " << mvdHit->GetZ()
399  << " ) " << flush;
400  }
401  }
402  else if ( detId == kMVDHitsPixel ) {// FIXME: Use FairRootManager::Instance()->GetBranchId()
403  mvdHit = (PndSdsHit*)fMvdPixelHitArray->At(hitIndices[itr][ihit][1]);
404  trackCand->AddHit(FairRootManager::Instance()->GetBranchId("MVDHitsPixel"),hitIndices[itr][ihit][1],mvdHit->GetPosition().Mag());
405  nMvdPixH += 1;
406  if ( fVerbose > 3 || printInfo ) {
407  cout << " ( " << mvdHit->GetX()
408  << " " << mvdHit->GetY()
409  << " " << mvdHit->GetZ()
410  << " ) " << flush;
411  }
412  }
413  if ( fVerbose > 3 || printInfo ) {
414  cout << endl;
415  }
416  }
417 
418  trackCand->Sort();
419 
420  TVector3 pos(0.,0.,0.);
421  TVector3 mom;
422  mom.SetMagThetaPhi(TMath::Abs(meanMom[itr]),meanThe[itr]*TMath::DegToRad(),meanPhi[itr]*TMath::DegToRad());
423 
424  Int_t charge = -1;
425  if ( mom.Mag() < 0. ) charge = 1;
426 
427  FairTrackParP* firstPar = new FairTrackParP(pos,mom,
428  TVector3(0.5, 0.5, 0.5),
429  0.1*mom,
430  charge,
431  pos,
432  TVector3(1.,0.,0.),
433  TVector3(0.,1.,0.));
434  FairTrackParP* lastPar = new FairTrackParP(TVector3(nGemHits,nMvdStrH,nMvdPixH), //should be something else;)
435  mom,
436  TVector3(0.5, 0.5, 0.5),
437  0.1*mom,
438  charge,
439  pos,
440  TVector3(1.,0.,0.),
441  TVector3(0.,1.,0.));
442 
443  new((*fTrackArray)[nofCreatedTracks]) PndTrack(*firstPar, *lastPar, *trackCand);
444 
445  if ( fVerbose > 2 || printInfo )
446  cout << "TRACK " << nofCreatedTracks << " HAS momentum: (" << meanMom[itr] << " GeV/c, " << meanThe[itr] << " deg theta, " << meanPhi[itr] << " deg phi)" << endl;
447 
448  if ( fVerbose > 3 || printInfo ) {
449  mom.SetMagThetaPhi(TMath::Abs(meanMom[itr]),meanThe[itr]*TMath::DegToRad(),meanPhi[itr]*TMath::DegToRad());
450 
451  cout << endl;
452  cout << " momentum: (" << meanMom[itr] << " GeV/c, " << meanThe[itr] << " deg theta, " << meanPhi[itr] << " deg phi)" << endl;
453  cout << "---------------------------------------------------------" << endl;
454  }
455 
456  nofCreatedTracks++;
457  }
458  return nofCreatedTracks;
459 }
460 // ------------------------------------------------------------
461 
462 // --- Private method to remove clone track -------------------
463 void PndMvdGemTrackFinderOnHits::RemoveCloneTracks(Int_t nofRecoTracks) {
464  Bool_t printInfo = kFALSE;//TRUE;
465 
466  if ( fVerbose > 4 || printInfo )
467  cout << "Trying to remove clone tracks" << endl;
468 
469  const Int_t kNofRecoTracks = nofRecoTracks;
470 
471  const Int_t maxNofHits = 100;
472  Double_t meanMom [kNofRecoTracks]; // mean momentum of segments in track
473  Double_t meanPhi [kNofRecoTracks]; // mean phi angle of segments in track
474  Double_t meanThe [kNofRecoTracks]; // mean theta angle of segments in track
475  Int_t hitIndices[kNofRecoTracks][maxNofHits][2]; // detId, hitId for each track hit
476  Int_t nofHits [kNofRecoTracks]; // number of hits in track
477  Int_t nofTS [kNofRecoTracks]; // number of segments in track
478  Double_t trackCons [kNofRecoTracks]; // kind of chi/ndf
479  Int_t whatToDo [kNofRecoTracks]; // =itr don't do anything, =-1 remove, =n attach to track n
480 
481  for ( Int_t itr = 0 ; itr < kNofRecoTracks ; itr++ ) {
482  meanMom [itr] = 0.;
483  meanPhi [itr] = 0.;
484  meanThe [itr] = 0.;
485  nofHits [itr] = 0;
486  nofTS [itr] = 0;
487  trackCons[itr] = 0.;
488  }
489 
490  Int_t itr = -1;
491  for ( size_t its = 0 ; its < fTrackSegments.size() ; its++ ) {
492  TrackSegment iterTS = fTrackSegments[its];
493  if ( iterTS.recoTrackIndex == -1 ) continue;
494  itr = iterTS.recoTrackIndex;
495 
496  if ( nofTS[itr] != 0 ) {
497  if ( meanPhi[itr]/nofTS[itr] < 5. && iterTS.trackPhi > 355. ) { iterTS.trackPhi -= 360.; }
498  if ( meanPhi[itr]/nofTS[itr] > 355. && iterTS.trackPhi < 5. ) { iterTS.trackPhi += 360.; }
499  }
500 
501  meanMom [itr] += iterTS.trackMom;
502  meanPhi [itr] += iterTS.trackPhi;
503  meanThe [itr] += iterTS.trackTheta;
504  nofTS [itr] += 1;
505  }
506  for ( itr = 0 ; itr < kNofRecoTracks ; itr++ ) {
507  meanMom[itr] = meanMom[itr]/nofTS[itr];
508  meanPhi[itr] = meanPhi[itr]/nofTS[itr];
509  meanThe[itr] = meanThe[itr]/nofTS[itr];
510  }
511 
512  for ( size_t its = 0 ; its < fTrackSegments.size() ; its++ ) {
513  TrackSegment iterTS = fTrackSegments[its];
514  if ( iterTS.recoTrackIndex == -1 ) continue;
515  itr = iterTS.recoTrackIndex;
516 
517  trackCons[itr] += TMath::Abs((iterTS.trackMom -meanMom[itr])/meanMom[itr]);
518  trackCons[itr] += TMath::Abs( iterTS.trackPhi -meanPhi[itr]);
519  trackCons[itr] += TMath::Abs( iterTS.trackTheta-meanThe[itr]);
520  }
521 
522  for ( itr = 0 ; itr < kNofRecoTracks ; itr++ ) {
523  trackCons[itr] = trackCons[itr]/(3.*nofTS[itr]);
524 
525  if ( trackCons[itr] > 1. && nofTS[itr] < 6 ) {
526  whatToDo[itr] = -1;
527  if ( fVerbose > 4 )
528  cout << "TRACK " << itr << " with " << nofTS[itr] << " segments will be deleted" << endl;
529  continue;
530  }
531  whatToDo[itr] = itr;
532 
533  for ( Int_t itr2 = 0 ; itr2 < itr ; itr2++ ) {
534  Double_t mmm = 0.5*(TMath::Abs(meanMom[itr])+TMath::Abs(meanMom[itr2]));
535  Double_t similarity = (TMath::Abs(meanPhi[itr]-meanPhi[itr2])+
536  TMath::Abs(meanThe[itr]-meanThe[itr2])+
537  TMath::Abs(meanMom[itr]-meanMom[itr2])/mmm);
538  // cout << " --> " << meanPhi[itr ] << " " << meanThe[itr ] << " " << meanMom[itr ] << " \\ " << endl;
539  // cout << " --> " << meanPhi[itr2] << " " << meanThe[itr2] << " " << meanMom[itr2] << " : " << similarity << endl;
540  // if ( TMath::Abs(meanPhi[itr]-meanPhi[itr2]) < 0. &&
541  // TMath::Abs(meanThe[itr]-meanThe[itr2]) < 0. &&
542  // (TMath::Abs(meanMom[itr]-meanMom[itr2]) < 0.3*TMath::Abs(meanMom[itr2])||
543  // TMath::Abs(meanMom[itr]-meanMom[itr2]) < 0.3*TMath::Abs(meanMom[itr ])) ) {
544  if ( similarity < 2.35 )
545  { // 2.35 chosen basen on several events
546  if ( whatToDo[itr2] == -1 ) whatToDo[itr2] = itr2;
547  whatToDo[itr] = whatToDo[itr2];
548  break;
549  }
550  }
551 
552  if ( fVerbose > 4 || printInfo ) {
553  cout << " TRACK " << itr << " with " << nofTS[itr] << " segments:" << endl;
554  cout << " MEAN = " << meanMom[itr] << " " << meanPhi[itr] << " " << meanThe[itr] << endl;
555  cout << " CONS = " << trackCons[itr] << " ========> " << whatToDo[itr] << endl;
556  }
557  }
558 
559  for ( size_t its = 0 ; its < fTrackSegments.size() ; its++ ) {
560  TrackSegment iterTS = fTrackSegments[its];
561  if ( iterTS.recoTrackIndex == -1 ) continue;
562  iterTS .recoTrackIndex = whatToDo[iterTS.recoTrackIndex];
563  fTrackSegments[its] = iterTS;
564  if ( iterTS.recoTrackIndex == -1 ) continue;
565  itr = iterTS.recoTrackIndex;
566 
567  Bool_t hitInTrack[2] = {kFALSE,kFALSE};
568  //Bool_t staInTrack[2] = {kFALSE,kFALSE}; //[R.K. 01/2017] unused variable?
569  for ( Int_t ih1 = 0 ; ih1 < 2 ; ih1++ ) {
570  for ( Int_t ihit = 0 ; ihit < nofHits[itr] ; ihit++ )
571  if ( hitIndices[itr][ihit][0] == iterTS.detId [ih1] )
572  if ( hitIndices[itr][ihit][1] == iterTS.hitIndex[ih1] ) {
573  hitInTrack[ih1] = kTRUE;
574  }
575  if ( hitInTrack[ih1] == kFALSE ) {
576  hitIndices[itr][nofHits[itr]][0] = iterTS.detId [ih1];
577  hitIndices[itr][nofHits[itr]][1] = iterTS.hitIndex[ih1];
578  nofHits[itr]++;
579  }
580  }
581  }
582  for ( itr = 0 ; itr < kNofRecoTracks ; itr++ ) {
583  if ( fVerbose > 3 )
584  {
585  cout << "---Track " << itr << " has " << nofTS[itr] << " segments and " << nofHits[itr] << " hits:" << endl;
586  for ( Int_t ihit = 0 ; ihit < nofHits[itr] ; ihit++ )
587  cout << " #" << ihit << ". "
588  << hitIndices[itr][ihit][0] << "."
589  << hitIndices[itr][ihit][1] << " "
590  << endl;
591  }
592  }
593 }
594 // ------------------------------------------------------------
595 
596 
597 // --- Private method to match track segments -----------------
598 Int_t PndMvdGemTrackFinderOnHits::MatchTrackSegments() {
599  // const Int_t kMaxNofSegments = fDigiPar->GetNStations()-1;
600 
601  Bool_t printInfo = kFALSE;//TRUE;
602 
603  Int_t nofRecoTracks = 0;
604 
605  for ( size_t itrc = 0 ; itrc < fTrackSegments.size() ; itrc++ ) {
606  TrackSegment origTS = fTrackSegments[itrc];
607  if ( origTS.recoTrackIndex != -1 ) continue; // this track segment has already been used
608  if ( fVerbose > 4 || printInfo )
609  cout << "MATCH TRACK SEGMENTS TO SEGM " << itrc << " with params: " << origTS.trackMom << " " << origTS.trackPhi << " " << origTS.trackTheta << endl;
610  // if ( origTS.stationIndex[1] != origTS.stationIndex[0]+1 ) continue; // only consider segments of hits on neighbouring stations
611 
612  vector<Int_t> trackSegs;
613  trackSegs.push_back(itrc);
614 
615  for ( size_t itrc2 = 0 ; itrc2 < fTrackSegments.size() ; itrc2++ ) {
616  if ( itrc2 == itrc ) continue; // do not match segment with itself
617  TrackSegment matchTS = fTrackSegments[itrc2];
618  if ( origTS.detId[0] == matchTS.detId[0] && origTS.detId[1] == matchTS.detId[1] ) continue;
619  if ( matchTS.recoTrackIndex != -1 ) continue; // this track segment has already been used
620  // if ( matchTS.stationIndex[1] != matchTS.stationIndex[0]+1 ) continue; // only consider segments of hits on neighbouring stations
621 
622  // check if track matches with any of the track segments in trackSegs
623  Double_t meanMom = 0.;
624  Double_t meanPhi = 0.;
625  Double_t meanThe = 0.;
626  Bool_t matching = kFALSE;
627  for ( size_t its = 0 ; its < trackSegs.size() ; its++ ) {
628  TrackSegment iterTS = fTrackSegments[trackSegs[its]];
629 
630  if ( TMath::Abs(matchTS.trackPhi -iterTS.trackPhi) < 1.0 &&
631  TMath::Abs(matchTS.trackTheta-iterTS.trackTheta) < 0.2 &&
632  TMath::Abs(matchTS.trackMom -iterTS.trackMom) < 0.1*TMath::Abs(iterTS.trackMom) ) {
633  matching = kTRUE;
634  break;
635  }
636  }
637 
638 // cout << "mean:(" << meanMom << "," << meanPhi << "," << meanThe << ") "
639 // << "does" << (matching?"":" not") << " match with "
640 // << "(" << matchTS.trackMom << "," << matchTS.trackPhi << "," << matchTS.trackTheta << ")" << endl;
641  if ( !matching ) continue; // matchTS does not match to any track segment in trackSegs
642 
643  // cout << "segment belongs to track" << endl;
644  // check if there exists in trackSegs a segment that uses different hits on some station
645  Bool_t sameStationSegment = kFALSE;
646  for ( size_t its = 0 ; its < trackSegs.size() ; its++ ) {
647  TrackSegment iterTS = fTrackSegments[trackSegs[its]];
648  // different stations
649  if ( matchTS.detId[0] != iterTS.detId[0] || matchTS.detId[1] != iterTS.detId[1] )
650  continue;
651  sameStationSegment = kTRUE;
652  if ( matchTS.hitIndex[0] == iterTS.hitIndex[0] && matchTS.hitIndex[1] != iterTS.hitIndex[1] )
653  continue; // the same segment
654  if ( fVerbose > 4 || printInfo )
655  cout << "there are already " << trackSegs.size() << " segments in this track: " << endl;
656  meanMom = 0.;
657  meanPhi = 0.;
658  meanThe = 0.;
659  for ( size_t its2 = 0 ; its2 < trackSegs.size() ; its2++ ) {
660  if ( its == its2 ) continue; // the track segment in question should not go to mean
661  TrackSegment aveTS = fTrackSegments[trackSegs[its2]];
662  if ( fVerbose > 4 || printInfo )
663  cout << its2 << " (" << trackSegs[its2] << ") > "
664  << aveTS.trackMom << " " << aveTS.trackPhi << " " << aveTS.trackTheta
665  << " ---- size = " << trackSegs.size()-1 << endl;
666  meanMom += aveTS.trackMom/(trackSegs.size()-1);
667  meanPhi += aveTS.trackPhi/(trackSegs.size()-1);
668  meanThe += aveTS.trackTheta/(trackSegs.size()-1);
669  }
670  if ( fVerbose > 4 || printInfo ) {
671  cout << its << " (" << trackSegs[its] << ") > " << iterTS.trackMom << " " << iterTS.trackPhi << " " << iterTS.trackTheta << endl;
672  cout << " this track conflicts with track " << its << " and is: " << endl;
673  cout << "C (" << itrc2 << ") > " << matchTS.trackMom << " " << matchTS.trackPhi << " " << matchTS.trackTheta << endl;
674  cout << " should decide on mean of other track segments, which is: " << endl;
675  cout << "MEAN > " << meanMom << " " << meanPhi << " " << meanThe << endl;
676  }
677  if ( meanPhi < 5. ) {
678  if ( iterTS.trackPhi > 355. ) iterTS.trackPhi-=360.;
679  if ( matchTS.trackPhi > 355. ) matchTS.trackPhi-=360.;
680  }
681  if ( meanPhi > 355 ) {
682  if ( iterTS.trackPhi < 5. ) iterTS.trackPhi+=360.;
683  if ( matchTS.trackPhi < 5. ) matchTS.trackPhi+=360.;
684  }
685  // check if new segment is better than the old
686  if ( TMath::Abs(matchTS.trackPhi -meanPhi) < TMath::Abs(iterTS.trackPhi -meanPhi) &&
687  TMath::Abs(matchTS.trackTheta-meanThe) < TMath::Abs(iterTS.trackTheta-meanThe) &&
688  TMath::Abs(matchTS.trackMom -meanMom) < TMath::Abs(iterTS.trackMom -meanMom) ) {
689  trackSegs[its] = itrc2;
690  continue;
691  }
692  }
693 
694  if ( sameStationSegment ) continue;
695 
696  // cout << "trackSegs ( now with " << trackSegs.size() << " ) increased by track: " << itrc2 << endl;
697  trackSegs.push_back(itrc2);
698  }
699  if ( trackSegs.size() < 5 ) continue;
700 
701  if ( fVerbose > 4 || printInfo )
702  cout << "track " << nofRecoTracks << " segments: " << endl;
703  for ( size_t its = 0 ; its < trackSegs.size() ; its++ ) {
704  TrackSegment iterTS = fTrackSegments[trackSegs[its]];
705  if ( fVerbose > 4 || printInfo )
706  {
707  cout << iterTS.detId[0] << " " << iterTS.detId[1] << " > " << flush;
708  for ( Int_t ihit = 0 ; ihit < 2 ; ihit++ )
709  cout << setw(4) << iterTS.hitIndex[ihit] << " " << flush;
710  cout << setw(11) << iterTS.trackMom << " " << setw(11) << iterTS.trackPhi << " " << setw(11) << iterTS.trackTheta << endl;
711  }
712  iterTS.recoTrackIndex = nofRecoTracks;
713  fTrackSegments[trackSegs[its]] = iterTS;
714  }
715  if ( fVerbose > 4 || printInfo )
716  cout << " seems to belong to one track" << endl;
717  nofRecoTracks++;
718  }
719 
720  if ( fVerbose || printInfo )
721  cout << "********* " << nofRecoTracks << " track candidates have been found" << endl;
722 
723  return nofRecoTracks;
724 
725 }
726 // ------------------------------------------------------------
727 
728 // ---- Private function to find track segments on 2 stations --------------------
729 Int_t PndMvdGemTrackFinderOnHits::FindTrackSegments(Int_t stat1Id, Int_t stat2Id) {
730  Bool_t printInfo = kFALSE;//TRUE;//FALSE;
731 
732  PndGemStation* stat1 = (PndGemStation*)fDigiPar->GetStation(stat1Id);
733  PndGemStation* stat2 = (PndGemStation*)fDigiPar->GetStation(stat2Id);
734  Double_t zStation1 = stat1->GetZ();
735  Double_t zStation2 = stat2->GetZ();
736 
737  Double_t zDistRatio = zStation2/zStation1;
738  Double_t zDiffRatio = zStation1/(zStation2-zStation1);
739  Double_t zCuDiff = zStation2*zStation2*zStation2-zStation1*zStation1*zStation1;
740  Double_t zSqDiff = zStation2*zStation2-zStation1*zStation1;
741  Double_t zDiff = zStation2-zStation1;
742 
743  Double_t par0_mom = fParMat0[0]*zCuDiff+fParMat0[1]*zSqDiff+fParMat0[2]*zDiff;
744  Double_t par1_mom = fParMat1[0]*zCuDiff+fParMat1[1]*zSqDiff+fParMat1[2]*zDiff;
745 
746  Int_t nGemHits = fGemHitArray->GetEntriesFast();
747 
748  PndGemHit* gemHit;
749  PndGemHit* gemHit2;
750 
751  for(Int_t iHit = 0; iHit < nGemHits; iHit++){
752  // Get the pointer to Gem hit
753  gemHit = (PndGemHit*) fGemHitArray->At(iHit);
754  if ( gemHit->GetStationNr() != stat1Id+1 ) continue;
755  zStation1 = gemHit->GetZ();
756 
757  if ( fVerbose > 3 || printInfo )
758  cout << "LOOKING FOR TRACK SEGMENTS STARTING AT " << gemHit->GetX() << " " << gemHit->GetY() << " " << gemHit->GetZ() << "( " << zStation1 << ") ndigihits = " << gemHit->GetNDigiHits() << endl;
759  // if ( gemHit->GetZ() > zStation1 ) continue;
760  if ( fVerbose > 3 || printInfo )
761  cout << "possible segment " << fTrackSegments.size() << " starts here " << gemHit->GetX() << " " << gemHit->GetY() << " " << gemHit->GetZ() << endl;
762  Double_t radius = TMath::Sqrt(gemHit->GetX()*gemHit->GetX()+gemHit->GetY()*gemHit->GetY());
763  Double_t pangle = TMath::ACos(gemHit->GetX()/radius);
764  if ( gemHit->GetY() < 0 )
765  pangle = 2.*TMath::Pi() - pangle;
766  //Double_t theta = TMath::RadToDeg()*TMath::ATan(radius/zStation1);
767  Double_t theta = fParThetaA * radius / zStation1 + fParThetaB;
768  if ( fVerbose > 3 || printInfo )
769  cout << " -> with theta of " << theta << " (radius = " << radius << " and phi angle = " << pangle*TMath::RadToDeg() << ")" << endl;
770  for(Int_t iHit2 = 0; iHit2 < nGemHits; iHit2++){
771  gemHit2 = (PndGemHit*) fGemHitArray->At(iHit2);
772  if ( gemHit2->GetStationNr() != stat2Id+1 ) continue;
773  // cout << "mixing station " << stat1Id+1 << " and " << stat2Id+1 << " hits, z = " << gemHit->GetZ() << " vs " << gemHit2->GetZ() << endl;
774  // if ( TMath::Abs(gemHit2->GetZ()-(zStation2-.6)) > 0.3 ) continue;
775  // if ( gemHit2->GetNDigiHits() < 0 ) continue;
776 
777  zStation2 = gemHit2->GetZ();
778 
779  zDistRatio = zStation2/zStation1;
780  zDiffRatio = zStation1/(zStation2-zStation1);
781  zCuDiff = zStation2*zStation2*zStation2-zStation1*zStation1*zStation1;
782  zSqDiff = zStation2*zStation2-zStation1*zStation1;
783  zDiff = zStation2-zStation1;
784 
785  par0_mom = fParMat0[0]*zCuDiff+fParMat0[1]*zSqDiff+fParMat0[2]*zDiff;
786  par1_mom = fParMat1[0]*zCuDiff+fParMat1[1]*zSqDiff+fParMat1[2]*zDiff;
787 
788  if ( fVerbose > 3 || printInfo )
789  cout << "trying to match it with " << gemHit2->GetX() << " " << gemHit2->GetY() << " " << gemHit2->GetZ() << endl;
790  Double_t radius2 = TMath::Sqrt(gemHit2->GetX()*gemHit2->GetX()+gemHit2->GetY()*gemHit2->GetY());
791  Double_t pangle2 = TMath::ACos(gemHit2->GetX()/radius2);
792  if ( gemHit2->GetY() < 0 )
793  pangle2 = 2.*TMath::Pi() - pangle2;
794 
795  if ( pangle < TMath::Pi()/2. && pangle2 > TMath::Pi()*3./2. ) pangle2 -= TMath::Pi()*2.;
796  if ( pangle > TMath::Pi()*3./2. && pangle2 < TMath::Pi()/2. ) pangle2 += TMath::Pi()*2.;
797 
798  if ( TMath::Abs(pangle-pangle2)*TMath::RadToDeg() > 40 ) continue;
799 
800  if ( fVerbose > 3 || printInfo )
801  cout << " (radius = " << radius2 << " and phi angle = " << pangle2*TMath::RadToDeg() << ")" << endl;
802 
803  Double_t expectedRad2 = (fParRadPhi0 + fParRadPhi2*TMath::RadToDeg()*TMath::RadToDeg()*(pangle-pangle2)*(pangle-pangle2))*radius*zDistRatio;
804  Double_t expRadUncert = 0.05*radius*zDistRatio;
805 
806  if ( fVerbose > 3 || printInfo ) {
807  cout << " expRad = "
808  << (fParRadPhi0 + fParRadPhi2*TMath::RadToDeg()*TMath::RadToDeg()*(pangle-pangle2)*(pangle-pangle2))
809  << " * " << radius << " * " << zDistRatio << endl;
810  cout << " -> while expected radius was " << expectedRad2 << endl;
811  }
812 
813  if ( radius2>expectedRad2+expRadUncert || radius2<expectedRad2-expRadUncert ) continue;
814 
815  if ( fVerbose > 4 || printInfo )
816  cout << "STRONG CORRELATION FOR THIS HIT!!!" << endl;
817  // calculate phi and momentum basing on the pangle-pangle2;
818  Double_t trackMomentum = 666.;
819  if ( (pangle-pangle2) != 0 )
820  trackMomentum = (par0_mom+par1_mom*radius)/((pangle-pangle2)*TMath::RadToDeg());;
821  if ( fVerbose > 3 || printInfo )
822  cout << "calculated track momentum is " << trackMomentum << endl;
823  Double_t trackPhiAngle = pangle+(pangle-pangle2)*zDiffRatio;
824  if ( trackPhiAngle < 0. ) trackPhiAngle += TMath::Pi()*2.;
825  if ( trackPhiAngle > TMath::Pi()*2. ) trackPhiAngle -= TMath::Pi()*2.;
826  if ( fVerbose > 3 || printInfo )
827  cout << "calculated phi is " << trackPhiAngle*TMath::RadToDeg() << endl;
828 
829  theta = ( fParTheta0 + 1. / ( TMath::Abs(trackMomentum) + fParTheta1 * zStation1 + fParTheta2 ) ) / zStation1 * radius + fParTheta3;
830 
831  TrackSegment tempTS;
832  tempTS.detId[0] = gemHit ->GetDetectorID();
833  tempTS.detId[1] = gemHit2->GetDetectorID();
834  tempTS.hitIndex[0] = iHit;
835  tempTS.hitIndex[1] = iHit2;
836  tempTS.trackMom = trackMomentum;
837  tempTS.trackPhi = trackPhiAngle*TMath::RadToDeg();
838  tempTS.trackTheta = theta;
839  tempTS.recoTrackIndex = -1;
840 
841  if ( fVerbose > 2 || printInfo )
842  {
843  cout << " found segment (stat. " << stat1Id << " & " << stat2Id << "), hits "
844  // << iHit << ", " << gemHit->GetNDigiHits() << ", " << iHit2 << ", " << gemHit2->GetNDigiHits()
845  << iHit << "(" << tempTS.detId[0] << "), " << iHit2 << "(" << tempTS.detId[1] << "), "
846  << " >>> " << trackMomentum << " GeV, " << trackPhiAngle*TMath::RadToDeg() << " deg, " << theta << " deg." << endl;
847  }
848 
849  fTrackSegments.push_back(tempTS);
850 
851  }
852  }
853 
854  if ( fVerbose > 2 || printInfo )
855  cout << "===>>> " << fTrackSegments.size() << " track segments" << endl;
856  return fTrackSegments.size();
857 }
858 // ------------------------------------------------------------
859 
860 // ---- Private function to find track segments between mvd and gem --------------------
861 Int_t PndMvdGemTrackFinderOnHits::FindTrackSegments(Int_t stat2Id) {
862  Bool_t printInfo = kFALSE;
863 
864  Double_t zStation1;
865  Double_t zStation2;
866 
867  Double_t zDistRatio;
868  Double_t zDiffRatio;
869  Double_t zCuDiff;
870  Double_t zSqDiff;
871  Double_t zDiff;
872 
873  Double_t par0_mom;
874  Double_t par1_mom;
875 
876  Int_t nMvdPixelHits = fMvdPixelHitArray->GetEntriesFast();
877  Int_t nMvdStripHits = fMvdStripHitArray->GetEntriesFast();
878  Int_t nGemHits = fGemHitArray->GetEntriesFast();
879 
880  PndSdsHit* mvdHit1;
881  PndGemHit* gemHit2;
882 
883  for ( Int_t iMvdHit1 = 0 ; iMvdHit1 < nMvdPixelHits+nMvdStripHits ; iMvdHit1++ ) {
884  Int_t mvdHit1Nr = iMvdHit1;
885  if ( iMvdHit1 < nMvdPixelHits )
886  mvdHit1 = (PndSdsHit*)fMvdPixelHitArray->At(mvdHit1Nr);
887  else {
888  mvdHit1Nr = iMvdHit1-nMvdPixelHits;
889  mvdHit1 = (PndSdsHit*)fMvdStripHitArray->At(mvdHit1Nr);
890  }
891  zStation1 = mvdHit1->GetZ();
892  if ( zStation1 < 4. ) continue; // take only hits in forward part
893 
894  if ( fVerbose > 3 || printInfo )
895  cout << "LOOKING FOR TRACK SEGMENTS STARTING AT " << mvdHit1->GetX() << " " << mvdHit1->GetY() << " " << mvdHit1->GetZ() << "( " << zStation1 << ")" << endl;
896  Double_t radius = TMath::Sqrt(mvdHit1->GetX()*mvdHit1->GetX()+mvdHit1->GetY()*mvdHit1->GetY());
897  Double_t pangle = TMath::ACos(mvdHit1->GetX()/radius);
898  if ( mvdHit1->GetY() < 0 )
899  pangle = 2.*TMath::Pi() - pangle;
900  //Double_t theta = TMath::RadToDeg()*TMath::ATan(radius/zStation1);
901  Double_t theta = fParThetaA * radius / zStation1 + fParThetaB;
902  if ( fVerbose > 3 || printInfo )
903  cout << " -> with theta of " << theta << " (radius = " << radius << " and phi angle = " << pangle*TMath::RadToDeg() << ")" << endl;
904 
905  for(Int_t iHit2 = 0; iHit2 < nGemHits; iHit2++){
906  gemHit2 = (PndGemHit*) fGemHitArray->At(iHit2);
907  if ( gemHit2->GetStationNr() != stat2Id+1 ) continue;
908 
909  zStation2 = gemHit2->GetZ();
910 
911  zDistRatio = zStation2/zStation1;
912  zDiffRatio = zStation1/(zStation2-zStation1);
913  zCuDiff = zStation2*zStation2*zStation2-zStation1*zStation1*zStation1;
914  zSqDiff = zStation2*zStation2-zStation1*zStation1;
915  zDiff = zStation2-zStation1;
916 
917  par0_mom = fParMat0[0]*zCuDiff+fParMat0[1]*zSqDiff+fParMat0[2]*zDiff;
918  par1_mom = fParMat1[0]*zCuDiff+fParMat1[1]*zSqDiff+fParMat1[2]*zDiff;
919 
920  if ( fVerbose > 3 || printInfo )
921  cout << "trying to match it with " << gemHit2->GetX() << " " << gemHit2->GetY() << " " << gemHit2->GetZ() << endl;
922  Double_t radius2 = TMath::Sqrt(gemHit2->GetX()*gemHit2->GetX()+gemHit2->GetY()*gemHit2->GetY());
923  Double_t pangle2 = TMath::ACos(gemHit2->GetX()/radius2);
924  if ( gemHit2->GetY() < 0 )
925  pangle2 = 2.*TMath::Pi() - pangle2;
926 
927  if ( pangle < TMath::Pi()/2. && pangle2 > TMath::Pi()*3./2. ) pangle2 -= TMath::Pi()*2.;
928  if ( pangle > TMath::Pi()*3./2. && pangle2 < TMath::Pi()/2. ) pangle2 += TMath::Pi()*2.;
929 
930  if ( TMath::Abs(pangle-pangle2)*TMath::RadToDeg() > 40 ) continue;
931 
932  if ( fVerbose > 3 || printInfo )
933  cout << " (radius = " << radius2 << " and phi angle = " << pangle2*TMath::RadToDeg() << ")" << endl;
934 
935  // Double_t expectedRad2 = (fParRadPhi0 + fParRadPhi2*TMath::RadToDeg()*TMath::RadToDeg()*(pangle-pangle2)*(pangle-pangle2))*radius*zDistRatio;
936  Double_t expectedRad2 = radius*zDistRatio;
937  Double_t expRadUncert = 0.05*radius*zDistRatio;
938 
939  if ( fVerbose > 3 || printInfo ) {
940  cout << " expRad = "
941  << (fParRadPhi0 + fParRadPhi2*TMath::RadToDeg()*TMath::RadToDeg()*(pangle-pangle2)*(pangle-pangle2))
942  << " * " << radius << " * " << zDistRatio << endl;
943  cout << " -> while expected radius was " << expectedRad2 << endl;
944  }
945 
946  if ( radius2>expectedRad2+expRadUncert || radius2<expectedRad2-expRadUncert ) continue;
947 
948  if ( fVerbose > 4 || printInfo )
949  cout << "STRONG CORRELATION FOR THIS HIT!!!" << endl;
950  // calculate phi and momentum basing on the pangle-pangle2;
951  Double_t trackMomentum = 666.;
952  if ( (pangle-pangle2) != 0 )
953  trackMomentum = (par0_mom+par1_mom*radius)/((pangle-pangle2)*TMath::RadToDeg());;
954  if ( fVerbose > 3 || printInfo )
955  cout << "calculated track momentum is " << trackMomentum << endl;
956  Double_t trackPhiAngle = pangle+(pangle-pangle2)*zDiffRatio;
957  if ( trackPhiAngle < 0. ) trackPhiAngle += TMath::Pi()*2.;
958  if ( trackPhiAngle > TMath::Pi()*2. ) trackPhiAngle -= TMath::Pi()*2.;
959  if ( fVerbose > 3 || printInfo )
960  cout << "calculated phi is " << trackPhiAngle*TMath::RadToDeg() << endl;
961 
962  theta = ( fParTheta0 + 1. / ( TMath::Abs(trackMomentum) + fParTheta1 * zStation1 + fParTheta2 ) ) / zStation1 * radius + fParTheta3;
963 
964  TrackSegment tempTS;
965  tempTS.detId[0] = -(mvdHit1->GetSensorID()+(mvdHit1->GetDetectorID()<<27));
966  tempTS.detId[1] = gemHit2->GetDetectorID();
967  tempTS.hitIndex[0] = mvdHit1Nr;
968  tempTS.hitIndex[1] = iHit2;
969  tempTS.trackMom = trackMomentum;
970  tempTS.trackPhi = trackPhiAngle*TMath::RadToDeg();
971  tempTS.trackTheta = theta;
972  tempTS.recoTrackIndex = -1;
973 
974  if ( fVerbose > 2 || printInfo )
975  {
976  cout << " found segment "
977  // << iHit << ", " << gemHit->GetNDigiHits() << ", " << iHit2 << ", " << gemHit2->GetNDigiHits()
978  << mvdHit1Nr << "(" << tempTS.detId[0]
979  //<< " + " << mvdHit1->GetSensorID() << " = " << mvdHit1->GetSensorID()+(tempTS.detId[0]<<27)
980  << "), " << iHit2 << "(" << tempTS.detId[1] << "), "
981  << " >>> " << trackMomentum << " GeV, " << trackPhiAngle*TMath::RadToDeg() << " deg, " << theta << " deg." << endl;
982  }
983 
984  fTrackSegments.push_back(tempTS);
985 
986  }
987  }
988 
989  if ( fVerbose > 2 || printInfo )
990  cout << "===>>> " << fTrackSegments.size() << " track segments" << endl;
991  return fTrackSegments.size();
992 }
993 // ------------------------------------------------------------
994 
995 // ---- Private function to find track segments between mvd and gem --------------------
996 Int_t PndMvdGemTrackFinderOnHits::FindTrackSegments() {
997  Bool_t printInfo = kFALSE;//kTRUE;
998 
999  Double_t zStation1;
1000  Double_t zStation2;
1001 
1002  Double_t zDistRatio;
1003  Double_t zDiffRatio;
1004  Double_t zCuDiff;
1005  Double_t zSqDiff;
1006  Double_t zDiff;
1007 
1008  Double_t par0_mom;
1009  Double_t par1_mom;
1010 
1011  Int_t nMvdPixelHits = fMvdPixelHitArray->GetEntriesFast();
1012  Int_t nMvdStripHits = fMvdStripHitArray->GetEntriesFast();
1013  //Int_t nGemHits = fGemHitArray->GetEntriesFast(); //[R.K. 01/2017] unused variable?
1014 
1015  PndSdsHit* mvdHit1;
1016  PndSdsHit* mvdHit2;
1017 
1018  for ( Int_t iMvdHit1 = 0 ; iMvdHit1 < nMvdPixelHits+nMvdStripHits ; iMvdHit1++ ) {
1019  Int_t mvdHit1Nr = iMvdHit1;
1020  if ( iMvdHit1 < nMvdPixelHits )
1021  mvdHit1 = (PndSdsHit*)fMvdPixelHitArray->At(mvdHit1Nr);
1022  else {
1023  mvdHit1Nr = iMvdHit1-nMvdPixelHits;
1024  mvdHit1 = (PndSdsHit*)fMvdStripHitArray->At(mvdHit1Nr);
1025  }
1026  zStation1 = mvdHit1->GetZ();
1027  if ( zStation1 < 4. ) continue; // take only hits in forward part
1028 
1029  if ( fVerbose > 3 || printInfo )
1030  cout << "LOOKING FOR TRACK SEGMENTS STARTING AT " << mvdHit1->GetX() << " " << mvdHit1->GetY() << " " << mvdHit1->GetZ() << "( " << zStation1 << ")" << endl;
1031  Double_t radius = TMath::Sqrt(mvdHit1->GetX()*mvdHit1->GetX()+mvdHit1->GetY()*mvdHit1->GetY());
1032  Double_t pangle = TMath::ACos(mvdHit1->GetX()/radius);
1033  if ( mvdHit1->GetY() < 0 )
1034  pangle = 2.*TMath::Pi() - pangle;
1035  //Double_t theta = TMath::RadToDeg()*TMath::ATan(radius/zStation1);
1036  Double_t theta = fParThetaA * radius / zStation1 + fParThetaB;
1037  if ( fVerbose > 3 || printInfo )
1038  cout << " -> with theta of " << theta << " (radius = " << radius << " and phi angle = " << pangle*TMath::RadToDeg() << ")" << endl;
1039 
1040  for ( Int_t iMvdHit2 = iMvdHit1+1 ; iMvdHit2 < nMvdPixelHits+nMvdStripHits ; iMvdHit2++ ) {
1041  Int_t mvdHit2Nr = iMvdHit2;
1042  if ( iMvdHit2 < nMvdPixelHits )
1043  mvdHit2 = (PndSdsHit*)fMvdPixelHitArray->At(mvdHit2Nr);
1044  else {
1045  mvdHit2Nr = iMvdHit2-nMvdPixelHits;
1046  mvdHit2 = (PndSdsHit*)fMvdStripHitArray->At(mvdHit2Nr);
1047  }
1048 
1049  zStation2 = mvdHit2->GetZ();
1050  if ( zStation2 < 4. ) continue; // take only hits in forward part
1051 
1052  if ( TMath::Abs(zStation1-zStation2) < 2. ) continue;
1053 
1054  zStation2 = mvdHit2->GetZ();
1055 
1056  zDistRatio = zStation2/zStation1;
1057  zDiffRatio = zStation1/(zStation2-zStation1);
1058  zCuDiff = zStation2*zStation2*zStation2-zStation1*zStation1*zStation1;
1059  zSqDiff = zStation2*zStation2-zStation1*zStation1;
1060  zDiff = zStation2-zStation1;
1061 
1062  par0_mom = fParMat0[0]*zCuDiff+fParMat0[1]*zSqDiff+fParMat0[2]*zDiff;
1063  par1_mom = fParMat1[0]*zCuDiff+fParMat1[1]*zSqDiff+fParMat1[2]*zDiff;
1064 
1065  if ( fVerbose > 3 || printInfo )
1066  cout << "trying to match it with " << mvdHit2->GetX() << " " << mvdHit2->GetY() << " " << mvdHit2->GetZ() << endl;
1067  Double_t radius2 = TMath::Sqrt(mvdHit2->GetX()*mvdHit2->GetX()+mvdHit2->GetY()*mvdHit2->GetY());
1068  Double_t pangle2 = TMath::ACos(mvdHit2->GetX()/radius2);
1069  if ( mvdHit2->GetY() < 0 )
1070  pangle2 = 2.*TMath::Pi() - pangle2;
1071 
1072  if ( pangle < TMath::Pi()/2. && pangle2 > TMath::Pi()*3./2. ) pangle2 -= TMath::Pi()*2.;
1073  if ( pangle > TMath::Pi()*3./2. && pangle2 < TMath::Pi()/2. ) pangle2 += TMath::Pi()*2.;
1074 
1075  if ( TMath::Abs(pangle-pangle2)*TMath::RadToDeg() > 40 ) continue;
1076 
1077  if ( fVerbose > 3 || printInfo )
1078  cout << " (radius = " << radius2 << " and phi angle = " << pangle2*TMath::RadToDeg() << ")" << endl;
1079 
1080  Double_t expectedRad2 = (fParRadPhi0 + fParRadPhi2*TMath::RadToDeg()*TMath::RadToDeg()*(pangle-pangle2)*(pangle-pangle2))*radius*zDistRatio;
1081  Double_t expRadUncert = 0.05*radius*zDistRatio;
1082 
1083  if ( fVerbose > 3 || printInfo ) {
1084  cout << " expRad = "
1085  << (fParRadPhi0 + fParRadPhi2*TMath::RadToDeg()*TMath::RadToDeg()*(pangle-pangle2)*(pangle-pangle2))
1086  << " * " << radius << " * " << zDistRatio << endl;
1087  cout << " -> while expected radius was " << expectedRad2 << endl;
1088  }
1089 
1090  if ( radius2>expectedRad2+expRadUncert || radius2<expectedRad2-expRadUncert ) continue;
1091 
1092  if ( fVerbose > 4 || printInfo )
1093  cout << "STRONG CORRELATION FOR THIS HIT!!!" << endl;
1094  // calculate phi and momentum basing on the pangle-pangle2;
1095  Double_t trackMomentum = 666.;
1096  if ( (pangle-pangle2) != 0 )
1097  trackMomentum = 1.6*(par0_mom+par1_mom*radius)/((pangle-pangle2)*TMath::RadToDeg());
1098  // this 1.6:
1099  // 12.11.2010: first version of the file, parameters are adjusted to GEM detector,
1100  // for pure MVD the found momentum is too far off
1101  if ( fVerbose > 3 || printInfo )
1102  cout << "calculated track momentum is " << trackMomentum << endl;
1103  Double_t trackPhiAngle = pangle+(pangle-pangle2)*zDiffRatio;
1104  if ( trackPhiAngle < 0. ) trackPhiAngle += TMath::Pi()*2.;
1105  if ( trackPhiAngle > TMath::Pi()*2. ) trackPhiAngle -= TMath::Pi()*2.;
1106  if ( fVerbose > 3 || printInfo )
1107  cout << "calculated phi is " << trackPhiAngle*TMath::RadToDeg() << endl;
1108 
1109  theta = ( fParTheta0 + 1. / ( TMath::Abs(trackMomentum) + fParTheta1 * zStation1 + fParTheta2 ) ) / zStation1 * radius + fParTheta3;
1110 
1111  TrackSegment tempTS;
1112  tempTS.detId[0] = -(mvdHit1->GetSensorID()+(mvdHit1->GetDetectorID()<<27));
1113  tempTS.detId[1] = -(mvdHit2->GetSensorID()+(mvdHit2->GetDetectorID()<<27));
1114  tempTS.hitIndex[0] = mvdHit1Nr;
1115  tempTS.hitIndex[1] = mvdHit2Nr;
1116  tempTS.trackMom = trackMomentum;
1117  tempTS.trackPhi = trackPhiAngle*TMath::RadToDeg();
1118  tempTS.trackTheta = theta;
1119  tempTS.recoTrackIndex = -1;
1120 
1121  if ( fVerbose > 2 || printInfo )
1122  {
1123  cout << " found segment "
1124  // << iHit << ", " << gemHit->GetNDigiHits() << ", " << mvdHit2Nr << ", " << mvdHit2->GetNDigiHits()
1125  << mvdHit1Nr << "(" << tempTS.detId[0]
1126  //<< " + " << mvdHit1->GetSensorID() << " = " << mvdHit1->GetSensorID()+(tempTS.detId[0]<<27)
1127  << "), " << mvdHit2Nr << "(" << tempTS.detId[1] << "), "
1128  << " >>> " << trackMomentum << " GeV, " << trackPhiAngle*TMath::RadToDeg() << " deg, " << theta << " deg." << endl;
1129  }
1130 
1131  fTrackSegments.push_back(tempTS);
1132 
1133  }
1134  }
1135 
1136  if ( fVerbose > 2 || printInfo )
1137  cout << "===>>> " << fTrackSegments.size() << " track segments" << endl;
1138  return fTrackSegments.size();
1139 }
1140 // ------------------------------------------------------------
1141 
1142 // --- Private method to print track segments -----------------
1143 void PndMvdGemTrackFinderOnHits::PrintTrackSegments() {
1144  /* PndGemHit* gemHit;
1145 
1146  const Int_t kNofGemStations = fDigiPar->GetNStations();
1147 
1148  for ( Int_t itrc = 0 ; itrc < fTrackSegments.size() ; itrc++ ) {
1149  TrackSegment tempTS = fTrackSegments[itrc];
1150  cout << tempTS.stationIndex[0] << " " << tempTS.stationIndex[1] << " >> segment " << itrc << ": " << flush;
1151  for ( Int_t ihit = 0 ; ihit < 4 ; ihit++ ) {
1152  if ( tempTS.hitIndex[ihit] == -1 ) continue;
1153  gemHit = (PndGemHit*) fGemHitArray->At(tempTS.hitIndex[ihit]);
1154  cout << " <" << gemHit->GetX() << "/" << gemHit->GetY() << "> " << flush;
1155  cout << setw(3) << tempTS.hitIndex[ihit] << "(" << flush;
1156  if ( gemHit->GetRefIndex() == -1 ) { cout << "--/-) " << flush; continue;}
1157  cout << setw(2) << gemHit->GetRefIndex() << ") " << flush;
1158 
1159  }
1160  cout << setw(11) << tempTS.trackMom << " " << setw(11) << tempTS.trackPhi << " " << setw(11) << tempTS.trackTheta << endl;
1161  }*/
1162 }
1163 // ------------------------------------------------------------
1164 
1165 // --- Private method to print track segments -----------------
1166 void PndMvdGemTrackFinderOnHits::PrintMCTrackSegments() {
1167  /* PndGemHit* gemHit;
1168  PndGemMCPoint* mcPoint;
1169 
1170  const Int_t kNofMCTracks = fMCTrackArray->GetEntriesFast();
1171 
1172  vector<Int_t> nofFiredStations(kNofMCTracks,0);
1173  vector<TVector3> mcTrackMomentum(kNofMCTracks,0);
1174 
1175  const Int_t kNofGemPoints = fMCGemPointArray->GetEntriesFast();
1176 
1177  const Int_t kNofGemStations = fDigiPar->GetNStations();
1178  Int_t nofPointsPerStation[kNofMCTracks][kNofGemStations];
1179  for ( Int_t imct = 0 ; imct < kNofMCTracks ; imct++ )
1180  for ( Int_t is = 0 ; is < kNofGemStations ; is++ )
1181  nofPointsPerStation[imct][is] = 0;
1182 
1183  for ( Int_t imcp = 0 ; imcp < kNofGemPoints ; imcp++ ) {
1184  mcPoint = (PndGemMCPoint*)fMCGemPointArray->At(imcp);
1185 
1186  Int_t stationNr = fDigiPar->GetStationNr(mcPoint->GetSensorId());
1187 
1188  nofPointsPerStation[mcPoint->GetTrackID()][stationNr-1] += 1;
1189  }
1190 
1191  Int_t nofCGM = 0;
1192  for ( Int_t imct = 0 ; imct < kNofMCTracks ; imct++ ) {
1193  // fMCTrackNofCrossedGemStations[imct] = nofCGM;
1194  nofCGM = 0;
1195  for ( Int_t is = 0 ; is < kNofGemStations ; is++ )
1196  if ( nofPointsPerStation[imct][is] > 0 )
1197  nofCGM++;
1198  nofFiredStations[imct] = nofCGM;
1199  PndMCTrack* mcTrack = (PndMCTrack*) fMCTrackArray->At(imct);
1200  mcTrackMomentum[imct] = mcTrack->GetMomentum();
1201  }
1202 
1203  vector<Int_t> nofTrSegments(kNofMCTracks,0);
1204  vector<Int_t> nofTrMCId(kNofMCTracks,0);
1205  vector<Int_t> segmentMCId(fTrackSegments.size(),-1);
1206 
1207  for ( Int_t itrc = 0 ; itrc < fTrackSegments.size() ; itrc++ ) {
1208  for ( Int_t itr = 0 ; itr < kNofMCTracks ; itr++ ) nofTrMCId[itr] = 0;
1209  TrackSegment tempTS = fTrackSegments[itrc];
1210  cout << tempTS.stationIndex[0] << " " << tempTS.stationIndex[1] << " >> segment " << itrc << ": " << flush;
1211  for ( Int_t ihit = 0 ; ihit < 4 ; ihit++ ) {
1212  if ( tempTS.hitIndex[ihit] == -1 ) continue;
1213  gemHit = (PndGemHit*) fGemHitArray->At(tempTS.hitIndex[ihit]);
1214  cout << " <" << gemHit->GetX() << "/" << gemHit->GetY() << "> " << flush;
1215  cout << setw(3) << tempTS.hitIndex[ihit] << "(" << flush;
1216  if ( gemHit->GetRefIndex() == -1 ) { cout << "--/-) " << flush; continue;}
1217  mcPoint = (PndGemMCPoint*) fMCGemPointArray->At(gemHit->GetRefIndex());
1218  cout << setw(2) << gemHit->GetRefIndex() << "/" << mcPoint->GetTrackID() << ") " << flush;
1219 
1220  if ( ihit < 2 && nofTrMCId[mcPoint->GetTrackID()] < 1 )
1221  nofTrMCId[mcPoint->GetTrackID()] += 1;
1222  if ( ihit > 1 && nofTrMCId[mcPoint->GetTrackID()] < 2 )
1223  nofTrMCId[mcPoint->GetTrackID()] += 2;
1224  }
1225  cout << setw(11) << tempTS.trackMom << " " << setw(11) << tempTS.trackPhi << " " << setw(11) << tempTS.trackTheta << endl;
1226  Int_t bestMCId = -1;
1227  for ( Int_t itr = 0 ; itr < kNofMCTracks ; itr++ ) {
1228  if ( nofTrMCId[itr] != 3 ) continue;
1229  if ( bestMCId > -1 ) { bestMCId = -1; break; }
1230  bestMCId = itr;
1231  }
1232  cout << " " << flush;
1233  if ( bestMCId == -1 ) cout << " -- NO MATCHING MC -----------------" << endl;
1234  else {
1235  cout << "MC " << setw(3) << bestMCId << " TRUTH: "
1236  << setw(11) << mcTrackMomentum[bestMCId].Mag() << " "
1237  << setw(11) << TMath::RadToDeg()*mcTrackMomentum[bestMCId].Phi()+(mcTrackMomentum[bestMCId].Phi()>=0.?0:360.) << " "
1238  << setw(11) << TMath::RadToDeg()*mcTrackMomentum[bestMCId].Theta() << endl;
1239  nofTrSegments[bestMCId] += 1;
1240  }
1241  segmentMCId[itrc] = bestMCId;
1242  }
1243 
1244  Int_t expNofS = 0;
1245  Int_t fndNofS = 0;
1246  for ( Int_t imct = 0 ; imct < kNofMCTracks ; imct++ ) {
1247  if ( nofTrSegments[imct] == 0 || nofFiredStations[imct] == 0 ) continue;
1248  cout << " track " << imct << " fired " << nofFiredStations[imct] << " stations, and " << nofTrSegments[imct] << " segments were created:" << endl;
1249 
1250  for ( Int_t itrc = 0 ; itrc < fTrackSegments.size() ; itrc++ ) {
1251  if ( segmentMCId[itrc] != imct ) continue;
1252  TrackSegment tempTS = fTrackSegments[itrc];
1253  cout << " " << itrc << " " << setw(11) << tempTS.trackMom << " " << setw(11) << tempTS.trackPhi << " " << setw(11) << tempTS.trackTheta << endl;
1254  }
1255 
1256  if ( mcTrackMomentum[imct].Mag() < 0.5 ) continue;
1257  Int_t expSegms = 0;
1258  for ( Int_t is = 0 ; is < nofFiredStations[imct] ; is++ ) for ( Int_t is2 = is+1 ; is2 < nofFiredStations[imct] ; is2++ ) expSegms++;
1259  expNofS += expSegms;
1260  fndNofS += nofTrSegments[imct];
1261  }
1262  fNofExpectedTrackSegments += expNofS;
1263  fNofFoundTrackSegments += fndNofS;
1264 
1265  cout << "********* FOUND " << fndNofS << " OUT OF " << expNofS << " SEGMENTS IN THIS EVENT" << endl;
1266  cout << "********* FOUND " << fNofFoundTrackSegments << " OUT OF " << fNofExpectedTrackSegments << " SEGMENTS IN ALL EVENTS" << endl;
1267  */
1268 }
1269 // ------------------------------------------------------------
1270 
1271 // --- Private method to print track candidates ---------------
1272 void PndMvdGemTrackFinderOnHits::PrintTracks(Int_t ) { // nofRecoTracks // [R.K.03/2017] unused variable(s)
1273  /*
1274  PndGemHit* gemHit;
1275 
1276  const Int_t kNofStatDbl = 2*fDigiPar->GetNStations();
1277 
1278  for ( Int_t itr = 0 ; itr < nofRecoTracks ; itr++ ) {
1279  Double_t meanMom = 0.;
1280  Double_t meanPhi = 0.;
1281  Double_t meanThe = 0.;
1282 // meanMom = 0.;
1283 // meanPhi = 0.;
1284 // meanThe = 0.;
1285 
1286  vector<Int_t> hitIndices(kNofStatDbl,-1);
1287  cout << "===================== TRACK " << itr << " ======================" << endl;
1288  Int_t nofTS = 0;
1289  for ( Int_t its = 0 ; its < fTrackSegments.size() ; its++ ) {
1290  TrackSegment iterTS = fTrackSegments[its];
1291  if ( iterTS.recoTrackIndex != itr ) continue;
1292  // for ( Int_t ih = 0 ; ih < 4 ; ih++ ) cout << its << " > " << ih << " > " << iterTS.hitIndex[ih] << endl;
1293  for ( Int_t istat = 0 ; istat < 2 ; istat++ ) {
1294  hitIndices[2*iterTS.stationIndex[istat] ] = iterTS.hitIndex[2*istat ];
1295  hitIndices[2*iterTS.stationIndex[istat]+1] = iterTS.hitIndex[2*istat+1];
1296  }
1297  cout << "segm " << its << " " << iterTS.trackMom << " " << iterTS.trackPhi << " " << iterTS.trackTheta << " "
1298  << iterTS.stationIndex[0] << ": " << iterTS.hitIndex[0] << " " << iterTS.hitIndex[1] << " / "
1299  << iterTS.stationIndex[1] << ": " << iterTS.hitIndex[2] << " " << iterTS.hitIndex[3] << endl;
1300  meanMom += iterTS.trackMom;
1301  meanPhi += iterTS.trackPhi;
1302  meanThe += iterTS.trackTheta;
1303  nofTS++;
1304  }
1305  meanMom = meanMom/nofTS;
1306  meanPhi = meanPhi/nofTS;
1307  meanThe = meanThe/nofTS;
1308 
1309  if ( nofTS == 0 ) { cout << "THIS TRACK WAS REMOVED " << endl; continue; }
1310 
1311  for ( Int_t ihit = 0 ; ihit < kNofStatDbl ; ihit++ ) {
1312  cout << ihit << "/" << hitIndices[ihit] << "/" << flush;
1313  if ( hitIndices[ihit] == -1 ) { cout << "- - " << flush; continue; }
1314  gemHit = (PndGemHit*) fGemHitArray->At(hitIndices[ihit]);
1315  if ( gemHit->GetRefIndex() == -1 ) { cout << "- " << flush; continue;}
1316  cout << setw(2) << gemHit->GetRefIndex() << "_ " << flush;
1317  }
1318  cout << endl;
1319  }
1320  */
1321 }
1322 // ------------------------------------------------------------
1323 
1324 // --- Private method to print track candidates ---------------
1325 void PndMvdGemTrackFinderOnHits::PrintMCTracks(Int_t ) { // nofRecoTracks // [R.K.03/2017] unused variable(s)
1326  /*
1327  PndGemHit* gemHit;
1328  FairMCPoint* mcPoint;
1329 
1330  const Int_t kNofStatDbl = 2*fDigiPar->GetNStations();
1331 
1332  for ( Int_t itr = 0 ; itr < nofRecoTracks ; itr++ ) {
1333  vector<Int_t> nofTrMCId(500,0);
1334 
1335  Double_t meanMom = 0.;
1336  Double_t meanPhi = 0.;
1337  Double_t meanThe = 0.;
1338 // meanMom = 0.;
1339 // meanPhi = 0.;
1340 // meanThe = 0.;
1341 
1342  vector<Int_t> hitIndices(kNofStatDbl,-1);
1343  cout << "===================== TRACK " << itr << " ======================" << endl;
1344  Int_t nofTS = 0;
1345  for ( Int_t its = 0 ; its < fTrackSegments.size() ; its++ ) {
1346  TrackSegment iterTS = fTrackSegments[its];
1347  if ( iterTS.recoTrackIndex != itr ) continue;
1348  // for ( Int_t ih = 0 ; ih < 4 ; ih++ ) cout << its << " > " << ih << " > " << iterTS.hitIndex[ih] << endl;
1349  for ( Int_t istat = 0 ; istat < 2 ; istat++ ) {
1350  hitIndices[2*iterTS.stationIndex[istat] ] = iterTS.hitIndex[2*istat ];
1351  hitIndices[2*iterTS.stationIndex[istat]+1] = iterTS.hitIndex[2*istat+1];
1352  }
1353  cout << "segm " << its << " " << iterTS.trackMom << " " << iterTS.trackPhi << " " << iterTS.trackTheta << " "
1354  << iterTS.stationIndex[0] << ": " << iterTS.hitIndex[0] << " " << iterTS.hitIndex[1] << " / "
1355  << iterTS.stationIndex[1] << ": " << iterTS.hitIndex[2] << " " << iterTS.hitIndex[3] << endl;
1356  meanMom += iterTS.trackMom;
1357  meanPhi += iterTS.trackPhi;
1358  meanThe += iterTS.trackTheta;
1359  nofTS++;
1360  }
1361  meanMom = meanMom/nofTS;
1362  meanPhi = meanPhi/nofTS;
1363  meanThe = meanThe/nofTS;
1364 
1365  if ( nofTS == 0 ) { cout << "THIS TRACK WAS REMOVED " << endl; continue; }
1366 
1367  for ( Int_t ihit = 0 ; ihit < kNofStatDbl ; ihit++ ) {
1368  cout << ihit << "/" << hitIndices[ihit] << "/" << flush;
1369  if ( hitIndices[ihit] == -1 ) { cout << "- - " << flush; continue; }
1370  gemHit = (PndGemHit*) fGemHitArray->At(hitIndices[ihit]);
1371  if ( gemHit->GetRefIndex() == -1 ) { cout << "- - " << flush; continue;}
1372  mcPoint = (FairMCPoint*) fMCGemPointArray->At(gemHit->GetRefIndex());
1373  cout << setw(2) << gemHit->GetRefIndex() << " _" << mcPoint->GetTrackID() << "_ " << flush;
1374  nofTrMCId[mcPoint->GetTrackID()] += 1;
1375  }
1376  cout << endl;
1377  Int_t bestMCId = -1;
1378  Int_t largestNofTracks = 0;
1379  for ( Int_t itm = 0 ; itm < 500 ; itm++ ) {
1380  if ( nofTrMCId[itm] == largestNofTracks ) { bestMCId = -1; }
1381  if ( nofTrMCId[itm] > largestNofTracks ) { bestMCId = itm; largestNofTracks = nofTrMCId[itm]; }
1382  }
1383  cout << " " << flush;
1384  cout << setw(11) << meanMom << " " << setw(11) << meanPhi << " " << setw(11) << meanThe << endl;
1385  cout << " " << flush;
1386  if ( bestMCId == -1 ) cout << " -- NO MATCHING MC -----------------" << endl;
1387  else {
1388  PndMCTrack* mcTrack = (PndMCTrack*) fMCTrackArray->At(bestMCId);
1389  TVector3 mcmom = mcTrack->GetMomentum();
1390  cout << "MC " << bestMCId << " TRUTH: "
1391  << setw(11) << mcmom.Mag() << " "
1392  << setw(11) << TMath::RadToDeg()*mcmom.Phi()+(mcmom.Phi()>=0.?0:360.) << " "
1393  << setw(11) << TMath::RadToDeg()*mcmom.Theta() << endl;
1394  }
1395  }
1396  */
1397 }
1398 // ------------------------------------------------------------
1399 
1400 // ----- Private method Finish --------------------------------------------
1401 void PndMvdGemTrackFinderOnHits::Finish() {
1402  fTrackArray ->Delete();
1403  fTrackCandArray->Delete();
1404 
1405  cout << "-------------------- " << fName.Data() << " : Summary -----------------------" << endl;
1406  cout << " Events: " << setw(10) << fNofEvents << endl;
1407  cout << "--------------------------------------------------------------------------------" << endl;
1408 }
1409 // ------------------------------------------------------------
1410 
1411 
1412 ClassImp(PndMvdGemTrackFinderOnHits)
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
PndMvdGemTrackFinderOnHits::fTrackSegments
std::vector< TrackSegment > fTrackSegments
Definition: PndMvdGemTrackFinderOnHits.h:115
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
PndSdsHit::GetPosition
TVector3 GetPosition() const
Definition: PndSdsHit.h:93
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
PndMCTrack::GetMomentum
TVector3 GetMomentum() const
Definition: PndMCTrack.h:78
PndMvdGemTrackFinderOnHits::CreateTracks
Int_t CreateTracks(Int_t nofRecoTracks)
Definition: PndMvdGemTrackFinderOnHits.cxx:285
PndPersistencyTask::SetPersistency
void SetPersistency(Bool_t val=kTRUE)
Definition: PndPersistencyTask.h:31
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
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
PndMvdGemTrackFinderOnHits::RemoveCloneTracks
void RemoveCloneTracks(Int_t nofRecoTracks)
Definition: PndMvdGemTrackFinderOnHits.cxx:463
CAMath::Abs
static T Abs(const T &x)
Definition: PndCAMath.h:39
PndTrackCand::AddHit
void AddHit(UInt_t detId, UInt_t hitId, Double_t rho)
Definition: PndTrackCand.cxx:36
Double_t
Double_t
Definition: SimCompleteLinkDef.h:6
PndGemHit::GetNDigiHits
Int_t GetNDigiHits() const
Definition: PndGemHit.h:70
rtdb
FairRuntimeDb * rtdb
Definition: hit_dirc.C:66
PndMvdGemTrackFinderOnHits::Exec
virtual void Exec(Option_t *opt)
Definition: PndMvdGemTrackFinderOnHits.cxx:197
PndGemDigiPar::GetNStations
Int_t GetNStations()
Definition: PndGemDigiPar.h:45
PndGemHit::GetStationNr
Int_t GetStationNr() const
Definition: PndGemHit.h:81
PndGemHit::GetPosition
TVector3 GetPosition() const
Definition: PndGemHit.h:71
PndGemDigiPar::GetTrackFinderOnHits_ParRadPhi2
Double_t GetTrackFinderOnHits_ParRadPhi2()
Definition: PndGemDigiPar.h:73
ClassImp
ClassImp(PndAnaContFact)
iVerbose
Int_t iVerbose
Definition: Lars/runMvdDigi.C:7
PndSdsHit::GetSensorID
Int_t GetSensorID() const
Definition: PndSdsHit.h:90
PndMvdGemTrackFinderOnHits
OnHits track finding algorithm.
Definition: PndMvdGemTrackFinderOnHits.h:39
Pi
Double_t Pi
Definition: PndGeoHypGeCrystal.cxx:10
PndGemDigiPar::GetTrackFinderOnHits_ParThetaA
Double_t GetTrackFinderOnHits_ParThetaA()
Definition: PndGemDigiPar.h:64
PndMvdGemTrackFinderOnHits::PrintMCTracks
void PrintMCTracks(Int_t nofRecoTracks)
Definition: PndMvdGemTrackFinderOnHits.cxx:1325
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:59 for FairRoot/PandaRoot by   doxygen 1.8.5