diff --git a/algo/ca/core/data/CaWindowData.h b/algo/ca/core/data/CaWindowData.h
index d0c091d7e83a8bf90c111eb876a8b4e616b5046e..0af05c83fcbde59432d0127632d4e4b16783927d 100644
--- a/algo/ca/core/data/CaWindowData.h
+++ b/algo/ca/core/data/CaWindowData.h
@@ -10,8 +10,11 @@
#pragma once
#include "CaConstants.h"
+#include "CaField.h"
#include "CaGrid.h"
#include "CaHit.h"
+#include "CaIteration.h"
+#include "CaMeasurementXy.h"
#include "CaTrack.h"
#include <array>
@@ -129,6 +132,26 @@ namespace cbm::algo::ca
/// \param iSt Index of tracking station
[[gnu::always_inline]] const Vector<HitIndex_t>& FullDsHitIndices(int iSt) const { return fvFullDsHitIndices[iSt]; }
+ /// \brief Accesses current iteration
+ [[gnu::always_inline]] void SetCurrentIteration(const Iteration* ptr) { fpCurrentIteration = ptr; }
+
+ /// \brief Accesses current iteration
+ [[gnu::always_inline]] const Iteration* CurrentIteration() const { return fpCurrentIteration; }
+
+ fvec& TargX() { return fTargX; }
+ fvec& TargY() { return fTargY; }
+ fvec& TargZ() { return fTargZ; }
+
+ fvec TargX() const { return fTargX; }
+ fvec TargY() const { return fTargY; }
+ fvec TargZ() const { return fTargZ; }
+
+ ca::FieldValue<fvec>& TargB() { return fTargB; }
+ const ca::FieldValue<fvec>& TargB() const { return fTargB; }
+
+ ca::MeasurementXy<fvec>& TargetMeasurement() { return fTargetMeasurement; }
+ const ca::MeasurementXy<fvec>& TargetMeasurement() const { return fTargetMeasurement; }
+
private:
static constexpr int kMaxNofStations = constants::size::MaxNstations; ///< Alias to max number of stations
@@ -162,5 +185,15 @@ namespace cbm::algo::ca
/// \brief Map of hit indices from sub-TS to full dataset
std::array<Vector<HitIndex_t>, kMaxNofStations> fvFullDsHitIndices{"WindowData::fvFullDSHitIndex"};
+ /// \brief Current track-finder iteration
+ const Iteration* fpCurrentIteration = nullptr;
+
+ fvec fTargX{constants::Undef<fvec>}; ///< target position x coordinate for the current iteration (modifiable)
+ fvec fTargY{constants::Undef<fvec>}; ///< target position y coordinate for the current iteration (modifiable)
+ fvec fTargZ{constants::Undef<fvec>}; ///< target position z coordinate for the current iteration (modifiable)
+
+ ca::FieldValue<fvec> fTargB _fvecalignment{}; // field in the target point (modifiable, do not touch!!)
+ ca::MeasurementXy<fvec> fTargetMeasurement _fvecalignment{}; // target constraint [cm]
+
} _fvecalignment;
} // namespace cbm::algo::ca
diff --git a/algo/ca/core/pars/CaIteration.h b/algo/ca/core/pars/CaIteration.h
index 5512083087e1d470a892155ea3eb0921f6def191..0704d3c103950f4e9d21bb79e0fc688d2de71869 100644
--- a/algo/ca/core/pars/CaIteration.h
+++ b/algo/ca/core/pars/CaIteration.h
@@ -145,6 +145,11 @@ namespace cbm::algo::ca
void SetMaxStationGap(int nSkipped) { fMaxStationGap = nSkipped; }
/// \brief Sets correction for accounting overlaping and iff z
+ ///
+ /// the maxDZ correction is set in order to take into account overlaping and iff z.
+ /// The reason is that low momentum tracks are too curved and goes not from target direction.
+ /// That's why sort by hit_y/hit_z is not work idealy. If sort by y then it is max diff between
+ /// same station's modules (~0.4cm).
void SetMaxDZ(float input) { fMaxDZ = input; }
/// \brief Sets max considered q/p for tracks
diff --git a/algo/ca/core/tracking/CaFramework.h b/algo/ca/core/tracking/CaFramework.h
index 007484d946034ad6504b462f2d3cb8990944aee4..0d4c18ad3c2c7d663dcd41f37cf55ac677b70e5f 100644
--- a/algo/ca/core/tracking/CaFramework.h
+++ b/algo/ca/core/tracking/CaFramework.h
@@ -178,8 +178,6 @@ namespace cbm::algo::ca
void PrintHits();
- float GetMaxInvMom() const { return fMaxInvMom[0]; }
-
/// \brief Gets monitor data
const TrackingMonitorData& GetMonitorData() const { return fMonitorData; }
@@ -245,8 +243,6 @@ namespace cbm::algo::ca
/// --- data used during finding iterations
unsigned int fCurrentIterationIndex{0}; ///< index of the corrent iteration, needed for debug output
- const Iteration* fpCurrentIteration = nullptr; ///< pointer to the current CA track finder iteration
-
public:
ca::WindowData fWData{};
@@ -261,40 +257,8 @@ namespace cbm::algo::ca
/// ----- Different parameters of CATrackFinder -----
- Tindex fFirstCAstation{0}; //first station used in CA
-
- // fNFindIterations - set number of interation for trackfinding
- // itetation of finding:
-
- int fNFindIterations = -1; // TODO investigate kAllPrimJumpIter & kAllSecJumpIter
-
- float fTrackChi2Cut{10.f};
- float fTripletFinalChi2Cut{10.f};
- float fTripletChi2Cut{5.f}; // cut for selecting triplets before collecting tracks.per one DoF
- float fDoubletChi2Cut{5.f};
- float fTimeCut1{0.f}; // TODO: please, specify "1" and "2" (S.Zharko)
- float fTimeCut2{0.f};
-
- /// correction in order to take into account overlaping and iff z. if sort by y then it is max diff between same station's modules (~0.4cm)
- fvec fMaxDZ{constants::Undef<fvec>};
-
- /// parameters which are different for different iterations. Set in the begin of CA Track Finder
-
- float fPickGather{constants::Undef<fvec>}; ///< same for attaching additional hits to track
- float fTripletLinkChi2{
- constants::Undef<fvec>}; ///< (dp2/dp_error2 < fTripletLinkChi2) => triplets are neighbours
- fvec fMaxInvMom{constants::Undef<fvec>}; ///< max considered q/p for tracks
- fvec fMaxSlopePV{constants::Undef<fvec>}; ///< max slope (tx\ty) in prim vertex
- float fMaxSlope{constants::Undef<fvec>}; ///< max slope (tx\ty) in 3d hit position of a triplet
-
- fvec fTargX{constants::Undef<fvec>}; ///< target position x coordinate for the current iteration (modifiable)
- fvec fTargY{constants::Undef<fvec>}; ///< target position y coordinate for the current iteration (modifiable)
- fvec fTargZ{constants::Undef<fvec>}; ///< target position z coordinate for the current iteration (modifiable)
-
bool fIsTargetField{false}; ///< is the magnetic field present at the target
- ca::FieldValue<fvec> fTargB _fvecalignment{}; // field in the target point (modifiable, do not touch!!)
- ca::MeasurementXy<fvec> TargetMeasurement _fvecalignment{}; // target constraint [cm]
} _fvecalignment;
diff --git a/algo/ca/core/tracking/CaTrackExtender.cxx b/algo/ca/core/tracking/CaTrackExtender.cxx
index b44b8b28e6e824d9cb1d19a879cf2d25f8d33a0f..efd6418ea3f3e6398331b11907eb183654ac4918 100644
--- a/algo/ca/core/tracking/CaTrackExtender.cxx
+++ b/algo/ca/core/tracking/CaTrackExtender.cxx
@@ -192,10 +192,11 @@ void TrackExtender::FindMoreHits(ca::Branch& t, TrackParamV& Tout, const bool up
int ista = ista2 + 2 * step; // skip one station. if there would be hit it has to be found on previous stap
- if (ista2 == frAlgo.fFirstCAstation) ista = ista2 + step;
-
- const fvec pickGather2 = frAlgo.fPickGather * frAlgo.fPickGather;
+ if (ista2 == frAlgo.fWData.CurrentIteration()->GetFirstStationIndex()) ista = ista2 + step;
+ const float pickGather = frAlgo.fWData.CurrentIteration()->GetPickGather();
+ const fvec pickGather2 = pickGather * pickGather;
+ const fvec maxDZ = frAlgo.fWData.CurrentIteration()->GetMaxDZ();
for (; (ista < frAlgo.GetParameters().GetNstationsActive()) && (ista >= 0); ista += step) { // CHECKME why ista2?
const ca::Station<fvec>& sta = frAlgo.GetParameters().GetStation(ista);
@@ -209,8 +210,8 @@ void TrackExtender::FindMoreHits(ca::Branch& t, TrackParamV& Tout, const bool up
const auto& grid = frAlgo.fWData.Grid(ista);
ca::GridArea area(grid, tr.X()[0], tr.Y()[0],
- (sqrt(frAlgo.fPickGather * tr.C00()) + grid.GetMaxRangeX() + frAlgo.fMaxDZ * abs(tr.Tx()))[0],
- (sqrt(frAlgo.fPickGather * tr.C11()) + grid.GetMaxRangeY() + frAlgo.fMaxDZ * abs(tr.Ty()))[0]);
+ (sqrt(pickGather * tr.C00()) + grid.GetMaxRangeX() + maxDZ * abs(tr.Tx()))[0],
+ (sqrt(pickGather * tr.C11()) + grid.GetMaxRangeY() + maxDZ * abs(tr.Ty()))[0]);
if (frAlgo.GetParameters().DevIsIgnoreHitSearchAreas()) {
area.DoLoopOverEntireGrid();
diff --git a/algo/ca/core/tracking/CaTrackFinderWindow.cxx b/algo/ca/core/tracking/CaTrackFinderWindow.cxx
index 5fc39952d04a696b6804d6826df688d2a154ebed..0336d2b10c0489bd7eb64511d9d92003221cd7af 100644
--- a/algo/ca/core/tracking/CaTrackFinderWindow.cxx
+++ b/algo/ca/core/tracking/CaTrackFinderWindow.cxx
@@ -72,6 +72,7 @@ bool TrackFinderWindow::checkTripletMatch(const ca::Triplet& l, const ca::Triple
if (r.GetLSta() != l.GetMSta()) return false;
+ const float tripletLinkChi2 = frAlgo.fWData.CurrentIteration()->GetTripletLinkChi2();
if (r.IsMomentumFitted()) {
assert(l.IsMomentumFitted());
@@ -81,7 +82,7 @@ bool TrackFinderWindow::checkTripletMatch(const ca::Triplet& l, const ca::Triple
if (!std::isfinite(dqp)) return false;
if (!std::isfinite(Cqp)) return false;
- if (dqp * dqp > frAlgo.fTripletLinkChi2 * Cqp) {
+ if (dqp * dqp > tripletLinkChi2 * Cqp) {
return false; // bad neighbour // CHECKME why do we need recheck it?? (it really change result)
}
dchi2 = dqp * dqp / Cqp;
@@ -100,8 +101,8 @@ bool TrackFinderWindow::checkTripletMatch(const ca::Triplet& l, const ca::Triple
if (!std::isfinite(Ctx)) return false;
if (!std::isfinite(Cty)) return false;
- if (dty * dty > frAlgo.fTripletLinkChi2 * Cty) return false;
- if (dtx * dtx > frAlgo.fTripletLinkChi2 * Ctx) return false;
+ if (dty * dty > tripletLinkChi2 * Cty) return false;
+ if (dtx * dtx > tripletLinkChi2 * Ctx) return false;
//dchi2 = 0.5f * (dtx * dtx / Ctx + dty * dty / Cty);
dchi2 = 0.;
@@ -170,9 +171,6 @@ void TrackFinderWindow::ReadWindowData()
void TrackFinderWindow::CaTrackFinderSlice()
{
- frAlgo.fNFindIterations = frAlgo.GetParameters().GetNcaIterations();
-
-
/********************************/ /**
* CATrackFinder routine setting
***********************************/
@@ -249,16 +247,13 @@ void TrackFinderWindow::CaTrackFinderSlice()
// ---- Loop over Track Finder iterations ----------------------------------------------------------------//
- // TODO: replace with exception
- assert(frAlgo.fNFindIterations == (int) frAlgo.GetParameters().GetCAIterations().size());
-
for (frAlgo.fCurrentIterationIndex = 0;
frAlgo.fCurrentIterationIndex < frAlgo.GetParameters().GetCAIterations().size();
++frAlgo.fCurrentIterationIndex) {
// current iteration of the CA tracker
const auto& caIteration = frAlgo.GetParameters().GetCAIterations()[frAlgo.fCurrentIterationIndex];
- frAlgo.fpCurrentIteration = &caIteration;
+ frAlgo.fWData.SetCurrentIteration(&caIteration);
if (frAlgo.fCurrentIterationIndex > 0) {
for (int ista = 0; ista < frAlgo.GetParameters().GetNstationsActive(); ++ista) {
@@ -280,48 +275,31 @@ void TrackFinderWindow::CaTrackFinderSlice()
// #pragma omp task
{
// --- SET PARAMETERS FOR THE ITERATION ---
-
- frAlgo.fFirstCAstation = caIteration.GetFirstStationIndex();
- frAlgo.fTrackChi2Cut = caIteration.GetTrackChi2Cut();
- frAlgo.fDoubletChi2Cut = caIteration.GetDoubletChi2Cut(); //11.3449 * 2.f / 3.f; // prob = 0.1
- frAlgo.fTripletChi2Cut = caIteration.GetTripletChi2Cut(); //21.1075; // prob = 0.01%
- frAlgo.fTripletFinalChi2Cut = caIteration.GetTripletFinalChi2Cut();
- frAlgo.fPickGather = caIteration.GetPickGather(); //3.0;
- frAlgo.fTripletLinkChi2 = caIteration.GetTripletLinkChi2(); //5.0;
- frAlgo.fMaxInvMom = caIteration.GetMaxQp(); //1.0 / 0.5; // max considered q/p
- frAlgo.fMaxSlopePV = caIteration.GetMaxSlopePV(); //1.1;
- frAlgo.fMaxSlope = caIteration.GetMaxSlope(); //2.748; // corresponds to 70 grad
-
// define the target
- frAlgo.fTargX = frAlgo.GetParameters().GetTargetPositionX();
- frAlgo.fTargY = frAlgo.GetParameters().GetTargetPositionY();
- frAlgo.fTargZ = frAlgo.GetParameters().GetTargetPositionZ();
+ frAlgo.fWData.TargX() = frAlgo.GetParameters().GetTargetPositionX();
+ frAlgo.fWData.TargY() = frAlgo.GetParameters().GetTargetPositionY();
+ frAlgo.fWData.TargZ() = frAlgo.GetParameters().GetTargetPositionZ();
fscal SigmaTargetX = caIteration.GetTargetPosSigmaX();
fscal SigmaTargetY = caIteration.GetTargetPosSigmaY(); // target constraint [cm]
// Select magnetic field. For primary tracks - fVtxFieldValue, for secondary tracks - st.fieldSlice
if (caIteration.GetPrimaryFlag()) {
- frAlgo.fTargB = frAlgo.GetParameters().GetVertexFieldValue();
+ frAlgo.fWData.TargB() = frAlgo.GetParameters().GetVertexFieldValue();
}
else {
- frAlgo.fTargB = frAlgo.GetParameters().GetStation(0).fieldSlice.GetFieldValue(0, 0);
+ frAlgo.fWData.TargB() = frAlgo.GetParameters().GetStation(0).fieldSlice.GetFieldValue(0, 0);
} // NOTE: calculates field frAlgo.fTargB in the center of 0th station
- frAlgo.fIsTargetField = (fabs(frAlgo.fTargB.y[0]) > 0.001);
-
- frAlgo.TargetMeasurement.SetX(frAlgo.fTargX);
- frAlgo.TargetMeasurement.SetY(frAlgo.fTargY);
- frAlgo.TargetMeasurement.SetDx2(SigmaTargetX * SigmaTargetX);
- frAlgo.TargetMeasurement.SetDxy(0);
- frAlgo.TargetMeasurement.SetDy2(SigmaTargetY * SigmaTargetY);
- frAlgo.TargetMeasurement.SetNdfX(1);
- frAlgo.TargetMeasurement.SetNdfY(1);
+ frAlgo.fIsTargetField = (fabs(frAlgo.fWData.TargB().y[0]) > 0.001);
- /// Set correction in order to take into account overlaping and iff z.
- /// The reason is that low momentum tracks are too curved and goes not from target direction. That's why sort by hit_y/hit_z is not work idealy
- /// If sort by y then it is max diff between same station's modules (~0.4cm)
- frAlgo.fMaxDZ = caIteration.GetMaxDZ(); //0;
+ frAlgo.fWData.TargetMeasurement().SetX(frAlgo.fWData.TargX());
+ frAlgo.fWData.TargetMeasurement().SetY(frAlgo.fWData.TargY());
+ frAlgo.fWData.TargetMeasurement().SetDx2(SigmaTargetX * SigmaTargetX);
+ frAlgo.fWData.TargetMeasurement().SetDxy(0);
+ frAlgo.fWData.TargetMeasurement().SetDy2(SigmaTargetY * SigmaTargetY);
+ frAlgo.fWData.TargetMeasurement().SetNdfX(1);
+ frAlgo.fWData.TargetMeasurement().SetNdfY(1);
}
}
@@ -334,14 +312,15 @@ void TrackFinderWindow::CaTrackFinderSlice()
// indices of the two neighbouring station, taking into account allowed gaps
std::vector<std::pair<int, int>> staPattern;
- for (int gap = 0; gap <= frAlgo.fpCurrentIteration->GetMaxStationGap(); gap++) {
+ for (int gap = 0; gap <= frAlgo.fWData.CurrentIteration()->GetMaxStationGap(); gap++) {
for (int i = 0; i <= gap; i++) {
staPattern.push_back(std::make_pair(1 + i, 2 + gap));
}
}
- for (int istal = frAlgo.GetParameters().GetNstationsActive() - 2; istal >= frAlgo.fFirstCAstation;
- istal--) { // start with downstream chambers
+ const int iStFirst = frAlgo.fWData.CurrentIteration()->GetFirstStationIndex();
+ for (int istal = frAlgo.GetParameters().GetNstationsActive() - 2; istal >= iStFirst; istal--) {
+ // start with downstream chambers
const auto& grid = frAlgo.fWData.Grid(istal);
Tindex nGridEntriesL = grid.GetEntries().size();
for (Tindex iel = 0; iel < nGridEntriesL; ++iel) {
@@ -362,9 +341,8 @@ void TrackFinderWindow::CaTrackFinderSlice()
// search for neighbouring triplets
frAlgo.fMonitorData.StartTimer(ETimer::NeighboringTripletSearch);
- for (int istal = frAlgo.GetParameters().GetNstationsActive() - 2; istal >= frAlgo.fFirstCAstation;
- istal--) { // start with downstream chambers
-
+ for (int istal = frAlgo.GetParameters().GetNstationsActive() - 2; istal >= iStFirst; istal--) {
+ // start with downstream chambers
for (unsigned int it = 0; it < fvTriplets[istal].size(); ++it) {
ca::Triplet& tr = fvTriplets[istal][it];
@@ -380,7 +358,7 @@ void TrackFinderWindow::CaTrackFinderSlice()
if (nNeighbours > 0) {
assert((int) neighStation >= istal + 1
- && (int) neighStation <= istal + 1 + frAlgo.fpCurrentIteration->GetMaxStationGap());
+ && (int) neighStation <= istal + 1 + frAlgo.fWData.CurrentIteration()->GetMaxStationGap());
}
unsigned char level = 0;
@@ -422,7 +400,7 @@ void TrackFinderWindow::CaTrackFinderSlice()
// min level to start triplet. So min track length = min_level+3.
int min_level = std::min(caIteration.GetMinNhits(), caIteration.GetMinNhitsStation0()) - 3;
- if (frAlgo.fpCurrentIteration->GetTrackFromTripletsFlag()) {
+ if (frAlgo.fWData.CurrentIteration()->GetTrackFromTripletsFlag()) {
min_level = 0;
}
@@ -454,8 +432,8 @@ void TrackFinderWindow::CaTrackFinderSlice()
//== Loop over triplets with the required level, find and store track candidates
- for (int istaF = frAlgo.fFirstCAstation;
- istaF <= frAlgo.GetParameters().GetNstationsActive() - 3 - firstTripletLevel; ++istaF) {
+ for (int istaF = iStFirst; istaF <= frAlgo.GetParameters().GetNstationsActive() - 3 - firstTripletLevel;
+ ++istaF) {
if (--nlevel == 0) break; //TODO: SG: this is not needed
@@ -469,12 +447,12 @@ void TrackFinderWindow::CaTrackFinderSlice()
// skip track candidates that are too short
- int minNhits = frAlgo.fpCurrentIteration->GetMinNhits();
+ int minNhits = frAlgo.fWData.CurrentIteration()->GetMinNhits();
if (fstTripLHit.Station() == 0) {
- minNhits = frAlgo.fpCurrentIteration->GetMinNhitsStation0();
+ minNhits = frAlgo.fWData.CurrentIteration()->GetMinNhitsStation0();
}
- if (frAlgo.fpCurrentIteration->GetTrackFromTripletsFlag()) {
+ if (frAlgo.fWData.CurrentIteration()->GetTrackFromTripletsFlag()) {
minNhits = 0;
}
@@ -519,9 +497,9 @@ void TrackFinderWindow::CaTrackFinderSlice()
best_tr.SetChi2(best_tr.Chi2() / ndf);
if (frAlgo.GetParameters().GetGhostSuppression()) {
if (3 == best_tr.NofHits()) {
- if (!frAlgo.fpCurrentIteration->GetPrimaryFlag() && (istaF != 0))
+ if (!frAlgo.fWData.CurrentIteration()->GetPrimaryFlag() && (istaF != 0))
continue; // too /*short*/ non-MAPS track
- if (frAlgo.fpCurrentIteration->GetPrimaryFlag() && (best_tr.Chi2() > 5.0)) continue;
+ if (frAlgo.fWData.CurrentIteration()->GetPrimaryFlag() && (best_tr.Chi2() > 5.0)) continue;
}
}
fvTrackCandidates.push_back(best_tr);
@@ -647,7 +625,7 @@ void TrackFinderWindow::CaTrackFinderSlice()
if (!tr.IsAlive()) continue;
- if (frAlgo.fpCurrentIteration->GetExtendTracksFlag()) {
+ if (frAlgo.fWData.CurrentIteration()->GetExtendTracksFlag()) {
if (tr.NofHits() < frAlgo.GetParameters().GetNstationsActive()) {
frAlgo.fTrackExtender.ExtendBranch(tr);
}
@@ -739,7 +717,8 @@ void TrackFinderWindow::CAFindTrack(int ista, ca::Branch& best_tr, const ca::Tri
ndf = curr_tr.NofHits() * 2 - 4;
}
- if (curr_tr.Chi2() > frAlgo.fTrackChi2Cut * ndf) {
+
+ if (curr_tr.Chi2() > frAlgo.fWData.CurrentIteration()->GetTrackChi2Cut() * ndf) {
return;
}
@@ -809,7 +788,7 @@ void TrackFinderWindow::CAFindTrack(int ista, ca::Branch& best_tr, const ca::Tri
ndf = new_L; // TODO: SG: 2 * (new_L + 2) - 5
}
- if (new_chi2 > frAlgo.fTrackChi2Cut * ndf) continue;
+ if (new_chi2 > frAlgo.fWData.CurrentIteration()->GetTrackChi2Cut() * ndf) continue;
// add new hit
new_tr[ista] = curr_tr;
diff --git a/algo/ca/core/tracking/CaTrackFitter.cxx b/algo/ca/core/tracking/CaTrackFitter.cxx
index 7167e5a0950907e4b90456a8f8623473c6358453..d25098a1ecf0619cb2229211ee24bd66b7fa3f1e 100644
--- a/algo/ca/core/tracking/CaTrackFitter.cxx
+++ b/algo/ca/core/tracking/CaTrackFitter.cxx
@@ -279,7 +279,7 @@ void TrackFitter::FitCaTracks()
{
fitpv.SetMask(fmask::One());
- ca::MeasurementXy<fvec> vtxInfo = frAlgo.TargetMeasurement;
+ ca::MeasurementXy<fvec> vtxInfo = frAlgo.fWData.TargetMeasurement();
vtxInfo.SetDx2(1.e-8);
vtxInfo.SetDxy(0.);
vtxInfo.SetDy2(1.e-8);
diff --git a/algo/ca/core/tracking/CaTripletConstructor.cxx b/algo/ca/core/tracking/CaTripletConstructor.cxx
index a3535c533f38ca846b6f1c38d157dd1de15b029e..031f0297a1295700538681e08cd2713bca6abc40 100644
--- a/algo/ca/core/tracking/CaTripletConstructor.cxx
+++ b/algo/ca/core/tracking/CaTripletConstructor.cxx
@@ -79,7 +79,7 @@ void TripletConstructor::InitStations(int istal, int istam, int istar)
}
fFit.SetParticleMass(fAlgo->GetDefaultParticleMass());
- if (fAlgo->fpCurrentIteration->GetElectronFlag()) {
+ if (fAlgo->fWData.CurrentIteration()->GetElectronFlag()) {
fFit.SetParticleMass(constants::phys::ElectronMass);
}
fFit.SetMask(fmask::One());
@@ -104,7 +104,7 @@ void TripletConstructor::CreateTripletsForHit(int istal, int istam, int istar, c
/// Get the field approximation. Add the target to parameters estimation.
/// Propagaete to the middle station of a triplet.
- fFit.SetQp0(fAlgo->fMaxInvMom);
+ fFit.SetQp0(fAlgo->fWData.CurrentIteration()->GetMaxQp());
ca::FieldValue<fvec> lB, mB, cB, rB _fvecalignment;
ca::FieldValue<fvec> l_B_global, centB_global _fvecalignment;
@@ -117,10 +117,10 @@ void TripletConstructor::CreateTripletsForHit(int istal, int istam, int istar, c
fvec zl = hitl.Z();
fvec time = hitl.T();
fvec timeEr2 = hitl.dT2();
- const fvec dzli = 1.f / (zl - fAlgo->fTargZ);
+ const fvec dzli = 1.f / (zl - fAlgo->fWData.TargZ());
- const fvec tx = (xl - fAlgo->fTargX) * dzli;
- const fvec ty = (yl - fAlgo->fTargY) * dzli;
+ const fvec tx = (xl - fAlgo->fWData.TargX()) * dzli;
+ const fvec ty = (yl - fAlgo->fWData.TargY()) * dzli;
TrackParamV& T = fFit.Tr();
@@ -133,18 +133,20 @@ void TripletConstructor::CreateTripletsForHit(int istal, int istam, int istar, c
T.Time() = time;
T.Vi() = constants::phys::SpeedOfLightInv;
- fvec txErr2 = fAlgo->fMaxSlopePV * fAlgo->fMaxSlopePV / fvec(9.);
- fvec qpErr2 = fAlgo->fMaxInvMom * fAlgo->fMaxInvMom / fvec(9.);
+ const fvec maxSlopePV = fAlgo->fWData.CurrentIteration()->GetMaxSlopePV();
+ const fvec maxQp = fAlgo->fWData.CurrentIteration()->GetMaxQp();
+ fvec txErr2 = maxSlopePV * maxSlopePV / fvec(9.);
+ fvec qpErr2 = maxQp * maxQp / fvec(9.);
T.ResetErrors(1., 1., txErr2, txErr2, qpErr2, (staL().timeInfo ? timeEr2 : 1.e6), 1.e10);
- T.InitVelocityRange(1. / fAlgo->fpCurrentIteration->GetMaxQp());
+ T.InitVelocityRange(1. / fAlgo->fWData.CurrentIteration()->GetMaxQp());
T.C00() = hitl.dX2();
T.C10() = hitl.dXY();
T.C11() = hitl.dY2();
- T.Ndf() = (fAlgo->fpCurrentIteration->GetPrimaryFlag()) ? fvec(2.) : fvec(0.);
+ T.Ndf() = (fAlgo->fWData.CurrentIteration()->GetPrimaryFlag()) ? fvec(2.) : fvec(0.);
T.NdfTime() = (staL().timeInfo ? 0 : -1);
// 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
@@ -155,7 +157,7 @@ void TripletConstructor::CreateTripletsForHit(int istal, int istam, int istar, c
{
ca::FieldValue<fvec> B0 = fFld0Sta[0]->fieldSlice.GetFieldValueForLine(T);
ca::FieldValue<fvec> B1 = fFld0Sta[1]->fieldSlice.GetFieldValueForLine(T);
- fld0.Set(fAlgo->fTargB, fAlgo->fTargZ, B0, fFld0Sta[0]->fZ, B1, fFld0Sta[1]->fZ);
+ fld0.Set(fAlgo->fWData.TargB(), fAlgo->fWData.TargZ(), B0, fFld0Sta[0]->fZ, B1, fFld0Sta[1]->fZ);
}
{ // field, made by the left hit, the middle station and the right station
@@ -168,7 +170,7 @@ void TripletConstructor::CreateTripletsForHit(int istal, int istam, int istar, c
// add the target constraint
- fFit.FilterWithTargetAtLine(fAlgo->fTargZ, fAlgo->TargetMeasurement, fld0);
+ fFit.FilterWithTargetAtLine(fAlgo->fWData.TargZ(), fAlgo->fWData.TargetMeasurement(), fld0);
fFit.MultipleScattering(fAlgo->fParameters.GetMaterialThickness(fIstaL, T.GetX(), T.GetY()));
@@ -197,7 +199,8 @@ void TripletConstructor::FindDoublets()
}
}
- CollectHits(fTrL, fIstaM, fAlgo->fDoubletChi2Cut, iMC, fHitsM_2, fAlgo->fParameters.GetMaxDoubletsPerSinglet());
+ CollectHits(fTrL, fIstaM, fAlgo->fWData.CurrentIteration()->GetDoubletChi2Cut(), iMC, fHitsM_2,
+ fAlgo->fParameters.GetMaxDoubletsPerSinglet());
}
@@ -214,6 +217,7 @@ void TripletConstructor::FitDoublets()
bool isMomentumFitted = (fAlgo->fIsTargetField || (staL().fieldStatus != 0) || (staM().fieldStatus != 0));
+ double maxQp = fAlgo->fWData.CurrentIteration()->GetMaxQp();
for (unsigned int i2 = 0; i2 < hitsMtmp.size(); i2++) {
const ca::HitIndex_t indM = hitsMtmp[i2];
@@ -240,7 +244,7 @@ void TripletConstructor::FitDoublets()
fFit.FilterTime(t_2, dt2_2, staM().timeInfo);
- fFit.SetQp0(isMomentumFitted ? fFit.Tr().GetQp() : fAlgo->fMaxInvMom);
+ fFit.SetQp0(isMomentumFitted ? fFit.Tr().GetQp() : maxQp);
fFit.MultipleScattering(fAlgo->fParameters.GetMaterialThickness(fIstaM, T2.GetX(), T2.Y()));
@@ -324,6 +328,8 @@ void TripletConstructor::FindRightHit()
}
// ---- Add the middle hits to parameters estimation. Propagate to right station. ----
+ double maxSlope = fAlgo->fWData.CurrentIteration()->GetMaxSlope();
+ double tripletChi2Cut = fAlgo->fWData.CurrentIteration()->GetTripletChi2Cut();
for (unsigned int i2 = 0; i2 < fHitsM_2.size(); i2++) {
fFit.SetTrack(fTracks_2[i2]);
@@ -360,10 +366,10 @@ void TripletConstructor::FindRightHit()
if (T2.C00()[0] < 0 || T2.C11()[0] < 0 || T2.C22()[0] < 0 || T2.C33()[0] < 0 || T2.C55()[0] < 0) {
break;
}
- if (fabs(T2.Tx()[0]) > fAlgo->fMaxSlope) {
+ if (fabs(T2.Tx()[0]) > maxSlope) {
break;
}
- if (fabs(T2.Ty()[0]) > fAlgo->fMaxSlope) {
+ if (fabs(T2.Ty()[0]) > maxSlope) {
break;
}
if (fAlgo->fParameters.DevIsMatchTripletsViaMc()) {
@@ -391,7 +397,7 @@ void TripletConstructor::FindRightHit()
}
Vector<ca::HitIndex_t> collectedHits;
- CollectHits(T2, fIstaR, fAlgo->fTripletChi2Cut, iMC, collectedHits, fAlgo->fParameters.GetMaxTripletPerDoublets());
+ CollectHits(T2, fIstaR, tripletChi2Cut, iMC, collectedHits, fAlgo->fParameters.GetMaxTripletPerDoublets());
if (collectedHits.size() >= fAlgo->fParameters.GetMaxTripletPerDoublets()) {
LOG(debug) << "Ca: GetMaxTripletPerDoublets==" << fAlgo->fParameters.GetMaxTripletPerDoublets()
@@ -445,7 +451,7 @@ void TripletConstructor::FitTriplets()
ca::TrackFit fit;
fit.SetMask(fmask::One());
- if (fAlgo->fpCurrentIteration->GetElectronFlag()) {
+ if (fAlgo->fWData.CurrentIteration()->GetElectronFlag()) {
fit.SetParticleMass(constants::phys::ElectronMass);
}
else {
@@ -466,6 +472,7 @@ void TripletConstructor::FitTriplets()
fvec ndfTrackModel = 4; // straight line
ndfTrackModel += isMomentumFitted ? 1 : 0; // track with momentum
+ const fvec maxQp = fAlgo->fWData.CurrentIteration()->GetMaxQp();
for (int i3 = 0; i3 < n3; ++i3) {
// prepare data
@@ -512,7 +519,7 @@ void TripletConstructor::FitTriplets()
};
fld.Set(B[0], sta[0].fZ, B[1], sta[1].fZ, B[2], sta[2].fZ);
- fldTarget.Set(fAlgo->fTargB, fAlgo->fTargZ, B[0], sta[0].fZ, B[1], sta[1].fZ);
+ fldTarget.Set(fAlgo->fWData.TargB(), fAlgo->fWData.TargZ(), B[0], sta[0].fZ, B[1], sta[1].fZ);
TrackParamV& T = fit.Tr();
@@ -530,8 +537,8 @@ void TripletConstructor::FitTriplets()
// fit forward
{
fit.SetQp0(T.Qp());
- fit.Qp0()(fit.Qp0() > fAlgo->GetMaxInvMom()) = fAlgo->GetMaxInvMom();
- fit.Qp0()(fit.Qp0() < -fAlgo->GetMaxInvMom()) = -fAlgo->GetMaxInvMom();
+ fit.Qp0()(fit.Qp0() > maxQp) = maxQp;
+ fit.Qp0()(fit.Qp0() < -maxQp) = -maxQp;
T.Qp() = 0.;
T.Vi() = 0.;
@@ -550,7 +557,7 @@ void TripletConstructor::FitTriplets()
T.NdfTime() = sta[ih0].timeInfo ? 0 : -1;
// add the target constraint
- fit.FilterWithTargetAtLine(fAlgo->fTargZ, fAlgo->TargetMeasurement, fldTarget);
+ fit.FilterWithTargetAtLine(fAlgo->fWData.TargZ(), fAlgo->fWData.TargetMeasurement(), fldTarget);
for (int ih = 1; ih < NHits; ++ih) {
fit.Extrapolate(z[ih], fld);
@@ -567,8 +574,8 @@ void TripletConstructor::FitTriplets()
// fit backward
{
fit.SetQp0(T.Qp());
- fit.Qp0()(fit.Qp0() > fAlgo->GetMaxInvMom()) = fAlgo->GetMaxInvMom();
- fit.Qp0()(fit.Qp0() < -fAlgo->GetMaxInvMom()) = -fAlgo->GetMaxInvMom();
+ fit.Qp0()(fit.Qp0() > maxQp) = maxQp;
+ fit.Qp0()(fit.Qp0() < -maxQp) = -maxQp;
int ih0 = NHits - 1;
T.X() = x[ih0];
@@ -642,6 +649,7 @@ void TripletConstructor::StoreTriplets()
fTriplets.clear();
fTriplets.reserve(n3);
+ const float tripletFinalChi2Cut = fAlgo->fWData.CurrentIteration()->GetTripletFinalChi2Cut();
for (Tindex i3 = 0; i3 < n3; ++i3) {
TrackParamV& T3 = fTracks_3[i3];
@@ -658,8 +666,8 @@ void TripletConstructor::StoreTriplets()
CBMCA_DEBUG_ASSERT(ihitm < fAlgo->fWData.HitStartIndexOnStation(fIstaM + 1));
CBMCA_DEBUG_ASSERT(ihitr < fAlgo->fWData.HitStartIndexOnStation(fIstaR + 1));
- if (!fAlgo->fpCurrentIteration->GetTrackFromTripletsFlag()) {
- if (chi2 > fAlgo->fTripletFinalChi2Cut) {
+ if (!fAlgo->fWData.CurrentIteration()->GetTrackFromTripletsFlag()) {
+ if (chi2 > tripletFinalChi2Cut) {
continue;
}
}
@@ -706,13 +714,13 @@ void TripletConstructor::CollectHits(const TrackParamV& Tr, const int iSta, cons
const fscal time = T.Time()[0];
const auto& grid = fAlgo->fWData.Grid(iSta);
- ca::GridArea area(grid, T.X()[0], T.Y()[0],
- (sqrt(Pick_m22 * T.C00()) + grid.GetMaxRangeX() + fAlgo->fMaxDZ * abs(T.Tx()))[0],
- (sqrt(Pick_m22 * T.C11()) + grid.GetMaxRangeY() + fAlgo->fMaxDZ * abs(T.Ty()))[0]);
+ const fvec maxDZ = fAlgo->fWData.CurrentIteration()->GetMaxDZ();
+ ca::GridArea area(grid, T.X()[0], T.Y()[0], (sqrt(Pick_m22 * T.C00()) + grid.GetMaxRangeX() + maxDZ * abs(T.Tx()))[0],
+ (sqrt(Pick_m22 * T.C11()) + grid.GetMaxRangeY() + maxDZ * abs(T.Ty()))[0]);
if constexpr (fDebugCollectHits) {
LOG(info) << "search area: " << T.X()[0] << " " << T.Y()[0] << " "
- << (sqrt(Pick_m22 * T.C00()) + grid.GetMaxRangeX() + fAlgo->fMaxDZ * abs(T.Tx()))[0] << " "
- << (sqrt(Pick_m22 * T.C11()) + grid.GetMaxRangeY() + fAlgo->fMaxDZ * abs(T.Ty()))[0];
+ << (sqrt(Pick_m22 * T.C00()) + grid.GetMaxRangeX() + maxDZ * abs(T.Tx()))[0] << " "
+ << (sqrt(Pick_m22 * T.C11()) + grid.GetMaxRangeY() + maxDZ * abs(T.Ty()))[0];
}
if (fAlgo->fParameters.DevIsIgnoreHitSearchAreas()) {
area.DoLoopOverEntireGrid();
@@ -797,7 +805,7 @@ void TripletConstructor::CollectHits(const TrackParamV& Tr, const int iSta, cons
auto [chi2x, chi2u] = ca::TrackFit::GetChi2XChi2U(mxy, x, y, C00, C10, C11);
// TODO: adjust the cut, cut on chi2x & chi2u separately
- if (!fAlgo->fpCurrentIteration->GetTrackFromTripletsFlag()) {
+ if (!fAlgo->fWData.CurrentIteration()->GetTrackFromTripletsFlag()) {
if (chi2x[0] > chi2Cut) {
if constexpr (fDebugCollectHits) {
LOG(info) << " hit chi2X is too large";
diff --git a/reco/L1/L1Algo/utils/L1AlgoEfficiencyPerformance.h b/reco/L1/L1Algo/utils/L1AlgoEfficiencyPerformance.h
index 36f97483603db5660055fba1a64f51601d9f7123..7023d4dc0e69de43984b09a28dbf09de0a621b57 100644
--- a/reco/L1/L1Algo/utils/L1AlgoEfficiencyPerformance.h
+++ b/reco/L1/L1Algo/utils/L1AlgoEfficiencyPerformance.h
@@ -202,7 +202,14 @@ void L1AlgoEfficiencyPerformance<NHits>::FillMC()
trlet.mcTrackId = mcTrk.GetId();
trlet.mcMotherTrackId = mcTrk.GetMotherId();
trlet.p = mcTrk.GetP();
- if (mcTrk.GetP() > 1. / fL1->fpAlgo->GetMaxInvMom()) trlet.SetAsReconstructable();
+ if (fL1->fpAlgo->fWData.CurrentIteration()) {
+ if (mcTrk.GetP() > 1. / fL1->fpAlgo->fWData.CurrentIteration()->GetMaxQp()) {
+ trlet.SetAsReconstructable();
+ }
+ }
+ else {
+ LOG(error) << "L1AlgoEfficiencyPerformance::FillMC(): attempt to call fpCurrentIteration, but it is a nullptr";
+ }
mcTracklets.push_back(trlet);
} // for Iter = stations