From 61cf78a66de3d80a168caf4ba2f9b48b5de5d40f Mon Sep 17 00:00:00 2001
From: "se.gorbunov" <se.gorbunov@gsi.de>
Date: Wed, 15 Feb 2023 00:08:34 +0000
Subject: [PATCH] L1 fit: unify track guess routines
---
reco/L1/CbmL1Performance.cxx | 3 -
reco/L1/L1Algo/L1Algo.h | 12 --
reco/L1/L1Algo/L1Fit.cxx | 129 ++++++++++++++++++
reco/L1/L1Algo/L1Fit.h | 3 +
reco/L1/L1Algo/L1TrackFitter.cxx | 226 ++-----------------------------
5 files changed, 145 insertions(+), 228 deletions(-)
diff --git a/reco/L1/CbmL1Performance.cxx b/reco/L1/CbmL1Performance.cxx
index dad4eb24ed..4b9f07235f 100644
--- a/reco/L1/CbmL1Performance.cxx
+++ b/reco/L1/CbmL1Performance.cxx
@@ -54,9 +54,6 @@
#include <map>
#include <vector>
-#include <cmath>
-#include <math.h>
-
#include "L1Algo/L1Algo.h"
#include "L1Algo/L1Def.h"
#include "L1Algo/L1Fit.h"
diff --git a/reco/L1/L1Algo/L1Algo.h b/reco/L1/L1Algo/L1Algo.h
index 5a8d7177cf..006a21fc24 100644
--- a/reco/L1/L1Algo/L1Algo.h
+++ b/reco/L1/L1Algo/L1Algo.h
@@ -336,18 +336,6 @@ public:
int istal, int istam, int istar, L1TrackPar* T_1, L1FieldRegion* fld_1, L1HitIndex_t* hitsl_1, Tindex& n_2,
L1Vector<L1HitIndex_t>& i1_2, L1Vector<L1HitIndex_t>& hitsm_2);
-
- /// ------ Subroutines used by L1Algo::KFTrackFitter() ------
-
- 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, fmask* w_time, fvec& dt2_last);
-
- void FilterFirst(L1Fit& track, fvec& x, fvec& y, fvec& t, fvec& dt2, fvec& dx2, fvec& dy2, fvec& dxy);
-
-
#ifdef DRAW
L1AlgoDraw* draw {nullptr};
void DrawRecoTracksTime(const L1Vector<CbmL1Track>& tracks);
diff --git a/reco/L1/L1Algo/L1Fit.cxx b/reco/L1/L1Algo/L1Fit.cxx
index 1b17d6982f..3a99d903db 100644
--- a/reco/L1/L1Algo/L1Fit.cxx
+++ b/reco/L1/L1Algo/L1Fit.cxx
@@ -1725,4 +1725,133 @@ void L1Fit::FilterChi2(const L1UMeasurementInfo& info, cnst& x, cnst& y, cnst& C
chi2 += zeta * zeta / (du2 + HCH);
}
+
+void L1Fit::GuessTrack(const fvec& trackZ, const fvec hitX[], const fvec hitY[], const fvec hitZ[], const fvec hitT[],
+ const fvec By[], const fmask hitW[], const fmask hitWtime[], int NHits)
+{
+ // 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%).
+
+ const fvec c_light(0.000299792458), c_light_i(fvec(1.) / c_light);
+
+ fvec A0 = 0., A1 = 0., A2 = 0., A3 = 0., A4 = 0., A5 = 0.;
+ fvec a0 = 0., a1 = 0., a2 = 0.;
+ fvec b0 = 0., b1 = 0., b2 = 0.;
+
+ fvec time = 0.;
+ fmask isTimeSet = fmask::Zero();
+
+ fvec prevZ = 0.;
+ fvec sy = 0., Sy = 0.; // field integrals
+
+ for (int i = 0; i < NHits; i++) {
+
+ fvec w = iif(hitW[i], fvec::One(), fvec::Zero());
+
+ fmask setTime = (!isTimeSet) && hitWtime[i] && hitW[i];
+ time(setTime) = hitT[i];
+ isTimeSet = isTimeSet || setTime;
+
+ fvec x = hitX[i];
+ fvec y = hitY[i];
+ fvec z = hitZ[i] - trackZ;
+
+ {
+ fvec dZ = z - prevZ;
+ Sy(hitW[i]) += dZ * sy + fvec(0.5) * dZ * dZ * By[i];
+ sy(hitW[i]) += dZ * By[i];
+ prevZ(hitW[i]) = z;
+ }
+
+ fvec S = Sy;
+
+ fvec wz = w * z;
+ fvec wS = w * S;
+
+ A0 += w;
+ A1 += wz;
+ A2 += wz * z;
+ A3 += wS;
+ A4 += wS * z;
+ A5 += wS * S;
+
+ a0 += w * x;
+ a1 += wz * x;
+ a2 += wS * x;
+
+ b0 += w * y;
+ b1 += wz * y;
+ b2 += wS * y;
+ }
+
+ fvec m00 = A0;
+ fvec m01 = A1;
+ fvec m02 = A3;
+
+ fvec m11 = A2;
+ fvec m12 = A4;
+
+ fvec m22 = A5;
+
+ fvec m21 = m12;
+
+ // { m00 m01 m02 } ( a0 )
+ // { m01 m11 m12 } = x * ( a1 )
+ // { m02 m21 m22 } ( a2 )
+
+ { // triangulation row 0
+ m11 = m00 * m11 - m01 * m01;
+ m12 = m00 * m12 - m01 * m02;
+ a1 = m00 * a1 - m01 * a0;
+
+ m21 = m00 * m21 - m02 * m01;
+ m22 = m00 * m22 - m02 * m02;
+ a2 = m00 * a2 - m02 * a0;
+ }
+
+ { // triangulation step row 1
+ m22 = m11 * m22 - m21 * m12;
+ a2 = m11 * a2 - m21 * a1;
+ }
+
+ fvec L = 0.;
+ { // diagonalization row 2
+ L(abs(m22) > fvec(1.e-4)) = a2 / m22;
+ a1 -= L * m12;
+ a0 -= L * m02;
+ }
+
+ { // diagonalization row 1
+ fTr.tx = a1 / m11;
+ a0 -= fTr.tx * m01;
+ }
+
+ { // diagonalization row 0
+ fTr.x = a0 / m00;
+ }
+
+ fvec txtx1 = fvec(1.) + fTr.tx * fTr.tx;
+ L = L / txtx1;
+ fvec L1 = L * fTr.tx;
+
+ A1 = A1 + A3 * L1;
+ A2 = A2 + (A4 + A4 + A5 * L1) * L1;
+ b1 += b2 * L1;
+
+ // { A0 A1 } = x * ( b0 )
+ // { A1 A2 } ( b1 )
+
+ A2 = A0 * A2 - A1 * A1;
+ b1 = A0 * b1 - A1 * b0;
+
+ fTr.ty = b1 / A2;
+ fTr.y = (b0 - A1 * fTr.ty) / A0;
+
+ fTr.qp = -L * c_light_i / sqrt(txtx1 + fTr.ty * fTr.ty);
+ fTr.t = time;
+ fTr.z = trackZ;
+ fTr.vi = 0.;
+ fQp0 = fTr.qp;
+}
+
+
#undef cnst
diff --git a/reco/L1/L1Algo/L1Fit.h b/reco/L1/L1Algo/L1Fit.h
index 72fe29a978..6ca7a9b1f9 100644
--- a/reco/L1/L1Algo/L1Fit.h
+++ b/reco/L1/L1Algo/L1Fit.h
@@ -130,6 +130,9 @@ public:
static void FilterChi2(const L1UMeasurementInfo& info, cnst& x, cnst& y, cnst& C00, cnst& C10, cnst& C11, fvec& chi2,
cnst& u, cnst& du2);
+ void GuessTrack(const fvec& trackZ, const fvec hitX[], const fvec hitY[], const fvec hitZ[], const fvec hitT[],
+ const fvec By[], const fmask hitW[], const fmask hitWtime[], int NHits);
+
private:
fmask fMask {fmask::One()}; // mask of active elements in simd vectors
fvec fMaskF {fvec::One()}; // floating point representation of fMask
diff --git a/reco/L1/L1Algo/L1TrackFitter.cxx b/reco/L1/L1Algo/L1TrackFitter.cxx
index 4ed3a4b1dc..4676463254 100644
--- a/reco/L1/L1Algo/L1TrackFitter.cxx
+++ b/reco/L1/L1Algo/L1TrackFitter.cxx
@@ -75,8 +75,8 @@ void L1Algo::L1KFTrackFitter()
fvec time_first;
fvec dt2_first;
- fvec x_last, x_2last;
- fvec y_last, y_2last;
+ fvec x_last;
+ fvec y_last;
fvec d_xx_lst;
fvec d_yy_lst;
fvec d_xy_lst;
@@ -85,7 +85,7 @@ void L1Algo::L1KFTrackFitter()
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 By[L1Constants::size::kMaxNstations];
fmask w[L1Constants::size::kMaxNstations];
fmask w_time[L1Constants::size::kMaxNstations]; // !!!
@@ -95,7 +95,7 @@ void L1Algo::L1KFTrackFitter()
fvec fldZ1;
fvec fldZ2;
fvec z_start;
- fvec z_end, z_2last;
+ fvec z_end;
L1FieldValue fB[L1Constants::size::kMaxNstations], fB_temp _fvecalignment;
@@ -182,34 +182,16 @@ void L1Algo::L1KFTrackFitter()
time_last[iVec] = time[ista][iVec];
dt2_last[iVec] = dt2[ista][iVec];
}
- else if (ih == nHitsTrack - 2) {
- z_2last[iVec] = z[ista][iVec];
- x_2last[iVec] = x[ista][iVec];
- y_2last[iVec] = y[ista][iVec];
- }
}
- fscal prevZ = z_end[iVec];
- fscal cursy = 0., curSy = 0.;
for (int ih = nHitsTrack - 1; ih >= 0; ih--) {
const int ista = iSta[ih];
const L1Station& st = fParameters.GetStation(ista);
- const fscal curZ = z[ista][iVec];
- fscal dZ = curZ - prevZ;
- fscal By = st.fieldSlice.cy[0][0];
- curSy += dZ * cursy + dZ * dZ * By / 2.;
- cursy += dZ * By;
- Sy[ista][iVec] = curSy;
- prevZ = curZ;
+ By[ista][iVec] = st.fieldSlice.cy[0][0];
}
}
- 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 {
- GuessVec(fit, x, y, z, Sy, w, nStations, &z_end, time, w_time);
- }
+ fit.GuessTrack(z_end, x, y, z, time, By, w, w_time, nStations);
if (kGlobal == fTrackingMode || kMcbm == fTrackingMode) { tr.qp = fvec(1. / 1.1); }
@@ -233,7 +215,6 @@ void L1Algo::L1KFTrackFitter()
tr.t = time_last;
tr.NDF = -3.0;
-
fldZ1 = z[ista];
sta[ista].fieldSlice.GetFieldValue(tr.x, tr.y, fldB1);
@@ -324,7 +305,13 @@ void L1Algo::L1KFTrackFitter()
ista = 0;
- FilterFirst(fit, x_first, y_first, time_first, dt2_first, d_xx_fst, d_yy_fst, d_xy_fst);
+ tr.ResetErrors(d_xx_fst, d_yy_fst, 0.1, 0.1, 1., dt2_first, 1.e2);
+ tr.C10 = d_xy_fst;
+
+ tr.x = x_first;
+ tr.y = y_first;
+ tr.t = time_first;
+ tr.NDF = -3.0;
fit.SetQp0(tr.qp);
@@ -374,190 +361,3 @@ void L1Algo::L1KFTrackFitter()
} // iter
}
}
-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, fmask* /*w_time*/, fvec& dt2_last)
-{
- L1TrackPar& tr = track.Tr();
-
- fvec dzi = fvec(1.) / (z_end - z_2last);
-
- tr.x = x_last;
- tr.y = y_last;
- tr.tx = (x_last - x_2last) * dzi;
- tr.ty = (y_last - y_2last) * dzi;
- tr.t = time_last;
- tr.z = z_end;
-
- tr.ResetErrors(1., 1., 1., 1., 1., dt2_last, 1.e2);
- tr.NDF = fvec(2.);
-}
-
-
-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%).
-{
-
- fvec A0, A1 = 0., A2 = 0., A3 = 0., A4 = 0., A5 = 0., a0, a1 = 0., a2 = 0., b0, b1 = 0., b2 = 0.;
- fvec z0, x, y, z, S, w, wz, wS;
-
- int i = NHits - 1;
- if (zCur) z0 = *zCur;
- else
- z0 = zV[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 = iif(wV[i], fvec::One(), fvec::Zero());
- z = zV[i] - z0;
- S = Sy[i];
- wz = w * z;
- wS = w * S;
- A0 += w;
- A1 += wz;
- A2 += wz * z;
- A3 += wS;
- A4 += wS * z;
- A5 += wS * S;
- a0 += w * x;
- a1 += wz * x;
- a2 += wS * x;
- b0 += w * y;
- b1 += wz * y;
- b2 += wS * y;
- }
-
- fvec A3A3 = A3 * A3;
- fvec A3A4 = A3 * A4;
- fvec A1A5 = A1 * A5;
- fvec A2A5 = A2 * A5;
- fvec A4A4 = A4 * A4;
-
- fvec det = fvec::One() / (-A2 * A3A3 + A1 * (A3A4 + A3A4 - A1A5) + A0 * (A2A5 - A4A4));
- fvec Ai0 = (-A4A4 + A2A5);
- fvec Ai1 = (A3A4 - A1A5);
- fvec Ai2 = (-A3A3 + A0 * A5);
- fvec Ai3 = (-A2 * A3 + A1 * A4);
- fvec Ai4 = (A1 * A3 - A0 * A4);
- fvec Ai5 = (-A1 * A1 + A0 * A2);
-
- fvec L, L1;
- tr.x = (Ai0 * a0 + Ai1 * a1 + Ai3 * a2) * det;
- tr.tx = (Ai1 * a0 + Ai2 * a1 + Ai4 * a2) * det;
- fvec txtx1 = fvec(1.) + tr.tx * tr.tx;
- L = (Ai3 * a0 + Ai4 * a1 + Ai5 * a2) * det / txtx1;
- L1 = L * tr.tx;
- A1 = A1 + A3 * L1;
- A2 = A2 + (A4 + A4 + A5 * L1) * L1;
- b1 += b2 * L1;
- det = fvec::One() / (-A1 * A1 + A0 * A2);
-
- tr.y = (A2 * b0 - A1 * b1) * det;
- tr.ty = (-A1 * b0 + A0 * b1) * det;
- tr.qp = -L * c_light_i / sqrt(txtx1 + tr.ty * tr.ty);
- tr.z = z0;
-}
-
-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();
-
- fvec A0, A1 = 0., A2 = 0., A3 = 0., A4 = 0., A5 = 0., a0, a1 = 0., a2 = 0., b0, b1 = 0., b2 = 0.;
- fvec z0, x, y, z, S, w, wz, wS, time;
-
- time = fvec(0.);
-
- int i = NHits - 1;
- if (zCur) { z0 = *zCur; }
- else {
- z0 = zV[i];
- }
-
- w = iif(wV[i], fvec::One(), fvec::Zero());
- A0 = w;
- a0 = w * xV[i];
- b0 = w * yV[i];
-
- fvec nhits = 0;
- for (i = 0; i < NHits; i++) {
- x = xV[i];
- y = yV[i];
- w = iif(wV[i], fvec::One(), fvec::Zero());
- ;
- if (timeV) nhits(w_time[i]) += fvec::One();
- z = zV[i] - z0;
- S = Sy[i];
-
- wz = w * z;
- wS = w * S;
-
- A0 += w;
- A1 += wz;
- A2 += wz * z;
- A3 += wS;
- A4 += wS * z;
- A5 += wS * S;
- a0 += w * x;
- a1 += wz * x;
- a2 += wS * x;
- b0 += w * y;
- b1 += wz * y;
- b2 += wS * y;
- if (timeV) {
- 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.));
- }
- }
-
- fvec A3A3 = A3 * A3;
- fvec A3A4 = A3 * A4;
- fvec A1A5 = A1 * A5;
- fvec A2A5 = A2 * A5;
- fvec A4A4 = A4 * A4;
-
- fvec det = fvec::One() / (-A2 * A3A3 + A1 * (A3A4 + A3A4 - A1A5) + A0 * (A2A5 - A4A4));
- fvec Ai0 = (-A4A4 + A2A5);
- fvec Ai1 = (A3A4 - A1A5);
- fvec Ai2 = (-A3A3 + A0 * A5);
- fvec Ai3 = (-A2 * A3 + A1 * A4);
- fvec Ai4 = (A1 * A3 - A0 * A4);
- fvec Ai5 = (-A1 * A1 + A0 * A2);
-
- fvec L, L1;
- tr.x = (Ai0 * a0 + Ai1 * a1 + Ai3 * a2) * det;
- tr.tx = (Ai1 * a0 + Ai2 * a1 + Ai4 * a2) * det;
- fvec txtx1 = fvec(1.) + tr.tx * tr.tx;
- L = (Ai3 * a0 + Ai4 * a1 + Ai5 * a2) * det / txtx1;
- L1 = L * tr.tx;
- A1 = A1 + A3 * L1;
- A2 = A2 + (A4 + A4 + A5 * L1) * L1;
- b1 += b2 * L1;
- det = fvec::One() / (-A1 * A1 + A0 * A2);
- tr.y = (A2 * b0 - A1 * b1) * det;
- tr.ty = (-A1 * b0 + A0 * b1) * det;
- tr.qp = -L * c_light_i / sqrt(txtx1 + tr.ty * tr.ty);
- if (timeV) { tr.t = iif(nhits > fvec::Zero(), time / nhits, fvec::Zero()); }
-
- tr.z = z0;
-}
-
-
-void L1Algo::FilterFirst(L1Fit& track, fvec& x, fvec& y, fvec& t, fvec& dt2, fvec& d_xx, fvec& d_yy, fvec& d_xy)
-{
- L1TrackPar& tr = track.Tr();
-
- tr.ResetErrors(d_xx, d_yy, 0.1, 0.1, 1., dt2, 1.e2);
- tr.C10 = d_xy;
-
- tr.x = x;
- tr.y = y;
- tr.t = t;
- tr.NDF = -3.0;
-}
--
GitLab