From 3b1c1ab3e17c7b95b68883eb86134f14c51fc8be Mon Sep 17 00:00:00 2001
From: "se.gorbunov" <se.gorbunov@gsi.de>
Date: Sat, 11 Feb 2023 00:00:51 +0000
Subject: [PATCH] L1 fit: make L1Fit members private
---
reco/L1/CbmL1Performance.cxx | 33 ++--
reco/L1/L1Algo/L1BranchExtender.cxx | 31 ++--
reco/L1/L1Algo/L1CATrackFinder.cxx | 54 +++---
reco/L1/L1Algo/L1CloneMerger.cxx | 12 +-
reco/L1/L1Algo/L1Fit.h | 24 ++-
reco/L1/L1Algo/L1TrackFitter.cxx | 215 ++++++++++++-----------
reco/L1/ParticleFinder/CbmL1PFFitter.cxx | 49 +++---
7 files changed, 212 insertions(+), 206 deletions(-)
diff --git a/reco/L1/CbmL1Performance.cxx b/reco/L1/CbmL1Performance.cxx
index f18fb3dc8a..0c06f06425 100644
--- a/reco/L1/CbmL1Performance.cxx
+++ b/reco/L1/CbmL1Performance.cxx
@@ -1228,12 +1228,8 @@ void CbmL1::TrackFitPerformance()
CbmL1MCTrack mc = *(it->GetMCTracks()[0]);
- {
- L1TrackPar trPar(it->T, it->C);
- fit.fTr = trPar;
- }
- fit.fQp0 = fit.fTr.qp;
- const L1TrackPar& tr = fit.fTr;
+ fit.SetTrack(it->T, it->C);
+ const L1TrackPar& tr = fit.Tr();
L1FieldRegion fld _fvecalignment;
fld.SetUseOriginalField();
@@ -1338,12 +1334,8 @@ void CbmL1::TrackFitPerformance()
L1FieldRegion fld _fvecalignment;
fld.SetUseOriginalField();
- {
- L1TrackPar trPar(it->TLast, it->CLast);
- fit.fTr = trPar;
- }
- fit.fQp0 = fit.fTr.qp;
- const L1TrackPar& tr = fit.fTr;
+ fit.SetTrack(it->TLast, it->CLast);
+ const L1TrackPar& tr = fit.Tr();
CbmL1MCPoint& mcP = fvMCPoints[iMC];
fit.Extrapolate(mcP.zOut, fld, fvec::One());
@@ -1394,12 +1386,8 @@ void CbmL1::TrackFitPerformance()
{ // vertex
CbmL1MCTrack mc = *(it->GetMCTracks()[0]);
- {
- L1TrackPar trTmp(it->T, it->C);
- fit.fTr = trTmp;
- }
- fit.fQp0 = fit.fTr.qp;
- L1TrackPar& tr = fit.fTr;
+ fit.SetTrack(it->T, it->C);
+ L1TrackPar& tr = fit.Tr();
// if (mc.mother_ID != -1) { // secondary
if (!mc.IsPrimary()) { // secondary
@@ -1417,7 +1405,7 @@ void CbmL1::TrackFitPerformance()
&& (dir * (mc.z - fpAlgo->GetParameters()->GetStation(iSta).fZ[0]) > 0);
iSta += dir) {
// cout << iSta << " " << dir << endl;
- auto radThick = fpAlgo->GetParameters()->GetMaterialThickness(iSta, fit.fTr.x, fit.fTr.y);
+ auto radThick = fpAlgo->GetParameters()->GetMaterialThickness(iSta, fit.Tr().x, fit.Tr().y);
fit.AddMsInMaterial(radThick, fvec::One());
fit.EnergyLossCorrection(radThick, fvec::One(), fvec::One());
}
@@ -1471,10 +1459,9 @@ void CbmL1::TrackFitPerformance()
iSta += dir) {
fit.Extrapolate(fpAlgo->GetParameters()->GetStation(iSta).fZ, fld, fvec::One());
-
- fit.AddMsInMaterial(fpAlgo->GetParameters()->GetMaterialThickness(iSta, fit.fTr.x, fit.fTr.y), fvec::One());
- fit.EnergyLossCorrection(fpAlgo->GetParameters()->GetMaterialThickness(iSta, fit.fTr.x, fit.fTr.y),
- fvec::One(), fvec::One());
+ auto radThick = fpAlgo->GetParameters()->GetMaterialThickness(iSta, fit.Tr().x, fit.Tr().y);
+ fit.AddMsInMaterial(radThick, fvec::One());
+ fit.EnergyLossCorrection(radThick, fvec::One(), fvec::One());
}
fit.Extrapolate(mc.z, fld, fvec::One());
}
diff --git a/reco/L1/L1Algo/L1BranchExtender.cxx b/reco/L1/L1Algo/L1BranchExtender.cxx
index 28306b6535..aad4120a27 100644
--- a/reco/L1/L1Algo/L1BranchExtender.cxx
+++ b/reco/L1/L1Algo/L1BranchExtender.cxx
@@ -29,8 +29,9 @@ void L1Algo::BranchFitterFast(const L1Branch& t, L1TrackPar& Tout, const bool di
L1Fit fit;
fit.SetParticleMass(GetDefaultParticleMass());
- fit.fTr = Tout;
- L1TrackPar& T = fit.fTr;
+
+ fit.SetTrack(Tout);
+ L1TrackPar& T = fit.Tr();
// get hits of current track
const L1Vector<L1HitIndex_t>& hits = t.fHits; // array of indeses of hits of current track
@@ -77,7 +78,7 @@ void L1Algo::BranchFitterFast(const L1Branch& t, L1TrackPar& Tout, const bool di
T.ty = (y1 - y0) * dzi;
T.qp = qp0;
}
- fit.fQp0 = qp0;
+ fit.SetQp0(qp0);
T.z = z0;
T.t = hit0.t;
@@ -165,8 +166,8 @@ void L1Algo::FindMoreHits(L1Branch& t, L1TrackPar& Tout, const bool dir,
L1Fit fit;
fit.SetParticleMass(GetDefaultParticleMass());
- fit.fQp0 = qp0;
- fit.fTr = Tout;
+ fit.SetTrack(Tout);
+ fit.SetQp0(qp0);
const signed short int step = -2 * static_cast<int>(dir) + 1; // increment for station index
const int iFirstHit = (dir) ? 2 : t.NHits - 3;
@@ -223,12 +224,14 @@ void L1Algo::FindMoreHits(L1Branch& t, L1TrackPar& Tout, const bool dir,
fscal r2_best = 1e8; // best distance to hit
int iHit_best = -1; // index of the best hit
- const fscal iz = 1.f / (fit.fTr.z[0] - fParameters.GetTargetPositionZ()[0]);
+ L1TrackPar& tr = fit.Tr();
+
+ const fscal iz = 1.f / (tr.z[0] - fParameters.GetTargetPositionZ()[0]);
- L1HitAreaTime area(vGridTime[ista], fit.fTr.x[0] * iz, fit.fTr.y[0] * iz,
- (sqrt(fPickGather * fit.fTr.C00) + fMaxDx[ista] + fMaxDZ * abs(fit.fTr.tx))[0] * iz,
- (sqrt(fPickGather * fit.fTr.C11) + fMaxDy[ista] + fMaxDZ * abs(fit.fTr.ty))[0] * iz,
- fit.fTr.t[0], sqrt(fit.fTr.C55[0]));
+ L1HitAreaTime area(vGridTime[ista], tr.x[0] * iz, tr.y[0] * iz,
+ (sqrt(fPickGather * tr.C00) + fMaxDx[ista] + fMaxDZ * abs(tr.tx))[0] * iz,
+ (sqrt(fPickGather * tr.C11) + fMaxDy[ista] + fMaxDZ * abs(tr.ty))[0] * iz, tr.t[0],
+ sqrt(tr.C55[0]));
for (L1HitIndex_t ih = -1; true;) { // loop over the hits in the area
@@ -248,8 +251,8 @@ void L1Algo::FindMoreHits(L1Branch& t, L1TrackPar& Tout, const bool dir,
const L1Hit& hit = (*vHitsUnused)[globalInd];
if (sta.timeInfo) {
- fscal dt = hit.t - fit.fTr.t[0];
- if (dt * dt > (fit.fTr.C55[0] + hit.dt2) * 25) continue;
+ fscal dt = hit.t - tr.t[0];
+ if (dt * dt > (tr.C55[0] + hit.dt2) * 25) continue;
}
//if (GetFUsed((*fStripFlag)[hit.f] | (*fStripFlag)[hit.b])) continue; // if used
@@ -294,7 +297,7 @@ void L1Algo::FindMoreHits(L1Branch& t, L1TrackPar& Tout, const bool dir,
fit.Filter(sta.frontInfo, hit.u, hit.du2, fvec::One());
fit.Filter(sta.backInfo, hit.v, hit.dv2, fvec::One());
fit.FilterTime(hit.t, hit.dt2, fvec::One(), sta.timeInfo);
- auto radThick = fParameters.GetMaterialThickness(ista, fit.fTr.x, fit.fTr.y);
+ auto radThick = fParameters.GetMaterialThickness(ista, tr.x, tr.y);
fit.AddMsInMaterial(radThick, fvec::One());
fit.EnergyLossCorrection(radThick, dir ? fvec(1.f) : fvec(-1.f), fvec::One());
@@ -327,7 +330,7 @@ void L1Algo::FindMoreHits(L1Branch& t, L1TrackPar& Tout, const bool dir,
}
}
- Tout = fit.fTr;
+ Tout = fit.Tr();
} // void L1Algo::FindMoreHits
diff --git a/reco/L1/L1Algo/L1CATrackFinder.cxx b/reco/L1/L1Algo/L1CATrackFinder.cxx
index 9dd1c86d8a..63db51f673 100644
--- a/reco/L1/L1Algo/L1CATrackFinder.cxx
+++ b/reco/L1/L1Algo/L1CATrackFinder.cxx
@@ -214,7 +214,7 @@ inline void L1Algo::findSingletsStep1( /// input 1st stage of singlet search
L1Fit fit;
fit.SetParticleMass(GetDefaultParticleMass());
- fit.fQp0 = fMaxInvMom;
+ fit.SetQp0(fMaxInvMom);
if (fpCurrentIteration->GetElectronFlag()) { fit.SetParticleMass(L1Constants::phys::kElectronMass); }
else {
@@ -222,7 +222,6 @@ inline void L1Algo::findSingletsStep1( /// input 1st stage of singlet search
}
for (int i1_V = 0; i1_V < n1_V; i1_V++) {
- L1TrackPar& T = fit.fTr;
// field made by the left hit, the target and the station istac in-between.
// is used for extrapolation to the target and to the middle hit
@@ -259,6 +258,8 @@ inline void L1Algo::findSingletsStep1( /// input 1st stage of singlet search
fld0.Set(fTargB, fTargZ, b00, fld0Sta0.fZ, b01, fld0Sta1.fZ);
fld1.Set(b10, fld1Sta0.fZ, b11, fld1Sta1.fZ, b12, fld1Sta2.fZ);
+ L1TrackPar& T = fit.Tr();
+
T.chi2 = 0.;
T.NDF = (fpCurrentIteration->GetPrimaryFlag()) ? fvec(2.) : fvec(0.);
@@ -333,8 +334,9 @@ inline void L1Algo::findDoubletsStep0(
unsigned int Ndoublets = 0;
const Tindex i1_V = i1 / fvec::size();
const Tindex i1_4 = i1 % fvec::size();
- fit.fTr = T_1[i1_V];
- L1TrackPar& T1 = fit.fTr;
+
+ fit.SetTrack(T_1[i1_V]);
+ L1TrackPar& T1 = fit.Tr();
// assert(T1.IsEntryConsistent(true, i1_4));
// if (!T1.IsEntryConsistent(false, i1_4)) continue;
@@ -503,7 +505,7 @@ inline void L1Algo::findDoubletsStep0(
}
} // loop over the hits in the area
- T_1[i1_V] = fit.fTr;
+ T_1[i1_V] = fit.Tr();
} // for i1
}
@@ -542,7 +544,7 @@ inline void L1Algo::findTripletsStep0( // input
L1Fit fit;
fit.SetParticleMass(fDefaultMass);
- fit.fQp0 = fvec(0.);
+ fit.SetQp0(fvec(0.));
n3 = 0;
Tindex n3_V = 0, n3_4 = 0;
@@ -563,9 +565,9 @@ inline void L1Algo::findTripletsStep0( // input
// ---- Add the middle hits to parameters estimation. Propagate to right station. ----
for (Tindex i2 = 0; i2 < n2;) {
- L1TrackPar& T2 = fit.fTr;
- T2 = L1TrackPar_0;
- fit.fQp0 = fvec(0.);
+ fit.SetTrack(L1TrackPar_0);
+ fit.SetQp0(fvec(0.));
+ L1TrackPar& T2 = fit.Tr();
L1FieldRegion f2;
// pack the data
@@ -616,11 +618,11 @@ inline void L1Algo::findTripletsStep0( // input
fit.Filter(stam.backInfo, u_back_2, dv2_2, fvec::One());
fit.FilterTime(t_2, dt2_2, fvec::One(), stam.timeInfo);
- fit.fQp0 = isMomentumFitted ? fit.fTr.qp : fMaxInvMom;
+ fit.SetQp0(isMomentumFitted ? fit.Tr().qp : fMaxInvMom);
fit.AddMsInMaterial(fParameters.GetMaterialThickness(iStaM, T2.x, T2.y), fvec::One());
- fit.fQp0 = fit.fTr.qp;
+ fit.SetQp0(fit.Tr().qp);
// extrapolate to the right hit station
@@ -683,10 +685,8 @@ inline void L1Algo::findTripletsStep0( // input
L1Fit fit3;
fit3.SetParticleMass(fDefaultMass);
- fit3.fQp0 = T2.qp;
-
- L1TrackPar& T_cur = fit3.fTr;
- T_cur = T2;
+ fit3.SetTrack(T2);
+ L1TrackPar& T_cur = fit3.Tr();
fit3.ExtrapolateLine(zr, fvec::One());
@@ -811,18 +811,12 @@ inline void L1Algo::findTripletsStep1( // input
fit.SetParticleMass(fDefaultMass);
for (Tindex i3_V = 0; i3_V < n3_V; ++i3_V) {
-
- L1TrackPar& T3 = fit.fTr;
- T3 = T_3[i3_V];
- fit.fQp0 = fit.fTr.qp;
-
+ fit.SetTrack(T_3[i3_V]);
fit.ExtrapolateLine(z_Pos[i3_V], fvec::One());
-
fit.Filter(star.frontInfo, u_front_[i3_V], du2_3[i3_V], fvec::One());
fit.Filter(star.backInfo, u_back_[i3_V], dv2_3[i3_V], fvec::One());
-
if (kMcbm != fTrackingMode) { fit.FilterTime(t_3[i3_V], dt2_3[i3_V], fvec::One(), star.timeInfo); }
- T_3[i3_V] = T3;
+ T_3[i3_V] = fit.Tr();
}
}
@@ -916,7 +910,7 @@ inline void L1Algo::findTripletsStep2(Tindex n3, int istal, int istam, int istar
fld.Set(B[0], sta[0].fZ, B[1], sta[1].fZ, B[2], sta[2].fZ);
fldTarget.Set(fTargB, fTargZ, B[0], sta[0].fZ, B[1], sta[1].fZ);
- L1TrackPar& T = fit.fTr;
+ L1TrackPar& T = fit.Tr();
// initial parameters
{
@@ -930,9 +924,9 @@ inline void L1Algo::findTripletsStep2(Tindex n3, int istal, int istam, int istar
for (int iiter = 0; iiter < nIterations; ++iiter) {
// fit forward
{
- fit.fQp0 = T.qp;
- if (fit.fQp0[0] > GetMaxInvMom()) { fit.fQp0 = GetMaxInvMom(); }
- if (fit.fQp0[0] < -GetMaxInvMom()) { fit.fQp0 = -GetMaxInvMom(); }
+ fit.SetQp0(T.qp);
+ fit.Qp0()(fit.Qp0() > GetMaxInvMom()) = GetMaxInvMom();
+ fit.Qp0()(fit.Qp0() < -GetMaxInvMom()) = -GetMaxInvMom();
T.ResetErrors(200., 200., 1., 1., 100., 1.e4);
//T.ResetErrors(200., 200., 10., 10., 100., 1.e4);
@@ -968,9 +962,9 @@ inline void L1Algo::findTripletsStep2(Tindex n3, int istal, int istam, int istar
// fit backward
{
- fit.fQp0 = T.qp;
- if (fit.fQp0[0] > GetMaxInvMom()) { fit.fQp0 = GetMaxInvMom(); }
- if (fit.fQp0[0] < -GetMaxInvMom()) { fit.fQp0 = -GetMaxInvMom(); }
+ fit.SetQp0(T.qp);
+ fit.Qp0()(fit.Qp0() > GetMaxInvMom()) = GetMaxInvMom();
+ fit.Qp0()(fit.Qp0() < -GetMaxInvMom()) = -GetMaxInvMom();
T.ResetErrors(200., 200., 1., 1., 100., 1.e4);
T.NDF = ndf;
diff --git a/reco/L1/L1Algo/L1CloneMerger.cxx b/reco/L1/L1Algo/L1CloneMerger.cxx
index ed665648f9..ee511944f2 100644
--- a/reco/L1/L1Algo/L1CloneMerger.cxx
+++ b/reco/L1/L1Algo/L1CloneMerger.cxx
@@ -78,14 +78,14 @@ void L1CloneMerger::Exec(L1Vector<L1Track>& extTracks, L1Vector<L1HitIndex_t>& e
L1Fit fitB;
fitB.SetParticleMass(frAlgo.GetDefaultParticleMass());
- fitB.fQp0 = fvec(0.);
+ fitB.SetQp0(fvec(0.));
L1Fit fitF;
fitF.SetParticleMass(frAlgo.GetDefaultParticleMass());
- fitF.fQp0 = fvec(0.);
+ fitF.SetQp0(fvec(0.));
- L1TrackPar& Tb = fitB.fTr;
- L1TrackPar& Tf = fitF.fTr;
+ L1TrackPar& Tb = fitB.Tr();
+ L1TrackPar& Tf = fitF.Tr();
L1FieldValue fBm, fBb, fBf _fvecalignment;
L1FieldRegion fld _fvecalignment;
@@ -133,7 +133,7 @@ void L1CloneMerger::Exec(L1Vector<L1Track>& extTracks, L1Vector<L1HitIndex_t>& e
Tb.C54 = extTracks[iTr].CFirst[19];
Tb.C55 = extTracks[iTr].CFirst[20];
- fitB.fQp0 = fitB.fTr.qp;
+ fitB.SetQp0(fitB.Tr().qp);
unsigned short staf = lastStation[jTr];
@@ -166,7 +166,7 @@ void L1CloneMerger::Exec(L1Vector<L1Track>& extTracks, L1Vector<L1HitIndex_t>& e
Tf.C54 = extTracks[jTr].CLast[19];
Tf.C55 = extTracks[jTr].CLast[20];
- fitF.fQp0 = fitF.fTr.qp;
+ fitF.SetQp0(fitF.Tr().qp);
if (fabs(Tf.t[0] - Tb.t[0]) > 3 * sqrt(Tf.C55[0] + Tb.C55[0])) continue;
unsigned short stam;
diff --git a/reco/L1/L1Algo/L1Fit.h b/reco/L1/L1Algo/L1Fit.h
index c057cf94e1..48cad6a908 100644
--- a/reco/L1/L1Algo/L1Fit.h
+++ b/reco/L1/L1Algo/L1Fit.h
@@ -1,6 +1,12 @@
/* Copyright (C) 2017-2021 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
SPDX-License-Identifier: GPL-3.0-only
- Authors: Maksym Zyzak [committer], Sergey Gorbunov */
+ Authors: Sergey Gorbunov [committer], Maksym Zyzak */
+
+
+/// \file L1Fit.h
+/// \brief Kalman filter fitter for L1 tracking
+/// \since 10.02.2023
+/// \author S.Gorbunov
#ifndef L1Fit_h
#define L1Fit_h
@@ -11,6 +17,8 @@
#include "L1UMeasurementInfo.h"
#include "L1XYMeasurementInfo.h"
+/// Kalman filter fitter for L1 tracking
+///
class L1Fit {
public:
@@ -37,6 +45,12 @@ public:
fQp0 = fTr.qp;
}
+ void SetQp0(const fvec& qp0) { fQp0 = qp0; }
+
+ L1TrackPar& Tr() { return fTr; }
+
+ fvec& Qp0() { return fQp0; }
+
void SetOneEntry(const int i0, const L1Fit& T1, const int i1);
void Print(int i = -1);
@@ -107,12 +121,12 @@ public:
static void FilterChi2(const L1UMeasurementInfo& info, const fvec& x, const fvec& y, const fvec& C00, const fvec& C10,
const fvec& C11, fvec& chi2, const fvec& u, const fvec& du2);
-public:
+private:
L1TrackPar fTr {};
- fvec fQp0;
+ fvec fQp0 {0.};
- fvec fMass = 0.10565800; // muon mass
- fvec fMass2 = fMass * fMass; // mass squared
+ fvec fMass {0.10565800}; // muon mass
+ fvec fMass2 {fMass * fMass}; // mass squared
fvec fPipeRadThick {7.87e-3f}; // 0.7 mm Aluminium // TODO:
fvec fTargetRadThick {3.73e-2f * 2}; // 250 mum Gold // TODO:
diff --git a/reco/L1/L1Algo/L1TrackFitter.cxx b/reco/L1/L1Algo/L1TrackFitter.cxx
index c3bf044a6c..6029484410 100644
--- a/reco/L1/L1Algo/L1TrackFitter.cxx
+++ b/reco/L1/L1Algo/L1TrackFitter.cxx
@@ -42,7 +42,7 @@ void L1Algo::L1KFTrackFitter()
int nTracks_SIMD = fvec::size();
L1Fit fit; // fitting parametr coresponding to current track
- L1TrackPar& tr = fit.fTr;
+ L1TrackPar& tr = fit.Tr();
fit.SetParticleMass(GetDefaultParticleMass());
@@ -219,7 +219,7 @@ void L1Algo::L1KFTrackFitter()
for (int iter = 0; iter < 2; iter++) { // 1.5 iterations
- fit.fQp0 = tr.qp;
+ fit.SetQp0(tr.qp);
// fit backward
@@ -285,21 +285,21 @@ void L1Algo::L1KFTrackFitter()
if (kGlobal == fTrackingMode) {
for (int vtxIter = 0; vtxIter < 2; vtxIter++) {
- fitpv.fQp0 = fitpv.fTr.qp;
- fitpv.fTr = fit.fTr;
- fitpv.fTr.qp = fitpv.fQp0;
+ fitpv.SetQp0(fitpv.Tr().qp);
+ fitpv.Tr() = fit.Tr();
+ fitpv.Tr().qp = fitpv.Qp0();
fitpv.Extrapolate(fParameters.GetTargetPositionZ(), fld, fvec(1.));
fitpv.Filter(vtxInfoX, fParameters.GetTargetPositionX(), fvec(1.e-8), fvec::One());
fitpv.Filter(vtxInfoY, fParameters.GetTargetPositionY(), fvec(1.e-8), fvec::One());
}
}
else {
- fitpv.fQp0 = fitpv.fTr.qp;
+ fitpv.SetQp0(fitpv.Tr().qp);
fitpv.Extrapolate(fParameters.GetTargetPositionZ(), fld, fvec(1.));
}
}
- const L1TrackPar& Tf = fit.fTr;
+ const L1TrackPar& Tf = fit.Tr();
for (int iVec = 0; iVec < nTracks_SIMD; iVec++) {
t[iVec]->TFirst[0] = Tf.x[iVec];
@@ -307,7 +307,7 @@ void L1Algo::L1KFTrackFitter()
t[iVec]->TFirst[2] = Tf.tx[iVec];
t[iVec]->TFirst[3] = Tf.ty[iVec];
t[iVec]->TFirst[4] = Tf.qp[iVec];
- if (kGlobal == fTrackingMode) { t[iVec]->TFirst[4] = fitpv.fTr.qp[iVec]; }
+ if (kGlobal == fTrackingMode) { t[iVec]->TFirst[4] = fitpv.Tr().qp[iVec]; }
t[iVec]->TFirst[5] = Tf.z[iVec];
t[iVec]->TFirst[6] = Tf.t[iVec];
@@ -338,35 +338,35 @@ void L1Algo::L1KFTrackFitter()
}
for (int iVec = 0; iVec < nTracks_SIMD; iVec++) {
- t[iVec]->Tpv[0] = fitpv.fTr.x[iVec];
- t[iVec]->Tpv[1] = fitpv.fTr.y[iVec];
- t[iVec]->Tpv[2] = fitpv.fTr.tx[iVec];
- t[iVec]->Tpv[3] = fitpv.fTr.ty[iVec];
- t[iVec]->Tpv[4] = fitpv.fTr.qp[iVec];
- t[iVec]->Tpv[5] = fitpv.fTr.z[iVec];
- t[iVec]->Tpv[6] = fitpv.fTr.t[iVec];
-
- t[iVec]->Cpv[0] = fitpv.fTr.C00[iVec];
- t[iVec]->Cpv[1] = fitpv.fTr.C10[iVec];
- t[iVec]->Cpv[2] = fitpv.fTr.C11[iVec];
- t[iVec]->Cpv[3] = fitpv.fTr.C20[iVec];
- t[iVec]->Cpv[4] = fitpv.fTr.C21[iVec];
- t[iVec]->Cpv[5] = fitpv.fTr.C22[iVec];
- t[iVec]->Cpv[6] = fitpv.fTr.C30[iVec];
- t[iVec]->Cpv[7] = fitpv.fTr.C31[iVec];
- t[iVec]->Cpv[8] = fitpv.fTr.C32[iVec];
- t[iVec]->Cpv[9] = fitpv.fTr.C33[iVec];
- t[iVec]->Cpv[10] = fitpv.fTr.C40[iVec];
- t[iVec]->Cpv[11] = fitpv.fTr.C41[iVec];
- t[iVec]->Cpv[12] = fitpv.fTr.C42[iVec];
- t[iVec]->Cpv[13] = fitpv.fTr.C43[iVec];
- t[iVec]->Cpv[14] = fitpv.fTr.C44[iVec];
- t[iVec]->Cpv[15] = fitpv.fTr.C50[iVec];
- t[iVec]->Cpv[16] = fitpv.fTr.C51[iVec];
- t[iVec]->Cpv[17] = fitpv.fTr.C52[iVec];
- t[iVec]->Cpv[18] = fitpv.fTr.C53[iVec];
- t[iVec]->Cpv[19] = fitpv.fTr.C54[iVec];
- t[iVec]->Cpv[20] = fitpv.fTr.C55[iVec];
+ t[iVec]->Tpv[0] = fitpv.Tr().x[iVec];
+ t[iVec]->Tpv[1] = fitpv.Tr().y[iVec];
+ t[iVec]->Tpv[2] = fitpv.Tr().tx[iVec];
+ t[iVec]->Tpv[3] = fitpv.Tr().ty[iVec];
+ t[iVec]->Tpv[4] = fitpv.Tr().qp[iVec];
+ t[iVec]->Tpv[5] = fitpv.Tr().z[iVec];
+ t[iVec]->Tpv[6] = fitpv.Tr().t[iVec];
+
+ t[iVec]->Cpv[0] = fitpv.Tr().C00[iVec];
+ t[iVec]->Cpv[1] = fitpv.Tr().C10[iVec];
+ t[iVec]->Cpv[2] = fitpv.Tr().C11[iVec];
+ t[iVec]->Cpv[3] = fitpv.Tr().C20[iVec];
+ t[iVec]->Cpv[4] = fitpv.Tr().C21[iVec];
+ t[iVec]->Cpv[5] = fitpv.Tr().C22[iVec];
+ t[iVec]->Cpv[6] = fitpv.Tr().C30[iVec];
+ t[iVec]->Cpv[7] = fitpv.Tr().C31[iVec];
+ t[iVec]->Cpv[8] = fitpv.Tr().C32[iVec];
+ t[iVec]->Cpv[9] = fitpv.Tr().C33[iVec];
+ t[iVec]->Cpv[10] = fitpv.Tr().C40[iVec];
+ t[iVec]->Cpv[11] = fitpv.Tr().C41[iVec];
+ t[iVec]->Cpv[12] = fitpv.Tr().C42[iVec];
+ t[iVec]->Cpv[13] = fitpv.Tr().C43[iVec];
+ t[iVec]->Cpv[14] = fitpv.Tr().C44[iVec];
+ t[iVec]->Cpv[15] = fitpv.Tr().C50[iVec];
+ t[iVec]->Cpv[16] = fitpv.Tr().C51[iVec];
+ t[iVec]->Cpv[17] = fitpv.Tr().C52[iVec];
+ t[iVec]->Cpv[18] = fitpv.Tr().C53[iVec];
+ t[iVec]->Cpv[19] = fitpv.Tr().C54[iVec];
+ t[iVec]->Cpv[20] = fitpv.Tr().C55[iVec];
}
if (iter == 1) continue; // only 1.5 iterations
@@ -377,7 +377,7 @@ void L1Algo::L1KFTrackFitter()
FilterFirst(fit, x_first, y_first, time_first, dt2_first, d_xx_fst, d_yy_fst, d_xy_fst);
- fit.fQp0 = tr.qp;
+ fit.SetQp0(tr.qp);
fldZ1 = z[ista];
sta[ista].fieldSlice.GetFieldValue(tr.x, tr.y, fldB1);
@@ -420,7 +420,7 @@ void L1Algo::L1KFTrackFitter()
t[iVec]->TLast[2] = tr.tx[iVec];
t[iVec]->TLast[3] = tr.ty[iVec];
t[iVec]->TLast[4] = tr.qp[iVec];
- if (kGlobal == fTrackingMode) { t[iVec]->TLast[4] = fitpv.fTr.qp[iVec]; }
+ if (kGlobal == fTrackingMode) { t[iVec]->TLast[4] = fitpv.Tr().qp[iVec]; }
t[iVec]->TLast[5] = tr.z[iVec];
t[iVec]->TLast[6] = tr.t[iVec];
@@ -452,33 +452,36 @@ void L1Algo::L1KFTrackFitter()
} // iter
}
}
-void L1Algo::GuessVecNoField(L1Fit& t, fvec& x_last, fvec& x_2last, fvec& y_last, fvec& y_2last, fvec& z_end,
+void L1Algo::GuessVecNoField(L1Fit& track, fvec& x_last, fvec& x_2last, fvec& y_last, fvec& y_2last, fvec& z_end,
fvec& z_2last, fvec& time_last, fvec* /*w_time*/, fvec& dt2_last)
{
+ L1TrackPar& tr = track.Tr();
+
fvec dzi = fvec(1.) / (z_end - z_2last);
- t.fTr.x = x_last;
- t.fTr.y = y_last;
- t.fTr.tx = (x_last - x_2last) * dzi;
- t.fTr.ty = (y_last - y_2last) * dzi;
- t.fTr.t = time_last;
- t.fTr.z = z_end;
- t.fTr.chi2 = fvec(0.);
- t.fTr.NDF = fvec(2.);
-
- t.fTr.C20 = t.fTr.C21 = fvec(0.);
- t.fTr.C30 = t.fTr.C31 = t.fTr.C32 = fvec(0.);
- t.fTr.C40 = t.fTr.C41 = t.fTr.C42 = t.fTr.C43 = fvec(0.);
- t.fTr.C50 = t.fTr.C51 = t.fTr.C52 = t.fTr.C53 = t.fTr.C54 = fvec(0.);
- t.fTr.C22 = t.fTr.C33 = vINF;
- t.fTr.C44 = fvec(1.);
- t.fTr.C55 = fvec(dt2_last);
+ tr.x = x_last;
+ tr.y = y_last;
+ tr.tx = (x_last - x_2last) * dzi;
+ tr.ty = (y_last - y_2last) * dzi;
+ tr.t = time_last;
+ tr.z = z_end;
+ tr.chi2 = fvec(0.);
+ tr.NDF = fvec(2.);
+
+ tr.C20 = tr.C21 = fvec(0.);
+ tr.C30 = tr.C31 = tr.C32 = fvec(0.);
+ tr.C40 = tr.C41 = tr.C42 = tr.C43 = fvec(0.);
+ tr.C50 = tr.C51 = tr.C52 = tr.C53 = tr.C54 = fvec(0.);
+ tr.C22 = tr.C33 = vINF;
+ tr.C44 = fvec(1.);
+ tr.C55 = fvec(dt2_last);
}
-void L1Algo::GuessVec(L1TrackPar& t, fvec* xV, fvec* yV, fvec* zV, fvec* Sy, fvec* wV, int NHits, fvec* zCur)
+void L1Algo::GuessVec(L1TrackPar& tr, fvec* xV, fvec* yV, fvec* zV, fvec* Sy, fvec* wV, int NHits, fvec* zCur)
// gives nice initial approximation for x,y,tx,ty - almost same as KF fit. qp - is shifted by 4%, resid_ual - ~3.5% (KF fit resid_ual - 1%).
{
+
fvec A0, A1 = ZERO, A2 = ZERO, A3 = ZERO, A4 = ZERO, A5 = ZERO, a0, a1 = ZERO, a2 = ZERO, b0, b1 = ZERO, b2 = ZERO;
fvec z0, x, y, z, S, w, wz, wS;
@@ -527,26 +530,28 @@ void L1Algo::GuessVec(L1TrackPar& t, fvec* xV, fvec* yV, fvec* zV, fvec* Sy, fve
fvec Ai5 = (-A1 * A1 + A0 * A2);
fvec L, L1;
- t.x = (Ai0 * a0 + Ai1 * a1 + Ai3 * a2) * det;
- t.tx = (Ai1 * a0 + Ai2 * a1 + Ai4 * a2) * det;
- fvec txtx1 = fvec(1.) + t.tx * t.tx;
+ tr.x = (Ai0 * a0 + Ai1 * a1 + Ai3 * a2) * det;
+ tr.tx = (Ai1 * a0 + Ai2 * a1 + Ai4 * a2) * det;
+ fvec txtx1 = fvec(1.) + tr.tx * tr.tx;
L = (Ai3 * a0 + Ai4 * a1 + Ai5 * a2) * det / txtx1;
- L1 = L * t.tx;
+ L1 = L * tr.tx;
A1 = A1 + A3 * L1;
A2 = A2 + (A4 + A4 + A5 * L1) * L1;
b1 += b2 * L1;
det = fvec::One() / (-A1 * A1 + A0 * A2);
- t.y = (A2 * b0 - A1 * b1) * det;
- t.ty = (-A1 * b0 + A0 * b1) * det;
- t.qp = -L * c_light_i / sqrt(txtx1 + t.ty * t.ty);
- t.z = z0;
+ tr.y = (A2 * b0 - A1 * b1) * det;
+ tr.ty = (-A1 * b0 + A0 * b1) * det;
+ tr.qp = -L * c_light_i / sqrt(txtx1 + tr.ty * tr.ty);
+ tr.z = z0;
}
-void L1Algo::GuessVec(L1Fit& t, fvec* xV, fvec* yV, fvec* zV, fvec* Sy, fvec* wV, int NHits, fvec* zCur, fvec* timeV,
- fvec* w_time)
+void L1Algo::GuessVec(L1Fit& track, fvec* xV, fvec* yV, fvec* zV, fvec* Sy, fvec* wV, int NHits, fvec* zCur,
+ fvec* timeV, fvec* w_time)
// gives nice initial approximation for x,y,tx,ty - almost same as KF fit. qp - is shifted by 4%, resid_ual - ~3.5% (KF fit resid_ual - 1%).
{
+ L1TrackPar& tr = track.Tr();
+
fvec A0, A1 = ZERO, A2 = ZERO, A3 = ZERO, A4 = ZERO, A5 = ZERO, a0, a1 = ZERO, a2 = ZERO, b0, b1 = ZERO, b2 = ZERO;
fvec z0, x, y, z, S, w, wz, wS, time;
@@ -610,51 +615,53 @@ void L1Algo::GuessVec(L1Fit& t, fvec* xV, fvec* yV, fvec* zV, fvec* Sy, fvec* wV
fvec Ai5 = (-A1 * A1 + A0 * A2);
fvec L, L1;
- t.fTr.x = (Ai0 * a0 + Ai1 * a1 + Ai3 * a2) * det;
- t.fTr.tx = (Ai1 * a0 + Ai2 * a1 + Ai4 * a2) * det;
- fvec txtx1 = fvec(1.) + t.fTr.tx * t.fTr.tx;
+ tr.x = (Ai0 * a0 + Ai1 * a1 + Ai3 * a2) * det;
+ tr.tx = (Ai1 * a0 + Ai2 * a1 + Ai4 * a2) * det;
+ fvec txtx1 = fvec(1.) + tr.tx * tr.tx;
L = (Ai3 * a0 + Ai4 * a1 + Ai5 * a2) * det / txtx1;
- L1 = L * t.fTr.tx;
+ L1 = L * tr.tx;
A1 = A1 + A3 * L1;
A2 = A2 + (A4 + A4 + A5 * L1) * L1;
b1 += b2 * L1;
- det = fvec::One() / (-A1 * A1 + A0 * A2);
- t.fTr.y = (A2 * b0 - A1 * b1) * det;
- t.fTr.ty = (-A1 * b0 + A0 * b1) * det;
- t.fTr.qp = -L * c_light_i / sqrt(txtx1 + t.fTr.ty * t.fTr.ty);
- if (timeV) { t.fTr.t = iif(nhits > fvec::Zero(), time / nhits, fvec::Zero()); }
+ det = fvec::One() / (-A1 * A1 + A0 * A2);
+ tr.y = (A2 * b0 - A1 * b1) * det;
+ tr.ty = (-A1 * b0 + A0 * b1) * det;
+ tr.qp = -L * c_light_i / sqrt(txtx1 + tr.ty * tr.ty);
+ if (timeV) { tr.t = iif(nhits > fvec::Zero(), time / nhits, fvec::Zero()); }
- t.fTr.z = z0;
+ tr.z = z0;
}
void L1Algo::FilterFirst(L1Fit& track, fvec& x, fvec& y, fvec& t, fvec& dt2, fvec& d_xx, fvec& d_yy, fvec& d_xy)
{
- track.fTr.C00 = d_xx;
- track.fTr.C10 = d_xy;
- track.fTr.C11 = d_yy;
- track.fTr.C20 = ZERO;
- track.fTr.C21 = ZERO;
- track.fTr.C22 = vINF;
- track.fTr.C30 = ZERO;
- track.fTr.C31 = ZERO;
- track.fTr.C32 = ZERO;
- track.fTr.C33 = vINF;
- track.fTr.C40 = ZERO;
- track.fTr.C41 = ZERO;
- track.fTr.C42 = ZERO;
- track.fTr.C43 = ZERO;
- track.fTr.C44 = ONE;
- track.fTr.C50 = ZERO;
- track.fTr.C51 = ZERO;
- track.fTr.C52 = ZERO;
- track.fTr.C53 = ZERO;
- track.fTr.C54 = ZERO;
- track.fTr.C55 = dt2;
-
- track.fTr.x = x;
- track.fTr.y = y;
- track.fTr.t = t;
- track.fTr.NDF = -3.0;
- track.fTr.chi2 = ZERO;
+ L1TrackPar& tr = track.Tr();
+
+ tr.C00 = d_xx;
+ tr.C10 = d_xy;
+ tr.C11 = d_yy;
+ tr.C20 = ZERO;
+ tr.C21 = ZERO;
+ tr.C22 = vINF;
+ tr.C30 = ZERO;
+ tr.C31 = ZERO;
+ tr.C32 = ZERO;
+ tr.C33 = vINF;
+ tr.C40 = ZERO;
+ tr.C41 = ZERO;
+ tr.C42 = ZERO;
+ tr.C43 = ZERO;
+ tr.C44 = ONE;
+ tr.C50 = ZERO;
+ tr.C51 = ZERO;
+ tr.C52 = ZERO;
+ tr.C53 = ZERO;
+ tr.C54 = ZERO;
+ tr.C55 = dt2;
+
+ tr.x = x;
+ tr.y = y;
+ tr.t = t;
+ tr.NDF = -3.0;
+ tr.chi2 = ZERO;
}
diff --git a/reco/L1/ParticleFinder/CbmL1PFFitter.cxx b/reco/L1/ParticleFinder/CbmL1PFFitter.cxx
index ff5f029fe6..36bb7140a7 100644
--- a/reco/L1/ParticleFinder/CbmL1PFFitter.cxx
+++ b/reco/L1/ParticleFinder/CbmL1PFFitter.cxx
@@ -99,22 +99,23 @@ inline CbmL1PFFitter::PFFieldRegion::PFFieldRegion(const L1FieldRegion& fld, int
void FilterFirst(L1Fit& fit, fvec& x, fvec& y, fvec& dxx, fvec& dxy, fvec& dyy)
{
- L1TrackPar& tr = fit.fTr;
- tr.C00 = dxx;
- tr.C10 = dxy;
- tr.C11 = dyy;
- tr.C20 = ZERO;
- tr.C21 = ZERO;
- tr.C22 = vINF;
- tr.C30 = ZERO;
- tr.C31 = ZERO;
- tr.C32 = ZERO;
- tr.C33 = vINF;
- tr.C40 = ZERO;
- tr.C41 = ZERO;
- tr.C42 = ZERO;
- tr.C43 = ZERO;
- tr.C44 = ONE;
+ L1TrackPar& tr = fit.Tr();
+
+ tr.C00 = dxx;
+ tr.C10 = dxy;
+ tr.C11 = dyy;
+ tr.C20 = ZERO;
+ tr.C21 = ZERO;
+ tr.C22 = vINF;
+ tr.C30 = ZERO;
+ tr.C31 = ZERO;
+ tr.C32 = ZERO;
+ tr.C33 = vINF;
+ tr.C40 = ZERO;
+ tr.C41 = ZERO;
+ tr.C42 = ZERO;
+ tr.C43 = ZERO;
+ tr.C44 = ONE;
tr.x = x;
tr.y = y;
@@ -142,7 +143,7 @@ void CbmL1PFFitter::Fit(vector<CbmStsTrack>& Tracks, vector<int>& pidHypo)
L1Fit fit;
fit.SetParticleMass(CbmL1::Instance()->fpAlgo->GetDefaultParticleMass());
- L1TrackPar& T = fit.fTr; // fitting parametr coresponding to current track
+ L1TrackPar& T = fit.Tr(); // fitting parametr coresponding to current track
CbmStsTrack* tr[fvec::size()] {nullptr};
@@ -302,7 +303,7 @@ void CbmL1PFFitter::Fit(vector<CbmStsTrack>& Tracks, vector<int>& pidHypo)
i = 0;
FilterFirst(fit, x_first, y_first, fstC00, fstC10, fstC11);
- fit.fQp0 = fit.fTr.qp;
+ fit.SetQp0(fit.Tr().qp);
z1 = z[i];
sta[i].fieldSlice.GetFieldValue(T.x, T.y, b1);
@@ -324,7 +325,7 @@ void CbmL1PFFitter::Fit(vector<CbmStsTrack>& Tracks, vector<int>& pidHypo)
fvec wIn = iif(initialised, fvec::One(), fvec::Zero());
fit.Extrapolate(z[i], fld, w1);
- auto radThick = CbmL1::Instance()->fpAlgo->GetParameters()->GetMaterialThickness(i, fit.fTr.x, fit.fTr.y);
+ auto radThick = CbmL1::Instance()->fpAlgo->GetParameters()->GetMaterialThickness(i, fit.Tr().x, fit.Tr().y);
fit.Filter(sta[i].frontInfo, u[i], du2[i], w1);
fit.Filter(sta[i].backInfo, v[i], dv2[i], w1);
fit.FilterTime(t[i], dt2[i], w1, sta[i].timeInfo);
@@ -366,7 +367,7 @@ void CbmL1PFFitter::Fit(vector<CbmStsTrack>& Tracks, vector<int>& pidHypo)
}
//fit backward
- fit.fQp0 = T.qp;
+ fit.SetQp0(T.qp);
i = nHits - 1;
@@ -393,7 +394,7 @@ void CbmL1PFFitter::Fit(vector<CbmStsTrack>& Tracks, vector<int>& pidHypo)
fvec wIn = iif(initialised, fvec::One(), fvec::Zero());
fit.Extrapolate(z[i], fld, w1);
- auto radThick = CbmL1::Instance()->fpAlgo->GetParameters()->GetMaterialThickness(i, fit.fTr.x, fit.fTr.y);
+ auto radThick = CbmL1::Instance()->fpAlgo->GetParameters()->GetMaterialThickness(i, fit.Tr().x, fit.Tr().y);
fit.Filter(sta[i].frontInfo, u[i], du2[i], w1);
fit.Filter(sta[i].backInfo, v[i], dv2[i], w1);
fit.FilterTime(t[i], dt2[i], w1, sta[i].timeInfo);
@@ -454,7 +455,7 @@ void CbmL1PFFitter::GetChiToVertex(vector<CbmStsTrack>& Tracks, vector<PFFieldRe
int nTracks_SIMD = fvec::size();
L1Fit fit;
- L1TrackPar& T = fit.fTr; // fitting parametr coresponding to current track
+ L1TrackPar& T = fit.Tr(); // fitting parametr coresponding to current track
CbmStsTrack* tr[fvec::size()] {nullptr};
@@ -550,14 +551,14 @@ void CbmL1PFFitter::GetChiToVertex(vector<CbmStsTrack>& Tracks, vector<PFFieldRe
field.emplace_back(fld, i);
}
- fit.fQp0 = fit.fTr.qp;
+ fit.SetQp0(fit.Tr().qp);
for (int iSt = nStations - 4; iSt >= 0; iSt--) {
fvec w = iif(T.z > (zSta[iSt] + fvec(2.5)), fvec::One(), fvec::Zero());
fit.Extrapolate(zSta[iSt], fld, w);
- auto radThick = CbmL1::Instance()->fpAlgo->GetParameters()->GetMaterialThickness(iSt, fit.fTr.x, fit.fTr.y);
+ auto radThick = CbmL1::Instance()->fpAlgo->GetParameters()->GetMaterialThickness(iSt, fit.Tr().x, fit.Tr().y);
fit.AddMsInMaterial(radThick, w);
fit.EnergyLossCorrection(radThick, fvec::One(), w);
}
--
GitLab