CbmRoot
LxCA.cxx
Go to the documentation of this file.
1 #include "LxCA.h"
2 #include "Lx.h"
3 #include <iostream>
4 #pragma GCC diagnostic ignored "-Weffc++"
5 
6 #ifdef CLUSTER_MODE
7 #include "kdtree++/kdtree.hpp"
8 #endif //CLUSTER_MODE
9 
10 #include <sys/time.h>
11 #ifdef BEST_SIX_POINTS
12 #include <set>
13 #endif //BEST_SIX_POINTS
14 
15 #include "base/CbmLitToolFactory.h"
16 #include "utils/CbmLitConverterFairTrackParam.h"
17 #include <cmath>
18 
19 using namespace std;
20 
21 scaltype tx_limits[LXSTATIONS];
22 scaltype tx_limits_sq[LXSTATIONS];
23 scaltype ty_limits[LXSTATIONS];
24 scaltype ty_limits_sq[LXSTATIONS];
25 scaltype x_disp_left_limits[LXSTATIONS];
26 scaltype x_disp_left_limits_sq[LXSTATIONS];
27 scaltype y_disp_left_limits[LXSTATIONS];
28 scaltype y_disp_left_limits_sq[LXSTATIONS];
29 scaltype x_disp_right_limits[LXSTATIONS];
30 scaltype x_disp_right_limits_sq[LXSTATIONS];
31 scaltype y_disp_right_limits[LXSTATIONS];
32 scaltype y_disp_right_limits_sq[LXSTATIONS];
33 /*scaltype x_coords[LXSTATIONS][LXLAYERS][LXMAXPOINTSONSTATION];
34 scaltype tx_vals[LXSTATIONS - 1][LXMAXPOINTSONSTATION];
35 scaltype x_errs[LXSTATIONS][LXLAYERS][LXMAXPOINTSONSTATION];
36 scaltype y_coords[LXSTATIONS][LXLAYERS][LXMAXPOINTSONSTATION];
37 scaltype ty_vals[LXSTATIONS][LXMAXPOINTSONSTATION];
38 scaltype y_errs[LXSTATIONS][LXLAYERS][LXMAXPOINTSONSTATION];
39 scaltype z_coords[LXSTATIONS][LXLAYERS][LXMAXPOINTSONSTATION];
40 LxPoint* point_refs[LXSTATIONS][LXLAYERS][LXMAXPOINTSONSTATION];
41 bool use_points[LXSTATIONS - 1][LXMAXPOINTSONSTATION];
42 int points_counts[LXSTATIONS][LXLAYERS];*/
43 
44 void LxSpace::InitGlobalCAArrays() {
45  memset(x_coords,
46  0,
47  sizeof(scaltype) * LXSTATIONS * LXLAYERS * LXMAXPOINTSONSTATION);
48  memset(tx_vals,
49  0,
50  sizeof(scaltype) * (LXSTATIONS - 1) * LXLAYERS * LXMAXPOINTSONSTATION);
51  memset(
52  x_errs, 0, sizeof(scaltype) * LXSTATIONS * LXLAYERS * LXMAXPOINTSONSTATION);
53  memset(y_coords,
54  0,
55  sizeof(scaltype) * LXSTATIONS * LXLAYERS * LXMAXPOINTSONSTATION);
56  memset(ty_vals,
57  0,
58  sizeof(scaltype) * (LXSTATIONS - 1) * LXLAYERS * LXMAXPOINTSONSTATION);
59  memset(
60  y_errs, 0, sizeof(scaltype) * LXSTATIONS * LXLAYERS * LXMAXPOINTSONSTATION);
61  memset(z_coords,
62  0,
63  sizeof(scaltype) * LXSTATIONS * LXLAYERS * LXMAXPOINTSONSTATION);
64  memset(point_refs,
65  0,
66  sizeof(LxPoint*) * LXSTATIONS * LXLAYERS * LXMAXPOINTSONSTATION);
67  memset(points_counts, 0, sizeof(int) * LXSTATIONS * LXLAYERS);
68 
69  for (int i = 0; i < stationsInAlgo - 2; ++i)
70  memset(&use_points[i][0], 0, sizeof(use_points[i]));
71 
72  for (int i = stationsInAlgo - 2; i < LXSTATIONS - 1; ++i)
73  memset(&use_points[i][0], 1, sizeof(use_points[LXSTATIONS - 2]));
74 }
75 
76 //#define LXSIMDIZE
77 
78 #ifdef LXSIMDIZE
79 
80 #include "immintrin.h"
81 
82 #define vectype __m256
83 #define veclen 8
84 #define vec_load_ps _mm256_load_ps
85 #define vec_store_ps _mm256_store_ps
86 #define vec_set1_ps _mm256_set1_ps
87 #define vec_add_ps _mm256_add_ps
88 #define vec_sub_ps _mm256_sub_ps
89 #define vec_mul_ps _mm256_mul_ps
90 #define vec_div_ps _mm256_div_ps
91 #define vec_cmp_ps _mm256_cmp_ps
92 #define vec_or_ps _mm256_or_ps
93 #define vec_sqrt_ps _mm256_sqrt_ps
94 
95 void LxSpace::RefineMiddlePoints() {
96  for (int station_number = LXSTATIONS - 1;
97  station_number >= stationsInAlgo - 1;
98  --station_number) {
99  for (int i = 0; i < points_counts[station_number][1]; ++i) {
100  scaltype xC = x_coords[station_number][1][i];
101  scaltype yC = y_coords[station_number][1][i];
102  scaltype zC = z_coords[station_number][1][i];
103  scaltype tx = xC / zC;
104  scaltype ty = yC / zC;
105  scaltype xErrCSq =
106  x_errs[station_number][1][i] * x_errs[station_number][1][i];
107  scaltype yErrCSq =
108  y_errs[station_number][1][i] * y_errs[station_number][1][i];
109  bool validC = false;
110 
111  for (int j = 0; j < points_counts[station_number][2]; ++j) {
112  scaltype zR = z_coords[station_number][2][j];
113  scaltype deltaZ = zR - zC;
114  scaltype xR = x_coords[station_number][2][j];
115  scaltype xErrRSq =
116  x_errs[station_number][2][j] * x_errs[station_number][2][j];
117  scaltype x = xC + tx * deltaZ;
118  scaltype deltaX = x - xR;
119  scaltype deltaXSq = deltaX * deltaX;
120 
121  if (deltaXSq
122  > 16 * (xErrCSq + xErrRSq) + x_disp_right_limits_sq[station_number])
123  continue;
124 
125  scaltype yR = y_coords[station_number][2][j];
126  scaltype yErrRSq =
127  y_errs[station_number][2][j] * y_errs[station_number][2][j];
128  scaltype y = yC + ty * deltaZ;
129  scaltype deltaY = y - yR;
130  scaltype deltaYSq = deltaY * deltaY;
131 
132  if (deltaYSq
133  > 16 * (yErrCSq + yErrRSq) + y_disp_right_limits_sq[station_number])
134  continue;
135 
136  validC = true;
137  break;
138  }
139 
140  if (validC) continue;
141 
142  for (int j = 0; j < points_counts[station_number][0]; ++j) {
143  scaltype zL = z_coords[station_number][0][j];
144  scaltype deltaZ = zL - zC;
145  scaltype xL = x_coords[station_number][0][j];
146  scaltype xErrLSq =
147  x_errs[station_number][0][j] * x_errs[station_number][0][j];
148  scaltype x = xC + tx * deltaZ;
149  scaltype deltaX = x - xL;
150  scaltype deltaXSq = deltaX * deltaX;
151 
152  if (deltaXSq
153  > 16 * (xErrCSq + xErrLSq) + x_disp_left_limits_sq[station_number])
154  continue;
155 
156  scaltype yL = y_coords[station_number][0][j];
157  scaltype yErrLSq =
158  y_errs[station_number][0][j] * y_errs[station_number][0][j];
159  scaltype y = yC + ty * deltaZ;
160  scaltype deltaY = y - yL;
161  scaltype deltaYSq = deltaY * deltaY;
162 
163  if (deltaYSq
164  > 16 * (yErrCSq + yErrLSq) + y_disp_left_limits_sq[station_number])
165  continue;
166 
167  validC = true;
168  break;
169  }
170 
171  if (!validC) {
172  use_points[station_number - 1][i] = 0;
173  point_refs[station_number][1][i]->valid = false;
174  }
175  }
176  }
177 }
178 
179 //#include <omp.h>
180 
181 void LxSpace::CalcTangents(int station_number) {
182  int l_station_number = station_number - 1;
183 
184  if (points_counts[station_number][1] > 200) {
185  //omp_set_num_threads(4);
186  //#pragma omp parallel for
187  for (int i = 0; i < points_counts[station_number][1]; i += veclen) {
188  if (!*reinterpret_cast<unsigned long long*>(
189  &use_points[l_station_number][i]))
190  continue;
191 
192  vectype x = vec_load_ps(&x_coords[station_number][1][i]);
193  vectype y = vec_load_ps(&y_coords[station_number][1][i]);
194  vectype z = vec_load_ps(&z_coords[station_number][1][i]);
195  vectype tx = vec_div_ps(x, z);
196  vectype ty = vec_div_ps(y, z);
197  vec_store_ps(&tx_vals[l_station_number][i], tx);
198  vec_store_ps(&ty_vals[l_station_number][i], ty);
199  }
200  } else {
201  for (int i = 0; i < points_counts[station_number][1]; i += veclen) {
202  if (!*reinterpret_cast<unsigned long long*>(
203  &use_points[l_station_number][i]))
204  continue;
205 
206  vectype x = vec_load_ps(&x_coords[station_number][1][i]);
207  vectype y = vec_load_ps(&y_coords[station_number][1][i]);
208  vectype z = vec_load_ps(&z_coords[station_number][1][i]);
209  vectype tx = vec_div_ps(x, z);
210  vectype ty = vec_div_ps(y, z);
211  vec_store_ps(&tx_vals[l_station_number][i], tx);
212  vec_store_ps(&ty_vals[l_station_number][i], ty);
213  }
214  }
215 }
216 
217 void LxSpace::BuildRaysGlobal() {
218  vectype tx_coeff_sq = vec_set1_ps(errorTxCoeffSq);
219  vectype ty_coeff_sq = vec_set1_ps(errorTyCoeffSq);
220  scaltype res_buf[veclen];
221  scaltype tx_buf[veclen];
222  scaltype ty_buf[veclen];
223  scaltype dtx_buf[veclen];
224  scaltype dty_buf[veclen];
225  int iter;
226 
227  //#pragma omp parallel for
228  for (int i = LXSTATIONS - 1; i > 0; --i) {
229  vectype tx_limit_sq = vec_set1_ps(tx_limits_sq[i]);
230  vectype ty_limit_sq = vec_set1_ps(ty_limits_sq[i]);
231  CalcTangents(i);
232  int i2 = i - 1;
233  int i3 = i2 - 1;
234 
235  for (int j = 0; j < points_counts[i][1]; ++j) {
236  if (!use_points[i2][j]) continue;
237 
238  vectype r_x = vec_set1_ps(x_coords[i][1][j]);
239  vectype r_dx = vec_set1_ps(x_errs[i][1][j]);
240  vectype r_y = vec_set1_ps(y_coords[i][1][j]);
241  vectype r_dy = vec_set1_ps(y_errs[i][1][j]);
242  vectype r_z = vec_set1_ps(z_coords[i][1][j]);
243  vectype tx0 = vec_set1_ps(tx_vals[i2][j]);
244  vectype ty0 = vec_set1_ps(ty_vals[i2][j]);
245 
246  vectype r_dx_sq = vec_mul_ps(r_dx, r_dx);
247  vectype r_dy_sq = vec_mul_ps(r_dy, r_dy);
248 
249  //#pragma omp parallel for
250  for (int k = 0; k < points_counts[i2][1]; k += veclen) {
251  vectype l_x = vec_load_ps(&x_coords[i2][1][k]);
252  vectype l_dx = vec_load_ps(&x_errs[i2][1][k]);
253  vectype l_y = vec_load_ps(&y_coords[i2][1][k]);
254  vectype l_dy = vec_load_ps(&y_errs[i2][1][k]);
255  vectype l_z = vec_load_ps(&z_coords[i2][1][k]);
256 
257  vectype delta_x = vec_sub_ps(l_x, r_x);
258  vectype delta_y = vec_sub_ps(l_y, r_y);
259  vectype delta_z = vec_sub_ps(l_z, r_z);
260  vectype delta_z_sq = vec_mul_ps(delta_z, delta_z);
261 
262  vectype tx = vec_div_ps(delta_x, delta_z);
263  vectype ty = vec_div_ps(delta_y, delta_z);
264  vectype l_dx_sq = vec_mul_ps(l_dx, l_dx);
265  vectype l_dy_sq = vec_mul_ps(l_dy, l_dy);
266 
267  vectype dtx_sq = vec_add_ps(r_dx_sq, l_dx_sq);
268  dtx_sq = vec_div_ps(dtx_sq, delta_z_sq);
269  vectype dtx = vec_sqrt_ps(dtx_sq);
270  dtx_sq = vec_mul_ps(dtx_sq, tx_coeff_sq);
271  dtx_sq = vec_add_ps(dtx_sq, tx_limit_sq);
272  vectype dty_sq = vec_add_ps(r_dy_sq, l_dy_sq);
273  dty_sq = vec_div_ps(dty_sq, delta_z_sq);
274  vectype dty = vec_sqrt_ps(dty_sq);
275  dty_sq = vec_mul_ps(dty_sq, ty_coeff_sq);
276  dty_sq = vec_add_ps(dtx_sq, ty_limit_sq);
277 
278  vectype diff_tx = vec_sub_ps(tx, tx0);
279  vectype diff_ty = vec_sub_ps(ty, ty0);
280  vectype diff_tx_sq = vec_mul_ps(diff_tx, diff_tx);
281  vectype diff_ty_sq = vec_mul_ps(diff_ty, diff_ty);
282 
283  vectype res_tx = vec_cmp_ps(diff_tx_sq, dtx_sq, _CMP_GT_OS);
284  vectype res_ty = vec_cmp_ps(diff_ty_sq, dty_sq, _CMP_GT_OS);
285  vectype res = vec_or_ps(res_tx, res_ty);
286  vec_store_ps(res_buf, res);
287  vec_store_ps(tx_buf, tx);
288  vec_store_ps(ty_buf, ty);
289  vec_store_ps(dtx_buf, dtx);
290  vec_store_ps(dty_buf, dty);
291  iter = k;
292 
293  for (int l = 0; l < veclen; ++l) {
294  if (res_buf[l]) {
295  ++iter;
296  continue;
297  }
298 
299  if (iter >= points_counts[i2][1]) break;
300 
301  point_refs[i][1][j]->CreateRay(point_refs[i2][1][iter],
302  tx_buf[l],
303  ty_buf[l],
304  dtx_buf[l],
305  dty_buf[l]);
306 
307  if (i3 > -1) use_points[i3][iter] = true;
308 
309  ++iter;
310  }
311  }
312  }
313  }
314 }
315 
316 #else //LXSIMDIZE
317 
318 #define _mm_malloc(X, Y) malloc((X))
319 #define _mm_free free
320 
321 scaltype x_coords[LXSTATIONS][LXLAYERS][LXMAXPOINTSONSTATION];
322 scaltype tx_vals[LXSTATIONS - 1][LXMAXPOINTSONSTATION];
323 scaltype x_errs[LXSTATIONS][LXLAYERS][LXMAXPOINTSONSTATION];
324 scaltype y_coords[LXSTATIONS][LXLAYERS][LXMAXPOINTSONSTATION];
325 scaltype ty_vals[LXSTATIONS][LXMAXPOINTSONSTATION];
326 scaltype y_errs[LXSTATIONS][LXLAYERS][LXMAXPOINTSONSTATION];
327 scaltype z_coords[LXSTATIONS][LXLAYERS][LXMAXPOINTSONSTATION];
328 LxPoint* point_refs[LXSTATIONS][LXLAYERS][LXMAXPOINTSONSTATION];
329 bool use_points[LXSTATIONS - 1][LXMAXPOINTSONSTATION];
330 int points_counts[LXSTATIONS][LXLAYERS];
331 
332 static void CalcTangents(int station_number) {
333  int l_station_number = station_number - 1;
334 
335  for (int i = 0; i < points_counts[station_number][1]; ++i) {
336  if (!use_points[l_station_number][i]) continue;
337 
338  tx_vals[l_station_number][i] =
339  x_coords[station_number][1][i] / z_coords[station_number][1][i];
340  ty_vals[l_station_number][i] =
341  y_coords[station_number][1][i] / z_coords[station_number][1][i];
342  }
343 }
344 
345 void BuildRaysGlobal() {
346  for (int i = LXSTATIONS - 1; i > 0; --i) {
347  scaltype tx_limit_sq = tx_limits[i] * tx_limits[i];
348  scaltype ty_limit_sq = ty_limits[i] * ty_limits[i];
349  CalcTangents(i);
350  int i2 = i - 1;
351  int i3 = i2 - 1;
352 
353  for (int j = 0; j < points_counts[i][1]; ++j) {
354  if (!use_points[i2][j]) continue;
355 
356  scaltype r_x = x_coords[i][1][j];
357  scaltype r_y = y_coords[i][1][j];
358  scaltype r_z = z_coords[i][1][j];
359  scaltype r_dx_sq = x_errs[i][1][j] * x_errs[i][1][j];
360  scaltype r_dy_sq = y_errs[i][1][j] * y_errs[i][1][j];
361 
362  for (int k = 0; k < points_counts[i2][1]; ++k) {
363  scaltype delta_x = x_coords[i2][1][k] - r_x;
364  scaltype delta_y = y_coords[i2][1][k] - r_y;
365  scaltype delta_z = z_coords[i2][1][k] - r_z;
366  scaltype delta_z_sq = delta_z * delta_z;
367 
368  scaltype tx = delta_x / delta_z;
369  scaltype ty = delta_y / delta_z;
370  scaltype l_dx_sq = x_errs[i2][1][k] * x_errs[i2][1][k];
371  scaltype l_dy_sq = y_errs[i2][1][k] * y_errs[i2][1][k];
372 
373  scaltype dtx_sq = (r_dx_sq + l_dx_sq) / delta_z_sq;
374  scaltype dtx = sqrt(dtx_sq);
375  dtx_sq *= errorTxCoeffSq;
376  dtx_sq += tx_limit_sq;
377  scaltype dty_sq = (r_dy_sq + l_dy_sq) / delta_z_sq;
378  scaltype dty = sqrt(dty_sq);
379  dty_sq *= errorTyCoeffSq;
380  dty_sq += ty_limit_sq;
381 
382  scaltype diff_tx = tx - tx_vals[i2][j];
383  scaltype diff_ty = ty - ty_vals[i2][j];
384  scaltype diff_tx_sq = diff_tx * diff_tx;
385  scaltype diff_ty_sq = diff_ty * diff_ty;
386 
387  if ((diff_tx_sq > dtx_sq) || (diff_ty_sq > dty_sq)) continue;
388 
389  point_refs[i][1][j]->CreateRay(point_refs[i2][1][k], tx, ty, dtx, dty);
390 
391  if (i3 > -1) use_points[i3][k] = true;
392  }
393  }
394  }
395 }
396 
397 #endif //LXSIMDIZE
398 
400 // Global declarations
402 
403 #ifdef CLUSTER_MODE
404 static scaltype squareRoot2 = sqrt(2.0);
405 static scaltype errorXcoeff = 4.0;
406 static scaltype errorYcoeff = 4.0;
407 static scaltype errorTxCoeff = 4.0;
408 static scaltype errorTyCoeff = 4.0;
409 #else //CLUSTER_MODE
410 static scaltype errorXcoeff = 4.0;
411 static scaltype errorYcoeff = 4.0;
412 #endif //CLUSTER_MODE
413 
414 static scaltype errorTxBreakCoeff = 4.0;
415 static scaltype errorTyBreakCoeff = 4.0;
416 
417 // These ..Ext.. coefficients are applied to squares of uncertainties. So they can be bigger than the previous.
418 //static scaltype errorExtXcoeff = 16;
419 //static scaltype errorExtYcoeff = 16;
420 //static scaltype errorExtTxCoeff = 16;
421 //static scaltype errorExtTyCoeff = 16;
422 
423 //static LxFinder* finderInstance = 0;
424 
426 // LxTrackCandidate
428 
429 struct LxTrackCandidate {
430  LxRay* rays[LXSTATIONS - 1]; // Rays are stored left to right.
431  Int_t length;
432  scaltype chi2;
433 
434  LxTrackCandidate(LxRay** rr, Int_t len, scaltype c2) : length(len), chi2(c2) {
435  for (int i = 0; i < len; ++i)
436  rays[i] = rr[i];
437 
438  for (int i = len; i < LXSTATIONS - 1; ++i)
439  rays[i] = 0;
440  }
441 };
442 
444 // LxTrack
446 
447 #ifdef CLUSTER_MODE
448 LxTrack::LxTrack(LxTrackCandidate* tc)
449  : matched(false)
450  , length(tc->length)
451  , chi2(tc->chi2)
452  , mcTrack(0)
453  , aX(0)
454  , bX(0)
455  , aY(0)
456  , bY(0)
457  , restoredPoints(0)
458  , externalTrack(0)
459 #ifdef USE_KALMAN_FIT
460  , x(0)
461  , y(0)
462  , z(0)
463  , dx(0)
464  , dy(0)
465  , tx(0)
466  , ty(0)
467  , dtx(0)
468  , dty(0)
469 #endif //USE_KALMAN_FIT
470 {
471  memset(rays, 0, sizeof(rays));
472  memset(points, 0, sizeof(points));
473 
474  for (Int_t i = 0; i < length; ++i) {
475  LxRay* ray = tc->rays[i];
476  rays[i] = ray;
477  LxPoint* point = ray->source;
478  point->used = true;
479  point->track = this;
480 
481  for (list<LxPoint*>::iterator j = ray->clusterPoints.begin();
482  j != ray->clusterPoints.end();
483  ++j) {
484  LxPoint* p = *j;
485  p->track = this;
486  }
487  }
488 
489  LxRay* ray = tc->rays[0];
490  LxPoint* point = ray->end;
491  point->used = true;
492  point->track = this;
493 }
494 #else //CLUSTER_MODE
495 LxTrack::LxTrack(LxTrackCandidate* tc)
496  : externalTrack(NULL)
497 #ifdef LX_EXT_LINK_SOPH
498  , extTrackCandidates()
499 #else //LX_EXT_LINK_SOPH
500  , extLinkChi2()
501 #endif
502  , matched(false)
503  , mcTrack(NULL)
504 #ifdef CALC_LINK_WITH_STS_EFF
505  , mcTracks()
506 #endif //CALC_LINK_WITH_STS_EFF
507  , length(tc->length)
508  , rays()
509  , points()
510  , chi2(tc->chi2)
511  , aX(0)
512  , bX(0)
513  , aY(0)
514  , bY(0)
515  , restoredPoints(0)
516 #ifdef USE_KALMAN_FIT
517  , x(0)
518  , y(0)
519  , z(0)
520  , dx(0)
521  , dy(0)
522  , tx(0)
523  , ty(0)
524  , dtx(0)
525  , dty(0)
526 #endif //USE_KALMAN_FIT
527  , clone(false)
528  , distanceOk(false)
529  , oppCharged(false)
530  , triggering(false) {
531  memset(rays, 0, sizeof(rays));
532  memset(points, 0, sizeof(points));
533 
534  for (int i = 0; i < length; ++i) {
535  LxRay* ray = tc->rays[i];
536  rays[i] = ray;
537  LxPoint* point = ray->source;
538  point->used = true;
539  point->track = this;
540  points[(i + 1) * LXLAYERS + LXMIDDLE] = point;
541 
542  if (point->artificial) ++restoredPoints;
543 
544  LxStation* station = ray->station;
545 
546  for (int j = 0; j < LXLAYERS; ++j) {
547  if (LXMIDDLE == j) continue;
548 
549  LxLayer* layer = station->layers[j];
550  LxPoint* p = layer->PickNearestPoint(ray);
551 
552  if (0 != p && !p->used) {
553  p->used = true;
554  p->track = this;
555  points[(i + 1) * LXLAYERS + j] = p;
556  }
557  }
558  }
559 
560  LxRay* ray = tc->rays[0];
561  LxPoint* point = ray->end;
562  point->used = true;
563  point->track = this;
564  points[LXMIDDLE] = point;
565 
566  if (point->artificial) ++restoredPoints;
567 
568  LxStation* station = ray->station->space->stations[0];
569 
570  for (int i = 0; i < LXLAYERS; ++i) {
571  if (LXMIDDLE == i) continue;
572 
573  LxLayer* layer = station->layers[i];
574  scaltype diffZ = layer->zCoord - ray->end->z;
575  scaltype x = ray->end->x + ray->tx * diffZ;
576  scaltype y = ray->end->y + ray->ty * diffZ;
577  LxPoint* p = layer->PickNearestPoint(x, y);
578 
579  if (0 != p && !p->used) {
580  p->used = true;
581  p->track = this;
582  points[i] = p;
583  }
584  }
585 
586  //cout << "This track contains: " << restoredPoints << " restored points" << endl;
587 }
588 #endif //CLUSTER_MODE
589 
590 static inline void Ask() {
591  char symbol;
592 
593  cout << "ask>";
594 
595  do {
596  cin.get(symbol);
597  } while (symbol != '\n');
598 
599  cout << endl;
600 }
601 
602 #ifdef USE_KALMAN_FIT
603 struct KFParams {
604  scaltype coord, tg, C11, C12, C21, C22;
605 };
606 
607 void LxTrack::Fit() {
608  LxRay* firstRay = rays[LXSTATIONS - LXFIRSTSTATION - 2];
609  LxPoint* firstPoint = firstRay->source;
610  KFParams params[2] = {
611  {firstPoint->x, firstRay->tx, firstPoint->dx, 0, 0, 1.0},
612  {firstPoint->y, firstRay->ty, firstPoint->dy, 0, 0, 1.0}};
613 
614  LxPoint* prevPoint = firstPoint;
615  scaltype aChi2 = 0;
616 
617  for (Int_t i = LXSTATIONS - LXFIRSTSTATION - 2; i >= 0; --i) {
618  LxRay* ray = rays[i];
619  LxPoint* point = ray->end;
620  LxStation* station = point->layer->station;
621  scaltype m[2] = {point->x, point->y};
622  scaltype V[2] = {point->dx * point->dx, point->dy * point->dy};
623  KFParams pPrev[2] = {params[0], params[1]};
624  scaltype deltaZ = point->z - prevPoint->z;
625  scaltype deltaZ2 = deltaZ * deltaZ;
626 
627  for (Int_t k = 0; k <= 1; ++k) {
628  // Extrapolate.
629  params[k].coord += pPrev[k].tg * deltaZ; // params[k].tg is unchanged.
630 
631  // Filter.
632  params[k].C11 += pPrev[k].C12 * deltaZ + pPrev[k].C21 * deltaZ
633  + pPrev[k].C22 * deltaZ2 + station->MSNoise[k][0][0];
634  params[k].C12 += pPrev[k].C22 * deltaZ + station->MSNoise[k][0][1];
635  params[k].C21 += pPrev[k].C22 * deltaZ + station->MSNoise[k][1][0];
636  params[k].C22 += station->MSNoise[k][1][1];
637 
638  scaltype S = 1.0 / (V[k] + params[k].C11);
639  scaltype Kcoord = params[k].C11 * S;
640  scaltype Ktg = params[k].C21 * S;
641  scaltype dzeta = m[k] - params[k].coord;
642  params[k].coord += Kcoord * dzeta;
643  params[k].tg += Ktg * dzeta;
644  params[k].C21 -= params[k].C11 * Ktg;
645  params[k].C22 -= params[k].C12 * Ktg;
646  params[k].C11 *= 1.0 - Kcoord;
647  params[k].C12 *= 1.0 - Kcoord;
648  aChi2 += dzeta * S * dzeta;
649  }
650 
651  prevPoint = point;
652  }
653 
654  x = params[0].coord;
655  dx = sqrt(params[0].C11);
656  tx = params[0].tg;
657  dtx = sqrt(params[0].C22);
658  y = params[1].coord;
659  dy = sqrt(params[1].C11);
660  ty = params[1].tg;
661  dty = sqrt(params[1].C22);
662  z = rays[0]->end->z;
663 }
664 #else //USE_KALMAN_FIT
665 void LxTrack::Fit() {
666  scaltype sumZ = 0;
667  scaltype sumZ2 = 0;
668  scaltype sumX = 0;
669  scaltype sumZX = 0;
670  scaltype sumY = 0;
671  scaltype sumZY = 0;
672 
673  for (int i = 0; i < LXSTATIONS - LXFIRSTSTATION; ++i) {
674  LxPoint* point = points[LXLAYERS * i + LXMIDDLE];
675  scaltype z = point->z;
676  scaltype x = point->x;
677  scaltype y = point->y;
678 
679  sumZ += z;
680  sumZ2 += z * z;
681  sumX += x;
682  sumZX += z * x;
683  sumY += y;
684  sumZY += z * y;
685  }
686 
687  aX = (LXSTATIONS * sumZX - sumZ * sumX) / (LXSTATIONS * sumZ2 - sumZ * sumZ);
688  bX = (sumX - aX * sumZ) / LXSTATIONS;
689  aY = (LXSTATIONS * sumZY - sumZ * sumY) / (LXSTATIONS * sumZ2 - sumZ * sumZ);
690  bY = (sumY - aY * sumZ) / LXSTATIONS;
691 
692  //LxPoint* point = points[LXLAYERS * (LXSTATIONS - 1) + LXMIDDLE];
693  //cout << "bX = " << bX << " , aX = " << aX << " , X / Z = " << point->x / point->z << endl;
694  //cout << "bY = " << bY << " , aY = " << aY << " , Y / Z = " << point->y / point->z << endl;
695  //Ask();
696 }
697 #endif //USE_KALMAN_FIT
698 
699 #ifdef LX_EXT_LINK_SOPH
700 void LxTrack::Rebind() {
701  externalTrack = 0;
702 
703  for (list<pair<LxExtTrack*, scaltype>>::iterator i =
704  extTrackCandidates.begin();
705  i != extTrackCandidates.end();
706  ++i) {
707  LxExtTrack* extTrack = i->first;
708  scaltype aChi2 = i->second;
709 
710  if (0 == extTrack->recoTrack.first) {
711  extTrack->recoTrack.first = this;
712  extTrack->recoTrack.second = aChi2;
713  externalTrack = extTrack;
714  break;
715  } else if (aChi2 < extTrack->recoTrack.second) {
716  LxTrack* anotherTrack = extTrack->recoTrack.first;
717  extTrack->recoTrack.first = this;
718  extTrack->recoTrack.second = aChi2;
719  externalTrack = extTrack;
720  anotherTrack->Rebind();
721  break;
722  }
723  }
724 }
725 #endif // LX_EXT_LINK_SOPH
726 
728 // LxPoint
730 
731 LxPoint::~LxPoint() {
732  for (list<LxRay*>::iterator i = rays.begin(); i != rays.end(); ++i)
733  delete *i;
734 }
735 
736 void LxPoint::CreateRay(LxPoint* lPoint,
737  scaltype tx,
738  scaltype ty,
739  scaltype dtx,
740  scaltype dty) {
741  rays.push_back(new LxRay(this,
742  lPoint,
743  tx,
744  ty,
745  dtx,
746  dty
747 #ifdef CLUSTER_MODE
748  ,
749  4
750 #endif //CLUSTER_MODE
751  ));
752 }
753 
755 // LxRay
757 
758 LxRay::LxRay(LxPoint* s,
759  LxPoint* e
760 #ifdef CLUSTER_MODE
761  ,
762  Int_t l
763 #endif //CLUSTER_MODE
764  )
765  : tx((e->x - s->x) / (e->z - s->z))
766  , ty((e->y - s->y) / (e->z - s->z))
767  , dtx(sqrt(e->dx * e->dx + s->dx * s->dx) / (s->z - e->z))
768  , dty(sqrt(e->dy * e->dy + s->dy * s->dy) / (s->z - e->z))
769  , source(s)
770  , end(e)
771  , station(s->layer->station)
772  , neighbours()
773 #ifdef CLUSTER_MODE
774  , level(l)
775  , used(false)
776  , neighbourhood()
777  , clusterPoints()
778 #endif //CLUSTER_MODE
779 #ifdef USE_KALMAN
780  , kalman()
781 #endif
782 {
783 }
784 
785 LxRay::LxRay(LxPoint* s,
786  LxPoint* e,
787  scaltype Tx,
788  scaltype Ty,
789  scaltype Dtx,
790  scaltype Dty
791 #ifdef CLUSTER_MODE
792  ,
793  Int_t l
794 #endif //CLUSTER_MODE
795  )
796  : tx(Tx)
797  , ty(Ty)
798  , dtx(Dtx)
799  , dty(Dty)
800  , source(s)
801  , end(e)
802  , station(s->layer->station)
803  , neighbours()
804 #ifdef CLUSTER_MODE
805  , level(l)
806  , used(false)
807  , neighbourhood()
808  , clusterPoints()
809 #endif //CLUSTER_MODE
810 #ifdef USE_KALMAN
811  , kalman()
812 #endif
813 
814 {
815 }
816 
818 // LxLayer
820 
821 LxLayer::LxLayer(LxStation* st, int lNum)
822  : station(st), layerNumber(lNum), zCoord(0) {}
823 
824 LxLayer::~LxLayer() { Clear(); }
825 
826 void LxLayer::Clear() {
827  for (list<LxPoint*>::iterator i = points.begin(); i != points.end(); ++i)
828  delete *i;
829 
830  points.clear();
831 }
832 
833 LxPoint* LxLayer::PickNearestPoint(scaltype x, scaltype y) {
834  LxPoint* result = 0;
835  scaltype minR2 = 0;
836 
837  for (list<LxPoint*>::iterator i = points.begin(); i != points.end(); ++i) {
838  LxPoint* point = *i;
839 
840  if (point->used) continue;
841 
842  scaltype diffX = point->x - x;
843  scaltype diffY = point->y - y;
844 
845  if (diffX < 0) diffX = -diffX;
846 
847  if (diffX > errorXcoeff * point->dx) continue;
848 
849  if (diffY < 0) diffY = -diffY;
850 
851  if (diffY > errorYcoeff * point->dy) continue;
852 
853  scaltype r2 = diffX * diffX + diffY * diffY;
854 
855  if (0 == result || r2 < minR2) {
856  result = point;
857  minR2 = r2;
858  }
859  }
860 
861  return result;
862 }
863 
864 LxPoint* LxLayer::PickNearestPoint(LxRay* ray) {
865  LxPoint* point = ray->source;
866  scaltype diffZ = zCoord - point->z;
867  scaltype x = point->x + ray->tx * diffZ;
868  scaltype y = point->y + ray->ty * diffZ;
869  return PickNearestPoint(x, y);
870 }
871 
872 LxPoint* LxLayer::PickNearestPoint(scaltype x,
873  scaltype y,
874  scaltype deltaX,
875  scaltype deltaY) {
876  LxPoint* result = 0;
877  scaltype minR2 = 0;
878  scaltype xLBound = x - deltaX;
879  scaltype xUBound = x + deltaX;
880  scaltype yLBound = y - deltaY;
881  scaltype yUBound = y + deltaY;
882 
883  for (list<LxPoint*>::iterator i = points.begin(); i != points.end(); ++i) {
884  LxPoint* point = *i;
885 
886  if (point->x < xLBound || point->x > xUBound || point->y < yLBound
887  || point->y > yUBound)
888  continue;
889 
890  scaltype diffX = point->x - x;
891  scaltype diffY = point->y - y;
892  scaltype r2 = diffX * diffX + diffY * diffY;
893 
894  if (0 == result || r2 < minR2) {
895  result = point;
896  minR2 = r2;
897  }
898  }
899 
900  return result;
901 }
902 
903 bool LxLayer::HasPointInRange(scaltype x,
904  scaltype y,
905  scaltype deltaX,
906  scaltype deltaY) {
907  scaltype xLBound = x - deltaX;
908  scaltype xUBound = x + deltaX;
909  scaltype yLBound = y - deltaY;
910  scaltype yUBound = y + deltaY;
911 
912  for (list<LxPoint*>::iterator i = points.begin(); i != points.end(); ++i) {
913  LxPoint* point = *i;
914 
915  if (xLBound < point->x && point->x < xUBound && yLBound < point->y
916  && point->y < yUBound)
917  return true;
918  }
919 
920  return false;
921 }
922 
924 // LxStation
926 
927 // These arrays are filled with data obtained during Monte Carlo simulations.
928 //static scaltype txLimits[] = { 0., 0.13, 0.11, 0.11, 0.125, 0.17 };
929 //static scaltype tyLimits[] = { 0., 0.125, 0.07, 0.07, 0.1, 0.21 };
930 
931 //static scaltype txBreakLimits[] = { 0., 0.085, 0.09, 0.13, 0.21, 0. };
932 //static scaltype tyBreakLimits[] = { 0., 0.07, 0.085, 0.11, 0.23, 0. };
933 
934 //static scaltype txBreakSigmas[] = { 0., 0.01242, 0.0115, 0.0134, 0.02955, 0. };
935 //static scaltype tyBreakSigmas[] = { 0., 0.01473, 0.01174, 0.0136, 0.02827, 0. };
936 
937 static scaltype dispersions01XSmall[] = {0.6, 1.0, 0.8, 0.8, 1.2, 1.2};
938 static scaltype dispersions01XBig[] = {1.9, 1.8, 1.6, 1.6, 1.9, 2.9};
939 static scaltype dispersions01YSmall[] = {0.4, 0.5, 0.7, 0.7, 1.0, 1.0};
940 static scaltype dispersions01YBig[] = {1.2, 1.3, 1.3, 1.3, 1.6, 2.9};
941 
942 static scaltype dispersions02XSmall[] = {0.7, 1.2, 1.2, 1.2, 1.1, 1.4};
943 static scaltype dispersions02XBig[] = {2.0, 1.9, 1.6, 1.7, 1.8, 2.3};
944 static scaltype dispersions02YSmall[] = {0.2, 0.6, 1.0, 0.9, 0.9, 1.4};
945 static scaltype dispersions02YBig[] = {1.5, 1.3, 1.3, 1.4, 1.7, 2.2};
946 
947 #ifdef CLUSTER_MODE
948 struct KDRayWrap {
949  LxRay* ray;
950  scaltype data[4];
951  static bool destroyRays;
952 
953  KDRayWrap(scaltype x, scaltype y, LxRay* r) : ray(r) {
954  data[0] = x;
955  data[1] = y;
956  data[2] = ray->tx;
957  data[3] = ray->ty;
958  }
959 
960  KDRayWrap(scaltype x, scaltype y, scaltype tx, scaltype ty)
961  : ray(0) // This constructor is used when setting search-range bounds.
962  {
963  data[0] = x;
964  data[1] = y;
965  data[2] = tx;
966  data[3] = ty;
967  }
968 
969  ~KDRayWrap() {
970  //if (destroyRays)
971  //delete ray;
972  }
973 
974  // Stuff required by libkdtree++
975  typedef scaltype value_type;
976 
977  value_type operator[](size_t n) const { return data[n]; }
978 };
979 
980 bool KDRayWrap::destroyRays = false;
981 
982 typedef KDTree::KDTree<4, KDRayWrap> KDRaysStorageType;
983 #endif //CLUSTER_MODE
984 
985 LxStation::LxStation(LxSpace* sp, int stNum)
986  : space(sp)
987  , stationNumber(stNum)
988  , zCoord(0)
989  , txLimit(0)
990  , tyLimit(0)
991  , txBreakLimit(0)
992  , tyBreakLimit(0)
993  , txBreakSigma(0)
994  , tyBreakSigma(0)
995  , disp01XSmall(dispersions01XSmall[stNum])
996  , disp01XBig(dispersions01XBig[stNum])
997  , disp01YSmall(dispersions01YSmall[stNum])
998  , disp01YBig(dispersions01YBig[stNum])
999  , disp02XSmall(dispersions02XSmall[stNum])
1000  , disp02XBig(dispersions02XBig[stNum])
1001  , disp02YSmall(dispersions02YSmall[stNum])
1002  , disp02YBig(dispersions02YBig[stNum])
1003 #ifdef CLUSTER_MODE
1004  , raysHandle(0)
1005  , clusteredRaysHandle(0)
1006  , clusterXLimit(0)
1007  , clusterXLimit2(0)
1008  , clusterYLimit(0)
1009  , clusterYLimit2(0)
1010  , clusterTxLimit(0)
1011  , clusterTxLimit2(0)
1012  , clusterTyLimit(0)
1013  , clusterTyLimit2(0)
1014 #endif //CLUSTER_MODE
1015 {
1016  for (int i = 0; i < LXLAYERS; ++i)
1017  layers.push_back(new LxLayer(this, i));
1018 
1019 #ifdef CLUSTER_MODE
1020  raysHandle = new KDRaysStorageType;
1021  clusteredRaysHandle = new KDRaysStorageType;
1022 #endif //CLUSTER_MODE
1023 }
1024 
1025 LxStation::~LxStation() {
1026  Clear();
1027 #ifdef CLUSTER_MODE
1028  KDRaysStorageType* rays = static_cast<KDRaysStorageType*>(raysHandle);
1029  delete rays;
1030  KDRaysStorageType* clusteredRays =
1031  static_cast<KDRaysStorageType*>(clusteredRaysHandle);
1032  delete clusteredRays;
1033 #endif //CLUSTER_MODE
1034 }
1035 
1036 void LxStation::Clear() {
1037  for (vector<LxLayer*>::iterator i = layers.begin(); i != layers.end(); ++i)
1038  (*i)->Clear();
1039 
1040 #ifdef CLUSTER_MODE
1041  KDRaysStorageType* rays = static_cast<KDRaysStorageType*>(raysHandle);
1042  KDRayWrap::destroyRays = true;
1043  rays->clear();
1044  KDRaysStorageType* clusteredRays =
1045  static_cast<KDRaysStorageType*>(clusteredRaysHandle);
1046  clusteredRays->clear();
1047  KDRayWrap::destroyRays = false;
1048 
1049  for (list<LxPoint*>::iterator i = clusteredPoints.begin();
1050  i != clusteredPoints.end();
1051  ++i)
1052  delete *i;
1053 
1054  clusteredPoints.clear();
1055 
1056  for (Int_t i = 0; i < 2 * LXLAYERS; ++i) {
1057  for (vector<list<LxRay*>*>::iterator j = clusters[i].begin();
1058  j != clusters[i].end();
1059  ++j)
1060  delete *j;
1061 
1062  clusters[i].clear();
1063  }
1064 #endif //CLUSTER_MODE
1065 }
1066 
1067 void LxStation::RestoreMiddlePoints() {
1068  LxLayer* lLayer = layers[0];
1069  LxLayer* mLayer = layers[LXMIDDLE];
1070  LxLayer* rLayer = layers[2];
1071  list<LxPoint*>& mPoints = mLayer->points;
1072  list<LxPoint*>& rPoints = rLayer->points;
1073 
1074  // 1. Loop through the points of the middle layer and check if they have at least on corresponding point on either left or right layer.
1075  for (list<LxPoint*>::iterator i = mPoints.begin(); i != mPoints.end(); ++i) {
1076  LxPoint* point = *i;
1077  scaltype tx = point->x / point->z;
1078  scaltype ty = point->y / point->z;
1079  scaltype diffZ = lLayer->zCoord - point->z;
1080  scaltype x = point->x + tx * diffZ;
1081  scaltype y = point->y + ty * diffZ;
1082 
1083  if (!lLayer->HasPointInRange(x, y, disp01XBig, disp01YBig)) {
1084  diffZ = rLayer->zCoord - point->z;
1085  x = point->x + tx * diffZ;
1086  y = point->y + tx * diffZ;
1087 
1088  if (!rLayer->HasPointInRange(x, y, disp01XBig, disp01YBig))
1089  point->valid = false;
1090  }
1091  }
1092 
1093  // 2. Loop through the points of the right station. If there are corresponding points in the middle station -- it is OK.
1094  // If there are no -- check corresponding points on the left station. If such points exist, choose the point nearest
1095  // to the predicted value among them and reconstruct a middle point as a point of intersection of the segment between
1096  // the left and corresponding right points and the middle layer.
1097  for (list<LxPoint*>::iterator i = rPoints.begin(); i != rPoints.end(); ++i) {
1098  LxPoint* rPoint = *i;
1099  scaltype tx = rPoint->x / rPoint->z;
1100  scaltype ty = rPoint->y / rPoint->z;
1101  scaltype diffZ = mLayer->zCoord - rPoint->z;
1102  scaltype x = rPoint->x + tx * diffZ;
1103  scaltype y = rPoint->y + ty * diffZ;
1104 
1105  if (mLayer->HasPointInRange(x, y, disp01XBig, disp01YBig)) continue;
1106 
1107  diffZ = lLayer->zCoord - rPoint->z;
1108  x = rPoint->x + tx * diffZ;
1109  y = rPoint->y + ty * diffZ;
1110 
1111  LxPoint* lPoint =
1112  lLayer->PickNearestPoint(x, y, disp02XSmall, disp02YSmall);
1113 
1114  if (0 == lPoint) continue;
1115 
1116  x = (lPoint->x + rPoint->x) / 2;
1117  y = (lPoint->y + rPoint->y) / 2;
1118  mLayer->AddPoint(
1119  -1, x, y, mLayer->zCoord, rPoint->dx, rPoint->dy, rPoint->dz, true);
1120  }
1121 }
1122 
1123 void LxStation::BuildRays() {
1124  if (stationNumber < LXFIRSTSTATION + 1) return;
1125 
1126  LxLayer* rLayer = layers[LXMIDDLE];
1127  LxStation* lStation =
1128  space
1129  ->stations[stationNumber
1130  - 1]; // Pointer to 'left' station. 'This' station is 'right'.
1131  LxLayer* lLayer = lStation->layers[LXMIDDLE];
1132 
1133  for (list<LxPoint*>::iterator rIter = rLayer->points.begin();
1134  rIter != rLayer->points.end();
1135  ++rIter) {
1136  LxPoint* rPoint = *rIter;
1137 
1138  if (!rPoint->valid) continue;
1139 
1140  scaltype tx1 = rPoint->x / rPoint->z;
1141  scaltype ty1 = rPoint->y / rPoint->z;
1142 
1143  for (list<LxPoint*>::iterator lIter = lLayer->points.begin();
1144  lIter != lLayer->points.end();
1145  ++lIter) {
1146  LxPoint* lPoint = *lIter;
1147 
1148  if (!lPoint->valid) continue;
1149 
1150  scaltype diffZ = lPoint->z - rPoint->z;
1151  scaltype tx = (lPoint->x - rPoint->x) / diffZ;
1152  scaltype ty = (lPoint->y - rPoint->y) / diffZ;
1153  scaltype dtx =
1154  -sqrt(lPoint->dx * lPoint->dx + rPoint->dx * rPoint->dx) / diffZ;
1155  scaltype dty =
1156  -sqrt(lPoint->dy * lPoint->dy + rPoint->dy * rPoint->dy) / diffZ;
1157  scaltype diffTx = tx - tx1;
1158  scaltype diffTy = ty - ty1;
1159 
1160  if (diffTx < 0) diffTx = -diffTx;
1161 
1162  if (diffTy < 0) diffTy = -diffTy;
1163 
1164  if (diffTx > txLimit + dtx * errorTxCoeff
1165  || diffTy > tyLimit + dty * errorTyCoeff)
1166  continue;
1167 
1168  rPoint->CreateRay(lPoint, tx, ty, dtx, dty);
1169  }
1170  }
1171 }
1172 
1173 #ifdef CLUSTER_MODE
1174 static void AddRayData(LxRay* ray,
1175  scaltype& lX,
1176  scaltype& lY,
1177  scaltype& lDx2,
1178  scaltype& lDy2,
1179  scaltype& rX,
1180  scaltype& rY,
1181  scaltype& rDx2,
1182  scaltype& rDy2) {
1183  LxPoint* lPoint = ray->end;
1184  scaltype deltaZ = lPoint->layer->station->zCoord - lPoint->z;
1185  lX += lPoint->x + ray->tx * deltaZ;
1186  lY += lPoint->y + ray->ty * deltaZ;
1187  lDx2 += lPoint->dx * lPoint->dx;
1188  lDy2 += lPoint->dy * lPoint->dy;
1189 
1190  LxPoint* rPoint = ray->source;
1191  deltaZ = rPoint->layer->station->zCoord - rPoint->z;
1192  rX += rPoint->x + ray->tx * deltaZ;
1193  rY += rPoint->y + ray->ty * deltaZ;
1194  rDx2 += rPoint->dx * rPoint->dx;
1195  rDy2 += rPoint->dy * rPoint->dy;
1196 }
1197 
1198 void LxStation::InsertClusterRay(Int_t levels,
1199  Int_t cardinality,
1200  LxRay* clusterRay) {
1201  if (levels < (LXLAYERS - 1) * (LXLAYERS - 1)) return;
1202 
1203  levels -= (LXLAYERS - 1) * (LXLAYERS - 1);
1204 
1205  while (clusters[levels].size() < cardinality)
1206  clusters[levels].push_back(new list<LxRay*>);
1207 
1208  clusters[levels][cardinality - 1]->push_back(clusterRay);
1209 }
1210 
1211 #ifdef BEST_SIX_POINTS
1212 struct LxLess {
1213  bool operator()(pair<LxPoint*, LxPoint*> lVal,
1214  pair<LxPoint*, LxPoint*> rVal) const {
1215  if (lVal.first < rVal.first) return true;
1216 
1217  return lVal.second < rVal.second;
1218  }
1219 };
1220 
1221 struct LxRaysCandidates {
1222  LxRay* data[LXLAYERS * LXLAYERS];
1223  scaltype chi2;
1224 
1225  LxRaysCandidates() : chi2(-1.0) { memset(data, 0, sizeof(data)); }
1226 };
1227 #endif //BEST_SIX_POINTS
1228 
1229 void LxStation::BuildRays2() {
1230  if (stationNumber < LXFIRSTSTATION + 1) return;
1231 
1232  KDRaysStorageType* rays = static_cast<KDRaysStorageType*>(raysHandle);
1233  LxStation* lStation =
1234  space
1235  ->stations[stationNumber
1236  - 1]; // Pointer to 'left' station. 'This' station is 'right'.
1237 
1238  for (Int_t i = 0; i < LXLAYERS; ++i) {
1239  LxLayer* rLayer = layers[i];
1240 
1241  for (Int_t j = 0; j < LXLAYERS; ++j) {
1242  LxLayer* lLayer = lStation->layers[j];
1243 
1244  for (list<LxPoint*>::iterator rIter = rLayer->points.begin();
1245  rIter != rLayer->points.end();
1246  ++rIter) {
1247  LxPoint* rPoint = *rIter;
1248  scaltype tx1 = rPoint->x / rPoint->z;
1249  scaltype ty1 = rPoint->y / rPoint->z;
1250  scaltype zShift = zCoord - rPoint->z;
1251 
1252  for (list<LxPoint*>::iterator lIter = lLayer->points.begin();
1253  lIter != lLayer->points.end();
1254  ++lIter) {
1255  LxPoint* lPoint = *lIter;
1256  scaltype diffZ = lPoint->z - rPoint->z;
1257  scaltype tx = (lPoint->x - rPoint->x) / diffZ;
1258  scaltype ty = (lPoint->y - rPoint->y) / diffZ;
1259  scaltype dtx =
1260  -sqrt(lPoint->dx * lPoint->dx + rPoint->dx * rPoint->dx) / diffZ;
1261  scaltype dty =
1262  -sqrt(lPoint->dy * lPoint->dy + rPoint->dy * rPoint->dy) / diffZ;
1263  scaltype diffTx = abs(tx - tx1);
1264  scaltype diffTy = abs(ty - ty1);
1265 
1266  if (diffTx > txLimit + dtx * errorTxCoeff
1267  || diffTy > tyLimit + dty * errorTyCoeff)
1268  continue;
1269 
1270  scaltype x = rPoint->x + tx * zShift;
1271  scaltype y = rPoint->y + ty * zShift;
1272  LxRay* ray = new LxRay(rPoint, lPoint, LXLAYERS * i + j);
1273  KDRayWrap rayWrap(x, y, ray);
1274  rays->insert(rayWrap);
1275  rPoint->rays.push_back(ray);
1276  }
1277  }
1278  }
1279  }
1280 
1281  timeval bTime, eTime;
1282  Int_t exeDuration;
1283  static Int_t maxExeDuration = 0;
1284  static Int_t goodClusters = 0;
1285  gettimeofday(&bTime, 0);
1286 
1287  for (KDRaysStorageType::iterator i = rays->begin(); i != rays->end(); ++i) {
1288  KDRayWrap& wrap = const_cast<KDRayWrap&>(*i);
1289  LxRay* ray = wrap.ray;
1290  KDTree::_Region<
1291  4,
1292  KDRayWrap,
1293  scaltype,
1294  KDTree::_Bracket_accessor<KDRayWrap>,
1295  std::less<KDTree::_Bracket_accessor<KDRayWrap>::result_type>>
1296  range(wrap);
1297  scaltype limitX =
1298  errorXcoeff
1299  * sqrt(clusterXLimit2 + 2.0 * ray->source->dx * ray->source->dx);
1300  scaltype limitY =
1301  errorYcoeff
1302  * sqrt(clusterYLimit2 + 2.0 * ray->source->dy * ray->source->dy);
1303  scaltype limitTx =
1304  errorTxCoeff * sqrt(clusterTxLimit2 + 2.0 * ray->dtx * ray->dtx);
1305  scaltype limitTy =
1306  errorTyCoeff * sqrt(clusterTyLimit2 + 2.0 * ray->dty * ray->dty);
1307  KDRayWrap boundsWrap(wrap.data[0] - limitX,
1308  wrap.data[1] - limitY,
1309  wrap.data[2] - limitTx,
1310  wrap.data[3] - limitTy);
1311  range.set_low_bound(boundsWrap, 0);
1312  range.set_low_bound(boundsWrap, 1);
1313  range.set_low_bound(boundsWrap, 2);
1314  range.set_low_bound(boundsWrap, 3);
1315  boundsWrap.data[0] = wrap.data[0] + limitX;
1316  boundsWrap.data[1] = wrap.data[1] + limitY;
1317  boundsWrap.data[2] = wrap.data[2] + limitTx;
1318  boundsWrap.data[3] = wrap.data[3] + limitTy;
1319  range.set_high_bound(boundsWrap, 0);
1320  range.set_high_bound(boundsWrap, 1);
1321  range.set_high_bound(boundsWrap, 2);
1322  range.set_high_bound(boundsWrap, 3);
1323  list<KDRayWrap> neighbours;
1324  rays->find_within_range(range,
1325  back_insert_iterator<list<KDRayWrap>>(neighbours));
1326  Int_t level_occupancy[LXLAYERS * LXLAYERS] = {};
1327 
1328  for (list<KDRayWrap>::iterator j = neighbours.begin();
1329  j != neighbours.end();
1330  ++j) {
1331  KDRayWrap& w = *j;
1332  LxRay* r = w.ray;
1333  level_occupancy[r->level] = 1;
1334 
1335  if (ray != r) ray->neighbourhood.push_back(r);
1336  }
1337 
1338  Int_t levels = 0;
1339 
1340  for (Int_t j = 0; j < LXLAYERS * LXLAYERS; ++j)
1341  levels += level_occupancy[j];
1342 
1343  InsertClusterRay(levels, ray->neighbourhood.size() + 1, ray);
1344  } // for (KDRaysStorageType::iterator i = rays->begin(); i != rays->end(); ++i)
1345 
1346  KDRaysStorageType* clusteredRays =
1347  static_cast<KDRaysStorageType*>(clusteredRaysHandle);
1348 
1349  for (Int_t i = 2 * LXLAYERS - 1; i >= 0; --i) {
1350 #ifdef DENSE_CLUSTERS_FIRST
1351  for (Int_t j = clusters[i].size() - 1; j >= 0; --j)
1352 #else //DENSE_CLUSTERS_FIRST
1353  for (Int_t j = 0; j < clusters[i].size(); ++j)
1354 #endif //DENSE_CLUSTERS_FIRST
1355  {
1356  list<LxRay*>* clusterRays = clusters[i][j];
1357 
1358  for (list<LxRay*>::iterator k = clusterRays->begin();
1359  k != clusterRays->end();
1360  ++k) {
1361  LxRay* ray = *k;
1362 
1363  if (0 == ray) continue;
1364 
1365 #ifdef BEST_SIX_POINTS
1366  if (ray->source->used || ray->end->used)
1367  //if (ray->used)
1368  continue;
1369 #else //BEST_SIX_POINTS
1370  if (ray->used) continue;
1371 #endif //BEST_SIX_POINTS
1372 
1373  Int_t levels = 0;
1374  Int_t level_occupancy[LXLAYERS * LXLAYERS] = {};
1375  scaltype lX = 0;
1376  scaltype lY = 0;
1377  scaltype lDx2 = 0;
1378  scaltype lDy2 = 0;
1379  scaltype lZ = ray->end->layer->station->zCoord;
1380  scaltype rX = 0;
1381  scaltype rY = 0;
1382  scaltype rDx2 = 0;
1383  scaltype rDy2 = 0;
1384  scaltype rZ = ray->source->layer->station->zCoord;
1385  scaltype deltaZ = rZ - lZ;
1386  AddRayData(ray, lX, lY, lDx2, lDy2, rX, rY, rDx2, rDy2);
1387  level_occupancy[ray->level] = 1;
1388  Int_t numRays = 1;
1389 
1390 #ifdef BEST_RAYS_ONLY
1391  pair<LxRay*, scaltype> bestNeighbours
1392  [LXLAYERS * LXLAYERS]; // Neigbours and chi2 to the cluster owner.
1393 
1394  for (Int_t l = 0; l < LXLAYERS * LXLAYERS; ++l) {
1395  bestNeighbours[l].first = 0;
1396  bestNeighbours[l].second = 0;
1397  }
1398 
1399 #ifdef BEST_SIX_POINTS
1400  set<LxPoint*> leftBestPoints[LXLAYERS];
1401  set<LxPoint*> rightBestPoints[LXLAYERS];
1402  map<pair<LxPoint*, LxPoint*>, LxRay*, LxLess>
1403  bestRaysMap[LXLAYERS * LXLAYERS];
1404  LxPoint* rayLeftPoint = ray->end;
1405  Int_t rayLeftInd = rayLeftPoint->layer->layerNumber;
1406  leftBestPoints[rayLeftInd].insert(rayLeftPoint);
1407  LxPoint* rayRightPoint = ray->source;
1408  Int_t rayRightInd = rayRightPoint->layer->layerNumber;
1409  rightBestPoints[rayRightInd].insert(rayRightPoint);
1410  bestRaysMap[ray->level][make_pair(rayLeftPoint, rayRightPoint)] = ray;
1411 
1412  for (list<LxRay*>::iterator l = ray->neighbourhood.begin();
1413  l != ray->neighbourhood.end();
1414  ++l) {
1415  LxRay* r = *l;
1416 
1417  /*if (r->used)
1418  continue;
1419 
1420  if (r->level == ray->level)
1421  continue;*/
1422 
1423  LxPoint* rLeftPoint = r->end;
1424  Int_t rLeftInd = rLeftPoint->layer->layerNumber;
1425 
1426  if (rLeftPoint->used) continue;
1427 
1428  LxPoint* rRightPoint = r->source;
1429  Int_t rRightInd = rRightPoint->layer->layerNumber;
1430 
1431  if (rRightPoint->used) continue;
1432 
1433  if (rLeftInd == rayLeftInd && rRightInd == rayRightInd) continue;
1434 
1435  leftBestPoints[rLeftInd].insert(rLeftPoint);
1436  rightBestPoints[rRightInd].insert(rRightPoint);
1437  bestRaysMap[r->level][make_pair(rLeftPoint, rRightPoint)] = r;
1438  }
1439 #endif //BEST_SIX_POINTS
1440 
1441 #endif //BEST_RAYS_ONLY
1442 
1443 #ifdef BEST_SIX_POINTS
1444  LxRaysCandidates bestRaysCandidates;
1445 
1446  for (set<LxPoint*>::iterator r0 = rightBestPoints[0].begin();
1447  r0 != rightBestPoints[0].end();
1448  ++r0) {
1449  for (set<LxPoint*>::iterator r1 = rightBestPoints[1].begin();
1450  r1 != rightBestPoints[1].end();
1451  ++r1) {
1452  for (set<LxPoint*>::iterator r2 = rightBestPoints[2].begin();
1453  r2 != rightBestPoints[2].end();
1454  ++r2) {
1455  for (set<LxPoint*>::iterator l0 = leftBestPoints[0].begin();
1456  l0 != leftBestPoints[0].end();
1457  ++l0) {
1458  for (set<LxPoint*>::iterator l1 = leftBestPoints[1].begin();
1459  l1 != leftBestPoints[1].end();
1460  ++l1) {
1461  for (set<LxPoint*>::iterator l2 = leftBestPoints[2].begin();
1462  l2 != leftBestPoints[2].end();
1463  ++l2) {
1464  LxRaysCandidates cand;
1465 
1466  for (Int_t rInd = 0; rInd < LXLAYERS; ++rInd) {
1467  LxPoint* rPoint = 0;
1468 
1469  switch (rInd) {
1470  case 0: rPoint = *r0; break;
1471 
1472  case 1: rPoint = *r1; break;
1473 
1474  default: rPoint = *r2;
1475  }
1476 
1477  for (Int_t lInd = 0; lInd < LXLAYERS; ++lInd) {
1478  LxPoint* lPoint = 0;
1479 
1480  switch (lInd) {
1481  case 0: lPoint = *l0; break;
1482 
1483  case 1: lPoint = *l1; break;
1484 
1485  default: lPoint = *l2;
1486  }
1487 
1488  Int_t level = rPoint->layer->layerNumber * LXLAYERS
1489  + lPoint->layer->layerNumber;
1490  map<pair<LxPoint*, LxPoint*>, LxRay*, LxLess>::iterator
1491  rIter =
1492  bestRaysMap[level].find(make_pair(lPoint, rPoint));
1493 
1494  if (rIter == bestRaysMap[level].end()) continue;
1495 
1496  if (level == ray->level) continue;
1497 
1498  cand.data[level] = rIter->second;
1499  }
1500  }
1501 
1502  Int_t candLevels = 0;
1503  cand.chi2 = 0;
1504 
1505  for (Int_t l = 0; l < LXLAYERS * LXLAYERS; ++l) {
1506  if (0 == cand.data[l]) continue;
1507 
1508  ++candLevels;
1509  LxRay* r = cand.data[l];
1510  scaltype diffTx = r->tx - ray->tx;
1511  scaltype diffTy = r->ty - ray->ty;
1512  scaltype diffZ = ray->source->z - r->source->z;
1513  scaltype diffX =
1514  r->source->x + r->tx * diffZ - ray->source->x;
1515  scaltype diffY =
1516  r->source->y + r->ty * diffZ - ray->source->y;
1517  cand.chi2 +=
1518  diffTx * diffTx / (ray->dtx * ray->dtx)
1519  + diffTy * diffTy / (ray->dty * ray->dty)
1520  + diffX * diffX / (ray->source->dx * ray->source->dx)
1521  + diffY * diffY / (ray->source->dy * ray->source->dy);
1522  }
1523 
1524  if (candLevels < LXLAYERS * LXLAYERS - 1) continue;
1525 
1526  if (bestRaysCandidates.chi2 < 0
1527  || bestRaysCandidates.chi2 > cand.chi2)
1528  bestRaysCandidates = cand;
1529  }
1530  }
1531  }
1532  }
1533  }
1534  }
1535 
1536  if (bestRaysCandidates.chi2 < 0) continue;
1537 
1538  levels = LXLAYERS * LXLAYERS;
1539 
1540  for (Int_t l = 0; l < LXLAYERS * LXLAYERS; ++l) {
1541  bestNeighbours[l].first = bestRaysCandidates.data[l];
1542  bestNeighbours[l].second = 0;
1543  }
1544 
1545 #else //BEST_SIX_POINTS
1546  for (list<LxRay*>::iterator l = ray->neighbourhood.begin();
1547  l != ray->neighbourhood.end();
1548  ++l) {
1549  LxRay* r = *l;
1550 
1551  if (r->used) continue;
1552 
1553 #ifdef BEST_RAYS_ONLY
1554  if (r->level == ray->level) continue;
1555 #endif //BEST_RAYS_ONLY
1556 
1557  scaltype diffTx = r->tx - ray->tx;
1558  scaltype diffTy = r->ty - ray->ty;
1559  scaltype diffZ = ray->source->z - r->source->z;
1560  scaltype diffX = r->source->x + r->tx * diffZ - ray->source->x;
1561  scaltype diffY = r->source->y + r->ty * diffZ - ray->source->y;
1562  scaltype aChi2 =
1563  diffTx * diffTx / (ray->dtx * ray->dtx)
1564  + diffTy * diffTy / (ray->dty * ray->dty)
1565  + diffX * diffX / (ray->source->dx * ray->source->dx)
1566  + diffY * diffY / (ray->source->dy * ray->source->dy);
1567 
1568 #ifdef BEST_RAYS_ONLY
1569  if (0 == bestNeighbours[r->level].first
1570  || aChi2 < bestNeighbours[r->level].second) {
1571  bestNeighbours[r->level].first = r;
1572  bestNeighbours[r->level].second = aChi2;
1573  }
1574 #else //BEST_RAYS_ONLY
1575  AddRayData(r, lX, lY, lDx2, lDy2, rX, rY, rDx2, rDy2);
1576 #endif //BEST_RAYS_ONLY
1577 
1578  level_occupancy[r->level] = 1;
1579  ++numRays;
1580  }
1581 
1582  for (Int_t l = 0; l < LXLAYERS * LXLAYERS; ++l)
1583  levels += level_occupancy[l];
1584 #endif //BEST_SIX_POINTS
1585 
1586  if (levels < LXLAYERS * LXLAYERS /*(LXLAYERS - 1) * (LXLAYERS - 1)*/)
1587  continue;
1588  else if (levels < i + (LXLAYERS - 1) * (LXLAYERS - 1)
1589 #if defined(DENSE_CLUSTERS_FIRST) && !defined(BEST_RAYS_ONLY)
1590  || numRays < j + 1
1591 #endif //defined(DENSE_CLUSTERS_FIRST) && !defined(BEST_RAYS_ONLY)
1592  ) {
1593  InsertClusterRay(
1594  levels,
1595  numRays,
1596  ray); // It should be safe to leave ray also in its previous place.
1597  continue;
1598  }
1599 
1600 #ifdef BEST_RAYS_ONLY
1601  numRays = 1;
1602 
1603  for (Int_t l = 0; l < LXLAYERS * LXLAYERS; ++l) {
1604  LxRay* r = bestNeighbours[l].first;
1605 
1606  if (0 == r) continue;
1607 
1608  AddRayData(r, lX, lY, lDx2, lDy2, rX, rY, rDx2, rDy2);
1609 #ifdef BEST_SIX_POINTS
1610  r->source->used = true;
1611  r->end->used = true;
1612 #else //BEST_SIX_POINTS
1613  r->used = true;
1614 #endif //BEST_SIX_POINTS
1615  ++numRays;
1616  }
1617 #endif //BEST_RAYS_ONLY
1618 
1619  lX /= numRays;
1620  lY /= numRays;
1621  scaltype lDx = ray->end->dx; //sqrt(lDx2) / numRays;
1622  scaltype lDy = ray->end->dy; //sqrt(lDy2) / numRays;
1623  rX /= numRays;
1624  rY /= numRays;
1625  scaltype rDx = ray->source->dx; //sqrt(rDx2) / numRays;
1626  scaltype rDy = ray->source->dy; //sqrt(rDy2) / numRays;
1627 
1628  LxPoint* lPoint =
1629  new LxPoint(lX, lY, lZ, lDx, lDy, 0, ray->end->layer, -1);
1630  clusteredPoints.push_back(lPoint);
1631  LxPoint* rPoint =
1632  new LxPoint(rX, rY, rZ, rDx, rDy, 0, ray->source->layer, -1);
1633  clusteredPoints.push_back(rPoint);
1634  LxRay* clRay = new LxRay(rPoint, lPoint, LXLAYERS + LXMIDDLE);
1635  clRay->clusterPoints.push_back(ray->source);
1636  clRay->clusterPoints.push_back(ray->end);
1637 #ifdef REMEMBER_CLUSTERED_RAYS_IN_POINTS
1638  ray->source->leftClusteredRay = clRay;
1639  ray->end->rightClusteredRay = clRay;
1640 #endif // REMEMBER_CLUSTERED_RAYS_IN_POINTS
1641 #ifdef BEST_SIX_POINTS
1642  ray->source->used = true;
1643  ray->end->used = true;
1644  //ray->used = true;
1645 #else //BEST_SIX_POINTS
1646  ray->used = true;
1647 #endif //BEST_SIX_POINTS
1648 
1649 #ifdef BEST_RAYS_ONLY
1650  for (Int_t l = 0; l < LXLAYERS * LXLAYERS; ++l) {
1651  LxRay* r = bestNeighbours[l].first;
1652 
1653  if (0 == r) continue;
1654 
1655  clRay->clusterPoints.push_back(r->source);
1656  clRay->clusterPoints.push_back(r->end);
1657 #ifdef REMEMBER_CLUSTERED_RAYS_IN_POINTS
1658  ray->source->leftClusteredRay = clRay;
1659  ray->end->rightClusteredRay = clRay;
1660 #endif // REMEMBER_CLUSTERED_RAYS_IN_POINTS
1661  }
1662 
1663 #ifdef REMOVE_SUBCLUSTER
1664  for (list<LxRay*>::iterator l = ray->neighbourhood.begin();
1665  l != ray->neighbourhood.end();
1666  ++l) {
1667  LxRay* r = *l;
1668 
1669 #ifdef BEST_SIX_POINTS
1670  //if (r->source->used || r->end->used)
1671  //continue;
1672 #else //BEST_SIX_POINTS
1673  if (r->used) continue;
1674 #endif //BEST_SIX_POINTS
1675 
1676  /*scaltype diffTx = r->tx - ray->tx;
1677 
1678  if (diffTx < 0)
1679  diffTx = -diffTx;
1680 
1681  if (diffTx > errorTxCoeff * ray->dtx * 2)
1682  continue;
1683 
1684  scaltype diffTy = r->ty - ray->ty;
1685 
1686  if (diffTy < 0)
1687  diffTy = -diffTy;
1688 
1689  if (diffTy > errorTyCoeff * ray->dty * 2)
1690  continue;
1691 
1692  scaltype diffZ = ray->source->z - r->source->z;
1693  scaltype diffX = ray->source->x - r->source->x - r->tx * diffZ;
1694 
1695  if (diffX < 0)
1696  diffX = -diffX;
1697 
1698  if (diffX > errorXcoeff * ray->source->dx * 2)
1699  continue;
1700 
1701  scaltype diffY = ray->source->y - r->source->y - r->ty * diffZ;
1702 
1703  if (diffY < 0)
1704  diffY = -diffY;
1705 
1706  if (diffY > errorYcoeff * ray->source->dy * 2)
1707  continue;*/
1708 
1709  clRay->clusterPoints.push_back(r->source);
1710  clRay->clusterPoints.push_back(r->end);
1711 #ifdef BEST_SIX_POINTS
1712  //r->source->used = true;
1713  //r->end->used = true;
1714 #else //BEST_SIX_POINTS
1715  r->used = true;
1716 #endif //BEST_SIX_POINTS
1717  }
1718 #endif //REMOVE_SUBCLUSTER
1719 
1720 #else //BEST_RAYS_ONLY
1721  for (list<LxRay*>::iterator l = ray->neighbourhood.begin();
1722  l != ray->neighbourhood.end();
1723  ++l) {
1724  LxRay* r = *l;
1725 
1726  if (r->used) continue;
1727 
1728  clRay->clusterPoints.push_back(r->source);
1729  clRay->clusterPoints.push_back(r->end);
1730  r->used = true;
1731  }
1732 #endif //BEST_RAYS_ONLY
1733 
1734  KDRayWrap clRayWrap(rX, rY, clRay);
1735  clusteredRays->insert(clRayWrap);
1736  } // for (list<LxRay*>::iterator k = rays->begin(); k != rays->end(); ++k)
1737  } // for (vector<list<LxRay*>*>::reverse_iterator j = clusters[i].rbegin(); j != clusters[i].rend(); ++j)
1738  } // for (Int_t i = 2 * LXLAYERS - 1; i >= 0; --i)
1739 
1740  gettimeofday(&eTime, 0);
1741  exeDuration =
1742  (eTime.tv_sec - bTime.tv_sec) * 1000000 + eTime.tv_usec - bTime.tv_usec;
1743 
1744  if (exeDuration > maxExeDuration) maxExeDuration = exeDuration;
1745 
1746  cout << "Max execution duration was: " << maxExeDuration << endl;
1747 }
1748 #endif //CLUSTER_MODE
1749 
1750 void LxStation::ConnectNeighbours() {
1751  if (stationNumber < LXFIRSTSTATION + 2) return;
1752 
1753  LxLayer* layer = layers[LXMIDDLE];
1754  LxStation* lStation = space->stations[stationNumber - 1];
1755 
1756  for (list<LxPoint*>::iterator i = layer->points.begin();
1757  i != layer->points.end();
1758  ++i) {
1759  LxPoint* rPoint = *i;
1760 
1761  if (!rPoint->valid) continue;
1762 
1763  for (list<LxRay*>::iterator j = rPoint->rays.begin();
1764  j != rPoint->rays.end();
1765  ++j) {
1766  LxRay* rRay = *j;
1767  LxPoint* ePoint = rRay->end;
1768 
1769  if (!ePoint->valid) continue;
1770 
1771  for (list<LxRay*>::iterator k = ePoint->rays.begin();
1772  k != ePoint->rays.end();
1773  ++k) {
1774  LxRay* lRay = *k;
1775  scaltype diffTx = lRay->tx - rRay->tx;
1776  scaltype diffTy = lRay->ty - rRay->ty;
1777 
1778  if (diffTx < 0) diffTx = -diffTx;
1779 
1780  if (diffTy < 0) diffTy = -diffTy;
1781 
1782  scaltype dtxb = sqrt(lRay->dtx * lRay->dtx + rRay->dtx * rRay->dtx);
1783  scaltype dtyb = sqrt(lRay->dty * lRay->dty + rRay->dty * rRay->dty);
1784 
1785  if (diffTx > lStation->txBreakLimit + dtxb * errorTxBreakCoeff
1786  || diffTy > lStation->tyBreakLimit + dtyb * errorTyBreakCoeff)
1787  continue;
1788 
1789  rRay->neighbours.push_back(lRay);
1790  }
1791  }
1792  }
1793 }
1794 
1796 // LxSpace
1798 
1799 LxSpace::LxSpace()
1800  : muchStsBreakX(0)
1801  , muchStsBreakY(0)
1802  , muchStsBreakTx(0)
1803  , muchStsBreakTy(0)
1804  , stationsInAlgo(LXSTATIONS) {
1805  for (int i = 0; i < LXSTATIONS; ++i)
1806  stations.push_back(new LxStation(this, i));
1807 
1808  x_coords = reinterpret_cast<scal_coords*>(_mm_malloc(
1809  sizeof(scaltype) * LXSTATIONS * LXLAYERS * LXMAXPOINTSONSTATION, 32));
1810  tx_vals = reinterpret_cast<scal_tans*>(
1811  _mm_malloc(sizeof(scaltype) * (LXSTATIONS - 1) * LXMAXPOINTSONSTATION, 32));
1812  x_errs = reinterpret_cast<scal_coords*>(_mm_malloc(
1813  sizeof(scaltype) * LXSTATIONS * LXLAYERS * LXMAXPOINTSONSTATION, 32));
1814  y_coords = reinterpret_cast<scal_coords*>(_mm_malloc(
1815  sizeof(scaltype) * LXSTATIONS * LXLAYERS * LXMAXPOINTSONSTATION, 32));
1816  ty_vals = reinterpret_cast<scal_tans*>(
1817  _mm_malloc(sizeof(scaltype) * (LXSTATIONS - 1) * LXMAXPOINTSONSTATION, 32));
1818  y_errs = reinterpret_cast<scal_coords*>(_mm_malloc(
1819  sizeof(scaltype) * LXSTATIONS * LXLAYERS * LXMAXPOINTSONSTATION, 32));
1820  z_coords = reinterpret_cast<scal_coords*>(_mm_malloc(
1821  sizeof(scaltype) * LXSTATIONS * LXLAYERS * LXMAXPOINTSONSTATION, 32));
1822 }
1823 
1824 LxSpace::~LxSpace() {
1825  _mm_free(x_coords);
1826  _mm_free(tx_vals);
1827  _mm_free(x_errs);
1828  _mm_free(y_coords);
1829  _mm_free(ty_vals);
1830  _mm_free(y_errs);
1831  _mm_free(z_coords);
1832 }
1833 
1834 void LxSpace::Clear() {
1835  for (vector<LxStation*>::iterator i = stations.begin(); i != stations.end();
1836  ++i)
1837  (*i)->Clear();
1838 
1839  for (list<LxTrack*>::iterator i = tracks.begin(); i != tracks.end(); ++i)
1840  delete *i;
1841 
1842  tracks.clear();
1843  extTracks.clear();
1844 }
1845 
1846 void LxSpace::RestoreMiddlePoints() {
1847  for (int i = LXFIRSTSTATION; i < LXSTATIONS; ++i)
1848  stations[i]->RestoreMiddlePoints();
1849 }
1850 
1851 void LxSpace::BuildRays() {
1852  for (int i = LXFIRSTSTATION + 1; i < LXSTATIONS; ++i)
1853  stations[i]->BuildRays();
1854 }
1855 
1856 #ifdef CLUSTER_MODE
1857 void LxSpace::BuildRays2() {
1858  for (Int_t i = LXSTATIONS - 1; i > LXFIRSTSTATION; --i) {
1859  stations[i]->BuildRays2();
1860 
1861  if (i <= LXFIRSTSTATION + 1) continue;
1862 
1863  for (Int_t j = 0; j < LXLAYERS; ++j) {
1864  LxLayer* layer = stations[i - 1]->layers[j];
1865 
1866  for (list<LxPoint*>::iterator k = layer->points.begin();
1867  k != layer->points.end();
1868  ++k) {
1869  LxPoint* point = *k;
1870  point->used = false;
1871  }
1872  }
1873  }
1874 }
1875 
1876 void LxSpace::ConnectNeighbours2() {
1877  for (Int_t i = LXSTATIONS - 1; i > LXFIRSTSTATION + 1; --i) {
1878  KDRaysStorageType* leftRays =
1879  static_cast<KDRaysStorageType*>(stations[i - 1]->clusteredRaysHandle);
1880  KDRaysStorageType* rightRays =
1881  static_cast<KDRaysStorageType*>(stations[i]->clusteredRaysHandle);
1882  scaltype txSigma = stations[i - 1]->txBreakSigma;
1883  scaltype txSigma2 = txSigma * txSigma;
1884  scaltype tySigma = stations[i - 1]->tyBreakSigma;
1885  scaltype tySigma2 = tySigma * tySigma;
1886  scaltype deltaZ = stations[i - 1]->zCoord - stations[i]->zCoord;
1887  scaltype deltaZ2 = deltaZ * deltaZ;
1888  scaltype deltaZ23 = 1.0; //deltaZ2 / 3.0;
1889 
1890  for (KDRaysStorageType::iterator j = rightRays->begin();
1891  j != rightRays->end();
1892  ++j) {
1893  KDRayWrap& wrap = const_cast<KDRayWrap&>(*j);
1894  LxRay* ray = wrap.ray;
1895  KDTree::_Region<
1896  4,
1897  KDRayWrap,
1898  scaltype,
1899  KDTree::_Bracket_accessor<KDRayWrap>,
1900  std::less<KDTree::_Bracket_accessor<KDRayWrap>::result_type>>
1901  range(wrap);
1902  LxPoint* point = ray->end;
1903  scaltype x = wrap.data[0] + wrap.data[2] * deltaZ;
1904  scaltype y = wrap.data[1] + wrap.data[3] * deltaZ;
1905  scaltype tx = wrap.data[2];
1906  scaltype ty = wrap.data[3];
1907  scaltype xRange =
1908  4 * sqrt(2 * point->dx * point->dx + txSigma2 * deltaZ23);
1909  scaltype yRange =
1910  4 * sqrt(2 * point->dy * point->dy + tySigma2 * deltaZ23);
1911  scaltype txRange = 4 * sqrt(2 * ray->dtx * ray->dtx + txSigma2);
1912  scaltype tyRange = 4 * sqrt(2 * ray->dty * ray->dty + tySigma2);
1913  KDRayWrap boundsWrap(x - xRange, y - yRange, tx - txRange, ty - tyRange);
1914  range.set_low_bound(boundsWrap, 0);
1915  range.set_low_bound(boundsWrap, 1);
1916  range.set_low_bound(boundsWrap, 2);
1917  range.set_low_bound(boundsWrap, 3);
1918  boundsWrap.data[0] = x + xRange;
1919  boundsWrap.data[1] = y + yRange;
1920  boundsWrap.data[2] = tx + txRange;
1921  boundsWrap.data[3] = ty + tyRange;
1922  range.set_high_bound(boundsWrap, 0);
1923  range.set_high_bound(boundsWrap, 1);
1924  range.set_high_bound(boundsWrap, 2);
1925  range.set_high_bound(boundsWrap, 3);
1926  list<KDRayWrap> neighbours;
1927  leftRays->find_within_range(
1928  range, back_insert_iterator<list<KDRayWrap>>(neighbours));
1929 
1930  for (list<KDRayWrap>::iterator k = neighbours.begin();
1931  k != neighbours.end();
1932  ++k) {
1933  KDRayWrap& w = *k;
1934  LxRay* r = w.ray;
1935  ray->neighbours.push_back(r);
1936  }
1937  }
1938  }
1939 }
1940 
1941 void LxSpace::BuildCandidates2(LxRay* ray,
1942  LxRay** rays,
1943  list<LxTrackCandidate*>& candidates,
1944  scaltype chi2) {
1945  int level = ray->station->stationNumber - LXFIRSTSTATION - 1;
1946  rays[level] = ray;
1947 
1948  if (0 == level) {
1949  LxTrackCandidate* tc =
1950  new LxTrackCandidate(rays, LXSTATIONS - LXFIRSTSTATION - 1, chi2);
1951  candidates.push_back(tc);
1952  return;
1953  }
1954 
1955  LxPoint* point = ray->end;
1956 
1957  for (list<LxRay*>::iterator i = ray->neighbours.begin();
1958  i != ray->neighbours.end();
1959  ++i) {
1960  LxRay* lRay = *i;
1961 
1962  if (lRay->source->used) continue;
1963 
1964  LxPoint* lPoint = lRay->source;
1965  LxStation* station = lRay->station;
1966  scaltype dx = point->x - lPoint->x;
1967  scaltype dx2 = dx * dx;
1968  scaltype sigmaX2 = point->dx * point->dx + lPoint->dx * lPoint->dx;
1969  scaltype dy = point->y - lPoint->y;
1970  scaltype dy2 = dy * dy;
1971  scaltype sigmaY2 = point->dy * point->dy + lPoint->dy * lPoint->dy;
1972  scaltype dtx = ray->tx - lRay->tx;
1973  scaltype dtx2 = dtx * dtx;
1974  scaltype sigmaTx2 = ray->dtx * ray->dtx + lRay->dtx * lRay->dtx
1975  + station->txBreakSigma * station->txBreakSigma;
1976  scaltype dty = ray->ty - lRay->ty;
1977  scaltype dty2 = dty * dty;
1978  scaltype sigmaTy2 = ray->dty * ray->dty + lRay->dty * lRay->dty
1979  + station->tyBreakSigma * station->tyBreakSigma;
1980 
1981  // Continue process of track building.
1982  BuildCandidates2(lRay,
1983  rays,
1984  candidates,
1985  chi2 + dx2 / sigmaX2 + dy2 / sigmaY2 + dtx2 / sigmaTx2
1986  + dty2 / sigmaTy2);
1987  }
1988 }
1989 
1990 void LxSpace::Reconstruct2() {
1991  LxStation* startStation = stations[LXSTATIONS - 1];
1992  KDRaysStorageType* startRays =
1993  static_cast<KDRaysStorageType*>(startStation->clusteredRaysHandle);
1994 
1995  for (KDRaysStorageType::iterator i = startRays->begin();
1996  i != startRays->end();
1997  ++i) {
1998  KDRayWrap& wrap = const_cast<KDRayWrap&>(*i);
1999  LxRay* ray = wrap.ray;
2000  LxRay* rays[LXSTATIONS - LXFIRSTSTATION - 1];
2001  list<LxTrackCandidate*> candidates;
2002  scaltype chi2 = 0;
2003  BuildCandidates2(ray, rays, candidates, chi2);
2004  LxTrackCandidate* bestCandidate = 0;
2005 
2006  for (list<LxTrackCandidate*>::iterator j = candidates.begin();
2007  j != candidates.end();
2008  ++j) {
2009  LxTrackCandidate* candidate = *j;
2010 
2011  if (0 == bestCandidate || candidate->chi2 < bestCandidate->chi2)
2012  bestCandidate = candidate;
2013  }
2014 
2015  if (0 != bestCandidate) {
2016  LxTrack* track = new LxTrack(bestCandidate);
2017  tracks.push_back(track);
2018  }
2019 
2020  for (list<LxTrackCandidate*>::iterator j = candidates.begin();
2021  j != candidates.end();
2022  ++j)
2023  delete *j;
2024  } // for (KDRaysStorageType::iterator i = rays->begin(); i != rays->end(); ++i)
2025 
2026  cout << "LxSpace::Reconstruct() found: " << tracks.size() << " tracks"
2027  << endl;
2028 }
2029 #endif //CLUSTER_MODE
2030 
2031 void LxSpace::ConnectNeighbours() {
2032  for (int i = LXFIRSTSTATION + 2; i < LXSTATIONS; ++i)
2033  stations[i]->ConnectNeighbours();
2034 }
2035 
2036 void LxSpace::BuildCandidates(int endStNum,
2037  LxRay* ray,
2038  LxRay** rays,
2039  list<LxTrackCandidate*>& candidates,
2040  scaltype chi2) {
2041  int level = ray->station->stationNumber - 1;
2042  rays[level] = ray;
2043 
2044  if (0 == level) {
2045 #ifdef USE_KALMAN
2046  LxTrackCandidate* tc =
2047  new LxTrackCandidate(rays, LXSTATIONS - 1, ray->kalman.chi2);
2048 #else // USE_KALMAN
2049  LxTrackCandidate* tc = new LxTrackCandidate(rays, endStNum, chi2);
2050 #endif //USE_KALMAN
2051  candidates.push_back(tc);
2052  return;
2053  }
2054 
2055  for (list<LxRay*>::iterator i = ray->neighbours.begin();
2056  i != ray->neighbours.end();
2057  ++i) {
2058  LxRay* lRay = *i;
2059 
2060  if (lRay->source->used) continue;
2061 
2062  LxStation* station = lRay->station;
2063  //scaltype dtx = ray->tx - lRay->tx;
2064  //dtx /= station->txBreakSigma;
2065  //scaltype dty = ray->ty - lRay->ty;
2066  //dty /= station->tyBreakSigma;
2067 
2068  scaltype dtx = ray->tx - lRay->tx;
2069  scaltype dtx2 = dtx * dtx;
2070  scaltype sigmaTx2 = ray->dtx * ray->dtx + lRay->dtx * lRay->dtx
2071  + station->txBreakSigma * station->txBreakSigma;
2072  scaltype dty = ray->ty - lRay->ty;
2073  scaltype dty2 = dty * dty;
2074  scaltype sigmaTy2 = ray->dty * ray->dty + lRay->dty * lRay->dty
2075  + station->tyBreakSigma * station->tyBreakSigma;
2076 
2077 #ifdef USE_KALMAN
2078  LxKalmanParams& kPrev = ray->kalman;
2079  LxKalmanParams& kalman = lRay->kalman;
2080  kalman.tx = kPrev.tx;
2081  kalman.ty = kPrev.ty;
2082  kalman.C11 = kPrev.C11 + station->txBreakSigma * station->txBreakSigma;
2083  kalman.C22 = kPrev.C22 + station->tyBreakSigma * station->tyBreakSigma;
2084  scaltype S11 = 1 / (kalman.C11 + lRay->dtx * lRay->dtx);
2085  scaltype S22 = 1 / (kalman.C22 + lRay->dty * lRay->dty);
2086  scaltype K11 = kalman.C11 * S11;
2087  scaltype K22 = kalman.C22 * S22;
2088  scaltype zetaTx = lRay->tx - kalman.tx;
2089  scaltype zetaTy = lRay->ty - kalman.ty;
2090  kalman.tx += K11 * zetaTx;
2091  kalman.ty += K22 * zetaTy;
2092  kalman.C11 = (1.0 - K11) * kalman.C11;
2093  kalman.C22 = (1.0 - K22) * kalman.C22;
2094  kalman.chi2 = kPrev.chi2 + zetaTx * S11 * zetaTx + zetaTy * S22 * zetaTy;
2095 #endif //USE_KALMAN
2096 
2097  // Continue process of track building.
2098  //BuildCandidates(lRay, rays, candidates, chi2 + dtx * dtx + dty * dty);
2099  BuildCandidates(endStNum,
2100  lRay,
2101  rays,
2102  candidates,
2103  chi2 + dtx2 / sigmaTx2 + dty2 / sigmaTy2);
2104  }
2105 }
2106 
2107 void LxSpace::Reconstruct() {
2108  for (int endStNum = LXSTATIONS - 1; endStNum >= stationsInAlgo - 1;
2109  --endStNum) {
2110  LxStation* startStation = stations[endStNum];
2111  LxLayer* startLayer = startStation->layers[LXMIDDLE];
2112  list<LxPoint*>& startPoints = startLayer->points;
2113 
2114  for (list<LxPoint*>::iterator i = startPoints.begin();
2115  i != startPoints.end();
2116  ++i) {
2117  LxPoint* point = *i;
2118 
2119  if (point->used) continue;
2120 
2121  if (!point->valid) continue;
2122 
2123  LxRay* rays[endStNum];
2124  list<LxTrackCandidate*> candidates;
2125  list<LxRay*>& startRays = point->rays;
2126 
2127  for (list<LxRay*>::iterator j = startRays.begin(); j != startRays.end();
2128  ++j) {
2129  LxRay* ray = *j;
2130  scaltype chi2 = 0;
2131 
2132 #ifdef USE_KALMAN
2133  LxKalmanParams& kalman = ray->kalman;
2134  kalman.tx = ray->tx;
2135  kalman.ty = ray->ty;
2136  kalman.C11 = ray->dtx * ray->dtx;
2137  kalman.C22 = ray->dty * ray->dty;
2138  kalman.chi2 = 0;
2139 #endif //USE_KALMAN
2140 
2141  BuildCandidates(endStNum, ray, rays, candidates, chi2);
2142  }
2143 
2144  LxTrackCandidate* bestCandidate = 0;
2145 
2146  for (list<LxTrackCandidate*>::iterator j = candidates.begin();
2147  j != candidates.end();
2148  ++j) {
2149  LxTrackCandidate* candidate = *j;
2150 
2151  if (0 == bestCandidate || candidate->chi2 < bestCandidate->chi2)
2152  bestCandidate = candidate;
2153  }
2154 
2155  if (0 != bestCandidate) {
2156  LxTrack* track = new LxTrack(bestCandidate);
2157  tracks.push_back(track);
2158  }
2159 
2160  for (list<LxTrackCandidate*>::iterator j = candidates.begin();
2161  j != candidates.end();
2162  ++j)
2163  delete *j;
2164  }
2165  } // for (int stNum = LXSTATIONS - 1; stNum >= stationsInAlgo - 1; --stNum)
2166 
2167  for (list<LxTrack*>::iterator i = tracks.begin(); i != tracks.end(); ++i) {
2168  LxTrack* track = *i;
2169  //cout << "LxSpace::Reconstruct(): found track with length = " << track->length << endl << "With points:";
2170 
2171  for (int j = 0; j < track->length * LXLAYERS; ++j) {
2172  LxPoint* point = track->points[j];
2173 
2174  //if (point)
2175  //cout << " " << point->z;
2176  //else
2177  //cout << " *";
2178  }
2179 
2180  //cout << endl;
2181  }
2182 
2183  RemoveClones();
2184 }
2185 
2186 void LxSpace::RemoveClones() {
2187  for (list<LxTrack*>::iterator i = tracks.begin(); i != tracks.end(); ++i) {
2188  LxTrack* firstTrack = *i;
2189  list<LxTrack*>::iterator i2 = i;
2190  ++i2;
2191 
2192  if (i2 == tracks.end()) break;
2193 
2194  LxTrack* secondTrack = *i2;
2195  Int_t neighbourPoints = 0;
2196  int minLength = firstTrack->length < secondTrack->length
2197  ? firstTrack->length
2198  : secondTrack->length;
2199 
2200  for (Int_t j = 0; j < minLength * LXLAYERS; ++j) {
2201  LxPoint* point1 = firstTrack->points[j];
2202  LxPoint* point2 = secondTrack->points[j];
2203 
2204  if (0 == point1 || 0 == point2) continue;
2205 
2206  scaltype dx = point1->dx > point2->dx ? point1->dx : point2->dx;
2207  scaltype dy = point1->dy > point2->dy ? point1->dy : point2->dy;
2208 
2209  if (abs(point2->x - point1->x) < 5.0 * dx
2210  && abs(point2->y - point1->y) < 5.0 * dy)
2211  ++neighbourPoints;
2212  }
2213 
2214  if (neighbourPoints < minLength * LXLAYERS / 2) continue;
2215 
2216  if (firstTrack->length > secondTrack->length)
2217  secondTrack->clone = true;
2218  else if (secondTrack->length > firstTrack->length)
2219  firstTrack->clone = true;
2220  else if (firstTrack->chi2 < secondTrack->chi2)
2221  secondTrack->clone = true;
2222  else
2223  firstTrack->clone = true;
2224  }
2225 }
2226 
2227 void LxSpace::FitTracks() {
2228  for (list<LxTrack*>::iterator i = tracks.begin(); i != tracks.end(); ++i)
2229  (*i)->Fit();
2230 }
2231 
2232 void LxSpace::JoinExtTracks() {
2233  /*scaltype sigmaX = 0.8548;//0.89;//1.202;
2234  scaltype sigmaX2 = sigmaX * sigmaX;
2235  scaltype sigmaY = 0.6233;//1.061;
2236  scaltype sigmaY2 = sigmaY * sigmaY;
2237  scaltype sigmaTx = 0.02349;//0.02426;
2238  scaltype sigmaTx2 = sigmaTx * sigmaTx;
2239  scaltype sigmaTy = 0.007941;//0.01082;
2240  scaltype sigmaTy2 = sigmaTy * sigmaTy;
2241  scaltype deltaMuPlus = -0.01883;
2242  scaltype sigmaTxMuPlus = 0.01105;
2243  scaltype sigmaTxMuPlus2 = sigmaTxMuPlus * sigmaTxMuPlus;
2244  scaltype deltaMuMinus = 0.020;
2245  scaltype sigmaTxMuMinus = 0.0118;
2246  scaltype sigmaTxMuMinus2 = sigmaTxMuMinus * sigmaTxMuMinus;
2247  scaltype covXTx = 0.155612;
2248  scaltype covYTy = 0.157198;*/
2249  //scaltype cutCoeff = 5.0;
2250 #ifdef USE_OLD_STS_LINKING_RULE
2251  scaltype sigmaX2 = muchStsBreakX * muchStsBreakX;
2252  scaltype sigmaY2 = muchStsBreakY * muchStsBreakY;
2253  scaltype sigmaTx2 = muchStsBreakTx * muchStsBreakTx;
2254  scaltype sigmaTy2 = muchStsBreakTy * muchStsBreakTy;
2255 #endif //USE_OLD_STS_LINKING_RULE
2256 
2257  TrackPropagatorPtr fPropagator =
2258  CbmLitToolFactory::CreateTrackPropagator("lit");
2259 
2260  for (list<LxTrack*>::iterator i = tracks.begin(); i != tracks.end(); ++i) {
2261  LxTrack* track = *i;
2262 
2263  if (track->clone) continue;
2264 
2265 #ifdef USE_KALMAN_FIT
2266  scaltype x = track->x;
2267  scaltype y = track->y;
2268  scaltype z = track->z;
2269  scaltype tx0 = track->tx;
2270  scaltype ty = track->ty;
2271  scaltype dxMuch = track->dx;
2272  scaltype dyMuch = track->dy;
2273  scaltype dtxMuch = track->dtx;
2274  scaltype dtyMuch = track->dty;
2275 #else //USE_KALMAN_FIT
2276  LxRay* ray = track->rays[0];
2277  LxPoint* point = ray->end;
2278  scaltype x = point->x;
2279  scaltype y = point->y;
2280  scaltype z = point->z;
2281  //scaltype tx0 = ray->tx;
2282  scaltype tx = ray->tx;
2283  scaltype ty = ray->ty;
2284  scaltype dxMuch = point->dx;
2285  scaltype dyMuch = point->dy;
2286  scaltype dtxMuch = ray->dtx;
2287  scaltype dtyMuch = ray->dty;
2288 #endif //USE_KALMAN_FIT
2289 
2290  // External track are already filtered by P and Pt.
2291  for (list<LxExtTrack>::iterator j = extTracks.begin(); j != extTracks.end();
2292  ++j) {
2293  LxExtTrack* extTrack = &(*j);
2294  const FairTrackParam* firstParam = extTrack->track->GetParamFirst();
2295  const FairTrackParam* lastParam = extTrack->track->GetParamLast();
2296 
2297 #ifndef USE_OLD_STS_LINKING_RULE
2298  CbmLitTrackParam litLastParam;
2299  CbmLitConverterFairTrackParam::FairTrackParamToCbmLitTrackParam(
2300  lastParam, &litLastParam);
2301 
2302  if (kLITERROR
2303  == fPropagator->Propagate(&litLastParam, stations[0]->zCoord, 13))
2304  continue;
2305 
2306  scaltype deltaX = abs(litLastParam.GetX() - x);
2307  scaltype deltaY = abs(litLastParam.GetY() - y);
2308  scaltype deltaTx = abs(litLastParam.GetTx() - tx0);
2309  scaltype deltaTy = abs(litLastParam.GetTy() - ty);
2310  scaltype sigmaX2 = dxMuch * dxMuch + litLastParam.GetCovariance(0);
2311  scaltype sigmaX = sqrt(sigmaX2);
2312  scaltype sigmaY2 = dyMuch * dyMuch + litLastParam.GetCovariance(5);
2313  scaltype sigmaY = sqrt(sigmaY2);
2314  scaltype sigmaTx2 = dtxMuch * dtxMuch + litLastParam.GetCovariance(9);
2315  scaltype sigmaTx = sqrt(sigmaTx2);
2316  scaltype sigmaTy2 = dtyMuch * dtyMuch + litLastParam.GetCovariance(12);
2317  scaltype sigmaTy = sqrt(sigmaTy2);
2318 
2319  //if (deltaX > cutCoeff * sigmaX || deltaY > cutCoeff * sigmaY ||
2320  //deltaTx > cutCoeff * sigmaTx || deltaTy > cutCoeff * sigmaTy)
2321  //{
2322  //continue;
2323  //}
2324 
2325  scaltype chi2 = deltaX * deltaX / sigmaX2 + deltaY * deltaY / sigmaY2
2326  + deltaTx * deltaTx / sigmaTx2
2327  + deltaTy * deltaTy / sigmaTy2;
2328 #else //USE_OLD_STS_LINKING_RULE
2329 
2330  //if ((tx0 - lastParam->GetTx()) * (lastParam->GetTx() - firstParam->GetTx()) < 0)
2331  //continue;
2332 
2333  /*scaltype muchCharge = tx0 - firstParam->GetTx();
2334  scaltype tx = muchCharge > 0 ? tx0 + deltaMuPlus : tx0 + deltaMuMinus;
2335 
2336  if (muchCharge > 0)
2337  {
2338  sigmaTx = sigmaTxMuPlus;
2339  sigmaTx2 = sigmaTxMuPlus2;
2340  }
2341  else
2342  {
2343  sigmaTx = sigmaTxMuMinus;
2344  sigmaTx2 = sigmaTxMuMinus2;
2345  }*/
2346 
2347  scaltype deltaZ = lastParam->GetZ() - z;
2348  scaltype deltaX = abs(lastParam->GetX() - x - tx * deltaZ);
2349  scaltype dxSts2 = lastParam->GetCovariance(0, 0);
2350  scaltype dySts2 = lastParam->GetCovariance(1, 1);
2351  scaltype dtxSts2 = lastParam->GetCovariance(2, 2);
2352  scaltype dtySts2 = lastParam->GetCovariance(3, 3);
2353  //scaltype sigmaXMeas2 = dxMuch * dxMuch + dxSts2 - covXTx * deltaZ;// deltaZ is negative.
2354  scaltype sigmaXMeas2 =
2355  dxMuch * dxMuch + dxSts2 + dtxMuch * dtxMuch * deltaZ * deltaZ;
2356  //scaltype sigmaXMeas = sqrt(sigmaXMeas2);
2357  //scaltype sigmaYMeas2 = dyMuch * dyMuch + dySts2 - covYTy * deltaZ;// deltaZ is negative.
2358  scaltype sigmaYMeas2 =
2359  dyMuch * dyMuch + dySts2 + dtyMuch * dtyMuch * deltaZ * deltaZ;
2360  //scaltype sigmaYMeas = sqrt(sigmaYMeas2);
2361  scaltype sigmaTxMeas2 = dtxMuch * dtxMuch + dtxSts2;
2362  //scaltype sigmaTxMeas = sqrt(sigmaTxMeas2);
2363  scaltype sigmaTyMeas2 = dtyMuch * dtyMuch + dtySts2;
2364  //scaltype sigmaTyMeas = sqrt(sigmaTyMeas2);
2365 
2366  //if (deltaX > cutCoeff * sqrt(sigmaX2 + sigmaXMeas2))
2367  //continue;
2368 
2369  scaltype deltaY = abs(lastParam->GetY() - y - ty * deltaZ);
2370 
2371  //if (deltaY > cutCoeff * sqrt(sigmaY2 + sigmaYMeas2))
2372  //continue;
2373 
2374  scaltype deltaTx = abs(lastParam->GetTx() - tx);
2375 
2376  //if (deltaTx > cutCoeff * sqrt(sigmaTx2 + sigmaTxMeas2))
2377  //continue;
2378 
2379  scaltype deltaTy = abs(lastParam->GetTy() - ty);
2380 
2381  //if (deltaTy > cutCoeff * sqrt(sigmaTy2 + sigmaTyMeas2))
2382  //continue;
2383 
2384  // Take the charge sign into account.
2385  //scaltype stsCharge = firstParam->GetQp();
2386  //scaltype angMomInv = muchCharge / stsCharge;
2387  //scaltype dAmi = abs(sqrt(dtxMuch * dtxMuch + firstParam->GetCovariance(2, 2)) / stsCharge);
2388 
2389  // Check if the MUCH track projection to XZ plane angle fit the STS track momentum.
2390  //if (0.18 - dAmi > angMomInv || 0.52 + dAmi < angMomInv)
2391  //if (0.26 - dAmi > angMomInv || 0.44 + dAmi < angMomInv)
2392  //continue;
2393 
2394  scaltype chi2 = deltaX * deltaX / (sigmaX2 + sigmaXMeas2)
2395  + deltaY * deltaY / (sigmaY2 + sigmaYMeas2)
2396  + deltaTx * deltaTx / (sigmaTx2 + sigmaTxMeas2)
2397  + deltaTy * deltaTy / (sigmaTy2 + sigmaTyMeas2);
2398 #endif //USE_OLD_STS_LINKING_RULE
2399 
2400 #ifdef LX_EXT_LINK_SOPH
2401  list<pair<LxExtTrack*, scaltype>>::iterator k =
2402  track->extTrackCandidates.begin();
2403 
2404  for (; k != track->extTrackCandidates.end() && chi2 >= k->second; ++k)
2405  ;
2406 
2407  pair<LxExtTrack*, scaltype> linkDesc(extTrack, chi2);
2408  track->extTrackCandidates.insert(k, linkDesc);
2409 #else //LX_EXT_LINK_SOPH
2410  if (0 == track->externalTrack || track->extLinkChi2 > chi2) {
2411  track->externalTrack = extTrack;
2412  track->extLinkChi2 = chi2;
2413  }
2414 #endif //LX_EXT_LINK_SOPH
2415  } // for (list<LxExtTrack>::iterator j = extTracks.begin(); j != extTracks.end(); ++j)
2416 
2417 #ifdef LX_EXT_LINK_SOPH
2418  for (list<pair<LxExtTrack*, scaltype>>::iterator j =
2419  track->extTrackCandidates.begin();
2420  j != track->extTrackCandidates.end();
2421  ++j) {
2422  LxExtTrack* extTrack = j->first;
2423  scaltype chi2 = j->second;
2424 
2425  if (0 == extTrack->recoTrack.first) {
2426  extTrack->recoTrack.first = track;
2427  extTrack->recoTrack.second = chi2;
2428  track->externalTrack = extTrack;
2429  break;
2430  } else if (chi2 < extTrack->recoTrack.second) {
2431  LxTrack* anotherTrack = extTrack->recoTrack.first;
2432  extTrack->recoTrack.first = track;
2433  extTrack->recoTrack.second = chi2;
2434  track->externalTrack = extTrack;
2435  anotherTrack->Rebind();
2436  break;
2437  }
2438  }
2439 #endif // LX_EXT_LINK_SOPH
2440  } // for (list<LxTrack*>::iterator i = tracks.begin(); i != tracks.end(); ++i)
2441 }
2442 
2443 void LxSpace::CheckArray(scaltype xs[LXSTATIONS][LXLAYERS],
2444  scaltype ys[LXSTATIONS][LXLAYERS],
2445  scaltype zs[LXSTATIONS][LXLAYERS],
2446  scaltype xDisp2Limits[LXSTATIONS],
2447  scaltype yDisp2Limits[LXSTATIONS],
2448  scaltype tx2Limits[LXSTATIONS],
2449  scaltype ty2Limits[LXSTATIONS],
2450  scaltype txBreak2Limits[LXSTATIONS],
2451  scaltype tyBreak2Limits[LXSTATIONS]) {
2452  for (int i = 0; i < LXSTATIONS; ++i) {
2453  scaltype diffZ = zs[i][0] - zs[i][1];
2454  scaltype tx = xs[i][1] / zs[i][1];
2455  scaltype ty = ys[i][1] / zs[i][1];
2456  scaltype dispXL = xs[i][0] - tx * diffZ - xs[i][1];
2457  scaltype dispYL = ys[i][0] - ty * diffZ - ys[i][1];
2458 
2459  if (dispXL < 0) dispXL = -dispXL;
2460 
2461  if (dispYL < 0) dispYL = -dispYL;
2462 
2463  diffZ = zs[i][2] - zs[i][1];
2464  scaltype dispXR = xs[i][2] - tx * diffZ - xs[i][1];
2465  scaltype dispYR = ys[i][2] - ty * diffZ - ys[i][1];
2466 
2467  if (dispXR < 0) dispXR = -dispXR;
2468 
2469  if (dispYR < 0) dispYR = -dispYR;
2470 
2471  scaltype dispX = dispXL < dispXR ? dispXL : dispXR;
2472  scaltype dispY = dispYL < dispYR ? dispYL : dispYR;
2473 
2474  if (stations[i]->disp01XBig - dispX < xDisp2Limits[i])
2475  xDisp2Limits[i] = stations[i]->disp01XBig - dispX;
2476 
2477  if (stations[i]->disp01YBig - dispY < yDisp2Limits[i])
2478  yDisp2Limits[i] = stations[i]->disp01YBig - dispY;
2479 
2480  if (i > 0) {
2481  diffZ = zs[i][1] - zs[i - 1][1];
2482  scaltype tx1 = xs[i][1] / zs[i][1];
2483  tx = (xs[i][1] - xs[i - 1][1]) / diffZ;
2484  scaltype dtx = tx - tx1;
2485 
2486  if (dtx < 0) dtx = -dtx;
2487 
2488  scaltype ty1 = ys[i][1] / zs[i][1];
2489  ty = (ys[i][1] - ys[i - 1][1]) / diffZ;
2490  scaltype dty = ty - ty1;
2491 
2492  if (dty < 0) dty = -dty;
2493 
2494  if (stations[i]->txLimit - dtx < tx2Limits[i])
2495  tx2Limits[i] = stations[i]->txLimit - dtx;
2496 
2497  if (stations[i]->tyLimit - dty < ty2Limits[i])
2498  ty2Limits[i] = stations[i]->tyLimit - dty;
2499 
2500  if (i < LXSTATIONS - 1) {
2501  diffZ = zs[i + 1][1] - zs[i][1];
2502  scaltype tx2 = (xs[i + 1][1] - xs[i][1]) / diffZ;
2503  scaltype ty2 = (ys[i + 1][1] - ys[i][1]) / diffZ;
2504  scaltype txBreak = tx2 - tx;
2505  scaltype tyBreak = ty2 - ty;
2506 
2507  if (txBreak < 0) txBreak = -txBreak;
2508 
2509  if (tyBreak < 0) tyBreak = -tyBreak;
2510 
2511  if (stations[i]->txBreakLimit - txBreak < txBreak2Limits[i])
2512  txBreak2Limits[i] = stations[i]->txBreakLimit - txBreak;
2513 
2514  if (stations[i]->tyBreakLimit - tyBreak < tyBreak2Limits[i])
2515  tyBreak2Limits[i] = stations[i]->tyBreakLimit - tyBreak;
2516  }
2517  }
2518  }
2519 }
2520 
2521 void LxSpace::CheckArray(scaltype xs[LXSTATIONS][LXLAYERS],
2522  scaltype ys[LXSTATIONS][LXLAYERS],
2523  scaltype zs[LXSTATIONS][LXLAYERS],
2524  list<LxPoint*> pts[LXSTATIONS][LXLAYERS],
2525  int level,
2526  scaltype xDisp2Limits[LXSTATIONS],
2527  scaltype yDisp2Limits[LXSTATIONS],
2528  scaltype tx2Limits[LXSTATIONS],
2529  scaltype ty2Limits[LXSTATIONS],
2530  scaltype txBreak2Limits[LXSTATIONS],
2531  scaltype tyBreak2Limits[LXSTATIONS]) {
2532  if (LXSTATIONS * LXLAYERS == level) {
2533  CheckArray(xs,
2534  ys,
2535  zs,
2536  xDisp2Limits,
2537  yDisp2Limits,
2538  tx2Limits,
2539  ty2Limits,
2540  txBreak2Limits,
2541  tyBreak2Limits);
2542  return;
2543  }
2544 
2545  int stNum = level / LXLAYERS;
2546  int lyNum = level % LXLAYERS;
2547  list<LxPoint*>& points = pts[stNum][lyNum];
2548 
2549  for (list<LxPoint*>::iterator i = points.begin(); i != points.end(); ++i) {
2550  LxPoint* point = *i;
2551 
2552  xs[stNum][lyNum] = point->x;
2553  ys[stNum][lyNum] = point->y;
2554  zs[stNum][lyNum] = point->z;
2555 
2556  CheckArray(xs,
2557  ys,
2558  zs,
2559  pts,
2560  level + 1,
2561  xDisp2Limits,
2562  yDisp2Limits,
2563  tx2Limits,
2564  ty2Limits,
2565  txBreak2Limits,
2566  tyBreak2Limits);
2567  }
2568 }
2569 
2570 void LxSpace::CheckArray(ostream& out, LxMCTrack& track) {
2571  scaltype xs[LXSTATIONS][LXLAYERS];
2572  scaltype ys[LXSTATIONS][LXLAYERS];
2573  scaltype zs[LXSTATIONS][LXLAYERS];
2574  list<LxPoint*> pts[LXSTATIONS][LXLAYERS];
2575  int inits[LXSTATIONS][LXLAYERS];
2576  scaltype xDisp2Limits[LXSTATIONS];
2577  scaltype yDisp2Limits[LXSTATIONS];
2578  scaltype tx2Limits[LXSTATIONS];
2579  scaltype ty2Limits[LXSTATIONS];
2580  scaltype txBreak2Limits[LXSTATIONS];
2581  scaltype tyBreak2Limits[LXSTATIONS];
2582  bool busyHits[LXSTATIONS];
2583 
2584  for (int i = 0; i < LXSTATIONS; ++i) {
2585  for (int j = 0; j < LXLAYERS; ++j)
2586  inits[i][j] = 0;
2587 
2588  xDisp2Limits[i] = stations[i]->disp01XBig;
2589  yDisp2Limits[i] = stations[i]->disp01YBig;
2590  tx2Limits[i] = stations[i]->txLimit;
2591  ty2Limits[i] = stations[i]->tyLimit;
2592  txBreak2Limits[i] = stations[i]->txBreakLimit;
2593  tyBreak2Limits[i] = stations[i]->tyBreakLimit;
2594  busyHits[i] = false;
2595  }
2596 
2597  for (vector<LxMCPoint*>::iterator i = track.Points.begin();
2598  i != track.Points.end();
2599  ++i) {
2600  LxMCPoint* point = *i;
2601 
2602  for (list<LxPoint*>::iterator j = point->lxPoints.begin();
2603  j != point->lxPoints.end();
2604  ++j) {
2605  LxPoint* lxPoint = *j;
2606  pts[point->stationNumber][point->layerNumber].push_back(lxPoint);
2607 
2608  if (lxPoint->used) busyHits[point->stationNumber] = true;
2609  }
2610 
2611  inits[point->stationNumber][point->layerNumber] += point->lxPoints.size();
2612  }
2613 
2614  bool ini = true;
2615 
2616  for (int i = 0; i < LXSTATIONS; ++i) {
2617  for (int j = 0; j < LXLAYERS; ++j)
2618  ini = ini && 0 < inits[i][j];
2619  }
2620 
2621  if (!ini) {
2622  out << "CheckArray: track does not contain all the points" << endl;
2623  return;
2624  }
2625 
2626  CheckArray(xs,
2627  ys,
2628  zs,
2629  pts,
2630  0,
2631  xDisp2Limits,
2632  yDisp2Limits,
2633  tx2Limits,
2634  ty2Limits,
2635  txBreak2Limits,
2636  tyBreak2Limits);
2637 
2638  for (int i = 0; i < LXSTATIONS; ++i) {
2639  out << "CheckArray on station " << i << ": ";
2640  out << "dispX to limit: " << xDisp2Limits[i];
2641  out << "; dispY to limit: " << yDisp2Limits[i];
2642 
2643  if (i > 0) {
2644  out << "; Tx to limit: " << tx2Limits[i];
2645  out << "; Ty to limit: " << ty2Limits[i];
2646 
2647  if (i < LXSTATIONS - 1) {
2648  out << "; tx break to limit: " << txBreak2Limits[i];
2649  out << "; ty break to limit: " << tyBreak2Limits[i];
2650  }
2651  }
2652 
2653  if (busyHits[i]) out << "; have busy hits";
2654 
2655  out << endl;
2656  }
2657 
2658  out << endl;
2659 }
y_disp_right_limits_sq
scaltype y_disp_right_limits_sq[LXSTATIONS]
Definition: LxCA.cxx:32
x_coords
scaltype x_coords[LXSTATIONS][LXLAYERS][LXMAXPOINTSONSTATION]
Definition: LxCA.cxx:321
CbmLitToolFactory::CreateTrackPropagator
static TrackPropagatorPtr CreateTrackPropagator(const string &name)
Create track propagation tool by name.
Definition: CbmLitToolFactory.cxx:58
LxTrack
Definition: LxCA.h:268
LxMCPoint
Definition: Simple/LxMC.h:16
LxStation::layers
std::vector< LxLayer * > layers
Definition: LxCA.h:190
LxLayer::AddPoint
LxPoint * AddPoint(int hitId, scaltype x, scaltype y, scaltype z, scaltype dx, scaltype dy, scaltype dz, bool isArtificial=false)
Definition: LxCA.h:161
dispersions01XSmall
static scaltype dispersions01XSmall[]
Definition: LxCA.cxx:937
LxTrack::bX
scaltype bX
Definition: LxCA.h:285
LxStation::RestoreMiddlePoints
void RestoreMiddlePoints()
Definition: LxCA.cxx:1067
LxPoint::artificial
bool artificial
Definition: LxCA.h:57
CbmTrack::GetParamLast
const FairTrackParam * GetParamLast() const
Definition: CbmTrack.h:62
fPropagator
TrackPropagatorPtr fPropagator
Definition: CbmGlobalTrackingTof.cxx:19
y_disp_left_limits
scaltype y_disp_left_limits[LXSTATIONS]
Definition: LxCA.cxx:27
LxSpace::stationsInAlgo
Int_t stationsInAlgo
Definition: LxCA.h:333
LxPoint::track
LxTrack * track
Definition: LxCA.h:58
LxTrack::LxTrack
LxTrack(LxTrackCandidate *tc)
Definition: LxCA.cxx:495
Ask
static void Ask()
Definition: LxCA.cxx:590
ty_limits
scaltype ty_limits[LXSTATIONS]
Definition: LxCA.cxx:23
LxSpace::BuildCandidates
void BuildCandidates(int endStNum, LxRay *ray, LxRay **rays, std::list< LxTrackCandidate * > &candidates, scaltype chi2)
dispersions02XBig
static scaltype dispersions02XBig[]
Definition: LxCA.cxx:943
errorTxBreakCoeff
static scaltype errorTxBreakCoeff
Definition: LxCA.cxx:414
muchStsBreakTy
static TH1F * muchStsBreakTy
Definition: Simple/LxTrackAna.cxx:42
LxStation::Clear
void Clear()
Definition: LxCA.cxx:1036
CbmLitTrackParam::GetX
litfloat GetX() const
Definition: CbmLitTrackParam.h:53
scaltype
#define scaltype
Definition: CbmGlobalTrackingDefs.h:17
LxStation::txLimit
scaltype txLimit
Definition: LxCA.h:208
memset
void memset(T *dest, T i, size_t num)
Definition: L1Grid.cxx:21
sqrt
friend F32vec4 sqrt(const F32vec4 &a)
Definition: L1/vectors/P4_F32vec4.h:41
LxRay::LxRay
LxRay(LxPoint *s, LxPoint *e)
Definition: LxCA.cxx:758
LxStation::~LxStation
~LxStation()
Definition: LxCA.cxx:1025
TrackPropagatorPtr
boost::shared_ptr< CbmLitTrackPropagator > TrackPropagatorPtr
Definition: CbmTofPtrTypes.h:23
LxTrack::length
int length
Definition: LxCA.h:280
LxLayer::Clear
void Clear()
Definition: LxCA.cxx:826
errorTyCoeff
#define errorTyCoeff
Definition: LxCA.h:16
CbmLitTrackParam
Data class for track parameters.
Definition: CbmLitTrackParam.h:29
x_errs
scaltype x_errs[LXSTATIONS][LXLAYERS][LXMAXPOINTSONSTATION]
Definition: LxCA.cxx:323
LxSpace::RestoreMiddlePoints
void RestoreMiddlePoints()
Definition: LxCA.cxx:1846
LxSpace::z_coords
scal_coords * z_coords
Definition: LxCA.h:316
LxLayer::~LxLayer
~LxLayer()
Definition: LxCA.cxx:824
errorXcoeff
static scaltype errorXcoeff
Definition: LxCA.cxx:410
tx_vals
scaltype tx_vals[LXSTATIONS - 1][LXMAXPOINTSONSTATION]
Definition: LxCA.cxx:322
LxSpace::tracks
std::list< LxTrack * > tracks
Definition: LxCA.h:326
i
int i
Definition: L1/vectors/P4_F32vec4.h:25
muchStsBreakY
static TH1F * muchStsBreakY
Definition: Simple/LxTrackAna.cxx:40
CbmLitTrackParam::GetY
litfloat GetY() const
Definition: CbmLitTrackParam.h:54
LxSpace::~LxSpace
~LxSpace()
Definition: LxCA.cxx:1824
LxTrackCandidate
Definition: LxCA.cxx:429
LxTrack::rays
LxRay * rays[LXSTATIONS - 1]
Definition: LxCA.h:281
CbmLitToolFactory.h
Tool factory for creation of littrack algorithms.
dispersions01YBig
static scaltype dispersions01YBig[]
Definition: LxCA.cxx:940
LXMAXPOINTSONSTATION
#define LXMAXPOINTSONSTATION
Definition: Simple/LxSettings.h:13
LxStation::disp02XSmall
scaltype disp02XSmall
Definition: LxCA.h:218
_mm_free
#define _mm_free
Definition: LxCA.cxx:319
LxRay::station
LxStation * station
Definition: LxCA.h:120
LxMCTrack::Points
std::vector< LxMCPoint * > Points
Definition: Simple/LxMC.h:31
LxSpace::extTracks
std::list< LxExtTrack > extTracks
Definition: LxCA.h:327
CbmLitTrackParam::GetCovariance
litfloat GetCovariance(int index) const
Definition: CbmLitTrackParam.h:60
LxExtTrack::recoTrack
std::pair< LxTrack *, Double_t > recoTrack
Definition: LxCATriplets.h:253
LxMCPoint::layerNumber
Int_t layerNumber
Definition: Simple/LxMC.h:18
LxSpace::RemoveClones
void RemoveClones()
Definition: LxCA.cxx:2186
LXLAYERS
#define LXLAYERS
Definition: Simple/LxSettings.h:8
BuildRaysGlobal
void BuildRaysGlobal()
Definition: LxCA.cxx:345
LxLayer::PickNearestPoint
LxPoint * PickNearestPoint(scaltype x, scaltype y)
Definition: LxCA.cxx:833
x_disp_left_limits
scaltype x_disp_left_limits[LXSTATIONS]
Definition: LxCA.cxx:25
LxPoint::y
scaltype y
Definition: LxCA.h:53
LxSpace
Definition: LxCA.h:309
LxTrack::chi2
scaltype chi2
Definition: LxCA.h:283
LXMIDDLE
#define LXMIDDLE
Definition: Simple/LxSettings.h:10
LxStation::stationNumber
int stationNumber
Definition: LxCA.h:206
LxSpace::ty_vals
scal_tans * ty_vals
Definition: LxCA.h:314
LxStation::disp01YBig
scaltype disp01YBig
Definition: LxCA.h:217
LxStation::txBreakSigma
scaltype txBreakSigma
Definition: LxCA.h:212
LxStation::tyBreakLimit
scaltype tyBreakLimit
Definition: LxCA.h:211
LxRay::dty
scaltype dty
Definition: LxCA.h:117
y_disp_right_limits
scaltype y_disp_right_limits[LXSTATIONS]
Definition: LxCA.cxx:31
LxPoint
Definition: LxCA.h:52
dispersions02YSmall
static scaltype dispersions02YSmall[]
Definition: LxCA.cxx:944
LxPoint::dx
scaltype dx
Definition: LxCA.h:54
LxTrackCandidate::rays
LxRay * rays[LXSTATIONS - 1]
Definition: LxCA.cxx:430
ty_limits_sq
scaltype ty_limits_sq[LXSTATIONS]
Definition: LxCA.cxx:24
LxStation::txBreakLimit
scaltype txBreakLimit
Definition: LxCA.h:210
LxExtTrack
Definition: LxCA.h:249
LxLayer::HasPointInRange
bool HasPointInRange(scaltype x, scaltype y, scaltype deltaX, scaltype deltaY)
Definition: LxCA.cxx:903
LxMCPoint::lxPoints
std::list< LxPoint * > lxPoints
Definition: Simple/LxMC.h:19
x_disp_right_limits_sq
scaltype x_disp_right_limits_sq[LXSTATIONS]
Definition: LxCA.cxx:30
LxSpace::InitGlobalCAArrays
void InitGlobalCAArrays()
Definition: LxCA.cxx:44
LxPoint::used
bool used
Definition: LxCA.h:55
_mm_malloc
#define _mm_malloc(X, Y)
Definition: LxCA.cxx:318
LxSpace::muchStsBreakY
scaltype muchStsBreakY
Definition: LxCA.h:329
LxMCPoint::stationNumber
Int_t stationNumber
Definition: Simple/LxMC.h:18
LxStation::LxStation
LxStation(LxSpace *sp, int stNum)
Definition: LxCA.cxx:985
LxLayer::points
std::list< LxPoint * > points
Definition: LxCA.h:153
LxRay::ty
scaltype ty
Definition: LxCA.h:116
LxRay::neighbours
std::list< LxRay * > neighbours
Definition: LxCA.h:121
dispersions02XSmall
static scaltype dispersions02XSmall[]
Definition: LxCA.cxx:942
LxPoint::rays
std::list< LxRay * > rays
Definition: LxCA.h:59
LxTrack::Rebind
void Rebind()
Definition: LxCATriplets.cxx:602
mcTracks
static vector< vector< QAMCTrack > > mcTracks
Definition: CbmTofHitFinderTBQA.cxx:112
z_coords
scaltype z_coords[LXSTATIONS][LXLAYERS][LXMAXPOINTSONSTATION]
Definition: LxCA.cxx:327
LxRay::tx
scaltype tx
Definition: LxCA.h:116
scal_coords
scaltype scal_coords[LXLAYERS][LXMAXPOINTSONSTATION]
Definition: LxCA.h:306
LxLayer::zCoord
scaltype zCoord
Definition: LxCA.h:156
dispersions01YSmall
static scaltype dispersions01YSmall[]
Definition: LxCA.cxx:939
LxTrack::extLinkChi2
scaltype extLinkChi2
Definition: LxCA.h:273
CALC_LINK_WITH_STS_EFF
#define CALC_LINK_WITH_STS_EFF
Definition: riplet/LxSettings.h:30
LxSpace::y_coords
scal_coords * y_coords
Definition: LxCA.h:313
LxPoint::z
scaltype z
Definition: LxCA.h:53
CbmLitConverterFairTrackParam.h
LxSpace::CheckArray
void CheckArray(scaltype xs[LXSTATIONS][LXLAYERS], scaltype ys[LXSTATIONS][LXLAYERS], scaltype zs[LXSTATIONS][LXLAYERS], scaltype xDisp2Limits[LXSTATIONS], scaltype yDisp2Limits[LXSTATIONS], scaltype tx2Limits[LXSTATIONS], scaltype ty2Limits[LXSTATIONS], scaltype txBreak2Limits[LXSTATIONS], scaltype tyBreak2Limits[LXSTATIONS])
Definition: LxCA.cxx:2443
LxSpace::muchStsBreakTy
scaltype muchStsBreakTy
Definition: LxCA.h:331
LxStation::disp01XBig
scaltype disp01XBig
Definition: LxCA.h:215
LxSpace::stations
std::vector< LxStation * > stations
Definition: LxCA.h:325
kLITERROR
@ kLITERROR
Definition: CbmLitEnums.h:25
LxSpace::CalcTangents
void CalcTangents(int station_number)
Lx.h
LxStation::ConnectNeighbours
void ConnectNeighbours()
Definition: LxCA.cxx:1750
dispersions02YBig
static scaltype dispersions02YBig[]
Definition: LxCA.cxx:945
LxSpace::BuildRaysGlobal
void BuildRaysGlobal()
errorTxCoeffSq
#define errorTxCoeffSq
Definition: LxCA.h:15
LxTrack::extTrackCandidates
std::list< std::pair< LxExtTrack *, Double_t > > extTrackCandidates
Definition: LxCATriplets.h:267
LxSpace::Reconstruct
void Reconstruct()
Definition: LxCA.cxx:2107
LXFIRSTSTATION
#define LXFIRSTSTATION
Definition: Simple/LxSettings.h:11
errorTyCoeffSq
#define errorTyCoeffSq
Definition: LxCA.h:17
LxStation::tyLimit
scaltype tyLimit
Definition: LxCA.h:209
LxTrackCandidate::length
Int_t length
Definition: LxCA.cxx:431
errorTxCoeff
#define errorTxCoeff
Definition: LxCA.h:14
LxTrack::aY
scaltype aY
Definition: LxCA.h:286
CbmTrack::GetParamFirst
const FairTrackParam * GetParamFirst() const
Definition: CbmTrack.h:61
LxLayer::layerNumber
int layerNumber
Definition: LxCA.h:155
errorYcoeff
static scaltype errorYcoeff
Definition: LxCA.cxx:411
errorTyBreakCoeff
static scaltype errorTyBreakCoeff
Definition: LxCA.cxx:415
LxExtTrack::track
CbmStsTrack * track
Definition: LxCA.h:250
use_points
bool use_points[LXSTATIONS - 1][LXMAXPOINTSONSTATION]
Definition: LxCA.cxx:329
dispersions01XBig
static scaltype dispersions01XBig[]
Definition: LxCA.cxx:938
LxSpace::FitTracks
void FitTracks()
Definition: LxCA.cxx:2227
LxRay::end
LxPoint * end
Definition: LxCA.h:119
LxPoint::valid
bool valid
Definition: LxCA.h:56
LxSpace::JoinExtTracks
void JoinExtTracks()
Definition: LxCA.cxx:2232
LxSpace::RefineMiddlePoints
void RefineMiddlePoints()
y_coords
scaltype y_coords[LXSTATIONS][LXLAYERS][LXMAXPOINTSONSTATION]
Definition: LxCA.cxx:324
y_disp_left_limits_sq
scaltype y_disp_left_limits_sq[LXSTATIONS]
Definition: LxCA.cxx:28
x_disp_left_limits_sq
scaltype x_disp_left_limits_sq[LXSTATIONS]
Definition: LxCA.cxx:26
LxSpace::ConnectNeighbours
void ConnectNeighbours()
Definition: LxCA.cxx:2031
x_disp_right_limits
scaltype x_disp_right_limits[LXSTATIONS]
Definition: LxCA.cxx:29
y_errs
scaltype y_errs[LXSTATIONS][LXLAYERS][LXMAXPOINTSONSTATION]
Definition: LxCA.cxx:326
points
TClonesArray * points
Definition: Analyze_matching.h:18
LxStation::zCoord
scaltype zCoord
Definition: LxCA.h:207
LxTrack::points
LxPoint * points[LXSTATIONS *LXLAYERS]
Definition: LxCA.h:282
LxStation
Definition: LxCA.h:189
LxCA.h
LxPoint::~LxPoint
~LxPoint()
Definition: LxCA.cxx:731
tx_limits_sq
scaltype tx_limits_sq[LXSTATIONS]
Definition: LxCA.cxx:22
LxTrack::clone
bool clone
Definition: LxCA.h:292
x
Double_t x
Definition: CbmMvdSensorDigiToHitTask.cxx:68
LxRay::source
LxPoint * source
Definition: LxCA.h:118
LxTrackCandidate::chi2
scaltype chi2
Definition: LxCA.cxx:432
LxTrack::restoredPoints
int restoredPoints
Definition: LxCA.h:288
LxStation::BuildRays
void BuildRays()
Definition: LxCA.cxx:1123
LxPoint::CreateRay
void CreateRay(LxPoint *lPoint, scaltype tx, scaltype ty, scaltype dtx, scaltype dty)
Definition: LxCA.cxx:736
m
__m128 m
Definition: L1/vectors/P4_F32vec4.h:26
y
Double_t y
Definition: CbmMvdSensorDigiToHitTask.cxx:68
tx_limits
scaltype tx_limits[LXSTATIONS]
Definition: LxCA.cxx:21
points_counts
int points_counts[LXSTATIONS][LXLAYERS]
Definition: LxCA.cxx:330
CbmLitTrackParam::GetTy
litfloat GetTy() const
Definition: CbmLitTrackParam.h:57
LxPoint::dz
scaltype dz
Definition: LxCA.h:54
LxStation::space
LxSpace * space
Definition: LxCA.h:205
LxTrack::bY
scaltype bY
Definition: LxCA.h:287
LxMCTrack
Definition: Simple/LxMC.h:27
LxStation::tyBreakSigma
scaltype tyBreakSigma
Definition: LxCA.h:213
LxSpace::muchStsBreakTx
scaltype muchStsBreakTx
Definition: LxCA.h:330
muchStsBreakX
static TH1F * muchStsBreakX
Definition: Simple/LxTrackAna.cxx:39
LxSpace::BuildRays
void BuildRays()
Definition: LxCA.cxx:1851
LxTrack::aX
scaltype aX
Definition: LxCA.h:284
scal_tans
scaltype scal_tans[LXMAXPOINTSONSTATION]
Definition: LxCA.h:307
LxPoint::layer
LxLayer * layer
Definition: LxCA.h:60
LxSpace::x_errs
scal_coords * x_errs
Definition: LxCA.h:312
LxSpace::y_errs
scal_coords * y_errs
Definition: LxCA.h:315
LxSpace::x_coords
scal_coords * x_coords
Definition: LxCA.h:310
CalcTangents
static void CalcTangents(int station_number)
Definition: LxCA.cxx:332
LxSpace::LxSpace
LxSpace()
Definition: LxCA.cxx:1799
point_refs
LxPoint * point_refs[LXSTATIONS][LXLAYERS][LXMAXPOINTSONSTATION]
Definition: LxCA.cxx:328
LxStation::disp02YSmall
scaltype disp02YSmall
Definition: LxCA.h:220
LxSpace::Clear
void Clear()
Definition: LxCA.cxx:1834
LXSTATIONS
#define LXSTATIONS
Definition: Simple/LxSettings.h:9
LxPoint::dy
scaltype dy
Definition: LxCA.h:54
LxLayer
Definition: LxCA.h:152
LxSpace::tx_vals
scal_tans * tx_vals
Definition: LxCA.h:311
CbmLitTrackParam::GetTx
litfloat GetTx() const
Definition: CbmLitTrackParam.h:56
LxLayer::LxLayer
LxLayer(LxStation *st, int lNum)
Definition: LxCA.cxx:821
LxLayer::station
LxStation * station
Definition: LxCA.h:154
LxPoint::x
scaltype x
Definition: LxCA.h:53
LxTrack::Fit
void Fit()
Definition: LxCA.cxx:665
muchStsBreakTx
static TH1F * muchStsBreakTx
Definition: Simple/LxTrackAna.cxx:41
LxRay::dtx
scaltype dtx
Definition: LxCA.h:117
CbmLitConverterFairTrackParam::FairTrackParamToCbmLitTrackParam
static void FairTrackParamToCbmLitTrackParam(const FairTrackParam *par, CbmLitTrackParam *litPar)
Definition: CbmLitConverterFairTrackParam.h:40
ty_vals
scaltype ty_vals[LXSTATIONS][LXMAXPOINTSONSTATION]
Definition: LxCA.cxx:325
operator[]
float & operator[](int i)
Definition: L1/vectors/P4_F32vec4.h:3
LxTrackCandidate::LxTrackCandidate
LxTrackCandidate(LxRay **rr, Int_t len, scaltype c2)
Definition: LxCA.cxx:434
LxRay
Definition: LxCA.h:115
LxSpace::muchStsBreakX
scaltype muchStsBreakX
Definition: LxCA.h:328
LxTrack::externalTrack
LxExtTrack * externalTrack
Definition: LxCA.h:269