From 7230f2d34d66fa579930afe7d0569e2bd8eb0d32 Mon Sep 17 00:00:00 2001
From: "se.gorbunov" <se.gorbunov@gsi.de>
Date: Tue, 13 Dec 2022 12:40:21 +0000
Subject: [PATCH] L1: fit: moving to L1TrackParFit utility..
---
reco/L1/L1Algo/L1Algo.h | 2 +
reco/L1/L1Algo/L1CATrackFinder.cxx | 140 ++++++++++-------------------
reco/L1/L1Algo/L1TrackFitter.cxx | 4 +-
reco/L1/L1Algo/L1TrackParFit.cxx | 43 +++++----
reco/L1/L1Algo/L1TrackParFit.h | 4 +-
5 files changed, 76 insertions(+), 117 deletions(-)
diff --git a/reco/L1/L1Algo/L1Algo.h b/reco/L1/L1Algo/L1Algo.h
index abb91bb991..8aa411d922 100644
--- a/reco/L1/L1Algo/L1Algo.h
+++ b/reco/L1/L1Algo/L1Algo.h
@@ -548,6 +548,8 @@ private:
fvec fTargY {L1Utils::kNaN}; ///< target position y coordinate for the current iteration (modifiable)
fvec fTargZ {L1Utils::kNaN}; ///< target position z coordinate for the current iteration (modifiable)
+ bool fIsTargetField {false}; ///< is the magnetic field present at the target
+
L1FieldValue fTargB _fvecalignment {}; // field in the target point (modifiable, do not touch!!)
L1XYMeasurementInfo TargetXYInfo _fvecalignment {}; // target constraint [cm]
diff --git a/reco/L1/L1Algo/L1CATrackFinder.cxx b/reco/L1/L1Algo/L1CATrackFinder.cxx
index 8e909cf558..dfe5662112 100644
--- a/reco/L1/L1Algo/L1CATrackFinder.cxx
+++ b/reco/L1/L1Algo/L1CATrackFinder.cxx
@@ -213,21 +213,16 @@ inline void L1Algo::findSingletsStep1( /// input 1st stage of singlet search
const L1Station& fld1Sta1 = fParameters.GetStation(fld1Ista1);
const L1Station& fld1Sta2 = fParameters.GetStation(fld1Ista2);
- L1Fit fitOld;
L1TrackParFit fit;
fit.fQp0 = fvec(0.);
- if (fpCurrentIteration->GetElectronFlag()) {
- fitOld.SetParticleMass(L1Constants::phys::kElectronMass);
- fit.SetParticleMass(L1Constants::phys::kElectronMass);
- }
+ if (fpCurrentIteration->GetElectronFlag()) { fit.SetParticleMass(L1Constants::phys::kElectronMass); }
else {
- fitOld.SetParticleMass(fDefaultMass);
fit.SetParticleMass(fDefaultMass);
}
for (int i1_V = 0; i1_V < n1_V; i1_V++) {
- L1TrackPar& T = fit1.fTr;
+ 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
@@ -295,41 +290,17 @@ inline void L1Algo::findSingletsStep1( /// input 1st stage of singlet search
std::tie(T.C00, T.C10, T.C11) = stal.FormXYCovarianceMatrix(du2_l[i1_V], dv2_l[i1_V]);
- //assert(T.IsConsistent(true, -1));
-
- // add the target
- {
- fvec eX, eY, J04, J14;
- fvec dz;
- dz = fTargZ - zl;
- L1ExtrapolateJXY0(T.tx, T.ty, dz, fld0, eX, eY, J04, J14);
- fvec J[6];
- J[0] = dz;
- J[1] = 0;
- J[2] = J04;
- J[3] = 0;
- J[4] = dz;
- J[5] = J14;
- L1FilterVtx(T, fTargX, fTargY, TargetXYInfo, eX, eY, J);
- }
+ fit.AddTargetToLine(fTargX, fTargY, fTargZ, TargetXYInfo, fld0);
}
- //assert(T.IsConsistent(true, -1));
-
- if (kMcbm == fTrackingMode) {
- fitOld.L1AddThickMaterial(T, fParameters.GetMaterialThickness(istal, T.x, T.y), fMaxInvMom, fvec::One(),
- stal.fZthick, 1);
- }
- else {
- fitOld.L1AddMaterial(T, fParameters.GetMaterialThickness(istal, T.x, T.y), fMaxInvMom, fvec::One());
- }
+ fit.AddMsInMaterial(fParameters.GetMaterialThickness(istal, T.x, T.y), fMaxInvMom, fvec::One());
//if ((istam >= fNstationsBeforePipe) && (istal <= fNstationsBeforePipe - 1)) {
- //fitOld.L1AddPipeMaterial(T, fMaxInvMom, fvec::One());
+ //fit.L1AddPipeMaterial(T, fMaxInvMom, fvec::One());
//}
fvec dz = stam.fZ - zl;
- L1ExtrapolateTime(T, dz, stam.timeInfo);
+ L1ExtrapolateTime(T, dz, fvec::One());
// extrapolate to the middle hit
L1Extrapolate0(T, stam.fZ, fld0);
@@ -487,13 +458,14 @@ inline void L1Algo::findDoubletsStep0(
/// Add the middle hits to parameters estimation. Propagate to right station.
/// Find the triplets(right hit). Reformat data in the portion of triplets.
inline void L1Algo::findTripletsStep0( // input
- L1HitPoint* vHits_r, int /*iStaL*/, int iStaM, int iStaR, L1HitPoint* vHits_m, L1TrackPar* T_1, L1FieldRegion* fld_1,
+ L1HitPoint* vHits_r, int iStaL, int iStaM, int iStaR, L1HitPoint* vHits_m, L1TrackPar* T_1, L1FieldRegion* fld_1,
L1HitIndex_t* hitsl_1, Tindex n2, L1Vector<L1HitIndex_t>& hitsm_2, L1Vector<L1HitIndex_t>& i1_2,
// output
Tindex& n3, L1Vector<L1TrackPar>& T_3, L1Vector<L1HitIndex_t>& hitsl_3, L1Vector<L1HitIndex_t>& hitsm_3,
L1Vector<L1HitIndex_t>& hitsr_3, L1Vector<fvec>& u_front_3, L1Vector<fvec>& u_back_3, L1Vector<fvec>& z_Pos_3,
L1Vector<fvec>& du2_3, L1Vector<fvec>& dv2_3, L1Vector<fvec>& t_3, L1Vector<fvec>& dt2_3)
{
+ const L1Station& stal = fParameters.GetStation(iStaL);
const L1Station& stam = fParameters.GetStation(iStaM);
const L1Station& star = fParameters.GetStation(iStaR);
@@ -511,10 +483,12 @@ inline void L1Algo::findTripletsStep0( // input
L1TrackPar_0.C55 = 1.f;
*/
- L1Fit fitOld;
+ bool isMomentumFitted = (fIsTargetField || (stal.fieldStatus != 0) || (stam.fieldStatus != 0));
+ //bool isTimeFitted = ((stal.timeInfo != 0) || (stam.timeInfo != 0));
+
L1TrackParFit fit;
- fitOld.SetParticleMass(fDefaultMass);
fit.SetParticleMass(fDefaultMass);
+ fit.fQp0 = fvec(0.);
n3 = 0;
Tindex n3_V = 0, n3_4 = 0;
@@ -578,42 +552,25 @@ inline void L1Algo::findTripletsStep0( // input
hitsm_2_tmp[n2_4] = hitsm_2[i2];
} // n2_4
- fvec dz = zPos_2 - T2.z;
-
- L1ExtrapolateTime(T2, dz, stam.timeInfo);
+ // add the middle hit
- // assert(T2.IsConsistent(true, n2_4));
-
- // add middle hit
- L1ExtrapolateLine(T2, zPos_2);
-
- // assert(T2.IsConsistent(true, n2_4));
-
- // L1TrackPar tStore1 = T2;
-
- L1Filter(T2, stam.frontInfo, u_front_2, du2_2, fvec::One());
- L1Filter(T2, stam.backInfo, u_back_2, dv2_2, fvec::One());
+ fit.ExtrapolateLine(zPos_2, fvec::One());
+ 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);
FilterTime(T2, t_2, dt2_2, stam.timeInfo);
- if (kMcbm == fTrackingMode) {
- fitOld.L1AddThickMaterial(T2, fParameters.GetMaterialThickness(iStaM, T2.x, T2.y), fMaxInvMom, fvec::One(),
- stam.fZthick, 1);
- }
- else if (kGlobal == fTrackingMode) {
- fitOld.L1AddMaterial(T2, fParameters.GetMaterialThickness(iStaM, T2.x, T2.y), fMaxInvMom, fvec::One());
- }
- else {
- fitOld.L1AddMaterial(T2, fParameters.GetMaterialThickness(iStaM, T2.x, T2.y), T2.qp, fvec::One());
- }
-
- //if ((iStaR >= fNstationsBeforePipe) && (iStaM <= fNstationsBeforePipe - 1)) { fitOld.L1AddPipeMaterial(T2, T2.qp, 1); }
+ fit.AddMsInMaterial(fParameters.GetMaterialThickness(iStaM, T2.x, T2.y), isMomentumFitted ? T2.qp : fMaxInvMom,
+ fvec::One());
- fvec dz2 = star.fZ - T2.z;
- L1ExtrapolateTime(T2, dz2, stam.timeInfo);
+ //if ((iStaR >= fNstationsBeforePipe) && (iStaM <= fNstationsBeforePipe - 1)) { fit.L1AddPipeMaterial(T2, T2.qp, 1); }
// extrapolate to the right hit station
+ //fit.Extrapolate(star.fZ, T2.qp, f2,fvec::One());
+ fvec dz2 = star.fZ - T2.z;
+ L1ExtrapolateTime(T2, dz2, fvec::One());
L1Extrapolate(T2, star.fZ, T2.qp, f2);
// assert(T2.IsConsistent(true, n2_4));
@@ -668,12 +625,14 @@ inline void L1Algo::findTripletsStep0( // input
auto [xr, yr] = star.ConvUVtoXY<fscal>(hitr.U(), hitr.V());
- L1TrackPar T_cur = T2;
+ L1TrackParFit fit3;
+ fit3.SetParticleMass(fDefaultMass);
+ fit3.fQp0 = T2.qp;
- fvec dz3 = zr - T_cur.z;
- L1ExtrapolateTime(T_cur, dz3, star.timeInfo);
+ L1TrackPar& T_cur = fit3.fTr;
+ T_cur = T2;
- L1ExtrapolateLine(T_cur, zr);
+ fit3.ExtrapolateLine(zr, fvec::One());
if ((star.timeInfo) && (stam.timeInfo))
if (fabs(T_cur.t[i2_4] - hitr.T()) > sqrt(T_cur.C55[i2_4] + hitr.dT2()) * 5) continue;
@@ -787,21 +746,22 @@ inline void L1Algo::findTripletsStep1( // input
// L1TrackPar *T_3
L1Vector<L1TrackPar>& T_3)
{
+ L1TrackParFit fit;
+ fit.SetParticleMass(fDefaultMass);
for (Tindex i3_V = 0; i3_V < n3_V; ++i3_V) {
- fvec dz = z_Pos[i3_V] - T_3[i3_V].z;
-
- L1TrackPar& T3 = T_3[i3_V];
+ L1TrackPar& T3 = fit.fTr;
+ T3 = T_3[i3_V];
+ fit.fQp0 = fit.fTr.qp;
- L1ExtrapolateTime(T3, dz, star.timeInfo);
+ fit.ExtrapolateLine(z_Pos[i3_V], fvec::One());
- L1ExtrapolateLine(T3, z_Pos[i3_V]);
+ 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());
- L1Filter(T3, star.frontInfo, u_front_[i3_V], du2_3[i3_V], fvec::One());
- L1Filter(T3, star.backInfo, u_back_[i3_V], dv2_3[i3_V], fvec::One());
-
- if (kMcbm != fTrackingMode) { FilterTime(T3, t_3[i3_V], dt2_3[i3_V], star.timeInfo); }
+ if (kMcbm != fTrackingMode) { fit.FilterTime(t_3[i3_V], dt2_3[i3_V], fvec::One(), star.timeInfo); }
+ T_3[i3_V] = T3;
}
}
@@ -929,24 +889,12 @@ inline void L1Algo::findTripletsStep2(Tindex n3, int istal, int istam, int istar
fit.Filter(sta[ih0].backInfo, v[ih0], dv2[ih0], fvec::One());
fit.FilterTime(t[ih0], dt2[ih0], fvec::One(), sta[ih0].timeInfo);
- { // add the target constraint
- fvec eX, eY, J04, J14;
- fvec dz;
- dz = fTargZ - T.z;
- L1ExtrapolateJXY0(T.tx, T.ty, dz, fldTarget, eX, eY, J04, J14);
- fvec J[6];
- J[0] = dz;
- J[1] = 0;
- J[2] = J04;
- J[3] = 0;
- J[4] = dz;
- J[5] = J14;
- L1FilterVtx(T, fTargX, fTargY, TargetXYInfo, eX, eY, J);
- }
+ // add the target constraint
+ fit.AddTargetToLine(fTargX, fTargY, fTargZ, TargetXYInfo, fldTarget);
for (int ih = 1; ih < NHits; ++ih) {
fit.Extrapolate(z[ih], qp0, fld, fvec::One());
- fit.AddMaterial(fParameters.GetMaterialThickness(ista[ih], T.x, T.y), qp0, fvec::One());
+ fit.AddMsInMaterial(fParameters.GetMaterialThickness(ista[ih], T.x, T.y), qp0, fvec::One());
//if (ista[ih] == fNstationsBeforePipe) { fit.AddPipeMaterial(qp0, fvec::One()); }
fit.Filter(sta[ih].frontInfo, u[ih], du2[ih], fvec::One());
fit.Filter(sta[ih].backInfo, v[ih], dv2[ih], fvec::One());
@@ -979,7 +927,7 @@ inline void L1Algo::findTripletsStep2(Tindex n3, int istal, int istam, int istar
for (int ih = NHits - 2; ih >= 0; --ih) {
fit.Extrapolate(z[ih], qp0, fld, fvec::One());
- fit.AddMaterial(fParameters.GetMaterialThickness(ista[ih], T.x, T.y), qp0, fvec::One());
+ fit.AddMsInMaterial(fParameters.GetMaterialThickness(ista[ih], T.x, T.y), qp0, fvec::One());
//if (ista[ih] == fNstationsBeforePipe - 1) { fit.AddPipeMaterial(qp0, fvec::One()); }
fit.Filter(sta[ih].frontInfo, u[ih], du2[ih], fvec::One());
fit.Filter(sta[ih].backInfo, v[ih], dv2[ih], fvec::One());
@@ -1759,6 +1707,8 @@ void L1Algo::CATrackFinder()
fParameters.GetStation(0).fieldSlice.GetFieldValue(0, 0, fTargB);
} // NOTE: calculates field fTargB in the center of 0th station
+ fIsTargetField = (fabs(fTargB.y[0]) > 0.001);
+
TargetXYInfo.C00 = SigmaTargetX * SigmaTargetX;
TargetXYInfo.C10 = 0;
TargetXYInfo.C11 = SigmaTargetY * SigmaTargetY;
diff --git a/reco/L1/L1Algo/L1TrackFitter.cxx b/reco/L1/L1Algo/L1TrackFitter.cxx
index 77329ba97a..0c17271136 100644
--- a/reco/L1/L1Algo/L1TrackFitter.cxx
+++ b/reco/L1/L1Algo/L1TrackFitter.cxx
@@ -533,7 +533,7 @@ void L1Algo::L1KFTrackFitter()
//fit.AddPipeMaterial(qp01, wExtr);
//fit.EnergyLossCorrection(fit.fPipeRadThick, qp01, fvec(1.f), wExtr);
//}
- fit.AddMaterial(fParameters.GetMaterialThickness(ista, tr.x, tr.y), qp01, wExtr);
+ fit.AddMsInMaterial(fParameters.GetMaterialThickness(ista, tr.x, tr.y), qp01, wExtr);
fit.EnergyLossCorrection(fParameters.GetMaterialThickness(ista, tr.x, tr.y), qp01, fvec(1.f), wExtr);
fit.Filter(sta[ista].frontInfo, u[ista], du2[ista], w1);
@@ -683,7 +683,7 @@ void L1Algo::L1KFTrackFitter()
//fit.AddPipeMaterial(qp01, wExtr);
//fit.EnergyLossCorrection(fit.fPipeRadThick, qp01, fvec(-1.f), wExtr);
//}
- fit.AddMaterial(fParameters.GetMaterialThickness(ista, tr.x, tr.y), qp01, wExtr);
+ fit.AddMsInMaterial(fParameters.GetMaterialThickness(ista, tr.x, tr.y), qp01, wExtr);
fit.EnergyLossCorrection(fParameters.GetMaterialThickness(ista, tr.x, tr.y), qp01, fvec(-1.f), wExtr);
fit.Filter(sta[ista].frontInfo, u[ista], du2[ista], w1);
diff --git a/reco/L1/L1Algo/L1TrackParFit.cxx b/reco/L1/L1Algo/L1TrackParFit.cxx
index 719980e2d3..ef3ba7252e 100644
--- a/reco/L1/L1Algo/L1TrackParFit.cxx
+++ b/reco/L1/L1Algo/L1TrackParFit.cxx
@@ -932,12 +932,35 @@ void L1TrackParFit::ExtrapolateLine(fvec z_out, const fvec& w)
{
// extrapolate the track assuming fQp0 == 0
//
+ // it is a copy of a sequence two routines
+ // L1ExtrapolateTime() and L1ExtrapolateLine()
+ // TODO: make it right
cnst c_light(29.9792458);
fvec dz = (z_out - fTr.z);
dz.setZero(w <= fvec::Zero());
L1TrackPar& T = fTr;
+
+ // extrapolate time
+ 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);
+
+ fvec c52 = T.C52;
+ fvec c53 = T.C53;
+
+ T.C50 += k1 * T.C20 + k2 * T.C30;
+ T.C51 += k1 * T.C21 + k2 * T.C31;
+ T.C52 += k1 * T.C22 + k2 * T.C32;
+ T.C53 += k1 * T.C32 + k2 * T.C33;
+ T.C54 += k1 * T.C42 + k2 * T.C43;
+ T.C55 += k1 * (T.C52 + c52) + k2 * (T.C53 + c53);
+
+ // extrapolate trajectory
+
T.x += T.tx * dz;
T.y += T.ty * dz;
T.z += dz;
@@ -960,26 +983,10 @@ void L1TrackParFit::ExtrapolateLine(fvec z_out, const fvec& w)
T.C40 += dz * T.C42;
T.C41 += dz * T.C43;
-
- 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);
-
- fvec c52 = T.C52;
- fvec c53 = T.C53;
-
- T.C50 += k1 * T.C20 + k2 * T.C30;
- T.C51 += k1 * T.C21 + k2 * T.C31;
- T.C52 += k1 * T.C22 + k2 * T.C32;
- T.C53 += k1 * T.C32 + k2 * T.C33;
- T.C54 += k1 * T.C42 + k2 * T.C43;
- T.C55 += k1 * (T.C52 + c52) + k2 * (T.C53 + c53);
}
-void L1TrackParFit::AddMaterial(const fvec& radThick, fvec qp0, fvec w)
+void L1TrackParFit::AddMsInMaterial(const fvec& radThick, fvec qp0, fvec w)
{
cnst ONE = 1.;
@@ -1004,7 +1011,7 @@ void L1TrackParFit::AddMaterial(const fvec& radThick, fvec qp0, fvec w)
fTr.C33 += w * (ONE + tyty) * a;
}
-void L1TrackParFit::AddThickMaterial(fvec radThick, fvec qp0, fvec w, fvec thickness, bool fDownstream)
+void L1TrackParFit::AddMsInThickMaterial(fvec radThick, fvec qp0, fvec w, fvec thickness, bool fDownstream)
{
cnst ONE = 1.;
diff --git a/reco/L1/L1Algo/L1TrackParFit.h b/reco/L1/L1Algo/L1TrackParFit.h
index 2bbce0de61..217351c5f2 100644
--- a/reco/L1/L1Algo/L1TrackParFit.h
+++ b/reco/L1/L1Algo/L1TrackParFit.h
@@ -68,8 +68,8 @@ public:
void EnergyLossCorrectionCarbon(const fvec& radThick, fvec& qp0, fvec direction, fvec w);
void EnergyLossCorrectionAl(const fvec& radThick, fvec& qp0, fvec direction, fvec w);
- void AddMaterial(const fvec& radThick, fvec qp0, fvec w);
- void AddThickMaterial(fvec radThick, fvec qp0, fvec w, fvec thickness, bool fDownstream);
+ void AddMsInMaterial(const fvec& radThick, fvec qp0, fvec w);
+ void AddMsInThickMaterial(fvec radThick, fvec qp0, fvec w, fvec thickness, bool fDownstream);
void GetExtrapolatedXYline(const fvec& z, const L1FieldRegion& F, fvec& extrX, fvec& extrY, fvec Jx[6],
fvec Jy[6]) const;
--
GitLab