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