BmnRoot 26.10.0-18-g4e3414f9
Software framework of the BM@N experiment
Loading...
Searching...
No Matches
L1CATrackFinder.cxx
Go to the documentation of this file.
1/*
2 *=====================================================
3 *
4 * CBM Level 1 Reconstruction
5 *
6 * Authors: I.Kisel, S.Gorbunov, I.Kulakov, M.Zyzak
7 * Documentation: V.Akishina
8 *
9 * e-mail : i.kulakov@gsi.de
10 *
11 *=====================================================
12 *
13 * Finds tracks using the Cellular Automaton algorithm
14 *
15 */
16
17#include "L1AddMaterial.h"
18#include "L1Algo.h"
19#include "L1Branch.h"
20#include "L1Extrapolation.h"
21#include "L1Filtration.h"
22#include "L1Grid.h"
23#include "L1HitArea.h"
24#include "L1HitPoint.h"
25#include "L1HitsSortHelper.h"
26#include "L1Portion.h"
27#include "L1Timer.h"
28#include "L1Track.h"
29#include "L1TrackPar.h"
30#include "TMath.h"
31#include "TStopwatch.h"
32
33#ifdef DRAW
34#include "utils/L1AlgoDraw.h"
35#endif
36#ifdef PULLS
37#include "utils/L1AlgoPulls.h"
38#endif
39#ifdef TRIP_PERFORMANCE
40#include "utils/L1AlgoEfficiencyPerformance.h"
41#endif
42#ifdef DOUB_PERFORMANCE
43#include "utils/L1AlgoEfficiencyPerformance.h"
44#endif
45
46#ifdef TBB
47#include "L1AlgoTBB.h"
48#endif
49
50#include <algorithm>
51#include <iostream>
52#include <list>
53#include <map>
54#include <stdio.h>
55
56// using namespace std;
57using std::cout;
58using std::endl;
59using std::vector;
60
61/// Prepare the portion of data of left hits of a triplet: all hits except the last and the second last station will be
62/// procesesed in the algorythm, the data is orginesed in order to be used by SIMD
63inline void L1Algo::f10( // input
64 int start_lh,
65 int n1,
66 L1HitPoint* vStsHits_l,
67 // output
68 fvec* u_front,
69 fvec* u_back,
70 fvec* zPos,
71 THitI* hitsl_1,
72 fvec* du_front,
73 fvec* du_back // AZ
74)
75{
76 const int end_lh = start_lh + n1;
77 for (int ilh = start_lh, i1 = 0; ilh < end_lh; ilh++, i1++) {
78 int i1_V = i1 / fvecLen;
79 int i1_4 = i1 % fvecLen;
80 L1HitPoint& hitl = vStsHits_l[ilh];
81
82 hitsl_1[i1] = ilh;
83 u_front[i1_V][i1_4] = hitl.U();
84 u_back[i1_V][i1_4] = hitl.V();
85 zPos[i1_V][i1_4] = hitl.Z();
86 if (du_front)
87 du_front[i1_V][i1_4] = hitl.DU(); // AZ
88 if (du_back)
89 du_back[i1_V][i1_4] = hitl.DV(); // AZ
90 }
91}
92
93/// Get the field approximation. Add the target to parameters estimation. Propagaete to the middle station of a triplet.
94inline void L1Algo::f11( /// input 1st stage of singlet search
95 int istal,
96 int istam, /// indexes of left and middle stations of a triplet
97 int n1_V,
98
99 fvec* u_front,
100 fvec* u_back,
101 fvec* zPos,
102 // output
103 // L1TrackPar *T_1,
104 L1TrackPar* T_1,
105 L1FieldRegion* fld_1,
106 fvec* du_front,
107 fvec* du_back // AZ
108)
109{
110 L1Station& stal = vStations[istal];
111 L1Station& stam = vStations[istam];
112 fvec zstal = stal.z;
113 fvec zstam = stam.z;
114
115 L1FieldRegion fld0; // by left hit, target and "center" (station in-between). Used here for extrapolation to
116 // target and to middle hit
117 L1FieldValue centB, l_B _fvecalignment; // field for singlet creation
118 L1FieldValue m_B _fvecalignment; // field for the next extrapolations
119
120 for (int i1_V = 0; i1_V < n1_V; i1_V++) {
121 L1TrackPar& T = T_1[i1_V];
122 L1FieldRegion& fld1 =
123 fld_1[i1_V]; // by middle hit, left hit and "center" . Will be used for extrapolation to right hit
124
125 // get the magnetic field along the track.
126 // (suppose track is straight line with origin in the target)
127 fvec& u = u_front[i1_V];
128 fvec& v = u_back[i1_V];
129 fvec xl, yl; // left(1-st) hit coor
130 fvec zl = zPos[i1_V];
131
132 fvec dzli = 1. / (zl - targZ);
133 StripsToCoor(u, v, xl, yl, stal);
134
135 fvec tx = (xl - targX) * dzli;
136 fvec ty = (yl - targY) * dzli;
137
138 // estimate field for singlet creation
139 int istac = istal / 2; // "center" station
140 // if (istal >= NMvdStations) // to make it in the same way for both with and w\o mvd
141 // istac = (istal-NMvdStations)/2+NMvdStations;
142 L1Station& stac = vStations[istac];
143 fvec zstac = stac.z;
144 stac.fieldSlice.GetFieldValue(targX + tx * (zstac - targZ), targY + ty * (zstac - targZ), centB);
145 stal.fieldSlice.GetFieldValue(targX + tx * (zstal - targZ), targY + ty * (zstal - targZ), l_B);
146 if (istac != istal) {
147 fld0.Set(l_B, stal.z, centB, stac.z, targB, targZ);
148 } else {
149 fld0.Set(l_B, zstal, targB, targZ);
150 }
151 // estimate field for the next extrapolations
152 stam.fieldSlice.GetFieldValue(targX + tx * (zstam - targZ), targY + ty * (zstam - targZ), m_B);
153#define USE_3HITS
154#ifndef USE_3HITS
155 if (istac != istal) {
156 fld1.Set(m_B, zstam, l_B, zstal, centB, zstac);
157 } else {
158 fld1.Set(m_B, zstam, l_B, zstal, targB, targZ);
159 }
160#else // if USE_3HITS // the best now
161 L1FieldValue r_B _fvecalignment;
162 L1Station& star = vStations[istam + 1];
163 fvec zstar = star.z;
164 star.fieldSlice.GetFieldValue(targX + tx * (zstar - targZ), targY + ty * (zstar - targZ), r_B);
165 fld1.Set(r_B, zstar, m_B, zstam, l_B, zstal);
166#endif // USE_3HITS
167
168 T.chi2 = 0.;
169 T.NDF = 2.;
170 if ((isec == kAllSecIter) || (isec == kAllSecEIter) || (isec == kAllSecJumpIter))
171 T.NDF = 0;
172 T.tx = tx;
173 T.ty = ty;
174 T.qp = 0.;
175 T.C20 = T.C21 = 0;
176 T.C30 = T.C31 = T.C32 = 0;
177 T.C40 = T.C41 = T.C42 = T.C43 = 0;
178 T.C22 = T.C33 = MaxSlope * MaxSlope / 9.;
179 T.C44 = MaxInvMom / 3. * MaxInvMom / 3.;
180
181// #define BEGIN_FROM_TARGET
182#ifndef BEGIN_FROM_TARGET // the best now
183
184 T.x = xl;
185 T.y = yl;
186 T.z = zl;
187 T.C00 = stal.XYInfo.C00;
188 T.C10 = stal.XYInfo.C10;
189 T.C11 = stal.XYInfo.C11;
190 // AZ - compute track covariance from hits
191 if (du_front) {
192 fvec& du = du_front[i1_V];
193 fvec& dv = du_back[i1_V];
194 //*
195 fvec c_f = stal.frontInfo.cos_phi;
196 fvec s_f = stal.frontInfo.sin_phi;
197 fvec c_b = stal.backInfo.cos_phi;
198 fvec s_b = stal.backInfo.sin_phi;
199 fvec det = c_f * s_b - s_f * c_b;
200 det *= det;
201 // cout << " 1: " << T.C00 << " " << T.C10 << " " << T.C11 << endl;
202 // T.C00 = ( s_b*s_b*stal.frontInfo.sigma2 + s_f*s_f*stal.backInfo.sigma2 )/det;
203 // T.C10 =-( s_b*c_b*stal.frontInfo.sigma2 + s_f*c_f*stal.backInfo.sigma2 )/det;
204 // T.C11 = ( c_b*c_b*stal.frontInfo.sigma2 + c_f*c_f*stal.backInfo.sigma2 )/det;
205 // cout << " 2: " << T.C00 << " " << T.C10 << " " << T.C11 << endl;
206 T.C00 = (s_b * s_b * du * du + s_f * s_f * dv * dv) / det;
207 T.C10 = -(s_b * c_b * du * du + s_f * c_f * dv * dv) / det;
208 T.C11 = (c_b * c_b * du * du + c_f * c_f * dv * dv) / det;
209 //*/
210 }
211 // AZ
212
213 // add the target
214 {
215 fvec eX, eY, J04, J14;
216 fvec dz;
217 dz = targZ - zl;
218 L1ExtrapolateJXY0(T.tx, T.ty, dz, fld0, eX, eY, J04, J14);
219 fvec J[6];
220 J[0] = dz;
221 J[1] = 0;
222 J[2] = J04;
223 J[3] = 0;
224 J[4] = dz;
225 J[5] = J14;
226 L1FilterVtx(T, targX, targY, TargetXYInfo, eX, eY, J);
227 }
228#else // TODO: doesn't work. Why?
229
230 T.x = targX;
231 T.y = targY;
232
233 T.z = targZ;
234 T.C00 = TargetXYInfo.C00;
235 T.C10 = TargetXYInfo.C10;
236 T.C11 = TargetXYInfo.C11;
237
238 // extrapolate to left hit
239 L1Extrapolate0(T, zl, fld0);
240 for (int ista = 0; ista <= istal - 1; ista++) {
241 if ((isec != kAllPrimEIter) && (isec != kAllSecEIter)) {
242#ifdef USE_RL_TABLE
243 L1AddMaterial(T, fRadThick[ista].GetRadThick(T.x, T.y), MaxInvMom);
244#else
245 L1AddMaterial(T, vStations[ista].materialInfo, MaxInvMom);
246#endif
247 if (ista == NMvdStations - 1)
248 L1AddPipeMaterial(T, MaxInvMom);
249 } else {
250#ifdef USE_RL_TABLE
251 L1AddMaterial(T, fRadThick[ista].GetRadThick(T.x, T.y), MaxInvMom, 1, 0.000511f * 0.000511f);
252#else
253 L1AddMaterial(T, vStations[ista].materialInfo, MaxInvMom, 1, 0.000511f * 0.000511f);
254#endif
255 if (ista == NMvdStations - 1)
256 L1AddPipeMaterial(T, MaxInvMom, 1, 0.000511f * 0.000511f);
257 }
258 }
259 // add left hit
260 L1Filter(T, stal.frontInfo, u);
261 L1Filter(T, stal.backInfo, v);
262#endif
263
264 if ((isec != kAllPrimEIter) && (isec != kAllSecEIter)) {
265#ifdef USE_RL_TABLE
266 L1AddMaterial(T, fRadThick[istal].GetRadThick(T.x, T.y), MaxInvMom);
267#else
268 L1AddMaterial(T, stal.materialInfo, MaxInvMom);
269#endif
270 if ((istam >= NMvdStations) && (istal <= NMvdStations - 1))
271 L1AddPipeMaterial(T, MaxInvMom);
272 } else {
273#ifdef USE_RL_TABLE
274 L1AddMaterial(T, fRadThick[istal].GetRadThick(T.x, T.y), MaxInvMom, 1, 0.000511f * 0.000511f);
275#else
276 L1AddMaterial(T, stal.materialInfo, MaxInvMom, 1, 0.000511f * 0.000511f);
277#endif
278 if ((istam >= NMvdStations) && (istal <= NMvdStations - 1))
279 L1AddPipeMaterial(T, MaxInvMom, 1, 0.000511f * 0.000511f);
280 }
281 L1Extrapolate0(T, zstam, fld0); // TODO: fld1 doesn't work!
282 // L1Extrapolate( T, zstam, T.qp, fld1 );
283
284 } // i1_V
285}
286
287/// Find the doublets. Reformat data in the portion of doublets.
288inline void L1Algo::f20( // input
289 int n1,
290 L1Station& stam,
291 L1HitPoint* vStsHits_l,
292 L1HitPoint* vStsHits_m,
293 int NHits_m,
294 L1TrackPar* T_1,
295 THitI* hitsl_1,
296 // output
297 int& n2,
298 vector<THitI>& i1_2,
299 int& start_mhit,
300#ifdef DOUB_PERFORMANCE
301 vector<THitI>& hitsl_2,
302#endif // DOUB_PERFORMANCE
303 vector<THitI>& hitsm_2,
304 vector<bool>& lmDuplets)
305{
306 n2 = 0; // number of doublet
307 for (int i1 = 0; i1 < n1; i1++) { // for each singlet
308 const int i1_V = i1 / fvecLen;
309 const int i1_4 = i1 % fvecLen;
310 const L1TrackPar& T1 = T_1[i1_V];
311
312 const int n2Saved = n2;
313
314 const THitI nl = vStsHits_l[hitsl_1[i1]].N();
315 fvec Pick_m22 =
316 (DOUBLET_CHI2_CUT
317 - T1.chi2); // if make it bigger the found hits will be rejected later because of the chi2 cut.
318
319 // -- collect possible doublets --
320
321 const fscal iz = 1 / T1.z[i1_4];
322 L1HitArea area(vGrid[&stam - vStations], T1.x[i1_4] * iz, T1.y[i1_4] * iz,
323 (sqrt(Pick_m22 * (T1.C00 + stam.XYInfo.C00)) + MaxDZ * fabs(T1.tx))[i1_4] * iz,
324 (sqrt(Pick_m22 * (T1.C11 + stam.XYInfo.C11)) + MaxDZ * fabs(T1.ty))[i1_4] * iz);
325 THitI imh = 0;
326 while (area.GetNext(imh)) {
327 const L1HitPoint& hitm = vStsHits_m[imh];
328 // check y-boundaries
329
330 // - check whether hit belong to the window ( track position +\- errors ) -
331 const fscal& zm = hitm.Z();
332
333 // check lower boundary
334 fvec y, C11;
335 L1ExtrapolateYC11Line(T1, zm, y, C11);
336 const fscal dy_est2 = (Pick_m22 * fabs(C11 + stam.XYInfo.C11))[i1_4];
337
338 const fscal dy = hitm.Y() - y[i1_4];
339 const fscal dy2 = dy * dy;
340 if (dy2 > dy_est2 && dy < 0)
341 continue;
342
343 // check upper boundary
344 const unsigned short int nm = hitm.N();
345 if ((nl != nm))
346 continue;
347
348 // check x-boundaries
349 fvec x, C00;
350 L1ExtrapolateXC00Line(T1, zm, x, C00);
351 const fscal dx_est2 = (Pick_m22 * fabs(C00 + stam.XYInfo.C00))[i1_4];
352 const fscal dx = hitm.X() - x[i1_4];
353 if (dx * dx > dx_est2)
354 continue;
355
356 // check chi2
357 fvec C10;
358 L1ExtrapolateC10Line(T1, zm, C10);
359 fvec chi2 = T1.chi2;
360 fvec du2 = hitm.DU(); // AZ
361 du2 *= du2; // AZ
362 // L1FilterChi2XYC00C10C11( stam.frontInfo, x, y, C00, C10, C11, chi2, hitm.U() );
363 L1FilterChi2XYC00C10C11(stam.frontInfo, x, y, C00, C10, C11, chi2, hitm.U(), &du2); // AZ
364#ifdef DO_NOT_SELECT_TRIPLETS
365 if (isec != TRACKS_FROM_TRIPLETS_ITERATION)
366#endif
367 if (chi2[i1_4] > DOUBLET_CHI2_CUT)
368 continue;
369 fvec dv2 = hitm.DV(); // AZ
370 dv2 *= dv2; // AZ
371 // L1FilterChi2 ( stam.backInfo, x, y, C00, C10, C11, chi2, hitm.V() );
372 L1FilterChi2(stam.backInfo, x, y, C00, C10, C11, chi2, hitm.V(), &dv2); // AZ
373#ifdef DO_NOT_SELECT_TRIPLETS
374 if (isec != TRACKS_FROM_TRIPLETS_ITERATION)
375#endif
376 if (chi2[i1_4] > DOUBLET_CHI2_CUT)
377 continue;
378
379 i1_2.push_back(i1);
380#ifdef DOUB_PERFORMANCE
381 hitsl_2.push_back(hitsl_1[i1]);
382#endif // DOUB_PERFORMANCE
383 hitsm_2.push_back(imh);
384
385 n2++;
386 }
387
388 lmDuplets[hitsl_1[i1]] = (n2Saved < n2);
389 } // for i1
390 // cout << " singl-doubl " << n1 << " " << n2 << endl; //AZ-200223
391 // for (int jj = 0; jj < hitsm_2.size(); ++jj) cout << hitsm_2[jj] << " "; cout << endl; //AZ-200223
392}
393
394/// Add the middle hits to parameters estimation. Propagate to right station.
395/// Find the triplets(right hit). Reformat data in the portion of triplets.
396inline void L1Algo::f30( // input
397 L1HitPoint* vStsHits_r,
398 L1Station& stam,
399 L1Station& star,
400
401 int istam,
402 int istar,
403 L1HitPoint* vStsHits_m,
404 L1TrackPar* T_1,
405 L1FieldRegion* fld_1,
406 THitI* hitsl_1,
407
408 int n2,
409 vector<THitI>& hitsm_2,
410 vector<THitI>& i1_2,
411
412 const vector<bool>& mrDuplets,
413 // output
414 int& n3,
415 nsL1::vector<L1TrackPar>::TSimd& T_3,
416 vector<THitI>& hitsl_3,
417 vector<THitI>& hitsm_3,
418 vector<THitI>& hitsr_3,
419 nsL1::vector<fvec>::TSimd& u_front_3,
420 nsL1::vector<fvec>::TSimd& u_back_3,
421 nsL1::vector<fvec>::TSimd& z_Pos_3,
422 nsL1::vector<fvec>::TSimd* du_front_3,
423 nsL1::vector<fvec>::TSimd* du_back_3 // AZ
424)
425{
426 THitI hitsl_2[fvecLen];
427 THitI hitsm_2_tmp[fvecLen];
428 fvec fvec_0;
429 // fvec fvec_1 = 1.; //AZ
430 L1TrackPar L1TrackPar_0;
431
432 int n3_V = 0, n3_4 = 0;
433
434 T_3.push_back(L1TrackPar_0);
435 u_front_3.push_back(fvec_0);
436 u_back_3.push_back(fvec_0);
437 z_Pos_3.push_back(fvec_0);
438 // AZ
439 if (du_front_3) {
440 du_front_3->push_back(fvec_0);
441 du_back_3->push_back(fvec_0);
442 // du_front_3->push_back(fvec_1);
443 // du_back_3->push_back(fvec_1);
444 }
445 // AZ
446
447 // ---- Add the middle hits to parameters estimation. Propagate to right station. ----
448 if (istar < NStations) {
449
450 for (int i2 = 0; i2 < n2;) {
451
452 L1TrackPar T2;
453 L1FieldRegion f2;
454 // pack the data
455 fvec u_front_2;
456 fvec u_back_2;
457 fvec zPos_2;
458 fvec du_front_2, du_back_2; // AZ
459 // fvec du_front_2 = 1.0, du_back_2 = 1.0; //AZ
460 int n2_4 = 0;
461 for (; n2_4 < fvecLen && i2 < n2; n2_4++, i2++) {
462
463 if (!mrDuplets[hitsm_2[i2]]) {
464 n2_4--;
465 // cout << " here: " << n2_4 << endl; //AZ-200223
466 continue;
467 }
468
469 const int i1 = i1_2[i2];
470 const int i1_V = i1 / fvecLen;
471 const int i1_4 = i1 % fvecLen;
472
473 const L1TrackPar& T1 = T_1[i1_V];
474 const L1FieldRegion& f1 = fld_1[i1_V];
475 T2.SetOneEntry(n2_4, T1, i1_4);
476 f2.SetOneEntry(n2_4, f1, i1_4);
477
478 const int imh = hitsm_2[i2];
479 const L1HitPoint& hitm = vStsHits_m[imh];
480 u_front_2[n2_4] = hitm.U();
481 u_back_2[n2_4] = hitm.V();
482 zPos_2[n2_4] = hitm.Z();
483 du_front_2[n2_4] = hitm.DU() * hitm.DU(); // AZ
484 du_back_2[n2_4] = hitm.DV() * hitm.DV(); // AZ
485
486 hitsl_2[n2_4] = hitsl_1[i1];
487 hitsm_2_tmp[n2_4] = hitsm_2[i2];
488 } // n2_4
489
490 // add middle hit
491 L1ExtrapolateLine(T2, zPos_2);
492 // L1Filter( T2, stam.frontInfo, u_front_2 );
493 // L1Filter( T2, stam.backInfo, u_back_2 );
494 fvec ww = 1.; // AZ
495 L1Filter(T2, stam.frontInfo, u_front_2, ww, &du_front_2); // AZ
496 L1Filter(T2, stam.backInfo, u_back_2, ww, &du_back_2); // AZ
497
498 if ((isec != kAllPrimEIter) && (isec != kAllSecEIter)) {
499#ifdef USE_RL_TABLE
500 L1AddMaterial(T2, fRadThick[istam].GetRadThick(T2.x, T2.y), T2.qp);
501#else
502 L1AddMaterial(T2, stam.materialInfo, T2.qp);
503#endif
504 if ((istar >= NMvdStations) && (istam <= NMvdStations - 1))
505 L1AddPipeMaterial(T2, T2.qp);
506 } else {
507#ifdef USE_RL_TABLE
508 L1AddMaterial(T2, fRadThick[istam].GetRadThick(T2.x, T2.y), T2.qp, 1, 0.000511f * 0.000511f);
509#else
510 L1AddMaterial(T2, stam.materialInfo, T2.qp, 1, 0.000511f * 0.000511f);
511#endif
512 if ((istar >= NMvdStations) && (istam <= NMvdStations - 1))
513 L1AddPipeMaterial(T2, T2.qp, 1, 0.000511f * 0.000511f);
514 }
515 // extrapolate to the right hit station
516 L1Extrapolate(T2, star.z, T2.qp, f2);
517
518 // ---- Find the triplets(right hit). Reformat data in the portion of triplets. ----
519
520 for (int i2_4 = 0; i2_4 < n2_4; i2_4++) {
521
522 // THitI nm = vStsHits_m[hitsm_2[i2]].n;
523 if (T2.C00[i2_4] < 0 || T2.C11[i2_4] < 0 || T2.C22[i2_4] < 0 || T2.C33[i2_4] < 0 || T2.C44[i2_4] < 0)
524 continue;
525
526 fvec Pick_r22 =
527 (TRIPLET_CHI2_CUT - T2.chi2); // if make it bigger the found hits will be rejected later because
528 // of the chi2 cut. find first possible hit
529
530#ifdef DO_NOT_SELECT_TRIPLETS
531 if (isec == TRACKS_FROM_TRIPLETS_ITERATION)
532 Pick_r22 = Pick_r2 + 1;
533#endif
534
535 const fscal iz = 1 / T2.z[i2_4];
536 L1HitArea area(vGrid[&star - vStations], T2.x[i2_4] * iz, T2.y[i2_4] * iz,
537 (sqrt(Pick_r22 * (T2.C00 + stam.XYInfo.C00)) + MaxDZ * fabs(T2.tx))[i2_4] * iz,
538 (sqrt(Pick_r22 * (T2.C11 + stam.XYInfo.C11)) + MaxDZ * fabs(T2.ty))[i2_4] * iz);
539 THitI irh = 0;
540 while (area.GetNext(irh)) {
541
542 const L1HitPoint& hitr = vStsHits_r[irh];
543
544 const fscal zr = hitr.Z();
545 const fscal yr = hitr.Y();
546
547 // - check whether hit belong to the window ( track position +\- errors ) -
548 // check lower boundary
549 fvec y, C11;
550 L1ExtrapolateYC11Line(T2, zr, y, C11);
551 const fscal dy_est2 =
552 (Pick_r22 * (fabs(C11 + star.XYInfo.C11)))[i2_4]; // TODO for FastPrim dx < dy - other sort is
553 // optimal. But not for doublets
554 const fscal dy = yr - y[i2_4];
555 const fscal dy2 = dy * dy;
556 if (dy2 > dy_est2 && dy < 0)
557 continue; // if (yr < y_minus_new[i2_4]) continue;
558
559 // check upper boundary
560 if (dy2 > dy_est2)
561 continue; // if (yr > y_plus_new [i2_4] ) continue;
562
563 // check x-boundaries
564 fvec x, C00;
565 L1ExtrapolateXC00Line(T2, zr, x, C00);
566 const fscal dx_est2 = (Pick_r22 * (fabs(C00 + star.XYInfo.C00)))[i2_4];
567 const fscal dx = hitr.X() - x[i2_4];
568 if (dx * dx > dx_est2)
569 continue;
570
571 // check chi2 // not effective
572 fvec C10;
573 L1ExtrapolateC10Line(T2, zr, C10);
574 fvec chi2 = T2.chi2;
575 fvec du2 = hitr.DU(), dv2 = hitr.DV(); // AZ
576 du2 *= du2; // AZ
577 dv2 *= dv2; // AZ
578 // L1FilterChi2XYC00C10C11( star.frontInfo, x, y, C00, C10, C11, chi2, hitr.U() );
579 // L1FilterChi2 ( star.backInfo, x, y, C00, C10, C11, chi2, hitr.V() );
580 L1FilterChi2XYC00C10C11(star.frontInfo, x, y, C00, C10, C11, chi2, hitr.U(), &du2); // AZ
581 L1FilterChi2(star.backInfo, x, y, C00, C10, C11, chi2, hitr.V(), &dv2); // AZ
582#ifdef DO_NOT_SELECT_TRIPLETS
583 if (isec != TRACKS_FROM_TRIPLETS_ITERATION)
584#endif
585 if (chi2[i2_4] > TRIPLET_CHI2_CUT || C00[i2_4] < 0 || C11[i2_4] < 0)
586 continue; // chi2_triplet < CHI2_CUT
587
588 // pack triplet
589 L1TrackPar& T3 = T_3[n3_V];
590
591 hitsl_3.push_back(hitsl_2[i2_4]);
592 hitsm_3.push_back(hitsm_2_tmp[i2_4]);
593 hitsr_3.push_back(irh);
594
595 T3.SetOneEntry(n3_4, T2, i2_4);
596 u_front_3[n3_V][n3_4] = hitr.U();
597 u_back_3[n3_V][n3_4] = hitr.V();
598 z_Pos_3[n3_V][n3_4] = zr;
599 // AZ
600 if (du_front_3) {
601 (*du_front_3)[n3_V][n3_4] = hitr.DU() * hitr.DU();
602 (*du_back_3)[n3_V][n3_4] = hitr.DV() * hitr.DV();
603 }
604 // AZ
605
606 n3++;
607 n3_V = n3 / fvecLen;
608 n3_4 = n3 % fvecLen;
609
610 if (!n3_4) {
611 T_3.push_back(L1TrackPar_0);
612 u_front_3.push_back(fvec_0);
613 u_back_3.push_back(fvec_0);
614 z_Pos_3.push_back(fvec_0);
615 // AZ
616 if (du_front_3) {
617 du_front_3->push_back(fvec_0);
618 du_back_3->push_back(fvec_0);
619 // du_front_3->push_back(fvec_1);
620 // du_back_3->push_back(fvec_1);
621 }
622 // AZ
623 }
624 }
625 } // i2_4
626 } // i2_V
627 } // if istar
628}
629
630/// Add the right hits to parameters estimation.
631inline void L1Algo::f31( // input
632 int n3_V,
633 L1Station& star,
634 nsL1::vector<fvec>::TSimd& u_front,
635 nsL1::vector<fvec>::TSimd& u_back,
636 nsL1::vector<fvec>::TSimd& z_Pos,
637 // output
638 // L1TrackPar *T_3
639 nsL1::vector<L1TrackPar>::TSimd& T_3,
640 nsL1::vector<fvec>::TSimd* du_front2,
641 nsL1::vector<fvec>::TSimd* du_back2 // AZ
642)
643{
644 fvec ww = 1.; // AZ
645 for (int i3_V = 0; i3_V < n3_V; i3_V++) {
646 L1ExtrapolateLine(T_3[i3_V], z_Pos[i3_V]);
647 // AZ
648 if (du_front2) {
649 L1Filter(T_3[i3_V], star.frontInfo, u_front[i3_V], ww, &(*du_front2)[i3_V]); // AZ
650 L1Filter(T_3[i3_V], star.backInfo, u_back[i3_V], ww, &(*du_back2)[i3_V]); // AZ
651 } else { // AZ
652 L1Filter(T_3[i3_V], star.frontInfo, u_front[i3_V]); // 2.1/100 sec
653 L1Filter(T_3[i3_V], star.backInfo, u_back[i3_V]); // 2.0/100 sec
654 }
655 };
656}
657
658/// Refit Triplets.
659inline void L1Algo::f32( // input // TODO not updated after gaps introduction
660 int n3,
661 int istal,
662 nsL1::vector<L1TrackPar>::TSimd& T_3,
663 vector<THitI>& hitsl_3,
664 vector<THitI>& hitsm_3,
665 vector<THitI>& hitsr_3,
666 int nIterations)
667{
668 // AZ-220223 This function is not used
669 cout << " !!! L1Algo::f32 should not be called !!! " << endl;
670 exit(9); // AZ-220223
671
672 const int NHits = 3; // triplets
673
674 // prepare data
675 int ista[NHits] = {istal, istal + 1, istal + 2};
676
677 L1Station sta[3];
678 for (int is = 0; is < NHits; is++) {
679 sta[is] = vStations[ista[is]];
680 };
681
682 for (int i3 = 0; i3 < n3; i3++) {
683 int i3_V = i3 / fvecLen;
684 int i3_4 = i3 % fvecLen;
685
686 L1TrackPar& T3 = T_3[i3_V];
687 fscal qp0 = T3.qp[i3_4];
688
689 // prepare data
690 THitI ihit[NHits] = {RealIHit[hitsl_3[i3] + StsHitsUnusedStartIndex[ista[0]]],
691 RealIHit[hitsm_3[i3] + StsHitsUnusedStartIndex[ista[1]]],
692 RealIHit[hitsr_3[i3] + StsHitsUnusedStartIndex[ista[2]]]};
693
694 fscal u[NHits], v[NHits], x[NHits], y[NHits], z[NHits];
695 fvec du2[NHits], dv2[NHits], ww = 1.; // AZ
696 for (int ih = 0; ih < NHits; ih++) {
697 L1StsHit& hit = vStsHits[ihit[ih]];
698 u[ih] = vStsStrips[hit.f];
699 v[ih] = vStsStripsB[hit.b];
700 du2[ih] = vStsStrips[hit.f].errf() * vStsStrips[hit.f].errf(); // AZ
701 dv2[ih] = vStsStripsB[hit.b].errf() * vStsStripsB[hit.b].errf(); // AZ
702 StripsToCoor(u[ih], v[ih], x[ih], y[ih], sta[ih]);
703 z[ih] = vStsZPos[hit.iz];
704 };
705
706 // initialize parameters
707 L1TrackPar T;
708
709 const fvec vINF = .1;
710 T.x = x[0];
711 T.y = y[0];
712
713 fvec dz01 = 1. / (z[1] - z[0]);
714 T.tx = (x[1] - x[0]) * dz01;
715 T.ty = (y[1] - y[0]) * dz01;
716
717 T.qp = qp0;
718 T.z = z[0];
719 T.chi2 = 0.;
720 T.NDF = 2.;
721 T.C00 = sta[0].XYInfo.C00;
722 T.C10 = sta[0].XYInfo.C10;
723 T.C11 = sta[0].XYInfo.C11;
724
725 T.C20 = T.C21 = 0;
726 T.C30 = T.C31 = T.C32 = 0;
727 T.C40 = T.C41 = T.C42 = T.C43 = 0;
728 T.C22 = T.C33 = vINF;
729 T.C44 = 1.;
730
731 // find field along the track
732 L1FieldValue B[3] _fvecalignment;
733 L1FieldRegion fld _fvecalignment;
734
735 fvec tx[3] = {(x[1] - x[0]) / (z[1] - z[0]), (x[2] - x[0]) / (z[2] - z[0]), (x[2] - x[1]) / (z[2] - z[1])};
736 fvec ty[3] = {(y[1] - y[0]) / (z[1] - z[0]), (y[2] - y[0]) / (z[2] - z[0]), (y[2] - y[1]) / (z[2] - z[1])};
737 for (int ih = 0; ih < NHits; ih++) {
738 fvec dz = (sta[ih].z - z[ih]);
739 sta[ih].fieldSlice.GetFieldValue(x[ih] + tx[ih] * dz, y[ih] + ty[ih] * dz, B[ih]);
740 };
741 fld.Set(B[0], sta[0].z, B[1], sta[1].z, B[2], sta[2].z);
742
743 // fit
744 for (int ih = 1; ih < NHits; ih++) {
745 L1Extrapolate(T, z[ih], T.qp, fld);
746
747 if ((isec != kAllPrimEIter) && (isec != kAllSecEIter)) {
748#ifdef USE_RL_TABLE
749 L1AddMaterial(T, fRadThick[ista[ih]].GetRadThick(T.x, T.y), T.qp);
750#else
751 L1AddMaterial(T, sta[ih].materialInfo, T.qp);
752#endif
753
754 if (ista[ih] == NMvdStations - 1)
755 L1AddPipeMaterial(T, T.qp);
756 } else {
757#ifdef USE_RL_TABLE
758 L1AddMaterial(T, fRadThick[ista[ih]].GetRadThick(T.x, T.y), T.qp, 1, 0.000511f * 0.000511f);
759#else
760 L1AddMaterial(T, sta[ih].materialInfo, T.qp, 1, 0.000511f * 0.000511f);
761#endif
762
763 if (ista[ih] == NMvdStations - 1)
764 L1AddPipeMaterial(T, T.qp, 1, 0.000511f * 0.000511f);
765 }
766
767 // L1Filter( T, sta[ih].frontInfo, u[ih] );
768 // L1Filter( T, sta[ih].backInfo, v[ih] );
769 L1Filter(T, sta[ih].frontInfo, u[ih], ww, &du2[ih]); // AZ
770 L1Filter(T, sta[ih].backInfo, v[ih], ww, &dv2[ih]); // AZ
771 }
772
773 // repeat several times in order to improve precision
774 for (int iiter = 0; iiter < nIterations; iiter++) {
775 // fit backward
776 // keep tx,ty,q/p
777 int ih = NHits - 1;
778 T.x = x[ih];
779 T.y = y[ih];
780 T.z = z[ih];
781 T.chi2 = 0.;
782 T.NDF = 2.;
783 T.C00 = sta[ih].XYInfo.C00;
784 T.C10 = sta[ih].XYInfo.C10;
785 T.C11 = sta[ih].XYInfo.C11;
786 // AZ
787 fvec c_f = sta[ih].frontInfo.cos_phi;
788 fvec s_f = sta[ih].frontInfo.sin_phi;
789 fvec c_b = sta[ih].backInfo.cos_phi;
790 fvec s_b = sta[ih].backInfo.sin_phi;
791 fvec det = c_f * s_b - s_f * c_b;
792 det *= det;
793 fvec du22 = du2[ih], dv22 = dv2[ih];
794 // du2 *= du2;
795 // dv2 *= dv2;
796 T.C00 = (s_b * s_b * du22 + s_f * s_f * dv22) / det;
797 T.C10 = -(s_b * c_b * du22 + s_f * c_f * dv22) / det;
798 T.C11 = (c_b * c_b * du22 + c_f * c_f * dv22) / det;
799 // AZ
800
801 T.C20 = T.C21 = 0;
802 T.C30 = T.C31 = T.C32 = 0;
803 T.C40 = T.C41 = T.C42 = T.C43 = 0;
804 T.C22 = T.C33 = vINF;
805 T.C44 = 1.;
806
807 // L1Filter( T, sta[ih].frontInfo, u[ih] );
808 // L1Filter( T, sta[ih].backInfo, v[ih] );
809 for (ih = NHits - 2; ih >= 0; ih--) {
810 L1Extrapolate(T, z[ih], T.qp, fld);
811
812 if ((isec != kAllPrimEIter) && (isec != kAllSecEIter)) {
813#ifdef USE_RL_TABLE
814 L1AddMaterial(T, fRadThick[ista[ih]].GetRadThick(T.x, T.y), T.qp);
815#else
816 L1AddMaterial(T, sta[ih].materialInfo, T.qp);
817#endif
818 if (ista[ih] == NMvdStations)
819 L1AddPipeMaterial(T, T.qp);
820 } else {
821#ifdef USE_RL_TABLE
822 L1AddMaterial(T, fRadThick[ista[ih]].GetRadThick(T.x, T.y), T.qp, 1, 0.000511f * 0.000511f);
823#else
824 L1AddMaterial(T, sta[ih].materialInfo, T.qp, 1, 0.000511f * 0.000511f);
825#endif
826 if (ista[ih] == NMvdStations)
827 L1AddPipeMaterial(T, T.qp, 1, 0.000511f * 0.000511f);
828 }
829
830 // L1Filter( T, sta[ih].frontInfo, u[ih] );
831 // L1Filter( T, sta[ih].backInfo, v[ih] );
832 L1Filter(T, sta[ih].frontInfo, u[ih], ww, &du2[ih]); // AZ
833 L1Filter(T, sta[ih].backInfo, v[ih], ww, &dv2[ih]); // AZ
834 }
835 // fit forward
836 ih = 0;
837 T.x = x[ih];
838 T.y = y[ih];
839 T.z = z[ih];
840 T.chi2 = 0.;
841 T.NDF = 2.;
842 T.C00 = sta[ih].XYInfo.C00;
843 T.C10 = sta[ih].XYInfo.C10;
844 T.C11 = sta[ih].XYInfo.C11;
845 // AZ
846 c_f = sta[ih].frontInfo.cos_phi;
847 s_f = sta[ih].frontInfo.sin_phi;
848 c_b = sta[ih].backInfo.cos_phi;
849 s_b = sta[ih].backInfo.sin_phi;
850 det = c_f * s_b - s_f * c_b;
851 det *= det;
852 du22 = du2[ih];
853 dv22 = dv2[ih];
854 // du2 *= du2;
855 // dv2 *= dv2;
856 T.C00 = (s_b * s_b * du22 + s_f * s_f * dv22) / det;
857 T.C10 = -(s_b * c_b * du22 + s_f * c_f * dv22) / det;
858 T.C11 = (c_b * c_b * du22 + c_f * c_f * dv22) / det;
859 // AZ
860
861 T.C20 = T.C21 = 0;
862 T.C30 = T.C31 = T.C32 = 0;
863 T.C40 = T.C41 = T.C42 = T.C43 = 0;
864 T.C22 = T.C33 = vINF;
865 T.C44 = 1.;
866
867 // L1Filter( T, sta[ih].frontInfo, u[ih] );
868 // L1Filter( T, sta[ih].backInfo, v[ih] );
869 for (ih = 1; ih < NHits; ih++) {
870 L1Extrapolate(T, z[ih], T.qp, fld);
871
872 if ((isec != kAllPrimEIter) && (isec != kAllSecEIter)) {
873#ifdef USE_RL_TABLE
874 L1AddMaterial(T, fRadThick[ista[ih]].GetRadThick(T.x, T.y), T.qp);
875#else
876 L1AddMaterial(T, sta[ih].materialInfo, T.qp);
877#endif
878 if (ista[ih] == NMvdStations + 1)
879 L1AddPipeMaterial(T, T.qp);
880 } else {
881#ifdef USE_RL_TABLE
882 L1AddMaterial(T, fRadThick[ista[ih]].GetRadThick(T.x, T.y), T.qp, 1, 0.000511f * 0.000511f);
883#else
884 L1AddMaterial(T, sta[ih].materialInfo, T.qp, 1, 0.000511f * 0.000511f);
885#endif
886 if (ista[ih] == NMvdStations + 1)
887 L1AddPipeMaterial(T, T.qp, 1, 0.000511f * 0.000511f);
888 }
889
890 // L1Filter( T, sta[ih].frontInfo, u[ih] );
891 // L1Filter( T, sta[ih].backInfo, v[ih] );
892 L1Filter(T, sta[ih].frontInfo, u[ih], ww, &du2[ih]); // AZ
893 L1Filter(T, sta[ih].backInfo, v[ih], ww, &dv2[ih]); // AZ
894 }
895 } // for iiter
896
897 T3.SetOneEntry(i3_4, T, i3_4);
898 } // i3
899} // f32
900
901/// Select triplets. Save them into vTriplets.
902inline void L1Algo::f4( // input
903 int n3,
904 int istal,
905 int istam,
906 int istar,
907 nsL1::vector<L1TrackPar>::TSimd& T_3,
908
909 vector<THitI>& hitsl_3,
910 vector<THitI>& hitsm_3,
911 vector<THitI>& hitsr_3,
912 // output
913 unsigned int& nstaltriplets,
914 std::vector<L1Triplet>& vTriplets_part,
915 unsigned int* TripStartIndexH,
916 unsigned int* TripStopIndexH
917 // #ifdef XXX
918 // ,unsigned int &stat_n_trip
919 // #endif
920)
921{
922
923 for (int i3 = 0; i3 < n3; i3++) {
924 int i3_V = i3 / fvecLen;
925 int i3_4 = i3 % fvecLen;
926 L1TrackPar& T3 = T_3[i3_V];
927
928 // select
929 fscal chi2 = T3.chi2[i3_4];
930#ifdef DO_NOT_SELECT_TRIPLETS
931 if (isec != TRACKS_FROM_TRIPLETS_ITERATION)
932#endif
933 if (!finite(chi2) || chi2 < 0 || chi2 > TRIPLET_CHI2_CUT)
934 continue;
935 // if ( isec == kAllSecIter && T3.qp[i3_4] +sqrt( T3.C44[i3_4] ) < 1./10. ) continue; // why does it cut so much
936 // of ExtraSec ?
937
938 // prepare data
939 fscal MaxInvMomS = MaxInvMom[0];
940 fscal qp = MaxInvMomS + T3.qp[i3_4];
941 if (qp < 0)
942 qp = 0;
943 if (qp > MaxInvMomS * 2)
944 qp = MaxInvMomS * 2;
945 fscal Cqp = 5. * sqrt(fabs(T3.C44[i3_4]));
946
947 fscal scale = 255 / (MaxInvMom[0] * 2);
948
949 qp = (static_cast<unsigned int>(qp * scale)) % 256;
950 Cqp = (static_cast<unsigned int>(Cqp * scale)) % 256;
951 Cqp += 1;
952
953 if (Cqp < 0)
954 Cqp = 0;
955 if (Cqp > 20)
956 Cqp = 20;
957 qp = static_cast<unsigned char>(qp);
958
959 THitI ihitl = hitsl_3[i3] + StsHitsUnusedStartIndex[istal];
960 THitI ihitm = hitsm_3[i3] + StsHitsUnusedStartIndex[istam];
961 THitI ihitr = hitsr_3[i3] + StsHitsUnusedStartIndex[istar];
962 L1_ASSERT(ihitl < StsHitsUnusedStopIndex[istal], ihitl << " < " << StsHitsUnusedStopIndex[istal]);
963 L1_ASSERT(ihitm < StsHitsUnusedStopIndex[istam], ihitm << " < " << StsHitsUnusedStopIndex[istam]);
964 L1_ASSERT(ihitr < StsHitsUnusedStopIndex[istar], ihitr << " < " << StsHitsUnusedStopIndex[istar]);
965
966 // save
967 L1Triplet trip;
968 trip.Set(ihitl, ihitm, ihitr, istal, istam, istar, 0, static_cast<unsigned char>(qp), chi2); // 1.2/1000 sec
969
970 trip.Cqp = static_cast<unsigned char>(Cqp);
971
972 // count statistics
973
974 vTriplets_part.push_back(trip); // 0.67/1000 sec
975
976 if (TripStopIndexH[ihitl] == 0)
977 TripStartIndexH[ihitl] = nstaltriplets;
978 TripStopIndexH[ihitl] = ++nstaltriplets; // < 0.001/100000
979
980 // #ifdef XXX
981 // stat_n_trip++;
982 // #endif
983 } // i3
984} // f4()
985
986/// Find neighbours of triplets. Calculate level of triplets.
987inline void L1Algo::f5( // input
988 // output
989 unsigned int* TripStartIndexH,
990 unsigned int* TripStopIndexH,
991 int* nlevel)
992{
993// cout << "Start Finding levels of triplets " << endl;
994#ifdef TRACKS_FROM_TRIPLETS
995 if (isec != TRACKS_FROM_TRIPLETS_ITERATION)
996#endif
997
998 for (int istal = NStations - 4; istal >= FIRSTCASTATION; istal--) {
999 for (int tripType = 0; tripType < 3; tripType++)
1000 { // tT = 0 - 123triplet, tT = 1 - 124triplet, tT = 2 - 134triplet
1001 if ((((isec != kFastPrimJumpIter) && (isec != kAllPrimJumpIter) && (isec != kAllSecJumpIter))
1002 && (tripType != 0))
1003 || (((isec == kFastPrimJumpIter) || (isec == kAllPrimJumpIter) || (isec == kAllSecJumpIter))
1004 && (tripType != 0) && (istal == NStations - 4)))
1005 continue;
1006
1007 int istam = istal + 1;
1008 int istar = istal + 2;
1009 switch (tripType) {
1010 case 1:
1011 istar++;
1012 break;
1013 case 2:
1014 istam++;
1015 istar++;
1016 break;
1017 }
1018
1019 unsigned int offset_m = TripStartIndex[istam];
1020
1021 for (int it = TripStartIndex[istal]; it < TripStopIndex[istal]; it++) {
1022 L1Triplet* trip = &(vTriplets[it]);
1023 if (istam != trip->GetMSta())
1024 continue;
1025 if (istar != trip->GetRSta())
1026 continue;
1027
1028 unsigned char level = 0;
1029 float chi2 = trip->GetChi2();
1030 unsigned char qp = trip->GetQp();
1031
1032 THitI ihitl = trip->GetLHit();
1033 THitI ihitm = trip->GetMHit();
1034 THitI ihitr = trip->GetRHit();
1035 L1_ASSERT(ihitl < StsHitsUnusedStopIndex[istal], ihitl << " < " << StsHitsUnusedStopIndex[istal]);
1036 L1_ASSERT(ihitm < StsHitsUnusedStopIndex[istam], ihitm << " < " << StsHitsUnusedStopIndex[istam]);
1037 L1_ASSERT(ihitr < StsHitsUnusedStopIndex[istar], ihitr << " < " << StsHitsUnusedStopIndex[istar]);
1038
1039 // neighbours should have 2 common hits
1040 unsigned int first_triplet;
1041 unsigned int last_triplet;
1042 unsigned int iN = TripStartIndexH[ihitm]; // first posible neighbour
1043 unsigned int lastN = TripStopIndexH[ihitm]; // last posible neighbour
1044
1045 // // find first triplet with 2 common hits
1046 // for (; (iN < lastN); ++iN){
1047 // unsigned int jhitm;
1048 // jhitm = vTriplets[iN].GetMHit();
1049 // if (jhitm == ihitr) break;
1050 // };
1051
1052 if (iN < lastN) { // if exist at least one trip
1053 L1_ASSERT(iN >= offset_m, iN << " >= " << offset_m);
1054
1055 first_triplet = iN;
1056
1057 // // find last triplets with 2 common hits
1058 // for (; iN < lastN; ++iN){
1059 // if (vTriplets[iN].GetLHit() != ihitm) break;
1060 // unsigned int jhitm;
1061 // jhitm = vTriplets[iN].GetMHit();
1062 // if (jhitm != ihitr) break;
1063 // };
1064 // last_triplet = iN-1;
1065 last_triplet = lastN;
1066
1067 vector<unsigned int> neighCands; // save neighbour candidates
1068 neighCands.reserve(8); // ~average is 1-2 for central, up to 5
1069 for (iN = first_triplet; iN <= last_triplet; ++iN) {
1070
1071 if (vTriplets[iN].GetMSta() != istar)
1072 continue; // neighbours should have 2 common hits
1073 if (vTriplets[iN].GetMHit() != ihitr)
1074 continue;
1075
1076 L1Triplet* tripn = &vTriplets[iN];
1077
1078 fscal qp2 = tripn->GetQp();
1079 fscal Cqp1 = trip->Cqp;
1080 fscal Cqp2 = tripn->Cqp;
1081 if (fabs(qp - qp2) > PickNeighbour * (Cqp1 + Cqp2))
1082 continue; // neighbours should have same qp
1083
1084 // calculate level
1085 unsigned char jlevel = tripn->GetLevel();
1086 if (level <= jlevel)
1087 level = jlevel + 1;
1088 if (level == jlevel + 1)
1089 neighCands.push_back(iN);
1090 }
1091 for (unsigned int in = 0; in < neighCands.size(); in++) {
1092 const int nLevel = vTriplets[neighCands[in]].GetLevel();
1093 if (level == nLevel + 1)
1094 trip->neighbours.push_back(neighCands[in] - offset_m);
1095 }
1096
1097 }; // if at least one triplet exists
1098
1099 // save information
1100 trip->Set(ihitl, ihitm, ihitr, istal, istam, istar, level, qp, chi2);
1101
1102 nlevel[level]++;
1103
1104 } // vTriplets
1105 } // tripType
1106 } // istal
1107}
1108
1109/// ------------------- doublets on station ----------------------
1110
1111inline void
1112 L1Algo::DupletsStaPort( /// creates duplets: input: @istal - start station number, @istam - last station number,
1113 /// @ip - index of portion, @&n_g1 - number of elements in portion, @*portionStopIndex
1114 int istal,
1115 int istam,
1116 unsigned int ip,
1117 vector<int>& n_g1,
1118 unsigned int* portionStopIndex,
1119
1120 // output
1121 L1TrackPar* T_1, /// @*T_1 - singlets perameters, @*fld_1 - field aproximation, @*hitsl_1- left hits of
1122 /// triplets, @&lmDuplets - existance of a doublet starting from the left hit,
1123 L1FieldRegion* fld_1, /// @&n_2 - number of douplets,@&i1_2 - index of 1st hit in portion indexed by doublet
1124 /// index, @&hitsm_2 - index of middle hit in hits array indexed by doublet index
1125 THitI* hitsl_1,
1126
1127 vector<bool>& lmDuplets,
1128
1129 int& n_2,
1130 vector<THitI>& i1_2,
1131 vector<THitI>& hitsm_2)
1132{
1133 if (istam < NStations) {
1134 L1Station& stam = vStations[istam];
1135
1136 // prepare data
1137 L1HitPoint* vStsHits_l = &((*vStsHitPointsUnused)[0]) + StsHitsUnusedStartIndex[istal];
1138 L1HitPoint* vStsHits_m = &((*vStsHitPointsUnused)[0]) + StsHitsUnusedStartIndex[istam];
1139
1140 int NHits_m = StsHitsUnusedStopIndex[istam] - StsHitsUnusedStartIndex[istam];
1141
1142 int start_mhit = 0; // hit on the middle stantion to start find new doublets
1143
1144 fvec u_front[Portion / fvecLen], u_back[Portion / fvecLen];
1145 fvec du_front[Portion / fvecLen], du_back[Portion / fvecLen]; // AZ
1146 fvec zPos[Portion / fvecLen];
1147
1148 /// prepare the portion of left hits data
1149
1150 int n1 = n_g1[ip];
1151
1152 f10( // input
1153 (ip - portionStopIndex[istal + 1]) * Portion, n1, vStsHits_l,
1154 // output
1155 u_front, u_back, zPos, hitsl_1, du_front, du_back // AZ
1156 );
1157
1158 for (int i = 0; i < n1; i++)
1159 L1_ASSERT(hitsl_1[i] < StsHitsUnusedStopIndex[istal] - StsHitsUnusedStartIndex[istal],
1160 hitsl_1[i] << " < " << StsHitsUnusedStopIndex[istal] - StsHitsUnusedStartIndex[istal]);
1161
1162 int n1_V = (n1 + fvecLen - 1) / fvecLen;
1163
1164 /// Get the field approximation. Add the target to parameters estimation. Propagaete to middle station.
1165
1166 f11(istal, istam, n1_V,
1167
1168 u_front, u_back, zPos,
1169 // output
1170 T_1, fld_1, du_front, du_back // AZ
1171 );
1172
1173 /// Find the doublets. Reformat data in the portion of doublets.
1174
1175#ifdef DOUB_PERFORMANCE
1176 vector<THitI> hitsl_2;
1177#endif // DOUB_PERFORMANCE
1178
1179 f20( // input
1180 n1, stam, vStsHits_l, vStsHits_m, NHits_m, T_1, hitsl_1,
1181 // output
1182 n_2, i1_2, start_mhit,
1183#ifdef DOUB_PERFORMANCE
1184 hitsl_2,
1185#endif // DOUB_PERFORMANCE
1186 hitsm_2, lmDuplets);
1187 // cout << " dub: " << istam << " " << lmDuplets.size() << endl; //AZ-200223
1188 // for (int jj = 0; jj < lmDuplets.size(); ++jj) cout << " " << lmDuplets[jj]; cout << endl; //AZ-200223
1189
1190 for (unsigned int i = 0; i < hitsm_2.size(); i++)
1191 L1_ASSERT(hitsm_2[i] < StsHitsUnusedStopIndex[istam] - StsHitsUnusedStartIndex[istam],
1192 hitsm_2[i] << " " << StsHitsUnusedStopIndex[istam] - StsHitsUnusedStartIndex[istam]);
1193
1194#ifdef DOUB_PERFORMANCE
1195 THitI* RealIHitL = &(RealIHit[StsHitsUnusedStartIndex[istal]]);
1196 THitI* RealIHitM = &(RealIHit[StsHitsUnusedStartIndex[istam]]);
1197 for (int i = 0; i < n2; i++) {
1198 // int i_4 = i%4;
1199 // int i_V = i/4;
1200 THitI iHits[2] = {RealIHitL[hitsl_2[i]], RealIHitM[hitsm_2[i]]};
1201 fL1Eff_doublets->AddOne(iHits);
1202 }
1203#endif // DOUB_PERFORMANCE
1204 }
1205}
1206
1207/// ------------------- Triplets on station ----------------------
1208
1209inline void L1Algo::TripletsStaPort( /// creates triplets: input: @istal - start station number, @istam - middle
1210 /// station number, @istar - last station number, @ip - index of portion, @&n_g1
1211 /// - numer of elements in portion, @*portionStopIndex
1212 int istal,
1213 int istam,
1214 int istar,
1215 ///@nstaltriplets - , @*portionStopIndex, @*T_1 - track parameters for singlets, @*fld_1 - field approximation for
1216 /// singlets, @&n_2 - number of doublets in portion
1217 /// @&n_2 - number of douplets,@&i1_2 - index of 1st hit in portion indexed by doublet index, @&hitsm_2 - index of
1218 /// middle hit in hits array indexed by doublet index
1219 unsigned int& nstaltriplets,
1220 L1TrackPar* T_1,
1221 L1FieldRegion* fld_1,
1222 THitI* hitsl_1,
1223
1224 int& n_2,
1225 unsigned int* portionStopIndex,
1226 vector<THitI>& i1_2,
1227 vector<THitI>& hitsm_2,
1228
1229 const vector<bool>& mrDuplets,
1230 /// output: @*vTriplets_part - array of triplets, @*TripStartIndexH, @*TripStopIndexH - start/stop index of a
1231 /// triplet in the array
1232
1233 // output
1234 std::vector<L1Triplet>* vTriplets_part,
1235 unsigned int* TripStartIndexH,
1236 unsigned int* TripStopIndexH)
1237{
1238
1239 if (istar < NStations) {
1240 // prepare data
1241 L1Station& stam = vStations[istam];
1242 L1Station& star = vStations[istar];
1243
1244 L1HitPoint* vStsHits_m = &((*vStsHitPointsUnused)[0]) + StsHitsUnusedStartIndex[istam];
1245
1246 L1HitPoint* vStsHits_r = 0;
1247 vStsHits_r = &((*vStsHitPointsUnused)[0]) + StsHitsUnusedStartIndex[istar];
1248
1249 int n3 = 0, n3_V;
1250
1251 nsL1::vector<L1TrackPar>::TSimd T_3;
1252
1253 vector<THitI> hitsl_3, hitsm_3, hitsr_3;
1254
1255 /// Add the middle hits to parameters estimation. Propagate to right station.
1256
1257 nsL1::vector<fvec>::TSimd u_front3, u_back3, z_pos3;
1258 nsL1::vector<fvec>::TSimd du_front3, du_back3; // AZ
1259
1260 T_3.reserve(MaxPortionTriplets / fvecLen);
1261 hitsl_3.reserve(MaxPortionTriplets);
1262 hitsm_3.reserve(MaxPortionTriplets);
1263 hitsr_3.reserve(MaxPortionTriplets);
1264 u_front3.reserve(MaxPortionTriplets / fvecLen);
1265 u_back3.reserve(MaxPortionTriplets / fvecLen);
1266 z_pos3.reserve(MaxPortionTriplets / fvecLen);
1267 du_front3.reserve(MaxPortionTriplets / fvecLen); // AZ
1268 du_back3.reserve(MaxPortionTriplets / fvecLen); // AZ
1269
1270 /// Find the triplets(right hit). Reformat data in the portion of triplets.
1271
1272 f30( // input
1273 vStsHits_r, stam, star,
1274
1275 istam, istar, vStsHits_m, T_1, fld_1, hitsl_1,
1276
1277 n_2, hitsm_2, i1_2,
1278
1279 mrDuplets,
1280 // output
1281 n3, T_3, hitsl_3, hitsm_3, hitsr_3, u_front3, u_back3, z_pos3, &du_front3, &du_back3 // AZ
1282 );
1283 // cout << " n3 " << n3 << endl; //AZ-200223
1284 n3_V = (n3 + fvecLen - 1) / fvecLen;
1285
1286 for (unsigned int i = 0; i < hitsl_3.size(); i++)
1287 L1_assert(hitsl_3[i] < StsHitsUnusedStopIndex[istal] - StsHitsUnusedStartIndex[istal]);
1288 for (unsigned int i = 0; i < hitsm_3.size(); i++)
1289 L1_assert(hitsm_3[i] < StsHitsUnusedStopIndex[istam] - StsHitsUnusedStartIndex[istam]);
1290 for (unsigned int i = 0; i < hitsr_3.size(); i++)
1291 L1_assert(hitsr_3[i] < StsHitsUnusedStopIndex[istar] - StsHitsUnusedStartIndex[istar]);
1292
1293 // if (n3 >= MaxPortionTriplets) cout << "isec: " << isec << " stantion: " << istal << " portion number:
1294 // " << ip << " CATrackFinder: Warning: Too many Triplets created in portion" << endl;
1295
1296 /// Add the right hits to parameters estimation.
1297 f31( // input
1298 n3_V, star, u_front3, u_back3, z_pos3,
1299 // output
1300 T_3, &du_front3, &du_back3 // AZ
1301 );
1302
1303 /// refit
1304 // f32( n3, istal, _RealIHit, T_3, hitsl_3, hitsm_3, hitsr_3, 0 );
1305
1306#ifdef TRIP_PERFORMANCE
1307 THitI* RealIHitL = &(RealIHit[StsHitsUnusedStartIndex[istal]]);
1308 THitI* RealIHitM = &(RealIHit[StsHitsUnusedStartIndex[istam]]);
1309 THitI* RealIHitR = &(RealIHit[StsHitsUnusedStartIndex[istar]]);
1310 for (int i = 0; i < n3; i++) {
1311 int i_4 = i % 4;
1312 int i_V = i / 4;
1313 THitI iHits[3] = {RealIHitL[hitsl_3[i]], RealIHitM[hitsm_3[i]], RealIHitR[hitsr_3[i]]};
1314#ifdef PULLS
1315 if (fL1Eff_triplets->AddOne(iHits))
1316 fL1Pulls->AddOne(T_3[i_V], i_4, iHits[2]);
1317#else
1318 fL1Eff_triplets->AddOne(iHits);
1319#endif
1320 if (T_3[i_V].chi2[i_4] < TRIPLET_CHI2_CUT)
1321 fL1Eff_triplets2->AddOne(iHits);
1322 }
1323#endif // TRIP_PERFORMANCE
1324
1325 /// Fill Triplets.
1326
1327 f4( // input
1328 n3, istal, istam, istar, T_3, hitsl_3, hitsm_3, hitsr_3,
1329 // output
1330 nstaltriplets, vTriplets_part[istal], TripStartIndexH, TripStopIndexH);
1331 }
1332}
1333
1334///**********************************************************************************************************************
1335///**********************************************************************************************************************
1336///**********************************************************************************************************************
1337///* *
1338///* *
1339///* *
1340///* CATrackFinder procedure *
1341///* *
1342///* *
1343///* *
1344///**********************************************************************************************************************
1345///**********************************************************************************************************************
1346///**********************************************************************************************************************
1347
1348void L1Algo::CATrackFinder()
1349{
1350
1351 // cout << "Start TrackFinder" << endl;
1352
1353 // cout<<" total STS hits "<<vStsHits.size()<<endl;
1354 // cout<<" total STS hits, strips, back strips "
1355 //<<vStsHits.size()<<" "<<vStsStrips.size()<<" "<<vStsStripsB.size()<<endl;
1356
1357#ifdef PULLS
1358 static L1AlgoPulls* l1Pulls_ = new L1AlgoPulls();
1359 fL1Pulls = l1Pulls_;
1360 fL1Pulls->Init();
1361#endif
1362#ifdef TRIP_PERFORMANCE
1363 static L1AlgoEfficiencyPerformance<3>* l1Eff_triplets_ = new L1AlgoEfficiencyPerformance<3>();
1364 fL1Eff_triplets = l1Eff_triplets_;
1365 fL1Eff_triplets->Init();
1366 static L1AlgoEfficiencyPerformance<3>* l1Eff_triplets2_ = new L1AlgoEfficiencyPerformance<3>();
1367 fL1Eff_triplets2 = l1Eff_triplets2_;
1368 fL1Eff_triplets2->Init();
1369#endif
1370#ifdef DOUB_PERFORMANCE
1371 static L1AlgoEfficiencyPerformance<2>* l1Eff_doublets_ = new L1AlgoEfficiencyPerformance<2>();
1372 fL1Eff_doublets = l1Eff_doublets_;
1373 fL1Eff_doublets->Init();
1374#endif
1375
1376#ifdef DRAW
1377 static L1AlgoDraw draw;
1378 draw.InitL1Draw(this);
1379#endif
1380
1381 cout.precision(6);
1382
1383 TStopwatch c_time; // for performance time
1384#if defined(XXX) || defined(COUNTERS)
1385 static unsigned int stat_N = 0; // number of events
1386 stat_N++;
1387#endif
1388
1389#ifdef XXX
1390 TStopwatch c_timerG; // global
1391 TStopwatch c_timerI; // for iterations
1392
1393 L1CATFIterTimerInfo gti; // global
1394 gti.Add("init ");
1395 gti.Add("iterts");
1396 gti.Add("merge ");
1397
1398 L1CATFTimerInfo ti;
1399 ti.SetNIter(fNFindIterations); // for iterations
1400 ti.Add("init ");
1401 // ti.Add("dblte1");
1402 // ti.Add("dblte2");
1403 ti.Add("tripl1");
1404 // ti.Add("tripl2");
1405 // ti.Add("tripl3");
1406 ti.Add("cpTrls");
1407 ti.Add("nghbrs");
1408 ti.Add("tracks");
1409 ti.Add("finish");
1410
1411 static L1CATFIterTimerInfo stat_gti = gti;
1412 static L1CATFTimerInfo stat_ti = ti;
1413
1414#endif
1415
1416#ifdef COUNTERS
1417 static unsigned int stat_nStartHits = 0;
1418 static unsigned int stat_nHits[fNFindIterations] = {0};
1419 static unsigned int stat_nSinglets[fNFindIterations] = {0};
1420 // static unsigned int stat_nDoublets[fNFindIterations] = {0};
1421 static unsigned int stat_nTriplets[fNFindIterations] = {0};
1422#endif
1423
1424 c_time.Start();
1425
1426#ifdef XXX
1427 c_timerG.Start();
1428#endif
1429
1430 vTracks.clear();
1431 vRecoHits.clear();
1432
1433 // set all hits as usable
1434 for (vector<unsigned char>::iterator is = vSFlag.begin(); is != vSFlag.end(); ++is)
1435 SetFUnUsed(*is);
1436 for (vector<unsigned char>::iterator is = vSFlagB.begin(); is != vSFlagB.end(); ++is)
1437 SetFUnUsed(*is);
1438
1439#ifdef XXX
1440 unsigned int
1441 // stat_max_n_trip = 0,
1442 stat_max_trip_size = 0,
1443 // stat_max_n_dup = 0,
1444 stat_max_BranchSize = 0, stat_max_n_branches = 0;
1445#endif
1446
1447#ifdef DRAW
1448 // draw.DrawInfo();
1449 draw.ClearVeiw();
1450 draw.DrawInputHits();
1451 draw.DrawTarget();
1452 draw.DrawMCTracks();
1453 draw.SaveCanvas("MC_");
1454 draw.DrawAsk();
1455#endif
1456
1457 vector<L1StsHit>
1458 vStsDontUsedHits_A; /// arrays for copy from vStsHits only hits, which aren't used in created tracks
1459 vector<L1StsHit> vStsDontUsedHits_B;
1460
1461 vector<L1HitPoint> vStsDontUsedHitsxy_A; /// arrays for contain prepared information
1462 vector<L1HitPoint> vStsDontUsedHitsxy_B;
1463
1464 vStsHitsUnused = &vStsDontUsedHits_B; /// array of hits uses on current iteration
1465 std::vector<L1StsHit>* vStsHitsUnused_buf = &vStsDontUsedHits_A; /// buffer for copy
1466
1467 vStsHitPointsUnused = &vStsDontUsedHitsxy_B; /// array of info for hits uses on current iteration
1468 std::vector<L1HitPoint>* vStsHitPointsUnused_buf = &vStsDontUsedHitsxy_A;
1469
1470 vStsHitsUnused_buf->clear();
1471 vStsHitPointsUnused_buf->clear();
1472 vector<THitI> RealIHit_v; /// index of hit in vStsHits indexed by index of hit in vStsHitsUnused;
1473 RealIHit_v.resize(vStsHits.size());
1474 RealIHit = &(RealIHit_v[0]);
1475#ifdef AZ
1476 vStsHits2Unus.resize(vStsHits.size()); // AZ-230223
1477 // vAllHitPoints.resize(vStsHits.size()); //AZ-240223
1478#endif
1479
1480 /// fill the arrays for the first iteration
1481 int ihN = 0;
1482 for (int ista = 0; ista < NStations; ista++) {
1483 StsHitsUnusedStartIndex[ista] = ihN;
1484 for (THitI ih = StsHitsStartIndex[ista]; ih < StsHitsStopIndex[ista]; ih++) {
1485 L1StsHit& hit = vStsHits[ih];
1486#ifndef AZ
1487 const L1HitPoint p = CreateHitPoint(hit, ista);
1488 // if (!( int(p.y/p.x*10000)%100 > 75 )) continue;
1489 vStsDontUsedHitsxy_B.push_back(p);
1490#else
1491 vStsHits2Unus[ih] = ihN; // AZ-230223
1492 // vAllHitPoints[ihN] = p; //AZ-240223
1493 const L1HitPoint& p1 = vAllHitPoints[ihN]; // AZ-240223
1494 /*cout << " HP: " << ihN << " " << p.X() << " " << p.Y() << " " << p.Z() << " "
1495 << p1.X() << " " << p1.Y() << " " << p1.Z() << endl; //AZ-240223
1496 cout << vStsStrips[hit.f].f << " " << vStsStripsB[hit.b].f << " " << p1.U() << " "
1497 << p1.V() << endl; //AZ-240223
1498 L1HitPoint p2 (p1.X(), p1.Y(), p1.Z(), p1.V(), p1.U(), p1.DV(), p1.DU(), hit.n); //AZ-240223*/
1499 vStsDontUsedHitsxy_B.push_back(p1); // AZ-240223
1500#endif
1501 RealIHit[ihN++] = ih;
1502 vStsDontUsedHits_B.push_back(hit);
1503 // AZ-240223 vStsDontUsedHitsxy_B.push_back( p );
1504 }
1505 StsHitsUnusedStopIndex[ista] = ihN;
1506 }
1507
1508#ifndef L1_NO_ASSERT
1509 for (int istal = NStations - 1; istal >= 0; istal--) {
1510 L1_ASSERT(StsHitsStopIndex[istal] <= (vStsHits).size(), StsHitsStopIndex[istal] << " <= " << (vStsHits).size());
1511 L1_ASSERT(StsHitsUnusedStopIndex[istal] <= (*vStsHitsUnused).size(),
1512 StsHitsUnusedStopIndex[istal] << " <= " << (*vStsHitsUnused).size());
1513 }
1514#endif // L1_NO_ASSERT
1515
1516 // create Grid
1517 const float hitDensity = sqrt(vStsHitsUnused->size()); /// hitDensity - square root of number of hits, helps to
1518 /// keep the number of hits per bin (sell of grid) the same
1519 float yStep = 0.5 / hitDensity; // empirics. 0.01*sqrt(2374) ~= 0.5
1520 if (yStep > 0.3)
1521 yStep = 0.3;
1522 const float xStep = yStep * 3;
1523
1524 for (int iS = 0; iS < NStations; iS++) {
1525 L1Grid& grid = vGrid[iS];
1526 grid.Create(-1, 1, -0.6, 0.6, xStep, yStep);
1527 }
1528
1529 // sort hits by grid and y/z
1530 L1HitsSortHelper sh(*vStsHitsUnused, *vStsHitPointsUnused, RealIHit_v, vGrid, StsHitsUnusedStartIndex,
1531 StsHitsUnusedStopIndex, NStations);
1532 /// sort hits by grid and y/z
1533 /// input: array of unused hits, unused points, RealIHit_v - array of initial hit indexes, grid,
1534 /// StsHitsUnusedStartIndex, StsHitsUnusedStopIndex - indexes of 1st and last hits on station, NStations - number
1535 /// of stations
1536 sh.Sort();
1537#ifdef AZ
1538 for (int j = 0; j < int(RealIHit_v.size()); ++j)
1539 vStsHits2Unus[RealIHit_v[j]] = j; // AZ-230223
1540#endif
1541
1542#ifdef XXX
1543 c_timerG.Stop();
1544 gti["init "] = c_timerG;
1545 c_timerG.Start(1);
1546#endif
1547
1548#ifdef COUNTERS
1549 cout << " Begin " << endl;
1550 stat_nStartHits += vStsHitsUnused->size();
1551 cout << " NStartHits = " << stat_nStartHits / stat_N << endl;
1552#endif // COUNTERS
1553 // iterations of finding:
1554 // kFastPrimIter = 0, // primary fast track
1555 // kFastPrimJumpIter, // primary fast tracks with gaps. can be dissabled by macro
1556 // kAllPrimIter, // primary all track
1557 // kAllPrimJumpIter, // primary tracks with gaps. can be dissabled by macro
1558 // kAllSecIter // secondary all track
1559 // cout << " AAA " << fNFindIterations << " " << NStations << endl; //AZ-179223
1560
1561 for (isec = 0; isec < fNFindIterations; isec++) { // all finder
1562
1563#ifdef COUNTERS
1564 unsigned int nSinglets = 0;
1565 // unsigned int nDoublets = 0;
1566 // unsigned int nTriplets = 0;
1567#endif
1568
1569#ifdef XXX
1570 // cout << " Begin of iteration " << isec << endl;
1571 TStopwatch c_timer;
1572#endif
1573
1574 for (int iS = 0; iS < NStations; iS++) {
1575 L1Grid& grid = vGrid[iS];
1576 grid.Fill(&((*vStsHitPointsUnused)[StsHitsUnusedStartIndex[iS]]),
1577 StsHitsUnusedStopIndex[iS] - StsHitsUnusedStartIndex[iS]);
1578 }
1579 /// input: 1st hit on station, number of hits on station
1580
1581 // --- SET PARAMETERS FOR THE ITERATION ---
1582
1583 // first station used in the triplets finding
1584 FIRSTCASTATION = 0;
1585 // if ( (isec == kAllSecIter) || (isec == kAllSecEIter) || (isec == kAllSecJumpIter) )
1586 // FIRSTCASTATION = 2;
1587
1588 DOUBLET_CHI2_CUT = 2.706; // prob = 0.1
1589 TRIPLET_CHI2_CUT = 5;
1590 switch (isec) {
1591 case kFastPrimIter:
1592 TRIPLET_CHI2_CUT = 7.815; // prob = 0.05
1593 break;
1594 case kAllPrimIter:
1595 case kAllPrimEIter:
1596 TRIPLET_CHI2_CUT = 7.815; // prob = 0.05
1597 break;
1598 case kAllPrimJumpIter:
1599 TRIPLET_CHI2_CUT = 6.252; // prob = 0.1
1600 break;
1601 case kAllSecIter:
1602 case kAllSecEIter:
1603 TRIPLET_CHI2_CUT = 6.252; // 2.706; // prob = 0.1
1604 break;
1605 }
1606
1607 Pick_gather = 3.0; /// coefficient for size of region for attach new hits to the created track
1608 if ((isec == kAllPrimIter) || (isec == kAllPrimEIter) || (isec == kAllPrimJumpIter) || (isec == kAllSecIter)
1609 || (isec == kAllSecEIter) || (isec == kAllSecJumpIter))
1610 Pick_gather = 4.0;
1611 Pick_gather = 7.0; // AZ
1612
1613 PickNeighbour = 1.0; // (PickNeighbour < dp/dp_error) => triplets are neighbours
1614 // if ( (isec == kFastPrimIter) )
1615 // PickNeighbour = 0.5; // TODO understand why works with 0.2
1616 PickNeighbour = 6.0; // AZ
1617
1618 MaxInvMom = 1.0 / 0.5; // max considered q/p
1619 if ((isec == kAllPrimJumpIter) || (isec == kAllSecIter) || (isec == kAllSecJumpIter))
1620 MaxInvMom = 1.0 / 0.1;
1621 if ((isec == kAllPrimIter) || (isec == kAllPrimEIter) || (isec == kAllSecEIter))
1622 MaxInvMom = 1. / 0.05;
1623
1624 MaxSlope = 1.1;
1625 if ( // (isec == kAllPrimIter) || (isec == kAllPrimEIter) || (isec == kAllPrimJumpIter) ||
1626 (isec == kAllSecIter) || (isec == kAllSecEIter) || (isec == kAllSecJumpIter))
1627 MaxSlope = 1.5;
1628
1629 // define the target
1630 targX = 0;
1631 targY = 0;
1632 targZ = 0; // suppose, what target will be at (0,0,0)
1633
1634 float SigmaTargetX = 0, SigmaTargetY = 0; // target constraint [cm]
1635 if ((isec == kFastPrimIter) || (isec == kFastPrimIter2) || (isec == kFastPrimJumpIter) || (isec == kAllPrimIter)
1636 || (isec == kAllPrimEIter) || (isec == kAllPrimJumpIter))
1637 { // target
1638 targB = vtxFieldValue;
1639 if ((isec == kFastPrimIter) || (isec == kAllPrimIter) || (isec == kAllPrimEIter))
1640 // AZ SigmaTargetX = SigmaTargetY = 0.01; // target
1641 SigmaTargetX = SigmaTargetY = 1.0; // target
1642 else
1643 // AZ SigmaTargetX = SigmaTargetY = 0.1;
1644 SigmaTargetX = SigmaTargetY = 5.0;
1645 }
1646 if ((isec == kAllSecIter) || (isec == kAllSecEIter) || (isec == kAllSecJumpIter))
1647 { // use outer radius of the 1st station as a constraint
1648 L1Station& st = vStations[0];
1649 SigmaTargetX = SigmaTargetY = 10; // st.Rmax[0];
1650 targZ = 0.; //-1.;
1651 st.fieldSlice.GetFieldValue(0, 0, targB);
1652 }
1653 // AZ SigmaTargetX = SigmaTargetY = 1.0; //AZ
1654 targZ = 0.0; // AZ
1655
1656 TargetXYInfo.C00 = SigmaTargetX * SigmaTargetX;
1657 TargetXYInfo.C10 = 0;
1658 TargetXYInfo.C11 = SigmaTargetY * SigmaTargetY;
1659
1660 /// Set correction in order to take into account overlaping and iff z.
1661 /// The reason is that low momentum tracks are too curved and goes not from target direction. That's why sort by
1662 /// hit_y/hit_z is not work idealy If sort by y then it is max diff between same station's modules (~0.4cm)
1663 MaxDZ = 0;
1664 if ((isec == kAllPrimIter) || (isec == kAllPrimEIter) || (isec == kAllPrimJumpIter) || (isec == kAllSecIter)
1665 || (isec == kAllSecEIter) || (isec == kAllSecJumpIter))
1666 MaxDZ = 0.1;
1667 MaxDZ = 2; // AZ
1668
1669 if (NStations > MaxNStations)
1670 cout << " CATrackFinder: Error: Too many Stantions" << endl;
1671
1672 // static unsigned int TripStartIndexH[MaxArrSize], TripStopIndexH[MaxArrSize]; // index region for
1673 // triplets, which begins from hit, indexed by index of lhit.
1674 vector<unsigned int> TripStartIndexH_v, TripStopIndexH_v;
1675 vector<unsigned int> TripStartIndexHG124_v, TripStopIndexHG124_v;
1676 vector<unsigned int> TripStartIndexHG134_v, TripStopIndexHG134_v;
1677
1678 TripStartIndexH_v.resize(vStsHitsUnused->size());
1679 TripStopIndexH_v.resize(vStsHitsUnused->size());
1680 TripStartIndexHG124_v.resize(vStsHitsUnused->size());
1681 TripStopIndexHG124_v.resize(vStsHitsUnused->size());
1682 TripStartIndexHG134_v.resize(vStsHitsUnused->size());
1683 TripStopIndexHG134_v.resize(vStsHitsUnused->size());
1684 unsigned int *TripStartIndexH = &(TripStartIndexH_v[0]), *TripStopIndexH = &(TripStopIndexH_v[0]);
1685 unsigned int *TripStartIndexHG124 = &(TripStartIndexHG124_v[0]),
1686 *TripStopIndexHG124 = &(TripStopIndexHG124_v[0]);
1687 unsigned int *TripStartIndexHG134 = &(TripStartIndexHG134_v[0]),
1688 *TripStopIndexHG134 = &(TripStopIndexHG134_v[0]);
1689 // cout << " BBB " << isec << " " << vStsHitsUnused->size() << endl; //AZ-170223
1690
1691 for (unsigned int i = 0; i < vStsHitsUnused->size(); i++) {
1692 TripStartIndexH[i] = 0; // actually don't need this one, will be initialized at f4(..)
1693 TripStopIndexH[i] = 0;
1694 TripStartIndexHG124[i] = 0;
1695 TripStopIndexHG124[i] = 0;
1696 TripStartIndexHG134[i] = 0;
1697 TripStopIndexHG134[i] = 0;
1698 }
1699
1700 vector<L1Triplet> vTriplets_part[MaxNStations]; // temporary arrays for parallelizing save triplets
1701 vector<L1Triplet> vTripletsG124_part[MaxNStations]; // gaped triplets. type 124
1702 vector<L1Triplet> vTripletsG134_part[MaxNStations]; // gaped triplets. type 134
1703
1704 { // constr of tripletss
1705
1706 // global arrays for keep data of doublets on all stations
1707 vector<int> n_g1;
1708 n_g1.reserve(MaxNPortion);
1709 unsigned int portionStopIndex[MaxNStations]; // number of portion on each station
1710
1711#ifndef L1_NO_ASSERT
1712 for (int istal = NStations - 1; istal >= 0; istal--) {
1713 L1_ASSERT(StsHitsUnusedStopIndex[istal] >= StsHitsUnusedStartIndex[istal],
1714 StsHitsUnusedStopIndex[istal] << " >= " << StsHitsUnusedStartIndex[istal]);
1715 L1_ASSERT(StsHitsUnusedStopIndex[istal] <= (*vStsHitsUnused).size(),
1716 StsHitsUnusedStopIndex[istal] << " <= " << (*vStsHitsUnused).size());
1717 }
1718#endif // L1_NO_ASSERT
1719 // unsigned int nPortions = 0; /// number of portions
1720 {
1721 /// possible left hits of triplets are splited in portions of 16 (4 SIMDs) to use memory faster
1722 portionStopIndex[NStations - 1] = 0;
1723 unsigned int ip = 0; // index of curent portion
1724
1725 for (int istal = NStations - 2; istal >= FIRSTCASTATION; istal--) { // start downstream chambers
1726 int NHits_l = StsHitsUnusedStopIndex[istal] - StsHitsUnusedStartIndex[istal];
1727 // cout << " nhits-sta " << NHits_l << " " << istal << endl; //AZ-190223
1728
1729 int NHits_l_P = NHits_l / Portion;
1730
1731 for (int ipp = 0; ipp < NHits_l_P; ipp++) {
1732 // n_g1[ip++] = Portion;
1733 n_g1.push_back(Portion);
1734#ifdef COUNTERS
1735 nSinglets += Portion;
1736#endif
1737 ip++;
1738 } // start_lh - portion of left hits
1739
1740 // n_g1[ip++] = NHits_l - NHits_l_P*Portion;
1741 n_g1.push_back(NHits_l - NHits_l_P * Portion);
1742#ifdef COUNTERS
1743 nSinglets += n_g1.back();
1744#endif
1745 ip++;
1746 portionStopIndex[istal] = ip;
1747 } // lstations
1748 // nPortions = ip;
1749 }
1750
1751#ifdef COUNTERS
1752 stat_nSinglets[isec] += nSinglets;
1753#endif
1754
1755 /// CREATE doublets
1756#ifdef XXX
1757 c_timer.Stop();
1758 ti[isec]["init "] = c_timer;
1759 c_timer.Start(1);
1760#endif
1761
1762 const unsigned int vPortion = Portion / fvecLen;
1763 L1TrackPar T_1[vPortion];
1764 L1FieldRegion fld_1[vPortion];
1765 THitI hitsl_1[Portion];
1766 L1TrackPar TG_1[vPortion];
1767 L1FieldRegion fldG_1[vPortion];
1768 THitI hitslG_1[Portion];
1769 vector<bool> lmDuplets[MaxNStations]; // is exist a doublet started from indexed by left hit
1770 vector<bool> lmDupletsG[MaxNStations]; // is exist a doublet started from indexed by left hit
1771 for (int istal = NStations - 2; istal >= FIRSTCASTATION; istal--) { // //start downstream chambers
1772 unsigned int nstaltriplets = 0; // find triplets begin from this stantion
1773 unsigned int nstaltripletsG124 = 0; // find triplets begin from this stantion
1774 unsigned int nstaltripletsG134 = 0; // find triplets begin from this stantion
1775
1776 int NHitsSta = StsHitsStopIndex[istal] - StsHitsStartIndex[istal];
1777 lmDuplets[istal].resize(NHitsSta);
1778 lmDupletsG[istal].resize(NHitsSta);
1779 for (unsigned int ip = portionStopIndex[istal + 1]; ip < portionStopIndex[istal]; ip++) {
1780 // nsL1::vector<L1TrackPar>::TSimd T_1; // estimation of parameters
1781 // nsL1::vector<L1FieldRegion>::TSimd fld_1; // magnetic field
1782 // vector<THitI> hitsl_1; // left hits indexed by number of singlets in portion(i2)
1783 vector<THitI> hitsm_2; /// middle hits indexed by number of doublets in portion(i1)
1784 vector<THitI>
1785 i1_2; /// index in portion of singlets(i1) indexed by index in portion of doublets(i2)
1786 int n_2; /// number of doublets in portion
1787
1788 hitsm_2.reserve(MaxPortionDoublets);
1789 i1_2.reserve(MaxPortionDoublets);
1790
1791 DupletsStaPort( // input
1792 istal, istal + 1, ip, n_g1, portionStopIndex,
1793
1794 // output
1795 T_1, fld_1, hitsl_1,
1796
1797 lmDuplets[istal],
1798
1799 n_2, i1_2, hitsm_2);
1800 // cout << " doubl: " << istal << " " << lmDuplets[istal].size() << endl; //AZ-200223
1801 TripletsStaPort( // input
1802 istal, istal + 1, istal + 2, nstaltriplets, T_1, fld_1, hitsl_1,
1803
1804 n_2, portionStopIndex, i1_2, hitsm_2,
1805
1806 lmDuplets[istal + 1],
1807 // output
1808 vTriplets_part, TripStartIndexH, TripStopIndexH);
1809 // cout << " tripl: " << istal << " " << vTriplets_part[istal].size() << endl; //AZ-200223
1810 if ((isec == kFastPrimJumpIter) || (isec == kAllPrimJumpIter) || (isec == kAllSecJumpIter)) {
1811 vector<THitI> hitsmG_2;
1812 vector<THitI> i1G_2;
1813 int nG_2;
1814 hitsmG_2.reserve(MaxPortionDoublets);
1815 i1G_2.reserve(MaxPortionDoublets);
1816
1817 DupletsStaPort( // input
1818 istal, istal + 2, ip, n_g1, portionStopIndex,
1819 // output
1820 TG_1, fldG_1, hitslG_1,
1821
1822 lmDupletsG[istal],
1823
1824 nG_2, i1G_2, hitsmG_2);
1825 TripletsStaPort( // input
1826 istal, istal + 1, istal + 3, nstaltripletsG124, T_1, fld_1, hitsl_1,
1827
1828 n_2, portionStopIndex, i1_2, hitsm_2,
1829
1830 lmDupletsG[istal + 1],
1831 // output
1832 vTripletsG124_part, TripStartIndexHG124, TripStopIndexHG124);
1833 TripletsStaPort( // input
1834 istal, istal + 2, istal + 3, nstaltripletsG134, TG_1, fldG_1, hitslG_1,
1835
1836 nG_2, portionStopIndex, i1G_2, hitsmG_2,
1837
1838 lmDuplets[istal + 2],
1839 // output
1840 vTripletsG134_part, TripStartIndexHG134, TripStopIndexHG134);
1841 }
1842 } // l-stations
1843 }
1844#ifdef XXX
1845 c_timer.Stop();
1846 ti[isec]["tripl1"] = c_timer;
1847#endif
1848 } // constr of triplets
1849 // #ifdef XXX
1850 // c_timer.Stop();
1851 // ti[isec]["tripls"] = c_timer;
1852 // #endif
1853
1854 // #ifdef XXX
1855 // cout<<"isec: " << isec << " n hits, dup, tr_c, trip: "
1856 // <<istat_nhits<<" "<<istat_nduplets<<" "<<istat_n_triplet_c<<" "<<istat_n_triplets<<endl;
1857 // #endif
1858
1859 /// save triplets in vTriplets
1860// #ifdef XXX
1861// cout << " copy triplets in vTriplets " << endl;
1862// #endif
1863#ifdef XXX
1864 c_timer.Start(1);
1865#endif
1866
1867 vTriplets.clear();
1868 int tSize = 0;
1869 for (int istal = NStations - 2; istal >= FIRSTCASTATION; istal--) {
1870 tSize += vTriplets_part[istal].size();
1871 tSize += vTripletsG124_part[istal].size();
1872 tSize += vTripletsG134_part[istal].size();
1873 }
1874 vTriplets.reserve(tSize); // somehow resize is slower...
1875
1876 for (int istal = NStations - 2; istal >= FIRSTCASTATION; istal--) {
1877 TripStartIndex[istal] = vTriplets.size();
1878
1879 for (unsigned int i = 0; i < vTriplets_part[istal].size(); i++) {
1880 vTriplets.push_back(vTriplets_part[istal][i]);
1881 }
1882 for (unsigned int i = 0; i < vTripletsG124_part[istal].size(); i++)
1883 { // it's initialized, so for all iterations it's ok
1884 vTriplets.push_back(vTripletsG124_part[istal][i]); /// save triplets in vTriplets with gaps
1885 }
1886 for (unsigned int i = 0; i < vTripletsG134_part[istal].size(); i++) {
1887 vTriplets.push_back(vTripletsG134_part[istal][i]); /// save triplets in vTriplets with gaps
1888 }
1889
1890 for (THitI i = StsHitsUnusedStartIndex[istal]; i < StsHitsUnusedStopIndex[istal]; i++) {
1891 TripStartIndexH[i] += TripStartIndex[istal];
1892 TripStopIndexH[i] += TripStartIndex[istal] + TripStopIndexHG124[i] - TripStartIndexHG124[i]
1893 + TripStopIndexHG134[i] - TripStartIndexHG134[i];
1894 }
1895
1896 TripStopIndex[istal] = vTriplets.size();
1897 /*cout << " triplets: " << istal << " " << vTriplets_part[istal].size() << " " <<
1898 vTripletsG124_part[istal].size() << " "
1899 << vTripletsG134_part[istal].size() << endl;*/ //AZ-200223
1900 }
1901
1902#ifdef COUNTERS
1903 stat_nTriplets[isec] += vTriplets.size();
1904#endif
1905
1906// #ifdef DRAW
1907// {
1908// draw.ClearVeiw();
1909// // draw.DrawInfo();
1910// draw.DrawRestHits(StsHitsUnusedStartIndex, StsHitsUnusedStopIndex, RealIHit);
1911// draw.DrawTriplets(vTriplets,RealIHit);
1912// TString name = "Trip_";
1913// name += isec;
1914// name += "_";
1915// draw.SaveCanvas(name);
1916// draw.DrawAsk();
1917// }
1918// #endif
1919
1920// #ifdef XXX
1921// { // calculate occupied memory size // TODO
1922
1923// // int sizeD, sizeN, sizeT, sizeF, sizeH, sizeTrip;
1924// // sizeD = sizeof(THitI)*MaxArrSize*MaxNStations*2; // Duplets_*,
1925// // sizeN = sizeof(int)*MaxNPortion*2; // n_g*
1926// // sizeT = sizeof(L1TrackPar)*MaxNPortion*Portion/fvecLen; // T_g1
1927// // sizeF = sizeof(L1FieldRegion)*MaxNPortion*Portion/fvecLen; // fld_g1
1928// // sizeH = sizeof(THitI)*MaxNPortion*(Portion + MaxPortionDoublets*2); // hitsl_g1, hitsm_g2,i1_g2
1929// // sizeTrip = sizeof(L1Triplet)*vTriplets.size();
1930
1931// int sizeD2 = 0, sizeN2, sizeT2, sizeF2, sizeH2, sizeTrip2;
1932// int ndoublets = 0;
1933// for(int i = 0; i < NStations; i++){
1934// sizeD2 += Duplets_hits[i].size();
1935// sizeD2 += Duplets_start[i].size();
1936// ndoublets += Duplets_hits[i].size();
1937// sizeD2*=sizeof(THitI);
1938// }
1939
1940// cout << "number of doublets: " << ndoublets << endl;
1941// cout << "number of triplets: " << vTriplets.size() << endl;
1942// /*
1943// sizeN2 = ( n_g1.size() + n_g2.size() )*sizeof(int);
1944// sizeT2 = T_g1.CalcSize();
1945// sizeF2 = fld_g1.CalcSize();
1946// sizeH2 = ( hitsl_g1.CalcSize() + hitsm_g2.CalcSize() + i1_g2.CalcSize());
1947
1948// cout << sizeD/1000 << " D " << sizeD2/1000 << " [1000elem]" << endl;
1949// cout << sizeN/1000 << " N " << sizeN2/1000 << " [1000elem]" << endl;
1950// cout << sizeT/1000 << " T " << sizeT2/1000 << " [1000elem]" << endl;
1951// cout << sizeF/1000 << " F " << sizeF2/1000 << " [1000elem]" << endl;
1952// cout << sizeH/1000 << " H " << sizeH2/1000 << " [1000elem]" << endl;
1953
1954// long int size = sizeD2 + sizeN2 + sizeT2+ sizeF2 + sizeH2 + sizeTrip;
1955// cout << "Size: " << size/1000 << " KB" << endl;*/
1956// }
1957// #endif
1958#ifdef XXX
1959 c_timer.Stop();
1960 ti[isec]["cpTrls"] = c_timer;
1961#endif
1962
1963 /// Find neighbours of triplets.
1964// #ifdef XXX
1965// cout << " Find neighbours of triplets. " << endl;
1966// #endif
1967#ifdef XXX
1968 c_timer.Start(1);
1969#endif
1970
1971 int nlevels[MaxNStations]; // number of triplets with some number of neighbours.
1972 for (int il = 0; il < NStations; ++il)
1973 nlevels[il] = 0;
1974
1975 f5( // input
1976 // output
1977 TripStartIndexH, TripStopIndexH, nlevels);
1978#ifdef XXX
1979 c_timer.Stop();
1980 ti[isec]["nghbrs"] = c_timer;
1981#endif
1982
1983 ///====================================================================
1984 ///= =
1985 ///= Collect track candidates. CREATE TRACKS =
1986 ///= =
1987 ///====================================================================
1988
1989// #ifdef XXX
1990// cout<<"---- Collect track candidates. ----"<<endl;
1991// #endif
1992#ifdef XXX
1993 c_timer.Start(1);
1994#endif
1995 int min_level = 0; // min level for start triplet. So min track length = min_level+3.
1996 // if (isec == kFastPrimJumpIter) min_level = 1;
1997 if ((isec == kAllSecIter) || (isec == kAllSecEIter) || (isec == kAllSecJumpIter))
1998 min_level = 1; // only the long low momentum tracks
1999#ifdef TRACKS_FROM_TRIPLETS
2000 if (isec == TRACKS_FROM_TRIPLETS_ITERATION)
2001 min_level = 0;
2002#endif
2003 // if (isec == -1) min_level = NStations-3 - 3; //only the longest tracks
2004
2005 // collect consequtive: the longest tracks, shorter, more shorter and so on
2006 for (int ilev = NStations - 3; ilev >= min_level; ilev--) { // choose length
2007 vector<L1Branch> vtrackcandidate;
2008 vtrackcandidate.clear();
2009
2010 // how many levels to check
2011 int nlevel = (NStations - 2) - ilev + 1;
2012 int ntracks = 0;
2013
2014 for (int istaF = FIRSTCASTATION; istaF <= NStations - 3 - ilev; istaF++) { // choose first track-station
2015
2016 // if( ( (isec == kAllSecIter) || (isec == kAllSecEIter) || (isec == kAllSecJumpIter) ) &&
2017 // ( istaF >= NStations-6) ) break; // ghost supression !!! too strong for Lambda. According to
2018 // I.Vassiliev NStations-4 is good
2019
2020 int trip_first = TripStartIndex[istaF];
2021 int trip_end = TripStopIndex[istaF];
2022 for (int itrip = trip_first; itrip < trip_end; itrip++) {
2023 L1Triplet* first_trip = &vTriplets[itrip];
2024
2025 // ghost supression !!!
2026#ifndef FIND_GAPED_TRACKS
2027 if (/*(isec == kFastPrimIter) ||*/ (isec == kAllPrimIter) || (isec == kAllPrimEIter)
2028 || (isec == kAllSecIter) || (isec == kAllSecEIter) || (isec == kAllSecJumpIter))
2029 {
2030#else
2031 if ((isec == kFastPrimIter) || (isec == kFastPrimIter2) || (isec == kFastPrimJumpIter)
2032 || (isec == kAllPrimIter) || (isec == kAllPrimEIter) || (isec == kAllPrimJumpIter)
2033 || (isec == kAllSecIter) || (isec == kAllSecEIter) || (isec == kAllSecJumpIter))
2034 {
2035#endif
2036#ifdef TRACKS_FROM_TRIPLETS
2037 if (isec != TRACKS_FROM_TRIPLETS_ITERATION)
2038#endif
2039 {
2040 if (isec != kFastPrimIter && isec != kAllPrimIter && isec != kAllPrimEIter
2041 && isec != kAllSecEIter)
2042 if (first_trip->GetLevel() == 0)
2043 continue; // ghost suppression // find track with 3 hits only if it was created
2044 // from a chain of triplets, but not from only one triplet
2045 if ((ilev == 0) && (GetFStation(vSFlag[(*vStsHitsUnused)[first_trip->GetLHit()].f]) != 0))
2046 continue; // ghost supression // collect only MAPS tracks-triplets
2047 }
2048 if (first_trip->GetLevel() < ilev)
2049 continue; // try only triplets, which can start track with ilev+3 length. w\o it have more
2050 // ghosts, but efficiency either
2051 }
2052
2053 if (GetFUsed(vSFlag[(*vStsHitsUnused)[first_trip->GetLHit()].f]
2054 | vSFlagB[(*vStsHitsUnused)[first_trip->GetLHit()].b]))
2055 continue;
2056
2057 // OK, search the best track candidate in the branches starting from the first_trip
2058
2059 L1Triplet* curr_trip = first_trip;
2060 L1Branch curr_tr;
2061 curr_tr.StsHits.push_back(RealIHit[first_trip->GetLHit()]);
2062 unsigned char curr_L = 1;
2063 fscal curr_chi2 = first_trip->GetChi2();
2064
2065 L1Branch best_tr = curr_tr;
2066 fscal best_chi2 = curr_chi2;
2067 unsigned char best_L = curr_L;
2068 int NCalls = 1;
2069#ifdef TRACKS_FROM_TRIPLETS
2070 best_tr.StsHits.push_back(RealIHit[first_trip->GetMHit()]);
2071 best_tr.StsHits.push_back(RealIHit[first_trip->GetRHit()]);
2072 best_L += 2;
2073#else
2074 CAFindTrack(istaF, best_tr, best_L, best_chi2, curr_trip, curr_tr, curr_L, curr_chi2, NCalls);
2075 // if (best_L >= 3){
2076 // BranchExtender(best_tr);
2077 // best_L = best_tr.StsHits.size();
2078 // }
2079#endif
2080
2081 // if( ( (isec == kAllSecIter) || (isec == kAllSecEIter) || (isec == kAllSecJumpIter) ) &&
2082 // (GetFStation((*vStsHitsUnused)[best_tr.StsHits[0]].f ) >= 4) ) break; // ghost supression
2083
2084 // {
2085 // if ( best_L < 3 ) continue;
2086 // BranchExtender(best_tr);
2087 // best_L = best_tr.StsHits.size();
2088 // }
2089
2090 if (best_L < ilev + 2)
2091 continue; // lose maximum one hit
2092 if (best_L < min_level + 3)
2093 continue; // should find all hits for min_level
2094
2095 int ndf = best_L * 2 - 5;
2096 best_chi2 = best_chi2 / ndf; // normalize
2097
2098 // if (best_chi2 > TRACK_CHI2_CUT) { // never works because of way neibour triplets and chi2 defined
2099 // !
2100 // cout << " A " << endl;
2101 // continue;
2102 // }
2103
2104 // BranchExtender(best_tr);
2105 // best_L = best_tr.StsHits.size();
2106
2107// if( (isec == kAllPrimIter) || (isec == kAllPrimEIter) || (isec == kAllPrimJumpIter) || (isec == kAllSecIter)
2108// || (isec == kAllSecEIter) || (isec == kAllSecJumpIter) )
2109#ifndef TRACKS_FROM_TRIPLETS
2110 if (fGhostSuppression) {
2111 if (best_L == 3) {
2112 // if( isec == kAllSecIter ) continue; // too /*short*/ secondary track
2113 if (((isec == kAllSecIter) || (isec == kAllSecEIter) || (isec == kAllSecJumpIter))
2114 && (istaF != 0))
2115 continue; // too /*short*/ non-MAPS track
2116 if ((isec != kAllSecIter) && (isec != kAllSecEIter) && (isec != kAllSecJumpIter)
2117 && (best_chi2 > 5.0))
2118 continue;
2119 }
2120 }
2121#endif
2122 // BranchExtender(best_tr);
2123 // best_L = best_tr.StsHits.size();
2124
2125 // store candidate
2126 best_tr.Set(istaF, best_L, best_chi2, first_trip->GetQpOrig(MaxInvMom[0]));
2127 vtrackcandidate.push_back(best_tr);
2128
2129 ntracks++;
2130 } // itrip
2131 // cout<<" level "<<ilev<<", station "<<ista<<" ok"<<endl;
2132 } // istaF
2133
2134 // cout <<"Level "<< ilev <<" track candidates "<< ntracks << endl;
2135 // cout << " total number of track candidates " << vtrackcandidate.size() <<endl;
2136
2137 if (--nlevel == 0)
2138 break;
2139 // -- make competition between tracks of the same length --
2140 // sort track candidates
2141 vector<L1Branch*> vptrackcandidate; // vptrackcandidate - array of pointers to vtrackcandidate
2142 vptrackcandidate.clear();
2143
2144 for (vector<L1Branch>::iterator trIt = vtrackcandidate.begin(); trIt != vtrackcandidate.end(); ++trIt) {
2145
2146 vptrackcandidate.push_back(&(*trIt));
2147 }
2148 // if (isec <= 1)
2149 sort(vptrackcandidate.begin(), vptrackcandidate.end(), L1Branch::comparePChi2);
2150 // else
2151 // sort(vptrackcandidate.begin(), vptrackcandidate.end(), L1Branch::comparePMomentum);
2152
2153 // choose and save tracks
2154 // int nUsed[10] = {0,0,0,0,0, 0,0,0,0,0};
2155 ntracks = 0;
2156 for (vector<L1Branch*>::iterator trIt = vptrackcandidate.begin(); trIt != vptrackcandidate.end(); ++trIt) {
2157 L1Branch* tr = *trIt;
2158#ifdef EXTEND_TRACKS
2159 BranchExtender(*tr);
2160#endif // EXTEND_TRACKS
2161
2162 // check if some hits have been used already
2163 int nused = 0;
2164 for (vector<THitI>::iterator phIt = tr->StsHits.begin(); phIt != tr->StsHits.end(); ++phIt) {
2165 L1StsHit& h = vStsHits[*phIt];
2166 if (GetFUsed(vSFlag[h.f] | vSFlagB[h.b]))
2167 nused++;
2168 }
2169#ifdef TRACKS_FROM_TRIPLETS
2170 if (isec != TRACKS_FROM_TRIPLETS_ITERATION)
2171#endif
2172 if (nused > 0)
2173 continue; // don't allow tracks have shared hits. Track will be shorter, leave it for the next
2174 // iteration
2175
2176 // if ( (isec == kFastPrimIter || isec == kAllPrimIter) && (nused > 2) ||
2177 // (isec != kFastPrimIter) && (isec != kAllPrimIter) && (nused > 0) ){/*nUsed[nused+1]++;*/
2178 // continue;}
2179
2180 //=======================================================
2181 // gather MAPS hits using the Kalman filter
2182 // was here, skipped for a moment
2183 //=======================================================
2184
2185 // if(isec <= kAllSecIter){ // check chi2. gives almost nothing
2186 // L1TrackPar T;
2187 // BranchFitterFast(*tr, T, 0);
2188 // if (T.chi2[0]/T.NDF[0] > TRACK_CHI2_CUT) continue;
2189 // }
2190
2191 // BranchExtender(*tr);
2192
2193 // store track
2194 for (vector<THitI>::iterator phitIt = tr->StsHits.begin(); phitIt != tr->StsHits.end(); ++phitIt) {
2195 L1StsHit& h = vStsHits[*phitIt];
2196#ifdef TRACKS_FROM_TRIPLETS
2197 if (isec != TRACKS_FROM_TRIPLETS_ITERATION)
2198#endif
2199 {
2200 SetFUsed(vSFlag[h.f]);
2201 SetFUsed(vSFlagB[h.b]);
2202 }
2203#ifdef TRACKS_FROM_TRIPLETS
2204 if (isec == TRACKS_FROM_TRIPLETS_ITERATION)
2205#endif
2206 vRecoHits.push_back(*phitIt);
2207 }
2208
2209 L1Track t;
2210 t.NHits = tr->StsHits.size();
2211 t.Momentum = tr->Momentum;
2212 for (int i = 0; i < 6; i++)
2213 t.TFirst[i] = t.TLast[i] = 0; // CHECKME don't need this
2214 for (int i = 0; i < 15; i++)
2215 t.CFirst[i] = t.CLast[i] = 10;
2216 t.chi2 = tr->chi2;
2217 t.NDF = t.NHits * 2 - 5;
2218
2219#ifdef TRACKS_FROM_TRIPLETS
2220 t.NDF = 3; // 3 hits + target
2221 if (isec == kAllSecIter || (isec == kAllSecEIter))
2222 t.NDF = 1; // 3 hits
2223 if (isec == TRACKS_FROM_TRIPLETS_ITERATION)
2224#endif
2225 vTracks.push_back(t);
2226 ntracks++;
2227
2228 } // i_trackCandidate
2229
2230// for (int iu = 0; iu < 10; iu++) cout << iu+1 << " " << nUsed[iu] << endl;
2231#ifdef XXX
2232 unsigned int Bsize = 0, nb = 0;
2233 for (vector<L1Branch*>::iterator trIt = vptrackcandidate.begin(); trIt != vptrackcandidate.end(); ++trIt) {
2234 L1Branch* tr = *trIt;
2235 Bsize += sizeof(L1Branch) + sizeof(THitI) * tr->StsHits.size();
2236 nb++;
2237 }
2238
2239 if (Bsize > stat_max_BranchSize)
2240 stat_max_BranchSize = Bsize;
2241 if (nb > stat_max_n_branches)
2242 stat_max_n_branches = nb;
2243#endif
2244
2245 // cout<< " track candidates stored: "<< ntracks <<endl;
2246 } // ilev
2247
2248#ifdef XXX
2249 c_timer.Stop();
2250 ti[isec]["tracks"] = c_timer;
2251 c_timer.Start(1);
2252#endif
2253 // renew hits array
2254 // #ifdef XXX
2255 // cout << " renew hits arrays " << endl;
2256 // #endif
2257
2258 int StsHitsUnusedStartIndex_temp;
2259 int ista = 0;
2260 int nDontUsedHits = 0;
2261 const int HSize = vStsHitsUnused->size();
2262 vStsHitsUnused_buf->clear();
2263 vStsHitPointsUnused_buf->clear();
2264 vStsHitsUnused_buf->reserve(HSize); // in fact a memory is already reserved after clear
2265 vStsHitPointsUnused_buf->reserve(HSize);
2266 for (; ista < NStations; ista++) {
2267
2268 StsHitsUnusedStartIndex_temp = StsHitsUnusedStartIndex[ista];
2269 StsHitsUnusedStartIndex[ista] = nDontUsedHits;
2270 for (THitI ih = StsHitsUnusedStartIndex_temp; ih < StsHitsUnusedStopIndex[ista]; ih++) {
2271 const L1StsHit& hit = (*vStsHitsUnused)[ih];
2272 if (GetFUsed(vSFlag[hit.f] | vSFlagB[hit.b])) {
2273 continue;
2274 } // if used
2275 vStsHitsUnused_buf->push_back(hit);
2276 vStsHitPointsUnused_buf->push_back((*vStsHitPointsUnused)[ih]);
2277
2278 RealIHit[nDontUsedHits] = RealIHit[ih];
2279#ifdef AZ
2280 vStsHits2Unus[RealIHit[ih]] = nDontUsedHits; // AZ-230223
2281#endif
2282 nDontUsedHits++;
2283 }
2284 StsHitsUnusedStopIndex[ista] = nDontUsedHits;
2285 }
2286 std::vector<L1StsHit>* vStsHitsUnused_temp = vStsHitsUnused;
2287 vStsHitsUnused = vStsHitsUnused_buf;
2288 vStsHitsUnused_buf = vStsHitsUnused_temp;
2289 std::vector<L1HitPoint>* vStsHitsUnused_temp2 = vStsHitPointsUnused;
2290 vStsHitPointsUnused = vStsHitPointsUnused_buf;
2291 vStsHitPointsUnused_buf = vStsHitsUnused_temp2;
2292
2293 vStsHitsUnused_buf->clear();
2294 vStsHitPointsUnused_buf->clear();
2295// cout << " clear hits array end " << endl;
2296
2297// cout<<"End of collecting track candidates "<<endl;
2298#ifdef XXX
2299 c_timer.Stop();
2300 ti[isec]["finish"] = c_timer;
2301#endif
2302#ifdef XXX
2303 // if( stat_max_n_trip<stat_n_trip ) stat_max_n_trip = vTriplets.size();
2304 unsigned int tsize = vTriplets.size() * sizeof(L1Triplet);
2305 if (stat_max_trip_size < tsize)
2306 stat_max_trip_size = tsize;
2307#endif
2308// #ifdef DRAW
2309// draw.ClearVeiw();
2310// // draw.DrawInfo();
2311// draw.DrawRestHits(StsHitsUnusedStartIndex, StsHitsUnusedStopIndex, RealIHit);
2312// draw.DrawRecoTracks();
2313// draw.SaveCanvas("Reco_"+isec+"_");
2314// draw.DrawAsk();
2315// #endif
2316
2317// #ifdef PULLS
2318// fL1Pulls->Build(1);
2319// #endif
2320#ifdef COUNTERS
2321 stat_nHits[isec] += vStsHitsUnused->size();
2322
2323 cout << "iter = " << isec << endl;
2324 cout << " NSinglets = " << stat_nSinglets[isec] / stat_N << endl;
2325 // cout << " NDoublets = " << stat_nDoublets[isec]/stat_N << endl;
2326 cout << " NTriplets = " << stat_nTriplets[isec] / stat_N << endl;
2327 cout << " NHitsUnused = " << stat_nHits[isec] / stat_N << endl;
2328#endif // COUNTERS
2329
2330 } // for (int isec
2331
2332#ifdef XXX
2333 c_timerG.Stop();
2334 gti["iterts"] = c_timerG;
2335 c_timerG.Start(1);
2336#endif
2337
2338#ifndef TRACKS_FROM_TRIPLETS
2339#ifdef MERGE_CLONES
2340 CAMergeClones();
2341#endif
2342#endif
2343
2344#ifdef XXX
2345 c_timerG.Stop();
2346 gti["merge "] = c_timerG;
2347#endif
2348
2349 // cout<< "Total number of tracks found " << vTracks.size() <<endl;
2350
2351 //==================================
2352
2353 c_time.Stop();
2354
2355 // cout << "End TrackFinder" << endl;
2356 // CATime = (double(c_time.CpuTime()));
2357 CATime = (double(c_time.RealTime()));
2358
2359#ifdef XXX
2360
2361 cout << endl << " --- Timers, ms --- " << endl;
2362 ti.Calc();
2363 stat_ti += ti;
2364 L1CATFTimerInfo tmp_ti = stat_ti / 0.001 / stat_N; // ms
2365
2366 tmp_ti.PrintReal();
2367 stat_gti += gti;
2368 L1CATFIterTimerInfo tmp_gti = stat_gti / 0.001 / stat_N; // ms
2369 tmp_gti.PrintReal(1);
2370
2371#if 0
2372 static long int NTimes =0, NHits=0, HitSize =0, NStrips=0, StripSize =0, NStripsB=0, StripSizeB =0,
2373 NDup=0, DupSize=0, NTrip=0, TripSize=0, NBranches=0, BranchSize=0, NTracks=0, TrackSize=0 ;
2374
2375 NTimes++;
2376 NHits += vStsHitsUnused->size();
2377 HitSize += vStsHitsUnused->size()*sizeof(L1StsHit);
2378 NStrips+= vStsStrips.size();
2379 StripSize += vStsStrips.size()*sizeof(fscal) + vSFlag.size()*sizeof(unsigned char);
2380 NStripsB+= vSFlagB.size();
2381 StripSizeB += vStsStripsB.size()*sizeof(fscal) + vSFlagB.size()*sizeof(unsigned char);
2382 NDup += stat_max_n_dup;
2383 DupSize += stat_max_n_dup*sizeof(/*short*/ int);
2384 NTrip += stat_max_n_trip;
2385 TripSize += stat_max_trip_size;
2386
2387 NBranches += stat_max_n_branches;
2388 BranchSize += stat_max_BranchSize;
2389 NTracks += vTracks.size();
2390 TrackSize += sizeof(L1Track)*vTracks.size() + sizeof(THitI)*vRecoHits.size();
2391 int k = 1024*NTimes;
2392
2393 cout<<"L1 Event size: \n"
2394 <<HitSize/k<<"kB for "<<NHits/NTimes<<" hits, \n"
2395 <<StripSize/k<<"kB for "<<NStrips/NTimes<<" strips, \n"
2396 <<StripSizeB/k<<"kB for "<<NStripsB/NTimes<<" stripsB, \n"
2397 <<DupSize/k<<"kB for "<<NDup/NTimes<<" doublets, \n"
2398 <<TripSize/k<<"kB for "<<NTrip/NTimes<<" triplets\n"
2399 <<BranchSize/k<<"kB for "<<NBranches/NTimes<<" branches, \n"
2400 <<TrackSize/k<<"kB for "<<NTracks/NTimes<<" tracks. "<<endl;
2401 cout<<" L1 total event size = "
2402 <<( HitSize + StripSize + StripSizeB + DupSize + TripSize + BranchSize + TrackSize )/k
2403 <<" Kb"<<endl;
2404#endif // 0
2405 // cout << "===> ... CA Track Finder " << endl << endl;
2406#endif
2407
2408#ifdef DRAW
2409 draw.ClearVeiw();
2410 draw.DrawInputHits();
2411 // draw.DrawInfo();
2412 draw.DrawRecoTracks();
2413 draw.SaveCanvas("Reco_");
2414 draw.DrawAsk();
2415#endif
2416#ifdef PULLS
2417 static int iEvee = 0;
2418 iEvee++;
2419 if (iEvee % 1 == 0)
2420 fL1Pulls->Build(1);
2421#endif
2422#ifdef DOUB_PERFORMANCE
2423 fL1Eff_doublets->CalculateEff();
2424 fL1Eff_doublets->Print("Doublets performance.", 1);
2425#endif
2426#ifdef TRIP_PERFORMANCE
2427 fL1Eff_triplets->CalculateEff();
2428 fL1Eff_triplets->Print("Triplet performance", 1);
2429// fL1Eff_triplets2->CalculateEff();
2430// fL1Eff_triplets2->Print("Triplet performance. After chi2 cut");
2431#endif
2432}
2433
2434/** *************************************************************
2435 * *
2436 * The routine performs recursive search for tracks *
2437 * *
2438 * I. Kisel 06.03.05 *
2439 * I.Kulakov 2012 *
2440 * *
2441 ****************************************************************/
2442
2443void L1Algo::CAFindTrack(int ista,
2444 L1Branch& best_tr,
2445 unsigned char& best_L,
2446 fscal& best_chi2,
2447 const L1Triplet* curr_trip,
2448 L1Branch& curr_tr,
2449 unsigned char& curr_L,
2450 fscal& curr_chi2,
2451 int& NCalls)
2452/// recursive search for tracks
2453/// input: @ista - station index, @&best_tr - best track for the privious call, @&best_L -
2454/// output: @&NCalls - number of function calls
2455{
2456 if (curr_trip->GetLevel() == 0) { // the end of the track -> check and store
2457
2458 // -- finish with current track
2459 // add rest of hits
2460 THitI ihitm = curr_trip->GetMHit();
2461 THitI ihitr = curr_trip->GetRHit();
2462 if (!GetFUsed(vSFlag[(*vStsHitsUnused)[ihitm].f] | vSFlagB[(*vStsHitsUnused)[ihitm].b])) {
2463 curr_tr.StsHits.push_back(RealIHit[ihitm]);
2464 curr_L++;
2465 }
2466 if (!GetFUsed(vSFlag[(*vStsHitsUnused)[ihitr].f] | vSFlagB[(*vStsHitsUnused)[ihitr].b])) {
2467 curr_tr.StsHits.push_back(RealIHit[ihitr]);
2468 curr_L++;
2469 }
2470 if (curr_chi2 > TRACK_CHI2_CUT * (curr_L * 2 - 5))
2471 return;
2472
2473 // // try to find more hits
2474 // #ifdef EXTEND_TRACKS
2475 // // if ( isec == kFastPrimJumpIter || isec == kAllPrimJumpIter || isec == kAllSecJumpIter )
2476 // if (curr_L >= 3){
2477 // //curr_chi2 = BranchExtender(curr_tr);
2478 // BranchExtender(curr_tr);
2479 // curr_L = curr_tr.StsHits.size();
2480 // // if( 2*curr_chi2 > (2*(curr_L*2-5) + 1) * 4*4 ) return;
2481 // }
2482 // #endif // EXTEND_TRACKS
2483
2484 // -- select the best
2485 if ((curr_L > best_L) || ((curr_L == best_L) && (curr_chi2 < best_chi2))) {
2486 best_tr = curr_tr;
2487 best_chi2 = curr_chi2;
2488 best_L = curr_L;
2489 }
2490
2491 return;
2492 }
2493
2494 // try to extend. try all possible triplets
2495 int offset = TripStartIndex[curr_trip->GetMSta()];
2496 int NNeighs = curr_trip->neighbours.size();
2497 for (int in = 0; in < NNeighs; in++) {
2498 int index = curr_trip->neighbours[in] + offset;
2499 L1Triplet* new_trip = &vTriplets[index];
2500
2501 // check new triplet
2502 const fscal qp1 = curr_trip->GetQp();
2503 const fscal qp2 = new_trip->GetQp();
2504 fscal dqp = fabs(qp1 - qp2);
2505 fscal Cqp = curr_trip->Cqp;
2506 Cqp += new_trip->Cqp;
2507 if (dqp > PickNeighbour * Cqp)
2508 continue; // bad neighbour // CHECKME why do we need recheck it?? (it really change result)
2509
2510 if (GetFUsed(vSFlag[(*vStsHitsUnused)[new_trip->GetLHit()].f]
2511 | vSFlagB[(*vStsHitsUnused)[new_trip->GetLHit()].b]))
2512 {
2513 // no used hits allowed -> compare and store track
2514 if ((curr_L > best_L) || ((curr_L == best_L) && (curr_chi2 < best_chi2))) {
2515 best_tr = curr_tr;
2516 best_chi2 = curr_chi2;
2517 best_L = curr_L;
2518 }
2519 } else { // add new triplet to the current track
2520
2521 // restore current track
2522 // save current tracklet
2523 L1Branch new_tr = curr_tr;
2524 unsigned char new_L = curr_L;
2525 fscal new_chi2 = curr_chi2;
2526
2527 // add new hit
2528 new_tr.StsHits.push_back(RealIHit[new_trip->GetLHit()]);
2529 new_L += 1;
2530 dqp = dqp / Cqp * 5.; // CHECKME: understand 5, why no sqrt(5)?
2531 new_chi2 += dqp * dqp;
2532 // dqp = dqp/Cqp; // IKu
2533 // new_chi2 += dqp*dqp*5.;
2534
2535 if (new_chi2 > TRACK_CHI2_CUT * new_L)
2536 continue;
2537
2538 const int new_ista = ista + new_trip->GetMSta() - new_trip->GetLSta();
2539 CAFindTrack(new_ista, best_tr, best_L, best_chi2, new_trip, new_tr, new_L, new_chi2, NCalls);
2540 NCalls++;
2541 } // add triplet to track
2542 } // for neighbours
2543
2544 return;
2545}
L1_ASSERT
#define L1_ASSERT(v, msg)
Definition CbmL1Def.h:42
L1_assert
#define L1_assert(v)
Definition CbmL1Def.h:49
L1AddPipeMaterial
void L1AddPipeMaterial(L1TrackPar &T, fvec qp0, fvec w=1, fvec mass2=0.1395679f *0.1395679f)
Definition L1AddMaterial.h:91
L1AddMaterial
void L1AddMaterial(L1TrackPar &T, fvec radThick, fvec qp0, fvec w=1, fvec mass2=0.1395679f *0.1395679f)
Definition L1AddMaterial.h:11
L1Extrapolate0
void L1Extrapolate0(L1TrackPar &T, fvec z_out, L1FieldRegion &F)
Definition L1Extrapolation.h:279
L1ExtrapolateC10Line
void L1ExtrapolateC10Line(const L1TrackPar &T, fvec z_out, fvec &C10)
Definition L1Extrapolation.h:416
L1ExtrapolateLine
void L1ExtrapolateLine(L1TrackPar &T, fvec z_out)
Definition L1Extrapolation.h:374
L1ExtrapolateJXY0
void L1ExtrapolateJXY0(fvec &tx, fvec &ty, fvec dz, L1FieldRegion &F, fvec &extrDx, fvec &extrDy, fvec &J04, fvec &J14)
Definition L1Extrapolation.h:573
L1ExtrapolateYC11Line
void L1ExtrapolateYC11Line(const L1TrackPar &T, fvec z_out, fvec &y, fvec &C11)
Definition L1Extrapolation.h:409
L1Extrapolate
void L1Extrapolate(L1TrackPar &T, fvec z_out, fvec qp0, L1FieldRegion &F, fvec *w=0)
Definition L1Extrapolation.h:11
L1ExtrapolateXC00Line
void L1ExtrapolateXC00Line(const L1TrackPar &T, fvec z_out, fvec &x, fvec &C00)
Definition L1Extrapolation.h:402
L1Filter
void L1Filter(L1TrackPar &T, L1UMeasurementInfo &info, fvec u, fvec w=1., fvec *du2=0)
Definition L1Filtration.h:13
L1FilterChi2XYC00C10C11
void L1FilterChi2XYC00C10C11(const L1UMeasurementInfo &info, fvec &x, fvec &y, fvec &C00, fvec &C10, fvec &C11, fvec &chi2, const fvec &u, const fvec *du2=0)
Definition L1Filtration.h:103
L1FilterChi2
void L1FilterChi2(const L1UMeasurementInfo &info, const fvec &x, const fvec &y, const fvec &C00, const fvec &C10, const fvec &C11, fvec &chi2, const fvec &u, const fvec *du2=0)
Definition L1Filtration.h:75
L1FilterVtx
void L1FilterVtx(L1TrackPar &T, fvec x, fvec y, L1XYMeasurementInfo &info, fvec extrDx, fvec extrDy, fvec J[])
Definition L1Filtration.h:143
THitI
unsigned int THitI
Definition L1StsHit.h:6
fscal
float fscal
Definition P4_F32vec4.h:232
fabs
friend F32vec4 fabs(const F32vec4 &a)
Definition P4_F32vec4.h:52
v
__m128 v
Definition P4_F32vec4.h:1
fvecLen
const int fvecLen
Definition P4_F32vec4.h:233
i
int i
Definition P4_F32vec4.h:22
fvec
F32vec4 fvec
Definition P4_F32vec4.h:231
L1AlgoDraw::SaveCanvas
void SaveCanvas(TString name)
Definition L1AlgoDraw.h:909
L1AlgoDraw::DrawTarget
void DrawTarget()
Definition L1AlgoDraw.h:580
L1AlgoDraw::ClearVeiw
void ClearVeiw()
Definition L1AlgoDraw.h:884
L1AlgoDraw::DrawRecoTracks
void DrawRecoTracks()
Definition L1AlgoDraw.h:322
L1AlgoDraw::InitL1Draw
void InitL1Draw(L1Algo *algo_)
Definition L1AlgoDraw.h:149
L1AlgoDraw::DrawMCTracks
void DrawMCTracks()
Definition L1AlgoDraw.h:173
L1AlgoDraw::DrawAsk
void DrawAsk()
Definition L1AlgoDraw.h:868
L1AlgoDraw::DrawInputHits
void DrawInputHits()
Definition L1AlgoDraw.h:634
L1Algo::StripsToCoor
void StripsToCoor(const fscal &u, const fscal &v, fscal &x, fscal &y, const L1Station &sta) const
convert strip positions to coordinates
Definition L1Algo.cxx:185
L1Algo::vStsHits
vector< L1StsHit > vStsHits
Definition L1Algo.h:162
L1Algo::vTracks
vector< L1Track > vTracks
--— Output data --—
Definition L1Algo.h:182
L1Algo::CreateHitPoint
L1HitPoint CreateHitPoint(const L1StsHit &hit, char ista)
full the hit point by hit information.
Definition L1Algo.cxx:209
L1Algo::CAMergeClones
void CAMergeClones()
Definition L1CAMergeClones.cxx:205
L1Algo::RealIHit
THitI * RealIHit
Definition L1Algo.h:174
L1Algo::targY
fvec targY
Definition L1Algo.h:534
L1Algo::f32
void f32(int n3, int istal, nsL1::vector< L1TrackPar >::TSimd &T_3, vector< THitI > &hitsl_3, vector< THitI > &hitsm_3, vector< THitI > &hitsr_3, int nIterations=0)
Refit Triplets.
Definition L1CATrackFinder.cxx:659
L1Algo::DupletsStaPort
void DupletsStaPort(int istal, int istam, unsigned int ip, vector< int > &n_g1, unsigned int *portionStopIndex, L1TrackPar *T_1, L1FieldRegion *fld_1, THitI *hitsl_1, vector< bool > &lmDuplets, int &n_2, vector< THitI > &i1_2, vector< THitI > &hitsm_2)
Find doublets on station.
Definition L1CATrackFinder.cxx:1112
L1Algo::SetFUnUsed
void SetFUnUsed(unsigned char &flag)
Definition L1Algo.h:278
L1Algo::targX
fvec targX
Definition L1Algo.h:534
L1Algo::_fvecalignment
L1Station vStations[MaxNStations] _fvecalignment
Definition L1Algo.h:155
L1Algo::kAllSecEIter
@ kAllSecEIter
Definition L1Algo.h:488
L1Algo::kAllSecIter
@ kAllSecIter
Definition L1Algo.h:486
L1Algo::kFastPrimIter
@ kFastPrimIter
Definition L1Algo.h:483
L1Algo::kAllSecJumpIter
@ kAllSecJumpIter
Definition L1Algo.h:492
L1Algo::kAllPrimIter
@ kAllPrimIter
Definition L1Algo.h:484
L1Algo::kFastPrimJumpIter
@ kFastPrimJumpIter
Definition L1Algo.h:490
L1Algo::kAllPrimEIter
@ kAllPrimEIter
Definition L1Algo.h:487
L1Algo::kAllPrimJumpIter
@ kAllPrimJumpIter
Definition L1Algo.h:485
L1Algo::kFastPrimIter2
@ kFastPrimIter2
Definition L1Algo.h:491
L1Algo::NStations
int NStations
Definition L1Algo.h:153
L1Algo::StsHitsStopIndex
THitI StsHitsStopIndex[MaxNStations+1]
Definition L1Algo.h:167
L1Algo::vStsZPos
vector< fscal > vStsZPos
Definition L1Algo.h:161
L1Algo::CATime
double CATime
Definition L1Algo.h:185
L1Algo::vStsStripsB
vector< L1Strip > vStsStripsB
Definition L1Algo.h:160
L1Algo::vAllHitPoints
vector< L1HitPoint > vAllHitPoints
Definition L1Algo.h:178
L1Algo::TRACK_CHI2_CUT
const float TRACK_CHI2_CUT
Definition L1Algo.h:522
L1Algo::GetFStation
unsigned char GetFStation(unsigned char flag)
– Flags routines –
Definition L1Algo.h:271
L1Algo::vSFlagB
vector< unsigned char > vSFlagB
Definition L1Algo.h:165
L1Algo::TripletsStaPort
void TripletsStaPort(int istal, int istam, int istar, unsigned int &nstaltriplets, L1TrackPar *T_1, L1FieldRegion *fld_1, THitI *hitsl_1, int &n_2, unsigned int *portionStopIndex, vector< THitI > &i1_2, vector< THitI > &hitsm_2, const vector< bool > &mrDuplets, vector< L1Triplet > *vTriplets_part, unsigned int *TripStartIndexH, unsigned int *TripStopIndexH)
Find triplets on station.
Definition L1CATrackFinder.cxx:1209
L1Algo::TRIPLET_CHI2_CUT
float TRIPLET_CHI2_CUT
Definition L1Algo.h:523
L1Algo::StsHitsUnusedStartIndex
THitI StsHitsUnusedStartIndex[MaxNStations+1]
Definition L1Algo.h:175
L1Algo::MaxDZ
fvec MaxDZ
Definition L1Algo.h:526
L1Algo::GetFUsed
bool GetFUsed(unsigned char flag)
Definition L1Algo.h:272
L1Algo::vSFlag
vector< unsigned char > vSFlag
Definition L1Algo.h:164
L1Algo::StsHitsStartIndex
THitI StsHitsStartIndex[MaxNStations+1]
Definition L1Algo.h:166
L1Algo::vRecoHits
vector< THitI > vRecoHits
Definition L1Algo.h:183
L1Algo::targZ
fvec targZ
Definition L1Algo.h:534
L1Algo::DOUBLET_CHI2_CUT
float DOUBLET_CHI2_CUT
Definition L1Algo.h:524
L1Algo::BranchExtender
fscal BranchExtender(L1Branch &t)
Find additional hits for existing track return chi2.
Definition L1TrackExtender.cxx:401
L1Algo::TripStopIndex
int TripStopIndex[MaxNStations - 2]
Definition L1Algo.h:562
L1Algo::f30
void f30(L1HitPoint *vStsHits_r, L1Station &stam, L1Station &star, int istam, int istar, L1HitPoint *vStsHits_m, L1TrackPar *T_1, L1FieldRegion *fld_1, THitI *hitsl_1, int n2, vector< THitI > &hitsm_2, vector< THitI > &i1_2, const vector< bool > &mrDuplets, int &n3, nsL1::vector< L1TrackPar >::TSimd &T_3, vector< THitI > &hitsl_3, vector< THitI > &hitsm_3, vector< THitI > &hitsr_3, nsL1::vector< fvec >::TSimd &u_front_3, nsL1::vector< fvec >::TSimd &u_back_3, nsL1::vector< fvec >::TSimd &z_Pos_3, nsL1::vector< fvec >::TSimd *du_front_3=0, nsL1::vector< fvec >::TSimd *du_back_3=0)
Add the middle hits to parameters estimation.
Definition L1CATrackFinder.cxx:396
L1Algo::fGhostSuppression
int fGhostSuppression
Definition L1Algo.h:564
L1Algo::NMvdStations
int NMvdStations
Definition L1Algo.h:154
L1Algo::MaxSlope
fvec MaxSlope
Definition L1Algo.h:533
L1Algo::f31
void f31(int n3_V, L1Station &star, nsL1::vector< fvec >::TSimd &u_front_3, nsL1::vector< fvec >::TSimd &u_back_3, nsL1::vector< fvec >::TSimd &z_Pos_3, nsL1::vector< L1TrackPar >::TSimd &T_3, nsL1::vector< fvec >::TSimd *du_front=0, nsL1::vector< fvec >::TSimd *du_back=0)
Add the right hits to parameters estimation.
Definition L1CATrackFinder.cxx:631
L1Algo::vStsHitPointsUnused
std::vector< L1HitPoint > * vStsHitPointsUnused
Definition L1Algo.h:173
L1Algo::Pick_gather
fvec Pick_gather
parameters which are different for different iterations. Set in the begin of CAL1TrackFinder
Definition L1Algo.h:530
L1Algo::f4
void f4(int n3, int istal, int istam, int istar, nsL1::vector< L1TrackPar >::TSimd &T_3, vector< THitI > &hitsl_3, vector< THitI > &hitsm_3, vector< THitI > &hitsr_3, unsigned int &nstaltriplets, vector< L1Triplet > &vTriplets_part, unsigned int *TripStartIndexH, unsigned int *TripStopIndexH)
Select triplets. Save them into vTriplets.
Definition L1CATrackFinder.cxx:902
L1Algo::MaxPortionDoublets
@ MaxPortionDoublets
Definition L1Algo.h:547
L1Algo::Portion
@ Portion
Definition L1Algo.h:545
L1Algo::MaxPortionTriplets
@ MaxPortionTriplets
Definition L1Algo.h:548
L1Algo::MaxNPortion
@ MaxNPortion
Definition L1Algo.h:549
L1Algo::isec
int isec
— data used during finding iterations
Definition L1Algo.h:171
L1Algo::TripStartIndex
int TripStartIndex[MaxNStations - 2]
Definition L1Algo.h:560
L1Algo::vStsHitsUnused
vector< L1StsHit > * vStsHitsUnused
Definition L1Algo.h:172
L1Algo::PickNeighbour
float PickNeighbour
Definition L1Algo.h:531
L1Algo::fNFindIterations
@ fNFindIterations
Definition L1Algo.h:502
L1Algo::vStsStrips
vector< L1Strip > vStsStrips
Definition L1Algo.h:159
L1Algo::vTriplets
vector< L1Triplet > vTriplets
Definition L1Algo.h:559
L1Algo::MaxInvMom
fvec MaxInvMom
Definition L1Algo.h:532
L1Algo::f11
void f11(int istal, int istam, int n1_V, fvec *u_front, fvec *u_back, fvec *zPos, L1TrackPar *T_1, L1FieldRegion *fld_1, fvec *du_front=0, fvec *du_back=0)
Get the field approximation. Add the target to parameters estimation. Propagate to middle station.
Definition L1CATrackFinder.cxx:94
L1Algo::vStsHits2Unus
vector< int > vStsHits2Unus
Definition L1Algo.h:177
L1Algo::FIRSTCASTATION
int FIRSTCASTATION
================================= DATA PART =================================
Definition L1Algo.h:476
L1Algo::CATrackFinder
void CATrackFinder()
The main procedure - find tracks.
Definition L1CATrackFinder.cxx:1348
L1Algo::CAFindTrack
void CAFindTrack(int ista, L1Branch &best_tr, unsigned char &best_L, fscal &best_chi2, const L1Triplet *curr_trip, L1Branch &curr_tr, unsigned char &curr_L, fscal &curr_chi2, int &NCalls)
--— Subroutines used by L1Algo::CATrackFinder() ---—
Definition L1CATrackFinder.cxx:2443
L1Algo::StsHitsUnusedStopIndex
THitI StsHitsUnusedStopIndex[MaxNStations+1]
Definition L1Algo.h:175
L1Algo::vGrid
L1Grid vGrid[MaxNStations]
Definition L1Algo.h:163
L1Algo::f20
void f20(int n1, L1Station &stam, L1HitPoint *vStsHits_l, L1HitPoint *vStsHits_m, int NHits_m, L1TrackPar *T_1, THitI *hitsl_1, int &n2, vector< THitI > &i1_2, int &start_mhit, vector< THitI > &hitsm_2, vector< bool > &lmDuplets)
Find the doublets. Reformat data in the portion of doublets.
Definition L1CATrackFinder.cxx:288
L1Algo::f5
void f5(unsigned int *TripStartIndexH, unsigned int *TripStopIndexH, int *nlevel)
Find neighbours of triplets. Calculate level of triplets.
Definition L1CATrackFinder.cxx:987
L1Algo::MaxNStations
@ MaxNStations
Definition L1Algo.h:151
L1Algo::f10
void f10(int start_lh, int n1, L1HitPoint *vStsHits_l, fvec *u_front, fvec *u_back, fvec *zPos, THitI *hitsl_1, fvec *du_front=0, fvec *du_back=0)
Prepare the portion of left hits data.
Definition L1CATrackFinder.cxx:63
L1Algo::SetFUsed
void SetFUsed(unsigned char &flag)
Definition L1Algo.h:276
L1Algo::fRadThick
vector< L1Material > fRadThick
Definition L1Algo.h:157
L1CATFIterTimerInfo::PrintReal
void PrintReal(int f=0)
Definition L1Timer.h:101
L1CATFIterTimerInfo::Add
void Add(string name)
Definition L1Timer.h:78
L1CATFTimerInfo::Add
void Add(string name)
Definition L1Timer.h:136
L1CATFTimerInfo::PrintReal
void PrintReal()
Definition L1Timer.h:169
L1CATFTimerInfo::SetNIter
void SetNIter(int n)
Definition L1Timer.h:134
L1FieldRegion::SetOneEntry
void SetOneEntry(const int i0, const L1FieldRegion &f1, const int i1)
Definition L1Field.h:193
L1FieldRegion::Set
void Set(const L1FieldValue &B0, const fvec B0z, const L1FieldValue &B1, const fvec B1z, const L1FieldValue &B2, const fvec B2z)
Definition L1Field.h:131
L1FieldSlice::GetFieldValue
void GetFieldValue(const fvec &x, const fvec &y, L1FieldValue &B) const
Definition L1Field.h:47
L1Grid
2-dimensional grid of pointers.
Definition L1Grid.h:49
L1Grid::Create
void Create(float yMin, float yMax, float zMin, float zMax, float sy, float sz)
Definition L1Grid.h:116
L1Grid::Fill
void Fill(const L1HitPoint *points, THitI n)
Definition L1Grid.h:138
L1HitArea::GetNext
bool GetNext(THitI &i)
look up the next hit in the requested area.
Definition L1HitArea.h:86
L1Station::XYInfo
L1XYMeasurementInfo XYInfo
Definition L1Station.h:31
L1Station::backInfo
L1UMeasurementInfo backInfo
Definition L1Station.h:30
L1Station::fieldSlice
L1FieldSlice fieldSlice
Definition L1Station.h:29
L1Station::materialInfo
L1MaterialInfo materialInfo
Definition L1Station.h:28
L1Station::frontInfo
L1UMeasurementInfo frontInfo
Definition L1Station.h:30
L1StsHit::iz
TZPosI iz
Definition L1StsHit.h:13
L1StsHit::f
TStripI f
Definition L1StsHit.h:12
L1StsHit::b
TStripI b
Definition L1StsHit.h:12
L1TrackPar::SetOneEntry
void SetOneEntry(const int i0, const L1TrackPar &T1, const int i1)
Definition L1TrackPar.h:106
L1Track::TFirst
double TFirst[6]
Definition L1Track.h:25
L1Track::chi2
double chi2
Definition L1Track.h:25
L1Track::Momentum
float Momentum
Definition L1Track.h:24
L1Track::NHits
unsigned char NHits
Definition L1Track.h:23
L1Track::NDF
int NDF
Definition L1Track.h:26
L1Track::CFirst
double CFirst[15]
Definition L1Track.h:25
L1Track::CLast
double CLast[15]
Definition L1Track.h:25
L1Track::TLast
double TLast[6]
Definition L1Track.h:25
L1Triplet::Set
void Set(unsigned int iHitL, unsigned int iHitM, unsigned int iHitR, unsigned int iStaL, unsigned int iStaM, unsigned int iStaR, unsigned char Level, unsigned char Qp, float Chi2)
Definition L1Triplet.h:37
L1Triplet::GetRSta
int GetRSta() const
Definition L1Triplet.h:95
L1Triplet::GetQp
unsigned char GetQp() const
Definition L1Triplet.h:80
L1Triplet::GetLSta
int GetLSta() const
Definition L1Triplet.h:91
L1Triplet::GetLHit
THitI GetLHit() const
Definition L1Triplet.h:63
L1Triplet::GetMHit
THitI GetMHit() const
Definition L1Triplet.h:64
L1Triplet::GetMSta
int GetMSta() const
Definition L1Triplet.h:93
L1Triplet::GetRHit
THitI GetRHit() const
Definition L1Triplet.h:65
L1Triplet::GetQpOrig
float GetQpOrig(float MaxInvMom)
Definition L1Triplet.h:89
L1Triplet::GetLevel
unsigned char GetLevel() const
Definition L1Triplet.h:74
L1Triplet::neighbours
std::vector< unsigned int > neighbours
Definition L1Triplet.h:20
L1Triplet::Cqp
unsigned char Cqp
Definition L1Triplet.h:18
L1Triplet::GetChi2
float GetChi2() const
Definition L1Triplet.h:82
TStopwatch::RealTime
Double_t RealTime()
TStopwatch::Start
void Start(Bool_t reset=kTRUE)
TStopwatch::Stop
void Stop()
sqrt
RVec< PromoteType< T > > sqrt(const RVec< T > &v)
L1Branch::Momentum
float Momentum
Definition L1Branch.h:38
L1Branch::chi2
float chi2
Definition L1Branch.h:38
L1Branch::StsHits
std::vector< THitI > StsHits
Definition L1Branch.h:39
L1Branch::Set
void Set(unsigned char iStation, unsigned char Length, float Chi2, float Qp)
Definition L1Branch.h:41
L1Branch::comparePChi2
static bool comparePChi2(const L1Branch *a, const L1Branch *b)
Definition L1Branch.h:57
L1HitPoint
contain strips positions and coordinates of hit
Definition L1HitPoint.h:9
L1HitPoint::V
fscal V() const
Definition L1HitPoint.h:43
L1HitPoint::DU
fscal DU() const
Definition L1HitPoint.h:44
L1HitPoint::U
fscal U() const
Definition L1HitPoint.h:42
L1HitPoint::Z
fscal Z() const
Definition L1HitPoint.h:41
L1HitPoint::N
unsigned short int N() const
Definition L1HitPoint.h:47
L1HitPoint::Y
fscal Y() const
Definition L1HitPoint.h:40
L1HitPoint::X
fscal X() const
Definition L1HitPoint.h:39
L1HitPoint::DV
fscal DV() const
Definition L1HitPoint.h:45
nsL1::vector::TSimd
std::vector< T > TSimd
Definition PSEUDO_F32vec1.h:148