CbmRoot
Loading...
Searching...
No Matches
CbmL1PFFitter.cxx
Go to the documentation of this file.
1/* Copyright (C) 2011-2021 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
2 SPDX-License-Identifier: GPL-3.0-only
3 Authors: Maksym Zyzak [committer] */
4
5/*
6 *=====================================================
7 *
8 * CBM Level 1 Reconstruction
9 *
10 * Authors: M.Zyzak
11 *
12 * e-mail :
13 *
14 *=====================================================
15 *
16 * SIMD Fitter
17 *
18 */
19
20#include "CbmL1PFFitter.h"
21
22#include "CbmL1.h"
23#include "CbmStsAddress.h"
24#include "CbmStsHit.h"
25#include "CbmStsSetup.h"
26#include "CbmStsTrack.h"
27#include "TClonesArray.h"
28
29//ca::Framework tools
30#include "CaFramework.h"
31#include "CaSimd.h"
32#include "CaStation.h"
33#include "CaToolsField.h"
34#include "CbmKFVertex.h"
35#include "FairRootManager.h"
36#include "KFParticleDatabase.h"
37#include "KfFieldRegion.h"
38#include "KfTrackKalmanFilter.h"
39#include "KfTrackParam.h"
40#include "TDatabasePDG.h"
41
42using namespace cbm::algo::ca;
43using cbm::algo::kf::TrackParamD;
44
45using std::vector;
46using namespace std;
47
48namespace NS_L1TrackFitter
49{
50 const fvec c_light(0.000299792458), c_light_i = fvec(1.) / c_light;
51 const fvec ZERO = fvec(0.);
52 const fvec ONE = fvec(1.);
53 const fvec vINF = fvec(0.1);
54} // namespace NS_L1TrackFitter
55using namespace NS_L1TrackFitter;
56
57
58inline void CbmL1PFFitter::PFFieldRegion::setFromL1FieldRegion(const kf::FieldRegion<fvec>& fld, int ind)
59{
60 assert(fld.GetFieldMode() == kf::EFieldMode::Intrpl);
61 const auto& coeff = fld.GetIntrplField()->GetCoefficients();
62
63 int i = 0;
64 for (int j = 0; j < 3; j++) {
65 for (int k = 0; k < 3; k++) {
66 fP[i] = coeff[j][k][ind];
67 i++;
68 }
69 }
70 fP[9] = fld.GetIntrplField()->GetZfirst()[ind];
71}
72
73inline void CbmL1PFFitter::PFFieldRegion::getL1FieldRegion(kf::FieldRegion<fvec>& fld, int ind)
74{
75 assert(fld.GetFieldMode() == kf::EFieldMode::Intrpl);
76
77 auto coeff = fld.GetIntrplField()->GetCoefficients();
78 auto z = fld.GetIntrplField()->GetZfirst();
79
80 int i = 0;
81 for (int j = 0; j < 3; j++) {
82 for (int k = 0; k < 3; k++) {
83 coeff[j][k][ind] = fP[i];
84 i++;
85 }
86 }
87 z[ind] = fP[9];
88
89 kf::detail::FieldRegionBase<fvec, kf::EFieldMode::Intrpl> fldInterpolated;
90 fldInterpolated.Set(coeff, fP[9]);
91 fld = kf::FieldRegion<fvec>(kf::EFieldType::Normal, fldInterpolated);
92}
93
94inline CbmL1PFFitter::PFFieldRegion::PFFieldRegion(const kf::FieldRegion<fvec>& fld, int i)
95{
96 setFromL1FieldRegion(fld, i);
97}
98
99
100CbmL1PFFitter::CbmL1PFFitter() {}
101
102CbmL1PFFitter::~CbmL1PFFitter() {}
103
104inline void CbmL1PFFitter::Initialize()
105{
106 if (fIsInitialised) {
107 return;
108 }
109
110 fNmvdStationsActive = 0;
111 fNstsStationsActive = 0;
112 fMvdHitArray = nullptr;
113 fStsHitArray = nullptr;
114
115 FairRootManager* manager = FairRootManager::Instance();
116
117 if (!manager) {
118 LOG(fatal) << "CbmL1PFFitter: no FairRootManager";
119 }
120
121 if (!CbmL1::Instance() || !CbmL1::Instance()->fpAlgo) {
122 LOG(fatal) << "CbmL1PFFitter: no CbmL1 task initialised ";
123 }
124
125 using cbm::algo::ca::EDetectorID;
126 // TODO: Remove CbmL1::Instance with itterations:
127 // 1) replace material with active setup
128 // 2) replace field with active setup
129 // 3) provide proper initialization of the setup
130 fNmvdStationsActive = CbmL1::Instance()->fpAlgo->GetParameters().GetNstationsActive(EDetectorID::kMvd);
131 fNstsStationsActive = CbmL1::Instance()->fpAlgo->GetParameters().GetNstationsActive(EDetectorID::kSts);
132
133 if (fNmvdStationsActive > 0) {
134 fMvdHitArray = static_cast<TClonesArray*>(manager->GetObject("MvdHit"));
135 }
136 if (fNstsStationsActive > 0) {
137 fStsHitArray = static_cast<TClonesArray*>(manager->GetObject("StsHit"));
138 }
139
140 cbm::ca::tools::SetOriginalCbmField();
141
142 fIsInitialised = true;
143}
144
145inline int CbmL1PFFitter::GetMvdStationIndex(const CbmMvdHit* hit)
146{
147 using cbm::algo::ca::EDetectorID;
148 return CbmL1::Instance()->fpAlgo->GetParameters().GetStationIndexActive(hit->GetStationNr(), EDetectorID::kMvd);
149}
150
151inline int CbmL1PFFitter::GetStsStationIndex(const CbmStsHit* hit)
152{
153 using cbm::algo::ca::EDetectorID;
154 return CbmL1::Instance()->fpAlgo->GetParameters().GetStationIndexActive(
155 CbmStsSetup::Instance()->GetStationNumber(hit->GetAddress()), EDetectorID::kSts);
156}
157
158
159void FilterFirst(kf::TrackKalmanFilter<fvec>& fit, kf::MeasurementXy<fvec>& mxy, fvec& t, fvec& dt2)
160{
161 TrackParamV& tr = fit.Tr();
162 tr.ResetErrors(mxy.Dx2(), mxy.Dy2(), 1., 1., 1., dt2, 1.e2);
163 tr.C10() = mxy.Dxy();
164 tr.X() = mxy.X();
165 tr.Y() = mxy.Y();
166 tr.Time() = t;
167 tr.Vi() = 0.;
168 tr.Ndf() = mxy.NdfX() + mxy.NdfY() - fvec(5.);
169 tr.NdfTime() = -2.;
170}
171
172void CbmL1PFFitter::Fit(std::vector<CbmStsTrack>& Tracks, const std::vector<CbmMvdHit>& vMvdHits,
173 const std::vector<CbmStsHit>& vStsHits, const std::vector<int>& pidHypo)
174{
175 // TODO: (!) replace CbmL1::Instance()->fpAlgo->GetParameters() with the cbm::ca::ParametersHandler::Instance()->Get()
176 // everywhere outside cbm::algo
177 const auto& activeTrackingSetup = CbmL1::Instance()->fpAlgo->GetParameters().GetActiveSetup();
178 Initialize();
179
180 kf::FieldValue<fvec> b0, b1, b2 _fvecalignment;
181 kf::FieldRegion<fvec> fld;
182 // fld.SetUseOriginalField();
183
184 static int nHits = CbmL1::Instance()->fpAlgo->GetParameters().GetNstationsActive();
185 int nTracks_SIMD = fvec::size();
186
187 kf::TrackKalmanFilter<fvec> fit;
188 fit.SetParticleMass(CbmL1::Instance()->fpAlgo->GetDefaultParticleMass());
189
190 TrackParamV& T = fit.Tr(); // fitting parametr coresponding to current track
191
192 CbmStsTrack* tr[fvec::size()]{nullptr};
193
194 int ista;
195 const ca::Station<fvec>* sta = CbmL1::Instance()->fpAlgo->GetParameters().GetStations().begin();
196
197 fvec x[nHits];
198 fvec y[nHits];
199 fvec z[nHits];
200 fvec t[nHits];
201 kf::MeasurementXy<fvec> mxy[nHits];
202 fvec dt2[nHits];
203 fmask w[nHits];
204
205 // fvec y_temp;
206 fvec x_first, y_first, t_first, x_last, y_last, t_last;
207 kf::MeasurementXy<fvec> mxy_first, mxy_last;
208 fvec dt2_first, dt2_last;
209
210 fvec z0, z1, z2, dz, z_start, z_end;
211 kf::FieldValue<fvec> fB[nHits];
212 kf::FieldValue<fvec> fB_temp _fvecalignment;
213
214 unsigned short N_vTracks = Tracks.size();
215
216
217 for (unsigned short itrack = 0; itrack < N_vTracks; itrack++) {
218 Tracks[itrack].SetPidHypo(pidHypo[itrack]);
219 }
220
221 fvec mass = 0.000511f;
222
223 for (unsigned short itrack = 0; itrack < N_vTracks; itrack += fvec::size()) {
224
225 T.Time() = 0.;
226 T.Vi() = 0.;
227 T.ResetErrors(1.e2, 1.e2, 1., 1., 1., 1.e6, 1.e2);
228
229 if (N_vTracks - itrack < static_cast<unsigned short>(fvec::size())) {
230 nTracks_SIMD = N_vTracks - itrack;
231 }
232 for (int iVec = 0; iVec < nTracks_SIMD; iVec++) {
233 tr[iVec] = &Tracks[itrack + iVec]; // current track
234 T.X()[iVec] = tr[iVec]->GetParamFirst()->GetX();
235 T.Y()[iVec] = tr[iVec]->GetParamFirst()->GetY();
236 T.Tx()[iVec] = tr[iVec]->GetParamFirst()->GetTx();
237 T.Ty()[iVec] = tr[iVec]->GetParamFirst()->GetTy();
238 T.Qp()[iVec] = tr[iVec]->GetParamFirst()->GetQp();
239 T.Time()[iVec] = 0.;
240 T.Vi()[iVec] = 0.;
241 T.Z()[iVec] = tr[iVec]->GetParamFirst()->GetZ();
242
243 for (int i = 0; i < 5; i++) {
244 for (int j = 0; j <= i; j++) {
245 T.C(i, j)[iVec] = tr[iVec]->GetParamFirst()->GetCovariance(i, j);
246 }
247 }
248
249 int pid = pidHypo[itrack + iVec];
250 if (pid == -1) {
251 pid = 211;
252 }
253 // mass[i] = TDatabasePDG::Instance()->GetParticle(pid)->Mass();
254 mass[iVec] = KFParticleDatabase::Instance()->GetMass(pid);
255 }
256
257 fit.SetParticleMass(mass);
258
259 // get hits of current track
260 for (int i = 0; i < nHits; i++) {
261 w[i] = fmask::Zero();
262 z[i] = sta[i].fZ;
263 x[i] = 0.;
264 y[i] = 0.;
265 t[i] = 0.;
266 mxy[i].SetX(0.);
267 mxy[i].SetY(0.);
268 mxy[i].SetDx2(1.);
269 mxy[i].SetDy2(1.);
270 mxy[i].SetDxy(0.);
271 mxy[i].SetNdfX(1.);
272 mxy[i].SetNdfY(1.);
273 dt2[i] = 1.;
274 fB[i].Set(0., 0., 0.);
275 }
276
277 for (int iVec = 0; iVec < nTracks_SIMD; iVec++) {
278 int nHitsTrackMvd = tr[iVec]->GetNofMvdHits();
279 int nHitsTrackSts = tr[iVec]->GetNofStsHits();
280 int nHitsTrack = nHitsTrackMvd + nHitsTrackSts;
281 for (int i = 0; i < nHitsTrack; i++) {
282
283 const CbmPixelHit* hit;
284 if (i < nHitsTrackMvd) {
285 int hitIndex = tr[iVec]->GetMvdHitIndex(i);
286 const CbmMvdHit* mvdHit = &(vMvdHits[hitIndex]);
287 ista = GetMvdStationIndex(mvdHit);
288 if (ista < 0) {
289 continue;
290 }
291 hit = mvdHit;
292 }
293 else {
294 int hitIndex = tr[iVec]->GetStsHitIndex(i - nHitsTrackMvd);
295 const CbmStsHit* stsHit = &(vStsHits[hitIndex]);
296 ista = GetStsStationIndex(stsHit);
297 if (ista < 0) {
298 continue;
299 }
300 hit = stsHit;
301 }
302
303 x[ista][iVec] = hit->GetX();
304 y[ista][iVec] = hit->GetY();
305 z[ista][iVec] = hit->GetZ();
306 t[ista][iVec] = hit->GetTime();
307
308 mxy[ista].X()[iVec] = hit->GetX();
309 mxy[ista].Y()[iVec] = hit->GetY();
310 mxy[ista].Dx2()[iVec] = hit->GetDx() * hit->GetDx();
311 mxy[ista].Dy2()[iVec] = hit->GetDy() * hit->GetDy();
312 mxy[ista].Dxy()[iVec] = hit->GetDxy();
313 mxy[ista].NdfX()[iVec] = 1.;
314 mxy[ista].NdfY()[iVec] = 1.;
315 dt2[ista][iVec] = hit->GetTimeError() * hit->GetTimeError();
316
317 w[ista][iVec] = true;
318
319
320 fB_temp = sta[ista].fieldSlice.GetFieldValue(x[ista], y[ista]);
321 fB[ista].SetSimdEntry(fB_temp.GetBx()[iVec], fB_temp.GetBy()[iVec], fB_temp.GetBz()[iVec], iVec);
322 if (i == 0) {
323 z_start[iVec] = z[ista][iVec];
324 x_first[iVec] = x[ista][iVec];
325 y_first[iVec] = y[ista][iVec];
326 t_first[iVec] = t[ista][iVec];
327 mxy_first.X()[iVec] = mxy[ista].X()[iVec];
328 mxy_first.Y()[iVec] = mxy[ista].Y()[iVec];
329 mxy_first.Dx2()[iVec] = mxy[ista].Dx2()[iVec];
330 mxy_first.Dy2()[iVec] = mxy[ista].Dy2()[iVec];
331 mxy_first.Dxy()[iVec] = mxy[ista].Dxy()[iVec];
332 mxy_first.NdfX()[iVec] = mxy[ista].NdfX()[iVec];
333 mxy_first.NdfY()[iVec] = mxy[ista].NdfY()[iVec];
334 dt2_first[iVec] = dt2[ista][iVec];
335 }
336 if (i == nHitsTrack - 1) {
337 z_end[iVec] = z[ista][iVec];
338 x_last[iVec] = x[ista][iVec];
339 y_last[iVec] = y[ista][iVec];
340 t_last[iVec] = t[ista][iVec];
341 mxy_last.X()[iVec] = mxy[ista].X()[iVec];
342 mxy_last.Y()[iVec] = mxy[ista].Y()[iVec];
343 mxy_last.Dx2()[iVec] = mxy[ista].Dx2()[iVec];
344 mxy_last.Dy2()[iVec] = mxy[ista].Dy2()[iVec];
345 mxy_last.Dxy()[iVec] = mxy[ista].Dxy()[iVec];
346 mxy_last.NdfX()[iVec] = mxy[ista].NdfX()[iVec];
347 mxy_last.NdfY()[iVec] = mxy[ista].NdfY()[iVec];
348 dt2_last[iVec] = dt2[ista][iVec];
349 }
350 }
351 }
352
353 // fit forward
354
355 int i = 0;
356
357 FilterFirst(fit, mxy_first, t_first, dt2_first);
358 fit.SetQp0(fit.Tr().Qp());
359
360 z1 = z[i];
361 b1 = sta[i].fieldSlice.GetFieldValue(T.X(), T.Y());
362 b1.SetSimdEntries(fB[i], w[i]);
363 z2 = z[i + 2];
364 dz = z2 - z1;
365 b2 = sta[i].fieldSlice.GetFieldValue(T.X() + T.Tx() * dz, T.Y() + T.Ty() * dz);
366 b2.SetSimdEntries(fB[i + 2], w[i + 2]);
367 fld.Set(b2, z2, b1, z1, b0, z0);
368 for (++i; i < nHits; i++) {
369 z0 = z[i];
370 dz = (z1 - z0);
371 b0 = sta[i].fieldSlice.GetFieldValue(T.X() - T.Tx() * dz, T.Y() - T.Ty() * dz);
372 b0.SetSimdEntries(fB[i], w[i]);
373 fld.Set(b0, z0, b1, z1, b2, z2);
374
375 fmask initialised = (z[i] <= z_end) & (z_start < z[i]);
376
377 fit.SetMask(initialised);
378 fit.Extrapolate(z[i], fld);
379 auto radThick = activeTrackingSetup.GetMaterial(i).GetThicknessX0(fit.Tr().X(), fit.Tr().Y());
380 fit.MultipleScattering(radThick);
381 fit.EnergyLossCorrection(radThick, kf::FitDirection::kDownstream);
382
383 fit.SetMask(initialised && w[i]);
384 fit.FilterXY(mxy[i]);
385 fit.FilterTime(t[i], dt2[i], fmask(sta[i].timeInfo));
386
387 b2 = b1;
388 z2 = z1;
389 b1 = b0;
390 z1 = z0;
391 }
392
393 TrackParamV Tout = T;
394 for (int iVec = 0; iVec < nTracks_SIMD; iVec++) {
395 FairTrackParam par;
396 par.SetX(T.X()[iVec]);
397 par.SetY(T.Y()[iVec]);
398 par.SetTx(T.Tx()[iVec]);
399 par.SetTy(T.Ty()[iVec]);
400 par.SetQp(T.Qp()[iVec]);
401 par.SetZ(T.Z()[iVec]);
402
403 for (int k = 0; k < 5; k++) {
404 for (int j = 0; j <= k; j++) {
405 par.SetCovariance(k, j, Tout.C(k, j)[iVec]);
406 }
407 }
408
409 tr[iVec]->SetParamLast(&par);
410 }
411
412 // fit backward
413
414 fit.SetQp0(T.Qp());
415
416 i = nHits - 1;
417
418 FilterFirst(fit, mxy_last, t_last, dt2_last);
419
420 z1 = z[i];
421 b1 = sta[i].fieldSlice.GetFieldValue(T.X(), T.Y());
422 b1.SetSimdEntries(fB[i], w[i]);
423
424 z2 = z[i - 2];
425 dz = z2 - z1;
426 b2 = sta[i].fieldSlice.GetFieldValue(T.X() + T.Tx() * dz, T.Y() + T.Ty() * dz);
427 b2.SetSimdEntries(fB[i - 2], w[i - 2]);
428 fld.Set(b2, z2, b1, z1, b0, z0);
429 for (--i; i >= 0; i--) {
430 z0 = z[i];
431 dz = (z1 - z0);
432 b0 = sta[i].fieldSlice.GetFieldValue(T.X() - T.Tx() * dz, T.Y() - T.Ty() * dz);
433 b0.SetSimdEntries(fB[i], w[i]);
434 fld.Set(b0, z0, b1, z1, b2, z2);
435
436 fmask initialised = (z[i] < z_end) & (z_start <= z[i]);
437
438 fit.SetMask(initialised);
439 fit.Extrapolate(z[i], fld);
440 auto radThick = activeTrackingSetup.GetMaterial(i).GetThicknessX0(fit.Tr().X(), fit.Tr().Y());
441 fit.MultipleScattering(radThick);
442 fit.EnergyLossCorrection(radThick, kf::FitDirection::kUpstream);
443
444 fit.SetMask(initialised && w[i]);
445 fit.FilterXY(mxy[i]);
446 fit.FilterTime(t[i], dt2[i], fmask(sta[i].timeInfo));
447
448 b2 = b1;
449 z2 = z1;
450 b1 = b0;
451 z1 = z0;
452 }
453
454 for (int iVec = 0; iVec < nTracks_SIMD; iVec++) {
455 FairTrackParam par;
456 par.SetX(T.X()[iVec]);
457 par.SetY(T.Y()[iVec]);
458 par.SetTx(T.Tx()[iVec]);
459 par.SetTy(T.Ty()[iVec]);
460 par.SetQp(T.Qp()[iVec]);
461 par.SetZ(T.Z()[iVec]);
462
463 for (int k = 0; k < 5; k++) {
464 for (int j = 0; j <= k; j++) {
465 par.SetCovariance(k, j, T.C(k, j)[iVec]);
466 }
467 }
468
469 tr[iVec]->SetParamFirst(&par);
470
471 tr[iVec]->SetChiSq(T.ChiSq()[iVec]);
472 tr[iVec]->SetNDF(static_cast<int>(T.Ndf()[iVec]));
473 }
474 }
475}
476
477void CbmL1PFFitter::Fit(vector<CbmStsTrack>& Tracks, const vector<int>& pidHypo)
478{
479
480 Initialize();
481
482 std::vector<CbmMvdHit> vMvdHits;
483 std::vector<CbmStsHit> vStsHits;
484
485 if (fMvdHitArray) {
486 for (int ih = 0; ih < fMvdHitArray->GetEntriesFast(); ih++) {
487 vMvdHits.push_back(*dynamic_cast<const CbmMvdHit*>(fMvdHitArray->At(ih)));
488 }
489 }
490
491 if (fStsHitArray) {
492 for (int ih = 0; ih < fStsHitArray->GetEntriesFast(); ih++) {
493 vStsHits.push_back(*dynamic_cast<const CbmStsHit*>(fStsHitArray->At(ih)));
494 }
495 }
496
497 Fit(Tracks, vMvdHits, vStsHits, pidHypo);
498}
499
500
501void CbmL1PFFitter::GetChiToVertex(vector<CbmStsTrack>& Tracks, vector<PFFieldRegion>& field, vector<float>& chiToVtx,
502 CbmKFVertex& primVtx, float chiPrim)
503{
504 // TODO: (!) replace CbmL1::Instance()->fpAlgo->GetParameters() with the cbm::ca::ParametersHandler::Instance()->Get()
505 // everywhere outside cbm::algo
506 const auto& activeTrackingSetup = CbmL1::Instance()->fpAlgo->GetParameters().GetActiveSetup();
507 Initialize();
508
509 chiToVtx.reserve(Tracks.size());
510
511 int nTracks_SIMD = fvec::size();
512
513 kf::TrackKalmanFilter<fvec> fit;
514 TrackParamV& T = fit.Tr(); // fitting parametr coresponding to current track
515
516 CbmStsTrack* tr[fvec::size()]{nullptr};
517
518 int nStations = fNmvdStationsActive + fNstsStationsActive;
519
520 const ca::Station<fvec>* sta = CbmL1::Instance()->fpAlgo->GetParameters().GetStations().begin();
521 fvec* zSta = new fvec[nStations];
522 for (int iSta = 0; iSta < nStations; iSta++) {
523 zSta[iSta] = sta[iSta].fZ;
524 }
525
526 field.reserve(Tracks.size());
527
528 kf::FieldRegion<fvec> fld _fvecalignment;
529 kf::FieldValue<fvec> fB[3], fB_temp _fvecalignment;
530 fvec zField[3];
531
532 unsigned short N_vTracks = Tracks.size();
533 int ista;
534 for (unsigned short itrack = 0; itrack < N_vTracks; itrack += fvec::size()) {
535 if (N_vTracks - itrack < static_cast<unsigned short>(fvec::size())) {
536 nTracks_SIMD = N_vTracks - itrack;
537 }
538
539 fvec mass2;
540 for (int iVec = 0; iVec < nTracks_SIMD; iVec++) {
541 tr[iVec] = &Tracks[itrack + iVec]; // current track
542 T.X()[iVec] = tr[iVec]->GetParamFirst()->GetX();
543 T.Y()[iVec] = tr[iVec]->GetParamFirst()->GetY();
544 T.Tx()[iVec] = tr[iVec]->GetParamFirst()->GetTx();
545 T.Ty()[iVec] = tr[iVec]->GetParamFirst()->GetTy();
546 T.Qp()[iVec] = tr[iVec]->GetParamFirst()->GetQp();
547 T.Z()[iVec] = tr[iVec]->GetParamFirst()->GetZ();
548
549 for (int i = 0; i < 5; i++) {
550 for (int j = 0; j <= i; j++) {
551 T.C(i, j)[iVec] = tr[iVec]->GetParamFirst()->GetCovariance(i, j);
552 }
553 }
554
555 // float mass = TDatabasePDG::Instance()->GetParticle(tr[iVec]->GetPidHypo())->Mass();
556 const float mass = KFParticleDatabase::Instance()->GetMass(tr[iVec]->GetPidHypo());
557 mass2[iVec] = mass * mass;
558
559 int nHitsTrackMvd = tr[iVec]->GetNofMvdHits();
560 for (int iH = 0; iH < 2; iH++) {
561 float posx = 0.f, posy = 0.f, posz = 0.f;
562
563 if (iH < nHitsTrackMvd) {
564 if (!fMvdHitArray) {
565 continue;
566 }
567 int hitIndex = tr[iVec]->GetMvdHitIndex(iH);
568 const CbmMvdHit* hit = dynamic_cast<const CbmMvdHit*>(fMvdHitArray->At(hitIndex));
569
570 posx = hit->GetX();
571 posy = hit->GetY();
572 posz = hit->GetZ();
573 ista = GetMvdStationIndex(hit);
574 if (ista < 0) {
575 continue;
576 }
577 }
578 else {
579 if (!fStsHitArray) {
580 continue;
581 }
582 int hitIndex = tr[iVec]->GetStsHitIndex(iH - nHitsTrackMvd);
583 const CbmStsHit* hit = dynamic_cast<const CbmStsHit*>(fStsHitArray->At(hitIndex));
584
585 posx = hit->GetX();
586 posy = hit->GetY();
587 posz = hit->GetZ();
588 ista = GetStsStationIndex(hit);
589 if (ista < 0) {
590 continue;
591 }
592 }
593
594 fB_temp = sta[ista].fieldSlice.GetFieldValue(posx, posy);
595 fB[iH + 1].SetSimdEntry(fB_temp.GetBx()[iVec], fB_temp.GetBy()[iVec], fB_temp.GetBz()[iVec], iVec);
596 zField[iH + 1][iVec] = posz;
597 }
598 }
599
600 fB[0] = CbmL1::Instance()->fpAlgo->GetParameters().GetVertexFieldValue();
601 zField[0] = CbmL1::Instance()->fpAlgo->GetParameters().GetTargetPositionZ();
602 fld.Set(fB[2], zField[2], fB[1], zField[1], fB[0], zField[0]);
603 for (int i = 0; i < nTracks_SIMD; i++) {
604 field.emplace_back(fld, i);
605 }
606
607 fit.SetQp0(fit.Tr().Qp());
608
609 for (int iSt = nStations - 4; iSt >= 0; iSt--) {
610 fit.SetMask(T.Z() > zSta[iSt] + fvec(2.5));
611 fit.Extrapolate(zSta[iSt], fld);
612 auto radThick = activeTrackingSetup.GetMaterial(iSt).GetThicknessX0(fit.Tr().X(), fit.Tr().Y());
613 fit.MultipleScattering(radThick);
614 fit.EnergyLossCorrection(radThick, kf::FitDirection::kUpstream);
615 }
616 fit.SetMask(fmask::One());
617 fit.Extrapolate(primVtx.GetRefZ(), fld);
618
619 // TODO: get it from parameters
620 constexpr float targetRadThick = 3.73e-2f * 2; // 250 mum Gold
621
622 fit.MultipleScattering(targetRadThick);
623 fit.EnergyLossCorrection(targetRadThick, kf::FitDirection::kUpstream);
624
625 Double_t Cv[3] = {primVtx.GetCovMatrix()[0], primVtx.GetCovMatrix()[1], primVtx.GetCovMatrix()[2]};
626
627 fvec dx = T.X() - fvec(primVtx.GetRefX());
628 fvec dy = T.Y() - fvec(primVtx.GetRefY());
629 fvec c[3] = {T.C00(), T.C10(), T.C11()};
630 c[0] += fvec(Cv[0]);
631 c[1] += fvec(Cv[1]);
632 c[2] += fvec(Cv[2]);
633 fvec d = c[0] * c[2] - c[1] * c[1];
634 fvec chi = sqrt(kf::utils::fabs(fvec(0.5) * (dx * dx * c[0] - fvec(2.) * dx * dy * c[1] + dy * dy * c[2]) / d));
635 chi.setZero(kf::utils::fabs(d) < fvec(1.e-20));
636
637 for (int iVec = 0; iVec < nTracks_SIMD; iVec++) {
638 chiToVtx.push_back(chi[iVec]);
639 }
640
641 if (chiPrim > 0) {
642 for (int iVec = 0; iVec < nTracks_SIMD; iVec++) {
643 if (chi[iVec] < chiPrim) {
644 FairTrackParam par;
645 par.SetX(T.X()[iVec]);
646 par.SetY(T.Y()[iVec]);
647 par.SetTx(T.Tx()[iVec]);
648 par.SetTy(T.Ty()[iVec]);
649 par.SetQp(T.Qp()[iVec]);
650 par.SetZ(T.Z()[iVec]);
651
652 for (int i = 0; i < 5; i++) {
653 for (int j = 0; j <= i; j++) {
654 par.SetCovariance(i, j, T.C(i, j)[iVec]);
655 }
656 }
657
658 tr[iVec]->SetParamFirst(&par);
659 }
660 }
661 }
662 }
663 delete[] zSta;
664}
665
666void CbmL1PFFitter::CalculateFieldRegion(vector<CbmStsTrack>& Tracks, vector<PFFieldRegion>& field)
667{
668 Initialize();
669
670 field.reserve(Tracks.size());
671
672 kf::FieldRegion<fvec> fld _fvecalignment;
673
674 int nTracks_SIMD = fvec::size();
675 TrackParamV T; // fitting parametr coresponding to current track
676
677 CbmStsTrack* tr[fvec::size()];
678
679 int ista;
680 const ca::Station<fvec>* sta = CbmL1::Instance()->fpAlgo->GetParameters().GetStations().begin();
681 kf::FieldValue<fvec> fB[3], fB_temp _fvecalignment;
682 fvec zField[3];
683
684 unsigned short N_vTracks = Tracks.size();
685
686 for (unsigned short itrack = 0; itrack < N_vTracks; itrack += fvec::size()) {
687 if (N_vTracks - itrack < static_cast<unsigned short>(fvec::size())) {
688 nTracks_SIMD = N_vTracks - itrack;
689 }
690
691 for (int i = 0; i < nTracks_SIMD; i++) {
692 tr[i] = &Tracks[itrack + i]; // current track
693 }
694
695 for (int iVec = 0; iVec < nTracks_SIMD; iVec++) {
696 int nHitsTrackMvd = tr[iVec]->GetNofMvdHits();
697 for (int iH = 0; iH < 2; iH++) {
698 float posx = 0.f, posy = 0.f, posz = 0.f;
699
700 if (iH < nHitsTrackMvd) {
701 if (!fMvdHitArray) {
702 continue;
703 }
704 int hitIndex = tr[iVec]->GetMvdHitIndex(iH);
705 const CbmMvdHit* hit = dynamic_cast<const CbmMvdHit*>(fMvdHitArray->At(hitIndex));
706
707 posx = hit->GetX();
708 posy = hit->GetY();
709 posz = hit->GetZ();
710 ista = GetMvdStationIndex(hit);
711 if (ista < 0) {
712 continue;
713 }
714 }
715 else {
716 if (!fStsHitArray) {
717 continue;
718 }
719 int hitIndex = tr[iVec]->GetStsHitIndex(iH - nHitsTrackMvd);
720 const CbmStsHit* hit = dynamic_cast<const CbmStsHit*>(fStsHitArray->At(hitIndex));
721
722 posx = hit->GetX();
723 posy = hit->GetY();
724 posz = hit->GetZ();
725 ista = GetStsStationIndex(hit);
726 if (ista < 0) {
727 continue;
728 }
729 }
730
731 fB_temp = sta[ista].fieldSlice.GetFieldValue(posx, posy);
732 fB[iH + 1].SetSimdEntry(fB_temp.GetBx()[iVec], fB_temp.GetBy()[iVec], fB_temp.GetBz()[iVec], iVec);
733 zField[iH + 1][iVec] = posz;
734 }
735 }
736
737 fB[0] = CbmL1::Instance()->fpAlgo->GetParameters().GetVertexFieldValue();
738 zField[0] = CbmL1::Instance()->fpAlgo->GetParameters().GetTargetPositionZ();
739 fld.Set(fB[2], zField[2], fB[1], zField[1], fB[0], zField[0]);
740 for (int i = 0; i < nTracks_SIMD; i++) {
741 field.emplace_back(fld, i);
742 }
743 }
744}
745
746void CbmL1PFFitter::CalculateFieldRegionAtLastPoint(vector<CbmStsTrack>& Tracks, vector<PFFieldRegion>& field)
747{
748 Initialize();
749
750 field.reserve(Tracks.size());
751
752 kf::FieldRegion<fvec> fld _fvecalignment;
753
754 int nTracks_SIMD = fvec::size();
755 TrackParamV T; // fitting parametr coresponding to current track
756
757 CbmStsTrack* tr[fvec::size()];
758
759 int ista;
760 const ca::Station<fvec>* sta = CbmL1::Instance()->fpAlgo->GetParameters().GetStations().begin();
761 kf::FieldValue<fvec> fB[3], fB_temp _fvecalignment;
762 fvec zField[3];
763
764 unsigned short N_vTracks = Tracks.size();
765
766 for (unsigned short itrack = 0; itrack < N_vTracks; itrack += fvec::size()) {
767 if (N_vTracks - itrack < static_cast<unsigned short>(fvec::size())) {
768 nTracks_SIMD = N_vTracks - itrack;
769 }
770
771 for (int i = 0; i < nTracks_SIMD; i++) {
772 tr[i] = &Tracks[itrack + i]; // current track
773 }
774
775 for (int iVec = 0; iVec < nTracks_SIMD; iVec++) {
776 int nHitsTrackMvd = tr[iVec]->GetNofMvdHits();
777 int nHits = tr[iVec]->GetTotalNofHits();
778 for (int iH = 0; iH < 3; iH++) {
779 float posx = 0.f, posy = 0.f, posz = 0.f;
780
781 int hitNumber = nHits - iH - 1;
782 if (hitNumber < nHitsTrackMvd) {
783 if (!fMvdHitArray) {
784 continue;
785 }
786 int hitIndex = tr[iVec]->GetMvdHitIndex(hitNumber);
787 const CbmMvdHit* hit = dynamic_cast<const CbmMvdHit*>(fMvdHitArray->At(hitIndex));
788
789 posx = hit->GetX();
790 posy = hit->GetY();
791 posz = hit->GetZ();
792 ista = GetMvdStationIndex(hit);
793 if (ista < 0) {
794 continue;
795 }
796 }
797 else {
798 if (!fStsHitArray) {
799 continue;
800 }
801 int hitIndex = tr[iVec]->GetStsHitIndex(hitNumber - nHitsTrackMvd);
802 const CbmStsHit* hit = dynamic_cast<const CbmStsHit*>(fStsHitArray->At(hitIndex));
803
804 posx = hit->GetX();
805 posy = hit->GetY();
806 posz = hit->GetZ();
807 ista = GetStsStationIndex(hit);
808 if (ista < 0) {
809 continue;
810 }
811 }
812
813 fB_temp = sta[ista].fieldSlice.GetFieldValue(posx, posy);
814
815 fB[iH].SetSimdEntry(fB_temp.GetBx()[iVec], fB_temp.GetBy()[iVec], fB_temp.GetBz()[iVec], iVec);
816 zField[iH][iVec] = posz;
817 }
818 }
819
820 fld.Set(fB[0], zField[0], fB[1], zField[1], fB[2], zField[2]);
821 for (int i = 0; i < nTracks_SIMD; i++) {
822 field.emplace_back(fld, i);
823 }
824 }
825}
CaToolsField.h
Tracking Field class (header)
FilterFirst
void FilterFirst(kf::TrackKalmanFilter< fvec > &fit, kf::MeasurementXy< fvec > &mxy, fvec &t, fvec &dt2)
Definition CbmL1PFFitter.cxx:159
CbmStsHit.h
Data class for a reconstructed hit in the STS.
CbmStsTrack.h
Data class for STS tracks.
KfFieldRegion.h
Magnetic flux density interpolation along the track vs. z-coordinate (header)
_fvecalignment
#define _fvecalignment
Definition KfSimdPseudo.h:281
sqrt
friend fvec sqrt(const fvec &a)
Definition KfSimdPseudo.h:118
KfTrackKalmanFilter.h
Track fit utilities for the CA tracking based on the Kalman filter.
CbmHit::GetTimeError
double GetTimeError() const
Definition CbmHit.h:77
CbmHit::GetTime
double GetTime() const
Definition CbmHit.h:76
CbmHit::GetAddress
int32_t GetAddress() const
Definition CbmHit.h:74
CbmHit::GetZ
double GetZ() const
Definition CbmHit.h:71
CbmKFVertex::GetRefZ
Double_t & GetRefZ()
Definition CbmKFVertex.h:27
CbmKFVertex::GetRefX
Double_t & GetRefX()
Definition CbmKFVertex.h:25
CbmKFVertex::GetCovMatrix
Double_t * GetCovMatrix()
Definition CbmKFVertex.h:28
CbmKFVertex::GetRefY
Double_t & GetRefY()
Definition CbmKFVertex.h:26
CbmL1PFFitter::GetMvdStationIndex
int GetMvdStationIndex(const CbmMvdHit *h)
Definition CbmL1PFFitter.cxx:145
CbmL1PFFitter::fStsHitArray
TClonesArray * fStsHitArray
Definition CbmL1PFFitter.h:68
CbmL1PFFitter::fMvdHitArray
TClonesArray * fMvdHitArray
Definition CbmL1PFFitter.h:67
CbmL1PFFitter::GetChiToVertex
void GetChiToVertex(std::vector< CbmStsTrack > &Tracks, std::vector< PFFieldRegion > &field, std::vector< float > &chiToVtx, CbmKFVertex &primVtx, float chiPrim=-1)
Definition CbmL1PFFitter.cxx:501
CbmL1PFFitter::Fit
void Fit(std::vector< CbmStsTrack > &Tracks, const std::vector< CbmMvdHit > &vMvdHits, const std::vector< CbmStsHit > &vStsHits, const std::vector< int > &pidHypo)
Definition CbmL1PFFitter.cxx:172
CbmL1PFFitter::GetStsStationIndex
int GetStsStationIndex(const CbmStsHit *h)
Definition CbmL1PFFitter.cxx:151
CbmL1PFFitter::CalculateFieldRegion
void CalculateFieldRegion(std::vector< CbmStsTrack > &Tracks, std::vector< PFFieldRegion > &Field)
Definition CbmL1PFFitter.cxx:666
CbmL1PFFitter::CalculateFieldRegionAtLastPoint
void CalculateFieldRegionAtLastPoint(std::vector< CbmStsTrack > &Tracks, std::vector< PFFieldRegion > &field)
Definition CbmL1PFFitter.cxx:746
CbmL1::Instance
static CbmL1 * Instance()
Pointer to CbmL1 instance.
Definition CbmL1.h:171
CbmL1::fpAlgo
ca::Framework * fpAlgo
Pointer to the L1 track finder algorithm.
Definition CbmL1.h:425
CbmMvdHit::GetStationNr
virtual int32_t GetStationNr() const
Definition CbmMvdHit.h:61
CbmPixelHit::GetDy
double GetDy() const
Definition CbmPixelHit.h:76
CbmPixelHit::GetDx
double GetDx() const
Definition CbmPixelHit.h:75
CbmPixelHit::GetY
double GetY() const
Definition CbmPixelHit.h:74
CbmPixelHit::GetX
double GetX() const
Definition CbmPixelHit.h:73
CbmPixelHit::GetDxy
double GetDxy() const
Definition CbmPixelHit.h:77
CbmStsHit
data class for a reconstructed 3-d hit in the STS
Definition CbmStsHit.h:35
CbmStsSetup::Instance
static CbmStsSetup * Instance()
Definition CbmStsSetup.cxx:303
CbmStsTrack::GetNofMvdHits
int32_t GetNofMvdHits() const
Definition CbmStsTrack.h:87
CbmStsTrack::GetMvdHitIndex
int32_t GetMvdHitIndex(int32_t iHit) const
Definition CbmStsTrack.h:75
CbmStsTrack::GetStsHitIndex
int32_t GetStsHitIndex(int32_t iHit) const
Definition CbmStsTrack.h:101
CbmStsTrack::GetTotalNofHits
virtual int32_t GetTotalNofHits() const
Definition CbmStsTrack.h:81
cbm::algo::ca::Framework::GetParameters
const Parameters< fvec > & GetParameters() const
Gets a pointer to the Framework parameters object.
Definition CaFramework.h:87
cbm::algo::ca::Station::fZ
DataT fZ
z position of station [cm]
Definition CaStation.h:29
cbm::algo::ca::Station::fieldSlice
kf::FieldSlice< DataT > fieldSlice
Magnetic field near the station.
Definition CaStation.h:33
cbm::algo::kf::FieldRegion
Magnetic field region, corresponding to a hit triplet.
Definition KfFieldRegion.h:257
cbm::algo::kf::FieldRegion::GetFieldMode
EFieldMode GetFieldMode() const
Gets field mode.
Definition KfFieldRegion.h:345
cbm::algo::kf::FieldRegion::GetIntrplField
const auto & GetIntrplField() const
Definition KfFieldRegion.h:385
cbm::algo::kf::FieldRegion::Set
void Set(const FieldValue< T > &b0, const T &z0, const FieldValue< T > &b1, const T &z1, const FieldValue< T > &b2, const T &z2)
Interpolates the field by three nodal points with the second order polynomial.
Definition KfFieldRegion.h:358
cbm::algo::kf::FieldSlice::GetFieldValue
constexpr FieldValue< T > GetFieldValue(const T &x, const T &y) const
Gets field value at a point on the transverse plane.
Definition KfFieldSlice.h:73
cbm::algo::kf::FieldValue
Magnetic flux density vector.
Definition KfFieldValue.h:31
cbm::algo::kf::FieldValue::SetSimdEntries
void SetSimdEntries(const FieldValue &other, const kf::utils::masktype< T > &mask)
Combines the current magnetic field value with another one using a mask.
Definition KfFieldValue.h:80
cbm::algo::kf::FieldValue::Set
void Set(const I &bx, const I &by, const I &bz)
Constructor from components.
Definition KfFieldValue.h:70
cbm::algo::kf::FieldValue::SetSimdEntry
void SetSimdEntry(double bx, double by, double bz, size_t i)
Sets magnetic flux density components to the field function.
Definition KfFieldValue.h:132
cbm::algo::kf::MeasurementXy
The class describes a 2D - measurement (x, y) in XY coordinate system.
Definition KfMeasurementXy.h:28
cbm::algo::kf::TrackKalmanFilter::Extrapolate
void Extrapolate(DataT z, const kf::FieldRegion< DataT > &F)
Definition KfTrackKalmanFilter.cxx:601
cbm::algo::kf::TrackKalmanFilter::FilterXY
void FilterXY(const kf::MeasurementXy< DataT > &m, bool skipUnmeasuredCoordinates=false)
filter the track with the XY measurement
Definition KfTrackKalmanFilter.cxx:189
cbm::algo::kf::TrackKalmanFilter::Tr
kf::TrackParam< DataT > & Tr()
Definition KfTrackKalmanFilter.h:70
cbm::algo::kf::TrackKalmanFilter::MultipleScattering
void MultipleScattering(DataT radThick, DataT tx0, DataT ty0, DataT qp0)
apply multiple scattering correction to the track with the given Qp0
Definition KfTrackKalmanFilter.cxx:980
cbm::algo::kf::TrackKalmanFilter::EnergyLossCorrection
void EnergyLossCorrection(DataT radThick, FitDirection direction)
Definition KfTrackKalmanFilter.cxx:1057
cbm::algo::kf::TrackKalmanFilter::FilterTime
void FilterTime(DataT t, DataT dt2, const DataTmask &m)
filter the track with the time measurement
Definition KfTrackKalmanFilter.cxx:105
cbm::algo::kf::TrackKalmanFilter::SetParticleMass
void SetParticleMass(DataT mass)
set particle mass for the fit
Definition KfTrackKalmanFilter.h:84
cbm::algo::kf::TrackParamBase::Ndf
T Ndf() const
Gets NDF of track fit model.
Definition KfTrackParam.h:111
cbm::algo::kf::TrackParamBase::Vi
T Vi() const
Gets inverse velocity [ns/cm].
Definition KfTrackParam.h:105
cbm::algo::kf::TrackParamBase::NdfTime
T NdfTime() const
Gets NDF of time measurements.
Definition KfTrackParam.h:117
cbm::algo::kf::TrackParamBase::ResetErrors
void ResetErrors(T c00, T c11, T c22, T c33, T c44, T c55, T c66)
Resets variances of track parameters and chi2, ndf values.
Definition KfTrackParam.h:744
cbm::algo::kf::TrackParamBase::X
T X() const
Gets x position [cm].
Definition KfTrackParam.h:87
cbm::algo::kf::TrackParamBase::Y
T Y() const
Gets y position [cm].
Definition KfTrackParam.h:90
cbm::algo::kf::TrackParamBase::Qp
T Qp() const
Gets charge over momentum [ec/GeV].
Definition KfTrackParam.h:99
cbm::algo::kf::TrackParamBase::Time
T Time() const
Gets time [ns].
Definition KfTrackParam.h:102
cbm::algo::kf::TrackParamBase::C
T C(int i, int j) const
Get covariance matrix element.
Definition KfTrackParam.h:123
cbm::algo::kf::detail::FieldRegionBase
Base class of the FieldRegion class (primary template)
Definition KfFieldRegion.h:81
cbm::algo::kf::fmask::Zero
static fmask Zero()
Definition KfSimdPseudo.h:62
cbm::algo::kf::fmask::One
static fmask One()
Definition KfSimdPseudo.h:60
cbm::algo::kf::fvec::size
static constexpr size_t size()
Definition KfSimdPseudo.h:123
cbm::algo::kf::fvec::setZero
void setZero(fmask m)
Definition KfSimdPseudo.h:161
NS_L1TrackFitter
Definition CbmL1PFFitter.cxx:49
NS_L1TrackFitter::c_light
const fvec c_light(0.000299792458)
NS_L1TrackFitter::ZERO
const fvec ZERO
Definition CbmL1PFFitter.cxx:51
NS_L1TrackFitter::vINF
const fvec vINF
Definition CbmL1PFFitter.cxx:53
NS_L1TrackFitter::ONE
const fvec ONE
Definition CbmL1PFFitter.cxx:52
NS_L1TrackFitter::c_light_i
const fvec c_light_i
Definition CbmL1PFFitter.cxx:50
cbm::algo::ca
TODO: SZh 8.11.2022: add selection of parameterisation.
Definition CaBranch.h:14
cbm::algo::ca::fmask
kf::fmask fmask
Definition CaSimd.h:15
cbm::algo::ca::EDetectorID
EDetectorID
Enumeration for the tracking detector subsystems in CBM-CA.
Definition CbmDefs.h:176
cbm::algo::ca::fvec
kf::fvec fvec
Definition CaSimd.h:13
cbm::ca::tools::SetOriginalCbmField
void SetOriginalCbmField()
pass the original magnetic field to L1Algo
Definition CaToolsField.h:21
std
Hash for CbmL1LinkKey.
Definition algo/base/Options.cxx:21
CbmL1PFFitter::PFFieldRegion::getL1FieldRegion
void getL1FieldRegion(cbm::algo::kf::FieldRegion< cbm::algo::kf::fvec > &, int i)
Definition CbmL1PFFitter.cxx:73
CbmL1PFFitter::PFFieldRegion::setFromL1FieldRegion
void setFromL1FieldRegion(const cbm::algo::kf::FieldRegion< cbm::algo::kf::fvec > &, int i)
Definition CbmL1PFFitter.cxx:58