CbmRoot
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 "CbmStsAddress.h"
20 #include "CbmStsSetup.h"
21 #include "CbmStsTrack.h"
22 #include "TClonesArray.h"
23 
24 //L1Algo tools
25 #include "CbmL1Track.h"
26 #include "L1AddMaterial.h"
27 #include "L1Algo.h"
28 #include "L1Extrapolation.h"
29 #include "L1Filtration.h"
30 #include "L1MaterialInfo.h"
31 #include "L1Station.h"
32 #include "L1TrackPar.h"
33 
34 #include "FairRootManager.h"
35 #include "TDatabasePDG.h"
36 
37 #include "CbmKFVertex.h"
38 #include "KFParticleDatabase.h"
39 
40 using std::vector;
41 
42 
43 namespace NS_L1TrackFitter {
44  const fvec c_light = 0.000299792458, c_light_i = 1. / c_light;
45  const fvec ZERO = 0.;
46  const fvec ONE = 1.;
47  const fvec vINF = 0.1;
48 } // namespace NS_L1TrackFitter
49 using namespace NS_L1TrackFitter;
50 
51 CbmL1PFFitter::CbmL1PFFitter() {}
52 
53 CbmL1PFFitter::~CbmL1PFFitter() {}
54 
55 void CbmL1PFFitter::FilterFirst(L1TrackPar& track,
56  fvec& x,
57  fvec& y,
58  L1Station& st) {
59  track.C00 = st.XYInfo.C00;
60  track.C10 = st.XYInfo.C10;
61  track.C11 = st.XYInfo.C11;
62  track.C20 = ZERO;
63  track.C21 = ZERO;
64  track.C22 = vINF;
65  track.C30 = ZERO;
66  track.C31 = ZERO;
67  track.C32 = ZERO;
68  track.C33 = vINF;
69  track.C40 = ZERO;
70  track.C41 = ZERO;
71  track.C42 = ZERO;
72  track.C43 = ZERO;
73  track.C44 = ONE;
74 
75  track.x = x;
76  track.y = y;
77  track.NDF = -3.0;
78  track.chi2 = ZERO;
79 }
80 
81 void CbmL1PFFitter::Fit(vector<CbmStsTrack>& Tracks, vector<int>& pidHypo) {
82 
83  L1FieldValue fB0, fB1, fB2 _fvecalignment;
84  L1FieldRegion fld _fvecalignment;
85 
86 
87  FairRootManager* fManger = FairRootManager::Instance();
88  TClonesArray* listStsHits = (TClonesArray*) fManger->GetObject("StsHit");
89  int NMvdStations = CbmL1::Instance()->algo->NMvdStations;
90  TClonesArray* listMvdHits = 0;
91  if (NMvdStations > 0.)
92  listMvdHits = (TClonesArray*) fManger->GetObject("MvdHit");
93 
94  static int nHits = CbmL1::Instance()->algo->NStations;
95  int iVec = 0, i = 0;
96  int nTracks_SIMD = fvecLen;
97  L1TrackPar T; // fitting parametr coresponding to current track
98 
99  CbmStsTrack* t[fvecLen];
100 
101  int ista;
102  L1Station* sta = CbmL1::Instance()->algo->vStations;
103  L1Station staFirst, staLast;
104  fvec* x = new fvec[nHits];
105  fvec* u = new fvec[nHits];
106  fvec* v = new fvec[nHits];
107  fvec* y = new fvec[nHits];
108  fvec* z = new fvec[nHits];
109  fvec* w = new fvec[nHits];
110  fvec y_temp;
111  fvec x_first, y_first, x_last, y_last;
112  fvec fz0, fz1, fz2, dz, z_start, z_end;
113  L1FieldValue* fB = new L1FieldValue[nHits];
114  L1FieldValue fB_temp _fvecalignment;
115 
116  unsigned short N_vTracks = Tracks.size();
117 
118 
119  for (unsigned short itrack = 0; itrack < N_vTracks; itrack++) {
120  Tracks[itrack].SetPidHypo(pidHypo[itrack]);
121  }
122 
123  fvec mass, mass2;
124 
125  for (unsigned short itrack = 0; itrack < N_vTracks; itrack += fvecLen) {
126 
127  if (N_vTracks - itrack < static_cast<unsigned short>(fvecLen))
128  nTracks_SIMD = N_vTracks - itrack;
129  for (i = 0; i < nTracks_SIMD; i++) {
130  t[i] = &Tracks[itrack + i]; // current track
131  T.x[i] = t[i]->GetParamFirst()->GetX();
132  T.y[i] = t[i]->GetParamFirst()->GetY();
133  T.tx[i] = t[i]->GetParamFirst()->GetTx();
134  T.ty[i] = t[i]->GetParamFirst()->GetTy();
135  T.qp[i] = t[i]->GetParamFirst()->GetQp();
136  T.z[i] = t[i]->GetParamFirst()->GetZ();
137  T.C00[i] = t[i]->GetParamFirst()->GetCovariance(0, 0);
138  T.C10[i] = t[i]->GetParamFirst()->GetCovariance(1, 0);
139  T.C11[i] = t[i]->GetParamFirst()->GetCovariance(1, 1);
140  T.C20[i] = t[i]->GetParamFirst()->GetCovariance(2, 0);
141  T.C21[i] = t[i]->GetParamFirst()->GetCovariance(2, 1);
142  T.C22[i] = t[i]->GetParamFirst()->GetCovariance(2, 2);
143  T.C30[i] = t[i]->GetParamFirst()->GetCovariance(3, 0);
144  T.C31[i] = t[i]->GetParamFirst()->GetCovariance(3, 1);
145  T.C32[i] = t[i]->GetParamFirst()->GetCovariance(3, 2);
146  T.C33[i] = t[i]->GetParamFirst()->GetCovariance(3, 3);
147  T.C40[i] = t[i]->GetParamFirst()->GetCovariance(4, 0);
148  T.C41[i] = t[i]->GetParamFirst()->GetCovariance(4, 1);
149  T.C42[i] = t[i]->GetParamFirst()->GetCovariance(4, 1);
150  T.C43[i] = t[i]->GetParamFirst()->GetCovariance(4, 3);
151  T.C44[i] = t[i]->GetParamFirst()->GetCovariance(4, 4);
152 
153  int pid = pidHypo[itrack + i];
154  if (pid == -1) pid = 211;
155  // mass[i] = TDatabasePDG::Instance()->GetParticle(pid)->Mass();
156  mass[i] = KFParticleDatabase::Instance()->GetMass(pid);
157  }
158  mass2 = mass * mass;
159  // get hits of current track
160  for (i = 0; i < nHits; i++) {
161  w[i] = ZERO;
162  z[i] = sta[i].z;
163  }
164 
165  for (iVec = 0; iVec < nTracks_SIMD; iVec++) {
166  int nHitsTrackMvd = t[iVec]->GetNofMvdHits();
167  int nHitsTrackSts = t[iVec]->GetNofStsHits();
168  int nHitsTrack = nHitsTrackMvd + nHitsTrackSts;
169  for (i = 0; i < nHitsTrack; i++) {
170  float posx = 0.f, posy = 0.f, posz = 0.f;
171 
172  if (i < nHitsTrackMvd) {
173  if (!listMvdHits) continue;
174  int hitIndex = t[iVec]->GetMvdHitIndex(i);
175  CbmMvdHit* hit =
176  L1_DYNAMIC_CAST<CbmMvdHit*>(listMvdHits->At(hitIndex));
177 
178  posx = hit->GetX();
179  posy = hit->GetY();
180  posz = hit->GetZ();
181  ista = hit->GetStationNr();
182  } else {
183  if (!listStsHits) continue;
184  int hitIndex = t[iVec]->GetHitIndex(i - nHitsTrackMvd);
185  CbmStsHit* hit =
186  L1_DYNAMIC_CAST<CbmStsHit*>(listStsHits->At(hitIndex));
187 
188  posx = hit->GetX();
189  posy = hit->GetY();
190  posz = hit->GetZ();
191  ista = CbmStsSetup::Instance()->GetStationNumber(hit->GetAddress())
192  + NMvdStations; //hit->GetStationNr() - 1 + NMvdStations;
193  }
194  w[ista][iVec] = 1.f;
195 
196  x[ista][iVec] = posx;
197  y[ista][iVec] = posy;
198  u[ista][iVec] = posx * sta[ista].frontInfo.cos_phi[0]
199  + posy * sta[ista].frontInfo.sin_phi[0];
200  v[ista][iVec] = posx * sta[ista].backInfo.cos_phi[0]
201  + posy * sta[ista].backInfo.sin_phi[0];
202  z[ista][iVec] = posz;
203  sta[ista].fieldSlice.GetFieldValue(x[ista], y[ista], fB_temp);
204  fB[ista].x[iVec] = fB_temp.x[iVec];
205  fB[ista].y[iVec] = fB_temp.y[iVec];
206  fB[ista].z[iVec] = fB_temp.z[iVec];
207  if (i == 0) {
208  z_start[iVec] = posz;
209  x_first[iVec] = x[ista][iVec];
210  y_first[iVec] = y[ista][iVec];
211  staFirst.XYInfo.C00[iVec] = sta[ista].XYInfo.C00[iVec];
212  staFirst.XYInfo.C10[iVec] = sta[ista].XYInfo.C10[iVec];
213  staFirst.XYInfo.C11[iVec] = sta[ista].XYInfo.C11[iVec];
214  }
215  if (i == nHitsTrack - 1) {
216  z_end[iVec] = posz;
217  x_last[iVec] = x[ista][iVec];
218  y_last[iVec] = y[ista][iVec];
219  staLast.XYInfo.C00[iVec] = sta[ista].XYInfo.C00[iVec];
220  staLast.XYInfo.C10[iVec] = sta[ista].XYInfo.C10[iVec];
221  staLast.XYInfo.C11[iVec] = sta[ista].XYInfo.C11[iVec];
222  }
223  }
224  }
225  // fit forward
226  i = 0;
227  FilterFirst(T, x_first, y_first, staFirst);
228  fvec qp0 = T.qp;
229  fz1 = z[i];
230  sta[i].fieldSlice.GetFieldValue(T.x, T.y, fB1);
231  fB1.Combine(fB[i], w[i]);
232  fz2 = z[i + 2];
233  dz = fz2 - fz1;
234  sta[i].fieldSlice.GetFieldValue(T.x + T.tx * dz, T.y + T.ty * dz, fB2);
235  fB2.Combine(fB[i + 2], w[i + 2]);
236  fld.Set(fB2, fz2, fB1, fz1, fB0, fz0);
237  for (++i; i < nHits; i++) {
238  fz0 = z[i];
239  dz = (fz1 - fz0);
240  sta[i].fieldSlice.GetFieldValue(T.x - T.tx * dz, T.y - T.ty * dz, fB0);
241  fB0.Combine(fB[i], w[i]);
242  fld.Set(fB0, fz0, fB1, fz1, fB2, fz2);
243 
244  fvec initialised = fvec(z[i] <= z_end) & fvec(z_start < z[i]);
245  fvec w1 = (w[i] & (initialised));
246  fvec wIn = (ONE & (initialised));
247 
248  L1Extrapolate(T, z[i], qp0, fld, &w1);
249  if (i == NMvdStations) {
250  L1AddPipeMaterial(T, qp0, wIn);
251  EnergyLossCorrection(T, mass2, PipeRadThick, qp0, fvec(-1.f), wIn);
252  }
253 #ifdef USE_RL_TABLE
254  L1AddMaterial(T,
255  CbmL1::Instance()->algo->fRadThick[i].GetRadThick(T.x, T.y),
256  qp0,
257  wIn,
258  mass2);
259  EnergyLossCorrection(
260  T,
261  mass2,
262  CbmL1::Instance()->algo->fRadThick[i].GetRadThick(T.x, T.y),
263  qp0,
264  -1,
265  wIn);
266 #else
267  L1AddMaterial(T, sta[i].materialInfo, qp0, wIn, mass2);
268 #endif
269  L1Filter(T, sta[i].frontInfo, u[i], w1);
270  L1Filter(T, sta[i].backInfo, v[i], w1);
271 
272  fB2 = fB1;
273  fz2 = fz1;
274  fB1 = fB0;
275  fz1 = fz0;
276  }
277 
278  L1TrackPar Tout = T;
279  for (iVec = 0; iVec < nTracks_SIMD; iVec++) {
280  FairTrackParam par;
281  par.SetX(T.x[iVec]);
282  par.SetY(T.y[iVec]);
283  par.SetTx(T.tx[iVec]);
284  par.SetTy(T.ty[iVec]);
285  par.SetQp(T.qp[iVec]);
286  par.SetZ(T.z[iVec]);
287 
288  par.SetCovariance(0, 0, Tout.C00[iVec]);
289  par.SetCovariance(1, 0, Tout.C10[iVec]);
290  par.SetCovariance(1, 1, Tout.C11[iVec]);
291  par.SetCovariance(2, 0, Tout.C20[iVec]);
292  par.SetCovariance(2, 1, Tout.C21[iVec]);
293  par.SetCovariance(2, 2, Tout.C22[iVec]);
294  par.SetCovariance(3, 0, Tout.C30[iVec]);
295  par.SetCovariance(3, 1, Tout.C31[iVec]);
296  par.SetCovariance(3, 2, Tout.C32[iVec]);
297  par.SetCovariance(3, 3, Tout.C33[iVec]);
298  par.SetCovariance(4, 0, Tout.C40[iVec]);
299  par.SetCovariance(4, 1, Tout.C41[iVec]);
300  par.SetCovariance(4, 2, Tout.C42[iVec]);
301  par.SetCovariance(4, 3, Tout.C43[iVec]);
302  par.SetCovariance(4, 4, Tout.C44[iVec]);
303  t[iVec]->SetParamLast(&par);
304  }
305 
306  //fit backward
307  qp0 = T.qp;
308 
309  i = nHits - 1;
310 
311  FilterFirst(T, x_last, y_last, staLast);
312 
313  fz1 = z[i];
314  sta[i].fieldSlice.GetFieldValue(T.x, T.y, fB1);
315  fB1.Combine(fB[i], w[i]);
316 
317  fz2 = z[i - 2];
318  dz = fz2 - fz1;
319  sta[i].fieldSlice.GetFieldValue(T.x + T.tx * dz, T.y + T.ty * dz, fB2);
320  fB2.Combine(fB[i - 2], w[i - 2]);
321  fld.Set(fB2, fz2, fB1, fz1, fB0, fz0);
322  for (--i; i >= 0; i--) {
323  fz0 = z[i];
324  dz = (fz1 - fz0);
325  sta[i].fieldSlice.GetFieldValue(T.x - T.tx * dz, T.y - T.ty * dz, fB0);
326  fB0.Combine(fB[i], w[i]);
327  fld.Set(fB0, fz0, fB1, fz1, fB2, fz2);
328 
329  fvec initialised = fvec(z[i] < z_end) & fvec(z_start <= z[i]);
330  fvec w1 = (w[i] & (initialised));
331  fvec wIn = (ONE & (initialised));
332 
333  L1Extrapolate(T, z[i], qp0, fld, &w1);
334  if (i == NMvdStations - 1) {
335  L1AddPipeMaterial(T, qp0, wIn);
336  EnergyLossCorrection(T, mass2, PipeRadThick, qp0, fvec(1.f), wIn);
337  }
338 #ifdef USE_RL_TABLE
339  L1AddMaterial(T,
340  CbmL1::Instance()->algo->fRadThick[i].GetRadThick(T.x, T.y),
341  qp0,
342  wIn,
343  mass2);
344  EnergyLossCorrection(
345  T,
346  mass2,
347  CbmL1::Instance()->algo->fRadThick[i].GetRadThick(T.x, T.y),
348  qp0,
349  1,
350  wIn);
351 #else
352  L1AddMaterial(T, sta[i].materialInfo, qp0, wIn, mass2);
353 #endif
354  L1Filter(T, sta[i].frontInfo, u[i], w1);
355  L1Filter(T, sta[i].backInfo, v[i], w1);
356 
357  fB2 = fB1;
358  fz2 = fz1;
359  fB1 = fB0;
360  fz1 = fz0;
361  }
362 
363  for (iVec = 0; iVec < nTracks_SIMD; iVec++) {
364  FairTrackParam par;
365  par.SetX(T.x[iVec]);
366  par.SetY(T.y[iVec]);
367  par.SetTx(T.tx[iVec]);
368  par.SetTy(T.ty[iVec]);
369  par.SetQp(T.qp[iVec]);
370  par.SetZ(T.z[iVec]);
371 
372  par.SetCovariance(0, 0, T.C00[iVec]);
373  par.SetCovariance(1, 0, T.C10[iVec]);
374  par.SetCovariance(1, 1, T.C11[iVec]);
375  par.SetCovariance(2, 0, T.C20[iVec]);
376  par.SetCovariance(2, 1, T.C21[iVec]);
377  par.SetCovariance(2, 2, T.C22[iVec]);
378  par.SetCovariance(3, 0, T.C30[iVec]);
379  par.SetCovariance(3, 1, T.C31[iVec]);
380  par.SetCovariance(3, 2, T.C32[iVec]);
381  par.SetCovariance(3, 3, T.C33[iVec]);
382  par.SetCovariance(4, 0, T.C40[iVec]);
383  par.SetCovariance(4, 1, T.C41[iVec]);
384  par.SetCovariance(4, 2, T.C42[iVec]);
385  par.SetCovariance(4, 3, T.C43[iVec]);
386  par.SetCovariance(4, 4, T.C44[iVec]);
387  t[iVec]->SetParamFirst(&par);
388 
389  t[iVec]->SetChiSq(T.chi2[iVec]);
390  t[iVec]->SetNDF(static_cast<int>(T.NDF[iVec]));
391  }
392  }
393 
394  delete[] x;
395  delete[] u;
396  delete[] v;
397  delete[] y;
398  delete[] z;
399  delete[] w;
400  delete[] fB;
401 }
402 
403 void CbmL1PFFitter::GetChiToVertex(vector<CbmStsTrack>& Tracks,
404  vector<L1FieldRegion>& field,
405  vector<float>& chiToVtx,
406  CbmKFVertex& primVtx,
407  float chiPrim) {
408  chiToVtx.reserve(Tracks.size());
409 
410  int nTracks_SIMD = fvecLen;
411  L1TrackPar T; // fitting parametr coresponding to current track
412 
413  CbmStsTrack* t[fvecLen];
414 
415  int nStations = CbmL1::Instance()->algo->NStations;
416  int NMvdStations = CbmL1::Instance()->algo->NMvdStations;
417  L1Station* sta = CbmL1::Instance()->algo->vStations;
418  fvec* zSta = new fvec[nStations];
419  for (int iSta = 0; iSta < nStations; iSta++)
420  zSta[iSta] = sta[iSta].z;
421 
422  field.reserve(Tracks.size());
423 
424  L1FieldRegion fld _fvecalignment;
425  L1FieldValue fB[3], fB_temp _fvecalignment;
426  fvec zField[3];
427 
428  FairRootManager* fManger = FairRootManager::Instance();
429  TClonesArray* listStsHits = (TClonesArray*) fManger->GetObject("StsHit");
430  TClonesArray* listMvdHits = 0;
431  if (NMvdStations > 0.)
432  listMvdHits = (TClonesArray*) fManger->GetObject("MvdHit");
433 
434  unsigned short N_vTracks = Tracks.size();
435  int ista;
436  for (unsigned short itrack = 0; itrack < N_vTracks; itrack += fvecLen) {
437  if (N_vTracks - itrack < static_cast<unsigned short>(fvecLen))
438  nTracks_SIMD = N_vTracks - itrack;
439 
440  fvec mass2;
441  for (int iVec = 0; iVec < nTracks_SIMD; iVec++) {
442  t[iVec] = &Tracks[itrack + iVec]; // current track
443  T.x[iVec] = t[iVec]->GetParamFirst()->GetX();
444  T.y[iVec] = t[iVec]->GetParamFirst()->GetY();
445  T.tx[iVec] = t[iVec]->GetParamFirst()->GetTx();
446  T.ty[iVec] = t[iVec]->GetParamFirst()->GetTy();
447  T.qp[iVec] = t[iVec]->GetParamFirst()->GetQp();
448  T.z[iVec] = t[iVec]->GetParamFirst()->GetZ();
449  T.C00[iVec] = t[iVec]->GetParamFirst()->GetCovariance(0, 0);
450  T.C10[iVec] = t[iVec]->GetParamFirst()->GetCovariance(1, 0);
451  T.C11[iVec] = t[iVec]->GetParamFirst()->GetCovariance(1, 1);
452  T.C20[iVec] = t[iVec]->GetParamFirst()->GetCovariance(2, 0);
453  T.C21[iVec] = t[iVec]->GetParamFirst()->GetCovariance(2, 1);
454  T.C22[iVec] = t[iVec]->GetParamFirst()->GetCovariance(2, 2);
455  T.C30[iVec] = t[iVec]->GetParamFirst()->GetCovariance(3, 0);
456  T.C31[iVec] = t[iVec]->GetParamFirst()->GetCovariance(3, 1);
457  T.C32[iVec] = t[iVec]->GetParamFirst()->GetCovariance(3, 2);
458  T.C33[iVec] = t[iVec]->GetParamFirst()->GetCovariance(3, 3);
459  T.C40[iVec] = t[iVec]->GetParamFirst()->GetCovariance(4, 0);
460  T.C41[iVec] = t[iVec]->GetParamFirst()->GetCovariance(4, 1);
461  T.C42[iVec] = t[iVec]->GetParamFirst()->GetCovariance(4, 1);
462  T.C43[iVec] = t[iVec]->GetParamFirst()->GetCovariance(4, 3);
463  T.C44[iVec] = t[iVec]->GetParamFirst()->GetCovariance(4, 4);
464  // float mass = TDatabasePDG::Instance()->GetParticle(t[iVec]->GetPidHypo())->Mass();
465  const float mass =
466  KFParticleDatabase::Instance()->GetMass(t[iVec]->GetPidHypo());
467  mass2[iVec] = mass * mass;
468 
469  int nHitsTrackMvd = t[iVec]->GetNofMvdHits();
470  for (int iH = 0; iH < 2; iH++) {
471  float posx = 0.f, posy = 0.f, posz = 0.f;
472 
473  if (iH < nHitsTrackMvd) {
474  if (!listMvdHits) continue;
475  int hitIndex = t[iVec]->GetMvdHitIndex(iH);
476  CbmMvdHit* hit =
477  L1_DYNAMIC_CAST<CbmMvdHit*>(listMvdHits->At(hitIndex));
478 
479  posx = hit->GetX();
480  posy = hit->GetY();
481  posz = hit->GetZ();
482  ista = hit->GetStationNr();
483  } else {
484  if (!listStsHits) continue;
485  int hitIndex = t[iVec]->GetHitIndex(iH - nHitsTrackMvd);
486  CbmStsHit* hit =
487  L1_DYNAMIC_CAST<CbmStsHit*>(listStsHits->At(hitIndex));
488 
489  posx = hit->GetX();
490  posy = hit->GetY();
491  posz = hit->GetZ();
492  ista = CbmStsSetup::Instance()->GetStationNumber(hit->GetAddress())
493  + NMvdStations; //hit->GetStationNr()-1+NMvdStations;
494  }
495 
496  sta[ista].fieldSlice.GetFieldValue(posx, posy, fB_temp);
497  fB[iH + 1].x[iVec] = fB_temp.x[iVec];
498  fB[iH + 1].y[iVec] = fB_temp.y[iVec];
499  fB[iH + 1].z[iVec] = fB_temp.z[iVec];
500  zField[iH + 1][iVec] = posz;
501  }
502  }
503 
504  fB[0] = CbmL1::Instance()->algo->GetVtxFieldValue();
505  zField[0] = 0;
506  fld.Set(fB[2], zField[2], fB[1], zField[1], fB[0], zField[0]);
507  field.push_back(fld);
508 
509  for (int iSt = nStations - 4; iSt >= 0; iSt--) {
510  fvec w = ONE;
511  fvec initialized = fvec(T.z > (zSta[iSt] + 2.5));
512  w = fvec(w & initialized);
513 
514  L1Extrapolate(T, zSta[iSt], T.qp, fld, &w);
515  if (iSt == NMvdStations - 1) {
516  L1AddPipeMaterial(T, T.qp, w, mass2);
517  EnergyLossCorrection(T, mass2, PipeRadThick, T.qp, fvec(1.f), w);
518  }
519  L1AddMaterial(
520  T,
521  CbmL1::Instance()->algo->fRadThick[iSt].GetRadThick(T.x, T.y),
522  T.qp,
523  w,
524  mass2);
525  EnergyLossCorrection(
526  T,
527  mass2,
528  CbmL1::Instance()->algo->fRadThick[iSt].GetRadThick(T.x, T.y),
529  T.qp,
530  fvec(1.f),
531  w);
532  }
533  if (NMvdStations <= 0) {
534  L1AddPipeMaterial(T, T.qp, ONE, mass2);
535  EnergyLossCorrection(T, mass2, PipeRadThick, T.qp, fvec(1.f), ONE);
536  }
537  L1Extrapolate(T, primVtx.GetRefZ(), T.qp, fld);
538  L1AddTargetMaterial(T, T.qp);
539 
540  Double_t Cv[3] = {primVtx.GetCovMatrix()[0],
541  primVtx.GetCovMatrix()[1],
542  primVtx.GetCovMatrix()[2]};
543 
544  fvec dx = T.x - primVtx.GetRefX();
545  fvec dy = T.y - primVtx.GetRefY();
546  fvec c[3] = {T.C00, T.C10, T.C11};
547  c[0] += Cv[0];
548  c[1] += Cv[1];
549  c[2] += Cv[2];
550  fvec d = c[0] * c[2] - c[1] * c[1];
551  fvec chi = sqrt(
552  fabs(0.5 * (dx * dx * c[0] - 2 * dx * dy * c[1] + dy * dy * c[2]) / d));
553  fvec isNull = fvec(fabs(d) < 1.e-20);
554  chi = fvec(fvec(!isNull) & chi) + fvec(isNull & fvec(0));
555 
556  for (int iVec = 0; iVec < nTracks_SIMD; iVec++)
557  chiToVtx.push_back(chi[iVec]);
558 
559  if (chiPrim > 0) {
560  for (int iVec = 0; iVec < nTracks_SIMD; iVec++) {
561  if (chi[iVec] < chiPrim) {
562  FairTrackParam par;
563  par.SetX(T.x[iVec]);
564  par.SetY(T.y[iVec]);
565  par.SetTx(T.tx[iVec]);
566  par.SetTy(T.ty[iVec]);
567  par.SetQp(T.qp[iVec]);
568  par.SetZ(T.z[iVec]);
569 
570  par.SetCovariance(0, 0, T.C00[iVec]);
571  par.SetCovariance(1, 0, T.C10[iVec]);
572  par.SetCovariance(1, 1, T.C11[iVec]);
573  par.SetCovariance(2, 0, T.C20[iVec]);
574  par.SetCovariance(2, 1, T.C21[iVec]);
575  par.SetCovariance(2, 2, T.C22[iVec]);
576  par.SetCovariance(3, 0, T.C30[iVec]);
577  par.SetCovariance(3, 1, T.C31[iVec]);
578  par.SetCovariance(3, 2, T.C32[iVec]);
579  par.SetCovariance(3, 3, T.C33[iVec]);
580  par.SetCovariance(4, 0, T.C40[iVec]);
581  par.SetCovariance(4, 1, T.C41[iVec]);
582  par.SetCovariance(4, 2, T.C42[iVec]);
583  par.SetCovariance(4, 3, T.C43[iVec]);
584  par.SetCovariance(4, 4, T.C44[iVec]);
585  t[iVec]->SetParamFirst(&par);
586  }
587  }
588  }
589  }
590  delete[] zSta;
591 }
592 
593 void CbmL1PFFitter::CalculateFieldRegion(vector<CbmStsTrack>& Tracks,
594  vector<L1FieldRegion>& field) {
595  field.reserve(Tracks.size());
596 
597  L1FieldRegion fld _fvecalignment;
598 
599  FairRootManager* fManger = FairRootManager::Instance();
600  TClonesArray* listStsHits = (TClonesArray*) fManger->GetObject("StsHit");
601  TClonesArray* listMvdHits = 0;
602  int NMvdStations = CbmL1::Instance()->algo->NMvdStations;
603  if (NMvdStations > 0.)
604  listMvdHits = (TClonesArray*) fManger->GetObject("MvdHit");
605 
606  int nTracks_SIMD = fvecLen;
607  L1TrackPar T; // fitting parametr coresponding to current track
608 
609  CbmStsTrack* t[fvecLen];
610 
611  int ista;
612  L1Station* sta = CbmL1::Instance()->algo->vStations;
613  L1FieldValue fB[3], fB_temp _fvecalignment;
614  fvec zField[3];
615 
616  unsigned short N_vTracks = Tracks.size();
617 
618  for (unsigned short itrack = 0; itrack < N_vTracks; itrack += fvecLen) {
619  if (N_vTracks - itrack < static_cast<unsigned short>(fvecLen))
620  nTracks_SIMD = N_vTracks - itrack;
621 
622  for (int i = 0; i < nTracks_SIMD; i++)
623  t[i] = &Tracks[itrack + i]; // current track
624 
625  for (int iVec = 0; iVec < nTracks_SIMD; iVec++) {
626  int nHitsTrackMvd = t[iVec]->GetNofMvdHits();
627  for (int iH = 0; iH < 2; iH++) {
628  float posx = 0.f, posy = 0.f, posz = 0.f;
629 
630  if (iH < nHitsTrackMvd) {
631  if (!listMvdHits) continue;
632  int hitIndex = t[iVec]->GetMvdHitIndex(iH);
633  CbmMvdHit* hit =
634  L1_DYNAMIC_CAST<CbmMvdHit*>(listMvdHits->At(hitIndex));
635 
636  posx = hit->GetX();
637  posy = hit->GetY();
638  posz = hit->GetZ();
639  ista = hit->GetStationNr();
640  } else {
641  if (!listStsHits) continue;
642  int hitIndex = t[iVec]->GetHitIndex(iH - nHitsTrackMvd);
643  CbmStsHit* hit =
644  L1_DYNAMIC_CAST<CbmStsHit*>(listStsHits->At(hitIndex));
645 
646  posx = hit->GetX();
647  posy = hit->GetY();
648  posz = hit->GetZ();
649  ista = CbmStsSetup::Instance()->GetStationNumber(hit->GetAddress())
650  + NMvdStations; //hit->GetStationNr()-1+NMvdStations;
651  }
652 
653  sta[ista].fieldSlice.GetFieldValue(posx, posy, fB_temp);
654  fB[iH + 1].x[iVec] = fB_temp.x[iVec];
655  fB[iH + 1].y[iVec] = fB_temp.y[iVec];
656  fB[iH + 1].z[iVec] = fB_temp.z[iVec];
657  zField[iH + 1][iVec] = posz;
658  }
659  }
660 
661  fB[0] = CbmL1::Instance()->algo->GetVtxFieldValue();
662  zField[0] = 0;
663  fld.Set(fB[2], zField[2], fB[1], zField[1], fB[0], zField[0]);
664  field.push_back(fld);
665  }
666 }
667 
668 void CbmL1PFFitter::CalculateFieldRegionAtLastPoint(
669  vector<CbmStsTrack>& Tracks,
670  vector<L1FieldRegion>& field) {
671  field.reserve(Tracks.size());
672 
673  L1FieldRegion fld _fvecalignment;
674 
675  FairRootManager* fManger = FairRootManager::Instance();
676  TClonesArray* listStsHits = (TClonesArray*) fManger->GetObject("StsHit");
677  TClonesArray* listMvdHits = 0;
678  int NMvdStations = CbmL1::Instance()->algo->NMvdStations;
679  if (NMvdStations > 0.)
680  listMvdHits = (TClonesArray*) fManger->GetObject("MvdHit");
681 
682  int nTracks_SIMD = fvecLen;
683  L1TrackPar T; // fitting parametr coresponding to current track
684 
685  CbmStsTrack* t[fvecLen];
686 
687  int ista;
688  L1Station* sta = CbmL1::Instance()->algo->vStations;
689  L1FieldValue fB[3], fB_temp _fvecalignment;
690  fvec zField[3];
691 
692  unsigned short N_vTracks = Tracks.size();
693 
694  for (unsigned short itrack = 0; itrack < N_vTracks; itrack += fvecLen) {
695  if (N_vTracks - itrack < static_cast<unsigned short>(fvecLen))
696  nTracks_SIMD = N_vTracks - itrack;
697 
698  for (int i = 0; i < nTracks_SIMD; i++)
699  t[i] = &Tracks[itrack + i]; // current track
700 
701  for (int iVec = 0; iVec < nTracks_SIMD; iVec++) {
702  int nHitsTrackMvd = t[iVec]->GetNofMvdHits();
703  int nHits = t[iVec]->GetNofHits();
704  for (int iH = 0; iH < 3; iH++) {
705  float posx = 0.f, posy = 0.f, posz = 0.f;
706 
707  int hitNumber = nHits - iH - 1;
708  if (hitNumber < nHitsTrackMvd) {
709  if (!listMvdHits) continue;
710  int hitIndex = t[iVec]->GetMvdHitIndex(hitNumber);
711  CbmMvdHit* hit =
712  L1_DYNAMIC_CAST<CbmMvdHit*>(listMvdHits->At(hitIndex));
713 
714  posx = hit->GetX();
715  posy = hit->GetY();
716  posz = hit->GetZ();
717  ista = hit->GetStationNr();
718  } else {
719  if (!listStsHits) continue;
720  int hitIndex = t[iVec]->GetHitIndex(hitNumber - nHitsTrackMvd);
721  CbmStsHit* hit =
722  L1_DYNAMIC_CAST<CbmStsHit*>(listStsHits->At(hitIndex));
723 
724  posx = hit->GetX();
725  posy = hit->GetY();
726  posz = hit->GetZ();
727  ista = CbmStsSetup::Instance()->GetStationNumber(hit->GetAddress())
728  + NMvdStations; //hit->GetStationNr()-1+NMvdStations;
729  }
730 
731  sta[ista].fieldSlice.GetFieldValue(posx, posy, fB_temp);
732 
733  fB[iH].x[iVec] = fB_temp.x[iVec];
734  fB[iH].y[iVec] = fB_temp.y[iVec];
735  fB[iH].z[iVec] = fB_temp.z[iVec];
736  zField[iH][iVec] = posz;
737  }
738  }
739 
740  fld.Set(fB[0], zField[0], fB[1], zField[1], fB[2], zField[2]);
741  field.push_back(fld);
742  }
743 }
CbmHit::GetZ
Double_t GetZ() const
Definition: CbmHit.h:70
CbmL1PFFitter::~CbmL1PFFitter
~CbmL1PFFitter()
Definition: CbmL1PFFitter.cxx:53
L1TrackPar::C10
fvec C10
Definition: L1TrackPar.h:9
L1AddTargetMaterial
void L1AddTargetMaterial(L1TrackPar &T, fvec qp0, fvec w=1, fvec mass2=0.10565f *0.10565f)
Definition: L1AddMaterial.h:771
L1TrackPar::qp
fvec qp
Definition: L1TrackPar.h:9
CbmL1PFFitter::CbmL1PFFitter
CbmL1PFFitter()
Definition: CbmL1PFFitter.cxx:51
f
float f
Definition: L1/vectors/P4_F32vec4.h:24
L1Algo.h
L1TrackPar::C20
fvec C20
Definition: L1TrackPar.h:9
_fvecalignment
#define _fvecalignment
Definition: L1/vectors/P4_F32vec4.h:254
F32vec4
Definition: L1/vectors/P4_F32vec4.h:47
CbmPixelHit::GetX
Double_t GetX() const
Definition: CbmPixelHit.h:83
L1Station.h
L1XYMeasurementInfo::C11
fvec C11
Definition: L1XYMeasurementInfo.h:11
CbmL1::algo
L1Algo * algo
Definition: CbmL1.h:119
L1XYMeasurementInfo::C10
fvec C10
Definition: L1XYMeasurementInfo.h:11
CbmPixelHit::GetY
Double_t GetY() const
Definition: CbmPixelHit.h:84
L1TrackPar::C41
fvec C41
Definition: L1TrackPar.h:10
CbmStsSetup.h
L1Station
Definition: L1Station.h:9
L1Station::fieldSlice
L1FieldSlice fieldSlice
Definition: L1Station.h:30
sqrt
friend F32vec4 sqrt(const F32vec4 &a)
Definition: L1/vectors/P4_F32vec4.h:41
L1TrackPar::C30
fvec C30
Definition: L1TrackPar.h:9
L1TrackPar::C31
fvec C31
Definition: L1TrackPar.h:9
CbmL1PFFitter::Fit
void Fit(std::vector< CbmStsTrack > &Tracks, std::vector< int > &pidHypo)
Definition: CbmL1PFFitter.cxx:81
L1Algo::GetVtxFieldValue
const L1FieldValue & GetVtxFieldValue() const
Definition: L1Algo.h:453
fvec
F32vec4 fvec
Definition: L1/vectors/P4_F32vec4.h:249
L1FieldValue::z
fvec z
Definition: L1Field.h:17
L1Extrapolate
void L1Extrapolate(L1TrackPar &T, fvec z_out, fvec qp0, const L1FieldRegion &F, fvec *w=0)
Definition: L1Extrapolation.h:314
CbmKFVertex::GetRefX
Double_t & GetRefX()
Definition: CbmKFVertex.h:23
i
int i
Definition: L1/vectors/P4_F32vec4.h:25
CbmMvdHit::GetStationNr
virtual Int_t GetStationNr() const
Definition: CbmMvdHit.h:61
L1TrackPar::C42
fvec C42
Definition: L1TrackPar.h:10
CbmStsSetup::Instance
static CbmStsSetup * Instance()
Definition: CbmStsSetup.cxx:293
L1Station::backInfo
L1UMeasurementInfo backInfo
Definition: L1Station.h:31
L1TrackPar::z
fvec z
Definition: L1TrackPar.h:9
CbmL1.h
CbmTrack::SetParamLast
void SetParamLast(const FairTrackParam *par)
Definition: CbmTrack.h:76
CbmTrack::SetNDF
void SetNDF(Int_t ndf)
Definition: CbmTrack.h:71
L1UMeasurementInfo::sin_phi
fvec sin_phi
Definition: L1UMeasurementInfo.h:12
L1TrackPar::ty
fvec ty
Definition: L1TrackPar.h:9
CbmL1PFFitter::CalculateFieldRegion
void CalculateFieldRegion(std::vector< CbmStsTrack > &Tracks, std::vector< L1FieldRegion > &Field)
Definition: CbmL1PFFitter.cxx:593
L1TrackPar::C32
fvec C32
Definition: L1TrackPar.h:9
L1AddMaterial
void L1AddMaterial(L1TrackPar &T, fvec radThick, fvec qp0, fvec w=1, fvec mass2=0.10565f *0.10565f)
Definition: L1AddMaterial.h:559
CbmStsTrack::GetNofHits
virtual Int_t GetNofHits() const
Definition: CbmStsTrack.h:76
CbmMvdHit
Definition: CbmMvdHit.h:29
L1TrackPar::y
fvec y
Definition: L1TrackPar.h:9
CbmL1PFFitter::CalculateFieldRegionAtLastPoint
void CalculateFieldRegionAtLastPoint(std::vector< CbmStsTrack > &Tracks, std::vector< L1FieldRegion > &field)
Definition: CbmL1PFFitter.cxx:668
CbmStsTrack::GetNofMvdHits
Int_t GetNofMvdHits() const
Definition: CbmStsTrack.h:84
CbmStsHit
data class for a reconstructed 3-d hit in the STS
Definition: CbmStsHit.h:31
L1FieldValue::y
fvec y
Definition: L1Field.h:17
L1Algo::NStations
int NStations
Definition: L1Algo.h:333
L1TrackPar::C44
fvec C44
Definition: L1TrackPar.h:10
CbmStsTrack.h
Data class for STS tracks.
L1FieldRegion
Definition: L1Field.h:85
fvecLen
const int fvecLen
Definition: L1/vectors/P4_F32vec4.h:251
EnergyLossCorrection
void EnergyLossCorrection(L1TrackPar &T, const fvec &mass2, const fvec &radThick, fvec &qp0, fvec direction, fvec w=1, fvec=0)
Definition: L1AddMaterial.h:235
L1Station::z
fvec z
Definition: L1Station.h:28
L1AddPipeMaterial
void L1AddPipeMaterial(L1TrackPar &T, fvec qp0, fvec w=1, fvec mass2=0.10565f *0.10565f)
Definition: L1AddMaterial.h:737
L1TrackPar::chi2
fvec chi2
Definition: L1TrackPar.h:10
CbmKFVertex::GetRefZ
Double_t & GetRefZ()
Definition: CbmKFVertex.h:25
CbmL1::Instance
static CbmL1 * Instance()
Definition: CbmL1.h:129
CbmL1PFFitter::GetChiToVertex
void GetChiToVertex(std::vector< CbmStsTrack > &Tracks, std::vector< L1FieldRegion > &field, std::vector< float > &chiToVtx, CbmKFVertex &primVtx, float chiPrim=-1)
Definition: CbmL1PFFitter.cxx:403
L1Filter
void L1Filter(L1TrackPar &T, L1UMeasurementInfo &info, fvec u, fvec w=1.)
Definition: L1Filtration.h:80
L1TrackPar::NDF
fvec NDF
Definition: L1TrackPar.h:10
d
double d
Definition: P4_F64vec2.h:24
CbmStsSetup::GetStationNumber
Int_t GetStationNumber(Int_t address)
Definition: CbmStsSetup.cxx:187
L1TrackPar::tx
fvec tx
Definition: L1TrackPar.h:9
CbmTrack::GetHitIndex
Int_t GetHitIndex(Int_t iHit) const
Definition: CbmTrack.h:54
CbmHit::GetAddress
Int_t GetAddress() const
Definition: CbmHit.h:73
L1TrackPar.h
NS_L1TrackFitter::c_light_i
const fvec c_light_i
Definition: L1TrackFitter.cxx:31
CbmL1PFFitter.h
L1TrackPar::C21
fvec C21
Definition: L1TrackPar.h:9
CbmTrack::SetChiSq
void SetChiSq(Double_t chiSq)
Definition: CbmTrack.h:70
CbmL1Track.h
L1TrackPar::x
fvec x
Definition: L1TrackPar.h:9
L1TrackPar::C00
fvec C00
Definition: L1TrackPar.h:9
CbmL1PFFitter::FilterFirst
void FilterFirst(L1TrackPar &track, fvec &x, fvec &y, L1Station &st)
Definition: CbmL1PFFitter.cxx:55
L1FieldSlice::GetFieldValue
void GetFieldValue(const fvec &x, const fvec &y, L1FieldValue &B) const
Definition: L1Field.h:41
L1Extrapolation.h
CbmTrack::GetParamFirst
const FairTrackParam * GetParamFirst() const
Definition: CbmTrack.h:61
L1TrackPar::C33
fvec C33
Definition: L1TrackPar.h:9
L1Algo::NMvdStations
int NMvdStations
Definition: L1Algo.h:334
v
__m128 v
Definition: L1/vectors/P4_F32vec4.h:1
CbmKFVertex::GetCovMatrix
Double_t * GetCovMatrix()
Definition: CbmKFVertex.h:26
L1AddMaterial.h
x
Double_t x
Definition: CbmMvdSensorDigiToHitTask.cxx:68
L1TrackPar
Definition: L1TrackPar.h:6
fabs
friend F32vec4 fabs(const F32vec4 &a)
Definition: L1/vectors/P4_F32vec4.h:60
y
Double_t y
Definition: CbmMvdSensorDigiToHitTask.cxx:68
L1TrackPar::C40
fvec C40
Definition: L1TrackPar.h:10
NS_L1TrackFitter::vINF
const fvec vINF
Definition: L1TrackFitter.cxx:34
CbmTrack::SetParamFirst
void SetParamFirst(const FairTrackParam *par)
Definition: CbmTrack.h:75
L1TrackPar::C11
fvec C11
Definition: L1TrackPar.h:9
CbmStsTrack::GetMvdHitIndex
Int_t GetMvdHitIndex(Int_t iHit) const
Definition: CbmStsTrack.h:70
CbmKFVertex::GetRefY
Double_t & GetRefY()
Definition: CbmKFVertex.h:24
CbmStsAddress.h
L1UMeasurementInfo::cos_phi
fvec cos_phi
Definition: L1UMeasurementInfo.h:12
CbmStsTrack
Definition: CbmStsTrack.h:37
L1MaterialInfo.h
PipeRadThick
const fvec PipeRadThick
Definition: L1AddMaterial.h:11
L1TrackPar::C43
fvec C43
Definition: L1TrackPar.h:10
NS_L1TrackFitter
Definition: L1TrackFitter.cxx:30
L1FieldValue
Definition: L1Field.h:11
L1XYMeasurementInfo::C00
fvec C00
Definition: L1XYMeasurementInfo.h:11
NS_L1TrackFitter::ZERO
const fvec ZERO
Definition: L1TrackFitter.cxx:32
CbmKFVertex.h
L1FieldValue::Combine
void Combine(L1FieldValue &B, fvec w)
Definition: L1Field.h:19
NS_L1TrackFitter::ONE
const fvec ONE
Definition: L1TrackFitter.cxx:33
L1Filtration.h
L1Station::frontInfo
L1UMeasurementInfo frontInfo
Definition: L1Station.h:31
CbmKFVertex
Definition: CbmKFVertex.h:6
L1FieldValue::x
fvec x
Definition: L1Field.h:15
L1TrackPar::C22
fvec C22
Definition: L1TrackPar.h:9
CbmStsTrack::GetNofStsHits
Int_t GetNofStsHits() const
Definition: CbmStsTrack.h:90
NS_L1TrackFitter::c_light
const fvec c_light
Definition: L1TrackFitter.cxx:31
L1Station::XYInfo
L1XYMeasurementInfo XYInfo
Definition: L1Station.h:32