From ae29ad9045d3d6746c42764c6ac96271136486cd Mon Sep 17 00:00:00 2001
From: "se.gorbunov" <se.gorbunov@gsi.de>
Date: Tue, 8 Nov 2022 17:59:55 +0000
Subject: [PATCH] L1: store squared hit errors
---
reco/L1/CbmL1ReadEvent.cxx | 13 +-
reco/L1/L1Algo/L1Algo.cxx | 4 +-
reco/L1/L1Algo/L1Algo.h | 18 +--
reco/L1/L1Algo/L1BranchExtender.cxx | 18 +--
reco/L1/L1Algo/L1CATrackFinder.cxx | 182 ++++++++++++++--------------
reco/L1/L1Algo/L1Filtration.h | 4 +-
reco/L1/L1Algo/L1Hit.h | 12 +-
reco/L1/L1Algo/L1HitPoint.h | 152 +++++------------------
reco/L1/L1Algo/L1Station.h | 17 ++-
reco/L1/L1Algo/L1TrackFitter.cxx | 101 +++++++--------
reco/L1/L1Algo/L1TrackParFit.cxx | 4 +-
reco/L1/L1Algo/L1TrackParFit.h | 2 +-
12 files changed, 214 insertions(+), 313 deletions(-)
diff --git a/reco/L1/CbmL1ReadEvent.cxx b/reco/L1/CbmL1ReadEvent.cxx
index 8c6e07cc59..c4695c66df 100644
--- a/reco/L1/CbmL1ReadEvent.cxx
+++ b/reco/L1/CbmL1ReadEvent.cxx
@@ -1261,17 +1261,18 @@ void CbmL1::ReadEvent(float& TsStart, float& TsLength, float& /*TsOverlap*/, int
assert(th.iStripB >= 0 || th.iStripB < NStrips);
L1Hit h;
+ h.iSt = th.iStation;
h.f = th.iStripF;
h.b = th.iStripB;
h.ID = th.id;
- h.t = th.time;
- h.dt = th.dt;
- h.du = th.du;
- h.dv = th.dv;
+ h.z = th.z;
h.u = th.u;
h.v = th.v;
- h.z = th.z;
- h.iSt = th.iStation;
+ h.t = th.time;
+ h.dt2 = th.dt * th.dt;
+ h.du2 = th.du * th.du;
+ h.dv2 = th.dv * th.dv;
+
// save hit
fvExternalHits.push_back(CbmL1Hit(iHit, th.ExtIndex, th.Det));
diff --git a/reco/L1/L1Algo/L1Algo.cxx b/reco/L1/L1Algo/L1Algo.cxx
index 71fc7a1a9b..4699109dc3 100644
--- a/reco/L1/L1Algo/L1Algo.cxx
+++ b/reco/L1/L1Algo/L1Algo.cxx
@@ -180,12 +180,12 @@ L1HitPoint L1Algo::CreateHitPoint(const L1Hit& hit)
{
/// full the hit point by hit information: takes hit as input (2 strips)
/// and creates hit_point with all coordinates (x,y,z,u,v,t);
- return L1HitPoint(hit.z, hit.u, hit.v, hit.du, hit.dv, hit.t, hit.dt);
+ return L1HitPoint(hit.z, hit.u, hit.v, hit.t, hit.du2, hit.dv2, hit.dt2);
}
void L1Algo::CreateHitPoint(const L1Hit& hit, L1HitPoint& point)
{
- point.Set(hit.z, hit.u, hit.v, hit.du, hit.dv, hit.t, hit.dt);
+ point.Set(hit.z, hit.u, hit.v, hit.t, hit.du2, hit.dv2, hit.dt2);
}
int L1Algo::GetMcTrackIdForHit(int iHit) const
diff --git a/reco/L1/L1Algo/L1Algo.h b/reco/L1/L1Algo/L1Algo.h
index f30ce0aaed..f00d23276e 100644
--- a/reco/L1/L1Algo/L1Algo.h
+++ b/reco/L1/L1Algo/L1Algo.h
@@ -264,16 +264,16 @@ public:
void findSingletsStep0( // input
Tindex start_lh, Tindex n1_l, L1HitPoint* Hits_l,
// output
- fvec* u_front_l, fvec* u_back_l, fvec* zPos_l, L1HitIndex_t* hitsl, fvec* HitTime_l, fvec* HitTimeEr, fvec* d_u,
- fvec* d_v);
+ fvec* u_front_l, fvec* u_back_l, fvec* zPos_l, L1HitIndex_t* hitsl, fvec* t_l, fvec* dt2_l, fvec* du2_l,
+ fvec* dv2_l);
/// Get the field approximation. Add the target to parameters estimation. Propagate to middle station.
void findSingletsStep1( // input
int istal, int istam, Tindex n1_V,
- fvec* u_front_l, fvec* u_back_l, fvec* zPos_l, fvec* HitTime_l, fvec* HitTimeEr,
+ fvec* u_front_l, fvec* u_back_l, fvec* zPos_l, fvec* t_l, fvec* dt2_l,
// output
- L1TrackPar* T_1, L1FieldRegion* fld_1, fvec* d_u, fvec* d_v);
+ L1TrackPar* T_1, L1FieldRegion* fld_1, fvec* du2_l, fvec* dv2_l);
/// Find the doublets. Reformat data in the portion of doublets.
void findDoubletsStep0( // input
@@ -298,12 +298,12 @@ public:
// 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>& du_, L1Vector<fvec>& dv_, L1Vector<fvec>& timeR, L1Vector<fvec>& timeER);
+ L1Vector<fvec>& du2_3, L1Vector<fvec>& dv2_3, L1Vector<fvec>& t_3, L1Vector<fvec>& dt2_3);
/// Add the right hits to parameters estimation.
void findTripletsStep1( // input
Tindex n3_V, const L1Station& star, L1Vector<fvec>& u_front_3, L1Vector<fvec>& u_back_3, L1Vector<fvec>& z_Pos_3,
- L1Vector<fvec>& du_, L1Vector<fvec>& dv_, L1Vector<fvec>& timeR, L1Vector<fvec>& timeER,
+ L1Vector<fvec>& du2_3, L1Vector<fvec>& dv2_3, L1Vector<fvec>& t_3, L1Vector<fvec>& dt2_3,
// output
L1Vector<L1TrackPar>& T_3);
@@ -346,11 +346,11 @@ public:
void FilterFirst(L1TrackPar& track, fvec& x, fvec& y, L1Station& st);
void FilterFirst(L1TrackParFit& track, fvec& x, fvec& y, fvec& t, L1Station& st);
- void FilterFirst(L1TrackParFit& track, fvec& x, fvec& y, fvec& t, fvec& t_er, L1Station& st);
+ void FilterFirst(L1TrackParFit& track, fvec& x, fvec& y, fvec& t, fvec& dt2, L1Station& st);
- void FilterFirst(L1TrackParFit& track, fvec& x, fvec& y, fvec& t, fvec& t_er, L1Station& st, fvec& dx, fvec& dy,
+ void FilterFirst(L1TrackParFit& track, fvec& x, fvec& y, fvec& t, fvec& dt2, L1Station& st, fvec& dx2, fvec& dy2,
fvec& dxy);
- void FilterFirstL(L1TrackParFit& track, fvec& x, fvec& y, fvec& t, fvec& t_er, L1Station& st, fvec& dx, fvec& dy,
+ void FilterFirstL(L1TrackParFit& track, fvec& x, fvec& y, fvec& t, fvec& dt2, L1Station& st, fvec& dx2, fvec& dy2,
fvec& dxy);
diff --git a/reco/L1/L1Algo/L1BranchExtender.cxx b/reco/L1/L1Algo/L1BranchExtender.cxx
index 73bed78e94..3ed23c7194 100644
--- a/reco/L1/L1Algo/L1BranchExtender.cxx
+++ b/reco/L1/L1Algo/L1BranchExtender.cxx
@@ -88,7 +88,7 @@ void L1Algo::BranchFitterFast(const L1Branch& t, L1TrackPar& T, const bool dir,
T.C10 = sta0.XYInfo.C10;
T.C11 = sta0.XYInfo.C11;
- if (fUseHitErrors) { std::tie(T.C00, T.C10, T.C11) = sta0.FormXYCovarianceMatrix(hit0.du, hit0.dv); }
+ if (fUseHitErrors) { std::tie(T.C00, T.C10, T.C11) = sta0.FormXYCovarianceMatrix(hit0.du2, hit0.dv2); }
T.C20 = T.C21 = 0;
T.C30 = T.C31 = T.C32 = 0;
@@ -96,7 +96,7 @@ void L1Algo::BranchFitterFast(const L1Branch& t, L1TrackPar& T, const bool dir,
T.C50 = T.C51 = T.C52 = T.C53 = T.C54 = 0;
T.C22 = T.C33 = vINF;
T.C44 = 1.;
- T.C55 = hit0.dt * hit0.dt;
+ T.C55 = hit0.dt2;
L1FieldValue fldB0, fldB1, fldB2 _fvecalignment;
L1FieldRegion fld _fvecalignment;
@@ -141,15 +141,15 @@ void L1Algo::BranchFitterFast(const L1Branch& t, L1TrackPar& T, const bool dir,
L1UMeasurementInfo info = sta.frontInfo;
- if (fUseHitErrors) { info.sigma2 = hit.du * hit.du; }
+ if (fUseHitErrors) { info.sigma2 = hit.du2; }
L1Filter(T, info, u);
info = sta.backInfo;
- if (fUseHitErrors) { info.sigma2 = hit.dv * hit.dv; }
+ if (fUseHitErrors) { info.sigma2 = hit.dv2; }
L1Filter(T, info, v);
- FilterTime(T, hit.t, hit.dt);
+ FilterTime(T, hit.t, hit.dt2);
fldB0 = fldB1;
fldB1 = fldB2;
@@ -269,7 +269,7 @@ void L1Algo::FindMoreHits(L1Branch& t, L1TrackPar& T, const bool dir,
ih += HitsUnusedStartIndex[ista];
const L1Hit& hit = (*vHitsUnused)[ih];
//TODO: bug, it should be hit.dt*hit.dt
- if (fabs(hit.t - T.t[0]) > sqrt(T.C55[0] + hit.dt * hit.dt) * 5) continue;
+ if (fabs(hit.t - T.t[0]) > sqrt(T.C55[0] + hit.dt2) * 5) continue;
//if (GetFUsed((*fStripFlag)[hit.f] | (*fStripFlag)[hit.b])) continue; // if used
//L1_SHOW(fvHitKeyFlags.size());
@@ -320,15 +320,15 @@ void L1Algo::FindMoreHits(L1Branch& t, L1TrackPar& T, const bool dir,
L1UMeasurementInfo info = sta.frontInfo;
- if (fUseHitErrors) { info.sigma2 = hit.du * hit.du; }
+ if (fUseHitErrors) { info.sigma2 = hit.du2; }
L1Filter(T, info, u);
info = sta.backInfo;
- if (fUseHitErrors) { info.sigma2 = hit.dv * hit.dv; }
+ if (fUseHitErrors) { info.sigma2 = hit.dv2; }
L1Filter(T, info, v);
- FilterTime(T, hit.t, hit.dt);
+ FilterTime(T, hit.t, hit.dt2);
fldB0 = fldB1;
fldB1 = fldB2;
diff --git a/reco/L1/L1Algo/L1CATrackFinder.cxx b/reco/L1/L1Algo/L1CATrackFinder.cxx
index 5facd7a142..3ac5105e2f 100644
--- a/reco/L1/L1Algo/L1CATrackFinder.cxx
+++ b/reco/L1/L1Algo/L1CATrackFinder.cxx
@@ -128,8 +128,7 @@ bool L1Algo::checkTripletMatch(const L1Triplet& l, const L1Triplet& r, fscal& dc
inline void L1Algo::findSingletsStep0( // input
Tindex start_lh, Tindex n1_l, L1HitPoint* Hits_l,
// output
- fvec* u_front_l, fvec* u_back_l, fvec* zPos_l, L1HitIndex_t* hitsl, fvec* HitTime_l, fvec* HitTimeEr, fvec* d_u,
- fvec* d_v)
+ fvec* u_front_l, fvec* u_back_l, fvec* zPos_l, L1HitIndex_t* hitsl, fvec* t_l, fvec* dt2_l, fvec* du2_l, fvec* dv2_l)
{
/// Prepare the portion of data of left hits of a triplet:
@@ -141,12 +140,12 @@ inline void L1Algo::findSingletsStep0( // input
if (lastV >= 0) {
// set some positive errors to unfilled part of vectors in order to avoid nans
L1HitPoint& hitl = Hits_l[0];
- d_u[lastV] = hitl.dU();
- d_v[lastV] = hitl.dV();
- HitTimeEr[lastV] = hitl.timeEr;
+ du2_l[lastV] = hitl.dU2();
+ dv2_l[lastV] = hitl.dV2();
+ dt2_l[lastV] = hitl.dT2();
u_front_l[lastV] = hitl.U();
u_back_l[lastV] = hitl.V();
- HitTime_l[lastV] = hitl.time;
+ t_l[lastV] = hitl.T();
zPos_l[lastV] = hitl.Z();
}
@@ -156,16 +155,16 @@ inline void L1Algo::findSingletsStep0( // input
L1HitPoint& hitl = Hits_l[ilh];
- HitTime_l[i1_V][i1_4] = hitl.time;
- HitTimeEr[i1_V][i1_4] = hitl.timeEr;
+ t_l[i1_V][i1_4] = hitl.T();
+ dt2_l[i1_V][i1_4] = hitl.dT2();
hitsl[i1] = ilh;
u_front_l[i1_V][i1_4] = hitl.U();
u_back_l[i1_V][i1_4] = hitl.V();
if (fUseHitErrors) {
- d_u[i1_V][i1_4] = hitl.dU();
- d_v[i1_V][i1_4] = hitl.dV();
+ du2_l[i1_V][i1_4] = hitl.dU2();
+ dv2_l[i1_V][i1_4] = hitl.dV2();
}
zPos_l[i1_V][i1_4] = hitl.Z();
@@ -177,10 +176,10 @@ inline void L1Algo::findSingletsStep1( /// input 1st stage of singlet search
int istal,
int istam, /// indexes of left and middle stations of a triplet
Tindex n1_V, ///
- fvec* u_front_l, fvec* u_back_l, fvec* zPos_l, fvec* HitTime_l, fvec* HitTimeEr,
+ fvec* u_front_l, fvec* u_back_l, fvec* zPos_l, fvec* t_l, fvec* dt2_l,
// output
// L1TrackPar *T_1,
- L1TrackPar* T_1, L1FieldRegion* fld_1, fvec* d_u, fvec* d_v)
+ L1TrackPar* T_1, L1FieldRegion* fld_1, fvec* du2_l, fvec* dv2_l)
{
/// Get the field approximation. Add the target to parameters estimation.
@@ -244,8 +243,8 @@ inline void L1Algo::findSingletsStep1( /// input 1st stage of singlet search
fvec& v = u_back_l[i1_V];
auto [xl, yl] = stal.ConvUVtoXY<fvec>(u, v);
fvec zl = zPos_l[i1_V];
- fvec& time = HitTime_l[i1_V];
- fvec& timeEr = HitTimeEr[i1_V];
+ fvec& time = t_l[i1_V];
+ fvec& timeEr2 = dt2_l[i1_V];
const fvec dzli = 1.f / (zl - fTargZ);
const fvec tx = (xl - fTargX) * dzli;
@@ -283,7 +282,7 @@ inline void L1Algo::findSingletsStep1( /// input 1st stage of singlet search
T.C44 = fMaxInvMom / fvec(3.) * fMaxInvMom / fvec(3.);
- T.C55 = timeEr * timeEr;
+ T.C55 = timeEr2;
{ // add the target constraint
T.x = xl;
@@ -293,7 +292,7 @@ inline void L1Algo::findSingletsStep1( /// input 1st stage of singlet search
T.C10 = stal.XYInfo.C10;
T.C11 = stal.XYInfo.C11;
- if (fUseHitErrors) { std::tie(T.C00, T.C10, T.C11) = stal.FormXYCovarianceMatrix(d_u[i1_V], d_v[i1_V]); }
+ if (fUseHitErrors) { std::tie(T.C00, T.C10, T.C11) = stal.FormXYCovarianceMatrix(du2_l[i1_V], dv2_l[i1_V]); }
//assert(T.IsConsistent(true, -1));
@@ -406,14 +405,14 @@ inline void L1Algo::findDoubletsStep0(
// Pick_m22 is not used, search for mean squared, 2nd version
// -- collect possible doublets --
- const fscal iz = 1.f / (T1.z[i1_4] - fParameters.GetTargetPositionZ()[0]);
- const fscal timeError = T1.C55[i1_4];
- const fscal time = T1.t[i1_4];
+ const fscal iz = 1.f / (T1.z[i1_4] - fParameters.GetTargetPositionZ()[0]);
+ const fscal timeError2 = T1.C55[i1_4];
+ const fscal time = T1.t[i1_4];
L1HitAreaTime areaTime(vGridTime[iStaM], T1.x[i1_4] * iz, T1.y[i1_4] * iz,
(sqrt(Pick_m22 * (T1.C00 + stam.XYInfo.C00)) + fMaxDZ * abs(T1.tx))[i1_4] * iz,
(sqrt(Pick_m22 * (T1.C11 + stam.XYInfo.C11)) + fMaxDZ * abs(T1.ty))[i1_4] * iz, time,
- sqrt(timeError) * 5);
+ sqrt(timeError2) * 5);
for (L1HitIndex_t imh = -1; true;) { // loop over the hits in the area
if (fParameters.DevIsIgnoreHitSearchAreas()) {
@@ -435,8 +434,8 @@ inline void L1Algo::findDoubletsStep0(
// check y-boundaries
//TODO: move hardcoded cuts to parameters
if ((stam.timeInfo) && (stal.timeInfo)) {
- if (fabs(time - hitm.time) > sqrt(timeError + hitm.timeEr * hitm.timeEr) * 5) continue;
- if (fabs(time - hitm.time) > 30) continue;
+ if (fabs(time - hitm.T()) > sqrt(timeError2 + hitm.dT2()) * 5) continue;
+ if (fabs(time - hitm.T()) > 30) continue;
}
// - check whether hit belong to the window ( track position +\- errors ) -
@@ -458,7 +457,7 @@ inline void L1Algo::findDoubletsStep0(
fscal dy_est2 = Pick_m22[i1_4] * fabs(C11[i1_4] + stam.XYInfo.C11[i1_4]);
/// Covariation matrix of the hit
- auto [dxxScalMhit, dxyScalMhit, dyyScalMhit] = stam.FormXYCovarianceMatrix(hitm.dU(), hitm.dV());
+ auto [dxxScalMhit, dxyScalMhit, dyyScalMhit] = stam.FormXYCovarianceMatrix(hitm.dU2(), hitm.dV2());
if (fUseHitErrors) { dy_est2 = Pick_m22[i1_4] * fabs(C11[i1_4] + dyyScalMhit); }
@@ -487,7 +486,7 @@ inline void L1Algo::findDoubletsStep0(
L1UMeasurementInfo info = stam.frontInfo;
- if (fUseHitErrors) info.sigma2 = hitm.dU() * hitm.dU();
+ if (fUseHitErrors) info.sigma2 = hitm.dU2();
L1FilterChi2XYC00C10C11(info, x, y, C00, C10, C11, chi2, hitm.U());
@@ -498,11 +497,11 @@ inline void L1Algo::findDoubletsStep0(
info = stam.backInfo;
- if (fUseHitErrors) info.sigma2 = hitm.dV() * hitm.dV();
+ if (fUseHitErrors) info.sigma2 = hitm.dV2();
L1FilterChi2(info, x, y, C00, C10, C11, chi2, hitm.V());
- // FilterTime(T1, hitm.time, hitm.timeEr);
+ // FilterTime(T1, hitm.T(), hitm.dT2());
if (!fpCurrentIteration->GetTrackFromTripletsFlag()) {
if (chi2[i1_4] > fDoubletChi2Cut) continue;
@@ -543,7 +542,7 @@ inline void L1Algo::findTripletsStep0( // input
// 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>& dv_, L1Vector<fvec>& du_, L1Vector<fvec>& timeR, L1Vector<fvec>& timeER)
+ L1Vector<fvec>& dv2_3, L1Vector<fvec>& du2_3, L1Vector<fvec>& t_3, L1Vector<fvec>& dt2_3)
{
int iStaM = &stam - fParameters.GetStations().begin();
int iStaR = &star - fParameters.GetStations().begin();
@@ -572,10 +571,10 @@ inline void L1Algo::findTripletsStep0( // input
u_front_3.reset(1, fvec::Zero());
u_back_3.reset(1, fvec::Zero());
z_Pos_3.reset(1, fvec::Zero());
- du_.reset(1, fvec::One());
- dv_.reset(1, fvec::One());
- timeR.reset(1, fvec::Zero());
- timeER.reset(1, fvec::One());
+ du2_3.reset(1, fvec::One());
+ dv2_3.reset(1, fvec::One());
+ t_3.reset(1, fvec::Zero());
+ dt2_3.reset(1, fvec::One());
assert(istar < fParameters.GetNstationsActive()); //TODO SG!!! check if it is needed
@@ -589,11 +588,11 @@ inline void L1Algo::findTripletsStep0( // input
// pack the data
fvec u_front_2 = 0.f;
fvec u_back_2 = 0.f;
- fvec du2 = 1.f;
- fvec dv2 = 1.f;
+ fvec du2_2 = 1.f;
+ fvec dv2_2 = 1.f;
fvec zPos_2 = 0.f;
- fvec timeM = 0.f;
- fvec timeMEr = 1.f;
+ fvec t_2 = 0.f;
+ fvec dt2_2 = 1.f;
size_t n2_4 = 0;
for (; n2_4 < fvec::size() && i2 < n2; i2++, n2_4++) {
@@ -616,11 +615,11 @@ inline void L1Algo::findTripletsStep0( // input
u_front_2[n2_4] = hitm.U();
u_back_2[n2_4] = hitm.V();
zPos_2[n2_4] = hitm.Z();
- timeM[n2_4] = hitm.time;
- timeMEr[n2_4] = hitm.timeEr;
+ t_2[n2_4] = hitm.T();
+ dt2_2[n2_4] = hitm.dT2();
// num[n2_4] = hitm.track;
- du2[n2_4] = hitm.dU();
- dv2[n2_4] = hitm.dV();
+ du2_2[n2_4] = hitm.dU2();
+ dv2_2[n2_4] = hitm.dV2();
hitsl_2[n2_4] = hitsl_1[i1];
hitsm_2_tmp[n2_4] = hitsm_2[i2];
@@ -641,7 +640,7 @@ inline void L1Algo::findTripletsStep0( // input
L1UMeasurementInfo info = stam.frontInfo;
- if (fUseHitErrors) info.sigma2 = du2 * du2;
+ if (fUseHitErrors) info.sigma2 = du2_2;
// TODO: SG: L1FilterNoField is wrong.
// TODO: If the field was present before,
@@ -677,7 +676,7 @@ inline void L1Algo::findTripletsStep0( // input
*/
info = stam.backInfo;
- if (fUseHitErrors) info.sigma2 = dv2 * dv2;
+ if (fUseHitErrors) info.sigma2 = dv2_2;
if (istam < fNfieldStations) { L1Filter(T2, info, u_back_2); }
else {
@@ -686,7 +685,7 @@ inline void L1Algo::findTripletsStep0( // input
// assert(T2.IsConsistent(true, n2_4));
- FilterTime(T2, timeM, timeMEr, stam.timeInfo);
+ FilterTime(T2, t_2, dt2_2, stam.timeInfo);
// assert(T2.IsConsistent(true, n2_4));
@@ -735,9 +734,9 @@ inline void L1Algo::findTripletsStep0( // input
if (fabs(T2.tx[i2_4]) > fMaxSlope) continue;
if (fabs(T2.ty[i2_4]) > fMaxSlope) continue;
- const fvec Pick_r22 = fTripletChi2Cut - T2.chi2 + (!fpCurrentIteration->GetTrackFromTripletsFlag() ? 0 : 1);
- const fscal timeError = T2.C55[i2_4];
- const fscal time = T2.t[i2_4];
+ const fvec Pick_r22 = fTripletChi2Cut - T2.chi2 + (!fpCurrentIteration->GetTrackFromTripletsFlag() ? 0 : 1);
+ const fscal timeError2 = T2.C55[i2_4];
+ const fscal time = T2.t[i2_4];
// find first possible hit
@@ -745,7 +744,7 @@ inline void L1Algo::findTripletsStep0( // input
L1HitAreaTime area(vGridTime[&star - fParameters.GetStations().begin()], T2.x[i2_4] * iz, T2.y[i2_4] * iz,
(sqrt(Pick_r22 * (T2.C00 + stam.XYInfo.C00)) + fMaxDZ * abs(T2.tx))[i2_4] * iz,
(sqrt(Pick_r22 * (T2.C11 + stam.XYInfo.C11)) + fMaxDZ * abs(T2.ty))[i2_4] * iz, time,
- sqrt(timeError) * 5);
+ sqrt(timeError2) * 5);
L1HitIndex_t irh = 0;
L1HitIndex_t doubletNtriplets = 0;
@@ -783,11 +782,11 @@ inline void L1Algo::findTripletsStep0( // input
L1ExtrapolateLine(T_cur, zr);
if ((star.timeInfo) && (stam.timeInfo))
- if (fabs(T_cur.t[i2_4] - hitr.time) > sqrt(T_cur.C55[i2_4] + hitr.timeEr) * 5) continue;
+ if (fabs(T_cur.t[i2_4] - hitr.T()) > sqrt(T_cur.C55[i2_4] + sqrt(hitr.dT2())) * 5) continue;
// TODO: SG: hardcoded cut of 30 ns
if ((star.timeInfo) && (stam.timeInfo))
- if (fabs(T_cur.t[i2_4] - hitr.time) > 30) continue;
+ if (fabs(T_cur.t[i2_4] - hitr.T()) > 30) continue;
// - check whether hit belong to the window ( track position +\- errors ) -
// check lower boundary
@@ -801,7 +800,7 @@ inline void L1Algo::findTripletsStep0( // input
+ star.XYInfo.C11[i2_4]))); // TODO for FastPrim dx < dy - other sort is optimal. But not for doublets
/// Covariation matrix of the hit
- auto [dxxScalRhit, dxyScalRhit, dyyScalRhit] = star.FormXYCovarianceMatrix(hitr.dU(), hitr.dV());
+ auto [dxxScalRhit, dxyScalRhit, dyyScalRhit] = star.FormXYCovarianceMatrix(hitr.dU2(), hitr.dV2());
if (fUseHitErrors) { dy_est2 = (Pick_r22[i2_4] * (fabs(C11[i2_4] + dyyScalRhit))); }
@@ -829,16 +828,16 @@ inline void L1Algo::findTripletsStep0( // input
info = star.frontInfo;
- if (fUseHitErrors) info.sigma2 = hitr.dU() * hitr.dU();
+ if (fUseHitErrors) info.sigma2 = hitr.dU2();
L1FilterChi2XYC00C10C11(info, x, y, C00, C10, C11, chi2, hitr.U());
info = star.backInfo;
- if (fUseHitErrors) info.sigma2 = hitr.dV() * hitr.dV();
+ if (fUseHitErrors) info.sigma2 = hitr.dV2();
L1FilterChi2(info, x, y, C00, C10, C11, chi2, hitr.V());
- FilterTime(T_cur, hitr.time, hitr.timeEr, star.timeInfo);
+ FilterTime(T_cur, hitr.T(), hitr.dT2(), star.timeInfo);
if (!fpCurrentIteration->GetTrackFromTripletsFlag()) {
@@ -877,11 +876,11 @@ inline void L1Algo::findTripletsStep0( // input
T3.SetOneEntry(n3_4, T2, i2_4);
u_front_3[n3_V][n3_4] = hitr.U();
u_back_3[n3_V][n3_4] = hitr.V();
- du_[n3_V][n3_4] = hitr.dU();
- dv_[n3_V][n3_4] = hitr.dV();
+ du2_3[n3_V][n3_4] = hitr.dU2();
+ dv2_3[n3_V][n3_4] = hitr.dV2();
z_Pos_3[n3_V][n3_4] = zr;
- timeR[n3_V][n3_4] = hitr.time;
- timeER[n3_V][n3_4] = hitr.timeEr;
+ t_3[n3_V][n3_4] = hitr.T();
+ dt2_3[n3_V][n3_4] = hitr.dT2();
n3++;
n3_V = n3 / fvec::size();
@@ -892,10 +891,10 @@ inline void L1Algo::findTripletsStep0( // input
u_front_3.push_back(fvec::Zero());
u_back_3.push_back(fvec::Zero());
z_Pos_3.push_back(fvec::Zero());
- du_.push_back(fvec::Zero());
- dv_.push_back(fvec::Zero());
- timeR.push_back(fvec::Zero());
- timeER.push_back(fvec::Zero());
+ du2_3.push_back(fvec::One());
+ dv2_3.push_back(fvec::One());
+ t_3.push_back(fvec::Zero());
+ dt2_3.push_back(fvec::One());
}
} // search area
@@ -906,7 +905,7 @@ inline void L1Algo::findTripletsStep0( // input
/// Add the right hits to parameters estimation.
inline void L1Algo::findTripletsStep1( // input
Tindex n3_V, const L1Station& star, L1Vector<fvec>& u_front_, L1Vector<fvec>& u_back_, L1Vector<fvec>& z_Pos,
- L1Vector<fvec>& dv_, L1Vector<fvec>& du_, L1Vector<fvec>& timeR, L1Vector<fvec>& timeER,
+ L1Vector<fvec>& dv2_3, L1Vector<fvec>& du2_3, L1Vector<fvec>& t_3, L1Vector<fvec>& dt2_3,
// output
// L1TrackPar *T_3
L1Vector<L1TrackPar>& T_3)
@@ -928,7 +927,7 @@ inline void L1Algo::findTripletsStep1( // input
L1UMeasurementInfo info = star.frontInfo;
- if (fUseHitErrors) info.sigma2 = du_[i3_V] * du_[i3_V];
+ if (fUseHitErrors) info.sigma2 = du2_3[i3_V];
bool noField = (&star - fParameters.GetStations().begin() >= fNfieldStations);
@@ -941,7 +940,7 @@ inline void L1Algo::findTripletsStep1( // input
info = star.backInfo;
- if (fUseHitErrors) info.sigma2 = dv_[i3_V] * dv_[i3_V];
+ if (fUseHitErrors) info.sigma2 = dv2_3[i3_V];
if (noField) { L1FilterNoField(T3, info, u_back_[i3_V]); }
else {
@@ -950,9 +949,7 @@ inline void L1Algo::findTripletsStep1( // input
// assert(T3.IsConsistent(true, -1));
- if (kGlobal != fTrackingMode && kMcbm != fTrackingMode) {
- FilterTime(T3, timeR[i3_V], timeER[i3_V], star.timeInfo);
- }
+ if (kGlobal != fTrackingMode && kMcbm != fTrackingMode) { FilterTime(T3, t_3[i3_V], dt2_3[i3_V], star.timeInfo); }
// assert(T3.IsConsistent(true, -1));
// FilterTime(T_3[i3_V], timeR[i3_V], timeER[i3_V]);
@@ -1018,7 +1015,7 @@ inline void L1Algo::findTripletsStep2(Tindex n3, int istal, int istam, int istar
}
fscal u[NHits], v[NHits], t[NHits], x[NHits], y[NHits], z[NHits];
- fscal du[NHits], dv[NHits], dt[NHits], du2[NHits], dv2[NHits], dt2[NHits];
+ fscal du2[NHits], dv2[NHits], dt2[NHits];
for (int ih = 0; ih < NHits; ++ih) {
const L1Hit& hit = fInputData.GetHit(ihit[ih]);
@@ -1027,13 +1024,10 @@ inline void L1Algo::findTripletsStep2(Tindex n3, int istal, int istam, int istar
t[ih] = hit.t;
std::tie(x[ih], y[ih]) = sta[ih].ConvUVtoXY<fscal>(u[ih], v[ih]);
z[ih] = hit.z;
- du[ih] = hit.du;
- dv[ih] = hit.dv;
- dt[ih] = hit.dt;
- du2[ih] = hit.du * hit.du;
- dv2[ih] = hit.dv * hit.dv;
- dt2[ih] = hit.dt * hit.dt;
+ du2[ih] = hit.du2;
+ dv2[ih] = hit.dv2;
+ dt2[ih] = hit.dt2;
};
// find the field along the track
@@ -1080,11 +1074,11 @@ inline void L1Algo::findTripletsStep2(Tindex n3, int istal, int istam, int istar
T.z = z[ih0];
T.t = t[ih0];
- //std::tie(T.C00, T.C10, T.C11) = sta[ih0].FormXYCovarianceMatrix(du[ih0], dv[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], dt[ih0], fvec::One(), sta[ih0].timeInfo);
+ fit.FilterTime(t[ih0], dt2[ih0], fvec::One(), sta[ih0].timeInfo);
{ // add the target constraint
fvec eX, eY, J04, J14;
@@ -1112,7 +1106,7 @@ inline void L1Algo::findTripletsStep2(Tindex n3, int istal, int istam, int istar
//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());
- fit.FilterTime(t[ih], dt[ih], fvec::One(), sta[ih].timeInfo);
+ fit.FilterTime(t[ih], dt2[ih], fvec::One(), sta[ih].timeInfo);
}
}
@@ -1134,7 +1128,7 @@ inline void L1Algo::findTripletsStep2(Tindex n3, int istal, int istam, int istar
T.t = t[ih0];
T.C55 = dt2[ih0];
- //std::tie(T.C00, T.C10, T.C11) = sta[ih0].FormXYCovarianceMatrix(du[ih0], dv[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], dt[ih0], fvec::One(), sta[ih0].timeInfo);
@@ -1422,18 +1416,18 @@ inline void L1Algo::CreatePortionOfDoublets(
fvec u_front[L1Constants::size::kSingletPortionSizeVec];
fvec u_back[L1Constants::size::kSingletPortionSizeVec];
- fvec dv0[L1Constants::size::kSingletPortionSizeVec];
- fvec du0[L1Constants::size::kSingletPortionSizeVec];
+ fvec dv2[L1Constants::size::kSingletPortionSizeVec];
+ fvec du2[L1Constants::size::kSingletPortionSizeVec];
fvec zPos[L1Constants::size::kSingletPortionSizeVec];
- fvec HitTime[L1Constants::size::kSingletPortionSizeVec];
- fvec HitTimeEr[L1Constants::size::kSingletPortionSizeVec];
+ fvec t[L1Constants::size::kSingletPortionSizeVec];
+ fvec dt2[L1Constants::size::kSingletPortionSizeVec];
/// prepare the portion of left hits data
findSingletsStep0( // input
iSingletPortion * L1Constants::size::kSingletPortionSize, singletPortionSize, vHits_l,
// output
- u_front, u_back, zPos, hitsl_1, HitTime, HitTimeEr, du0, dv0);
+ u_front, u_back, zPos, hitsl_1, t, dt2, du2, dv2);
for (Tindex i = 0; i < singletPortionSize; ++i)
L1_ASSERT(hitsl_1[i] < HitsUnusedStopIndex[istal] - HitsUnusedStartIndex[istal],
@@ -1443,9 +1437,9 @@ inline void L1Algo::CreatePortionOfDoublets(
/// Get the field approximation. Add the target to parameters estimation. Propagaete to middle station.
- findSingletsStep1(istal, istam, portionSize_V, u_front, u_back, zPos, HitTime, HitTimeEr,
+ findSingletsStep1(istal, istam, portionSize_V, u_front, u_back, zPos, t, dt2,
// output
- T_1, fld_1, du0, dv0);
+ T_1, fld_1, du2, dv2);
/// Find the doublets. Reformat data in the portion of doublets.
@@ -1540,10 +1534,10 @@ inline void L1Algo::CreatePortionOfTriplets(
L1Vector<fvec>& u_front3 = fTripletHitR_Ufront[Thread];
L1Vector<fvec>& u_back3 = fTripletHitR_Uback[Thread];
L1Vector<fvec>& z_pos3 = fTripletHitR_Z[Thread];
- L1Vector<fvec>& timeR = fTripletHitR_Time[Thread];
- L1Vector<fvec>& timeER = fTripletHitR_TimeErr[Thread];
- L1Vector<fvec>& du3 = fTripletHitR_dUfront[Thread];
- L1Vector<fvec>& dv3 = fTripletHitR_dUback[Thread];
+ L1Vector<fvec>& t_3 = fTripletHitR_Time[Thread];
+ L1Vector<fvec>& dt2_3 = fTripletHitR_TimeErr[Thread];
+ L1Vector<fvec>& du2_3 = fTripletHitR_dUfront[Thread];
+ L1Vector<fvec>& dv2_3 = fTripletHitR_dUback[Thread];
T_3.clear();
hitsl_3.clear();
@@ -1552,10 +1546,10 @@ inline void L1Algo::CreatePortionOfTriplets(
u_front3.clear();
u_back3.clear();
z_pos3.clear();
- du3.clear();
- dv3.clear();
- timeR.clear();
- timeER.clear();
+ du2_3.clear();
+ dv2_3.clear();
+ t_3.clear();
+ dt2_3.clear();
/// Find the triplets(right hit). Reformat data in the portion of triplets.
@@ -1567,7 +1561,7 @@ inline void L1Algo::CreatePortionOfTriplets(
n_2, hitsm_2, i1_2,
// output
- n3, T_3, hitsl_3, hitsm_3, hitsr_3, u_front3, u_back3, z_pos3, du3, dv3, timeR, timeER);
+ n3, T_3, hitsl_3, hitsm_3, hitsr_3, u_front3, u_back3, z_pos3, du2_3, dv2_3, t_3, dt2_3);
for (Tindex i = 0; i < static_cast<Tindex>(hitsl_3.size()); ++i)
@@ -1582,7 +1576,7 @@ inline void L1Algo::CreatePortionOfTriplets(
Tindex n3_V = (n3 + fvec::size() - 1) / fvec::size();
findTripletsStep1( // input
- n3_V, star, u_front3, u_back3, z_pos3, du3, dv3, timeR, timeER,
+ n3_V, star, u_front3, u_back3, z_pos3, du2_3, dv2_3, t_3, dt2_3,
// output
T_3);
diff --git a/reco/L1/L1Algo/L1Filtration.h b/reco/L1/L1Algo/L1Filtration.h
index f9975ef476..9797c173d8 100644
--- a/reco/L1/L1Algo/L1Filtration.h
+++ b/reco/L1/L1Algo/L1Filtration.h
@@ -14,7 +14,7 @@
#define cnst const fvec
-inline void FilterTime(L1TrackPar& T, fvec t, fvec dt, fvec timeInfo = fvec::One(), fvec w = fvec::One())
+inline void FilterTime(L1TrackPar& T, fvec t, fvec dt2, fvec timeInfo = fvec::One(), fvec w = fvec::One())
{
// filter track with a time measurement
@@ -30,8 +30,6 @@ inline void FilterTime(L1TrackPar& T, fvec t, fvec dt, fvec timeInfo = fvec::One
w.setZero(timeInfo <= fvec::Zero());
- fvec dt2 = dt * dt;
-
// when dt0 is much smaller than current time error,
// set track time exactly to the measurement value without filtering
// it helps to keep the initial time errors reasonably small
diff --git a/reco/L1/L1Algo/L1Hit.h b/reco/L1/L1Algo/L1Hit.h
index acdea6e3c4..983b0462a9 100644
--- a/reco/L1/L1Algo/L1Hit.h
+++ b/reco/L1/L1Algo/L1Hit.h
@@ -38,9 +38,9 @@ public:
float v = 0.f; ///< measured V coordinate [cm]
float t = 0.f; ///< measured time [ns]
float z = 0.f; ///< fixed Z coordinate [cm]
- float du = 0.f; ///< measured uncertainty of U coordinate [cm]
- float dv = 0.f; ///< measured uncertainty of V coordinate [cm]
- float dt = 0.f; ///< measured uncertainty of time [ns]
+ float du2 = 0.f; ///< measured uncertainty of U coordinate [cm]
+ float dv2 = 0.f; ///< measured uncertainty of V coordinate [cm]
+ float dt2 = 0.f; ///< measured uncertainty of time [ns]
int ID = 0; ///< index of hit before hits sorting
int iSt = 0; ///< index of station in the active stations array
// TODO: Test speed penalty of using iSt index
@@ -56,9 +56,9 @@ private:
ar& v;
ar& t;
ar& z;
- ar& du;
- ar& dv;
- ar& dt;
+ ar& du2;
+ ar& dv2;
+ ar& dt2;
ar& ID;
ar& iSt;
}
diff --git a/reco/L1/L1Algo/L1HitPoint.h b/reco/L1/L1Algo/L1HitPoint.h
index ceb27e5c91..5fd6b8111f 100644
--- a/reco/L1/L1Algo/L1HitPoint.h
+++ b/reco/L1/L1Algo/L1HitPoint.h
@@ -5,140 +5,50 @@
#ifndef _L1HitPoint_h_
#define _L1HitPoint_h_
-/// contain strips positions and coordinates of hit
-#if 1
+/// contain strips positions and coordinates of the hit
+
struct L1HitPoint {
- L1HitPoint()
- : // x(0.f)
- // , y(0.f)
- // , dx(0.f)
- // , dy(0.f)
- // , dxy(0.f)
- u(0.f)
- , v(0.f)
- , du(0.f)
- , dv(0.f)
- , z(0.f)
- , time(0.f)
- , timeEr(2.9f) {};
- L1HitPoint(
- // fscal dx_,
- // fscal dy_,
- // fscal dxy_,
- fscal du_, fscal dv_, fscal z_, fscal v_, fscal u_, float time_, float /*timeEv1_*/ = 0,
- float timeEr_ = 2.9f)
- : // : x(x_)
- // , y(y_)
- // , dx(dx_)
- // , dy(dy_)
- // , dxy(dxy_)
- u(u_)
- , v(v_)
- , du(du_)
- , dv(dv_)
- , z(z_)
- , time(time_)
- , timeEr(timeEr_) {};
- // L1HitPoint(fscal x_, fscal y_, fscal z_, fscal v_, fscal u_, fscal time_, unsigned short int n_ = 0):
- // x(x_), y(y_), z(z_), u(u_), v(v_), time(time_){};
+ L1HitPoint() = default;
- // fscal Xs() const { return X() / Z(); }
- // fscal Ys() const { return Y() / Z(); } // value to sort hits by
+ L1HitPoint(fscal z_, fscal v_, fscal u_, fscal t_, fscal du2_, fscal dv2_, fscal dt2_)
+ : z(z_)
+ , u(u_)
+ , v(v_)
+ , t(t_)
+ , du2(du2_)
+ , dv2(dv2_)
+ , dt2(dt2_) {};
- // fscal X() const { return x; }
- // fscal Y() const { return y; }
- // fscal dX() const { return dx; }
- // fscal dY() const { return dy; }
- // fscal dXY() const { return dxy; }
- fscal dU() const { return du; }
- fscal dV() const { return dv; }
fscal Z() const { return z; }
fscal U() const { return u; }
fscal V() const { return v; }
- // unsigned short int N() const { return n; }
- // int GetSortIndex() const { return SortIndex; }
- // void SetX(fscal X1) { x = X1; }
- // void SetY(fscal Y1) { y = Y1; }
- void SetZ(fscal Z1) { z = Z1; }
- void SetU(fscal U1) { u = U1; }
- void SetV(fscal V1) { v = V1; }
-
- void Set(
- // const float& dx1,
- // const float& dy1,
- // const float& xy1,
- const float& z1, const fscal& u1, const fscal& v1, const fscal& du1, const fscal& dv1, const float& time1,
- float timeEr1)
- {
- // x = x1;
- // y = y1;
- // dx = dx1;
- // dy = dy1;
- du = du1;
- // dxy = xy1;
- dv = dv1;
- z = z1;
- u = u1;
- v = v1;
- time = time1;
- timeEr = timeEr1;
- }
-
-private:
- // x, y,
- float u, v, du, dv;
- float z; // TODO: may be we should use iz
- // x\u, v - front and back strips; x, y, z - coordinates of hits
-
-public:
- float time, timeEr;
- // int track;
-};
-#else
-
-static const float R = 60;
-static const float shortPackingConstant = 2 * R / 65535.f;
-//TODO: change the Z conversion. Hit Z maybe negative in the new setup.
-static const float MZ = 110;
-static const float shortPackingConstantZ = MZ / 65535.f;
+ fscal T() const { return t; }
-/// contain strips positions and coordinates of hit
-struct L1HitPoint {
- L1HitPoint() {};
- L1HitPoint(fscal x_, fscal y_, fscal z_, fscal v_, fscal u_, unsigned short int n_ = 0)
- : x(f2s(x_))
- , y(f2s(y_))
- , z(f2sZ(z_))
- , u(f2s(u_))
- , v(f2s(v_))
- , n(n_) {};
+ fscal dU2() const { return du2; }
+ fscal dV2() const { return dv2; }
+ fscal dT2() const { return dt2; }
- fscal Xs() const { return X() / Z(); }
- fscal Ys() const { return Y() / Z(); } // value to sort hits by
+ void SetZ(fscal z_) { z = z_; }
+ void SetU(fscal u_) { u = u_; }
+ void SetV(fscal v_) { v = v_; }
+ void SetT(fscal t_) { t = t_; }
- fscal X() const { return s2f(x); }
- fscal Y() const { return s2f(y); }
- fscal Z() const { return s2fZ(z); }
- fscal U() const { return s2f(u); }
- fscal V() const { return s2f(v); }
-
- unsigned short int N() const { return n; }
+ void Set(const fscal& z_, const fscal& u_, const fscal& v_, const fscal& t_, const fscal& du2_, const fscal& dv2_,
+ const fscal& dt2_)
+ {
+ z = z_;
+ u = u_;
+ v = v_;
+ t = t_;
+ du2 = du2_;
+ dv2 = dv2_;
+ dt2 = dt2_;
+ }
private:
- //unsigned short int
- unsigned short int f2s(float f) const { return (f + R) / shortPackingConstant; }
- float s2f(unsigned short int f) const { return (float(f) + 0.5) * shortPackingConstant - R; }
-
- unsigned short int f2sZ(float f) const { return (f) / shortPackingConstantZ; }
- float s2fZ(unsigned short int f) const { return (float(f) + 0.5) * shortPackingConstantZ; }
-
- unsigned short int x, y;
- unsigned short int z; // TODO: may be we should use iz
- unsigned short int u,
- v; // x\u, v - front and back strips; x, y, z - coordinates of hits
+ fscal z {0.}, u {0.}, v {0.}, t {0.}, du2 {0.}, dv2 {0.}, dt2 {0.};
};
-#endif
#endif
diff --git a/reco/L1/L1Algo/L1Station.h b/reco/L1/L1Algo/L1Station.h
index 7e176b6ef0..c0103b1b7f 100644
--- a/reco/L1/L1Algo/L1Station.h
+++ b/reco/L1/L1Algo/L1Station.h
@@ -95,7 +95,7 @@ public:
/// \param dv V coordinate uncertainty [length unit]
/// \return Covariance matrix of hit position in Cartesian coordinates: [dxx, dxy, dyy] [length unit squared]
template<typename T, std::enable_if_t<std::is_floating_point<T>::value || std::is_same<T, fvec>::value, bool> = true>
- std::tuple<T, T, T> FormXYCovarianceMatrix(T du, T dv) const;
+ std::tuple<T, T, T> FormXYCovarianceMatrix(T du2, T dv2) const;
} _fvecalignment;
@@ -133,20 +133,19 @@ std::pair<T, T> L1Station::ConvUVtoXY(T u, T v) const
// ---------------------------------------------------------------------------------------------------------------------
//
template<typename T, std::enable_if_t<std::is_floating_point<T>::value || std::is_same<T, fvec>::value, bool>>
-std::tuple<T, T, T> L1Station::FormXYCovarianceMatrix(T du, T dv) const
+std::tuple<T, T, T> L1Station::FormXYCovarianceMatrix(T du2, T dv2) const
{
if constexpr (std::is_same<T, fvec>::value) {
- return std::make_tuple(
- (xInfo.cos_phi * du) * (xInfo.cos_phi * du) + (xInfo.sin_phi * dv) * (xInfo.sin_phi * dv), // dx2
- (xInfo.cos_phi * du) * (yInfo.cos_phi * du) + (xInfo.sin_phi * dv) * (yInfo.sin_phi * dv), // dxy
- (yInfo.cos_phi * du) * (yInfo.cos_phi * du) + (yInfo.sin_phi * dv) * (yInfo.sin_phi * dv) // dy2
+ return std::make_tuple(xInfo.cos_phi * xInfo.cos_phi * du2 + xInfo.sin_phi * xInfo.sin_phi * dv2, // dx2
+ xInfo.cos_phi * yInfo.cos_phi * du2 + xInfo.sin_phi * yInfo.sin_phi * dv2, // dxy
+ yInfo.cos_phi * yInfo.cos_phi * du2 + yInfo.sin_phi * yInfo.sin_phi * dv2 // dy2
);
}
else {
return std::make_tuple(
- (xInfo.cos_phi[0] * du) * (xInfo.cos_phi[0] * du) + (xInfo.sin_phi[0] * dv) * (xInfo.sin_phi[0] * dv), // dx2
- (xInfo.cos_phi[0] * du) * (yInfo.cos_phi[0] * du) + (xInfo.sin_phi[0] * dv) * (yInfo.sin_phi[0] * dv), // dxy
- (yInfo.cos_phi[0] * du) * (yInfo.cos_phi[0] * du) + (yInfo.sin_phi[0] * dv) * (yInfo.sin_phi[0] * dv) // dy2
+ xInfo.cos_phi[0] * xInfo.cos_phi[0] * du2 + xInfo.sin_phi[0] * xInfo.sin_phi[0] * dv2, // dx2
+ xInfo.cos_phi[0] * yInfo.cos_phi[0] * du2 + xInfo.sin_phi[0] * yInfo.sin_phi[0] * dv2, // dxy
+ yInfo.cos_phi[0] * yInfo.cos_phi[0] * du2 + yInfo.sin_phi[0] * yInfo.sin_phi[0] * dv2 // dy2
);
}
}
diff --git a/reco/L1/L1Algo/L1TrackFitter.cxx b/reco/L1/L1Algo/L1TrackFitter.cxx
index 8a9a07ccf3..cfeb883c74 100644
--- a/reco/L1/L1Algo/L1TrackFitter.cxx
+++ b/reco/L1/L1Algo/L1TrackFitter.cxx
@@ -347,8 +347,8 @@ void L1Algo::L1KFTrackFitter()
// NOTE: u- and v-axes are axes, orthogonal to front and back strips of the station, respectively.
fvec u[L1Constants::size::kMaxNstations]; // Hit position along the u-axis [cm]
fvec v[L1Constants::size::kMaxNstations]; // Hit position along the v-axis [cm]
- fvec d_u[L1Constants::size::kMaxNstations]; // Hit position uncertainty along the u-axis [cm]
- fvec d_v[L1Constants::size::kMaxNstations]; // Hit position uncertainty along the v-axis [cm]
+ fvec du2[L1Constants::size::kMaxNstations]; // Hit position uncertainty along the u-axis [cm]
+ fvec dv2[L1Constants::size::kMaxNstations]; // Hit position uncertainty along the v-axis [cm]
fvec x[L1Constants::size::kMaxNstations]; // Hit position along the x-axis [cm]
fvec y[L1Constants::size::kMaxNstations]; // Hit position along the y-axis [cm]
@@ -358,8 +358,8 @@ void L1Algo::L1KFTrackFitter()
fvec z[L1Constants::size::kMaxNstations]; // Hit position along the z-axis (precised) [cm]
- fvec time[L1Constants::size::kMaxNstations]; // Hit time [ns]
- fvec timeEr[L1Constants::size::kMaxNstations]; // Hit time uncertainty [ns]
+ fvec time[L1Constants::size::kMaxNstations]; // Hit time [ns]
+ fvec dt2[L1Constants::size::kMaxNstations]; // Hit time uncertainty [ns] squared
fvec x_first;
fvec y_first;
@@ -368,7 +368,7 @@ void L1Algo::L1KFTrackFitter()
fvec d_xy_fst;
fvec time_first;
- fvec time_er_first;
+ fvec dt2_first;
fvec x_last;
fvec y_last;
@@ -377,8 +377,8 @@ void L1Algo::L1KFTrackFitter()
fvec d_xy_lst;
fvec time_last;
- fvec time_er_last;
- // fvec time_er_lst; /// TODO: Why are there two different variables for the time error on the last station?
+ fvec dt2_last;
+ // 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];
@@ -436,16 +436,16 @@ void L1Algo::L1KFTrackFitter()
u[ista][iVec] = hit.u;
v[ista][iVec] = hit.v;
- d_u[ista][iVec] = hit.du;
- d_v[ista][iVec] = hit.dv;
+ du2[ista][iVec] = hit.du2;
+ dv2[ista][iVec] = hit.dv2;
std::tie(x_temp, y_temp) = sta[ista].ConvUVtoXY<fvec>(u[ista], v[ista]);
x[ista][iVec] = x_temp[iVec];
y[ista][iVec] = y_temp[iVec];
time[ista][iVec] = hit.t;
- timeEr[ista][iVec] = hit.dt;
+ dt2[ista][iVec] = hit.dt2;
z[ista][iVec] = hit.z;
sta[ista].fieldSlice.GetFieldValue(x[ista], y[ista], fB_temp);
- std::tie(d_xx[ista], d_xy[ista], d_yy[ista]) = sta[ista].FormXYCovarianceMatrix(d_u[ista], d_v[ista]);
+ std::tie(d_xx[ista], d_xy[ista], d_yy[ista]) = sta[ista].FormXYCovarianceMatrix(du2[ista], dv2[ista]);
fB[ista].x[iVec] = fB_temp.x[iVec];
fB[ista].y[iVec] = fB_temp.y[iVec];
@@ -455,7 +455,7 @@ void L1Algo::L1KFTrackFitter()
x_first[iVec] = x[ista][iVec];
y_first[iVec] = y[ista][iVec];
time_first[iVec] = time[ista][iVec];
- time_er_first[iVec] = timeEr[ista][iVec];
+ dt2_first[iVec] = dt2[ista][iVec];
d_xx_fst[iVec] = d_xx[ista][iVec];
d_yy_fst[iVec] = d_yy[ista][iVec];
d_xy_fst[iVec] = d_xy[ista][iVec];
@@ -471,7 +471,7 @@ void L1Algo::L1KFTrackFitter()
d_yy_lst[iVec] = d_yy[ista][iVec];
d_xy_lst[iVec] = d_xy[ista][iVec];
time_last[iVec] = time[ista][iVec];
- time_er_last[iVec] = timeEr[ista][iVec];
+ dt2_last[iVec] = dt2[ista][iVec];
staLast.XYInfo.C00[iVec] = sta[ista].XYInfo.C00[iVec];
staLast.XYInfo.C10[iVec] = sta[ista].XYInfo.C10[iVec];
staLast.XYInfo.C11[iVec] = sta[ista].XYInfo.C11[iVec];
@@ -505,10 +505,10 @@ void L1Algo::L1KFTrackFitter()
int ista = nStations - 1;
- time_last = iif(w_time[ista] > fvec::Zero(), time_last, fvec::Zero());
- time_er_last = iif(w_time[ista] > fvec::Zero(), time_er_last, fvec(100.));
+ time_last = iif(w_time[ista] > fvec::Zero(), time_last, fvec::Zero());
+ dt2_last = iif(w_time[ista] > fvec::Zero(), dt2_last, fvec(100.));
- FilterFirst(fit, x_last, y_last, time_last, time_er_last, staLast, d_xx_lst, d_yy_lst, d_xy_lst);
+ FilterFirst(fit, x_last, y_last, time_last, dt2_last, staLast, d_xx_lst, d_yy_lst, d_xy_lst);
fldZ1 = z[ista];
@@ -551,9 +551,9 @@ void L1Algo::L1KFTrackFitter()
fit.AddMaterial(sta[ista].materialInfo, qp01, wExtr);
}
- fit.Filter(sta[ista].frontInfo, u[ista], d_u[ista] * d_u[ista], w1);
- fit.Filter(sta[ista].backInfo, v[ista], d_v[ista] * d_v[ista], w1);
- fit.FilterTime(time[ista], timeEr[ista], w1_time, sta[ista].timeInfo);
+ 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);
fldB2 = fldB1;
fldZ2 = fldZ1;
@@ -668,7 +668,7 @@ void L1Algo::L1KFTrackFitter()
ista = 0;
- FilterFirst(fit, x_first, y_first, time_first, time_er_first, staFirst, d_xx_fst, d_yy_fst, d_xy_fst);
+ FilterFirst(fit, x_first, y_first, time_first, dt2_first, staFirst, d_xx_fst, d_yy_fst, d_xy_fst);
qp01 = tr.qp;
@@ -708,9 +708,9 @@ void L1Algo::L1KFTrackFitter()
fit.AddMaterial(sta[ista].materialInfo, qp01, wExtr);
}
- fit.Filter(sta[ista].frontInfo, u[ista], d_u[ista] * d_u[ista], w1);
- fit.Filter(sta[ista].backInfo, v[ista], d_v[ista] * d_v[ista], w1);
- fit.FilterTime(time[ista], timeEr[ista], w1_time, sta[ista].timeInfo);
+ 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);
fldB2 = fldB1;
fldZ2 = fldZ1;
@@ -791,8 +791,8 @@ void L1Algo::L1KFTrackFitterMuch()
fvec u[L1Constants::size::kMaxNstations];
fvec v[L1Constants::size::kMaxNstations];
- fvec d_u[L1Constants::size::kMaxNstations];
- fvec d_v[L1Constants::size::kMaxNstations];
+ fvec du2[L1Constants::size::kMaxNstations];
+ fvec dv2[L1Constants::size::kMaxNstations];
fvec x[L1Constants::size::kMaxNstations];
fvec y[L1Constants::size::kMaxNstations];
@@ -803,7 +803,7 @@ void L1Algo::L1KFTrackFitterMuch()
fvec z[L1Constants::size::kMaxNstations];
fvec time[L1Constants::size::kMaxNstations];
- fvec timeEr[L1Constants::size::kMaxNstations];
+ fvec dt2[L1Constants::size::kMaxNstations];
fvec x_first;
fvec y_first;
@@ -812,7 +812,7 @@ void L1Algo::L1KFTrackFitterMuch()
fvec d_xy_fst;
fvec time_first;
- fvec time_er_fst;
+ fvec dt2_fst;
fvec x_last;
fvec y_last;
@@ -821,7 +821,7 @@ void L1Algo::L1KFTrackFitterMuch()
fvec d_xy_lst;
fvec time_last;
- fvec time_er_lst;
+ fvec dt2_lst;
fvec dz; /// !!!
fvec Sy[L1Constants::size::kMaxNstations];
@@ -883,10 +883,10 @@ void L1Algo::L1KFTrackFitterMuch()
x[ista][iVec] = x_temp[iVec];
y[ista][iVec] = y_temp[iVec];
time[ista][iVec] = hit.t;
- timeEr[ista][iVec] = hit.dt;
- d_u[ista][iVec] = hit.du;
- d_v[ista][iVec] = hit.dv;
- std::tie(d_xx[ista], d_xy[ista], d_yy[ista]) = sta[ista].FormXYCovarianceMatrix(d_u[ista], d_v[ista]);
+ dt2[ista][iVec] = hit.dt2;
+ du2[ista][iVec] = hit.du2;
+ dv2[ista][iVec] = hit.dv2;
+ std::tie(d_xx[ista], d_xy[ista], d_yy[ista]) = sta[ista].FormXYCovarianceMatrix(du2[ista], dv2[ista]);
// mom[ista][iVec] = hit.p;
z[ista][iVec] = hit.z;
@@ -899,7 +899,7 @@ void L1Algo::L1KFTrackFitterMuch()
x_first[iVec] = x[ista][iVec];
y_first[iVec] = y[ista][iVec];
time_first[iVec] = time[ista][iVec];
- time_er_fst[iVec] = timeEr[ista][iVec];
+ dt2_fst[iVec] = dt2[ista][iVec];
d_xx_fst[iVec] = d_xx[ista][iVec];
d_yy_fst[iVec] = d_yy[ista][iVec];
d_xy_fst[iVec] = d_xy[ista][iVec];
@@ -912,7 +912,7 @@ void L1Algo::L1KFTrackFitterMuch()
x_last[iVec] = x[ista][iVec];
y_last[iVec] = y[ista][iVec];
time_last[iVec] = time[ista][iVec];
- time_er_lst[iVec] = timeEr[ista][iVec];
+ dt2_lst[iVec] = dt2[ista][iVec];
d_xx_lst[iVec] = d_xx[ista][iVec];
d_yy_lst[iVec] = d_yy[ista][iVec];
d_xy_lst[iVec] = d_xy[ista][iVec];
@@ -962,7 +962,7 @@ void L1Algo::L1KFTrackFitterMuch()
FilterFirst(T, x_first, y_first, staFirst);
- FilterFirst(T1, x_first, y_first, time_first, time_er_fst, staFirst, d_xx_fst, d_yy_fst, d_xy_fst);
+ FilterFirst(T1, x_first, y_first, time_first, dt2_fst, staFirst, d_xx_fst, d_yy_fst, d_xy_fst);
fldZ1 = z[i];
@@ -1029,7 +1029,7 @@ void L1Algo::L1KFTrackFitterMuch()
L1Filter(T, info, v[i], w1);
T1.Filter(info, v[i], d_yy[i], w1);
- T1.FilterTime(time[i], timeEr[i], w1, sta[i].timeInfo);
+ T1.FilterTime(time[i], dt2[i], w1, sta[i].timeInfo);
}
if (i >= 8) {
@@ -1107,7 +1107,7 @@ void L1Algo::L1KFTrackFitterMuch()
info.sigma2 = d_yy[i];
L1Filter(T, info, v[i], w1);
T1.Filter(info, v[i], d_yy[i], w1);
- T1.FilterTime(time[i], timeEr[i], w1, sta[i].timeInfo);
+ T1.FilterTime(time[i], dt2[i], w1, sta[i].timeInfo);
}
}
// fit.L1AddHalfMaterial( T, sta[i].materialInfo, qp0 );
@@ -1156,7 +1156,7 @@ void L1Algo::L1KFTrackFitterMuch()
FilterFirst(T, x_last, y_last, staLast);
- FilterFirstL(T1, x_last, y_last, time_last, time_er_lst, staLast, d_xx_lst, d_yy_lst, d_xy_lst);
+ FilterFirstL(T1, x_last, y_last, time_last, dt2_lst, staLast, d_xx_lst, d_yy_lst, d_xy_lst);
qp0 = T.qp;
qp01 = tr.qp;
@@ -1250,7 +1250,7 @@ void L1Algo::L1KFTrackFitterMuch()
L1Filter(T, info, v[i], w1);
T1.Filter(info, v[i], d_yy[i], w1);
- T1.FilterTime(time[i], timeEr[i], w1, sta[i].timeInfo);
+ T1.FilterTime(time[i], dt2[i], w1, sta[i].timeInfo);
}
if (i < fNfieldStations - 1) {
@@ -1294,13 +1294,13 @@ void L1Algo::L1KFTrackFitterMuch()
}
L1UMeasurementInfo info = sta[i].frontInfo;
- // info.sigma2 = d_u[i] * d_u[i];
+ info.sigma2 = du2[i];
T1.Filter(info, u[i], info.sigma2, w1);
- info = sta[i].backInfo;
- // info.sigma2 = d_v[i] * d_v[i];
+ info = sta[i].backInfo;
+ info.sigma2 = dv2[i];
T1.Filter(info, v[i], info.sigma2, w1);
- T1.FilterTime(time[i], timeEr[i], w1, sta[i].timeInfo);
+ T1.FilterTime(time[i], dt2[i], w1, sta[i].timeInfo);
fldB2 = fldB1;
@@ -1598,7 +1598,8 @@ void L1Algo::FilterFirst(L1TrackParFit& track, fvec& x, fvec& y, fvec& t, L1Stat
track.fTr.chi2 = ZERO;
}
-void L1Algo::FilterFirst(L1TrackParFit& track, fvec& x, fvec& y, fvec& t, fvec& dt, L1Station& st, fvec& /*d_xx*/,
+
+void L1Algo::FilterFirst(L1TrackParFit& track, fvec& x, fvec& y, fvec& t, fvec& dt2, L1Station& st, fvec& /*d_xx*/,
fvec& /*d_yy*/, fvec& /*d_xy*/)
{
track.fTr.C00 = st.XYInfo.C00;
@@ -1621,7 +1622,7 @@ void L1Algo::FilterFirst(L1TrackParFit& track, fvec& x, fvec& y, fvec& t, fvec&
track.fTr.C52 = ZERO;
track.fTr.C53 = ZERO;
track.fTr.C54 = ZERO;
- track.fTr.C55 = dt * dt;
+ track.fTr.C55 = dt2;
track.fTr.x = x;
track.fTr.y = y;
@@ -1630,7 +1631,8 @@ void L1Algo::FilterFirst(L1TrackParFit& track, fvec& x, fvec& y, fvec& t, fvec&
track.fTr.chi2 = ZERO;
}
-void L1Algo::FilterFirst(L1TrackParFit& track, fvec& x, fvec& y, fvec& /*t*/, fvec& dt, L1Station& st)
+
+void L1Algo::FilterFirst(L1TrackParFit& track, fvec& x, fvec& y, fvec& /*t*/, fvec& dt2, L1Station& st)
{
track.fTr.C00 = st.XYInfo.C00;
track.fTr.C10 = st.XYInfo.C10;
@@ -1652,7 +1654,7 @@ void L1Algo::FilterFirst(L1TrackParFit& track, fvec& x, fvec& y, fvec& /*t*/, fv
track.fTr.C52 = ZERO;
track.fTr.C53 = ZERO;
track.fTr.C54 = ZERO;
- track.fTr.C55 = dt * dt;
+ track.fTr.C55 = dt2;
track.fTr.x = x;
track.fTr.y = y;
@@ -1660,8 +1662,7 @@ void L1Algo::FilterFirst(L1TrackParFit& track, fvec& x, fvec& y, fvec& /*t*/, fv
track.fTr.chi2 = ZERO;
}
-
-void L1Algo::FilterFirstL(L1TrackParFit& track, fvec& x, fvec& y, fvec& /*t*/, fvec& dt, L1Station& /*st*/, fvec& d_xx,
+void L1Algo::FilterFirstL(L1TrackParFit& track, fvec& x, fvec& y, fvec& /*t*/, fvec& dt2, L1Station& /*st*/, fvec& d_xx,
fvec& d_yy, fvec& d_xy)
{
// initialize covariance matrix
@@ -1688,7 +1689,7 @@ void L1Algo::FilterFirstL(L1TrackParFit& track, fvec& x, fvec& y, fvec& /*t*/, f
track.fTr.C52 = ZERO;
track.fTr.C53 = ZERO;
track.fTr.C54 = ZERO;
- track.fTr.C55 = dt * dt;
+ track.fTr.C55 = dt2;
track.fTr.x = x;
track.fTr.y = y;
diff --git a/reco/L1/L1Algo/L1TrackParFit.cxx b/reco/L1/L1Algo/L1TrackParFit.cxx
index fd310acdd4..974c6e10ff 100644
--- a/reco/L1/L1Algo/L1TrackParFit.cxx
+++ b/reco/L1/L1Algo/L1TrackParFit.cxx
@@ -144,7 +144,7 @@ void L1TrackParFit::FilterNoP(L1UMeasurementInfo& info, fvec u, fvec w)
fTr.C55 -= K5 * F5;
}
-void L1TrackParFit::FilterTime(fvec t, fvec dt, fvec w, fvec timeInfo)
+void L1TrackParFit::FilterTime(fvec t, fvec dt2, fvec w, fvec timeInfo)
{
// filter track with a time measurement
@@ -160,8 +160,6 @@ void L1TrackParFit::FilterTime(fvec t, fvec dt, fvec w, fvec timeInfo)
w.setZero(timeInfo <= fvec::Zero());
- fvec dt2 = dt * dt;
-
// when dt0 is much smaller than current time error,
// set track time exactly to the measurement value without filtering
// it helps to keep the initial time errors reasonably small
diff --git a/reco/L1/L1Algo/L1TrackParFit.h b/reco/L1/L1Algo/L1TrackParFit.h
index cb965a9b47..2a800cf265 100644
--- a/reco/L1/L1Algo/L1TrackParFit.h
+++ b/reco/L1/L1Algo/L1TrackParFit.h
@@ -49,7 +49,7 @@ public:
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 t0, fvec dt0, fvec w = 1., fvec timeInfo = 1.);
+ void FilterTime(fvec t, fvec dt2, fvec w, fvec timeInfo);
void FilterNoP(L1UMeasurementInfo& info, fvec u, fvec w = 1.);
void Extrapolate(fvec z_out, fvec qp0, const L1FieldRegion& F, const fvec& w);
--
GitLab