FairRoot/PandaRoot
PndBarrelTrackFinder.cxx
Go to the documentation of this file.
1 //* $Id: */
2 
3 // -------------------------------------------------------------------------
4 // ----- PndBarrelTrackFinder source file -----
5 // ----- Created 05/12/2010 by R. Karabowicz -----
6 // -------------------------------------------------------------------------
7 
8 // Includes from GEM
9 #include "PndBarrelTrackFinder.h"
10 
11 // Includes from base
12 #include "FairRootManager.h"
13 #include "FairRunAna.h"
14 #include "FairRuntimeDb.h"
15 #include "FairLink.h"
16 
17 // Includes from ROOT
18 #include "TClonesArray.h"
19 #include "TObjArray.h"
20 #include "TMath.h"
21 #include "TGeoManager.h"
22 #include "TGeoNode.h"
23 #include "TRandom.h"
24 
25 #include "FairTrackParam.h"
26 #include "FairHit.h"
27 #include "PndTrack.h"
28 #include "PndTrackCand.h"
29 #include "PndTrackCandHit.h"
30 
31 #include "PndDetectorList.h"
32 #include "PndGemHit.h"
33 #include "PndSdsHit.h"
34 #include "PndSttHit.h"
35 #include "PndSttMapCreator.h"
36 #include "PndSttTube.h"
37 #include "PndSttTubeMap.h"
38 
39 #include <iostream>
40 #include <iomanip>
41 #include <map>
42 
43 #include "TCanvas.h" // DRAW
44 #include "TMarker.h" // DRAW
45 #include "TEllipse.h" // DRAW
46 #include "TSystem.h" // DRAW
47 #include "TLine.h" // DRAW
48 #include "TLatex.h" // DRAW
49 #include "TH1F.h" // DRAW
50 
51 
52 using std::cout;
53 using std::cerr;
54 using std::flush;
55 using std::endl;
56 using std::pair;
57 using std::setw;
58 using std::left;
59 using std::right;
60 using std::fixed;
61 using std::setprecision;
62 using std::map;
63 
64 
65 
66 // ----- Default constructor ------------------------------------------
67 PndBarrelTrackFinder::PndBarrelTrackFinder() : PndPersistencyTask("Barrel Track Finder", 1), fGeoType(0) {
68  for ( Int_t idet = 0 ; idet < 4 ; idet++ ) {
69  fIncludeDet [idet] = kFALSE;
70  fHitArray [idet] = NULL;
71  fDetName [idet] = "";
72  fDetType [idet] = 0;
73  fHitArrayName[idet] = "";
74  }
75  fBarrelTrackArray = NULL;
76  fBarrelTrackCandArray = NULL;
77 
78  fTNofEvents = 0;
79  fTNofTracks = 0;
80  fTracksVector.clear();
81 
82  fHitDetId.clear();
83  fHitDetNo.clear();
84  fHitVector.clear();
85  fHitVectDI.clear();
86  fHitVectHN.clear();
87 
88  fMaximalDist = 0.01;
89  fReasonableDist = 0.2;
90 
91  fMaximalZ_PD = 2.;
92  fReasonableZ_PD = 5.;
93 
94  fMaximalRadDiff = 0.1;
95  fMaximalPhiDiff = 0.1;
96 
97  fTubeArray = NULL;
98  fSttParameters = NULL;
99 
100  fDrawDetails = kFALSE;
101  fDrawHistos = kFALSE;
102 
103  Reset();
104 }
105 // -------------------------------------------------------------------------
106 
107 
108 
109 // ----- Standard constructor ------------------------------------------
110 PndBarrelTrackFinder::PndBarrelTrackFinder(Int_t iVerbose)
111  : PndPersistencyTask("Barrel Track Finder", iVerbose) {
112  for ( Int_t idet = 0 ; idet < 4 ; idet++ ) {
113  fIncludeDet[idet] = kFALSE;
114  fHitArray [idet] = NULL;
115  fDetName [idet] = "";
116  fDetType [idet] = 0;
117  fHitArrayName[idet] = "";
118  }
119  fBarrelTrackArray = NULL;
120  fBarrelTrackCandArray = NULL;
121 
122  fTNofEvents = 0;
123  fTNofTracks = 0;
124  fTracksVector.clear();
125 
126  fHitDetId.clear();
127  fHitDetNo.clear();
128  fHitVector.clear();
129  fHitVectDI.clear();
130  fHitVectHN.clear();
131 
132  fMaximalDist = 0.01;
133  fReasonableDist = 0.2;
134 
135  fMaximalZ_PD = 2.;
136  fReasonableZ_PD = 5.;
137 
138  fMaximalRadDiff = 0.1;
139  fMaximalPhiDiff = 0.1;
140 
141  fTubeArray = NULL;
142  fSttParameters = NULL;
143 
144  Reset();
145 }
146 // -------------------------------------------------------------------------
147 
148 
149 
150 // ----- Constructor with name -----------------------------------------
151 PndBarrelTrackFinder::PndBarrelTrackFinder(const char* name, Int_t iVerbose)
152  : PndPersistencyTask(name, iVerbose) {
153  for ( Int_t idet = 0 ; idet < 4 ; idet++ ) {
154  fIncludeDet[idet] = kFALSE;
155  fHitArray [idet] = NULL;
156  fDetName [idet] = "";
157  fDetType [idet] = 0;
158  fHitArrayName[idet] = "";
159  }
160  fBarrelTrackArray = NULL;
161  fBarrelTrackCandArray = NULL;
162 
163  fTNofEvents = 0;
164  fTNofTracks = 0;
165  fTracksVector.clear();
166 
167  fHitDetId.clear();
168  fHitDetNo.clear();
169  fHitVector.clear();
170  fHitVectDI.clear();
171  fHitVectHN.clear();
172 
173  fMaximalDist = 0.01;
174  fReasonableDist = 0.2;
175  fMaximalZ_PD = 2.;
176  fReasonableZ_PD = 5.;
177 
178  fMaximalRadDiff = 0.1;
179  fMaximalPhiDiff = 0.1;
180 
181  fTubeArray = NULL;
182  fSttParameters = NULL;
183 
184  SetPersistency(kTRUE);
185  Reset();
186 }
187 // -------------------------------------------------------------------------
188 
189 
190 
191 // ----- Destructor ----------------------------------------------------
192 PndBarrelTrackFinder::~PndBarrelTrackFinder() {
193  fBarrelTrackArray ->Delete();
194  fBarrelTrackCandArray->Delete();
195  Reset();
196 }
197 // -------------------------------------------------------------------------
198 
199 // -------------------------------------------------------------------------
200 void PndBarrelTrackFinder::UseMvdSttGem(const Bool_t useMvd, const Bool_t useStt, const Bool_t useGem) {
201  fIncludeDet[0] = useMvd;
202  fIncludeDet[1] = useMvd;
203  fIncludeDet[2] = useStt;
204  fIncludeDet[3] = useGem;
205 };
206 // -------------------------------------------------------------------------
207 
208 // ----- Public method Exec --------------------------------------------
209 void PndBarrelTrackFinder::Exec(Option_t*) {
210  Reset();
211  fBarrelTrackArray ->Delete();
212  fBarrelTrackCandArray->Delete();
213 
214  if ( fVerbose > 0 )
215  cout << "=============== EVENT " << fTNofEvents << " =================" << endl;
216 
217  //Int_t nofCreatedTracks = 0; //[R.K. 01/2017] unused variable?
218 
219  Int_t nofHits = 0;
220  Int_t firstDH[4];
221 
222  // putting hits into one common array
223  for ( Int_t idet = 0 ; idet < 4 ; idet++ ) {
224  firstDH[idet] = nofHits;
225  if ( fIncludeDet[idet] == kTRUE ) {
226  nofHits += fHitArray[idet]->GetEntriesFast();
227  // cout << fDetName[idet].Data() << " hits from " << firstDH[idet] << " to " << nofHits << " ( " << fHitArray[idet]->GetEntriesFast() << " hits )" << endl;
228  for ( Int_t ihit = firstDH[idet] ; ihit < nofHits ; ihit++ ) { // hit mixing could go here
229  fHitDetId.push_back(idet);
230  fHitDetNo.push_back(ihit-firstDH[idet]);
231  }
232  }
233  }
234  std::vector<Bool_t> hitUsed(fHitDetId.size(),kFALSE);
235 
236  FairHit* detHit;
237  // DRAW
238 
239  if ( fCanvas ) {
240  fCanvas->Clear();
241  Double_t minX = -60.;
242  Double_t maxX = +60.;
243  Double_t minY = -60.;
244  Double_t maxY = +60.;
245 
246  fCanvas->Range(minX,minY,maxX,maxY);
247  fCanvas->SetFillColor(0);
248  fCanvas->SetBorderSize(0);
249  fCanvas->SetBorderMode(0);
250  fCanvas->SetFrameFillColor(0);
251  fCanvas->SetGrid(0,0);
252  fCanvas->SetLogy(0);
253 
254  // middle cross
255  TLine* horLft = new TLine(-5., 0.,-2., 0.); horLft->SetLineWidth(2); horLft->Draw();
256  TLine* horRgt = new TLine( 2., 0., 5., 0.); horRgt->SetLineWidth(2); horRgt->Draw();
257  TLine* verTop = new TLine( 0., 2., 0., 5.); verTop->SetLineWidth(2); verTop->Draw();
258  TLine* verBot = new TLine( 0.,-5., 0.,-2.); verBot->SetLineWidth(2); verBot->Draw();
259  for ( Int_t xline = -100 ; xline <= 100 ; xline+=10 ) {
260  TLine* verGrid = new TLine(xline,-100,xline,100);
261  verGrid->SetLineStyle(2);
262  verGrid->Draw();
263  }
264  for ( Int_t yline = -100 ; yline <= 100 ; yline+=10 ) {
265  TLine* horGrid = new TLine(-100,yline,100,yline);
266  horGrid->SetLineStyle(2);
267  horGrid->Draw();
268  }
269 
270  for ( Int_t iell = 0 ; iell < 1000 ; iell++ ) {
271  fEllipse[iell]->SetX1(100);
272  fEllipse[iell]->SetY1(100);
273  fEllipse[iell]->SetR1(0);
274  fEllipse[iell]->SetR2(0);
275  fEllipse[iell]->Draw();
276  fStageDesc->Draw();
277  }
278 
279  for ( Int_t ihit = 0 ; ihit < nofHits ; ihit++ ) {
280 
281  if ( fHitDetId[ihit] == 0 || fHitDetId[ihit] == 1 ) {
282  detHit = (PndSdsHit*)fHitArray[fHitDetId[ihit]]->At(fHitDetNo[ihit]);
283  TMarker* marker = new TMarker(detHit->GetX(),detHit->GetY(),20);
284  marker->SetMarkerColor(3);
285  marker->Draw();
286  }
287  else if ( fHitDetId[ihit] == 2 ) {
288  detHit = (PndSttHit*)fHitArray[fHitDetId[ihit]]->At(fHitDetNo[ihit]);
289  Int_t iTube = ((PndSttHit*)detHit)->GetTubeID();
290  PndSttTube *sttTube = nullptr;
291  if (fGeoType == 1)
292  sttTube = (PndSttTube*) fTubeArray->At(iTube);
293  else if (fGeoType == 2){
294  sttTube = PndSttTubeMap::Instance()->GetTube(iTube);
295  }
296 
297  if ( sttTube->GetWireDirection().Z() < 1. )
298  continue; // don't draw skewed hits (how else to put them?)
299 
300  TEllipse* ellipse = new TEllipse(detHit->GetX(),detHit->GetY(),((PndSttHit*)detHit)->GetIsochrone());
301  ellipse->SetFillStyle(0);
302  ellipse->SetLineWidth(2);
303  ellipse->SetLineColor(2);
304  ellipse->Draw();
305  }
306  else if ( fHitDetId[ihit] == 3 ) {
307  detHit = (PndGemHit*)fHitArray[fHitDetId[ihit]]->At(fHitDetNo[ihit]);
308  TMarker* marker = new TMarker(detHit->GetX(),detHit->GetY(),20);
309  marker->SetMarkerColor(6);
310  marker->Draw();
311  }
312  }
313  fStageDesc->SetText(-50.,50.,"HITS");
314  fStageDesc->Paint();
315  fCanvas->Modified();
316  fCanvas->Update();
317  gSystem->Sleep(1000);
318  }
319 
320  Int_t hitN = 0;
321  // looping over hits as they are in the common array
322 
323  for ( Int_t ihit = 0 ; ihit < nofHits ; ihit++ ) {
324 
325  hitN = (hitN+gRandom->Integer(nofHits))%nofHits;
326  while ( 1==1 ) {
327  if ( hitUsed[hitN] == kFALSE ) break;
328  hitN = (hitN+1)%nofHits;
329  }
330  hitUsed[hitN] = kTRUE;
331  // cout << " for hit " << ihit << " would like to take hit " << hitN << " >> " << fHitDetId[hitN] << "." << fHitDetNo[hitN] << endl;
332 
333  // hitN = (firstDH[2]+ihit)%nofHits; // hit mixing could go here
334 
335  if ( fHitDetId[hitN] == 0 || fHitDetId[hitN] == 1 )
336  detHit = (PndSdsHit*)fHitArray[fHitDetId[hitN]]->At(fHitDetNo[hitN]);
337  else if ( fHitDetId[hitN] == 2 )
338  detHit = (PndSttHit*)fHitArray[fHitDetId[hitN]]->At(fHitDetNo[hitN]);
339  else if ( fHitDetId[hitN] == 3 )
340  detHit = (PndGemHit*)fHitArray[fHitDetId[hitN]]->At(fHitDetNo[hitN]);
341 
342  if ( fVerbose > 3 ) {
343  cout << "--------------------------------------------------------" << endl;
344  cout << ihit << " (" << hitN << ") > " << fHitDetId[hitN] << "." << fHitDetNo[hitN] << " in " << fDetName[fHitDetId[hitN]].Data() << " at "
345  << detHit->GetX() << " " << detHit->GetY() << " " << detHit->GetZ()
346  << " error " << detHit->GetDx() << " " << detHit->GetDy() << " " << detHit->GetDz() << endl;
347  }
348  if ( fVerbose > 3 )
349  PrintTracks();
350 
351  Bool_t sttHit = kFALSE;
352  Bool_t skewedSttHit = kFALSE;
353  if ( fHitDetId[hitN] == 2 ) {
354  sttHit = kTRUE;
355  Int_t iTube = ((PndSttHit*)detHit)->GetTubeID();
356  PndSttTube *sttTube = nullptr;
357  if (fGeoType == 1)
358  sttTube = (PndSttTube*) fTubeArray->At(iTube);
359  else if (fGeoType == 2){
360  sttTube = PndSttTubeMap::Instance()->GetTube(iTube);
361  }
362  if (sttTube == nullptr)
363  std::cout << "-E- PndBarrelTrackFinder::Exec no sttTube " << iTube << std::endl;
364  if ( sttTube->GetWireDirection().Z() < 1. ) {
365  skewedSttHit = kTRUE;
366  if ( fVerbose > 3 )
367  cout << "SKEWED HIT" << endl;
368  }
369  }
370 
371  Bool_t hitBelongsToKnownTrack = kFALSE;
372  for ( size_t itr = 0 ; itr < fTracksVector.size() ; itr++ ) {
373  if ( !sttHit ) { if ( MatchHitWithZInfoTT (detHit,fHitDetId[hitN],fHitDetNo[hitN],itr) ) hitBelongsToKnownTrack = kTRUE; }
374  else
375  if ( skewedSttHit ) { if ( MatchSkewedSttHitTT (detHit,fHitDetId[hitN],fHitDetNo[hitN],itr) ) hitBelongsToKnownTrack = kTRUE; }
376  else { if ( MatchParallelSttHitTT(detHit,fHitDetId[hitN],fHitDetNo[hitN],itr) ) hitBelongsToKnownTrack = kTRUE; }
377  }
378  if ( hitBelongsToKnownTrack ) continue;
379 
380  if ( skewedSttHit ) {
381  AddHitToPreviousHits (detHit,fHitDetId[hitN],fHitDetNo[hitN]);
382  continue;
383  }
384 
385  Bool_t hitMatchedToHit = kFALSE;
386  if ( fVerbose > 3 )
387  cout << "TRYING TO MATCH TO PREVIOUS " << fHitVector.size() << " HITS" << endl;
388  for ( Int_t ihitP = fHitVector.size()-1 ; ihitP >= 0 ; ihitP-- ) {
389  if ( MatchHitToHit(detHit,fHitDetId[hitN],fHitDetNo[hitN],ihitP) ) hitMatchedToHit = kTRUE;
390  }
391  if ( hitMatchedToHit ) continue;
392 
393  AddHitToPreviousHits (detHit,fHitDetId[hitN],fHitDetNo[hitN]);
394  }
395  if ( fCanvas ) {
396  DrawTracks();
397  fStageDesc->SetText(-50.,50.,"TRACKS FOUND");
398  fStageDesc->Paint();
399  fCanvas->Modified();
400  fCanvas->Update();
401  gSystem->Sleep(1000);
402  }
403 // PrintTracks();
404 // cout << "REMOVING TRACKS!" << endl;
405  RemoveShortTracks();
406  // PrintTracks();
407  if ( fVerbose > 0 )
408  PrintTracks();
409  if ( fCanvas ) {
410  DrawTracks();
411  fStageDesc->SetText(-50.,50.,"SHORT TRACKS REMOVED");
412  fStageDesc->Paint();
413  fCanvas->Modified();
414  fCanvas->Update();
415  gSystem->Sleep(1000);
416  }
417 
418  // deal with left previous hits
419  for ( Int_t iuh = fHitVector.size()-1 ; iuh >= 0 ; iuh-- ) {
420  if ( fVerbose > 3 )
421  cout << "PREVIOUS HIT " << iuh << " out of " << fHitVector.size() << endl;
422  Bool_t sttHit = kFALSE;
423  Bool_t skewedSttHit = kFALSE;
424  if ( fHitVectDI[iuh] == 2 ) {
425  sttHit = kTRUE;
426  Int_t iTube = ((PndSttHit*)fHitVector[iuh])->GetTubeID();
427  PndSttTube *sttTube = nullptr;
428  if (fGeoType == 1)
429  sttTube = (PndSttTube*) fTubeArray->At(iTube);
430  else if (fGeoType == 2){
431  sttTube = PndSttTubeMap::Instance()->GetTube(iTube);
432  }
433  if ( sttTube->GetWireDirection().Z() < 1. ) {
434  skewedSttHit = kTRUE;
435  }
436  }
437 
438  Bool_t hitBelongsToKnownTrack = kFALSE;
439  for ( size_t itr = 0 ; itr < fTracksVector.size() ; itr++ ) {
440  if ( HitBelongsToTrack(fHitVectDI[iuh],fHitVectHN[iuh],itr) ) continue;
441  if ( !sttHit ) { if ( MatchHitWithZInfoTT (fHitVector[iuh],fHitVectDI[iuh],fHitVectHN[iuh],itr) ) hitBelongsToKnownTrack = kTRUE; }
442  else
443  if ( skewedSttHit ) { if ( MatchSkewedSttHitTT (fHitVector[iuh],fHitVectDI[iuh],fHitVectHN[iuh],itr) ) hitBelongsToKnownTrack = kTRUE; }
444  else { if ( MatchParallelSttHitTT(fHitVector[iuh],fHitVectDI[iuh],fHitVectHN[iuh],itr) ) hitBelongsToKnownTrack = kTRUE; }
445  }
446  if ( hitBelongsToKnownTrack ) {
447  RemoveHitFromPreviousHits(iuh);
448  continue;
449  }
450 
451  }
452 
453  // cout << "==============================================================" << endl;
454  // cout << "==============================================================" << endl;
455  // cout << "there are still " << fHitVector.size() << " unused hits: " << flush;
456  // for ( Int_t iuh = 0 ; iuh < fHitVector.size() ; iuh++ )
457  // cout << fHitVectDI[iuh] << "." << fHitVectHN[iuh] << " " << flush;
458  // cout << endl;
459  // cout << "==============================================================" << endl;
460  // cout << "==============================================================" << endl;
461  // cout << "FOUND " << fTracksVector.size() << " tracks:" << endl;
462  if ( fVerbose > 0 )
463  PrintTracks();
464  if ( fCanvas ) {
465  DrawTracks();
466  fStageDesc->SetText(-50.,50.,"TRACKS REORGANIZED");
467  fStageDesc->Paint();
468  fCanvas->Modified();
469  fCanvas->Update();
470  gSystem->Sleep(1000);
471  }
472  // cout << "==============================================================" << endl;
473  // cout << "REMOVE CLONE TRACKS" << endl;
474  RemoveCloneTracks();
475  if ( fVerbose > 0 )
476  PrintTracks();
477  if ( fCanvas ) {
478  DrawTracks();
479  fStageDesc->SetText(-50.,50.,"CLONE TRACKS REMOVED");
480  fStageDesc->Paint();
481  fCanvas->Modified();
482  fCanvas->Update();
483  gSystem->Sleep(5000);
484  }
485 
486  WriteTracks();
487 
488  fTNofEvents++;
489  fTNofTracks += fTracksVector.size();
490 
491  fTracksVector.clear();
492  fHitDetId.clear();
493  fHitDetNo.clear();
494  fHitVector.clear();
495  fHitVectDI.clear();
496  fHitVectHN.clear();
497 
498  return;
499 }
500 // -------------------------------------------------------------------------
501 
502 // ----- Private method MatchHitToTrack -------------------------------
503 Bool_t PndBarrelTrackFinder::MatchHitToTrack (FairHit* thisHit, Int_t detId, Int_t hitNo, Int_t trackNo) {
504  if ( thisHit ) cout << detId << "." << hitNo << " to " << trackNo << endl;
505  return kFALSE;
506 }
507 // -------------------------------------------------------------------------
508 
509 // ----- Private method HitBelongsToTrack -------------------------------
510 Bool_t PndBarrelTrackFinder::HitBelongsToTrack(Int_t detId, Int_t hitNo, Int_t trackNo) {
511  for ( size_t ihitP = 0 ; ihitP < fTracksVector[trackNo].trackHitD.size() ; ihitP++ ) {
512  if ( fTracksVector[trackNo].trackHitD[ihitP] == detId &&
513  fTracksVector[trackNo].trackHitN[ihitP] == hitNo )
514  return kTRUE;
515  }
516  return kFALSE;
517 }
518 // -------------------------------------------------------------------------
519 
520 // ----- Private method MatchSkewedSttHitTT -------------------------------
521 Bool_t PndBarrelTrackFinder::MatchSkewedSttHitTT (FairHit* thisHit, Int_t detId, Int_t hitNo, Int_t trackNo) {
522  Bool_t printInfo = kFALSE;//kTRUE;
523 // if ( fVerbose > 4 )
524 // cout << "skewed hit" << endl;
525  if ( TMath::Abs(fTracksVector[trackNo].meanR) < 1. ) return kFALSE;
526 
527  Int_t iTube = ((PndSttHit*)thisHit)->GetTubeID();
528  PndSttTube *sttTube = nullptr;
529  if (fGeoType == 1)
530  sttTube = (PndSttTube*) fTubeArray->At(iTube);
531  else if (fGeoType == 2){
532  sttTube = PndSttTubeMap::Instance()->GetTube(iTube);
533  }
534 
535  Double_t ra = TMath::ATan(TMath::Sqrt(sttTube->GetWireDirection().X()*sttTube->GetWireDirection().X()+
536  sttTube->GetWireDirection().Y()*sttTube->GetWireDirection().Y())/
537  sttTube->GetWireDirection().Z());
538  Double_t a = 0.5*((PndSttHit*)thisHit)->GetIsochrone()*(1.+1./TMath::Cos(ra));
539 
540  Double_t tubePar[8];
541  tubePar[0] = sttTube->GetPosition().X();
542  tubePar[1] = sttTube->GetPosition().Y();
543  tubePar[2] = sttTube->GetPosition().Z();
544  tubePar[3] = sttTube->GetHalfLength();
545  tubePar[4] = sttTube->GetWireDirection().X();
546  tubePar[5] = sttTube->GetWireDirection().Y();
547  tubePar[6] = sttTube->GetWireDirection().Z();
548  tubePar[7] = a;
549 
550  if ( fVerbose > 4 || printInfo )
551  cout << "matching to track " << trackNo << " with " << fTracksVector[trackNo].trackHits.size() << " hits" << endl;
552 
553  Bool_t trackFits = kFALSE;
554 
555  Double_t meanPhi = CalcPhi(fTracksVector[trackNo].meanX,fTracksVector[trackNo].meanY);
556  Double_t meanCirc[4] = {fTracksVector[trackNo].meanX,
557  fTracksVector[trackNo].meanY,
558  fTracksVector[trackNo].meanR,
559  meanPhi};
560  Double_t intReg[4];
561 
562  Int_t nofReg = FindInterestingRegions(meanCirc,tubePar,intReg);
563  if ( nofReg==0 ) return kFALSE;
564 
565  if ( fVerbose > 4 || printInfo ) {
566  cout << " --> connects to " << trackNo << " with ("
567  << meanPhi << " "
568  << meanCirc[2] << " "
569  << fTracksVector[trackNo].meanZ_P << ") at z_p = " << flush;
570  for ( Int_t ireg = 0 ; ireg < nofReg ; ireg++ )
571  cout << intReg[ireg] << " " << flush;
572  cout << endl;
573  }
574 
575  // there is no mean z_p, attach all we can
576  if ( fTracksVector[trackNo].meanZ_P > 1.e6 ) {
577  for ( Int_t ireg = 0 ; ireg < nofReg ; ireg++ ) {
578  TrackParameter tp0;
579  tp0.x = meanCirc[0];
580  tp0.y = meanCirc[1];
581  tp0.r = meanCirc[2];
582  tp0.z_p = intReg[ireg];
583  tp0.n = 400/nofReg;//fTracksVector[trackNo].trackHits.size();
584  fTracksVector[trackNo].trackPars.push_back(tp0);
585  if ( fVerbose > 4 || printInfo )
586  cout << "creating another trackPar to track " << trackNo << " with z_p -> " << intReg[ireg] << endl;
587  }
588  fTracksVector[trackNo].trackHits.push_back(thisHit);
589  fTracksVector[trackNo].trackHitD.push_back(detId);
590  fTracksVector[trackNo].trackHitN.push_back(hitNo);
591 
592  ExtractMeanZ_PFromTrack(trackNo);
593  }
594  // there is already mean z_p !!!, check if any is good
595  else {
596  // cout << "should compare if any of new found z_p: " << flush;
597  for ( Int_t ireg = 0 ; ireg < nofReg ; ireg++ ) {
598  //cout << intReg[ireg] << " " << flush;
599  if ( TMath::Abs(intReg[ireg]-fTracksVector[trackNo].meanZ_P) < fReasonableZ_PD )
600  trackFits = kTRUE;
601  }
602  // cout << " is close enough to mean: " << fTracksVector[trackNo].meanZ_P << endl;
603  if ( !trackFits ) return kFALSE;
604 
605  fTracksVector[trackNo].trackHits.push_back(thisHit);
606  fTracksVector[trackNo].trackHitD.push_back(detId);
607  fTracksVector[trackNo].trackHitN.push_back(hitNo);
608  }
609 
610  return kTRUE;
611 }
612 // -------------------------------------------------------------------------
613 
614 // ----- Private method MatchParallelSttHitTT -------------------------------
615 Bool_t PndBarrelTrackFinder::MatchParallelSttHitTT (FairHit* thisHit, Int_t detId, Int_t hitNo, Int_t trackNo) {
616  Bool_t hitBelongsToTrack = kFALSE;
617 
618  Double_t sH1[3] = {thisHit->GetX(),thisHit->GetY(),0.};
619  sH1[2] = ((PndSttHit*)thisHit)->GetIsochrone();
620  Double_t sH2[3]; // in the loop over track previous hits
621  Double_t sH3[3] = {0.,0.,0.}; // set to origin now, should be another hit (double for loop over hits)
622  Double_t sH4[3]; // another loop over hits
623  Double_t circPar[3];
624 
625  Double_t thisHitPhi = CalcPhi(sH1[0],sH1[1]);
626  Double_t thisTrackPhi = CalcPhi(fTracksVector[trackNo].trackHits[0]->GetX(),fTracksVector[trackNo].trackHits[0]->GetY());
627  Double_t phiHTDiffr = TMath::Abs(thisHitPhi-thisTrackPhi);
628  if ( phiHTDiffr > TMath::Pi()/4. && phiHTDiffr < TMath::Pi()*7./4. ) return kFALSE; // checking if hit and track match approximately, should make it radius dependent
629 
630  if ( TMath::Abs(fTracksVector[trackNo].meanR) > 1. ) {
631  circPar[0] = fTracksVector[trackNo].meanX;
632  circPar[1] = fTracksVector[trackNo].meanY;
633  circPar[2] = fTracksVector[trackNo].meanR;
634 
635  if ( FindCircDist(circPar,sH1) < fReasonableDist )
636  hitBelongsToTrack = kTRUE;
637  if ( !hitBelongsToTrack ) return kFALSE;
638 
639  // putting thisHit into track
640  fTracksVector[trackNo].trackHits.push_back(thisHit);
641  fTracksVector[trackNo].trackHitD.push_back(detId);
642  fTracksVector[trackNo].trackHitN.push_back(hitNo);
643 
644  if ( fVerbose > 3 )
645  cout << "this stt hit belongs to track " << trackNo << " (cause dist = " << FindCircDist(circPar,sH1) << ")" << endl;
646  }
647  else {
648  // track has no mean information yet
649  for ( size_t ipar = 0 ; ipar < fTracksVector[trackNo].trackPars.size() ; ipar++ ) {
650  if ( fTracksVector[trackNo].trackPars[ipar].n < 0.9*fTracksVector[trackNo].trackHits.size() ) continue;
651 
652  circPar[0] = fTracksVector[trackNo].trackPars[ipar].x;
653  circPar[1] = fTracksVector[trackNo].trackPars[ipar].y;
654  circPar[2] = fTracksVector[trackNo].trackPars[ipar].r;
655 
656  if ( fVerbose > 4 ) {
657  cout << " ..... to parameter (" << trackNo << "." << ipar << ") " << CalcPhi(circPar[0],circPar[1]) << " " << circPar[2] << ": distance is " << FindCircDist(sH1,circPar) << flush;
658  cout << endl;
659  }
660 
661  if ( FindCircDist(sH1,circPar) < fMaximalDist ) {
662  hitBelongsToTrack = kTRUE;
663 
664  fTracksVector[trackNo].trackPars[ipar].n += 1; // because thisHit is close to parameter ipar
665 
666  if ( fVerbose > 4 )
667  cout << " par: " << CalcPhi(circPar[0],circPar[1]) << " | " << circPar[2] << " of track " << trackNo
668  << " matches to this hit with dist " << FindCircDist(sH1,circPar) << endl;
669  }
670  }
671  if ( !hitBelongsToTrack ) return kFALSE;
672 
673  // calculating track parameters with previous hits belonging to the track
674  for ( size_t ihitP = 0 ; ihitP < fTracksVector[trackNo].trackHitD.size() ; ihitP++ ) {
675  sH2[0] = fTracksVector[trackNo].trackHits[ihitP]->GetX();
676  sH2[1] = fTracksVector[trackNo].trackHits[ihitP]->GetY();
677  sH2[2] = 0.;
678  if ( fTracksVector[trackNo].trackHitD[ihitP] == 2 )
679  sH2[2] = ((PndSttHit*)fTracksVector[trackNo].trackHits[ihitP])->GetIsochrone();
680 
681  Int_t stepCno = 1;
682  if ( sH2[2] < 1.e-5 ) stepCno += 2;
683 
684  for ( Int_t icirc = 0 ; icirc < 4 ; icirc += stepCno ) {
685  if ( FindCircPar(sH1,sH2,sH3,icirc,circPar) ) {
686  if ( circPar[2] < 10. || circPar[2] > 100000. ) continue; // for now remove too egzotic tracks
687 
688  Double_t newParZ_P = 2.e6;
689  if ( sH2[2] < 1.e-5 ) newParZ_P = CalcZ_P(circPar,sH2[0],sH2[1],fTracksVector[trackNo].trackHits[ihitP]->GetZ());
690 
691  Int_t nofCloseHits = 2;
692  for ( size_t ihitF = 0 ; ihitF < fTracksVector[trackNo].trackHitD.size() ; ihitF++ ) {
693  if ( ihitF == ihitP ) continue;
694  sH4[0] = fTracksVector[trackNo].trackHits[ihitF]->GetX();
695  sH4[1] = fTracksVector[trackNo].trackHits[ihitF]->GetY();
696  sH4[2] = 0.;
697  if ( fTracksVector[trackNo].trackHitD[ihitF] == 2 )
698  sH4[2] = ((PndSttHit*)fTracksVector[trackNo].trackHits[ihitF])->GetIsochrone();
699 
700  if ( FindCircDist(sH4,circPar) < fMaximalDist || // if the hit is very very close to the new parameter in xy plane
701  ( newParZ_P < 1.e6 && FindCircDist(sH4,circPar) < fReasonableDist && // if the hit is close to the new parameter in xy plane and pz space
702  TMath::Abs(CalcZ_P(circPar,sH4[0],sH4[1],fTracksVector[trackNo].trackHits[ihitF]->GetZ())-newParZ_P) < fMaximalZ_PD ) ) {
703  nofCloseHits += 1;
704  }
705  }
706 
707  if ( fVerbose > 5 )
708  cout << "there are " << nofCloseHits << " with limit set to " << 0.6*fTracksVector[trackNo].trackHits.size() << " (0.6*" << fTracksVector[trackNo].trackHits.size() << ")" << endl;
709  if ( nofCloseHits < 0.6*fTracksVector[trackNo].trackHits.size() ) continue; // the track param has not enough close hits among previous track hits
710 
711  TrackParameter tp0;
712  tp0.x = circPar[0];
713  tp0.y = circPar[1];
714  tp0.r = circPar[2];
715  tp0.z_p = newParZ_P;
716  tp0.n = nofCloseHits;
717 
718  fTracksVector[trackNo].trackPars.push_back(tp0);
719 
720  }
721  }
722  }
723 
724  // putting thisHit into track
725  fTracksVector[trackNo].trackHits.push_back(thisHit);
726  fTracksVector[trackNo].trackHitD.push_back(detId);
727  fTracksVector[trackNo].trackHitN.push_back(hitNo);
728 
729  if ( fVerbose > 5 ) {
730  cout << "HIT DOES " << (hitBelongsToTrack?"":"NOT ") << "BELONG TO TRACK " << trackNo << " WITH "
731  << fTracksVector[trackNo].trackHits.size() << " HITS, " << fTracksVector[trackNo].trackPars.size() << " PARS: phi,rad,z_p,n: " << endl;
732  for ( size_t ipar = 0 ; ipar < fTracksVector[trackNo].trackPars.size() ; ipar++ )
733  cout << CalcPhi(fTracksVector[trackNo].trackPars[ipar].x,fTracksVector[trackNo].trackPars[ipar].y) << " "
734  << fTracksVector[trackNo].trackPars[ipar].r << " "
735  << fTracksVector[trackNo].trackPars[ipar].z_p << " "
736  << fTracksVector[trackNo].trackPars[ipar].n << " | " << flush;
737  cout << endl;
738  }
739 
740  ExtractMeanRPhiFromTrack(trackNo);
741  }
742 
743  return kTRUE;
744 }
745 // -------------------------------------------------------------------------
746 
747 // ----- Private method MatchHitWithZInfoTT -------------------------------
748 Bool_t PndBarrelTrackFinder::MatchHitWithZInfoTT (FairHit* thisHit, Int_t detId, Int_t hitNo, Int_t trackNo) {
749  Bool_t hitBelongsToTrack = kFALSE;
750 
751  Double_t sH1[3] = {thisHit->GetX(),thisHit->GetY(),0.};
752  Double_t sH2[3]; // in the loop over track previous hits
753  Double_t sH3[3] = {0.,0.,0.}; // set to origin now, should be another hit (double for loop over hits)
754  Double_t sH4[3]; // another loop over hits
755  Double_t circPar[3];
756 
757  Double_t thisHitPhi = CalcPhi(sH1[0],sH1[1]);
758  Double_t thisTrackPhi = CalcPhi(fTracksVector[trackNo].trackHits[0]->GetX(),fTracksVector[trackNo].trackHits[0]->GetY());
759  Double_t phiHTDiffr = TMath::Abs(thisHitPhi-thisTrackPhi);
760  if ( phiHTDiffr > TMath::Pi()/4. && phiHTDiffr < TMath::Pi()*7./4. ) return kFALSE; // checking if hit and track match approximately, should make it radius dependent
761 
762  if ( TMath::Abs(fTracksVector[trackNo].meanR) > 1. ) {
763  circPar[0] = fTracksVector[trackNo].meanX;
764  circPar[1] = fTracksVector[trackNo].meanY;
765  circPar[2] = fTracksVector[trackNo].meanR;
766 
767  if ( fVerbose > 3 )
768  cout << "hit " << detId << "." << hitNo << " distance to track " << trackNo << " is " << FindCircDist(circPar,sH1) << " with errors ( "
769  << ((PndGemHit*)thisHit)->GetDr() << " , " << ((PndGemHit*)thisHit)->GetDp() << " ) " << endl;
770 
771  if ( FindCircDist(circPar,sH1) < fReasonableDist )
772  hitBelongsToTrack = kTRUE;
773  if ( !hitBelongsToTrack ) return kFALSE;
774 
775  if ( fTracksVector[trackNo].meanZ_P > 1.e6 ) {
776  TrackParameter tp0;
777  tp0.x = circPar[0];
778  tp0.y = circPar[1];
779  tp0.r = circPar[2];
780  tp0.z_p = CalcZ_P(circPar,sH1[0],sH1[1],thisHit->GetZ());
781  tp0.n = 400;
782  fTracksVector[trackNo].trackPars.push_back(tp0);
783 
784  // putting thisHit into track
785  fTracksVector[trackNo].trackHits.push_back(thisHit);
786  fTracksVector[trackNo].trackHitD.push_back(detId);
787  fTracksVector[trackNo].trackHitN.push_back(hitNo);
788  ExtractMeanZ_PFromTrack(trackNo);
789  }
790  else {
791  // cout << "should compare if new found z_p " << CalcZ_P(circPar,sH1[0],sH1[1],thisHit->GetZ()) << " is close enough to mean: " << fTracksVector[trackNo].meanZ_P << endl;
792 
793  if ( TMath::Abs(CalcZ_P(circPar,sH1[0],sH1[1],thisHit->GetZ())-fTracksVector[trackNo].meanZ_P) < fReasonableZ_PD ) {
794  // putting thisHit into track
795  fTracksVector[trackNo].trackHits.push_back(thisHit);
796  fTracksVector[trackNo].trackHitD.push_back(detId);
797  fTracksVector[trackNo].trackHitN.push_back(hitNo);
798  }
799  else
800  return kFALSE;
801  }
802  if ( fVerbose > 3 )
803  cout << " This hit belongs to track " << trackNo << " with " << CalcZ_P(circPar,sH1[0],sH1[1],thisHit->GetZ()) << endl;
804  }
805  else {
806  // track has no mean information yet
807  for ( size_t ipar = 0 ; ipar < fTracksVector[trackNo].trackPars.size() ; ipar++ ) {
808  if ( fTracksVector[trackNo].trackPars[ipar].n < 0.6*fTracksVector[trackNo].trackHits.size() ) continue;
809 
810  circPar[0] = fTracksVector[trackNo].trackPars[ipar].x;
811  circPar[1] = fTracksVector[trackNo].trackPars[ipar].y;
812  circPar[2] = fTracksVector[trackNo].trackPars[ipar].r;
813 
814  if ( fVerbose > 4 ) {
815  cout << " ..... to parameter (" << trackNo << "." << ipar << ") " << CalcPhi(circPar[0],circPar[1]) << " " << circPar[2] << ": distance is " << FindCircDist(sH1,circPar) << flush;
816  cout << " new z_p is " << CalcZ_P(circPar,sH1[0],sH1[1],thisHit->GetZ()) << " vs " << fTracksVector[trackNo].trackPars[ipar].z_p << " old" << flush;
817  cout << endl;
818  }
819 
820  if ( ( FindCircDist(sH1,circPar) < fMaximalDist &&
821  TMath::Abs(CalcZ_P(circPar,sH1[0],sH1[1],thisHit->GetZ())-fTracksVector[trackNo].trackPars[ipar].z_p) < 3.*fReasonableZ_PD ) ||
822  ( FindCircDist(sH1,circPar) < 3.*fReasonableDist &&
823  TMath::Abs(CalcZ_P(circPar,sH1[0],sH1[1],thisHit->GetZ())-fTracksVector[trackNo].trackPars[ipar].z_p) < fMaximalZ_PD ) ) {
824  hitBelongsToTrack = kTRUE;
825 
826  fTracksVector[trackNo].trackPars[ipar].n += 1; // because thisHit is close to parameter ipar
827 
828  if ( fVerbose > 4 )
829  cout << " par: " << CalcPhi(circPar[0],circPar[1]) << " | " << circPar[2] << " of track " << trackNo
830  << " matches to this hit with dist " << FindCircDist(sH1,circPar) << endl;
831  }
832  }
833  if ( !hitBelongsToTrack ) return kFALSE;
834 
835  // calculating track parameters with previous hits belonging to the track
836  for ( size_t ihitP = 0 ; ihitP < fTracksVector[trackNo].trackHitD.size() ; ihitP++ ) {
837  sH2[0] = fTracksVector[trackNo].trackHits[ihitP]->GetX();
838  sH2[1] = fTracksVector[trackNo].trackHits[ihitP]->GetY();
839  sH2[2] = 0.;
840  if ( fTracksVector[trackNo].trackHitD[ihitP] == 2 )
841  sH2[2] = ((PndSttHit*)fTracksVector[trackNo].trackHits[ihitP])->GetIsochrone();
842 
843  Int_t stepCno = 2;
844  if ( sH2[2] < 1.e-5 ) stepCno += 2;
845 
846  for ( Int_t icirc = 0 ; icirc < 4 ; icirc += stepCno ) {
847  if ( FindCircPar(sH1,sH2,sH3,icirc,circPar) ) {
848  if ( circPar[2] < 10. || circPar[2] > 100000. ) continue; // for now remove too egzotic tracks
849 
850  Double_t newParZ_P = CalcZ_P(circPar,sH1[0],sH1[1],thisHit->GetZ());
851 
852  Int_t nofCloseHits = 2;
853  for ( size_t ihitF = 0 ; ihitF < fTracksVector[trackNo].trackHitD.size() ; ihitF++ ) {
854  if ( ihitF == ihitP ) continue;
855  sH4[0] = fTracksVector[trackNo].trackHits[ihitF]->GetX();
856  sH4[1] = fTracksVector[trackNo].trackHits[ihitF]->GetY();
857  sH4[2] = 0.;
858  if ( fTracksVector[trackNo].trackHitD[ihitF] == 2 )
859  sH4[2] = ((PndSttHit*)fTracksVector[trackNo].trackHits[ihitF])->GetIsochrone();
860 
861 // if ( FindCircDist(sH4,circPar) < fMaximalDist || // if the hit is very very close to the new parameter in xy plane
862 // ( FindCircDist(sH4,circPar) < fReasonableDist && // if the hit is close to the new parameter in xy plane and pz space
863 // TMath::Abs(CalcZ_P(circPar,sH4[0],sH4[1],fTracksVector[trackNo].trackHits[ihitF]->GetZ())-newParZ_P) < fMaximalZ_PD ) ) {
864  if ( ( FindCircDist(sH4,circPar) < fMaximalDist &&
865  TMath::Abs(CalcZ_P(circPar,sH4[0],sH4[1],fTracksVector[trackNo].trackHits[ihitF]->GetZ())-newParZ_P) < 3.*fReasonableZ_PD ) ||
866  ( FindCircDist(sH4,circPar) < 3.*fReasonableDist &&
867  TMath::Abs(CalcZ_P(circPar,sH4[0],sH4[1],fTracksVector[trackNo].trackHits[ihitF]->GetZ())-newParZ_P) < fMaximalZ_PD ) ) {
868  nofCloseHits += 1;
869  }
870  }
871 
872  if ( fVerbose > 5 )
873  cout << "there are " << nofCloseHits << " with limit set to " << 0.6*fTracksVector[trackNo].trackHits.size() << " (0.6*" << fTracksVector[trackNo].trackHits.size() << ")" << endl;
874  if ( nofCloseHits < 0.6*fTracksVector[trackNo].trackHits.size() ) continue; // the track param has not enough close hits among previous track hits
875 
876  TrackParameter tp0;
877  tp0.x = circPar[0];
878  tp0.y = circPar[1];
879  tp0.r = circPar[2];
880  tp0.z_p = newParZ_P;
881  tp0.n = nofCloseHits;
882 
883  fTracksVector[trackNo].trackPars.push_back(tp0);
884 
885  }
886  }
887  }
888 
889  // putting thisHit into track
890  fTracksVector[trackNo].trackHits.push_back(thisHit);
891  fTracksVector[trackNo].trackHitD.push_back(detId);
892  fTracksVector[trackNo].trackHitN.push_back(hitNo);
893 
894  if ( fVerbose > 5 ) {
895  cout << "HIT DOES " << (hitBelongsToTrack?"":"NOT ") << "BELONG TO TRACK " << trackNo << " WITH "
896  << fTracksVector[trackNo].trackHits.size() << " HITS, " << fTracksVector[trackNo].trackPars.size() << " PARS: phi,rad,z_p,n: " << endl;
897  for ( size_t ipar = 0 ; ipar < fTracksVector[trackNo].trackPars.size() ; ipar++ )
898  cout << CalcPhi(fTracksVector[trackNo].trackPars[ipar].x,fTracksVector[trackNo].trackPars[ipar].y) << " "
899  << fTracksVector[trackNo].trackPars[ipar].r << " "
900  << fTracksVector[trackNo].trackPars[ipar].z_p << " "
901  << fTracksVector[trackNo].trackPars[ipar].n << " | " << flush;
902  cout << endl;
903  }
904 
905  ExtractMeanRPhiFromTrack(trackNo);
906  }
907 
908  return kTRUE;
909 }
910 // -------------------------------------------------------------------------
911 
912 // ----- Private method MatchHitToHit -------------------------------
913 Bool_t PndBarrelTrackFinder::MatchHitToHit (FairHit* thisHit, Int_t detId, Int_t hitNo, Int_t prevHNo) {
914  Double_t sH1[3] = {thisHit->GetX(),thisHit->GetY(),0.};
915  if ( detId == 2 )
916  sH1[2] = ((PndSttHit*)thisHit)->GetIsochrone();
917 
918  Double_t sH2[3] = {fHitVector[prevHNo]->GetX(),fHitVector[prevHNo]->GetY(),0.};
919  if ( fHitVectDI[prevHNo] == 2 )
920  sH2[2] = ((PndSttHit*)fHitVector[prevHNo])->GetIsochrone();
921 
922  Double_t sH3[3] = {0.,0.,0.}; // set to origin now, should be another hit (double for loop over hits)
923  Double_t circPar[3];
924 
925  // cout << fHitVector[prevHNo]->GetZ() << " " << flush;
926  if ( fHitVectDI[prevHNo] == 2 ) {
927  Int_t prevTube = ((PndSttHit*)fHitVector[prevHNo])->GetTubeID();
928  PndSttTube *sttTube = nullptr;
929  if (fGeoType == 1)
930  sttTube = (PndSttTube*) fTubeArray->At(prevTube);
931  else if (fGeoType == 2){
932  sttTube = PndSttTubeMap::Instance()->GetTube(prevTube);
933  }
934  if ( sttTube->GetWireDirection().Z() < 1. )
935  return kFALSE;
936  sH2[2] = ((PndSttHit*)fHitVector[prevHNo])->GetIsochrone();
937  }
938  else {
939  if ( detId == fHitVectDI[prevHNo] )
940  if ( TMath::Abs(thisHit->GetZ()-fHitVector[prevHNo]->GetZ()) < 3. )
941  return kFALSE;
942  if ( detId != 2 ) {
943  Double_t thisRad = TMath::Sqrt(sH1[0]*sH1[0]+sH1[1]*sH1[1]);
944  Double_t prevRad = TMath::Sqrt(sH2[0]*sH2[0]+sH2[1]*sH2[1]);
945  if ( fVerbose > 3 ) {
946  cout << "matching r z " << thisRad << " " << thisHit->GetZ() << " (" << sH1[0] << " , " << sH1[1] << ") " << endl;
947  cout << " with r z " << prevRad << " " << fHitVector[prevHNo]->GetZ() << " (" << sH2[0] << " , " << sH2[1] << ") " << endl;
948  }
949  }
950  }
951 
952 
953  Bool_t hitsMatched = kFALSE;
954 
955  Double_t thisHitPhi = CalcPhi(sH1[0],sH1[1]);
956  Double_t prevHitPhi = CalcPhi(sH2[0],sH2[1]);
957 
958  Double_t phiDiffr = TMath::Abs(thisHitPhi-prevHitPhi);
959  if ( phiDiffr > TMath::Pi()/10. && phiDiffr < TMath::Pi()*19./10. ) return kFALSE;
960 
961  Int_t stepCno = 1;
962  if ( sH1[2] < 1.e-5 ) stepCno = 2;
963  if ( sH2[2] < 1.e-5 ) stepCno += 2;
964 // cout << "stepCno = " << stepCno << " cause " << sH1[2] << " and " << sH2[2] << endl;
965 
966  Int_t fTN = -1;
967  for ( Int_t icirc = 0 ; icirc < 4 ; icirc += stepCno ) {
968  // cout << icirc << endl;
969  if ( FindCircPar(sH1,sH2,sH3,icirc,circPar) ) {
970  // cout << " ---> " << circPar[2] << endl;
971  if ( circPar[2] < 10. || circPar[2] > 100000. ) continue;
972 
973  TrackParameter tp0;
974  tp0.x = circPar[0];
975  tp0.y = circPar[1];
976  tp0.r = circPar[2];
977  tp0.z_p = 2.e6;
978  // temporary recognition of non-STT hit
979  // should take only parameters created using hit with VALID z information
980  if ( sH1[2] < 1.e-5 && sH2[2] < 1.e-5 ) {
981  if ( TMath::Abs(CalcZ_P(circPar,sH1[0],sH1[1],thisHit->GetZ()) - CalcZ_P(circPar,sH2[0],sH2[1],fHitVector[prevHNo]->GetZ())) > 3.*fReasonableZ_PD ) {
982  continue;
983  }
984  tp0.z_p = (CalcZ_P(circPar,sH1[0],sH1[1],thisHit->GetZ()) + CalcZ_P(circPar,sH2[0],sH2[1],fHitVector[prevHNo]->GetZ()))/2.;
985  }
986  else {
987  if ( sH1[2] < 1.e-5 ) tp0.z_p = CalcZ_P(circPar,sH1[0],sH1[1],thisHit->GetZ());
988  else if ( sH2[2] < 1.e-5 ) tp0.z_p = CalcZ_P(circPar,sH2[0],sH2[1],fHitVector[prevHNo]->GetZ());
989  }
990 
991  if ( fVerbose > 3 )
992  if ( fHitVectDI[prevHNo] != 2 && detId != 2 ) {
993  cout << " combining hit " << fHitVectDI[prevHNo] << "." << fHitVectHN[prevHNo]
994  << " with hit " << detId << "." << hitNo << " --> z_p(1) = "
995  << CalcZ_P(circPar,sH2[0],sH2[1],fHitVector[prevHNo]->GetZ()) << " --> z_p(2) = "
996  << CalcZ_P(circPar,sH1[0],sH1[1],thisHit->GetZ()) << endl;
997  }
998  // cout << "and new z_p = " << tp0.z_p << endl;
999 
1000  tp0.n = 2;
1001 
1002  if ( fTN == -1 ) {
1003  TrackBasis tb0;
1004  tb0.trackPars.push_back(tp0);
1005  tb0.trackHits.push_back(fHitVector[prevHNo]);
1006  tb0.trackHitD.push_back(fHitVectDI[prevHNo]);
1007  tb0.trackHitN.push_back(fHitVectHN[prevHNo]);
1008  tb0.trackHits.push_back(thisHit);
1009  tb0.trackHitD.push_back(detId);
1010  tb0.trackHitN.push_back(hitNo);
1011  tb0.meanX = 0.;
1012  tb0.meanY = 0.;
1013  tb0.meanR = 0.;
1014  tb0.meanZ_P = 2.e6;
1015 
1016  fTN = fTracksVector.size();
1017  fTracksVector.push_back(tb0);
1018 
1019  hitsMatched = kTRUE;
1020  // cout << "creating track " << fTracksVector.size() << " from hits "
1021  // << fHitVectDI[prevHNo] << "." << fHitVectHN[prevHNo] << " and "
1022  // << fHitDetId[hitN] << "." << fHitDetNo[hitN] << endl;
1023  if ( fVerbose > 3 ) {
1024  cout << "!!! created track " << fTracksVector.size()-1 << " with " << fHitVectDI[prevHNo] << "." << fHitVectHN[prevHNo] << endl;
1025  // if ( ihit == 54 ) {
1026  cout << " >> with " << CalcPhi(tp0.x,tp0.y) << " (" << tp0.x << "," << tp0.y << ") " << tp0.r << " " << tp0.z_p << endl; //"(" << CalcPhi(tp0.x-sH1[0],tp0.y-sH1[1]) << " on " << detHit->GetZ() << ")" << endl;
1027  // Double_t tmpx = sH1[0]-tp0.x;
1028  // Double_t tmpy = sH1[1]-tp0.y;
1029  // cout << " should be: " << sH1[0] << " " << sH1[1] << " ---> " << tmpx << " " << tmpy << endl;
1030  // Double_t rotA = TMath::Pi()-CalcPhi(tp0.x,tp0.y);
1031  // cout << "rotate by " << rotA*TMath::RadToDeg() << " deg" << endl;
1032  // Double_t tmpx2 = tmpx*TMath::Cos(rotA)-tmpy*TMath::Sin(rotA);
1033  // Double_t tmpy2 = tmpx*TMath::Sin(rotA)+tmpy*TMath::Cos(rotA);
1034  // cout << " and then : " << tmpx << " " << tmpy << " ---> " << tmpx2 << " " << tmpy2 << endl;
1035  // cout << " and finally: " << CalcPhi(tmpx2,tmpy2) << endl;
1036  }
1037  }
1038  else {
1039  fTracksVector[fTN].trackPars.push_back(tp0);
1040  if ( fVerbose > 3 ) {
1041  cout << " >> with " << CalcPhi(tp0.x,tp0.y) << " (" << tp0.x << "," << tp0.y << ") " << tp0.r << " " << tp0.z_p << "(" << CalcPhi(tp0.x-sH1[0],tp0.y-sH1[1]) << " on " << thisHit->GetZ() << ")" << endl;
1042  }
1043  }
1044  }
1045  }
1046 
1047  if ( hitsMatched ) RemoveHitFromPreviousHits(prevHNo);
1048  return hitsMatched;
1049 }
1050 // -------------------------------------------------------------------------
1051 
1052 // ----- Private method ExtractMeanRPhiFromTrack -------------------------------
1053 Bool_t PndBarrelTrackFinder::ExtractMeanRPhiFromTrack (Int_t trackNo) {
1054  if ( fVerbose > 3 )
1055  cout << "extracting mean rphi from track " << trackNo << " with " << fTracksVector[trackNo].trackHits.size() << " hits." << endl;
1056  if ( fTracksVector[trackNo].trackHits.size() < 4 ) return kFALSE;
1057  if ( TMath::Abs(fTracksVector[trackNo].meanR) > 0.1 ) return kTRUE;
1058 
1059  Double_t meanPhi = 0.;
1060  Double_t meanRad = 0.;
1061  Double_t goodPNR = 0.;
1062  for ( size_t ipar = 0 ; ipar < fTracksVector[trackNo].trackPars.size() ; ipar++ ) {
1063  if ( fTracksVector[trackNo].trackPars[ipar].n < 0.9*fTracksVector[trackNo].trackHitD.size() ) continue;
1064  meanPhi += CalcPhi(fTracksVector[trackNo].trackPars[ipar].x,fTracksVector[trackNo].trackPars[ipar].y);
1065  meanRad += fTracksVector[trackNo].trackPars[ipar].r;
1066  goodPNR += 1.;
1067  }
1068  if ( goodPNR < 0.5 ) return kFALSE;
1069  meanRad = meanRad/goodPNR;
1070  meanPhi = meanPhi/goodPNR;
1071  if ( fVerbose > 3 )
1072  cout << "mean rad phi is " << meanRad << " " << meanPhi << endl;
1073 
1074  Double_t sH1[3] = {meanRad*TMath::Cos(meanPhi),
1075  meanRad*TMath::Sin(meanPhi),
1076  meanRad};
1077  Double_t sH2[3] = {0.,0.,0.};
1078 
1079  Int_t nofCloseHits = 0;
1080  for ( size_t ihit = 0 ; ihit < fTracksVector[trackNo].trackHits.size() ; ihit++ ) {
1081  sH2[0] = fTracksVector[trackNo].trackHits[ihit]->GetX();
1082  sH2[1] = fTracksVector[trackNo].trackHits[ihit]->GetY();
1083  sH2[2] = 0.;
1084  if ( fTracksVector[trackNo].trackHitD[ihit] == 2 ) {
1085  sH2[2] = ((PndSttHit*)fTracksVector[trackNo].trackHits[ihit])->GetIsochrone();
1086  }
1087 
1088  if ( FindCircDist(sH1,sH2) < fReasonableDist ) {
1089  nofCloseHits++;
1090  }
1091  }
1092  if ( fVerbose > 3 )
1093  cout << "there are " << nofCloseHits << " close hits." << endl;
1094  if ( nofCloseHits < 0.7*fTracksVector[trackNo].trackHitD.size() ) return kFALSE;
1095 
1096  fTracksVector[trackNo].meanX = sH1[0];
1097  fTracksVector[trackNo].meanY = sH1[1];
1098  fTracksVector[trackNo].meanR = sH1[2];
1099 
1100  for ( size_t ipar = 0 ; ipar < fTracksVector[trackNo].trackPars.size() ; ipar++ ) {
1101  if ( fTracksVector[trackNo].trackPars[ipar].n < 0.9*fTracksVector[trackNo].trackHitD.size() ) continue;
1102  if ( fTracksVector[trackNo].trackPars[ipar].z_p > 1.e6 ) continue;
1103  fTracksVector[trackNo].trackPars[ipar].n = 400;
1104  }
1105 
1106  ExtractMeanZ_PFromTrack(trackNo);
1107 
1108  return kTRUE;
1109 }
1110 // -------------------------------------------------------------------------
1111 
1112 // ----- Private method ExtractMeanZ_PFromTrack -------------------------------
1113 Bool_t PndBarrelTrackFinder::ExtractMeanZ_PFromTrack (Int_t trackNo) {
1114  if ( fVerbose > 3 )
1115  cout << "extracting mean Z_P from track " << trackNo << " with " << fTracksVector[trackNo].trackHits.size() << " hits." << endl;
1116  if ( fTracksVector[trackNo].meanZ_P < 1.e6 ) return kTRUE;
1117 
1118  // count hits with z information
1119  Int_t nofZHits = 0;
1120  Double_t circPar[4] = {fTracksVector[trackNo].meanX,
1121  fTracksVector[trackNo].meanY,
1122  fTracksVector[trackNo].meanR,
1123  CalcPhi(fTracksVector[trackNo].meanX,fTracksVector[trackNo].meanY)};
1124 
1125  std::vector<Double_t> trackHitsZ_P;
1126 
1127  for ( size_t ihit = 0 ; ihit < fTracksVector[trackNo].trackHits.size() ; ihit++ ) {
1128  if ( fTracksVector[trackNo].trackHitD[ihit] != 2 ) {
1129  nofZHits++;
1130  trackHitsZ_P.push_back(CalcZ_P(circPar,
1131  fTracksVector[trackNo].trackHits[ihit]->GetX(),
1132  fTracksVector[trackNo].trackHits[ihit]->GetY(),
1133  fTracksVector[trackNo].trackHits[ihit]->GetZ()));
1134  }
1135  else {
1136  Int_t tubeId = ((PndSttHit*)fTracksVector[trackNo].trackHits[ihit])->GetTubeID();
1137  PndSttTube *sttTube = nullptr;
1138  if (fGeoType == 1)
1139  sttTube = (PndSttTube*) fTubeArray->At(tubeId);
1140  else if (fGeoType == 2){
1141  sttTube = PndSttTubeMap::Instance()->GetTube(tubeId);
1142  }
1143 
1144  if ( sttTube->GetWireDirection().Z() < 1. ) {
1145  nofZHits++;
1146 
1147  Double_t ra = TMath::ATan(TMath::Sqrt(sttTube->GetWireDirection().X()*sttTube->GetWireDirection().X()+
1148  sttTube->GetWireDirection().Y()*sttTube->GetWireDirection().Y())/
1149  sttTube->GetWireDirection().Z());
1150  Double_t a = 0.5*((PndSttHit*)fTracksVector[trackNo].trackHits[ihit])->GetIsochrone()*(1.+1./TMath::Cos(ra));
1151 
1152  Double_t tubePar[8];
1153  tubePar[0] = sttTube->GetPosition().X();
1154  tubePar[1] = sttTube->GetPosition().Y();
1155  tubePar[2] = sttTube->GetPosition().Z();
1156  tubePar[3] = sttTube->GetHalfLength();
1157  tubePar[4] = sttTube->GetWireDirection().X();
1158  tubePar[5] = sttTube->GetWireDirection().Y();
1159  tubePar[6] = sttTube->GetWireDirection().Z();
1160  tubePar[7] = a;
1161 
1162  Double_t intReg[4];
1163 
1164  Int_t nofReg = FindInterestingRegions(circPar,tubePar,intReg);
1165  for ( Int_t ireg = 0 ; ireg < nofReg ; ireg++ ) {
1166  trackHitsZ_P.push_back(intReg[ireg]);
1167  }
1168  }
1169  }
1170  }
1171  if ( nofZHits < 4 ) return kFALSE;
1172 
1173  if ( fVerbose > 3 )
1174  cout << "trying to etract mean z_p. There are " << nofZHits << " hits with Z information, parameters z_p: " << flush;
1175 
1176  Int_t nofClosePars = 0;
1177  Int_t highestNPars = 2;
1178  Int_t highestNPPar = -1;
1179  Double_t meanZ_P = 0.;
1180  Double_t sumOfDist = 0.;
1181  Double_t smallestSOD = 1000.;
1182  Double_t bestMeanZ_P = 0.;
1183  if ( fVerbose > 3 )
1184  cout << " got " << trackHitsZ_P.size() << " z_p parameters: " << flush;
1185  for ( size_t ipar = 0 ; ipar < trackHitsZ_P.size() ; ipar++ ) {
1186 
1187  meanZ_P = 0.;
1188  sumOfDist = 0.;
1189  nofClosePars = 0;
1190 
1191  for ( size_t ipar2 = 0 ; ipar2 < trackHitsZ_P.size() ; ipar2++ ) {
1192  if ( TMath::Abs(trackHitsZ_P[ipar]-trackHitsZ_P[ipar2]) > fReasonableZ_PD ) continue;
1193  meanZ_P += trackHitsZ_P[ipar2];
1194  sumOfDist += TMath::Abs(trackHitsZ_P[ipar]-trackHitsZ_P[ipar2]);
1195  nofClosePars++;
1196  }
1197 
1198  if ( nofClosePars >= highestNPars ) {
1199  if ( nofClosePars > highestNPars ) smallestSOD = 1000.;
1200  if ( sumOfDist < smallestSOD ) {
1201  highestNPars = nofClosePars;
1202  highestNPPar = ipar;
1203  smallestSOD = sumOfDist;
1204  bestMeanZ_P = meanZ_P/((Double_t)(nofClosePars));
1205  }
1206  }
1207 
1208  if ( fVerbose > 3 )
1209  cout << trackHitsZ_P[ipar] << "(" << ipar << "/" << nofClosePars << "/" << sumOfDist << ") " << endl;
1210  }
1211  if ( fVerbose > 3 ) {
1212  cout << "best parameter number " << highestNPPar << " with " << highestNPars << " close parameters, sumOfDist = " << smallestSOD << " and meanZ_P = " << bestMeanZ_P << endl;
1213  }
1214 
1215  if ( highestNPars < 3 )
1216  return kFALSE;
1217 
1218  fTracksVector[trackNo].meanZ_P = bestMeanZ_P;
1219 
1220  return kTRUE;
1221 }
1222 // -------------------------------------------------------------------------
1223 
1224 // ----- Private method PrintTracks -------------------------------
1225 void PndBarrelTrackFinder::PrintTracks() {
1226  if ( fTracksVector.size() > 0 ) {
1227  cout << " --------- " << flush;
1228  for ( Int_t idet = 0 ; idet < 4 ; idet++ ) {
1229  if ( fIncludeDet[idet] == kTRUE ) {
1230  cout << "\033[" << 91+idet << "m"
1231  << fDetName[idet].Data()
1232  << "\033[0m " << flush;
1233  if ( idet == 2 )
1234  cout << "\033[33m(skewed)\033[0m " << flush;
1235 
1236  // cout << "\033[" << (skewedSttHit?33:91+fTracksVector[itr].trackHitD[ihit]) << "m"
1237  // << fTracksVector[itr].trackHitN[ihit] << " "
1238  // << "\033[0m" << flush;
1239  }
1240  }
1241  cout << "-------------" << endl;
1242  }
1243 
1244  for ( size_t itr = 0 ; itr < fTracksVector.size() ; itr++ ) {
1245  cout << " | track " << setw(2) << itr << ": " << flush;
1246 
1247  if ( TMath::Abs(fTracksVector[itr].meanR) > 0.1 ) {
1248  for ( size_t ihit = 0 ; ihit < fTracksVector[itr].trackHitD.size() ; ihit++ ) {
1249  Bool_t skewedSttHit = kFALSE;
1250  if ( fTracksVector[itr].trackHitD[ihit] == 2 ) {
1251  Int_t iTube = ((PndSttHit*) fTracksVector[itr].trackHits[ihit])->GetTubeID();
1252  PndSttTube *sttTube = nullptr;
1253  if (fGeoType == 1)
1254  sttTube = (PndSttTube*) fTubeArray->At(iTube);
1255  else if (fGeoType == 2){
1256  sttTube = PndSttTubeMap::Instance()->GetTube(iTube);
1257  }
1258  if ( sttTube->GetWireDirection().Z() < 1. ) {
1259  skewedSttHit = kTRUE;
1260  }
1261  }
1262 
1263  cout << "\033[" << (skewedSttHit?33:91+fTracksVector[itr].trackHitD[ihit]) << "m"
1264  << fTracksVector[itr].trackHitN[ihit] << " "
1265  << "\033[0m" << flush;
1266  }
1267  cout << endl;
1268 
1269  Double_t thisPhi = CalcPhi(fTracksVector[itr].meanX,
1270  fTracksVector[itr].meanY);
1271  Double_t thisRad = fTracksVector[itr].meanR;
1272  Double_t thisZ_P = fTracksVector[itr].meanZ_P;
1273  if ( TMath::IsNaN(thisZ_P) ) thisZ_P = 2.e6;
1274 
1275  Double_t calcPt = thisRad/165.974;
1276  Double_t calcPz = thisZ_P*0.0058997;
1277  Double_t calcP = TMath::Sqrt(calcPt*calcPt+calcPz*calcPz);
1278  Double_t calcThe = TMath::ACos(calcPz/calcP);
1279 
1280  Double_t calcPhi = thisPhi-TMath::Pi()/2.;
1281  //Double_t calcChg = -1.; //[R.K. 01/2017] unused variable?
1282  if ( calcPhi < 0. ) calcPhi += TMath::Pi()*2.;
1283  Double_t firpPhi = CalcPhi(fTracksVector[itr].trackHits[0]->GetX(),fTracksVector[itr].trackHits[0]->GetY());
1284  if ( TMath::Abs(calcPhi-firpPhi) > TMath::Pi()/2. && TMath::Abs(calcPhi-firpPhi) < 3.*TMath::Pi()/2. ) {
1285  calcPhi = thisPhi+TMath::Pi()/2.;
1286  if ( calcPhi > 2.*TMath::Pi() ) calcPhi -= TMath::Pi()*2.;
1287  //calcChg = 1.; //[R.K. 01/2017] unused variable?
1288  }
1289  if ( calcPhi > TMath::Pi() )
1290  calcPhi -= 2.*TMath::Pi();
1291 
1292  TVector3 trackMomentum;
1293  trackMomentum.SetMagThetaPhi(calcP,calcThe,calcPhi);
1294 
1295  cout << " | "
1296  << "\033[96m"
1297  << setw(7) << setprecision(6) << thisPhi << " "
1298  << setw(7) << setprecision(6) << thisRad << " "
1299  << setw(7) << setprecision(6) << thisZ_P << " "
1300  << setw(7) << setprecision(6) << calcPt << " "
1301  << setw(7) << setprecision(6) << calcPz << " "
1302  << setw(7) << setprecision(6) << calcP << " "
1303  << setw(7) << setprecision(6) << calcPhi << " "
1304  << setw(7) << setprecision(6) << calcThe << "\033[0m" << endl;
1305  }
1306 
1307  else {
1308  for ( size_t ihit = 0 ; ihit < fTracksVector[itr].trackHitD.size() ; ihit++ )
1309  cout << setw(1) << fTracksVector[itr].trackHitD[ihit] << "." << fTracksVector[itr].trackHitN[ihit] << " " << flush;
1310  cout << endl;
1311 
1312  for ( size_t ipar = 0 ; ipar < fTracksVector[itr].trackPars.size() ; ipar++ ) {
1313  cout << " | /" << setw(2) << ipar << " "
1314  << setw(6) << setprecision(6) << CalcPhi(fTracksVector[itr].trackPars[ipar].x,fTracksVector[itr].trackPars[ipar].y) << " "
1315  << setw(6) << setprecision(6) << fTracksVector[itr].trackPars[ipar].r << " "
1316  << setw(6) << setprecision(6) << fTracksVector[itr].trackPars[ipar].z_p << " "
1317  << setw(6) << setprecision(6) << fTracksVector[itr].trackPars[ipar].n << endl;
1318  }
1319  }
1320 
1321  cout << " ------------------------------------------------------------" << endl;
1322  }
1323 }
1324 // -------------------------------------------------------------------------
1325 
1326 // ----- Private method PrintTracks -------------------------------
1327 void PndBarrelTrackFinder::DrawTracks() {
1328  for ( Int_t iell = 0 ; iell < 1000 ; iell++ ) {
1329  fEllipse[iell]->SetX1(100);
1330  fEllipse[iell]->SetY1(100);
1331  fEllipse[iell]->SetR1(0);
1332  fEllipse[iell]->SetR2(0);
1333  fEllipse[iell]->Paint();
1334  }
1335 
1336  for ( size_t itr = 0 ; itr < fTracksVector.size() ; itr++ ) {
1337  if ( TMath::Abs(fTracksVector[itr].meanR) > 0.1 ) {
1338  for ( size_t ihit = 0 ; ihit < fTracksVector[itr].trackHitD.size() ; ihit++ ) {
1339  //Bool_t skewedSttHit = kFALSE; //[R.K. 01/2017] unused variable?
1340  if ( fTracksVector[itr].trackHitD[ihit] == 2 ) {
1341  Int_t iTube = ((PndSttHit*) fTracksVector[itr].trackHits[ihit])->GetTubeID();
1342  PndSttTube *sttTube = nullptr;
1343  if (fGeoType == 1)
1344  sttTube = (PndSttTube*) fTubeArray->At(iTube);
1345  else if (fGeoType == 2){
1346  sttTube = PndSttTubeMap::Instance()->GetTube(iTube);
1347  }
1348  if ( sttTube->GetWireDirection().Z() < 1. ) {
1349  //skewedSttHit = kTRUE; //[R.K. 01/2017] unused variable?
1350  }
1351  }
1352  }
1353 
1354  //Double_t thisPhi = CalcPhi(fTracksVector[itr].meanX,
1355  //fTracksVector[itr].meanY); //[R.K. 01/2017] unused variable?
1356  //Double_t thisRad = fTracksVector[itr].meanR; //[R.K. 01/2017] unused variable?
1357 
1358  fEllipse[itr]->SetX1(fTracksVector[itr].meanX);
1359  fEllipse[itr]->SetY1(fTracksVector[itr].meanY);
1360  fEllipse[itr]->SetR1(fTracksVector[itr].meanR);
1361  fEllipse[itr]->SetR2(fTracksVector[itr].meanR);
1362  fEllipse[itr]->SetFillStyle(0);
1363  fEllipse[itr]->SetLineWidth(2);
1364  fEllipse[itr]->SetLineColor(51+2*itr);
1365  fEllipse[itr]->Paint();
1366  }
1367  }
1368 }
1369 // -------------------------------------------------------------------------
1370 
1371 // ----- Private method CleanTracks -------------------------------
1372 Int_t PndBarrelTrackFinder::CleanTracks() {
1373  return -1;
1374 }
1375 // -------------------------------------------------------------------------
1376 
1377 // ----- Private method RemoveShortTracks -------------------------------
1378 void PndBarrelTrackFinder::RemoveShortTracks() {
1379  for ( Int_t itr = fTracksVector.size()-1 ; itr >= 0 ; itr-- ) {
1380  if ( TMath::Abs(fTracksVector[itr].meanZ_P) < 1.e6 && TMath::Abs(fTracksVector[itr].meanR) > 1. ) continue;
1381  for ( size_t ihit = 0 ; ihit < fTracksVector[itr].trackHits.size() ; ihit++ ) {
1382  AddHitToPreviousHits (fTracksVector[itr].trackHits[ihit],
1383  fTracksVector[itr].trackHitD[ihit],
1384  fTracksVector[itr].trackHitN[ihit]);
1385  }
1386  for ( size_t intr = itr ; intr < fTracksVector.size()-1 ; intr++ ) {
1387  fTracksVector[intr] = fTracksVector[intr+1];
1388  }
1389  fTracksVector.pop_back();
1390  }
1391 }
1392 // -------------------------------------------------------------------------
1393 
1394 // ----- Private method RemoveCloneTracks -------------------------------
1395 void PndBarrelTrackFinder::RemoveCloneTracks() {
1396  Bool_t printInfo = kFALSE;
1397  for ( Int_t itr1 = fTracksVector.size()-1 ; itr1 >= 0 ; itr1-- ) {
1398  for ( Int_t itr2 = fTracksVector.size()-1 ; itr2 > itr1 ; itr2-- ) {
1399  fClonesXDiff ->Fill( fTracksVector[itr1].meanX -fTracksVector[itr2].meanX);
1400  fClonesYDiff ->Fill( fTracksVector[itr1].meanY -fTracksVector[itr2].meanY);
1401  fClonesRDiff ->Fill( fTracksVector[itr1].meanR -fTracksVector[itr2].meanR);
1402  fClonesZ_PDiff ->Fill( fTracksVector[itr1].meanZ_P-fTracksVector[itr2].meanZ_P);
1403  fClonesXDiffRel ->Fill((fTracksVector[itr1].meanX -fTracksVector[itr2].meanX )/fTracksVector[itr2].meanR);
1404  fClonesYDiffRel ->Fill((fTracksVector[itr1].meanY -fTracksVector[itr2].meanY )/fTracksVector[itr2].meanR);
1405  fClonesRDiffRel ->Fill((fTracksVector[itr1].meanR -fTracksVector[itr2].meanR )/fTracksVector[itr2].meanR);
1406  fClonesZ_PDiffRel->Fill((fTracksVector[itr1].meanZ_P-fTracksVector[itr2].meanZ_P)/fTracksVector[itr2].meanZ_P);
1407  }
1408  }
1409 
1410  for ( Int_t itr1 = fTracksVector.size()-1 ; itr1 >= 0 ; itr1-- ) {
1411  Double_t tr1Phi = CalcPhi(fTracksVector[itr1].meanX,fTracksVector[itr1].meanY);
1412  Int_t nofTr1Hits = fTracksVector[itr1].trackHitD.size();
1413  for ( Int_t itr2 = fTracksVector.size()-1 ; itr2 > itr1 ; itr2-- ) {
1414  Double_t tr2Phi = CalcPhi(fTracksVector[itr2].meanX,fTracksVector[itr2].meanY);
1415 
1416 // if ( (mcMag-reMag) < 0.3*mcMag &&
1417 // (mcMag-reMag) > -0.3*mcMag &&
1418 // TMath::Abs(mcPhi-rePhi) < 5. &&
1419 // TMath::Abs(mcThe-reThe) < 5. ) {
1420 
1421 // if ( TMath::Abs(fTracksVector[itr1].meanR -fTracksVector[itr2].meanR ) < 0.1*fTracksVector[itr1].meanR &&
1422 // TMath::Abs(fTracksVector[itr1].meanZ_P-fTracksVector[itr2].meanZ_P) < fReasonableZ_PD+0.05*TMath::Abs(fTracksVector[itr1].meanZ_P) &&
1423 // TMath::Abs(tr1Phi-tr2Phi) < 0.1*TMath::RadToDeg() ) {
1424  if ( TMath::Abs(fTracksVector[itr1].meanX -fTracksVector[itr2].meanX ) < 0.2*TMath::Abs(fTracksVector[itr1].meanR) &&
1425  TMath::Abs(fTracksVector[itr1].meanY -fTracksVector[itr2].meanY ) < 0.2*TMath::Abs(fTracksVector[itr1].meanR) &&
1426  TMath::Abs(fTracksVector[itr1].meanR -fTracksVector[itr2].meanR ) < 0.3*fTracksVector[itr1].meanR &&
1427  TMath::Abs(fTracksVector[itr1].meanZ_P-fTracksVector[itr2].meanZ_P) < 0.4*TMath::Abs(fTracksVector[itr1].meanZ_P) &&
1428  TMath::Abs(tr1Phi-tr2Phi) < 0.2*TMath::RadToDeg() ) {
1429  if ( fVerbose > 0 || printInfo ) {
1430  cout << "seems that tracks " << itr1 << " (" << fTracksVector[itr1].meanX << "," << fTracksVector[itr1].meanY << "," << fTracksVector[itr1].meanR << "," << fTracksVector[itr1].meanZ_P << " )" << endl;
1431  cout << " and " << itr2 << " (" << fTracksVector[itr2].meanX << "," << fTracksVector[itr2].meanY << "," << fTracksVector[itr2].meanR << "," << fTracksVector[itr2].meanZ_P << " ) ARE THE SAME!" << endl;
1432  }
1433 
1434  Int_t nofTr2Hits = fTracksVector[itr2].trackHitD.size();
1435  for ( Int_t ith2 = 0 ; ith2 < nofTr2Hits ; ith2++ ) {
1436  Bool_t hitExists = kFALSE;
1437  for ( Int_t ith1 = 0 ; ith1 < nofTr1Hits ; ith1++ ) {
1438  if ( fTracksVector[itr1].trackHitD[ith1] == fTracksVector[itr2].trackHitD[ith2] &&
1439  fTracksVector[itr1].trackHitN[ith1] == fTracksVector[itr2].trackHitN[ith2] ) {
1440  hitExists = kTRUE;
1441  break;
1442  }
1443  }
1444  if ( hitExists ) continue;
1445  fTracksVector[itr1].trackHits.push_back(fTracksVector[itr2].trackHits[ith2]);
1446  fTracksVector[itr1].trackHitD.push_back(fTracksVector[itr2].trackHitD[ith2]);
1447  fTracksVector[itr1].trackHitN.push_back(fTracksVector[itr2].trackHitN[ith2]);
1448  // cout << endl;
1449  }
1450  for ( size_t intr = itr2 ; intr < fTracksVector.size()-1 ; intr++ ) {
1451  fTracksVector[intr] = fTracksVector[intr+1];
1452  }
1453  fTracksVector.pop_back();
1454  }
1455  }
1456  }
1457 }
1458 // -------------------------------------------------------------------------
1459 
1460 // ----- Private method WriteTracks -------------------------------
1461 Int_t PndBarrelTrackFinder::WriteTracks() {
1462  Int_t nofCreatedTracks = 0;
1463 
1464  for ( size_t itr = 0 ; itr < fTracksVector.size() ; itr++ ) {
1465  Double_t thisPhi = CalcPhi(fTracksVector[itr].meanX,
1466  fTracksVector[itr].meanY);
1467  Double_t thisRad = fTracksVector[itr].meanR;
1468  Double_t thisZ_P = fTracksVector[itr].meanZ_P;
1469  if ( TMath::IsNaN(thisZ_P) ) thisZ_P = 2.e6;
1470 
1471  // for ( Int_t ihit = 0 ; ihit < fTracksVector[itr].trackHitD.size() ; ihit++ )
1472  // cout << fTracksVector[itr].trackHitD[ihit] << "." << fTracksVector[itr].trackHitN[ihit] << " " << flush;
1473  // cout << endl;
1474  Double_t calcPt = thisRad/165.974;
1475  Double_t calcPz = thisZ_P*0.0058997;
1476  Double_t calcP = TMath::Sqrt(calcPt*calcPt+calcPz*calcPz);
1477  Double_t calcThe = TMath::ACos(calcPz/calcP);
1478 
1479  Double_t calcPhi = thisPhi-TMath::Pi()/2.;
1480  Int_t calcChg = -1;
1481  if ( calcPhi < 0. ) calcPhi += TMath::Pi()*2.;
1482  Double_t firpPhi = CalcPhi(fTracksVector[itr].trackHits[0]->GetX(),fTracksVector[itr].trackHits[0]->GetY());
1483  if ( TMath::Abs(calcPhi-firpPhi) > TMath::Pi()/2. && TMath::Abs(calcPhi-firpPhi) < 3.*TMath::Pi()/2. ) {
1484  calcPhi = thisPhi+TMath::Pi()/2.;
1485  if ( calcPhi > 2.*TMath::Pi() ) calcPhi -= TMath::Pi()*2.;
1486  calcChg = 1;
1487  }
1488  if ( calcPhi > TMath::Pi() )
1489  calcPhi -= 2.*TMath::Pi();
1490 
1491  // cout << "F " << thisPhi << " " << thisRad << " " << thisZ_P << " " << calcPt << " " << calcPhi << " " << calcPz << " " << calcP << " " << calcPhi << " " << calcThe << endl;//flush;
1492 
1493  // cout << " MTP: " << calcP << " " << calcPhi << " " << calcThe << " and charge = " << calcChg << flush;
1494 
1495  TVector3 trackPosition(0.,0.,0.);
1496  TVector3 trackMomentum;
1497  trackMomentum.SetMagThetaPhi(calcP,calcThe,calcPhi);
1498 
1499  PndTrackCand* trackCand = new((*fBarrelTrackCandArray)[nofCreatedTracks]) PndTrackCand();
1500 
1501  std::vector<Int_t> trackHitsPerDet(4,0);
1502  for ( size_t ihit = 0 ; ihit < fTracksVector[itr].trackHitD.size() ; ihit++ ) {
1503  trackHitsPerDet[fTracksVector[itr].trackHitD[ihit]] += 1;
1504  Double_t tempPos = TMath::Sqrt(fTracksVector[itr].trackHits[ihit]->GetX()*fTracksVector[itr].trackHits[ihit]->GetX()+
1505  fTracksVector[itr].trackHits[ihit]->GetY()*fTracksVector[itr].trackHits[ihit]->GetY()+
1506  fTracksVector[itr].trackHits[ihit]->GetZ()*fTracksVector[itr].trackHits[ihit]->GetZ());
1507 // cout << ">>> " << fDetName[fTracksVector[itr].trackHitD[ihit]] << " hit will be written as "
1508 // << fDetType[fTracksVector[itr].trackHitD[ihit]] << "." << fTracksVector[itr].trackHitN[ihit] << endl;
1509 
1510  trackCand->AddHit(FairRootManager::Instance()->GetBranchId(fHitArrayName[fTracksVector[itr].trackHitD[ihit]].Data()),
1511  fTracksVector[itr].trackHitN[ihit],
1512  tempPos);
1513  }
1514  trackCand->Sort();
1515 
1516 // cout << "track hits in :" << flush;
1517 // for ( Int_t itch = 0 ; itch < trackCand->GetNHits() ; itch++ ) {
1518 // cout << trackCand->GetSortedHit(itch).GetDetId() << " " << flush;
1519 // }
1520 // cout << endl;
1521 
1522  FairTrackParP firstPar(trackPosition,trackMomentum,
1523  TVector3(0.5, 0.5, 0.5),
1524  0.1*trackMomentum,
1525  calcChg,
1526  trackPosition,
1527  TVector3(1.,0.,0.),
1528  TVector3(0.,1.,0.));
1529  TVector3 tempVect1(trackHitsPerDet[0],trackHitsPerDet[1],trackHitsPerDet[2]);
1530  TVector3 tempVect2(trackHitsPerDet[3],-1.,-1.);
1531 
1532  // FairTrackParP* lastPar = new FairTrackParP(trackPosition,trackMomentum,
1533  FairTrackParP lastPar (tempVect1,tempVect2,
1534  TVector3(0.5, 0.5, 0.5),
1535  0.1*trackMomentum,
1536  calcChg,
1537  trackPosition,
1538  TVector3(1.,0.,0.),
1539  TVector3(0.,1.,0.));
1540 
1541  new((*fBarrelTrackArray)[nofCreatedTracks]) PndTrack(firstPar,lastPar,*trackCand);
1542 
1543  nofCreatedTracks++;
1544 
1545  // cout << endl;
1546  }
1547 
1548  return nofCreatedTracks;
1549 }
1550 // -------------------------------------------------------------------------
1551 
1552 // ----- Private method RemoveHitFromPreviousHits -------------------------------
1553 void PndBarrelTrackFinder::RemoveHitFromPreviousHits(Int_t hitNo) {
1554  // cout << endl << "There are " << fTracksVector.size() << " segments in the event" << endl;
1555  for ( size_t irh = hitNo ; irh < fHitVector.size()-1 ; irh++ ) {
1556  fHitVector[irh] = fHitVector[irh+1];
1557  fHitVectDI[irh] = fHitVectDI[irh+1];
1558  fHitVectHN[irh] = fHitVectHN[irh+1];
1559  }
1560  fHitVector.pop_back();
1561  fHitVectDI.pop_back();
1562  fHitVectHN.pop_back();
1563 }
1564 // -------------------------------------------------------------------------
1565 
1566 // ----- Private method AddHitToPreviousHits -------------------------------
1567 void PndBarrelTrackFinder::AddHitToPreviousHits (FairHit* thisHit, Int_t detId, Int_t hitNo) {
1568  if ( fVerbose > 3 )
1569  cout << "!!! ADDING HIT " << detId << "." << hitNo << endl;
1570  fHitVector.push_back(thisHit);
1571  fHitVectDI.push_back(detId);
1572  fHitVectHN.push_back(hitNo);
1573 }
1574 // -------------------------------------------------------------------------
1575 
1576 
1577 // ----- Private method FindCircPar -------------------------------
1578 // the method finds a circle cl tangent to the three given circles c1, c2, c3 (circle given by x,y,r)
1579 // parameter cno (from 0 to 7) is a bit mask on whether a given circle is inside or outside circle cl
1580 Bool_t PndBarrelTrackFinder::FindCircPar(Double_t* c1, Double_t* c2, Double_t* c3, Int_t cno, Double_t* cl) {
1581  if ( cno < 0 || cno > 7 ) return kFALSE;
1582 // cout << "c1 : " << c1[0] << "," << c1[1] << " - " << c1[2] << endl;
1583 // cout << "c2 : " << c2[0] << "," << c2[1] << " - " << c2[2] << endl;
1584 // cout << "c3 : " << c3[0] << "," << c3[1] << " - " << c3[2] << endl;
1585  Double_t s1 = -1+2*(cno%2); cno/=2;
1586  Double_t s2 = -1+2*(cno%2); cno/=2;
1587  Double_t s3 = -1+2*(cno%2);
1588 
1589  Double_t a2 = 2.*(c1[0]-c2[0]);
1590  Double_t b2 = 2.*(c1[1]-c2[1]);
1591  Double_t e2 = 2.*(s1*c1[2]-s2*c2[2]);
1592  Double_t d2 = (c1[0]*c1[0]+c1[1]*c1[1]-c1[2]*c1[2]) - (c2[0]*c2[0]+c2[1]*c2[1]-c2[2]*c2[2]);
1593 
1594  Double_t a3 = 2.*(c1[0]-c3[0]);
1595  Double_t b3 = 2.*(c1[1]-c3[1]);
1596  Double_t e3 = 2.*(s1*c1[2]-s3*c3[2]);
1597  Double_t d3 = (c1[0]*c1[0]+c1[1]*c1[1]-c1[2]*c1[2]) - (c3[0]*c3[0]+c3[1]*c3[1]-c3[2]*c3[2]);
1598  // cout << "a2 " << a2 << " a3 " << a3 << " b2 " << b2 << " b3 " << b3 << " /// a2*b3-a3*b2 = " << a2*b3-a3*b2 << endl;
1599  if ( TMath::Abs(a2*b3-a3*b2) < 0.001 ) return kFALSE;
1600 
1601  Double_t ax = (d2*b3-d3*b2)/(a2*b3-a3*b2);
1602  Double_t bx = (e2*b3-e3*b2)/(a2*b3-a3*b2);
1603  Double_t ay = (a2*d3-a3*d2)/(a2*b3-a3*b2);
1604  Double_t by = (a2*e3-a3*e2)/(a2*b3-a3*b2);
1605 
1606  Double_t sa = bx*bx+by*by-1.;
1607  Double_t sb = 2.*bx*(ax-c1[0])+2.*by*(ay-c1[1])+2.*s1*c1[2];
1608  Double_t sc = (ax-c1[0])*(ax-c1[0])+(ay-c1[1])*(ay-c1[1])-c1[2]*c1[2];
1609  if ( sa == 0 ) return kFALSE;
1610  Double_t dt = sb*sb-4.*sa*sc;
1611  if ( dt < 0 ) return kFALSE;
1612 
1613  cl[2] = (-sb-TMath::Sqrt(dt))/(2.*sa);
1614  cl[0] = ax+bx*cl[2];
1615  cl[1] = ay+by*cl[2];
1616  cl[2] = TMath::Abs(cl[2]);
1617  if ( TMath::IsNaN(cl[2]) ) {
1618  return kFALSE;
1619  }
1620  return kTRUE;
1621 }
1622 // -------------------------------------------------------------------------
1623 
1624 // ----- Private method FindCircDist -------------------------------
1625  Double_t PndBarrelTrackFinder::FindCircDist(Double_t* c1, Double_t* c2) {
1626  Double_t circDist = TMath::Sqrt((c1[0]-c2[0])*(c1[0]-c2[0])+
1627  (c1[1]-c2[1])*(c1[1]-c2[1]));
1628  if ( circDist >= c1[2]+c2[2] ) { return circDist-c1[2]-c2[2];}
1629  if ( circDist <= TMath::Abs(c1[2]-c2[2]) ) { return TMath::Abs(c1[2]-c2[2])-circDist;}
1630  if ( c1[2] > c2[2] ) { return c2[2]-TMath::Abs(circDist-c1[2]);}
1631  return c1[2]-TMath::Abs(circDist-c2[2]);
1632 }
1633 // -------------------------------------------------------------------------
1634 
1635 // ----- Private method CalcPhi ----------------------------------------
1636 Int_t PndBarrelTrackFinder::FindInterestingRegions(Double_t* circ, Double_t* tube, Double_t* reg) {
1637  Double_t ax = tube[4]/tube[6];
1638  Double_t ay = tube[5]/tube[6];
1639  Double_t a = ax*ax+ay*ay;
1640  Double_t bx = (tube[0]-ax*tube[2]-circ[0]);
1641  Double_t by = (tube[1]-ay*tube[2]-circ[1]);
1642  Double_t b = 2.*ax*bx+2.*ay*by;
1643  Double_t c = bx*bx+by*by;
1644 
1645  Double_t delta[2];
1646  delta[0] = b*b-4.*a*(c-(circ[2]+tube[7])*(circ[2]+tube[7]));
1647  delta[1] = b*b-4.*a*(c-(circ[2]-tube[7])*(circ[2]-tube[7]));
1648 
1649  Int_t nofReg = 0;
1650 
1651  for ( Int_t idel = 0 ; idel < 2 ; idel++ ) {
1652  if ( delta[idel] < 0 ) continue;
1653  Double_t z[2];
1654  z[0] = (-b-TMath::Sqrt(delta[idel]))/2./a;
1655  z[1] = (-b+TMath::Sqrt(delta[idel]))/2./a;
1656  // cout << " tube crosses circle in " << z[0] << " or " << z[1] << endl;
1657 
1658  for ( Int_t iz = 0 ; iz < 2 ; iz++ ) {
1659  Double_t z_p = 2.e6;
1660  if ( TMath::Abs(z[iz]-tube[2]) > tube[3] ) continue;
1661  // cout << "checking z = " << z[iz] << endl;
1662  Double_t x = tube[0] + ax * ( z[iz] - tube[2] );
1663  Double_t y = tube[1] + ay * ( z[iz] - tube[2] );
1664  Double_t xd = x - circ[0];
1665  Double_t yd = y - circ[1];
1666  Double_t rd = TMath::Sqrt(xd*xd+yd*yd);
1667  Double_t pointPhi = TMath::ACos(xd/rd);
1668  if ( yd < 0 )
1669  pointPhi = 2.*TMath::Pi() - pointPhi;
1670  pointPhi += TMath::Pi();
1671  if ( pointPhi > 2.*TMath::Pi() ) pointPhi -= 2.*TMath::Pi();
1672  Double_t phiDiffr = (pointPhi-circ[3]);
1673  if ( phiDiffr > TMath::Pi() ) phiDiffr -= 2.*TMath::Pi();
1674  if ( phiDiffr < -TMath::Pi() ) phiDiffr += 2.*TMath::Pi();
1675  z_p = z[iz]/TMath::Abs(phiDiffr);
1676  // cout << "the z_p = " << z_p << " cause dz = " << z[iz] << " and dphi = " << phiDiffr << endl;
1677  // cout << "got " << z_p << " while from other calc. it is: " << CalcZ_P(circ,x,y,z[iz]) << endl;
1678  reg[nofReg] = z_p;
1679  nofReg++;
1680  }
1681  }
1682  return nofReg;
1683 }
1684 // -------------------------------------------------------------------------
1685 
1686 // ----- Private method CalcPhi ----------------------------------------
1687 Double_t PndBarrelTrackFinder::CalcPhi(Double_t x, Double_t y) {
1688  Double_t r = TMath::Sqrt(x*x+y*y);
1689  Double_t p = TMath::ACos(x/r);
1690  if ( y < 0 )
1691  p = 2.*TMath::Pi() - p;
1692  return p;
1693 }
1694 // -------------------------------------------------------------------------
1695 
1696 // ----- Private method CalcZ_P ----------------------------------------
1697 Double_t PndBarrelTrackFinder::CalcZ_P(Double_t* circ, Double_t hx, Double_t hy, Double_t hz) {
1698  Double_t tmpx = hx-circ[0];
1699  Double_t tmpy = hy-circ[1];
1700  Double_t rotA = TMath::Pi()-CalcPhi(circ[0],circ[1]);
1701  Double_t tmpx2 = tmpx*TMath::Cos(rotA)-tmpy*TMath::Sin(rotA);
1702  Double_t tmpy2 = tmpx*TMath::Sin(rotA)+tmpy*TMath::Cos(rotA);
1703  Double_t phiT = CalcPhi(tmpx2,tmpy2);
1704  if ( phiT > TMath::Pi() ) phiT = -(2.*TMath::Pi()-phiT);
1705  return hz/TMath::Abs(phiT);
1706  //Double_t phiT = CalcPhi(hx-circ[0],hy-circ[1]);
1707  // return phiT/hz;
1708 }
1709 // -------------------------------------------------------------------------
1710 
1711 // ----- Private method SetParContainers -------------------------------
1712 void PndBarrelTrackFinder::SetParContainers() {
1713 
1714  // Get run and runtime database
1715  FairRunAna* run = FairRunAna::Instance();
1716  if ( ! run ) Fatal("SetParContainers", "No analysis run");
1717 
1718  FairRuntimeDb* db = run->GetRuntimeDb();
1719  if ( ! db ) Fatal("SetParContainers", "No runtime database");
1720 
1721  fSttParameters = (PndGeoSttPar*) db->getContainer("PndGeoSttPar");
1722 
1723  cout << "-I- PndBarrelTrackFinder::SetParContainers()." << endl;
1724 }
1725 // -------------------------------------------------------------------------
1726 
1727 
1728 
1729 // ----- Private method Init -------------------------------------------
1730 InitStatus PndBarrelTrackFinder::Init() {
1731  // fTrackSegments.clear();
1732 
1733  // Get and check FairRootManager
1734  FairRootManager* ioman = FairRootManager::Instance();
1735  if( !ioman ) {
1736  cout << "-E- "<< GetName() <<"::Init: "
1737  << "RootManager not instantised!" << endl;
1738  return kERROR;
1739  }
1740 
1741  // Get the pointer to the singleton FairRunAna object
1742  FairRunAna* ana = FairRunAna::Instance();
1743  if(NULL == ana) {
1744  cout << "-E- "<< GetName() <<"::Init :"
1745  <<" no FairRunAna object!" << endl;
1746  return kERROR;
1747  }
1748  // Get the pointer to run-time data base
1749  FairRuntimeDb* rtdb = ana->GetRuntimeDb();
1750  if(NULL == rtdb) {
1751  cout << "-E- "<< GetName() <<"::Init :"
1752  <<" no runtime database!" << endl;
1753  return kERROR;
1754  }
1755 
1756  fDetName[0] = "MVD Pixel";
1757  fDetName[1] = "MVD Strip";
1758  fDetName[2] = "STT";
1759  fDetName[3] = "GEM";
1760 
1761  fHitArrayName[0] = "MVDHitsPixel";
1762  fHitArrayName[1] = "MVDHitsStrip";
1763  fHitArrayName[2] = "STTHit";
1764  fHitArrayName[3] = "GEMHit";
1765 
1766  fDetType[0] = kMVDHitsPixel;// FIXME: Use FairRootManager::Instance()->GetBranchId()
1767  fDetType[1] = kMVDHitsStrip;// FIXME: Use FairRootManager::Instance()->GetBranchId()
1768  fDetType[2] = kSttHit;
1769  fDetType[3] = kGemHit;
1770 
1771 
1772  for ( Int_t idet = 0 ; idet < 4 ; idet++ ) {
1773  if ( fIncludeDet[idet] == kTRUE ) {
1774  fHitArray[idet] = (TClonesArray*) ioman->GetObject(fHitArrayName[idet].Data());
1775 
1776  if ( !fHitArray[idet] ) {
1777  cout << "-E- " << GetName() << "::Init: No " << fHitArrayName[idet].Data() << " array!" << endl;
1778  return kERROR;
1779  }
1780  std::cout << "-I- " << GetName() << ": " << fHitArrayName[idet].Data() << " array found" << std::endl;
1781  }
1782  }
1783 
1784  // Create and register PndTrack and PndTrackCand arrays
1785  fBarrelTrackArray = new TClonesArray("PndTrack",100);
1786  ioman->Register("BarrelTrack", "Barrel Track", fBarrelTrackArray, GetPersistency());
1787 
1788  fBarrelTrackCandArray = new TClonesArray("PndTrackCand",100);
1789  ioman->Register("BarrelTrackCand", "Barrel TrackCand", fBarrelTrackCandArray, GetPersistency());
1790 
1791  std::cout << "-I- " << GetName() << ": Initialization successfull" << std::endl;
1792  std::cout << "-I- " << GetName() << ": Looking for tracks in " << flush;
1793  for ( Int_t idet = 0 ; idet < 4 ; idet++ ) {
1794  if ( fIncludeDet[idet] == kFALSE ) continue;
1795  cout << fDetName[idet].Data() << ", ";
1796  }
1797  cout << "\b\b. " << endl;
1798  cout << "================================================================================" << endl;
1799 
1800  fGeoType = fSttParameters->GetGeometryType();
1801  if ( fGeoType == 1){
1802  PndSttMapCreator *mapper = new PndSttMapCreator(fSttParameters);
1803  fTubeArray = mapper->FillTubeArray();
1804  }
1805 // cout << "****************************************************" << endl;
1806 // cout << "** " << fTubeArray->GetEntriesFast() << " TUBES ***********************************" << endl;
1807 // cout << "****************************************************" << endl;
1808 
1809  // DRAW
1810  fCanvas = NULL;
1811  if ( fDrawDetails )
1812  fCanvas = new TCanvas("AnimatedCanvas","AnimatedCanvas",10,10,800,800);
1813 
1814  if ( fCanvas ) {
1815  for ( Int_t iell = 0 ; iell < 1000 ; iell++ ) {
1816  fEllipse[iell] = new TEllipse(100,100,0);
1817  }
1818  fStageDesc = new TLatex(0,0,"Init");
1819  }
1820 
1821  fClonesXDiff = new TH1F("fClonesXDiff" ,"fClonesXDiff" ,20000,-100.,100.);
1822  fClonesYDiff = new TH1F("fClonesYDiff" ,"fClonesYDiff" ,20000,-100.,100.);
1823  fClonesRDiff = new TH1F("fClonesRDiff" ,"fClonesRDiff" ,20000,-100.,100.);
1824  fClonesZ_PDiff = new TH1F("fClonesZ_PDiff" ,"fClonesZ_PDiff" ,20000,-100.,100.);
1825  fClonesXDiffRel = new TH1F("fClonesXDiffRel" ,"fClonesXDiffRel" ,20000,-100.,100.);
1826  fClonesYDiffRel = new TH1F("fClonesYDiffRel" ,"fClonesYDiffRel" ,20000,-100.,100.);
1827  fClonesRDiffRel = new TH1F("fClonesRDiffRel" ,"fClonesRDiffRel" ,20000,-100.,100.);
1828  fClonesZ_PDiffRel = new TH1F("fClonesZ_PDiffRel","fClonesZ_PDiffRel",20000,-100.,100.);
1829 
1830  return kSUCCESS;
1831 }
1832 // -------------------------------------------------------------------------
1833 
1834 
1835 
1836 // ----- Private method ReInit -----------------------------------------
1837 InitStatus PndBarrelTrackFinder::ReInit() {
1838 
1839  return kSUCCESS;
1840 
1841 }
1842 // -------------------------------------------------------------------------
1843 
1844 
1845 
1846 // ----- Private method Reset ------------------------------------------
1847 void PndBarrelTrackFinder::Reset() {
1848 }
1849 // -------------------------------------------------------------------------
1850 
1851 // ----- Public method Finish ------------------------------------------
1852 void PndBarrelTrackFinder::Finish() {
1853  fBarrelTrackArray ->Clear();
1854  fBarrelTrackCandArray->Clear();
1855 
1856  if ( fDrawHistos ) {
1857  new TCanvas("xdiff","xdiff",10,10,1000,800);
1858  fClonesXDiff ->Draw();
1859  new TCanvas("ydiff","ydiff",10,10,1000,800);
1860  fClonesYDiff ->Draw();
1861  new TCanvas("rdiff","rdiff",10,10,1000,800);
1862  fClonesRDiff ->Draw();
1863  new TCanvas("zdiff","`diff",10,10,1000,800);
1864  fClonesZ_PDiff ->Draw();
1865  new TCanvas("xdiffR","xdiffR",10,10,1000,800);
1866  fClonesXDiffRel ->Draw();
1867  new TCanvas("ydiffR","ydiffR",10,10,1000,800);
1868  fClonesYDiffRel ->Draw();
1869  new TCanvas("rdiffR","rdiffR",10,10,1000,800);
1870  fClonesRDiffRel ->Draw();
1871  new TCanvas("zdiffR","zdiffR",10,10,1000,800);
1872  fClonesZ_PDiffRel->Draw();
1873  }
1874 
1875  cout << "-------------------- " << fName.Data() << " : Summary -----------------------" << endl;
1876  cout << " Events: " << setw(10) << fTNofEvents << endl;
1877  cout << " Tracks: " << setw(10) << fTNofTracks << " ( " << (Double_t)fTNofTracks/((Double_t)fTNofEvents) << " per event )" << endl;
1878  cout << "--------------------------------------------------------------------------------" << endl;
1879 }
1880 // -------------------------------------------------------------------------
1881 
1882 
1883 
1884 ClassImp(PndBarrelTrackFinder)
1885 
fVerbose
int fVerbose
Definition: poormantracks.C:24
z_p
Double_t z_p
Definition: anal_point_fwendcap.C:21
run
Int_t run
Definition: autocutx.C:47
PndBarrelTrackFinder::FindInterestingRegions
Int_t FindInterestingRegions(Double_t *circ, Double_t *tube, Double_t *reg)
Definition: PndBarrelTrackFinder.cxx:1636
r
double r
Definition: RiemannTest.C:14
PndBarrelTrackFinder::HitBelongsToTrack
Bool_t HitBelongsToTrack(Int_t detId, Int_t hitNo, Int_t trackNo)
Definition: PndBarrelTrackFinder.cxx:510
b
TTree * b
Definition: GammaSpectraAnalysis.C:24
PndSttTube::GetHalfLength
Double_t GetHalfLength()
Definition: PndSttTube.cxx:99
PndBarrelTrackFinder::fHitVectHN
std::vector< Int_t > fHitVectHN
Definition: PndBarrelTrackFinder.h:117
map
PndTransMap * map
Definition: sim_emc_apd.C:99
CAMath::Sqrt
static T Sqrt(const T &x)
Definition: PndCAMath.h:37
PndBarrelTrackFinder::RemoveHitFromPreviousHits
void RemoveHitFromPreviousHits(Int_t hitNo)
Definition: PndBarrelTrackFinder.cxx:1553
TrackBasis::trackHits
std::vector< FairHit * > trackHits
Definition: PndBarrelTrackFinder.h:55
CAMath::Sin
static T Sin(const T &x)
Definition: PndCAMath.h:42
n
int n
Definition: HypStatDecay.cxx:21
PndBarrelTrackFinder::fTNofEvents
Int_t fTNofEvents
event counter
Definition: PndBarrelTrackFinder.h:130
PndBarrelTrackFinder::ExtractMeanRPhiFromTrack
Bool_t ExtractMeanRPhiFromTrack(Int_t trackNo)
Definition: PndBarrelTrackFinder.cxx:1053
c2
c2
Definition: plot_dirc.C:39
PndBarrelTrackFinder::fBarrelTrackCandArray
TClonesArray * fBarrelTrackCandArray
Definition: PndBarrelTrackFinder.h:127
PndPersistencyTask::SetPersistency
void SetPersistency(Bool_t val=kTRUE)
Definition: PndPersistencyTask.h:31
PndBarrelTrackFinder::fHitArray
TClonesArray * fHitArray[4]
Definition: PndBarrelTrackFinder.h:120
PndSttTubeMap::GetTube
PndSttTube * GetTube(int tubeId)
Definition: PndSttTubeMap.h:22
PndBarrelTrackFinder::fHitVectDI
std::vector< Int_t > fHitVectDI
Definition: PndBarrelTrackFinder.h:116
PndGeoSttPar::GetGeometryType
Int_t GetGeometryType()
Definition: PndGeoSttPar.h:28
CAMath::Cos
static T Cos(const T &x)
Definition: PndCAMath.h:43
PndBarrelTrackFinder::CalcZ_P
Double_t CalcZ_P(Double_t *circ, Double_t hx, Double_t hy, Double_t hz)
Definition: PndBarrelTrackFinder.cxx:1697
p
Double_t p
Definition: anasim.C:58
TrackBasis::trackHitD
std::vector< Int_t > trackHitD
Definition: PndBarrelTrackFinder.h:56
At
cout<< "blue = Monte Carlo "<< endl;cout<< "red = Helix Hit "<< endl;cout<< "green = Center Of Tubes "<< endl;for(Int_t k=0;k< track->GetEntriesFast();k++){PndSttTrack *stttrack=(PndSttTrack *) track-> At(k)
Definition: checkhelixhit.C:56
PndSttTubeMap::Instance
static PndSttTubeMap * Instance()
Definition: PndSttTubeMap.cxx:24
PndBarrelTrackFinder::fHitDetNo
std::vector< Int_t > fHitDetNo
Definition: PndBarrelTrackFinder.h:113
CAMath::Abs
static T Abs(const T &x)
Definition: PndCAMath.h:39
b2
Double_t b2
Definition: createRootGeoFileBwEnd_FwEndCapGeo.C:22
TrackBasis::trackHitN
std::vector< Int_t > trackHitN
Definition: PndBarrelTrackFinder.h:57
a
Int_t a
Definition: anaLmdDigi.C:126
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
PndSttTube::GetPosition
TVector3 GetPosition()
Definition: PndSttTube.cxx:87
rtdb
FairRuntimeDb * rtdb
Definition: hit_dirc.C:66
PndSttMapCreator::FillTubeArray
TClonesArray * FillTubeArray()
Definition: PndSttMapCreator.cxx:128
PndBarrelTrackFinder::fEllipse
TEllipse * fEllipse[1000]
Definition: PndBarrelTrackFinder.h:98
PndBarrelTrackFinder::fHitDetId
std::vector< Int_t > fHitDetId
Definition: PndBarrelTrackFinder.h:112
z
Double_t z
Definition: createRootFscGeometryFile.C:146
PndBarrelTrackFinder::MatchHitWithZInfoTT
Bool_t MatchHitWithZInfoTT(FairHit *thisHit, Int_t detId, Int_t hitNo, Int_t trackNo)
Definition: PndBarrelTrackFinder.cxx:748
c1
c1
Definition: plot_dirc.C:35
PndBarrelTrackFinder::MatchSkewedSttHitTT
Bool_t MatchSkewedSttHitTT(FairHit *thisHit, Int_t detId, Int_t hitNo, Int_t trackNo)
Definition: PndBarrelTrackFinder.cxx:521
name
TString name
Definition: createRootGeoFile.C:82
c
GeV c
Definition: ReadMesDecCorrRes_PiT.C:28
c3
c3
Definition: plot_dirc.C:50
PndBarrelTrackFinder::fHitVector
std::vector< FairHit * > fHitVector
Definition: PndBarrelTrackFinder.h:115
tmpx
Double_t tmpx
Definition: anaLmdDigi.C:65
x
Double_t x
Definition: createRootFscGeometryFile.C:159
PndBarrelTrackFinder::UseMvdSttGem
void UseMvdSttGem(const Bool_t useMvd, const Bool_t useStt, const Bool_t useGem)
Definition: PndBarrelTrackFinder.cxx:200
PndBarrelTrackFinder::FindCircDist
Double_t FindCircDist(Double_t *c1, Double_t *c2)
Definition: PndBarrelTrackFinder.cxx:1625
PndBarrelTrackFinder::MatchHitToTrack
Bool_t MatchHitToTrack(FairHit *thisHit, Int_t detId, Int_t hitNo, Int_t trackNo)
Definition: PndBarrelTrackFinder.cxx:503
PndBarrelTrackFinder::Exec
virtual void Exec(Option_t *opt)
Definition: PndBarrelTrackFinder.cxx:209
PndBarrelTrackFinder::fTracksVector
std::vector< TrackBasis > fTracksVector
Definition: PndBarrelTrackFinder.h:110
ClassImp
ClassImp(PndAnaContFact)
iVerbose
Int_t iVerbose
Definition: Lars/runMvdDigi.C:7
PndBarrelTrackFinder::MatchHitToHit
Bool_t MatchHitToHit(FairHit *thisHit, Int_t detId, Int_t hitNo, Int_t prevHNo)
Definition: PndBarrelTrackFinder.cxx:913
y
Double_t y
Definition: createRootFscGeometryFile.C:159
tmpy
Double_t tmpy
Definition: anaLmdDigi.C:65
PndBarrelTrackFinder::fBarrelTrackArray
TClonesArray * fBarrelTrackArray
Output array of PndTrackCands.
Definition: PndBarrelTrackFinder.h:126
PndBarrelTrackFinder::AddHitToPreviousHits
void AddHitToPreviousHits(FairHit *thisHit, Int_t detId, Int_t hitNo)
Definition: PndBarrelTrackFinder.cxx:1567
TrackBasis::trackPars
std::vector< TrackParameter > trackPars
Definition: PndBarrelTrackFinder.h:54
Pi
Double_t Pi
Definition: PndGeoHypGeCrystal.cxx:10
PndBarrelTrackFinder::CalcPhi
Double_t CalcPhi(Double_t x, Double_t y)
Definition: PndBarrelTrackFinder.cxx:1687
PndBarrelTrackFinder::FindCircPar
Bool_t FindCircPar(Double_t *c1, Double_t *c2, Double_t *c3, Int_t cno, Double_t *cl)
Definition: PndBarrelTrackFinder.cxx:1580
PndBarrelTrackFinder::MatchParallelSttHitTT
Bool_t MatchParallelSttHitTT(FairHit *thisHit, Int_t detId, Int_t hitNo, Int_t trackNo)
Definition: PndBarrelTrackFinder.cxx:615
PndBarrelTrackFinder::ExtractMeanZ_PFromTrack
Bool_t ExtractMeanZ_PFromTrack(Int_t trackNo)
Definition: PndBarrelTrackFinder.cxx:1113
PndBarrelTrackFinder::ReInit
virtual InitStatus ReInit()
Definition: PndBarrelTrackFinder.cxx:1837
PndSttTube::GetWireDirection
TVector3 GetWireDirection()
Definition: PndSttTube.cxx:107
PndBarrelTrackFinder::Init
virtual InitStatus Init()
Definition: PndBarrelTrackFinder.cxx:1730
Generated on Sun Aug 18 2019 03:01:17 for FairRoot/PandaRoot by   doxygen 1.8.8