CbmRoot
CbmL1CATrdTrackFinderSA.cxx
Go to the documentation of this file.
1 // -----------------------------------------------------------------------
2 // ----- CbmL1CATrdTrackFinderSA -----
3 // ----- Created 2/12/2006 by A. Bubak & M. Krauze -----
4 // -----------------------------------------------------------------------
5 
6 #include "CbmL1CATrdTrackFinderSA.h"
7 
8 #include "CbmL1Def.h"
9 
10 #include "CbmL1.h"
11 #include "CbmL1TrdTracklet.h"
12 #include "CbmL1TrdTracklet4.h"
13 
14 #include "CbmGeoTrdPar.h"
15 #include "CbmKF.h"
16 #include "CbmKFTrack.h"
17 #include "CbmKFTrdHit.h"
18 #include "CbmMCTrack.h"
19 #include "CbmStsTrack.h"
20 #include "CbmTrdHit.h"
21 #include "CbmTrdPoint.h"
22 #include "CbmTrdTrack.h"
23 #include "CbmTrdTrackFitterKF.h"
24 #include "CbmVertex.h"
25 #include "FairBaseParSet.h"
26 #include "FairDetector.h"
27 #include "FairGeoNode.h"
28 #include "FairMCPoint.h"
29 #include "FairRootManager.h"
30 #include "FairRunAna.h"
31 #include "FairRuntimeDb.h"
32 
33 #include "TClonesArray.h"
34 #include "TH1F.h"
35 #include "TH2F.h"
36 #include "TLinearFitter.h"
37 
38 #include <cmath>
39 #include <iostream>
40 #include <list>
41 #include <map>
42 #include <vector>
43 
44 using std::cout;
45 using std::endl;
46 using std::fabs;
47 using std::map;
48 using std::pair;
49 using std::set;
50 using std::vector;
51 
52 //#include "L1CATrdDraw.h"
53 
54 ClassImp(CbmL1CATrdTrackFinderSA);
55 
56 
57 CbmL1CATrdTrackFinderSA* CbmL1CATrdTrackFinderSA::fInstance = 0;
58 
59 // ----------------------- Default constructor ---------------------------
60 CbmL1CATrdTrackFinderSA::CbmL1CATrdTrackFinderSA()
61  : geoLayer()
62  , fNTrdHits(0)
63  , fNoTrdStations(0)
64  , fNoTrdPerStation(0)
65  , planeHits()
66  , fivTrdHitArr()
67  , fArrayTrdHit(new TClonesArray("CbmTrdHit"))
68  , fArrayTrdTrack(new TClonesArray("CbmTrdTrack"))
69  , TrdPar(0)
70  , fEvents(0)
71  , fNofEvents(0)
72  , fMCTrackArray(0)
73  , fMCPointArray(0)
74  , trdTrackFitterKF(new CbmTrdTrackFitterKF(1, 11))
75  , fVerbose(1)
76  ,
77 
78  fvTempArray()
79  , fvFoundTracks()
80  , tempTrack()
81  ,
82 
83  iStation1()
84  , iStation2()
85  , fImapSt1()
86  , fImapSt2()
87  , itTrackletsLeft()
88  , itTrackletsRight()
89  ,
90 
91  iHitMap1(0)
92  , iHitMap2(0)
93  , iHitMapY1(0)
94  , iHitMapY2(0)
95  , fTrd13_Z(0)
96  , fTrd14_Z(0)
97  , fTrd21_Z(0)
98  , fTrd24_Z(0)
99  , fTrd31_Z(0)
100  , fKfTrack(0)
101  , createSegments()
102  , findNeighbour()
103  , createSPs()
104  , createSPs_SL()
105  , sortSPs()
106  , doFind()
107  , sortHits()
108  , tagTracks()
109  , createTracks()
110  , selectTracks()
111  , delTime()
112  , secondLoopTime()
113  , refittingKF()
114  , thirdLoopTime()
115  , totCreateSegments(0)
116  , totFindNeighbour(0)
117  , totCreateSPs(0)
118  , totCreateSPs_SL(0)
119  , totSortSPs(0)
120  , totDoFind(0)
121  , totSortHits(0)
122  , totTagTracks(0)
123  , totCreateTracks(0)
124  , totDelTime(0)
125  , totSecondLoopTime(0)
126  , totThirdLoopTime(0)
127  , totSelectTracks(0)
128  , totRefittingKF(0)
129  , fh_chi2hit(0)
130  , fh_chi2hit_plane(0)
131  , fDistLongX(0)
132  , fDistLongY(0)
133  , fDistShortX(0)
134  , fDistShortY(0)
135  , fDistLongBX(0)
136  , fDistLongBY(0)
137  , fDistShortBX(0)
138  , fDistShortBY(0)
139  , fDistY12(0)
140  , fMomDistLongPrimaryX(0)
141  , fMomDistLongPrimaryY(0)
142  , fMomDistLongExtraX(0)
143  , fMomDistLongExtraY(0)
144  , fMomDistExtrapolPrimaryX(0)
145  , fMomDistExtrapolPrimaryY(0)
146  , fMomDistExtrapolExtraX(0)
147  , fMomDistExtrapolExtraY(0)
148  , fMomDistShortPrimaryX(0)
149  , fMomDistShortPrimaryY(0)
150  , fMomDistShortExtraX(0)
151  , fMomDistShortExtraY(0)
152  , fDistY(0)
153  , fDistX(0)
154  , fPlane1Ydens(0)
155  , fPlane5Ydens(0)
156  , fPlane9Ydens(0)
157  , fSPlength(0)
158  , fSPlengthMC(0)
159  , fYat0(0)
160  , fYat0MC(0)
161  , fNoEvTime(0)
162  , fh_SP_xDiff_MC(0)
163  , fh_SP_yDiff_MC(0)
164  , fh_SP_xDiff_nMC(0)
165  , fh_SP_yDiff_nMC(0)
166  , fUsedHitsPerPlane(0)
167  , fUnUsedHitsPerPlane(0)
168  , fRUsedHits()
169  , fRUnUsedHits()
170  , fTotHits() {
171  if (!fInstance) fInstance = this;
172 }
173 // -----------------------------------------------------------------------
174 
175 // --------------------------- Destructor --------------------------------
176 CbmL1CATrdTrackFinderSA::~CbmL1CATrdTrackFinderSA() {
177  delete fArrayTrdHit;
178  delete fArrayTrdTrack;
179  delete trdTrackFitterKF;
180  // delete fKfTrack;
181  if (fInstance == this) fInstance = 0;
182 }
183 // -----------------------------------------------------------------------
184 
185 
186 // ------------------------- Initialisation ------------------------------
187 void CbmL1CATrdTrackFinderSA::Init() {
188  // Create histogramms
189  CreateHistogramms();
190  // Activate data branches
191  DataBranches();
192 
193  FairRootManager* ioman = FairRootManager::Instance();
194  if (!ioman) {
195  cout << "-E- CbmL1CATrdTrackFinderSA::Init: "
196  << "RootManager not instantised!" << endl;
197  return;
198  }
199 
200  // Get MCTrack array
201  fMCTrackArray = L1_DYNAMIC_CAST<TClonesArray*>(ioman->GetObject("MCTrack"));
202  if (!fMCTrackArray) {
203  cout << "-E- CbmL1CATrdTrackFinderSA::Init: No MCTrack array!" << endl;
204  return;
205  }
206 
207  // Get MCPoint array
208  fMCPointArray = L1_DYNAMIC_CAST<TClonesArray*>(ioman->GetObject("TrdPoint"));
209  if (!fMCPointArray) {
210  cout << "-E- CbmL1CATrdTrackFinderSA::Init: No MCPoint array!" << endl;
211  return;
212  }
213 
214  // Determine the TRD layout
215  TrdLayout();
216 
217  // Init KF Track fitter
218  trdTrackFitterKF->Init();
219 }
220 // -----------------------------------------------------------------------
221 
222 
223 // -------------------- Algorithm implementation -------------------------
224 Int_t CbmL1CATrdTrackFinderSA::DoFind(TClonesArray* hitArray,
225  TClonesArray* trackArray) {
226 
227 #ifdef DRAW
228  InitL1Draw();
229 #endif
230 
231  fArrayTrdHit = hitArray;
232  fArrayTrdTrack = trackArray;
233 
234  fNofEvents++;
235 
236  doFind.Start();
237 
238  // Initialize control counters
239  // Int_t nNoTrdPoint = 0;
240  Int_t nNoTrdHit = 0;
241  Int_t trdPlaneID = 0;
242  Int_t ptIndex = 0;
243  Int_t noHitsPerTracklet = 4;
244 
245  // Create pointers to TRD hit and TrdPoint
246  CbmTrdHit* pHit = NULL;
247  CbmTrdPoint* pMCpt = NULL;
248  // CbmMCTrack* pMCtr = NULL;
249 
250  //--- Function: CreateSpacePoint --------------------------------------
251  Double_t // b
252  sigmaA_FL = 2, //3 3 3 2.5 2.5
253  sigmaB_FL = 2, //3 3 3 2.5 2.5
254  sigmaA_SL = 3, //4 4 4 3 3
255  sigmaB_SL = 3, //4 4 4 3 3
256  sigmaA_TL = 4, //4 4 4 4 4
257  sigmaB_TL = 4; //4 4 4 4 4
258  //--- Function: CreateSegments; Data from MC, look inside function ----
259  Double_t dY_FL = 0.3, //0.5 0.5 0.5 0.4 0.4
260  dX_FL = 0.3, //0.5 0.5 0.5 0.4 0.4
261  dY_SL = 0.5, //0.7 0.6 0.8 0.6 0.6
262  dX_SL = 0.5, //0.8 0.8 0.8 0.6 0.6
263  dY_TL = 0.7, //0.7 0.6 0.8 0.7 0.8
264  dX_TL = 0.7; //0.8 0.8 0.8 0.7 0.8
265 
266  //--- Function: FindNeighbour ------------------------------------------
267  //distance from Y-propagated point, around which we look for neighbours
268  Double_t distPropLongY_FL = 2.5, //3 3 3 3 2
269  distPropLongX_FL = 2.5, //3 3 3 3 2
270  distPropLongY_SL = 3, //4 4 4 3 3
271  distPropLongX_SL = 3, //8 6 8 3 3
272  distPropLongY_TL = 4, //4 4 4 3 4
273  distPropLongX_TL = 4; //8 6 8 3 4
274 
275  //----------------------------------------------------------------------
276  Bool_t competition = true;
277  Bool_t removeUsedHits = true;
278  Bool_t secondLoop = true;
279  Bool_t thirdLoop = true;
280  Int_t segment = 0; //0 = long tracks, 1 = segments only
281  //Int_t nrLoop = 1;
282 
283  //---- Simulate detector inefficiency -----------------------------------
284  Int_t accept = 100;
285  set<Int_t> globalSetRejectedHits;
286  globalSetRejectedHits.clear();
287  //----------------------------------------------------------------------
288 
289 
290  sortHits.Start();
291  //--- hit sorting procedure ------------------------------------------------
292  Int_t nTrdHit = hitArray->GetEntriesFast();
293  for (Int_t iHit = 0; iHit < nTrdHit; iHit++) {
294  // Loop over hits. Get corresponding MCPoint and MCTrack index
295  pHit = L1_DYNAMIC_CAST<CbmTrdHit*>(hitArray->At(iHit));
296 
297  if (!Rejection(accept)) {
298  // Simulate detector inefficiency
299  globalSetRejectedHits.insert(iHit);
300  continue;
301  }
302  if (!pHit) {
303  cout << "-E- CbmL1CATrdTrackFinderSA::DoFind: Empty slot "
304  << "in HitArray at position " << iHit << endl;
305  nNoTrdHit++;
306  continue;
307  }
308 
309  ptIndex = pHit->GetRefId();
310  //if (ptIndex < 0) continue; // fake or background hit
311  pMCpt = L1_DYNAMIC_CAST<CbmTrdPoint*>(fMCPointArray->At(ptIndex));
312  //if (!pMCpt) {
313  // nNoTrdPoint++;
314  // continue;
315  //}
316  // pMCtr = L1_DYNAMIC_CAST<CbmMCTrack*>(fMCTrackArray->At(pMCpt->GetTrackID()));
317  //if ( ! pMCtr ) continue;
318 
319  trdPlaneID = pHit->GetPlaneId();
320  Int_t trdPlaneIDN = trdPlaneID - 1;
321 
322  planeHits.mcTrackID = pMCpt->GetTrackID();
323  planeHits.hitIndex = iHit;
324  planeHits.X = pHit->GetX();
325  planeHits.Y = pHit->GetY();
326  planeHits.Z = pHit->GetZ();
327  planeHits.DX = pHit->GetDx();
328  planeHits.DY = pHit->GetDy();
329  planeHits.planeID = trdPlaneID;
330  fvTrdHitArr[trdPlaneIDN].push_back(planeHits);
331 
332  // cout << "geoLayer: " <<geoLayer.Y[trdPlaneID] <<", "<< geoLayer.Z[trdPlaneID] << endl;
333  //ptIndex = pHit->GetRefIndex();
334  //pMCpt = (CbmTrdPoint*) (fMCPointArray->At(ptIndex));
335 
336  //planeHits.hitIndex = ptIndex;
337  //planeHits.X = pMCpt->GetX();
338  //planeHits.Y = pMCpt->GetY();
339  //planeHits.Z = pMCpt->GetZ();
340  //planeHits.DX = 0;
341  //planeHits.DY = 0;;
342  //fvTrdPointArr[trdPlaneIDN].push_back(planeHits);
343  }
344 
345  for (Int_t i = 0; i < 12; i++) {
346  sort(fvTrdHitArr[i].begin(), fvTrdHitArr[i].end(), CompareY);
347  }
348 
349  /*
350  //-----------------------------------------------------------
351  vector <LayerWithHits>::iterator ivT;
352  LayerWithHits pl;
353  for(ivT = fvTrdHitArr[0].begin();
354  ivT != fvTrdHitArr[0].end();
355  ivT++){
356  pl = (*ivT);
357  cout <<"(X:Y) "<<pl.X<<"\t"<< pl.Y << endl;
358  fPlane1Ydens->Fill(pl.Y);
359  }
360 
361  for(ivT = fvTrdHitArr[5].begin();
362  ivT != fvTrdHitArr[5].end();
363  ivT++){
364  pl = (*ivT);
365  // cout << pl.Y << endl;
366  fPlane5Ydens->Fill(pl.Y);
367  }
368  for(ivT = fvTrdHitArr[9].begin();
369  ivT != fvTrdHitArr[9].end();
370  ivT++){
371  pl = (*ivT);
372  // cout << pl.Y << endl;
373  fPlane9Ydens->Fill(pl.Y);
374  }
375  //-----------------------------------------------------------
376  */
377  sortHits.Stop();
378 
379  if (fVerbose > 2) {
380  cout << "-I- CbmL1CATrdTrackFinderSA::DoFind " << endl;
381  for (Int_t i = 0; i < 12; i++) {
382  cout << " Size of vector " << i << ": " << fvTrdHitArr[i].size() << endl;
383  }
384  }
385 
386  //create vectors that will hold the tracklet (4-hit) objects
387  vector<CbmL1TrdTracklet4*> clTracklets[3];
388  vector<CbmL1TrdTracklet4*> clTrackletsNew[3];
389 
390  //create vectors that will hold the mini-tracklet (2-hit) objects - Space Points
391  vector<CbmL1TrdTracklet*> clSpacePoints[6];
392 
393  if (noHitsPerTracklet == 4) { //noHitsPerTracklet == 4
394  //------ creation of mini tracklet 1-2 -------------------------------------
395  cout << endl
396  << "-I- CbmL1CATrdTrackFinderSA::DoFind " << endl
397  << "--------------------------------------------------" << endl
398  << " ### Start creation of Space Points" << endl
399  << "--------------------------------------------------" << endl;
400 
401  createSPs.Start();
402  for (Int_t i = 0, j = 0; i < 6; i++, j = j + 2) {
403  CreateSpacePoints(fvTrdHitArr[j],
404  fvTrdHitArr[j + 1],
405  clSpacePoints[i],
406  sigmaA_FL,
407  sigmaB_FL);
408  }
409  createSPs.Stop();
410 
411  //------- Sorting spacepoints -------------------------------------------------
412  sortSPs.Start();
413  for (Int_t i = 0; i < 6; i++) {
414  sort(clSpacePoints[i].begin(),
415  clSpacePoints[i].end(),
416  CbmL1TrdTracklet::Compare1);
417  }
418  sortSPs.Stop();
419 
420  // --------------- Creation tracklet from 4 space points ----------------
421  cout << endl
422  << "-I- CbmL1CATrdTrackFinderSA::DoFind " << endl
423  << "--------------------------------------------------" << endl
424  << " ### Start creation tracklets" << endl
425  << "--------------------------------------------------" << endl;
426  createSegments.Start();
427  //joining SP together to create 4-hit tracklets
428  for (Int_t i = 0, j = 0; i < 3; i++, j = j + 2) {
429  CreateSegments(
430  clSpacePoints[j], clSpacePoints[j + 1], clTracklets[i], dX_FL, dY_FL);
431  }
432  createSegments.Stop();
433 
434  //--- TESTING -------------------------------
435  /* How many segments have 4 hits belonging to the same track (has the same trackID)*/
436  /*
437  vector<CbmL1TrdTracklet4*>::iterator itclTrackletsRight;
438  vector<CbmL1TrdTracklet4*>::iterator iTrackletT;
439  Int_t iSumSegment = 0;
440  CbmTrdHit *trdHit;
441  CbmTrdPoint* trdPoint;
442 
443  for(int i = 0; i < 3; i++){
444  Double_t clSegmentSize = clTracklets[i].size();
445  cout <<"Tracklet size: "<< clSegmentSize << endl;
446  iSumSegment = 0;
447 
448  for(iTrackletT = clTracklets[i].begin();
449  iTrackletT != clTracklets[i].end();
450  iTrackletT++){
451 
452  Int_t k = 0;
453  Int_t Ind0 = (*iTrackletT)->GetInd(0);
454  trdHit = L1_DYNAMIC_CAST<CbmTrdHit*>( hitArray->At(Ind0);
455  trdPoint = (CbmTrdPoint*) fMCPointArray->At(trdHit->GetRefIndex());
456  Double_t trackID1 = trdPoint->GetTrackID();
457 
458  for(int iw = 1; iw < 4; iw++){
459  Int_t Ind1 = (*iTrackletT)->GetInd(iw);
460  trdHit = L1_DYNAMIC_CAST<CbmTrdHit*>( hitArray->At(Ind1);
461  trdPoint = (CbmTrdPoint*) fMCPointArray->At(trdHit->GetRefIndex());
462  Double_t trackID2 = trdPoint->GetTrackID();
463 
464  if(trackID1 == trackID2) k++;
465  }
466  if(k == 3){
467  iSumSegment++;
468  }
469  }
470 
471  Double_t sumSegmentWith4 = (iSumSegment*100)/clSegmentSize;
472  //printf("(size %f, %i, %f)\n", clSegmentSize, iSumSegment, sumSegmentWith4);
473  }
474  //---------------------------------------------
475  */
476 
477  //---------- BEGIN: Find friend ------------------------------------------------
478  findNeighbour.Start();
479 
480  cout << "--- Finding neighbour 14-58" << endl;
481  FindNeighbour(
482  clTracklets[0], clTracklets[1], distPropLongX_FL, distPropLongY_FL);
483 
484  cout << "--- Finding neighbour 58-912" << endl;
485  FindNeighbour(
486  clTracklets[1], clTracklets[2], distPropLongX_FL, distPropLongY_FL);
487 
488  findNeighbour.Stop();
489 
490  } //end of noHitsPerTracklet == 4
491 
492  //###################################################################################################
493 
494 
495  cout << endl
496  << "-I- CbmL1CATrdTrackFinderSA::DoFind " << endl
497  << "--------------------------------------------------" << endl
498  << " ### (FL) in Event No. " << fNofEvents << " ###" << endl
499  << "--------------------------------------------------" << endl;
500 
501  if (fVerbose > 2) {
502  cout << "size of segment vector 14 = " << clTracklets[0].size() << endl
503  << "size of segment vector 58 = " << clTracklets[1].size() << endl
504  << "size of segment vector 912 = " << clTracklets[2].size() << endl;
505  }
506 
507  if (fVerbose > 1)
508  cout << "-I- CbmL1CATrdTrackFinderSA::DoFind " << endl
509  << "Tracklet finding phase completed." << endl
510  << "Now constructing tracks from tracklets..." << endl;
511 
512 
513  tagTracks.Start();
514 
515  if (noHitsPerTracklet == 4) { //noHitsPerTracklet == 4
516  //BEGIN: Tagging procedure ------------------------------------------------------
517  cout << endl
518  << "-I- CbmL1CATrdTrackFinderSA::DoFind " << endl
519  << "--------------------------------------------------" << endl
520  << " ### Starting tagging procedure" << endl
521  << "--------------------------------------------------" << endl;
522 
523  TagSegments(clTracklets[1], clTracklets[2], 1);
524  cout << "--- Tagging 58-921 done." << endl;
525 
526  TagSegments(clTracklets[0], clTracklets[1], 0);
527  cout << "--- Tagging 14-58 done." << endl;
528  } //end of noHitsPerTracklet == 4
529 
530  tagTracks.Stop();
531  //-------------------------------------------------------------------------
532 
533 
534  //--------------------------------------------------------------------------------------
535  // Counting number of segments which have a value of 4, 3, 2, ... etc.
536  CbmL1TrdTracklet4* clSegRight; //segment nearer to the target
537 
538  map<Int_t, Int_t> segValues14;
539  map<Int_t, Int_t> segValues58;
540  map<Int_t, Int_t> segValues912;
541  Int_t segValue;
542 
543  vector<CbmL1TrdTracklet4*>::iterator itclTrackletsRight;
544  //for tracklets in the 1st station
545  for (itclTrackletsRight = clTracklets[0].begin();
546  itclTrackletsRight != clTracklets[0].end();
547  itclTrackletsRight++) {
548 
549  clSegRight = *itclTrackletsRight;
550  segValue = clSegRight->GetVal();
551  segValues14[segValue]++;
552  if (segValue == 3) clTrackletsNew[0].push_back(clSegRight);
553  }
554 
555  //for tracklets in the 2nd station
556  for (itclTrackletsRight = clTracklets[1].begin();
557  itclTrackletsRight != clTracklets[1].end();
558  itclTrackletsRight++) {
559 
560  clSegRight = *itclTrackletsRight;
561  segValue = clSegRight->GetVal();
562  segValues58[segValue]++;
563  if (segValue == 2) clTrackletsNew[1].push_back(clSegRight);
564  }
565 
566  //for tracklets in the 3rd station
567  for (itclTrackletsRight = clTracklets[2].begin();
568  itclTrackletsRight != clTracklets[2].end();
569  itclTrackletsRight++) {
570 
571  clSegRight = *itclTrackletsRight;
572  segValue = clSegRight->GetVal();
573  segValues912[segValue]++;
574  if (segValue == 1) clTrackletsNew[2].push_back(clSegRight);
575  }
576 
577  //--- For test only ----------------------------------------------------
578  if (fVerbose > 0) {
579  map<Int_t, Int_t>::iterator mIt;
580 
581  cout << "--- Map no. 14 has: " << endl;
582  for (mIt = segValues14.begin(); mIt != segValues14.end(); mIt++) {
583  cout << mIt->first << " segment has a count of " << mIt->second << endl;
584  }
585 
586  cout << "--- Map no. 58 has: " << endl;
587  for (mIt = segValues58.begin(); mIt != segValues58.end(); mIt++) {
588  cout << mIt->first << " segment has a count of " << mIt->second << endl;
589  }
590 
591  cout << "--- Map no. 912 has: " << endl;
592  for (mIt = segValues912.begin(); mIt != segValues912.end(); mIt++) {
593  cout << mIt->first << " segment has a count of " << mIt->second << endl;
594  }
595 
596  cout << endl;
597  }
598  //-------------------------------------------------------------------------
599 
600  //vector<CbmTrdTrack*> vTrdTrackCand;
601  //vector<CbmTrdTrack*> vTempTrdTrackCand;
602  //vector<CbmTrdTrack*>::iterator ivTempTrdTrackCand;
603  //vector<CbmTrdTrack*> trackCand1;
604 
605  cout << "--------------------------------------------------" << endl
606  << " ### Starting creating the tracks " << endl
607  << "--------------------------------------------------" << endl;
608 
609  createTracks.Start();
610 
611  set<Int_t> globalSetUsedHits;
612  globalSetUsedHits.clear();
613 
614  if (segment == 0) { //segment == 0 -> we create long tracks
615 
616  if (fVerbose > 2) {
617  cout << "Outer area: " << endl
618  << "--- size of fArrayTrdTrack = "
619  << fArrayTrdTrack->GetEntriesFast() << endl
620  << "--- size of globalSetUsedHits = " << globalSetUsedHits.size()
621  << endl;
622  }
623 
624  //create long tracks from 3 sorts of tracklets - station 1,2,3
625  CreateTracks(clTracklets[0],
626  clTracklets[1],
627  clTracklets[2],
628  globalSetUsedHits,
629  removeUsedHits,
630  competition,
631  1);
632 
633  createTracks.Stop();
634 
635  //clearing and removing objects that are not needed anymore
636  for (Int_t i = 0; i < 3; i++) {
637  DeleteTracklets(clTracklets[i]);
638  clTracklets[i].clear();
639  }
640  //----------------------------------------------------
641 
642  for (Int_t i = 0; i < 6; i++) {
643  DeleteTracklets(clSpacePoints[i]);
644  clSpacePoints[i].clear();
645  }
646 
647  //--- [Second Loop] ------------------------------------------------------
648  secondLoopTime.Start();
649 
650  cout << endl
651  << "-I- CbmL1CATrdTrackFinderSA::DoFind (Second loop)" << endl
652  << "--------------------------------------------------" << endl
653  << " ### (SL) in Event No. " << fNofEvents << " ###" << endl
654  << "--------------------------------------------------" << endl;
655 
656  if (secondLoop) {
657 
658  for (Int_t i = 0; i < 12; i++) {
659  fvTrdHitArr[i].clear();
660  }
661 
662  set<Int_t>::iterator iglobalSetUsedHits;
663  set<Int_t>::iterator iglobalSetRejectedHits;
664 
665  for (Int_t iHit = 0; iHit < nTrdHit; iHit++) {
666  pHit = L1_DYNAMIC_CAST<CbmTrdHit*>(hitArray->At(iHit));
667 
668  // if(!Rejection(accept)) continue;
669 
670  if (!pHit) {
671  cout << "-E- CbmL1CATrdTrackFinderSA::DoFind: Empty slot "
672  << "in HitArray at position " << iHit << endl;
673  nNoTrdHit++;
674  continue;
675  }
676 
677  //---- Simulate detector inefficiency -----------------------------------
678  iglobalSetRejectedHits = globalSetRejectedHits.find(iHit);
679  if (iglobalSetRejectedHits != globalSetRejectedHits.end()) continue;
680  //----------------------------------------------------------------------
681 
682  iglobalSetUsedHits = globalSetUsedHits.find(iHit);
683  if (iglobalSetUsedHits != globalSetUsedHits.end()) continue;
684 
685  trdPlaneID = pHit->GetPlaneId();
686  Int_t trdPlaneIDN = trdPlaneID - 1;
687 
688  planeHits.mcTrackID = pMCpt->GetTrackID();
689  planeHits.hitIndex = iHit;
690  planeHits.X = pHit->GetX();
691  planeHits.Y = pHit->GetY();
692  planeHits.Z = pHit->GetZ();
693  planeHits.DX = pHit->GetDx();
694  planeHits.DY = pHit->GetDy();
695  fvTrdHitArr[trdPlaneIDN].push_back(planeHits);
696  }
697 
698  for (Int_t i = 0; i < 12; i++) {
699  sort(fvTrdHitArr[i].begin(), fvTrdHitArr[i].end(), CompareY);
700  }
701 
702  //for(Int_t i = 0; i < 12; i++){
703  // cout <<" Size of vector "<< i <<": "<< fvTrdHitArr[i].size() << endl;
704  //}
705  //--------------------------------------------------------------------------------------
706 
707  createSPs_SL.Start();
708 
709  //performing second step - make all previous tasks but with less restrictive conditions
710  cout << "[Second Loop]::Creating space points" << endl;
711 
712  for (Int_t i = 0, j = 0; i < 6; i++, j = j + 2) {
713  CreateSpacePoints(fvTrdHitArr[j],
714  fvTrdHitArr[j + 1],
715  clSpacePoints[i],
716  sigmaA_SL,
717  sigmaB_SL);
718  }
719  createSPs_SL.Stop();
720 
721  for (Int_t i = 0; i < 6; i++) {
722  sort(clSpacePoints[i].begin(),
723  clSpacePoints[i].end(),
724  CbmL1TrdTracklet::Compare1);
725  }
726 
727  cout << "[Second Loop]::Creating tracklets" << endl;
728  for (Int_t i = 0, j = 0; i < 3; i++, j = j + 2) {
729  CreateSegments(
730  clSpacePoints[j], clSpacePoints[j + 1], clTracklets[i], dX_SL, dY_SL);
731  }
732 
733  if (fVerbose > 2) {
734  cout << "--- size of segment vector 14 = " << clTracklets[0].size()
735  << endl
736  << "--- size of segment vector 58 = " << clTracklets[1].size()
737  << endl
738  << "--- size of segment vector 912 = " << clTracklets[2].size()
739  << endl;
740  }
741 
742  cout << "[Second Loop]::Finding friends 14-58" << endl;
743  FindNeighbour(
744  clTracklets[0], clTracklets[1], distPropLongX_SL, distPropLongY_SL);
745 
746  cout << "[Second Loop]::Finding friends 58-912" << endl;
747  FindNeighbour(
748  clTracklets[1], clTracklets[2], distPropLongX_SL, distPropLongY_SL);
749 
750  cout << "[Second Loop]::Tagging segments 58-912" << endl;
751  TagSegments(clTracklets[1], clTracklets[2], 1);
752 
753  //--------------------------------------------------------------------------------------------
754  cout << "[Second Loop]::Tagging segments 14-58" << endl;
755  TagSegments(clTracklets[0], clTracklets[1], 0);
756 
757  //--------------------------------------------------------------------------------------------
758  cout << "[Second Loop]::Creating tracks" << endl;
759  CreateTracks(clTracklets[0],
760  clTracklets[1],
761  clTracklets[2],
762  globalSetUsedHits,
763  removeUsedHits,
764  competition,
765  2);
766  //competition -> selects the best track candidate from each branch
767  }
768 
769  for (Int_t i = 0; i < 3; i++) {
770  DeleteTracklets(clTracklets[i]);
771  clTracklets[i].clear();
772  }
773  for (Int_t i = 0; i < 6; i++) {
774  DeleteTracklets(clSpacePoints[i]);
775  clSpacePoints[i].clear();
776  }
777 
778  secondLoopTime.Stop();
779 
780  thirdLoopTime.Start();
781 
782  cout << endl
783  << "-I- CbmL1CATrdTrackFinderSA::DoFind (Third loop)" << endl
784  << "--------------------------------------------------" << endl
785  << " ### (TL) in Event No. " << fNofEvents << " ###" << endl
786  << "--------------------------------------------------" << endl;
787 
788  if (thirdLoop) {
789 
790  for (Int_t i = 0; i < 12; i++) {
791  fvTrdHitArr[i].clear();
792  }
793 
794  set<Int_t>::iterator iglobalSetUsedHits;
795  set<Int_t>::iterator iglobalSetRejectedHits;
796 
797  for (Int_t iHit = 0; iHit < nTrdHit; iHit++) {
798  pHit = L1_DYNAMIC_CAST<CbmTrdHit*>(hitArray->At(iHit));
799 
800  // if(!Rejection(accept)) continue;
801 
802  if (!pHit) {
803  cout << "-E- CbmL1CATrdTrackFinderSA::DoFind: Empty slot "
804  << "in HitArray at position " << iHit << endl;
805  nNoTrdHit++;
806  continue;
807  }
808 
809  //---- Simulate detector inefficiency -----------------------------------
810  iglobalSetRejectedHits = globalSetRejectedHits.find(iHit);
811  if (iglobalSetRejectedHits != globalSetRejectedHits.end()) continue;
812  //----------------------------------------------------------------------
813 
814  iglobalSetUsedHits = globalSetUsedHits.find(iHit);
815  if (iglobalSetUsedHits != globalSetUsedHits.end()) continue;
816 
817  trdPlaneID = pHit->GetPlaneId();
818  Int_t trdPlaneIDN = trdPlaneID - 1;
819 
820  planeHits.mcTrackID = pMCpt->GetTrackID();
821  planeHits.hitIndex = iHit;
822  planeHits.X = pHit->GetX();
823  planeHits.Y = pHit->GetY();
824  planeHits.Z = pHit->GetZ();
825  planeHits.DX = pHit->GetDx();
826  planeHits.DY = pHit->GetDy();
827  fvTrdHitArr[trdPlaneIDN].push_back(planeHits);
828  }
829 
830  for (Int_t i = 0; i < 12; i++) {
831  sort(fvTrdHitArr[i].begin(), fvTrdHitArr[i].end(), CompareY);
832  }
833 
834  for (Int_t i = 0, j = 0; i < 6; i++, j = j + 2) {
835  CreateSpacePoints(fvTrdHitArr[j],
836  fvTrdHitArr[j + 1],
837  clSpacePoints[i],
838  sigmaA_TL,
839  sigmaB_TL);
840  }
841 
842  for (Int_t i = 0; i < 6; i++) {
843  sort(clSpacePoints[i].begin(),
844  clSpacePoints[i].end(),
845  CbmL1TrdTracklet::Compare1);
846  }
847 
848  for (Int_t i = 0, j = 0; i < 3; i++, j = j + 2) {
849  CreateSegments(
850  clSpacePoints[j], clSpacePoints[j + 1], clTracklets[i], dX_TL, dY_TL);
851  }
852 
853  FindNeighbour(
854  clTracklets[0], clTracklets[1], distPropLongX_TL, distPropLongY_TL);
855  FindNeighbour(
856  clTracklets[1], clTracklets[2], distPropLongX_TL, distPropLongY_TL);
857  TagSegments(clTracklets[1], clTracklets[2], 1);
858 
859  TagSegments(clTracklets[0], clTracklets[1], 0);
860 
861  CreateTracks(clTracklets[0],
862  clTracklets[1],
863  clTracklets[2],
864  globalSetUsedHits,
865  removeUsedHits,
866  competition,
867  3);
868  }
869 
870  thirdLoopTime.Stop();
871 
872  } //end segment == 0
873 
874 
875  if (segment == 1) { //segment == 1 -> we create only 4-hit, 1-segment tracks!
876  CreateAndManageSegments(clTracklets[0], clTracklets[1], clTracklets[2]);
877  }
878 
879 
880  //------ Used and Unused hists in the TRD's planes ------------------------
881  //------ Histogramming
882  set<Int_t>::iterator iglSetUsedHits;
883  nTrdHit = hitArray->GetEntriesFast();
884 
885  for (Int_t iHit = 0; iHit < nTrdHit; iHit++) {
886  pHit = L1_DYNAMIC_CAST<CbmTrdHit*>(hitArray->At(iHit));
887 
888  trdPlaneID = pHit->GetPlaneId();
889  fTotHits[trdPlaneID]++;
890 
891  iglSetUsedHits = globalSetUsedHits.find(iHit);
892  if (iglSetUsedHits != globalSetUsedHits.end()) {
893  fRUsedHits[trdPlaneID]++;
894  //fUsedHitsPerPlane->Fill(trdPlaneID-1);
895  } else {
896  fRUnUsedHits[trdPlaneID]++;
897  //fUnUsedHitsPerPlane->Fill(trdPlaneID-1);
898  }
899  }
900  //---------------------------------------------------------------------------
901 
902 
903  cout << endl
904  << "-I- CbmL1CATrdTrackFinderSA::DoFind " << endl
905  << "--------------------------------------------------" << endl
906  << " ### Summary" << endl
907  << "--------------------------------------------------" << endl
908  << "--- Number of found tracks: " << fArrayTrdTrack->GetEntriesFast()
909  << endl
910  << endl;
911 
912  delTime.Start();
913 
914  cout << "\n-I- CbmL1CATrdTrackFinderSA::DoFind " << endl
915  << "--------------------------------------------------" << endl
916  << " ### Delete:Clear" << endl
917  << "--------------------------------------------------" << endl;
918 
919  for (Int_t i = 0; i < 12; i++) {
920  fvTrdHitArr[i].clear();
921  }
922 
923  //----------------------------------------------------
924  for (Int_t i = 0; i < 3; i++) {
925  DeleteTracklets(clTracklets[i]);
926  clTracklets[i].clear();
927  }
928 
929  //----------------------------------------------------
930  for (Int_t i = 0; i < 6; i++) {
931  DeleteTracklets(clSpacePoints[i]);
932  clSpacePoints[i].clear();
933  }
934  //----------------------------------------------------
935 
936  //segValues14.clear();
937  //segValues58.clear();
938  //segValues912.clear();
939 
940  delTime.Stop();
941  doFind.Stop();
942 
943 
944  //-------------------------------------------------------
945 #ifdef DRAW
946 
947  fNTrdHits = hitArray->GetEntriesFast();
948 
949  for (Int_t iHit = 0; iHit < fNTrdHits; iHit++) {
950  pHit = L1_DYNAMIC_CAST<CbmTrdHit*>( hitArray->At(iHit);
951 
952  ptIndex = pHit->GetRefIndex();
953  pMCpt = (CbmTrdPoint*) (fMCPointArray->At(ptIndex));
954 // pMCtr = (CbmMCTrack*) fMCTrackArray->At(pMCpt->GetTrackID());
955 
956  trdPlaneID = pHit->GetPlaneID();
957  Int_t trdPlaneIDN = trdPlaneID-1;
958 
959  planeHits.mcTrackID = pMCpt->GetTrackID();
960  planeHits.hitIndex = iHit;
961  planeHits.X = pHit->GetX();
962  planeHits.Y = pHit->GetY();
963  planeHits.Z = pHit->GetZ();
964  planeHits.DX = pHit->GetDx();
965  planeHits.DY = pHit->GetDy();
966  planeHits.planeID = trdPlaneID;
967  fvTrdHitArr[trdPlaneIDN].push_back(planeHits);
968  }
969 
970  Draw2();
971  DrawAsk();
972 
973  for (Int_t i = 0; i < 12; i++) {
974  fvTrdHitArr[i].clear();
975  }
976 
977 #endif
978 
979 
980  //-------------------------------------------------------
981 
982 
983  //-----------------------------------------------------------------
984  Double_t rtime;
985  Double_t ctime;
986  Double_t qtime = Double_t(fNofEvents);
987 
988  cout << endl
989  << "\n-I- CbmL1CATrdTrackFinderSA::DoFind " << endl
990  << "--------------------------------------------------" << endl
991  << " ### Time" << endl
992  << "--------------------------------------------------" << endl;
993  cout << " Do find: ";
994 
995  rtime = doFind.RealTime();
996  ctime = doFind.CpuTime();
997  totDoFind += ctime;
998  printf("RealTime=%4.2f s, CpuTime=%4.2f s, TotCpu=%4.2f s, %4.2f/event\n",
999  rtime,
1000  ctime,
1001  totDoFind,
1002  totDoFind / qtime);
1003 
1004  cout << " Sort Hits: ";
1005  rtime = sortHits.RealTime();
1006  ctime = sortHits.CpuTime();
1007  totSortHits += ctime;
1008  printf("RealTime=%4.2f s, CpuTime=%4.2f s, TotCpu=%4.2f s, %4.2f/event\n",
1009  rtime,
1010  ctime,
1011  totSortHits,
1012  totSortHits / qtime);
1013 
1014  cout << " Create SPs: ";
1015  rtime = createSPs.RealTime();
1016  ctime = createSPs.CpuTime();
1017  totCreateSPs += ctime;
1018  printf("RealTime=%4.2f s, CpuTime=%4.2f s, TotCpu=%4.2f s, %4.2f/event\n",
1019  rtime,
1020  ctime,
1021  totCreateSPs,
1022  totCreateSPs / qtime);
1023 
1024  cout << " Sort SPs: ";
1025  rtime = sortSPs.RealTime();
1026  ctime = sortSPs.CpuTime();
1027  totSortSPs += ctime;
1028  printf("RealTime=%4.2f s, CpuTime=%4.2f s, TotCpu=%4.2f s, %4.2f/event\n",
1029  rtime,
1030  ctime,
1031  totSortSPs,
1032  totSortSPs / qtime);
1033 
1034  cout << " Create segments: ";
1035  rtime = createSegments.RealTime();
1036  ctime = createSegments.CpuTime();
1037  totCreateSegments += ctime;
1038  printf("RealTime=%4.2f s, CpuTime=%4.2f s, TotCpu=%4.2f s, %4.2f/event\n",
1039  rtime,
1040  ctime,
1041  totCreateSegments,
1042  totCreateSegments / qtime);
1043 
1044  cout << " Find friend: ";
1045  rtime = findNeighbour.RealTime();
1046  ctime = findNeighbour.CpuTime();
1047  totFindNeighbour += ctime;
1048  printf("RealTime=%4.2f s, CpuTime=%4.2f s, TotCpu=%4.2f s, %4.2f/event\n",
1049  rtime,
1050  ctime,
1051  totFindNeighbour,
1052  totFindNeighbour / qtime);
1053 
1054  cout << " Tag Tracklets: ";
1055  rtime = tagTracks.RealTime();
1056  ctime = tagTracks.CpuTime();
1057  totTagTracks += ctime;
1058  printf("RealTime=%4.2f s, CpuTime=%4.2f s, TotCpu=%4.2f s, %4.2f/event\n",
1059  rtime,
1060  ctime,
1061  totTagTracks,
1062  totTagTracks / qtime);
1063 
1064  cout << " Create Tracks: ";
1065  rtime = createTracks.RealTime();
1066  ctime = createTracks.CpuTime();
1067  totCreateTracks += ctime;
1068  printf("RealTime=%4.2f s, CpuTime=%4.2f s, TotCpu=%4.2f s, %4.2f/event\n",
1069  rtime,
1070  ctime,
1071  totCreateTracks,
1072  totCreateTracks / qtime);
1073 
1074  cout << "Refitting Tracks: ";
1075  rtime = refittingKF.RealTime();
1076  ctime = refittingKF.CpuTime();
1077  totRefittingKF += ctime;
1078  printf("RealTime=%4.2f s, CpuTime=%4.2f s, TotCpu=%4.2f s, %4.2f/event\n",
1079  rtime,
1080  ctime,
1081  totRefittingKF,
1082  totRefittingKF / qtime);
1083 
1084  cout << " Second Loop: ";
1085  rtime = secondLoopTime.RealTime();
1086  ctime = secondLoopTime.CpuTime();
1087  totSecondLoopTime += ctime;
1088  printf("RealTime=%4.2f s, CpuTime=%4.2f s, TotCpu=%4.2f s, %4.2f/event\n",
1089  rtime,
1090  ctime,
1091  totSecondLoopTime,
1092  totSecondLoopTime / qtime);
1093 
1094  cout << " Third Loop: ";
1095  rtime = thirdLoopTime.RealTime();
1096  ctime = thirdLoopTime.CpuTime();
1097  totThirdLoopTime += ctime;
1098  printf("RealTime=%4.2f s, CpuTime=%4.2f s, TotCpu=%4.2f s, %4.2f/event\n",
1099  rtime,
1100  ctime,
1101  totThirdLoopTime,
1102  totThirdLoopTime / qtime);
1103 
1104  cout << "Deleting Objects: ";
1105  rtime = delTime.RealTime();
1106  ctime = delTime.CpuTime();
1107  totDelTime += ctime;
1108  printf("RealTime=%4.2f s, CpuTime=%4.2f s, TotCpu=%4.2f s, %4.2f/event\n",
1109  rtime,
1110  ctime,
1111  totDelTime,
1112  totDelTime / qtime);
1113 
1114  cout << "[Second Loop]" << endl;
1115  cout << " Create SPs: ";
1116  rtime = createSPs_SL.RealTime();
1117  ctime = createSPs_SL.CpuTime();
1118  totCreateSPs_SL += ctime;
1119  printf("RealTime=%4.2f s, CpuTime=%4.2f s, TotCpu=%4.2f s, %4.2f/event\n",
1120  rtime,
1121  ctime,
1122  totCreateSPs_SL,
1123  totCreateSPs_SL / qtime);
1124 
1125 
1126  cout << endl
1127  << "-I- CbmL1CATrdTrackFinderSA::DoFind " << endl
1128  << "--------------------------------------------------" << endl
1129  << " ### END: Track constructing phase completed" << endl
1130  << "--------------------------------------------------" << endl;
1131 
1132  fNoEvTime->Fill(fArrayTrdTrack->GetEntriesFast(), doFind.RealTime());
1133 
1134  return 1;
1135 }
1136 // -----------------------------------------------------------------------
1137 
1138 void CbmL1CATrdTrackFinderSA::DeleteTracklets(vector<CbmL1TrdTracklet4*> vect) {
1139 
1140  vector<CbmL1TrdTracklet4*>::iterator it;
1141  for (it = vect.begin(); it != vect.end(); it++) {
1142  delete (*it);
1143  }
1144 }
1145 
1146 // -----------------------------------------------------------------------
1147 
1148 void CbmL1CATrdTrackFinderSA::DeleteTracklets(vector<CbmL1TrdTracklet*> vect) {
1149  vector<CbmL1TrdTracklet*>::iterator it;
1150  for (it = vect.begin(); it != vect.end(); it++) {
1151  delete (*it);
1152  }
1153 }
1154 
1155 
1156 // ---------------------- Create histogramms -----------------------------
1157 void CbmL1CATrdTrackFinderSA::CreateHistogramms() {
1158 
1159  fUsedHitsPerPlane =
1160  new TH1F("h_UsedHitsPerPlane", "h_UsedHitsPerPlane", 13, 0, 13);
1161  fUnUsedHitsPerPlane =
1162  new TH1F("h_UnUsedHitsPerPlane", "h_UnUsedHitsPerPlane", 13, 0, 13);
1163 
1164 
1165  // Normalized distance to hit
1166  fh_chi2hit =
1167  new TH1F("h_chi2hit", "Normalized distance to hit", 500, 0., 50.);
1168  fh_chi2hit_plane = new TH2F("h_chi2hit_plane",
1169  "Normalized distance to hit vs. plane number",
1170  500,
1171  0.,
1172  50.,
1173  12,
1174  0.,
1175  12.);
1176 
1177  fDistLongX = new TH1F("DistLongX", "DistLongX", 100, -200, 200);
1178  fDistLongY = new TH1F("DistLongY", "DistLongY", 100, -200, 200);
1179  fDistShortX = new TH1F("DistShortX", "DistShortX", 100, -10, 10);
1180  fDistShortY = new TH1F("DistShortY", "DistShortY", 100, -10, 10);
1181 
1182  fDistLongBX = new TH1F("DistLongBX", "DistLongBX", 100, -200, 200);
1183  fDistLongBY = new TH1F("DistLongBY", "DistLongBY", 100, -200, 200);
1184  fDistShortBX = new TH1F("DistShortBX", "DistShortBX", 100, -200, 200);
1185  fDistShortBY = new TH1F("DistShortBY", "DistShortBY", 100, -200, 200);
1186 
1187  fDistY12 = new TH1F("DistY12", "DistY12", 100, -1, 1);
1188 
1189  fMomDistLongPrimaryX = new TH2F(
1190  "MomDistLongPrimaryX", "MomDistLongPrimaryX", 100, 0, 10, 100, -30, 30);
1191  fMomDistLongPrimaryY = new TH2F(
1192  "MomDistLongPrimaryY", "MomDistLongPrimaryY", 100, 0, 10, 100, -30, 30);
1193  fMomDistExtrapolPrimaryX = new TH2F("MomDistExtrapolPrimaryX",
1194  "MomDistExtrapolPrimaryX",
1195  100,
1196  0,
1197  10,
1198  200,
1199  -20,
1200  20);
1201  fMomDistExtrapolPrimaryY = new TH2F("MomDistExtrapolPrimaryY",
1202  "MomDistExtrapolPrimaryY",
1203  100,
1204  0,
1205  10,
1206  200,
1207  -20,
1208  20);
1209 
1210  fMomDistLongExtraX = new TH2F(
1211  "MomDistLongExtraX", "MomDistLongExtraX", 100, 0, 10, 100, -30, 30);
1212  fMomDistLongExtraY = new TH2F(
1213  "MomDistLongExtraY", "MomDistLongExtraY", 100, 0, 10, 100, -30, 30);
1214  fMomDistExtrapolExtraX = new TH2F(
1215  "MomDistExtrapolExtraX", "MomDistExtrapolExtraX", 100, 0, 10, 200, -20, 20);
1216  fMomDistExtrapolExtraY = new TH2F(
1217  "MomDistExtrapolExtraY", "MomDistExtrapolExtraY", 100, 0, 10, 200, -20, 20);
1218 
1219  fMomDistShortPrimaryX = new TH2F(
1220  "MomDistShortPrimaryX", "MomDistShortPrimaryX", 100, 0, 10, 100, -5, 5);
1221  fMomDistShortPrimaryY = new TH2F(
1222  "MomDistShortPrimaryY", "MomDistShortPrimaryY", 100, 0, 10, 100, -5, 5);
1223  fMomDistShortExtraX = new TH2F(
1224  "MomDistShortExtraX", "MomDistShortExtraX", 100, 0, 10, 100, -5, 5);
1225  fMomDistShortExtraY = new TH2F(
1226  "MomDistShortExtraY", "MomDistShortExtraY", 100, 0, 10, 100, -5, 5);
1227 
1228  fDistY = new TH1F("DistY", "DistY", 1000, -10, 10);
1229  fDistX = new TH1F("DistX", "DistX", 1000, -10, 10);
1230 
1231  fPlane1Ydens = new TH1F("Plane1Ydens", "Plane1Ydens", 1000, -1000, 1000);
1232  fPlane5Ydens = new TH1F("Plane5Ydens", "Plane5Ydens", 1000, -1000, 1000);
1233  fPlane9Ydens = new TH1F("Plane9Ydens", "Plane9Ydens", 1000, -1000, 1000);
1234 
1235  fSPlengthMC = new TH1F("SPlengthMC", "SPlengthMC", 200, 0, 20);
1236  fSPlength = new TH1F("SPlength", "SPlength", 200, 0, 20);
1237 
1238  fYat0 = new TH1F("Yat0", "Yat0", 500, -500, 500);
1239  fYat0MC = new TH1F("Yat0MC", "Yat0MC", 500, -500, 500);
1240 
1241  fNoEvTime = new TH2F("NoEvTime", "NoEvTime", 1000, 0, 1000, 1000, 0, 5);
1242 
1243  fh_SP_xDiff_MC = new TH1F("fh_SP_xDiff_MC", "fh_SP_xDiff_MC", 400, -20, 20);
1244  fh_SP_yDiff_MC = new TH1F("fh_SP_yDiff_MC", "fh_SP_yDiff_MC", 400, -20, 20);
1245 
1246  fh_SP_xDiff_nMC =
1247  new TH1F("fh_SP_xDiff_nMC", "fh_SP_xDiff_nMC", 400, -20, 20);
1248  fh_SP_yDiff_nMC =
1249  new TH1F("fh_SP_yDiff_nMC", "fh_SP_yDiff_nMC", 400, -20, 20);
1250 }
1251 // -----------------------------------------------------------------------
1252 
1253 
1254 // -------------------- Activates data branches --------------------------
1255 void CbmL1CATrdTrackFinderSA::DataBranches() {
1256 
1257  // Get pointer to the ROOT manager
1258  FairRootManager* rootMgr = FairRootManager::Instance();
1259  if (NULL == rootMgr) {
1260  cout << "-E- CbmL1CATrdTrackFinderSA::DataBranches : "
1261  << "ROOT manager is not instantiated" << endl;
1262  return;
1263  }
1264 
1265  fArrayTrdHit = L1_DYNAMIC_CAST<TClonesArray*>(rootMgr->GetObject("TrdHit"));
1266  if (!fArrayTrdHit) {
1267  cout << "-E- CbmL1CATrdTrackFinderSA::Init: No TrdHit array!" << endl;
1268  return;
1269  }
1270 }
1271 // -----------------------------------------------------------------------
1272 
1273 
1274 // ------------------- Determines the TRD layout -------------------------
1275 void CbmL1CATrdTrackFinderSA::TrdLayout() {
1276 
1277  // Get the pointer to the singleton FairRunAna object
1278  FairRunAna* ana = FairRunAna::Instance();
1279  if (NULL == ana) {
1280  cout << "-E- CbmL1CATrdTrackFinderSA::TrdLayout :"
1281  << " no FairRunAna object!" << endl;
1282  return;
1283  }
1284  // Get the pointer to run-time data base
1285  FairRuntimeDb* rtdb = ana->GetRuntimeDb();
1286  if (NULL == rtdb) {
1287  cout << "-E- CbmL1CATrdTrackFinderSA::TrdLayout :"
1288  << " no runtime database!" << endl;
1289  return;
1290  }
1291  // Get the pointer to container of base parameters
1292  FairBaseParSet* baseParSet =
1293  L1_DYNAMIC_CAST<FairBaseParSet*>(rtdb->getContainer("FairBaseParSet"));
1294  if (NULL == baseParSet) {
1295  cout << "-E- CbmL1CATrdTrackFinderSA::TrdLayout :"
1296  << " no container of base parameters!" << endl;
1297  return;
1298  }
1299 
1300  TrdPar = L1_DYNAMIC_CAST<CbmGeoTrdPar*>(rtdb->findContainer("CbmGeoTrdPar"));
1301  TObjArray* Nodes = TrdPar->GetGeoSensitiveNodes();
1302  for (Int_t i = 0; i < Nodes->GetEntries(); i++) {
1303  FairGeoNode* node = dynamic_cast<FairGeoNode*>(Nodes->At(i));
1304  //FairGeoNode *node = (FairGeoNode*) Nodes->At(i);
1305  if (!node) continue;
1306 
1307  TString name = node->getName();
1308  //TString Sname = node->getShapePointer()->GetName();
1309  FairGeoVector nodeV = node->getLabTransform()->getTranslation(); //in cm
1310 
1311  // FairGeoBasicShape *gShapePointer = (FairGeoBasicShape *)node->getShapePointer();
1312  //gShapePointer->printParam();
1313 
1314  // node->print();
1315 
1316  //Int_t id = node->getMCid();
1317  //cout <<"name: "<< name <<"\tid: "<< id << endl;
1318  cout << " name: " << name << "\t(z): " << nodeV.Z() << endl;
1319 
1320  // if(name.Contains("trd13")) fTrd13_Z = nodeV.Z();
1321  // if(name.Contains("trd14")) fTrd14_Z = nodeV.Z();
1322  // if(name.Contains("trd21")) fTrd21_Z = nodeV.Z();
1323  // if(name.Contains("trd24")) fTrd24_Z = nodeV.Z();
1324  // if(name.Contains("trd31")) fTrd31_Z = nodeV.Z();
1325 
1326  if (name == "trd1gas#3") fTrd13_Z = Int_t(nodeV.Z());
1327  if (name == "trd1gas#4") fTrd14_Z = Int_t(nodeV.Z());
1328  if (name == "trd2gas#1") fTrd21_Z = Int_t(nodeV.Z());
1329  if (name == "trd2gas#4") fTrd24_Z = Int_t(nodeV.Z());
1330  if (name == "trd3gas#1") fTrd31_Z = Int_t(nodeV.Z());
1331 
1332  // if(name == "trd13") fTrd13_Z = nodeV.Z();
1333  // if(name == "trd14") fTrd14_Z = nodeV.Z();
1334  // if(name == "trd21") fTrd21_Z = nodeV.Z();
1335  // if(name == "trd24") fTrd24_Z = nodeV.Z();
1336  // if(name == "trd31") fTrd31_Z = nodeV.Z();
1337 
1338  if (name == "trd1gas#1") geoLayer.Z[0] = nodeV.Z();
1339  if (name == "trd1gas#2") geoLayer.Z[1] = nodeV.Z();
1340  if (name == "trd1gas#3") geoLayer.Z[2] = nodeV.Z();
1341  if (name == "trd1gas#4") geoLayer.Z[3] = nodeV.Z();
1342  if (name == "trd2gas#1") geoLayer.Z[4] = nodeV.Z();
1343  if (name == "trd2gas#2") geoLayer.Z[5] = nodeV.Z();
1344  if (name == "trd2gas#3") geoLayer.Z[6] = nodeV.Z();
1345  if (name == "trd2gas#4") geoLayer.Z[7] = nodeV.Z();
1346  if (name == "trd3gas#1") geoLayer.Z[8] = nodeV.Z();
1347  if (name == "trd3gas#2") geoLayer.Z[9] = nodeV.Z();
1348  if (name == "trd3gas#3") geoLayer.Z[10] = nodeV.Z();
1349  if (name == "trd3gas#4") geoLayer.Z[11] = nodeV.Z();
1350  }
1351 
1352 #ifdef DRAW
1353  for (int i = 0; i < 4; i++) {
1354  geoLayer.Y[i] = 273.1;
1355  }
1356  for (int i = 4; i < 8; i++) {
1357  geoLayer.Y[i] = 396.0;
1358  }
1359  for (int i = 8; i < 12; i++) {
1360  geoLayer.Y[i] = 518.9;
1361  }
1362 
1363  Int_t scaleZ = 300;
1364  Int_t scaleZA = 0;
1365  for (int i = 0; i < 12; i++) {
1366  if (i == 4) scaleZA = 1000;
1367  if (i == 8) scaleZA *= 2;
1368 
1369  geoLayer.scale[i] += ((i + 1) * scaleZ) + scaleZA;
1370  geoLayer.Z[i] += geoLayer.scale[i];
1371  }
1372 #endif
1373 
1374  /*
1375  TrdPar = (CbmGeoTrdPar*) (rtdb->findContainer("CbmGeoTrdPar"));
1376  //TObjArray *Nodes = TrdPar->GetGeoSensitiveNodes();
1377  Nodes = TrdPar->GetGeoSensitiveNodes();
1378  for( Int_t i=0;i<Nodes->GetEntries(); i++) {
1379  FairGeoNode *node = dynamic_cast<FairGeoNode*> (Nodes->At(i));
1380  if ( !node ) continue;
1381 
1382  TString name = node->getName();
1383  FairGeoVector nodeV = node->getLabTransform()->getTranslation(); //in cm
1384 
1385  Int_t uid = i+1;
1386  cout <<"name, uid, Z: "<< name <<", "<< uid <<", "<< nodeV.Z() << endl;
1387 
1388  if(uid == 3) fTrd13_Z = nodeV.Z();
1389  if(uid == 4) fTrd14_Z = nodeV.Z();
1390  if(uid == 5) fTrd21_Z = nodeV.Z();
1391  if(uid == 8) fTrd24_Z = nodeV.Z();
1392  if(uid == 9) fTrd31_Z = nodeV.Z();
1393  }
1394  */
1395 
1396  // Get the pointer to detector list
1397  TObjArray* detList = baseParSet->GetDetList();
1398  if (NULL == detList) {
1399  cout << "-E- CbmL1CATrdTrackFinderSA::TrdLayout :"
1400  << " no detector list!" << endl;
1401  return;
1402  }
1403  // Find TRD detector
1404  FairDetector* trd =
1405  L1_DYNAMIC_CAST<FairDetector*>(detList->FindObject("TRD"));
1406  if (NULL == trd) {
1407  cout << "-E- CbmL1CATrdTrackFinderSA::TrdLayout :"
1408  << " no TRD detector!" << endl;
1409  return;
1410  }
1411  // Determine the geometry version
1412  TString name = trd->GetGeometryFileName();
1413  if (name.Contains("9")) {
1414  cout << "-I- CbmL1CATrdTrackFinderSA::TrdLayout :"
1415  << " TRD geometry : 3x3." << endl;
1416  fNoTrdStations = 3;
1417  fNoTrdPerStation = 3;
1418  } else if (name.Contains("12")) {
1419  cout << "-I- CbmL1CATrdTrackFinderSA::TrdLayout :"
1420  << " TRD geometry : 3x4." << endl;
1421  fNoTrdStations = 3;
1422  fNoTrdPerStation = 4;
1423  } else if (name.Contains("6x2")) {
1424  cout << "-I- CbmL1CATrdTrackFinderSA::TrdLayout :"
1425  << " TRD geometry : 6x2." << endl;
1426  fNoTrdStations = 6;
1427  fNoTrdPerStation = 2;
1428  }
1429  fNoTrdStations = 3;
1430  fNoTrdPerStation = 4;
1431 }
1432 // -----------------------------------------------------------------------
1433 
1434 
1435 // ------------------------- Write histogramms ---------------------------
1436 void CbmL1CATrdTrackFinderSA::WriteHistogramms() {
1437 
1438  fh_chi2hit->Write();
1439  fh_chi2hit_plane->Write();
1440 
1441  fDistLongBX->Write();
1442  fDistLongBY->Write();
1443  fDistShortBX->Write();
1444  fDistShortBY->Write();
1445 
1446  fDistY12->Write();
1447 
1448  //----------------------------------------------------------
1449  fMomDistLongPrimaryX->Write();
1450  fMomDistLongPrimaryY->Write();
1451  fMomDistLongExtraX->Write();
1452  fMomDistLongExtraY->Write();
1453 
1454  fMomDistExtrapolPrimaryX->Write();
1455  fMomDistExtrapolPrimaryY->Write();
1456  fMomDistExtrapolExtraX->Write();
1457  fMomDistExtrapolExtraY->Write();
1458 
1459  fMomDistShortPrimaryX->Write();
1460  fMomDistShortPrimaryY->Write();
1461  fMomDistShortExtraX->Write();
1462  fMomDistShortExtraY->Write();
1463 
1464  fDistLongX->Write();
1465  fDistLongY->Write();
1466  fDistShortX->Write();
1467  fDistShortY->Write();
1468 
1469  fDistY->Write();
1470  fDistX->Write();
1471 
1472  fPlane1Ydens->Write();
1473  fPlane5Ydens->Write();
1474  fPlane9Ydens->Write();
1475 
1476  fSPlength->Write();
1477  fSPlengthMC->Write();
1478 
1479  fYat0->Write();
1480  fYat0MC->Write();
1481 
1482  fNoEvTime->Write();
1483 
1484  fh_SP_xDiff_MC->Write();
1485  fh_SP_yDiff_MC->Write();
1486 
1487  fh_SP_xDiff_nMC->Write();
1488  fh_SP_yDiff_nMC->Write();
1489 
1490 
1491  //---------------------------------------------------------------------
1492  map<Int_t, Int_t>::iterator iUsedHits;
1493  map<Int_t, Int_t>::iterator iUnUsedHits;
1494  map<Int_t, Int_t>::iterator iTotHits;
1495  Double_t nContent;
1496 
1497  for (iUsedHits = fRUsedHits.begin(); iUsedHits != fRUsedHits.end();
1498  iUsedHits++) {
1499 
1500  iTotHits = fTotHits.find(iUsedHits->first);
1501  nContent = (iUsedHits->second * 100) / iTotHits->second;
1502  // cout << iUsedHits->first <<", "
1503  // << iUsedHits->second <<", "
1504  // << iTotHits->second <<", "
1505  // << nContent << endl;
1506  fUsedHitsPerPlane->SetBinContent((iUsedHits->first), nContent);
1507  }
1508  for (iUnUsedHits = fRUnUsedHits.begin(); iUnUsedHits != fRUnUsedHits.end();
1509  iUnUsedHits++) {
1510 
1511  iTotHits = fTotHits.find(iUnUsedHits->first);
1512  nContent = (iUnUsedHits->second * 100) / iTotHits->second;
1513  fUnUsedHitsPerPlane->SetBinContent((iUnUsedHits->first), nContent);
1514  }
1515  fUsedHitsPerPlane->Write();
1516  fUnUsedHitsPerPlane->Write();
1517  //---------------------------------------------------------------------
1518 }
1519 // -----------------------------------------------------------------------
1520 
1521 
1522 // -----------------------------------------------------------------------
1523 Double_t CbmL1CATrdTrackFinderSA::FitLinear(CbmTrdTrack* tr,
1524  Int_t /*var*/ = 1) {
1525  //fit using a least square method
1526 
1527 
1528  // cout << ">>> Trying to fit a track..." << endl;
1529  Int_t noHits = tr->GetNofHits();
1530  // cout << "No Hits in the track : " << noHits << endl;
1531  Int_t iHit;
1532  CbmTrdHit* hit;
1533 
1534  Double_t C1[12];
1535  Double_t C2[12];
1536  Double_t Z[12];
1537 
1538  Double_t n = 12;
1539 
1540  for (int i = 0; i < 12; i++) {
1541  C1[i] = 0;
1542  C2[i] = 0;
1543  Z[i] = 0;
1544  }
1545 
1546  // Int_t ind = 0;
1547 
1548  for (int i = 0; i < noHits; i++) {
1549 
1550  iHit = tr->GetHitIndex(i);
1551  hit = L1_DYNAMIC_CAST<CbmTrdHit*>(fArrayTrdHit->At(iHit));
1552 
1553  C1[i] = hit->GetX();
1554  C2[i] = hit->GetY();
1555  Z[i] = hit->GetZ();
1556  }
1557 
1558  Double_t sumXZ = 0, sumYZ = 0, sumX = 0, sumY = 0, sumZx = 0, sumZy = 0,
1559  sumX2 = 0, sumY2 = 0, sumZx2 = 0, sumZy2 = 0;
1560 
1561  for (int i = 0; i < 12; i += 2) {
1562  // cout << "C1[" << i << "]=" << C1[i] << endl;
1563  // cout << "C2[" << i << "]=" << C2[i] << endl;
1564 
1565  // if(!(i % 2))
1566  {
1567  sumXZ += C1[i] * Z[i];
1568 
1569  sumX += C1[i];
1570  sumX2 += C1[i] * C1[i];
1571 
1572  sumZx += Z[i];
1573  sumZx2 += Z[i] * Z[i];
1574  }
1575  }
1576 
1577  for (int i = 1; i < 12; i += 2) {
1578  // cout << "C1[" << i << "]=" << C1[i] << endl;
1579  // cout << "C2[" << i << "]=" << C2[i] << endl;
1580 
1581  // if(!((1i % 2))
1582  {
1583  sumYZ += C2[i] * Z[i];
1584 
1585  sumY += C2[i];
1586  sumY2 += C2[i] * C2[i];
1587 
1588  sumZy += Z[i];
1589  sumZy2 += Z[i] * Z[i];
1590  }
1591  }
1592 
1593 
1594  cout << "";
1595  // Double_t a = (n*sumXZ - sumX*sumZ)/(n*sumX2 - sumX*sumX);
1596  // Double_t b = (sumZ - a*sumX)/n;
1597  Double_t r1 =
1598  (n * sumXZ - sumX * sumZx)
1599  / sqrt((n * sumX2 - sumX * sumX) * (n * sumZx2 - sumZx * sumZx));
1600  // Double_t Sa = sqrt((n*(sumY2 -a*sumXZ-b*sumY))/((n-2)*(n*sumX2-sumX*sumX)));
1601  // Double_t Sb = Sa*sqrt(sumX2/n);
1602 
1603 
1604  Double_t r2 =
1605  (n * sumYZ - sumY * sumZy)
1606  / sqrt((n * sumY2 - sumY * sumY) * (n * sumZy2 - sumZy * sumZy));
1607 
1608 
1609  // cout << "Linear coefficients: a= "
1610  // << a << ", Sa= " << Sa << ", b= " << b << ", Sb= " << Sb << endl
1611  // << ", r=" << r << ", chi2=" << tr->GetChi2() << endl;
1612 
1613 
1614  return sqrt(r1 * r1 + r2 * r2);
1615 }
1616 // -----------------------------------------------------------------------
1617 
1618 
1619 // -----------------------------------------------------------------------
1620 Double_t CbmL1CATrdTrackFinderSA::DistTwoTrackletsY(Int_t iIndexFirst,
1621  Int_t iIndexSecond,
1622  Double_t zed) {
1623  //returns an extrapolated coordinate
1624  CbmTrdHit *pHitFirst, *pHitSecond;
1625  Double_t dist;
1626 
1627  Double_t Y1 = 0;
1628  Double_t Z1 = 0;
1629 
1630  // Get hits from first and second plane in station
1631  if (iIndexFirst != -1) {
1632  pHitFirst = L1_DYNAMIC_CAST<CbmTrdHit*>(fArrayTrdHit->At(iIndexFirst));
1633  Y1 = pHitFirst->GetY();
1634  Z1 = pHitFirst->GetZ();
1635  }
1636 
1637  pHitSecond = L1_DYNAMIC_CAST<CbmTrdHit*>(fArrayTrdHit->At(iIndexSecond));
1638 
1639  // Get position Y & Z of hits
1640 
1641  Double_t Y4 = pHitSecond->GetY();
1642  Double_t Z4 = pHitSecond->GetZ();
1643 
1644  Double_t a = (Y4 - Y1) / (Z4 - Z1);
1645  Double_t b = Y1 - a * Z1;
1646  Double_t Y = a * zed + b; // - Y1;
1647  dist = Y;
1648  // cout << "From the coords (y1,z1) = ("
1649  // << Y1 << ", " << Z1 << endl
1650  // << "From the coords (y2,z2) = ("
1651  // << Y4 << ", " << Z4 << endl
1652  // << " we get the value = "
1653  // << dist
1654  // << endl;
1655  return dist;
1656 }
1657 // -----------------------------------------------------------------------
1658 
1659 
1660 // -----------------------------------------------------------------------
1661 
1662 Double_t CbmL1CATrdTrackFinderSA::DistTwoTrackletsX(Int_t iIndexFirst,
1663  Int_t iIndexSecond,
1664  Double_t zed) {
1665  CbmTrdHit *pHitFirst, *pHitSecond;
1666  Double_t dist;
1667 
1668  Double_t X1 = 0;
1669  Double_t Z1 = 0;
1670 
1671  // Get hits from first and second plane in station
1672  if (iIndexFirst != -1) {
1673  pHitFirst = L1_DYNAMIC_CAST<CbmTrdHit*>(fArrayTrdHit->At(iIndexFirst));
1674  X1 = pHitFirst->GetX();
1675  Z1 = pHitFirst->GetZ();
1676  }
1677 
1678  pHitSecond = L1_DYNAMIC_CAST<CbmTrdHit*>(fArrayTrdHit->At(iIndexSecond));
1679 
1680 
1681  Double_t X4 = pHitSecond->GetX();
1682  Double_t Z4 = pHitSecond->GetZ();
1683 
1684  Double_t a = (X4 - X1) / (Z4 - Z1);
1685  Double_t b = X1 - a * Z1;
1686  // Double_t X = (dz)/(dx)*(zed - Z1);
1687  Double_t X = a * zed + b;
1688  // = dy/dz * (zed + Z1) - X1;
1689  dist = X;
1690  // cout << "From the coords (x1,z1) = ("
1691  // << X1 << ", " << Z1 << endl
1692  // << "From the coords (x4,z4) = ("
1693  // << X4 << ", " << Z4 << endl
1694  // << "and at zed = " << zed << endl
1695  // << " we get the value = "
1696  // << dist
1697  // << endl;
1698  return dist;
1699 }
1700 // -----------------------------------------------------------------------
1701 
1702 // -----------------------------------------------------------------------
1703 void CbmL1CATrdTrackFinderSA::FindNeighbour(vector<CbmL1TrdTracklet4*>& v1,
1704  vector<CbmL1TrdTracklet4*>& v2,
1705  Double_t dY,
1706  Double_t dX) {
1707 
1708  vector<CbmL1TrdTracklet4*>::iterator iv1;
1709  CbmL1TrdTracklet4* tr1;
1710 
1711  Double_t extY = 0,
1712  extX = 0, //two extrapolated coords from the first tracklet
1713  mesY = 0, mesX = 0, //two measured coords from the second tracklet
1714  aminY, amaxY, aminX, amaxX;
1715  Int_t indexA, indexB;
1716  Int_t Left = 0, Right, Mid = 0;
1717 
1718  Int_t il = 0;
1719  for (iv1 = v1.begin(); iv1 != v1.end(); iv1++) {
1720  il++;
1721  tr1 = *iv1;
1722  extY = tr1->GetExt(0);
1723  extX = tr1->GetExt(1);
1724  indexA = tr1->GetIndex();
1725  amaxY = extY + dY;
1726  aminY = extY - dY;
1727  amaxX = extX + dX;
1728  aminX = extX - dX;
1729 
1730  Left = 0;
1731  Right = v2.size();
1732  Mid = 0;
1733 
1734  // ----------- looking by bisection ------------------
1735  while (Left < Right) {
1736  Mid = (Left + Right) / 2;
1737  mesY = v2[Mid]->GetCoord(0);
1738 
1739  if (amaxY < mesY) {
1740  Left = Mid + 1;
1741  } else {
1742  Right = Mid - 1;
1743  }
1744 
1745  /*
1746  cout << il << ": Size: "<< v1.size()
1747  <<" mesY: "<< v2[Mid]->GetCoord(0)
1748  <<" extY: "<< extY
1749  <<" Mid: "<< Mid
1750  << endl;
1751  */
1752  }
1753  //----------------------------------------------------
1754 
1755  Int_t size = v2.size();
1756  for (Int_t i = Mid; i < size; i++) {
1757  mesY = v2[i]->GetCoord(0);
1758  mesX = v2[i]->GetCoord(1);
1759  //cout <<"(mesX,MesY) "<< mesX <<"\t"<< mesY << endl;
1760  indexB = v2[i]->GetIndex();
1761  //cout <<"indexB "<< indexB << endl;
1762 
1763  if (mesY > amaxY) continue;
1764  if (mesY < aminY) break;
1765 
1766  // cout <<"----------- Setting vAccost [Right|Left]"<< endl;
1767 
1768  // cout << aminX <<", "<< mesX <<", "<< amaxX << endl;
1769 
1770  if (aminX < mesX && mesX < amaxX) {
1771  tr1->vAccostRight.push_back(indexB);
1772  v2[i]->vAccostLeft.push_back(indexA);
1773  // cout << aminX <<", "<< mesX <<", "<< amaxX << endl;
1774  }
1775  }
1776  }
1777 }
1778 // -----------------------------------------------------------------------
1779 
1780 
1781 // -----------------------------------------------------------------------
1782 Bool_t CbmL1CATrdTrackFinderSA::OverlapsHitsXY(Int_t posA, Int_t posB) {
1783  Bool_t overlap = false;
1784 
1785  CbmTrdHit *pHitA, *pHitB;
1786 
1787  Double_t hitPosA_X, hitPosA_Y, hitPosB_X, hitPosB_Y,
1788  dMul1 = 3, //sigma multiplier for calculating the range of overlapping hits
1789  dMul2 = 3; //sigma multiplier for calculating the range of overlapping hits
1790 
1791  Double_t hitPosA_DX, hitPosA_DY, hitPosB_DX, hitPosB_DY;
1792 
1793  pHitA = L1_DYNAMIC_CAST<CbmTrdHit*>(fArrayTrdHit->At(posA));
1794  hitPosA_X = pHitA->GetX();
1795  hitPosA_Y = pHitA->GetY();
1796  hitPosA_DX = pHitA->GetDx();
1797  hitPosA_DY = pHitA->GetDy();
1798 
1799  pHitB = L1_DYNAMIC_CAST<CbmTrdHit*>(fArrayTrdHit->At(posB));
1800  hitPosB_X = pHitB->GetX();
1801  hitPosB_Y = pHitB->GetY();
1802  hitPosB_DX = pHitB->GetDx();
1803  hitPosB_DY = pHitB->GetDy();
1804 
1805  // cout
1806  // << "Overlap function...\n"
1807  // << "X1=" << hitPosA_X << ", dX1=" << hitPosA_DX << ", Y1=" << hitPosA_Y << ", dY1=" << hitPosA_DY
1808  // << "X2=" << hitPosB_X << ", dX2=" << hitPosB_DX << ", Y2=" << hitPosB_Y << ", dY2=" << hitPosB_DY
1809  // << endl;
1810 
1811 
1812  if (((hitPosA_X + dMul1 * hitPosA_DX) >= (hitPosB_X - dMul2 * hitPosB_DX))
1813  && ((hitPosA_X - dMul1 * hitPosA_DX) <= (hitPosB_X + dMul2 * hitPosB_DX))
1814  && ((hitPosA_Y + dMul2 * hitPosA_DY) >= (hitPosB_Y - dMul1 * hitPosB_DY))
1815  && ((hitPosA_Y - dMul2 * hitPosA_DY)
1816  <= (hitPosB_Y + dMul1 * hitPosB_DY))) {
1817  overlap = true;
1818  }
1819 
1820  // CbmTrdPoint *pt1 = (CbmTrdPoint*)fMCPointArray->At(posA);
1821  // CbmTrdPoint *pt2 = (CbmTrdPoint*)fMCPointArray->At(posB);
1822 
1823  // if(overlap && (pt1->GetTrackID() == pt2->GetTrackID())) return true;
1824  // else return false;
1825 
1826  return overlap;
1827 }
1828 
1829 
1830 // -----------------------------------------------------------------------
1831 void CbmL1CATrdTrackFinderSA::TagSegments(
1832  vector<CbmL1TrdTracklet4*>& clTrackletsA,
1833  vector<CbmL1TrdTracklet4*>& clTrackletsB,
1834  Int_t /*noCombSegments*/) {
1835  //asigning numbers to each segment, 4- highest, 0 - lowest
1836 
1837  // cout << "TagSegments: engaging... " << endl;
1838 
1839  // vector<CbmL1TrdTracklet4*>::iterator itclTrackletsLeft;
1840  vector<CbmL1TrdTracklet4*>::iterator itclTrackletsRight;
1841 
1842  CbmL1TrdTracklet4* clSegRight;
1843 
1844  vector<Int_t> vLeft;
1845  vector<Int_t> vRight;
1846 
1847  vector<Int_t>::iterator ivLeft;
1848  // vector <Int_t>::iterator ivRight;
1849 
1850  Int_t valLeft = 0, valRight = 0;
1851 
1852  // CbmL1TrdTracklet4* clSegLeft; //segment farther from the target
1853  // CbmL1TrdTracklet4* clSegRight; //segment nearer to the target
1854 
1855 
1856  // cout << "Engaging Right loop " << endl;
1857  for (itclTrackletsRight = clTrackletsB.begin();
1858  itclTrackletsRight != clTrackletsB.end();
1859  itclTrackletsRight++) {
1860 
1861  clSegRight = *itclTrackletsRight;
1862  //cout << "---> X of tracklet = " << clSegRight->GetCoord(0) << endl;
1863  vLeft = clSegRight->vAccostLeft;
1864  //cout << "Size of Right vector = " << clSegRight->vAccostRight.size() << endl;
1865  //cout << "Size of Left vector = " << clSegRight->vAccostLeft.size() << endl;
1866  // iHitRightA = clSegRight->GetInd(0);
1867  // iHitRightB = clSegRight->GetInd(1);
1868 
1869  // cout << "Tracklets 912: iHitRight: " << iHitRight << "\t" << endl;
1870  if (vLeft.size() > 0) {
1871  // cout << "Processing Comrade vector of size " << vLeft.size() << endl;
1872  valRight = clSegRight->GetVal();
1873  if (valRight == 0) {
1874  clSegRight->SetVal(1);
1875  valRight++;
1876  }
1877  for (ivLeft = vLeft.begin(); ivLeft != vLeft.end(); ivLeft++) {
1878  valLeft = (clTrackletsA[*ivLeft])->GetVal();
1879  if (valLeft <= valRight) clTrackletsA[*ivLeft]->SetVal(valRight + 1);
1880  }
1881  }
1882 
1883  // if(noCombSegments == 1){
1884  // if((*itclTrackletsRight)->GetVal() != 2) clTrackletsB.erase(itclTrackletsRight);
1885  //}
1886 
1887  //cout << "clSegRight = " << clSegRight->GetVal() << endl;
1888  //if(isAlone == true) clSegRight->SetIsAlone(true);
1889  }
1890 }
1891 // -----------------------------------------------------------------------
1892 
1893 
1894 // -----------------------------------------------------------------------
1895 void CbmL1CATrdTrackFinderSA::CreateTracks(
1896  vector<CbmL1TrdTracklet4*> clTracklets14,
1897  vector<CbmL1TrdTracklet4*> clTracklets58,
1898  vector<CbmL1TrdTracklet4*> clTracklets912,
1899  set<Int_t>& globalSetUsedHits,
1900  Bool_t removeUsedHits,
1901  Bool_t competition,
1902  Int_t /*nrLoop*/) {
1903  //create long tracks from previously selected segments
1904 
1905  if (fVerbose > 2) {
1906  cout << "Inner area: " << endl
1907  << "--- size of fArrayTrdTrack = "
1908  << fArrayTrdTrack->GetEntriesFast() << endl
1909  << "--- size of globalSetUsedHits = " << globalSetUsedHits.size()
1910  << endl;
1911  }
1912 
1913  CbmL1TrdTracklet4 *clTrdSeg1, *clTrdSeg2;
1914 
1915  vector<CbmL1TrdTracklet4*>::iterator itclSeg3;
1916 
1917  Bool_t isEmpty = true;
1918 
1919  Int_t licz = 0, counter = 0;
1920 
1921  CbmTrdTrack* tr1 = NULL;
1922  // CbmTrdTrack* tr2 = NULL;
1923  // CbmTrdHit* hit1 = NULL;
1924  //CbmTrdHit* hit2 = NULL;
1925 
1926  vector<CbmTrdTrack*> vTempTrdTrackCand;
1927  vector<CbmTrdTrack*>::iterator ivTempTrdTrackCand;
1928 
1929  set<Int_t>::iterator iglobalSetUsedHits;
1930 
1931  vector<CbmTrdTrack*>::iterator ivTrdTrackCand;
1932 
1933  Int_t iFakeTrack = 0;
1934  Int_t iHit = 0;
1935  Bool_t bTrueTrack = true;
1936 
1937  FairTrackParam* trParam;
1938  trParam = new FairTrackParam();
1939  trParam->SetQp(0.);
1940 
1941  CbmVertex vtx;
1942  TVector3 vzero(0, 0, 0);
1943  vtx.Position(vzero);
1944 
1945  Int_t iSecondLoop = 1;
1946  Int_t trMax =
1947  1; //number of tracks from one tree accepted as tracks candidates
1948  Int_t noHitsAccepted =
1949  0; //number of used hits above which the tracks is recognized as fake; 0 = all tracks are fake.
1950 
1951  vector<TempTrackStruct> auxTrack;
1952  auxTrack.clear();
1953 
1954  Int_t trdInd;
1955  vector<TempTrackStruct>::iterator ikol;
1956 
1957  cout << "--- Engaging main loop..." << endl;
1958 
1959  if (fVerbose > 2)
1960  cout << "--- No of outer iterations to go: " << clTracklets14.size()
1961  << endl;
1962  // Int_t vSize = Int_t(clTracklets14.size());
1963  //for(int i = 0; i < vSize; i++) {
1964  // if( (clTracklets14[i])->GetVal() == 3) licz++;
1965  // }
1966  // cout << "* No of 3 in vector: " << licz << endl;
1967 
1968  Int_t trackArrayInd = fArrayTrdTrack->GetEntriesFast();
1969  Int_t trlsNo[] = {0, 0, 0};
1970 
1971  for (itclSeg3 = clTracklets14.begin(); itclSeg3 != clTracklets14.end();
1972  itclSeg3++) { //loop over tracklets with value of 3
1973  if ((*itclSeg3)->GetVal() != 3) {
1974  trlsNo[0]++;
1975  continue;
1976  }
1977  for (int i = 0; i < 12; i++)
1978  tempTrack.M[i] = 0;
1979 
1980  tempTrack.M[0] = (*itclSeg3)->GetInd(0);
1981  tempTrack.M[1] = (*itclSeg3)->GetInd(1);
1982  tempTrack.M[2] = (*itclSeg3)->GetInd(2);
1983  tempTrack.M[3] = (*itclSeg3)->GetInd(3);
1984 
1985  Int_t iInd2 = 0;
1986  Int_t size2 = Int_t((*itclSeg3)->vAccostRight.size());
1987  for (Int_t iSeg2 = 0; iSeg2 < size2;
1988  iSeg2++) { //loop over clTracklets with value of 2
1989  iInd2 = (*itclSeg3)->vAccostRight[iSeg2];
1990  clTrdSeg2 = clTracklets58[iInd2];
1991  if (clTrdSeg2->GetVal() != 2) continue;
1992 
1993  tempTrack.M[4] = clTrdSeg2->GetInd(0);
1994  tempTrack.M[5] = clTrdSeg2->GetInd(1);
1995  tempTrack.M[6] = clTrdSeg2->GetInd(2);
1996  tempTrack.M[7] = clTrdSeg2->GetInd(3);
1997 
1998  Int_t iInd1 = 0;
1999  Int_t size1 = Int_t(clTrdSeg2->vAccostRight.size());
2000  for (Int_t iSeg1 = 0; iSeg1 < size1;
2001  iSeg1++) { //loop over clTracklets with value of 1
2002  iInd1 = clTrdSeg2->vAccostRight[iSeg1];
2003  clTrdSeg1 = clTracklets912[iInd1];
2004  if (clTrdSeg1->GetVal() != 1) continue;
2005 
2006  tempTrack.M[8] = clTrdSeg1->GetInd(0);
2007  tempTrack.M[9] = clTrdSeg1->GetInd(1);
2008  tempTrack.M[10] = clTrdSeg1->GetInd(2);
2009  tempTrack.M[11] = clTrdSeg1->GetInd(3);
2010 
2011  counter++;
2012  isEmpty = false;
2013 
2014  /*
2015  //-------------------------------------
2016  tr2 = new CbmTrdTrack();
2017  tr2->SetParamFirst(*trParam);
2018  tr2->SetParamLast(*trParam);
2019 
2020  for(Int_t we = 0; we < 12; we++) {
2021  hit1 = L1_DYNAMIC_CAST<CbmTrdHit*>( fArrayTrdHit->At(tempTrack.M[we]);
2022  tr2->AddHit(tempTrack.M[we], hit1);
2023  }
2024  tr2->SortHits();
2025  FitKF(tr2);
2026  // trdTrackFitterKF->DoFit(tr2);
2027  tempTrack.Chi2 = tr2->GetChi2();
2028  delete tr2;
2029  //------------------------------------- */
2030 
2031  //tempTrack.Chi2 = FitTLinearFitter(tempTrack.M);
2032  //tempTrack.Chi2 = FitLinear(tempTrack.M, 1);
2033  //tempTrack.Chi2 = FitLSM(tempTrack.M);
2034  tempTrack.Chi2 = Fit(tempTrack.M);
2035 
2036  auxTrack.push_back(tempTrack);
2037 
2038  //if((counter%1000) == 0) {
2039  //printf("\rTracks: %i",counter);
2040  //fflush(stdout);
2041  //}
2042  } //end::loop 1
2043  } //end::loop 2
2044 
2045  //if(iSecondLoop%2 == 0)
2046  {
2047  iSecondLoop = 0;
2048  sort(auxTrack.begin(), auxTrack.end(), CompareChi2);
2049  //if(nrLoop == 1) trMax = 1000;
2050  //if(nrLoop == 2) trMax = 1000;
2051  // if(nrLoop == 1) cout <<"auxTrack.size(): "<< auxTrack.size() << endl;
2052  Int_t li = 0;
2053  //CbmTrdTrack *T3;
2054  //Double_t chiMax = 1000;
2055 
2056  for (ikol = auxTrack.begin(); ikol != auxTrack.end(); ikol++) {
2057  if (li >= trMax) break;
2058 
2059  tr1 = new CbmTrdTrack();
2060 
2061  for (Int_t we = 0; we < 12; we++) {
2062  trdInd = (*ikol).M[we];
2063  // hit1 = L1_DYNAMIC_CAST<CbmTrdHit*>( fArrayTrdHit->At(trdInd) );
2064 
2065  tr1->AddHit(trdInd, kTRDHIT);
2066  }
2067  // tr1->SortHits();
2068  //tr1->SetChi2((*ikol).Chi2);
2069  //cout << tempTrack.Chi2 << endl;
2070 
2071  /*// Set reliable Q/p value. Based on MC values ----------------
2072  CbmTrdPoint *pMCpt;
2073  CbmMCTrack *mcTr;
2074  Int_t ptIndex;
2075  Double_t mom;
2076  trdInd = tr1->GetTrdHitIndex(0);
2077  hit1 = L1_DYNAMIC_CAST<CbmTrdHit*>( fArrayTrdHit->At(trdInd);
2078  ptIndex = hit1->GetRefIndex();
2079  pMCpt = (CbmTrdPoint*) (fMCPointArray->At(ptIndex));
2080 
2081  mcTr = (CbmMCTrack*)fMCTrackArray->At(pMCpt->GetTrackID());
2082  mom = mcTr->GetMomentum().Mag();
2083  Int_t pdgCode = mcTr->GetPdgCode();
2084  Int_t charge = PdgToCharge(pdgCode);
2085  trParam->SetQp(charge/mom);
2086  //trParam->SetQp(1.);
2087  //cout <<"(pdg, Qp) "<< pdgCode <<", "<< charge/mom << endl;
2088  //-------------------------------------------------------------
2089  */
2090 
2091  tr1->SetParamFirst(trParam);
2092  tr1->SetParamLast(trParam);
2093 
2094  //FitKF(tr1);
2095  // if(tr1->GetChi2() < chiMax){
2096  // T3 = tr1;
2097  // chiMax = tr1->GetChi2();
2098  //}
2099  //if((*ikol).Chi2 < 3)
2100  vTempTrdTrackCand.push_back(tr1);
2101  li++;
2102  }
2103  //vTempTrdTrackCand.push_back(T3);
2104  auxTrack.clear();
2105  // cout << "------------------"<< endl;
2106  }
2107  //cout << iSecondLoop << endl;
2108  iSecondLoop++;
2109  } //loop over 3s
2110 
2111  //cout << clTracklets14.size() <<"\t"
2112  // << trlsNo[0] << endl;
2113 
2114  // if(nrLoop == 2) competition = false;
2115  if (competition) {
2116  refittingKF.Start();
2117  cout << "\n--- Refiting by KF..." << endl;
2118  /* Refit using KF fitter */
2119  //Int_t iKlops = 0;
2120  //cout << "vTempTrdTrackCand.size(): " << vTempTrdTrackCand.size() << endl;
2121  for (ivTempTrdTrackCand = vTempTrdTrackCand.begin();
2122  ivTempTrdTrackCand != vTempTrdTrackCand.end();
2123  ivTempTrdTrackCand++) {
2124 
2125  //FitKF(*ivTempTrdTrackCand);
2126  //trdTrackFitterKF->DoFit(*ivTempTrdTrackCand, &vtx);
2127  trdTrackFitterKF->DoFit(*ivTempTrdTrackCand);
2128  // Double_t mX = (*ivTempTrdTrackCand)->GetParamLast()->GetX();
2129  // Double_t mY = (*ivTempTrdTrackCand)->GetParamLast()->GetY();
2130  // Double_t mXY = sqrt(pow(mX,2)+pow(mY,2));
2131  // Double_t mChi2 = (*ivTempTrdTrackCand)->GetChi2();
2132  // (*ivTempTrdTrackCand)->SetChi2(mChi2*mX2);
2133  //cout <<"Chi2: "<< (*ivTempTrdTrackCand)->GetChi2() << endl;
2134  //iKlops++;
2135  //cout << iKlops <<": GetNofTrdHits(): " << (*ivTempTrdTrackCand)->GetNofTrdHits() << endl;
2136  }
2137  refittingKF.Stop();
2138 
2139  sort(
2140  vTempTrdTrackCand.begin(), vTempTrdTrackCand.end(), CompareChi2TrdTrack);
2141  }
2142 
2143  CbmTrdTrack* trCand;
2144  Int_t firstHitSunk = 0;
2145  for (ivTrdTrackCand = vTempTrdTrackCand.begin();
2146  ivTrdTrackCand != vTempTrdTrackCand.end();
2147  ivTrdTrackCand++) { //Loop over track candidates; mark used hits
2148  trCand = (*ivTrdTrackCand);
2149  //cout << "Track no " << trackInd++
2150  // << " has " << (*ivTrdTrackCand)->GetNofTrdHits()
2151  // << " hits and chi2= " << (*ivTrdTrackCand)->GetChi2() << endl;
2152 
2153  // if(nrLoop == 2) removeUsedHits = false;
2154 
2155  if (removeUsedHits) { //RemoveHits
2156  Int_t noHitsA = trCand->GetNofHits();
2157 
2158  bTrueTrack = true;
2159  firstHitSunk = 0;
2160 
2161  for (int i = 0; i < noHitsA; i++) {
2162  iHit = trCand->GetHitIndex(i);
2163  iglobalSetUsedHits = globalSetUsedHits.find(iHit);
2164 
2165  if (iglobalSetUsedHits != globalSetUsedHits.end()) {
2166  if (firstHitSunk == noHitsAccepted) {
2167  bTrueTrack = false;
2168  iFakeTrack++;
2169  break;
2170  }
2171  firstHitSunk++;
2172  }
2173  }
2174 
2175  // Int_t k = 0;
2176  // for(int i = 0; i < noHitsA; i++) {
2177  // if(i != 11) k = SA[i]+SA[i+1];
2178  // // if(i != 10) k = SA[i]+SA[i+1]+SA[i+2];
2179  // if(k == noHitsAccepted){
2180  // bTrueTrack = false;
2181  // iFakeTrack++;
2182  // break;
2183  // }
2184  // }
2185 
2186  //if(firstHitSunk != 0) cout << firstHitSunk << endl;
2187  if (bTrueTrack) {
2188  new ((*fArrayTrdTrack)[trackArrayInd]) CbmTrdTrack(*trCand);
2189  trackArrayInd++;
2190  for (int i = 0; i < noHitsA; i++) {
2191  iHit = trCand->GetHitIndex(i);
2192  globalSetUsedHits.insert(iHit);
2193  }
2194  delete trCand;
2195  }
2196 
2197  } else { // end of RemoveHits
2198  new ((*fArrayTrdTrack)[trackArrayInd]) CbmTrdTrack(*trCand);
2199  trackArrayInd++;
2200  delete trCand;
2201  }
2202  }
2203 
2204  cout << "\n--- Finding tracks finished" << endl;
2205  cout << ":::Track candidates: " << vTempTrdTrackCand.size() << endl
2206  << ":::Fake tracks: " << iFakeTrack << endl
2207  << ":::fArrayTrdTrack: " << fArrayTrdTrack->GetEntriesFast() << endl;
2208 
2209  vTempTrdTrackCand.clear();
2210  // delete trParam;
2211 
2212  if (!isEmpty) {
2213  licz++;
2214  isEmpty = true;
2215  }
2216 }
2217 // -----------------------------------------------------------------------
2218 
2219 
2220 // -----------------------------------------------------------------------
2221 void CbmL1CATrdTrackFinderSA::CreateSegments(
2222  vector<CbmL1TrdTracklet*> clSpacePointsAB,
2223  vector<CbmL1TrdTracklet*> clSpacePointsCD,
2224  vector<CbmL1TrdTracklet4*>& clTrackletsAD,
2225  Double_t dX,
2226  Double_t dY) {
2227  //--- method to create 4-hit tracklets from Space Points (smaller 2-hit tracklets) ---
2228  vector<CbmL1TrdTracklet*>::iterator iSpacePtA;
2229  vector<CbmL1TrdTracklet*>::iterator iSpacePtB;
2230  vector<CbmL1TrdTracklet*>::iterator iSpacePtStart;
2231 
2232  Int_t iIndex = 0, //index
2233  indSpacePtA, //index of 1st miniTracklet
2234  indSpacePtB, //index of 2nd miniTracklet
2235  indSpacePtC, //index of 3rd miniTracklet
2236  indSpacePtD, //index of 4th miniTracklet
2237  noXSurv = 0, //number of tracklets that survived the X distance cut
2238  noYSurv = 0, //number of tracklets that survived the Y distance cut
2239  noXYSurv = 0, //number of tracklets that survived both X & Y distance cuts
2240  noAllPairs = 0, iSegAcc14 = 0;
2241 
2242  CbmL1TrdTracklet4* clTr;
2243 
2244  CbmL1TrdTracklet *trLeft, *trRight;
2245 
2246  iSpacePtStart = clSpacePointsCD.begin();
2247 
2248  //---------- BEGIN: Creation tracklet 12-34 ---------------------------------
2249  for (iSpacePtA = clSpacePointsAB.begin(); iSpacePtA != clSpacePointsAB.end();
2250  iSpacePtA++) {
2251  trLeft = *iSpacePtA;
2252  indSpacePtA = trLeft->GetIndLeft();
2253  indSpacePtB = trLeft->GetIndRight();
2254 
2255  // cout << "GetCoord1: " << trLeft->GetCoord1() << endl;
2256  //cout << "--------------------" << clSpacePointsCD.size() << endl;
2257 
2258 
2259  //--- extrapolate --------------------------------------------
2260  Double_t y1 = trLeft->GetCoord(0), x1 = trLeft->GetCoord(1),
2261  y1z = trLeft->GetCoord(2), x1z = trLeft->GetCoord(3), x2 = 0,
2262  y2 = 0, distBetweenLayer = x1z - y1z,
2263  y2z = y1z + distBetweenLayer * 2, x2z = x1z + distBetweenLayer * 2;
2264 
2265  y2 = (y1 / y1z) * y2z;
2266  x2 = (x1 / x1z) * x2z;
2267  //--- end: extrapolate --------------------------------------
2268 
2269  Bool_t bFirst = true;
2270 
2271  //for(y = 0; y < trsize; y++){//y
2272  //trRight = clSpacePointsCD[y];
2273  //cout << trRight->GetCoord(0) << endl;
2274  // Bool_t bFirst = true;
2275  for (iSpacePtB = iSpacePtStart; iSpacePtB != clSpacePointsCD.end();
2276  iSpacePtB++) {
2277  trRight = *iSpacePtB;
2278 
2279  indSpacePtC = trRight->GetIndLeft();
2280  indSpacePtD = trRight->GetIndRight();
2281  noAllPairs++;
2282 
2283  Double_t yB = trRight->GetCoord(0), xB = trRight->GetCoord(1),
2284  yBz = trRight->GetCoord(2), xBz = trRight->GetCoord(3);
2285 
2286  /*
2287  //--- for estimation dx and dy ----------------------------
2288  trdHit1 = L1_DYNAMIC_CAST<CbmTrdHit*>( fArrayTrdHit->At(indSpacePtA);
2289  Int_t refInd1 = trdHit1->GetRefIndex();
2290  trdHit1 = L1_DYNAMIC_CAST<CbmTrdHit*>( fArrayTrdHit->At(indSpacePtB);
2291  Int_t refInd2 = trdHit1->GetRefIndex();
2292 
2293  trdHit1 = L1_DYNAMIC_CAST<CbmTrdHit*>( fArrayTrdHit->At(indSpacePtC);
2294  Int_t refInd3 = trdHit1->GetRefIndex();
2295  trdHit1 = L1_DYNAMIC_CAST<CbmTrdHit*>( fArrayTrdHit->At(indSpacePtD);
2296  Int_t refInd4 = trdHit1->GetRefIndex();
2297 
2298  CbmTrdPoint
2299  *pPointAy,
2300  *pPointBx,
2301  *pPointCy,
2302  *pPointDx;
2303 
2304  pPointAy = (CbmTrdPoint*) fMCPointArray->At(refInd1);
2305  pPointBx = (CbmTrdPoint*) fMCPointArray->At(refInd2);
2306  pPointCy = (CbmTrdPoint*) fMCPointArray->At(refInd3);
2307  pPointDx = (CbmTrdPoint*) fMCPointArray->At(refInd4);
2308 
2309  Int_t trIDAy = pPointAy->GetTrackID();
2310  Int_t trIDBx = pPointBx->GetTrackID();
2311  Int_t trIDCy = pPointCy->GetTrackID();
2312  Int_t trIDDx = pPointDx->GetTrackID();
2313 
2314  CbmMCTrack *pMCtrack;
2315  Int_t noTRDPoints;
2316 
2317  Double_t
2318  x1 = pPointBx->GetX(),
2319  y1 = pPointAy->GetY(),
2320  z1y = pPointAy->GetZ(),
2321  z1x = pPointBx->GetZ(),
2322  x2 = 0,
2323  y2 = 0,
2324  z2y = z1y+distBetweenLayer*2,
2325  z2x = z1x+distBetweenLayer*2;
2326 
2327  y2 = (y1/z1y)*z2y;
2328  x2 = (x1/z1x)*z2x;
2329 
2330  if(trIDAy == trIDCy){
2331  pMCtrack = (CbmMCTrack*) fMCTrackArray->At(trIDAy);
2332  noTRDPoints = pMCtrack->GetTrdPoints();
2333 
2334  if(noTRDPoints == 12){
2335  Double_t fPosY = pPointCy->GetY();
2336  fDistY->Fill(fabs(fPosY-y2));
2337  }
2338  }
2339 
2340  if(trIDBx == trIDDx){
2341  pMCtrack = (CbmMCTrack*) fMCTrackArray->At(trIDBx);
2342  noTRDPoints = pMCtrack->GetTrdPoints();
2343 
2344  if(noTRDPoints == 12){
2345  Double_t fPosX = pPointDx->GetX();
2346  fDistX->Fill(fabs(fPosX-x2));
2347  }
2348  }
2349  //--- end: for estimation dx and dy -------------------------------
2350  */
2351 
2352  if (trRight->GetCoord(0) > y2 + dY) {
2353  continue;
2354  } else {
2355  if (bFirst) {
2356  bFirst = false;
2357  iSpacePtStart = iSpacePtB;
2358  }
2359  }
2360  if (trRight->GetCoord(0) < y2 - dY) break;
2361 
2362 
2363  if ((trRight->GetCoord(1) < x2 + dX)
2364  && (trRight->GetCoord(1) > x2 - dX)) {
2365 
2366  iSegAcc14++; //we have one more tracklet
2367  //we create a new tracklet object to add it into a vector
2368  clTr = new CbmL1TrdTracklet4();
2369 
2370  clTr->SetCoord(0, y1);
2371  clTr->SetCoord(1, x1);
2372  clTr->SetCoord(2, yB);
2373  clTr->SetCoord(3, xB);
2374 
2375  clTr->M[0] = y1z;
2376  clTr->M[1] = x1z;
2377  clTr->M[2] = yBz;
2378  clTr->M[3] = xBz;
2379 
2380  //zed is a z-axis value to which we extrapolate
2381  Double_t zed = 0;
2382 
2383  //cout << trdHit1->GetZ() <<" : "<< trRight->GetCoord(3) << endl;
2384  if (Int_t(trRight->GetCoord(3) + 1) == fTrd14_Z) {
2385  zed = fTrd21_Z;
2386  //cout <<"fTrd21_Z: "<< fTrd21_Z << endl;
2387  }
2388  //cout << trRight->GetCoord(3) <<"\t"<< fTrd14_Z << endl;
2389  if (Int_t(trRight->GetCoord(3) + 1) == fTrd24_Z) {
2390  zed = fTrd31_Z;
2391  //cout <<"fTrd31_Z: "<< fTrd31_Z << endl;
2392  }
2393  //cout << trRight->GetCoord(3) <<"\t"<< fTrd24_Z << endl;
2394 
2395  //if(trRight->GetPlanesID(1) == 4) zed = fTrd21_Z;
2396  //if(trRight->GetPlanesID(1) == 8) zed = fTrd31_Z;
2397  // if(trRight->GetPlanesID(1) == 4) zed = 970;
2398  //if(trRight->GetPlanesID(1) == 8) zed = fTrd31_Z;
2399 
2400  /* //Extrapolated by CbmTrdTrackFitterKF::Extrapolate insteas fo DistTwoTrackletsY
2401  //Extrapolate && fo DistTwoTrackletsY give the same results
2402  CbmTrdTrack *trdTrack = new CbmTrdTrack();
2403  CbmTrdHit *hit1;
2404 
2405  hit1 = L1_DYNAMIC_CAST<CbmTrdHit*>( fArrayTrdHit->At(indSpacePtA);
2406  trdTrack->AddHit(indSpacePtA, hit1);
2407  hit1 = L1_DYNAMIC_CAST<CbmTrdHit*>( fArrayTrdHit->At(indSpacePtB);
2408  trdTrack->AddHit(indSpacePtB, hit1);
2409  hit1 = L1_DYNAMIC_CAST<CbmTrdHit*>( fArrayTrdHit->At(indSpacePtC);
2410  trdTrack->AddHit(indSpacePtC, hit1);
2411  hit1 = L1_DYNAMIC_CAST<CbmTrdHit*>( fArrayTrdHit->At(indSpacePtD);
2412  trdTrack->AddHit(indSpacePtD, hit1);
2413 
2414  trdTrack->SortHits();
2415 
2416  FairTrackParam *trParam;
2417  trParam = new FairTrackParam();
2418  trParam->SetQp(0.);
2419  trdTrack->SetParamFirst(*trParam);
2420  trdTrack->SetParamLast(*trParam);
2421 
2422  trdTrackFitterKF->DoFit(trdTrack);
2423 
2424  //FairTrackParam *e_track = new FairTrackParam();
2425  trdTrackFitterKF->Extrapolate(trdTrack, zed, trParam);
2426 
2427  clTr->SetExt(0, trParam->GetY());
2428  clTr->SetExt(1, trParam->GetX());
2429 
2430  cout <<"(MY, KF) "
2431  << DistTwoTrackletsY(indSpacePtB, indSpacePtD, zed) <<", "
2432  << trParam->GetX() << endl;
2433  */
2434 
2435  //set the extrapolated coordinates - Y
2436  clTr->SetExt(0, DistTwoTrackletsY(indSpacePtA, indSpacePtC, zed));
2437  //set the extrapolated coordinates - X
2438  clTr->SetExt(
2439  1,
2440  DistTwoTrackletsX(indSpacePtB, indSpacePtD, zed + distBetweenLayer));
2441  //set the indices of all 4 point in tracklet
2442  clTr->SetInd(0, indSpacePtA);
2443  clTr->SetInd(1, indSpacePtB);
2444  clTr->SetInd(2, indSpacePtC);
2445  clTr->SetInd(3, indSpacePtD);
2446  //initial value is always 0; it is modified if necessary
2447  clTr->SetVal(0);
2448 
2449  //position of the tracklet in tracklets vector
2450  clTr->SetIndex(iIndex);
2451  iIndex++;
2452 
2453  //add a tracklet to a vector
2454  clTrackletsAD.push_back(clTr);
2455  noXYSurv++;
2456  }
2457  }
2458  }
2459  if (fVerbose > 2) {
2460  cout << " ----------- Tracklet 12-34 ------------------" << endl;
2461  cout << " Number of X survivors: " << noXSurv << endl;
2462  cout << " Number of Y survivors: " << noYSurv << endl;
2463  cout << "Number of XY survivors: " << noXYSurv << endl;
2464  cout << " Number of all pairs: " << noAllPairs << endl;
2465  }
2466 }
2467 // -----------------------------------------------------------------------
2468 
2469 void CbmL1CATrdTrackFinderSA::CreateSpacePoints(
2470  vector<LayerWithHits> vTrdHitArrayA,
2471  vector<LayerWithHits> vTrdHitArrayB,
2472  vector<CbmL1TrdTracklet*>& clSpacePointsAB,
2473  Double_t sigmaA,
2474  Double_t sigmaB) {
2475  //method to create Space Points (SP) from two individual detector hist
2476 
2477  vector<LayerWithHits>::iterator itA;
2478  vector<LayerWithHits>::iterator itB;
2479  vector<LayerWithHits>::iterator itStart;
2480 
2481  Int_t noOverlapsAB = 0,
2482  // noMCOverlapsAB = 0,
2483  iInd = 0;
2484 
2485  CbmL1TrdTracklet* clSpacePt;
2486 
2487  Int_t indA, indB;
2488 
2489  Double_t A_X, A_Y, A_Z, A_DX, A_DY;
2490 
2491  Double_t B_X, B_Y, B_Z, B_DX, B_DY;
2492 
2493  Int_t A_planeID, B_planeID;
2494 
2495  Int_t A_mcTrID, B_mcTrID;
2496 
2497  // Double_t SPlength;
2498 
2499 
2500  /*
2501  68,0%; <= 1sigma
2502  95,5%; <= 2sigma
2503  99,7%, <= 3sigma
2504  */
2505 
2506  itStart = vTrdHitArrayB.begin();
2507 
2508  for (itA = vTrdHitArrayA.begin(); itA != vTrdHitArrayA.end(); itA++) {
2509 
2510  indA = (*itA).hitIndex;
2511  A_X = (*itA).X;
2512  A_Y = (*itA).Y;
2513  A_Z = (*itA).Z;
2514  A_DX = (*itA).DX;
2515  A_DY = (*itA).DY;
2516  A_planeID = (*itA).planeID;
2517  A_mcTrID = (*itA).mcTrackID;
2518 
2519  /*//--- extrapolate --------------------------------------------
2520  Double_t distBetweenLayer = fTrd14_Z-fTrd13_Z;
2521  cout << "distBetweenLayer: " << distBetweenLayer << endl;
2522 
2523  Double_t
2524  x1 = A_X,
2525  y1 = A_Y,
2526  z1y = A_Z,
2527  z1x = A_Z,
2528  x2 = 0,
2529  y2 = 0,
2530  z2y = z1y+distBetweenLayer, // 12 is a distance beetwen plane in TRD
2531  z2x = z1x+distBetweenLayer; // 12 is a distance beetwen plane in TRD
2532 
2533  y2 = (y1/z1y)*z2y;
2534  x2 = (x1/z1x)*z2x;
2535  //--- end: extrapolate --------------------------------------*/
2536 
2537  Bool_t bFirst = true;
2538  for (itB = itStart; itB != vTrdHitArrayB.end(); itB++) {
2539 
2540  indB = (*itB).hitIndex;
2541  B_X = (*itB).X;
2542  B_Y = (*itB).Y;
2543  B_Z = (*itB).Z;
2544  B_DX = (*itB).DX;
2545  B_DY = (*itB).DY;
2546  B_planeID = (*itB).planeID;
2547  B_mcTrID = (*itB).mcTrackID;
2548 
2549  if (B_Y + sigmaB * B_DY < A_Y - sigmaA * A_DY) {
2550  continue;
2551  } else {
2552  /* Finding the bottom level to begin the second loop.
2553  The value of bottom level is taken from previous "second loop".*/
2554  if (bFirst) {
2555  itStart = itB;
2556  bFirst = false;
2557  }
2558  }
2559 
2560  if (B_Y - sigmaB * B_DY > A_Y + sigmaA * A_DY) break;
2561 
2562  // cout <<"(A_X, A_DX, sigmaA*A_DX) "<< A_X <<", "<< A_DX <<", "<< sigmaA*A_DX <<", "<< endl;
2563  // cout <<"(A_Y, A_DY, sigmaA*A_DY) "<< A_Y <<", "<< A_DY <<", "<< sigmaA*A_DY <<", "<< endl;
2564  // cout <<"(B_X, B_DX, sigmaB*B_DX) "<< B_X <<", "<< B_DX <<", "<< sigmaB*B_DX <<", "<< endl;
2565  // cout <<"(B_Y, B_DY, sigmaB*B_DY) "<< B_Y <<", "<< B_DY <<", "<< sigmaB*B_DY <<", "<< endl;
2566  // cout <<"-----------------------------------------"<< endl;
2567  // if (B_Y < y2 - sigmaB*A_DY) continue;
2568  // if (B_Y > y2 + sigmaB*A_DY) break;
2569  // cout <<"[Extrapol, Real, Diff] "
2570  // << y2 <<", "<< B_Y <<", "<< B_Y-y2 << endl;
2571 
2572  if ((B_X - sigmaB * B_DX < A_X + sigmaA * A_DX)
2573  && (B_X + sigmaB * B_DX > A_X - sigmaA * A_DX)) {
2574 
2575  noOverlapsAB++;
2576  clSpacePt = new CbmL1TrdTracklet();
2577 
2578  clSpacePt->SetIndLeft(indA);
2579  clSpacePt->SetIndRight(indB);
2580  clSpacePt->SetVal(0);
2581  clSpacePt->SetIndex(iInd);
2582  clSpacePt->SetCoord(0, A_Y);
2583  clSpacePt->SetCoord(1, B_X);
2584  clSpacePt->SetCoord(2, A_Z);
2585  clSpacePt->SetCoord(3, B_Z);
2586  clSpacePt->SetPlanesID(A_planeID, B_planeID);
2587 
2588  // //----------------------------------------
2589  // TVector3 t1, t2, t3;
2590  // t1.SetXYZ(A_X, A_Y, A_Z);
2591  // t2.SetXYZ(B_X, B_Y, B_Z);
2592  // t3 = t2 - t1;
2593  // SPlength = t3.Mag();
2594 
2595  // Double_t
2596  // x = 0,
2597  // y = 0,
2598  // z = 0,
2599  // t = 0;
2600 
2601  // t = (z - A_Z)/(B_Z - A_Z);
2602  // x = A_X + t*(B_X-A_X);
2603  // y = A_Y + t*(B_Y-A_Y);
2604 
2605  // //SPlength = sqrt(pow(B_Y-A_Y,2)+pow(B_Y-A_Y,2));
2606  // if(Station == 1){
2607  // if(A_mcTrID == B_mcTrID){
2608  // noMCOverlapsAB++;
2609  // fSPlengthMC->Fill(SPlength);
2610  // fYat0MC->Fill(x);
2611  // }else{
2612  // fSPlength->Fill(SPlength);
2613  // fYat0->Fill(x);
2614  // }
2615  // }
2616  iInd++;
2617 
2618  clSpacePointsAB.push_back(clSpacePt);
2619 
2620  if (A_mcTrID == B_mcTrID) {
2621  fh_SP_xDiff_MC->Fill(A_X - B_X);
2622  fh_SP_yDiff_MC->Fill(A_Y - B_Y);
2623  } else {
2624  fh_SP_xDiff_nMC->Fill(A_X - B_X);
2625  fh_SP_yDiff_nMC->Fill(A_Y - B_Y);
2626  }
2627  }
2628  }
2629  }
2630  iInd = 0;
2631  if (fVerbose > 1) cout << "--- No Space Points: " << noOverlapsAB << endl;
2632  // cout <<"--- (No SP, MCSP) "
2633  // << noOverlapsAB <<", "
2634  // << noMCOverlapsAB <<endl;
2635 }
2636 
2637 
2638 // -----------------------------------------------------------------------
2639 void CbmL1CATrdTrackFinderSA::CreateAndManageSegments(
2640  vector<CbmL1TrdTracklet4*> clTracklets14,
2641  vector<CbmL1TrdTracklet4*> clTracklets58,
2642  vector<CbmL1TrdTracklet4*> clTracklets912) {
2643  //CreateAndManageSegments with no long track creation, 1 track = 1 segment
2644 
2645  vector<CbmL1TrdTracklet4*>::iterator itclSeg;
2646 
2647  CbmTrdTrack* tr1 = NULL;
2648  // CbmTrdHit* hit1 = NULL;
2649  // CbmTrdHit* hit2 = NULL;
2650  // CbmTrdHit* hit3 = NULL;
2651  // CbmTrdHit* hit4 = NULL;
2652 
2653  CbmL1TrdTracklet4* clTrdSeg;
2654 
2655  Int_t trackArrayInd = 0;
2656 
2657  Int_t indTrdHit1, indTrdHit2, indTrdHit3, indTrdHit4;
2658 
2659  for (itclSeg = clTracklets14.begin(); itclSeg != clTracklets14.end();
2660  itclSeg++) {
2661 
2662  clTrdSeg = *itclSeg;
2663  tr1 = new CbmTrdTrack();
2664 
2665  indTrdHit1 = clTrdSeg->GetInd(0);
2666  // hit1 = L1_DYNAMIC_CAST<CbmTrdHit*>( fArrayTrdHit->At(indTrdHit1) );
2667  tr1->AddHit(indTrdHit1, kTRDHIT);
2668 
2669  indTrdHit2 = clTrdSeg->GetInd(1);
2670  // hit2 = L1_DYNAMIC_CAST<CbmTrdHit*>( fArrayTrdHit->At(indTrdHit2) );
2671  tr1->AddHit(indTrdHit2, kTRDHIT);
2672 
2673  indTrdHit3 = clTrdSeg->GetInd(2);
2674  // hit3 = L1_DYNAMIC_CAST<CbmTrdHit*>( fArrayTrdHit->At(indTrdHit3) );
2675  tr1->AddHit(indTrdHit3, kTRDHIT);
2676 
2677  indTrdHit4 = clTrdSeg->GetInd(3);
2678  // hit4 = L1_DYNAMIC_CAST<CbmTrdHit*>( fArrayTrdHit->At(indTrdHit4) );
2679  tr1->AddHit(indTrdHit4, kTRDHIT);
2680 
2681  // tr1->SortHits();
2682 
2683  new ((*fArrayTrdTrack)[trackArrayInd]) CbmTrdTrack(*tr1);
2684 
2685  delete tr1;
2686  trackArrayInd++;
2687  }
2688 
2689  for (itclSeg = clTracklets58.begin(); itclSeg != clTracklets58.end();
2690  itclSeg++) {
2691 
2692  clTrdSeg = *itclSeg;
2693  tr1 = new CbmTrdTrack();
2694 
2695  indTrdHit1 = clTrdSeg->GetInd(0);
2696  // hit1 = L1_DYNAMIC_CAST<CbmTrdHit*>( fArrayTrdHit->At(indTrdHit1) );
2697  tr1->AddHit(indTrdHit1, kTRDHIT);
2698 
2699  indTrdHit2 = clTrdSeg->GetInd(1);
2700  // hit2 = L1_DYNAMIC_CAST<CbmTrdHit*>( fArrayTrdHit->At(indTrdHit2) );
2701  tr1->AddHit(indTrdHit2, kTRDHIT);
2702 
2703  indTrdHit3 = clTrdSeg->GetInd(2);
2704  // hit3 = L1_DYNAMIC_CAST<CbmTrdHit*>( fArrayTrdHit->At(indTrdHit3) );
2705  tr1->AddHit(indTrdHit3, kTRDHIT);
2706 
2707  indTrdHit4 = clTrdSeg->GetInd(3);
2708  // hit4 = L1_DYNAMIC_CAST<CbmTrdHit*>( fArrayTrdHit->At(indTrdHit4) );
2709  tr1->AddHit(indTrdHit4, kTRDHIT);
2710 
2711  // tr1->SortHits();
2712 
2713  new ((*fArrayTrdTrack)[trackArrayInd]) CbmTrdTrack(*tr1);
2714 
2715  delete tr1;
2716  trackArrayInd++;
2717  }
2718 
2719  for (itclSeg = clTracklets912.begin(); itclSeg != clTracklets912.end();
2720  itclSeg++) {
2721 
2722  clTrdSeg = *itclSeg;
2723  tr1 = new CbmTrdTrack();
2724 
2725  indTrdHit1 = clTrdSeg->GetInd(0);
2726  // hit1 = L1_DYNAMIC_CAST<CbmTrdHit*>( fArrayTrdHit->At(indTrdHit1) );
2727  tr1->AddHit(indTrdHit1, kTRDHIT);
2728 
2729  indTrdHit2 = clTrdSeg->GetInd(1);
2730  // hit2 = L1_DYNAMIC_CAST<CbmTrdHit*>( fArrayTrdHit->At(indTrdHit2) );
2731  tr1->AddHit(indTrdHit2, kTRDHIT);
2732 
2733  indTrdHit3 = clTrdSeg->GetInd(2);
2734  // hit3 = L1_DYNAMIC_CAST<CbmTrdHit*>( fArrayTrdHit->At(indTrdHit3) );
2735  tr1->AddHit(indTrdHit3, kTRDHIT);
2736 
2737  indTrdHit4 = clTrdSeg->GetInd(3);
2738  // hit4 = L1_DYNAMIC_CAST<CbmTrdHit*>( fArrayTrdHit->At(indTrdHit4) );
2739  tr1->AddHit(indTrdHit4, kTRDHIT);
2740 
2741  // tr1->SortHits();
2742 
2743  new ((*fArrayTrdTrack)[trackArrayInd]) CbmTrdTrack(*tr1);
2744 
2745  delete tr1;
2746  trackArrayInd++;
2747  }
2748 }
2749 // -----------------------------------------------------------------------
2750 
2751 
2752 // -----------------------------------------------------------------------
2753 Double_t CbmL1CATrdTrackFinderSA::Fit(
2754  CbmTrdTrack*
2755  tr) { //fits a track with a straight line and caluculates each point's deviation
2756 
2757  Double_t chi2 = 0;
2758 
2759  Int_t iHit;
2760  CbmTrdHit* pHit[12];
2761 
2762  Double_t y1, y2, z1, z2, yA[12], zA[12], yB[12], zB[12];
2763 
2764  Double_t dist1 = 0;
2765  Double_t dist2 = 0;
2766  // Double_t sum1 = 0;
2767  // Double_t sum2 = 0;
2768  Double_t yS1 = 0;
2769  Double_t yS2 = 0;
2770 
2771  Int_t noHits = tr->GetNofHits();
2772 
2773  for (int i = 0; i < noHits; i++) {
2774  iHit = tr->GetHitIndex(i);
2775 
2776  pHit[i] = L1_DYNAMIC_CAST<CbmTrdHit*>(fArrayTrdHit->At(iHit));
2777  }
2778 
2779  //first and last X(orY) pointa in the track
2780  yA[0] = pHit[0]->GetY();
2781  zA[0] = pHit[0]->GetZ();
2782  yB[10] = pHit[10]->GetY();
2783  zB[10] = pHit[10]->GetZ();
2784 
2785  //first and last Y(orX) pointa in the track
2786  yA[1] = pHit[1]->GetX();
2787  zA[1] = pHit[1]->GetZ();
2788  yB[11] = pHit[11]->GetX();
2789  zB[11] = pHit[11]->GetZ();
2790 
2791  Int_t t1 = 0;
2792  Int_t t2 = 1;
2793 
2794  for (int i = 0; i < 4; i++) {
2795 
2796  //for X(orY) coordinates
2797  t1 += 2;
2798  y1 = fabs(pHit[t1]->GetY());
2799  z1 = pHit[t1]->GetZ();
2800 
2801  yS1 = fabs((((yB[10] - yA[0]) * (z1 - zA[0])) / (zB[10] - zA[0])) + yA[0]);
2802  dist1 = fabs(yS1 - y1);
2803 
2804  //for Y(orX) coordinates
2805  t2 += 2;
2806  y2 = fabs(pHit[t2]->GetX());
2807  z2 = pHit[t2]->GetZ();
2808 
2809  yS2 = fabs((((yB[11] - yA[1]) * (z2 - zA[1])) / (zB[11] - zA[1])) + yA[1]);
2810  dist2 = fabs(yS2 - y2);
2811 
2812  chi2 += (dist1 + dist2) / 2;
2813  }
2814 
2815  // cout << "chi2: " << chi2 << endl;
2816 
2817 
2818  // for(int i = 0; i < 4; i++){
2819 
2820  // //for X(orY) coordinates
2821  // t1+=2;
2822  // y1 = fabs(pHit[t1]->Y());
2823  // z1 = pHit[t1]->Z();
2824 
2825  // yS1 = fabs((((yB[10]-yA[0])*(z1-zA[0]))/(zB[10]-zA[0]))+yA[0]);
2826  // dist1 = fabs(yS1-y1);
2827  // sum1+=dist1;
2828 
2829  // //for Y(orX) coordinates
2830  // t2+=2;
2831  // y2 = fabs(pHit[t2]->X());
2832  // z2 = pHit[t2]->Z();
2833 
2834  // yS2 = fabs((((yB[11]-yA[1])*(z2-zA[1]))/(zB[11]-zA[1]))+yA[1]);
2835  // dist2 = fabs(yS2-y2);
2836  // sum2+=dist2;
2837  // }
2838  // chi2=sum1+sum2;
2839 
2840 
2841  // Calculate length of the track ---------------
2842  // Double_t
2843  // length0 = 0,
2844  // length1 = 0,
2845  // length2 = 0;
2846 
2847  // t1 = 0;
2848  // t2 = 1;
2849 
2850  // Double_t
2851  // x = 0,
2852  // y = 0,
2853  // z = 0;
2854 
2855  // for(int i = 0; i < 4; i++){
2856  // y = fabs(pHit[t1+2]->Y())-fabs(pHit[t1]->Y());
2857  // z = (pHit[t1+2]->Z())-(pHit[t1]->Z());
2858  // length1 += sqrt(pow(y,2) + pow(z,2));
2859  // t1+=2;
2860 
2861  // x = fabs(pHit[t2+2]->X())-fabs(pHit[t2]->X());
2862  // z = (pHit[t2+2]->Z())-(pHit[t2]->Z());
2863  // length2 += sqrt(pow(x,2) + pow(z,2));
2864  // t2+=2;
2865  // }
2866 
2867  // width0 = width1 + width2;
2868  // // cout << width0 << endl;
2869 
2870  return chi2;
2871  //return length0;
2872 }
2873 // -----------------------------------------------------------------------
2874 
2875 Double_t CbmL1CATrdTrackFinderSA::Fit(
2876  Int_t M
2877  []) { //fits a track with a straight line and caluculates each point's deviation
2878 
2879  Double_t chi2 = 0;
2880 
2881  Int_t iHit;
2882  CbmTrdHit* pHit[12];
2883 
2884  Double_t y1, y2, z1, z2, yA[12], zA[12], yB[12], zB[12];
2885 
2886  Double_t dist1 = 0;
2887  Double_t dist2 = 0;
2888  Double_t yS1 = 0;
2889  Double_t yS2 = 0;
2890 
2891  Int_t noHits = 12;
2892 
2893  for (int i = 0; i < noHits; i++) {
2894  iHit = M[i];
2895  pHit[i] = L1_DYNAMIC_CAST<CbmTrdHit*>(fArrayTrdHit->At(iHit));
2896  }
2897 
2898  //first and last X(orY) pointa in the track
2899  yA[0] = pHit[0]->GetY();
2900  zA[0] = pHit[0]->GetZ();
2901  yB[10] = pHit[10]->GetY();
2902  zB[10] = pHit[10]->GetZ();
2903 
2904  //first and last Y(orX) pointa in the track
2905  yA[1] = pHit[1]->GetX();
2906  zA[1] = pHit[1]->GetZ();
2907  yB[11] = pHit[11]->GetX();
2908  zB[11] = pHit[11]->GetZ();
2909 
2910  Int_t t1 = 0;
2911  Int_t t2 = 1;
2912 
2913  for (int i = 0; i < 4; i++) {
2914 
2915  //for X(orY) coordinates
2916  t1 += 2;
2917  y1 = fabs(pHit[t1]->GetY());
2918  z1 = pHit[t1]->GetZ();
2919 
2920  yS1 = fabs((((yB[10] - yA[0]) * (z1 - zA[0])) / (zB[10] - zA[0])) + yA[0]);
2921  dist1 = fabs(yS1 - y1);
2922 
2923  //for Y(orX) coordinates
2924  t2 += 2;
2925  y2 = fabs(pHit[t2]->GetX());
2926  z2 = pHit[t2]->GetZ();
2927 
2928  yS2 = fabs((((yB[11] - yA[1]) * (z2 - zA[1])) / (zB[11] - zA[1])) + yA[1]);
2929  dist2 = fabs(yS2 - y2);
2930 
2931  //chi2 += (dist2+dist1)/2;
2932  chi2 += sqrt(pow(dist2, 2) + pow(dist1, 2));
2933  }
2934  return chi2;
2935 }
2936 // -----------------------------------------------------------------------
2937 
2938 // -----------------------------------------------------------------------
2939 Double_t CbmL1CATrdTrackFinderSA::FitLinear(Int_t M[], Int_t /*var*/ = 1) {
2940  //fit using a least square method
2941 
2942 
2943  Int_t noHits = 12;
2944 
2945  Int_t iHit;
2946  CbmTrdHit* hit;
2947 
2948  Double_t C1[12];
2949  Double_t C2[12];
2950  Double_t Z[12];
2951 
2952  Double_t n = 12;
2953 
2954  for (int i = 0; i < 12; i++) {
2955  C1[i] = 0;
2956  C2[i] = 0;
2957  Z[i] = 0;
2958  }
2959 
2960  // Int_t ind = 0;
2961 
2962  for (int i = 0; i < noHits; i++) {
2963 
2964  iHit = M[i];
2965  hit = L1_DYNAMIC_CAST<CbmTrdHit*>(fArrayTrdHit->At(iHit));
2966 
2967  C1[i] = hit->GetX();
2968  C2[i] = hit->GetY();
2969  Z[i] = hit->GetZ();
2970  }
2971 
2972  Double_t sumXZ = 0, sumYZ = 0, sumX = 0, sumY = 0, sumZx = 0, sumZy = 0,
2973  sumX2 = 0, sumY2 = 0, sumZx2 = 0, sumZy2 = 0;
2974 
2975  for (int i = 0; i < 12; i += 2) {
2976  // cout << "C1[" << i << "]=" << C1[i] << endl;
2977  // cout << "C2[" << i << "]=" << C2[i] << endl;
2978  // if(!(i % 2))
2979  {
2980  sumXZ += C1[i] * Z[i];
2981 
2982  sumX += C1[i];
2983  sumX2 += C1[i] * C1[i];
2984 
2985  sumZx += Z[i];
2986  sumZx2 += Z[i] * Z[i];
2987  }
2988  }
2989 
2990  for (int i = 1; i < 12; i += 2) {
2991  // cout << "C1[" << i << "]=" << C1[i] << endl;
2992  // cout << "C2[" << i << "]=" << C2[i] << endl;
2993  // if(!((1i % 2))
2994  {
2995  sumYZ += C2[i] * Z[i];
2996 
2997  sumY += C2[i];
2998  sumY2 += C2[i] * C2[i];
2999 
3000  sumZy += Z[i];
3001  sumZy2 += Z[i] * Z[i];
3002  }
3003  }
3004 
3005 
3006  cout << "";
3007  // Double_t a = (n*sumXZ - sumX*sumZ)/(n*sumX2 - sumX*sumX);
3008  // Double_t b = (sumZ - a*sumX)/n;
3009  Double_t r1 =
3010  (n * sumXZ - sumX * sumZx)
3011  / sqrt((n * sumX2 - sumX * sumX) * (n * sumZx2 - sumZx * sumZx));
3012  // Double_t Sa = sqrt((n*(sumY2 -a*sumXZ-b*sumY))/((n-2)*(n*sumX2-sumX*sumX)));
3013  // Double_t Sb = Sa*sqrt(sumX2/n);
3014 
3015 
3016  Double_t r2 =
3017  (n * sumYZ - sumY * sumZy)
3018  / sqrt((n * sumY2 - sumY * sumY) * (n * sumZy2 - sumZy * sumZy));
3019 
3020 
3021  // cout << "Linear coefficients: a= "
3022  // << a << ", Sa= " << Sa << ", b= " << b << ", Sb= " << Sb << endl
3023  // << ", r=" << r << ", chi2=" << tr->GetChi2() << endl;
3024 
3025 
3026  return sqrt(r1 * r1 + r2 * r2);
3027 }
3028 
3029 Double_t CbmL1CATrdTrackFinderSA::FitLSM(Int_t M[]) {
3030 
3031  Double_t chi2 = 0;
3032  Int_t iHit;
3033  CbmTrdHit* pHit[12];
3034 
3035  Int_t noHits = 12;
3036  for (int i = 0; i < noHits; i++) {
3037  iHit = M[i];
3038  pHit[i] = L1_DYNAMIC_CAST<CbmTrdHit*>(fArrayTrdHit->At(iHit));
3039  }
3040 
3041  Double_t xAv = 0, yAv = 0, zAvx = 0, zAvy = 0, sumXiXav = 0, sumYiYav = 0,
3042  sumXiXav_ZiZAv = 0, sumYiYav_ZiZAv = 0, sumZiZav_x = 0,
3043  sumZiZav_y = 0, bY = 0, bX = 0;
3044 
3045  for (int i = 0; i < 12; i += 2) {
3046  yAv += pHit[i]->GetY();
3047  zAvy += pHit[i]->GetZ();
3048  }
3049  for (int i = 1; i < 12; i += 2) {
3050  xAv += pHit[i]->GetX();
3051  zAvx += pHit[i]->GetZ();
3052  }
3053 
3054  for (int i = 0; i < 12; i += 2) {
3055  sumYiYav_ZiZAv += ((pHit[i]->GetY()) - yAv) * ((pHit[i]->GetZ()) - zAvy);
3056  sumYiYav += pow((pHit[i]->GetY()) - yAv, 2);
3057  sumZiZav_y += pow((pHit[i]->GetZ()) - zAvy, 2);
3058  }
3059  for (int i = 1; i < 12; i += 2) {
3060  sumXiXav_ZiZAv += ((pHit[i]->GetX()) - xAv) * ((pHit[i]->GetZ()) - zAvx);
3061  sumXiXav += pow((pHit[i]->GetX()) - xAv, 2);
3062  sumZiZav_x += pow((pHit[i]->GetZ()) - zAvx, 2);
3063  }
3064 
3065  bY = sumYiYav_ZiZAv / sumYiYav;
3066  bX = sumXiXav_ZiZAv / sumXiXav;
3067 
3068  chi2 = sqrt((sumZiZav_y - bY * sumYiYav_ZiZAv) / (4))
3069  + sqrt((sumZiZav_x - bX * sumXiXav_ZiZAv) / (4));
3070 
3071  return chi2;
3072 }
3073 
3074 
3075 Double_t CbmL1CATrdTrackFinderSA::FitKF(CbmTrdTrack* pTrack) {
3076 
3077  // Declare variables outside the loop
3078  CbmTrdHit* pHit = NULL;
3079  CbmKFTrdHit* pKFHit = NULL;
3080  Int_t hitIndex = 0;
3081  // Int_t materialIndex = 0;
3082  Double_t eLoss = 0.;
3083 
3084  // Loop over TRD hits of the track
3085  Int_t nTracks = pTrack->GetNofHits();
3086  for (Int_t iHit = 0; iHit < nTracks; iHit++) {
3087  // Get current hit index
3088  hitIndex = pTrack->GetHitIndex(iHit);
3089  //Get the pointer to the CbmTrdHit
3090  pHit = L1_DYNAMIC_CAST<CbmTrdHit*>(fArrayTrdHit->At(hitIndex));
3091  // Accumulate TR energy loss
3092  eLoss += pHit->GetELoss();
3093  // Create CbmKFTrdHit
3094  pKFHit = new CbmKFTrdHit();
3095  pKFHit->Create(pHit);
3096  // materialIndex = pKFHit->MaterialIndex;
3097  // Add to the KFTrack
3098  fKfTrack->fHits.push_back(pKFHit);
3099  } // Loop over TRD hits
3100 
3101  fKfTrack->GetRefChi2() = 0.;
3102  fKfTrack->GetRefNDF() = 0;
3103  fKfTrack->SetTrackParam(
3104  *(const_cast<FairTrackParam*>(pTrack->GetParamLast())));
3105 
3106  fKfTrack->Fit();
3107  // Store parameters at first layer
3108  fKfTrack->GetTrackParam(
3109  *(const_cast<FairTrackParam*>(pTrack->GetParamFirst())));
3110  // Store parameters at last layer
3111  fKfTrack->GetTrackParam(
3112  *(const_cast<FairTrackParam*>(pTrack->GetParamLast())));
3113  // Store chi2 of fit
3114  pTrack->SetChiSq(fKfTrack->GetRefChi2());
3115  pTrack->SetNDF(fKfTrack->GetRefNDF());
3116 
3117  // Store accumulated TR energy loss
3118  pTrack->SetELoss(eLoss);
3119 
3120  // Delete CbmKFTrdHit objects
3121  vector<CbmKFHit*>::iterator it;
3122  for (it = fKfTrack->fHits.begin(); it != fKfTrack->fHits.end(); it++) {
3123  pKFHit = L1_DYNAMIC_CAST<CbmKFTrdHit*>(*it);
3124  delete pKFHit;
3125  }
3126  fKfTrack->fHits.clear();
3127 
3128  return 0;
3129 }
3130 
3131 Double_t CbmL1CATrdTrackFinderSA::FitTLinearFitter(Int_t M[]) {
3132 
3133  Double_t chi2 = 0;
3134  Double_t chi2x = 0;
3135  Double_t chi2y = 0;
3136  const Int_t nrPnts = 12;
3137  Int_t k = 0;
3138  Int_t w = 0;
3139  Int_t j = 6;
3140  Double_t ax[6];
3141  Double_t ay[6];
3142  Double_t az1[6];
3143  // Double_t az2[6];
3144 
3145  CbmTrdHit* trdHit;
3146  for (int i = 0; i < nrPnts; i++) {
3147  trdHit = L1_DYNAMIC_CAST<CbmTrdHit*>(fArrayTrdHit->At(M[i]));
3148  k = i + 1;
3149  if ((k % 2) == 0) {
3150  ax[w] = trdHit->GetX();
3151  az1[w] = trdHit->GetZ();
3152  w++;
3153  } else {
3154  ay[w] = trdHit->GetY();
3155  // az2[w] = trdHit->GetZ();
3156  }
3157  }
3158 
3159  TLinearFitter* lf;
3160  lf = new TLinearFitter(2);
3161  lf->SetFormula("hyp2");
3162  lf->StoreData(0);
3163 
3164  TVectorD x;
3165  x.Use(j, ax);
3166  TVectorD z1;
3167  z1.Use(j, az1);
3168  lf->AssignData(j, 1, ax, az1);
3169  lf->Eval();
3170  chi2x = lf->GetChisquare();
3171 
3172  delete lf;
3173  //--------------------------------
3174  lf = new TLinearFitter(2);
3175  lf->SetFormula("hyp2");
3176  lf->StoreData(0);
3177 
3178  TVectorD y;
3179  x.Use(j, ay);
3180  TVectorD z2;
3181  z2.Use(j, az1);
3182  lf->AssignData(j, 1, ay, az1);
3183  lf->Eval();
3184  chi2y = lf->GetChisquare();
3185 
3186  chi2 = chi2x + chi2y;
3187 
3188  return chi2;
3189 }
3190 
3191 Int_t CbmL1CATrdTrackFinderSA::PdgToCharge(Int_t pdgCode) {
3192  Int_t value;
3193  switch (pdgCode) { //switch
3194  case 11: //electron
3195  {
3196  value = -1;
3197  break;
3198  }
3199  case -11: //positron
3200  {
3201  value = +1;
3202  break;
3203  }
3204  case 13: //muon-
3205  {
3206  value = -1;
3207  break;
3208  }
3209  case -13: //muon+
3210  {
3211  value = +1;
3212  break;
3213  }
3214  case 22: {
3215  value = 0; //gamma
3216  break;
3217  }
3218  case 211: {
3219  value = +1; //pi+
3220  break;
3221  }
3222  case -211: {
3223  value = -1; //pi-
3224  break;
3225  }
3226  case 111: {
3227  value = 0; //pi0
3228  break;
3229  }
3230  case 213: {
3231  value = +1; //rho+
3232  break;
3233  }
3234  case -213: {
3235  value = -1; //rho-
3236  break;
3237  }
3238  case 113: {
3239  value = 0; //rho0
3240  break;
3241  }
3242  case 221: {
3243  value = 0; //eta
3244  break;
3245  }
3246  case 323: {
3247  value = 0; //omega
3248  break;
3249  }
3250  case 333: {
3251  value = 0; //phi
3252  break;
3253  }
3254  case 321: {
3255  value = +1; //kaon+
3256  break;
3257  }
3258  case -321: {
3259  value = -1; //kaon-
3260  break;
3261  }
3262 
3263  case 443: {
3264  value = 0; //jpsi
3265  break;
3266  }
3267  case 2112: {
3268  value = 0; //neutron
3269  break;
3270  }
3271  case 2212: {
3272  value = +1; //proton
3273  break;
3274  }
3275  case -2212: {
3276  value = 0; //anti-proton (?!)
3277  break;
3278  }
3279  case 3222: {
3280  value = +1; //sigma+
3281  break;
3282  }
3283  case 3112: {
3284  value = -1; //sigma-
3285  break;
3286  }
3287  case 3212: {
3288  value = 0; //sigma0
3289  break;
3290  }
3291  default: {
3292  value = -1000;
3293  break;
3294  }
3295  }
3296  return value;
3297 }
3298 
3299 Bool_t CbmL1CATrdTrackFinderSA::Rejection(Double_t Procent, Int_t num) {
3300 
3301  Bool_t accept = true;
3302  Int_t random;
3303 
3304  random = (rand() % num);
3305 
3306  if (random >= Procent) accept = false;
3307 
3308  return accept;
3309 }
CbmL1TrdTracklet::SetCoord
void SetCoord(Int_t ind, Double_t val)
Definition: CbmL1TrdTracklet.h:29
CbmL1CATrdTrackFinderSA::totSortSPs
Double_t totSortSPs
Definition: CbmL1CATrdTrackFinderSA.h:160
CbmL1CATrdTrackFinderSA::totSecondLoopTime
Double_t totSecondLoopTime
Definition: CbmL1CATrdTrackFinderSA.h:161
CbmL1CATrdTrackFinderSA::DeleteTracklets
void DeleteTracklets(std::vector< CbmL1TrdTracklet4 * > vect)
Definition: CbmL1CATrdTrackFinderSA.cxx:1138
CbmHit::GetZ
Double_t GetZ() const
Definition: CbmHit.h:70
CbmL1TrdTracklet::GetIndRight
Int_t GetIndRight()
Definition: CbmL1TrdTracklet.h:17
CbmL1CATrdTrackFinderSA::fDistLongX
TH1F * fDistLongX
Definition: CbmL1CATrdTrackFinderSA.h:250
CbmL1CATrdTrackFinderSA::fh_SP_yDiff_nMC
TH1F * fh_SP_yDiff_nMC
Definition: CbmL1CATrdTrackFinderSA.h:296
CbmL1CATrdTrackFinderSA::TrdPar
CbmGeoTrdPar * TrdPar
Definition: CbmL1CATrdTrackFinderSA.h:96
CbmL1CATrdTrackFinderSA::findNeighbour
TStopwatch findNeighbour
Definition: CbmL1CATrdTrackFinderSA.h:155
CbmL1CATrdTrackFinderSA::fMomDistShortPrimaryX
TH2F * fMomDistShortPrimaryX
Definition: CbmL1CATrdTrackFinderSA.h:272
CbmL1CATrdTrackFinderSA::fMomDistExtrapolExtraY
TH2F * fMomDistExtrapolExtraY
Definition: CbmL1CATrdTrackFinderSA.h:270
CbmVertex.h
CbmTrack::GetParamLast
const FairTrackParam * GetParamLast() const
Definition: CbmTrack.h:62
CbmL1TrdTracklet4::SetIndex
void SetIndex(Int_t index)
Definition: CbmL1TrdTracklet4.h:45
CbmL1CATrdTrackFinderSA::fDistShortBX
TH1F * fDistShortBX
Definition: CbmL1CATrdTrackFinderSA.h:257
CbmL1CATrdTrackFinderSA::fVerbose
Int_t fVerbose
Definition: CbmL1CATrdTrackFinderSA.h:110
CbmL1CATrdTrackFinderSA::totRefittingKF
Double_t totRefittingKF
Definition: CbmL1CATrdTrackFinderSA.h:162
CbmPixelHit::GetX
Double_t GetX() const
Definition: CbmPixelHit.h:83
CbmKF.h
CbmL1CATrdTrackFinderSA::fKfTrack
CbmKFTrack * fKfTrack
Definition: CbmL1CATrdTrackFinderSA.h:152
CbmL1CATrdTrackFinderSA::Fit
Double_t Fit(CbmTrdTrack *tr)
Definition: CbmL1CATrdTrackFinderSA.cxx:2753
CbmTrdTrackFitterKF::DoFit
Int_t DoFit(CbmTrdTrack *pTrack)
Definition: CbmTrdTrackFitterKF.cxx:83
CbmTrack::GetNofHits
virtual Int_t GetNofHits() const
Definition: CbmTrack.h:53
CbmL1CATrdTrackFinderSA::Layer::Z
Double_t Z[12]
Definition: CbmL1CATrdTrackFinderSA.h:61
CbmPixelHit::GetY
Double_t GetY() const
Definition: CbmPixelHit.h:84
CbmL1CATrdTrackFinderSA::fh_SP_yDiff_MC
TH1F * fh_SP_yDiff_MC
Definition: CbmL1CATrdTrackFinderSA.h:293
sqrt
friend F32vec4 sqrt(const F32vec4 &a)
Definition: L1/vectors/P4_F32vec4.h:41
CbmPixelHit::GetDx
Double_t GetDx() const
Definition: CbmPixelHit.h:85
CbmTrdHit::GetELoss
Double_t GetELoss() const
Definition: CbmTrdHit.h:79
CbmL1CATrdTrackFinderSA::totTagTracks
Double_t totTagTracks
Definition: CbmL1CATrdTrackFinderSA.h:160
CbmTrdTrackFitterKF
Definition: CbmTrdTrackFitterKF.h:15
CbmL1CATrdTrackFinderSA::WriteHistogramms
void WriteHistogramms()
Definition: CbmL1CATrdTrackFinderSA.cxx:1436
CbmL1CATrdTrackFinderSA::fDistShortY
TH1F * fDistShortY
Definition: CbmL1CATrdTrackFinderSA.h:253
CbmL1CATrdTrackFinderSA::fNoTrdPerStation
Int_t fNoTrdPerStation
Definition: CbmL1CATrdTrackFinderSA.h:69
CbmL1TrdTracklet4::vAccostLeft
std::vector< Int_t > vAccostLeft
Definition: CbmL1TrdTracklet4.h:56
CbmL1TrdTracklet.h
CbmKFTrdHit.h
rootMgr
static FairRootManager * rootMgr
Definition: CbmDeviceHitBuilderTof.cxx:72
CbmL1CATrdTrackFinderSA::fMomDistLongExtraY
TH2F * fMomDistLongExtraY
Definition: CbmL1CATrdTrackFinderSA.h:265
CbmL1CATrdTrackFinderSA::fNTrdHits
Int_t fNTrdHits
Definition: CbmL1CATrdTrackFinderSA.h:65
CbmL1CATrdTrackFinderSA::CbmL1CATrdTrackFinderSA
CbmL1CATrdTrackFinderSA()
Definition: CbmL1CATrdTrackFinderSA.cxx:60
CbmKFTrack::GetTrackParam
void GetTrackParam(FairTrackParam &track)
Definition: CbmKFTrack.cxx:45
CbmL1CATrdTrackFinderSA::fPlane5Ydens
TH1F * fPlane5Ydens
Definition: CbmL1CATrdTrackFinderSA.h:281
i
int i
Definition: L1/vectors/P4_F32vec4.h:25
CbmL1CATrdTrackFinderSA::sortHits
TStopwatch sortHits
Definition: CbmL1CATrdTrackFinderSA.h:156
CbmL1CATrdTrackFinderSA::DistTwoTrackletsX
Double_t DistTwoTrackletsX(Int_t iIndexFirst, Int_t iIndexSecond, Double_t zed)
Definition: CbmL1CATrdTrackFinderSA.cxx:1662
CbmTrdHit
data class for a reconstructed Energy-4D measurement in the TRD
Definition: CbmTrdHit.h:35
CbmL1TrdTracklet4::vAccostRight
std::vector< Int_t > vAccostRight
Definition: CbmL1TrdTracklet4.h:57
CbmL1TrdTracklet4::GetIndex
Int_t GetIndex()
Definition: CbmL1TrdTracklet4.h:30
CbmTrdTrackFitterKF.h
CbmL1CATrdTrackFinderSA::fMomDistExtrapolExtraX
TH2F * fMomDistExtrapolExtraX
Definition: CbmL1CATrdTrackFinderSA.h:269
CbmL1TrdTracklet4
Definition: CbmL1TrdTracklet4.h:7
CbmL1CATrdTrackFinderSA::fNofEvents
Int_t fNofEvents
Definition: CbmL1CATrdTrackFinderSA.h:100
CbmL1CATrdTrackFinderSA::Init
void Init()
Definition: CbmL1CATrdTrackFinderSA.cxx:187
CbmL1.h
CbmL1CATrdTrackFinderSA::fArrayTrdTrack
TClonesArray * fArrayTrdTrack
Definition: CbmL1CATrdTrackFinderSA.h:93
CbmHit::GetRefId
Int_t GetRefId() const
Definition: CbmHit.h:72
CbmL1CATrdTrackFinderSA::totSortHits
Double_t totSortHits
Definition: CbmL1CATrdTrackFinderSA.h:160
CbmL1CATrdTrackFinderSA::fMomDistShortExtraY
TH2F * fMomDistShortExtraY
Definition: CbmL1CATrdTrackFinderSA.h:275
CbmL1CATrdTrackFinderSA::CreateAndManageSegments
void CreateAndManageSegments(std::vector< CbmL1TrdTracklet4 * > clTracklets14, std::vector< CbmL1TrdTracklet4 * > clTracklets58, std::vector< CbmL1TrdTracklet4 * > clTracklets912)
Definition: CbmL1CATrdTrackFinderSA.cxx:2639
CbmL1CATrdTrackFinderSA::planeHits
struct CbmL1CATrdTrackFinderSA::LayerWithHits planeHits
CbmL1CATrdTrackFinderSA::fDistShortBY
TH1F * fDistShortBY
Definition: CbmL1CATrdTrackFinderSA.h:258
CbmL1CATrdTrackFinderSA::CreateTracks
void CreateTracks(std::vector< CbmL1TrdTracklet4 * > clTracklets14, std::vector< CbmL1TrdTracklet4 * > clTracklets58, std::vector< CbmL1TrdTracklet4 * > clTracklets912, std::set< Int_t > &setUsedHits, Bool_t removeUsedHits, Bool_t competition=true, Int_t nrLoop=0)
Definition: CbmL1CATrdTrackFinderSA.cxx:1895
CbmL1CATrdTrackFinderSA::fMomDistShortExtraX
TH2F * fMomDistShortExtraX
Definition: CbmL1CATrdTrackFinderSA.h:274
CbmL1CATrdTrackFinderSA::OverlapsHitsXY
Bool_t OverlapsHitsXY(Int_t posA, Int_t posB)
Definition: CbmL1CATrdTrackFinderSA.cxx:1782
CbmTrack::SetParamLast
void SetParamLast(const FairTrackParam *par)
Definition: CbmTrack.h:76
CbmTrack::SetNDF
void SetNDF(Int_t ndf)
Definition: CbmTrack.h:71
CbmL1TrdTracklet4::GetExt
Double_t GetExt(Int_t ind)
Definition: CbmL1TrdTracklet4.h:48
CbmL1CATrdTrackFinderSA::fNoEvTime
TH2F * fNoEvTime
Definition: CbmL1CATrdTrackFinderSA.h:290
CbmPixelHit::GetDy
Double_t GetDy() const
Definition: CbmPixelHit.h:86
CbmL1CATrdTrackFinderSA::DataBranches
void DataBranches()
Definition: CbmL1CATrdTrackFinderSA.cxx:1255
CbmL1CATrdTrackFinderSA::fYat0MC
TH1F * fYat0MC
Definition: CbmL1CATrdTrackFinderSA.h:288
CbmL1CATrdTrackFinderSA::fInstance
static CbmL1CATrdTrackFinderSA * fInstance
Definition: CbmL1CATrdTrackFinderSA.h:49
CbmL1CATrdTrackFinderSA::fMCPointArray
TClonesArray * fMCPointArray
Definition: CbmL1CATrdTrackFinderSA.h:104
CbmL1CATrdTrackFinderSA::doFind
TStopwatch doFind
Definition: CbmL1CATrdTrackFinderSA.h:156
CbmL1CATrdTrackFinderSA::fDistShortX
TH1F * fDistShortX
Definition: CbmL1CATrdTrackFinderSA.h:252
CbmL1CATrdTrackFinderSA::FitLSM
Double_t FitLSM(Int_t M[])
Definition: CbmL1CATrdTrackFinderSA.cxx:3029
CbmL1TrdTracklet4::SetExt
void SetExt(Int_t ind, Double_t val)
Definition: CbmL1TrdTracklet4.h:53
ClassImp
ClassImp(CbmL1CATrdTrackFinderSA)
CbmL1CATrdTrackFinderSA::tagTracks
TStopwatch tagTracks
Definition: CbmL1CATrdTrackFinderSA.h:156
CbmL1TrdTracklet::SetIndLeft
void SetIndLeft(Int_t indLeft)
Definition: CbmL1TrdTracklet.h:24
CbmL1Def.h
CbmL1CATrdTrackFinderSA::fMomDistShortPrimaryY
TH2F * fMomDistShortPrimaryY
Definition: CbmL1CATrdTrackFinderSA.h:273
CbmL1CATrdTrackFinderSA::refittingKF
TStopwatch refittingKF
Definition: CbmL1CATrdTrackFinderSA.h:157
CbmL1CATrdTrackFinderSA::CompareChi2TrdTrack
static Bool_t CompareChi2TrdTrack(const CbmTrdTrack *a, const CbmTrdTrack *b)
Definition: CbmL1CATrdTrackFinderSA.h:306
CbmL1CATrdTrackFinderSA::delTime
TStopwatch delTime
Definition: CbmL1CATrdTrackFinderSA.h:156
CbmKFTrack::GetRefChi2
Double_t & GetRefChi2()
array[15] of covariance matrix
Definition: CbmKFTrack.h:60
CbmL1CATrdTrackFinderSA.h
CbmL1CATrdTrackFinderSA::createSPs_SL
TStopwatch createSPs_SL
Definition: CbmL1CATrdTrackFinderSA.h:155
CbmL1CATrdTrackFinderSA::CreateSegments
void CreateSegments(std::vector< CbmL1TrdTracklet * > clSpacePointsAB, std::vector< CbmL1TrdTracklet * > clSpacePointsCD, std::vector< CbmL1TrdTracklet4 * > &clTrackletsAD, Double_t dX, Double_t dY)
Definition: CbmL1CATrdTrackFinderSA.cxx:2221
CbmL1CATrdTrackFinderSA::LayerWithHits::planeID
Int_t planeID
Definition: CbmL1CATrdTrackFinderSA.h:82
CbmL1CATrdTrackFinderSA::CreateSpacePoints
void CreateSpacePoints(std::vector< LayerWithHits > vTrdHitArrayA, std::vector< LayerWithHits > vTrdHitArrayB, std::vector< CbmL1TrdTracklet * > &clSpacePointsAB, Double_t sigmaA, Double_t sigmaB)
Definition: CbmL1CATrdTrackFinderSA.cxx:2469
CbmL1CATrdTrackFinderSA::LayerWithHits::DX
Double_t DX
Definition: CbmL1CATrdTrackFinderSA.h:80
CbmStsTrack.h
Data class for STS tracks.
CbmKFTrack::GetRefNDF
Int_t & GetRefNDF()
Chi^2 after fit.
Definition: CbmKFTrack.h:61
z2
Double_t z2[nSects2]
Definition: pipe_v16a_mvdsts100.h:11
CbmL1CATrdTrackFinderSA::CompareY
static Bool_t CompareY(LayerWithHits A, LayerWithHits B)
Definition: CbmL1CATrdTrackFinderSA.h:120
CbmL1CATrdTrackFinderSA::TrdLayout
void TrdLayout()
Definition: CbmL1CATrdTrackFinderSA.cxx:1275
CbmL1CATrdTrackFinderSA::LayerWithHits::X
Double_t X
Definition: CbmL1CATrdTrackFinderSA.h:77
CbmL1CATrdTrackFinderSA::LayerWithHits::mcTrackID
Int_t mcTrackID
Definition: CbmL1CATrdTrackFinderSA.h:76
CbmL1CATrdTrackFinderSA::DistTwoTrackletsY
Double_t DistTwoTrackletsY(Int_t iIndexFirst, Int_t iIndexSecond, Double_t zed)
Definition: CbmL1CATrdTrackFinderSA.cxx:1620
CbmL1CATrdTrackFinderSA::sortSPs
TStopwatch sortSPs
Definition: CbmL1CATrdTrackFinderSA.h:155
CbmL1CATrdTrackFinderSA::tempTrack
struct CbmL1CATrdTrackFinderSA::TempTrackStruct tempTrack
CbmVertex::Position
void Position(TVector3 &pos) const
Definition: CbmVertex.h:74
CbmL1CATrdTrackFinderSA::fDistY12
TH1F * fDistY12
Definition: CbmL1CATrdTrackFinderSA.h:260
CbmL1CATrdTrackFinderSA::fArrayTrdHit
TClonesArray * fArrayTrdHit
Definition: CbmL1CATrdTrackFinderSA.h:92
CbmL1CATrdTrackFinderSA::fUnUsedHitsPerPlane
TH1F * fUnUsedHitsPerPlane
Definition: CbmL1CATrdTrackFinderSA.h:298
CbmTrack::GetHitIndex
Int_t GetHitIndex(Int_t iHit) const
Definition: CbmTrack.h:54
CbmL1CATrdTrackFinderSA::fvTrdHitArr
std::vector< LayerWithHits > fvTrdHitArr[12]
Definition: CbmL1CATrdTrackFinderSA.h:86
CbmL1CATrdTrackFinderSA::DoFind
Int_t DoFind(TClonesArray *hitArray, TClonesArray *trackArray)
Definition: CbmL1CATrdTrackFinderSA.cxx:224
CbmL1CATrdTrackFinderSA::fMomDistExtrapolPrimaryY
TH2F * fMomDistExtrapolPrimaryY
Definition: CbmL1CATrdTrackFinderSA.h:268
CbmVertex
Definition: CbmVertex.h:26
CbmL1CATrdTrackFinderSA::fh_chi2hit_plane
TH2F * fh_chi2hit_plane
Definition: CbmL1CATrdTrackFinderSA.h:247
CbmL1TrdTracklet::GetCoord
Double_t GetCoord(Int_t ind)
Definition: CbmL1TrdTracklet.h:21
CbmL1CATrdTrackFinderSA::totCreateTracks
Double_t totCreateTracks
Definition: CbmL1CATrdTrackFinderSA.h:160
CbmL1CATrdTrackFinderSA::fTrd21_Z
Double_t fTrd21_Z
Definition: CbmL1CATrdTrackFinderSA.h:150
CbmL1CATrdTrackFinderSA::fNoTrdStations
Int_t fNoTrdStations
Definition: CbmL1CATrdTrackFinderSA.h:68
CbmL1CATrdTrackFinderSA::TempTrackStruct::M
Int_t M[12]
Definition: CbmL1CATrdTrackFinderSA.h:128
CbmKFTrdHit::Create
void Create(CbmTrdHit *hit)
Definition: CbmKFTrdHit.cxx:23
CbmTrack::SetChiSq
void SetChiSq(Double_t chiSq)
Definition: CbmTrack.h:70
CbmL1CATrdTrackFinderSA::totDelTime
Double_t totDelTime
Definition: CbmL1CATrdTrackFinderSA.h:161
CbmTrdHit.h
Class for hits in TRD detector.
CbmL1CATrdTrackFinderSA::createSPs
TStopwatch createSPs
Definition: CbmL1CATrdTrackFinderSA.h:155
CbmL1CATrdTrackFinderSA::geoLayer
struct CbmL1CATrdTrackFinderSA::Layer geoLayer
CbmKFTrack.h
CbmL1CATrdTrackFinderSA::fTrd24_Z
Double_t fTrd24_Z
Definition: CbmL1CATrdTrackFinderSA.h:150
CbmL1CATrdTrackFinderSA::fMCTrackArray
TClonesArray * fMCTrackArray
Definition: CbmL1CATrdTrackFinderSA.h:103
CbmL1CATrdTrackFinderSA::fh_chi2hit
TH1F * fh_chi2hit
Definition: CbmL1CATrdTrackFinderSA.h:243
CbmTrdTrackFitterKF::Init
void Init()
Definition: CbmTrdTrackFitterKF.cxx:60
CbmL1CATrdTrackFinderSA::createTracks
TStopwatch createTracks
Definition: CbmL1CATrdTrackFinderSA.h:156
CbmL1CATrdTrackFinderSA::fTotHits
std::map< Int_t, Int_t > fTotHits
Definition: CbmL1CATrdTrackFinderSA.h:303
CbmL1CATrdTrackFinderSA::fMomDistLongExtraX
TH2F * fMomDistLongExtraX
Definition: CbmL1CATrdTrackFinderSA.h:264
CbmTrack::GetParamFirst
const FairTrackParam * GetParamFirst() const
Definition: CbmTrack.h:61
CbmL1CATrdTrackFinderSA::Layer::scale
Double_t scale[12]
Definition: CbmL1CATrdTrackFinderSA.h:62
kTRDHIT
@ kTRDHIT
Definition: CbmHit.h:25
CbmL1CATrdTrackFinderSA::trdTrackFitterKF
CbmTrdTrackFitterKF * trdTrackFitterKF
Definition: CbmL1CATrdTrackFinderSA.h:107
CbmL1TrdTracklet4::GetVal
Int_t GetVal()
Definition: CbmL1TrdTracklet4.h:25
CbmTrdPoint
Definition: CbmTrdPoint.h:23
CbmL1TrdTracklet::SetVal
void SetVal(Int_t segVal)
Definition: CbmL1TrdTracklet.h:26
CbmL1CATrdTrackFinderSA::fYat0
TH1F * fYat0
Definition: CbmL1CATrdTrackFinderSA.h:287
CbmL1CATrdTrackFinderSA::fDistLongBY
TH1F * fDistLongBY
Definition: CbmL1CATrdTrackFinderSA.h:256
CbmL1CATrdTrackFinderSA::fPlane1Ydens
TH1F * fPlane1Ydens
Definition: CbmL1CATrdTrackFinderSA.h:280
CbmTrdTrack
Definition: CbmTrdTrack.h:22
CbmL1CATrdTrackFinderSA::fMomDistLongPrimaryY
TH2F * fMomDistLongPrimaryY
Definition: CbmL1CATrdTrackFinderSA.h:263
CbmL1TrdTracklet::SetIndRight
void SetIndRight(Int_t indRight)
Definition: CbmL1TrdTracklet.h:25
CbmL1CATrdTrackFinderSA::LayerWithHits::DY
Double_t DY
Definition: CbmL1CATrdTrackFinderSA.h:81
CbmL1CATrdTrackFinderSA::fDistLongY
TH1F * fDistLongY
Definition: CbmL1CATrdTrackFinderSA.h:251
CbmTrdTrack::SetELoss
void SetELoss(Double_t eLoss)
Definition: CbmTrdTrack.h:43
counter
int counter
Definition: CbmMvdSensorDigiToHitTask.cxx:72
CbmL1CATrdTrackFinderSA::fRUnUsedHits
std::map< Int_t, Int_t > fRUnUsedHits
Definition: CbmL1CATrdTrackFinderSA.h:302
CbmL1CATrdTrackFinderSA::LayerWithHits::Y
Double_t Y
Definition: CbmL1CATrdTrackFinderSA.h:78
CbmL1CATrdTrackFinderSA::PdgToCharge
Int_t PdgToCharge(Int_t pdgCode)
Definition: CbmL1CATrdTrackFinderSA.cxx:3191
CbmL1CATrdTrackFinderSA::fTrd31_Z
Double_t fTrd31_Z
Definition: CbmL1CATrdTrackFinderSA.h:150
CbmL1CATrdTrackFinderSA::CompareChi2
static Bool_t CompareChi2(TempTrackStruct A, TempTrackStruct B)
Definition: CbmL1CATrdTrackFinderSA.h:131
CbmL1CATrdTrackFinderSA::fSPlengthMC
TH1F * fSPlengthMC
Definition: CbmL1CATrdTrackFinderSA.h:285
CbmMCTrack.h
CbmL1CATrdTrackFinderSA::TempTrackStruct::Chi2
Double_t Chi2
Definition: CbmL1CATrdTrackFinderSA.h:126
CbmL1CATrdTrackFinderSA::FitTLinearFitter
Double_t FitTLinearFitter(Int_t M[])
Definition: CbmL1CATrdTrackFinderSA.cxx:3131
CbmL1CATrdTrackFinderSA::fDistY
TH1F * fDistY
Definition: CbmL1CATrdTrackFinderSA.h:277
CbmL1TrdTracklet4::SetVal
void SetVal(Int_t segVal)
Definition: CbmL1TrdTracklet4.h:40
CbmL1CATrdTrackFinderSA::createSegments
TStopwatch createSegments
Definition: CbmL1CATrdTrackFinderSA.h:155
CbmL1CATrdTrackFinderSA::CreateHistogramms
void CreateHistogramms()
Definition: CbmL1CATrdTrackFinderSA.cxx:1157
CbmL1CATrdTrackFinderSA::fUsedHitsPerPlane
TH1F * fUsedHitsPerPlane
Definition: CbmL1CATrdTrackFinderSA.h:298
x
Double_t x
Definition: CbmMvdSensorDigiToHitTask.cxx:68
CbmL1CATrdTrackFinderSA::FitKF
Double_t FitKF(CbmTrdTrack *pTrack)
Definition: CbmL1CATrdTrackFinderSA.cxx:3075
fabs
friend F32vec4 fabs(const F32vec4 &a)
Definition: L1/vectors/P4_F32vec4.h:60
y
Double_t y
Definition: CbmMvdSensorDigiToHitTask.cxx:68
CbmL1CATrdTrackFinderSA::totCreateSPs_SL
Double_t totCreateSPs_SL
Definition: CbmL1CATrdTrackFinderSA.h:159
CbmTrdPoint.h
CbmTrack::SetParamFirst
void SetParamFirst(const FairTrackParam *par)
Definition: CbmTrack.h:75
CbmL1TrdTracklet::Compare1
static Bool_t Compare1(CbmL1TrdTracklet *tr1, CbmL1TrdTracklet *tr2)
Definition: CbmL1TrdTracklet.h:39
CbmL1CATrdTrackFinderSA::~CbmL1CATrdTrackFinderSA
virtual ~CbmL1CATrdTrackFinderSA()
Definition: CbmL1CATrdTrackFinderSA.cxx:176
CbmL1TrdTracklet::SetPlanesID
void SetPlanesID(Int_t A_ID, Int_t B_ID)
Definition: CbmL1TrdTracklet.h:30
CbmL1CATrdTrackFinderSA::fMomDistExtrapolPrimaryX
TH2F * fMomDistExtrapolPrimaryX
Definition: CbmL1CATrdTrackFinderSA.h:267
CbmL1CATrdTrackFinderSA::totCreateSPs
Double_t totCreateSPs
Definition: CbmL1CATrdTrackFinderSA.h:159
CbmL1CATrdTrackFinderSA::totFindNeighbour
Double_t totFindNeighbour
Definition: CbmL1CATrdTrackFinderSA.h:159
CbmL1TrdTracklet::GetIndLeft
Int_t GetIndLeft()
Definition: CbmL1TrdTracklet.h:16
CbmTrdTrack.h
CbmL1CATrdTrackFinderSA::fDistLongBX
TH1F * fDistLongBX
Definition: CbmL1CATrdTrackFinderSA.h:255
z1
Double_t z1[nSects1]
Definition: pipe_v16a_mvdsts100.h:6
CbmL1CATrdTrackFinderSA::fRUsedHits
std::map< Int_t, Int_t > fRUsedHits
Definition: CbmL1CATrdTrackFinderSA.h:301
CbmKFTrdHit
Definition: CbmKFTrdHit.h:23
CbmL1CATrdTrackFinderSA::thirdLoopTime
TStopwatch thirdLoopTime
Definition: CbmL1CATrdTrackFinderSA.h:157
CbmL1CATrdTrackFinderSA::totCreateSegments
Double_t totCreateSegments
Definition: CbmL1CATrdTrackFinderSA.h:159
CbmL1CATrdTrackFinderSA::LayerWithHits::Z
Double_t Z
Definition: CbmL1CATrdTrackFinderSA.h:79
CbmL1CATrdTrackFinderSA::Rejection
Bool_t Rejection(Double_t Procent, Int_t num=100)
Definition: CbmL1CATrdTrackFinderSA.cxx:3299
CbmL1CATrdTrackFinderSA
Definition: CbmL1CATrdTrackFinderSA.h:35
CbmL1CATrdTrackFinderSA::Layer::Y
Double_t Y[12]
Definition: CbmL1CATrdTrackFinderSA.h:60
CbmL1TrdTracklet4::GetInd
Int_t GetInd(Int_t i)
Definition: CbmL1TrdTracklet4.h:19
CbmL1CATrdTrackFinderSA::totThirdLoopTime
Double_t totThirdLoopTime
Definition: CbmL1CATrdTrackFinderSA.h:161
CbmL1CATrdTrackFinderSA::totDoFind
Double_t totDoFind
Definition: CbmL1CATrdTrackFinderSA.h:160
CbmL1TrdTracklet::SetIndex
void SetIndex(Int_t index)
Definition: CbmL1TrdTracklet.h:28
CbmKFTrackInterface::Fit
Int_t Fit(Bool_t downstream=1)
Definition: CbmKFTrackInterface.cxx:101
CbmTrack::AddHit
void AddHit(Int_t index, HitType type)
Definition: CbmTrack.cxx:75
CbmL1CATrdTrackFinderSA::fSPlength
TH1F * fSPlength
Definition: CbmL1CATrdTrackFinderSA.h:284
CbmL1CATrdTrackFinderSA::fh_SP_xDiff_MC
TH1F * fh_SP_xDiff_MC
Definition: CbmL1CATrdTrackFinderSA.h:292
CbmKFTrack::SetTrackParam
void SetTrackParam(const FairTrackParam &track)
Definition: CbmKFTrack.cxx:34
CbmL1CATrdTrackFinderSA::LayerWithHits::hitIndex
Int_t hitIndex
Definition: CbmL1CATrdTrackFinderSA.h:75
CbmL1CATrdTrackFinderSA::secondLoopTime
TStopwatch secondLoopTime
Definition: CbmL1CATrdTrackFinderSA.h:157
CbmL1CATrdTrackFinderSA::fTrd14_Z
Double_t fTrd14_Z
Definition: CbmL1CATrdTrackFinderSA.h:150
CbmL1CATrdTrackFinderSA::fDistX
TH1F * fDistX
Definition: CbmL1CATrdTrackFinderSA.h:278
CbmL1CATrdTrackFinderSA::fTrd13_Z
Double_t fTrd13_Z
Definition: CbmL1CATrdTrackFinderSA.h:150
CbmL1TrdTracklet4.h
CbmL1CATrdTrackFinderSA::FindNeighbour
void FindNeighbour(std::vector< CbmL1TrdTracklet4 * > &v1, std::vector< CbmL1TrdTracklet4 * > &v2, Double_t dY, Double_t dX)
Definition: CbmL1CATrdTrackFinderSA.cxx:1703
CbmKFTrack::fHits
std::vector< CbmKFHit * > fHits
Definition: CbmKFTrack.h:29
CbmL1CATrdTrackFinderSA::fPlane9Ydens
TH1F * fPlane9Ydens
Definition: CbmL1CATrdTrackFinderSA.h:282
CbmL1TrdTracklet4::M
Double_t M[4]
Definition: CbmL1TrdTracklet4.h:73
CbmL1TrdTracklet
Definition: CbmL1TrdTracklet.h:7
CbmL1TrdTracklet4::SetCoord
void SetCoord(Int_t ind, Double_t val)
Definition: CbmL1TrdTracklet4.h:37
CbmL1CATrdTrackFinderSA::FitLinear
Double_t FitLinear(CbmTrdTrack *tr, Int_t var)
Definition: CbmL1CATrdTrackFinderSA.cxx:1523
CbmL1CATrdTrackFinderSA::fh_SP_xDiff_nMC
TH1F * fh_SP_xDiff_nMC
Definition: CbmL1CATrdTrackFinderSA.h:295
CbmL1CATrdTrackFinderSA::fMomDistLongPrimaryX
TH2F * fMomDistLongPrimaryX
Definition: CbmL1CATrdTrackFinderSA.h:262
CbmL1CATrdTrackFinderSA::TagSegments
void TagSegments(std::vector< CbmL1TrdTracklet4 * > &clTrackletsA, std::vector< CbmL1TrdTracklet4 * > &clTrackletsB, Int_t noCombSegments=0)
Definition: CbmL1CATrdTrackFinderSA.cxx:1831
CbmL1TrdTracklet4::SetInd
void SetInd(Int_t ind, Int_t val)
Definition: CbmL1TrdTracklet4.h:34
CbmTrdHit::GetPlaneId
Int_t GetPlaneId() const
Inherited from CbmBaseHit.
Definition: CbmTrdHit.h:73