CbmRoot
CbmGlobalTrackingTof.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 "CbmGlobalTrackingTof.h"
8 #include "CbmPixelHit.h"
9 #include "CbmTofAddress.h"
10 #include "CbmTofHit.h"
11 #include "TGeoBBox.h"
12 #include "TGeoManager.h"
13 #include "TMath.h"
14 #include "base/CbmLitToolFactory.h"
15 #include "utils/CbmLitConverter.h"
16 #include <iostream>
17 #include <utility>
18 
19 TrackPropagatorPtr fPropagator;
20 TrackPropagatorPtr fLinePropagator;
21 TrackUpdatePtr fFilter;
22 
23 using std::cout;
24 using std::endl;
25 using std::list;
26 using std::map;
27 using std::pair;
28 using std::set;
29 using std::swap;
30 
31 /* (VF) Unused function
32 static bool GaussSolve(double coeffs[3][4], double result[3])
33 {
34  int indices[3];
35 
36  for (int i = 0; i < 3; ++i)
37  indices[i] = i;
38 
39  for (int i = 2; i > 0; --i)
40  {
41  if (0 == coeffs[i][i + 1])
42  {
43  for (int j = i; j > 0; --j)
44  {
45  if (0 != coeffs[i][j])
46  {
47  for (int k = 0; k < 3; ++k)
48  swap(coeffs[k][j], coeffs[k][i + 1]);
49 
50  swap(indices[j - 1], indices[i]);
51  break;
52  }
53  }
54  }
55 
56  double den = coeffs[i][i + 1];
57 
58  if (0 == den)
59  return false;
60 
61  for (int j = 0; j < i; ++j)
62  {
63  for (int k = 0; k <= i; ++k)
64  coeffs[j][k] += -coeffs[j][i + 1] * coeffs[i][k] / den;
65  }
66  }
67 
68  double solution[3];
69 
70  for (int i = 0; i < 3; ++i)
71  {
72  double s = coeffs[i][0];
73 
74  for (int j = 0; j < i; ++j)
75  s += coeffs[i][j + 1] * solution[j];
76 
77  if (0 == coeffs[i][i + 1])
78  return false;
79 
80  solution[i] = -s / coeffs[i][i + 1];
81  }
82 
83  for (int i = 0; i < 3; ++i)
84  result[indices[i]] = solution[i];
85 
86  return true;
87 }
88 */
89 
90 struct Point {
91  double x;
92  double y;
93  double z;
94 };
95 
96 struct Line {
97  double x0;
98  double y0;
99  double z0;
100 
101  double cosX;
102  double cosY;
103  double cosZ;
104 };
105 
106 struct Segment {
107  Point a;
108  Point b;
109 };
110 
111 struct Plane {
112  double x0;
113  double y0;
114  double z0;
115 
116  double cosX1;
117  double cosY1;
118  double cosZ1;
119 
120  double cosX2;
121  double cosY2;
122  double cosZ2;
123 };
124 
125 /* (VF) Unused function
126 static bool Intersect(const Plane& plane, const Line& line, double result[3])
127 {
128  double coeffs[3][4] = {
129  { plane.x0 - line.x0, plane.cosX1, plane.cosX2, -line.cosX },
130  { plane.y0 - line.y0, plane.cosY1, plane.cosY2, -line.cosY },
131  { plane.z0 - line.z0, plane.cosZ1, plane.cosZ2, -line.cosZ }
132  };
133 
134  return GaussSolve(coeffs, result);
135 }
136 */
137 
138 struct Rectangle {
139  Plane plane;
140  double s1;
141  double s2;
142 };
143 
144 /* (VF) Unused function
145 static bool Intersect(const Rectangle& rect, const Line& line, double result[3])
146 {
147  if (!Intersect(rect.plane, line, result))
148  return false;
149 
150  return 0 <= result[0] && result[0] <= rect.s1 && 0 <= result[1] && result[1] <= rect.s2;
151 }
152 */
153 
154 /* (VF) Unused function
155 static bool Intersect(const Rectangle& rect, const Line& line)
156 {
157  double result[3];
158 
159  if (!Intersect(rect.plane, line, result))
160  return false;
161 
162  return 0 <= result[0] && result[0] <= rect.s1 && 0 <= result[1] && result[1] <= rect.s2;
163 }
164 */
165 
166 /* (VF) Unused function
167 static bool Intersect(const Rectangle& rect, const Segment& segment)
168 {
169  double len = TMath::Sqrt((segment.b.x - segment.a.x) * (segment.b.x - segment.a.x) +
170  (segment.b.y - segment.a.y) * (segment.b.y - segment.a.y) +
171  (segment.b.z - segment.a.z) * (segment.b.z - segment.a.z));
172  double cosX = (segment.b.x - segment.a.x) / len;
173  double cosY = (segment.b.y - segment.a.y) / len;
174  double cosZ = (segment.b.z - segment.a.z) / len;
175  Line line = { segment.a.x, segment.a.y, segment.a.z, cosX, cosY, cosZ };
176  double result[3];
177 
178  if (!Intersect(rect.plane, line, result))
179  return false;
180 
181  return 0 <= result[2] && result[2] <= len;
182 }
183 */
184 
185 struct Direction {
186  double cosX;
187  double cosY;
188  double cosZ;
189 };
190 
191 struct TBin {
192  list<Int_t> fHitInds;
193 
194  TBin() : fHitInds() {}
195 
196  void Clear() { fHitInds.clear(); }
197 };
198 
199 #ifdef CBM_GLOBALTB_TOF_3D_CUBOIDS
200 struct Cuboid {
201  double width;
202  double height;
203  double thickness;
204  double minX;
205  double maxX;
206  double minY;
207  double maxY;
208  double minZ;
209  double maxZ;
210  Direction dirWidth;
211  Direction dirHeight;
212  Direction dirThickness;
213  Point vertices[2][2][2];
214  Segment edges[12];
215  Rectangle faces[3][2];
216 
217  //Plane plane1;
218  //Plane plane2;
219 
220  int fNofTBins;
221  TBin* fTBins;
222 
223  Int_t fFullModId;
224 
225  void Clear() {
226  for (int i = 0; i < fNofTBins; ++i)
227  fTBins[i].Clear();
228  }
229 
230  void Calc(int nofTBins = 0) {
231  fNofTBins = nofTBins;
232 
233  if (fNofTBins) fTBins = new TBin[fNofTBins];
234 
235  dirWidth = {(vertices[0][0][1].x - vertices[0][0][0].x) / width,
236  (vertices[0][0][1].y - vertices[0][0][0].y) / width,
237  (vertices[0][0][1].z - vertices[0][0][0].z) / width};
238  dirHeight = {(vertices[0][1][0].x - vertices[0][0][0].x) / height,
239  (vertices[0][1][0].y - vertices[0][0][0].y) / height,
240  (vertices[0][1][0].z - vertices[0][0][0].z) / height};
241  dirThickness = {(vertices[1][0][0].x - vertices[0][0][0].x) / thickness,
242  (vertices[1][0][0].y - vertices[0][0][0].y) / thickness,
243  (vertices[1][0][0].z - vertices[0][0][0].z) / thickness};
244 
245  minX = 1000000;
246  maxX = -1000000;
247  minY = 1000000;
248  maxY = -1000000;
249  minZ = 1000000;
250  maxZ = -1000000;
251 
252  int dims[3];
253  int ind = 0;
254 
255  for (int i = 0; i < 2; ++i) {
256  dims[0] = i;
257 
258  for (int j = 0; j < 2; ++j) {
259  dims[1] = j;
260 
261  for (int k = 0; k < 2; ++k) {
262  dims[2] = k;
263 
264  for (int l = 0; l < 3; ++l) {
265  if (0 == dims[l]) {
266  int dims2[3] = {dims[0], dims[1], dims[2]};
267  dims2[l] = 1;
268  edges[ind++] = {{vertices[dims[0]][dims[1]][dims[2]].x,
269  vertices[dims[0]][dims[1]][dims[2]].y,
270  vertices[dims[0]][dims[1]][dims[2]].z},
271  {vertices[dims2[0]][dims2[1]][dims2[2]].x,
272  vertices[dims2[0]][dims2[1]][dims2[2]].y,
273  vertices[dims2[0]][dims2[1]][dims2[2]].z}};
274  }
275  }
276 
277  if (minX > vertices[i][j][k].x) minX = vertices[i][j][k].x;
278 
279  if (maxX < vertices[i][j][k].x) maxX = vertices[i][j][k].x;
280 
281  if (minY > vertices[i][j][k].y) minY = vertices[i][j][k].y;
282 
283  if (maxY < vertices[i][j][k].y) maxY = vertices[i][j][k].y;
284 
285  if (minZ > vertices[i][j][k].z) minZ = vertices[i][j][k].z;
286 
287  if (maxZ < vertices[i][j][k].z) maxZ = vertices[i][j][k].z;
288  }
289  }
290  }
291 
292  for (int i = 0; i < 3; ++i) {
293  int j = (i + 1) % 3;
294  int k = (i + 2) % 3;
295  int dims0[3] = {0, 0, 0};
296  int dims1[3] = {0, 0, 0};
297  int dims2[3] = {0, 0, 0};
298  dims1[j] = 1;
299  dims2[k] = 1;
300  double len1 =
301  TMath::Sqrt((vertices[dims1[0]][dims1[1]][dims1[2]].x
302  - vertices[dims0[0]][dims0[1]][dims0[2]].x)
303  * (vertices[dims1[0]][dims1[1]][dims1[2]].x
304  - vertices[dims0[0]][dims0[1]][dims0[2]].x)
305  + (vertices[dims1[0]][dims1[1]][dims1[2]].y
306  - vertices[dims0[0]][dims0[1]][dims0[2]].y)
307  * (vertices[dims1[0]][dims1[1]][dims1[2]].y
308  - vertices[dims0[0]][dims0[1]][dims0[2]].y)
309  + (vertices[dims1[0]][dims1[1]][dims1[2]].z
310  - vertices[dims0[0]][dims0[1]][dims0[2]].z)
311  * (vertices[dims1[0]][dims1[1]][dims1[2]].z
312  - vertices[dims0[0]][dims0[1]][dims0[2]].z));
313  double cosX1 = (vertices[dims1[0]][dims1[1]][dims1[2]].x
314  - vertices[dims0[0]][dims0[1]][dims0[2]].x)
315  / len1;
316  double cosY1 = (vertices[dims1[0]][dims1[1]][dims1[2]].y
317  - vertices[dims0[0]][dims0[1]][dims0[2]].y)
318  / len1;
319  double cosZ1 = (vertices[dims1[0]][dims1[1]][dims1[2]].z
320  - vertices[dims0[0]][dims0[1]][dims0[2]].z)
321  / len1;
322 
323  double len2 =
324  TMath::Sqrt((vertices[dims2[0]][dims2[1]][dims2[2]].x
325  - vertices[dims0[0]][dims0[1]][dims0[2]].x)
326  * (vertices[dims2[0]][dims2[1]][dims2[2]].x
327  - vertices[dims0[0]][dims0[1]][dims0[2]].x)
328  + (vertices[dims2[0]][dims2[1]][dims2[2]].y
329  - vertices[dims0[0]][dims0[1]][dims0[2]].y)
330  * (vertices[dims2[0]][dims2[1]][dims2[2]].y
331  - vertices[dims0[0]][dims0[1]][dims0[2]].y)
332  + (vertices[dims2[0]][dims2[1]][dims2[2]].z
333  - vertices[dims0[0]][dims0[1]][dims0[2]].z)
334  * (vertices[dims2[0]][dims2[1]][dims2[2]].z
335  - vertices[dims0[0]][dims0[1]][dims0[2]].z));
336  double cosX2 = (vertices[dims2[0]][dims2[1]][dims2[2]].x
337  - vertices[dims0[0]][dims0[1]][dims0[2]].x)
338  / len2;
339  double cosY2 = (vertices[dims2[0]][dims2[1]][dims2[2]].y
340  - vertices[dims0[0]][dims0[1]][dims0[2]].y)
341  / len2;
342  double cosZ2 = (vertices[dims2[0]][dims2[1]][dims2[2]].z
343  - vertices[dims0[0]][dims0[1]][dims0[2]].z)
344  / len2;
345 
346  faces[i][0] = {{vertices[dims0[0]][dims0[1]][dims0[2]].x,
347  vertices[dims0[0]][dims0[1]][dims0[2]].y,
348  vertices[dims0[0]][dims0[1]][dims0[2]].z,
349  cosX1,
350  cosY1,
351  cosZ1,
352  cosX2,
353  cosY2,
354  cosZ2},
355  len1,
356  len2};
357 
358  dims0[i] = 1;
359  dims1[i] = 1;
360  dims2[i] = 1;
361  len1 = TMath::Sqrt((vertices[dims1[0]][dims1[1]][dims1[2]].x
362  - vertices[dims0[0]][dims0[1]][dims0[2]].x)
363  * (vertices[dims1[0]][dims1[1]][dims1[2]].x
364  - vertices[dims0[0]][dims0[1]][dims0[2]].x)
365  + (vertices[dims1[0]][dims1[1]][dims1[2]].y
366  - vertices[dims0[0]][dims0[1]][dims0[2]].y)
367  * (vertices[dims1[0]][dims1[1]][dims1[2]].y
368  - vertices[dims0[0]][dims0[1]][dims0[2]].y)
369  + (vertices[dims1[0]][dims1[1]][dims1[2]].z
370  - vertices[dims0[0]][dims0[1]][dims0[2]].z)
371  * (vertices[dims1[0]][dims1[1]][dims1[2]].z
372  - vertices[dims0[0]][dims0[1]][dims0[2]].z));
373  cosX1 = (vertices[dims1[0]][dims1[1]][dims1[2]].x
374  - vertices[dims0[0]][dims0[1]][dims0[2]].x)
375  / len1;
376  cosY1 = (vertices[dims1[0]][dims1[1]][dims1[2]].y
377  - vertices[dims0[0]][dims0[1]][dims0[2]].y)
378  / len1;
379  cosZ1 = (vertices[dims1[0]][dims1[1]][dims1[2]].z
380  - vertices[dims0[0]][dims0[1]][dims0[2]].z)
381  / len1;
382 
383  len2 = TMath::Sqrt((vertices[dims2[0]][dims2[1]][dims2[2]].x
384  - vertices[dims0[0]][dims0[1]][dims0[2]].x)
385  * (vertices[dims2[0]][dims2[1]][dims2[2]].x
386  - vertices[dims0[0]][dims0[1]][dims0[2]].x)
387  + (vertices[dims2[0]][dims2[1]][dims2[2]].y
388  - vertices[dims0[0]][dims0[1]][dims0[2]].y)
389  * (vertices[dims2[0]][dims2[1]][dims2[2]].y
390  - vertices[dims0[0]][dims0[1]][dims0[2]].y)
391  + (vertices[dims2[0]][dims2[1]][dims2[2]].z
392  - vertices[dims0[0]][dims0[1]][dims0[2]].z)
393  * (vertices[dims2[0]][dims2[1]][dims2[2]].z
394  - vertices[dims0[0]][dims0[1]][dims0[2]].z));
395  cosX2 = (vertices[dims2[0]][dims2[1]][dims2[2]].x
396  - vertices[dims0[0]][dims0[1]][dims0[2]].x)
397  / len2;
398  cosY2 = (vertices[dims2[0]][dims2[1]][dims2[2]].y
399  - vertices[dims0[0]][dims0[1]][dims0[2]].y)
400  / len2;
401  cosZ2 = (vertices[dims2[0]][dims2[1]][dims2[2]].z
402  - vertices[dims0[0]][dims0[1]][dims0[2]].z)
403  / len2;
404 
405  faces[i][1] = {{vertices[dims0[0]][dims0[1]][dims0[2]].x,
406  vertices[dims0[0]][dims0[1]][dims0[2]].y,
407  vertices[dims0[0]][dims0[1]][dims0[2]].z,
408  cosX1,
409  cosY1,
410  cosZ1,
411  cosX2,
412  cosY2,
413  cosZ2},
414  len1,
415  len2};
416  }
417  }
418 };
419 
420 static bool Inside(const Cuboid& cuboid, const Point& point) {
421  const Point& O = cuboid.vertices[0][0][0];
422  double projX = (point.x - O.x) * cuboid.dirWidth.cosX
423  + (point.y - O.y) * cuboid.dirWidth.cosY
424  + (point.z - O.z) * cuboid.dirWidth.cosZ;
425 
426  if (0 > projX || projX > cuboid.width) return false;
427 
428  double projY = (point.x - O.x) * cuboid.dirHeight.cosX
429  + (point.y - O.y) * cuboid.dirHeight.cosY
430  + (point.z - O.z) * cuboid.dirHeight.cosZ;
431 
432  if (0 > projY || projY > cuboid.height) return false;
433 
434  double projZ = (point.x - O.x) * cuboid.dirThickness.cosX
435  + (point.y - O.y) * cuboid.dirThickness.cosY
436  + (point.z - O.z) * cuboid.dirThickness.cosZ;
437 
438  if (0 > projZ || projZ > cuboid.thickness) return false;
439 
440  return true;
441 }
442 
443 static bool Intersect(const Cuboid& cuboid1, const Cuboid& cuboid2) {
444  // Check if cuboid2 is inside cuboid1
445  for (int i = 0; i < 2; ++i) {
446  for (int j = 0; j < 2; ++j) {
447  for (int k = 0; k < 2; ++k) {
448  if (Inside(cuboid1, cuboid2.vertices[i][j][k])) return true;
449  }
450  }
451  }
452 
453  // Check if cuboid1 is inside cuboid2
454  for (int i = 0; i < 2; ++i) {
455  for (int j = 0; j < 2; ++j) {
456  for (int k = 0; k < 2; ++k) {
457  if (Inside(cuboid2, cuboid1.vertices[i][j][k])) return true;
458  }
459  }
460  }
461 
462  // Check if one edge of the cuboid2 intersects one face of the cuboid1
463  for (int i = 0; i < 3; ++i) {
464  for (int j = 0; j < 2; ++j) {
465  const Rectangle& face = cuboid1.faces[i][j];
466 
467  for (int k = 0; k < 12; ++k) {
468  const Segment& edge = cuboid2.edges[k];
469 
470  if (Intersect(face, edge)) return true;
471  }
472  }
473  }
474 
475  // Check if one edge of the cuboid1 intersects one face of the cuboid2
476  for (int i = 0; i < 3; ++i) {
477  for (int j = 0; j < 2; ++j) {
478  const Rectangle& face = cuboid2.faces[i][j];
479 
480  for (int k = 0; k < 12; ++k) {
481  const Segment& edge = cuboid1.edges[k];
482 
483  if (Intersect(face, edge)) return true;
484  }
485  }
486  }
487 
488  return false;
489 }
490 #endif //CBM_GLOBALTB_TOF_3D_CUBOIDS
491 
492 struct XBin {
493 #ifdef CBM_GLOBALTB_TOF_3D_CUBOIDS
494  list<Cuboid*> fCuboids;
495 
496  void Clear() {
497  for (list<Cuboid*>::iterator i = fCuboids.begin(); i != fCuboids.end();
498  ++i) {
499  Cuboid* cuboid = *i;
500  cuboid->Clear();
501  }
502  }
503 #else //CBM_GLOBALTB_TOF_3D_CUBOIDS
504  int fNofTBins;
505  TBin* fTBins;
506  explicit XBin(int nofTBins)
507  : fNofTBins(nofTBins), fTBins(new TBin[fNofTBins]) {}
508 
509  void Clear() {
510  for (int i = 0; i < fNofTBins; ++i) {
511  fTBins[i].Clear();
512  }
513  }
514 #endif //CBM_GLOBALTB_TOF_3D_CUBOIDS
515 };
516 
517 struct YBin {
518  int fNofXBins;
519  XBin* fXBins;
520 
521  YBin(int nofXBins, int nofTBins)
522  : fNofXBins(nofXBins)
523  , fXBins(
524  reinterpret_cast<XBin*>(new unsigned char[fNofXBins * sizeof(XBin)])) {
525  for (int i = 0; i < fNofXBins; ++i)
526 #ifdef CBM_GLOBALTB_TOF_3D_CUBOIDS
527  new (&fXBins[i]) XBin;
528 #else //CBM_GLOBALTB_TOF_3D_CUBOIDS
529  new (&fXBins[i]) XBin(nofTBins);
530 #endif //CBM_GLOBALTB_TOF_3D_CUBOIDS
531  }
532 
533  void Clear() {
534  for (int i = 0; i < fNofXBins; ++i)
535  fXBins[i].Clear();
536  }
537 };
538 
539 struct ZBin {
540  int fNofYBins;
541  YBin* fYBins;
542 
543  ZBin(int nofYBins, int nofXBins, int nofTBins)
544  : fNofYBins(nofYBins)
545  , fYBins(
546  reinterpret_cast<YBin*>(new unsigned char[fNofYBins * sizeof(YBin)])) {
547  for (int i = 0; i < fNofYBins; ++i)
548  new (&fYBins[i]) YBin(nofXBins, nofTBins);
549  }
550 
551  void Clear() {
552  for (int i = 0; i < fNofYBins; ++i)
553  fYBins[i].Clear();
554  }
555 };
556 
557 CbmGlobalTrackingTofGeometry::CbmGlobalTrackingTofGeometry()
558  : fC(0)
559  , fPDG(211)
560  , fChi2Cut(50)
561  ,
562 #ifdef CBM_GLOBALTB_TOF_3D_CUBOIDS
563  fCuboids()
564  ,
565 #endif //CBM_GLOBALTB_TOF_3D_CUBOIDS
566  fNofTBins(1000)
567  , fNofXBins(60)
568  , fNofYBins(60)
569  , fNofZBins(5)
570  , fZBins(0)
571  , fTBinSize(100)
572  , fMinX(1000000)
573  , fMaxX(-1000000)
574  , fXBinSize(0)
575  , fMinY(1000000)
576  , fMaxY(-1000000)
577  , fYBinSize(0)
578  , fMinZ(1000000)
579  , fMaxZ(-1000000)
580  , fZBinSize(0)
581  , fStartTime(0)
582  , fEndTime(0)
583  , fTofHits(0) //, fPropagator()
584 {
585 }
586 
587 CbmGlobalTrackingTofGeometry::~CbmGlobalTrackingTofGeometry() {
588 #ifdef CBM_GLOBALTB_TOF_3D_CUBOIDS
589  for (list<Cuboid*>::iterator i = fCuboids.begin(); i != fCuboids.end(); ++i)
590  delete *i;
591 #endif //CBM_GLOBALTB_TOF_3D_CUBOIDS
592 }
593 
594 static void
595 FindGeoChild(TGeoNode* node, const char* name, list<TGeoNode*>& results) {
596  Int_t nofChildren = node->GetNdaughters();
597 
598  for (Int_t i = 0; i < nofChildren; ++i) {
599  TGeoNode* child = node->GetDaughter(i);
600  TString childName(child->GetName());
601 
602  if (childName.Contains(name, TString::kIgnoreCase))
603  results.push_back(child);
604  }
605 }
606 
607 bool CbmGlobalTrackingTofGeometry::Read() {
608  fC = 1.e-7 * TMath::C();
609  fPropagator = CbmLitToolFactory::CreateTrackPropagator("lit");
610  fLinePropagator = CbmLitToolFactory::CreateTrackPropagator("line");
611  fFilter = CbmLitToolFactory::CreateTrackUpdate("kalman");
612  TGeoNavigator* pNavigator = gGeoManager->GetCurrentNavigator();
613  gGeoManager->cd("/cave_1");
614  list<TGeoNode*> tofNodes;
615  FindGeoChild(gGeoManager->GetCurrentNode(), "tof", tofNodes);
616 
617  if (tofNodes.empty()) return false;
618 
619  TGeoNode* tofNode = tofNodes.front();
620  pNavigator->CdDown(tofNode);
621 
622  list<TGeoNode*> tofModules;
623  FindGeoChild(tofNode, "module", tofModules);
624 
625  for (list<TGeoNode*>::iterator i = tofModules.begin(); i != tofModules.end();
626  ++i) {
627  TGeoNode* tofModule = *i;
628  const char* modName = tofModule->GetName();
629  const char* firstUnd = strchr(modName, '_');
630  const char* lastUnd = strrchr(modName, '_');
631 
632  if (0 == firstUnd || 0 == lastUnd) continue;
633 
634  //int modType = atoi(firstUnd + 1); (VF) unused
635  //int modNum = atoi(lastUnd + 1); (VF) unused
636 
637  pNavigator->CdDown(tofModule);
638  list<TGeoNode*> tofGasBoxes;
639  FindGeoChild(tofModule, "gas_box", tofGasBoxes);
640 
641  for (list<TGeoNode*>::iterator j = tofGasBoxes.begin();
642  j != tofGasBoxes.end();
643  ++j) {
644  TGeoNode* tofGasBox = *j;
645  pNavigator->CdDown(tofGasBox);
646  list<TGeoNode*> tofModuleCounters;
647  FindGeoChild(tofGasBox, "counter", tofModuleCounters);
648 
649  for (list<TGeoNode*>::iterator k = tofModuleCounters.begin();
650  k != tofModuleCounters.end();
651  ++k) {
652  TGeoNode* tofModuleCounter = *k;
653  const char* counterName = tofModuleCounter->GetName();
654 
655  const char* lastUnd2 = strrchr(counterName, '_');
656 
657  if (0 == lastUnd2) continue;
658 
659  //int counterNum = atoi(lastUnd2 + 1); (VF) unused
660 
661  pNavigator->CdDown(tofModuleCounter);
662  TGeoVolume* tofModuleCounterVol = gGeoManager->GetCurrentVolume();
663  const TGeoBBox* tofModuleCounterShape =
664  static_cast<const TGeoBBox*>(tofModuleCounterVol->GetShape());
665  Double_t halfwidth = tofModuleCounterShape->GetDX();
666  Double_t halfheight = tofModuleCounterShape->GetDY();
667  Double_t halfthickness = tofModuleCounterShape->GetDZ();
668 
669 #ifdef CBM_GLOBALTB_TOF_3D_CUBOIDS
670  Cuboid* cuboid =
671  new Cuboid({2 * halfwidth, 2 * halfheight, 2 * halfthickness});
672  cuboid->fFullModId = CbmTofAddress::GetModFullId(
673  CbmTofAddress::GetUniqueAddress(modNum, counterNum, 0, 0, modType));
674 #endif //CBM_GLOBALTB_TOF_3D_CUBOIDS
675 
676  /*Int_t modId0 = CbmTofAddress::GetModFullId(CbmTofAddress::GetUniqueAddress(modNum, counterNum, 0, 0, modType));
677  Int_t modId1 = CbmTofAddress::GetModFullId(CbmTofAddress::GetUniqueAddress(modNum, counterNum, 1, 0, modType));
678  Int_t modId2 = CbmTofAddress::GetModFullId(CbmTofAddress::GetUniqueAddress(modNum, counterNum, 2, 0, modType));
679  Int_t modId3 = CbmTofAddress::GetModFullId(CbmTofAddress::GetUniqueAddress(modNum, counterNum, 3, 0, modType));
680  Int_t modId0_1 = CbmTofAddress::GetModFullId(CbmTofAddress::GetUniqueAddress(modNum, counterNum, 0, 1, modType));
681  Int_t modId1_1 = CbmTofAddress::GetModFullId(CbmTofAddress::GetUniqueAddress(modNum, counterNum, 1, 1, modType));
682  Int_t modId2_1 = CbmTofAddress::GetModFullId(CbmTofAddress::GetUniqueAddress(modNum, counterNum, 2, 1, modType));
683  Int_t modId3_1 = CbmTofAddress::GetModFullId(CbmTofAddress::GetUniqueAddress(modNum, counterNum, 3, 1, modType));*/
684 
685  for (int t = 0; t < 2; ++t) // t means thickness
686  {
687  for (int h = 0; h < 2; ++h) // h means height
688  {
689  for (int w = 0; w < 2; ++w) // w means width
690  {
691  Double_t localCoords[3] = {w > 0 ? halfwidth : -halfwidth,
692  h > 0 ? halfheight : -halfheight,
693  t > 0 ? halfthickness
694  : -halfthickness};
695  Double_t globalCoords[3];
696  gGeoManager->LocalToMaster(localCoords, globalCoords);
697 #ifdef CBM_GLOBALTB_TOF_3D_CUBOIDS
698  cuboid->vertices[t][h][w] = {
699  globalCoords[0], globalCoords[1], globalCoords[2]};
700 #endif //CBM_GLOBALTB_TOF_3D_CUBOIDS
701 
702  if (fMinX > globalCoords[0]) fMinX = globalCoords[0];
703 
704  if (fMaxX < globalCoords[0]) fMaxX = globalCoords[0];
705 
706  if (fMinY > globalCoords[1]) fMinY = globalCoords[1];
707 
708  if (fMaxY < globalCoords[1]) fMaxY = globalCoords[1];
709 
710  if (fMinZ > globalCoords[2]) fMinZ = globalCoords[2];
711 
712  if (fMaxZ < globalCoords[2]) fMaxZ = globalCoords[2];
713  }
714  }
715  }
716 
717 #ifdef CBM_GLOBALTB_TOF_3D_CUBOIDS
718  cuboid->Calc(fNofTBins);
719  fCuboids.push_back(cuboid);
720 #endif //CBM_GLOBALTB_TOF_3D_CUBOIDS
721  pNavigator->CdUp();
722  }
723 
724  pNavigator->CdUp();
725  }
726 
727  pNavigator->CdUp();
728  }
729 
730  pNavigator->CdUp();
731 
732  fXBinSize = (fMaxX - fMinX) / fNofXBins;
733  fYBinSize = (fMaxY - fMinY) / fNofYBins;
734  fZBinSize = (fMaxZ - fMinZ) / fNofZBins;
735 
736  // Create search bins
737  fZBins = reinterpret_cast<ZBin*>(new unsigned char[fNofZBins * sizeof(ZBin)]);
738 
739  for (int i = 0; i < fNofZBins; ++i)
740  new (&fZBins[i]) ZBin(fNofYBins, fNofXBins, fNofTBins);
741 
742 #ifdef CBM_GLOBALTB_TOF_3D_CUBOIDS
743  // Create references to the existing cuboids
744  for (list<Cuboid*>::iterator iter = fCuboids.begin(); iter != fCuboids.end();
745  ++iter) {
746  Cuboid* cuboid = *iter;
747  int zIndMin = GetZInd(cuboid->minZ);
748  int zIndMax = GetZInd(cuboid->maxZ);
749  int yIndMin = GetYInd(cuboid->minY);
750  int yIndMax = GetYInd(cuboid->maxY);
751  int xIndMin = GetXInd(cuboid->minX);
752  int xIndMax = GetXInd(cuboid->maxX);
753 
754  for (int zInd = zIndMin; zInd <= zIndMax; ++zInd) {
755  ZBin& zBin = fZBins[zInd];
756 
757  for (int yInd = yIndMin; yInd <= yIndMax; ++yInd) {
758  YBin& yBin = zBin.fYBins[yInd];
759 
760  for (int xInd = xIndMin; xInd <= xIndMax; ++xInd) {
761  XBin& xBin = yBin.fXBins[xInd];
762  Cuboid cuboid2 = {fXBinSize, fYBinSize, fZBinSize};
763 
764  for (int i = 0; i < 2; ++i) {
765  for (int j = 0; j < 2; ++j) {
766  for (int k = 0; k < 2; ++k)
767  cuboid2.vertices[i][j][k] = {fMinX + (xInd + k) * fXBinSize,
768  fMinY + (yInd + j) * fYBinSize,
769  fMinZ + (zInd + i) * fZBinSize};
770  }
771  }
772 
773  cuboid2.Calc();
774 
775  if (Intersect(*cuboid, cuboid2)) xBin.fCuboids.push_back(cuboid);
776  }
777  }
778  }
779  } // for (list<Cuboid*>::iterator i = fCuboids.begin(); i != fCuboids.end(); ++i)
780 #endif //CBM_GLOBALTB_TOF_3D_CUBOIDS
781 
782  return true;
783 }
784 
785 void CbmGlobalTrackingTofGeometry::Clear() {
786  for (int i = 0; i < fNofZBins; ++i)
787  fZBins[i].Clear();
788 }
789 
790 int globalNofHits = 0;
791 int globalNofHitsT = 0;
792 int globalNofHitsM = 0;
793 
794 void CbmGlobalTrackingTofGeometry::Prepare(timetype startTime) {
795  Clear();
796  fStartTime = startTime;
797  fEndTime = fStartTime + fNofTBins * fTBinSize;
798 
799  int nofHits = fTofHits->GetEntriesFast();
800  globalNofHits += nofHits;
801 
802  for (int i = 0; i < nofHits; ++i) {
803  const CbmTofHit* hit = static_cast<const CbmTofHit*>(fTofHits->At(i));
804  scaltype z = hit->GetZ();
805 
806  if (z < fMinZ || z > fMaxZ) continue;
807 
808  int zInd = GetZInd(z);
809  ZBin& zBin = fZBins[zInd];
810  scaltype y = hit->GetY();
811 
812  if (y < fMinY || y > fMaxY) continue;
813 
814  int yInd = GetYInd(y);
815  YBin& yBin = zBin.fYBins[yInd];
816  scaltype x = hit->GetX();
817 
818  if (x < fMinX || x > fMaxX) continue;
819 
820  int xInd = GetXInd(x);
821  XBin& xBin = yBin.fXBins[xInd];
822 
823 #ifdef CBM_GLOBALTB_TOF_3D_CUBOIDS
824  if (xBin.fCuboids.empty()) continue;
825 #endif //CBM_GLOBALTB_TOF_3D_CUBOIDS
826 
827  timetype t = hit->GetTime();
828 
829  if (t < fStartTime || t > fEndTime) continue;
830 
831  ++globalNofHitsT;
832 
833  int tInd = GetTInd(t);
834 
835 #ifdef CBM_GLOBALTB_TOF_3D_CUBOIDS
836  Int_t address = hit->GetAddress();
837  Int_t moduleId = CbmTofAddress::GetModFullId(address);
838 
839  for (list<Cuboid*>::iterator j = xBin.fCuboids.begin();
840  j != xBin.fCuboids.end();
841  ++j) {
842  Cuboid* mrpc = *j;
843 
844  if (mrpc->fFullModId == moduleId) {
845  mrpc->fTBins[tInd].fHitInds.push_back(i);
846  ++globalNofHitsM;
847  }
848  }
849 #else //CBM_GLOBALTB_TOF_3D_CUBOIDS
850  TBin& tBin = xBin.fTBins[tInd];
851  tBin.fHitInds.push_back(i);
852 #endif //CBM_GLOBALTB_TOF_3D_CUBOIDS
853  }
854 }
855 
856 int nofMRPCIntersections = 0;
857 int nofToFIntersections = 0;
858 int nofMRPCIntersectionsT = 0;
859 
860 void CbmGlobalTrackingTofGeometry::Find(FairTrackParam& trackParams,
861  timetype trackTime,
862  timetype errT,
863  Int_t& tofHitInd,
864  Double_t& length)
865 //void CbmGlobalTrackingTofGeometry::Find(scaltype x0, scaltype errX, scaltype y0, scaltype errY, scaltype z0, scaltype t0, scaltype errT,
866 //scaltype tx, scaltype errTx, scaltype ty, scaltype errTy, Int_t& tofHitInd)
867 //void CbmGlobalTrackingTofGeometry::Find(scaltype x0, scaltype errXSq, scaltype y0, scaltype errYSq, scaltype z0, scaltype t0, scaltype errT,
868 //scaltype tx, scaltype errTx, scaltype ty, scaltype errTy, Int_t& tofHitInd)
869 {
870  tofHitInd = -1;
871  Double_t stsTrackLength = length;
872  length = 0;
873  //double x0 = trackParams.GetX();
874 
875  //if (x0 < fMinX || x0 > fMaxX)
876  //return;
877 
878  //double y0 = trackParams.GetY();
879 
880  //if (y0 < fMinY || y0 > fMaxY)
881  //return;
882 
883  //double z0 = trackParams.GetZ();
884  //double t0 = fStartTime;//trackParams.GetTime();
885  //double tx = trackParams.GetTx();
886  //double ty = trackParams.GetTy();
887 
888  //double deltaZ1 = fMinZ - z0;
889  //double x1 = x0 + tx * deltaZ1;
890  //double y1 = y0 + ty * deltaZ1;
891  double z1 = fMinZ;
892  CbmTrackParam cbmTrackParams;
893  cbmTrackParams.Set(trackParams, trackTime, errT);
894  CbmLitTrackParam litTrackParams;
895  CbmLitConverterFairTrackParam::FairTrackParamToCbmLitTrackParam(
896  &cbmTrackParams, &litTrackParams);
897 
898  if (fPropagator->Propagate(&litTrackParams, z1, fPDG, 0, &length)
899  == kLITERROR) {
900  length = 0;
901  return;
902  }
903 
904  double x1 = litTrackParams.GetX();
905  double y1 = litTrackParams.GetY();
906 
907  //if (fMinX > x1 || x1 > fMaxX || fMinY > y1 || y1 > fMaxY)
908  //return;
909 
910  // First check if track exits from the ToF XYZ manifold from the back face.
911  double deltaZMax = fMaxZ - z1;
912  double tx = litTrackParams.GetTx();
913  double xMax = x1 + tx * deltaZMax;
914  double ty = litTrackParams.GetTy();
915  double yMax = y1 + ty * deltaZMax;
916  double normLen = TMath::Sqrt(1 + tx * tx + ty * ty);
917  //double t1 = trackTime + z1 * normLen / fC;
918  double t1 = trackTime + (stsTrackLength + length) / fC;
919  //double t1 = litTrackParams.GetTime();
920  double zMax;
921 
922  map<int, map<int, map<int, double>>> inds;
923 
924  /*for (int i = 0; i < 25; ++i)
925  {
926  double val = litTrackParams.GetCovariance(i);
927  int qq = 0;
928  }*/
929 
930  double deltaX = 4 * TMath::Sqrt(litTrackParams.GetCovariance(0));
931  double deltaY = 4 * TMath::Sqrt(litTrackParams.GetCovariance(6));
932  double deltaT = 4 * errT;
933  //double deltaT = 4 * litTrackParams.GetTimeError();
934  //Find(x1 - deltaX, y1 - deltaY, z1, tx, ty, inds);
935  //Find(x1 + deltaX, y1 - deltaY, z1, tx, ty, inds);
936  //Find(x1 - deltaX, y1 + deltaY, z1, tx, ty, inds);
937  //Find(x1 + deltaX, y1 + deltaY, z1, tx, ty, inds);
938 
939 
940  if (fMinX <= xMax && xMax <= fMaxX && fMinY <= yMax && yMax <= fMaxY)
941  zMax = fMaxZ;
942  else if (ty > 0 && (fMaxY - y1) / ty < fMaxZ - z1
943  && fMinX <= x1 + tx * (fMaxY - y1) / ty
944  && x1 + tx * (fMaxY - y1) / ty <= fMaxX) {
945  zMax = z1 + (fMaxY - y1) / ty;
946  yMax = fMaxY;
947  xMax = x1 + tx * (fMaxY - y1) / ty;
948  } else if (ty < 0 && (fMinY - y1) / ty < fMaxZ - z1
949  && fMinX <= x1 + tx * (fMinY - y1) / ty
950  && x1 + tx * (fMinY - y1) / ty <= fMaxX) {
951  zMax = z1 + (fMinY - y1) / ty;
952  yMax = fMinY;
953  xMax = x1 + tx * (fMinY - y1) / ty;
954  } else if (tx > 0 && (fMaxX - x1) / tx < fMaxZ - z1
955  && fMinY <= y1 + ty * (fMaxX - x1) / tx
956  && y1 + ty * (fMaxX - x1) / tx <= fMaxY) {
957  zMax = z1 + (fMaxX - x1) / tx;
958  yMax = y1 + ty * (fMaxX - x1) / tx;
959  xMax = fMaxX;
960  } else if (tx < 0 && (fMinX - x1) / tx < fMaxZ - z1
961  && fMinY <= y1 + ty * (fMinX - x1) / tx
962  && y1 + ty * (fMinX - x1) / tx <= fMaxY) {
963  zMax = z1 + (fMinX - x1) / tx;
964  yMax = y1 + ty * (fMinX - x1) / tx;
965  xMax = fMinX;
966  } else
967  return;
968 
969  ++nofToFIntersections;
970 
971 #ifdef CBM_GLOBALTB_TOF_3D_CUBOIDS
972  set<const Cuboid*> cuboidSet;
973 #else //CBM_GLOBALTB_TOF_3D_CUBOIDS
974  double minChi2 = std::numeric_limits<double>::max();
975  const CbmTofHit* mergingHit = 0;
976  double deltaLength = 0;
977  Line line = {x1, y1, z1, tx / normLen, ty / normLen, 1 / normLen};
978 #endif //CBM_GLOBALTB_TOF_3D_CUBOIDS
979 
980  set<const TBin*> neighbourhood;
981 
982  int zMinInd = GetZInd(fMinZ);
983  int zMaxInd = GetZInd(zMax);
984 
985  for (int zInd = zMinInd; zInd <= zMaxInd; ++zInd) {
986  const ZBin& zBin = fZBins[zInd];
987  double startZ = fMinZ + zInd * fZBinSize;
988  double stopZ = startZ + fZBinSize;
989  int yIndDelta = ty < 0 ? -1 : 1;
990  double startY = y1 + (startZ - z1) * ty; // - yIndDelta * deltaY;
991  double stopY = y1 + (stopZ - z1) * ty; // + yIndDelta * deltaY;
992  int yStartInd = GetYInd(startY);
993  int yStopInd = GetYInd(stopY);
994 
995  for (int yInd = yStartInd; true; yInd += yIndDelta) {
996  //const YBin& yBin = zBin.fYBins[yInd]; (VF) unused
997  double startZy;
998  double stopZy;
999 
1000  if (0 == ty) {
1001  startZy = startZ;
1002  stopZy = stopZ;
1003  } else if (0 > ty) {
1004  startZy = z1 + (fMinY + (yInd + 1) * fYBinSize - y1) / ty;
1005  stopZy = z1 + (fMinY + yInd * fYBinSize - y1) / ty;
1006  } else {
1007  startZy = z1 + (fMinY + yInd * fYBinSize - y1) / ty;
1008  stopZy = z1 + (fMinY + (yInd + 1) * fYBinSize - y1) / ty;
1009  }
1010 
1011  if (startZy < startZ) startZy = startZ;
1012 
1013  if (stopZy > stopZ) stopZy = stopZ;
1014 
1015  //int xIndDelta = tx < 0 ? -1 : 1; (VF) unused
1016  double startX = x1 + (startZy - z1) * tx; // - xIndDelta * deltaX;
1017  double stopX = x1 + (stopZy - z1) * tx; // + xIndDelta * deltaX;
1018  //int xStartInd = GetXInd(startX); (VF) unused
1019  //int xStopInd = GetXInd(stopX); (VF) unused
1020  //#ifndef CBM_GLOBALTB_TOF_3D_CUBOIDS
1021  double extLen = ((startZy + stopZy) / 2 - z1) * normLen;
1022  double extT = t1 + extLen / fC;
1023  //int tInd = (extT - fStartTime) / fTBinSize;
1024  int tLowInd = (extT - deltaT - fStartTime) / fTBinSize;
1025  int tHighInd = (extT + deltaT - fStartTime) / fTBinSize;
1026 
1027  if (tHighInd >= 0 || tLowInd < fNofTBins) {
1028  if (tLowInd < 0) tLowInd = 0;
1029 
1030  if (tHighInd >= fNofTBins) tHighInd = fNofTBins - 1;
1031 
1032  double startY2 = startY < stopY ? startY - deltaY : stopY - deltaY;
1033  double stopY2 = startY < stopY ? stopY + deltaY : startY + deltaY;
1034  int yStartInd2 = GetYInd(startY2);
1035  int yStopInd2 = GetYInd(stopY2);
1036  double startX2 = startX < stopX ? startX - deltaX : stopX - deltaX;
1037  double stopX2 = startX < stopX ? stopX + deltaX : startX + deltaX;
1038  int xStartInd2 = GetXInd(startX2);
1039  int xStopInd2 = GetXInd(stopX2);
1040 
1041  for (int yInd2 = yStartInd2; yInd2 <= yStopInd2; ++yInd2) {
1042  const YBin& yBin2 = zBin.fYBins[yInd2];
1043 
1044  for (int xInd2 = xStartInd2; xInd2 <= xStopInd2; ++xInd2) {
1045  const XBin& xBin2 = yBin2.fXBins[xInd2];
1046 
1047  for (int tInd2 = tLowInd; tInd2 <= tHighInd; ++tInd2)
1048  neighbourhood.insert(&xBin2.fTBins[tInd2]);
1049  }
1050  }
1051  }
1052 
1053  if (yInd == yStopInd) break;
1054  }
1055  }
1056 
1057  for (set<const TBin*>::const_iterator tBinIter = neighbourhood.begin();
1058  tBinIter != neighbourhood.end();
1059  ++tBinIter) {
1060  const TBin* tBin = *tBinIter;
1061 
1062  for (list<Int_t>::const_iterator hitIndIter = tBin->fHitInds.begin();
1063  hitIndIter != tBin->fHitInds.end();
1064  ++hitIndIter) {
1065  Int_t hitInd = *hitIndIter;
1066  const CbmTofHit* hit =
1067  static_cast<const CbmTofHit*>(fTofHits->At(hitInd));
1068  double L01Sq = (hit->GetX() - x1) * (hit->GetX() - x1)
1069  + (hit->GetY() - y1) * (hit->GetY() - y1)
1070  + (hit->GetZ() - z1) * (hit->GetZ() - z1);
1071  //double L01 = TMath::Sqrt(L01Sq);
1072  double L02 = (hit->GetX() - x1) * line.cosX
1073  + (hit->GetY() - y1) * line.cosY
1074  + (hit->GetZ() - z1) * line.cosZ;
1075  double extT2 = t1 + L02 / fC;
1076  /*double x = x1 + L02 * line.cosX;
1077  double y = y1 + L02 * line.cosY;
1078  double z = z1 + L02 * line.cosZ;
1079  double extT2 = t1 + TMath::Sqrt(x * x + y * y + z * z) / fC;*/
1080  double chi2 =
1081  (L01Sq - L02 * L02)
1082  / (litTrackParams.GetCovariance(0) + hit->GetDx() * hit->GetDx()
1083  + litTrackParams.GetCovariance(5) + hit->GetDy() * hit->GetDy())
1084  + (hit->GetTime() - extT2) * (hit->GetTime() - extT2)
1085  / (/*litTrackParams.GetTimeError() * litTrackParams.GetTimeError()*/
1086  errT * errT
1087  + hit->GetTimeError() * hit->GetTimeError());
1088 
1089  if (chi2 < minChi2) {
1090  tofHitInd = hitInd;
1091  mergingHit = hit;
1092  minChi2 = chi2;
1093  deltaLength = L02;
1094  }
1095  }
1096  }
1097 
1098  if (minChi2 > fChi2Cut) {
1099  tofHitInd = -1;
1100  length = 0;
1101  } else {
1102  length += deltaLength;
1103  fLinePropagator->Propagate(&litTrackParams, mergingHit->GetZ(), fPDG);
1104  CbmLitPixelHit litHit;
1105  CbmLitConverter::CbmPixelHitToCbmLitPixelHit(
1106  mergingHit, tofHitInd, &litHit);
1107  double mergingChi2 = 0;
1108  fFilter->Update(&litTrackParams, &litHit, mergingChi2);
1109  CbmLitConverterFairTrackParam::CbmLitTrackParamToFairTrackParam(
1110  &litTrackParams, &trackParams);
1111  }
1112 
1113 #if 0
1114  if (tInd >= 0 && tInd < fNofTBins)
1115  {
1116 //#endif//CBM_GLOBALTB_TOF_3D_CUBOIDS
1117  for (int xInd = xStartInd; true; xInd += xIndDelta)
1118  {
1119  const XBin& xBin = yBin.fXBins[xInd];
1120 
1121 #ifdef CBM_GLOBALTB_TOF_3D_CUBOIDS
1122  for (list<Cuboid*>::const_iterator i = xBin.fCuboids.begin(); i != xBin.fCuboids.end(); ++i)
1123  cuboidSet.insert(*i);
1124 #else //CBM_GLOBALTB_TOF_3D_CUBOIDS
1125  const TBin& tBin = xBin.fTBins[tInd];
1126 
1127  for (list<Int_t>::const_iterator hitIndIter = tBin.fHitInds.begin(); hitIndIter != tBin.fHitInds.end(); ++hitIndIter)
1128  {
1129  Int_t hitInd = *hitIndIter;
1130  const CbmTofHit* hit = static_cast<const CbmTofHit*> (fTofHits->At(hitInd));
1131  double L01Sq = (hit->GetX() - x1) * (hit->GetX() - x1) + (hit->GetY() - y1) * (hit->GetY() - y1) + (hit->GetZ() - z1) * (hit->GetZ() - z1);
1132  //double L01 = TMath::Sqrt(L01Sq);
1133  double L02 = (hit->GetX() - x1) * line.cosX + (hit->GetY() - y1) * line.cosY + (hit->GetZ() - z1) * line.cosZ;
1134  double extT2 = t1 + L02 / fC;
1135  double chi2 = (L01Sq - L02 * L02) / (litTrackParams.GetCovariance(0) + hit->GetDx() * hit->GetDx() + litTrackParams.GetCovariance(6) + hit->GetDy() * hit->GetDy()) +
1136  (hit->GetTime() - extT2) * (hit->GetTime() - extT2) / (litTrackParams.GetTimeError() * litTrackParams.GetTimeError() + hit->GetTimeError() * hit->GetTimeError());
1137 
1138  if (chi2 < minChi2)
1139  {
1140  tofHitInd = hitInd;
1141  minChi2 = chi2;
1142  }
1143  }
1144 #endif //CBM_GLOBALTB_TOF_3D_CUBOIDS
1145 
1146  if (xInd == xStopInd)
1147  break;
1148  }
1149 #ifndef CBM_GLOBALTB_TOF_3D_CUBOIDS
1150  }
1151 #endif //CBM_GLOBALTB_TOF_3D_CUBOIDS
1152 
1153  if (yInd == yStopInd)
1154  break;
1155  }
1156  }
1157 #endif // 0
1158 
1159  /*double minChi2 = std::numeric_limits<double>::max();
1160  Line line = { x1, y1, z1, tx / normLen, ty / normLen, 1 / normLen };
1161 
1162  for (map<int, map<int, map<int, double> > >::const_iterator i = inds.begin(); i != inds.end(); ++i)
1163  {
1164  const map<int, map<int, double> >& yInds = i->second;
1165 
1166  for (map<int, map<int, double> >::const_iterator j = yInds.begin(); j != yInds.end(); ++j)
1167  {
1168  const map<int, double>& xInds = j->second;
1169 
1170  for (map<int, double>::const_iterator k = xInds.begin(); k != xInds.end(); ++k)
1171  {
1172  const XBin& xBin = fZBins[i->first].fYBins[j->first].fXBins[k->first];
1173  double extT = t1 + k->second / fC;
1174  int tInd = (extT - fStartTime) / fTBinSize;
1175  const TBin& tBin = xBin.fTBins[tInd];
1176 
1177  for (list<Int_t>::const_iterator hitIndIter = tBin.fHitInds.begin(); hitIndIter != tBin.fHitInds.end(); ++hitIndIter)
1178  {
1179  Int_t hitInd = *hitIndIter;
1180  const CbmTofHit* hit = static_cast<const CbmTofHit*> (fTofHits->At(hitInd));
1181  double L01Sq = (hit->GetX() - x1) * (hit->GetX() - x1) + (hit->GetY() - y1) * (hit->GetY() - y1) + (hit->GetZ() - z1) * (hit->GetZ() - z1);
1182  //double L01 = TMath::Sqrt(L01Sq);
1183  double L02 = (hit->GetX() - x1) * line.cosX + (hit->GetY() - y1) * line.cosY + (hit->GetZ() - z1) * line.cosZ;
1184  double extT2 = t1 + L02 / fC;
1185  double chi2 = (L01Sq - L02 * L02) / (litTrackParams.GetCovariance(0) + hit->GetDx() * hit->GetDx() + litTrackParams.GetCovariance(6) + hit->GetDy() * hit->GetDy()) +
1186  (hit->GetTime() - extT2) * (hit->GetTime() - extT2) / (errT * errT + hit->GetTimeError() * hit->GetTimeError());
1187 
1188  if (chi2 < minChi2) {
1189  tofHitInd = hitInd;
1190  minChi2 = chi2;
1191  }
1192  }
1193  }
1194  }
1195  }*/
1196 
1197 #ifdef CBM_GLOBALTB_TOF_3D_CUBOIDS
1198  if (cuboidSet.empty()) return;
1199 
1200  bool inMRPC = false;
1201  bool inMRPCT = false;
1202 
1203  double minChi2 = std::numeric_limits<double>::max();
1204  double normLen = TMath::Sqrt(1 + tx * tx + ty * ty);
1205  Line line = {x1, y1, z1, tx / normLen, ty / normLen, 1 / normLen};
1206 
1207  for (set<const Cuboid*>::const_iterator i = cuboidSet.begin();
1208  i != cuboidSet.end();
1209  ++i) {
1210  const Cuboid* cuboid = *i;
1211  const Rectangle& face1 = cuboid->faces[0][0];
1212  const Rectangle& face2 = cuboid->faces[0][1];
1213  double result[3];
1214 
1215  if (Intersect(face1, line, result) || Intersect(face2, line, result)) {
1216  inMRPC = true;
1217  double extT = t0 + result[2] / fC;
1218  int tInd = (extT - fStartTime) / fTBinSize;
1219 
1220  if (tInd < 0 || tInd >= fNofTBins) continue;
1221 
1222  inMRPCT = true;
1223 
1224  //const TBin& tBin = xBin.f
1225  double extX = x1 + line.cosX * result[2];
1226  double extY = y1 + line.cosY * result[2];
1227  double extZ = z1 + line.cosZ * result[2];
1228  TBin& tBin = cuboid->fTBins[tInd];
1229 
1230  for (list<Int_t>::const_iterator hitIndIter = tBin.fHitInds.begin();
1231  hitIndIter != tBin.fHitInds.end();
1232  ++hitIndIter) {
1233  Int_t hitInd = *hitIndIter;
1234  const CbmTofHit* hit =
1235  static_cast<const CbmTofHit*>(fTofHits->At(hitInd));
1236  double chi2 =
1237  (hit->GetX() - extX) * (hit->GetX() - extX)
1238  / (errX * errX + hit->GetDx() * hit->GetDx())
1239  + (hit->GetY() - extY) * (hit->GetY() - extY)
1240  / (errY * errY + hit->GetDy() * hit->GetDy())
1241  + (hit->GetTime() - extT) * (hit->GetTime() - extT)
1242  / (litTrackParams.GetTimeError() * litTrackParams.GetTimeError()
1243  + hit->GetTimeError() * hit->GetTimeError());
1244 
1245  if (chi2 < minChi2) {
1246  tofHitInd = hitInd;
1247  minChi2 = chi2;
1248  }
1249  }
1250  }
1251  }
1252 
1253  if (inMRPC) ++nofMRPCIntersections;
1254 
1255  if (inMRPCT) ++nofMRPCIntersectionsT;
1256 #endif //CBM_GLOBALTB_TOF_3D_CUBOIDS
1257 }
XBin::fNofTBins
int fNofTBins
Definition: CbmGlobalTrackingTof.cxx:504
CbmLitToolFactory::CreateTrackPropagator
static TrackPropagatorPtr CreateTrackPropagator(const string &name)
Create track propagation tool by name.
Definition: CbmLitToolFactory.cxx:58
CbmHit::GetZ
Double_t GetZ() const
Definition: CbmHit.h:70
Plane::z0
double z0
Definition: CbmGlobalTrackingTof.cxx:114
CbmGlobalTrackingTofGeometry::fC
scaltype fC
Definition: CbmGlobalTrackingTof.h:120
fPropagator
TrackPropagatorPtr fPropagator
Definition: CbmGlobalTrackingTof.cxx:19
YBin::Clear
void Clear()
Definition: CbmGlobalTrackingTof.cxx:533
Line::cosX
double cosX
Definition: CbmGlobalTrackingTof.cxx:101
Direction
Definition: CbmGlobalTrackingTof.cxx:185
CbmGlobalTrackingTofGeometry::GetZInd
int GetZInd(scaltype z) const
Definition: CbmGlobalTrackingTof.h:81
Rectangle::s1
double s1
Definition: CbmGlobalTrackingTof.cxx:140
CbmGlobalTrackingTofGeometry::fMaxY
scaltype fMaxY
Definition: CbmGlobalTrackingTof.h:136
CbmGlobalTrackingTofGeometry::fNofTBins
int fNofTBins
Definition: CbmGlobalTrackingTof.h:126
CbmPixelHit::GetX
Double_t GetX() const
Definition: CbmPixelHit.h:83
CbmGlobalTrackingTof.h
CbmLitTrackParam::GetX
litfloat GetX() const
Definition: CbmLitTrackParam.h:53
scaltype
#define scaltype
Definition: CbmGlobalTrackingDefs.h:17
CbmPixelHit::GetY
Double_t GetY() const
Definition: CbmPixelHit.h:84
CbmGlobalTrackingTofGeometry::fEndTime
timetype fEndTime
Definition: CbmGlobalTrackingTof.h:142
Rectangle::s2
double s2
Definition: CbmGlobalTrackingTof.cxx:141
CbmPixelHit::GetDx
Double_t GetDx() const
Definition: CbmPixelHit.h:85
TrackPropagatorPtr
boost::shared_ptr< CbmLitTrackPropagator > TrackPropagatorPtr
Definition: CbmTofPtrTypes.h:23
CbmLitTrackParam
Data class for track parameters.
Definition: CbmLitTrackParam.h:29
i
int i
Definition: L1/vectors/P4_F32vec4.h:25
Line::y0
double y0
Definition: CbmGlobalTrackingTof.cxx:98
CbmGlobalTrackingTofGeometry::fMinX
scaltype fMinX
Definition: CbmGlobalTrackingTof.h:132
TrackUpdatePtr
boost::shared_ptr< CbmLitTrackUpdate > TrackUpdatePtr
Definition: CbmTofPtrTypes.h:26
CbmLitTrackParam::GetY
litfloat GetY() const
Definition: CbmLitTrackParam.h:54
CbmLitToolFactory.h
Tool factory for creation of littrack algorithms.
TBin::fHitInds
list< Int_t > fHitInds
Definition: CbmGlobalTrackingTof.cxx:192
nofToFIntersections
int nofToFIntersections
Definition: CbmGlobalTrackingTof.cxx:857
CbmLitTrackParam::GetCovariance
litfloat GetCovariance(int index) const
Definition: CbmLitTrackParam.h:60
CbmGlobalTrackingTofGeometry::CbmGlobalTrackingTofGeometry
CbmGlobalTrackingTofGeometry()
Definition: CbmGlobalTrackingTof.cxx:557
fLinePropagator
TrackPropagatorPtr fLinePropagator
Definition: CbmGlobalTrackingTof.cxx:20
Plane::cosX1
double cosX1
Definition: CbmGlobalTrackingTof.cxx:116
Rectangle::plane
Plane plane
Definition: CbmGlobalTrackingTof.cxx:139
fFilter
TrackUpdatePtr fFilter
Definition: CbmGlobalTrackingTof.cxx:21
CbmGlobalTrackingTofGeometry::fZBinSize
scaltype fZBinSize
Definition: CbmGlobalTrackingTof.h:140
ZBin::fNofYBins
int fNofYBins
Definition: CbmGlobalTrackingTof.cxx:540
globalNofHitsT
int globalNofHitsT
Definition: CbmGlobalTrackingTof.cxx:791
YBin::YBin
YBin(int nofXBins, int nofTBins)
Definition: CbmGlobalTrackingTof.cxx:521
CbmTrackParam::Set
void Set(const FairTrackParam &ftp, Double_t time=0., Double_t timeError=0.)
Definition: CbmTrackParam.cxx:12
CbmGlobalTrackingTofGeometry::~CbmGlobalTrackingTofGeometry
~CbmGlobalTrackingTofGeometry()
Definition: CbmGlobalTrackingTof.cxx:587
CbmPixelHit::GetDy
Double_t GetDy() const
Definition: CbmPixelHit.h:86
CbmGlobalTrackingTofGeometry::GetYInd
int GetYInd(scaltype y) const
Definition: CbmGlobalTrackingTof.h:70
CbmHit::GetTimeError
Double_t GetTimeError() const
Definition: CbmHit.h:76
CbmLitConverterFairTrackParam::CbmLitTrackParamToFairTrackParam
static void CbmLitTrackParamToFairTrackParam(const CbmLitTrackParam *litPar, FairTrackParam *par)
Definition: CbmLitConverterFairTrackParam.h:110
CbmGlobalTrackingTofGeometry::fMinZ
scaltype fMinZ
Definition: CbmGlobalTrackingTof.h:138
Direction::cosX
double cosX
Definition: CbmGlobalTrackingTof.cxx:186
FindGeoChild
static void FindGeoChild(TGeoNode *node, const char *name, list< TGeoNode * > &results)
Definition: CbmGlobalTrackingTof.cxx:595
Line::x0
double x0
Definition: CbmGlobalTrackingTof.cxx:97
Plane::cosX2
double cosX2
Definition: CbmGlobalTrackingTof.cxx:120
CbmGlobalTrackingTofGeometry::fMinY
scaltype fMinY
Definition: CbmGlobalTrackingTof.h:135
ZBin::fYBins
YBin * fYBins
Definition: CbmGlobalTrackingTof.cxx:541
XBin::Clear
void Clear()
Definition: CbmGlobalTrackingTof.cxx:509
CbmGlobalTrackingTofGeometry::GetTInd
int GetTInd(timetype t) const
Definition: CbmGlobalTrackingTof.h:92
h
Data class with information on a STS local track.
Line
Definition: CbmGlobalTrackingTof.cxx:96
Direction::cosY
double cosY
Definition: CbmGlobalTrackingTof.cxx:187
CbmGlobalTrackingTofGeometry::fPDG
Int_t fPDG
Definition: CbmGlobalTrackingTof.h:121
CbmGlobalTrackingTofGeometry::fTofHits
TClonesArray * fTofHits
Definition: CbmGlobalTrackingTof.h:143
CbmGlobalTrackingTofGeometry::fNofZBins
int fNofZBins
Definition: CbmGlobalTrackingTof.h:129
Point::z
double z
Definition: CbmGlobalTrackingTof.cxx:93
CbmGlobalTrackingTofGeometry::fMaxX
scaltype fMaxX
Definition: CbmGlobalTrackingTof.h:133
CbmHit::GetTime
Double_t GetTime() const
Definition: CbmHit.h:75
Plane::x0
double x0
Definition: CbmGlobalTrackingTof.cxx:112
CbmGlobalTrackingTofGeometry::fMaxZ
scaltype fMaxZ
Definition: CbmGlobalTrackingTof.h:139
CbmGlobalTrackingTofGeometry::Read
bool Read()
Definition: CbmGlobalTrackingTof.cxx:607
CbmHit::GetAddress
Int_t GetAddress() const
Definition: CbmHit.h:73
CbmGlobalTrackingTofGeometry::fChi2Cut
Double_t fChi2Cut
Definition: CbmGlobalTrackingTof.h:122
nofMRPCIntersectionsT
int nofMRPCIntersectionsT
Definition: CbmGlobalTrackingTof.cxx:858
Point
Definition: CbmGlobalTrackingTof.cxx:90
kLITERROR
@ kLITERROR
Definition: CbmLitEnums.h:25
Plane::cosY1
double cosY1
Definition: CbmGlobalTrackingTof.cxx:117
CbmTofAddress.h
CbmLitPixelHit
Base data class for pixel hits.
Definition: CbmLitPixelHit.h:22
YBin::fXBins
XBin * fXBins
Definition: CbmGlobalTrackingTof.cxx:519
Segment::b
Point b
Definition: CbmGlobalTrackingTof.cxx:108
TBin::Clear
void Clear()
Definition: CbmGlobalTrackingTof.cxx:196
Plane::y0
double y0
Definition: CbmGlobalTrackingTof.cxx:113
CbmPixelHit.h
CbmGlobalTrackingTofGeometry::Find
void Find(FairTrackParam &trackParams, timetype trackTime, timetype errT, Int_t &tofHitInd, Double_t &length)
Definition: CbmGlobalTrackingTof.cxx:860
Segment::a
Point a
Definition: CbmGlobalTrackingTof.cxx:107
CbmGlobalTrackingTofGeometry::fXBinSize
scaltype fXBinSize
Definition: CbmGlobalTrackingTof.h:134
Direction::cosZ
double cosZ
Definition: CbmGlobalTrackingTof.cxx:188
globalNofHitsM
int globalNofHitsM
Definition: CbmGlobalTrackingTof.cxx:792
Rectangle
Definition: CbmGlobalTrackingTof.cxx:138
ZBin::Clear
void Clear()
Definition: CbmGlobalTrackingTof.cxx:551
TBin::TBin
TBin()
Definition: CbmGlobalTrackingTof.cxx:194
YBin::fNofXBins
int fNofXBins
Definition: CbmGlobalTrackingTof.cxx:518
timetype
#define timetype
Definition: CbmGlobalTrackingDefs.h:18
CbmLitToolFactory::CreateTrackUpdate
static TrackUpdatePtr CreateTrackUpdate(const string &name)
Create track update tool by name.
Definition: CbmLitToolFactory.cxx:73
Segment
Definition: CbmGlobalTrackingTof.cxx:106
x
Double_t x
Definition: CbmMvdSensorDigiToHitTask.cxx:68
CbmGlobalTrackingTofGeometry::fTBinSize
scaltype fTBinSize
Definition: CbmGlobalTrackingTof.h:131
CbmTrackParam
Definition: CbmTrackParam.h:22
CbmGlobalTrackingTofGeometry::fStartTime
timetype fStartTime
Definition: CbmGlobalTrackingTof.h:141
Line::cosZ
double cosZ
Definition: CbmGlobalTrackingTof.cxx:103
y
Double_t y
Definition: CbmMvdSensorDigiToHitTask.cxx:68
CbmLitConverter.h
CbmGlobalTrackingTofGeometry::Prepare
void Prepare(timetype startTime)
Definition: CbmGlobalTrackingTof.cxx:794
XBin::XBin
XBin(int nofTBins)
Definition: CbmGlobalTrackingTof.cxx:506
Plane::cosY2
double cosY2
Definition: CbmGlobalTrackingTof.cxx:121
TBin
Definition: CbmGlobalTrackingTof.cxx:191
Line::cosY
double cosY
Definition: CbmGlobalTrackingTof.cxx:102
globalNofHits
int globalNofHits
Definition: CbmGlobalTrackingTof.cxx:790
CbmGlobalTrackingTofGeometry::fNofXBins
int fNofXBins
Definition: CbmGlobalTrackingTof.h:127
CbmGlobalTrackingTofGeometry::fNofYBins
int fNofYBins
Definition: CbmGlobalTrackingTof.h:128
CbmLitTrackParam::GetTy
litfloat GetTy() const
Definition: CbmLitTrackParam.h:57
nofMRPCIntersections
int nofMRPCIntersections
Definition: CbmGlobalTrackingTof.cxx:856
Point::y
double y
Definition: CbmGlobalTrackingTof.cxx:92
CbmTofHit
Definition: core/data/tof/CbmTofHit.h:26
Plane::cosZ2
double cosZ2
Definition: CbmGlobalTrackingTof.cxx:122
Plane
Definition: CbmGlobalTrackingTof.cxx:111
Point::x
double x
Definition: CbmGlobalTrackingTof.cxx:91
XBin::fTBins
TBin * fTBins
Definition: CbmGlobalTrackingTof.cxx:505
CbmTofAddress::GetUniqueAddress
static UInt_t GetUniqueAddress(UInt_t Sm, UInt_t Rpc, UInt_t Channel, UInt_t Side=0, UInt_t SmType=0)
Definition: CbmTofAddress.h:124
ZBin::ZBin
ZBin(int nofYBins, int nofXBins, int nofTBins)
Definition: CbmGlobalTrackingTof.cxx:543
CbmGlobalTrackingTofGeometry::fZBins
ZBin * fZBins
Definition: CbmGlobalTrackingTof.h:130
z1
Double_t z1[nSects1]
Definition: pipe_v16a_mvdsts100.h:6
Line::z0
double z0
Definition: CbmGlobalTrackingTof.cxx:99
max
friend F32vec4 max(const F32vec4 &a, const F32vec4 &b)
Definition: L1/vectors/P4_F32vec4.h:36
XBin
Definition: CbmGlobalTrackingTof.cxx:492
CbmGlobalTrackingTofGeometry::GetXInd
int GetXInd(scaltype x) const
Definition: CbmGlobalTrackingTof.h:59
YBin
Definition: CbmGlobalTrackingTof.cxx:517
CbmLitConverter::CbmPixelHitToCbmLitPixelHit
static void CbmPixelHitToCbmLitPixelHit(const CbmPixelHit *hit, Int_t index, CbmLitPixelHit *litHit)
Definition: CbmLitConverter.h:46
CbmLitTrackParam::GetTx
litfloat GetTx() const
Definition: CbmLitTrackParam.h:56
CbmTofAddress::GetModFullId
static Int_t GetModFullId(UInt_t address)
Definition: CbmTofAddress.h:111
CbmGlobalTrackingTofGeometry::fYBinSize
scaltype fYBinSize
Definition: CbmGlobalTrackingTof.h:137
ZBin
Definition: CbmGlobalTrackingTof.cxx:539
CbmLitConverterFairTrackParam::FairTrackParamToCbmLitTrackParam
static void FairTrackParamToCbmLitTrackParam(const FairTrackParam *par, CbmLitTrackParam *litPar)
Definition: CbmLitConverterFairTrackParam.h:40
CbmGlobalTrackingTofGeometry::Clear
void Clear()
Definition: CbmGlobalTrackingTof.cxx:785
Plane::cosZ1
double cosZ1
Definition: CbmGlobalTrackingTof.cxx:118