From 33d4ad32036d4bd68bd982510e7e5e83d0ab636b Mon Sep 17 00:00:00 2001
From: "se.gorbunov" <se.gorbunov@gsi.de>
Date: Sat, 11 Feb 2023 00:57:24 +0000
Subject: [PATCH] L1 fit: improve L1Fit interface
---
reco/L1/CbmL1Performance.cxx | 19 +-
reco/L1/L1Algo/L1Algo.h | 10 +-
reco/L1/L1Algo/L1BranchExtender.cxx | 29 +-
reco/L1/L1Algo/L1CATrackFinder.cxx | 71 +--
reco/L1/L1Algo/L1CloneMerger.cxx | 6 +-
reco/L1/L1Algo/L1Field.h | 15 +-
reco/L1/L1Algo/L1Fit.cxx | 557 ++++++++++++-----------
reco/L1/L1Algo/L1Fit.h | 96 ++--
reco/L1/L1Algo/L1TrackFitter.cxx | 88 ++--
reco/L1/ParticleFinder/CbmL1PFFitter.cxx | 58 +--
10 files changed, 491 insertions(+), 458 deletions(-)
diff --git a/reco/L1/CbmL1Performance.cxx b/reco/L1/CbmL1Performance.cxx
index 0c06f06425..956f31f2f0 100644
--- a/reco/L1/CbmL1Performance.cxx
+++ b/reco/L1/CbmL1Performance.cxx
@@ -1103,6 +1103,7 @@ void CbmL1::TrackFitPerformance()
L1Fit fit;
fit.SetParticleMass(fpAlgo->GetDefaultParticleMass());
+ fit.SetMask(fmask::One());
if (first_call) {
first_call = 0;
@@ -1236,7 +1237,7 @@ void CbmL1::TrackFitPerformance()
CbmL1MCPoint& mcP = fvMCPoints[mc.Points[0]];
- fit.Extrapolate(mcP.zIn, fld, fvec::One());
+ fit.Extrapolate(mcP.zIn, fld);
double dx = tr.x[0] - mcP.xIn;
double dy = tr.y[0] - mcP.yIn;
@@ -1338,7 +1339,7 @@ void CbmL1::TrackFitPerformance()
const L1TrackPar& tr = fit.Tr();
CbmL1MCPoint& mcP = fvMCPoints[iMC];
- fit.Extrapolate(mcP.zOut, fld, fvec::One());
+ fit.Extrapolate(mcP.zOut, fld);
h_fitL[0]->Fill((tr.x[0] - mcP.xOut) * 1.e4);
h_fitL[1]->Fill((tr.y[0] - mcP.yOut) * 1.e4);
@@ -1395,7 +1396,7 @@ void CbmL1::TrackFitPerformance()
{ // extrapolate to vertex
L1FieldRegion fld _fvecalignment;
fld.SetUseOriginalField();
- fit.Extrapolate(mc.z, fld, fvec::One());
+ fit.Extrapolate(mc.z, fld);
// add material
const int fSta = fvHitStore[it->Hits[0]].iStation;
const int dir = int((mc.z - fpAlgo->GetParameters()->GetStation(fSta).fZ[0])
@@ -1406,8 +1407,8 @@ void CbmL1::TrackFitPerformance()
iSta += dir) {
// cout << iSta << " " << dir << endl;
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.AddMsInMaterial(radThick);
+ fit.EnergyLossCorrection(radThick, fvec::One());
}
}
if (mc.z != tr.z[0]) continue;
@@ -1458,12 +1459,12 @@ void CbmL1::TrackFitPerformance()
&& (dir * (mc.z - fpAlgo->GetParameters()->GetStation(iSta).fZ[0]) > 0);
iSta += dir) {
- fit.Extrapolate(fpAlgo->GetParameters()->GetStation(iSta).fZ, fld, fvec::One());
+ fit.Extrapolate(fpAlgo->GetParameters()->GetStation(iSta).fZ, fld);
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.AddMsInMaterial(radThick);
+ fit.EnergyLossCorrection(radThick, fvec::One());
}
- fit.Extrapolate(mc.z, fld, fvec::One());
+ fit.Extrapolate(mc.z, fld);
}
if (mc.z != tr.z[0]) continue;
diff --git a/reco/L1/L1Algo/L1Algo.h b/reco/L1/L1Algo/L1Algo.h
index 7d6a8b9682..5a8d7177cf 100644
--- a/reco/L1/L1Algo/L1Algo.h
+++ b/reco/L1/L1Algo/L1Algo.h
@@ -339,11 +339,11 @@ public:
/// ------ Subroutines used by L1Algo::KFTrackFitter() ------
- void GuessVec(L1TrackPar& t, fvec* xV, fvec* yV, fvec* zV, fvec* Sy, fvec* wV, int NHits, fvec* zCur = 0);
- void GuessVec(L1Fit& t, fvec* xV, fvec* yV, fvec* zV, fvec* Sy, fvec* wV, int NHits, fvec* zCur = 0, fvec* timeV = 0,
- fvec* w_time = 0);
+ void GuessVec(L1TrackPar& t, fvec* xV, fvec* yV, fvec* zV, fvec* Sy, fmask* wV, int NHits, fvec* zCur = 0);
+ void GuessVec(L1Fit& t, fvec* xV, fvec* yV, fvec* zV, fvec* Sy, fmask* wV, int NHits, fvec* zCur = 0, fvec* timeV = 0,
+ fmask* w_time = 0);
void GuessVecNoField(L1Fit& t, 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);
+ fvec& time_last, fmask* w_time, fvec& dt2_last);
void FilterFirst(L1Fit& track, fvec& x, fvec& y, fvec& t, fvec& dt2, fvec& dx2, fvec& dy2, fvec& dxy);
@@ -371,7 +371,7 @@ public:
/// Track fitting procedures
- void L1KFTrackFitter(); // version from SIMD-KF benchmark
+ void L1KFTrackFitter(); // version from SIMD-KF benchmark
float GetMaxInvMom() const { return fMaxInvMom[0]; }
diff --git a/reco/L1/L1Algo/L1BranchExtender.cxx b/reco/L1/L1Algo/L1BranchExtender.cxx
index aad4120a27..64beea3659 100644
--- a/reco/L1/L1Algo/L1BranchExtender.cxx
+++ b/reco/L1/L1Algo/L1BranchExtender.cxx
@@ -29,7 +29,7 @@ void L1Algo::BranchFitterFast(const L1Branch& t, L1TrackPar& Tout, const bool di
L1Fit fit;
fit.SetParticleMass(GetDefaultParticleMass());
-
+ fit.SetMask(fmask::One());
fit.SetTrack(Tout);
L1TrackPar& T = fit.Tr();
@@ -120,13 +120,13 @@ void L1Algo::BranchFitterFast(const L1Branch& t, L1TrackPar& Tout, const bool di
const L1Station& sta = fParameters.GetStation(ista);
- fit.Extrapolate(hit.z, fld, fvec::One());
- 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);
+ fit.Extrapolate(hit.z, fld);
+ fit.Filter(sta.frontInfo, hit.u, hit.du2);
+ fit.Filter(sta.backInfo, hit.v, hit.dv2);
+ fit.FilterTime(hit.t, hit.dt2, sta.timeInfo);
auto radThick = fParameters.GetMaterialThickness(ista, T.x, T.y);
- fit.AddMsInMaterial(radThick, fvec::One());
- fit.EnergyLossCorrection(radThick, fvec(-1.f), fvec::One());
+ fit.AddMsInMaterial(radThick);
+ fit.EnergyLossCorrection(radThick, fvec(-1.f));
fldB0 = fldB1;
fldB1 = fldB2;
@@ -166,6 +166,7 @@ void L1Algo::FindMoreHits(L1Branch& t, L1TrackPar& Tout, const bool dir,
L1Fit fit;
fit.SetParticleMass(GetDefaultParticleMass());
+ fit.SetMask(fmask::One());
fit.SetTrack(Tout);
fit.SetQp0(qp0);
@@ -219,7 +220,7 @@ void L1Algo::FindMoreHits(L1Branch& t, L1TrackPar& Tout, const bool dir,
const L1Station& sta = fParameters.GetStation(ista);
- fit.Extrapolate(sta.fZ, fld, fvec::One());
+ fit.Extrapolate(sta.fZ, fld);
fscal r2_best = 1e8; // best distance to hit
int iHit_best = -1; // index of the best hit
@@ -293,13 +294,13 @@ void L1Algo::FindMoreHits(L1Branch& t, L1TrackPar& Tout, const bool dir,
auto [x, y] = sta.ConvUVtoXY<fvec>(hit.u, hit.v);
- fit.Extrapolate(hit.z, fld, fvec::One());
- 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);
+ fit.Extrapolate(hit.z, fld);
+ fit.Filter(sta.frontInfo, hit.u, hit.du2);
+ fit.Filter(sta.backInfo, hit.v, hit.dv2);
+ fit.FilterTime(hit.t, hit.dt2, sta.timeInfo);
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());
+ fit.AddMsInMaterial(radThick);
+ fit.EnergyLossCorrection(radThick, dir ? fvec(1.f) : fvec(-1.f));
fldB0 = fldB1;
fldB1 = fldB2;
diff --git a/reco/L1/L1Algo/L1CATrackFinder.cxx b/reco/L1/L1Algo/L1CATrackFinder.cxx
index 63db51f673..3f941d0a96 100644
--- a/reco/L1/L1Algo/L1CATrackFinder.cxx
+++ b/reco/L1/L1Algo/L1CATrackFinder.cxx
@@ -213,7 +213,7 @@ inline void L1Algo::findSingletsStep1( /// input 1st stage of singlet search
L1Fit fit;
fit.SetParticleMass(GetDefaultParticleMass());
-
+ fit.SetMask(fmask::One());
fit.SetQp0(fMaxInvMom);
if (fpCurrentIteration->GetElectronFlag()) { fit.SetParticleMass(L1Constants::phys::kElectronMass); }
@@ -295,11 +295,11 @@ inline void L1Algo::findSingletsStep1( /// input 1st stage of singlet search
fit.AddTargetToLine(fTargX, fTargY, fTargZ, TargetXYInfo, fld0);
}
- fit.AddMsInMaterial(fParameters.GetMaterialThickness(istal, T.x, T.y), fvec::One());
+ fit.AddMsInMaterial(fParameters.GetMaterialThickness(istal, T.x, T.y));
// extrapolate to the middle hit
- fit.ExtrapolateLine(stam.fZ, fld1, fvec::One());
+ fit.ExtrapolateLine(stam.fZ, fld1);
T_1[i1_V] = T;
} // i1_V
@@ -544,6 +544,7 @@ inline void L1Algo::findTripletsStep0( // input
L1Fit fit;
fit.SetParticleMass(fDefaultMass);
+ fit.SetMask(fmask::One());
fit.SetQp0(fvec(0.));
n3 = 0;
@@ -612,21 +613,21 @@ inline void L1Algo::findTripletsStep0( // input
// add the middle hit
- fit.ExtrapolateLine(zPos_2, fvec::One());
+ fit.ExtrapolateLine(zPos_2);
- fit.Filter(stam.frontInfo, u_front_2, du2_2, fvec::One());
- fit.Filter(stam.backInfo, u_back_2, dv2_2, fvec::One());
- fit.FilterTime(t_2, dt2_2, fvec::One(), stam.timeInfo);
+ fit.Filter(stam.frontInfo, u_front_2, du2_2);
+ fit.Filter(stam.backInfo, u_back_2, dv2_2);
+ fit.FilterTime(t_2, dt2_2, stam.timeInfo);
fit.SetQp0(isMomentumFitted ? fit.Tr().qp : fMaxInvMom);
- fit.AddMsInMaterial(fParameters.GetMaterialThickness(iStaM, T2.x, T2.y), fvec::One());
+ fit.AddMsInMaterial(fParameters.GetMaterialThickness(iStaM, T2.x, T2.y));
fit.SetQp0(fit.Tr().qp);
// extrapolate to the right hit station
- fit.Extrapolate(star.fZ, f2, fvec::One());
+ fit.Extrapolate(star.fZ, f2);
// assert(T2.IsConsistent(true, n2_4));
@@ -685,10 +686,11 @@ inline void L1Algo::findTripletsStep0( // input
L1Fit fit3;
fit3.SetParticleMass(fDefaultMass);
+ fit3.SetMask(fmask::One());
fit3.SetTrack(T2);
L1TrackPar& T_cur = fit3.Tr();
- fit3.ExtrapolateLine(zr, fvec::One());
+ fit3.ExtrapolateLine(zr);
if ((star.timeInfo) && (stam.timeInfo))
if (fabs(T_cur.t[i2_4] - hitr.T()) > sqrt(T_cur.C55[i2_4] + hitr.dT2()) * 5) continue;
@@ -729,7 +731,7 @@ inline void L1Algo::findTripletsStep0( // input
L1Fit::FilterChi2(star.backInfo, x, y, C00, C10, C11, chi2, hitr.V(), hitr.dV2());
- fit3.FilterTime(hitr.T(), hitr.dT2(), fvec::One(), star.timeInfo);
+ fit3.FilterTime(hitr.T(), hitr.dT2(), star.timeInfo);
if (!fpCurrentIteration->GetTrackFromTripletsFlag()) {
if (chi2[i2_4] > fTripletChi2Cut || C00[i2_4] < 0 || C11[i2_4] < 0 || T_cur.C55[i2_4] < 0) {
@@ -809,13 +811,13 @@ inline void L1Algo::findTripletsStep1( // input
{
L1Fit fit;
fit.SetParticleMass(fDefaultMass);
-
+ fit.SetMask(fmask::One());
for (Tindex i3_V = 0; i3_V < n3_V; ++i3_V) {
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); }
+ fit.ExtrapolateLine(z_Pos[i3_V]);
+ fit.Filter(star.frontInfo, u_front_[i3_V], du2_3[i3_V]);
+ fit.Filter(star.backInfo, u_back_[i3_V], dv2_3[i3_V]);
+ if (kMcbm != fTrackingMode) { fit.FilterTime(t_3[i3_V], dt2_3[i3_V], star.timeInfo); }
T_3[i3_V] = fit.Tr();
}
}
@@ -833,6 +835,7 @@ inline void L1Algo::findTripletsStep2(Tindex n3, int istal, int istam, int istar
ca::tools::Debugger::Instance().AddNtuple("triplets", "ev:mc:sta:p:vx:vy:vz:chi2:ndf");
L1Fit fit;
+ fit.SetMask(fmask::One());
if (fpCurrentIteration->GetElectronFlag()) { fit.SetParticleMass(L1Constants::phys::kElectronMass); }
else {
@@ -940,21 +943,21 @@ inline void L1Algo::findTripletsStep2(Tindex n3, int istal, int istam, int istar
//std::tie(T.C00, T.C10, T.C11) = sta[ih0].FormXYCovarianceMatrix(du2[ih0], dv2[ih0]);
- fit.Filter(sta[ih0].frontInfo, u[ih0], du2[ih0], fvec::One());
- fit.Filter(sta[ih0].backInfo, v[ih0], dv2[ih0], fvec::One());
- fit.FilterTime(t[ih0], dt2[ih0], fvec::One(), sta[ih0].timeInfo);
+ fit.Filter(sta[ih0].frontInfo, u[ih0], du2[ih0]);
+ fit.Filter(sta[ih0].backInfo, v[ih0], dv2[ih0]);
+ fit.FilterTime(t[ih0], dt2[ih0], sta[ih0].timeInfo);
// add the target constraint
fit.AddTargetToLine(fTargX, fTargY, fTargZ, TargetXYInfo, fldTarget);
for (int ih = 1; ih < NHits; ++ih) {
- fit.Extrapolate(z[ih], fld, fvec::One());
+ fit.Extrapolate(z[ih], fld);
auto radThick = fParameters.GetMaterialThickness(ista[ih], T.x, T.y);
- fit.AddMsInMaterial(radThick, fvec::One());
- fit.EnergyLossCorrection(radThick, fvec(-1.f), fvec::One());
- fit.Filter(sta[ih].frontInfo, u[ih], du2[ih], fvec::One());
- fit.Filter(sta[ih].backInfo, v[ih], dv2[ih], fvec::One());
- fit.FilterTime(t[ih], dt2[ih], fvec::One(), sta[ih].timeInfo);
+ fit.AddMsInMaterial(radThick);
+ fit.EnergyLossCorrection(radThick, fvec(-1.f));
+ fit.Filter(sta[ih].frontInfo, u[ih], du2[ih]);
+ fit.Filter(sta[ih].backInfo, v[ih], dv2[ih]);
+ fit.FilterTime(t[ih], dt2[ih], sta[ih].timeInfo);
}
}
@@ -977,18 +980,18 @@ inline void L1Algo::findTripletsStep2(Tindex n3, int istal, int istam, int istar
T.C55 = dt2[ih0];
//std::tie(T.C00, T.C10, T.C11) = sta[ih0].FormXYCovarianceMatrix(du2[ih0], dv2[ih0]);
- fit.Filter(sta[ih0].frontInfo, u[ih0], du2[ih0], fvec::One());
- fit.Filter(sta[ih0].backInfo, v[ih0], dv2[ih0], fvec::One());
- fit.FilterTime(t[ih0], dt2[ih0], fvec::One(), sta[ih0].timeInfo);
+ fit.Filter(sta[ih0].frontInfo, u[ih0], du2[ih0]);
+ fit.Filter(sta[ih0].backInfo, v[ih0], dv2[ih0]);
+ fit.FilterTime(t[ih0], dt2[ih0], sta[ih0].timeInfo);
for (int ih = NHits - 2; ih >= 0; --ih) {
- fit.Extrapolate(z[ih], fld, fvec::One());
+ fit.Extrapolate(z[ih], fld);
auto radThick = fParameters.GetMaterialThickness(ista[ih], T.x, T.y);
- fit.AddMsInMaterial(radThick, fvec::One());
- fit.EnergyLossCorrection(radThick, fvec(1.f), fvec::One());
- fit.Filter(sta[ih].frontInfo, u[ih], du2[ih], fvec::One());
- fit.Filter(sta[ih].backInfo, v[ih], dv2[ih], fvec::One());
- fit.FilterTime(t[ih], dt2[ih], fvec::One(), sta[ih].timeInfo);
+ fit.AddMsInMaterial(radThick);
+ fit.EnergyLossCorrection(radThick, fvec(1.f));
+ fit.Filter(sta[ih].frontInfo, u[ih], du2[ih]);
+ fit.Filter(sta[ih].backInfo, v[ih], dv2[ih]);
+ fit.FilterTime(t[ih], dt2[ih], sta[ih].timeInfo);
}
}
} // for iiter
diff --git a/reco/L1/L1Algo/L1CloneMerger.cxx b/reco/L1/L1Algo/L1CloneMerger.cxx
index ee511944f2..80a1f062b3 100644
--- a/reco/L1/L1Algo/L1CloneMerger.cxx
+++ b/reco/L1/L1Algo/L1CloneMerger.cxx
@@ -78,10 +78,12 @@ void L1CloneMerger::Exec(L1Vector<L1Track>& extTracks, L1Vector<L1HitIndex_t>& e
L1Fit fitB;
fitB.SetParticleMass(frAlgo.GetDefaultParticleMass());
+ fitB.SetMask(fmask::One());
fitB.SetQp0(fvec(0.));
L1Fit fitF;
fitF.SetParticleMass(frAlgo.GetDefaultParticleMass());
+ fitF.SetMask(fmask::One());
fitF.SetQp0(fvec(0.));
L1TrackPar& Tb = fitB.Tr();
@@ -188,8 +190,8 @@ void L1CloneMerger::Exec(L1Vector<L1Track>& extTracks, L1Vector<L1HitIndex_t>& e
fvec zMiddle = fvec(0.5) * (Tb.z + Tf.z);
- fitF.Extrapolate(zMiddle, fld, fvec::One());
- fitB.Extrapolate(zMiddle, fld, fvec::One());
+ fitF.Extrapolate(zMiddle, fld);
+ fitB.Extrapolate(zMiddle, fld);
fvec Chi2Tracks(0.);
FilterTracks(&(Tf.x), &(Tf.C00), &(Tb.x), &(Tb.C00), 0, 0, &Chi2Tracks);
diff --git a/reco/L1/L1Algo/L1Field.h b/reco/L1/L1Algo/L1Field.h
index 55efd242b9..086dca9c1c 100644
--- a/reco/L1/L1Algo/L1Field.h
+++ b/reco/L1/L1Algo/L1Field.h
@@ -19,8 +19,10 @@ public:
/// Combines the magnetic field with another field value object using weight
/// \param B other field value to combine with
- /// \param w weight from 0 to 1 // TODO: Do we need any checks here? (S.Zharko)
- void Combine(L1FieldValue& B, fvec w); // TODO: Shouldn't the B parameter be const? (S.Zharko)
+ /// \param w weight from 0 to 1
+ void Combine(const L1FieldValue& B, const fvec& w);
+
+ void Combine(const L1FieldValue& B, const fmask& w);
/// Consistency checker
void CheckConsistency() const;
@@ -43,13 +45,20 @@ public:
}
} _fvecalignment;
-[[gnu::always_inline]] inline void L1FieldValue::Combine(L1FieldValue& B, fvec w)
+[[gnu::always_inline]] inline void L1FieldValue::Combine(const L1FieldValue& B, const fvec& w)
{
x += w * (B.x - x);
y += w * (B.y - y);
z += w * (B.z - z);
}
+[[gnu::always_inline]] inline void L1FieldValue::Combine(const L1FieldValue& B, const fmask& w)
+{
+ x(w) = B.x;
+ y(w) = B.y;
+ z(w) = B.z;
+}
+
/// Class represents a set of magnetic field approximation coefficients
///
// TODO: Crosscheck the default content (S.Zharko)
diff --git a/reco/L1/L1Algo/L1Fit.cxx b/reco/L1/L1Algo/L1Fit.cxx
index cd79f57cb4..0639c696bf 100644
--- a/reco/L1/L1Algo/L1Fit.cxx
+++ b/reco/L1/L1Algo/L1Fit.cxx
@@ -6,7 +6,7 @@
#define cnst const fvec
-void L1Fit::Filter(const L1UMeasurementInfo& info, const fvec& u, const fvec& sigma2, const fvec& w)
+void L1Fit::Filter(const L1UMeasurementInfo& info, cnst& u, cnst& sigma2)
{
fvec zeta, HCH;
fvec F0, F1, F2, F3, F4, F5;
@@ -26,17 +26,17 @@ void L1Fit::Filter(const L1UMeasurementInfo& info, const fvec& u, const fvec& si
F5 = info.cos_phi * fTr.C50 + info.sin_phi * fTr.C51;
const fmask maskDoFilter = (HCH < sigma2 * 16.f);
- //const fvec maskDoFilter = _f32vec4_true;
+ //cnst maskDoFilter = _f32vec4_true;
// correction to HCH is needed for the case when sigma2 is so small
// with respect to HCH that it disappears due to the roundoff error
//
- fvec wi = w / (sigma2 + fvec(1.0000001) * HCH);
- fvec zetawi = w * zeta / (iif(maskDoFilter, sigma2, fvec::Zero()) + HCH);
+ fvec wi = fMaskF / (sigma2 + fvec(1.0000001) * HCH);
+ fvec zetawi = fMaskF * zeta / (iif(maskDoFilter, sigma2, fvec::Zero()) + HCH);
// fTr.chi2 += iif( maskDoFilter, zeta * zetawi, fvec::Zero() );
fTr.chi2 += zeta * zeta * wi;
- fTr.NDF += w;
+ fTr.NDF += fMaskF;
K1 = F1 * wi;
K2 = F2 * wi;
@@ -75,7 +75,7 @@ void L1Fit::Filter(const L1UMeasurementInfo& info, const fvec& u, const fvec& si
}
-void L1Fit::FilterTime(fvec t, fvec dt2, fvec w, fvec timeInfo)
+void L1Fit::FilterTime(cnst& t, cnst& dt2, cnst& timeInfo)
{
// filter track with a time measurement
@@ -89,6 +89,8 @@ void L1Fit::FilterTime(fvec t, fvec dt2, fvec w, fvec timeInfo)
fvec HCH = fTr.C55;
+ fvec w = fMaskF;
+
w.setZero(timeInfo <= fvec::Zero());
// when dt0 is much smaller than current time error,
@@ -142,7 +144,7 @@ void L1Fit::FilterTime(fvec t, fvec dt2, fvec w, fvec timeInfo)
}
-void L1Fit::FilterXY(const L1XYMeasurementInfo& info, fvec x, fvec y)
+void L1Fit::FilterXY(const L1XYMeasurementInfo& info, cnst& x, cnst& y)
{
cnst TWO(2.);
@@ -222,8 +224,8 @@ void L1Fit::FilterXY(const L1XYMeasurementInfo& info, fvec x, fvec y)
fTr.C55 -= K50 * F50 + K51 * F51;
}
-void L1Fit::FilterExtrapolatedXY(const fvec& x, const fvec& y, const L1XYMeasurementInfo& info, const fvec& extrX,
- const fvec& extrY, const fvec Jx[6], const fvec Jy[6])
+void L1Fit::FilterExtrapolatedXY(cnst& x, cnst& y, const L1XYMeasurementInfo& info, cnst& extrX, cnst& extrY,
+ cnst Jx[6], cnst Jy[6])
{
// add a 2-D measurenent (x,y) at some z, that differs from fTr.z
// extrX, extrY are extrapolated track parameters at z, Jx, Jy are derivatives of the extrapolation
@@ -304,22 +306,22 @@ void L1Fit::FilterExtrapolatedXY(const fvec& x, const fvec& y, const L1XYMeasure
void L1Fit::Extrapolate // extrapolates track parameters and returns jacobian for extrapolation of CovMatrix
- (fvec z_out, // extrapolate to this z position
- const L1FieldRegion& F, const fvec& w)
+ (cnst& z_out, // extrapolate to this z position
+ const L1FieldRegion& F)
{
// use Q/p linearisation at fQp0
fvec sgn = iif(fTr.z < z_out, fvec(1.), fvec(-1.));
- while (!(w * abs(z_out - fTr.z) <= fvec(1.e-6)).isFull()) {
+ while (!(fMaskF * abs(z_out - fTr.z) <= fvec(1.e-6)).isFull()) {
fvec zNew = fTr.z + sgn * fvec(50.); // max. 50 cm step
zNew(sgn * (z_out - zNew) <= fvec(0.)) = z_out;
- ExtrapolateStep(zNew, F, w);
+ ExtrapolateStep(zNew, F);
}
}
void L1Fit::ExtrapolateStep // extrapolates track parameters and returns jacobian for extrapolation of CovMatrix
- (fvec z_out, // extrapolate to this z position
- const L1FieldRegion& F, const fvec& w)
+ (cnst& zOut, // extrapolate to this z position
+ const L1FieldRegion& F)
{
// use Q/p linearisation at fQp0
@@ -357,13 +359,13 @@ void L1Fit::ExtrapolateStep // extrapolates track parameters and returns jacobi
//----------------------------------------------------------------
- z_out = iif((fvec(0.f) < w), z_out, fTr.z);
+ cnst zMasked = iif(fMask, zOut, fTr.z);
- fvec qp_in = fTr.qp;
- const fvec h = (z_out - fTr.z);
+ fvec qp_in = fTr.qp;
+ cnst h = (zMasked - fTr.z);
- const fvec a[4] = {0., h * fvec(0.5), h * fvec(0.5), h};
- const fvec c[4] = {h / fvec(6.), h / fvec(3.), h / fvec(3.), h / fvec(6.)};
+ cnst a[4] = {0., h * fvec(0.5), h * fvec(0.5), h};
+ cnst c[4] = {h / fvec(6.), h / fvec(3.), h / fvec(3.), h / fvec(6.)};
fvec f0[4]; // ( dx/dz )
fvec f1[4]; // ( dy/dz )
@@ -393,9 +395,9 @@ void L1Fit::ExtrapolateStep // extrapolates track parameters and returns jacobi
}
fvec B[3];
- const fvec& Bx = B[0];
- const fvec& By = B[1];
- const fvec& Bz = B[2];
+ cnst& Bx = B[0];
+ cnst& By = B[1];
+ cnst& Bz = B[2];
F.Get(r[0], r[1], z, B);
@@ -444,7 +446,7 @@ void L1Fit::ExtrapolateStep // extrapolates track parameters and returns jacobi
fTr.tx = r0[2] + (c[0] * k[1][2] + c[1] * k[2][2] + c[2] * k[3][2] + c[3] * k[4][2]);
fTr.ty = r0[3] + (c[0] * k[1][3] + c[1] * k[2][3] + c[2] * k[3][3] + c[3] * k[4][3]);
fTr.t = r0[4] + (c[0] * k[1][4] + c[1] * k[2][4] + c[2] * k[3][4] + c[3] * k[4][4]);
- fTr.z = z_out;
+ fTr.z = zMasked;
}
@@ -658,8 +660,8 @@ void L1Fit::ExtrapolateStep // extrapolates track parameters and returns jacobi
void
L1Fit::ExtrapolateStepAnalytic // extrapolates track parameters and returns jacobian for extrapolation of CovMatrix
- (fvec z_out, // extrapolate to this z position
- const L1FieldRegion& F, const fvec& w)
+ (cnst& z_out, // extrapolate to this z position
+ const L1FieldRegion& F)
{
//
// Part of the analytic extrapolation formula with error (c_light*B*dz)^4/4!
@@ -673,72 +675,72 @@ void
cnst c1(1.), c2(2.), c3(3.), c4(4.), c6(6.), c9(9.), c15(15.), c18(18.), c45(45.), c2i(1. / 2.), c3i(1. / 3.),
c6i(1. / 6.), c12i(1. / 12.);
- const fvec qp = fTr.qp;
- fvec dz = (z_out - fTr.z);
+ cnst qp = fTr.qp;
+ fvec dz = (z_out - fTr.z);
- dz.setZero(w <= fvec::Zero());
+ dz.setZero(!fMask);
- const fvec dz2 = dz * dz;
- const fvec dz3 = dz2 * dz;
+ cnst dz2 = dz * dz;
+ cnst dz3 = dz2 * dz;
// construct coefficients
- const fvec x = fTr.tx;
- const fvec y = fTr.ty;
- const fvec xx = x * x;
- const fvec xy = x * y;
- const fvec yy = y * y;
- const fvec y2 = y * c2;
- const fvec x2 = x * c2;
- const fvec x4 = x * c4;
- const fvec xx31 = xx * c3 + c1;
- const fvec xx159 = xx * c15 + c9;
-
- const fvec Ay = -xx - c1;
- const fvec Ayy = x * (xx * c3 + c3);
- const fvec Ayz = -c2 * xy;
- const fvec Ayyy = -(c15 * xx * xx + c18 * xx + c3);
-
- const fvec Ayy_x = c3 * xx31;
- const fvec Ayyy_x = -x4 * xx159;
-
- const fvec Bx = yy + c1;
- const fvec Byy = y * xx31;
- const fvec Byz = c2 * xx + c1;
- const fvec Byyy = -xy * xx159;
-
- const fvec Byy_x = c6 * xy;
- const fvec Byyy_x = -y * (c45 * xx + c9);
- const fvec Byyy_y = -x * xx159;
+ cnst x = fTr.tx;
+ cnst y = fTr.ty;
+ cnst xx = x * x;
+ cnst xy = x * y;
+ cnst yy = y * y;
+ cnst y2 = y * c2;
+ cnst x2 = x * c2;
+ cnst x4 = x * c4;
+ cnst xx31 = xx * c3 + c1;
+ cnst xx159 = xx * c15 + c9;
+
+ cnst Ay = -xx - c1;
+ cnst Ayy = x * (xx * c3 + c3);
+ cnst Ayz = -c2 * xy;
+ cnst Ayyy = -(c15 * xx * xx + c18 * xx + c3);
+
+ cnst Ayy_x = c3 * xx31;
+ cnst Ayyy_x = -x4 * xx159;
+
+ cnst Bx = yy + c1;
+ cnst Byy = y * xx31;
+ cnst Byz = c2 * xx + c1;
+ cnst Byyy = -xy * xx159;
+
+ cnst Byy_x = c6 * xy;
+ cnst Byyy_x = -y * (c45 * xx + c9);
+ cnst Byyy_y = -x * xx159;
// end of coefficients calculation
- const fvec t2 = c1 + xx + yy;
- const fvec t = sqrt(t2);
- const fvec h = fQp0 * c_light;
- const fvec ht = h * t;
+ cnst t2 = c1 + xx + yy;
+ cnst t = sqrt(t2);
+ cnst h = fQp0 * c_light;
+ cnst ht = h * t;
// get field integrals
- const fvec ddz = fTr.z - F.z0;
- const fvec Fx0 = F.cx0 + F.cx1 * ddz + F.cx2 * ddz * ddz;
- const fvec Fx1 = (F.cx1 + c2 * F.cx2 * ddz) * dz;
- const fvec Fx2 = F.cx2 * dz2;
- const fvec Fy0 = F.cy0 + F.cy1 * ddz + F.cy2 * ddz * ddz;
- const fvec Fy1 = (F.cy1 + c2 * F.cy2 * ddz) * dz;
- const fvec Fy2 = F.cy2 * dz2;
- const fvec Fz0 = F.cz0 + F.cz1 * ddz + F.cz2 * ddz * ddz;
- const fvec Fz1 = (F.cz1 + c2 * F.cz2 * ddz) * dz;
- const fvec Fz2 = F.cz2 * dz2;
+ cnst ddz = fTr.z - F.z0;
+ cnst Fx0 = F.cx0 + F.cx1 * ddz + F.cx2 * ddz * ddz;
+ cnst Fx1 = (F.cx1 + c2 * F.cx2 * ddz) * dz;
+ cnst Fx2 = F.cx2 * dz2;
+ cnst Fy0 = F.cy0 + F.cy1 * ddz + F.cy2 * ddz * ddz;
+ cnst Fy1 = (F.cy1 + c2 * F.cy2 * ddz) * dz;
+ cnst Fy2 = F.cy2 * dz2;
+ cnst Fz0 = F.cz0 + F.cz1 * ddz + F.cz2 * ddz * ddz;
+ cnst Fz1 = (F.cz1 + c2 * F.cz2 * ddz) * dz;
+ cnst Fz2 = F.cz2 * dz2;
//
- const fvec sx = (Fx0 + Fx1 * c2i + Fx2 * c3i);
- const fvec sy = (Fy0 + Fy1 * c2i + Fy2 * c3i);
- const fvec sz = (Fz0 + Fz1 * c2i + Fz2 * c3i);
+ cnst sx = (Fx0 + Fx1 * c2i + Fx2 * c3i);
+ cnst sy = (Fy0 + Fy1 * c2i + Fy2 * c3i);
+ cnst sz = (Fz0 + Fz1 * c2i + Fz2 * c3i);
- const fvec Sx = (Fx0 * c2i + Fx1 * c6i + Fx2 * c12i);
- const fvec Sy = (Fy0 * c2i + Fy1 * c6i + Fy2 * c12i);
- const fvec Sz = (Fz0 * c2i + Fz1 * c6i + Fz2 * c12i);
+ cnst Sx = (Fx0 * c2i + Fx1 * c6i + Fx2 * c12i);
+ cnst Sy = (Fy0 * c2i + Fy1 * c6i + Fy2 * c12i);
+ cnst Sz = (Fz0 * c2i + Fz1 * c6i + Fz2 * c12i);
fvec syz;
{
@@ -758,8 +760,8 @@ void
+ Fy2 * (c20 * Fz0 + c21 * Fz1 + c22 * Fz2);
}
- const fvec syy = sy * sy * c2i;
- const fvec syyy = syy * sy * c3i;
+ cnst syy = sy * sy * c2i;
+ cnst syyy = syy * sy * c3i;
fvec Syy;
{
@@ -773,71 +775,71 @@ void
constexpr double d = 181440.;
cnst c000(7560 / d), c001(9 * 1008 / d), c002(5 * 1008 / d), c011(21 * 180 / d), c012(24 * 180 / d),
c022(7 * 180 / d), c111(540 / d), c112(945 / d), c122(560 / d), c222(112 / d);
- const fvec Fy22 = Fy2 * Fy2;
- Syyy = Fy0 * (Fy0 * (c000 * Fy0 + c001 * Fy1 + c002 * Fy2) + Fy1 * (c011 * Fy1 + c012 * Fy2) + c022 * Fy22)
+ cnst Fy22 = Fy2 * Fy2;
+ Syyy = Fy0 * (Fy0 * (c000 * Fy0 + c001 * Fy1 + c002 * Fy2) + Fy1 * (c011 * Fy1 + c012 * Fy2) + c022 * Fy22)
+ Fy1 * (Fy1 * (c111 * Fy1 + c112 * Fy2) + c122 * Fy22) + c222 * Fy22 * Fy2;
}
- const fvec sA1 = sx * xy + sy * Ay + sz * y;
- const fvec sA1_x = sx * y - sy * x2;
- const fvec sA1_y = sx * x + sz;
-
- const fvec sB1 = sx * Bx - sy * xy - sz * x;
- const fvec sB1_x = -sy * y - sz;
- const fvec sB1_y = sx * y2 - sy * x;
-
- const fvec SA1 = Sx * xy + Sy * Ay + Sz * y;
- const fvec SA1_x = Sx * y - Sy * x2;
- const fvec SA1_y = Sx * x + Sz;
-
- const fvec SB1 = Sx * Bx - Sy * xy - Sz * x;
- const fvec SB1_x = -Sy * y - Sz;
- const fvec SB1_y = Sx * y2 - Sy * x;
-
-
- const fvec sA2 = syy * Ayy + syz * Ayz;
- const fvec sA2_x = syy * Ayy_x - syz * y2;
- const fvec sA2_y = -syz * x2;
- const fvec sB2 = syy * Byy + syz * Byz;
- const fvec sB2_x = syy * Byy_x + syz * x4;
- const fvec sB2_y = syy * xx31;
-
- const fvec SA2 = Syy * Ayy + Syz * Ayz;
- const fvec SA2_x = Syy * Ayy_x - Syz * y2;
- const fvec SA2_y = -Syz * x2;
- const fvec SB2 = Syy * Byy + Syz * Byz;
- const fvec SB2_x = Syy * Byy_x + Syz * x4;
- const fvec SB2_y = Syy * xx31;
-
- const fvec sA3 = syyy * Ayyy;
- const fvec sA3_x = syyy * Ayyy_x;
- const fvec sB3 = syyy * Byyy;
- const fvec sB3_x = syyy * Byyy_x;
- const fvec sB3_y = syyy * Byyy_y;
-
-
- const fvec SA3 = Syyy * Ayyy;
- const fvec SA3_x = Syyy * Ayyy_x;
- const fvec SB3 = Syyy * Byyy;
- const fvec SB3_x = Syyy * Byyy_x;
- const fvec SB3_y = Syyy * Byyy_y;
-
- const fvec ht1 = ht * dz;
- const fvec ht2 = ht * ht * dz2;
- const fvec ht3 = ht * ht * ht * dz3;
- const fvec ht1sA1 = ht1 * sA1;
- const fvec ht1sB1 = ht1 * sB1;
- const fvec ht1SA1 = ht1 * SA1;
- const fvec ht1SB1 = ht1 * SB1;
- const fvec ht2sA2 = ht2 * sA2;
- const fvec ht2SA2 = ht2 * SA2;
- const fvec ht2sB2 = ht2 * sB2;
- const fvec ht2SB2 = ht2 * SB2;
- const fvec ht3sA3 = ht3 * sA3;
- const fvec ht3sB3 = ht3 * sB3;
- const fvec ht3SA3 = ht3 * SA3;
- const fvec ht3SB3 = ht3 * SB3;
+ cnst sA1 = sx * xy + sy * Ay + sz * y;
+ cnst sA1_x = sx * y - sy * x2;
+ cnst sA1_y = sx * x + sz;
+
+ cnst sB1 = sx * Bx - sy * xy - sz * x;
+ cnst sB1_x = -sy * y - sz;
+ cnst sB1_y = sx * y2 - sy * x;
+
+ cnst SA1 = Sx * xy + Sy * Ay + Sz * y;
+ cnst SA1_x = Sx * y - Sy * x2;
+ cnst SA1_y = Sx * x + Sz;
+
+ cnst SB1 = Sx * Bx - Sy * xy - Sz * x;
+ cnst SB1_x = -Sy * y - Sz;
+ cnst SB1_y = Sx * y2 - Sy * x;
+
+
+ cnst sA2 = syy * Ayy + syz * Ayz;
+ cnst sA2_x = syy * Ayy_x - syz * y2;
+ cnst sA2_y = -syz * x2;
+ cnst sB2 = syy * Byy + syz * Byz;
+ cnst sB2_x = syy * Byy_x + syz * x4;
+ cnst sB2_y = syy * xx31;
+
+ cnst SA2 = Syy * Ayy + Syz * Ayz;
+ cnst SA2_x = Syy * Ayy_x - Syz * y2;
+ cnst SA2_y = -Syz * x2;
+ cnst SB2 = Syy * Byy + Syz * Byz;
+ cnst SB2_x = Syy * Byy_x + Syz * x4;
+ cnst SB2_y = Syy * xx31;
+
+ cnst sA3 = syyy * Ayyy;
+ cnst sA3_x = syyy * Ayyy_x;
+ cnst sB3 = syyy * Byyy;
+ cnst sB3_x = syyy * Byyy_x;
+ cnst sB3_y = syyy * Byyy_y;
+
+
+ cnst SA3 = Syyy * Ayyy;
+ cnst SA3_x = Syyy * Ayyy_x;
+ cnst SB3 = Syyy * Byyy;
+ cnst SB3_x = Syyy * Byyy_x;
+ cnst SB3_y = Syyy * Byyy_y;
+
+ cnst ht1 = ht * dz;
+ cnst ht2 = ht * ht * dz2;
+ cnst ht3 = ht * ht * ht * dz3;
+ cnst ht1sA1 = ht1 * sA1;
+ cnst ht1sB1 = ht1 * sB1;
+ cnst ht1SA1 = ht1 * SA1;
+ cnst ht1SB1 = ht1 * SB1;
+ cnst ht2sA2 = ht2 * sA2;
+ cnst ht2SA2 = ht2 * SA2;
+ cnst ht2sB2 = ht2 * sB2;
+ cnst ht2SB2 = ht2 * SB2;
+ cnst ht3sA3 = ht3 * sA3;
+ cnst ht3sB3 = ht3 * sB3;
+ cnst ht3SA3 = ht3 * SA3;
+ cnst ht3SB3 = ht3 * SB3;
fTr.x += (x + ht1SA1 + ht2SA2 + ht3SA3) * dz;
fTr.y += (y + ht1SB1 + ht2SB2 + ht3SB3) * dz;
@@ -845,32 +847,32 @@ void
fTr.ty += ht1sB1 + ht2sB2 + ht3sB3;
fTr.z += dz;
- const fvec ctdz = c_light * t * dz;
- const fvec ctdz2 = c_light * t * dz2;
+ cnst ctdz = c_light * t * dz;
+ cnst ctdz2 = c_light * t * dz2;
- const fvec dqp = qp - fQp0;
- const fvec t2i = c1 / t2;
- const fvec xt2i = x * t2i;
- const fvec yt2i = y * t2i;
- const fvec tmp0 = ht1SA1 + c2 * ht2SA2 + c3 * ht3SA3;
- const fvec tmp1 = ht1SB1 + c2 * ht2SB2 + c3 * ht3SB3;
- const fvec tmp2 = ht1sA1 + c2 * ht2sA2 + c3 * ht3sA3;
- const fvec tmp3 = ht1sB1 + c2 * ht2sB2 + c3 * ht3sB3;
+ cnst dqp = qp - fQp0;
+ cnst t2i = c1 / t2;
+ cnst xt2i = x * t2i;
+ cnst yt2i = y * t2i;
+ cnst tmp0 = ht1SA1 + c2 * ht2SA2 + c3 * ht3SA3;
+ cnst tmp1 = ht1SB1 + c2 * ht2SB2 + c3 * ht3SB3;
+ cnst tmp2 = ht1sA1 + c2 * ht2sA2 + c3 * ht3sA3;
+ cnst tmp3 = ht1sB1 + c2 * ht2sB2 + c3 * ht3sB3;
- const fvec j02 = dz * (c1 + xt2i * tmp0 + ht1 * SA1_x + ht2 * SA2_x + ht3 * SA3_x);
- const fvec j12 = dz * (xt2i * tmp1 + ht1 * SB1_x + ht2 * SB2_x + ht3 * SB3_x);
- const fvec j22 = c1 + xt2i * tmp2 + ht1 * sA1_x + ht2 * sA2_x + ht3 * sA3_x;
- const fvec j32 = xt2i * tmp3 + ht1 * sB1_x + ht2 * sB2_x + ht3 * sB3_x;
+ cnst j02 = dz * (c1 + xt2i * tmp0 + ht1 * SA1_x + ht2 * SA2_x + ht3 * SA3_x);
+ cnst j12 = dz * (xt2i * tmp1 + ht1 * SB1_x + ht2 * SB2_x + ht3 * SB3_x);
+ cnst j22 = c1 + xt2i * tmp2 + ht1 * sA1_x + ht2 * sA2_x + ht3 * sA3_x;
+ cnst j32 = xt2i * tmp3 + ht1 * sB1_x + ht2 * sB2_x + ht3 * sB3_x;
- const fvec j03 = dz * (yt2i * tmp0 + ht1 * SA1_y + ht2 * SA2_y);
- const fvec j13 = dz * (c1 + yt2i * tmp1 + ht1 * SB1_y + ht2 * SB2_y + ht3 * SB3_y);
- const fvec j23 = yt2i * tmp2 + ht1 * sA1_y + ht2 * sA2_y;
- const fvec j33 = c1 + yt2i * tmp3 + ht1 * sB1_y + ht2 * sB2_y + ht3 * sB3_y;
+ cnst j03 = dz * (yt2i * tmp0 + ht1 * SA1_y + ht2 * SA2_y);
+ cnst j13 = dz * (c1 + yt2i * tmp1 + ht1 * SB1_y + ht2 * SB2_y + ht3 * SB3_y);
+ cnst j23 = yt2i * tmp2 + ht1 * sA1_y + ht2 * sA2_y;
+ cnst j33 = c1 + yt2i * tmp3 + ht1 * sB1_y + ht2 * sB2_y + ht3 * sB3_y;
- const fvec j04 = ctdz2 * (SA1 + c2 * ht1 * SA2 + c3 * ht2 * SA3);
- const fvec j14 = ctdz2 * (SB1 + c2 * ht1 * SB2 + c3 * ht2 * SB3);
- const fvec j24 = ctdz * (sA1 + c2 * ht1 * sA2 + c3 * ht2 * sA3);
- const fvec j34 = ctdz * (sB1 + c2 * ht1 * sB2 + c3 * ht2 * sB3);
+ cnst j04 = ctdz2 * (SA1 + c2 * ht1 * SA2 + c3 * ht2 * SA3);
+ cnst j14 = ctdz2 * (SB1 + c2 * ht1 * SB2 + c3 * ht2 * SB3);
+ cnst j24 = ctdz * (sA1 + c2 * ht1 * sA2 + c3 * ht2 * sA3);
+ cnst j34 = ctdz * (sB1 + c2 * ht1 * sB2 + c3 * ht2 * sB3);
// extrapolate inverse momentum
@@ -882,27 +884,27 @@ void
// covariance matrix transport
- const fvec c42 = fTr.C42, c43 = fTr.C43;
+ cnst c42 = fTr.C42, c43 = fTr.C43;
- const fvec cj00 = fTr.C00 + fTr.C20 * j02 + fTr.C30 * j03 + fTr.C40 * j04;
- // const fvec cj10 = fTr.C10 + fTr.C21*j02 + fTr.C31*j03 + fTr.C41*j04;
- const fvec cj20 = fTr.C20 + fTr.C22 * j02 + fTr.C32 * j03 + c42 * j04;
- const fvec cj30 = fTr.C30 + fTr.C32 * j02 + fTr.C33 * j03 + c43 * j04;
+ cnst cj00 = fTr.C00 + fTr.C20 * j02 + fTr.C30 * j03 + fTr.C40 * j04;
+ // cnst cj10 = fTr.C10 + fTr.C21*j02 + fTr.C31*j03 + fTr.C41*j04;
+ cnst cj20 = fTr.C20 + fTr.C22 * j02 + fTr.C32 * j03 + c42 * j04;
+ cnst cj30 = fTr.C30 + fTr.C32 * j02 + fTr.C33 * j03 + c43 * j04;
- const fvec cj01 = fTr.C10 + fTr.C20 * j12 + fTr.C30 * j13 + fTr.C40 * j14;
- const fvec cj11 = fTr.C11 + fTr.C21 * j12 + fTr.C31 * j13 + fTr.C41 * j14;
- const fvec cj21 = fTr.C21 + fTr.C22 * j12 + fTr.C32 * j13 + c42 * j14;
- const fvec cj31 = fTr.C31 + fTr.C32 * j12 + fTr.C33 * j13 + c43 * j14;
+ cnst cj01 = fTr.C10 + fTr.C20 * j12 + fTr.C30 * j13 + fTr.C40 * j14;
+ cnst cj11 = fTr.C11 + fTr.C21 * j12 + fTr.C31 * j13 + fTr.C41 * j14;
+ cnst cj21 = fTr.C21 + fTr.C22 * j12 + fTr.C32 * j13 + c42 * j14;
+ cnst cj31 = fTr.C31 + fTr.C32 * j12 + fTr.C33 * j13 + c43 * j14;
- // const fvec cj02 = fTr.C20*j22 + fTr.C30*j23 + fTr.C40*j24;
- // const fvec cj12 = fTr.C21*j22 + fTr.C31*j23 + fTr.C41*j24;
- const fvec cj22 = fTr.C22 * j22 + fTr.C32 * j23 + c42 * j24;
- const fvec cj32 = fTr.C32 * j22 + fTr.C33 * j23 + c43 * j24;
+ // cnst cj02 = fTr.C20*j22 + fTr.C30*j23 + fTr.C40*j24;
+ // cnst cj12 = fTr.C21*j22 + fTr.C31*j23 + fTr.C41*j24;
+ cnst cj22 = fTr.C22 * j22 + fTr.C32 * j23 + c42 * j24;
+ cnst cj32 = fTr.C32 * j22 + fTr.C33 * j23 + c43 * j24;
- // const fvec cj03 = fTr.C20*j32 + fTr.C30*j33 + fTr.C40*j34;
- // const fvec cj13 = fTr.C21*j32 + fTr.C31*j33 + fTr.C41*j34;
- const fvec cj23 = fTr.C22 * j32 + fTr.C32 * j33 + c42 * j34;
- const fvec cj33 = fTr.C32 * j32 + fTr.C33 * j33 + c43 * j34;
+ // cnst cj03 = fTr.C20*j32 + fTr.C30*j33 + fTr.C40*j34;
+ // cnst cj13 = fTr.C21*j32 + fTr.C31*j33 + fTr.C41*j34;
+ cnst cj23 = fTr.C22 * j32 + fTr.C32 * j33 + c42 * j34;
+ cnst cj33 = fTr.C32 * j32 + fTr.C33 * j33 + c43 * j34;
fTr.C40 += c42 * j02 + c43 * j03 + fTr.C44 * j04; // cj40
fTr.C41 += c42 * j12 + c43 * j13 + fTr.C44 * j14; // cj41
@@ -923,18 +925,18 @@ void
//cout<<"Extrapolation ok"<<endl;
}
-void L1Fit::ExtrapolateLine(fvec z_out, const L1FieldRegion& F, const fvec& w)
+void L1Fit::ExtrapolateLine(cnst& z_out, const L1FieldRegion& F)
{
// extrapolate the track assuming fQp0 == 0
// TODO: write special simplified procedure
//
auto qp0 = fQp0;
fQp0 = fvec(0.);
- Extrapolate(z_out, F, w);
+ Extrapolate(z_out, F);
fQp0 = qp0;
}
-void L1Fit::ExtrapolateLine(fvec z_out, const fvec& w)
+void L1Fit::ExtrapolateLine(cnst& z_out)
{
// extrapolate the track assuming fQp0 == 0
//
@@ -944,7 +946,7 @@ void L1Fit::ExtrapolateLine(fvec z_out, const fvec& w)
cnst c_light(29.9792458);
fvec dz = (z_out - fTr.z);
- dz.setZero(w <= fvec::Zero());
+ dz.setZero(!fMask);
L1TrackPar& T = fTr;
@@ -952,8 +954,8 @@ void L1Fit::ExtrapolateLine(fvec z_out, const fvec& w)
fvec L = sqrt(T.tx * T.tx + T.ty * T.ty + fvec(1.));
T.t += dz * L / c_light;
- const fvec k1 = dz * T.tx / (c_light * L);
- const fvec k2 = dz * T.ty / (c_light * L);
+ cnst k1 = dz * T.tx / (c_light * L);
+ cnst k2 = dz * T.ty / (c_light * L);
fvec c52 = T.C52;
fvec c53 = T.C53;
@@ -971,17 +973,17 @@ void L1Fit::ExtrapolateLine(fvec z_out, const fvec& w)
T.y += T.ty * dz;
T.z += dz;
- const fvec dzC32_in = dz * T.C32;
+ cnst dzC32_in = dz * T.C32;
T.C21 += dzC32_in;
T.C10 += dz * (T.C21 + T.C30);
- const fvec C20_in = T.C20;
+ cnst C20_in = T.C20;
T.C20 += dz * T.C22;
T.C00 += dz * (T.C20 + C20_in);
- const fvec C31_in = T.C31;
+ cnst C31_in = T.C31;
T.C31 += dz * T.C33;
T.C11 += dz * (T.C31 + C31_in);
@@ -992,7 +994,7 @@ void L1Fit::ExtrapolateLine(fvec z_out, const fvec& w)
}
-void L1Fit::AddMsInMaterial(const fvec& radThick, fvec w)
+void L1Fit::AddMsInMaterial(cnst& radThick)
{
cnst ONE = 1.;
@@ -1012,12 +1014,12 @@ void L1Fit::AddMsInMaterial(const fvec& radThick, fvec w)
//fvec a = ( (ONE+mass2*qp0*qp0t)*radThick*s0*s0 );
fvec a = ((t + fMass2 * fQp0 * qp0t) * radThick * s0 * s0);
- fTr.C22 += w * txtx1 * a;
- fTr.C32 += w * tx * ty * a;
- fTr.C33 += w * (ONE + tyty) * a;
+ fTr.C22(fMask) += txtx1 * a;
+ fTr.C32(fMask) += tx * ty * a;
+ fTr.C33(fMask) += (ONE + tyty) * a;
}
-void L1Fit::AddMsInThickMaterial(fvec radThick, fvec w, fvec thickness, bool fDownstream)
+void L1Fit::AddMsInThickMaterial(cnst& radThick, cnst& thickness, bool fDownstream)
{
cnst ONE = 1.;
@@ -1042,39 +1044,39 @@ void L1Fit::AddMsInThickMaterial(fvec radThick, fvec w, fvec thickness, bool fDo
fvec T23 = (thickness * thickness) / fvec(3.0);
fvec T2 = thickness / fvec(2.0);
- fTr.C00 += w * txtx1 * a * T23;
- fTr.C10 += w * tx * ty * a * T23;
- fTr.C20 += w * txtx1 * a * D * T2;
- fTr.C30 += w * tx * ty * a * D * T2;
+ fTr.C00(fMask) += txtx1 * a * T23;
+ fTr.C10(fMask) += tx * ty * a * T23;
+ fTr.C20(fMask) += txtx1 * a * D * T2;
+ fTr.C30(fMask) += tx * ty * a * D * T2;
- fTr.C11 += w * (ONE + tyty) * a * T23;
- fTr.C21 += w * tx * ty * a * D * T2;
- fTr.C31 += w * (ONE + tyty) * a * D * T2;
+ fTr.C11(fMask) += (ONE + tyty) * a * T23;
+ fTr.C21(fMask) += tx * ty * a * D * T2;
+ fTr.C31(fMask) += (ONE + tyty) * a * D * T2;
- fTr.C22 += w * txtx1 * a;
- fTr.C32 += w * tx * ty * a;
- fTr.C33 += w * (ONE + tyty) * a;
+ fTr.C22(fMask) += txtx1 * a;
+ fTr.C32(fMask) += tx * ty * a;
+ fTr.C33(fMask) += (ONE + tyty) * a;
}
-void L1Fit::EnergyLossCorrection(const fvec& radThick, fvec upstreamDirection, fvec w)
+void L1Fit::EnergyLossCorrection(cnst& radThick, cnst& upstreamDirection)
{
- const fvec qp2cut(1. / (10. * 10.)); // 10 GeV cut
- const fvec qp02 = max(fQp0 * fQp0, qp2cut);
- const fvec p2 = fvec(1.) / qp02;
- const fvec E2 = fMass2 + p2;
+ cnst qp2cut(1. / (10. * 10.)); // 10 GeV cut
+ cnst qp02 = max(fQp0 * fQp0, qp2cut);
+ cnst p2 = fvec(1.) / qp02;
+ cnst E2 = fMass2 + p2;
- const fvec bethe = ApproximateBetheBloch(p2 / fMass2);
+ cnst bethe = ApproximateBetheBloch(p2 / fMass2);
fvec tr = sqrt(fvec(1.f) + fTr.tx * fTr.tx + fTr.ty * fTr.ty);
- const fvec dE = bethe * radThick * tr * 2.33f * 9.34961f;
+ cnst dE = bethe * radThick * tr * 2.33f * 9.34961f;
- const fvec ECorrected = sqrt(E2) + upstreamDirection * dE;
- const fvec E2Corrected = ECorrected * ECorrected;
+ cnst ECorrected = sqrt(E2) + upstreamDirection * dE;
+ cnst E2Corrected = ECorrected * ECorrected;
fvec corr = sqrt(p2 / (E2Corrected - fMass2));
- fmask ok = (corr == corr) & (fvec::Zero() < w);
+ fmask ok = (corr == corr) && fMask;
corr = iif(ok, corr, fvec::One());
fQp0 *= corr;
@@ -1088,28 +1090,28 @@ void L1Fit::EnergyLossCorrection(const fvec& radThick, fvec upstreamDirection, f
}
template<int atomicZ>
-void L1Fit::EnergyLossCorrection(float atomicA, float rho, float radLen, const fvec& radThick, fvec direction, fvec w)
+void L1Fit::EnergyLossCorrection(float atomicA, float rho, float radLen, cnst& radThick, cnst& direction)
{
- const fvec qp2cut(1. / (10. * 10.)); // 10 GeV cut
- const fvec qp02 = max(fQp0 * fQp0, qp2cut);
- const fvec p2 = fvec(1.) / qp02;
- const fvec E2 = fMass2 + p2;
+ cnst qp2cut(1. / (10. * 10.)); // 10 GeV cut
+ cnst qp02 = max(fQp0 * fQp0, qp2cut);
+ cnst p2 = fvec(1.) / qp02;
+ cnst E2 = fMass2 + p2;
fvec i;
if (atomicZ < 13) i = (12. * atomicZ + 7.) * 1.e-9;
else
i = (9.76 * atomicZ + 58.8 * std::pow(atomicZ, -0.19)) * 1.e-9;
- const fvec bethe = ApproximateBetheBloch(p2 / fMass2, rho, 0.20, 3.00, i, atomicZ / atomicA);
+ cnst bethe = ApproximateBetheBloch(p2 / fMass2, rho, 0.20, 3.00, i, atomicZ / atomicA);
fvec tr = sqrt(fvec(1.f) + fTr.tx * fTr.tx + fTr.ty * fTr.ty);
fvec dE = bethe * radThick * tr * radLen * rho;
- const fvec E2Corrected = (sqrt(E2) + direction * dE) * (sqrt(E2) + direction * dE);
- fvec corr = sqrt(p2 / (E2Corrected - fMass2));
- fmask ok = (corr == corr) & (fvec::Zero() < w);
- corr = iif(ok, corr, fvec::One());
+ cnst E2Corrected = (sqrt(E2) + direction * dE) * (sqrt(E2) + direction * dE);
+ fvec corr = sqrt(p2 / (E2Corrected - fMass2));
+ fmask ok = (corr == corr) && fMask;
+ corr = iif(ok, corr, fvec::One());
fQp0 *= corr;
fTr.qp *= corr;
@@ -1151,35 +1153,35 @@ void L1Fit::EnergyLossCorrection(float atomicA, float rho, float radLen, const f
}
-void L1Fit::EnergyLossCorrectionIron(const fvec& radThick, fvec direction, fvec w)
+void L1Fit::EnergyLossCorrectionIron(cnst& radThick, cnst& direction)
{
constexpr int atomicZ = 26;
constexpr float atomicA = 55.845f;
constexpr float rho = 7.87;
constexpr float radLen = 1.758f;
- EnergyLossCorrection<atomicZ>(atomicA, rho, radLen, radThick, direction, w);
+ EnergyLossCorrection<atomicZ>(atomicA, rho, radLen, radThick, direction);
}
-void L1Fit::EnergyLossCorrectionCarbon(const fvec& radThick, fvec direction, fvec w)
+void L1Fit::EnergyLossCorrectionCarbon(cnst& radThick, cnst& direction)
{
constexpr int atomicZ = 6;
constexpr float atomicA = 12.011f;
constexpr float rho = 2.265;
constexpr float radLen = 18.76f;
- EnergyLossCorrection<atomicZ>(atomicA, rho, radLen, radThick, direction, w);
+ EnergyLossCorrection<atomicZ>(atomicA, rho, radLen, radThick, direction);
}
-void L1Fit::EnergyLossCorrectionAl(const fvec& radThick, fvec direction, fvec w)
+void L1Fit::EnergyLossCorrectionAl(cnst& radThick, cnst& direction)
{
constexpr int atomicZ = 13;
constexpr float atomicA = 26.981f;
constexpr float rho = 2.70f;
constexpr float radLen = 2.265f;
- EnergyLossCorrection<atomicZ>(atomicA, rho, radLen, radThick, direction, w);
+ EnergyLossCorrection<atomicZ>(atomicA, rho, radLen, radThick, direction);
}
-void L1Fit::GetExtrapolatedXYline(const fvec& z, const L1FieldRegion& F, fvec& extrX, fvec& extrY, fvec Jx[6],
+void L1Fit::GetExtrapolatedXYline(cnst& z, const L1FieldRegion& F, fvec& extrX, fvec& extrY, fvec Jx[6],
fvec Jy[6]) const
{
// extrapolate track assuming it is straight (qp==0)
@@ -1187,9 +1189,9 @@ void L1Fit::GetExtrapolatedXYline(const fvec& z, const L1FieldRegion& F, fvec& e
cnst c_light(0.000299792458), c1(1.), c2i(0.5), c6i(1. / 6.), c12i(1. / 12.);
- const fvec& tx = fTr.tx;
- const fvec& ty = fTr.ty;
- fvec dz = z - fTr.z;
+ cnst& tx = fTr.tx;
+ cnst& ty = fTr.ty;
+ fvec dz = z - fTr.z;
fvec dz2 = dz * dz;
fvec dzc6i = dz * c6i;
@@ -1227,8 +1229,8 @@ void L1Fit::GetExtrapolatedXYline(const fvec& z, const L1FieldRegion& F, fvec& e
}
-void L1Fit::AddTargetToLine(const fvec& targX, const fvec& targY, const fvec& targZ,
- const L1XYMeasurementInfo& targXYInfo, const L1FieldRegion& F)
+void L1Fit::AddTargetToLine(cnst& targX, cnst& targY, cnst& targZ, const L1XYMeasurementInfo& targXYInfo,
+ const L1FieldRegion& F)
{
// Add the target constraint to a straight line track
@@ -1237,7 +1239,7 @@ void L1Fit::AddTargetToLine(const fvec& targX, const fvec& targY, const fvec& ta
FilterExtrapolatedXY(targX, targY, targXYInfo, eX, eY, Jx, Jy);
}
-fvec L1Fit::ApproximateBetheBloch(const fvec& bg2)
+fvec L1Fit::ApproximateBetheBloch(cnst& bg2)
{
//
// This is the parameterization of the Bethe-Bloch formula inspired by Geant.
@@ -1253,38 +1255,37 @@ fvec L1Fit::ApproximateBetheBloch(const fvec& bg2)
// The returned value is in [GeV/(g/cm^2)].
//
- const fvec kp0 = 2.33f;
- const fvec kp1 = 0.20f;
- const fvec kp2 = 3.00f;
- const fvec kp3 = 173e-9f;
- const fvec kp4 = 0.49848f;
+ cnst kp0 = 2.33f;
+ cnst kp1 = 0.20f;
+ cnst kp2 = 3.00f;
+ cnst kp3 = 173e-9f;
+ cnst kp4 = 0.49848f;
constexpr float mK = 0.307075e-3f; // [GeV*cm^2/g]
constexpr float _2me = 1.022e-3f; // [GeV/c^2]
- const fvec rho = kp0;
- const fvec x0 = kp1 * 2.303f;
- const fvec x1 = kp2 * 2.303f;
- const fvec mI = kp3;
- const fvec mZA = kp4;
- const fvec maxT = _2me * bg2; // neglecting the electron mass
+ cnst rho = kp0;
+ cnst x0 = kp1 * 2.303f;
+ cnst x1 = kp2 * 2.303f;
+ cnst mI = kp3;
+ cnst mZA = kp4;
+ cnst maxT = _2me * bg2; // neglecting the electron mass
//*** Density effect
fvec d2(0.f);
- const fvec x = 0.5f * log(bg2);
- const fvec lhwI = log(28.816f * 1e-9f * sqrt(rho * mZA) / mI);
+ cnst x = 0.5f * log(bg2);
+ cnst lhwI = log(28.816f * 1e-9f * sqrt(rho * mZA) / mI);
- fmask init = x > x1;
- d2 = iif(init, lhwI + x - 0.5f, fvec::Zero());
- const fvec r = (x1 - x) / (x1 - x0);
- init = (x > x0) & (x1 > x);
- d2 = iif(init, lhwI + x - 0.5f + (0.5f - lhwI - x0) * r * r * r, d2);
+ fmask init = x > x1;
+ d2 = iif(init, lhwI + x - 0.5f, fvec::Zero());
+ cnst r = (x1 - x) / (x1 - x0);
+ init = (x > x0) & (x1 > x);
+ d2 = iif(init, lhwI + x - 0.5f + (0.5f - lhwI - x0) * r * r * r, d2);
return mK * mZA * (fvec(1.f) + bg2) / bg2
* (0.5f * log(_2me * bg2 * maxT / (mI * mI)) - bg2 / (fvec(1.f) + bg2) - d2);
}
-fvec L1Fit::ApproximateBetheBloch(const fvec& bg2, const fvec& kp0, const fvec& kp1, const fvec& kp2, const fvec& kp3,
- const fvec& kp4)
+fvec L1Fit::ApproximateBetheBloch(cnst& bg2, cnst& kp0, cnst& kp1, cnst& kp2, cnst& kp3, cnst& kp4)
{
//
// This is the parameterization of the Bethe-Bloch formula inspired by Geant.
@@ -1300,31 +1301,31 @@ fvec L1Fit::ApproximateBetheBloch(const fvec& bg2, const fvec& kp0, const fvec&
// The returned value is in [GeV/(g/cm^2)].
//
- // const fvec &kp0 = 2.33f;
- // const fvec &kp1 = 0.20f;
- // const fvec &kp2 = 3.00f;
- // const fvec &kp3 = 173e-9f;
- // const fvec &kp4 = 0.49848f;
+ // cnst &kp0 = 2.33f;
+ // cnst &kp1 = 0.20f;
+ // cnst &kp2 = 3.00f;
+ // cnst &kp3 = 173e-9f;
+ // cnst &kp4 = 0.49848f;
constexpr float mK = 0.307075e-3f; // [GeV*cm^2/g]
constexpr float _2me = 1.022e-3f; // [GeV/c^2]
- const fvec& rho = kp0;
- const fvec x0 = kp1 * 2.303f;
- const fvec x1 = kp2 * 2.303f;
- const fvec& mI = kp3;
- const fvec& mZA = kp4;
- const fvec maxT = _2me * bg2; // neglecting the electron mass
+ cnst& rho = kp0;
+ cnst x0 = kp1 * 2.303f;
+ cnst x1 = kp2 * 2.303f;
+ cnst& mI = kp3;
+ cnst& mZA = kp4;
+ cnst maxT = _2me * bg2; // neglecting the electron mass
//*** Density effect
fvec d2(0.f);
- const fvec x = 0.5f * log(bg2);
- const fvec lhwI = log(28.816f * 1e-9f * sqrt(rho * mZA) / mI);
+ cnst x = 0.5f * log(bg2);
+ cnst lhwI = log(28.816f * 1e-9f * sqrt(rho * mZA) / mI);
- fmask init = x > x1;
- d2 = iif(init, lhwI + x - 0.5f, fvec::Zero());
- const fvec r = (x1 - x) / (x1 - x0);
- init = (x > x0) & (x1 > x);
- d2 = iif(init, lhwI + x - 0.5f + (0.5f - lhwI - x0) * r * r * r, d2);
+ fmask init = x > x1;
+ d2 = iif(init, lhwI + x - 0.5f, fvec::Zero());
+ cnst r = (x1 - x) / (x1 - x0);
+ init = (x > x0) & (x1 > x);
+ d2 = iif(init, lhwI + x - 0.5f + (0.5f - lhwI - x0) * r * r * r, d2);
return mK * mZA * (fvec(1.f) + bg2) / bg2
* (0.5f * log(_2me * bg2 * maxT / (mI * mI)) - bg2 / (fvec(1.f) + bg2) - d2);
@@ -1332,7 +1333,7 @@ fvec L1Fit::ApproximateBetheBloch(const fvec& bg2, const fvec& kp0, const fvec&
void L1Fit::FilterChi2XYC00C10C11(const L1UMeasurementInfo& info, fvec& x, fvec& y, fvec& C00, fvec& C10, fvec& C11,
- fvec& chi2, const fvec& u, const fvec& du2)
+ fvec& chi2, cnst& u, cnst& du2)
{
fvec wi, zeta, zetawi, HCH;
fvec F0, F1;
@@ -1361,8 +1362,8 @@ void L1Fit::FilterChi2XYC00C10C11(const L1UMeasurementInfo& info, fvec& x, fvec&
}
-void L1Fit::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)
+void L1Fit::FilterChi2(const L1UMeasurementInfo& info, cnst& x, cnst& y, cnst& C00, cnst& C10, cnst& C11, fvec& chi2,
+ cnst& u, cnst& du2)
{
fvec zeta, HCH;
fvec F0, F1;
diff --git a/reco/L1/L1Algo/L1Fit.h b/reco/L1/L1Algo/L1Fit.h
index 48cad6a908..4e95f91e7f 100644
--- a/reco/L1/L1Algo/L1Fit.h
+++ b/reco/L1/L1Algo/L1Fit.h
@@ -17,6 +17,8 @@
#include "L1UMeasurementInfo.h"
#include "L1XYMeasurementInfo.h"
+#define cnst const fvec
+
/// Kalman filter fitter for L1 tracking
///
class L1Fit {
@@ -32,6 +34,12 @@ public:
SetTrack(tr, C);
}
+ void SetMask(const fmask& m)
+ {
+ fMask = m;
+ fMaskF = iif(m, fvec::One(), fvec::Zero());
+ }
+
template<typename T>
void SetTrack(const T* tr, const T* C)
{
@@ -45,7 +53,7 @@ public:
fQp0 = fTr.qp;
}
- void SetQp0(const fvec& qp0) { fQp0 = qp0; }
+ void SetQp0(cnst& qp0) { fQp0 = qp0; }
L1TrackPar& Tr() { return fTr; }
@@ -73,55 +81,56 @@ public:
/// get the particle mass squared
fvec GetParticleMass2() const { return fMass2; }
- void Filter(const L1UMeasurementInfo& info, const fvec& u, const fvec& sigma2, const fvec& w);
- void FilterXY(const L1XYMeasurementInfo& info, fvec x, fvec y);
- void FilterTime(fvec t, fvec dt2, fvec w, fvec timeInfo);
- void FilterExtrapolatedXY(const fvec& x, const fvec& y, const L1XYMeasurementInfo& info, const fvec& extrX,
- const fvec& extrY, const fvec Jx[6], const fvec Jy[6]);
+ void Filter(const L1UMeasurementInfo& info, cnst& u, cnst& sigma2);
+ void FilterXY(const L1XYMeasurementInfo& info, cnst& x, cnst& y);
+ void FilterTime(cnst& t, cnst& dt2, cnst& timeInfo);
+ void FilterExtrapolatedXY(cnst& x, cnst& y, const L1XYMeasurementInfo& info, cnst& extrX, cnst& extrY, cnst Jx[6],
+ cnst Jy[6]);
- void Extrapolate(fvec z_out, const L1FieldRegion& F, const fvec& w);
+ void Extrapolate(cnst& z_out, const L1FieldRegion& F);
- void ExtrapolateStep(fvec z_out, const L1FieldRegion& F, const fvec& w);
- void ExtrapolateStepAnalytic(fvec z_out, const L1FieldRegion& F, const fvec& w);
- void ExtrapolateLine(fvec z_out, const fvec& w);
- void ExtrapolateLine(fvec z_out, const L1FieldRegion& F, const fvec& w);
+ void ExtrapolateStep(cnst& z_out, const L1FieldRegion& F);
+ void ExtrapolateStepAnalytic(cnst& z_out, const L1FieldRegion& F);
+ void ExtrapolateLine(cnst& z_out);
+ void ExtrapolateLine(cnst& z_out, const L1FieldRegion& F);
- void EnergyLossCorrection(const fvec& radThick, fvec upstreamDirection, fvec w);
+ void EnergyLossCorrection(cnst& radThick, cnst& upstreamDirection);
template<int atomicZ>
- void EnergyLossCorrection(float atomicA, float rho, float radLen, const fvec& radThick, fvec direction, fvec w);
+ void EnergyLossCorrection(float atomicA, float rho, float radLen, cnst& radThick, cnst& direction);
- void EnergyLossCorrectionIron(const fvec& radThick, fvec direction, fvec w);
- void EnergyLossCorrectionCarbon(const fvec& radThick, fvec direction, fvec w);
- void EnergyLossCorrectionAl(const fvec& radThick, fvec direction, fvec w);
+ void EnergyLossCorrectionIron(cnst& radThick, cnst& direction);
+ void EnergyLossCorrectionCarbon(cnst& radThick, cnst& direction);
+ void EnergyLossCorrectionAl(cnst& radThick, cnst& direction);
- void AddMsInMaterial(const fvec& radThick, fvec w);
+ void AddMsInMaterial(cnst& radThick);
- void AddMsInThickMaterial(fvec radThick, fvec w, fvec thickness, bool fDownstream);
+ void AddMsInThickMaterial(cnst& radThick, cnst& thickness, bool fDownstream);
- void GetExtrapolatedXYline(const fvec& z, const L1FieldRegion& F, fvec& extrX, fvec& extrY, fvec Jx[6],
- fvec Jy[6]) const;
+ void GetExtrapolatedXYline(cnst& z, const L1FieldRegion& F, fvec& extrX, fvec& extrY, fvec Jx[6], fvec Jy[6]) const;
- void ExtrapolateXC00Line(fvec z_out, fvec& x, fvec& C00) const;
- void ExtrapolateYC11Line(fvec z_out, fvec& y, fvec& C11) const;
- void ExtrapolateC10Line(fvec z_out, fvec& C10) const;
+ void ExtrapolateXC00Line(cnst& z_out, fvec& x, fvec& C00) const;
+ void ExtrapolateYC11Line(cnst& z_out, fvec& y, fvec& C11) const;
+ void ExtrapolateC10Line(cnst& z_out, fvec& C10) const;
- void AddTargetToLine(const fvec& targX, const fvec& targY, const fvec& targZ, const L1XYMeasurementInfo& targXYInfo,
+ void AddTargetToLine(cnst& targX, cnst& targY, cnst& targZ, const L1XYMeasurementInfo& targXYInfo,
const L1FieldRegion& F);
- static fvec ApproximateBetheBloch(const fvec& bg2);
+ static fvec ApproximateBetheBloch(cnst& bg2);
- static fvec ApproximateBetheBloch(const fvec& bg2, const fvec& kp0, const fvec& kp1, const fvec& kp2, const fvec& kp3,
- const fvec& kp4);
+ static fvec ApproximateBetheBloch(cnst& bg2, cnst& kp0, cnst& kp1, cnst& kp2, cnst& kp3, cnst& kp4);
static void FilterChi2XYC00C10C11(const L1UMeasurementInfo& info, fvec& x, fvec& y, fvec& C00, fvec& C10, fvec& C11,
- fvec& chi2, const fvec& u, const fvec& du2);
+ fvec& chi2, cnst& u, cnst& du2);
- 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);
+ static void FilterChi2(const L1UMeasurementInfo& info, cnst& x, cnst& y, cnst& C00, cnst& C10, cnst& C11, fvec& chi2,
+ cnst& u, cnst& du2);
private:
+ fmask fMask {fmask::One()}; // mask of active elements in simd vectors
+ fvec fMaskF {fvec::One()}; // floating point representation of fMask
+
L1TrackPar fTr {};
fvec fQp0 {0.};
@@ -130,7 +139,6 @@ private:
fvec fPipeRadThick {7.87e-3f}; // 0.7 mm Aluminium // TODO:
fvec fTargetRadThick {3.73e-2f * 2}; // 250 mum Gold // TODO:
-
} _fvecalignment;
// =============================================================================================
@@ -144,23 +152,27 @@ inline void L1Fit::SetOneEntry(const int i0, const L1Fit& T1, const int i1)
fQp0[i0] = T1.fQp0[i1];
}
-inline void L1Fit::ExtrapolateXC00Line(fvec z_out, fvec& x, fvec& C00) const
+inline void L1Fit::ExtrapolateXC00Line(cnst& z_out, fvec& x, fvec& C00) const
{
- const fvec dz = (z_out - fTr.z);
- x = fTr.x + fTr.tx * dz;
- C00 = fTr.C00 + dz * (2 * fTr.C20 + dz * fTr.C22);
+ cnst dz = (z_out - fTr.z);
+ x = fTr.x + fTr.tx * dz;
+ C00 = fTr.C00 + dz * (2 * fTr.C20 + dz * fTr.C22);
}
-inline void L1Fit::ExtrapolateYC11Line(fvec z_out, fvec& y, fvec& C11) const
+inline void L1Fit::ExtrapolateYC11Line(cnst& z_out, fvec& y, fvec& C11) const
{
- const fvec dz = (z_out - fTr.z);
- y = fTr.y + fTr.ty * dz;
- C11 = fTr.C11 + dz * (2 * fTr.C31 + dz * fTr.C33);
+ cnst dz = (z_out - fTr.z);
+ y = fTr.y + fTr.ty * dz;
+ C11 = fTr.C11 + dz * (2 * fTr.C31 + dz * fTr.C33);
}
-inline void L1Fit::ExtrapolateC10Line(fvec z_out, fvec& C10) const
+inline void L1Fit::ExtrapolateC10Line(cnst& z_out, fvec& C10) const
{
- const fvec dz = (z_out - fTr.z);
- C10 = fTr.C10 + dz * (fTr.C21 + fTr.C30 + dz * fTr.C32);
+ cnst dz = (z_out - fTr.z);
+ C10 = fTr.C10 + dz * (fTr.C21 + fTr.C30 + dz * fTr.C32);
}
+
+#undef cnst
+
+
#endif
diff --git a/reco/L1/L1Algo/L1TrackFitter.cxx b/reco/L1/L1Algo/L1TrackFitter.cxx
index 6029484410..1bc56c0b0d 100644
--- a/reco/L1/L1Algo/L1TrackFitter.cxx
+++ b/reco/L1/L1Algo/L1TrackFitter.cxx
@@ -88,8 +88,9 @@ void L1Algo::L1KFTrackFitter()
// fvec dt2_lst; /// TODO: Why are there two different variables for the time error on the last station?
fvec Sy[L1Constants::size::kMaxNstations];
- fvec w[L1Constants::size::kMaxNstations];
- fvec w_time[L1Constants::size::kMaxNstations]; // !!!
+ fmask w[L1Constants::size::kMaxNstations];
+ fmask w_time[L1Constants::size::kMaxNstations]; // !!!
+
fvec y_temp;
fvec x_temp;
fvec fldZ0;
@@ -123,8 +124,8 @@ void L1Algo::L1KFTrackFitter()
// get hits of current track
for (int ista = 0; ista < nStations; ista++) {
- w[ista] = ZERO;
- w_time[ista] = ZERO;
+ w[ista] = fmask::Zero();
+ w_time[ista] = fmask::Zero();
z[ista] = ZSta[ista];
}
@@ -143,8 +144,8 @@ void L1Algo::L1KFTrackFitter()
//if (sta[ista].fieldStatus) { isFieldPresent[iVec] = true; }
iSta[ih] = ista;
- w[ista][iVec] = 1.;
- if (sta[ista].timeInfo) { w_time[ista][iVec] = 1.; }
+ w[ista][iVec] = true;
+ if (sta[ista].timeInfo) { w_time[ista][iVec] = true; }
u[ista][iVec] = hit.u;
v[ista][iVec] = hit.v;
@@ -208,10 +209,9 @@ void L1Algo::L1KFTrackFitter()
if (kGlobal == fTrackingMode || kMcbm == fTrackingMode) {
GuessVecNoField(fit, x_last, x_2last, y_last, y_2last, z_end, z_2last, time_last, w_time, dt2_last);
}
-
- else
-
+ else {
GuessVec(fit, x, y, z, Sy, w, nStations, &z_end, time, w_time);
+ }
if (kGlobal == fTrackingMode || kMcbm == fTrackingMode) { tr.qp = fvec(1. / 1.1); }
@@ -225,8 +225,8 @@ void L1Algo::L1KFTrackFitter()
int ista = nStations - 1;
- time_last = iif(w_time[ista] > fvec::Zero(), time_last, fvec::Zero());
- dt2_last = iif(w_time[ista] > fvec::Zero(), dt2_last, fvec(1.e6));
+ time_last = iif(w_time[ista], time_last, fvec::Zero());
+ dt2_last = iif(w_time[ista], dt2_last, fvec(1.e6));
FilterFirst(fit, x_last, y_last, time_last, dt2_last, d_xx_lst, d_yy_lst, d_xy_lst);
@@ -251,18 +251,18 @@ void L1Algo::L1KFTrackFitter()
fld.Set(fldB0, fldZ0, fldB1, fldZ1, fldB2, fldZ2);
fmask initialised = (z[ista] < z_end) & (z_start <= z[ista]);
- fvec w1 = iif(initialised, w[ista], fvec::Zero());
- fvec wExtr = iif(initialised, fvec::One(), fvec::Zero());
- fvec w1_time = iif(initialised, w_time[ista], fvec::Zero());
fld1 = fld;
- fit.Extrapolate(z[ista], fld1, wExtr);
- fit.AddMsInMaterial(fParameters.GetMaterialThickness(ista, tr.x, tr.y), wExtr);
- fit.EnergyLossCorrection(fParameters.GetMaterialThickness(ista, tr.x, tr.y), fvec(1.f), wExtr);
- fit.Filter(sta[ista].frontInfo, u[ista], du2[ista], w1);
- fit.Filter(sta[ista].backInfo, v[ista], dv2[ista], w1);
- fit.FilterTime(time[ista], dt2[ista], w1_time, sta[ista].timeInfo);
+ fit.SetMask(initialised);
+ fit.Extrapolate(z[ista], fld1);
+ fit.AddMsInMaterial(fParameters.GetMaterialThickness(ista, tr.x, tr.y));
+ fit.EnergyLossCorrection(fParameters.GetMaterialThickness(ista, tr.x, tr.y), fvec(1.f));
+
+ fit.SetMask(initialised && w[ista]);
+ fit.Filter(sta[ista].frontInfo, u[ista], du2[ista]);
+ fit.Filter(sta[ista].backInfo, v[ista], dv2[ista]);
+ fit.FilterTime(time[ista], dt2[ista], sta[ista].timeInfo);
fldB2 = fldB1;
@@ -275,6 +275,8 @@ void L1Algo::L1KFTrackFitter()
L1Fit fitpv = fit;
{
+ fitpv.SetMask(fmask::One());
+
L1UMeasurementInfo vtxInfoX;
vtxInfoX.cos_phi = 1.;
vtxInfoX.sin_phi = 0.;
@@ -288,14 +290,14 @@ void L1Algo::L1KFTrackFitter()
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());
+ fitpv.Extrapolate(fParameters.GetTargetPositionZ(), fld);
+ fitpv.Filter(vtxInfoX, fParameters.GetTargetPositionX(), fvec(1.e-8));
+ fitpv.Filter(vtxInfoY, fParameters.GetTargetPositionY(), fvec(1.e-8));
}
}
else {
fitpv.SetQp0(fitpv.Tr().qp);
- fitpv.Extrapolate(fParameters.GetTargetPositionZ(), fld, fvec(1.));
+ fitpv.Extrapolate(fParameters.GetTargetPositionZ(), fld);
}
}
@@ -397,16 +399,15 @@ void L1Algo::L1KFTrackFitter()
fld.Set(fldB0, fldZ0, fldB1, fldZ1, fldB2, fldZ2);
fmask initialised = (z[ista] <= z_end) & (z_start < z[ista]);
- fvec w1 = iif(initialised, w[ista], fvec::Zero());
- fvec w1_time = iif(initialised, w_time[ista], fvec::Zero());
- fvec wExtr = iif(initialised, fvec::One(), fvec::Zero());
- fit.Extrapolate(z[ista], fld, w1);
- fit.AddMsInMaterial(fParameters.GetMaterialThickness(ista, tr.x, tr.y), wExtr);
- fit.EnergyLossCorrection(fParameters.GetMaterialThickness(ista, tr.x, tr.y), fvec(-1.f), wExtr);
- fit.Filter(sta[ista].frontInfo, u[ista], du2[ista], w1);
- fit.Filter(sta[ista].backInfo, v[ista], dv2[ista], w1);
- fit.FilterTime(time[ista], dt2[ista], w1_time, sta[ista].timeInfo);
+ fit.SetMask(initialised);
+ fit.Extrapolate(z[ista], fld);
+ fit.AddMsInMaterial(fParameters.GetMaterialThickness(ista, tr.x, tr.y));
+ fit.EnergyLossCorrection(fParameters.GetMaterialThickness(ista, tr.x, tr.y), fvec(-1.f));
+ fit.SetMask(initialised && w[ista]);
+ fit.Filter(sta[ista].frontInfo, u[ista], du2[ista]);
+ fit.Filter(sta[ista].backInfo, v[ista], dv2[ista]);
+ fit.FilterTime(time[ista], dt2[ista], sta[ista].timeInfo);
fldB2 = fldB1;
fldZ2 = fldZ1;
@@ -453,7 +454,7 @@ void L1Algo::L1KFTrackFitter()
}
}
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)
+ fvec& z_2last, fvec& time_last, fmask* /*w_time*/, fvec& dt2_last)
{
L1TrackPar& tr = track.Tr();
@@ -478,7 +479,7 @@ void L1Algo::GuessVecNoField(L1Fit& track, fvec& x_last, fvec& x_2last, fvec& y_
}
-void L1Algo::GuessVec(L1TrackPar& tr, 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, fmask* 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%).
{
@@ -489,14 +490,14 @@ void L1Algo::GuessVec(L1TrackPar& tr, fvec* xV, fvec* yV, fvec* zV, fvec* Sy, fv
if (zCur) z0 = *zCur;
else
z0 = zV[i];
- w = wV[i];
+ w = iif(wV[i], fvec::One(), fvec::Zero());
A0 = w;
a0 = w * xV[i];
b0 = w * yV[i];
for (i = 0; i < NHits; i++) {
x = xV[i];
y = yV[i];
- w = wV[i];
+ w = iif(wV[i], fvec::One(), fvec::Zero());
z = zV[i] - z0;
S = Sy[i];
wz = w * z;
@@ -546,8 +547,8 @@ void L1Algo::GuessVec(L1TrackPar& tr, fvec* xV, fvec* yV, fvec* zV, fvec* Sy, fv
tr.z = z0;
}
-void L1Algo::GuessVec(L1Fit& track, 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, fmask* wV, int NHits, fvec* zCur,
+ fvec* timeV, fmask* 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();
@@ -563,7 +564,7 @@ void L1Algo::GuessVec(L1Fit& track, fvec* xV, fvec* yV, fvec* zV, fvec* Sy, fvec
z0 = zV[i];
}
- w = wV[i];
+ w = iif(wV[i], fvec::One(), fvec::Zero());
A0 = w;
a0 = w * xV[i];
b0 = w * yV[i];
@@ -572,8 +573,9 @@ void L1Algo::GuessVec(L1Fit& track, fvec* xV, fvec* yV, fvec* zV, fvec* Sy, fvec
for (i = 0; i < NHits; i++) {
x = xV[i];
y = yV[i];
- w = wV[i];
- if (timeV) nhits = nhits + w_time[i];
+ w = iif(wV[i], fvec::One(), fvec::Zero());
+ ;
+ if (timeV) nhits(w_time[i]) += fvec::One();
z = zV[i] - z0;
S = Sy[i];
@@ -596,7 +598,7 @@ void L1Algo::GuessVec(L1Fit& track, fvec* xV, fvec* yV, fvec* zV, fvec* Sy, fvec
fvec dx = xV[i] - fParameters.GetTargetPositionX();
fvec dy = yV[i] - fParameters.GetTargetPositionY();
fvec dz = zV[i] - fParameters.GetTargetPositionZ();
- time += w_time[i] * (timeV[i] - sqrt(dx * dx + dy * dy + dz * dz) / fvec(30.));
+ time(w_time[i]) += (timeV[i] - sqrt(dx * dx + dy * dy + dz * dz) / fvec(30.));
}
}
diff --git a/reco/L1/ParticleFinder/CbmL1PFFitter.cxx b/reco/L1/ParticleFinder/CbmL1PFFitter.cxx
index 36bb7140a7..ab9b8c0aad 100644
--- a/reco/L1/ParticleFinder/CbmL1PFFitter.cxx
+++ b/reco/L1/ParticleFinder/CbmL1PFFitter.cxx
@@ -159,7 +159,8 @@ void CbmL1PFFitter::Fit(vector<CbmStsTrack>& Tracks, vector<int>& pidHypo)
fvec* dt2 = new fvec[nHits];
fvec* y = new fvec[nHits];
fvec* z = new fvec[nHits];
- fvec* w = new fvec[nHits];
+ fmask* w = new fmask[nHits];
+
// fvec y_temp;
fvec x_first, y_first, x_last, y_last;
fvec fstC00, fstC10, fstC11;
@@ -214,7 +215,7 @@ void CbmL1PFFitter::Fit(vector<CbmStsTrack>& Tracks, vector<int>& pidHypo)
// get hits of current track
for (i = 0; i < nHits; i++) {
- w[i] = ZERO;
+ w[i] = fmask::Zero();
z[i] = sta[i].fZ;
}
@@ -262,7 +263,7 @@ void CbmL1PFFitter::Fit(vector<CbmStsTrack>& Tracks, vector<int>& pidHypo)
dv2[ista][iVec] = hit->GetDv() * hit->GetDv();
dt2[ista][iVec] = hit->GetTimeError() * hit->GetTimeError();
}
- w[ista][iVec] = 1.f;
+ w[ista][iVec] = true;
x[ista][iVec] = posx;
y[ista][iVec] = posy;
@@ -321,16 +322,17 @@ void CbmL1PFFitter::Fit(vector<CbmStsTrack>& Tracks, vector<int>& pidHypo)
fld.Set(b0, z0, b1, z1, b2, z2);
fmask initialised = (z[i] <= z_end) & (z_start < z[i]);
- fvec w1 = iif(initialised, w[i], fvec::Zero());
- fvec wIn = iif(initialised, fvec::One(), fvec::Zero());
- fit.Extrapolate(z[i], fld, w1);
+ fit.SetMask(initialised);
+ fit.Extrapolate(z[i], fld);
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);
- fit.AddMsInMaterial(radThick, wIn);
- fit.EnergyLossCorrection(radThick, -fvec::One(), wIn);
+ fit.AddMsInMaterial(radThick);
+ fit.EnergyLossCorrection(radThick, -fvec::One());
+
+ fit.SetMask(initialised && w[i]);
+ fit.Filter(sta[i].frontInfo, u[i], du2[i]);
+ fit.Filter(sta[i].backInfo, v[i], dv2[i]);
+ fit.FilterTime(t[i], dt2[i], sta[i].timeInfo);
b2 = b1;
z2 = z1;
@@ -390,16 +392,17 @@ void CbmL1PFFitter::Fit(vector<CbmStsTrack>& Tracks, vector<int>& pidHypo)
fld.Set(b0, z0, b1, z1, b2, z2);
fmask initialised = (z[i] < z_end) & (z_start <= z[i]);
- fvec w1 = iif(initialised, w[i], fvec::Zero());
- fvec wIn = iif(initialised, fvec::One(), fvec::Zero());
- fit.Extrapolate(z[i], fld, w1);
+ fit.SetMask(initialised);
+ fit.Extrapolate(z[i], fld);
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);
- fit.AddMsInMaterial(radThick, wIn);
- fit.EnergyLossCorrection(radThick, fvec::One(), wIn);
+ fit.AddMsInMaterial(radThick);
+ fit.EnergyLossCorrection(radThick, fvec::One());
+
+ fit.SetMask(initialised && w[i]);
+ fit.Filter(sta[i].frontInfo, u[i], du2[i]);
+ fit.Filter(sta[i].backInfo, v[i], dv2[i]);
+ fit.FilterTime(t[i], dt2[i], sta[i].timeInfo);
b2 = b1;
z2 = z1;
@@ -554,17 +557,16 @@ void CbmL1PFFitter::GetChiToVertex(vector<CbmStsTrack>& Tracks, vector<PFFieldRe
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);
+ fit.SetMask(T.z > zSta[iSt] + fvec(2.5));
+ fit.Extrapolate(zSta[iSt], fld);
auto radThick = CbmL1::Instance()->fpAlgo->GetParameters()->GetMaterialThickness(iSt, fit.Tr().x, fit.Tr().y);
- fit.AddMsInMaterial(radThick, w);
- fit.EnergyLossCorrection(radThick, fvec::One(), w);
+ fit.AddMsInMaterial(radThick);
+ fit.EnergyLossCorrection(radThick, fvec::One());
}
- fit.Extrapolate(primVtx.GetRefZ(), fld, fvec::One());
- fit.AddMsInMaterial(fit.GetTargetRadThick(), fvec::One());
- fit.EnergyLossCorrection(fit.GetTargetRadThick(), fvec::One(), fvec::One());
+ fit.SetMask(fmask::One());
+ fit.Extrapolate(primVtx.GetRefZ(), fld);
+ fit.AddMsInMaterial(fit.GetTargetRadThick());
+ fit.EnergyLossCorrection(fit.GetTargetRadThick(), fvec::One());
Double_t Cv[3] = {primVtx.GetCovMatrix()[0], primVtx.GetCovMatrix()[1], primVtx.GetCovMatrix()[2]};
--
GitLab