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