From 17ed65c48a45eb17e5d21e8e881c2938ab9608d5 Mon Sep 17 00:00:00 2001
From: "se.gorbunov" <se.gorbunov@gsi.de>
Date: Thu, 9 Feb 2023 16:23:11 +0000
Subject: [PATCH] L1 fit: unify fit routines
---
reco/L1/CbmL1.cxx | 2 +-
reco/L1/CbmL1Hit.h | 2 +-
reco/L1/CbmL1Performance.cxx | 71 +++---
reco/L1/L1Algo/L1Algo.h | 2 +-
reco/L1/L1Algo/L1BranchExtender.cxx | 95 ++++----
reco/L1/L1Algo/L1CATrackFinder.cxx | 88 +++----
reco/L1/L1Algo/L1CloneMerger.cxx | 21 +-
reco/L1/L1Algo/L1TrackFitter.cxx | 282 -----------------------
reco/L1/L1Algo/L1TrackPar.h | 112 ++++-----
reco/L1/L1Algo/L1TrackParFit.cxx | 159 ++++++++++++-
reco/L1/L1Algo/L1TrackParFit.h | 72 +++++-
reco/L1/ParticleFinder/CbmL1PFFitter.cxx | 282 +++++++++++------------
12 files changed, 544 insertions(+), 644 deletions(-)
diff --git a/reco/L1/CbmL1.cxx b/reco/L1/CbmL1.cxx
index eabdfb260f..fa44b00e76 100644
--- a/reco/L1/CbmL1.cxx
+++ b/reco/L1/CbmL1.cxx
@@ -926,7 +926,7 @@ void CbmL1::Reconstruct(CbmEvent* event)
else {
fvSortedStsHitsIndexes.clear();
fvSortedStsHitsIndexes.reserve(nStsHits);
- for (unsigned int i = 0; i < nStsHits; i++)
+ for (int i = 0; i < nStsHits; i++)
fvSortedStsHitsIndexes.push_back(i);
}
// -----------------------------------------------------------------------
diff --git a/reco/L1/CbmL1Hit.h b/reco/L1/CbmL1Hit.h
index 4f7e645b46..64b274f706 100644
--- a/reco/L1/CbmL1Hit.h
+++ b/reco/L1/CbmL1Hit.h
@@ -18,7 +18,7 @@ public:
int hitId = 0; ///< index of L1Hit in fInputData::fvHits array. Should be equal to index of this in L1->vHits
int extIndex = 0; ///< index of hit in the TClonesArray array
- int Det = 0; ///< station index
+ int Det = 0; ///< detector ID (mvd/sts/etc)
float x = 0.f; ///< measured X coordinate
float y = 0.f; ///< measured Y coordinate
float t = 0.f; ///< measured time
diff --git a/reco/L1/CbmL1Performance.cxx b/reco/L1/CbmL1Performance.cxx
index 5b38e7d06c..08af527039 100644
--- a/reco/L1/CbmL1Performance.cxx
+++ b/reco/L1/CbmL1Performance.cxx
@@ -56,9 +56,7 @@
#include "L1Algo/L1Algo.h"
#include "L1Algo/L1Def.h"
-#include "L1Algo/L1Extrapolation.h" // for vertex pulls
-#include "L1Algo/L1Fit.h" // for vertex pulls
-
+#include "L1Algo/L1TrackParFit.h"
using std::cout;
using std::endl;
@@ -1103,10 +1101,9 @@ void CbmL1::TrackFitPerformance()
static bool first_call = 1;
- L1Fit fit;
+ L1TrackParFit fit;
fit.SetParticleMass(fpAlgo->GetDefaultParticleMass());
-
if (first_call) {
first_call = 0;
@@ -1231,16 +1228,19 @@ void CbmL1::TrackFitPerformance()
CbmL1MCTrack mc = *(it->GetMCTracks()[0]);
- L1TrackPar trPar(it->T, it->C);
+ {
+ L1TrackPar trPar(it->T, it->C);
+ fit.fTr = trPar;
+ }
+ fit.fQp0 = fit.fTr.qp;
+ const L1TrackPar& tr = fit.fTr;
L1FieldRegion fld _fvecalignment;
fld.SetUseOriginalField();
CbmL1MCPoint& mcP = fvMCPoints[mc.Points[0]];
- L1Extrapolate(trPar, mcP.zIn, trPar.qp, fld);
-
- const L1TrackPar& tr = trPar;
+ fit.Extrapolate(mcP.zIn, fit.fQp0, fld, fvec::One());
double dx = tr.x[0] - mcP.xIn;
double dy = tr.y[0] - mcP.yIn;
@@ -1334,14 +1334,19 @@ void CbmL1::TrackFitPerformance()
const int last_station = fvHitStore[it->Hits.back()].iStation;
CbmL1MCTrack mc = *(it->GetMCTracks()[0]);
- L1TrackPar trPar(it->TLast, it->CLast);
+
L1FieldRegion fld _fvecalignment;
fld.SetUseOriginalField();
- CbmL1MCPoint& mcP = fvMCPoints[iMC];
- L1Extrapolate(trPar, mcP.zOut, trPar.qp, fld);
+ {
+ L1TrackPar trPar(it->TLast, it->CLast);
+ fit.fTr = trPar;
+ }
+ fit.fQp0 = fit.fTr.qp;
+ const L1TrackPar& tr = fit.fTr;
- const L1TrackPar& tr = trPar;
+ CbmL1MCPoint& mcP = fvMCPoints[iMC];
+ fit.Extrapolate(mcP.zOut, fit.fQp0, fld, fvec::One());
h_fitL[0]->Fill((tr.x[0] - mcP.xOut) * 1.e4);
h_fitL[1]->Fill((tr.y[0] - mcP.yOut) * 1.e4);
@@ -1389,7 +1394,12 @@ void CbmL1::TrackFitPerformance()
{ // vertex
CbmL1MCTrack mc = *(it->GetMCTracks()[0]);
- L1TrackPar trPar(it->T, it->C);
+ {
+ L1TrackPar trTmp(it->T, it->C);
+ fit.fTr = trTmp;
+ }
+ fit.fQp0 = fit.fTr.qp;
+ L1TrackPar& tr = fit.fTr;
// if (mc.mother_ID != -1) { // secondary
if (!mc.IsPrimary()) { // secondary
@@ -1397,7 +1407,7 @@ void CbmL1::TrackFitPerformance()
{ // extrapolate to vertex
L1FieldRegion fld _fvecalignment;
fld.SetUseOriginalField();
- L1Extrapolate(trPar, mc.z, trPar.qp, fld);
+ fit.Extrapolate(mc.z, fit.fQp0, fld, fvec::One());
// add material
const int fSta = fvHitStore[it->Hits[0]].iStation;
const int dir = int((mc.z - fpAlgo->GetParameters()->GetStation(fSta).fZ[0])
@@ -1407,15 +1417,11 @@ void CbmL1::TrackFitPerformance()
&& (dir * (mc.z - fpAlgo->GetParameters()->GetStation(iSta).fZ[0]) > 0);
iSta += dir) {
// cout << iSta << " " << dir << endl;
- fit.L1AddMaterial(trPar, fpAlgo->GetParameters()->GetMaterialThickness(iSta, trPar.x, trPar.y), trPar.qp,
- fvec::One());
-
- if (iSta + dir == fNMvdStations - 1) fit.L1AddPipeMaterial(trPar, trPar.qp, 1);
+ fit.AddMsInMaterial(fpAlgo->GetParameters()->GetMaterialThickness(iSta, fit.fTr.x, fit.fTr.y), fit.fQp0,
+ fvec::One());
}
}
- if (mc.z != trPar.z[0]) continue;
-
- const L1TrackPar& tr = trPar;
+ if (mc.z != tr.z[0]) continue;
// static int good = 0;
// static int bad = 0;
@@ -1463,21 +1469,16 @@ void CbmL1::TrackFitPerformance()
&& (dir * (mc.z - fpAlgo->GetParameters()->GetStation(iSta).fZ[0]) > 0);
iSta += dir) {
- L1Extrapolate(trPar, fpAlgo->GetParameters()->GetStation(iSta).fZ[0], trPar.qp, fld);
- fit.L1AddMaterial(trPar, fpAlgo->GetParameters()->GetMaterialThickness(iSta, trPar.x, trPar.y), trPar.qp,
- 1);
- fit.EnergyLossCorrection(trPar, fpAlgo->GetParameters()->GetMaterialThickness(iSta, trPar.x, trPar.y),
- trPar.qp, fvec(1.), fvec(1.));
- if (iSta + dir == fNMvdStations - 1) {
- fit.L1AddPipeMaterial(trPar, trPar.qp, 1);
- fit.EnergyLossCorrection(trPar, fit.PipeRadThick, trPar.qp, fvec(1.f), fvec(1.f));
- }
+ fit.Extrapolate(fpAlgo->GetParameters()->GetStation(iSta).fZ, fit.fQp0, fld, fvec::One());
+
+ fit.AddMsInMaterial(fpAlgo->GetParameters()->GetMaterialThickness(iSta, fit.fTr.x, fit.fTr.y), fit.fQp0,
+ fvec::One());
+ fit.EnergyLossCorrection(fpAlgo->GetParameters()->GetMaterialThickness(iSta, fit.fTr.x, fit.fTr.y),
+ fit.fQp0, fvec::One(), fvec::One());
}
- L1Extrapolate(trPar, mc.z, trPar.qp, fld);
+ fit.Extrapolate(mc.z, fit.fQp0, fld, fvec::One());
}
- if (mc.z != trPar.z[0]) continue;
-
- const L1TrackPar& tr = trPar;
+ if (mc.z != tr.z[0]) continue;
double dx = tr.x[0] - mc.x;
double dy = tr.y[0] - mc.y;
diff --git a/reco/L1/L1Algo/L1Algo.h b/reco/L1/L1Algo/L1Algo.h
index 798ecf952c..3aa9afa64d 100644
--- a/reco/L1/L1Algo/L1Algo.h
+++ b/reco/L1/L1Algo/L1Algo.h
@@ -370,7 +370,7 @@ public:
void CATrackFinder();
/// Track fitting procedures
- void KFTrackFitter_simple(); // version, which use procedured used during the reconstruction
+
void L1KFTrackFitter(); // version from SIMD-KF benchmark
diff --git a/reco/L1/L1Algo/L1BranchExtender.cxx b/reco/L1/L1Algo/L1BranchExtender.cxx
index 07e8185119..770db8283a 100644
--- a/reco/L1/L1Algo/L1BranchExtender.cxx
+++ b/reco/L1/L1Algo/L1BranchExtender.cxx
@@ -6,9 +6,6 @@
#include "L1Algo.h"
#include "L1Branch.h"
-#include "L1Extrapolation.h"
-#include "L1Filtration.h"
-#include "L1Fit.h"
#include "L1HitArea.h"
#include "L1HitPoint.h"
#include "L1Track.h"
@@ -25,12 +22,15 @@ using std::endl;
/// dir - 0 - forward, 1 - backward
/// qp0 - momentum for extrapolation
/// initialize - should be params ititialized. 1 - yes.
-void L1Algo::BranchFitterFast(const L1Branch& t, L1TrackPar& T, const bool dir, const fvec qp0, const bool initParams)
+void L1Algo::BranchFitterFast(const L1Branch& t, L1TrackPar& Tout, const bool dir, const fvec qp0,
+ const bool initParams)
{
L1_assert(t.NHits >= 3);
- L1Fit fit;
+ L1TrackParFit fit;
fit.SetParticleMass(GetDefaultParticleMass());
+ fit.fTr = Tout;
+ L1TrackPar& T = fit.fTr;
// get hits of current track
const L1Vector<L1HitIndex_t>& hits = t.fHits; // array of indeses of hits of current track
@@ -77,6 +77,7 @@ void L1Algo::BranchFitterFast(const L1Branch& t, L1TrackPar& T, const bool dir,
T.ty = (y1 - y0) * dzi;
T.qp = qp0;
}
+ fit.fQp0 = qp0;
T.z = z0;
T.t = hit0.t;
@@ -93,7 +94,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.dt2;
+ T.C55 = sta0.timeInfo ? hit0.dt2 : vINF;
L1FieldValue fldB0, fldB1, fldB2 _fvecalignment;
L1FieldRegion fld _fvecalignment;
@@ -118,38 +119,24 @@ void L1Algo::BranchFitterFast(const L1Branch& t, L1TrackPar& T, const bool dir,
const L1Station& sta = fParameters.GetStation(ista);
- float z_sta = hit.z;
-
- fvec dz = z_sta - T.z;
-
- if (kMcbm == fTrackingMode) { L1ExtrapolateLine(T, hit.z); }
- else {
- L1Extrapolate(T, hit.z, qp0, fld);
- }
- L1ExtrapolateTime(T, dz);
- fit.L1AddMaterial(T, fParameters.GetMaterialThickness(ista, T.x, T.y), qp0, fvec::One());
- //if ((step * ista <= step * (fNstationsBeforePipe + (step + 1) / 2 - 1))
- // && (step * ista_prev >= step * (fNstationsBeforePipe + (step + 1) / 2 - 1 - step)))
- //fit.L1AddPipeMaterial(T, qp0, 1);
-
- fvec u = hit.u;
- fvec v = hit.v;
-
- L1Filter(T, sta.frontInfo, u, hit.du2, fvec::One());
- L1Filter(T, sta.backInfo, v, hit.dv2, fvec::One());
- FilterTime(T, hit.t, hit.dt2);
+ fit.Extrapolate(hit.z, fit.fQp0, 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.AddMsInMaterial(fParameters.GetMaterialThickness(ista, T.x, T.y), fit.fQp0, fvec::One());
fldB0 = fldB1;
fldB1 = fldB2;
fldZ0 = fldZ1;
fldZ1 = fldZ2;
- auto [x, y] = sta.ConvUVtoXY<fvec>(u, v);
+ auto [x, y] = sta.ConvUVtoXY<fvec>(hit.u, hit.v);
sta.fieldSlice.GetFieldValue(x, y, fldB2);
fldZ2 = sta.fZ;
fld.Set(fldB2, fldZ2, fldB1, fldZ1, fldB0, fldZ0);
} // i
+ Tout = T;
} // void L1Algo::BranchFitterFast
/// like BranchFitterFast but more precise
@@ -168,13 +155,16 @@ void L1Algo::BranchFitter(const L1Branch& t, L1TrackPar& T, const bool dir, cons
/// dir - 0 - forward, 1 - backward
/// qp0 - momentum for extrapolation
/// initialize - should be params ititialized. 1 - yes.
-void L1Algo::FindMoreHits(L1Branch& t, L1TrackPar& T, const bool dir,
+void L1Algo::FindMoreHits(L1Branch& t, L1TrackPar& Tout, const bool dir,
const fvec qp0) // TODO take into account pipe
{
L1Vector<L1HitIndex_t> newHits {"L1TrackExtender::newHits"};
newHits.reserve(fParameters.GetNstationsActive());
- L1Fit fit;
+
+ L1TrackParFit fit;
fit.SetParticleMass(GetDefaultParticleMass());
+ fit.fQp0 = qp0;
+ fit.fTr = Tout;
const signed short int step = -2 * static_cast<int>(dir) + 1; // increment for station index
const int iFirstHit = (dir) ? 2 : t.NHits - 3;
@@ -226,22 +216,17 @@ void L1Algo::FindMoreHits(L1Branch& t, L1TrackPar& T, const bool dir,
const L1Station& sta = fParameters.GetStation(ista);
- fvec dz = sta.fZ - T.z;
-
- if (kMcbm == fTrackingMode) { L1ExtrapolateLine(T, sta.fZ); }
- else {
- L1Extrapolate(T, sta.fZ, qp0, fld);
- }
- L1ExtrapolateTime(T, dz);
+ fit.Extrapolate(sta.fZ, fit.fQp0, fld, fvec::One());
fscal r2_best = 1e8; // best distance to hit
int iHit_best = -1; // index of the best hit
- const fscal iz = 1.f / (T.z[0] - fParameters.GetTargetPositionZ()[0]);
+ const fscal iz = 1.f / (fit.fTr.z[0] - fParameters.GetTargetPositionZ()[0]);
- L1HitAreaTime area(vGridTime[ista], T.x[0] * iz, T.y[0] * iz,
- (sqrt(fPickGather * T.C00) + fMaxDx[ista] + fMaxDZ * abs(T.tx))[0] * iz,
- (sqrt(fPickGather * T.C11) + fMaxDy[ista] + fMaxDZ * abs(T.ty))[0] * iz, T.t[0], sqrt(T.C55[0]));
+ L1HitAreaTime area(vGridTime[ista], fit.fTr.x[0] * iz, fit.fTr.y[0] * iz,
+ (sqrt(fPickGather * fit.fTr.C00) + fMaxDx[ista] + fMaxDZ * abs(fit.fTr.tx))[0] * iz,
+ (sqrt(fPickGather * fit.fTr.C11) + fMaxDy[ista] + fMaxDZ * abs(fit.fTr.ty))[0] * iz,
+ fit.fTr.t[0], sqrt(fit.fTr.C55[0]));
for (L1HitIndex_t ih = -1; true;) { // loop over the hits in the area
@@ -260,7 +245,10 @@ void L1Algo::FindMoreHits(L1Branch& t, L1TrackPar& T, const bool dir,
const L1Hit& hit = (*vHitsUnused)[globalInd];
- if (fabs(hit.t - T.t[0]) > sqrt(T.C55[0] + hit.dt2) * 5) continue;
+ if (sta.timeInfo) {
+ fscal dt = hit.t - fit.fTr.t[0];
+ if (dt * dt > (fit.fTr.C55[0] + hit.dt2) * 25) continue;
+ }
//if (GetFUsed((*fStripFlag)[hit.f] | (*fStripFlag)[hit.b])) continue; // if used
//L1_SHOW(fvHitKeyFlags.size());
@@ -272,14 +260,14 @@ void L1Algo::FindMoreHits(L1Branch& t, L1TrackPar& T, const bool dir,
GetHitCoor(hit, xh, yh, zh, sta); // faster
fvec y, C11;
- L1ExtrapolateYC11Line(T, zh, y, C11);
+ fit.ExtrapolateYC11Line(zh, y, C11);
// fscal dym_est = ( fPickGather * sqrt(fabs(C11[0]+sta.XYInfo.C11[0])) );
// fscal y_minus_new = y[0] - dym_est;
// if (yh < y_minus_new) continue; // CHECKME take into account overlaping?
fvec x, C00;
- L1ExtrapolateXC00Line(T, zh, x, C00);
+ fit.ExtrapolateXC00Line(zh, x, C00);
fscal d_x = xh - x[0];
fscal d_y = yh - y[0];
@@ -297,21 +285,14 @@ void L1Algo::FindMoreHits(L1Branch& t, L1TrackPar& T, const bool dir,
newHits.push_back((*RealIHitP)[iHit_best]);
const L1Hit& hit = (*vHitsUnused)[iHit_best];
- fvec u = hit.u;
- fvec v = hit.v;
- fvec z = hit.z;
- auto [x, y] = sta.ConvUVtoXY<fvec>(u, v);
- fvec dz1 = z - T.z;
+ auto [x, y] = sta.ConvUVtoXY<fvec>(hit.u, hit.v);
- L1ExtrapolateTime(T, dz1);
-
- L1ExtrapolateLine(T, z);
- fit.L1AddMaterial(T, fParameters.GetMaterialThickness(ista, T.x, T.y), qp0, fvec::One());
-
- L1Filter(T, sta.frontInfo, u, hit.du2, fvec::One());
- L1Filter(T, sta.backInfo, v, hit.dv2, fvec::One());
- FilterTime(T, hit.t, hit.dt2);
+ fit.Extrapolate(hit.z, fit.fQp0, 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.AddMsInMaterial(fParameters.GetMaterialThickness(ista, fit.fTr.x, fit.fTr.y), fit.fQp0, fvec::One());
fldB0 = fldB1;
fldB1 = fldB2;
@@ -342,6 +323,8 @@ void L1Algo::FindMoreHits(L1Branch& t, L1TrackPar& T, const bool dir,
}
}
+ Tout = fit.fTr;
+
} // void L1Algo::FindMoreHits
/// Try to extrapolate and find additional hits on other stations
diff --git a/reco/L1/L1Algo/L1CATrackFinder.cxx b/reco/L1/L1Algo/L1CATrackFinder.cxx
index eada2696ef..868a7354de 100644
--- a/reco/L1/L1Algo/L1CATrackFinder.cxx
+++ b/reco/L1/L1Algo/L1CATrackFinder.cxx
@@ -25,15 +25,13 @@
#include "L1Algo.h"
#include "L1Assert.h"
#include "L1Branch.h"
-#include "L1Extrapolation.h"
-#include "L1Filtration.h"
-#include "L1Fit.h"
#include "L1Grid.h"
#include "L1HitArea.h"
#include "L1HitPoint.h"
#include "L1Portion.h"
#include "L1Track.h"
#include "L1TrackPar.h"
+#include "L1TrackParFit.h"
#ifdef _OPENMP
#include "omp.h"
@@ -214,6 +212,8 @@ inline void L1Algo::findSingletsStep1( /// input 1st stage of singlet search
const L1Station& fld1Sta2 = fParameters.GetStation(fld1Ista2);
L1TrackParFit fit;
+ fit.SetParticleMass(GetDefaultParticleMass());
+
fit.fQp0 = fvec(0.);
if (fpCurrentIteration->GetElectronFlag()) { fit.SetParticleMass(L1Constants::phys::kElectronMass); }
@@ -264,29 +264,31 @@ inline void L1Algo::findSingletsStep1( /// input 1st stage of singlet search
// TODO: iteration parameter: "Starting NDF of track parameters"
// NDF = number of track parameters (6: x, y, tx, ty, qp, time) - number of measured parameters (3: x, y, time) on station or (2: x, y) on target
- // Alternative: Iteration can find tracks starting from target or from station: -> use a FLAG
- T.tx = tx;
- T.ty = ty;
- T.t = time;
+ T.x = xl;
+ T.y = yl;
+ T.z = zl;
+ T.tx = tx;
+ T.ty = ty;
+ T.qp = fvec(0.);
+ T.t = time;
+ fit.fQp0 = fvec(0.);
- T.qp = fvec(0.);
T.C20 = T.C21 = fvec(0.);
T.C30 = T.C31 = T.C32 = fvec(0.);
T.C40 = T.C41 = T.C42 = T.C43 = fvec(0.);
T.C50 = T.C51 = T.C52 = T.C53 = T.C54 = fvec(0.);
+
T.C22 = T.C33 = fMaxSlopePV * fMaxSlopePV / fvec(9.);
+
if (kGlobal == fTrackingMode || kMcbm == fTrackingMode) T.C22 = T.C33 = fvec(10.);
- // TODO: Why 9. and 10.?
+ // TODO: Why 9. and 10.?
T.C44 = fMaxInvMom / fvec(3.) * fMaxInvMom / fvec(3.);
T.C55 = timeEr2;
{ // add the target constraint
- T.x = xl;
- T.y = yl;
- T.z = zl;
std::tie(T.C00, T.C10, T.C11) = stal.FormXYCovarianceMatrix(du2_l[i1_V], dv2_l[i1_V]);
@@ -295,15 +297,10 @@ inline void L1Algo::findSingletsStep1( /// input 1st stage of singlet search
fit.AddMsInMaterial(fParameters.GetMaterialThickness(istal, T.x, T.y), fMaxInvMom, fvec::One());
- //if ((istam >= fNstationsBeforePipe) && (istal <= fNstationsBeforePipe - 1)) {
- //fit.L1AddPipeMaterial(T, fMaxInvMom, fvec::One());
- //}
+ // extrapolate to the middle hit
- fvec dz = stam.fZ - zl;
- L1ExtrapolateTime(T, dz, fvec::One());
+ fit.ExtrapolateLine(stam.fZ, fld1, fvec::One());
- // extrapolate to the middle hit
- L1Extrapolate0(T, stam.fZ, fld0);
T_1[i1_V] = T;
} // i1_V
}
@@ -329,12 +326,16 @@ inline void L1Algo::findDoubletsStep0(
n2 = 0; // number of doublets
+ L1TrackParFit fit;
+ fit.SetParticleMass(GetDefaultParticleMass());
+
for (Tindex i1 = 0; i1 < n1; ++i1) // for each singlet
{
unsigned int Ndoublets = 0;
const Tindex i1_V = i1 / fvec::size();
const Tindex i1_4 = i1 % fvec::size();
- L1TrackPar& T1 = T_1[i1_V];
+ fit.fTr = T_1[i1_V];
+ L1TrackPar& T1 = fit.fTr;
// assert(T1.IsEntryConsistent(true, i1_4));
// if (!T1.IsEntryConsistent(false, i1_4)) continue;
@@ -397,7 +398,8 @@ inline void L1Algo::findDoubletsStep0(
// if ( dt*dt > dt_est2 && dt < 0 ) continue;
fvec y, C11;
- L1ExtrapolateYC11Line(T1, zm, y, C11);
+ fit.ExtrapolateYC11Line(zm, y, C11);
+ //L1ExtrapolateYC11Line(T1, zm, y, C11);
/// Covariation matrix of the hit
auto [dxxScalMhit, dxyScalMhit, dyyScalMhit] = stam.FormXYCovarianceMatrix(hitm.dU2(), hitm.dV2());
@@ -412,7 +414,8 @@ inline void L1Algo::findDoubletsStep0(
// check x-boundaries
fvec x, C00;
- L1ExtrapolateXC00Line(T1, zm, x, C00);
+ fit.ExtrapolateXC00Line(zm, x, C00);
+ //L1ExtrapolateXC00Line(T1, zm, x, C00);
fscal dx_est2 = Pick_m22[i1_4] * fabs(C00[i1_4] + dxxScalMhit);
@@ -422,16 +425,18 @@ inline void L1Algo::findDoubletsStep0(
// check chi2
fvec C10;
- L1ExtrapolateC10Line(T1, zm, C10);
+ fit.ExtrapolateC10Line(zm, C10);
+ //L1ExtrapolateC10Line(T1, zm, C10);
+
fvec chi2 = T1.chi2;
- L1FilterChi2XYC00C10C11(stam.frontInfo, x, y, C00, C10, C11, chi2, hitm.U(), hitm.dU2());
+ L1TrackParFit::FilterChi2XYC00C10C11(stam.frontInfo, x, y, C00, C10, C11, chi2, hitm.U(), hitm.dU2());
if (!fpCurrentIteration->GetTrackFromTripletsFlag()) {
if (chi2[i1_4] > fDoubletChi2Cut) continue;
}
- L1FilterChi2(stam.backInfo, x, y, C00, C10, C11, chi2, hitm.V(), hitm.dV2());
+ L1TrackParFit::FilterChi2(stam.backInfo, x, y, C00, C10, C11, chi2, hitm.V(), hitm.dV2());
// FilterTime(T1, hitm.T(), hitm.dT2());
@@ -501,6 +506,8 @@ inline void L1Algo::findDoubletsStep0(
}
} // loop over the hits in the area
+ T_1[i1_V] = fit.fTr;
+
} // for i1
}
@@ -561,6 +568,8 @@ inline void L1Algo::findTripletsStep0( // input
for (Tindex i2 = 0; i2 < n2;) {
L1TrackPar& T2 = fit.fTr;
T2 = L1TrackPar_0;
+ fit.fQp0 = fvec(0.);
+
L1FieldRegion f2;
// pack the data
fvec u_front_2 = 0.f;
@@ -608,20 +617,19 @@ inline void L1Algo::findTripletsStep0( // input
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);
+ fit.FilterTime(t_2, dt2_2, fvec::One(), stam.timeInfo);
+
+ fit.fQp0 = fit.fTr.qp;
- fit.AddMsInMaterial(fParameters.GetMaterialThickness(iStaM, T2.x, T2.y), isMomentumFitted ? T2.qp : fMaxInvMom,
+ fit.AddMsInMaterial(fParameters.GetMaterialThickness(iStaM, T2.x, T2.y), isMomentumFitted ? fit.fQp0 : fMaxInvMom,
fvec::One());
+
//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);
+ fit.Extrapolate(star.fZ, fit.fQp0, f2, fvec::One());
// assert(T2.IsConsistent(true, n2_4));
@@ -697,7 +705,8 @@ inline void L1Algo::findTripletsStep0( // input
// - check whether hit belong to the window ( track position +\- errors ) -
// check lower boundary
fvec y, C11;
- L1ExtrapolateYC11Line(T2, zr, y, C11);
+ fit3.ExtrapolateYC11Line(zr, y, C11);
+ //L1ExtrapolateYC11Line(T2, zr, y, C11);
/// Covariation matrix of the hit
auto [dxxScalRhit, dxyScalRhit, dyyScalRhit] = star.FormXYCovarianceMatrix(hitr.dU2(), hitr.dV2());
@@ -710,8 +719,7 @@ inline void L1Algo::findTripletsStep0( // input
if (dY2 > dy_est2) continue; // if (yr > y_plus_new [i2_4] ) continue;
// check x-boundaries
fvec x, C00;
-
- L1ExtrapolateXC00Line(T2, zr, x, C00);
+ fit3.ExtrapolateXC00Line(zr, x, C00);
fscal dx_est2 = (Pick_r22[i2_4] * (fabs(C00[i2_4] + dxxScalRhit)));
@@ -719,15 +727,15 @@ inline void L1Algo::findTripletsStep0( // input
if (dX * dX > dx_est2) continue;
// check chi2 // not effective
fvec C10;
- L1ExtrapolateC10Line(T2, zr, C10);
- fvec chi2 = T2.chi2;
+ fit3.ExtrapolateC10Line(zr, C10);
- L1FilterChi2XYC00C10C11(star.frontInfo, x, y, C00, C10, C11, chi2, hitr.U(), hitr.dU2());
+ fvec chi2 = T2.chi2;
- L1FilterChi2(star.backInfo, x, y, C00, C10, C11, chi2, hitr.V(), hitr.dV2());
+ L1TrackParFit::FilterChi2XYC00C10C11(star.frontInfo, x, y, C00, C10, C11, chi2, hitr.U(), hitr.dU2());
- FilterTime(T_cur, hitr.T(), hitr.dT2(), star.timeInfo);
+ L1TrackParFit::FilterChi2(star.backInfo, x, y, C00, C10, C11, chi2, hitr.V(), hitr.dV2());
+ fit3.FilterTime(hitr.T(), hitr.dT2(), fvec::One(), star.timeInfo);
if (!fpCurrentIteration->GetTrackFromTripletsFlag()) {
if (chi2[i2_4] > fTripletChi2Cut || C00[i2_4] < 0 || C11[i2_4] < 0 || T_cur.C55[i2_4] < 0) {
diff --git a/reco/L1/L1Algo/L1CloneMerger.cxx b/reco/L1/L1Algo/L1CloneMerger.cxx
index 9d2d28b4ad..219410cae2 100644
--- a/reco/L1/L1Algo/L1CloneMerger.cxx
+++ b/reco/L1/L1Algo/L1CloneMerger.cxx
@@ -12,9 +12,9 @@
#include <iostream>
#include "L1Algo.h"
-#include "L1Extrapolation.h"
#include "L1Parameters.h"
#include "L1Track.h"
+#include "L1TrackParFit.h"
// ---------------------------------------------------------------------------------------------------------------------
//
@@ -76,8 +76,16 @@ void L1CloneMerger::Exec(L1Vector<L1Track>& extTracks, L1Vector<L1HitIndex_t>& e
isDownstreamNeighbour[iTr] = false;
}
- L1TrackPar Tb;
- L1TrackPar Tf;
+ L1TrackParFit fitB;
+ fitB.SetParticleMass(frAlgo.GetDefaultParticleMass());
+ fitB.fQp0 = fvec(0.);
+
+ L1TrackParFit fitF;
+ fitF.SetParticleMass(frAlgo.GetDefaultParticleMass());
+ fitF.fQp0 = fvec(0.);
+
+ L1TrackPar& Tb = fitB.fTr;
+ L1TrackPar& Tf = fitF.fTr;
L1FieldValue fBm, fBb, fBf _fvecalignment;
L1FieldRegion fld _fvecalignment;
@@ -125,6 +133,8 @@ void L1CloneMerger::Exec(L1Vector<L1Track>& extTracks, L1Vector<L1HitIndex_t>& e
Tb.C54 = extTracks[iTr].CFirst[19];
Tb.C55 = extTracks[iTr].CFirst[20];
+ fitB.fQp0 = fitB.fTr.qp;
+
unsigned short staf = lastStation[jTr];
Tf.x = extTracks[jTr].TLast[0];
@@ -156,6 +166,7 @@ void L1CloneMerger::Exec(L1Vector<L1Track>& extTracks, L1Vector<L1HitIndex_t>& e
Tf.C54 = extTracks[jTr].CLast[19];
Tf.C55 = extTracks[jTr].CLast[20];
+ fitF.fQp0 = fitF.fTr.qp;
if (fabs(Tf.t[0] - Tb.t[0]) > 3 * sqrt(Tf.C55[0] + Tb.C55[0])) continue;
unsigned short stam;
@@ -177,8 +188,8 @@ void L1CloneMerger::Exec(L1Vector<L1Track>& extTracks, L1Vector<L1HitIndex_t>& e
fvec zMiddle = fvec(0.5) * (Tb.z + Tf.z);
- L1Extrapolate(Tf, zMiddle, Tf.qp, fld);
- L1Extrapolate(Tb, zMiddle, Tb.qp, fld);
+ fitF.Extrapolate(zMiddle, fitF.fQp0, fld, fvec::One());
+ fitB.Extrapolate(zMiddle, fitB.fQp0, fld, fvec::One());
fvec Chi2Tracks(0.);
FilterTracks(&(Tf.x), &(Tf.C00), &(Tb.x), &(Tb.C00), 0, 0, &Chi2Tracks);
diff --git a/reco/L1/L1Algo/L1TrackFitter.cxx b/reco/L1/L1Algo/L1TrackFitter.cxx
index e1b4f6bbed..3926955a1f 100644
--- a/reco/L1/L1Algo/L1TrackFitter.cxx
+++ b/reco/L1/L1Algo/L1TrackFitter.cxx
@@ -6,9 +6,6 @@
#include <vector>
#include "L1Algo.h"
-#include "L1Extrapolation.h"
-#include "L1Filtration.h" // for KFTrackFitter_simple
-#include "L1Fit.h"
#include "L1TrackPar.h"
#include "L1TrackParFit.h"
@@ -26,285 +23,6 @@ namespace NS_L1TrackFitter
} // namespace NS_L1TrackFitter
using namespace NS_L1TrackFitter;
-/// Fit reconstracted track like it fitted d_uring the reconstruction.
-void L1Algo::KFTrackFitter_simple() // TODO: Add pipe.
-{
- // cout << " Start KF Track Fitter " << endl;
- int start_hit = 0; // for interation in fRecoHits[]
- for (unsigned int itrack = 0; itrack < fTracks.size(); itrack++) {
- L1Track& t = fTracks[itrack]; // current track
-
- // get hits of current track
- std::vector<unsigned short int> hits; // array of indeses of hits of current track
- hits.clear();
- int nHits = t.NHits;
- for (int i = 0; i < nHits; i++) {
- hits.push_back(fRecoHits[start_hit++]);
- }
-
- L1TrackPar T; // fitting parametr coresponding to current track
-
- L1Fit fit;
- fit.SetParticleMass(GetDefaultParticleMass());
-
- fvec qp0(0.25);
- //fvec qp0 = 2./t.Momentum;
- for (int iter = 0; iter < 3; iter++) {
- //cout<<" Back 1"<<endl;
- { // fit backward
- const L1Hit& hit0 = fInputData.GetHit(hits[nHits - 1]);
- const L1Hit& hit1 = fInputData.GetHit(hits[nHits - 2]);
- const L1Hit& hit2 = fInputData.GetHit(hits[nHits - 3]);
-
- int ista0 = hit0.iSt;
- int ista1 = hit1.iSt;
- int ista2 = hit2.iSt;
-
-
- //cout<<"back: ista012="<<ista0<<" "<<ista1<<" "<<ista2<<endl;
- const L1Station& sta0 = fParameters.GetStation(ista0);
- const L1Station& sta1 = fParameters.GetStation(ista1);
- const L1Station& sta2 = fParameters.GetStation(ista2);
-
- fvec u0 = hit0.u;
- fvec v0 = hit0.v;
- auto [x0, y0] = sta0.ConvUVtoXY<fvec>(u0, v0);
- fvec z0 = hit0.z;
-
- fvec u1 = hit1.u;
- fvec v1 = hit1.v;
- auto [x1, y1] = sta1.ConvUVtoXY<fvec>(u1, v1);
- fvec z1 = hit1.z;
-
- fvec u2 = hit2.u;
- fvec v2 = hit2.v;
- auto [x2, y2] = sta1.ConvUVtoXY<fvec>(u2, v2);
- // fvec z2 = hit2.z;
-
- fvec dzi = fvec(1.) / (z1 - z0);
-
- // const fvec vINF = .1;
- T.x = x0;
- T.y = y0;
- if (iter == 0) {
- T.tx = (x1 - x0) * dzi;
- T.ty = (y1 - y0) * dzi;
- }
-
- T.qp = qp0;
- T.z = z0;
- T.chi2 = fvec(0.);
- T.NDF = fvec(2.);
-
- std::tie(T.C00, T.C10, T.C11) = sta0.FormXYCovarianceMatrix(hit0.du2, hit0.dv2);
-
- T.C20 = T.C21 = fvec(0.);
- T.C30 = T.C31 = T.C32 = fvec(0.);
- T.C40 = T.C41 = T.C42 = T.C43 = fvec(0.);
- T.C22 = T.C33 = vINF;
- T.C44 = fvec(1.);
-
- // static L1FieldValue fldB0, fldB1, fldB2 _fvecalignment;
- // static L1FieldRegion fld _fvecalignment;
- L1FieldValue fldB0, fldB1, fldB2 _fvecalignment;
- L1FieldRegion fld _fvecalignment;
- fvec fldZ0 = sta1.fZ; // suppose field is smoth
- fvec fldZ1 = sta2.fZ;
- fvec fldZ2 = sta0.fZ;
-
- sta1.fieldSlice.GetFieldValue(x1, y1, fldB0);
- sta2.fieldSlice.GetFieldValue(x2, y2, fldB1);
- sta0.fieldSlice.GetFieldValue(x0, y0, fldB2);
-
- fld.Set(fldB2, fldZ2, fldB1, fldZ1, fldB0, fldZ0);
-
- int ista = ista2;
- //cout<<"\nfit, iter=:"<<iter<<endl;
- for (int i = nHits - 2; i >= 0; i--) {
- // if( fabs(T.qp[0])>2. ) break; // iklm. Don't know it need for
- const L1Hit& hit = fInputData.GetHit(hits[i]);
- ista = hit.iSt;
-
- const L1Station& sta = fParameters.GetStation(ista);
-
- // L1Extrapolate( T, hit.z, qp0, fld );
- L1ExtrapolateLine(T, hit.z);
- // T.L1Extrapolate( sta.z, qp0, fld );
- // L1Extrapolate( T, hit.z, qp0, fld );
- fit.L1AddMaterial(T, fParameters.GetMaterialThickness(i, T.x, T.y), qp0, ONE);
-
- // if (ista == fNstationsBeforePipe - 1) fit.L1AddPipeMaterial( T, qp0, 1);
-
- fvec u = hit.u;
- fvec v = hit.v;
- auto [x, y] = sta.ConvUVtoXY<fvec>(u, v);
- L1Filter(T, sta.frontInfo, u, hit.du2, fvec::One());
- L1Filter(T, sta.backInfo, v, hit.dv2, fvec::One());
- fldB0 = fldB1;
- fldB1 = fldB2;
- fldZ0 = fldZ1;
- fldZ1 = fldZ2;
- sta.fieldSlice.GetFieldValue(x, y, fldB2);
-
- fldZ2 = sta.fZ;
- fld.Set(fldB2, fldZ2, fldB1, fldZ1, fldB0, fldZ0);
- } // i
-
- // write received parametres in track
- t.TFirst[0] = T.x[0];
- t.TFirst[1] = T.y[0];
- t.TFirst[2] = T.tx[0];
- t.TFirst[3] = T.ty[0];
- t.TFirst[4] = T.qp[0];
- t.TFirst[5] = T.z[0];
-
- t.CFirst[0] = T.C00[0];
- t.CFirst[1] = T.C10[0];
- t.CFirst[2] = T.C11[0];
- t.CFirst[3] = T.C20[0];
- t.CFirst[4] = T.C21[0];
- t.CFirst[5] = T.C22[0];
- t.CFirst[6] = T.C30[0];
- t.CFirst[7] = T.C31[0];
- t.CFirst[8] = T.C32[0];
- t.CFirst[9] = T.C33[0];
- t.CFirst[10] = T.C40[0];
- t.CFirst[11] = T.C41[0];
- t.CFirst[12] = T.C42[0];
- t.CFirst[13] = T.C43[0];
- t.CFirst[14] = T.C44[0];
-
- t.chi2 = T.chi2[0];
- t.NDF = static_cast<int>(T.NDF[0]);
- qp0 = T.qp;
- } // fit backward
-
- // fit forward
- {
- //T.qp = first_trip->GetQpOrig(MaxInvMom);
- const L1Hit& hit0 = fInputData.GetHit(hits[0]);
- const L1Hit& hit1 = fInputData.GetHit(hits[1]);
- const L1Hit& hit2 = fInputData.GetHit(hits[2]);
-
- int ista0 = hit0.iSt;
- int ista1 = hit1.iSt;
- int ista2 = hit2.iSt;
-
- const L1Station& sta0 = fParameters.GetStation(ista0);
- const L1Station& sta1 = fParameters.GetStation(ista1);
- const L1Station& sta2 = fParameters.GetStation(ista2);
-
- fvec u0 = hit0.u;
- fvec v0 = hit0.v;
- auto [x0, y0] = sta0.ConvUVtoXY<fvec>(u0, v0);
- fvec z0 = hit0.z;
-
- fvec u1 = hit1.u;
- fvec v1 = hit1.v;
- auto [x1, y1] = sta1.ConvUVtoXY<fvec>(u1, v1);
- // fvec z1 = hit1.z;
-
- fvec u2 = hit2.u;
- fvec v2 = hit2.v;
- auto [x2, y2] = sta2.ConvUVtoXY<fvec>(u2, v2);
- // fvec z2 = hit2.z;
-
- // fvec dzi = 1./(z1-z0);
-
- //fvec qp0 = first_trip->GetQpOrig(MaxInvMom);
-
- // const fvec vINF = .1;
- T.chi2 = fvec(0.);
- T.NDF = fvec(2.);
- T.x = x0;
- T.y = y0;
- // T.tx = (x1-x0)*dzi;
- // T.ty = (y1-y0)*dzi;
- // qp0 = 0;
- T.qp = qp0;
- T.z = z0;
-
- 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;
- T.C40 = T.C41 = T.C42 = T.C43 = 0;
- T.C22 = T.C33 = vINF;
- T.C44 = 1.;
-
- // static L1FieldValue fldB0, fldB1, fldB2 _fvecalignment;
- // static L1FieldRegion fld _fvecalignment;
- L1FieldValue fldB0, fldB1, fldB2 _fvecalignment;
- L1FieldRegion fld _fvecalignment;
- fvec fldZ0 = sta1.fZ;
- fvec fldZ1 = sta2.fZ;
- fvec fldZ2 = sta0.fZ;
-
- sta1.fieldSlice.GetFieldValue(x1, y1, fldB0);
- sta2.fieldSlice.GetFieldValue(x2, y2, fldB1);
- sta0.fieldSlice.GetFieldValue(x0, y0, fldB2);
-
- fld.Set(fldB2, fldZ2, fldB1, fldZ1, fldB0, fldZ0);
- int ista = ista2;
-
- for (int i = 1; i < nHits; i++) {
- const L1Hit& hit = fInputData.GetHit(hits[i]);
- ista = hit.iSt;
- const L1Station& sta = fParameters.GetStation(ista);
- fvec u = hit.u;
- fvec v = hit.v;
- auto [x, y] = sta.ConvUVtoXY<fvec>(u, v);
-
- // L1Extrapolate( T, hit.z, qp0, fld );
- L1ExtrapolateLine(T, hit.z);
- // T.L1Extrapolate( sta.z, qp0, fld );
- // L1Extrapolate( T, hit.z, qp0, fld );
- fit.L1AddMaterial(T, fParameters.GetMaterialThickness(i, T.x, T.y), qp0, ONE);
-
- // if (ista == fNstationsBeforePipe) fit.L1AddPipeMaterial( T, qp0, 1);
- L1Filter(T, sta.frontInfo, u, hit.du2, fvec::One());
- L1Filter(T, sta.backInfo, v, hit.dv2, fvec::One());
-
- fldB0 = fldB1;
- fldB1 = fldB2;
- fldZ0 = fldZ1;
- fldZ1 = fldZ2;
- sta.fieldSlice.GetFieldValue(x, y, fldB2);
- fldZ2 = sta.fZ;
- fld.Set(fldB2, fldZ2, fldB1, fldZ1, fldB0, fldZ0);
- }
-
- // write received parametres in track
- t.TLast[0] = T.x[0];
- t.TLast[1] = T.y[0];
- t.TLast[2] = T.tx[0];
- t.TLast[3] = T.ty[0];
- t.TLast[4] = T.qp[0];
- t.TLast[5] = T.z[0];
-
- t.CLast[0] = T.C00[0];
- t.CLast[1] = T.C10[0];
- t.CLast[2] = T.C11[0];
- t.CLast[3] = T.C20[0];
- t.CLast[4] = T.C21[0];
- t.CLast[5] = T.C22[0];
- t.CLast[6] = T.C30[0];
- t.CLast[7] = T.C31[0];
- t.CLast[8] = T.C32[0];
- t.CLast[9] = T.C33[0];
- t.CLast[10] = T.C40[0];
- t.CLast[11] = T.C41[0];
- t.CLast[12] = T.C42[0];
- t.CLast[13] = T.C43[0];
- t.CLast[14] = T.C44[0];
-
- t.chi2 += T.chi2[0];
- t.NDF += static_cast<int>(T.NDF[0]);
- }
- qp0 = T.qp;
- }
- } // for(int itrack
-}
void L1Algo::L1KFTrackFitter()
{
diff --git a/reco/L1/L1Algo/L1TrackPar.h b/reco/L1/L1Algo/L1TrackPar.h
index 4f129149e1..261944c4bc 100644
--- a/reco/L1/L1Algo/L1TrackPar.h
+++ b/reco/L1/L1Algo/L1TrackPar.h
@@ -13,72 +13,52 @@
class L1TrackPar {
public:
- fvec x, y, tx, ty, qp, z, t, C00, C10, C11, C20, C21, C22, C30, C31, C32, C33, C40, C41, C42, C43, C44, C50, C51, C52,
- C53, C54, C55, chi2, NDF;
-
- L1TrackPar()
- : x(0)
- , y(0)
- , tx(0)
- , ty(0)
- , qp(0)
- , z(0)
- , t(0)
- , C00(0)
- , C10(0)
- , C11(0)
- , C20(0)
- , C21(0)
- , C22(0)
- , C30(0)
- , C31(0)
- , C32(0)
- , C33(0)
- , C40(0)
- , C41(0)
- , C42(0)
- , C43(0)
- , C44(0)
- , C50(0)
- , C51(0)
- , C52(0)
- , C53(0)
- , C54(0)
- , C55(0)
- , chi2(0)
- , NDF(0) {};
-
- L1TrackPar(double* T, double* C)
- : x(T[0])
- , y(T[1])
- , tx(T[2])
- , ty(T[3])
- , qp(T[4])
- , z(T[5])
- , t(T[6])
- , C00(C[0])
- , C10(C[1])
- , C11(C[2])
- , C20(C[3])
- , C21(C[4])
- , C22(C[5])
- , C30(C[6])
- , C31(C[7])
- , C32(C[8])
- , C33(C[9])
- , C40(C[10])
- , C41(C[11])
- , C42(C[12])
- , C43(C[13])
- , C44(C[14])
- , C50(C[15])
- , C51(C[16])
- , C52(C[17])
- , C53(C[18])
- , C54(C[19])
- , C55(C[20])
- , chi2(0)
- , NDF(0) {};
+ fvec x {0.}, y {0.}, tx {0.}, ty {0.}, qp {0.}, z {0.}, t {0.};
+
+ fvec C00 {0.}, C10 {0.}, C11 {0.}, C20 {0.}, C21 {0.}, C22 {0.}, C30 {0.}, C31 {0.}, C32 {0.}, C33 {0.}, C40 {0.},
+ C41 {0.}, C42 {0.}, C43 {0.}, C44 {0.}, C50 {0.}, C51 {0.}, C52 {0.}, C53 {0.}, C54 {0.}, C55 {0.}, chi2 {0.},
+ NDF {0.};
+
+ L1TrackPar() = default;
+
+ L1TrackPar(const double* T, const double* C) { Set(T, C); }
+
+ void Set(const double* T, const double* C)
+ {
+ x = T[0];
+ y = T[1];
+ tx = T[2];
+ ty = T[3];
+ qp = T[4];
+ z = T[5];
+ t = T[6];
+
+ chi2 = 0.;
+ NDF = 0.;
+
+ C00 = C[0];
+ C10 = C[1];
+ C11 = C[2];
+ C20 = C[3];
+ C21 = C[4];
+ C22 = C[5];
+ C30 = C[6];
+ C31 = C[7];
+ C32 = C[8];
+ C33 = C[9];
+ C40 = C[10];
+ C41 = C[11];
+ C42 = C[12];
+ C43 = C[13];
+ C44 = C[14];
+ C50 = C[15];
+ C51 = C[16];
+ C52 = C[17];
+ C53 = C[18];
+ C54 = C[19];
+ C55 = C[20];
+ }
+
void SetOneEntry(const int i0, const L1TrackPar& T1, const int i1);
diff --git a/reco/L1/L1Algo/L1TrackParFit.cxx b/reco/L1/L1Algo/L1TrackParFit.cxx
index ef3ba7252e..817b036a2a 100644
--- a/reco/L1/L1Algo/L1TrackParFit.cxx
+++ b/reco/L1/L1Algo/L1TrackParFit.cxx
@@ -4,11 +4,6 @@
#include "L1TrackParFit.h"
-#include "L1Extrapolation.h"
-#include "L1Filtration.h"
-#include "L1Fit.h"
-
-
#define cnst const fvec
void L1TrackParFit::Filter(const L1UMeasurementInfo& info, const fvec& u, const fvec& sigma2, const fvec& w)
@@ -927,6 +922,15 @@ void L1TrackParFit::
//cout<<"Extrapolation ok"<<endl;
}
+void L1TrackParFit::ExtrapolateLine(fvec z_out, const L1FieldRegion& F, const fvec& w)
+{
+ // extrapolate the track assuming fQp0 == 0
+ //
+ auto qp0 = fQp0;
+ fQp0 = fvec(0.);
+ Extrapolate(z_out, fQp0, F, w);
+ fQp0 = qp0;
+}
void L1TrackParFit::ExtrapolateLine(fvec z_out, const fvec& w)
{
@@ -1058,7 +1062,7 @@ void L1TrackParFit::EnergyLossCorrection(const fvec& radThick, fvec& qp0, fvec d
const fvec p2 = fvec(1.) / qp02;
const fvec E2 = fMass2 + p2;
- const fvec bethe = L1Fit::ApproximateBetheBloch(p2 / fMass2);
+ const fvec bethe = ApproximateBetheBloch(p2 / fMass2);
fvec tr = sqrt(fvec(1.f) + fTr.tx * fTr.tx + fTr.ty * fTr.ty);
@@ -1093,7 +1097,7 @@ void L1TrackParFit::EnergyLossCorrection(float atomicA, float rho, float radLen,
else
i = (9.76 * atomicZ + 58.8 * std::pow(atomicZ, -0.19)) * 1.e-9;
- const fvec bethe = L1Fit::ApproximateBetheBloch(p2 / fMass2, rho, 0.20, 3.00, i, atomicZ / atomicA);
+ const fvec 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);
@@ -1230,4 +1234,145 @@ void L1TrackParFit::AddTargetToLine(const fvec& targX, const fvec& targY, const
FilterExtrapolatedXY(targX, targY, targXYInfo, eX, eY, Jx, Jy);
}
+fvec L1TrackParFit::ApproximateBetheBloch(const fvec& bg2)
+{
+ //
+ // This is the parameterization of the Bethe-Bloch formula inspired by Geant.
+ //
+ // bg2 - (beta*gamma)^2
+ // kp0 - density [g/cm^3]
+ // kp1 - density effect first junction point
+ // kp2 - density effect second junction point
+ // kp3 - mean excitation energy [GeV]
+ // kp4 - mean Z/A
+ //
+ // The default values for the kp* parameters are for silicon.
+ // 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;
+
+ 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
+
+ //*** 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);
+
+ 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);
+
+ return mK * mZA * (fvec(1.f) + bg2) / bg2
+ * (0.5f * log(_2me * bg2 * maxT / (mI * mI)) - bg2 / (fvec(1.f) + bg2) - d2);
+}
+
+fvec L1TrackParFit::ApproximateBetheBloch(const fvec& bg2, const fvec& kp0, const fvec& kp1, const fvec& kp2,
+ const fvec& kp3, const fvec& kp4)
+{
+ //
+ // This is the parameterization of the Bethe-Bloch formula inspired by Geant.
+ //
+ // bg2 - (beta*gamma)^2
+ // kp0 - density [g/cm^3]
+ // kp1 - density effect first junction point
+ // kp2 - density effect second junction point
+ // kp3 - mean excitation energy [GeV]
+ // kp4 - mean Z/A
+ //
+ // The default values for the kp* parameters are for silicon.
+ // 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;
+
+ 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
+
+ //*** 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);
+
+ 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);
+
+ return mK * mZA * (fvec(1.f) + bg2) / bg2
+ * (0.5f * log(_2me * bg2 * maxT / (mI * mI)) - bg2 / (fvec(1.f) + bg2) - d2);
+}
+
+
+void L1TrackParFit::FilterChi2XYC00C10C11(const L1UMeasurementInfo& info, fvec& x, fvec& y, fvec& C00, fvec& C10,
+ fvec& C11, fvec& chi2, const fvec& u, const fvec& du2)
+{
+ fvec wi, zeta, zetawi, HCH;
+ fvec F0, F1;
+ fvec K1;
+
+ zeta = info.cos_phi * x + info.sin_phi * y - u;
+
+ // F = CH'
+ F0 = info.cos_phi * C00 + info.sin_phi * C10;
+ F1 = info.cos_phi * C10 + info.sin_phi * C11;
+
+ HCH = (F0 * info.cos_phi + F1 * info.sin_phi);
+
+ wi = fvec(1.) / (du2 + HCH);
+ zetawi = zeta * wi;
+ chi2 += zeta * zetawi;
+
+ K1 = F1 * wi;
+
+ x -= F0 * zetawi;
+ y -= F1 * zetawi;
+
+ C00 -= F0 * F0 * wi;
+ C10 -= K1 * F0;
+ C11 -= K1 * F1;
+}
+
+
+void L1TrackParFit::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)
+{
+ fvec zeta, HCH;
+ fvec F0, F1;
+
+ zeta = info.cos_phi * x + info.sin_phi * y - u;
+
+ // F = CH'
+ F0 = info.cos_phi * C00 + info.sin_phi * C10;
+ F1 = info.cos_phi * C10 + info.sin_phi * C11;
+
+ HCH = (F0 * info.cos_phi + F1 * info.sin_phi);
+
+ chi2 += zeta * zeta / (du2 + HCH);
+}
+
#undef cnst
diff --git a/reco/L1/L1Algo/L1TrackParFit.h b/reco/L1/L1Algo/L1TrackParFit.h
index 217351c5f2..267f6dad92 100644
--- a/reco/L1/L1Algo/L1TrackParFit.h
+++ b/reco/L1/L1Algo/L1TrackParFit.h
@@ -14,18 +14,23 @@
class L1TrackParFit {
public:
- L1TrackPar fTr {};
- fvec fQp0;
+ L1TrackParFit() = default;
- fvec fMass = 0.10565800; // muon mass
- fvec fMass2 = fMass * fMass; // mass squared
+ L1TrackParFit(const L1TrackPar& t) { SetTrack(t); }
- fvec fPipeRadThick {7.87e-3f}; // 0.7 mm Aluminium // TODO:
- fvec fTargetRadThick {3.73e-2f * 2}; // 250 mum Gold // TODO:
+ L1TrackParFit(const double* T, const double* C) { SetTrack(T, C); }
- L1TrackParFit() = default;
+ void SetTrack(const double* T, const double* C)
+ {
+ fTr.Set(T, C);
+ fQp0 = fTr.qp;
+ }
- L1TrackParFit(double* T, double* C) : fTr(T, C) {}
+ void SetTrack(const L1TrackPar& t)
+ {
+ fTr = t;
+ fQp0 = fTr.qp;
+ }
void SetOneEntry(const int i0, const L1TrackParFit& T1, const int i1);
@@ -33,8 +38,11 @@ public:
//Fit functionality
+ /// get target radiation thickness
+ fvec GetTargetRadThick() const { return fTargetRadThick; }
+
/// set particle mass for the fit
- void SetParticleMass(float mass)
+ void SetParticleMass(fvec mass)
{
fMass = mass;
fMass2 = mass * mass;
@@ -56,6 +64,7 @@ public:
void ExtrapolateStep(fvec z_out, fvec qp0, const L1FieldRegion& F, const fvec& w);
void ExtrapolateStepAnalytic(fvec z_out, fvec qp0, 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 EnergyLossCorrection(const fvec& radThick, fvec& qp0, fvec direction, fvec w);
@@ -74,9 +83,35 @@ public:
void GetExtrapolatedXYline(const fvec& 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 AddTargetToLine(const fvec& targX, const fvec& targY, const fvec& targZ, const L1XYMeasurementInfo& targXYInfo,
const L1FieldRegion& F);
+ static fvec ApproximateBetheBloch(const fvec& bg2);
+
+ static fvec ApproximateBetheBloch(const fvec& bg2, const fvec& kp0, const fvec& kp1, const fvec& kp2, const fvec& kp3,
+ const fvec& 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);
+
+ 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);
+
+public:
+ L1TrackPar fTr {};
+ fvec fQp0;
+
+ fvec fMass = 0.10565800; // muon mass
+ fvec fMass2 = fMass * fMass; // mass squared
+
+ fvec fPipeRadThick {7.87e-3f}; // 0.7 mm Aluminium // TODO:
+ fvec fTargetRadThick {3.73e-2f * 2}; // 250 mum Gold // TODO:
+
} _fvecalignment;
// =============================================================================================
@@ -90,4 +125,23 @@ inline void L1TrackParFit::SetOneEntry(const int i0, const L1TrackParFit& T1, co
fQp0[i0] = T1.fQp0[i1];
}
+inline void L1TrackParFit::ExtrapolateXC00Line(fvec 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);
+}
+
+inline void L1TrackParFit::ExtrapolateYC11Line(fvec 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);
+}
+
+inline void L1TrackParFit::ExtrapolateC10Line(fvec z_out, fvec& C10) const
+{
+ const fvec dz = (z_out - fTr.z);
+ C10 = fTr.C10 + dz * (fTr.C21 + fTr.C30 + dz * fTr.C32);
+}
#endif
diff --git a/reco/L1/ParticleFinder/CbmL1PFFitter.cxx b/reco/L1/ParticleFinder/CbmL1PFFitter.cxx
index c6aefbad11..b5f16dc2d2 100644
--- a/reco/L1/ParticleFinder/CbmL1PFFitter.cxx
+++ b/reco/L1/ParticleFinder/CbmL1PFFitter.cxx
@@ -38,13 +38,12 @@
#include "KFParticleDatabase.h"
#include "L1Algo.h" // Also defines L1Parameters
#include "L1Def.h"
-#include "L1Extrapolation.h"
#include "L1Field.h"
-#include "L1Filtration.h"
-#include "L1Fit.h"
-#include "L1Material.h"
#include "L1Station.h"
#include "L1TrackPar.h"
+#include "L1TrackParFit.h"
+
+//typedef L1Fit1 L1Fit;
using std::vector;
@@ -98,28 +97,29 @@ inline void CbmL1PFFitter::PFFieldRegion::getL1FieldRegion(L1FieldRegion& fld, i
inline CbmL1PFFitter::PFFieldRegion::PFFieldRegion(const L1FieldRegion& fld, int i) { setFromL1FieldRegion(fld, i); }
-void FilterFirst(L1TrackPar& track, fvec& x, fvec& y, fvec& dxx, fvec& dxy, fvec& dyy)
+void FilterFirst(L1TrackParFit& fit, fvec& x, fvec& y, fvec& dxx, fvec& dxy, fvec& dyy)
{
- track.C00 = dxx;
- track.C10 = dxy;
- track.C11 = dyy;
- track.C20 = ZERO;
- track.C21 = ZERO;
- track.C22 = vINF;
- track.C30 = ZERO;
- track.C31 = ZERO;
- track.C32 = ZERO;
- track.C33 = vINF;
- track.C40 = ZERO;
- track.C41 = ZERO;
- track.C42 = ZERO;
- track.C43 = ZERO;
- track.C44 = ONE;
-
- track.x = x;
- track.y = y;
- track.NDF = -3.0;
- track.chi2 = ZERO;
+ L1TrackPar& tr = fit.fTr;
+ tr.C00 = dxx;
+ tr.C10 = dxy;
+ tr.C11 = dyy;
+ tr.C20 = ZERO;
+ tr.C21 = ZERO;
+ tr.C22 = vINF;
+ tr.C30 = ZERO;
+ tr.C31 = ZERO;
+ tr.C32 = ZERO;
+ tr.C33 = vINF;
+ tr.C40 = ZERO;
+ tr.C41 = ZERO;
+ tr.C42 = ZERO;
+ tr.C43 = ZERO;
+ tr.C44 = ONE;
+
+ tr.x = x;
+ tr.y = y;
+ tr.NDF = -3.0;
+ tr.chi2 = ZERO;
}
void CbmL1PFFitter::Fit(vector<CbmStsTrack>& Tracks, vector<int>& pidHypo)
@@ -138,12 +138,13 @@ void CbmL1PFFitter::Fit(vector<CbmStsTrack>& Tracks, vector<int>& pidHypo)
static int nHits = CbmL1::Instance()->fpAlgo->GetParameters()->GetNstationsActive();
int iVec = 0, i = 0;
int nTracks_SIMD = fvec::size();
- L1TrackPar T; // fitting parametr coresponding to current track
- L1Fit fit;
- fit.SetParticleMass(0.000511f); // muon
+ L1TrackParFit fit;
+ fit.SetParticleMass(CbmL1::Instance()->fpAlgo->GetDefaultParticleMass());
- CbmStsTrack* t[fvec::size()] {nullptr};
+ L1TrackPar& T = fit.fTr; // fitting parametr coresponding to current track
+
+ CbmStsTrack* tr[fvec::size()] {nullptr};
int ista;
const L1Station* sta = CbmL1::Instance()->fpAlgo->GetParameters()->GetStations().begin();
@@ -151,8 +152,10 @@ void CbmL1PFFitter::Fit(vector<CbmStsTrack>& Tracks, vector<int>& pidHypo)
fvec* x = new fvec[nHits];
fvec* u = new fvec[nHits];
fvec* v = new fvec[nHits];
+ fvec* t = new fvec[nHits];
fvec* du2 = new fvec[nHits];
fvec* dv2 = new fvec[nHits];
+ fvec* dt2 = new fvec[nHits];
fvec* y = new fvec[nHits];
fvec* z = new fvec[nHits];
fvec* w = new fvec[nHits];
@@ -178,28 +181,28 @@ void CbmL1PFFitter::Fit(vector<CbmStsTrack>& Tracks, vector<int>& pidHypo)
if (N_vTracks - itrack < static_cast<unsigned short>(fvec::size())) nTracks_SIMD = N_vTracks - itrack;
for (i = 0; i < nTracks_SIMD; i++) {
- t[i] = &Tracks[itrack + i]; // current track
- T.x[i] = t[i]->GetParamFirst()->GetX();
- T.y[i] = t[i]->GetParamFirst()->GetY();
- T.tx[i] = t[i]->GetParamFirst()->GetTx();
- T.ty[i] = t[i]->GetParamFirst()->GetTy();
- T.qp[i] = t[i]->GetParamFirst()->GetQp();
- T.z[i] = t[i]->GetParamFirst()->GetZ();
- T.C00[i] = t[i]->GetParamFirst()->GetCovariance(0, 0);
- T.C10[i] = t[i]->GetParamFirst()->GetCovariance(1, 0);
- T.C11[i] = t[i]->GetParamFirst()->GetCovariance(1, 1);
- T.C20[i] = t[i]->GetParamFirst()->GetCovariance(2, 0);
- T.C21[i] = t[i]->GetParamFirst()->GetCovariance(2, 1);
- T.C22[i] = t[i]->GetParamFirst()->GetCovariance(2, 2);
- T.C30[i] = t[i]->GetParamFirst()->GetCovariance(3, 0);
- T.C31[i] = t[i]->GetParamFirst()->GetCovariance(3, 1);
- T.C32[i] = t[i]->GetParamFirst()->GetCovariance(3, 2);
- T.C33[i] = t[i]->GetParamFirst()->GetCovariance(3, 3);
- T.C40[i] = t[i]->GetParamFirst()->GetCovariance(4, 0);
- T.C41[i] = t[i]->GetParamFirst()->GetCovariance(4, 1);
- T.C42[i] = t[i]->GetParamFirst()->GetCovariance(4, 1);
- T.C43[i] = t[i]->GetParamFirst()->GetCovariance(4, 3);
- T.C44[i] = t[i]->GetParamFirst()->GetCovariance(4, 4);
+ tr[i] = &Tracks[itrack + i]; // current track
+ T.x[i] = tr[i]->GetParamFirst()->GetX();
+ T.y[i] = tr[i]->GetParamFirst()->GetY();
+ T.tx[i] = tr[i]->GetParamFirst()->GetTx();
+ T.ty[i] = tr[i]->GetParamFirst()->GetTy();
+ T.qp[i] = tr[i]->GetParamFirst()->GetQp();
+ T.z[i] = tr[i]->GetParamFirst()->GetZ();
+ T.C00[i] = tr[i]->GetParamFirst()->GetCovariance(0, 0);
+ T.C10[i] = tr[i]->GetParamFirst()->GetCovariance(1, 0);
+ T.C11[i] = tr[i]->GetParamFirst()->GetCovariance(1, 1);
+ T.C20[i] = tr[i]->GetParamFirst()->GetCovariance(2, 0);
+ T.C21[i] = tr[i]->GetParamFirst()->GetCovariance(2, 1);
+ T.C22[i] = tr[i]->GetParamFirst()->GetCovariance(2, 2);
+ T.C30[i] = tr[i]->GetParamFirst()->GetCovariance(3, 0);
+ T.C31[i] = tr[i]->GetParamFirst()->GetCovariance(3, 1);
+ T.C32[i] = tr[i]->GetParamFirst()->GetCovariance(3, 2);
+ T.C33[i] = tr[i]->GetParamFirst()->GetCovariance(3, 3);
+ T.C40[i] = tr[i]->GetParamFirst()->GetCovariance(4, 0);
+ T.C41[i] = tr[i]->GetParamFirst()->GetCovariance(4, 1);
+ T.C42[i] = tr[i]->GetParamFirst()->GetCovariance(4, 1);
+ T.C43[i] = tr[i]->GetParamFirst()->GetCovariance(4, 3);
+ T.C44[i] = tr[i]->GetParamFirst()->GetCovariance(4, 4);
int pid = pidHypo[itrack + i];
if (pid == -1) pid = 211;
@@ -215,35 +218,40 @@ void CbmL1PFFitter::Fit(vector<CbmStsTrack>& Tracks, vector<int>& pidHypo)
}
for (iVec = 0; iVec < nTracks_SIMD; iVec++) {
- int nHitsTrackMvd = t[iVec]->GetNofMvdHits();
- int nHitsTrackSts = t[iVec]->GetNofStsHits();
+ int nHitsTrackMvd = tr[iVec]->GetNofMvdHits();
+ int nHitsTrackSts = tr[iVec]->GetNofStsHits();
int nHitsTrack = nHitsTrackMvd + nHitsTrackSts;
for (i = 0; i < nHitsTrack; i++) {
- float posx = 0.f, posy = 0.f, posz = 0.f;
+ float posx = 0.f, posy = 0.f, posz = 0.f, time = 0.;
if (i < nHitsTrackMvd) {
if (!listMvdHits) continue;
- int hitIndex = t[iVec]->GetMvdHitIndex(i);
+ int hitIndex = tr[iVec]->GetMvdHitIndex(i);
CbmMvdHit* hit = L1_DYNAMIC_CAST<CbmMvdHit*>(listMvdHits->At(hitIndex));
posx = hit->GetX();
posy = hit->GetY();
posz = hit->GetZ();
+ time = hit->GetTime();
+
// ista = hit->GetStationNr();
ista =
CbmL1::Instance()->fpAlgo->GetParameters()->GetStationIndexActive(hit->GetStationNr(), L1DetectorID::kMvd);
if (ista == -1) continue;
du2[ista][iVec] = hit->GetDx() * hit->GetDx();
dv2[ista][iVec] = hit->GetDy() * hit->GetDy();
+ dt2[ista][iVec] = hit->GetTimeError() * hit->GetTimeError();
}
else {
if (!listStsHits) continue;
- int hitIndex = t[iVec]->GetHitIndex(i - nHitsTrackMvd);
+ int hitIndex = tr[iVec]->GetHitIndex(i - nHitsTrackMvd);
CbmStsHit* hit = L1_DYNAMIC_CAST<CbmStsHit*>(listStsHits->At(hitIndex));
posx = hit->GetX();
posy = hit->GetY();
posz = hit->GetZ();
+ time = hit->GetTime();
+
// ista = CbmStsSetup::Instance()->GetStationNumber(hit->GetAddress())
// + NMvdStations; //hit->GetStationNr() - 1 + NMvdStations;
ista = CbmL1::Instance()->fpAlgo->GetParameters()->GetStationIndexActive(
@@ -251,6 +259,7 @@ void CbmL1PFFitter::Fit(vector<CbmStsTrack>& Tracks, vector<int>& pidHypo)
if (ista == -1) continue;
du2[ista][iVec] = hit->GetDu() * hit->GetDu();
dv2[ista][iVec] = hit->GetDv() * hit->GetDv();
+ dt2[ista][iVec] = hit->GetTimeError() * hit->GetTimeError();
}
w[ista][iVec] = 1.f;
@@ -259,6 +268,7 @@ void CbmL1PFFitter::Fit(vector<CbmStsTrack>& Tracks, vector<int>& pidHypo)
u[ista][iVec] = posx * sta[ista].frontInfo.cos_phi[0] + posy * sta[ista].frontInfo.sin_phi[0];
v[ista][iVec] = posx * sta[ista].backInfo.cos_phi[0] + posy * sta[ista].backInfo.sin_phi[0];
z[ista][iVec] = posz;
+ t[ista][iVec] = time;
sta[ista].fieldSlice.GetFieldValue(x[ista], y[ista], fB_temp);
fB[ista].x[iVec] = fB_temp.x[iVec];
@@ -291,10 +301,10 @@ void CbmL1PFFitter::Fit(vector<CbmStsTrack>& Tracks, vector<int>& pidHypo)
// fit forward
i = 0;
- FilterFirst(T, x_first, y_first, fstC00, fstC10, fstC11);
+ FilterFirst(fit, x_first, y_first, fstC00, fstC10, fstC11);
+ fit.fQp0 = fit.fTr.qp;
- fvec qp0 = T.qp;
- z1 = z[i];
+ z1 = z[i];
sta[i].fieldSlice.GetFieldValue(T.x, T.y, b1);
b1.Combine(fB[i], w[i]);
z2 = z[i + 2];
@@ -313,16 +323,13 @@ void CbmL1PFFitter::Fit(vector<CbmStsTrack>& Tracks, vector<int>& pidHypo)
fvec w1 = iif(initialised, w[i], fvec::Zero());
fvec wIn = iif(initialised, fvec::One(), fvec::Zero());
- L1Extrapolate(T, z[i], qp0, fld, &w1);
- if (i == NMvdStations) { // TODO: How a hit can be also a station? (S.Zharko)
- fit.L1AddPipeMaterial(T, qp0, wIn);
- fit.EnergyLossCorrection(T, fit.PipeRadThick, qp0, fvec(-1.f), wIn);
- }
- fit.L1AddMaterial(T, CbmL1::Instance()->fpAlgo->GetParameters()->GetMaterialThickness(i, T.x, T.y), qp0, wIn);
- fit.EnergyLossCorrection(T, CbmL1::Instance()->fpAlgo->GetParameters()->GetMaterialThickness(i, T.x, T.y), qp0,
- -1, wIn);
- L1Filter(T, sta[i].frontInfo, u[i], du2[i], w1);
- L1Filter(T, sta[i].backInfo, v[i], dv2[i], w1);
+ fit.Extrapolate(z[i], fit.fQp0, fld, w1);
+ auto radThick = CbmL1::Instance()->fpAlgo->GetParameters()->GetMaterialThickness(i, fit.fTr.x, fit.fTr.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, fit.fQp0, wIn);
+ fit.EnergyLossCorrection(radThick, fit.fQp0, -fvec::One(), wIn);
b2 = b1;
z2 = z1;
@@ -355,15 +362,15 @@ void CbmL1PFFitter::Fit(vector<CbmStsTrack>& Tracks, vector<int>& pidHypo)
par.SetCovariance(4, 2, Tout.C42[iVec]);
par.SetCovariance(4, 3, Tout.C43[iVec]);
par.SetCovariance(4, 4, Tout.C44[iVec]);
- t[iVec]->SetParamLast(&par);
+ tr[iVec]->SetParamLast(&par);
}
//fit backward
- qp0 = T.qp;
+ fit.fQp0 = T.qp;
i = nHits - 1;
- FilterFirst(T, x_last, y_last, lstC00, lstC10, lstC11);
+ FilterFirst(fit, x_last, y_last, lstC00, lstC10, lstC11);
z1 = z[i];
sta[i].fieldSlice.GetFieldValue(T.x, T.y, b1);
@@ -385,16 +392,13 @@ void CbmL1PFFitter::Fit(vector<CbmStsTrack>& Tracks, vector<int>& pidHypo)
fvec w1 = iif(initialised, w[i], fvec::Zero());
fvec wIn = iif(initialised, fvec::One(), fvec::Zero());
- L1Extrapolate(T, z[i], qp0, fld, &w1);
- if (i == NMvdStations - 1) {
- fit.L1AddPipeMaterial(T, qp0, wIn);
- fit.EnergyLossCorrection(T, fit.PipeRadThick, qp0, fvec(1.f), wIn);
- }
- fit.L1AddMaterial(T, CbmL1::Instance()->fpAlgo->GetParameters()->GetMaterialThickness(i, T.x, T.y), qp0, wIn);
- fit.EnergyLossCorrection(T, CbmL1::Instance()->fpAlgo->GetParameters()->GetMaterialThickness(i, T.x, T.y), qp0, 1,
- wIn);
- L1Filter(T, sta[i].frontInfo, u[i], du2[i], w1);
- L1Filter(T, sta[i].backInfo, v[i], dv2[i], w1);
+ fit.Extrapolate(z[i], fit.fQp0, fld, w1);
+ auto radThick = CbmL1::Instance()->fpAlgo->GetParameters()->GetMaterialThickness(i, fit.fTr.x, fit.fTr.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, fit.fQp0, wIn);
+ fit.EnergyLossCorrection(radThick, fit.fQp0, fvec::One(), wIn);
b2 = b1;
z2 = z1;
@@ -426,10 +430,10 @@ void CbmL1PFFitter::Fit(vector<CbmStsTrack>& Tracks, vector<int>& pidHypo)
par.SetCovariance(4, 2, T.C42[iVec]);
par.SetCovariance(4, 3, T.C43[iVec]);
par.SetCovariance(4, 4, T.C44[iVec]);
- t[iVec]->SetParamFirst(&par);
+ tr[iVec]->SetParamFirst(&par);
- t[iVec]->SetChiSq(T.chi2[iVec]);
- t[iVec]->SetNDF(static_cast<int>(T.NDF[iVec]));
+ tr[iVec]->SetChiSq(T.chi2[iVec]);
+ tr[iVec]->SetNDF(static_cast<int>(T.NDF[iVec]));
}
}
@@ -448,10 +452,11 @@ void CbmL1PFFitter::GetChiToVertex(vector<CbmStsTrack>& Tracks, vector<PFFieldRe
chiToVtx.reserve(Tracks.size());
int nTracks_SIMD = fvec::size();
- L1TrackPar T; // fitting parametr coresponding to current track
- L1Fit fit;
- CbmStsTrack* t[fvec::size()] {nullptr};
+ L1TrackParFit fit;
+ L1TrackPar& T = fit.fTr; // fitting parametr coresponding to current track
+
+ CbmStsTrack* tr[fvec::size()] {nullptr};
int nStations = CbmL1::Instance()->fpAlgo->GetParameters()->GetNstationsActive();
int NMvdStations = CbmL1::Instance()->fpAlgo->GetNstationsBeforePipe();
@@ -478,39 +483,39 @@ void CbmL1PFFitter::GetChiToVertex(vector<CbmStsTrack>& Tracks, vector<PFFieldRe
fvec mass2;
for (int iVec = 0; iVec < nTracks_SIMD; iVec++) {
- t[iVec] = &Tracks[itrack + iVec]; // current track
- T.x[iVec] = t[iVec]->GetParamFirst()->GetX();
- T.y[iVec] = t[iVec]->GetParamFirst()->GetY();
- T.tx[iVec] = t[iVec]->GetParamFirst()->GetTx();
- T.ty[iVec] = t[iVec]->GetParamFirst()->GetTy();
- T.qp[iVec] = t[iVec]->GetParamFirst()->GetQp();
- T.z[iVec] = t[iVec]->GetParamFirst()->GetZ();
- T.C00[iVec] = t[iVec]->GetParamFirst()->GetCovariance(0, 0);
- T.C10[iVec] = t[iVec]->GetParamFirst()->GetCovariance(1, 0);
- T.C11[iVec] = t[iVec]->GetParamFirst()->GetCovariance(1, 1);
- T.C20[iVec] = t[iVec]->GetParamFirst()->GetCovariance(2, 0);
- T.C21[iVec] = t[iVec]->GetParamFirst()->GetCovariance(2, 1);
- T.C22[iVec] = t[iVec]->GetParamFirst()->GetCovariance(2, 2);
- T.C30[iVec] = t[iVec]->GetParamFirst()->GetCovariance(3, 0);
- T.C31[iVec] = t[iVec]->GetParamFirst()->GetCovariance(3, 1);
- T.C32[iVec] = t[iVec]->GetParamFirst()->GetCovariance(3, 2);
- T.C33[iVec] = t[iVec]->GetParamFirst()->GetCovariance(3, 3);
- T.C40[iVec] = t[iVec]->GetParamFirst()->GetCovariance(4, 0);
- T.C41[iVec] = t[iVec]->GetParamFirst()->GetCovariance(4, 1);
- T.C42[iVec] = t[iVec]->GetParamFirst()->GetCovariance(4, 1);
- T.C43[iVec] = t[iVec]->GetParamFirst()->GetCovariance(4, 3);
- T.C44[iVec] = t[iVec]->GetParamFirst()->GetCovariance(4, 4);
- // float mass = TDatabasePDG::Instance()->GetParticle(t[iVec]->GetPidHypo())->Mass();
- const float mass = KFParticleDatabase::Instance()->GetMass(t[iVec]->GetPidHypo());
+ tr[iVec] = &Tracks[itrack + iVec]; // current track
+ T.x[iVec] = tr[iVec]->GetParamFirst()->GetX();
+ T.y[iVec] = tr[iVec]->GetParamFirst()->GetY();
+ T.tx[iVec] = tr[iVec]->GetParamFirst()->GetTx();
+ T.ty[iVec] = tr[iVec]->GetParamFirst()->GetTy();
+ T.qp[iVec] = tr[iVec]->GetParamFirst()->GetQp();
+ T.z[iVec] = tr[iVec]->GetParamFirst()->GetZ();
+ T.C00[iVec] = tr[iVec]->GetParamFirst()->GetCovariance(0, 0);
+ T.C10[iVec] = tr[iVec]->GetParamFirst()->GetCovariance(1, 0);
+ T.C11[iVec] = tr[iVec]->GetParamFirst()->GetCovariance(1, 1);
+ T.C20[iVec] = tr[iVec]->GetParamFirst()->GetCovariance(2, 0);
+ T.C21[iVec] = tr[iVec]->GetParamFirst()->GetCovariance(2, 1);
+ T.C22[iVec] = tr[iVec]->GetParamFirst()->GetCovariance(2, 2);
+ T.C30[iVec] = tr[iVec]->GetParamFirst()->GetCovariance(3, 0);
+ T.C31[iVec] = tr[iVec]->GetParamFirst()->GetCovariance(3, 1);
+ T.C32[iVec] = tr[iVec]->GetParamFirst()->GetCovariance(3, 2);
+ T.C33[iVec] = tr[iVec]->GetParamFirst()->GetCovariance(3, 3);
+ T.C40[iVec] = tr[iVec]->GetParamFirst()->GetCovariance(4, 0);
+ T.C41[iVec] = tr[iVec]->GetParamFirst()->GetCovariance(4, 1);
+ T.C42[iVec] = tr[iVec]->GetParamFirst()->GetCovariance(4, 1);
+ T.C43[iVec] = tr[iVec]->GetParamFirst()->GetCovariance(4, 3);
+ T.C44[iVec] = tr[iVec]->GetParamFirst()->GetCovariance(4, 4);
+ // float mass = TDatabasePDG::Instance()->GetParticle(tr[iVec]->GetPidHypo())->Mass();
+ const float mass = KFParticleDatabase::Instance()->GetMass(tr[iVec]->GetPidHypo());
mass2[iVec] = mass * mass;
- int nHitsTrackMvd = t[iVec]->GetNofMvdHits();
+ int nHitsTrackMvd = tr[iVec]->GetNofMvdHits();
for (int iH = 0; iH < 2; iH++) {
float posx = 0.f, posy = 0.f, posz = 0.f;
if (iH < nHitsTrackMvd) {
if (!listMvdHits) continue;
- int hitIndex = t[iVec]->GetMvdHitIndex(iH);
+ int hitIndex = tr[iVec]->GetMvdHitIndex(iH);
CbmMvdHit* hit = L1_DYNAMIC_CAST<CbmMvdHit*>(listMvdHits->At(hitIndex));
posx = hit->GetX();
@@ -520,7 +525,7 @@ void CbmL1PFFitter::GetChiToVertex(vector<CbmStsTrack>& Tracks, vector<PFFieldRe
}
else {
if (!listStsHits) continue;
- int hitIndex = t[iVec]->GetHitIndex(iH - nHitsTrackMvd);
+ int hitIndex = tr[iVec]->GetHitIndex(iH - nHitsTrackMvd);
CbmStsHit* hit = L1_DYNAMIC_CAST<CbmStsHit*>(listStsHits->At(hitIndex));
posx = hit->GetX();
@@ -545,25 +550,20 @@ void CbmL1PFFitter::GetChiToVertex(vector<CbmStsTrack>& Tracks, vector<PFFieldRe
field.emplace_back(fld, i);
}
+ fit.fQp0 = fit.fTr.qp;
+
for (int iSt = nStations - 4; iSt >= 0; iSt--) {
fvec w = iif(T.z > (zSta[iSt] + fvec(2.5)), fvec::One(), fvec::Zero());
- L1Extrapolate(T, zSta[iSt], T.qp, fld, &w);
- if (iSt == NMvdStations - 1) {
- fit.L1AddPipeMaterial(T, T.qp, w);
- fit.EnergyLossCorrection(T, fit.PipeRadThick, T.qp, fvec(1.f), w);
- }
- fit.L1AddMaterial(T, CbmL1::Instance()->fpAlgo->GetParameters()->GetMaterialThickness(iSt, T.x, T.y), T.qp, w);
- fit.EnergyLossCorrection(T, CbmL1::Instance()->fpAlgo->GetParameters()->GetMaterialThickness(iSt, T.x, T.y), T.qp,
- fvec(1.f), w);
- }
- if (NMvdStations <= 0) {
- fit.L1AddPipeMaterial(T, T.qp, ONE);
- fit.EnergyLossCorrection(T, fit.PipeRadThick, T.qp, fvec(1.f), ONE);
+ fit.Extrapolate(zSta[iSt], fit.fQp0, fld, w);
+ auto radThick = CbmL1::Instance()->fpAlgo->GetParameters()->GetMaterialThickness(iSt, fit.fTr.x, fit.fTr.y);
+ fit.AddMsInMaterial(radThick, fit.fQp0, w);
+ fit.EnergyLossCorrection(radThick, fit.fQp0, fvec::One(), w);
}
- L1Extrapolate(T, primVtx.GetRefZ(), T.qp, fld);
- fit.L1AddTargetMaterial(T, T.qp, 1);
+ fit.Extrapolate(primVtx.GetRefZ(), fit.fQp0, fld, fvec::One());
+ fit.AddMsInMaterial(fit.GetTargetRadThick(), fit.fQp0, fvec::One());
+ fit.EnergyLossCorrection(fit.GetTargetRadThick(), fit.fQp0, fvec::One(), fvec::One());
Double_t Cv[3] = {primVtx.GetCovMatrix()[0], primVtx.GetCovMatrix()[1], primVtx.GetCovMatrix()[2]};
@@ -607,7 +607,7 @@ void CbmL1PFFitter::GetChiToVertex(vector<CbmStsTrack>& Tracks, vector<PFFieldRe
par.SetCovariance(4, 2, T.C42[iVec]);
par.SetCovariance(4, 3, T.C43[iVec]);
par.SetCovariance(4, 4, T.C44[iVec]);
- t[iVec]->SetParamFirst(&par);
+ tr[iVec]->SetParamFirst(&par);
}
}
}
@@ -630,7 +630,7 @@ void CbmL1PFFitter::CalculateFieldRegion(vector<CbmStsTrack>& Tracks, vector<PFF
int nTracks_SIMD = fvec::size();
L1TrackPar T; // fitting parametr coresponding to current track
- CbmStsTrack* t[fvec::size()];
+ CbmStsTrack* tr[fvec::size()];
int ista;
const L1Station* sta = CbmL1::Instance()->fpAlgo->GetParameters()->GetStations().begin();
@@ -643,16 +643,16 @@ void CbmL1PFFitter::CalculateFieldRegion(vector<CbmStsTrack>& Tracks, vector<PFF
if (N_vTracks - itrack < static_cast<unsigned short>(fvec::size())) nTracks_SIMD = N_vTracks - itrack;
for (int i = 0; i < nTracks_SIMD; i++)
- t[i] = &Tracks[itrack + i]; // current track
+ tr[i] = &Tracks[itrack + i]; // current track
for (int iVec = 0; iVec < nTracks_SIMD; iVec++) {
- int nHitsTrackMvd = t[iVec]->GetNofMvdHits();
+ int nHitsTrackMvd = tr[iVec]->GetNofMvdHits();
for (int iH = 0; iH < 2; iH++) {
float posx = 0.f, posy = 0.f, posz = 0.f;
if (iH < nHitsTrackMvd) {
if (!listMvdHits) continue;
- int hitIndex = t[iVec]->GetMvdHitIndex(iH);
+ int hitIndex = tr[iVec]->GetMvdHitIndex(iH);
CbmMvdHit* hit = L1_DYNAMIC_CAST<CbmMvdHit*>(listMvdHits->At(hitIndex));
posx = hit->GetX();
@@ -662,7 +662,7 @@ void CbmL1PFFitter::CalculateFieldRegion(vector<CbmStsTrack>& Tracks, vector<PFF
}
else {
if (!listStsHits) continue;
- int hitIndex = t[iVec]->GetHitIndex(iH - nHitsTrackMvd);
+ int hitIndex = tr[iVec]->GetHitIndex(iH - nHitsTrackMvd);
CbmStsHit* hit = L1_DYNAMIC_CAST<CbmStsHit*>(listStsHits->At(hitIndex));
posx = hit->GetX();
@@ -704,7 +704,7 @@ void CbmL1PFFitter::CalculateFieldRegionAtLastPoint(vector<CbmStsTrack>& Tracks,
int nTracks_SIMD = fvec::size();
L1TrackPar T; // fitting parametr coresponding to current track
- CbmStsTrack* t[fvec::size()];
+ CbmStsTrack* tr[fvec::size()];
int ista;
const L1Station* sta = CbmL1::Instance()->fpAlgo->GetParameters()->GetStations().begin();
@@ -717,19 +717,19 @@ void CbmL1PFFitter::CalculateFieldRegionAtLastPoint(vector<CbmStsTrack>& Tracks,
if (N_vTracks - itrack < static_cast<unsigned short>(fvec::size())) nTracks_SIMD = N_vTracks - itrack;
for (int i = 0; i < nTracks_SIMD; i++) {
- t[i] = &Tracks[itrack + i]; // current track
+ tr[i] = &Tracks[itrack + i]; // current track
}
for (int iVec = 0; iVec < nTracks_SIMD; iVec++) {
- int nHitsTrackMvd = t[iVec]->GetNofMvdHits();
- int nHits = t[iVec]->GetNofHits();
+ int nHitsTrackMvd = tr[iVec]->GetNofMvdHits();
+ int nHits = tr[iVec]->GetNofHits();
for (int iH = 0; iH < 3; iH++) {
float posx = 0.f, posy = 0.f, posz = 0.f;
int hitNumber = nHits - iH - 1;
if (hitNumber < nHitsTrackMvd) {
if (!listMvdHits) continue;
- int hitIndex = t[iVec]->GetMvdHitIndex(hitNumber);
+ int hitIndex = tr[iVec]->GetMvdHitIndex(hitNumber);
CbmMvdHit* hit = L1_DYNAMIC_CAST<CbmMvdHit*>(listMvdHits->At(hitIndex));
posx = hit->GetX();
@@ -739,7 +739,7 @@ void CbmL1PFFitter::CalculateFieldRegionAtLastPoint(vector<CbmStsTrack>& Tracks,
}
else {
if (!listStsHits) continue;
- int hitIndex = t[iVec]->GetHitIndex(hitNumber - nHitsTrackMvd);
+ int hitIndex = tr[iVec]->GetHitIndex(hitNumber - nHitsTrackMvd);
CbmStsHit* hit = L1_DYNAMIC_CAST<CbmStsHit*>(listStsHits->At(hitIndex));
posx = hit->GetX();
--
GitLab