From eb563b4942b749c0a961f52bc665e2b5929eb258 Mon Sep 17 00:00:00 2001
From: "se.gorbunov" <se.gorbunov@gsi.de>
Date: Thu, 18 Aug 2022 16:22:23 +0000
Subject: [PATCH] L1: update of filtration
---
reco/L1/CbmL1.cxx | 9 +++
reco/L1/L1Algo/L1CATrackFinder.cxx | 46 +++++++++--
reco/L1/L1Algo/L1Filtration.h | 122 +++++++++++++++--------------
reco/L1/L1Algo/L1TrackPar.h | 111 ++++++++++++++++++--------
4 files changed, 194 insertions(+), 94 deletions(-)
diff --git a/reco/L1/CbmL1.cxx b/reco/L1/CbmL1.cxx
index 08f1d78d78..bbd90eb1ab 100644
--- a/reco/L1/CbmL1.cxx
+++ b/reco/L1/CbmL1.cxx
@@ -800,7 +800,16 @@ InitStatus CbmL1::Init()
trd2dIter2.SetMinLevelTripletStart(1);
trd2dIter2.SetPrimaryFlag(true);
+ // Initialize CA track finder iterations sequence
+
fpInitManager->SetCAIterationsNumberCrosscheck(1);
+ /*
+ fpInitManager->SetCAIterationsNumberCrosscheck(5);
+ fpInitManager->PushBackCAIteration(trackingIterFastPrim);
+ fpInitManager->PushBackCAIteration(trackingIterAllPrim);
+ fpInitManager->PushBackCAIteration(trackingIterAllPrimJump);
+ fpInitManager->PushBackCAIteration(trackingIterAllSec);
+ */
fpInitManager->PushBackCAIteration(trd2dIter2);
}
else {
diff --git a/reco/L1/L1Algo/L1CATrackFinder.cxx b/reco/L1/L1Algo/L1CATrackFinder.cxx
index 1e69f08ed7..0bd214f415 100644
--- a/reco/L1/L1Algo/L1CATrackFinder.cxx
+++ b/reco/L1/L1Algo/L1CATrackFinder.cxx
@@ -694,6 +694,8 @@ inline void L1Algo::findTripletsStep0( // input
// assert(T2.IsConsistent(true, n2_4));
+ // L1TrackPar tStore1 = T2;
+
L1UMeasurementInfo info = stam.frontInfo;
if (fUseHitErrors) info.sigma2 = du2 * du2;
@@ -708,7 +710,28 @@ inline void L1Algo::findTripletsStep0( // input
L1FilterNoField(T2, info, u_front_2);
}
- // assert(T2.IsConsistent(true, n2_4));
+ // L1TrackPar tStore2 = T2;
+
+ //assert(T2.IsConsistent(true, n2_4));
+
+ // a good place to debug the fit
+
+ /*
+ if (tStore1.IsConsistent(false, n2_4) && !tStore2.IsConsistent(false, n2_4)) {
+ cout << " i2 " << i2 << " dx2 " << dx2 << endl;
+ cout << "\n\n before filtration: \n\n" << endl;
+ tStore1.Print();
+ tStore1.IsConsistent(true, n2_4);
+ cout << "\n\n after filtration: \n\n" << endl;
+ tStore2.Print();
+ tStore2.IsConsistent(true, n2_4);
+ cout << "\n\n " << n2_4 << " vector elements are filled." << endl;
+ cout << "measurement: cos " << info.cos_phi[0] << " sin " << info.sin_phi[0] << " u " << u_front_2 << " s2 "
+ << sqrt(info.sigma2) << endl;
+ cout << "\n\n" << endl;
+ exit(0);
+ }
+ */
info = stam.backInfo;
if (fUseHitErrors) info.sigma2 = dv2 * dv2;
@@ -946,7 +969,7 @@ inline void L1Algo::findTripletsStep0( // input
dv_.resize(n3 / fvecLen);
timeR.resize(n3 / fvecLen);
timeER.resize(n3 / fvecLen);
- cout << "GetMaxTripletPerDoublets==" << fParameters.GetMaxTripletPerDoublets()
+ cout << "L1: GetMaxTripletPerDoublets==" << fParameters.GetMaxTripletPerDoublets()
<< " reached in findTripletsStep0()" << endl;
//assert(0);
break;
@@ -1259,13 +1282,17 @@ inline void L1Algo::findTripletsStep3( // input
ihitl_prev = 1 + hitsl_3[i3];
#ifdef DO_NOT_SELECT_TRIPLETS
- if (isec != TRACKS_FROM_TRIPLETS_ITERATION)
+ if (isec != TRACKS_FROM_TRIPLETS_ITERATION) {
+ if (chi2 > fTripletChi2Cut) { continue; }
+ }
+#else
+ if (chi2 > fTripletChi2Cut) { continue; }
#endif // DO_NOT_SELECT_TRIPLETS
- // assert(finite(chi2) && chi2 > 0);
+ // assert(finite(chi2) && chi2 > 0);
- if (!finite(chi2) || chi2 < 0) { continue; }
- if (chi2 > fTripletChi2Cut) { continue; }
+ if (!finite(chi2) || chi2 < 0) { continue; }
+ //if (!T3.IsEntryConsistent(false, i3_4)) { continue; }
fscal qp = T3.qp[i3_4];
@@ -2820,6 +2847,13 @@ inline void L1Algo::CAFindTrack(int ista, L1Branch& best_tr, unsigned char& best
fscal Cty = curr_trip->GetCty();
Cty += new_trip.GetCty();
+ // it shouldn't happen, but happens sometimes
+
+ if (!std::isfinite(dtx)) continue;
+ if (!std::isfinite(dty)) continue;
+ if (!std::isfinite(Ctx)) continue;
+ if (!std::isfinite(Cty)) continue;
+
assert(std::isfinite(dtx));
assert(std::isfinite(dty));
assert(std::isfinite(Ctx));
diff --git a/reco/L1/L1Algo/L1Filtration.h b/reco/L1/L1Algo/L1Filtration.h
index 2310f934fc..6e1a095e47 100644
--- a/reco/L1/L1Algo/L1Filtration.h
+++ b/reco/L1/L1Algo/L1Filtration.h
@@ -14,42 +14,45 @@
#define cnst const fvec
-inline void FilterTime(L1TrackPar& T, fvec t0, fvec dt0, fvec timeInfo = 1., fvec w = 1.)
+inline void FilterTime(L1TrackPar& T, fvec t, fvec dt, fvec timeInfo = 1., fvec w = 1.)
{
- fvec wi, zeta, zetawi, HCH;
- fvec F0, F1, F2, F3, F4, F5;
- fvec K1, K2, K3, K4, K5;
-
- zeta = T.t - t0;
+ // filter track with a time measurement
// F = CH'
- F0 = T.C50;
- F1 = T.C51;
+ fvec F0 = T.C50;
+ fvec F1 = T.C51;
+ fvec F2 = T.C52;
+ fvec F3 = T.C53;
+ fvec F4 = T.C54;
+ fvec F5 = T.C55;
- HCH = T.C55;
+ fvec HCH = T.C55;
- F2 = T.C52;
- F3 = T.C53;
- F4 = T.C54;
- F5 = T.C55;
+ w = w & (timeInfo > 0.f);
-#if 1 // use mask
- const fvec mask = (timeInfo > 0);
- wi = mask & w / (dt0 * dt0 + HCH);
- zetawi = zeta * wi;
- T.chi2 += mask & (zeta * zetawi);
-#else
- wi = w / (dt0 * dt0 + HCH);
- zetawi = zeta * wi;
- T.chi2 += zeta * zetawi;
-#endif // 0
+ 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
+ // the calculations in the covariance matrix are not affected
+
+ const fvec maskDoFilter = (HCH < dt2 * 16.f);
+ //const fvec maskDoFilter = _f32vec4_true;
+
+ fvec wi = w / (dt2 + 1.0000001f * HCH);
+ fvec zeta = T.t - t;
+ fvec zetawi = zeta / ((maskDoFilter & dt2) + HCH);
+
+ //T.chi2 += maskDoFilter & (zeta * zetawi);
+ T.chi2 += zeta * zeta * wi;
T.NDF += w;
- K1 = F1 * wi;
- K2 = F2 * wi;
- K3 = F3 * wi;
- K4 = F4 * wi;
- K5 = F5 * wi;
+ fvec K1 = F1 * wi;
+ fvec K2 = F2 * wi;
+ fvec K3 = F3 * wi;
+ fvec K4 = F4 * wi;
+ fvec K5 = F5 * wi;
T.x -= F0 * zetawi;
T.y -= F1 * zetawi;
@@ -81,42 +84,46 @@ inline void FilterTime(L1TrackPar& T, fvec t0, fvec dt0, fvec timeInfo = 1., fve
T.C55 -= K5 * F5;
}
+
inline void L1Filter(L1TrackPar& T, const L1UMeasurementInfo& info, fvec u, fvec w = 1.)
{
- fvec wi, zeta, zetawi, HCH;
- fvec F0, F1, F2, F3, F4, F5;
- fvec K1, K2, K3, K4, K5;
+ // filter the track T with an 1-D measurement u
- zeta = info.cos_phi * T.x + info.sin_phi * T.y - u;
+ fvec zeta = info.cos_phi * T.x + info.sin_phi * T.y - u;
// F = CH'
- F0 = info.cos_phi * T.C00 + info.sin_phi * T.C10;
- F1 = info.cos_phi * T.C10 + info.sin_phi * T.C11;
-
- HCH = (F0 * info.cos_phi + F1 * info.sin_phi);
-
- F2 = info.cos_phi * T.C20 + info.sin_phi * T.C21;
- F3 = info.cos_phi * T.C30 + info.sin_phi * T.C31;
- F4 = info.cos_phi * T.C40 + info.sin_phi * T.C41;
- F5 = info.cos_phi * T.C50 + info.sin_phi * T.C51;
-
-#if 0 // use mask
- const fvec mask = (HCH < info.sigma2 * 16.);
- wi = w/( (mask & info.sigma2) +HCH );
- zetawi = zeta *wi;
- T.chi2 += mask & (zeta * zetawi);
-#else
- wi = w / (info.sigma2 + HCH);
- zetawi = zeta * wi;
- T.chi2 += zeta * zetawi;
-#endif // 0
+ fvec F0 = info.cos_phi * T.C00 + info.sin_phi * T.C10;
+ fvec F1 = info.cos_phi * T.C10 + info.sin_phi * T.C11;
+ fvec F2 = info.cos_phi * T.C20 + info.sin_phi * T.C21;
+ fvec F3 = info.cos_phi * T.C30 + info.sin_phi * T.C31;
+ fvec F4 = info.cos_phi * T.C40 + info.sin_phi * T.C41;
+ fvec F5 = info.cos_phi * T.C50 + info.sin_phi * T.C51;
+
+ fvec HCH = (F0 * info.cos_phi + F1 * info.sin_phi);
+
+ // when the measurement error info.sigma2 is much smaller than the current track error HCH,
+ // move track exactly to the measurement without filtering
+ // it helps to keep the initial track errors reasonably small
+ // the calculations in the covariance matrix are not affected
+
+ //const fvec maskDoFilter = (HCH < info.sigma2 * 16.f);
+ const fvec 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 / (info.sigma2 + 1.0000001f * HCH);
+ fvec zetawi = zeta / ((maskDoFilter & info.sigma2) + HCH);
+
+ // T.chi2 += maskDoFilter & (zeta * zetawi);
+ T.chi2 += zeta * zeta * wi;
T.NDF += w;
- K1 = F1 * wi;
- K2 = F2 * wi;
- K3 = F3 * wi;
- K4 = F4 * wi;
- K5 = F5 * wi;
+ fvec K1 = F1 * wi;
+ fvec K2 = F2 * wi;
+ fvec K3 = F3 * wi;
+ fvec K4 = F4 * wi;
+ fvec K5 = F5 * wi;
T.x -= F0 * zetawi;
T.y -= F1 * zetawi;
@@ -148,6 +155,7 @@ inline void L1Filter(L1TrackPar& T, const L1UMeasurementInfo& info, fvec u, fvec
T.C55 -= K5 * F5;
}
+
inline void L1FilterNoField(L1TrackPar& T, const L1UMeasurementInfo& info, fvec u, fvec w = 1.)
{
fvec wi, zeta, zetawi, HCH;
diff --git a/reco/L1/L1Algo/L1TrackPar.h b/reco/L1/L1Algo/L1TrackPar.h
index 3361d53ae2..687a880c68 100644
--- a/reco/L1/L1Algo/L1TrackPar.h
+++ b/reco/L1/L1Algo/L1TrackPar.h
@@ -81,6 +81,13 @@ public:
void SetOneEntry(const int i0, const L1TrackPar& T1, const int i1);
+ fvec C(int i, int j) const
+ {
+ const fvec* c = &C00;
+ int ind = (j <= i) ? i * (1 + i) / 2 + j : j * (1 + j) / 2 + i;
+ return c[ind];
+ }
+
void Print(int i = -1) const;
void PrintCorrelations(int i = -1) const;
@@ -144,7 +151,10 @@ inline void L1TrackPar::Print(int i) const
inline void L1TrackPar::PrintCorrelations(int i) const
{
+
+ std::_Ios_Fmtflags flagSafe = std::cout.flags();
// std::cout.setf(std::ios::scientific, std::ios::floatfield);
+ std::cout << std::setprecision(6);
if (i == -1) {
fvec s0 = sqrt(C00);
@@ -192,6 +202,7 @@ inline void L1TrackPar::PrintCorrelations(int i) const
std::cout << " " << C50[i] / s5 / s0 << " " << C51[i] / s5 / s1 << " " << C52[i] / s5 / s2 << " "
<< C53[i] / s5 / s3 << " " << C54[i] / s5 / s4 << std::endl;
}
+ std::cout.flags(flagSafe);
}
inline void L1TrackPar::SetOneEntry(const int i0, const L1TrackPar& T1, const int i1)
@@ -229,48 +240,86 @@ inline void L1TrackPar::SetOneEntry(const int i0, const L1TrackPar& T1, const in
NDF[i0] = T1.NDF[i1];
} // SetOneEntry
-inline bool L1TrackPar::IsEntryConsistent(bool printWhenWrong, int i) const
+inline bool L1TrackPar::IsEntryConsistent(bool printWhenWrong, int k) const
{
bool ok = true;
- const fvec* v = &x;
- const fvec* v1 = &NDF;
- for (; v < v1; v++) {
- ok = ok && std::isfinite((*v)[i]);
+ // verify that all the numbers in the object are valid floats
+ const fvec* memberFirst = &x;
+ const fvec* memberLast = &NDF;
+ for (int i = 0; i < &memberLast - &memberFirst + 1; i++) {
+ if (!std::isfinite(memberFirst[i][k])) {
+ ok = false;
+ if (printWhenWrong) {
+ std::cout << " L1TrackPar member N " << i << ", vector entry " << k
+ << " is not a number: " << memberFirst[i][k];
+ }
+ }
}
- // verify diagonal elements
-
- ok = ok && (C00[i] > 0.f);
- ok = ok && (C11[i] > 0.f);
- ok = ok && (C22[i] > 0.f);
- ok = ok && (C33[i] > 0.f);
- ok = ok && (C44[i] > 0.f);
- ok = ok && (C55[i] > 0.f);
-
- // verify non-diagonal elements
- ok = ok && (C10[i] * C10[i] <= C11[i] * C00[i]);
- ok = ok && (C20[i] * C20[i] <= C22[i] * C00[i]);
- ok = ok && (C21[i] * C21[i] <= C22[i] * C11[i]);
+ // verify diagonal elements.
+ // Cii is a squared dispersion of i-th parameter, it must be positive
- ok = ok && (C30[i] * C30[i] <= C33[i] * C00[i]);
- ok = ok && (C31[i] * C31[i] <= C33[i] * C11[i]);
- ok = ok && (C32[i] * C32[i] <= C33[i] * C22[i]);
+ for (int i = 0; i < 6; i++) {
+ if (C(i, i)[k] <= 0.f) {
+ ok = false;
+ if (printWhenWrong) {
+ std::cout << " L1TrackPar: C[" << i << "," << i << "], vec entry " << k << " is not positive: " << C(i, i)[k]
+ << std::endl;
+ }
+ }
+ }
- ok = ok && (C40[i] * C40[i] <= C44[i] * C00[i]);
- ok = ok && (C41[i] * C41[i] <= C44[i] * C11[i]);
- ok = ok && (C42[i] * C42[i] <= C44[i] * C22[i]);
- ok = ok && (C43[i] * C43[i] <= C44[i] * C33[i]);
+ // verify non-diagonal elements.
+ // Cij/sqrt(Cii*Cjj) is a correlation between i-th and j-th parameter,
+ // it must belong to [-1,1]
+
+ for (int i = 1; i < 6; i++) {
+ for (int j = 0; j < i; j++) {
+ double tolerance = 1.0;
+ if (C(i, j)[k] * C(i, j)[k] > tolerance * (C(i, i)[k] * C(j, j)[k])) {
+ ok = false;
+ if (printWhenWrong) {
+ std::cout << " L1TrackPar: correlation [" << i << "," << j << "], vec entry " << k
+ << " is too large: " << C(i, i)[k] / sqrt(C(i, i)[k] * C(j, j)[k]) << std::endl;
+ }
+ }
+ }
+ }
- ok = ok && (C50[i] * C50[i] <= C55[i] * C00[i]);
- ok = ok && (C51[i] * C51[i] <= C55[i] * C11[i]);
- ok = ok && (C52[i] * C52[i] <= C55[i] * C22[i]);
- ok = ok && (C53[i] * C53[i] <= C55[i] * C33[i]);
- ok = ok && (C54[i] * C54[i] <= C55[i] * C44[i]);
+ // verify triplets of correlations
+ // Kxy * Kxy + Kxz * Kxz + Kyz * Kyz <= 1 + 2 * Kxy * Kxz * Kyz
+
+ for (int i = 2; i < 6; i++) {
+ for (int j = 1; j < i; j++) {
+ for (int m = 0; m < j; m++) {
+ double tolerance = 1.0;
+ double Cxx = C(i, i)[k];
+ double Cyy = C(j, j)[k];
+ double Czz = C(m, m)[k];
+ double Cxy = C(i, j)[k];
+ double Cxz = C(i, m)[k];
+ double Cyz = C(j, m)[k];
+ if (Cxx * Cyz * Cyz + Cyy * Cxz * Cxz + Czz * Cxy * Cxy
+ > tolerance * (Cxx * Cyy * Czz + 2. * Cxy * Cyz * Cxz)) {
+ ok = false;
+ if (printWhenWrong) {
+ double Kxy = Cxy / sqrt(Cxx * Cyy);
+ double Kxz = Cxz / sqrt(Cxx * Czz);
+ double Kyz = Cyz / sqrt(Cyy * Czz);
+ std::cout << " L1TrackPar: correlations between parametetrs " << i << ", " << j << ", " << m
+ << ", vec entry " << k << " are wrong: " << Kxy << " " << Kxz << " " << Kyz << std::endl;
+ std::cout << " inequation: " << Kxy * Kxy + Kxz * Kxz + Kyz * Kyz << " > " << 1 + 2 * Kxy * Kxz * Kyz
+ << std::endl;
+ }
+ }
+ }
+ }
+ }
if (!ok && printWhenWrong) {
std::cout << "L1TrackPar parameters are not consistent: " << std::endl;
- Print(i);
+ Print(k);
}
return ok;
}
--
GitLab