BmnRoot 26.10.0-18-g4e3414f9
Software framework of the BM@N experiment
Loading...
Searching...
No Matches
CbmL1PFFitter.cxx
Go to the documentation of this file.
1/*
2 *=====================================================
3 *
4 * CBM Level 1 Reconstruction
5 *
6 * Authors: M.Zyzak
7 *
8 * e-mail :
9 *
10 *=====================================================
11 *
12 * SIMD Fitter for CbmL1Track class
13 *
14 */
15
16#include "CbmL1PFFitter.h"
17
18#include "CbmL1.h"
19#include "CbmMvdHit.h"
20#include "CbmStsHit.h"
21#include "CbmStsTrack.h"
22#include "TClonesArray.h"
23
24// L1Algo tools
25#include "CbmKFParticleDatabase.h"
26#include "CbmKFVertex.h"
27#include "CbmL1Track.h"
28#include "FairRootManager.h"
29#include "L1Algo.h"
30#include "L1Extrapolation.h"
31#include "L1MaterialInfo.h"
32#include "L1Station.h"
33#include "L1TrackPar.h"
34#include "TDatabasePDG.h"
35
36using std::vector;
37
38namespace NS_L1TrackFitter
39{
40const fvec c_light = 0.000299792458, c_light_i = 1. / c_light;
41const fvec ZERO = 0.;
42const fvec ONE = 1.;
43const fvec vINF = 0.1;
44} // namespace NS_L1TrackFitter
45using namespace NS_L1TrackFitter;
46
47CbmL1PFFitter::CbmL1PFFitter() {}
48
49CbmL1PFFitter::~CbmL1PFFitter() {}
50
51inline void CbmL1PFFitter::AddMaterial(L1TrackPar& T, fvec radThick, fvec qp0, fvec& mass2, fvec& w)
52{
53 static const fvec ZERO = 0.0f, ONE = 1.;
54
55 fvec tx = T.tx;
56 fvec ty = T.ty;
57 fvec txtx = tx * tx;
58 fvec tyty = ty * ty;
59 fvec txtx1 = txtx + ONE;
60 fvec h = txtx + tyty;
61 fvec t = sqrt(txtx1 + tyty);
62 fvec h2 = h * h;
63 fvec qp0t = qp0 * t;
64
65 static const fvec c1 = 0.0136f, c2 = c1 * 0.038f, c3 = c2 * 0.5f, c4 = -c3 / 2.0f, c5 = c3 / 3.0f, c6 = -c3 / 4.0f;
66
67 fvec s0 = (c1 + c2 * log(radThick) + c3 * h + h2 * (c4 + c5 * h + c6 * h2)) * qp0t;
68 fvec a = ((ONE + mass2 * qp0 * qp0t) * radThick * s0 * s0);
69 // std::cout <<" a " << a << std::endl;
70 // a=0.000005;
71 fvec zero = ZERO;
72 fvec init = fvec(w > zero);
73 T.C22 += (init & (txtx1 * a));
74 T.C32 += (init & (tx * ty * a));
75 T.C33 += (init & ((ONE + tyty) * a));
76}
77
78inline void CbmL1PFFitter::AddPipeMaterial(L1TrackPar& T, fvec qp0, fvec& mass2, fvec& w)
79{
80 static const fvec ZERO = 0.0f, ONE = 1.f;
81
82 static const fscal RadThick = 0.0009f; // 0.5/18.76;
83 static const fscal logRadThick = log(RadThick);
84 fvec tx = T.tx;
85 fvec ty = T.ty;
86 fvec txtx = tx * tx;
87 fvec tyty = ty * ty;
88 fvec txtx1 = txtx + ONE;
89 fvec h = txtx + tyty;
90 fvec t = sqrt(txtx1 + tyty);
91 fvec h2 = h * h;
92 fvec qp0t = qp0 * t;
93
94 static const fvec c1 = 0.0136f, c2 = c1 * 0.038f, c3 = c2 * 0.5f, c4 = -c3 / 2.0f, c5 = c3 / 3.0f, c6 = -c3 / 4.0f;
95 fvec s0 = (c1 + c2 * (logRadThick) + c3 * h + h2 * (c4 + c5 * h + c6 * h2)) * qp0t;
96 fvec a = ((ONE + mass2 * qp0 * qp0t) * RadThick * s0 * s0);
97 fvec zero = ZERO;
98 fvec init = fvec(w > zero);
99 T.C22 += (init & (txtx1 * a));
100 T.C32 += (init & (tx * ty * a));
101 T.C33 += (init & ((ONE + tyty) * a));
102}
103
104void CbmL1PFFitter::FilterFirst(L1TrackPar& track, fvec& x, fvec& y, fvec& w, L1Station& st)
105{
106 fvec w1 = ONE - w;
107 fvec sigma21 = w * st.XYInfo.C00 + w1 * vINF;
108 fvec sigma22 = w * st.XYInfo.C10;
109 fvec sigma23 = w * st.XYInfo.C11 + w1 * vINF;
110 // initialize covariance matrix
111 track.C00 = sigma21;
112 track.C10 = sigma22;
113 track.C11 = sigma23;
114 /* track.C00= vINF;
115 track.C10= ZERO; track.C11= vINF;*/
116 track.C20 = ZERO;
117 track.C21 = ZERO;
118 track.C22 = vINF;
119 track.C30 = ZERO;
120 track.C31 = ZERO;
121 track.C32 = ZERO;
122 track.C33 = vINF;
123 track.C40 = ZERO;
124 track.C41 = ZERO;
125 track.C42 = ZERO;
126 track.C43 = ZERO;
127 track.C44 = vINF;
128
129 track.x = w * x + w1 * track.x;
130 track.y = w * y + w1 * track.y;
131 track.NDF = -3.0;
132 track.chi2 = ZERO;
133}
134
135void CbmL1PFFitter::FilterLast(L1TrackPar& track, fvec& x, fvec& y, fvec& w, L1Station& st)
136{
137 fvec w1 = ONE - w;
138 fvec sigma21 = w * st.XYInfo.C00 + w1 * vINF;
139 fvec sigma22 = w * st.XYInfo.C10;
140 fvec sigma23 = w * st.XYInfo.C11 + w1 * vINF;
141 // initialize covariance matrix
142 track.C00 = sigma21;
143 track.C10 = sigma22;
144 track.C11 = sigma23;
145 track.C20 = ZERO;
146 track.C21 = ZERO;
147 track.C22 = vINF;
148 track.C30 = ZERO;
149 track.C31 = ZERO;
150 track.C32 = ZERO;
151 track.C33 = vINF;
152 track.C40 = ZERO;
153 track.C41 = ZERO;
154 track.C42 = ZERO;
155 track.C43 = ZERO;
156 track.C44 = ONE;
157
158 track.x = w * x + w1 * track.x;
159 track.y = w * y + w1 * track.y;
160 track.NDF = -3.0;
161 track.chi2 = ZERO;
162}
163
164void CbmL1PFFitter::Filter(L1TrackPar& T, L1UMeasurementInfo& info, fvec& u, fvec& w)
165{
166 register fvec wi, zeta, zetawi, HCH;
167 register fvec F0, F1, F2, F3, F4;
168 register fvec K1, K2, K3, K4;
169
170 zeta = info.cos_phi * T.x + info.sin_phi * T.y - u;
171
172 // F = CH'
173 F0 = info.cos_phi * T.C00 + info.sin_phi * T.C10;
174 F1 = info.cos_phi * T.C10 + info.sin_phi * T.C11;
175
176 HCH = (F0 * info.cos_phi + F1 * info.sin_phi);
177
178 F2 = info.cos_phi * T.C20 + info.sin_phi * T.C21;
179 F3 = info.cos_phi * T.C30 + info.sin_phi * T.C31;
180 F4 = info.cos_phi * T.C40 + info.sin_phi * T.C41;
181
182 wi = 1. / (info.sigma2 + HCH);
183
184 wi *= w;
185
186 zetawi = zeta * wi;
187 T.chi2 += zeta * zetawi;
188 T.NDF += w;
189
190 K1 = F1 * wi;
191 K2 = F2 * wi;
192 K3 = F3 * wi;
193 K4 = F4 * wi;
194
195 T.x -= F0 * zetawi;
196 T.y -= F1 * zetawi;
197 T.tx -= F2 * zetawi;
198 T.ty -= F3 * zetawi;
199 T.qp -= F4 * zetawi;
200
201 T.C00 -= F0 * F0 * wi;
202 T.C10 -= K1 * F0;
203 T.C11 -= K1 * F1;
204 T.C20 -= K2 * F0;
205 T.C21 -= K2 * F1;
206 T.C22 -= K2 * F2;
207 T.C30 -= K3 * F0;
208 T.C31 -= K3 * F1;
209 T.C32 -= K3 * F2;
210 T.C33 -= K3 * F3;
211 T.C40 -= K4 * F0;
212 T.C41 -= K4 * F1;
213 T.C42 -= K4 * F2;
214 T.C43 -= K4 * F3;
215 T.C44 -= K4 * F4;
216}
217
218void CbmL1PFFitter::Fit(vector<CbmL1Track>& Tracks, fvec mass)
219{
220 // cout << " Start L1 Track Fitter " << endl;
221 // int start_hit = 0; // for interation in vRecoHits[]
222
223 fvec mass2 = mass * mass;
224
225 L1FieldValue fB0, fB1, fB2 _fvecalignment;
226 L1FieldRegion fld _fvecalignment;
227
228 TClonesArray* listStsHits = CbmL1::Instance()->listStsHits;
229 TClonesArray* listMvdHits = CbmL1::Instance()->listMvdHits;
230 int NMvdStations = CbmL1::Instance()->algo->NMvdStations;
231
232 static int nHits = CbmL1::Instance()->algo->NStations;
233 int iVec = 0, i = 0;
234 int nTracks_SIMD = fvecLen;
235 L1TrackPar T; // fitting parametr coresponding to current track
236
237 CbmL1Track* t[fvecLen];
238
239 int ista;
240 L1Station* sta = CbmL1::Instance()->algo->vStations;
241 fvec* x = new fvec[nHits];
242 fvec* u = new fvec[nHits];
243 fvec* v = new fvec[nHits];
244 fvec* y = new fvec[nHits];
245 fvec* z = new fvec[nHits];
246 fvec* w = new fvec[nHits];
247 fvec y_temp;
248 fvec fz0, fz1, fz2, dz, z_start, z_end;
249 L1FieldValue* fB = new L1FieldValue[nHits];
250 L1FieldValue fB_temp _fvecalignment;
251
252 unsigned short N_vTracks = Tracks.size();
253
254 for (unsigned short itrack = 0; itrack < N_vTracks; itrack++) {
255 Tracks[itrack].mass = mass[0];
256 }
257
258 for (unsigned short itrack = 0; itrack < N_vTracks; itrack += fvecLen) {
259 if (N_vTracks - itrack < static_cast<unsigned short>(fvecLen))
260 nTracks_SIMD = N_vTracks - itrack;
261
262 for (i = 0; i < nTracks_SIMD; i++) {
263 t[i] = &Tracks[itrack + i]; // current track
264 T.x[i] = t[i]->TLast[0];
265 T.y[i] = t[i]->TLast[1];
266 T.tx[i] = t[i]->TLast[2];
267 T.ty[i] = t[i]->TLast[3];
268 T.qp[i] = t[i]->TLast[4];
269 T.z[i] = t[i]->TLast[5];
270 T.C00[i] = t[i]->CLast[0];
271 T.C10[i] = t[i]->CLast[1];
272 T.C11[i] = t[i]->CLast[2];
273 T.C20[i] = t[i]->CLast[3];
274 T.C21[i] = t[i]->CLast[4];
275 T.C22[i] = t[i]->CLast[5];
276 T.C30[i] = t[i]->CLast[6];
277 T.C31[i] = t[i]->CLast[7];
278 T.C32[i] = t[i]->CLast[8];
279 T.C33[i] = t[i]->CLast[9];
280 T.C40[i] = t[i]->CLast[10];
281 T.C41[i] = t[i]->CLast[11];
282 T.C42[i] = t[i]->CLast[12];
283 T.C43[i] = t[i]->CLast[13];
284 T.C44[i] = t[i]->CLast[14];
285 }
286
287 // get hits of current track
288 for (i = 0; i < nHits; i++) {
289 w[i] = ZERO;
290 z[i] = sta[i].z;
291 }
292
293 for (iVec = 0; iVec < nTracks_SIMD; iVec++) {
294 int nHitsTrack = t[iVec]->GetNOfHits();
295 for (i = 0; i < nHitsTrack; i++) {
296 float posx = 0.f, posy = 0.f, posz = 0.f;
297 CbmL1HitStore& hs = CbmL1::Instance()->vHitStore[t[iVec]->StsHits[i]];
298 if (hs.ExtIndex < 0) {
299 CbmMvdHit* hit = L1_DYNAMIC_CAST<CbmMvdHit*>(listMvdHits->At(-hs.ExtIndex - 1));
300 posx = hit->GetX();
301 posy = hit->GetY();
302 posz = hit->GetZ();
303 ista = posz < 7.f ? 0 : 1;
304 } else {
305 CbmStsHit* hit = L1_DYNAMIC_CAST<CbmStsHit*>(listStsHits->At(hs.ExtIndex));
306 posx = hit->GetX();
307 posy = hit->GetY();
308 posz = hit->GetZ();
309 ista = hit->GetStationNr() - 1 + NMvdStations;
310 }
311 // std::cout << ista << std::endl;
312 w[ista][iVec] = 1.f;
313
314 x[ista][iVec] = posx;
315 y[ista][iVec] = posy;
316 u[ista][iVec] = posx * sta[ista].frontInfo.cos_phi[0] + posy * sta[ista].frontInfo.sin_phi[0];
317 v[ista][iVec] = posx * sta[ista].backInfo.cos_phi[0] + posy * sta[ista].backInfo.sin_phi[0];
318 z[ista][iVec] = posz;
319 sta[ista].fieldSlice.GetFieldValue(x[ista], y[ista], fB_temp);
320 fB[ista].x[iVec] = fB_temp.x[iVec];
321 fB[ista].y[iVec] = fB_temp.y[iVec];
322 fB[ista].z[iVec] = fB_temp.z[iVec];
323 if (i == 0)
324 z_start[iVec] = posz;
325 if (i == nHitsTrack - 1)
326 z_end[iVec] = posz;
327 }
328 }
329
330 // fit backward
331
332 fvec qp0 = T.qp;
333
334 i = nHits - 1;
335
336 FilterFirst(T, x[i], y[i], w[i], sta[i]);
337 AddMaterial(T, CbmL1::Instance()->algo->fRadThick[i].GetRadThick(T.x, T.y), qp0, mass2, w[i]);
338
339 fz1 = z[i];
340 sta[i].fieldSlice.GetFieldValue(T.x, T.y, fB1);
341 fB1.Combine(fB[i], w[i]);
342
343 fz2 = z[i - 2];
344 dz = fz2 - fz1;
345 sta[i].fieldSlice.GetFieldValue(T.x + T.tx * dz, T.y + T.ty * dz, fB2);
346 fB2.Combine(fB[i - 2], w[i - 2]);
347 fld.Set(fB2, fz2, fB1, fz1, fB0, fz0);
348 for (--i; i >= 0; i--) {
349 fz0 = z[i];
350 dz = (fz1 - fz0);
351 sta[i].fieldSlice.GetFieldValue(T.x - T.tx * dz, T.y - T.ty * dz, fB0);
352 fB0.Combine(fB[i], w[i]);
353 fld.Set(fB0, fz0, fB1, fz1, fB2, fz2);
354
355 // fvec one = ONE;
356 fvec zero = ZERO;
357 fvec initialised = fvec(z[i] > z_end) & fvec(z_start > z[i]);
358 fvec w1 = (zero & initialised) + (w[i] & (!initialised));
359
360 L1Extrapolate(T, z[i], qp0, fld, &w1);
361 if (i == NMvdStations - 1)
362 AddPipeMaterial(T, qp0, mass2, w1);
363 AddMaterial(T, CbmL1::Instance()->algo->fRadThick[i].GetRadThick(T.x, T.y), qp0, mass2, w1);
364 Filter(T, sta[i].frontInfo, u[i], w[i]);
365 Filter(T, sta[i].backInfo, v[i], w[i]);
366 fB2 = fB1;
367 fz2 = fz1;
368 fB1 = fB0;
369 fz1 = fz0;
370 }
371
372 for (iVec = 0; iVec < nTracks_SIMD; iVec++) {
373 t[iVec]->T[0] = T.x[iVec];
374 t[iVec]->T[1] = T.y[iVec];
375 t[iVec]->T[2] = T.tx[iVec];
376 t[iVec]->T[3] = T.ty[iVec];
377 t[iVec]->T[4] = T.qp[iVec];
378 // t[iVec]->TFirst[5] = z_start[iVec];
379 t[iVec]->T[5] = T.z[iVec];
380
381 t[iVec]->C[0] = T.C00[iVec];
382 t[iVec]->C[1] = T.C10[iVec];
383 t[iVec]->C[2] = T.C11[iVec];
384 t[iVec]->C[3] = T.C20[iVec];
385 t[iVec]->C[4] = T.C21[iVec];
386 t[iVec]->C[5] = T.C22[iVec];
387 t[iVec]->C[6] = T.C30[iVec];
388 t[iVec]->C[7] = T.C31[iVec];
389 t[iVec]->C[8] = T.C32[iVec];
390 t[iVec]->C[9] = T.C33[iVec];
391 t[iVec]->C[10] = T.C40[iVec];
392 t[iVec]->C[11] = T.C41[iVec];
393 t[iVec]->C[12] = T.C42[iVec];
394 t[iVec]->C[13] = T.C43[iVec];
395 t[iVec]->C[14] = T.C44[iVec];
396
397 t[iVec]->chi2 = T.chi2[iVec];
398 t[iVec]->NDF = static_cast<int>(T.NDF[iVec]);
399 }
400
401 // fit forward
402
403 i = 0;
404 FilterLast(T, x[i], y[i], w[i], sta[i]);
405 qp0 = T.qp;
406 AddMaterial(T, CbmL1::Instance()->algo->fRadThick[i].GetRadThick(T.x, T.y), qp0, mass2, w[i]);
407
408 fz1 = z[i];
409 sta[i].fieldSlice.GetFieldValue(T.x, T.y, fB1);
410 fB1.Combine(fB[i], w[i]);
411
412 fz2 = z[i + 2];
413 dz = fz2 - fz1;
414 sta[i].fieldSlice.GetFieldValue(T.x + T.tx * dz, T.y + T.ty * dz, fB2);
415 fB2.Combine(fB[i + 2], w[i + 2]);
416 fld.Set(fB2, fz2, fB1, fz1, fB0, fz0);
417
418 for (++i; i < nHits; i++) {
419 fz0 = z[i];
420 dz = (fz1 - fz0);
421 sta[i].fieldSlice.GetFieldValue(T.x - T.tx * dz, T.y - T.ty * dz, fB0);
422 fB0.Combine(fB[i], w[i]);
423 fld.Set(fB0, fz0, fB1, fz1, fB2, fz2);
424
425 // fvec one = ONE;
426 fvec zero = ZERO;
427 fvec initialised = fvec(z[i] > z_end) & fvec(z_start > z[i]);
428 fvec w1 = (zero & initialised) + (w[i] & (!initialised));
429
430 L1Extrapolate(T, z[i], qp0, fld, &w1);
431 if (i == NMvdStations)
432 AddPipeMaterial(T, qp0, mass2, w1);
433 AddMaterial(T, CbmL1::Instance()->algo->fRadThick[i].GetRadThick(T.x, T.y), qp0, mass2, w1);
434 Filter(T, sta[i].frontInfo, u[i], w[i]);
435 Filter(T, sta[i].backInfo, v[i], w[i]);
436
437 fB2 = fB1;
438 fz2 = fz1;
439 fB1 = fB0;
440 fz1 = fz0;
441 }
442
443 for (iVec = 0; iVec < nTracks_SIMD; iVec++) {
444 t[iVec]->TLast[0] = T.x[iVec];
445 t[iVec]->TLast[1] = T.y[iVec];
446 t[iVec]->TLast[2] = T.tx[iVec];
447 t[iVec]->TLast[3] = T.ty[iVec];
448 t[iVec]->TLast[4] = T.qp[iVec];
449 // t[iVec]->TLast[5] = z_end[iVec];
450 t[iVec]->TLast[5] = T.z[iVec];
451
452 t[iVec]->CLast[0] = T.C00[iVec];
453 t[iVec]->CLast[1] = T.C10[iVec];
454 t[iVec]->CLast[2] = T.C11[iVec];
455 t[iVec]->CLast[3] = T.C20[iVec];
456 t[iVec]->CLast[4] = T.C21[iVec];
457 t[iVec]->CLast[5] = T.C22[iVec];
458 t[iVec]->CLast[6] = T.C30[iVec];
459 t[iVec]->CLast[7] = T.C31[iVec];
460 t[iVec]->CLast[8] = T.C32[iVec];
461 t[iVec]->CLast[9] = T.C33[iVec];
462 t[iVec]->CLast[10] = T.C40[iVec];
463 t[iVec]->CLast[11] = T.C41[iVec];
464 t[iVec]->CLast[12] = T.C42[iVec];
465 t[iVec]->CLast[13] = T.C43[iVec];
466 t[iVec]->CLast[14] = T.C44[iVec];
467
468 t[iVec]->chi2 = T.chi2[iVec];
469 t[iVec]->NDF = static_cast<int>(T.NDF[iVec]);
470 }
471 }
472
473 delete[] x;
474 delete[] u;
475 delete[] v;
476 delete[] y;
477 delete[] z;
478 delete[] w;
479 delete[] fB;
480}
481
482void CbmL1PFFitter::GetChiToVertex(vector<CbmL1Track>& Tracks, vector<float>& chiToVtx, CbmKFVertex& primVtx)
483{
484 vector<L1FieldRegion> B;
485 CalculateFieldRegion(Tracks, B);
486
487 int nTracks_SIMD = fvecLen;
488 L1TrackPar T; // fitting parametr coresponding to current track
489
490 CbmL1Track* t[fvecLen];
491
492 int nStations = CbmL1::Instance()->algo->NStations;
493 int NMvdStations = CbmL1::Instance()->algo->NMvdStations;
494 L1Station* sta = CbmL1::Instance()->algo->vStations;
495 fvec* zSta = new fvec[nStations];
496 for (int iSta = 0; iSta < nStations; iSta++)
497 zSta[iSta] = sta[iSta].z;
498
499 unsigned short N_vTracks = Tracks.size();
500
501 for (unsigned short itrack = 0; itrack < N_vTracks; itrack += fvecLen) {
502 if (N_vTracks - itrack < static_cast<unsigned short>(fvecLen))
503 nTracks_SIMD = N_vTracks - itrack;
504
505 fvec mass2;
506 for (int iVec = 0; iVec < nTracks_SIMD; iVec++) {
507 t[iVec] = &Tracks[itrack + iVec]; // current track
508 T.x[iVec] = t[iVec]->T[0];
509 T.y[iVec] = t[iVec]->T[1];
510 T.tx[iVec] = t[iVec]->T[2];
511 T.ty[iVec] = t[iVec]->T[3];
512 T.qp[iVec] = t[iVec]->T[4];
513 T.z[iVec] = t[iVec]->T[5];
514 T.C00[iVec] = t[iVec]->C[0];
515 T.C10[iVec] = t[iVec]->C[1];
516 T.C11[iVec] = t[iVec]->C[2];
517 T.C20[iVec] = t[iVec]->C[3];
518 T.C21[iVec] = t[iVec]->C[4];
519 T.C22[iVec] = t[iVec]->C[5];
520 T.C30[iVec] = t[iVec]->C[6];
521 T.C31[iVec] = t[iVec]->C[7];
522 T.C32[iVec] = t[iVec]->C[8];
523 T.C33[iVec] = t[iVec]->C[9];
524 T.C40[iVec] = t[iVec]->C[10];
525 T.C41[iVec] = t[iVec]->C[11];
526 T.C42[iVec] = t[iVec]->C[12];
527 T.C43[iVec] = t[iVec]->C[13];
528 T.C44[iVec] = t[iVec]->C[14];
529 mass2[iVec] = t[iVec]->mass * t[iVec]->mass;
530 }
531
532 L1FieldRegion& fld = B[itrack / fvecLen];
533 for (int iSt = nStations - 4; iSt >= 0; iSt--) {
534 // fvec zero = ZERO;
535 fvec w = ONE;
536 fvec initialized = fvec(T.z > (zSta[iSt] + 2.5));
537 w = fvec(w & initialized);
538
539 L1Extrapolate(T, zSta[iSt], T.qp, fld, &w);
540 if (iSt == NMvdStations - 1)
541 AddPipeMaterial(T, T.qp, mass2, w);
542 AddMaterial(T, CbmL1::Instance()->algo->fRadThick[iSt].GetRadThick(T.x, T.y), T.qp, mass2, w);
543 }
544 fvec ONE = 1;
545 if (NMvdStations <= 0)
546 AddPipeMaterial(T, T.qp, mass2, ONE);
547 L1Extrapolate(T, primVtx.GetRefZ(), T.qp, fld);
548
549 Double_t Cv[3] = {primVtx.GetCovMatrix()[0], primVtx.GetCovMatrix()[1], primVtx.GetCovMatrix()[2]};
550
551 fvec dx = T.x - primVtx.GetRefX();
552 fvec dy = T.y - primVtx.GetRefY();
553 fvec c[3] = {T.C00, T.C10, T.C11};
554 c[0] += Cv[0];
555 c[1] += Cv[1];
556 c[2] += Cv[2];
557 fvec d = c[0] * c[2] - c[1] * c[1];
558 fvec chi = sqrt(fabs(0.5 * (dx * dx * c[0] - 2 * dx * dy * c[1] + dy * dy * c[2]) / d));
559 fvec isNull = fvec(fabs(d) < 1.e-20);
560 chi = fvec(fvec(!isNull) & chi) + fvec(isNull & fvec(0));
561
562 for (int iVec = 0; iVec < nTracks_SIMD; iVec++)
563 chiToVtx.push_back(chi[iVec]);
564 }
565 delete[] zSta;
566}
567
568void CbmL1PFFitter::CalculateFieldRegion(vector<CbmL1Track>& Tracks, vector<L1FieldRegion>& field)
569{
570 L1FieldRegion fld _fvecalignment;
571
572 TClonesArray* listStsHits = CbmL1::Instance()->listStsHits;
573 TClonesArray* listMvdHits = CbmL1::Instance()->listMvdHits;
574 int NMvdStations = CbmL1::Instance()->algo->NMvdStations;
575
576 int nTracks_SIMD = fvecLen;
577 L1TrackPar T; // fitting parametr coresponding to current track
578
579 CbmL1Track* t[fvecLen];
580
581 int ista;
582 L1Station* sta = CbmL1::Instance()->algo->vStations;
583 L1FieldValue fB[3], fB_temp _fvecalignment;
584 fvec zField[3];
585
586 unsigned short N_vTracks = Tracks.size();
587
588 for (unsigned short itrack = 0; itrack < N_vTracks; itrack += fvecLen) {
589 if (N_vTracks - itrack < static_cast<unsigned short>(fvecLen))
590 nTracks_SIMD = N_vTracks - itrack;
591
592 for (int i = 0; i < nTracks_SIMD; i++)
593 t[i] = &Tracks[itrack + i]; // current track
594
595 for (int iVec = 0; iVec < nTracks_SIMD; iVec++) {
596 for (int iH = 0; iH < 2; iH++) // it is assumed, that track has at least 3 hits
597 {
598 float posx = 0.f, posy = 0.f, posz = 0.f;
599 CbmL1HitStore& hs = CbmL1::Instance()->vHitStore[t[iVec]->StsHits[iH]];
600 if (hs.ExtIndex < 0) {
601 CbmMvdHit* hit = L1_DYNAMIC_CAST<CbmMvdHit*>(listMvdHits->At(-hs.ExtIndex - 1));
602 posx = hit->GetX();
603 posy = hit->GetY();
604 posz = hit->GetZ();
605 ista = posz < 7.f ? 0 : 1;
606 } else {
607 CbmStsHit* hit = L1_DYNAMIC_CAST<CbmStsHit*>(listStsHits->At(hs.ExtIndex));
608 posx = hit->GetX();
609 posy = hit->GetY();
610 posz = hit->GetZ();
611 ista = hit->GetStationNr() - 1 + NMvdStations;
612 }
613 sta[ista].fieldSlice.GetFieldValue(posx, posy, fB_temp);
614 fB[iH + 1].x[iVec] = fB_temp.x[iVec];
615 fB[iH + 1].y[iVec] = fB_temp.y[iVec];
616 fB[iH + 1].z[iVec] = fB_temp.z[iVec];
617 zField[iH + 1][iVec] = posz;
618 }
619 }
620
621 fB[0] = CbmL1::Instance()->algo->GetvtxFieldValue();
622 zField[0] = 0;
623 fld.Set(fB[2], zField[2], fB[1], zField[1], fB[0], zField[0]);
624 field.push_back(fld);
625 }
626}
627
628void CbmL1PFFitter::Fit(vector<CbmStsTrack>& Tracks, int pidHypo)
629{
630 // cout << " Start L1 Track Fitter " << endl;
631 // int start_hit = 0; // for interation in vRecoHits[]
632
633 fvec mass = TDatabasePDG::Instance()->GetParticle(pidHypo)->Mass();
634
635 fvec mass2 = mass * mass;
636
637 L1FieldValue fB0, fB1, fB2 _fvecalignment;
638 L1FieldRegion fld _fvecalignment;
639
640 FairRootManager* fManger = FairRootManager::Instance();
641 TClonesArray* listStsHits = (TClonesArray*)fManger->GetObject("StsHit");
642 int NMvdStations = CbmL1::Instance()->algo->NMvdStations;
643 TClonesArray* listMvdHits = 0;
644 if (NMvdStations > 0.)
645 listMvdHits = (TClonesArray*)fManger->GetObject("MvdHit");
646
647 static int nHits = CbmL1::Instance()->algo->NStations;
648 int iVec = 0, i = 0;
649 int nTracks_SIMD = fvecLen;
650 L1TrackPar T; // fitting parametr coresponding to current track
651
652 CbmStsTrack* t[fvecLen];
653
654 int ista;
655 L1Station* sta = CbmL1::Instance()->algo->vStations;
656 fvec* x = new fvec[nHits];
657 fvec* u = new fvec[nHits];
658 fvec* v = new fvec[nHits];
659 fvec* y = new fvec[nHits];
660 fvec* z = new fvec[nHits];
661 fvec* w = new fvec[nHits];
662 fvec y_temp;
663 fvec fz0, fz1, fz2, dz, z_start, z_end;
664 L1FieldValue* fB = new L1FieldValue[nHits];
665 L1FieldValue fB_temp _fvecalignment;
666
667 unsigned short N_vTracks = Tracks.size();
668
669 for (unsigned short itrack = 0; itrack < N_vTracks; itrack++) {
670 Tracks[itrack].SetPidHypo(pidHypo);
671 }
672
673 for (unsigned short itrack = 0; itrack < N_vTracks; itrack += fvecLen) {
674 if (N_vTracks - itrack < static_cast<unsigned short>(fvecLen))
675 nTracks_SIMD = N_vTracks - itrack;
676
677 for (i = 0; i < nTracks_SIMD; i++) {
678 t[i] = &Tracks[itrack + i]; // current track
679 T.x[i] = t[i]->GetParamLast()->GetX();
680 T.y[i] = t[i]->GetParamLast()->GetY();
681 T.tx[i] = t[i]->GetParamLast()->GetTx();
682 T.ty[i] = t[i]->GetParamLast()->GetTy();
683 T.qp[i] = t[i]->GetParamLast()->GetQp();
684 T.z[i] = t[i]->GetParamLast()->GetZ();
685 T.C00[i] = t[i]->GetParamLast()->GetCovariance(0, 0);
686 T.C10[i] = t[i]->GetParamLast()->GetCovariance(1, 0);
687 T.C11[i] = t[i]->GetParamLast()->GetCovariance(1, 1);
688 T.C20[i] = t[i]->GetParamLast()->GetCovariance(2, 0);
689 T.C21[i] = t[i]->GetParamLast()->GetCovariance(2, 1);
690 T.C22[i] = t[i]->GetParamLast()->GetCovariance(2, 2);
691 T.C30[i] = t[i]->GetParamLast()->GetCovariance(3, 0);
692 T.C31[i] = t[i]->GetParamLast()->GetCovariance(3, 1);
693 T.C32[i] = t[i]->GetParamLast()->GetCovariance(3, 2);
694 T.C33[i] = t[i]->GetParamLast()->GetCovariance(3, 3);
695 T.C40[i] = t[i]->GetParamLast()->GetCovariance(4, 0);
696 T.C41[i] = t[i]->GetParamLast()->GetCovariance(4, 1);
697 T.C42[i] = t[i]->GetParamLast()->GetCovariance(4, 1);
698 T.C43[i] = t[i]->GetParamLast()->GetCovariance(4, 3);
699 T.C44[i] = t[i]->GetParamLast()->GetCovariance(4, 4);
700 }
701
702 // get hits of current track
703 for (i = 0; i < nHits; i++) {
704 w[i] = ZERO;
705 z[i] = sta[i].z;
706 }
707
708 for (iVec = 0; iVec < nTracks_SIMD; iVec++) {
709 int nHitsTrackMvd = t[iVec]->GetNMvdHits();
710 int nHitsTrackSts = t[iVec]->GetNStsHits();
711 int nHitsTrack = nHitsTrackMvd + nHitsTrackSts;
712 for (i = 0; i < nHitsTrack; i++) {
713 float posx = 0.f, posy = 0.f, posz = 0.f;
714
715 if (i < nHitsTrackMvd) {
716 if (!listMvdHits)
717 continue;
718 int hitIndex = t[iVec]->GetMvdHitIndex(i);
719 CbmMvdHit* hit = L1_DYNAMIC_CAST<CbmMvdHit*>(listMvdHits->At(hitIndex));
720
721 posx = hit->GetX();
722 posy = hit->GetY();
723 posz = hit->GetZ();
724 ista = posz < 7.f ? 0 : 1;
725 } else {
726 if (!listStsHits)
727 continue;
728 int hitIndex = t[iVec]->GetStsHitIndex(i - nHitsTrackMvd);
729 CbmStsHit* hit = L1_DYNAMIC_CAST<CbmStsHit*>(listStsHits->At(hitIndex));
730
731 posx = hit->GetX();
732 posy = hit->GetY();
733 posz = hit->GetZ();
734 ista = hit->GetStationNr() - 1 + NMvdStations;
735 }
736
737 w[ista][iVec] = 1.f;
738
739 x[ista][iVec] = posx;
740 y[ista][iVec] = posy;
741 u[ista][iVec] = posx * sta[ista].frontInfo.cos_phi[0] + posy * sta[ista].frontInfo.sin_phi[0];
742 v[ista][iVec] = posx * sta[ista].backInfo.cos_phi[0] + posy * sta[ista].backInfo.sin_phi[0];
743 z[ista][iVec] = posz;
744 sta[ista].fieldSlice.GetFieldValue(x[ista], y[ista], fB_temp);
745 fB[ista].x[iVec] = fB_temp.x[iVec];
746 fB[ista].y[iVec] = fB_temp.y[iVec];
747 fB[ista].z[iVec] = fB_temp.z[iVec];
748 if (i == 0)
749 z_start[iVec] = posz;
750 if (i == nHitsTrack - 1)
751 z_end[iVec] = posz;
752 }
753 }
754
755 // fit backward
756
757 fvec qp0 = T.qp;
758
759 i = nHits - 1;
760
761 FilterFirst(T, x[i], y[i], w[i], sta[i]);
762 AddMaterial(T, CbmL1::Instance()->algo->fRadThick[i].GetRadThick(T.x, T.y), qp0, mass2, w[i]);
763
764 fz1 = z[i];
765 sta[i].fieldSlice.GetFieldValue(T.x, T.y, fB1);
766 fB1.Combine(fB[i], w[i]);
767
768 fz2 = z[i - 2];
769 dz = fz2 - fz1;
770 sta[i].fieldSlice.GetFieldValue(T.x + T.tx * dz, T.y + T.ty * dz, fB2);
771 fB2.Combine(fB[i - 2], w[i - 2]);
772 fld.Set(fB2, fz2, fB1, fz1, fB0, fz0);
773 for (--i; i >= 0; i--) {
774 fz0 = z[i];
775 dz = (fz1 - fz0);
776 sta[i].fieldSlice.GetFieldValue(T.x - T.tx * dz, T.y - T.ty * dz, fB0);
777 fB0.Combine(fB[i], w[i]);
778 fld.Set(fB0, fz0, fB1, fz1, fB2, fz2);
779
780 // fvec one = ONE;
781 fvec zero = ZERO;
782 fvec initialised = fvec(z[i] > z_end) & fvec(z_start > z[i]);
783 fvec w1 = (zero & initialised) + (w[i] & (!initialised));
784
785 L1Extrapolate(T, z[i], qp0, fld, &w1);
786 if (i == NMvdStations - 1)
787 AddPipeMaterial(T, qp0, mass2, w1);
788 AddMaterial(T, CbmL1::Instance()->algo->fRadThick[i].GetRadThick(T.x, T.y), qp0, mass2, w1);
789 Filter(T, sta[i].frontInfo, u[i], w[i]);
790 Filter(T, sta[i].backInfo, v[i], w[i]);
791 fB2 = fB1;
792 fz2 = fz1;
793 fB1 = fB0;
794 fz1 = fz0;
795 }
796
797 for (iVec = 0; iVec < nTracks_SIMD; iVec++) {
798 t[iVec]->GetParamFirst()->SetX(T.x[iVec]);
799 t[iVec]->GetParamFirst()->SetY(T.y[iVec]);
800 t[iVec]->GetParamFirst()->SetTx(T.tx[iVec]);
801 t[iVec]->GetParamFirst()->SetTy(T.ty[iVec]);
802 t[iVec]->GetParamFirst()->SetQp(T.qp[iVec]);
803 // t[iVec]->TFirst[5] = z_start[iVec];
804 t[iVec]->GetParamFirst()->SetZ(T.z[iVec]);
805
806 t[iVec]->GetParamFirst()->SetCovariance(0, 0, T.C00[iVec]);
807 t[iVec]->GetParamFirst()->SetCovariance(1, 0, T.C10[iVec]);
808 t[iVec]->GetParamFirst()->SetCovariance(1, 1, T.C11[iVec]);
809 t[iVec]->GetParamFirst()->SetCovariance(2, 0, T.C20[iVec]);
810 t[iVec]->GetParamFirst()->SetCovariance(2, 1, T.C21[iVec]);
811 t[iVec]->GetParamFirst()->SetCovariance(2, 2, T.C22[iVec]);
812 t[iVec]->GetParamFirst()->SetCovariance(3, 0, T.C30[iVec]);
813 t[iVec]->GetParamFirst()->SetCovariance(3, 1, T.C31[iVec]);
814 t[iVec]->GetParamFirst()->SetCovariance(3, 2, T.C32[iVec]);
815 t[iVec]->GetParamFirst()->SetCovariance(3, 3, T.C33[iVec]);
816 t[iVec]->GetParamFirst()->SetCovariance(4, 0, T.C40[iVec]);
817 t[iVec]->GetParamFirst()->SetCovariance(4, 1, T.C41[iVec]);
818 t[iVec]->GetParamFirst()->SetCovariance(4, 2, T.C42[iVec]);
819 t[iVec]->GetParamFirst()->SetCovariance(4, 3, T.C43[iVec]);
820 t[iVec]->GetParamFirst()->SetCovariance(4, 4, T.C44[iVec]);
821
822 t[iVec]->SetChi2(T.chi2[iVec]);
823 t[iVec]->SetNDF(static_cast<int>(T.NDF[iVec]));
824 }
825
826 // fit forward
827
828 i = 0;
829 FilterLast(T, x[i], y[i], w[i], sta[i]);
830 qp0 = T.qp;
831 AddMaterial(T, CbmL1::Instance()->algo->fRadThick[i].GetRadThick(T.x, T.y), qp0, mass2, w[i]);
832
833 fz1 = z[i];
834 sta[i].fieldSlice.GetFieldValue(T.x, T.y, fB1);
835 fB1.Combine(fB[i], w[i]);
836
837 fz2 = z[i + 2];
838 dz = fz2 - fz1;
839 sta[i].fieldSlice.GetFieldValue(T.x + T.tx * dz, T.y + T.ty * dz, fB2);
840 fB2.Combine(fB[i + 2], w[i + 2]);
841 fld.Set(fB2, fz2, fB1, fz1, fB0, fz0);
842
843 for (++i; i < nHits; i++) {
844 fz0 = z[i];
845 dz = (fz1 - fz0);
846 sta[i].fieldSlice.GetFieldValue(T.x - T.tx * dz, T.y - T.ty * dz, fB0);
847 fB0.Combine(fB[i], w[i]);
848 fld.Set(fB0, fz0, fB1, fz1, fB2, fz2);
849
850 // fvec one = ONE;
851 fvec zero = ZERO;
852 fvec initialised = fvec(z[i] > z_end) & fvec(z_start > z[i]);
853 fvec w1 = (zero & initialised) + (w[i] & (!initialised));
854
855 L1Extrapolate(T, z[i], qp0, fld, &w1);
856 if (i == NMvdStations)
857 AddPipeMaterial(T, qp0, mass2, w1);
858 AddMaterial(T, CbmL1::Instance()->algo->fRadThick[i].GetRadThick(T.x, T.y), qp0, mass2, w1);
859 Filter(T, sta[i].frontInfo, u[i], w[i]);
860 Filter(T, sta[i].backInfo, v[i], w[i]);
861
862 fB2 = fB1;
863 fz2 = fz1;
864 fB1 = fB0;
865 fz1 = fz0;
866 }
867
868 for (iVec = 0; iVec < nTracks_SIMD; iVec++) {
869 t[iVec]->GetParamLast()->SetX(T.x[iVec]);
870 t[iVec]->GetParamLast()->SetY(T.y[iVec]);
871 t[iVec]->GetParamLast()->SetTx(T.tx[iVec]);
872 t[iVec]->GetParamLast()->SetTy(T.ty[iVec]);
873 t[iVec]->GetParamLast()->SetQp(T.qp[iVec]);
874 // t[iVec]->TFirst[5] = z_start[iVec];
875 t[iVec]->GetParamLast()->SetZ(T.z[iVec]);
876
877 t[iVec]->GetParamLast()->SetCovariance(0, 0, T.C00[iVec]);
878 t[iVec]->GetParamLast()->SetCovariance(1, 0, T.C10[iVec]);
879 t[iVec]->GetParamLast()->SetCovariance(1, 1, T.C11[iVec]);
880 t[iVec]->GetParamLast()->SetCovariance(2, 0, T.C20[iVec]);
881 t[iVec]->GetParamLast()->SetCovariance(2, 1, T.C21[iVec]);
882 t[iVec]->GetParamLast()->SetCovariance(2, 2, T.C22[iVec]);
883 t[iVec]->GetParamLast()->SetCovariance(3, 0, T.C30[iVec]);
884 t[iVec]->GetParamLast()->SetCovariance(3, 1, T.C31[iVec]);
885 t[iVec]->GetParamLast()->SetCovariance(3, 2, T.C32[iVec]);
886 t[iVec]->GetParamLast()->SetCovariance(3, 3, T.C33[iVec]);
887 t[iVec]->GetParamLast()->SetCovariance(4, 0, T.C40[iVec]);
888 t[iVec]->GetParamLast()->SetCovariance(4, 1, T.C41[iVec]);
889 t[iVec]->GetParamLast()->SetCovariance(4, 2, T.C42[iVec]);
890 t[iVec]->GetParamLast()->SetCovariance(4, 3, T.C43[iVec]);
891 t[iVec]->GetParamLast()->SetCovariance(4, 4, T.C44[iVec]);
892 }
893 }
894
895 delete[] x;
896 delete[] u;
897 delete[] v;
898 delete[] y;
899 delete[] z;
900 delete[] w;
901 delete[] fB;
902}
903
904void CbmL1PFFitter::GetChiToVertex(vector<CbmStsTrack>& Tracks,
905 vector<L1FieldRegion>& field,
906 vector<float>& chiToVtx,
907 CbmKFVertex& primVtx,
908 float chiPrim)
909{
910 chiToVtx.reserve(Tracks.size());
911
912 int nTracks_SIMD = fvecLen;
913 L1TrackPar T; // fitting parametr coresponding to current track
914
915 CbmStsTrack* t[fvecLen];
916
917 int nStations = CbmL1::Instance()->algo->NStations;
918 int NMvdStations = CbmL1::Instance()->algo->NMvdStations;
919 L1Station* sta = CbmL1::Instance()->algo->vStations;
920 fvec* zSta = new fvec[nStations];
921 for (int iSta = 0; iSta < nStations; iSta++)
922 zSta[iSta] = sta[iSta].z;
923
924 field.reserve(Tracks.size());
925
926 L1FieldRegion fld _fvecalignment;
927 L1FieldValue fB[3], fB_temp _fvecalignment;
928 fvec zField[3];
929
930 FairRootManager* fManger = FairRootManager::Instance();
931 TClonesArray* listStsHits = (TClonesArray*)fManger->GetObject("StsHit");
932 TClonesArray* listMvdHits = 0;
933 if (NMvdStations > 0.)
934 listMvdHits = (TClonesArray*)fManger->GetObject("MvdHit");
935
936 unsigned short N_vTracks = Tracks.size();
937 int ista;
938 for (unsigned short itrack = 0; itrack < N_vTracks; itrack += fvecLen) {
939 if (N_vTracks - itrack < static_cast<unsigned short>(fvecLen))
940 nTracks_SIMD = N_vTracks - itrack;
941
942 fvec mass2;
943 for (int iVec = 0; iVec < nTracks_SIMD; iVec++) {
944 t[iVec] = &Tracks[itrack + iVec]; // current track
945 T.x[iVec] = t[iVec]->GetParamFirst()->GetX();
946 T.y[iVec] = t[iVec]->GetParamFirst()->GetY();
947 T.tx[iVec] = t[iVec]->GetParamFirst()->GetTx();
948 T.ty[iVec] = t[iVec]->GetParamFirst()->GetTy();
949 T.qp[iVec] = t[iVec]->GetParamFirst()->GetQp();
950 T.z[iVec] = t[iVec]->GetParamFirst()->GetZ();
951 T.C00[iVec] = t[iVec]->GetParamFirst()->GetCovariance(0, 0);
952 T.C10[iVec] = t[iVec]->GetParamFirst()->GetCovariance(1, 0);
953 T.C11[iVec] = t[iVec]->GetParamFirst()->GetCovariance(1, 1);
954 T.C20[iVec] = t[iVec]->GetParamFirst()->GetCovariance(2, 0);
955 T.C21[iVec] = t[iVec]->GetParamFirst()->GetCovariance(2, 1);
956 T.C22[iVec] = t[iVec]->GetParamFirst()->GetCovariance(2, 2);
957 T.C30[iVec] = t[iVec]->GetParamFirst()->GetCovariance(3, 0);
958 T.C31[iVec] = t[iVec]->GetParamFirst()->GetCovariance(3, 1);
959 T.C32[iVec] = t[iVec]->GetParamFirst()->GetCovariance(3, 2);
960 T.C33[iVec] = t[iVec]->GetParamFirst()->GetCovariance(3, 3);
961 T.C40[iVec] = t[iVec]->GetParamFirst()->GetCovariance(4, 0);
962 T.C41[iVec] = t[iVec]->GetParamFirst()->GetCovariance(4, 1);
963 T.C42[iVec] = t[iVec]->GetParamFirst()->GetCovariance(4, 1);
964 T.C43[iVec] = t[iVec]->GetParamFirst()->GetCovariance(4, 3);
965 T.C44[iVec] = t[iVec]->GetParamFirst()->GetCovariance(4, 4);
966 // float mass = TDatabasePDG::Instance()->GetParticle(t[iVec]->GetPidHypo())->Mass();
967 const float mass = CbmKFParticleDatabase::Instance()->GetMass(t[iVec]->GetPidHypo());
968 mass2[iVec] = mass * mass;
969
970 int nHitsTrackMvd = t[iVec]->GetNMvdHits();
971 for (int iH = 0; iH < 2; iH++) {
972 float posx = 0.f, posy = 0.f, posz = 0.f;
973
974 if (iH < nHitsTrackMvd) {
975 if (!listMvdHits)
976 continue;
977 int hitIndex = t[iVec]->GetMvdHitIndex(iH);
978 CbmMvdHit* hit = L1_DYNAMIC_CAST<CbmMvdHit*>(listMvdHits->At(hitIndex));
979
980 posx = hit->GetX();
981 posy = hit->GetY();
982 posz = hit->GetZ();
983 ista = posz < 7.f ? 0 : 1;
984 } else {
985 if (!listStsHits)
986 continue;
987 int hitIndex = t[iVec]->GetStsHitIndex(iH - nHitsTrackMvd);
988 CbmStsHit* hit = L1_DYNAMIC_CAST<CbmStsHit*>(listStsHits->At(hitIndex));
989
990 posx = hit->GetX();
991 posy = hit->GetY();
992 posz = hit->GetZ();
993 ista = hit->GetStationNr() - 1 + NMvdStations;
994 }
995
996 sta[ista].fieldSlice.GetFieldValue(posx, posy, fB_temp);
997 fB[iH + 1].x[iVec] = fB_temp.x[iVec];
998 fB[iH + 1].y[iVec] = fB_temp.y[iVec];
999 fB[iH + 1].z[iVec] = fB_temp.z[iVec];
1000 zField[iH + 1][iVec] = posz;
1001 }
1002 }
1003
1004 fB[0] = CbmL1::Instance()->algo->GetvtxFieldValue();
1005 zField[0] = 0;
1006 fld.Set(fB[2], zField[2], fB[1], zField[1], fB[0], zField[0]);
1007 field.push_back(fld);
1008
1009 for (int iSt = nStations - 4; iSt >= 0; iSt--) {
1010 // fvec zero = ZERO;
1011 fvec w = ONE;
1012 fvec initialized = fvec(T.z > (zSta[iSt] + 2.5));
1013 w = fvec(w & initialized);
1014
1015 L1Extrapolate(T, zSta[iSt], T.qp, fld, &w);
1016 if (iSt == NMvdStations - 1)
1017 AddPipeMaterial(T, T.qp, mass2, w);
1018 AddMaterial(T, CbmL1::Instance()->algo->fRadThick[iSt].GetRadThick(T.x, T.y), T.qp, mass2, w);
1019 }
1020 fvec ONE = 1;
1021 if (NMvdStations <= 0)
1022 AddPipeMaterial(T, T.qp, mass2, ONE);
1023 L1Extrapolate(T, primVtx.GetRefZ(), T.qp, fld);
1024
1025 Double_t Cv[3] = {primVtx.GetCovMatrix()[0], primVtx.GetCovMatrix()[1], primVtx.GetCovMatrix()[2]};
1026
1027 fvec dx = T.x - primVtx.GetRefX();
1028 fvec dy = T.y - primVtx.GetRefY();
1029 fvec c[3] = {T.C00, T.C10, T.C11};
1030 c[0] += Cv[0];
1031 c[1] += Cv[1];
1032 c[2] += Cv[2];
1033 fvec d = c[0] * c[2] - c[1] * c[1];
1034 fvec chi = sqrt(fabs(0.5 * (dx * dx * c[0] - 2 * dx * dy * c[1] + dy * dy * c[2]) / d));
1035 fvec isNull = fvec(fabs(d) < 1.e-20);
1036 chi = fvec(fvec(!isNull) & chi) + fvec(isNull & fvec(0));
1037
1038 for (int iVec = 0; iVec < nTracks_SIMD; iVec++)
1039 chiToVtx.push_back(chi[iVec]);
1040
1041 if (chiPrim > 0) {
1042 for (int iVec = 0; iVec < nTracks_SIMD; iVec++) {
1043 if (chi[iVec] < chiPrim) {
1044 t[iVec]->GetParamFirst()->SetX(T.x[iVec]);
1045 t[iVec]->GetParamFirst()->SetY(T.y[iVec]);
1046 t[iVec]->GetParamFirst()->SetTx(T.tx[iVec]);
1047 t[iVec]->GetParamFirst()->SetTy(T.ty[iVec]);
1048 t[iVec]->GetParamFirst()->SetQp(T.qp[iVec]);
1049 t[iVec]->GetParamFirst()->SetZ(T.z[iVec]);
1050
1051 t[iVec]->GetParamFirst()->SetCovariance(0, 0, T.C00[iVec]);
1052 t[iVec]->GetParamFirst()->SetCovariance(1, 0, T.C10[iVec]);
1053 t[iVec]->GetParamFirst()->SetCovariance(1, 1, T.C11[iVec]);
1054 t[iVec]->GetParamFirst()->SetCovariance(2, 0, T.C20[iVec]);
1055 t[iVec]->GetParamFirst()->SetCovariance(2, 1, T.C21[iVec]);
1056 t[iVec]->GetParamFirst()->SetCovariance(2, 2, T.C22[iVec]);
1057 t[iVec]->GetParamFirst()->SetCovariance(3, 0, T.C30[iVec]);
1058 t[iVec]->GetParamFirst()->SetCovariance(3, 1, T.C31[iVec]);
1059 t[iVec]->GetParamFirst()->SetCovariance(3, 2, T.C32[iVec]);
1060 t[iVec]->GetParamFirst()->SetCovariance(3, 3, T.C33[iVec]);
1061 t[iVec]->GetParamFirst()->SetCovariance(4, 0, T.C40[iVec]);
1062 t[iVec]->GetParamFirst()->SetCovariance(4, 1, T.C41[iVec]);
1063 t[iVec]->GetParamFirst()->SetCovariance(4, 2, T.C42[iVec]);
1064 t[iVec]->GetParamFirst()->SetCovariance(4, 3, T.C43[iVec]);
1065 t[iVec]->GetParamFirst()->SetCovariance(4, 4, T.C44[iVec]);
1066 }
1067 }
1068 }
1069 }
1070 delete[] zSta;
1071}
1072
1073void CbmL1PFFitter::CalculateFieldRegion(vector<CbmStsTrack>& Tracks, vector<L1FieldRegion>& field)
1074{
1075 field.reserve(Tracks.size());
1076
1077 L1FieldRegion fld _fvecalignment;
1078
1079 FairRootManager* fManger = FairRootManager::Instance();
1080 TClonesArray* listStsHits = (TClonesArray*)fManger->GetObject("StsHit");
1081 TClonesArray* listMvdHits = 0;
1082 int NMvdStations = CbmL1::Instance()->algo->NMvdStations;
1083 if (NMvdStations > 0.)
1084 listMvdHits = (TClonesArray*)fManger->GetObject("MvdHit");
1085
1086 int nTracks_SIMD = fvecLen;
1087 L1TrackPar T; // fitting parametr coresponding to current track
1088
1089 CbmStsTrack* t[fvecLen];
1090
1091 int ista;
1092 L1Station* sta = CbmL1::Instance()->algo->vStations;
1093 L1FieldValue fB[3], fB_temp _fvecalignment;
1094 fvec zField[3];
1095
1096 unsigned short N_vTracks = Tracks.size();
1097
1098 for (unsigned short itrack = 0; itrack < N_vTracks; itrack += fvecLen) {
1099 if (N_vTracks - itrack < static_cast<unsigned short>(fvecLen))
1100 nTracks_SIMD = N_vTracks - itrack;
1101
1102 for (int i = 0; i < nTracks_SIMD; i++)
1103 t[i] = &Tracks[itrack + i]; // current track
1104
1105 for (int iVec = 0; iVec < nTracks_SIMD; iVec++) {
1106 int nHitsTrackMvd = t[iVec]->GetNMvdHits();
1107 for (int iH = 0; iH < 2; iH++) {
1108 float posx = 0.f, posy = 0.f, posz = 0.f;
1109
1110 if (iH < nHitsTrackMvd) {
1111 if (!listMvdHits)
1112 continue;
1113 int hitIndex = t[iVec]->GetMvdHitIndex(iH);
1114 CbmMvdHit* hit = L1_DYNAMIC_CAST<CbmMvdHit*>(listMvdHits->At(hitIndex));
1115
1116 posx = hit->GetX();
1117 posy = hit->GetY();
1118 posz = hit->GetZ();
1119 ista = posz < 7.f ? 0 : 1;
1120 } else {
1121 if (!listStsHits)
1122 continue;
1123 int hitIndex = t[iVec]->GetStsHitIndex(iH - nHitsTrackMvd);
1124 CbmStsHit* hit = L1_DYNAMIC_CAST<CbmStsHit*>(listStsHits->At(hitIndex));
1125
1126 posx = hit->GetX();
1127 posy = hit->GetY();
1128 posz = hit->GetZ();
1129 ista = hit->GetStationNr() - 1 + NMvdStations;
1130 }
1131
1132 sta[ista].fieldSlice.GetFieldValue(posx, posy, fB_temp);
1133 fB[iH + 1].x[iVec] = fB_temp.x[iVec];
1134 fB[iH + 1].y[iVec] = fB_temp.y[iVec];
1135 fB[iH + 1].z[iVec] = fB_temp.z[iVec];
1136 zField[iH + 1][iVec] = posz;
1137 }
1138 }
1139
1140 fB[0] = CbmL1::Instance()->algo->GetvtxFieldValue();
1141 zField[0] = 0;
1142 fld.Set(fB[2], zField[2], fB[1], zField[1], fB[0], zField[0]);
1143 field.push_back(fld);
1144 }
1145}
L1_DYNAMIC_CAST
#define L1_DYNAMIC_CAST
Definition CbmL1Def.h:26
L1Extrapolate
void L1Extrapolate(L1TrackPar &T, fvec z_out, fvec qp0, L1FieldRegion &F, fvec *w=0)
Definition L1Extrapolation.h:11
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
_fvecalignment
#define _fvecalignment
Definition P4_F32vec4.h:236
i
int i
Definition P4_F32vec4.h:22
fvec
F32vec4 fvec
Definition P4_F32vec4.h:231
Double_t
double Double_t
CbmKFParticleDatabase::Instance
static CbmKFParticleDatabase * Instance()
Definition CbmKFParticleDatabase.h:23
CbmKFVertex::GetRefZ
Double_t & GetRefZ()
Definition CbmKFVertex.h:39
CbmKFVertex::GetRefX
Double_t & GetRefX()
I/O parameters for vertex fit routines.
Definition CbmKFVertex.h:37
CbmKFVertex::GetCovMatrix
Double_t * GetCovMatrix()
Definition CbmKFVertex.h:40
CbmKFVertex::GetRefY
Double_t & GetRefY()
Definition CbmKFVertex.h:38
CbmL1HitStore::ExtIndex
int ExtIndex
Definition CbmL1.h:42
CbmL1PFFitter::Filter
void Filter(L1TrackPar &T, L1UMeasurementInfo &info, fvec &u, fvec &w)
Definition CbmL1PFFitter.cxx:164
CbmL1PFFitter::AddPipeMaterial
void AddPipeMaterial(L1TrackPar &T, fvec qp0, fvec &mass2, fvec &w)
Definition CbmL1PFFitter.cxx:78
CbmL1PFFitter::GetChiToVertex
void GetChiToVertex(std::vector< CbmL1Track > &Tracks, std::vector< float > &chiToVtx, CbmKFVertex &primVtx)
CbmL1PFFitter::AddMaterial
void AddMaterial(L1TrackPar &T, fvec radThick, fvec qp0, fvec &mass2, fvec &w)
Definition CbmL1PFFitter.cxx:51
CbmL1PFFitter::FilterLast
void FilterLast(L1TrackPar &track, fvec &x, fvec &y, fvec &w, L1Station &st)
Definition CbmL1PFFitter.cxx:135
CbmL1PFFitter::CalculateFieldRegion
void CalculateFieldRegion(std::vector< CbmL1Track > &Tracks, std::vector< L1FieldRegion > &Field)
CbmL1PFFitter::Fit
void Fit(std::vector< CbmL1Track > &Tracks, fvec mass=0.1395679f)
CbmL1PFFitter::FilterFirst
void FilterFirst(L1TrackPar &track, fvec &x, fvec &y, fvec &w, L1Station &st)
Definition CbmL1PFFitter.cxx:104
CbmL1Track::TLast
double TLast[6]
Definition CbmL1Track.h:57
CbmL1Track::StsHits
vector< int > StsHits
Definition CbmL1Track.h:58
CbmL1Track::CLast
double CLast[15]
Definition CbmL1Track.h:57
CbmL1Track::GetNOfHits
int GetNOfHits()
Number of Degrees of Freedom after fit.
Definition CbmL1Track.h:43
CbmL1::algo
L1Algo * algo
Definition CbmL1.h:54
CbmL1::Instance
static CbmL1 * Instance()
reconstructed tracks
Definition CbmL1.h:59
CbmL1::vHitStore
vector< CbmL1HitStore > vHitStore
Definition CbmL1.h:95
CbmL1::listMvdHits
TClonesArray * listMvdHits
Mvd MC points.
Definition CbmL1.h:168
CbmL1::listStsHits
TClonesArray * listStsHits
Definition CbmL1.h:165
CbmStsHit::GetStationNr
virtual Int_t GetStationNr() const
Definition CbmStsHit.h:67
CbmStsTrack::SetNDF
void SetNDF(Int_t ndf)
Definition CbmStsTrack.h:87
CbmStsTrack::GetParamLast
FairTrackParam * GetParamLast()
Definition CbmStsTrack.h:70
CbmStsTrack::GetNMvdHits
Int_t GetNMvdHits() const
Definition CbmStsTrack.h:61
CbmStsTrack::GetMvdHitIndex
Int_t GetMvdHitIndex(Int_t iHit) const
Definition CbmStsTrack.h:63
CbmStsTrack::SetChi2
void SetChi2(Double_t chi2)
Definition CbmStsTrack.h:86
CbmStsTrack::GetNStsHits
Int_t GetNStsHits() const
Accessors.
Definition CbmStsTrack.h:60
CbmStsTrack::GetStsHitIndex
Int_t GetStsHitIndex(Int_t iHit) const
Definition CbmStsTrack.h:62
CbmStsTrack::GetParamFirst
FairTrackParam * GetParamFirst()
Definition CbmStsTrack.h:69
L1Algo::NStations
int NStations
Definition L1Algo.h:153
L1Algo::NMvdStations
int NMvdStations
Definition L1Algo.h:154
L1Algo::GetvtxFieldValue
L1FieldValue GetvtxFieldValue()
Definition L1Algo.h:189
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
L1FieldValue::Combine
void Combine(L1FieldValue &B, fvec w)
Definition L1Field.h:22
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::frontInfo
L1UMeasurementInfo frontInfo
Definition L1Station.h:30
TDatabasePDG::Instance
static TDatabasePDG * Instance()
TDatabasePDG::GetParticle
TParticlePDG * GetParticle(const char *name) const
TObjArray::At
TObject * At(Int_t idx) const override
TParticlePDG::Mass
Double_t Mass() const
log
RVec< PromoteType< T > > log(const RVec< T > &v)
sqrt
RVec< PromoteType< T > > sqrt(const RVec< T > &v)
NS_L1TrackFitter
Definition L1TrackFitter.cxx:28
NS_L1TrackFitter::ZERO
const fvec ZERO
Definition L1TrackFitter.cxx:30
NS_L1TrackFitter::vINF
const fvec vINF
Definition L1TrackFitter.cxx:32
NS_L1TrackFitter::ONE
const fvec ONE
Definition L1TrackFitter.cxx:31
NS_L1TrackFitter::c_light_i
const fvec c_light_i
Definition L1TrackFitter.cxx:29
NS_L1TrackFitter::c_light
const fvec c_light
Definition L1TrackFitter.cxx:29
ROOT::Math::Chebyshev::T
double T(double x)
CbmL1TrackPar::C
double C[15]
Definition CbmL1TrackPar.h:22