diff --git a/algo/ca/core/tracking/CaTrackExtender.cxx b/algo/ca/core/tracking/CaTrackExtender.cxx
index 7af7d187ecb36098eec904f015d1c69e734c6077..6e3ccf6656918097c4b26cce737695dc4552dd34 100644
--- a/algo/ca/core/tracking/CaTrackExtender.cxx
+++ b/algo/ca/core/tracking/CaTrackExtender.cxx
@@ -55,9 +55,9 @@ void TrackExtender::FitBranchFast(const ca::Branch& t, TrackParamV& Tout, const
const int iFirstHit = (direction == FitDirection::kUpstream) ? nHits - 1 : 0;
const int iLastHit = (direction == FitDirection::kUpstream) ? 0 : nHits - 1;
- const ca::Hit& hit0 = fInputData->GetHit(hits[iFirstHit]);
- const ca::Hit& hit1 = fInputData->GetHit(hits[iFirstHit + step]);
- const ca::Hit& hit2 = fInputData->GetHit(hits[iFirstHit + 2 * step]);
+ const ca::Hit& hit0 = frWData->Hit(hits[iFirstHit]);
+ const ca::Hit& hit1 = frWData->Hit(hits[iFirstHit + step]);
+ const ca::Hit& hit2 = frWData->Hit(hits[iFirstHit + 2 * step]);
int ista0 = hit0.Station();
int ista1 = hit1.Station();
@@ -113,7 +113,7 @@ void TrackExtender::FitBranchFast(const ca::Branch& t, TrackParamV& Tout, const
fld.Set(fldB2, fldZ2, fldB1, fldZ1, fldB0, fldZ0);
for (int i = iFirstHit + step; step * i <= step * iLastHit; i += step) {
- const ca::Hit& hit = fInputData->GetHit(hits[i]);
+ const ca::Hit& hit = frWData->Hit(hits[i]);
int ista = hit.Station();
const ca::Station<fvec>& sta = fParameters.GetStation(ista);
@@ -162,9 +162,9 @@ void TrackExtender::FindMoreHits(ca::Branch& t, TrackParamV& Tout, const FitDire
const int iFirstHit = (direction == FitDirection::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
- const ca::Hit& hit1 = fInputData->GetHit(t.Hits()[iFirstHit + step]);
- const ca::Hit& hit2 = fInputData->GetHit(t.Hits()[iFirstHit + 2 * step]);
+ const ca::Hit& hit0 = frWData->Hit(t.Hits()[iFirstHit]); // optimize
+ const ca::Hit& hit1 = frWData->Hit(t.Hits()[iFirstHit + step]);
+ const ca::Hit& hit2 = frWData->Hit(t.Hits()[iFirstHit + 2 * step]);
const int ista0 = hit0.Station();
const int ista1 = hit1.Station();
@@ -267,7 +267,7 @@ void TrackExtender::FindMoreHits(ca::Branch& t, TrackParamV& Tout, const FitDire
const ca::Hit& hit = frWData->Hit(iHit_best);
- newHits.push_back(hit.Id());
+ newHits.push_back(iHit_best);
fit.Extrapolate(hit.Z(), fld);
fit.FilterHit(hit, fmask(sta.timeInfo));
@@ -308,10 +308,9 @@ void TrackExtender::FindMoreHits(ca::Branch& t, TrackParamV& Tout, const FitDire
} // void ca::Framework::FindMoreHits
/// Try to extrapolate and find additional hits on other stations
-fscal TrackExtender::ExtendBranch(ca::Branch& t, const InputData& input, WindowData& wData)
+fscal TrackExtender::ExtendBranch(ca::Branch& t, WindowData& wData)
{
- fInputData = &input;
- frWData = &wData;
+ frWData = &wData;
// const unsigned int minNHits = 3;
TrackParamV T;
diff --git a/algo/ca/core/tracking/CaTrackExtender.h b/algo/ca/core/tracking/CaTrackExtender.h
index be2dfea9c050f8ef792db655c85acc7893d12fe3..48260c0b8432ffd142c0d41ede908dfcb67a118f 100644
--- a/algo/ca/core/tracking/CaTrackExtender.h
+++ b/algo/ca/core/tracking/CaTrackExtender.h
@@ -51,7 +51,7 @@ namespace cbm::algo::ca
/// Find additional hits for existing track both downstream and upstream
/// \return chi2
- fscal ExtendBranch(ca::Branch& t, const ca::InputData& input, WindowData& wData);
+ fscal ExtendBranch(ca::Branch& t, WindowData& wData);
private:
///-------------------------------
@@ -88,7 +88,6 @@ namespace cbm::algo::ca
/// Data members
const Parameters<fvec>& fParameters; ///< Object of Framework parameters class
- const InputData* fInputData; ///< Tracking input data
WindowData* frWData;
fscal fDefaultMass{constants::phys::MuonMass}; ///< mass of the propagated particle [GeV/c2]
};
diff --git a/algo/ca/core/tracking/CaTrackFinder.cxx b/algo/ca/core/tracking/CaTrackFinder.cxx
index 9326f1d97b0f37967d058d951f838fff8748e456..c3650fe62102f80a8da91a14d1f58b294f093ec0 100644
--- a/algo/ca/core/tracking/CaTrackFinder.cxx
+++ b/algo/ca/core/tracking/CaTrackFinder.cxx
@@ -114,9 +114,7 @@ TrackFinder::Output_t TrackFinder::FindTracks(const InputData& input, TimesliceH
fStatNhitsTotal = 0;
// calculate possible event time for the hits (fHitTimeInfo array)
- const int nDataStreams = input.GetNdataStreams();
-
- for (int iStream = 0; iStream < nDataStreams; ++iStream) {
+ for (int iStream = 0; iStream < input.GetNdataStreams(); ++iStream) {
fscal maxTimeBeforeHit = std::numeric_limits<fscal>::lowest();
const int nStreamHits = input.GetStreamNhits(iStream);
@@ -203,38 +201,35 @@ TrackFinder::Output_t TrackFinder::FindTracks(const InputData& input, TimesliceH
}
// int nWindowsThread = nWindows / fNofThreads;
- fscal windowDelta = fWindowLength - fWindowOverlap - fWindowMargin;
+ const fscal windowDelta = fWindowLength - fWindowOverlap - fWindowMargin;
//LOG(info) << "CA: estimated number of time windows: " << nWindows;
- std::vector<fscal> vWindowStartThread(fNofThreads), vWindowEndThread(fNofThreads);
-
- vWindowStartThread[0] = fStatTsStart;
+ std::vector<std::pair<fscal, fscal>> vWindowRangeThread(fNofThreads);
{ // Estimation of number of hits in time windows
//Timer time;
//time.Start();
- const int nSt = fParameters.GetNstationsActive();
+ const HitIndex_t nHitsTot = input.GetNhits();
+ const int nSt = fParameters.GetNstationsActive();
// Count number of hits per window and station
- const double a = fStatTsStart + fWindowOverlap + fWindowMargin;
- const double b = fWindowLength - fWindowOverlap - fWindowMargin;
- HitIndex_t nHitsTot = input.GetNhits();
+ const double a = fStatTsStart + fWindowOverlap + fWindowMargin;
+ const double b = fWindowLength - fWindowOverlap - fWindowMargin;
std::vector<HitIndex_t> nHitsWindowSta(nWindows * nSt, 0);
+
for (HitIndex_t iHit = 0; iHit < nHitsTot; ++iHit) {
- const auto& hit = input.GetHit(iHit);
- //const auto& timeInfo = fHitTimeInfo[iHit];
- size_t iWindow = static_cast<size_t>((hit.T() - a) / b);
+ const auto& hit = input.GetHit(iHit);
+ const size_t iWindow = static_cast<size_t>((hit.T() - a) / b);
if (iWindow >= static_cast<size_t>(nWindows)) {
LOG(error) << "ca: Hit out of range: iHit = " << iHit << ", time = " << hit.T() * 1.e-6
<< " ms, window = " << iWindow;
continue;
}
-
++nHitsWindowSta[hit.Station() + iWindow * nSt];
}
// Remove hits from the "monster events"
if (ca::TrackingMode::kMcbm == fTrackingMode) {
- auto maxNofHitsSta = static_cast<HitIndex_t>(50 * fWindowLength / 1.e3);
+ const auto maxNofHitsSta = static_cast<HitIndex_t>(50 * fWindowLength / 1.e3);
for (auto& content : nHitsWindowSta) {
if (content > maxNofHitsSta) {
content = 0;
@@ -251,32 +246,33 @@ TrackFinder::Output_t TrackFinder::FindTracks(const InputData& input, TimesliceH
}
// Get time range for threads
- HitIndex_t nHitsPerThread = nHitsCollected / fNofThreads;
- auto windowIt = nHitsWindow.begin();
- size_t iWbegin = 0;
+ const HitIndex_t nHitsPerThread = nHitsCollected / fNofThreads;
+ auto windowIt = nHitsWindow.begin();
+ vWindowRangeThread[0].first = fStatTsStart;
for (int iTh = 1; iTh < fNofThreads; ++iTh) {
- windowIt = std::lower_bound(windowIt, nHitsWindow.end(), iTh * nHitsPerThread);
- iWbegin = std::distance(nHitsWindow.begin(), windowIt) + 1;
- vWindowStartThread[iTh] = fStatTsStart + iWbegin * windowDelta;
- vWindowEndThread[iTh - 1] = vWindowStartThread[iTh];
+ windowIt = std::lower_bound(windowIt, nHitsWindow.end(), iTh * nHitsPerThread);
+ const size_t iWbegin = std::distance(nHitsWindow.begin(), windowIt) + 1;
+ vWindowRangeThread[iTh].first = fStatTsStart + iWbegin * windowDelta;
+ vWindowRangeThread[iTh - 1].second = vWindowRangeThread[iTh].first;
}
+ vWindowRangeThread[fNofThreads - 1].second = fStatTsEnd;
+
//time.Stop();
//LOG(info) << "Thread boarders estimation time: " << time.GetTotalMs() << " ms";
LOG(debug) << "Fraction of hits from monster events: " << static_cast<double>(nHitsTot - nHitsCollected) / nHitsTot;
}
-
// cut data into sub-timeslices and process them one by one
//for (int iThread = 0; iThread < fNofThreads; ++iThread) {
// vWindowStartThread[iThread] = fStatTsStart + iThread * nWindowsThread * windowDelta;
// vWindowEndThread[iThread] =
// vWindowStartThread[iThread] + nWindowsThread * windowDelta + fWindowOverlap + fWindowMargin;
//}
- vWindowEndThread[fNofThreads - 1] = fStatTsEnd;
for (int iThread = 0; iThread < fNofThreads; ++iThread) {
- double start = vWindowStartThread[iThread] * 1.e-6;
- double end = vWindowEndThread[iThread] * 1.e-6;
+ auto& entry = vWindowRangeThread[iThread];
+ double start = entry.first * 1.e-6;
+ double end = entry.second * 1.e-6;
LOG(trace) << "Thread: " << iThread << " from " << start << " ms to " << end << " ms (delta = " << end - start
<< " ms)";
}
@@ -287,8 +283,8 @@ TrackFinder::Output_t TrackFinder::FindTracks(const InputData& input, TimesliceH
fMonitorData.StopTimer(ETimer::PrepareTimeslice);
// Save tracks
if (fNofThreads == 1) {
- this->FindTracksThread(input, 0, std::ref(vWindowStartThread[0]), std::ref(vWindowEndThread[0]),
- std::ref(vStatNwindows[0]), std::ref(vStatNhitsProcessed[0]));
+ this->FindTracksThread(input, 0, std::ref(vWindowRangeThread[0]), std::ref(vStatNwindows[0]),
+ std::ref(vStatNhitsProcessed[0]));
fMonitorData.StartTimer(ETimer::StoreTracksFinal);
recoTracks = std::move(fvRecoTracks[0]);
recoHits = std::move(fvRecoHitIndices[0]);
@@ -299,8 +295,8 @@ TrackFinder::Output_t TrackFinder::FindTracks(const InputData& input, TimesliceH
vThreadList.reserve(fNofThreads);
for (int iTh = 0; iTh < fNofThreads; ++iTh) {
vThreadList.emplace_back(&TrackFinder::FindTracksThread, this, std::ref(input), iTh,
- std::ref(vWindowStartThread[iTh]), std::ref(vWindowEndThread[iTh]),
- std::ref(vStatNwindows[iTh]), std::ref(vStatNhitsProcessed[iTh]));
+ std::ref(vWindowRangeThread[iTh]), std::ref(vStatNwindows[iTh]),
+ std::ref(vStatNhitsProcessed[iTh]));
}
for (auto& th : vThreadList) {
if (th.joinable()) {
@@ -350,29 +346,32 @@ TrackFinder::Output_t TrackFinder::FindTracks(const InputData& input, TimesliceH
// ---------------------------------------------------------------------------------------------------------------------
//
-void TrackFinder::FindTracksThread(const InputData& input, int iThread, fscal& windowStart, fscal& windowEnd,
+void TrackFinder::FindTracksThread(const InputData& input, int iThread, std::pair<fscal, fscal>& windowRange,
int& statNwindows, int& statNhitsProcessed)
{
//std::stringstream filename;
//filename << "./dbg_caTrackFinder::FindTracksThread_" << iThread << ".txt";
//std::ofstream out(filename.str());
- auto& monitor = fvMonitorDataThread[iThread];
Timer timer;
timer.Start();
+
+ auto& monitor = fvMonitorDataThread[iThread];
monitor.StartTimer(ETimer::TrackingThread);
monitor.StartTimer(ETimer::PrepareThread);
// Init vectors
- auto& wData = fvWData[iThread];
+ auto& tracks = fvRecoTracks[iThread];
+ auto& hitIndices = fvRecoHitIndices[iThread];
+ auto& wData = fvWData[iThread];
{
const int nStations = fParameters.GetNstationsActive();
const size_t nHitsTot = input.GetNhits();
const size_t nHitsExpected = 2 * nHitsTot;
const size_t nTracksExpected = 2 * nHitsTot / nStations;
- fvRecoTracks[iThread].clear();
- fvRecoTracks[iThread].reserve(nTracksExpected / fNofThreads);
- fvRecoHitIndices[iThread].clear();
- fvRecoHitIndices[iThread].reserve(nHitsExpected / fNofThreads);
+ tracks.clear();
+ tracks.reserve(nTracksExpected / fNofThreads);
+ hitIndices.clear();
+ hitIndices.reserve(nHitsExpected / fNofThreads);
if (iThread != 0) {
wData.HitKeyFlags() = fvWData[0].HitKeyFlags();
}
@@ -382,25 +381,26 @@ void TrackFinder::FindTracksThread(const InputData& input, int iThread, fscal& w
}
}
- const int nDataStreams = input.GetNdataStreams();
- std::vector<HitIndex_t> sliceFirstHit(nDataStreams, 0);
- std::vector<HitIndex_t> sliceLastHit(nDataStreams, 0);
+ // Begin and end index of hit-range for streams
+ std::vector<std::pair<HitIndex_t, HitIndex_t>> streamHitRanges(input.GetNdataStreams(), {0, 0});
// Define first hit, skip all the hits, which are before the first window
- for (int iStream = 0; iStream < nDataStreams; ++iStream) {
- sliceFirstHit[iStream] = input.GetStreamStartIndex(iStream);
- sliceLastHit[iStream] = input.GetStreamStopIndex(iStream);
- for (HitIndex_t caHitId = sliceFirstHit[iStream]; caHitId < sliceLastHit[iStream]; ++caHitId) {
+ for (size_t iStream = 0; iStream < streamHitRanges.size(); ++iStream) {
+ auto& range = streamHitRanges[iStream];
+ range.first = input.GetStreamStartIndex(iStream);
+ range.second = input.GetStreamStopIndex(iStream);
+
+ for (HitIndex_t caHitId = range.first; caHitId < range.second; ++caHitId) {
const ca::Hit& h = input.GetHit(caHitId);
if (wData.IsHitKeyUsed(h.FrontKey()) || wData.IsHitKeyUsed(h.BackKey())) {
continue;
}
const CaHitTimeInfo& info = fHitTimeInfo[caHitId];
- if (info.fMaxTimeBeforeHit < windowStart) {
- sliceFirstHit[iStream] = caHitId + 1;
+ if (info.fMaxTimeBeforeHit < windowRange.first) {
+ range.first = caHitId + 1;
}
- if (info.fMinTimeAfterHit > windowEnd) {
- sliceLastHit[iStream] = caHitId;
+ if (info.fMinTimeAfterHit > windowRange.second) {
+ range.second = caHitId;
}
}
}
@@ -408,58 +408,60 @@ void TrackFinder::FindTracksThread(const InputData& input, int iThread, fscal& w
int statLastLogTimeChunk = -1;
// Track finder algorithm for the time window
- ca::TrackFinderWindow trackFinderWindow(fParameters, fDefaultMass, fTrackingMode, fvMonitorDataThread[iThread]);
+ ca::TrackFinderWindow trackFinderWindow(fParameters, fDefaultMass, fTrackingMode, monitor);
trackFinderWindow.InitTimeslice(input.GetNhitKeys());
monitor.StopTimer(ETimer::PrepareThread);
- bool areUntouchedDataLeft = true; // is the whole TS processed
-
- while (areUntouchedDataLeft) {
- fvMonitorDataThread[iThread].IncrementCounter(ECounter::SubTS);
+ while (true) {
+ monitor.IncrementCounter(ECounter::SubTS);
// select the sub-slice hits
for (int iS = 0; iS < fParameters.GetNstationsActive(); ++iS) {
wData.TsHitIndices(iS).clear();
}
- areUntouchedDataLeft = false;
+ bool areUntouchedDataLeft = false; // is the whole TS processed
// TODO: SG: skip empty regions and start the subslice with the earliest hit
statNwindows++;
//out << statNwindows << ' ';
- int statNwindowHits = 0;
+ monitor.StartTimer(ETimer::PrepareWindow);
- fvMonitorDataThread[iThread].StartTimer(ETimer::PrepareWindow);
- for (int iStream = 0; iStream < nDataStreams; ++iStream) {
- for (ca::HitIndex_t caHitId = sliceFirstHit[iStream]; caHitId < sliceLastHit[iStream]; ++caHitId) {
- const CaHitTimeInfo& info = fHitTimeInfo[caHitId];
- const ca::Hit& h = input.GetHit(caHitId);
+ for (auto& range : streamHitRanges) {
+ for (HitIndex_t caHitId = range.first; caHitId < range.second; ++caHitId) {
+ const ca::Hit& h = input.GetHit(caHitId);
if (wData.IsHitKeyUsed(h.FrontKey()) || wData.IsHitKeyUsed(h.BackKey())) {
// the hit is already reconstructed
continue;
}
- if (info.fEventTimeMin > windowStart + fWindowLength) { // the hit is too late for the sub slice
+ const CaHitTimeInfo& info = fHitTimeInfo[caHitId];
+ if (info.fEventTimeMax < windowRange.first) {
+ // the hit belongs to previous sub-slices
+ continue;
+ }
+ if (info.fMinTimeAfterHit > windowRange.first + fWindowLength) {
+ // this hit and all later hits are out of the sub-slice
areUntouchedDataLeft = true;
- if (info.fMinTimeAfterHit > windowStart + fWindowLength) {
- // this hit and all later hits are out of the sub-slice
- break;
- }
+ break;
+ }
+ if (info.fEventTimeMin > windowRange.first + fWindowLength) {
+ // the hit is too late for the sub slice
+ areUntouchedDataLeft = true;
+ continue;
+ }
+
+ // the hit belongs to the sub-slice
+ wData.TsHitIndices(h.Station()).push_back(caHitId);
+ if (info.fMaxTimeBeforeHit < windowRange.first + fWindowLength - fWindowOverlap) {
+ range.first = caHitId + 1; // this hit and all hits before are before the overlap
}
- else if (windowStart <= info.fEventTimeMax) { // the hit belongs to the sub-slice
- wData.TsHitIndices(h.Station()).push_back(caHitId);
- statNwindowHits++;
- if (info.fMaxTimeBeforeHit < windowStart + fWindowLength - fWindowOverlap) {
- // this hit and all hits before are before the overlap
- sliceFirstHit[iStream] = caHitId + 1;
- }
- } // else the hit has been alread processed in previous sub-slices
}
}
//out << statNwindowHits << ' ';
//if (statNwindowHits == 0) { // Empty window
- // fvMonitorDataThread[iThread].StopTimer(ETimer::PrepareWindow);
+ // monitor.StopTimer(ETimer::PrepareWindow);
// out << 0 << ' ' << 0 << ' ' << 0 << '\n';
// continue;
//}
@@ -470,54 +472,55 @@ void TrackFinder::FindTracksThread(const InputData& input, int iThread, fscal& w
for (int ista = 0; ista < fParameters.GetNstationsActive(); ++ista) {
int nHitsSta = static_cast<int>(wData.TsHitIndices(ista).size());
if (nHitsSta > maxStationHits) {
- statNwindowHits -= nHitsSta;
wData.TsHitIndices(ista).clear();
}
}
}
+ int statNwindowHits = 0;
+ for (int ista = 0; ista < fParameters.GetNstationsActive(); ++ista) {
+ statNwindowHits += wData.TsHitIndices(ista).size();
+ }
statNhitsProcessed += statNwindowHits;
// print the LOG for every 10 ms of data processed
if constexpr (0) {
- int currentChunk = (int) ((windowStart - fStatTsStart) / 10.e6);
+ int currentChunk = (int) ((windowRange.first - fStatTsStart) / 10.e6);
if (!areUntouchedDataLeft || currentChunk > statLastLogTimeChunk) {
statLastLogTimeChunk = currentChunk;
- double dataRead = 100. * (windowStart + fWindowLength - fStatTsStart) / (fStatTsEnd - fStatTsStart);
+ double dataRead = 100. * (windowRange.first + fWindowLength - fStatTsStart) / (fStatTsEnd - fStatTsStart);
if (dataRead > 100.) {
dataRead = 100.;
}
- LOG(debug) << "CA tracker process sliding window N " << statNwindows << ": time " << windowStart / 1.e6
+ LOG(debug) << "CA tracker process sliding window N " << statNwindows << ": time " << windowRange.first / 1.e6
<< " ms + " << fWindowLength / 1.e3 << " us) with " << statNwindowHits << " hits. "
<< " Processing " << dataRead << " % of the TS time and "
<< 100. * statNhitsProcessed / fStatNhitsTotal << " % of TS hits."
- << " Already reconstructed " << fvRecoTracks[iThread].size() << " tracks on thread #" << iThread;
+ << " Already reconstructed " << tracks.size() << " tracks on thread #" << iThread;
}
}
-
//out << statNwindowHits << ' ';
-
- fvMonitorDataThread[iThread].StopTimer(ETimer::PrepareWindow);
+ monitor.StopTimer(ETimer::PrepareWindow);
//Timer trackingInWindow; //DBG
//trackingInWindow.Start();
- fvMonitorDataThread[iThread].StartTimer(ETimer::TrackingWindow);
+ monitor.StartTimer(ETimer::TrackingWindow);
trackFinderWindow.CaTrackFinderSlice(input, wData);
- fvMonitorDataThread[iThread].StopTimer(ETimer::TrackingWindow);
+ monitor.StopTimer(ETimer::TrackingWindow);
//trackingInWindow.Stop();
//out << trackingInWindow.GetTotalMs() << ' ';
// save reconstructed tracks with no hits in the overlap region
- //if (windowStart > 13.23e6 && windowStart < 13.26e6) {
- windowStart += fWindowLength - fWindowOverlap;
+ //if (windowRange.first > 13.23e6 && windowRange.first < 13.26e6) {
+ windowRange.first += fWindowLength - fWindowOverlap;
// we should add hits from reconstructed but not stored tracks to the new sub-timeslice
// we do it in a simple way by extending the tsStartNew
// TODO: only add those hits from the region before tsStartNew that belong to the not stored tracks
//out << fvWData[iThread].RecoHitIndices().size() << ' ';
//out << fvWData[iThread].RecoTracks().size() << '\n';
- fvMonitorDataThread[iThread].StartTimer(ETimer::StoreTracksWindow);
+ monitor.StartTimer(ETimer::StoreTracksWindow);
auto trackFirstHit = wData.RecoHitIndices().begin();
for (const auto& track : wData.RecoTracks()) {
@@ -525,7 +528,7 @@ void TrackFinder::FindTracksThread(const InputData& input, int iThread, fscal& w
const bool isTrackCompletelyInOverlap =
std::all_of(trackFirstHit, trackFirstHit + track.fNofHits, [&](int caHitId) {
CaHitTimeInfo& info = fHitTimeInfo[caHitId];
- return info.fEventTimeMax >= windowStart;
+ return info.fEventTimeMax >= windowRange.first;
});
// Don't save tracks from the overlap region, since they might have additional hits in the next subslice.
@@ -539,27 +542,30 @@ void TrackFinder::FindTracksThread(const InputData& input, int iThread, fscal& w
wData.IsHitKeyUsed(h.BackKey()) = static_cast<int>(useFlag);
if (useFlag) {
- fvRecoHitIndices[iThread].push_back(caHitId);
+ hitIndices.push_back(caHitId);
}
}
if (useFlag) {
- fvRecoTracks[iThread].push_back(track);
+ tracks.push_back(track);
}
trackFirstHit += track.fNofHits;
} // sub-timeslice tracks
- fvMonitorDataThread[iThread].StopTimer(ETimer::StoreTracksWindow);
+ monitor.StopTimer(ETimer::StoreTracksWindow);
- if (windowStart > windowEnd) {
+ if (windowRange.first > windowRange.second) {
break;
}
else {
- windowStart -= fWindowMargin; // do 5 ns margin
+ windowRange.first -= fWindowMargin; // do 5 ns margin
}
- }
- fvMonitorDataThread[iThread].StopTimer(ETimer::TrackingThread);
+ if (!areUntouchedDataLeft) {
+ break;
+ }
+ } // while(true)
+ monitor.StopTimer(ETimer::TrackingThread);
//timer.Stop();
//LOG(info) << "CA: finishing tracking on thread " << iThread << " (time: " << timer.GetTotalMs() << " ms, "
// << "hits processed: " << statNhitsProcessed << ", "
- // << "hits used: " << fvRecoHitIndices[iThread].size() << ')';
+ // << "hits used: " << hitIndices.size() << ')';
}
diff --git a/algo/ca/core/tracking/CaTrackFinder.h b/algo/ca/core/tracking/CaTrackFinder.h
index 21f4d56551d639cdb07ba03abb43d56a06bf6169..5a1cd6d810617c0f11636b6a24256f093697c041 100644
--- a/algo/ca/core/tracking/CaTrackFinder.h
+++ b/algo/ca/core/tracking/CaTrackFinder.h
@@ -57,7 +57,7 @@ namespace cbm::algo::ca
private:
// -------------------------------
// Private methods
- void FindTracksThread(const InputData& input, int iThread, fscal& windowStart, fscal& windowEnd, int& statNwindows,
+ void FindTracksThread(const InputData& input, int iThread, std::pair<fscal, fscal>& windowRange, int& statNwindows,
int& statNhitsProcessed);
// bool checkTripletMatch(const ca::Triplet& l, const ca::Triplet& r, fscal& dchi2) const;
diff --git a/algo/ca/core/tracking/CaTrackFinderWindow.cxx b/algo/ca/core/tracking/CaTrackFinderWindow.cxx
index 6acaa5091c51db548a95d1a9733bfd06f889dd51..7b8e04d7dec6871537dc85f0a8e1c6de66476e53 100644
--- a/algo/ca/core/tracking/CaTrackFinderWindow.cxx
+++ b/algo/ca/core/tracking/CaTrackFinderWindow.cxx
@@ -44,7 +44,6 @@ using Track = cbm::algo::ca::Track;
using namespace cbm::algo::ca;
// ---------------------------------------------------------------------------------------------------------------------
-//
TrackFinderWindow::TrackFinderWindow(const ca::Parameters<fvec>& pars, const fscal mass, const ca::TrackingMode& mode,
ca::TrackingMonitorData& monitorData)
: fParameters(pars)
@@ -55,20 +54,9 @@ TrackFinderWindow::TrackFinderWindow(const ca::Parameters<fvec>& pars, const fsc
, fCloneMerger(pars, mass)
, fTrackFitter(pars, mass, mode)
{
- for (unsigned int i = 0; i < constants::size::MaxNstations; ++i) {
- fvTriplets[i].SetName(std::stringstream() << "TrackFinderWindow::fTriplets[" << i << "]");
- }
}
-
-// ---------------------------------------------------------------------------------------------------------------------
-//
-TrackFinderWindow::~TrackFinderWindow() {}
-
// ---------------------------------------------------------------------------------------------------------------------
-//
-
-
bool TrackFinderWindow::checkTripletMatch(const ca::Triplet& l, const ca::Triplet& r, fscal& dchi2,
WindowData& wData) const
{
@@ -129,7 +117,73 @@ bool TrackFinderWindow::checkTripletMatch(const ca::Triplet& l, const ca::Triple
// **************************************************************************************************
// ---------------------------------------------------------------------------------------------------------------------
-//
+void TrackFinderWindow::CaTrackFinderSlice(const ca::InputData& input, WindowData& wData)
+{
+ // Init windows
+ frMonitorData.StartTimer(ETimer::InitWindow);
+ ReadWindowData(input.GetHits(), wData);
+ frMonitorData.StopTimer(ETimer::InitWindow);
+
+ // Init grids
+ frMonitorData.StartTimer(ETimer::PrepareGrid);
+ PrepareGrid(input.GetHits(), wData);
+ frMonitorData.StopTimer(ETimer::PrepareGrid);
+
+ // Run CA iterations
+ frMonitorData.StartTimer(ETimer::FindTracks);
+ auto& caIterations = fParameters.GetCAIterations();
+ for (auto iter = caIterations.begin(); iter != caIterations.end(); ++iter) {
+
+ // ----- Prepare iteration
+ frMonitorData.StartTimer(ETimer::PrepareIteration);
+ PrepareCAIteration(*iter, wData, iter == caIterations.begin());
+ frMonitorData.StopTimer(ETimer::PrepareIteration);
+
+ // ----- Triplets construction -----
+ frMonitorData.StartTimer(ETimer::ConstructTriplets);
+ ConstructTriplets(wData);
+ frMonitorData.StopTimer(ETimer::ConstructTriplets);
+
+ // ----- Search for neighbouring triplets -----
+ frMonitorData.StartTimer(ETimer::SearchNeighbours);
+ SearchNeighbors(wData);
+ frMonitorData.StopTimer(ETimer::SearchNeighbours);
+
+ // ----- Collect track candidates and create tracks
+ frMonitorData.StartTimer(ETimer::CreateTracks);
+ CreateTracks(wData, *iter, std::get<2>(fTripletData));
+ frMonitorData.StopTimer(ETimer::CreateTracks);
+
+ // ----- Suppress strips of suppressed hits
+ frMonitorData.StartTimer(ETimer::SuppressHitKeys);
+ for (unsigned int ih = 0; ih < wData.Hits().size(); ih++) {
+ if (wData.IsHitSuppressed(ih)) {
+ const ca::Hit& hit = wData.Hit(ih);
+ wData.IsHitKeyUsed(hit.FrontKey()) = 1;
+ wData.IsHitKeyUsed(hit.BackKey()) = 1;
+ }
+ }
+ frMonitorData.StopTimer(ETimer::SuppressHitKeys);
+ } // ---- Loop over Track Finder iterations: END ----//
+ frMonitorData.StopTimer(ETimer::FindTracks);
+
+ // Fit tracks
+ frMonitorData.StartTimer(ETimer::FitTracks);
+ fTrackFitter.FitCaTracks(input, wData);
+ frMonitorData.StopTimer(ETimer::FitTracks);
+
+ // Merge clones
+ frMonitorData.StartTimer(ETimer::MergeClones);
+ fCloneMerger.Exec(input, wData);
+ frMonitorData.StopTimer(ETimer::MergeClones);
+
+ // Fit tracks
+ frMonitorData.StartTimer(ETimer::FitTracks);
+ fTrackFitter.FitCaTracks(input, wData);
+ frMonitorData.StopTimer(ETimer::FitTracks);
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
void TrackFinderWindow::ReadWindowData(const Vector<Hit>& hits, WindowData& wData)
{
int nHits = 0;
@@ -165,17 +219,9 @@ void TrackFinderWindow::ReadWindowData(const Vector<Hit>& hits, WindowData& wDat
wData.RecoHitIndices().reserve(2 * nHits);
}
-void TrackFinderWindow::CaTrackFinderSlice(const ca::InputData& input, WindowData& wData)
+// ---------------------------------------------------------------------------------------------------------------------
+void TrackFinderWindow::PrepareGrid(const Vector<Hit>& hits, WindowData& wData)
{
- /********************************/ /**
- * CATrackFinder routine setting
- ***********************************/
-
- frMonitorData.StartTimer(ETimer::InitWindow);
- ReadWindowData(input.GetHits(), wData);
- frMonitorData.StopTimer(ETimer::InitWindow);
-
- frMonitorData.StartTimer(ETimer::PrepareGrid);
for (int iS = 0; iS < fParameters.GetNstationsActive(); ++iS) {
fscal lasttime = std::numeric_limits<fscal>::infinity();
@@ -186,11 +232,10 @@ void TrackFinderWindow::CaTrackFinderSlice(const ca::InputData& input, WindowDat
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));
+ const ca::Hit& h = hits[wData.TsHitIndex(iS, ih)];
gridMinX = std::min(gridMinX, h.X());
gridMinY = std::min(gridMinY, h.Y());
-
gridMaxX = std::max(gridMaxX, h.X());
gridMaxY = std::max(gridMaxY, h.Y());
@@ -227,402 +272,372 @@ void TrackFinderWindow::CaTrackFinderSlice(const ca::InputData& input, WindowDat
wData.HitKeyFlags());
/* clang-format on */
}
- frMonitorData.StopTimer(ETimer::PrepareGrid);
-
- // ----- Find tracks -----
- //
- frMonitorData.StartTimer(ETimer::FindTracks);
- auto& caIterations = fParameters.GetCAIterations();
-
- for (auto iter = caIterations.begin(); iter != caIterations.end(); ++iter) {
+}
- // ----- Prepare iteration
- //
- frMonitorData.StartTimer(ETimer::PrepareIteration);
- const auto& caIteration = *iter; // Dereferencing iterator to get the value
- wData.SetCurrentIteration(&caIteration);
+// ---------------------------------------------------------------------------------------------------------------------
+void TrackFinderWindow::PrepareCAIteration(const ca::Iteration& caIteration, WindowData& wData, const bool isFirst)
+{
+ wData.SetCurrentIteration(&caIteration);
- // 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());
- }
+ // Check if it is not the first element
+ if (!isFirst) {
+ for (int ista = 0; ista < fParameters.GetNstationsActive(); ++ista) {
+ wData.Grid(ista).RemoveUsedHits(wData.Hits(), wData.HitKeyFlags());
}
+ }
- // --> 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);
- }
+ // --> frAlgo.fIsWindowHitSuppressed.reset(frAlgo.fWindowHits.size(), 0);
+ wData.ResetHitSuppressionFlags(); // TODO: ??? No effect?
- // #pragma omp task
- {
- // --- SET PARAMETERS FOR THE ITERATION ---
- // define the target
- const fscal SigmaTargetX = caIteration.GetTargetPosSigmaX();
- const fscal SigmaTargetY = caIteration.GetTargetPosSigmaY(); // target constraint [cm]
+ // --- 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);
- }
+ // 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);
+// ---------------------------------------------------------------------------------------------------------------------
+void TrackFinderWindow::ConstructTriplets(WindowData& wData)
+{
+ auto& [vHitFirstTriplet, vHitNofTriplets, vTriplets] = fTripletData;
+ vHitFirstTriplet.reset(wData.Hits().size(), 0); /// link hit -> first triplet { hit, *, *}
+ vHitNofTriplets.reset(wData.Hits().size(), 0); /// link hit ->n triplets { hit, *, *}
- // ----- Triplets construction -----
- //
- frMonitorData.StartTimer(ETimer::ConstructTriplets);
- ca::TripletConstructor constructor(fParameters, input, wData, fDefaultMass, fTrackingMode);
+ ca::TripletConstructor constructor(fParameters, wData, fDefaultMass, fTrackingMode);
- // indices of the two neighbouring station, taking into account allowed gaps
- std::vector<std::pair<int, int>> staPattern;
- for (int gap = 0; gap <= wData.CurrentIteration()->GetMaxStationGap(); gap++) {
- for (int i = 0; i <= gap; i++) {
- staPattern.push_back(std::make_pair(1 + i, 2 + gap));
- }
+ // prepare triplet storage
+ for (int j = 0; j < fParameters.GetNstationsActive(); j++) {
+ const size_t nHitsStation = wData.TsHitIndices(j).size();
+ vTriplets[j].clear();
+ vTriplets[j].reserve(2 * nHitsStation);
+ }
+
+ // indices of the two neighbouring station, taking into account allowed gaps
+ std::vector<std::pair<int, int>> staPattern;
+ for (int gap = 0; gap <= wData.CurrentIteration()->GetMaxStationGap(); gap++) {
+ for (int i = 0; i <= gap; i++) {
+ staPattern.push_back(std::make_pair(1 + i, 2 + gap));
}
+ }
- 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);
- 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(),
- constructor.GetTriplets().end());
- }
- vHitFirstTriplet[ihitl] = PackTripletId(istal, oldSize);
- vHitNofTriplets[ihitl] = fvTriplets[istal].size() - oldSize;
+ for (int istal = fParameters.GetNstationsActive() - 2; istal >= wData.CurrentIteration()->GetFirstStationIndex();
+ istal--) {
+ // start with downstream chambers
+ const auto& grid = wData.Grid(istal);
+ for (auto& entry : grid.GetEntries()) {
+ ca::HitIndex_t ihitl = entry.GetObjectId();
+ const size_t oldSize = vTriplets[istal].size();
+ for (auto& pattern : staPattern) {
+ constructor.CreateTripletsForHit(fvTriplets, istal, istal + pattern.first, istal + pattern.second, ihitl);
+ vTriplets[istal].insert(vTriplets[istal].end(), fvTriplets.begin(), fvTriplets.end());
}
- } // istal
- frMonitorData.StopTimer(ETimer::ConstructTriplets);
-
- // ----- Search for neighbouring triplets -----
- //
- frMonitorData.StartTimer(ETimer::SearchNeighbours);
+ vHitFirstTriplet[ihitl] = PackTripletId(istal, oldSize);
+ vHitNofTriplets[ihitl] = vTriplets[istal].size() - oldSize;
+ }
+ } // istal
+}
- for (int istal = fParameters.GetNstationsActive() - 2; istal >= iStFirst; istal--) {
- // start with downstream chambers
+// ---------------------------------------------------------------------------------------------------------------------
+void TrackFinderWindow::SearchNeighbors(WindowData& wData)
+{
+ auto& [vHitFirstTriplet, vHitNofTriplets, vTriplets] = fTripletData;
- for (ca::Triplet& tr : fvTriplets[istal]) {
- unsigned int nNeighbours = vHitNofTriplets[tr.GetMHit()];
- unsigned int neighLocation = vHitFirstTriplet[tr.GetMHit()];
- unsigned int neighStation = TripletId2Station(neighLocation);
- unsigned int neighTriplet = TripletId2Triplet(neighLocation);
+ for (int istal = fParameters.GetNstationsActive() - 2; istal >= wData.CurrentIteration()->GetFirstStationIndex();
+ istal--) {
+ // start with downstream chambers
- if (nNeighbours > 0) {
- assert((int) neighStation >= istal + 1
- && (int) neighStation <= istal + 1 + wData.CurrentIteration()->GetMaxStationGap());
- }
+ for (ca::Triplet& tr : vTriplets[istal]) {
+ unsigned int nNeighbours = vHitNofTriplets[tr.GetMHit()];
+ unsigned int neighLocation = vHitFirstTriplet[tr.GetMHit()];
+ unsigned int neighStation = TripletId2Station(neighLocation);
+ unsigned int neighTriplet = TripletId2Triplet(neighLocation);
- unsigned char level = 0;
+ if (nNeighbours > 0) {
+ assert((int) neighStation >= istal + 1
+ && (int) neighStation <= istal + 1 + wData.CurrentIteration()->GetMaxStationGap());
+ }
- for (unsigned int iN = 0; iN < nNeighbours; ++iN, ++neighTriplet, ++neighLocation) {
+ unsigned char level = 0;
- ca::Triplet& neighbour = fvTriplets[neighStation][neighTriplet];
+ for (unsigned int iN = 0; iN < nNeighbours; ++iN, ++neighTriplet, ++neighLocation) {
- fscal dchi2 = 0.;
- if (!checkTripletMatch(tr, neighbour, dchi2, wData)) continue;
+ ca::Triplet& neighbour = vTriplets[neighStation][neighTriplet];
- if (tr.GetFNeighbour() == 0) {
- tr.SetFNeighbour(neighLocation);
- }
- tr.SetNNeighbours(neighLocation - tr.GetFNeighbour() + 1);
+ fscal dchi2 = 0.;
+ if (!checkTripletMatch(tr, neighbour, dchi2, wData)) continue;
- level = std::max(level, static_cast<unsigned char>(neighbour.GetLevel() + 1));
+ if (tr.GetFNeighbour() == 0) {
+ tr.SetFNeighbour(neighLocation);
}
- tr.SetLevel(level);
- } // neighbour search
-
- frMonitorData.IncrementCounter(ECounter::Triplet, fvTriplets[istal].size());
+ tr.SetNNeighbours(neighLocation - tr.GetFNeighbour() + 1);
- } // istal
- frMonitorData.StopTimer(ETimer::SearchNeighbours);
+ level = std::max(level, static_cast<unsigned char>(neighbour.GetLevel() + 1));
+ }
+ tr.SetLevel(level);
+ }
+ frMonitorData.IncrementCounter(ECounter::Triplet, vTriplets[istal].size());
+ }
+}
- // ----- Collect track candidates and create tracks
- //
- frMonitorData.StartTimer(ETimer::CreateTracks);
+// ---------------------------------------------------------------------------------------------------------------------
+void TrackFinderWindow::CreateTracks(WindowData& wData, const ca::Iteration& caIteration, TripletArray_t& vTriplets)
+{
+ // min level to start triplet. So min track length = min_level+3.
+ const int min_level = wData.CurrentIteration()->GetTrackFromTripletsFlag()
+ ? 0
+ : std::min(caIteration.GetMinNhits(), caIteration.GetMinNhitsStation0()) - 3;
+
+ // collect consequtive: the longest tracks, shorter, more shorter and so on
+ for (int firstTripletLevel = fParameters.GetNstationsActive() - 3; firstTripletLevel >= min_level;
+ firstTripletLevel--) {
+ // choose length in triplets number - firstTripletLevel - the maximum possible triplet level among all triplets in the searched candidate
+ CreateTrackCandidates(wData, vTriplets, min_level, firstTripletLevel);
+ DoCompetitionLoop(wData);
+ SelectTracks(wData);
+ }
+}
- // min level to start triplet. So min track length = min_level+3.
- const int min_level = wData.CurrentIteration()->GetTrackFromTripletsFlag()
- ? 0
- : std::min(caIteration.GetMinNhits(), caIteration.GetMinNhitsStation0()) - 3;
+// ---------------------------------------------------------------------------------------------------------------------
+void TrackFinderWindow::CreateTrackCandidates(WindowData& wData, TripletArray_t& vTriplets, const int min_level,
+ const int firstTripletLevel)
+{
+ // how many levels to check
+ int nlevel = (fParameters.GetNstationsActive() - 2) - firstTripletLevel + 1;
- ca::Branch new_tr[constants::size::MaxNstations];
- ca::Branch best_tr;
+ const unsigned char min_best_l =
+ (firstTripletLevel > min_level) ? firstTripletLevel + 2 : min_level + 3; // loose maximum
- // collect consequtive: the longest tracks, shorter, more shorter and so on
- for (int firstTripletLevel = fParameters.GetNstationsActive() - 3; firstTripletLevel >= min_level;
- firstTripletLevel--) {
- // choose length in triplets number - firstTripletLevel - the maximum possible triplet level among all triplets in the searched candidate
+ // Uses persistent field to save memory allocations.
+ // fNewTr is only used here!
+ ca::Branch(&new_tr)[constants::size::MaxNstations] = fNewTr;
- // how many levels to check
- int nlevel = (fParameters.GetNstationsActive() - 2) - firstTripletLevel + 1;
+ fvTrackCandidates.clear();
+ fvTrackCandidates.reserve(wData.Hits().size() / 10);
- const unsigned char min_best_l =
- (firstTripletLevel > min_level) ? firstTripletLevel + 2 : min_level + 3; // loose maximum
+ for (const auto& h : wData.Hits()) {
+ fvHitKeyToTrack[h.FrontKey()] = -1;
+ fvHitKeyToTrack[h.BackKey()] = -1;
+ }
- fvTrackCandidates.clear();
- fvTrackCandidates.reserve(wData.Hits().size() / 10);
+ //== Loop over triplets with the required level, find and store track candidates
+ for (int istaF = wData.CurrentIteration()->GetFirstStationIndex();
+ istaF <= fParameters.GetNstationsActive() - 3 - firstTripletLevel; ++istaF) {
- for (const auto& h : wData.Hits()) {
- fvHitKeyToTrack[h.FrontKey()] = -1;
- fvHitKeyToTrack[h.BackKey()] = -1;
- }
+ if (--nlevel == 0) break; //TODO: SG: this is not needed
- //== Loop over triplets with the required level, find and store track candidates
+ for (ca::Triplet& first_trip : vTriplets[istaF]) {
- for (int istaF = iStFirst; istaF <= fParameters.GetNstationsActive() - 3 - firstTripletLevel; ++istaF) {
+ const auto& fstTripLHit = wData.Hit(first_trip.GetLHit());
+ if (wData.IsHitKeyUsed(fstTripLHit.FrontKey()) || wData.IsHitKeyUsed(fstTripLHit.BackKey())) {
+ continue;
+ }
- if (--nlevel == 0) break; //TODO: SG: this is not needed
+ // skip track candidates that are too short
- for (ca::Triplet& first_trip : fvTriplets[istaF]) {
+ int minNhits = wData.CurrentIteration()->GetMinNhits();
- const auto& fstTripLHit = wData.Hit(first_trip.GetLHit());
- if (wData.IsHitKeyUsed(fstTripLHit.FrontKey()) || wData.IsHitKeyUsed(fstTripLHit.BackKey())) {
- continue;
- }
+ if (fstTripLHit.Station() == 0) {
+ minNhits = wData.CurrentIteration()->GetMinNhitsStation0();
+ }
+ if (wData.CurrentIteration()->GetTrackFromTripletsFlag()) {
+ minNhits = 0;
+ }
- // skip track candidates that are too short
+ if (3 + first_trip.GetLevel() < minNhits) {
+ continue;
+ }
- int minNhits = wData.CurrentIteration()->GetMinNhits();
+ // Collect triplets, which can start a track with length equal to firstTipletLevel + 3. This cut suppresses
+ // ghost tracks, but does not affect the efficiency
+ if (first_trip.GetLevel() < firstTripletLevel) {
+ continue;
+ }
- if (fstTripLHit.Station() == 0) {
- minNhits = wData.CurrentIteration()->GetMinNhitsStation0();
- }
- if (wData.CurrentIteration()->GetTrackFromTripletsFlag()) {
- minNhits = 0;
- }
+ ca::Branch curr_tr;
+ curr_tr.AddHit(first_trip.GetLHit());
+ curr_tr.SetChi2(first_trip.GetChi2());
- if (3 + first_trip.GetLevel() < minNhits) {
- continue;
- }
+ ca::Branch best_tr = curr_tr;
- // Collect triplets, which can start a track with length equal to firstTipletLevel + 3. This cut suppresses
- // ghost tracks, but does not affect the efficiency
- if (first_trip.GetLevel() < firstTripletLevel) {
- continue;
- }
+ /// reqursive func to build a tree of possible track-candidates and choose the best
+ CAFindTrack(istaF, best_tr, &first_trip, curr_tr, min_best_l, new_tr, wData, vTriplets);
- ca::Branch curr_tr;
- curr_tr.AddHit(fstTripLHit.Id());
- curr_tr.SetChi2(first_trip.GetChi2());
+ if (best_tr.NofHits() < firstTripletLevel + 2) continue; // loose maximum one hit
- best_tr = curr_tr;
+ if (best_tr.NofHits() < min_level + 3) continue; // should find all hits for min_level
- /// reqursive func to build a tree of possible track-candidates and choose the best
- CAFindTrack(istaF, best_tr, &first_trip, curr_tr, min_best_l, new_tr, wData);
+ if (best_tr.NofHits() < minNhits) {
+ continue;
+ }
- if (best_tr.NofHits() < firstTripletLevel + 2) continue; // loose maximum one hit
+ int ndf = best_tr.NofHits() * 2 - 5;
- if (best_tr.NofHits() < min_level + 3) continue; // should find all hits for min_level
+ // TODO: automatize the NDF calculation
- if (best_tr.NofHits() < minNhits) {
- continue;
- }
+ if (ca::TrackingMode::kGlobal == fTrackingMode || ca::TrackingMode::kMcbm == fTrackingMode) {
+ ndf = best_tr.NofHits() * 2 - 4;
+ }
- int ndf = best_tr.NofHits() * 2 - 5;
+ best_tr.SetChi2(best_tr.Chi2() / ndf);
+ if (fParameters.GetGhostSuppression()) {
+ if (3 == best_tr.NofHits()) {
+ if (!wData.CurrentIteration()->GetPrimaryFlag() && (istaF != 0)) continue; // too /*short*/ non-MAPS track
+ if (wData.CurrentIteration()->GetPrimaryFlag() && (best_tr.Chi2() > 5.0)) continue;
+ }
+ }
+ fvTrackCandidates.push_back(best_tr);
+ ca::Branch& tr = fvTrackCandidates.back();
+ tr.SetStation(istaF);
+ tr.SetId(fvTrackCandidates.size() - 1);
+ tr.SetAlive(true);
+ if constexpr (fDebug) {
+ std::stringstream s;
+ s << "iter " << wData.CurrentIteration()->GetName() << ", track candidate " << fvTrackCandidates.size() - 1
+ << " found, L = " << best_tr.NofHits() << " chi2= " << best_tr.Chi2() << " hits: ";
+ for (auto hitIdLoc : tr.Hits()) {
+ const auto hitId = wData.Hit(hitIdLoc).Id();
+ s << hitId << " (mc " << ca::Framework::GetMcTrackIdForCaHit(hitId) << ") ";
+ }
+ LOG(info) << s.str();
+ }
+ } // itrip
+ } // istaF
- // TODO: automatize the NDF calculation
+ for (ca::Branch& tr : fvTrackCandidates) {
+ tr.SetAlive(false);
+ }
+}
- if (ca::TrackingMode::kGlobal == fTrackingMode || ca::TrackingMode::kMcbm == fTrackingMode) {
- ndf = best_tr.NofHits() * 2 - 4;
- }
+// ---------------------------------------------------------------------------------------------------------------------
+void TrackFinderWindow::DoCompetitionLoop(const WindowData& wData)
+{
+ for (int iComp = 0; (iComp < 100); ++iComp) {
- best_tr.SetChi2(best_tr.Chi2() / ndf);
- if (fParameters.GetGhostSuppression()) {
- if (3 == best_tr.NofHits()) {
- if (!wData.CurrentIteration()->GetPrimaryFlag() && (istaF != 0))
- continue; // too /*short*/ non-MAPS track
- if (wData.CurrentIteration()->GetPrimaryFlag() && (best_tr.Chi2() > 5.0)) continue;
- }
- }
- fvTrackCandidates.push_back(best_tr);
- ca::Branch& tr = fvTrackCandidates.back();
- tr.SetStation(istaF);
- tr.SetId(fvTrackCandidates.size() - 1);
- tr.SetAlive(true);
- if constexpr (fDebug) {
- std::stringstream s;
- 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) << ") ";
- }
- LOG(info) << s.str();
- }
- } // itrip
- } // istaF
+ bool repeatCompetition = false;
- for (ca::Branch& tr : fvTrackCandidates) {
- tr.SetAlive(false);
+ // == Loop over track candidates and mark their strips
+ for (ca::Branch& tr : fvTrackCandidates) {
+ if (tr.IsAlive()) {
+ continue;
}
+ for (auto& hitId : tr.Hits()) {
- bool repeatCompetition = true;
- int iComp = 0;
- for (iComp = 0; (iComp < 100) && repeatCompetition; ++iComp) {
-
- repeatCompetition = false;
-
- // == Loop over track candidates and mark their strips
- for (ca::Branch& tr : fvTrackCandidates) {
- if (tr.IsAlive()) {
- continue;
- }
- for (auto& hitId : tr.Hits()) {
-
- 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 {
- return false;
- }
- }
- else {
- strip = tr.Id();
- }
- return true;
- };
-
- const ca::Hit& h = input.GetHit(hitId);
- if (!updateStrip(fvHitKeyToTrack[h.FrontKey()])) { // front strip
- break;
+ 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();
}
- if (!updateStrip(fvHitKeyToTrack[h.BackKey()])) { // back strip
- break;
- }
- } // loop over hits
- } // itrack
-
- // == Check if some suppressed candidates are still alive
-
- for (ca::Branch& tr : fvTrackCandidates) {
- if (tr.IsAlive()) {
- continue;
- }
-
- tr.SetAlive(true);
- 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
- for (auto hitId : tr.Hits()) {
- const ca::Hit& h = input.GetHit(hitId);
- if (fvHitKeyToTrack[h.FrontKey()] == tr.Id()) {
- fvHitKeyToTrack[h.FrontKey()] = -1;
- }
- if (fvHitKeyToTrack[h.BackKey()] == tr.Id()) {
- fvHitKeyToTrack[h.BackKey()] = -1;
- }
+ else {
+ return false;
}
}
else {
- repeatCompetition = true;
+ strip = tr.Id();
}
- } // itrack
-
- } // competitions
+ return true;
+ };
- // LOG(info) << " N Competitions: " << iComp ;
+ const ca::Hit& h = wData.Hit(hitId);
+ if (!updateStrip(fvHitKeyToTrack[h.FrontKey()])) { // front strip
+ break;
+ }
+ if (!updateStrip(fvHitKeyToTrack[h.BackKey()])) { // back strip
+ break;
+ }
+ } // loop over hits
+ } // itrack
- // ==
+ // == Check if some suppressed candidates are still alive
- for (Tindex iCandidate = 0; iCandidate < (Tindex) fvTrackCandidates.size(); ++iCandidate) {
- ca::Branch& tr = fvTrackCandidates[iCandidate];
+ for (ca::Branch& tr : fvTrackCandidates) {
+ if (tr.IsAlive()) {
+ continue;
+ }
- if constexpr (fDebug) {
- LOG(info) << "iter " << wData.CurrentIteration()->GetName() << ", track candidate " << iCandidate
- << ": alive = " << tr.IsAlive();
- }
- if (!tr.IsAlive()) continue;
+ tr.SetAlive(true);
+ for (const auto& hitIndex : tr.Hits()) {
+ if (!tr.IsAlive()) break;
+ const ca::Hit& h = wData.Hit(hitIndex);
+ tr.SetAlive((fvHitKeyToTrack[h.FrontKey()] == tr.Id()) && (fvHitKeyToTrack[h.BackKey()] == tr.Id()));
+ }
- if (wData.CurrentIteration()->GetExtendTracksFlag()) {
- if (tr.NofHits() < fParameters.GetNstationsActive()) {
- fTrackExtender.ExtendBranch(tr, input, wData);
+ if (!tr.IsAlive()) { // release strips
+ for (auto hitId : tr.Hits()) {
+ const ca::Hit& h = wData.Hit(hitId);
+ if (fvHitKeyToTrack[h.FrontKey()] == tr.Id()) {
+ fvHitKeyToTrack[h.FrontKey()] = -1;
}
- }
-
- for (auto iHit : tr.Hits()) {
- const ca::Hit& hit = input.GetHit(iHit);
- /// used strips are marked
- wData.IsHitKeyUsed(hit.FrontKey()) = 1;
- wData.IsHitKeyUsed(hit.BackKey()) = 1;
- wData.RecoHitIndices().push_back(iHit);
- }
- Track t;
- t.fNofHits = tr.NofHits();
- wData.RecoTracks().push_back(t);
- if (0) { // SG debug
- std::stringstream s;
- s << "store track " << iCandidate << " chi2= " << tr.Chi2() << "\n";
- s << " hits: ";
- for (unsigned int ih = 0; ih < tr.Hits().size(); ih++) {
- s << ca::Framework::GetMcTrackIdForCaHit(tr.Hits()[ih]) << " ";
+ if (fvHitKeyToTrack[h.BackKey()] == tr.Id()) {
+ fvHitKeyToTrack[h.BackKey()] = -1;
}
- LOG(info) << s.str();
}
+ }
+ else {
+ repeatCompetition = true;
+ }
+ } // itrack
- } // tracks
+ if (!repeatCompetition) break;
+ } // competitions
+}
- } // firstTripletLevel
- frMonitorData.StopTimer(ETimer::CreateTracks);
+// ---------------------------------------------------------------------------------------------------------------------
+void TrackFinderWindow::SelectTracks(WindowData& wData)
+{
+ for (Tindex iCandidate = 0; iCandidate < (Tindex) fvTrackCandidates.size(); ++iCandidate) {
+ ca::Branch& tr = fvTrackCandidates[iCandidate];
- frMonitorData.StartTimer(ETimer::SuppressHitKeys);
- // suppress strips of suppressed hits
- for (unsigned int ih = 0; ih < wData.Hits().size(); ih++) {
- if (wData.IsHitSuppressed(ih)) {
- const ca::Hit& hit = wData.Hit(ih);
- wData.IsHitKeyUsed(hit.FrontKey()) = 1;
- wData.IsHitKeyUsed(hit.BackKey()) = 1;
- }
+ if constexpr (fDebug) {
+ LOG(info) << "iter " << wData.CurrentIteration()->GetName() << ", track candidate " << iCandidate
+ << ": alive = " << tr.IsAlive();
}
- frMonitorData.StopTimer(ETimer::SuppressHitKeys);
- } // ---- Loop over Track Finder iterations: END -----------------------------------------------------------//
- frMonitorData.StopTimer(ETimer::FindTracks);
+ if (!tr.IsAlive()) continue;
- // Fit tracks
- frMonitorData.StartTimer(ETimer::FitTracks);
- fTrackFitter.FitCaTracks(input, wData);
- frMonitorData.StopTimer(ETimer::FitTracks);
-
- // Merge clones
- frMonitorData.StartTimer(ETimer::MergeClones);
- fCloneMerger.Exec(input, wData);
- frMonitorData.StopTimer(ETimer::MergeClones);
+ if (wData.CurrentIteration()->GetExtendTracksFlag()) {
+ if (tr.NofHits() < fParameters.GetNstationsActive()) {
+ fTrackExtender.ExtendBranch(tr, wData);
+ }
+ }
- frMonitorData.StartTimer(ETimer::FitTracks);
- fTrackFitter.FitCaTracks(input, wData);
- frMonitorData.StopTimer(ETimer::FitTracks);
+ for (auto iHit : tr.Hits()) {
+ const ca::Hit& hit = wData.Hit(iHit);
+ /// used strips are marked
+ wData.IsHitKeyUsed(hit.FrontKey()) = 1;
+ wData.IsHitKeyUsed(hit.BackKey()) = 1;
+ wData.RecoHitIndices().push_back(hit.Id());
+ }
+ Track t;
+ t.fNofHits = tr.NofHits();
+ wData.RecoTracks().push_back(t);
+ if (0) { // SG debug
+ std::stringstream s;
+ s << "store track " << iCandidate << " chi2= " << tr.Chi2() << "\n";
+ s << " hits: ";
+ for (auto hitLoc : tr.Hits()) {
+ auto hitId = wData.Hit(hitLoc).Id();
+ s << ca::Framework::GetMcTrackIdForCaHit(hitId) << " ";
+ }
+ LOG(info) << s.str();
+ }
+ } // tracks
}
-
/** *************************************************************
* *
* The routine performs recursive search for tracks *
@@ -632,8 +647,10 @@ void TrackFinderWindow::CaTrackFinderSlice(const ca::InputData& input, WindowDat
* *
****************************************************************/
+// ---------------------------------------------------------------------------------------------------------------------
void TrackFinderWindow::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)
+ unsigned char min_best_l, ca::Branch* new_tr, WindowData& wData,
+ TripletArray_t& vTriplets)
/// recursive search for tracks
/// input: @ista - station index, @&best_tr - best track for the privious call
/// output: @&NCalls - number of function calls
@@ -647,10 +664,10 @@ void TrackFinderWindow::CAFindTrack(int ista, ca::Branch& best_tr, const ca::Tri
const auto& hitR = wData.Hit(curr_trip->GetRHit());
if (!(wData.IsHitKeyUsed(hitM.FrontKey()) || wData.IsHitKeyUsed(hitM.BackKey()))) {
- curr_tr.AddHit(hitM.Id());
+ curr_tr.AddHit(curr_trip->GetMHit());
}
if (!(wData.IsHitKeyUsed(hitR.FrontKey()) || wData.IsHitKeyUsed(hitR.BackKey()))) {
- curr_tr.AddHit(hitR.Id());
+ curr_tr.AddHit(curr_trip->GetRHit());
}
//if( curr_tr.NofHits() < min_best_l - 1 ) return; // suppouse that only one hit can be added by extender
@@ -681,7 +698,7 @@ void TrackFinderWindow::CAFindTrack(int ista, ca::Branch& best_tr, const ca::Tri
unsigned int Station = TripletId2Station(ID);
unsigned int Triplet = TripletId2Triplet(ID);
- const ca::Triplet& new_trip = fvTriplets[Station][Triplet];
+ const ca::Triplet& new_trip = vTriplets[Station][Triplet];
fscal dchi2 = 0.;
if (!checkTripletMatch(*curr_trip, new_trip, dchi2, wData)) continue;
@@ -732,12 +749,12 @@ void TrackFinderWindow::CAFindTrack(int ista, ca::Branch& best_tr, const ca::Tri
// add new hit
new_tr[ista] = curr_tr;
- new_tr[ista].AddHit(hitL.Id());
+ new_tr[ista].AddHit(new_trip.GetLHit());
const int new_ista = ista + new_trip.GetMSta() - new_trip.GetLSta();
new_tr[ista].SetChi2(new_chi2);
- CAFindTrack(new_ista, best_tr, &new_trip, new_tr[ista], min_best_l, new_tr, wData);
+ CAFindTrack(new_ista, best_tr, &new_trip, new_tr[ista], min_best_l, new_tr, wData, vTriplets);
} // add triplet to track
} // for neighbours
} // level = 0
diff --git a/algo/ca/core/tracking/CaTrackFinderWindow.h b/algo/ca/core/tracking/CaTrackFinderWindow.h
index 8983a527feda0851b203522f9f38aad5969490c3..eb7296cf63ec46abf92e6ef89c40c2c140a78947 100644
--- a/algo/ca/core/tracking/CaTrackFinderWindow.h
+++ b/algo/ca/core/tracking/CaTrackFinderWindow.h
@@ -35,7 +35,7 @@ namespace cbm::algo::ca
TrackFinderWindow(const ca::Parameters<fvec>& pars, const fscal mass, const ca::TrackingMode& mode,
ca::TrackingMonitorData& monitorData);
/// Destructor
- ~TrackFinderWindow();
+ ~TrackFinderWindow() = default;
/// Copy constructor
TrackFinderWindow(const TrackFinderWindow&) = default;
@@ -58,10 +58,35 @@ namespace cbm::algo::ca
///-------------------------------
/// Private methods
+ typedef std::array<Vector<ca::Triplet>, constants::size::MaxNstations> TripletArray_t;
+
+ // TripletData_t:
+ // 1 -> vHitFirstTriplet /// link hit -> first triplet { hit, *, *}
+ // 2 -> vHitNofTriplets /// link hit ->n triplets { hit, *, *}
+ // 3 -> vTriplets;
+ typedef std::tuple<Vector<int>, Vector<int>, TripletArray_t> TripletData_t;
+
+ void ConstructTriplets(WindowData& wData);
+
void ReadWindowData(const Vector<Hit>& hits, WindowData& wData);
+ void PrepareGrid(const Vector<Hit>& hits, WindowData& wData);
+
+ void PrepareCAIteration(const ca::Iteration& caIteration, WindowData& wData, const bool isFirst);
+
+ void CreateTracks(WindowData& wData, const ca::Iteration& caIteration, TripletArray_t& vTriplets);
+
+ void CreateTrackCandidates(WindowData& wData, TripletArray_t& vTriplets, const int min_level,
+ const int firstTripletLevel);
+
+ void DoCompetitionLoop(const WindowData& wData);
+
+ void SelectTracks(WindowData& wData);
+
+ void SearchNeighbors(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);
+ unsigned char min_best_l, ca::Branch* new_tr, WindowData& wData, TripletArray_t& vTriplets);
bool checkTripletMatch(const ca::Triplet& l, const ca::Triplet& r, fscal& dchi2, WindowData& wData) const;
@@ -86,16 +111,6 @@ namespace cbm::algo::ca
///-------------------------------
/// Data members
- /// \brief Created triplets vs station index
- std::array<Vector<ca::Triplet>, constants::size::MaxNstations> fvTriplets;
-
- /// \brief Track candidates created out of adjacent triplets before the final track selection.
- /// \note The candidates may share any amount of hits.
- Vector<ca::Branch> fvTrackCandidates{"TrackFinderWindow::fTrackCandidates"};
-
- /// \note Global array for a given thread
- Vector<int> fvHitKeyToTrack{"TrackFinderWindow::fvHitKeyToTrack"};
-
static constexpr bool fDebug = false; // print debug info
const Parameters<fvec>& fParameters; ///< Object of Framework parameters class
@@ -106,6 +121,23 @@ namespace cbm::algo::ca
TrackExtender fTrackExtender; ///< Object of the track extender algorithm
CloneMerger fCloneMerger; ///< Object of the clone merger algorithm
TrackFitter fTrackFitter; ///< Object of the track extender algorithm
+
+ /// \note Global array for a given thread
+ Vector<int> fvHitKeyToTrack{"TrackFinderWindow::fvHitKeyToTrack"};
+
+ // Persistent to avoid memory allocations.
+ // Only used in CreateTrackCandidates().
+ ca::Branch fNewTr[constants::size::MaxNstations];
+
+ // Track candidates created out of adjacent triplets before the final track selection.
+ // The candidates may share any amount of hits.
+ Vector<ca::Branch> fvTrackCandidates{"TrackFinderWindow::fvTrackCandidates"};
+
+ // Triplets and related meta data
+ TripletData_t fTripletData;
+
+ // Triplet temporary storage. Only used in ConstructTriplets().
+ Vector<ca::Triplet> fvTriplets;
};
// ********************************************
diff --git a/algo/ca/core/tracking/CaTrackFit.h b/algo/ca/core/tracking/CaTrackFit.h
index 3aea397d7ffbea2a56f9ab8bf7346899eb3d4c21..dfff3a899c45bda34770d41a064db07b9cf059b8 100644
--- a/algo/ca/core/tracking/CaTrackFit.h
+++ b/algo/ca/core/tracking/CaTrackFit.h
@@ -48,6 +48,8 @@ namespace cbm::algo::ca
TrackFit(const TrackParamBase<DataT>& t) { SetTrack(t); }
+ TrackFit(const DataTmask& m, bool fitV) : fMask(m), fDoFitVelocity(fitV) {}
+
template<typename T>
TrackFit(const TrackParamBase<T>& t)
{
diff --git a/algo/ca/core/tracking/CaTripletConstructor.cxx b/algo/ca/core/tracking/CaTripletConstructor.cxx
index 8fdc40b109bd0c86ae0ebe9fd4b6d75b38b11f9e..63d42c05323e5ec9404f819b84b4bbf08e61cd17 100644
--- a/algo/ca/core/tracking/CaTripletConstructor.cxx
+++ b/algo/ca/core/tracking/CaTripletConstructor.cxx
@@ -19,10 +19,9 @@
using namespace cbm::algo::ca;
-TripletConstructor::TripletConstructor(const ca::Parameters<fvec>& pars, const ca::InputData& input, WindowData& wData,
- const fscal mass, const ca::TrackingMode& mode)
+TripletConstructor::TripletConstructor(const ca::Parameters<fvec>& pars, WindowData& wData, const fscal mass,
+ const ca::TrackingMode& mode)
: fParameters(pars)
- , fInputData(input)
, frWData(wData)
, fDefaultMass(mass)
, fTrackingMode(mode)
@@ -37,9 +36,6 @@ 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);
}
@@ -87,64 +83,60 @@ bool TripletConstructor::InitStations(int istal, int istam, int istar)
}
-void TripletConstructor::CreateTripletsForHit(int istal, int istam, int istar, ca::HitIndex_t iHitL)
+void TripletConstructor::CreateTripletsForHit(Vector<ca::Triplet>& tripletsOut, int istal, int istam, int istar,
+ ca::HitIndex_t iHitL)
{
if (!InitStations(istal, istam, istar)) {
+ tripletsOut.clear();
return;
}
- if (frWData.CurrentIteration()->GetElectronFlag()) {
- fFit.SetParticleMass(constants::phys::ElectronMass);
- }
- else {
- fFit.SetParticleMass(fDefaultMass);
- }
+ ca::TrackFit<fvec> fit(fmask::One(), true);
+ fit.SetParticleMass(frWData.CurrentIteration()->GetElectronFlag() ? constants::phys::ElectronMass : fDefaultMass);
+ fit.SetQp0(frWData.CurrentIteration()->GetMaxQp());
fIhitL = iHitL;
const auto& hitl = frWData.Hit(fIhitL);
// fit a straight line through the target and the left hit.
+ TrackParamV& T = fit.Tr();
{
-
/// Get the field approximation. Add the target to parameters estimation.
/// Propagaete to the middle station of a triplet.
-
- fFit.SetQp0(frWData.CurrentIteration()->GetMaxQp());
-
- ca::FieldValue<fvec> lB, mB, cB, rB _fvecalignment;
- ca::FieldValue<fvec> l_B_global, centB_global _fvecalignment;
+ //ca::FieldValue<fvec> lB, mB, cB, rB _fvecalignment; currently not used
+ //ca::FieldValue<fvec> l_B_global, centB_global _fvecalignment; currently not used
// get the magnetic field along the track.
// (suppose track is straight line with origin in the target)
- 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() = hitl.T();
- T.Vi() = constants::phys::SpeedOfLightInv;
-
- const fvec maxSlopePV = frWData.CurrentIteration()->GetMaxSlopePV();
- const fvec maxQp = frWData.CurrentIteration()->GetMaxQp();
- const fvec txErr2 = maxSlopePV * maxSlopePV / fvec(9.);
- const fvec qpErr2 = maxQp * maxQp / fvec(9.);
-
- T.ResetErrors(1., 1., txErr2, txErr2, qpErr2, (fStaL->timeInfo ? hitl.dT2() : 1.e6), 1.e10);
-
- T.InitVelocityRange(1. / frWData.CurrentIteration()->GetMaxQp());
-
- T.C00() = hitl.dX2();
- T.C10() = hitl.dXY();
- T.C11() = hitl.dY2();
-
- T.Ndf() = (frWData.CurrentIteration()->GetPrimaryFlag()) ? fvec(2.) : fvec(0.);
- T.NdfTime() = (fStaL->timeInfo ? 0 : -1);
+ {
+ 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() = hitl.T();
+ T.Vi() = constants::phys::SpeedOfLightInv;
+
+ const fvec maxSlopePV = frWData.CurrentIteration()->GetMaxSlopePV();
+ const fvec maxQp = frWData.CurrentIteration()->GetMaxQp();
+ const fvec txErr2 = maxSlopePV * maxSlopePV / fvec(9.);
+ const fvec qpErr2 = maxQp * maxQp / fvec(9.);
+
+ T.ResetErrors(1., 1., txErr2, txErr2, qpErr2, (fStaL->timeInfo ? hitl.dT2() : 1.e6), 1.e10);
+ T.InitVelocityRange(1. / frWData.CurrentIteration()->GetMaxQp());
+
+ T.C00() = hitl.dX2();
+ T.C10() = hitl.dXY();
+ T.C11() = hitl.dY2();
+ T.Ndf() = (frWData.CurrentIteration()->GetPrimaryFlag()) ? fvec(2.) : fvec(0.);
+ 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
+ // 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
// field made by the left hit, the target and the station istac in-between.
// is used for extrapolation to the target and to the middle hit
@@ -164,66 +156,66 @@ void TripletConstructor::CreateTripletsForHit(int istal, int istam, int istar, c
}
// add the target constraint
-
- fFit.FilterWithTargetAtLine(fParameters.GetTargetPositionZ(), frWData.TargetMeasurement(), fld0);
-
- fFit.MultipleScattering(fParameters.GetMaterialThickness(fIstaL, T.GetX(), T.GetY()));
+ fit.FilterWithTargetAtLine(fParameters.GetTargetPositionZ(), frWData.TargetMeasurement(), fld0);
+ fit.MultipleScattering(fParameters.GetMaterialThickness(fIstaL, T.GetX(), T.GetY()));
// extrapolate to the middle hit
-
- fFit.ExtrapolateLine(fStaM->fZ, fFldL);
-
- fTrL = T;
+ fit.ExtrapolateLine(fStaM->fZ, fFldL);
}
- FindDoublets();
- FitDoublets();
- FindTriplets();
-}
-
-
-void TripletConstructor::FindDoublets()
-{
/// Find the doublets. Reformat data into portions of doublets.
- int iMC = -1;
- if (fParameters.DevIsMatchDoubletsViaMc()) {
- iMC = ca::Framework::GetMcTrackIdForWindowHit(fIhitL);
- if (iMC < 0) {
+ auto FindDoubletHits = [&]() {
+ const bool matchMc = fParameters.DevIsMatchDoubletsViaMc();
+ const int iMC = matchMc ? ca::Framework::GetMcTrackIdForWindowHit(fIhitL) : -1;
+ fDoubletData.second.clear();
+ if (iMC < 0 && matchMc) {
return;
}
+ CollectHits(fDoubletData.second, fit, fIstaM, frWData.CurrentIteration()->GetDoubletChi2Cut(), iMC,
+ fParameters.GetMaxDoubletsPerSinglet());
+ };
+
+ FindDoubletHits();
+ FindDoublets(fit);
+
+ //D.Smith 28.8.24: Moving this upward (before doublet finding) changes QA output slightly
+ if (fIstaR >= fParameters.GetNstationsActive()) {
+ tripletsOut.clear();
+ return;
}
- CollectHits(fTrL, fIstaM, frWData.CurrentIteration()->GetDoubletChi2Cut(), iMC, fHitsM_2,
- fParameters.GetMaxDoubletsPerSinglet());
+ FindTripletHits();
+ FindTriplets();
+ SelectTriplets(tripletsOut);
}
-void TripletConstructor::FitDoublets()
+void TripletConstructor::FindDoublets(ca::TrackFit<fvec>& fit)
{
-
// ---- Add the middle hits to parameters estimation ----
+ Vector<TrackParamV>& tracks = fDoubletData.first;
+ Vector<ca::HitIndex_t>& hitsM = fDoubletData.second;
- Vector<ca::HitIndex_t> hitsMtmp("TripletConstructor::hitsMtmp", fHitsM_2);
-
- fHitsM_2.clear();
- fTracks_2.clear();
- fTracks_2.reserve(hitsMtmp.size());
+ tracks.clear();
+ tracks.reserve(hitsM.size());
const bool isMomentumFitted = (fIsTargetField || (fStaL->fieldStatus != 0) || (fStaM->fieldStatus != 0));
- double maxQp = frWData.CurrentIteration()->GetMaxQp();
- for (unsigned int i2 = 0; i2 < hitsMtmp.size(); i2++) {
+ const TrackParamV Tr = fit.Tr(); // copy contents of fit
- const ca::HitIndex_t indM = hitsMtmp[i2];
+ auto it2 = hitsM.begin();
+ for (auto it = hitsM.begin(); it != hitsM.end(); it++) {
+
+ const ca::HitIndex_t indM = *it;
const ca::Hit& hitm = frWData.Hit(indM);
if (frWData.IsHitSuppressed(indM)) {
continue;
}
- TrackParamV& T2 = fFit.Tr();
- T2 = fTrL;
- fFit.SetQp0(fvec(0.f));
+ TrackParamV& T2 = fit.Tr();
+ T2 = Tr;
+ fit.SetQp0(fvec(0.f));
fvec z_2 = hitm.Z();
ca::MeasurementXy<fvec> m_2(hitm.X(), hitm.Y(), hitm.dX2(), hitm.dY2(), hitm.dXY(), fvec::One(), fvec::One());
@@ -231,13 +223,12 @@ void TripletConstructor::FitDoublets()
fvec dt2_2 = hitm.dT2();
// add the middle hit
-
- fFit.ExtrapolateLineNoField(z_2);
- fFit.FilterXY(m_2);
- fFit.FilterTime(t_2, dt2_2, fmask(fStaM->timeInfo));
- fFit.SetQp0(isMomentumFitted ? fFit.Tr().GetQp() : maxQp);
- fFit.MultipleScattering(fParameters.GetMaterialThickness(fIstaM, T2.GetX(), T2.Y()));
- fFit.SetQp0(fFit.Tr().Qp());
+ fit.ExtrapolateLineNoField(z_2);
+ fit.FilterXY(m_2);
+ fit.FilterTime(t_2, dt2_2, fmask(fStaM->timeInfo));
+ fit.SetQp0(isMomentumFitted ? fit.Tr().GetQp() : frWData.CurrentIteration()->GetMaxQp());
+ fit.MultipleScattering(fParameters.GetMaterialThickness(fIstaM, T2.GetX(), T2.Y()));
+ fit.SetQp0(fit.Tr().Qp());
// check if there are other hits close to the doublet on the same station
if (ca::kMcbm != fTrackingMode) {
@@ -247,9 +238,9 @@ void TripletConstructor::FitDoublets()
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++) {
+ for (auto itClone = it + 1; itClone != hitsM.end(); itClone++) {
- const int indClone = hitsMtmp[iClone];
+ const int indClone = *itClone;
const ca::Hit& hitClone = frWData.Hit(indClone);
const fscal dz = hitClone.Z() - T2.Z()[0];
@@ -283,61 +274,56 @@ void TripletConstructor::FitDoublets()
}
}
- fHitsM_2.push_back(indM);
- fTracks_2.push_back(T2);
- } // i2
+ tracks.push_back(T2);
+ *it2 = indM;
+ it2++;
+ } // it
+ hitsM.erase(it2, hitsM.end());
}
-void TripletConstructor::FindTriplets()
+void TripletConstructor::FindTripletHits()
{
- if (fIstaR >= fParameters.GetNstationsActive()) {
- return;
- }
- FindRightHit();
- FitTriplets();
- StoreTriplets();
-}
-
+ //auto& [tracks_2, hitsM_2] = doublets; TO DO: Reactivate when MacOS compiler bug is fixed.
+ Vector<TrackParamV>& tracks_2 = fDoubletData.first;
+ Vector<ca::HitIndex_t>& hitsM_2 = fDoubletData.second;
-void TripletConstructor::FindRightHit()
-{
+ Vector<ca::HitIndex_t>& hitsM_3 = std::get<1>(fTripletData);
+ Vector<ca::HitIndex_t>& hitsR_3 = std::get<2>(fTripletData);
/// Add the middle hits to parameters estimation. Propagate to right station.
/// Find the triplets(right hit). Reformat data in the portion of triplets.
-
- fFit.SetQp0(fvec(0.));
-
- if (fIstaR >= fParameters.GetNstationsActive()) return;
+ ca::TrackFit<fvec> fit(fmask::One(), true);
+ fit.SetParticleMass(frWData.CurrentIteration()->GetElectronFlag() ? constants::phys::ElectronMass : fDefaultMass);
+ fit.SetQp0(fvec(0.));
{
- int maxTriplets = fHitsM_2.size() * fParameters.GetMaxTripletPerDoublets();
- fHitsM_3.clear();
- fHitsR_3.clear();
- fHitsM_3.reserve(maxTriplets);
- fHitsR_3.reserve(maxTriplets);
+ const int maxTriplets = hitsM_2.size() * fParameters.GetMaxTripletPerDoublets();
+ hitsM_3.clear();
+ hitsR_3.clear();
+ hitsM_3.reserve(maxTriplets);
+ hitsR_3.reserve(maxTriplets);
}
// ---- Add the middle hits to parameters estimation. Propagate to right station. ----
const double maxSlope = frWData.CurrentIteration()->GetMaxSlope();
const double tripletChi2Cut = frWData.CurrentIteration()->GetTripletChi2Cut();
- for (unsigned int i2 = 0; i2 < fHitsM_2.size(); i2++) {
+ for (unsigned int i2 = 0; i2 < hitsM_2.size(); i2++) {
- fFit.SetTrack(fTracks_2[i2]);
- TrackParamV& T2 = fFit.Tr();
+ fit.SetTrack(tracks_2[i2]);
+ TrackParamV& T2 = fit.Tr();
// extrapolate to the right hit station
-
- fFit.Extrapolate(fStaR->fZ, fFldL);
+ fit.Extrapolate(fStaR->fZ, fFldL);
if constexpr (fDebugDublets) {
- ca::HitIndex_t iwhit[2] = {fIhitL, fHitsM_2[i2]};
+ ca::HitIndex_t iwhit[2] = {fIhitL, hitsM_2[i2]};
ca::HitIndex_t ihit[2] = {frWData.Hit(iwhit[0]).Id(), frWData.Hit(iwhit[1]).Id()};
const int ih0 = ihit[0];
const int ih1 = ihit[1];
- const ca::Hit& h0 = fInputData.GetHit(ih0);
- const ca::Hit& h1 = fInputData.GetHit(ih1);
+ const ca::Hit& h0 = frWData.Hit(iwhit[0]);
+ const ca::Hit& h1 = frWData.Hit(iwhit[1]);
LOG(info) << "\n====== Extrapolated Doublet : "
<< " iter " << frWData.CurrentIteration()->GetName() << " hits: {" << fIstaL << "/" << ih0 << " "
@@ -364,7 +350,7 @@ void TripletConstructor::FindRightHit()
return true;
}
if (fParameters.DevIsMatchTripletsViaMc()) {
- int indM = fHitsM_2[i2];
+ int indM = hitsM_2[i2];
iMC = ca::Framework::GetMcTrackIdForWindowHit(fIhitL);
if (iMC < 0 || iMC != ca::Framework::GetMcTrackIdForWindowHit(indM)) {
return true;
@@ -389,7 +375,7 @@ void TripletConstructor::FindRightHit()
}
Vector<ca::HitIndex_t> collectedHits;
- CollectHits(T2, fIstaR, tripletChi2Cut, iMC, collectedHits, fParameters.GetMaxTripletPerDoublets());
+ CollectHits(collectedHits, fit, fIstaR, tripletChi2Cut, iMC, fParameters.GetMaxTripletPerDoublets());
if (collectedHits.size() >= fParameters.GetMaxTripletPerDoublets()) {
// FU, 28.08.2024, Comment the following log lines since it spams the output
@@ -404,8 +390,8 @@ void TripletConstructor::FindRightHit()
}
for (ca::HitIndex_t& irh : collectedHits) {
if constexpr (fDebugDublets) {
- ca::HitIndex_t ihit = frWData.Hit(irh).Id();
- const ca::Hit& h = fInputData.GetHit(ihit);
+ const ca::Hit& h = frWData.Hit(irh);
+ ca::HitIndex_t ihit = h.Id();
LOG(info) << " hit " << ihit << " " << h.ToString();
}
if (frWData.IsHitSuppressed(irh)) {
@@ -415,8 +401,8 @@ void TripletConstructor::FindRightHit()
continue;
}
// pack the triplet
- fHitsM_3.push_back(fHitsM_2[i2]);
- fHitsR_3.push_back(irh);
+ hitsM_3.push_back(hitsM_2[i2]);
+ hitsR_3.push_back(irh);
} // search area
if constexpr (fDebugDublets) {
LOG(info) << "======== end of extrapolated doublet ==== \n";
@@ -424,16 +410,17 @@ void TripletConstructor::FindRightHit()
} // i2
}
-
-void TripletConstructor::FitTriplets()
+void TripletConstructor::FindTriplets()
{
constexpr int nIterations = 2;
- const Tindex n3 = fHitsM_3.size();
- assert(fHitsM_3.size() == fHitsR_3.size());
+ Vector<TrackParamV>& tracks = std::get<0>(fTripletData);
+ Vector<ca::HitIndex_t>& hitsM = std::get<1>(fTripletData);
+ Vector<ca::HitIndex_t>& hitsR = std::get<2>(fTripletData);
+ assert(hitsM.size() == hitsR.size());
- fTracks_3.clear();
- fTracks_3.reset(n3, TrackParamV());
+ tracks.clear();
+ tracks.reserve(hitsM.size());
/// Refit Triplets
if constexpr (fDebugTriplets) {
@@ -443,13 +430,7 @@ void TripletConstructor::FitTriplets()
ca::TrackFit<fvec> fit;
fit.SetMask(fmask::One());
-
- if (frWData.CurrentIteration()->GetElectronFlag()) {
- fit.SetParticleMass(constants::phys::ElectronMass);
- }
- else {
- fit.SetParticleMass(fDefaultMass);
- }
+ fit.SetParticleMass(frWData.CurrentIteration()->GetElectronFlag() ? constants::phys::ElectronMass : fDefaultMass);
// prepare data
const int NHits = 3; // triplets
@@ -461,12 +442,10 @@ void TripletConstructor::FitTriplets()
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) {
+ for (size_t i3 = 0; i3 < hitsM.size(); ++i3) {
// prepare data
- const ca::HitIndex_t iwhit[NHits] = {fIhitL, fHitsM_3[i3], fHitsR_3[i3]};
+ const ca::HitIndex_t iwhit[NHits] = {fIhitL, hitsM[i3], hitsR[i3]};
const ca::HitIndex_t ihit[NHits] = {frWData.Hit(iwhit[0]).Id(), frWData.Hit(iwhit[1]).Id(),
frWData.Hit(iwhit[2]).Id()};
@@ -476,6 +455,9 @@ void TripletConstructor::FitTriplets()
int mc2 = ca::Framework::GetMcTrackIdForCaHit(ihit[1]);
int mc3 = ca::Framework::GetMcTrackIdForCaHit(ihit[2]);
if ((mc1 != mc2) || (mc1 != mc3)) {
+ // D.S.: Added to preserve the ordering when switching from index-based
+ // access to push_back(). Discuss with SZ.
+ tracks.push_back(TrackParamV());
continue;
}
}
@@ -485,7 +467,7 @@ void TripletConstructor::FitTriplets()
ca::MeasurementXy<fvec> mxy[NHits];
for (int ih = 0; ih < NHits; ++ih) {
- const ca::Hit& hit = fInputData.GetHit(ihit[ih]);
+ const ca::Hit& hit = frWData.Hit(iwhit[ih]);
mxy[ih] = ca::MeasurementXy<fvec>(hit.X(), hit.Y(), hit.dX2(), hit.dY2(), hit.dXY(), fvec::One(), fvec::One());
x[ih] = hit.X();
@@ -503,6 +485,7 @@ void TripletConstructor::FitTriplets()
fvec tx[3] = {(x[1] - x[0]) / (z[1] - z[0]), (x[2] - x[0]) / (z[2] - z[0]), (x[2] - x[1]) / (z[2] - z[1])};
fvec ty[3] = {(y[1] - y[0]) / (z[1] - z[0]), (y[2] - y[0]) / (z[2] - z[0]), (y[2] - y[1]) / (z[2] - z[1])};
+
for (int ih = 0; ih < NHits; ++ih) {
fvec dz = (sta[ih].fZ - z[ih]);
B[ih] = sta[ih].fieldSlice.GetFieldValue(x[ih] + tx[ih] * dz, y[ih] + ty[ih] * dz);
@@ -526,6 +509,7 @@ void TripletConstructor::FitTriplets()
for (int iiter = 0; iiter < nIterations; ++iiter) {
auto fitTrack = [&](int startIdx, int endIdx, int step, FitDirection direction) {
+ const fvec maxQp = frWData.CurrentIteration()->GetMaxQp();
fit.SetQp0(T.Qp());
fit.Qp0()(fit.Qp0() > maxQp) = maxQp;
fit.Qp0()(fit.Qp0() < -maxQp) = -maxQp;
@@ -569,7 +553,7 @@ void TripletConstructor::FitTriplets()
fitTrack(NHits - 1, -1, -1, FitDirection::kUpstream);
} // for iiter
- fTracks_3[i3] = T;
+ tracks.push_back(T);
if constexpr (fDebugTriplets) {
int ih0 = ihit[0];
@@ -580,9 +564,9 @@ void TripletConstructor::FitTriplets()
int mc3 = ca::Framework::GetMcTrackIdForCaHit(ih2);
if (1 || (mc1 >= 0) && (mc1 == mc2) && (mc1 == mc3)) {
- const ca::Hit& h0 = fInputData.GetHit(ih0);
- const ca::Hit& h1 = fInputData.GetHit(ih1);
- const ca::Hit& h2 = fInputData.GetHit(ih2);
+ const ca::Hit& h0 = frWData.Hit(iwhit[0]);
+ const ca::Hit& h1 = frWData.Hit(iwhit[1]);
+ const ca::Hit& h2 = frWData.Hit(iwhit[2]);
//const CbmL1MCTrack& mctr = CbmL1::Instance()->GetMcTracks()[mc1];
LOG(info) << "== fitted triplet: "
<< " iter " << frWData.CurrentIteration()->GetName() << " hits: {" << fIstaL << "/" << ih0 << " "
@@ -599,41 +583,41 @@ void TripletConstructor::FitTriplets()
}
}
} //i3
-
-} // FitTriplets
+} // FindTriplets
-void TripletConstructor::StoreTriplets()
+void TripletConstructor::SelectTriplets(Vector<ca::Triplet>& tripletsOut)
{
- /// Selects good triplets and saves them into fTriplets.
+ /// Selects good triplets and saves them into tripletsOut.
/// Finds neighbouring triplets at the next station.
- const Tindex n3 = fHitsM_3.size();
+ Vector<TrackParamV>& tracks = std::get<0>(fTripletData);
+ Vector<ca::HitIndex_t>& hitsM = std::get<1>(fTripletData);
+ Vector<ca::HitIndex_t>& hitsR = std::get<2>(fTripletData);
bool isMomentumFitted = ((fStaL->fieldStatus != 0) || (fStaM->fieldStatus != 0) || (fStaR->fieldStatus != 0));
- fTriplets.clear();
- fTriplets.reserve(n3);
+ tripletsOut.clear();
+ tripletsOut.reserve(hitsM.size());
- const fscal tripletFinalChi2Cut = frWData.CurrentIteration()->GetTripletFinalChi2Cut();
- for (Tindex i3 = 0; i3 < n3; ++i3) {
+ for (size_t i3 = 0; i3 < hitsM.size(); ++i3) {
- TrackParamV& T3 = fTracks_3[i3];
+ TrackParamV& T3 = tracks[i3];
// TODO: SG: normalize chi2, separate cuts on time and space
const fscal chi2 = T3.GetChiSq()[0] + T3.GetChiSqTime()[0];
const ca::HitIndex_t ihitl = fIhitL;
- const ca::HitIndex_t ihitm = fHitsM_3[i3];
- const ca::HitIndex_t ihitr = fHitsR_3[i3];
+ const ca::HitIndex_t ihitm = hitsM[i3];
+ const ca::HitIndex_t ihitr = hitsR[i3];
CBMCA_DEBUG_ASSERT(ihitl < frWData.HitStartIndexOnStation(fIstaL + 1));
CBMCA_DEBUG_ASSERT(ihitm < frWData.HitStartIndexOnStation(fIstaM + 1));
CBMCA_DEBUG_ASSERT(ihitr < frWData.HitStartIndexOnStation(fIstaR + 1));
if (!frWData.CurrentIteration()->GetTrackFromTripletsFlag()) {
- if (chi2 > tripletFinalChi2Cut) {
+ if (chi2 > frWData.CurrentIteration()->GetTripletFinalChi2Cut()) {
continue;
}
}
@@ -651,24 +635,22 @@ void TripletConstructor::StoreTriplets()
// removing it leads to a higher ghost ratio
Cqp += 0.001;
- fTriplets.emplace_back(ihitl, ihitm, ihitr, fIstaL, fIstaM, fIstaR, 0, 0, 0, chi2, qp, Cqp, T3.Tx()[0], T3.C22()[0],
- T3.Ty()[0], T3.C33()[0], isMomentumFitted);
+ tripletsOut.emplace_back(ihitl, ihitm, ihitr, fIstaL, fIstaM, fIstaR, 0, 0, 0, chi2, qp, Cqp, T3.Tx()[0],
+ T3.C22()[0], T3.Ty()[0], T3.C33()[0], isMomentumFitted);
}
}
-void TripletConstructor::CollectHits(const TrackParamV& Tr, const int iSta, const double chi2Cut, const int iMC,
- Vector<ca::HitIndex_t>& collectedHits, int maxNhits)
+void TripletConstructor::CollectHits(Vector<ca::HitIndex_t>& collectedHits, ca::TrackFit<fvec>& fit, const int iSta,
+ const double chi2Cut, const int iMC, const int maxNhits)
{
/// Collect hits on a station
-
collectedHits.clear();
collectedHits.reserve(maxNhits);
const ca::Station<fvec>& sta = fParameters.GetStation(iSta);
- fFit.SetTrack(Tr);
- TrackParamV& T = fFit.Tr();
+ TrackParamV& T = fit.Tr();
//LOG(info) << T.chi2[0] ;
T.ChiSq() = 0.;
@@ -730,7 +712,7 @@ void TripletConstructor::CollectHits(const TrackParamV& Tr, const int iSta, cons
const fscal z = hit.Z();
// check y-boundaries
- const auto [y, C11] = fFit.ExtrapolateLineYdY2(z);
+ const auto [y, C11] = fit.ExtrapolateLineYdY2(z);
const fscal dy_est = sqrt(Pick_m22[0] * C11[0]) + hit.RangeY();
const fscal dY = hit.Y() - y[0];
@@ -742,7 +724,7 @@ void TripletConstructor::CollectHits(const TrackParamV& Tr, const int iSta, cons
}
// check x-boundaries
- const auto [x, C00] = fFit.ExtrapolateLineXdX2(z);
+ const auto [x, C00] = fit.ExtrapolateLineXdX2(z);
const fscal dx_est = sqrt(Pick_m22[0] * C00[0]) + hit.RangeX();
const fscal dX = hit.X() - x[0];
@@ -756,7 +738,7 @@ void TripletConstructor::CollectHits(const TrackParamV& Tr, const int iSta, cons
// check chi2
ca::MeasurementXy<fvec> mxy(hit.X(), hit.Y(), hit.dX2(), hit.dY2(), hit.dXY(), fvec::One(), fvec::One());
- const fvec C10 = fFit.ExtrapolateLineDxy(z);
+ const fvec C10 = fit.ExtrapolateLineDxy(z);
const auto [chi2x, chi2u] = ca::TrackFit<fvec>::GetChi2XChi2U(mxy, x, y, C00, C10, C11);
// TODO: adjust the cut, cut on chi2x & chi2u separately
diff --git a/algo/ca/core/tracking/CaTripletConstructor.h b/algo/ca/core/tracking/CaTripletConstructor.h
index af9bed635c144b58073a61ee1204f2eb0cd8237f..6282b0da73cec42920508f154000259ab532c2b0 100644
--- a/algo/ca/core/tracking/CaTripletConstructor.h
+++ b/algo/ca/core/tracking/CaTripletConstructor.h
@@ -30,8 +30,8 @@ namespace cbm::algo::ca
/// Constructor
/// \param nThreads Number of threads for multi-threaded mode
- TripletConstructor(const ca::Parameters<fvec>& pars, const ca::InputData& input, WindowData& wData,
- const fscal mass, const ca::TrackingMode& mode);
+ TripletConstructor(const ca::Parameters<fvec>& pars, WindowData& wData, const fscal mass,
+ const ca::TrackingMode& mode);
/// Copy constructor
TripletConstructor(const TripletConstructor&) = delete;
@@ -49,37 +49,32 @@ namespace cbm::algo::ca
~TripletConstructor() = default;
/// ------ FUNCTIONAL PART ------
-
- void CreateTripletsForHit(int istal, int istam, int istar, ca::HitIndex_t ihl);
-
- const Vector<ca::Triplet>& GetTriplets() const { return fTriplets; }
+ void CreateTripletsForHit(Vector<ca::Triplet>& tripletsOut, int istal, int istam, int istar, ca::HitIndex_t ihl);
private:
+ typedef std::pair<Vector<TrackParamV>, Vector<ca::HitIndex_t>> Doublet_t;
+ typedef std::tuple<Vector<TrackParamV>, Vector<ca::HitIndex_t>, Vector<ca::HitIndex_t>> Triplet_t;
+
bool InitStations(int istal, int istam, int istar);
/// Find the doublets. Reformat data in the portion of doublets.
- void FindDoublets();
- void FitDoublets();
-
- /// Find triplets on station
- void FindTriplets();
+ void FindDoublets(ca::TrackFit<fvec>& fit);
/// Add the middle hits to parameters estimation. Propagate to right station.
/// Find the triplets (right hit). Reformat data in the portion of triplets.
- void FindRightHit();
+ void FindTripletHits();
- /// Fit Triplets
- void FitTriplets();
+ /// Find triplets on station
+ void FindTriplets();
- /// Select triplets. Save them into vTriplets.
- void StoreTriplets();
+ /// Select good triplets.
+ void SelectTriplets(Vector<ca::Triplet>& tripletsOut);
- void CollectHits(const TrackParamV& Tr, const int iSta, const double chi2Cut, const int iMC,
- Vector<ca::HitIndex_t>& collectedHits, int maxNhits);
+ void CollectHits(Vector<ca::HitIndex_t>& collectedHits, ca::TrackFit<fvec>& fit, const int iSta,
+ const double chi2Cut, const int iMC, const int maxNhits);
private:
const Parameters<fvec>& fParameters; ///< Object of Framework parameters class
- const InputData& fInputData; ///< Tracking input data
WindowData& frWData;
fscal fDefaultMass{constants::phys::MuonMass}; ///< mass of the propagated particle [GeV/c2]
ca::TrackingMode fTrackingMode;
@@ -103,19 +98,13 @@ namespace cbm::algo::ca
fFld1Sta[3]; // three stations for approximating the field between the left and the right hit
ca::HitIndex_t fIhitL; ///< index of the left hit in fAlgo->fWindowHits
- TrackParamV fTrL;
ca::FieldRegion<fvec> fFldL;
- Vector<ca::HitIndex_t> fHitsM_2{"TripletConstructor::fHitsM_2"};
- Vector<TrackParamV> fTracks_2{"TripletConstructor::fTracks_2"};
-
- Vector<TrackParamV> fTracks_3{"TripletConstructor::fTracks_3"};
- Vector<ca::HitIndex_t> fHitsM_3{"TripletConstructor::fHitsM_3"};
- Vector<ca::HitIndex_t> fHitsR_3{"TripletConstructor::fHitsR_3"};
-
- Vector<ca::Triplet> fTriplets{"TripletConstructor::fTriplets"};
+ // Persistent storage for triplet tracks and hits
+ Triplet_t fTripletData;
- ca::TrackFit<fvec> fFit;
+ // Persistent storage for doublet tracks and hits
+ Doublet_t fDoubletData;
private:
static constexpr bool fDebugDublets = false; // print debug info for dublets
diff --git a/algo/ca/core/utils/CaVector.h b/algo/ca/core/utils/CaVector.h
index 16fee44c40f0ccb15a052db5a5a6bbd557b99b82..aabed81edf283f68c2e159fca1f97827c7d01c2d 100644
--- a/algo/ca/core/utils/CaVector.h
+++ b/algo/ca/core/utils/CaVector.h
@@ -265,6 +265,7 @@ namespace cbm::algo::ca
using Tbase::cend;
using Tbase::clear;
using Tbase::end;
+ using Tbase::erase;
using Tbase::insert; //TODO:: make it private
using Tbase::pop_back;
using Tbase::rbegin;