diff --git a/algo/ca/core/data/CaWindowData.h b/algo/ca/core/data/CaWindowData.h
index a9c8e4366712bea50ca6b52553b36d4a103a146c..b599ae6a03da1658eed297236b0208e01f3c5712 100644
--- a/algo/ca/core/data/CaWindowData.h
+++ b/algo/ca/core/data/CaWindowData.h
@@ -149,12 +149,6 @@ namespace cbm::algo::ca
/// \brief Accesses current iteration
[[gnu::always_inline]] const Iteration* CurrentIteration() const { return fpCurrentIteration; }
- /// \brief Accesses current iteration index (DEBUG)
- [[gnu::always_inline]] int CurrentIterationIndex() const { return fCurrentIterationIndex; }
-
- /// \brief Accesses current iteration index (DEBUG)
- [[gnu::always_inline]] void SetCurrentIterationIndex(int id) { fCurrentIterationIndex = id; }
-
/// \brief Accesses magnetic field in starting point (target or first station)
[[gnu::always_inline]] ca::FieldValue<fvec>& TargB() { return fTargB; }
@@ -203,7 +197,6 @@ namespace cbm::algo::ca
/// \brief Current track-finder iteration
const Iteration* fpCurrentIteration = nullptr;
- int fCurrentIterationIndex = -1; // TODO: Get rid of this field, it can cause bugs
ca::FieldValue<fvec> fTargB _fvecalignment{}; ///< field in the target point (modifiable, do not touch!!)
ca::MeasurementXy<fvec> fTargetMeasurement _fvecalignment{}; ///< target constraint
diff --git a/algo/ca/core/tracking/CaTrackFinder.cxx b/algo/ca/core/tracking/CaTrackFinder.cxx
index 6afe5987e10de282076714fff6eb41641c8afe0b..6c19e6c26ff65a208cc6c540cbd0cef68a53cafa 100644
--- a/algo/ca/core/tracking/CaTrackFinder.cxx
+++ b/algo/ca/core/tracking/CaTrackFinder.cxx
@@ -189,9 +189,8 @@ TrackFinder::Output_t TrackFinder::FindTracks(const InputData& input, TimesliceH
}
}
}
-
- fStatTsEnd = std::max(fStatTsEnd, fStatTsStart); // all hits belong to one sub-timeslice
- fStatTsEnd = std::min(fStatTsEnd, fStatTsStart + 500.e6f); // 500 ms maximal length of the TS
+ // all hits belong to one sub-timeslice; 500 ms maximal length of the TS
+ fStatTsEnd = std::clamp(fStatTsEnd, fStatTsStart, fStatTsStart + 500.e6f);
LOG(debug) << "CA tracker process time slice " << fStatTsStart * 1.e-6 << " -- " << fStatTsEnd * 1.e-6
<< " [ms] with " << fStatNhitsTotal << " hits";
@@ -213,7 +212,7 @@ TrackFinder::Output_t TrackFinder::FindTracks(const InputData& input, TimesliceH
{ // Estimation of number of hits in time windows
//Timer time;
//time.Start();
- int nSt = fParameters.GetNstationsActive();
+ const int nSt = fParameters.GetNstationsActive();
// Count number of hits per window and station
const double a = fStatTsStart + fWindowOverlap + fWindowMargin;
diff --git a/algo/ca/core/tracking/CaTrackFinderWindow.cxx b/algo/ca/core/tracking/CaTrackFinderWindow.cxx
index ff9b0a5b8ed3d8996c7d3e0b5237a357679eef33..4fdb1f2ab8469d09fa52f58a01584aaf33c55dfc 100644
--- a/algo/ca/core/tracking/CaTrackFinderWindow.cxx
+++ b/algo/ca/core/tracking/CaTrackFinderWindow.cxx
@@ -81,7 +81,6 @@ bool TrackFinderWindow::checkTripletMatch(const ca::Triplet& l, const ca::Triple
if (r.GetMSta() != l.GetRSta()) return false;
if (r.GetLSta() != l.GetMSta()) return false;
-
const float tripletLinkChi2 = wData.CurrentIteration()->GetTripletLinkChi2();
if (r.IsMomentumFitted()) {
assert(l.IsMomentumFitted());
@@ -132,7 +131,7 @@ bool TrackFinderWindow::checkTripletMatch(const ca::Triplet& l, const ca::Triple
// ---------------------------------------------------------------------------------------------------------------------
//
-void TrackFinderWindow::ReadWindowData(const ca::InputData& input, WindowData& wData)
+void TrackFinderWindow::ReadWindowData(const Vector<Hit>& hits, WindowData& wData)
{
int nHits = 0;
for (int iS = 0; iS < fParameters.GetNstationsActive(); iS++) {
@@ -146,7 +145,7 @@ void TrackFinderWindow::ReadWindowData(const ca::InputData& input, WindowData& w
for (int iS = 0; iS < fParameters.GetNstationsActive(); iS++) {
int iFstHit = wData.HitStartIndexOnStation(iS);
for (ca::HitIndex_t ih = 0; ih < wData.TsHitIndices(iS).size(); ++ih) {
- ca::Hit h = input.GetHit(wData.TsHitIndex(iS, ih));
+ ca::Hit h = hits[wData.TsHitIndex(iS, ih)]; /// D.S. 29.7.24: There is a copy operation here. Can be avoided?
h.SetId(wData.TsHitIndex(iS, ih));
wData.Hit(iFstHit + ih) = h;
}
@@ -160,20 +159,11 @@ void TrackFinderWindow::ReadWindowData(const ca::InputData& input, WindowData& w
LOG(info) << "===== ";
}
-
wData.RecoTracks().clear();
wData.RecoTracks().reserve(2 * nHits / fParameters.GetNstationsActive());
wData.RecoHitIndices().clear();
wData.RecoHitIndices().reserve(2 * nHits);
-
- fvTrackCandidates.clear();
- fvTrackCandidates.reserve(nHits / 10);
- for (unsigned int iS = 0; iS < constants::size::MaxNstations; iS++) {
- int nHitsStation = wData.TsHitIndices(iS).size();
- fvTriplets[iS].clear();
- fvTriplets[iS].reserve(2 * nHitsStation);
- }
}
void TrackFinderWindow::CaTrackFinderSlice(const ca::InputData& input, WindowData& wData)
@@ -183,66 +173,51 @@ void TrackFinderWindow::CaTrackFinderSlice(const ca::InputData& input, WindowDat
***********************************/
frMonitorData.StartTimer(ETimer::InitWindow);
- ReadWindowData(input, wData);
+ ReadWindowData(input.GetHits(), wData);
frMonitorData.StopTimer(ETimer::InitWindow);
frMonitorData.StartTimer(ETimer::PrepareGrid);
for (int iS = 0; iS < fParameters.GetNstationsActive(); ++iS) {
- bool timeUninitialised = 1;
- fscal lasttime = 0;
- fscal starttime = 0;
-
-
- fscal gridMinX = -0.1;
- fscal gridMaxX = 0.1;
- fscal gridMinY = -0.1;
- fscal gridMaxY = 0.1;
+ fscal lasttime = std::numeric_limits<fscal>::infinity();
+ fscal starttime = -std::numeric_limits<fscal>::infinity();
+ fscal gridMinX = -0.1;
+ fscal gridMaxX = 0.1;
+ fscal gridMinY = -0.1;
+ fscal gridMaxY = 0.1;
for (ca::HitIndex_t ih = 0; ih < wData.TsHitIndices(iS).size(); ++ih) {
const ca::Hit& h = input.GetHit(wData.TsHitIndex(iS, ih));
- if (h.X() < gridMinX) {
- gridMinX = h.X();
- }
- if (h.X() > gridMaxX) {
- gridMaxX = h.X();
- }
- if (h.Y() < gridMinY) {
- gridMinY = h.Y();
- }
- if (h.Y() > gridMaxY) {
- gridMaxY = h.Y();
- }
+ gridMinX = std::min(gridMinX, h.X());
+ gridMinY = std::min(gridMinY, h.Y());
+
+ gridMaxX = std::max(gridMaxX, h.X());
+ gridMaxY = std::max(gridMaxY, h.Y());
const fscal time = h.T();
assert(std::isfinite(time));
- if (timeUninitialised || lasttime < time) {
- lasttime = time;
- }
- if (timeUninitialised || starttime > time) {
- starttime = time;
- }
- timeUninitialised = false;
+ lasttime = std::min(lasttime, time);
+ starttime = std::max(starttime, time);
}
// TODO: changing the grid also changes the result. Investigate why it happens.
// TODO: find the optimal grid size
- int nSliceHits = wData.TsHitIndices(iS).size();
- fscal sizeY = gridMaxY - gridMinY;
- fscal sizeX = gridMaxX - gridMinX;
- int nBins2D = 1 + nSliceHits;
+ const int nSliceHits = wData.TsHitIndices(iS).size();
+ const fscal sizeY = gridMaxY - gridMinY;
+ const fscal sizeX = gridMaxX - gridMinX;
+ const int nBins2D = 1 + nSliceHits;
+
// TODO: SG: the coefficients should be removed
+ const fscal scale = fParameters.GetStation(iS).GetZ<float>() - fParameters.GetTargetPositionZ()[0];
+ const fscal maxScale = 0.3 * scale;
+ const fscal minScale = 0.01 * scale;
+
fscal yStep = 0.3 * sizeY / sqrt(nBins2D);
fscal xStep = 0.8 * sizeX / sqrt(nBins2D);
-
- fscal scale = fParameters.GetStation(iS).GetZ<float>() - fParameters.GetTargetPositionZ()[0];
-
- if (yStep > 0.3 * scale) yStep = 0.3 * scale;
- if (yStep < 0.01 * scale) yStep = 0.01 * scale;
- if (xStep > 0.3 * scale) xStep = 0.3 * scale;
- if (xStep < 0.01 * scale) xStep = 0.01 * scale;
+ yStep = std::clamp(yStep, minScale, maxScale);
+ xStep = std::clamp(xStep, minScale, maxScale);
auto& grid = wData.Grid(iS);
grid.BuildBins(gridMinX, gridMaxX, gridMinY, gridMaxY, xStep, yStep);
@@ -258,63 +233,61 @@ void TrackFinderWindow::CaTrackFinderSlice(const ca::InputData& input, WindowDat
// ----- Find tracks -----
//
frMonitorData.StartTimer(ETimer::FindTracks);
- // TODO: Get rid of fCurrentTierationIndex. If it is needed of logs, please used caIteration.GetName()
- for (int iIteration = 0; iIteration < fParameters.GetNcaIterations(); ++iIteration) {
+ auto& caIterations = fParameters.GetCAIterations();
+
+ for (auto iter = caIterations.begin(); iter != caIterations.end(); ++iter) {
// ----- Prepare iteration
//
frMonitorData.StartTimer(ETimer::PrepareIteration);
- const auto& caIteration = fParameters.GetCAIterations()[iIteration];
+ const auto& caIteration = *iter; // Dereferencing iterator to get the value
wData.SetCurrentIteration(&caIteration);
- wData.SetCurrentIterationIndex(iIteration);
- if (iIteration > 0) {
+ // Use iter to check if it is not the first element
+ if (iter != caIterations.begin()) {
for (int ista = 0; ista < fParameters.GetNstationsActive(); ++ista) {
wData.Grid(ista).RemoveUsedHits(wData.Hits(), wData.HitKeyFlags());
}
}
- int nHits = wData.Hits().size();
// --> frAlgo.fIsWindowHitSuppressed.reset(frAlgo.fWindowHits.size(), 0);
wData.ResetHitSuppressionFlags(); // TODO: ??? No effect?
+ //
for (int j = 0; j < fParameters.GetNstationsActive(); j++) {
+ const size_t nHitsStation = wData.TsHitIndices(j).size();
fvTriplets[j].clear();
+ fvTriplets[j].reserve(2 * nHitsStation);
}
- Vector<int> vHitFirstTriplet(nHits, 0); /// link hit -> first triplet { hit, *, *}
- Vector<int> vHitNofTriplets(nHits, 0); /// link hit ->n triplets { hit, *, *}
-
+ // #pragma omp task
{
- // #pragma omp task
- {
- // --- SET PARAMETERS FOR THE ITERATION ---
- // define the target
- 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()) {
- wData.TargB() = fParameters.GetVertexFieldValue();
- }
- else {
- wData.TargB() = fParameters.GetStation(0).fieldSlice.GetFieldValue(0, 0);
- } // NOTE: calculates field frAlgo.fTargB in the center of 0th station
-
- wData.TargetMeasurement().SetX(fParameters.GetTargetPositionX());
- wData.TargetMeasurement().SetY(fParameters.GetTargetPositionY());
- wData.TargetMeasurement().SetDx2(SigmaTargetX * SigmaTargetX);
- wData.TargetMeasurement().SetDxy(0);
- wData.TargetMeasurement().SetDy2(SigmaTargetY * SigmaTargetY);
- wData.TargetMeasurement().SetNdfX(1);
- wData.TargetMeasurement().SetNdfY(1);
+ // --- SET PARAMETERS FOR THE ITERATION ---
+ // define the target
+ const fscal SigmaTargetX = caIteration.GetTargetPosSigmaX();
+ const fscal SigmaTargetY = caIteration.GetTargetPosSigmaY(); // target constraint [cm]
+
+ // Select magnetic field. For primary tracks - fVtxFieldValue, for secondary tracks - st.fieldSlice
+ if (caIteration.GetPrimaryFlag()) {
+ wData.TargB() = fParameters.GetVertexFieldValue();
}
+ else {
+ wData.TargB() = fParameters.GetStation(0).fieldSlice.GetFieldValue(0, 0);
+ } // NOTE: calculates field frAlgo.fTargB in the center of 0th station
+
+ wData.TargetMeasurement().SetX(fParameters.GetTargetPositionX());
+ wData.TargetMeasurement().SetY(fParameters.GetTargetPositionY());
+ wData.TargetMeasurement().SetDx2(SigmaTargetX * SigmaTargetX);
+ wData.TargetMeasurement().SetDxy(0);
+ wData.TargetMeasurement().SetDy2(SigmaTargetY * SigmaTargetY);
+ wData.TargetMeasurement().SetNdfX(1);
+ wData.TargetMeasurement().SetNdfY(1);
}
+
frMonitorData.StopTimer(ETimer::PrepareIteration);
// ----- Triplets construction -----
//
frMonitorData.StartTimer(ETimer::ConstructTriplets);
-
ca::TripletConstructor constructor(fParameters, input, wData, fDefaultMass, fTrackingMode);
// indices of the two neighbouring station, taking into account allowed gaps
@@ -325,14 +298,16 @@ void TrackFinderWindow::CaTrackFinderSlice(const ca::InputData& input, WindowDat
}
}
+ Vector<int> vHitFirstTriplet(wData.Hits().size(), 0); /// link hit -> first triplet { hit, *, *}
+ Vector<int> vHitNofTriplets(wData.Hits().size(), 0); /// link hit ->n triplets { hit, *, *}
+
const int iStFirst = wData.CurrentIteration()->GetFirstStationIndex();
for (int istal = fParameters.GetNstationsActive() - 2; istal >= iStFirst; istal--) {
// start with downstream chambers
- const auto& grid = wData.Grid(istal);
- Tindex nGridEntriesL = grid.GetEntries().size();
- for (Tindex iel = 0; iel < nGridEntriesL; ++iel) {
- ca::HitIndex_t ihitl = grid.GetEntries()[iel].GetObjectId();
- auto oldSize = fvTriplets[istal].size();
+ const auto& grid = wData.Grid(istal);
+ for (auto& entry : grid.GetEntries()) {
+ ca::HitIndex_t ihitl = entry.GetObjectId();
+ const size_t oldSize = fvTriplets[istal].size();
for (auto& pattern : staPattern) {
constructor.CreateTripletsForHit(istal, istal + pattern.first, istal + pattern.second, ihitl);
fvTriplets[istal].insert(fvTriplets[istal].end(), constructor.GetTriplets().begin(),
@@ -350,13 +325,8 @@ void TrackFinderWindow::CaTrackFinderSlice(const ca::InputData& input, WindowDat
for (int istal = fParameters.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];
- tr.SetLevel(0);
- tr.SetFNeighbour(0);
- tr.SetNNeighbours(0);
+ for (ca::Triplet& tr : fvTriplets[istal]) {
unsigned int nNeighbours = vHitNofTriplets[tr.GetMHit()];
unsigned int neighLocation = vHitFirstTriplet[tr.GetMHit()];
unsigned int neighStation = TripletId2Station(neighLocation);
@@ -379,14 +349,9 @@ void TrackFinderWindow::CaTrackFinderSlice(const ca::InputData& input, WindowDat
if (tr.GetFNeighbour() == 0) {
tr.SetFNeighbour(neighLocation);
}
-
tr.SetNNeighbours(neighLocation - tr.GetFNeighbour() + 1);
- const unsigned char neighbourLevel = neighbour.GetLevel();
-
- if (level <= neighbourLevel) {
- level = neighbourLevel + 1;
- }
+ level = std::max(level, static_cast<unsigned char>(neighbour.GetLevel() + 1));
}
tr.SetLevel(level);
} // neighbour search
@@ -401,19 +366,13 @@ void TrackFinderWindow::CaTrackFinderSlice(const ca::InputData& input, WindowDat
frMonitorData.StartTimer(ETimer::CreateTracks);
// min level to start triplet. So min track length = min_level+3.
- int min_level = std::min(caIteration.GetMinNhits(), caIteration.GetMinNhitsStation0()) - 3;
+ const int min_level = wData.CurrentIteration()->GetTrackFromTripletsFlag()
+ ? 0
+ : std::min(caIteration.GetMinNhits(), caIteration.GetMinNhitsStation0()) - 3;
- if (wData.CurrentIteration()->GetTrackFromTripletsFlag()) {
- min_level = 0;
- }
-
- ca::Branch curr_tr;
ca::Branch new_tr[constants::size::MaxNstations];
ca::Branch best_tr;
- int ndf = 1;
-
-
// collect consequtive: the longest tracks, shorter, more shorter and so on
for (int firstTripletLevel = fParameters.GetNstationsActive() - 3; firstTripletLevel >= min_level;
firstTripletLevel--) {
@@ -425,7 +384,6 @@ void TrackFinderWindow::CaTrackFinderSlice(const ca::InputData& input, WindowDat
const unsigned char min_best_l =
(firstTripletLevel > min_level) ? firstTripletLevel + 2 : min_level + 3; // loose maximum
-
fvTrackCandidates.clear();
fvTrackCandidates.reserve(wData.Hits().size() / 10);
@@ -440,9 +398,8 @@ void TrackFinderWindow::CaTrackFinderSlice(const ca::InputData& input, WindowDat
if (--nlevel == 0) break; //TODO: SG: this is not needed
- for (Tindex itrip = 0; itrip < (Tindex) fvTriplets[istaF].size(); ++itrip) {
+ for (ca::Triplet& first_trip : fvTriplets[istaF]) {
- ca::Triplet& first_trip = (fvTriplets[istaF][itrip]);
const auto& fstTripLHit = wData.Hit(first_trip.GetLHit());
if (wData.IsHitKeyUsed(fstTripLHit.FrontKey()) || wData.IsHitKeyUsed(fstTripLHit.BackKey())) {
continue;
@@ -469,9 +426,7 @@ void TrackFinderWindow::CaTrackFinderSlice(const ca::InputData& input, WindowDat
continue;
}
- curr_tr.SetChi2(0.);
- curr_tr.SetStation(0);
- curr_tr.ResetHits();
+ ca::Branch curr_tr;
curr_tr.AddHit(fstTripLHit.Id());
curr_tr.SetChi2(first_trip.GetChi2());
@@ -488,7 +443,7 @@ void TrackFinderWindow::CaTrackFinderSlice(const ca::InputData& input, WindowDat
continue;
}
- ndf = best_tr.NofHits() * 2 - 5;
+ int ndf = best_tr.NofHits() * 2 - 5;
// TODO: automatize the NDF calculation
@@ -511,7 +466,7 @@ void TrackFinderWindow::CaTrackFinderSlice(const ca::InputData& input, WindowDat
tr.SetAlive(true);
if constexpr (fDebug) {
std::stringstream s;
- s << "iter " << wData.CurrentIterationIndex() << ", track candidate " << fvTrackCandidates.size() - 1
+ s << "iter " << wData.CurrentIteration()->GetName() << ", track candidate " << fvTrackCandidates.size() - 1
<< " found, L = " << best_tr.NofHits() << " chi2= " << best_tr.Chi2() << " hits: ";
for (auto hitId : tr.Hits()) {
s << hitId << " (mc " << ca::Framework::GetMcTrackIdForCaHit(hitId) << ") ";
@@ -529,54 +484,37 @@ void TrackFinderWindow::CaTrackFinderSlice(const ca::InputData& input, WindowDat
int iComp = 0;
for (iComp = 0; (iComp < 100) && repeatCompetition; ++iComp) {
- // == Loop over track candidates and mark their strips
-
repeatCompetition = false;
+ // == Loop over track candidates and mark their strips
for (ca::Branch& tr : fvTrackCandidates) {
-
if (tr.IsAlive()) {
continue;
}
for (auto& hitId : tr.Hits()) {
- const ca::Hit& h = input.GetHit(hitId);
- bool isAlive = true;
- { // front strip
- auto& stripF = fvHitKeyToTrack[h.FrontKey()];
- if ((stripF >= 0) && (stripF != tr.Id())) { // strip is used by other candidate
- const auto& other = fvTrackCandidates[stripF];
- if (!other.IsAlive() && tr.IsBetterThan(other)) {
- stripF = tr.Id();
- }
- else {
- isAlive = false;
- }
- }
- else {
- stripF = tr.Id();
- }
- if (!isAlive) {
- break;
- }
- }
- { // back strip
- auto& stripB = fvHitKeyToTrack[h.BackKey()];
- if ((stripB >= 0) && (stripB != tr.Id())) { // strip is used by other candidate
- const auto& other = fvTrackCandidates[stripB];
- if (!other.IsAlive() && tr.IsBetterThan(other)) {
- stripB = tr.Id();
+ auto updateStrip = [&](int& strip) {
+ if ((strip >= 0) && (strip != tr.Id())) { // strip is used by other candidate
+ const auto& other = fvTrackCandidates[strip];
+ if (!other.IsAlive() && tr.IsBetterThan(other)) { // D.S. 29.7.24: should this be || instead of && ?
+ strip = tr.Id();
}
else {
- isAlive = false;
+ return false;
}
}
else {
- stripB = tr.Id();
- }
- if (!isAlive) {
- break;
+ strip = tr.Id();
}
+ return true;
+ };
+
+ const ca::Hit& h = input.GetHit(hitId);
+ if (!updateStrip(fvHitKeyToTrack[h.FrontKey()])) { // front strip
+ break;
+ }
+ if (!updateStrip(fvHitKeyToTrack[h.BackKey()])) { // back strip
+ break;
}
} // loop over hits
} // itrack
@@ -589,10 +527,10 @@ void TrackFinderWindow::CaTrackFinderSlice(const ca::InputData& input, WindowDat
}
tr.SetAlive(true);
- for (int iHit = 0; tr.IsAlive() && (iHit < (int) tr.Hits().size()); ++iHit) {
- const ca::Hit& h = input.GetHit(tr.Hits()[iHit]);
- tr.SetAlive(tr.IsAlive() && (fvHitKeyToTrack[h.FrontKey()] == tr.Id())
- && (fvHitKeyToTrack[h.BackKey()] == tr.Id()));
+ for (const auto& hitIndex : tr.Hits()) {
+ if (!tr.IsAlive()) break;
+ const ca::Hit& h = input.GetHit(hitIndex);
+ tr.SetAlive((fvHitKeyToTrack[h.FrontKey()] == tr.Id()) && (fvHitKeyToTrack[h.BackKey()] == tr.Id()));
}
if (!tr.IsAlive()) { // release strips
@@ -621,10 +559,9 @@ void TrackFinderWindow::CaTrackFinderSlice(const ca::InputData& input, WindowDat
ca::Branch& tr = fvTrackCandidates[iCandidate];
if constexpr (fDebug) {
- LOG(info) << "iter " << wData.CurrentIterationIndex() << ", track candidate " << iCandidate
+ LOG(info) << "iter " << wData.CurrentIteration()->GetName() << ", track candidate " << iCandidate
<< ": alive = " << tr.IsAlive();
}
-
if (!tr.IsAlive()) continue;
if (wData.CurrentIteration()->GetExtendTracksFlag()) {
@@ -713,7 +650,6 @@ void TrackFinderWindow::CAFindTrack(int ista, ca::Branch& best_tr, const ca::Tri
if (!(wData.IsHitKeyUsed(hitM.FrontKey()) || wData.IsHitKeyUsed(hitM.BackKey()))) {
curr_tr.AddHit(hitM.Id());
}
-
if (!(wData.IsHitKeyUsed(hitR.FrontKey()) || wData.IsHitKeyUsed(hitR.BackKey()))) {
curr_tr.AddHit(hitR.Id());
}
@@ -725,13 +661,10 @@ void TrackFinderWindow::CAFindTrack(int ista, ca::Branch& best_tr, const ca::Tri
if (ca::TrackingMode::kGlobal == fTrackingMode || ca::TrackingMode::kMcbm == fTrackingMode) {
ndf = curr_tr.NofHits() * 2 - 4;
}
-
-
if (curr_tr.Chi2() > wData.CurrentIteration()->GetTrackChi2Cut() * ndf) {
return;
}
-
// -- select the best
if ((curr_tr.NofHits() > best_tr.NofHits())
|| ((curr_tr.NofHits() == best_tr.NofHits()) && (curr_tr.Chi2() < best_tr.Chi2()))) {
@@ -745,8 +678,7 @@ void TrackFinderWindow::CAFindTrack(int ista, ca::Branch& best_tr, const ca::Tri
for (Tindex in = 0; in < N_neighbour; in++) {
- unsigned int ID = curr_trip->GetFNeighbour() + in;
-
+ unsigned int ID = curr_trip->GetFNeighbour() + in;
unsigned int Station = TripletId2Station(ID);
unsigned int Triplet = TripletId2Triplet(ID);
@@ -768,7 +700,6 @@ void TrackFinderWindow::CAFindTrack(int ista, ca::Branch& best_tr, const ca::Tri
unsigned char new_L = curr_tr.NofHits() + 1;
fscal new_chi2 = curr_tr.Chi2() + dchi2;
-
if constexpr (0) { //SGtrd2d debug!!
int mc01 = ca::Framework::GetMcTrackIdForWindowHit(curr_trip->GetLHit());
int mc02 = ca::Framework::GetMcTrackIdForWindowHit(curr_trip->GetMHit());
diff --git a/algo/ca/core/tracking/CaTrackFinderWindow.h b/algo/ca/core/tracking/CaTrackFinderWindow.h
index 5e6aa450a0915e3d99989c575e1466b056a14c5f..0bb4b0e90577dfdf716ffc90ccddc5f36e62289b 100644
--- a/algo/ca/core/tracking/CaTrackFinderWindow.h
+++ b/algo/ca/core/tracking/CaTrackFinderWindow.h
@@ -64,7 +64,7 @@ namespace cbm::algo::ca
///-------------------------------
/// Private methods
- void ReadWindowData(const ca::InputData& input, WindowData& wData);
+ void ReadWindowData(const Vector<Hit>& hits, WindowData& wData);
void CAFindTrack(int ista, ca::Branch& best_tr, const ca::Triplet* curr_trip, ca::Branch& curr_tr,
unsigned char min_best_l, ca::Branch* new_tr, WindowData& wData);
diff --git a/algo/ca/core/tracking/CaTripletConstructor.cxx b/algo/ca/core/tracking/CaTripletConstructor.cxx
index e10dea41f48c182fb0ecb70c52240c0e1844d978..f07fe9653ba6d7c1c7b90c18d0ba2d62857672a5 100644
--- a/algo/ca/core/tracking/CaTripletConstructor.cxx
+++ b/algo/ca/core/tracking/CaTripletConstructor.cxx
@@ -37,38 +37,30 @@ TripletConstructor::TripletConstructor(const ca::Parameters<fvec>& pars, const c
- fParameters.GetStations().cbegin();
fIsTargetField = (fabs(frWData.TargB().y[0]) > 0.001);
+
+ fFit.SetMask(fmask::One());
+ fFit.SetDoFitVelocity(true);
}
-void TripletConstructor::InitStations(int istal, int istam, int istar)
+bool TripletConstructor::InitStations(int istal, int istam, int istar)
{
- fIsInitialized = false;
-
fIstaL = istal;
fIstaM = istam;
fIstaR = istar;
if (fIstaM >= fParameters.GetNstationsActive()) {
- return;
+ return false;
}
-
fStaL = &fParameters.GetStation(fIstaL);
fStaM = &fParameters.GetStation(fIstaM);
fStaR = &fParameters.GetStation(fIstaR);
-
{ // two stations for approximating the field between the target and the left hit
- int sta1 = fIstaL;
- if (sta1 >= fNfieldStations) {
- sta1 = fNfieldStations - 1;
- }
- if (sta1 < 1) {
- sta1 = 1;
- }
- int sta0 = sta1 / 2; // station between fIstaL and the target
+ const int sta1 = (fNfieldStations <= 1) ? 1 : std::clamp(fIstaL, 1, fNfieldStations - 1);
+ const int sta0 = sta1 / 2; // station between fIstaL and the target
assert(0 <= sta0 && sta0 < sta1 && sta1 < fParameters.GetNstationsActive());
-
fFld0Sta[0] = &fParameters.GetStation(sta0);
fFld0Sta[1] = &fParameters.GetStation(sta1);
}
@@ -81,12 +73,7 @@ void TripletConstructor::InitStations(int istal, int istam, int istar)
if (sta2 >= fParameters.GetNstationsActive()) {
sta2 = fIstaM;
sta1 = fIstaM - 1;
- if (sta0 >= sta1) {
- sta0 = sta1 - 1;
- }
- if (sta0 < 0) {
- sta0 = 0;
- }
+ sta0 = (sta1 <= 0) ? 0 : std::clamp(sta0, 0, sta1 - 1);
}
if (fParameters.GetNstationsActive() >= 3) {
assert(0 <= sta0 && sta0 < sta1 && sta1 < sta2 && sta2 < fParameters.GetNstationsActive());
@@ -96,23 +83,22 @@ void TripletConstructor::InitStations(int istal, int istam, int istar)
fFld1Sta[1] = &fParameters.GetStation(sta1);
fFld1Sta[2] = &fParameters.GetStation(sta2);
}
-
- fFit.SetParticleMass(fDefaultMass);
- if (frWData.CurrentIteration()->GetElectronFlag()) {
- fFit.SetParticleMass(constants::phys::ElectronMass);
- }
- fFit.SetMask(fmask::One());
- fFit.SetDoFitVelocity(true);
-
- fIsInitialized = true;
+ return true;
}
void TripletConstructor::CreateTripletsForHit(int istal, int istam, int istar, ca::HitIndex_t iHitL)
{
- InitStations(istal, istam, istar);
+ if (!InitStations(istal, istam, istar)) {
+ return;
+ }
- if (!fIsInitialized) return;
+ if (frWData.CurrentIteration()->GetElectronFlag()) {
+ fFit.SetParticleMass(constants::phys::ElectronMass);
+ }
+ else {
+ fFit.SetParticleMass(fDefaultMass);
+ }
fIhitL = iHitL;
const auto& hitl = frWData.Hit(fIhitL);
@@ -130,34 +116,24 @@ void TripletConstructor::CreateTripletsForHit(int istal, int istam, int istar, c
// get the magnetic field along the track.
// (suppose track is straight line with origin in the target)
-
- fvec xl = hitl.X();
- fvec yl = hitl.Y();
- fvec zl = hitl.Z();
- fvec time = hitl.T();
- fvec timeEr2 = hitl.dT2();
- const fvec dzli = 1.f / (zl - fParameters.GetTargetPositionZ());
-
- const fvec tx = (xl - fParameters.GetTargetPositionX()) * dzli;
- const fvec ty = (yl - fParameters.GetTargetPositionY()) * dzli;
-
- TrackParamV& T = fFit.Tr();
-
- T.X() = xl;
- T.Y() = yl;
- T.Z() = zl;
- T.Tx() = tx;
- T.Ty() = ty;
+ TrackParamV& T = fFit.Tr();
+ const fvec dzli = 1.f / (hitl.Z() - fParameters.GetTargetPositionZ());
+
+ T.X() = hitl.X();
+ T.Y() = hitl.Y();
+ T.Z() = hitl.Z();
+ T.Tx() = (hitl.X() - fParameters.GetTargetPositionX()) * dzli;
+ T.Ty() = (hitl.Y() - fParameters.GetTargetPositionY()) * dzli;
T.Qp() = fvec(0.);
- T.Time() = time;
+ T.Time() = hitl.T();
T.Vi() = constants::phys::SpeedOfLightInv;
const fvec maxSlopePV = frWData.CurrentIteration()->GetMaxSlopePV();
const fvec maxQp = frWData.CurrentIteration()->GetMaxQp();
- fvec txErr2 = maxSlopePV * maxSlopePV / fvec(9.);
- fvec qpErr2 = maxQp * maxQp / fvec(9.);
+ const fvec txErr2 = maxSlopePV * maxSlopePV / fvec(9.);
+ const fvec qpErr2 = maxQp * maxQp / fvec(9.);
- T.ResetErrors(1., 1., txErr2, txErr2, qpErr2, (staL().timeInfo ? timeEr2 : 1.e6), 1.e10);
+ T.ResetErrors(1., 1., txErr2, txErr2, qpErr2, (fStaL->timeInfo ? hitl.dT2() : 1.e6), 1.e10);
T.InitVelocityRange(1. / frWData.CurrentIteration()->GetMaxQp());
@@ -166,7 +142,7 @@ void TripletConstructor::CreateTripletsForHit(int istal, int istam, int istar, c
T.C11() = hitl.dY2();
T.Ndf() = (frWData.CurrentIteration()->GetPrimaryFlag()) ? fvec(2.) : fvec(0.);
- T.NdfTime() = (staL().timeInfo ? 0 : -1);
+ T.NdfTime() = (fStaL->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
@@ -195,7 +171,7 @@ void TripletConstructor::CreateTripletsForHit(int istal, int istam, int istar, c
// extrapolate to the middle hit
- fFit.ExtrapolateLine(staM().fZ, fFldL);
+ fFit.ExtrapolateLine(fStaM->fZ, fFldL);
fTrL = T;
}
@@ -233,7 +209,7 @@ void TripletConstructor::FitDoublets()
fTracks_2.clear();
fTracks_2.reserve(hitsMtmp.size());
- bool isMomentumFitted = (fIsTargetField || (staL().fieldStatus != 0) || (staM().fieldStatus != 0));
+ const bool isMomentumFitted = (fIsTargetField || (fStaL->fieldStatus != 0) || (fStaM->fieldStatus != 0));
double maxQp = frWData.CurrentIteration()->GetMaxQp();
for (unsigned int i2 = 0; i2 < hitsMtmp.size(); i2++) {
@@ -257,51 +233,46 @@ void TripletConstructor::FitDoublets()
// add the middle hit
fFit.ExtrapolateLineNoField(z_2);
-
fFit.FilterXY(m_2);
-
- fFit.FilterTime(t_2, dt2_2, staM().timeInfo);
-
+ fFit.FilterTime(t_2, dt2_2, fStaM->timeInfo);
fFit.SetQp0(isMomentumFitted ? fFit.Tr().GetQp() : maxQp);
-
fFit.MultipleScattering(fParameters.GetMaterialThickness(fIstaM, T2.GetX(), T2.Y()));
-
fFit.SetQp0(fFit.Tr().Qp());
// check if there are other hits close to the doublet on the same station
if (ca::kMcbm != fTrackingMode) {
// TODO: SG: adjust cuts, put them to parameters
- fscal tx = T2.Tx()[0];
- fscal ty = T2.Ty()[0];
- fscal tt = T2.Vi()[0] * sqrt(1. + tx * tx + ty * ty); // dt/dl * dl/dz
+ const fscal tx = T2.Tx()[0];
+ const fscal ty = T2.Ty()[0];
+ const fscal tt = T2.Vi()[0] * sqrt(1. + tx * tx + ty * ty); // dt/dl * dl/dz
for (unsigned int iClone = i2 + 1; iClone < hitsMtmp.size(); iClone++) {
- int indClone = hitsMtmp[iClone];
+ const int indClone = hitsMtmp[iClone];
const ca::Hit& hitClone = frWData.Hit(indClone);
- fscal dz = hitClone.Z() - T2.Z()[0];
+ const fscal dz = hitClone.Z() - T2.Z()[0];
- if ((staM().timeInfo) && (T2.NdfTime()[0] >= 0)) {
- fscal dt = T2.Time()[0] + tt * dz - hitClone.T();
+ if ((fStaM->timeInfo) && (T2.NdfTime()[0] >= 0)) {
+ const fscal dt = T2.Time()[0] + tt * dz - hitClone.T();
if (!(fabs(dt) <= 3.5 * sqrt(T2.C55()[0]) + hitClone.RangeT())) {
continue;
}
}
- fscal dx = T2.GetX()[0] + tx * dz - hitClone.X();
+ const fscal dx = T2.GetX()[0] + tx * dz - hitClone.X();
if (!(fabs(dx) <= 3.5 * sqrt(T2.C00()[0]) + hitClone.RangeX())) {
continue;
}
- fscal dy = T2.Y()[0] + ty * dz - hitClone.Y();
+ const fscal dy = T2.Y()[0] + ty * dz - hitClone.Y();
if (!(fabs(dy) <= 3.5 * sqrt(T2.C11()[0]) + hitClone.RangeY())) {
continue;
}
if (fParameters.DevIsSuppressOverlapHitsViaMc()) {
- int iMC = ca::Framework::GetMcTrackIdForWindowHit(fIhitL);
+ const int iMC = ca::Framework::GetMcTrackIdForWindowHit(fIhitL);
if ((iMC != ca::Framework::GetMcTrackIdForWindowHit(indM))
|| (iMC != ca::Framework::GetMcTrackIdForWindowHit(indClone))) {
continue;
@@ -348,8 +319,8 @@ void TripletConstructor::FindRightHit()
}
// ---- Add the middle hits to parameters estimation. Propagate to right station. ----
- double maxSlope = frWData.CurrentIteration()->GetMaxSlope();
- double tripletChi2Cut = frWData.CurrentIteration()->GetTripletChi2Cut();
+ const double maxSlope = frWData.CurrentIteration()->GetMaxSlope();
+ const double tripletChi2Cut = frWData.CurrentIteration()->GetTripletChi2Cut();
for (unsigned int i2 = 0; i2 < fHitsM_2.size(); i2++) {
fFit.SetTrack(fTracks_2[i2]);
@@ -357,14 +328,14 @@ void TripletConstructor::FindRightHit()
// extrapolate to the right hit station
- fFit.Extrapolate(staR().fZ, fFldL);
+ fFit.Extrapolate(fStaR->fZ, fFldL);
if constexpr (fDebugDublets) {
ca::HitIndex_t iwhit[2] = {fIhitL, fHitsM_2[i2]};
ca::HitIndex_t ihit[2] = {frWData.Hit(iwhit[0]).Id(), frWData.Hit(iwhit[1]).Id()};
- int ih0 = ihit[0];
- int ih1 = ihit[1];
+ const int ih0 = ihit[0];
+ const int ih1 = ihit[1];
const ca::Hit& h0 = fInputData.GetHit(ih0);
const ca::Hit& h1 = fInputData.GetHit(ih1);
@@ -377,30 +348,31 @@ void TripletConstructor::FindRightHit()
}
// ---- Find the triplets(right hit). Reformat data in the portion of triplets. ----
- bool isDoubletGood = false;
- int iMC = -1;
- do {
+ int iMC = -1;
+
+ auto rejectDoublet = [&]() -> bool {
if (ca::TrackingMode::kSts == fTrackingMode && (T2.C44()[0] < 0)) {
- break;
+ return true;
}
if (T2.C00()[0] < 0 || T2.C11()[0] < 0 || T2.C22()[0] < 0 || T2.C33()[0] < 0 || T2.C55()[0] < 0) {
- break;
+ return true;
}
if (fabs(T2.Tx()[0]) > maxSlope) {
- break;
+ return true;
}
if (fabs(T2.Ty()[0]) > maxSlope) {
- break;
+ return true;
}
if (fParameters.DevIsMatchTripletsViaMc()) {
int indM = fHitsM_2[i2];
iMC = ca::Framework::GetMcTrackIdForWindowHit(fIhitL);
if (iMC < 0 || iMC != ca::Framework::GetMcTrackIdForWindowHit(indM)) {
- break;
+ return true;
}
}
- isDoubletGood = true;
- } while (0);
+ return false;
+ };
+ const bool isDoubletGood = !rejectDoublet();
if constexpr (fDebugDublets) {
if (isDoubletGood) {
@@ -428,8 +400,7 @@ void TripletConstructor::FindRightHit()
if constexpr (fDebugDublets) {
LOG(info) << " collected " << collectedHits.size() << " hits on the right station ";
}
- for (unsigned int ih = 0; ih < collectedHits.size(); ih++) {
- ca::HitIndex_t irh = collectedHits[ih];
+ for (ca::HitIndex_t& irh : collectedHits) {
if constexpr (fDebugDublets) {
ca::HitIndex_t ihit = frWData.Hit(irh).Id();
const ca::Hit& h = fInputData.GetHit(ihit);
@@ -456,7 +427,7 @@ void TripletConstructor::FitTriplets()
{
constexpr int nIterations = 2;
- Tindex n3 = fHitsM_3.size();
+ const Tindex n3 = fHitsM_3.size();
assert(fHitsM_3.size() == fHitsR_3.size());
fTracks_3.clear();
@@ -477,28 +448,26 @@ void TripletConstructor::FitTriplets()
else {
fit.SetParticleMass(fDefaultMass);
}
- const int NHits = 3; // triplets
// prepare data
- int ista[NHits] = {fIstaL, fIstaM, fIstaR};
-
- ca::Station<fvec> sta[3];
- for (int is = 0; is < NHits; ++is) {
- sta[is] = fParameters.GetStation(ista[is]);
- };
+ const int NHits = 3; // triplets
+ const int ista[NHits] = {fIstaL, fIstaM, fIstaR};
- bool isMomentumFitted = ((staL().fieldStatus != 0) || (staM().fieldStatus != 0) || (staR().fieldStatus != 0));
+ const ca::Station<fvec> sta[NHits] = {fParameters.GetStation(ista[0]), fParameters.GetStation(ista[1]),
+ fParameters.GetStation(ista[2])};
- fvec ndfTrackModel = 4; // straight line
- ndfTrackModel += isMomentumFitted ? 1 : 0; // track with momentum
+ const bool isMomentumFitted = ((fStaL->fieldStatus != 0) || (fStaM->fieldStatus != 0) || (fStaR->fieldStatus != 0));
+ const fvec ndfTrackModel = 4 + (isMomentumFitted ? 1 : 0); // straight line or track with momentum
const fvec maxQp = frWData.CurrentIteration()->GetMaxQp();
+
for (int i3 = 0; i3 < n3; ++i3) {
// prepare data
- ca::HitIndex_t iwhit[NHits] = {fIhitL, fHitsM_3[i3], fHitsR_3[i3]};
+ const ca::HitIndex_t iwhit[NHits] = {fIhitL, fHitsM_3[i3], fHitsR_3[i3]};
- ca::HitIndex_t ihit[NHits] = {frWData.Hit(iwhit[0]).Id(), frWData.Hit(iwhit[1]).Id(), frWData.Hit(iwhit[2]).Id()};
+ const ca::HitIndex_t ihit[NHits] = {frWData.Hit(iwhit[0]).Id(), frWData.Hit(iwhit[1]).Id(),
+ frWData.Hit(iwhit[2]).Id()};
if (fParameters.DevIsMatchTripletsViaMc()) {
int mc1 = ca::Framework::GetMcTrackIdForCaHit(ihit[0]);
@@ -553,19 +522,19 @@ void TripletConstructor::FitTriplets()
// repeat several times in order to improve the precision
for (int iiter = 0; iiter < nIterations; ++iiter) {
- // fit forward
- {
+
+ auto fitTrack = [&](int startIdx, int endIdx, int step, float energyLossSign) {
fit.SetQp0(T.Qp());
fit.Qp0()(fit.Qp0() > maxQp) = maxQp;
fit.Qp0()(fit.Qp0() < -maxQp) = -maxQp;
- T.Qp() = 0.;
- T.Vi() = 0.;
- int ih0 = 0;
+ int ih0 = startIdx;
T.X() = x[ih0];
T.Y() = y[ih0];
T.Z() = z[ih0];
T.Time() = t[ih0];
+ T.Qp() = 0.;
+ T.Vi() = 0.;
T.ResetErrors(1., 1., 1., 1., 100., (sta[ih0].timeInfo ? dt2[ih0] : 1.e6), 1.e2);
T.C00() = mxy[ih0].Dx2();
@@ -575,52 +544,27 @@ void TripletConstructor::FitTriplets()
T.Ndf() = -ndfTrackModel + 2;
T.NdfTime() = sta[ih0].timeInfo ? 0 : -1;
- // add the target constraint
- fit.FilterWithTargetAtLine(fParameters.GetTargetPositionZ(), frWData.TargetMeasurement(), fldTarget);
+ if (startIdx == 0) { //only for the forward fit
+ fit.FilterWithTargetAtLine(fParameters.GetTargetPositionZ(), frWData.TargetMeasurement(), fldTarget);
+ }
- for (int ih = 1; ih < NHits; ++ih) {
+ for (int ih = startIdx + step; ih != endIdx; ih += step) {
fit.Extrapolate(z[ih], fld);
auto radThick = fParameters.GetMaterialThickness(ista[ih], T.X(), T.Y());
fit.MultipleScattering(radThick);
- fit.EnergyLossCorrection(radThick, fvec(-1.f));
+ fit.EnergyLossCorrection(radThick, fvec(energyLossSign));
fit.FilterXY(mxy[ih]);
fit.FilterTime(t[ih], dt2[ih], sta[ih].timeInfo);
}
- }
-
- if (iiter == nIterations - 1) break;
-
- // fit backward
- {
- fit.SetQp0(T.Qp());
- fit.Qp0()(fit.Qp0() > maxQp) = maxQp;
- fit.Qp0()(fit.Qp0() < -maxQp) = -maxQp;
-
- int ih0 = NHits - 1;
- T.X() = x[ih0];
- T.Y() = y[ih0];
- T.Z() = z[ih0];
- T.Time() = t[ih0];
- T.Qp() = 0.;
- T.Vi() = 0.;
+ };
- T.ResetErrors(1., 1., 1., 1., 100., (sta[ih0].timeInfo ? dt2[ih0] : 1.e6), 1.e2);
- T.C00() = mxy[ih0].Dx2();
- T.C10() = mxy[ih0].Dxy();
- T.C11() = mxy[ih0].Dy2();
+ // Fit forward
+ fitTrack(0, NHits, 1, -1.f);
- T.Ndf() = -ndfTrackModel + 2;
- T.NdfTime() = sta[ih0].timeInfo ? 0 : -1;
+ if (iiter == nIterations - 1) break;
- for (int ih = NHits - 2; ih >= 0; --ih) {
- fit.Extrapolate(z[ih], fld);
- auto radThick = fParameters.GetMaterialThickness(ista[ih], T.X(), T.Y());
- fit.MultipleScattering(radThick);
- fit.EnergyLossCorrection(radThick, fvec(1.f));
- fit.FilterXY(mxy[ih]);
- fit.FilterTime(t[ih], dt2[ih], sta[ih].timeInfo);
- }
- }
+ // Fit backward
+ fitTrack(NHits - 1, -1, -1, 1.f);
} // for iiter
fTracks_3[i3] = T;
@@ -639,11 +583,11 @@ void TripletConstructor::FitTriplets()
const ca::Hit& h2 = fInputData.GetHit(ih2);
//const CbmL1MCTrack& mctr = CbmL1::Instance()->GetMcTracks()[mc1];
LOG(info) << "== fitted triplet: "
- << " iter " << frWData.CurrentIterationIndex() << " hits: {" << fIstaL << "/" << ih0 << " " << fIstaM
- << "/" << ih1 << " " << fIstaR << "/" << ih2 << "} xyz: {" << h0.X() << " " << h0.Y() << " " << h0.Z()
- << "}, {" << h1.X() << " " << h1.Y() << " " << h1.Z() << "}, {" << h2.X() << ", " << h2.Y() << ", "
- << h2.Z() << "} chi2 " << T.GetChiSq()[0] << " ndf " << T.Ndf()[0] << " chi2time " << T.ChiSqTime()[0]
- << " ndfTime " << T.NdfTime()[0];
+ << " iter " << frWData.CurrentIteration()->GetName() << " hits: {" << fIstaL << "/" << ih0 << " "
+ << fIstaM << "/" << ih1 << " " << fIstaR << "/" << ih2 << "} xyz: {" << h0.X() << " " << h0.Y() << " "
+ << h0.Z() << "}, {" << h1.X() << " " << h1.Y() << " " << h1.Z() << "}, {" << h2.X() << ", " << h2.Y()
+ << ", " << h2.Z() << "} chi2 " << T.GetChiSq()[0] << " ndf " << T.Ndf()[0] << " chi2time "
+ << T.ChiSqTime()[0] << " ndfTime " << T.NdfTime()[0];
/*
cbm::ca::tools::Debugger::Instance().FillNtuple(
"triplets", mctr.iEvent, frAlgo.fCurrentIterationIndex, ih0, h0.X(), h0.Y(), h0.Z(), ih1, h1.X(), h1.Y(),
@@ -662,9 +606,9 @@ void TripletConstructor::StoreTriplets()
/// Selects good triplets and saves them into fTriplets.
/// Finds neighbouring triplets at the next station.
- Tindex n3 = fHitsM_3.size();
+ const Tindex n3 = fHitsM_3.size();
- bool isMomentumFitted = ((staL().fieldStatus != 0) || (staM().fieldStatus != 0) || (staR().fieldStatus != 0));
+ bool isMomentumFitted = ((fStaL->fieldStatus != 0) || (fStaM->fieldStatus != 0) || (fStaR->fieldStatus != 0));
fTriplets.clear();
fTriplets.reserve(n3);
@@ -676,7 +620,7 @@ void TripletConstructor::StoreTriplets()
// TODO: SG: normalize chi2, separate cuts on time and space
- fscal chi2 = T3.GetChiSq()[0] + T3.GetChiSqTime()[0];
+ const fscal chi2 = T3.GetChiSq()[0] + T3.GetChiSqTime()[0];
const ca::HitIndex_t ihitl = fIhitL;
const ca::HitIndex_t ihitm = fHitsM_3[i3];
@@ -691,13 +635,11 @@ void TripletConstructor::StoreTriplets()
continue;
}
}
-
// assert(std::isfinite(chi2) && chi2 > 0);
if (!std::isfinite(chi2) || chi2 < 0) {
continue;
}
-
//if (!T3.IsEntryConsistent(true, 0)) { continue; }
fscal qp = T3.Qp()[0];
@@ -727,9 +669,9 @@ void TripletConstructor::CollectHits(const TrackParamV& Tr, const int iSta, cons
TrackParamV& T = fFit.Tr();
//LOG(info) << T.chi2[0] ;
T.ChiSq() = 0.;
- // if make it bigger the found hits will be rejected later because of the chi2 cut.
- const fvec Pick_m22 = (fvec(chi2Cut) - T.GetChiSq());
+ // if make it bigger the found hits will be rejected later because of the chi2 cut.
+ const fvec Pick_m22 = (fvec(chi2Cut) - T.GetChiSq());
const fscal timeError2 = T.C55()[0];
const fscal time = T.Time()[0];
@@ -746,10 +688,8 @@ void TripletConstructor::CollectHits(const TrackParamV& Tr, const int iSta, cons
area.DoLoopOverEntireGrid();
}
- ca::HitIndex_t ih = 0;
-
- while (area.GetNextObjectId(ih) && ((int) collectedHits.size() < maxNhits)) { // loop over station hits
-
+ // loop over station hits (index incremented in condition)
+ for (ca::HitIndex_t ih = 0; area.GetNextObjectId(ih) && ((int) collectedHits.size() < maxNhits);) {
if (frWData.IsHitSuppressed(ih)) {
continue;
@@ -760,8 +700,8 @@ void TripletConstructor::CollectHits(const TrackParamV& Tr, const int iSta, cons
LOG(info) << "hit in the search area: " << hit.ToString();
LOG(info) << " check the hit.. ";
}
- if (iMC >= 0) { // match via MC
+ if (iMC >= 0) { // match via MC
if (iMC != ca::Framework::GetMcTrackIdForWindowHit(ih)) {
if constexpr (fDebugCollectHits) {
LOG(info) << " hit mc does not match";
@@ -787,14 +727,9 @@ void TripletConstructor::CollectHits(const TrackParamV& Tr, const int iSta, cons
const fscal z = hit.Z();
// check y-boundaries
- auto [y, C11] = fFit.ExtrapolateLineYdY2(z);
-
- fscal dy_est = sqrt(Pick_m22[0] * C11[0]) + hit.RangeY();
-
- fscal xm = hit.X();
- fscal ym = hit.Y();
-
- const fscal dY = ym - y[0];
+ const auto [y, C11] = fFit.ExtrapolateLineYdY2(z);
+ const fscal dy_est = sqrt(Pick_m22[0] * C11[0]) + hit.RangeY();
+ const fscal dY = hit.Y() - y[0];
if (fabs(dY) > dy_est) {
if constexpr (fDebugCollectHits) {
@@ -804,11 +739,9 @@ void TripletConstructor::CollectHits(const TrackParamV& Tr, const int iSta, cons
}
// check x-boundaries
- auto [x, C00] = fFit.ExtrapolateLineXdX2(z);
-
- fscal dx_est = sqrt(Pick_m22[0] * C00[0]) + hit.RangeX();
-
- const fscal dX = xm - x[0];
+ const auto [x, C00] = fFit.ExtrapolateLineXdX2(z);
+ const fscal dx_est = sqrt(Pick_m22[0] * C00[0]) + hit.RangeX();
+ const fscal dX = hit.X() - x[0];
if (fabs(dX) > dx_est) {
if constexpr (fDebugCollectHits) {
@@ -818,12 +751,10 @@ void TripletConstructor::CollectHits(const TrackParamV& Tr, const int iSta, cons
}
// check chi2
-
- fvec C10 = fFit.ExtrapolateLineDxy(z);
-
ca::MeasurementXy<fvec> mxy(hit.X(), hit.Y(), hit.dX2(), hit.dY2(), hit.dXY(), fvec::One(), fvec::One());
- auto [chi2x, chi2u] = ca::TrackFit::GetChi2XChi2U(mxy, x, y, C00, C10, C11);
+ const fvec C10 = fFit.ExtrapolateLineDxy(z);
+ const auto [chi2x, chi2u] = ca::TrackFit::GetChi2XChi2U(mxy, x, y, C00, C10, C11);
// TODO: adjust the cut, cut on chi2x & chi2u separately
if (!frWData.CurrentIteration()->GetTrackFromTripletsFlag()) {
diff --git a/algo/ca/core/tracking/CaTripletConstructor.h b/algo/ca/core/tracking/CaTripletConstructor.h
index 162e0c5279d8231381e32e469d3d7c9d1f31ca98..2c39596c4b1ac1cc4db9f909e8e5528eb6cbc341 100644
--- a/algo/ca/core/tracking/CaTripletConstructor.h
+++ b/algo/ca/core/tracking/CaTripletConstructor.h
@@ -52,15 +52,15 @@ namespace cbm::algo::ca
/// Destructor
~TripletConstructor() = default;
-
/// ------ FUNCTIONAL PART ------
- void InitStations(int istal, int istam, int istar);
-
void CreateTripletsForHit(int istal, int istam, int istar, ca::HitIndex_t ihl);
const Vector<ca::Triplet>& GetTriplets() const { return fTriplets; }
+ private:
+ bool InitStations(int istal, int istam, int istar);
+
/// Find the doublets. Reformat data in the portion of doublets.
void FindDoublets();
void FitDoublets();
@@ -81,14 +81,6 @@ namespace cbm::algo::ca
void CollectHits(const TrackParamV& Tr, const int iSta, const double chi2Cut, const int iMC,
Vector<ca::HitIndex_t>& collectedHits, int maxNhits);
- private:
- /// left station
- const ca::Station<fvec>& staL() { return *fStaL; }
- /// middle station
- const ca::Station<fvec>& staM() { return *fStaM; }
- /// right station
- const ca::Station<fvec>& staR() { return *fStaR; }
-
private:
const Parameters<fvec>& fParameters; ///< Object of Framework parameters class
const InputData& fInputData; ///< Tracking input data
@@ -101,8 +93,6 @@ namespace cbm::algo::ca
// Number of stations situated in field region (MVD + STS in CBM)
int fNfieldStations{0};
- bool fIsInitialized{false}; ///< is initialized;
-
int fIstaL{-1}; ///< left station index
int fIstaM{-1}; ///< middle station index
int fIstaR{-1}; ///< right station index