From 5f90ed3589c2235753375ed1425619c67d5979fc Mon Sep 17 00:00:00 2001
From: "se.gorbunov" <se.gorbunov@gsi.de>
Date: Fri, 23 Aug 2024 20:34:59 +0200
Subject: [PATCH] Ca: fit: use enum for upstream/downstream direction flag
---
algo/ca/core/tracking/CaTrackExtender.cxx | 42 +++++++++----------
algo/ca/core/tracking/CaTrackExtender.h | 8 ++--
algo/ca/core/tracking/CaTrackFit.cxx | 20 +++++----
algo/ca/core/tracking/CaTrackFit.h | 20 ++++++---
algo/ca/core/tracking/CaTrackFitter.cxx | 4 +-
.../ca/core/tracking/CaTripletConstructor.cxx | 12 +++---
reco/KF/CbmKfTrackFitter.cxx | 2 +-
reco/KF/CbmKfTrackFitter.h | 7 +---
reco/KF/ParticleFitter/CbmL1PFFitter.cxx | 12 +++---
reco/L1/CbmL1Performance.cxx | 4 +-
reco/L1/qa/CbmCaTrackTypeQa.cxx | 22 +++++-----
11 files changed, 80 insertions(+), 73 deletions(-)
diff --git a/algo/ca/core/tracking/CaTrackExtender.cxx b/algo/ca/core/tracking/CaTrackExtender.cxx
index 2b7e606765..40cb9e3ed0 100644
--- a/algo/ca/core/tracking/CaTrackExtender.cxx
+++ b/algo/ca/core/tracking/CaTrackExtender.cxx
@@ -37,7 +37,7 @@ TrackExtender::~TrackExtender() {}
// ---------------------------------------------------------------------------------------------------------------------
//
-void TrackExtender::FitBranchFast(const ca::Branch& t, TrackParamV& Tout, const bool upstream, const fvec qp0,
+void TrackExtender::FitBranchFast(const ca::Branch& t, TrackParamV& Tout, const FitDirection direction, const fvec qp0,
const bool initParams)
{
CBMCA_DEBUG_ASSERT(t.NHits >= 3);
@@ -52,9 +52,9 @@ void TrackExtender::FitBranchFast(const ca::Branch& t, TrackParamV& Tout, const
const Vector<ca::HitIndex_t>& hits = t.Hits(); // array of indeses of hits of current track
const int nHits = t.NofHits();
- const signed short int step = -2 * static_cast<int>(upstream) + 1; // increment for station index
- const int iFirstHit = (upstream) ? nHits - 1 : 0;
- const int iLastHit = (upstream) ? 0 : nHits - 1;
+ const signed short int step = (direction == kUpstream ? -1 : 1); // increment for station index
+ const int iFirstHit = (direction == kUpstream) ? nHits - 1 : 0;
+ const int iLastHit = (direction == kUpstream) ? 0 : nHits - 1;
const ca::Hit& hit0 = fInputData->GetHit(hits[iFirstHit]);
const ca::Hit& hit1 = fInputData->GetHit(hits[iFirstHit + step]);
@@ -122,7 +122,7 @@ void TrackExtender::FitBranchFast(const ca::Branch& t, TrackParamV& Tout, const
fit.FilterHit(hit, fmask(sta.timeInfo));
auto radThick = fParameters.GetMaterialThickness(ista, T.X(), T.Y());
fit.MultipleScattering(radThick);
- fit.EnergyLossCorrection(radThick, fvec(upstream ? 1. : -1.));
+ fit.EnergyLossCorrection(radThick, direction);
fldB0 = fldB1;
fldB1 = fldB2;
@@ -137,18 +137,18 @@ void TrackExtender::FitBranchFast(const ca::Branch& t, TrackParamV& Tout, const
} // void ca::Framework::BranchFitterFast
/// like BranchFitterFast but more precise
-void TrackExtender::FitBranch(const ca::Branch& t, TrackParamV& T, const bool upstream, const fvec qp0,
+void TrackExtender::FitBranch(const ca::Branch& t, TrackParamV& T, const FitDirection direction, const fvec qp0,
const bool initParams)
{
- FitBranchFast(t, T, upstream, qp0, initParams);
+ FitBranchFast(t, T, direction, qp0, initParams);
for (int i = 0; i < 1; i++) {
- FitBranchFast(t, T, !upstream, T.Qp(), false);
- FitBranchFast(t, T, upstream, T.Qp(), false);
+ FitBranchFast(t, T, !direction, T.Qp(), false);
+ FitBranchFast(t, T, direction, T.Qp(), false);
}
} // void ca::Framework::BranchFitter
-void TrackExtender::FindMoreHits(ca::Branch& t, TrackParamV& Tout, const bool upstream, const fvec qp0)
+void TrackExtender::FindMoreHits(ca::Branch& t, TrackParamV& Tout, const FitDirection direction, const fvec qp0)
{
Vector<ca::HitIndex_t> newHits{"ca::TrackExtender::newHits"};
newHits.reserve(fParameters.GetNstationsActive());
@@ -159,8 +159,8 @@ void TrackExtender::FindMoreHits(ca::Branch& t, TrackParamV& Tout, const bool up
fit.SetTrack(Tout);
fit.SetQp0(qp0);
- const signed short int step = -2 * static_cast<int>(upstream) + 1; // increment for station index
- const int iFirstHit = (upstream) ? 2 : t.NofHits() - 3;
+ const signed short int step = (direction == kUpstream) ? -1 : 1; // increment for station index
+ const int iFirstHit = (direction == kUpstream) ? 2 : t.NofHits() - 3;
// int ista = fInputData->GetHit(t.Hits[iFirstHit]).iSt + 2 * step; // current station. set to the end of track
const ca::Hit& hit0 = fInputData->GetHit(t.Hits()[iFirstHit]); // optimize
@@ -274,7 +274,7 @@ void TrackExtender::FindMoreHits(ca::Branch& t, TrackParamV& Tout, const bool up
fit.FilterHit(hit, fmask(sta.timeInfo));
auto radThick = fParameters.GetMaterialThickness(ista, tr.X(), tr.Y());
fit.MultipleScattering(radThick);
- fit.EnergyLossCorrection(radThick, upstream ? fvec(1.f) : fvec(-1.f));
+ fit.EnergyLossCorrection(radThick, direction);
fldB0 = fldB1;
fldB1 = fldB2;
@@ -290,7 +290,7 @@ void TrackExtender::FindMoreHits(ca::Branch& t, TrackParamV& Tout, const bool up
const unsigned int NNewHits = newHits.size();
t.RefHits().enlarge(NNewHits + NOldHits);
- if (upstream) { // backward
+ if (direction == kUpstream) {
for (int i = NOldHits - 1; i >= 0; i--) {
t.RefHits()[NNewHits + i] = t.RefHits()[i];
}
@@ -298,7 +298,7 @@ void TrackExtender::FindMoreHits(ca::Branch& t, TrackParamV& Tout, const bool up
t.RefHits()[i] = newHits[ii];
}
}
- else { // forward
+ else { // downstream
for (unsigned int i = 0; i < newHits.size(); i++) {
t.RefHits()[NOldHits + i] = newHits[i];
}
@@ -318,18 +318,14 @@ fscal TrackExtender::ExtendBranch(ca::Branch& t, const InputData& input, WindowD
TrackParamV T;
// downstream
- bool upstream = 0;
- FitBranch(t, T, upstream, 0.0);
-
- FindMoreHits(t, T, upstream, T.Qp());
+ FitBranch(t, T, kDownstream, 0.0);
+ FindMoreHits(t, T, kDownstream, T.Qp());
// upstream
- upstream = 1;
- FitBranchFast(t, T, upstream, T.Qp(), false);
-
- FindMoreHits(t, T, upstream, T.Qp());
+ FitBranchFast(t, T, kUpstream, T.Qp(), false);
+ FindMoreHits(t, T, kUpstream, T.Qp());
return T.GetChiSq()[0];
}
diff --git a/algo/ca/core/tracking/CaTrackExtender.h b/algo/ca/core/tracking/CaTrackExtender.h
index 5368e84439..2c9ae46e99 100644
--- a/algo/ca/core/tracking/CaTrackExtender.h
+++ b/algo/ca/core/tracking/CaTrackExtender.h
@@ -13,6 +13,7 @@
#include "CaHit.h"
#include "CaParameters.h"
#include "CaSimd.h"
+#include "CaTrackFit.h"
#include "CaVector.h"
#include "CaWindowData.h"
#include "KfTrackParam.h"
@@ -67,7 +68,7 @@ namespace cbm::algo::ca
/// \param T - track params
/// \param dir - 0 - forward, 1 - backward
/// \param qp0 - momentum value to linearize the extrapolation
- void FindMoreHits(ca::Branch& t, TrackParamV& T, const bool dir, const fvec qp0);
+ void FindMoreHits(ca::Branch& t, TrackParamV& T, const FitDirection direction, const fvec qp0);
/// Fits the branch. Does few passes over the hits.
/// \param t - track branch with hits
@@ -75,7 +76,8 @@ namespace cbm::algo::ca
/// \param dir - false - forward, true - backward
/// \param qp0 - momentum value to linearize the extrapolation
/// \param initParams - should params be ititialized. 1 - yes.
- void FitBranch(const ca::Branch& t, TrackParamV& T, const bool dir, const fvec qp0, const bool initParams = true);
+ void FitBranch(const ca::Branch& t, TrackParamV& T, const FitDirection direction, const fvec qp0,
+ const bool initParams = true);
/// Fits the branch. Does only one pass over the hits.
@@ -84,7 +86,7 @@ namespace cbm::algo::ca
/// \param dir - false - forward, true - backward
/// \param qp0 - momentum value to linearize the extrapolation
/// \param initParams - should params be ititialized. 1 - yes.
- void FitBranchFast(const ca::Branch& t, TrackParamV& T, const bool dir, const fvec qp0,
+ void FitBranchFast(const ca::Branch& t, TrackParamV& T, const FitDirection direction, const fvec qp0,
const bool initParams = true);
private:
diff --git a/algo/ca/core/tracking/CaTrackFit.cxx b/algo/ca/core/tracking/CaTrackFit.cxx
index 3666eef749..62a03a96b9 100644
--- a/algo/ca/core/tracking/CaTrackFit.cxx
+++ b/algo/ca/core/tracking/CaTrackFit.cxx
@@ -1054,7 +1054,7 @@ namespace cbm::algo::ca
template<typename DataT>
- void TrackFit<DataT>::EnergyLossCorrection(DataT radThick, DataT upstreamDirection)
+ void TrackFit<DataT>::EnergyLossCorrection(DataT radThick, FitDirection direction)
{
cnst qp2cut(1. / (10. * 10.)); // 10 GeV cut
cnst qp02 = utils::max(fQp0 * fQp0, qp2cut);
@@ -1065,9 +1065,11 @@ namespace cbm::algo::ca
DataT tr = sqrt(DataT(1.f) + fTr.Tx() * fTr.Tx() + fTr.Ty() * fTr.Ty());
- cnst dE = bethe * radThick * tr * 2.33f * 9.34961f;
+ DataT dE = bethe * radThick * tr * 2.33f * 9.34961f;
- cnst ECorrected = sqrt(E2) + upstreamDirection * dE;
+ if (direction == kDownstream) dE = -dE;
+
+ cnst ECorrected = sqrt(E2) + dE;
cnst E2Corrected = ECorrected * ECorrected;
DataT corr = sqrt(p2 / (E2Corrected - fMass2));
@@ -1087,7 +1089,7 @@ namespace cbm::algo::ca
template<typename DataT>
void TrackFit<DataT>::EnergyLossCorrection(int atomicZ, DataTscal atomicA, DataTscal rho, DataTscal radLen,
- DataT radThick, DataT upstreamDirection)
+ DataT radThick, FitDirection direction)
{
cnst qp2cut(1. / (10. * 10.)); // 10 GeV cut
cnst qp02 = utils::max(fQp0 * fQp0, qp2cut);
@@ -1106,7 +1108,11 @@ namespace cbm::algo::ca
DataT dE = bethe * radThick * tr * radLen * rho;
- cnst E2Corrected = (sqrt(E2) + upstreamDirection * dE) * (sqrt(E2) + upstreamDirection * dE);
+ if (direction == kDownstream) dE = -dE;
+
+ cnst ECorrected = (sqrt(E2) + dE);
+ cnst E2Corrected = ECorrected * ECorrected;
+
DataT corr = sqrt(p2 / (E2Corrected - fMass2));
DataTmask ok = (corr == corr) && fMask;
corr = utils::iif(ok, corr, DataT(1.));
@@ -1123,8 +1129,8 @@ namespace cbm::algo::ca
DataT STEP = radThick * tr * radLen;
DataT EMASS(0.511 * 1e-3); // GeV
- DataT BETA = P / sqrt(E2Corrected);
- DataT GAMMA = sqrt(E2Corrected) / fMass;
+ DataT BETA = P / ECorrected;
+ DataT GAMMA = ECorrected / fMass;
// Calculate xi factor (KeV).
DataT XI = (DataT(153.5) * Z * STEP * RHO) / (A * BETA * BETA);
diff --git a/algo/ca/core/tracking/CaTrackFit.h b/algo/ca/core/tracking/CaTrackFit.h
index a4f6b1f8cf..9b06b7d44e 100644
--- a/algo/ca/core/tracking/CaTrackFit.h
+++ b/algo/ca/core/tracking/CaTrackFit.h
@@ -32,15 +32,23 @@ namespace cbm::algo::ca
class Hit;
+ enum FitDirection
+ {
+ kUpstream,
+ kDownstream
+ };
+
+ inline FitDirection operator!(FitDirection d) { return d == kUpstream ? kDownstream : kUpstream; }
+
/// Track fit utilities for the CA tracking based on the Kalman Filter
///
template<typename DataT>
class TrackFit {
+ public:
using DataTscal = utils::scaltype<DataT>;
using DataTmask = utils::masktype<DataT>;
- public:
TrackFit() = default;
TrackFit(const TrackParamBase<DataT>& t) { SetTrack(t); }
@@ -98,7 +106,7 @@ namespace cbm::algo::ca
DataT GetParticleMass2() const { return fMass2; }
/// set max extrapolation step [cm]
- void SetMaxExtrapolationStep(DataTscal step) { fMaxExtraplationStep = DataT(step); }
+ void SetMaxExtrapolationStep(double step) { fMaxExtraplationStep = DataT(step); }
/// get the particle mass
DataT GetMaxExtrapolationStep() const { return fMaxExtraplationStep; }
@@ -142,8 +150,8 @@ namespace cbm::algo::ca
/// apply energy loss correction to the track
/// \param radThick - radiation length of the material
- /// \param upstreamDirection - direction of the track (1 or 0)
- void EnergyLossCorrection(DataT radThick, DataT upstreamDirection);
+ /// \param direction - direction of the track
+ void EnergyLossCorrection(DataT radThick, FitDirection direction);
/// apply energy loss correction to the track
/// more accurate formula using material atomic numbers
@@ -152,9 +160,9 @@ namespace cbm::algo::ca
/// \param rho - density of the material
/// \param radLen - radiation length of the material
/// \param radThick - radiation length of the material
- /// \param upstreamDirection - direction of the track (1 or 0)
+ /// \param direction - direction of the track
void EnergyLossCorrection(int atomicZ, DataTscal atomicA, DataTscal rho, DataTscal radLen, DataT radThick,
- DataT upstreamDirection);
+ FitDirection direction);
/// apply multiple scattering correction to the track with the given Qp0
diff --git a/algo/ca/core/tracking/CaTrackFitter.cxx b/algo/ca/core/tracking/CaTrackFitter.cxx
index b9dc278886..e9d5d47a35 100644
--- a/algo/ca/core/tracking/CaTrackFitter.cxx
+++ b/algo/ca/core/tracking/CaTrackFitter.cxx
@@ -265,7 +265,7 @@ void TrackFitter::FitCaTracks(const ca::InputData& input, WindowData& wData)
fit.SetMask(initialised);
fit.Extrapolate(z[ista], fld1);
fit.MultipleScattering(fParameters.GetMaterialThickness(ista, tr.X(), tr.Y()));
- fit.EnergyLossCorrection(fParameters.GetMaterialThickness(ista, tr.X(), tr.Y()), fvec(1.f));
+ fit.EnergyLossCorrection(fParameters.GetMaterialThickness(ista, tr.X(), tr.Y()), kUpstream);
fit.SetMask(initialised && w[ista]);
fit.FilterXY(mxy[ista]);
@@ -369,7 +369,7 @@ void TrackFitter::FitCaTracks(const ca::InputData& input, WindowData& wData)
fit.SetMask(initialised);
fit.Extrapolate(z[ista], fld);
fit.MultipleScattering(fParameters.GetMaterialThickness(ista, tr.X(), tr.Y()));
- fit.EnergyLossCorrection(fParameters.GetMaterialThickness(ista, tr.X(), tr.Y()), fvec(-1.f));
+ fit.EnergyLossCorrection(fParameters.GetMaterialThickness(ista, tr.X(), tr.Y()), kDownstream);
fit.SetMask(initialised && w[ista]);
fit.FilterXY(mxy[ista]);
fit.FilterTime(time[ista], dt2[ista], fmask(sta[ista].timeInfo));
diff --git a/algo/ca/core/tracking/CaTripletConstructor.cxx b/algo/ca/core/tracking/CaTripletConstructor.cxx
index 5a1227b582..d1b70880c4 100644
--- a/algo/ca/core/tracking/CaTripletConstructor.cxx
+++ b/algo/ca/core/tracking/CaTripletConstructor.cxx
@@ -524,7 +524,7 @@ void TripletConstructor::FitTriplets()
// repeat several times in order to improve the precision
for (int iiter = 0; iiter < nIterations; ++iiter) {
- auto fitTrack = [&](int startIdx, int endIdx, int step, fscal energyLossSign) {
+ auto fitTrack = [&](int startIdx, int endIdx, int step, FitDirection direction) {
fit.SetQp0(T.Qp());
fit.Qp0()(fit.Qp0() > maxQp) = maxQp;
fit.Qp0()(fit.Qp0() < -maxQp) = -maxQp;
@@ -553,19 +553,19 @@ void TripletConstructor::FitTriplets()
fit.Extrapolate(z[ih], fld);
auto radThick = fParameters.GetMaterialThickness(ista[ih], T.X(), T.Y());
fit.MultipleScattering(radThick);
- fit.EnergyLossCorrection(radThick, fvec(energyLossSign));
+ fit.EnergyLossCorrection(radThick, direction);
fit.FilterXY(mxy[ih]);
fit.FilterTime(t[ih], dt2[ih], fmask(sta[ih].timeInfo));
}
};
- // Fit forward
- fitTrack(0, NHits, 1, -1.f);
+ // Fit downstream
+ fitTrack(0, NHits, 1, kDownstream);
if (iiter == nIterations - 1) break;
- // Fit backward
- fitTrack(NHits - 1, -1, -1, 1.f);
+ // Fit upstream
+ fitTrack(NHits - 1, -1, -1, kUpstream);
} // for iiter
fTracks_3[i3] = T;
diff --git a/reco/KF/CbmKfTrackFitter.cxx b/reco/KF/CbmKfTrackFitter.cxx
index afe921dc99..38e50ada95 100644
--- a/reco/KF/CbmKfTrackFitter.cxx
+++ b/reco/KF/CbmKfTrackFitter.cxx
@@ -456,7 +456,7 @@ void CbmKfTrackFitter::FilterFirstMeasurement(const FitNode& n)
}
-void CbmKfTrackFitter::AddMaterialEffects(CbmKfTrackFitter::FitNode& n, CbmKfTrackFitter::Direction direction)
+void CbmKfTrackFitter::AddMaterialEffects(CbmKfTrackFitter::FitNode& n, ca::FitDirection direction)
{
// add material effects
if (n.fMaterialLayer < 0) {
diff --git a/reco/KF/CbmKfTrackFitter.h b/reco/KF/CbmKfTrackFitter.h
index 73786e3216..61add5eff4 100644
--- a/reco/KF/CbmKfTrackFitter.h
+++ b/reco/KF/CbmKfTrackFitter.h
@@ -137,12 +137,7 @@ class CbmKfTrackFitter {
private:
void FilterFirstMeasurement(const FitNode& n);
- enum Direction
- {
- kUpstream,
- kDownstream
- };
- void AddMaterialEffects(FitNode& n, Direction direction);
+ void AddMaterialEffects(FitNode& n, ca::FitDirection direction);
// combine two tracks
bool Smooth(kf::TrackParamBase<double>& t1, const kf::TrackParamBase<double>& t2);
diff --git a/reco/KF/ParticleFitter/CbmL1PFFitter.cxx b/reco/KF/ParticleFitter/CbmL1PFFitter.cxx
index ad8f858906..1c6f3c6c42 100644
--- a/reco/KF/ParticleFitter/CbmL1PFFitter.cxx
+++ b/reco/KF/ParticleFitter/CbmL1PFFitter.cxx
@@ -374,7 +374,7 @@ void CbmL1PFFitter::Fit(std::vector<CbmStsTrack>& Tracks, const std::vector<CbmM
fit.Extrapolate(z[i], fld);
auto radThick = CbmL1::Instance()->fpAlgo->GetParameters().GetMaterialThickness(i, fit.Tr().X(), fit.Tr().Y());
fit.MultipleScattering(radThick);
- fit.EnergyLossCorrection(radThick, -fvec::One());
+ fit.EnergyLossCorrection(radThick, kDownstream);
fit.SetMask(initialised && w[i]);
fit.FilterXY(mxy[i]);
@@ -435,7 +435,7 @@ void CbmL1PFFitter::Fit(std::vector<CbmStsTrack>& Tracks, const std::vector<CbmM
fit.Extrapolate(z[i], fld);
auto radThick = CbmL1::Instance()->fpAlgo->GetParameters().GetMaterialThickness(i, fit.Tr().X(), fit.Tr().Y());
fit.MultipleScattering(radThick);
- fit.EnergyLossCorrection(radThick, fvec::One());
+ fit.EnergyLossCorrection(radThick, kUpstream);
fit.SetMask(initialised && w[i]);
fit.FilterXY(mxy[i]);
@@ -606,7 +606,7 @@ void CbmL1PFFitter::GetChiToVertex(vector<CbmStsTrack>& Tracks, vector<PFFieldRe
fit.Extrapolate(zSta[iSt], fld);
auto radThick = CbmL1::Instance()->fpAlgo->GetParameters().GetMaterialThickness(iSt, fit.Tr().X(), fit.Tr().Y());
fit.MultipleScattering(radThick);
- fit.EnergyLossCorrection(radThick, fvec::One());
+ fit.EnergyLossCorrection(radThick, kUpstream);
}
fit.SetMask(fmask::One());
fit.Extrapolate(primVtx.GetRefZ(), fld);
@@ -615,7 +615,7 @@ void CbmL1PFFitter::GetChiToVertex(vector<CbmStsTrack>& Tracks, vector<PFFieldRe
constexpr float targetRadThick = 3.73e-2f * 2; // 250 mum Gold
fit.MultipleScattering(targetRadThick);
- fit.EnergyLossCorrection(targetRadThick, fvec::One());
+ fit.EnergyLossCorrection(targetRadThick, kUpstream);
Double_t Cv[3] = {primVtx.GetCovMatrix()[0], primVtx.GetCovMatrix()[1], primVtx.GetCovMatrix()[2]};
@@ -626,8 +626,8 @@ void CbmL1PFFitter::GetChiToVertex(vector<CbmStsTrack>& Tracks, vector<PFFieldRe
c[1] += fvec(Cv[1]);
c[2] += fvec(Cv[2]);
fvec d = c[0] * c[2] - c[1] * c[1];
- fvec chi = sqrt(abs(fvec(0.5) * (dx * dx * c[0] - fvec(2.) * dx * dy * c[1] + dy * dy * c[2]) / d));
- chi.setZero(abs(d) < fvec(1.e-20));
+ fvec chi = sqrt(utils::fabs(fvec(0.5) * (dx * dx * c[0] - fvec(2.) * dx * dy * c[1] + dy * dy * c[2]) / d));
+ chi.setZero(utils::fabs(d) < fvec(1.e-20));
for (int iVec = 0; iVec < nTracks_SIMD; iVec++) {
chiToVtx.push_back(chi[iVec]);
diff --git a/reco/L1/CbmL1Performance.cxx b/reco/L1/CbmL1Performance.cxx
index ff8367828d..7ad879550b 100644
--- a/reco/L1/CbmL1Performance.cxx
+++ b/reco/L1/CbmL1Performance.cxx
@@ -1310,7 +1310,7 @@ void CbmL1::TrackFitPerformance()
// LOG(info) << iSta << " " << dir;
auto radThick = fpAlgo->GetParameters().GetMaterialThickness(iSta, fit.Tr().GetX(), fit.Tr().GetY());
fit.MultipleScattering(radThick);
- fit.EnergyLossCorrection(radThick, fvec::One());
+ fit.EnergyLossCorrection(radThick, kUpstream);
}
}
if (mc.GetStartZ() != tr.GetZ()[0]) continue;
@@ -1380,7 +1380,7 @@ void CbmL1::TrackFitPerformance()
fit.Extrapolate(fpAlgo->GetParameters().GetStation(iSta).fZ, fld);
auto radThick = fpAlgo->GetParameters().GetMaterialThickness(iSta, fit.Tr().GetX(), fit.Tr().GetY());
fit.MultipleScattering(radThick);
- fit.EnergyLossCorrection(radThick, fvec::One());
+ fit.EnergyLossCorrection(radThick, kUpstream);
}
fit.Extrapolate(mc.GetStartZ(), fld);
}
diff --git a/reco/L1/qa/CbmCaTrackTypeQa.cxx b/reco/L1/qa/CbmCaTrackTypeQa.cxx
index a12ce6b99b..c552a74c6d 100644
--- a/reco/L1/qa/CbmCaTrackTypeQa.cxx
+++ b/reco/L1/qa/CbmCaTrackTypeQa.cxx
@@ -109,13 +109,13 @@ void TrackTypeQa::Init()
//
// ** Distributions of MC tracks **
//
- fph_mc_pMC = MakeQaObject<TH1F>("mc_pMC", "", kBinsP, kLoP, kUpP);
- fph_mc_etaMC = MakeQaObject<TH1F>("mc_etaMC", "", kBinsETA, kLoETA, kUpETA);
- fph_mc_yMC = MakeQaObject<TH1F>("mc_yMC", "", kBinsY, kLoY, kUpY);
- fph_mc_ptMC_yMC = MakeQaObject<TH2F>("mc_ptMC_yMC", "", kBinsY, kLoY, kUpY, kBinsPT, kLoPT, kUpPT);
- fph_mc_ptMC = MakeQaObject<TH1F>("mc_ptMC", "", kBinsPT, kLoPT, kUpPT);
- fph_mc_phiMC = MakeQaObject<TH1F>("mc_phiMC", "", kBinsPHI, kLoPHI, kUpPHI);
- fph_mc_thetaMC = MakeQaObject<TH1F>("mc_thetaMC", "", kBinsTHETA, kLoTHETA, kUpTHETA);
+ fph_mc_pMC = MakeQaObject<TH1F>("mc_pMC", "", kBinsP, kLoP, kUpP);
+ fph_mc_etaMC = MakeQaObject<TH1F>("mc_etaMC", "", kBinsETA, kLoETA, kUpETA);
+ fph_mc_yMC = MakeQaObject<TH1F>("mc_yMC", "", kBinsY, kLoY, kUpY);
+ fph_mc_ptMC_yMC = MakeQaObject<TH2F>("mc_ptMC_yMC", "", kBinsY, kLoY, kUpY, kBinsPT, kLoPT, kUpPT);
+ fph_mc_ptMC = MakeQaObject<TH1F>("mc_ptMC", "", kBinsPT, kLoPT, kUpPT);
+ fph_mc_phiMC = MakeQaObject<TH1F>("mc_phiMC", "", kBinsPHI, kLoPHI, kUpPHI);
+ fph_mc_thetaMC = MakeQaObject<TH1F>("mc_thetaMC", "", kBinsTHETA, kLoTHETA, kUpTHETA);
fph_mc_thetaMC_phiMC =
MakeQaObject<TH2F>("mc_thetaMC_phiMC", "", kBinsPHI, kLoPHI, kUpPHI, kBinsTHETA, kLoTHETA, kUpTHETA);
fph_mc_txMC = MakeQaObject<TH1F>("mc_txMC", "", kBinsTX, kLoTX, kUpTX);
@@ -290,18 +290,18 @@ void TrackTypeQa::FillRecoTrack(int iTrkReco, double weight)
MCPoint mcTrkVertex = mcTrack.GetVertexPoint();
TrackParamV trPar(recoTrack);
fTrackFit.SetTrack(trPar);
- fTrackFit.Extrapolate(mcTrkVertex.GetZ(), fFieldRegion);
- // Add material
+
+ // Add material effects
int iStFst = (*fpvHits)[recoTrack.GetFirstHitIndex()].GetStationId();
double dZ = mcTrkVertex.GetZ() - fpParameters->GetStation(iStFst).fZ[0];
- int direction = static_cast<int>(dZ / std::fabs(dZ));
+ int direction = (dZ > 0.) ? 1 : -1;
for (int iSt = iStFst; (iSt >= 0) && (iSt < fpParameters->GetNstationsActive())
&& (direction * (mcTrkVertex.GetZ() - fpParameters->GetStation(iSt).fZ[0]) > 0);
iSt += direction) {
fTrackFit.Extrapolate(fpParameters->GetStation(iSt).fZ, fFieldRegion);
auto radLength = fpParameters->GetMaterialThickness(iSt, fTrackFit.Tr().GetX(), fTrackFit.Tr().GetY());
fTrackFit.MultipleScattering(radLength);
- fTrackFit.EnergyLossCorrection(radLength, fvec::One());
+ fTrackFit.EnergyLossCorrection(radLength, (direction > 0) ? kDownstream : kUpstream);
}
fTrackFit.Extrapolate(mcTrkVertex.GetZ(), fFieldRegion);
const TrackParamV& trParExtr = fTrackFit.Tr();
--
GitLab