CbmRoot
LxTBMLTask.cxx
Go to the documentation of this file.
1 /*
2  * To change this license header, choose License Headers in Project Properties.
3  * To change this template file, choose Tools | Templates
4  * and open the template in the editor.
5  */
6 
7 #include "LxTBMLTask.h"
8 #include "CbmCluster.h"
9 #include "CbmMCDataManager.h"
10 #include "CbmMCTrack.h"
11 #include "CbmMatch.h"
12 #include "CbmMuchGeoScheme.h"
13 #include "CbmMuchPixelHit.h"
14 #include "CbmMuchPoint.h"
15 #include "FairLogger.h"
16 #include "LxTBMatEffs.h"
17 #include "Simple/LxCA.h"
18 #include "Simple/LxSettings.h"
19 #include "TClonesArray.h"
20 #include "TDatabasePDG.h"
21 #include "TGeoArb8.h"
22 #include "TGeoBoolNode.h"
23 #include "TGeoCompositeShape.h"
24 #include "TGeoManager.h"
25 #include "TH1F.h"
26 #include "TMath.h"
27 #include "TRandom3.h"
28 
29 #ifdef __MACH__
30 #include <mach/mach_time.h>
31 #include <sys/time.h>
32 #define CLOCK_REALTIME 0
33 #define CLOCK_MONOTONIC 0
34 inline int clock_gettime(int clk_id, struct timespec* t) {
35  mach_timebase_info_data_t timebase;
36  mach_timebase_info(&timebase);
37  uint64_t time;
38  time = mach_absolute_time();
39  double nseconds =
40  ((double) time * (double) timebase.numer) / ((double) timebase.denom);
41  double seconds =
42  ((double) time * (double) timebase.numer) / ((double) timebase.denom * 1e9);
43  t->tv_sec = seconds;
44  t->tv_nsec = nseconds;
45  return 0;
46 }
47 #else
48 #include <time.h>
49 #endif
50 
51 using namespace std;
52 
53 ClassImp(LxTBMLFinder)
54 
55  Double_t speedOfLight = 0;
56 //static scaltype magneticFieldCorrections[] = { 1.5, 1.0, 1.7, 1.1 };
57 static scaltype magneticFieldCorrections[] = {0, 0, 0, 0};
58 
59 #ifdef LXTB_DEBUG
60 struct DebugTrack2 {
61  list<LxTbBinnedPoint*> points[NOF_STATIONS][NOF_LAYERS];
62 };
63 
64 static map<Int_t, DebugTrack2> debugTracks2;
65 static int nofOverXR = 0;
66 static int nofOverYR = 0;
67 static int nofOverTR = 0;
68 
69 static int nofOverXL = 0;
70 static int nofOverYL = 0;
71 static int nofOverTL = 0;
72 
73 static int nofFoundR = 0;
74 static int nofNotFoundR = 0;
75 
76 vector<vector<LxTBMLFinder::TrackDataHolder>>* gMCTracks = 0;
77 #endif //LXTB_DEBUG
78 
79 struct LxTbMLStation {
80  struct Q {
81  scaltype Q11, Q12, Q21, Q22;
82  };
83 
84  int fStationNumber;
85  LxTbLayer* fLayers;
86  LxTbAbsorber fAbsorber;
87  scaltype fDeltaThetaX;
88  scaltype fThetaX;
89  scaltype fScatXRL;
90  scaltype fScatXLS;
91  scaltype fDeltaThetaY;
92  scaltype fThetaY;
93  scaltype fScatYRL;
94  scaltype fScatYLS;
95  Q qs[2];
96 
97  LxTbMLStation(int stationNumber, int nofxb, int nofyb, int noftb)
98  : fStationNumber(stationNumber)
99  , fLayers(reinterpret_cast<LxTbLayer*>(
100  new unsigned char[NOF_LAYERS * sizeof(LxTbLayer)]))
101  , fDeltaThetaX(0)
102  , fThetaX(0)
103  , fScatXRL(0)
104  , fScatXLS(0)
105  , fDeltaThetaY(0)
106  , fThetaY(0)
107  , fScatYRL(0)
108  , fScatYLS(0)
109  , fHandleMPoint(*this)
110  , fHandleRPoint(*this)
111  , fHandleLPoint(*this) {
112  for (int i = 0; i < NOF_LAYERS; ++i)
113  new (&fLayers[i]) LxTbLayer(nofxb, nofyb, noftb);
114  }
115 
116  void Init() {
117  for (int i = 0; i < NOF_LAYERS; ++i) {
118  //fLayers[i].xBinLength = (fLayers[i].maxX - fLayers[i].minX) / fLayers[i].nofXBins;
119  //fLayers[i].yBinLength = (fLayers[i].maxY - fLayers[i].minY) / fLayers[i].nofYXBins;
120  fLayers[i].Init();
121  }
122  }
123 
124  void Init2() {
125  fHandleMPoint.Init();
126  fHandleRPoint.Init();
127  fHandleLPoint.Init();
128  }
129 
130  void Clear() {
131  for (int i = 0; i < NOF_LAYERS; ++i)
132  fLayers[i].Clear();
133  }
134 
135  void SetMinT(timetype v) {
136  for (int i = 0; i < NOF_LAYERS; ++i)
137  fLayers[i].SetMinT(v);
138  }
139 
140  struct HandleMPoint {
141  LxTbMLStation& station;
142  timetype c;
143  scaltype deltaZr;
144  scaltype scatXRL;
145 
146  explicit HandleMPoint(LxTbMLStation& parent)
147  : station(parent), c(0), deltaZr(0), scatXRL(0) {}
148 
149  void Init() {
150  c = speedOfLight;
151  deltaZr = station.fLayers[2].z - station.fLayers[1].z;
152  scatXRL = sqrt(station.fScatXRL * station.fScatXRL
153  + magneticFieldCorrections[station.fStationNumber]
154  * magneticFieldCorrections[station.fStationNumber]);
155  }
156 
157  void operator()(LxTbBinnedPoint& point) {
158  scaltype txR = point.x / station.fLayers[1].z;
159  scaltype tyR = point.y / station.fLayers[1].z;
160  scaltype pXr = point.x + txR * deltaZr;
161  scaltype pYr = point.y + tyR * deltaZr;
162  scaltype trajLenR = sqrt(1 + txR * txR + tyR * tyR) * deltaZr;
163  timetype pTr = point.t + 1.e9 * trajLenR / c; // 1.e9 to convert to ns.
164  station.fHandleRPoint.mPoint = &point;
165  station.fHandleLPoint.mPoint = &point;
166 #ifdef LXTB_DEBUG
167  struct Debug {
168  list<LxTbBinnedPoint*> points;
169 
170  void operator()(LxTbBinnedPoint& point) { points.push_back(&point); }
171  };
172 
173  Debug debug;
174  IterateNeighbourhood(station.fLayers[2],
175  pXr,
176  point.dx,
177  scatXRL,
178  pYr,
179  point.dy,
180  station.fScatYRL,
181  pTr,
182  point.dt,
183  debug);
184 
185  if (LAST_STATION == station.fStationNumber) {
186  bool isSignal = false;
187  bool found = false;
188 
189  for (list<LxTbBinnedPoint::PointDesc>::const_iterator i =
190  point.mcRefs.begin();
191  i != point.mcRefs.end();
192  ++i) {
193  const LxTbBinnedPoint::PointDesc& pd = *i;
194  const vector<vector<LxTBMLFinder::TrackDataHolder>>& mcTracks =
195  *gMCTracks;
196 
197  if (mcTracks[pd.eventId][pd.trackId].isSignal) isSignal = true;
198 
199  for (list<LxTbBinnedPoint*>::iterator j =
200  debugTracks2[pd.trackId]
201  .points[point.stationNumber][2]
202  .begin();
203  j
204  != debugTracks2[pd.trackId].points[point.stationNumber][2].end();
205  ++j) {
206  LxTbBinnedPoint* a = *j;
207 
208  for (list<LxTbBinnedPoint*>::iterator k = debug.points.begin();
209  k != debug.points.end();
210  ++k) {
211  LxTbBinnedPoint* b = *k;
212 
213  if (a == b) found = true;
214  }
215  }
216  }
217 
218  if (isSignal) {
219  if (!found)
220  ++nofNotFoundR;
221  else
222  ++nofFoundR;
223  }
224 
225  struct Debug2 {
226  void operator()(LxTbBinnedPoint& point) {
227  cout << "Debug2: Point: (" << point.x << ", " << point.y << ", "
228  << point.t << ")" << endl;
229  }
230  };
231 
232  Debug2 debug2;
233  IterateLayer(station.fLayers[2], debug2);
234  }
235 #endif //LXTB_DEBUG
236  IterateNeighbourhood(station.fLayers[2],
237  pXr,
238  point.dx,
239  scatXRL,
240  pYr,
241  point.dy,
242  station.fScatYRL,
243  pTr,
244  point.dt,
245  station.fHandleRPoint);
246 #ifdef LXTB_DEBUG
247  int qq = 0;
248  qq += 10;
249 #endif //LXTB_DEBUG
250  }
251  };
252 
253  HandleMPoint fHandleMPoint;
254 
255  struct HandleRPoint {
256  LxTbMLStation& station;
257  timetype c;
258  scaltype deltaZl;
259  scaltype scatXLL;
260  LxTbBinnedPoint* mPoint;
261 
262  explicit HandleRPoint(LxTbMLStation& parent)
263  : station(parent), c(0), deltaZl(0), scatXLL(0), mPoint(0) {}
264 
265  void Init() {
266  c = speedOfLight;
267  deltaZl = station.fLayers[0].z - station.fLayers[1].z;
268  scatXLL = magneticFieldCorrections[station.fStationNumber];
269  }
270 
271  void operator()(LxTbBinnedPoint& point) {
272  scaltype txL =
273  (point.x - mPoint->x) / (station.fLayers[2].z - station.fLayers[1].z);
274  scaltype tyL =
275  (point.y - mPoint->y) / (station.fLayers[2].z - station.fLayers[1].z);
276  scaltype pXl = mPoint->x + txL * deltaZl;
277  scaltype pYl = mPoint->y + tyL * deltaZl;
278  scaltype trajLenL = sqrt(1 + txL * txL + tyL * tyL) * deltaZl;
279  timetype pTl =
280  mPoint->t
281  + 1.e9 * trajLenL
282  / c; // 1.e9 to convert to ns and trajLenL is negative.
283  station.fHandleLPoint.rPoint = &point;
284  IterateNeighbourhood(
285  station.fLayers[0],
286  pXl,
287  sqrt(mPoint->dx * mPoint->dx + point.dx * point.dx),
288  scatXLL,
289  pYl,
290  sqrt(mPoint->dy * mPoint->dy + point.dy * point.dy),
291  0,
292  pTl,
293  mPoint->dt /*sqrt(mPoint->dt * mPoint->dt + point.dt * point.dt) / 2*/,
294  station.fHandleLPoint);
295  }
296  };
297 
298  HandleRPoint fHandleRPoint;
299 
300  struct HandleLPoint {
301  LxTbMLStation& station;
302  scaltype deltaZ;
303  LxTbBinnedPoint* mPoint;
304  LxTbBinnedPoint* rPoint;
305 
306  explicit HandleLPoint(LxTbMLStation& parent)
307  : station(parent)
308  , deltaZ(parent.fLayers[2].z - parent.fLayers[0].z)
309  , mPoint(0)
310  , rPoint(0) {}
311 
312  void Init() { deltaZ = station.fLayers[0].z - station.fLayers[2].z; }
313 
314  void operator()(LxTbBinnedPoint& point) {
315  LxTbBinnedTriplet* triplet =
316  new LxTbBinnedTriplet(&point, rPoint, deltaZ);
317  mPoint->triplets.push_back(triplet);
318  }
319  };
320 
321  HandleLPoint fHandleLPoint;
322 
323  void Reconstruct() { IterateLayer<HandleMPoint>(fLayers[1], fHandleMPoint); }
324 };
325 
326 static void
327 FindGeoChild(TGeoNode* node, const char* name, list<TGeoNode*>& results) {
328  Int_t nofChildren = node->GetNdaughters();
329 
330  for (Int_t i = 0; i < nofChildren; ++i) {
331  TGeoNode* child = node->GetDaughter(i);
332  TString childName(child->GetName());
333 
334  if (childName.Contains(name, TString::kIgnoreCase))
335  results.push_back(child);
336  }
337 }
338 
339 struct SignalParticle {
340  const char* fName;
341  Int_t fPdgCode;
342  scaltype fMinEnergy;
343 };
344 
345 static SignalParticle particleDescs[] = {{"omega", 223, 1.5}, {"", -1, 0}};
346 
347 struct LxTbDetector {
348  LxTbMLStation* fStations;
349  list<LxTBMLFinder::Chain*> recoTracks;
350  SignalParticle* fSignalParticle;
351 
352  LxTbDetector(int nofxb, int nofyb, int noftb)
353  : fStations(reinterpret_cast<LxTbMLStation*>(
354  new unsigned char[NOF_STATIONS * sizeof(LxTbMLStation)]))
355  , fSignalParticle(&particleDescs[0])
356  , fHandleLastPoint(*this) {
357  for (int i = 0; i < NOF_STATIONS; ++i)
358  new (&fStations[i]) LxTbMLStation(i, nofxb, nofyb, noftb);
359  }
360 
361  void Init() {
362  speedOfLight =
363  100 * TMath::C(); // Multiply by 100 to express in centimeters.
364  gMuonMass = TDatabasePDG::Instance()->GetParticle(13)->Mass();
365  gElectronMass = TDatabasePDG::Instance()->GetParticle(11)->Mass();
366  HandleGeometry();
367 
368  for (int i = 0; i < NOF_STATIONS; ++i) {
369  LxTbMLStation& station = fStations[i];
370  station.Init();
371  }
372 
373  scaltype E = fSignalParticle->fMinEnergy; // GeV
374  scaltype E0 = E;
375 
376  for (int i = 0; i < NOF_STATIONS; ++i) {
377  LxTbMLStation& station = fStations[i];
378  scaltype L = station.fAbsorber.width; // / cos(3.14159265 * 15 / 180);
379  E -= EnergyLoss(E, L, &station.fAbsorber);
380  scaltype Escat = (E0 + E) / 2;
381  //scaltype Escat = E;
382  scaltype deltaTheta = CalcThetaPrj(Escat, L, &station.fAbsorber);
383  station.fDeltaThetaX = deltaTheta;
384  station.fDeltaThetaY = deltaTheta;
385 
386  if (i > 0) {
387  scaltype deltaZLS =
388  station.fLayers[1].z - fStations[i - 1].fLayers[1].z;
389  station.fScatXLS = station.fDeltaThetaX * deltaZLS;
390  station.fScatYLS = station.fDeltaThetaY * deltaZLS;
391  }
392 
393  scaltype q0XSq = station.fDeltaThetaX * station.fDeltaThetaX;
394  scaltype q0YSq = station.fDeltaThetaY * station.fDeltaThetaY;
395  station.qs[0] = {q0XSq * L * L / 3, q0XSq * L / 2, q0XSq * L / 2, q0XSq};
396  station.qs[1] = {q0YSq * L * L / 3, q0YSq * L / 2, q0YSq * L / 2, q0YSq};
397  E0 = E;
398  }
399 
400  scaltype Ls[NOF_STATIONS + 1];
401  Ls[0] = fStations[0].fAbsorber.zCoord;
402 
403  for (int i = 1; i < NOF_STATIONS; ++i)
404  Ls[i] = fStations[i].fAbsorber.zCoord - Ls[i - 1];
405 
406  Ls[NOF_STATIONS] = fStations[LAST_STATION].fLayers[1].z - Ls[LAST_STATION];
407 
408  for (int s = 0; s < NOF_STATIONS; ++s) {
409  LxTbMLStation& station = fStations[s];
410  scaltype L = station.fLayers[1].z;
411  int n = s + 1;
412  scaltype thetaXSq = 0;
413  scaltype thetaYSq = 0;
414 
415  for (int i = 1; i <= n; ++i) {
416  scaltype sumLi = 0;
417 
418  for (int j = 0; j < i; ++j)
419  sumLi += Ls[j];
420 
421  thetaXSq +=
422  sumLi * sumLi * fStations[i].fDeltaThetaX * fStations[i].fDeltaThetaX;
423  thetaYSq +=
424  sumLi * sumLi * fStations[i].fDeltaThetaY * fStations[i].fDeltaThetaY;
425  }
426 
427  station.fThetaX = sqrt(thetaXSq) / L;
428  station.fThetaY = sqrt(thetaYSq) / L;
429  scaltype deltaZRL = station.fLayers[2].z - station.fLayers[1].z;
430  station.fScatXRL = station.fThetaX * deltaZRL;
431  station.fScatYRL = station.fThetaY * deltaZRL;
432  }
433 
434  for (int i = 0; i < NOF_STATIONS; ++i) {
435  LxTbMLStation& station = fStations[i];
436  station.Init2();
437  }
438 
439  fHandleRPoint.Init();
440  }
441 
442  void Clear() {
443  recoTracks.clear();
444 
445  for (int i = 0; i < NOF_STATIONS; ++i)
446  fStations[i].Clear();
447  }
448 
449  void SetMinT(timetype v) {
450  for (int i = 0; i < NOF_STATIONS; ++i)
451  fStations[i].SetMinT(v);
452  }
453 
454  void HandleGeometry() {
455  Double_t localCoords[3] = {0., 0., 0.};
456  Double_t globalCoords[3];
457  TGeoNavigator* pNavigator = gGeoManager->GetCurrentNavigator();
458  gGeoManager->cd("/cave_1");
459  list<TGeoNode*> detectors;
460  FindGeoChild(gGeoManager->GetCurrentNode(), "much", detectors);
461 
462  for (list<TGeoNode*>::iterator i = detectors.begin(); i != detectors.end();
463  ++i) {
464  TGeoNode* detector = *i;
465  pNavigator->CdDown(detector);
466  list<TGeoNode*> stations;
467  FindGeoChild(detector, "station", stations);
468  int stationNumber = 0;
469 
470  for (list<TGeoNode*>::iterator j = stations.begin(); j != stations.end();
471  ++j) {
472  TGeoNode* station = *j;
473  pNavigator->CdDown(station);
474  list<TGeoNode*> layers;
475  FindGeoChild(station, "layer", layers);
476  int layerNumber = 0;
477 
478  for (list<TGeoNode*>::iterator k = layers.begin(); k != layers.end();
479  ++k) {
480  TGeoNode* layer = *k;
481  pNavigator->CdDown(layer);
482  gGeoManager->LocalToMaster(localCoords, globalCoords);
483 
484  fStations[stationNumber].fLayers[layerNumber].z = globalCoords[2];
485  fStations[stationNumber].fLayers[layerNumber].minX = 0;
486  fStations[stationNumber].fLayers[layerNumber].maxX = 0;
487  fStations[stationNumber].fLayers[layerNumber].xBinLength = 0;
488  fStations[stationNumber].fLayers[layerNumber].minY = 0;
489  fStations[stationNumber].fLayers[layerNumber].maxY = 0;
490  fStations[stationNumber].fLayers[layerNumber].yBinLength = 0;
491 
492  list<TGeoNode*> actives;
493  FindGeoChild(layer, "active", actives);
494 
495  for (list<TGeoNode*>::iterator l = actives.begin();
496  l != actives.end();
497  ++l) {
498  TGeoNode* active = *l;
499  pNavigator->CdDown(active);
500  TGeoCompositeShape* cs = dynamic_cast<TGeoCompositeShape*>(
501  active->GetVolume()->GetShape());
502  TGeoBoolNode* bn = cs->GetBoolNode();
503  TGeoTrap* trap = dynamic_cast<TGeoTrap*>(bn->GetLeftShape());
504 
505  if (0 != trap) {
506  Double_t* xy = trap->GetVertices();
507 
508  for (int m = 0; m < 4; ++m) {
509  Double_t localActiveCoords[3] = {xy[2 * m], xy[2 * m + 1], 0.};
510  Double_t globalActiveCoords[3];
511  gGeoManager->LocalToMaster(localActiveCoords,
512  globalActiveCoords);
513 
514  if (fStations[stationNumber].fLayers[layerNumber].minY
515  > globalActiveCoords[1])
516  fStations[stationNumber].fLayers[layerNumber].minY =
517  globalActiveCoords[1];
518 
519  if (fStations[stationNumber].fLayers[layerNumber].maxY
520  < globalActiveCoords[1])
521  fStations[stationNumber].fLayers[layerNumber].maxY =
522  globalActiveCoords[1];
523 
524  if (fStations[stationNumber].fLayers[layerNumber].minX
525  > globalActiveCoords[0])
526  fStations[stationNumber].fLayers[layerNumber].minX =
527  globalActiveCoords[0];
528 
529  if (fStations[stationNumber].fLayers[layerNumber].maxX
530  < globalActiveCoords[0])
531  fStations[stationNumber].fLayers[layerNumber].maxX =
532  globalActiveCoords[0];
533  }
534  }
535 
536  pNavigator->CdUp();
537  }
538 
539  pNavigator->CdUp();
540  ++layerNumber;
541  }
542 
543  ++stationNumber;
544  pNavigator->CdUp();
545  }
546 
547  int nofStations = stationNumber;
548 
549  list<TGeoNode*> absorbers;
550  FindGeoChild(detector, "absorber", absorbers);
551  int absorberNumber = 0;
552 
553  for (list<TGeoNode*>::iterator j = absorbers.begin();
554  j != absorbers.end();
555  ++j) {
556  TGeoNode* absorber = *j;
557  pNavigator->CdDown(absorber);
558  TGeoVolume* absVol = gGeoManager->GetCurrentVolume();
559  const TGeoBBox* absShape =
560  static_cast<const TGeoBBox*>(absVol->GetShape());
561 
562  if (absorberNumber < nofStations) {
563  gGeoManager->LocalToMaster(localCoords, globalCoords);
564  fStations[absorberNumber].fAbsorber.zCoord =
565  globalCoords[2] - absShape->GetDZ();
566  fStations[absorberNumber].fAbsorber.width = 2 * absShape->GetDZ();
567  fStations[absorberNumber].fAbsorber.radLength =
568  absVol->GetMaterial()->GetRadLen();
569  fStations[absorberNumber].fAbsorber.rho =
570  absVol->GetMaterial()->GetDensity();
571  fStations[absorberNumber].fAbsorber.Z = absVol->GetMaterial()->GetZ();
572  fStations[absorberNumber].fAbsorber.A = absVol->GetMaterial()->GetA();
573  }
574 
575  ++absorberNumber;
576  pNavigator->CdUp();
577  }
578 
579  pNavigator->CdUp();
580  }
581  } // void HandleGeometry()
582 
583  struct HandleRPoint {
584  LxTbMLStation* rStation;
585  LxTbLayer* lLayer;
586  scaltype deltaZ;
587  timetype c;
588 
589  HandleRPoint() : rStation(0), lLayer(0), deltaZ(0), c(0) {}
590 
591  void Init() { c = speedOfLight; }
592 
593  void operator()(LxTbBinnedPoint& point) {
594  for (list<LxTbBinnedTriplet*>::iterator i = point.triplets.begin();
595  i != point.triplets.end();
596  ++i) {
597  LxTbBinnedTriplet* triplet = *i;
598  scaltype x = point.x + triplet->tx * deltaZ;
599  scaltype y = point.y + triplet->ty * deltaZ;
600  scaltype trajLen =
601  sqrt(1 + triplet->tx * triplet->tx + triplet->ty * triplet->ty)
602  * deltaZ;
603  timetype t =
604  point.t
605  + 1.e9 * trajLen
606  / c; // 1.e9 to convert to ns and trajLenL is negative.
607  handleLPoint.rTriplet = triplet;
608  IterateNeighbourhood(
609  *lLayer,
610  x,
611  sqrt(point.dx * point.dx
612  + triplet->dtx * triplet->dtx * deltaZ * deltaZ),
613  rStation->fDeltaThetaX * deltaZ,
614  y,
615  sqrt(point.dy * point.dy
616  + triplet->dty * triplet->dty * deltaZ * deltaZ),
617  rStation->fDeltaThetaY * deltaZ,
618  t,
619  point.dt,
620  handleLPoint);
621  }
622  }
623 
624  struct HandleLPoint {
625  LxTbBinnedTriplet* rTriplet;
626  HandleLPoint() : rTriplet(0) {}
627 
628  void operator()(LxTbBinnedPoint& point) {
629  rTriplet->neighbours.push_back(&point);
630  }
631  };
632 
633  HandleLPoint handleLPoint;
634  };
635 
636  HandleRPoint fHandleRPoint;
637 
638  struct HandleLastPoint {
639  struct KFParamsCoord {
640  scaltype coord, tg, C11, C12, C21, C22;
641 
642  void Clear() {
643  coord = 0;
644  tg = 0;
645  C11 = 1.0;
646  C12 = 0;
647  C21 = 0;
648  C22 = 0;
649  }
650  };
651 
652  struct KFParams {
653  KFParamsCoord xParams;
654  KFParamsCoord yParams;
655  scaltype chi2;
656 
657  void Clear() {
658  xParams.Clear();
659  yParams.Clear();
660  chi2 = 0;
661  }
662  };
663 
664  LxTbDetector& detector;
665 
666  explicit HandleLastPoint(LxTbDetector& parent) : detector(parent) {}
667 
668  KFParamsCoord KFAddPointCoord(KFParamsCoord prevParam,
669  scaltype m,
670  scaltype V,
671  scaltype& chi2,
672  int stationNumber,
673  int layerNumber,
674  int coordNumber) {
675  KFParamsCoord param = prevParam;
676  const LxTbMLStation& station = detector.fStations[stationNumber];
677  const LxTbLayer& layer = station.fLayers[layerNumber];
678  const LxTbMLStation::Q& Q = station.qs[coordNumber];
679  scaltype deltaZ =
680  LAST_LAYER == layerNumber
681  ? LAST_STATION == stationNumber
682  ? 0
683  : layer.z - detector.fStations[stationNumber + 1].fLayers[0].z
684  : layer.z - station.fLayers[layerNumber + 1].z;
685  scaltype deltaZSq = deltaZ * deltaZ;
686 
687  // Extrapolate.
688  param.coord += prevParam.tg * deltaZ; // params[k].tg is unchanged.
689 
690  // Filter.
691  if (LAST_LAYER == layerNumber) {
692  param.C11 += prevParam.C12 * deltaZ + prevParam.C21 * deltaZ
693  + prevParam.C22 * deltaZSq + Q.Q11;
694  param.C12 += prevParam.C22 * deltaZ + Q.Q12;
695  param.C21 += prevParam.C22 * deltaZ + Q.Q21;
696  param.C22 += Q.Q22;
697  } else {
698  param.C11 += prevParam.C12 * deltaZ + prevParam.C21 * deltaZ
699  + prevParam.C22 * deltaZSq;
700  param.C12 += prevParam.C22 * deltaZ;
701  param.C21 += prevParam.C22 * deltaZ;
702  }
703 
704  scaltype S = 1.0 / (V + param.C11);
705  scaltype Kcoord = param.C11 * S;
706  scaltype Ktg = param.C21 * S;
707  scaltype dzeta = m - param.coord;
708  param.coord += Kcoord * dzeta;
709  param.tg += Ktg * dzeta;
710  param.C21 -= param.C11 * Ktg;
711  param.C22 -= param.C12 * Ktg;
712  param.C11 *= 1.0 - Kcoord;
713  param.C12 *= 1.0 - Kcoord;
714  chi2 += dzeta * S * dzeta;
715  return param;
716  }
717 
718  KFParams KFAddPoint(KFParams prevParam,
719  scaltype m[2],
720  scaltype V[2],
721  int stationNumber,
722  int layerNumber) {
723  KFParams param = {KFAddPointCoord(prevParam.xParams,
724  m[0],
725  V[0],
726  param.chi2,
727  stationNumber,
728  layerNumber,
729  0),
730  KFAddPointCoord(prevParam.yParams,
731  m[1],
732  V[1],
733  param.chi2,
734  stationNumber,
735  layerNumber,
736  1)};
737  return param;
738  }
739 
740  KFParams KFAddTriplet(
741  KFParams param,
742  LxTbBinnedPoint* trackCandidatePoints[NOF_STATIONS][NOF_LAYERS],
743  int level) {
744  LxTbBinnedPoint** points = trackCandidatePoints[level];
745 
746  for (int i = LAST_LAYER; i >= 0; --i) {
747  LxTbBinnedPoint* point = points[i];
748  scaltype m[2] = {point->x, point->y};
749  scaltype V[2] = {point->dx * point->dx, point->dy * point->dy};
750  param = KFAddPoint(param, m, V, level, i);
751  }
752 
753  return param;
754  }
755 
756  void HandleTriplet(
757  LxTbBinnedTriplet* triplet,
758  LxTbBinnedPoint* trackCandidatePoints[NOF_STATIONS][NOF_LAYERS],
759  list<LxTBMLFinder::Chain>& chains,
760  int level,
761  KFParams kfParams) {
762  trackCandidatePoints[level][0] = triplet->lPoint;
763  trackCandidatePoints[level][2] = triplet->rPoint;
764  kfParams = KFAddTriplet(kfParams, trackCandidatePoints, level);
765 
766  if (0 == level)
767  chains.push_back(
768  LxTBMLFinder::Chain(trackCandidatePoints, kfParams.chi2));
769  else {
770  for (list<LxTbBinnedPoint*>::iterator i = triplet->neighbours.begin();
771  i != triplet->neighbours.end();
772  ++i)
773  HandlePoint(*i, trackCandidatePoints, chains, level - 1, kfParams);
774  }
775  }
776 
777  void
778  HandlePoint(LxTbBinnedPoint* point,
779  LxTbBinnedPoint* trackCandidatePoints[NOF_STATIONS][NOF_LAYERS],
780  list<LxTBMLFinder::Chain>& chains,
781  int level,
782  KFParams kfParams) {
783  trackCandidatePoints[level][1] = point;
784 
785  for (list<LxTbBinnedTriplet*>::iterator i = point->triplets.begin();
786  i != point->triplets.end();
787  ++i) {
788  LxTbBinnedTriplet* triplet = *i;
789  HandleTriplet(triplet, trackCandidatePoints, chains, level, kfParams);
790  }
791  }
792 
793  void operator()(LxTbBinnedPoint& point) {
794  LxTbBinnedPoint* trackCandidatePoints[NOF_STATIONS][NOF_LAYERS];
795  list<LxTBMLFinder::Chain> chains;
796  KFParams kfParams = {{0, 0, 1.0, 0, 0, 1.0}, {0, 0, 1.0, 0, 0, 1.0}, 0};
797  HandlePoint(&point, trackCandidatePoints, chains, LAST_STATION, kfParams);
798  const LxTBMLFinder::Chain* bestChain = 0;
799  scaltype chi2 = 0;
800 
801  for (list<LxTBMLFinder::Chain>::const_iterator i = chains.begin();
802  i != chains.end();
803  ++i) {
804  const LxTBMLFinder::Chain& chain = *i;
805 
806  if (0 == bestChain || chain.chi2 < chi2) {
807  bestChain = &chain;
808  chi2 = chain.chi2;
809  }
810  }
811 
812  if (0 != bestChain)
813  detector.recoTracks.push_back(new LxTBMLFinder::Chain(*bestChain));
814  }
815  };
816 
817  HandleLastPoint fHandleLastPoint;
818 
819 #ifdef LXTB_DEBUG
820  struct Debug {
821  void operator()(LxTbBinnedPoint& point) {
822  cout << "Point Point Point!!!" << endl;
823 
824  for (list<LxTbBinnedTriplet*>::iterator i = point.triplets.begin();
825  i != point.triplets.end();
826  ++i)
827  cout << "Triplet Triplet Triplet!!!" << endl;
828  }
829  };
830 
831  Debug debug;
832 #endif //LXTB_DEBUG
833 
834  void Reconstruct() {
835  for (int i = LAST_STATION; i >= 0; --i) {
836  LxTbMLStation& rStation = fStations[i];
837  rStation.Reconstruct();
838 
839 #ifdef LXTB_DEBUG
840  if (LAST_STATION == i) {
841  cout << "Points and triplets dump for the station #:" << i << endl;
842  IterateLayer(rStation.fLayers[0], debug);
843  IterateLayer(rStation.fLayers[1], debug);
844  IterateLayer(rStation.fLayers[2], debug);
845  }
846 #endif //LXTB_DEBUG
847 
848  if (i > 0) {
849  LxTbLayer& rLayer = rStation.fLayers[1];
850  LxTbMLStation& lStation = fStations[i - 1];
851  LxTbLayer& lLayer = lStation.fLayers[1];
852  fHandleRPoint.rStation = &rStation;
853  fHandleRPoint.lLayer = &lLayer;
854  fHandleRPoint.deltaZ = lLayer.z - rLayer.z;
855  IterateLayer(rLayer, fHandleRPoint);
856  }
857  }
858 
859  LxTbMLStation& lastStation = fStations[LAST_STATION];
860  LxTbLayer& lastLayer = lastStation.fLayers[1];
861  IterateLayer(lastLayer, fHandleLastPoint);
862  }
863 };
864 
865 static LxTbMLStation* gStations;
866 
867 LxTBMLFinder::LxTBMLFinder()
868  : fReconstructor(0)
869  , fIsEvByEv(true)
870  , fNofXBins(20)
871  , fNofYBins(20)
872  , fNofTBins(fIsEvByEv ? 5 : 1000)
873  , fNEvents(1000) {}
874 
875 #ifdef LXTB_DEBUG
876 TH1F* deltaXRHisto[NOF_STATIONS];
877 TH1F* deltaYRHisto[NOF_STATIONS];
878 TH1F* deltaTRHisto[NOF_STATIONS];
879 
880 TH1F* deltaXLHisto[NOF_STATIONS];
881 TH1F* deltaYLHisto[NOF_STATIONS];
882 TH1F* deltaTLHisto[NOF_STATIONS];
883 #endif //LXTB_DEBUG
884 
885 InitStatus LxTBMLFinder::Init() {
886 #ifdef LXTB_DEBUG
887  gMCTracks = &fMCTracks;
888 #endif //LXTB_DEBUG
889  FairRootManager* ioman = FairRootManager::Instance();
890 
891  if (0 == ioman) LOG(fatal) << "No FairRootManager";
892 
893  int nofEventsInFile = ioman->CheckMaxEventNo();
894 
895  if (nofEventsInFile < fNEvents) fNEvents = nofEventsInFile;
896 
897  LxTbDetector* pReconstructor =
898  new LxTbDetector(fNofXBins, fNofYBins, fNofTBins);
899  fReconstructor = pReconstructor;
900  pReconstructor->Init();
901 
902  fMuchPixelHits = static_cast<TClonesArray*>(ioman->GetObject("MuchPixelHit"));
903 
904 #ifdef LXTB_QA
905  CbmMCDataManager* mcManager =
906  static_cast<CbmMCDataManager*>(ioman->GetObject("MCDataManager"));
907  fMuchMCPoints = mcManager->InitBranch("MuchPoint");
908  fMuchClusters = static_cast<TClonesArray*>(ioman->GetObject("MuchCluster"));
909  fMuchPixelDigiMatches =
910  static_cast<TClonesArray*>(ioman->GetObject("MuchDigiMatch"));
911  CbmMCDataArray* mcTracks = mcManager->InitBranch("MCTrack");
912 
913  for (int i = 0; i < fNEvents; ++i) {
914  Int_t evSize = mcTracks->Size(0, i);
915  fMCTracks.push_back(vector<TrackDataHolder>());
916 
917  if (0 >= evSize) continue;
918 
919  vector<TrackDataHolder>& evTracks = fMCTracks.back();
920  const CbmMCTrack* posTrack = 0;
921  const CbmMCTrack* negTrack = 0;
922 
923  for (int j = 0; j < evSize; ++j) {
924  evTracks.push_back(TrackDataHolder());
925  const CbmMCTrack* mcTrack =
926  static_cast<const CbmMCTrack*>(mcTracks->Get(0, i, j));
927 
928  if (mcTrack->GetPdgCode() == 13 || mcTrack->GetPdgCode() == -13) {
929  Double_t m = mcTrack->GetMass();
930  Int_t motherId = mcTrack->GetMotherId();
931 
932  if (motherId >= 0
933  && static_cast<const CbmMCTrack*>(mcTracks->Get(0, i, motherId))
934  ->GetPdgCode()
935  == pReconstructor->fSignalParticle->fPdgCode) {
936  //const CbmMCTrack* motherTrack = static_cast<const CbmMCTrack*> (mcTracks->Get(0, i, motherId));
937 
938  //if (fFinder->fSignalParticle->fPdgCode == motherTrack->GetPdgCode())
939  {
940  evTracks.back().isSignal = true;
941  evTracks.back().isPos = mcTrack->GetPdgCode() == -13;
942 
943  if (-13 == mcTrack->GetPdgCode())
944  posTrack = mcTrack;
945  else
946  negTrack = mcTrack;
947  }
948  }
949  }
950  } // for (int j = 0; j < evSize; ++j)
951  } // for (int i = 0; i < fNEvents; ++i)
952 
953  fEventTimes.resize(fNEvents);
954  fEventTimes[0] = 50;
955 
956  for (int i = 1; i < fNEvents; ++i)
957  fEventTimes[i] = fEventTimes[i - 1] + 100;
958 
959  for (int i = 0; i < fNEvents; ++i) {
960  Int_t evSize = fMuchMCPoints->Size(0, i);
961  fMuchPoints.push_back(vector<PointDataHolder>());
962 
963  if (0 >= evSize) continue;
964 
965  //++numEvents;
966  vector<PointDataHolder>& evPoints = fMuchPoints.back();
967 
968  for (int j = 0; j < evSize; ++j) {
969  const CbmMuchPoint* pMuchPt =
970  static_cast<const CbmMuchPoint*>(fMuchMCPoints->Get(0, i, j));
971  PointDataHolder muchPt;
972  muchPt.x = (pMuchPt->GetXIn() + pMuchPt->GetXOut()) / 2;
973  muchPt.y = (pMuchPt->GetYIn() + pMuchPt->GetYOut()) / 2;
974  muchPt.t = fEventTimes[i] + pMuchPt->GetTime();
975  muchPt.eventId = i;
976  muchPt.trackId = pMuchPt->GetTrackID();
977  muchPt.pointId = j;
978  muchPt.stationNumber =
979  CbmMuchGeoScheme::GetStationIndex(pMuchPt->GetDetectorId());
980  muchPt.layerNumber =
981  CbmMuchGeoScheme::GetLayerIndex(pMuchPt->GetDetectorId());
982  evPoints.push_back(muchPt);
983  fMCTracks[muchPt.eventId][muchPt.trackId]
984  .pointInds[muchPt.stationNumber][muchPt.layerNumber] = muchPt.pointId;
985  }
986  }
987 
988  for (vector<vector<TrackDataHolder>>::iterator i = fMCTracks.begin();
989  i != fMCTracks.end();
990  ++i) {
991  vector<TrackDataHolder>& evTracks = *i;
992 
993  for (vector<TrackDataHolder>::iterator j = evTracks.begin();
994  j != evTracks.end();
995  ++j) {
996  TrackDataHolder& track = *j;
997 
998  if (!track.isSignal) continue;
999 
1000  for (int k = 0; k < NOF_STATIONS; ++k) {
1001  for (int l = 0; l < NOF_LAYERS; ++l) {
1002  if (track.pointInds[k][l] < 0) {
1003  track.isSignal = false;
1004  break;
1005  }
1006  }
1007 
1008  if (!track.isSignal) break;
1009  }
1010  }
1011  }
1012 
1013 #ifdef LXTB_DEBUG
1014  for (int i = 0; i < NOF_STATIONS; ++i) {
1015  char buf[64];
1016  sprintf(buf, "deltaXRHisto_%d", i);
1017  deltaXRHisto[i] = new TH1F(buf, buf, 300, -15., 15.);
1018  sprintf(buf, "deltaYRHisto_%d", i);
1019  deltaYRHisto[i] = new TH1F(buf, buf, 300, -15., 15.);
1020  sprintf(buf, "deltaTRHisto_%d", i);
1021  deltaTRHisto[i] = new TH1F(buf, buf, 300, -15., 15.);
1022 
1023  sprintf(buf, "deltaXLHisto_%d", i);
1024  deltaXLHisto[i] = new TH1F(buf, buf, 300, -15., 15.);
1025  sprintf(buf, "deltaYLHisto_%d", i);
1026  deltaYLHisto[i] = new TH1F(buf, buf, 300, -15., 15.);
1027  sprintf(buf, "deltaTLHisto_%d", i);
1028  deltaTLHisto[i] = new TH1F(buf, buf, 300, -15., 15.);
1029  }
1030 #endif //LXTB_DEBUG
1031 #endif //LXTB_QA
1032 
1033  return kSUCCESS; //, kERROR, kFATAL
1034 }
1035 
1036 static Int_t currentEventN = 0;
1037 static unsigned long long tsStartTime = 0;
1038 
1039 #ifdef LXTB_QA
1040 static long fullDuration = 0;
1041 static int nofTriggerings = 0;
1042 #endif //LXTB_QA
1043 
1044 #ifdef LXTB_EMU_TS
1045 static Double_t min_ts_time = 100000;
1046 static Double_t max_ts_time = -100000;
1047 static list<LxTbBinnedPoint> ts_points;
1048 #endif //LXTB_EMU_TS
1049 
1050 #ifdef LXTB_DEBUG
1051 /*struct LxTbDebug
1052  {
1053  struct Triplet
1054  {
1055  Int_t left;
1056  Int_t middle;
1057  Int_t right;
1058 
1059  Triplet(Int_t l, Int_t m, Int_t r) : left(l), middle(m), right(r) {}
1060  };
1061 
1062  struct TrLess
1063  {
1064  bool operator()(const Triplet& a, const Triplet& b) const
1065  {
1066  if (a.left < b.left)
1067  return true;
1068  else if (a.middle < b.middle)
1069  return true;
1070  else if (a.right < b.right)
1071  return true;
1072 
1073  return false;
1074  }
1075  };
1076 
1077  map<Triplet, bool, TrLess> triplets[NOF_STATIONS];
1078  int stationNumber;
1079 
1080  explicit LxTbDebug(vector<LxTBMLFinder::TrackDataHolder>& mcTracks) : stationNumber(-1)
1081  {
1082  for (vector<LxTBMLFinder::TrackDataHolder>::const_iterator i = mcTracks.begin(); i != mcTracks.end(); ++i)
1083  {
1084  const LxTBMLFinder::TrackDataHolder& track = *i;
1085 
1086  if (!track.isSignal)
1087  continue;
1088 
1089  for (int j = 0; j < NOF_STATIONS; ++j)
1090  {
1091  Triplet trip(track.pointInds[j][0], track.pointInds[j][1], track.pointInds[j][2]);
1092  triplets[j][trip] = false;
1093  }
1094  }
1095  }
1096 
1097  void operator()(const LxTbBinnedPoint& point)
1098  {
1099  if (!point.use)
1100  return;
1101 
1102  for (list<LxTbBinnedTriplet*>::const_iterator i = point.triplets.begin(); i != point.triplets.end(); ++i)
1103  {
1104  LxTbBinnedTriplet* trip = *i;
1105 
1106  for (list<LxTbBinnedPoint::PointDesc>::const_iterator j = point.mcRefs.begin(); j != point.mcRefs.end(); ++j)
1107  {
1108  Int_t mId = j->pointId;
1109 
1110  for (list<LxTbBinnedPoint::PointDesc>::const_iterator k = trip->rPoint->mcRefs.begin(); k != trip->rPoint->mcRefs.end(); ++k)
1111  {
1112  Int_t rId = k->pointId;
1113 
1114  for (list<LxTbBinnedPoint::PointDesc>::const_iterator l = trip->lPoint->mcRefs.begin(); l != trip->lPoint->mcRefs.end(); ++l)
1115  {
1116  Int_t lId = l->pointId;
1117  Triplet mcTrip(lId, mId, rId);
1118  map<Triplet, bool, TrLess>::iterator iter = triplets[stationNumber].find(mcTrip);
1119 
1120  if (iter != triplets[stationNumber].end())
1121  {
1122  iter->second = true;
1123  cout << stationNumber << " " << iter->first.left << " " << iter->first.middle << " " << iter->first.right << " " << iter->second << endl;
1124  }
1125  }
1126  }
1127  }
1128  }
1129  }
1130  };*/
1131 #endif //LXTB_DEBUG
1132 
1133 void LxTBMLFinder::Exec(Option_t* opt) {
1134 #ifdef LXTB_DEBUG
1135  debugTracks2.clear();
1136 #endif //LXTB_DEBUG
1137  LxTbDetector* pReconstructor = static_cast<LxTbDetector*>(fReconstructor);
1138  pReconstructor->Clear();
1139  pReconstructor->SetMinT(tsStartTime);
1140 
1141  for (int i = 0; i < fMuchPixelHits->GetEntriesFast(); ++i) {
1142  const CbmMuchPixelHit* hit =
1143  static_cast<const CbmMuchPixelHit*>(fMuchPixelHits->At(i));
1144  Int_t stationNumber = CbmMuchGeoScheme::GetStationIndex(hit->GetAddress());
1145  Int_t layerNumber = CbmMuchGeoScheme::GetLayerIndex(hit->GetAddress());
1146  scaltype x = hit->GetX();
1147  scaltype y = hit->GetY();
1148  timetype t = hit->GetTime();
1149  scaltype dx = hit->GetDx();
1150  scaltype dy = hit->GetDy();
1151  timetype dt = 4; //hit->GetTimeError();
1152  LxTbBinnedPoint point(
1153  x, dx, y, dy, t, dt, i, LAST_STATION == stationNumber);
1154 #ifdef LXTB_QA
1155  point.isTrd = false;
1156  point.stationNumber = stationNumber;
1157  point.layerNumber = layerNumber;
1158  Int_t clusterId = hit->GetRefId();
1159  const CbmCluster* cluster =
1160  static_cast<const CbmCluster*>(fMuchClusters->At(clusterId));
1161  Int_t nDigis = cluster->GetNofDigis();
1162  double avT = 0;
1163 #ifdef LXTB_EMU_TS
1164  double avTErr = 0;
1165 #endif //LXTB_EMU_TS
1166  int nofT = 0;
1167 
1168  for (Int_t j = 0; j < nDigis; ++j) {
1169  const CbmMatch* digiMatch = static_cast<const CbmMatch*>(
1170  fMuchPixelDigiMatches->At(cluster->GetDigi(j)));
1171  Int_t nMCs = digiMatch->GetNofLinks();
1172 
1173  for (Int_t k = 0; k < nMCs; ++k) {
1174  const CbmLink& lnk = digiMatch->GetLink(k);
1175  Int_t eventId = fIsEvByEv ? currentEventN : lnk.GetEntry();
1176  Int_t pointId = lnk.GetIndex();
1177  const FairMCPoint* pMCPt = static_cast<const FairMCPoint*>(
1178  fMuchMCPoints->Get(0, eventId, pointId));
1179  Int_t trackId = pMCPt->GetTrackID();
1180  LxTbBinnedPoint::PointDesc ptDesc = {eventId, pointId, trackId};
1181  point.mcRefs.push_back(ptDesc);
1182  Double_t deltaT = fMuchPoints[eventId][pointId].t;
1183 #ifdef LXTB_EMU_TS
1184  deltaT += gRandom->Gaus(0, 4);
1185  avTErr += 4 * 4;
1186 #endif //LXTB_EMU_TS
1187  avT += deltaT;
1188  ++nofT;
1189  }
1190  }
1191 
1192  if (nofT > 0) {
1193  avT /= nofT;
1194 #ifdef LXTB_EMU_TS
1195  avTErr = TMath::Sqrt(avTErr);
1196  avTErr /= nofT;
1197  dt = avT;
1198 #endif //LXTB_EMU_TS
1199  }
1200 
1201  t = avT;
1202 #endif //LXTB_QA
1203  point.t = t;
1204  //point.t = tsStartTime + hit->GetTime();
1205  //point.dt = hit->GetTimeError();
1206 
1207 #ifdef LXTB_EMU_TS
1208  ts_points.push_back(point);
1209 
1210  if (min_ts_time > t) min_ts_time = t;
1211 
1212  if (max_ts_time < t) max_ts_time = t;
1213 #else //LXTB_EMU_TS
1214 
1215 #ifdef LXTB_DEBUG
1216  int qq = 0;
1217 
1218  if (LAST_STATION == stationNumber && LAST_LAYER == layerNumber) qq += 10;
1219 #endif //LXTB_DEBUG
1220  scaltype minY =
1221  pReconstructor->fStations[stationNumber].fLayers[layerNumber].minY;
1222  scaltype binSizeY =
1223  pReconstructor->fStations[stationNumber].fLayers[layerNumber].yBinLength;
1224  int lastYBin = (pReconstructor->fStations[stationNumber]
1225  .fLayers[layerNumber]
1226  .lastYBinNumber);
1227  scaltype minX =
1228  pReconstructor->fStations[stationNumber].fLayers[layerNumber].minX;
1229  scaltype binSizeX =
1230  pReconstructor->fStations[stationNumber].fLayers[layerNumber].xBinLength;
1231  int lastXBin = pReconstructor->fStations[stationNumber]
1232  .fLayers[layerNumber]
1233  .lastXBinNumber;
1234  int last_timebin = pReconstructor->fStations[stationNumber]
1235  .fLayers[layerNumber]
1236  .lastTimeBinNumber;
1237 
1238  int tInd =
1239  (t - pReconstructor->fStations[stationNumber].fLayers[layerNumber].minT)
1240  / TIMEBIN_LENGTH;
1241 
1242  if (tInd < 0)
1243  tInd = 0;
1244  else if (tInd > last_timebin)
1245  tInd = last_timebin;
1246 
1247  LxTbTYXBin& tyxBin = pReconstructor->fStations[stationNumber]
1248  .fLayers[layerNumber]
1249  .tyxBins[tInd];
1250  int yInd = (y - minY) / binSizeY;
1251 
1252  if (yInd < 0)
1253  yInd = 0;
1254  else if (yInd > lastYBin)
1255  yInd = lastYBin;
1256 
1257  LxTbYXBin& yxBin = tyxBin.yxBins[yInd];
1258  int xInd = (x - minX) / binSizeX;
1259 
1260  if (xInd < 0)
1261  xInd = 0;
1262  else if (xInd > lastXBin)
1263  xInd = lastXBin;
1264 
1265  LxTbXBin& xBin = yxBin.xBins[xInd];
1266  xBin.points.push_back(point);
1267 
1268  if (LAST_STATION == stationNumber) {
1269  xBin.use = true;
1270  yxBin.use = true;
1271  tyxBin.use = true;
1272 
1273  if (0 == layerNumber)
1274  xBin.use = true;
1275  else if (1 == layerNumber)
1276  xBin.use = true;
1277  else if (2 == layerNumber)
1278  xBin.use = true;
1279  }
1280 #endif //LXTB_EMU_TS
1281  } // for (int i = 0; i < fMuchPixelHits->GetEntriesFast(); ++i)
1282 
1283 #ifdef LXTB_DEBUG
1284  struct DebugTrack {
1285  list<Int_t> pointIds[NOF_STATIONS][NOF_LAYERS];
1286  };
1287 
1288  map<Int_t, DebugTrack> debugTracks;
1289 
1290  for (int i = 0; i < NOF_STATIONS; ++i) {
1291  for (int j = 0; j < NOF_LAYERS; ++j) {
1292  const LxTbLayer& layer = pReconstructor->fStations[i].fLayers[j];
1293 
1294  for (int k = 0; k < layer.nofTYXBins; ++k) {
1295  LxTbTYXBin& tyxBin = layer.tyxBins[k];
1296 
1297  for (int l = 0; l < layer.nofYXBins; ++l) {
1298  LxTbYXBin& yxBin = tyxBin.yxBins[l];
1299 
1300  for (int m = 0; m < layer.nofXBins; ++m) {
1301  LxTbXBin& xBin = yxBin.xBins[m];
1302 
1303  for (std::list<LxTbBinnedPoint>::iterator n = xBin.points.begin();
1304  n != xBin.points.end();
1305  ++n) {
1306  LxTbBinnedPoint& point = *n;
1307 
1308  for (list<LxTbBinnedPoint::PointDesc>::const_iterator o =
1309  point.mcRefs.begin();
1310  o != point.mcRefs.end();
1311  ++o) {
1312  const LxTbBinnedPoint::PointDesc& pd = *o;
1313  debugTracks[pd.trackId]
1314  .pointIds[point.stationNumber][point.layerNumber]
1315  .push_back(pd.pointId);
1316  debugTracks2[pd.trackId]
1317  .points[point.stationNumber][point.layerNumber]
1318  .push_back(&point);
1319  }
1320  }
1321  }
1322  }
1323  }
1324  }
1325  }
1326 #endif //LXTB_DEBUG
1327 
1328  timespec ts;
1329  clock_gettime(CLOCK_REALTIME, &ts);
1330  long beginTime = ts.tv_sec * 1000000000 + ts.tv_nsec;
1331 
1332  pReconstructor->Reconstruct();
1333 
1334  clock_gettime(CLOCK_REALTIME, &ts);
1335  long endTime = ts.tv_sec * 1000000000 + ts.tv_nsec;
1336  fullDuration += endTime - beginTime;
1337 
1338 #ifdef LXTB_DEBUG
1339  static int nofTriplesAll[NOF_STATIONS] = {0, 0, 0, 0};
1340  static int nofTriplesFound[NOF_STATIONS] = {0, 0, 0, 0};
1341 
1342  struct Debug {
1343  struct Triplet {
1344  Int_t left;
1345  Int_t middle;
1346  Int_t right;
1347 
1348  Triplet(Int_t l, Int_t m, Int_t r) : left(l), middle(m), right(r) {}
1349  };
1350 
1351  struct TrLess {
1352  bool operator()(const Triplet& a, const Triplet& b) const {
1353  if (a.left < b.left)
1354  return true;
1355  else if (a.middle < b.middle)
1356  return true;
1357  else if (a.right < b.right)
1358  return true;
1359 
1360  return false;
1361  }
1362  };
1363 
1364  map<Triplet, bool, TrLess> triplets[NOF_STATIONS];
1365  int stationNumber;
1366 
1367  explicit Debug(vector<TrackDataHolder>& mcTracks) : stationNumber(-1) {
1368  for (vector<TrackDataHolder>::const_iterator i = mcTracks.begin();
1369  i != mcTracks.end();
1370  ++i) {
1371  const TrackDataHolder& track = *i;
1372 
1373  if (!track.isSignal) continue;
1374 
1375  for (int j = 0; j < NOF_STATIONS; ++j) {
1376  Triplet trip(track.pointInds[j][0],
1377  track.pointInds[j][1],
1378  track.pointInds[j][2]);
1379  triplets[j][trip] = false;
1380  }
1381  }
1382  }
1383 
1384  void operator()(const LxTbBinnedPoint& point) {
1385  if (!point.use) return;
1386 
1387  for (list<LxTbBinnedTriplet*>::const_iterator i = point.triplets.begin();
1388  i != point.triplets.end();
1389  ++i) {
1390  LxTbBinnedTriplet* trip = *i;
1391 
1392  for (list<LxTbBinnedPoint::PointDesc>::const_iterator j =
1393  point.mcRefs.begin();
1394  j != point.mcRefs.end();
1395  ++j) {
1396  Int_t mId = j->pointId;
1397 
1398  for (list<LxTbBinnedPoint::PointDesc>::const_iterator k =
1399  trip->rPoint->mcRefs.begin();
1400  k != trip->rPoint->mcRefs.end();
1401  ++k) {
1402  Int_t rId = k->pointId;
1403 
1404  for (list<LxTbBinnedPoint::PointDesc>::const_iterator l =
1405  trip->lPoint->mcRefs.begin();
1406  l != trip->lPoint->mcRefs.end();
1407  ++l) {
1408  Int_t lId = l->pointId;
1409  Triplet mcTrip(lId, mId, rId);
1410  map<Triplet, bool, TrLess>::iterator iter =
1411  triplets[stationNumber].find(mcTrip);
1412 
1413  if (iter != triplets[stationNumber].end()) iter->second = true;
1414  }
1415  }
1416  }
1417  }
1418  }
1419  };
1420 
1421  Debug debug(fMCTracks[currentEventN]);
1422 
1423  /*for (int i = 0; i < fMuchPixelHits->GetEntriesFast(); ++i)
1424  {
1425  const CbmMuchPixelHit* hit = static_cast<const CbmMuchPixelHit*> (fMuchPixelHits->At(i));
1426  Int_t stationNumber = CbmMuchGeoScheme::GetStationIndex(hit->GetAddress());
1427  Int_t layerNumber = CbmMuchGeoScheme::GetLayerIndex(hit->GetAddress());
1428  Int_t clusterId = hit->GetRefId();
1429  const CbmCluster* cluster = static_cast<const CbmCluster*> (fMuchClusters->At(clusterId));
1430  Int_t nDigis = cluster->GetNofDigis();
1431 
1432  for (Int_t j = 0; j < nDigis; ++j)
1433  {
1434  const CbmMatch* digiMatch = static_cast<const CbmMatch*> (fMuchPixelDigiMatches->At(cluster->GetDigi(j)));
1435  Int_t nMCs = digiMatch->GetNofLinks();
1436 
1437  for (Int_t k = 0; k < nMCs; ++k)
1438  {
1439  const CbmLink& lnk = digiMatch->GetLink(k);
1440  Int_t pointId = lnk.GetIndex();
1441  const FairMCPoint* pMCPt = static_cast<const FairMCPoint*> (fMuchMCPoints->Get(0, currentEventN, pointId));
1442  Int_t trackId = pMCPt->GetTrackID();
1443  debugTracks[trackId].pointIds[stationNumber][layerNumber].push_back(pointId);
1444  }
1445  }
1446  }*/
1447 
1448  /*for (map<Int_t, DebugTrack>::const_iterator i = debugTracks.begin(); i != debugTracks.end(); ++i)
1449  {
1450  const DebugTrack& dt = i->second;
1451 
1452  for (int j = 0; j < NOF_STATIONS; ++j)
1453  {
1454  for (list<Int_t>::const_iterator k = dt.pointIds[j][1].begin(); k != dt.pointIds[j][1].end(); ++k)
1455  {
1456  Int_t mId = *k;
1457 
1458  for (list<Int_t>::const_iterator l = dt.pointIds[j][2].begin(); l != dt.pointIds[j][2].end(); ++l)
1459  {
1460  Int_t rId = *l;
1461 
1462  for (list<Int_t>::const_iterator m = dt.pointIds[j][0].begin(); m != dt.pointIds[j][0].end(); ++m)
1463  {
1464  Int_t lId = *m;
1465  Debug::Triplet mcTrip(lId, mId, rId);
1466  map<Debug::Triplet, bool, Debug::TrLess>::iterator iter = debug.triplets[j].find(mcTrip);
1467 
1468  if (iter != debug.triplets[j].end())
1469  iter->second = true;
1470  }
1471  }
1472  }
1473  }
1474  }*/
1475 
1476  for (map<Int_t, DebugTrack2>::const_iterator i = debugTracks2.begin();
1477  i != debugTracks2.end();
1478  ++i) {
1479  const DebugTrack2& dt = i->second;
1480 
1481  for (int j = 0; j < NOF_STATIONS; ++j) {
1482  scaltype deltaZr = pReconstructor->fStations[j].fLayers[2].z
1483  - pReconstructor->fStations[j].fLayers[1].z;
1484  scaltype deltaZl = pReconstructor->fStations[j].fLayers[0].z
1485  - pReconstructor->fStations[j].fLayers[1].z;
1486 
1487  for (list<LxTbBinnedPoint*>::const_iterator k = dt.points[j][1].begin();
1488  k != dt.points[j][1].end();
1489  ++k) {
1490  LxTbBinnedPoint* mPt = *k;
1491  scaltype rTx = mPt->x / pReconstructor->fStations[j].fLayers[1].z;
1492  scaltype rTy = mPt->y / pReconstructor->fStations[j].fLayers[1].z;
1493  scaltype rX = mPt->x + rTx * deltaZr;
1494  scaltype rY = mPt->y + rTy * deltaZr;
1495  scaltype trajLenR = sqrt(1 + rTx * rTx + rTy * rTy) * deltaZr;
1496  timetype rT =
1497  mPt->t + 1.e9 * trajLenR / speedOfLight; // 1.e9 to convert to ns.
1498 
1499  for (list<LxTbBinnedPoint*>::const_iterator l = dt.points[j][2].begin();
1500  l != dt.points[j][2].end();
1501  ++l) {
1502  LxTbBinnedPoint* rPt = *l;
1503  scaltype wXr =
1504  NOF_SIGMAS
1505  * sqrt(pReconstructor->fStations[j].fHandleMPoint.scatXRL
1506  * pReconstructor->fStations[j].fHandleMPoint.scatXRL
1507  + mPt->dx * mPt->dx + rPt->dx * rPt->dx);
1508  scaltype wYr = NOF_SIGMAS
1509  * sqrt(pReconstructor->fStations[j].fScatYRL
1510  * pReconstructor->fStations[j].fScatYRL
1511  + mPt->dy * mPt->dy + rPt->dy * rPt->dy);
1512  scaltype wTr =
1513  NOF_SIGMAS * sqrt(mPt->dt * mPt->dt + rPt->dt * rPt->dt);
1514  deltaXRHisto[j]->Fill(rPt->x - rX);
1515  deltaYRHisto[j]->Fill(rPt->y - rY);
1516  deltaTRHisto[j]->Fill(rPt->t - rT);
1517 
1518  if (LAST_STATION == j) {
1519  if (abs(rPt->x - rX) > wXr) ++nofOverXR;
1520 
1521  if (abs(rPt->y - rY) > wYr) ++nofOverYR;
1522 
1523  if (abs(rPt->t - rT) > wTr) ++nofOverTR;
1524  }
1525 
1526  for (list<LxTbBinnedPoint*>::const_iterator m =
1527  dt.points[j][0].begin();
1528  m != dt.points[j][0].end();
1529  ++m) {
1530  LxTbBinnedPoint* lPt = *m;
1531  scaltype wXl =
1532  NOF_SIGMAS
1533  * sqrt(pReconstructor->fStations[j].fHandleRPoint.scatXLL
1534  * pReconstructor->fStations[j].fHandleRPoint.scatXLL
1535  + mPt->dx * mPt->dx + rPt->dx * rPt->dx
1536  + lPt->dx * lPt->dx);
1537  scaltype wYl =
1538  NOF_SIGMAS
1539  * sqrt(mPt->dy * mPt->dy + rPt->dy * rPt->dy + lPt->dy * lPt->dy);
1540  scaltype wTl =
1541  NOF_SIGMAS * sqrt(mPt->dt * mPt->dt + lPt->dt * lPt->dt);
1542  scaltype lTx = (rPt->x - mPt->x)
1543  / (pReconstructor->fStations[j].fLayers[2].z
1544  - pReconstructor->fStations[j].fLayers[1].z);
1545  scaltype lTy = (rPt->y - mPt->y)
1546  / (pReconstructor->fStations[j].fLayers[2].z
1547  - pReconstructor->fStations[j].fLayers[1].z);
1548  scaltype lX = mPt->x + lTx * deltaZl;
1549  scaltype lY = mPt->y + lTy * deltaZl;
1550  scaltype trajLenL = sqrt(1 + lTx * lTx + lTy * lTy) * deltaZl;
1551  timetype lT =
1552  mPt->t
1553  + 1.e9 * trajLenL
1554  / speedOfLight; // 1.e9 to convert to ns and trajLenL is negative.
1555  deltaXLHisto[j]->Fill(lPt->x - lX);
1556  deltaYLHisto[j]->Fill(lPt->y - lY);
1557  deltaTLHisto[j]->Fill(lPt->t - lT);
1558 
1559  if (LAST_STATION == j) {
1560  if (abs(lPt->x - lX) > wXl) ++nofOverXL;
1561 
1562  if (abs(lPt->y - lY) > wYl) ++nofOverYL;
1563 
1564  if (abs(lPt->t - lT) > wTl) ++nofOverTL;
1565  }
1566  }
1567  }
1568  }
1569  }
1570  }
1571 
1572  cout << "nofOverXR = " << nofOverXR << endl;
1573  cout << "nofOverYR = " << nofOverYR << endl;
1574  cout << "nofOverTR = " << nofOverTR << endl;
1575 
1576  cout << "nofOverXL = " << nofOverXL << endl;
1577  cout << "nofOverYL = " << nofOverYL << endl;
1578  cout << "nofOverTL = " << nofOverTL << endl;
1579 
1580  cout << "nofFoundR = " << nofFoundR << endl;
1581  cout << "nofNotFoundR = " << nofNotFoundR << endl;
1582 
1583  for (int i = LAST_STATION; i >= 0; --i) {
1584  debug.stationNumber = i;
1585  IterateLayer(pReconstructor->fStations[i].fLayers[1], debug);
1586 
1587  nofTriplesAll[i] += debug.triplets[i].size();
1588 
1589  for (map<Debug::Triplet, bool, Debug::TrLess>::const_iterator j =
1590  debug.triplets[i].begin();
1591  j != debug.triplets[i].end();
1592  ++j) {
1593  if (j->second) ++nofTriplesFound[i];
1594  }
1595 
1596  double foundPerc = 100 * nofTriplesFound[i];
1597  foundPerc /= nofTriplesAll[i];
1598  cout << "For station #:" << i << " found triplets: " << foundPerc << " ( "
1599  << nofTriplesFound[i] << " / " << nofTriplesAll[i] << " )" << endl;
1600  }
1601 #endif //LXTB_DEBUG
1602 
1603  cout << "In event #" << currentEventN
1604  << " reconstructed: " << pReconstructor->recoTracks.size() << " tracks"
1605  << endl;
1606 
1607  static int nofSignals = 0;
1608  static int nofRecoSignals = 0;
1609  bool hasPos = false;
1610  bool hasNeg = false;
1611 
1612  for (vector<TrackDataHolder>::const_iterator i =
1613  fMCTracks[currentEventN].begin();
1614  i != fMCTracks[currentEventN].end();
1615  ++i) {
1616  const TrackDataHolder& dh = *i;
1617 
1618  if (!dh.isSignal) continue;
1619 
1620  if (dh.isPos)
1621  hasPos = true;
1622  else
1623  hasNeg = true;
1624  }
1625 
1626  if (hasPos && hasNeg) ++nofSignals;
1627 
1628  bool hasPos2 = false;
1629  bool hasNeg2 = false;
1630 
1631  for (list<Chain*>::const_iterator i = pReconstructor->recoTracks.begin();
1632  i != pReconstructor->recoTracks.end();
1633  ++i) {
1634  const Chain* track = *i;
1635 
1636  if ((track->points[0][2]->x - track->points[0][0]->x)
1637  / (pReconstructor->fStations[0].fLayers[2].z
1638  - pReconstructor->fStations[0].fLayers[0].z)
1639  - track->points[0][1]->x / pReconstructor->fStations[0].fLayers[1].z
1640  > 0)
1641  hasPos2 = true;
1642  else
1643  hasNeg2 = true;
1644  }
1645 
1646  if (hasPos2 && hasNeg2) {
1647  ++nofTriggerings;
1648 
1649  if (hasNeg && hasPos) ++nofRecoSignals;
1650  }
1651 
1652  cout << "The number of triggerings: " << nofTriggerings << endl;
1653  cout << "The triggering efficiency: "
1654  << (0 == nofSignals ? 100 : 100 * nofRecoSignals / nofSignals) << " % ( "
1655  << nofRecoSignals << " / " << nofSignals << " )" << endl;
1656 
1657  recoTracks.splice(recoTracks.end(), pReconstructor->recoTracks);
1658  ++currentEventN;
1659  tsStartTime += 100;
1660 }
1661 
1662 struct RecoTrackData {
1663  Int_t eventId;
1664  Int_t trackId;
1665 
1666  RecoTrackData(Int_t eId, Int_t tId) : eventId(eId), trackId(tId) {}
1667 };
1668 
1669 struct RTDLess {
1670  bool operator()(const RecoTrackData& x, const RecoTrackData& y) const {
1671  if (x.eventId < y.eventId) return true;
1672 
1673  return x.trackId < y.trackId;
1674  }
1675 };
1676 
1677 static void SaveHisto(TH1* histo) {
1678  TFile* curFile = TFile::CurrentFile();
1679  TString histoName = histo->GetName();
1680  histoName += ".root";
1681  TFile fh(histoName.Data(), "RECREATE");
1682  histo->Write();
1683  fh.Close();
1684  delete histo;
1685  TFile::CurrentFile() = curFile;
1686 }
1687 
1688 void LxTBMLFinder::Finish() {
1689 #ifdef LXTB_EMU_TS
1690  Double_t tCoeff = TIMEBIN_LENGTH * nof_timebins / (max_ts_time - min_ts_time);
1691 
1692  for (list<LxTbBinnedPoint>::iterator i = ts_points.begin();
1693  i != ts_points.end();
1694  ++i) {
1695  LxTbBinnedPoint& point = *i;
1696  point.t = (point.t - min_ts_time) * tCoeff;
1697  point.dt *= tCoeff;
1698 
1699  bool isTrd = point.isTrd;
1700  Int_t stationNumber = point.stationNumber;
1701  scaltype minY = (isTrd ? fFinder->trdStation.minY
1702  : fFinder->stations[stationNumber].minY);
1703  scaltype binSizeY = (isTrd ? fFinder->trdStation.binSizeY
1704  : fFinder->stations[stationNumber].binSizeY);
1705  int lastYBin = (isTrd ? fFinder->trdStation.lastYBin
1706  : fFinder->stations[stationNumber].lastYBin);
1707  scaltype minX = (isTrd ? fFinder->trdStation.minX
1708  : fFinder->stations[stationNumber].minX);
1709  scaltype binSizeX = (isTrd ? fFinder->trdStation.binSizeX
1710  : fFinder->stations[stationNumber].binSizeX);
1711  int lastXBin = (isTrd ? fFinder->trdStation.lastXBin
1712  : fFinder->stations[stationNumber].lastXBin);
1713 
1714  int tInd = (point.t - fFinder->minT) / CUR_TIMEBIN_LENGTH;
1715 
1716  if (tInd < 0)
1717  tInd = 0;
1718  else if (tInd > last_timebin)
1719  tInd = last_timebin;
1720 
1721  LxTbTYXBin& tyxBin =
1722  (isTrd ? fFinder->trdStation.tyxBinsArr[stationNumber][tInd]
1723  : fFinder->stations[stationNumber].tyxBins[tInd]);
1724  int yInd = (point.y - minY) / binSizeY;
1725 
1726  if (yInd < 0)
1727  yInd = 0;
1728  else if (yInd > lastYBin)
1729  yInd = lastYBin;
1730 
1731  LxTbYXBin& yxBin = tyxBin.yxBins[yInd];
1732  int xInd = (point.x - minX) / binSizeX;
1733 
1734  if (xInd < 0)
1735  xInd = 0;
1736  else if (xInd > lastXBin)
1737  xInd = lastXBin;
1738 
1739  LxTbXBin& xBin = yxBin.xBins[xInd];
1740  xBin.points.push_back(point);
1741 
1742  if ( CUR_LAST_STATION == stationNumber) {
1743  xBin.use = true;
1744  yxBin.use = true;
1745  tyxBin.use = true;
1746  }
1747  }
1748 
1749 #ifdef LXTB_TIE
1750  fDetector->SetTSBegin(0);
1751 #endif //LXTB_TIE
1752 
1753  fFinder->Reconstruct();
1754 
1755  for (list<timetype>::iterator i = shortSignalMCTimes.begin();
1756  i != shortSignalMCTimes.end();
1757  ++i) {
1758  timetype& v = *i;
1759  v = (v - min_ts_time) * tCoeff;
1760  }
1761 
1762  for (list<timetype>::iterator i = longSignalMCTimes.begin();
1763  i != longSignalMCTimes.end();
1764  ++i) {
1765  timetype& v = *i;
1766  v = (v - min_ts_time) * tCoeff;
1767  }
1768 
1769  for (int i = 0; i < fFinder->nofTrackBins; ++i) {
1770  list<LxTbBinnedFinder::Chain*>& recoTracksBin = fFinder->recoTracks[i];
1771 
1772  for (list<LxTbBinnedFinder::Chain*>::const_iterator j =
1773  recoTracksBin.begin();
1774  j != recoTracksBin.end();
1775  ++j)
1776  recoTracks.push_back(*j);
1777  }
1778 
1779  SpliceTriggerings(triggerTimes_trd0_sign0_dist0,
1780  fFinder->triggerTimes_trd0_sign0_dist0);
1781  SpliceTriggerings(triggerTimes_trd0_sign0_dist1,
1782  fFinder->triggerTimes_trd0_sign0_dist1);
1783  SpliceTriggerings(triggerTimes_trd0_sign1_dist0,
1784  fFinder->triggerTimes_trd0_sign1_dist0);
1785  SpliceTriggerings(triggerTimes_trd0_sign1_dist1,
1786  fFinder->triggerTimes_trd0_sign1_dist1);
1787  SpliceTriggerings(triggerTimes_trd1_sign0_dist0,
1788  fFinder->triggerTimes_trd1_sign0_dist0);
1789  SpliceTriggerings(triggerTimes_trd1_sign0_dist1,
1790  fFinder->triggerTimes_trd1_sign0_dist1);
1791  SpliceTriggerings(triggerTimes_trd1_sign1_dist0,
1792  fFinder->triggerTimes_trd1_sign1_dist0);
1793  SpliceTriggerings(triggerTimes_trd1_sign1_dist1,
1794  fFinder->triggerTimes_trd1_sign1_dist1);
1795 #endif //LXTB_EMU_TS
1796  cout << "LxTbBinnedFinder::Reconstruct() full duration was: " << fullDuration
1797  << endl;
1798 
1799  int nofRecoTracks = recoTracks.size();
1800  cout << "LxTbBinnedFinder::Reconstruct() the number of found tracks: "
1801  << nofRecoTracks << endl;
1802 
1803 #ifdef LXTB_QA
1804  static int nofSignalTracks = 0;
1805  static int nofRecoSignalTracks = 0;
1806  int eventN = 0;
1807 
1808  for (vector<vector<TrackDataHolder>>::iterator i = fMCTracks.begin();
1809  i != fMCTracks.end();
1810  ++i) {
1811  vector<TrackDataHolder>& evTracks = *i;
1812 
1813  for (vector<TrackDataHolder>::iterator j = evTracks.begin();
1814  j != evTracks.end();
1815  ++j) {
1816  TrackDataHolder& track = *j;
1817 
1818  if (!track.isSignal) continue;
1819 
1820  ++nofSignalTracks;
1821 
1822  int nofMatchPoints = 0;
1823 
1824  for (list<Chain*>::const_iterator k = recoTracks.begin();
1825  k != recoTracks.end();
1826  ++k) {
1827  const Chain* chain = *k;
1828 
1829  for (int l = 0; l < NOF_STATIONS; ++l) {
1830  for (int m = 0; m < NOF_LAYERS; ++m) {
1831  bool pointsMatched = false;
1832 
1833  for (list<LxTbBinnedPoint::PointDesc>::const_iterator n =
1834  chain->points[l][m]->mcRefs.begin();
1835  n != chain->points[l][m]->mcRefs.end();
1836  ++n) {
1837  if (n->eventId == eventN && n->pointId == track.pointInds[l][m]) {
1838  pointsMatched = true;
1839  break;
1840  }
1841  }
1842 
1843  if (pointsMatched) ++nofMatchPoints;
1844  }
1845  }
1846  }
1847 
1848  if (nofMatchPoints >= NOF_STATIONS * NOF_LAYERS * 0.7) {
1849  ++nofRecoSignalTracks;
1850  continue;
1851  }
1852  }
1853 
1854  ++eventN;
1855  }
1856 
1857  double eff =
1858  0 == nofSignalTracks ? 100 : 100.0 * nofRecoSignalTracks / nofSignalTracks;
1859  cout << "Reconstruction efficiency is: " << eff << "% [ "
1860  << nofRecoSignalTracks << " / " << nofSignalTracks << " ]" << endl;
1861 
1862  int nofRightRecoTracks = 0;
1863 
1864  for (list<Chain*>::const_iterator i = recoTracks.begin();
1865  i != recoTracks.end();
1866  ++i) {
1867  const Chain* chain = *i;
1868  map<RecoTrackData, int, RTDLess> nofTracks;
1869 
1870  for (int j = 0; j < NOF_STATIONS; ++j) {
1871  for (int k = 0; k < NOF_LAYERS; ++k) {
1872  int stMask = 1 << j * NOF_LAYERS + k;
1873 
1874  for (list<LxTbBinnedPoint::PointDesc>::const_iterator l =
1875  chain->points[j][k]->mcRefs.begin();
1876  l != chain->points[j][k]->mcRefs.end();
1877  ++l) {
1878  RecoTrackData st(l->eventId, l->trackId);
1879  map<RecoTrackData, int, RTDLess>::iterator nofTIter =
1880  nofTracks.find(st);
1881 
1882  if (nofTIter != nofTracks.end())
1883  nofTIter->second |= stMask;
1884  else
1885  nofTracks[st] = stMask;
1886  }
1887  }
1888  }
1889 
1890  int nofPoints = 0;
1891 
1892  for (map<RecoTrackData, int, RTDLess>::const_iterator j = nofTracks.begin();
1893  j != nofTracks.end();
1894  ++j) {
1895  int nofp = 0;
1896 
1897  for (int k = 0; k < NOF_STATIONS; ++k) {
1898  for (int l = 0; l < NOF_LAYERS; ++l) {
1899  if (j->second & (1 << k * NOF_LAYERS + l)) ++nofp;
1900  }
1901  }
1902 
1903  if (nofp > nofPoints) nofPoints = nofp;
1904  }
1905 
1906  if (nofPoints >= NOF_STATIONS * NOF_LAYERS * 0.7) ++nofRightRecoTracks;
1907  }
1908 
1909  eff =
1910  0 == recoTracks.size() ? 100 : 100.0 * nofRightRecoTracks / nofRecoTracks;
1911  cout << "Non ghosts are: " << eff << "% [ " << nofRightRecoTracks << " / "
1912  << nofRecoTracks << " ]" << endl;
1913 
1914  //cout << "Have: " << shortSignalMCTimes.size() << " short signaling events" << endl;
1915  //cout << "Have: " << longSignalMCTimes.size() << " long signaling events" << endl;
1916  cout << "Have: " << currentEventN << " events" << endl;
1917 
1918  ofstream nofEventsFile("nof_events.txt", ios_base::out | ios_base::trunc);
1919  nofEventsFile << currentEventN;
1920 
1921  ofstream nofTriggeringsFile("nof_triggerings.txt",
1922  ios_base::out | ios_base::trunc);
1923  nofTriggeringsFile << nofTriggerings;
1924 
1925  /*ofstream nofShortSignalsFile("nof_short_signals.txt", ios_base::out | ios_base::trunc);
1926  nofShortSignalsFile << shortSignalMCTimes.size();
1927 
1928  ofstream nofLongSignalsFile("nof_long_signals.txt", ios_base::out | ios_base::trunc);
1929  nofLongSignalsFile << longSignalMCTimes.size();
1930 
1931  PrintTrigger(triggerTimes_trd0_sign0_dist0, shortSignalMCTimes, "triggerTimes_trd0_sign0_dist0");
1932  PrintTrigger(triggerTimes_trd0_sign0_dist1, shortSignalMCTimes, "triggerTimes_trd0_sign0_dist1");
1933  PrintTrigger(triggerTimes_trd0_sign1_dist0, shortSignalMCTimes, "triggerTimes_trd0_sign1_dist0");
1934  PrintTrigger(triggerTimes_trd0_sign1_dist1, shortSignalMCTimes, "triggerTimes_trd0_sign1_dist1");
1935  PrintTrigger(triggerTimes_trd1_sign0_dist0, longSignalMCTimes, "triggerTimes_trd1_sign0_dist0");
1936  PrintTrigger(triggerTimes_trd1_sign0_dist1, longSignalMCTimes, "triggerTimes_trd1_sign0_dist1");
1937  PrintTrigger(triggerTimes_trd1_sign1_dist0, longSignalMCTimes, "triggerTimes_trd1_sign1_dist0");
1938  PrintTrigger(triggerTimes_trd1_sign1_dist1, longSignalMCTimes, "triggerTimes_trd1_sign1_dist1", true);
1939 
1940  Int_t nofTriggerDigis = 0;
1941 
1942  for (set<Int_t>::const_iterator i = fFinder->triggerEventNumber.begin(); i != fFinder->triggerEventNumber.end(); ++i)
1943  nofTriggerDigis += nof_ev_digis[*i];
1944 
1945  ofstream nofTriggerDigisFile("nof_trigger_digis.txt", ios_base::out | ios_base::trunc);
1946  nofTriggerDigisFile << nofTriggerDigis;
1947  ofstream nofDigisFile("nof_digis.txt", ios_base::out | ios_base::trunc);
1948  nofDigisFile << nof_digis;*/
1949 
1950 #ifdef LXTB_DEBUG
1951  for (int i = 0; i < NOF_STATIONS; ++i) {
1952  SaveHisto(deltaXRHisto[i]);
1953  SaveHisto(deltaYRHisto[i]);
1954  SaveHisto(deltaTRHisto[i]);
1955  SaveHisto(deltaXLHisto[i]);
1956  SaveHisto(deltaYLHisto[i]);
1957  SaveHisto(deltaTLHisto[i]);
1958  }
1959 #endif //LXTB_DEBUG
1960 #endif //LXTB_QA
1961 
1962  for (list<Chain*>::iterator i = recoTracks.begin(); i != recoTracks.end();
1963  ++i)
1964  delete *i;
1965 }
LxTbMLStation::HandleRPoint::Init
void Init()
Definition: LxTBMLTask.cxx:265
LxTbDetector::HandleRPoint::lLayer
LxTbLayer * lLayer
Definition: LxTBMLTask.cxx:585
LxTBMLFinder::PointDataHolder::stationNumber
Int_t stationNumber
Definition: LxTBMLTask.h:59
CbmMuchPoint::GetXOut
Double_t GetXOut() const
Definition: CbmMuchPoint.h:73
LxTBMLFinder::PointDataHolder::pointId
Int_t pointId
Definition: LxTBMLTask.h:58
LxTbBinnedPoint::PointDesc::pointId
Int_t pointId
Definition: LxTBBinned.h:96
CbmMCTrack::GetMotherId
Int_t GetMotherId() const
Definition: CbmMCTrack.h:71
SignalParticle::fName
const char * fName
Definition: LxTBMLTask.cxx:340
LxTbDetector::HandleRPoint::HandleLPoint
Definition: LxTBMLTask.cxx:624
gMuonMass
static scaltype gMuonMass
Definition: LxTBMatEffs.h:19
SignalParticle::fPdgCode
Int_t fPdgCode
Definition: LxTBMLTask.cxx:341
CbmMCTrack::GetMass
Double_t GetMass() const
Mass of the associated particle.
Definition: CbmMCTrack.cxx:114
CbmMatch
Definition: CbmMatch.h:22
CbmMCDataManager.h
LxTbMLStation
Definition: LxTBMLTask.cxx:79
LxTBMLFinder::TrackDataHolder::isPos
bool isPos
Definition: LxTBMLTask.h:66
LxTbMLStation::HandleLPoint::station
LxTbMLStation & station
Definition: LxTBMLTask.cxx:301
LxTbDetector::fHandleLastPoint
HandleLastPoint fHandleLastPoint
Definition: LxTBMLTask.cxx:817
LxTbDetector::HandleLastPoint::KFParams::yParams
KFParamsCoord yParams
Definition: LxTBMLTask.cxx:654
LxTbMLStation::Clear
void Clear()
Definition: LxTBMLTask.cxx:130
CUR_LAST_STATION
#define CUR_LAST_STATION
Definition: LxTBTask.h:33
CbmMuchPoint
Definition: CbmMuchPoint.h:21
LxTbBinnedPoint::use
bool use
Definition: LxTBBinned.h:85
LxTbMLStation::HandleRPoint::mPoint
LxTbBinnedPoint * mPoint
Definition: LxTBMLTask.cxx:260
NOF_LAYERS
#define NOF_LAYERS
Definition: LxTBBinned2.h:22
LxTbDetector::HandleLastPoint::KFAddPointCoord
KFParamsCoord KFAddPointCoord(KFParamsCoord prevParam, scaltype m, scaltype V, scaltype &chi2, int stationNumber, int layerNumber, int coordNumber)
Definition: LxTBMLTask.cxx:668
LxTbLayer::minT
timetype minT
Definition: LxTBBinned2.h:254
LAST_LAYER
#define LAST_LAYER
Definition: LxTBBinned2.h:23
CbmPixelHit::GetX
Double_t GetX() const
Definition: CbmPixelHit.h:83
triggerTimes_trd0_sign1_dist0
static list< pair< timetype, timetype > > triggerTimes_trd0_sign1_dist0
Definition: LxTBTask.cxx:399
CbmMatch::GetLink
const CbmLink & GetLink(Int_t i) const
Definition: CbmMatch.h:35
CbmMatch::GetNofLinks
Int_t GetNofLinks() const
Definition: CbmMatch.h:38
LxTbMLStation::Q::Q21
scaltype Q21
Definition: LxTBMLTask.cxx:81
LxTBMLFinder::fMCTracks
std::vector< std::vector< TrackDataHolder > > fMCTracks
Definition: LxTBMLTask.h:104
LxTbAbsorber::Z
scaltype Z
Definition: LxTBMatEffs.h:27
LxTbLayer::Init
void Init()
Definition: LxTBBinned2.h:232
LxTbTYXBin
Definition: LxTBBinned.h:180
LxTbDetector::recoTracks
list< LxTBMLFinder::Chain * > recoTracks
Definition: LxTBMLTask.cxx:349
LxTbDetector::HandleRPoint::handleLPoint
HandleLPoint handleLPoint
Definition: LxTBMLTask.cxx:633
LxTbBinnedPoint::mcRefs
std::list< PointDesc > mcRefs
Definition: LxTBBinned.h:100
currentEventN
static Int_t currentEventN
Definition: LxTBMLTask.cxx:1036
scaltype
#define scaltype
Definition: CbmGlobalTrackingDefs.h:17
LxTbMLStation::HandleLPoint
Definition: LxTBMLTask.cxx:300
triggerTimes_trd0_sign0_dist0
static list< pair< timetype, timetype > > triggerTimes_trd0_sign0_dist0
Definition: LxTBTask.cxx:397
CbmPixelHit::GetY
Double_t GetY() const
Definition: CbmPixelHit.h:84
LxTbMLStation::fScatXLS
scaltype fScatXLS
Definition: LxTBMLTask.cxx:90
LxTBMLFinder::fMuchPoints
std::vector< std::vector< PointDataHolder > > fMuchPoints
Definition: LxTBMLTask.h:103
triggerTimes_trd1_sign0_dist0
static list< pair< timetype, timetype > > triggerTimes_trd1_sign0_dist0
Definition: LxTBTask.cxx:401
LxTbMLStation::HandleMPoint::deltaZr
scaltype deltaZr
Definition: LxTBMLTask.cxx:143
sqrt
friend F32vec4 sqrt(const F32vec4 &a)
Definition: L1/vectors/P4_F32vec4.h:41
CbmPixelHit::GetDx
Double_t GetDx() const
Definition: CbmPixelHit.h:85
CbmMCDataArray::Size
Int_t Size(Int_t fileNumber, Int_t eventNumber)
Definition: CbmMCDataArray.cxx:133
LxTbLayer::yBinLength
scaltype yBinLength
Definition: LxTBBinned2.h:257
LxTbMLStation::Q::Q22
scaltype Q22
Definition: LxTBMLTask.cxx:81
LxTbMLStation::fThetaY
scaltype fThetaY
Definition: LxTBMLTask.cxx:92
RecoTrackData::eventId
Int_t eventId
Definition: LxTBMLTask.cxx:1663
LxTbMLStation::fDeltaThetaY
scaltype fDeltaThetaY
Definition: LxTBMLTask.cxx:91
LxTbAbsorber::zCoord
scaltype zCoord
Definition: LxTBMatEffs.h:23
LxTbTYXBin::yxBins
LxTbYXBin * yxBins
Definition: LxTBBinned.h:181
LxTbDetector::HandleRPoint::HandleLPoint::operator()
void operator()(LxTbBinnedPoint &point)
Definition: LxTBMLTask.cxx:628
triggerTimes_trd0_sign0_dist1
static list< pair< timetype, timetype > > triggerTimes_trd0_sign0_dist1
Definition: LxTBTask.cxx:398
CbmMCTrack::GetPdgCode
Int_t GetPdgCode() const
Definition: CbmMCTrack.h:70
i
int i
Definition: L1/vectors/P4_F32vec4.h:25
tsStartTime
static unsigned long long tsStartTime
Definition: LxTBMLTask.cxx:1037
LxTBMLFinder::Finish
void Finish()
Definition: LxTBMLTask.cxx:1688
shortSignalMCTimes
static list< timetype > shortSignalMCTimes
Definition: LxTBTask.cxx:394
LxTbMLStation::HandleRPoint::operator()
void operator()(LxTbBinnedPoint &point)
Definition: LxTBMLTask.cxx:271
max_ts_time
static Double_t max_ts_time
Definition: LxTBMLTask.cxx:1046
RecoTrackData
Definition: LxTBMLTask.cxx:1662
LxTbYXBin::use
bool use
Definition: LxTBBinned.h:164
NOF_SIGMAS
#define NOF_SIGMAS
Definition: LxTBBinned.h:52
CbmMCDataManager::InitBranch
CbmMCDataArray * InitBranch(const char *name)
Definition: CbmMCDataManager.cxx:106
LxTbLayer::lastXBinNumber
int lastXBinNumber
Definition: LxTBBinned2.h:249
gMCTracks
static vector< vector< MCTrackDesc > > gMCTracks
Definition: CbmBinnedTrackerQANew.cxx:36
CbmMuchPoint::GetDetectorId
Int_t GetDetectorId() const
Definition: CbmMuchPoint.h:69
nofTriggerings
static int nofTriggerings
Definition: LxTBMLTask.cxx:1041
CUR_TIMEBIN_LENGTH
#define CUR_TIMEBIN_LENGTH
Definition: LxTBTask.h:36
LxTBMLFinder::fMuchMCPoints
CbmMCDataArray * fMuchMCPoints
Definition: LxTBMLTask.h:96
LxTbDetector::HandleLastPoint::KFParams::chi2
scaltype chi2
Definition: LxTBMLTask.cxx:655
LxTbDetector::HandleLastPoint::KFParamsCoord
Definition: LxTBMLTask.cxx:639
TIMEBIN_LENGTH
#define TIMEBIN_LENGTH
Definition: LxTBBinned2.h:27
LxTbMLStation::Reconstruct
void Reconstruct()
Definition: LxTBMLTask.cxx:323
LxTbDetector::HandleLastPoint::KFAddTriplet
KFParams KFAddTriplet(KFParams param, LxTbBinnedPoint *trackCandidatePoints[NOF_STATIONS][NOF_LAYERS], int level)
Definition: LxTBMLTask.cxx:740
CbmMCDataArray
Access to a MC data branch for time-based analysis.
Definition: CbmMCDataArray.h:35
LxTbXBin::points
std::list< LxTbBinnedPoint > points
Definition: LxTBBinned.h:145
LxTBMLTask.h
CbmHit::GetRefId
Int_t GetRefId() const
Definition: CbmHit.h:72
LxTbMLStation::Q::Q11
scaltype Q11
Definition: LxTBMLTask.cxx:81
LxTbMLStation::HandleRPoint::HandleRPoint
HandleRPoint(LxTbMLStation &parent)
Definition: LxTBMLTask.cxx:262
triggerTimes_trd1_sign0_dist1
static list< pair< timetype, timetype > > triggerTimes_trd1_sign0_dist1
Definition: LxTBTask.cxx:402
LxTbMLStation::HandleRPoint::station
LxTbMLStation & station
Definition: LxTBMLTask.cxx:256
LxTbMLStation::HandleRPoint::c
timetype c
Definition: LxTBMLTask.cxx:257
LxTBMLFinder::fReconstructor
void * fReconstructor
Definition: LxTBMLTask.h:90
LxTbYXBin::xBins
LxTbXBin * xBins
Definition: LxTBBinned.h:162
LxTbMLStation::HandleLPoint::mPoint
LxTbBinnedPoint * mPoint
Definition: LxTBMLTask.cxx:303
LxTbDetector::HandleLastPoint::KFParams::xParams
KFParamsCoord xParams
Definition: LxTBMLTask.cxx:653
LxTbMLStation::HandleMPoint::station
LxTbMLStation & station
Definition: LxTBMLTask.cxx:141
LxTbDetector::HandleLastPoint::KFParamsCoord::Clear
void Clear()
Definition: LxTBMLTask.cxx:642
CbmMatch.h
LxTbBinnedPoint::isTrd
bool isTrd
Definition: LxTBBinned.h:91
LxTbMLStation::fScatXRL
scaltype fScatXRL
Definition: LxTBMLTask.cxx:89
LxTbMLStation::Init
void Init()
Definition: LxTBMLTask.cxx:116
LxTbLayer::z
scaltype z
Definition: LxTBBinned2.h:259
CbmPixelHit::GetDy
Double_t GetDy() const
Definition: CbmPixelHit.h:86
LxTbMLStation::qs
Q qs[2]
Definition: LxTBMLTask.cxx:95
LxTbDetector::HandleRPoint::operator()
void operator()(LxTbBinnedPoint &point)
Definition: LxTBMLTask.cxx:593
LxTbDetector::HandleRPoint::Init
void Init()
Definition: LxTBMLTask.cxx:591
LxTbLayer::tyxBins
LxTbTYXBin * tyxBins
Definition: LxTBBinned2.h:243
LxTBMLFinder::PointDataHolder::layerNumber
Int_t layerNumber
Definition: LxTBMLTask.h:60
LxTbBinnedPoint::dx
scaltype dx
Definition: LxTBBinned.h:80
triggerTimes_trd0_sign1_dist1
static list< pair< timetype, timetype > > triggerTimes_trd0_sign1_dist1
Definition: LxTBTask.cxx:400
LxTbMLStation::HandleLPoint::Init
void Init()
Definition: LxTBMLTask.cxx:312
LxTbMLStation::fScatYRL
scaltype fScatYRL
Definition: LxTBMLTask.cxx:93
CbmCluster::GetNofDigis
Int_t GetNofDigis() const
Number of digis in cluster.
Definition: CbmCluster.h:69
LxTbDetector::fSignalParticle
SignalParticle * fSignalParticle
Definition: LxTBMLTask.cxx:350
LxTbXBin
Definition: LxTBBinned.h:144
LxTBMLFinder::fMuchPixelDigiMatches
TClonesArray * fMuchPixelDigiMatches
Definition: LxTBMLTask.h:99
LxTBMLFinder::PointDataHolder
Definition: LxTBMLTask.h:51
CbmMuchPoint::GetYOut
Double_t GetYOut() const
Definition: CbmMuchPoint.h:74
LxTbAbsorber::radLength
scaltype radLength
Definition: LxTBMatEffs.h:25
LxTBMLFinder::recoTracks
std::list< Chain * > recoTracks
Definition: LxTBMLTask.h:100
LxTbMLStation::fHandleRPoint
HandleRPoint fHandleRPoint
Definition: LxTBMLTask.cxx:298
LxTBMLFinder::fMuchClusters
TClonesArray * fMuchClusters
Definition: LxTBMLTask.h:98
LxTbXBin::use
bool use
Definition: LxTBBinned.h:146
LxTBMLFinder::Chain::chi2
scaltype chi2
Definition: LxTBMLTask.h:25
LxTbDetector::Init
void Init()
Definition: LxTBMLTask.cxx:361
LxTbMLStation::fDeltaThetaX
scaltype fDeltaThetaX
Definition: LxTBMLTask.cxx:87
LxTbBinnedPoint::triplets
std::list< LxTbBinnedTriplet * > triplets
Definition: LxTBBinned2.h:41
LxTBMLFinder::PointDataHolder::eventId
Int_t eventId
Definition: LxTBMLTask.h:56
LxTbDetector::HandleLastPoint
Definition: LxTBMLTask.cxx:638
LxTBMLFinder::fIsEvByEv
bool fIsEvByEv
Definition: LxTBMLTask.h:91
RecoTrackData::RecoTrackData
RecoTrackData(Int_t eId, Int_t tId)
Definition: LxTBMLTask.cxx:1666
LxTbLayer::nofYXBins
int nofYXBins
Definition: LxTBBinned2.h:245
SaveHisto
static void SaveHisto(TH1 *histo)
Definition: LxTBMLTask.cxx:1677
LxTBMLFinder::TrackDataHolder
Definition: LxTBMLTask.h:63
LxTbMLStation::fHandleLPoint
HandleLPoint fHandleLPoint
Definition: LxTBMLTask.cxx:321
CbmMuchGeoScheme::GetLayerIndex
static Int_t GetLayerIndex(Int_t address)
Definition: CbmMuchGeoScheme.h:71
CbmMuchGeoScheme::GetStationIndex
static Int_t GetStationIndex(Int_t address)
Definition: CbmMuchGeoScheme.h:68
LxTbDetector::HandleLastPoint::KFParamsCoord::C12
scaltype C12
Definition: LxTBMLTask.cxx:640
LxTbMLStation::HandleRPoint::scatXLL
scaltype scatXLL
Definition: LxTBMLTask.cxx:259
ClassImp
ClassImp(LxTBMLFinder) Double_t speedOfLight=0
LxTbMLStation::HandleMPoint::Init
void Init()
Definition: LxTBMLTask.cxx:149
LxTbDetector::HandleLastPoint::HandleLastPoint
HandleLastPoint(LxTbDetector &parent)
Definition: LxTBMLTask.cxx:666
LxTBMLFinder
Definition: LxTBMLTask.h:21
CbmMuchPoint.h
LxTbMLStation::HandleRPoint
Definition: LxTBMLTask.cxx:255
LxTBMLFinder::PointDataHolder::y
Double_t y
Definition: LxTBMLTask.h:53
LxTbMLStation::HandleMPoint::HandleMPoint
HandleMPoint(LxTbMLStation &parent)
Definition: LxTBMLTask.cxx:146
LxTbLayer::minY
scaltype minY
Definition: LxTBBinned2.h:252
LxTbDetector::HandleRPoint
Definition: LxTBMLTask.cxx:583
RTDLess
Definition: LxTBMLTask.cxx:1669
LxTbBinnedTriplet::lPoint
LxTbBinnedPoint * lPoint
Definition: LxTBBinned2.h:108
mcTracks
static vector< vector< QAMCTrack > > mcTracks
Definition: CbmTofHitFinderTBQA.cxx:112
LxTbDetector::HandleLastPoint::HandleTriplet
void HandleTriplet(LxTbBinnedTriplet *triplet, LxTbBinnedPoint *trackCandidatePoints[NOF_STATIONS][NOF_LAYERS], list< LxTBMLFinder::Chain > &chains, int level, KFParams kfParams)
Definition: LxTBMLTask.cxx:756
LxTbDetector::HandleLastPoint::KFParamsCoord::coord
scaltype coord
Definition: LxTBMLTask.cxx:640
LxTbMLStation::fStationNumber
int fStationNumber
Definition: LxTBMLTask.cxx:84
CbmMCDataArray::Get
TObject * Get(const CbmLink *lnk)
Definition: CbmMCDataArray.h:47
CbmHit::GetTime
Double_t GetTime() const
Definition: CbmHit.h:75
ts_points
static list< LxTbBinnedPoint > ts_points
Definition: LxTBMLTask.cxx:1047
LxTbDetector::HandleLastPoint::KFParamsCoord::C21
scaltype C21
Definition: LxTBMLTask.cxx:640
LxTBMLFinder::Chain
Definition: LxTBMLTask.h:23
LxTBMLFinder::LxTBMLFinder
LxTBMLFinder()
Definition: LxTBMLTask.cxx:867
LxTbDetector::HandleLastPoint::KFParamsCoord::tg
scaltype tg
Definition: LxTBMLTask.cxx:640
LxTbLayer::lastTimeBinNumber
int lastTimeBinNumber
Definition: LxTBBinned2.h:247
CbmHit::GetAddress
Int_t GetAddress() const
Definition: CbmHit.h:73
LxTbMLStation::Init2
void Init2()
Definition: LxTBMLTask.cxx:124
LxTbMLStation::HandleMPoint::scatXRL
scaltype scatXRL
Definition: LxTBMLTask.cxx:144
LxTbDetector::HandleLastPoint::KFParamsCoord::C22
scaltype C22
Definition: LxTBMLTask.cxx:640
LxTBMLFinder::Init
InitStatus Init()
Definition: LxTBMLTask.cxx:885
LxTbBinnedPoint::PointDesc::trackId
Int_t trackId
Definition: LxTBBinned.h:97
LxSettings.h
LxTbDetector::HandleLastPoint::KFAddPoint
KFParams KFAddPoint(KFParams prevParam, scaltype m[2], scaltype V[2], int stationNumber, int layerNumber)
Definition: LxTBMLTask.cxx:718
LxTbMLStation::HandleMPoint::operator()
void operator()(LxTbBinnedPoint &point)
Definition: LxTBMLTask.cxx:157
fullDuration
static long fullDuration
Definition: LxTBMLTask.cxx:1040
FindGeoChild
static void FindGeoChild(TGeoNode *node, const char *name, list< TGeoNode * > &results)
Definition: LxTBMLTask.cxx:327
IterateNeighbourhood
void IterateNeighbourhood(LxTbLayer &layer, scaltype x, scaltype dx, scaltype scatX, scaltype y, scaltype dy, scaltype scatY, timetype t, timetype dt, HandlePoint &handlePoint)
Definition: LxTBBinned2.h:382
LxTbLayer::lastYBinNumber
int lastYBinNumber
Definition: LxTBBinned2.h:248
LxTBMLFinder::fNEvents
int fNEvents
Definition: LxTBMLTask.h:107
LxTbBinnedTriplet::neighbours
std::list< LxTbBinnedPoint * > neighbours
Definition: LxTBBinned2.h:114
magneticFieldCorrections
static scaltype magneticFieldCorrections[]
Definition: LxTBMLTask.cxx:57
LxTbMLStation::fThetaX
scaltype fThetaX
Definition: LxTBMLTask.cxx:88
LxTbBinnedPoint::y
scaltype y
Definition: LxTBBinned.h:81
LxTbLayer::nofTYXBins
int nofTYXBins
Definition: LxTBBinned2.h:244
LxTBMLFinder::fMuchPixelHits
TClonesArray * fMuchPixelHits
Definition: LxTBMLTask.h:97
LxTbBinnedPoint::PointDesc
Definition: LxTBBinned.h:94
gElectronMass
static scaltype gElectronMass
Definition: LxTBMatEffs.h:20
LxTbDetector::LxTbDetector
LxTbDetector(int nofxb, int nofyb, int noftb)
Definition: LxTBMLTask.cxx:352
SignalParticle
Definition: LxTBMLTask.cxx:339
LxTbDetector::HandleRPoint::deltaZ
scaltype deltaZ
Definition: LxTBMLTask.cxx:586
LxTbMLStation::fLayers
LxTbLayer * fLayers
Definition: LxTBMLTask.cxx:85
CbmMCDataManager
Task class creating and managing CbmMCDataArray objects.
Definition: CbmMCDataManager.h:27
particleDescs
static SignalParticle particleDescs[]
Definition: LxTBMLTask.cxx:345
LxTbDetector::SetMinT
void SetMinT(timetype v)
Definition: LxTBMLTask.cxx:449
LxTbBinnedPoint::PointDesc::eventId
Int_t eventId
Definition: LxTBBinned.h:95
LxTbBinnedPoint::dy
scaltype dy
Definition: LxTBBinned.h:82
LxTbDetector::HandleLastPoint::KFParams
Definition: LxTBMLTask.cxx:652
LxTbBinnedTriplet::ty
scaltype ty
Definition: LxTBBinned2.h:112
min_ts_time
static Double_t min_ts_time
Definition: LxTBMLTask.cxx:1045
LxTbDetector::HandleLastPoint::HandlePoint
void HandlePoint(LxTbBinnedPoint *point, LxTbBinnedPoint *trackCandidatePoints[NOF_STATIONS][NOF_LAYERS], list< LxTBMLFinder::Chain > &chains, int level, KFParams kfParams)
Definition: LxTBMLTask.cxx:778
LxTbMLStation::fScatYLS
scaltype fScatYLS
Definition: LxTBMLTask.cxx:94
LxTbMLStation::HandleLPoint::HandleLPoint
HandleLPoint(LxTbMLStation &parent)
Definition: LxTBMLTask.cxx:306
LxTbMLStation::Q::Q12
scaltype Q12
Definition: LxTBMLTask.cxx:81
LxTbTYXBin::use
bool use
Definition: LxTBBinned.h:183
LxTbDetector::HandleLastPoint::KFParamsCoord::C11
scaltype C11
Definition: LxTBMLTask.cxx:640
LAST_STATION
#define LAST_STATION
Definition: LxTBBinned2.h:25
CalcThetaPrj
static scaltype CalcThetaPrj(scaltype E, scaltype x, const LxTbAbsorber *mat)
Definition: LxTBMatEffs.h:125
LxTbMLStation::fAbsorber
LxTbAbsorber fAbsorber
Definition: LxTBMLTask.cxx:86
triggerTimes_trd1_sign1_dist0
static list< pair< timetype, timetype > > triggerTimes_trd1_sign1_dist0
Definition: LxTBTask.cxx:403
LxTbMLStation::HandleRPoint::deltaZl
scaltype deltaZl
Definition: LxTBMLTask.cxx:258
points
TClonesArray * points
Definition: Analyze_matching.h:18
CbmMCTrack.h
LxTbMLStation::Q
Definition: LxTBMLTask.cxx:80
LxTBMLFinder::PointDataHolder::t
Double_t t
Definition: LxTBMLTask.h:55
v
__m128 v
Definition: L1/vectors/P4_F32vec4.h:1
LxTbBinnedPoint::dt
timetype dt
Definition: LxTBBinned.h:84
NOF_STATIONS
#define NOF_STATIONS
Definition: LxTBBinned2.h:24
timetype
#define timetype
Definition: CbmGlobalTrackingDefs.h:18
CbmCluster
Base class for cluster objects.
Definition: CbmCluster.h:26
LxCA.h
LxTbYXBin
Definition: LxTBBinned.h:161
LxTbDetector::HandleGeometry
void HandleGeometry()
Definition: LxTBMLTask.cxx:454
LxTbAbsorber::width
scaltype width
Definition: LxTBMatEffs.h:24
LxTBMLFinder::PointDataHolder::trackId
Int_t trackId
Definition: LxTBMLTask.h:57
LxTbDetector::HandleRPoint::c
timetype c
Definition: LxTBMLTask.cxx:587
CbmMCTrack
Definition: CbmMCTrack.h:34
speedOfLight
Double_t speedOfLight
CbmMuchPoint::GetYIn
Double_t GetYIn() const
Definition: CbmMuchPoint.h:71
x
Double_t x
Definition: CbmMvdSensorDigiToHitTask.cxx:68
LxTbMLStation::HandleMPoint
Definition: LxTBMLTask.cxx:140
LxTbMLStation::fHandleMPoint
HandleMPoint fHandleMPoint
Definition: LxTBMLTask.cxx:253
LxTBMLFinder::fNofXBins
int fNofXBins
Definition: LxTBMLTask.h:92
LxTbAbsorber::A
scaltype A
Definition: LxTBMatEffs.h:28
m
__m128 m
Definition: L1/vectors/P4_F32vec4.h:26
CbmMuchPixelHit.h
Class for pixel hits in MUCH detector.
y
Double_t y
Definition: CbmMvdSensorDigiToHitTask.cxx:68
LxTbAbsorber
Definition: LxTBMatEffs.h:22
LxTbDetector::HandleLastPoint::detector
LxTbDetector & detector
Definition: LxTBMLTask.cxx:664
LxTbMLStation::HandleLPoint::rPoint
LxTbBinnedPoint * rPoint
Definition: LxTBMLTask.cxx:304
LxTbDetector::HandleRPoint::rStation
LxTbMLStation * rStation
Definition: LxTBMLTask.cxx:584
CbmMuchPoint::GetXIn
Double_t GetXIn() const
Definition: CbmMuchPoint.h:70
LxTbLayer::maxY
scaltype maxY
Definition: LxTBBinned2.h:253
longSignalMCTimes
static list< timetype > longSignalMCTimes
Definition: LxTBTask.cxx:395
LxTbDetector::fStations
LxTbMLStation * fStations
Definition: LxTBMLTask.cxx:348
LxTBMLFinder::TrackDataHolder::isSignal
bool isSignal
Definition: LxTBMLTask.h:65
RTDLess::operator()
bool operator()(const RecoTrackData &x, const RecoTrackData &y) const
Definition: LxTBMLTask.cxx:1670
LxTbMLStation::SetMinT
void SetMinT(timetype v)
Definition: LxTBMLTask.cxx:135
LxTbDetector::HandleRPoint::HandleLPoint::HandleLPoint
HandleLPoint()
Definition: LxTBMLTask.cxx:626
IterateLayer
void IterateLayer(LxTbLayer &layer, HandlePoint &handlePoint)
Definition: LxTBBinned2.h:351
LxTbMLStation::LxTbMLStation
LxTbMLStation(int stationNumber, int nofxb, int nofyb, int noftb)
Definition: LxTBMLTask.cxx:97
EnergyLoss
static scaltype EnergyLoss(scaltype E, scaltype L, const LxTbAbsorber *mat)
Definition: LxTBMatEffs.h:79
LxTbBinnedTriplet
Definition: LxTBBinned2.h:107
CbmMuchPixelHit
Definition: CbmMuchPixelHit.h:17
CbmCluster.h
Base class for cluster objects.
LxTbDetector::HandleLastPoint::operator()
void operator()(LxTbBinnedPoint &point)
Definition: LxTBMLTask.cxx:793
LxTbDetector::fHandleRPoint
HandleRPoint fHandleRPoint
Definition: LxTBMLTask.cxx:636
LxTbDetector::HandleRPoint::HandleLPoint::rTriplet
LxTbBinnedTriplet * rTriplet
Definition: LxTBMLTask.cxx:625
CbmMuchGeoScheme.h
LxTbDetector
Definition: LxTBMLTask.cxx:347
RecoTrackData::trackId
Int_t trackId
Definition: LxTBMLTask.cxx:1664
LxTbBinnedTriplet::dty
scaltype dty
Definition: LxTBBinned2.h:113
LxTbBinnedPoint::t
timetype t
Definition: LxTBBinned.h:83
LxTbMLStation::HandleMPoint::c
timetype c
Definition: LxTBMLTask.cxx:142
LxTbDetector::HandleRPoint::HandleRPoint
HandleRPoint()
Definition: LxTBMLTask.cxx:589
LxTBMLFinder::TrackDataHolder::pointInds
Int_t pointInds[NOF_STATIONS][NOF_LAYERS]
Definition: LxTBMLTask.h:64
LxTbBinnedPoint
Definition: LxTBBinned.h:78
LxTbDetector::HandleLastPoint::KFParams::Clear
void Clear()
Definition: LxTBMLTask.cxx:657
LxTbBinnedTriplet::tx
scaltype tx
Definition: LxTBBinned2.h:110
LxTbMLStation::HandleLPoint::deltaZ
scaltype deltaZ
Definition: LxTBMLTask.cxx:302
LxTBMatEffs.h
gStations
static LxTbMLStation * gStations
Definition: LxTBMLTask.cxx:865
LxTbBinnedTriplet::rPoint
LxTbBinnedPoint * rPoint
Definition: LxTBBinned2.h:109
LxTBMLFinder::fEventTimes
std::vector< Double_t > fEventTimes
Definition: LxTBMLTask.h:105
LxTbLayer::minX
scaltype minX
Definition: LxTBBinned2.h:250
CbmCluster::GetDigi
Int_t GetDigi(Int_t index) const
Get digi at position index.
Definition: CbmCluster.h:76
LxTbBinnedPoint::layerNumber
Int_t layerNumber
Definition: LxTBBinned2.h:48
LxTBMLFinder::Exec
void Exec(Option_t *opt)
Definition: LxTBMLTask.cxx:1133
LxTBMLFinder::Chain::points
LxTbBinnedPoint * points[NOF_STATIONS][NOF_LAYERS]
Definition: LxTBMLTask.h:24
LxTbMLStation::HandleLPoint::operator()
void operator()(LxTbBinnedPoint &point)
Definition: LxTBMLTask.cxx:314
LxTbBinnedPoint::x
scaltype x
Definition: LxTBBinned.h:79
LxTbDetector::Reconstruct
void Reconstruct()
Definition: LxTBMLTask.cxx:834
LxTBMLFinder::fNofYBins
int fNofYBins
Definition: LxTBMLTask.h:93
LxTBMLFinder::fNofTBins
int fNofTBins
Definition: LxTBMLTask.h:94
LxTbLayer::maxX
scaltype maxX
Definition: LxTBBinned2.h:251
LxTbBinnedPoint::stationNumber
Int_t stationNumber
Definition: LxTBBinned.h:92
LxTbLayer
Definition: LxTBBinned2.h:201
SpliceTriggerings
static void SpliceTriggerings(list< pair< timetype, timetype >> &out, LxTbBinnedFinder::TriggerTimeArray &in)
Definition: LxTBTask.cxx:1298
triggerTimes_trd1_sign1_dist1
static list< pair< timetype, timetype > > triggerTimes_trd1_sign1_dist1
Definition: LxTBTask.cxx:404
LxTbBinnedTriplet::dtx
scaltype dtx
Definition: LxTBBinned2.h:111
LxTbDetector::Clear
void Clear()
Definition: LxTBMLTask.cxx:442
SignalParticle::fMinEnergy
scaltype fMinEnergy
Definition: LxTBMLTask.cxx:342
LxTbLayer::nofXBins
int nofXBins
Definition: LxTBBinned2.h:246
LxTBMLFinder::PointDataHolder::x
Double_t x
Definition: LxTBMLTask.h:52
LxTbAbsorber::rho
scaltype rho
Definition: LxTBMatEffs.h:26
LxTbLayer::xBinLength
scaltype xBinLength
Definition: LxTBBinned2.h:256