diff --git a/algo/ca/core/data/CaWindowData.h b/algo/ca/core/data/CaWindowData.h index dc25c37a30f97c9afe275dae68979ee6198562cf..67269e9533f11712885b8ec3d3ffcf026bf56c0a 100644 --- a/algo/ca/core/data/CaWindowData.h +++ b/algo/ca/core/data/CaWindowData.h @@ -22,7 +22,7 @@ namespace cbm::algo::ca { /// \class WindowData - /// \brief Container for internal data, processed on a single sub-TS window + /// \brief Container for internal data, processed on a single time window class WindowData { public: /// \brief Constructor @@ -74,6 +74,17 @@ namespace cbm::algo::ca return fvHitStartIndexOnStation[iStation]; } + /// \brief Hit key flag: if this hit or cluster was already used + /// \param iKey Index of the key (index of the hit or the cluster) + [[gnu::always_inline]] unsigned char IsHitKeyUsed(HitKeyIndex_t iKey) const { return fvbHitKeyFlags[iKey]; } + + /// \brief Hit key flag: if this hit or cluster was already used + /// \param iKey Index of the key (index of the hit or the cluster) + [[gnu::always_inline]] unsigned char& IsHitKeyUsed(HitKeyIndex_t iKey) { return fvbHitKeyFlags[iKey]; } + + /// \brief Access to the hit key flags container + [[gnu::always_inline]] Vector<unsigned char>& HitKeyFlags() { return fvbHitKeyFlags; } + /// \brief Index of the first hit on the station /// \param iStation Index of the tracking station [[gnu::always_inline]] HitIndex_t HitStartIndexOnStation(int iStation) const @@ -162,8 +173,12 @@ namespace cbm::algo::ca /// hit.Id is replaced by the hit index in fInputData Vector<ca::Hit> fvHits{"WindowData::fHits"}; + /// \brief List of used hit keys + /// \note Global for all the time windows within one thread + Vector<unsigned char> fvbHitKeyFlags{"WindowData::fvbHitKeyFlags"}; + /// \brief Flag, if the hit is suppressed for tracking - Vector<unsigned char> fvbHitSuppressed{"WindowData::fbHitSuppressed"}; + Vector<unsigned char> fvbHitSuppressed{"WindowData::fvbHitSuppressed"}; /// \brief First hit index of the station std::array<HitIndex_t, kMaxNofStations + 1> fvHitStartIndexOnStation = {0}; diff --git a/algo/ca/core/tracking/CaFramework.h b/algo/ca/core/tracking/CaFramework.h index 9d9ec0488b9a6466dfb73b74effda4f4157a3289..74b85031a249962b7dc88a1ea7d273ce42b559c8 100644 --- a/algo/ca/core/tracking/CaFramework.h +++ b/algo/ca/core/tracking/CaFramework.h @@ -18,6 +18,7 @@ #include "CaTimer.h" #include "CaTrack.h" #include "CaTrackExtender.h" +#include "CaTrackExtraInfo.h" #include "CaTrackFinder.h" #include "CaTrackFinderWindow.h" #include "CaTrackFitter.h" @@ -207,6 +208,7 @@ namespace cbm::algo::ca int GetNofThreads() const { return fNofThreads; } + bool IfCollectTrackExtraInfo() const { return fbCollectTrackExtraInfo; } private: int fNstationsBeforePipe{0}; ///< number of stations before pipe (MVD stations in CBM) @@ -219,10 +221,9 @@ namespace cbm::algo::ca // *************************** Parameters<fvec> fParameters; ///< Object of Framework parameters class - InputData fInputData; ///< Tracking input data + InputData fInputData; ///< Tracking input data - Vector<unsigned char> fvHitKeyFlags{ - "Framework::fvHitKeyFlags"}; ///< List of key flags: has been this hit or cluster already used + bool fbCollectTrackExtraInfo = true; ///< Flag, if to fill the additional track info structure TrackingMonitorData fMonitorData{}; ///< Tracking monitor data (statistics per call) @@ -241,6 +242,7 @@ namespace cbm::algo::ca double fCaRecoTime{0.}; // time of the track finder + fitter Vector<Track> fRecoTracks{"Framework::fRecoTracks"}; ///< reconstructed tracks + Vector<TrackExtraInfo> fRecoTrackExtraInfo{"Framework::fRecoTrackExtraInfo"}; ///< Extra info on tracks Vector<ca::HitIndex_t> fRecoHits{"Framework::fRecoHits"}; ///< packed hits of reconstructed tracks fvec EventTime{0.f}; diff --git a/algo/ca/core/tracking/CaTrackExtender.cxx b/algo/ca/core/tracking/CaTrackExtender.cxx index 9dcfddbd2a4a2e9e2559b7aae586dd2201eec5d6..6ece464b2e2b50b9ad7963b9cea32f49c9cb1d95 100644 --- a/algo/ca/core/tracking/CaTrackExtender.cxx +++ b/algo/ca/core/tracking/CaTrackExtender.cxx @@ -233,7 +233,7 @@ void TrackExtender::FindMoreHits(ca::Branch& t, TrackParamV& Tout, const bool up //if (GetFUsed((*fStripFlag)[hit.f] | (*fStripFlag)[hit.b])) continue; // if used - if (frAlgo.fvHitKeyFlags[hit.FrontKey()] || frAlgo.fvHitKeyFlags[hit.BackKey()]) { + if (frWData.IsHitKeyUsed(hit.FrontKey()) || frWData.IsHitKeyUsed(hit.BackKey())) { continue; } diff --git a/algo/ca/core/tracking/CaTrackFinder.cxx b/algo/ca/core/tracking/CaTrackFinder.cxx index 2dfedb8f9ae5cb4e402104574a6f1a5e4ee297c8..96bc482b184e26a3a1fb912d4f68091775a2b1a1 100644 --- a/algo/ca/core/tracking/CaTrackFinder.cxx +++ b/algo/ca/core/tracking/CaTrackFinder.cxx @@ -56,24 +56,27 @@ void TrackFinder::FindTracks() auto timerStart = std::chrono::high_resolution_clock::now(); // ----- Reset data arrays ------------------------------------------------------------------------------------------- - frAlgo.fvHitKeyFlags.reset(frAlgo.fInputData.GetNhitKeys(), 0); - size_t nHitsExpected = 2 * frAlgo.fInputData.GetNhits(); size_t nTracksExpected = 2 * frAlgo.fInputData.GetNhits() / frAlgo.GetParameters().GetNstationsActive(); for (int iThread = 0; iThread < frAlgo.GetNofThreads(); ++iThread) { fvRecoTracks[iThread].clear(); fvRecoTracks[iThread].reserve(nTracksExpected / frAlgo.GetNofThreads()); + if (frAlgo.IfCollectTrackExtraInfo()) { + fvRecoTrackExtraInfo[iThread].clear(); + fvRecoTrackExtraInfo[iThread].reserve(nTracksExpected / frAlgo.GetNofThreads()); + } fvRecoHitIndices[iThread].clear(); fvRecoHitIndices[iThread].reserve(nHitsExpected / frAlgo.GetNofThreads()); frAlgo.fvTrackFinderWindow[iThread].InitTimeslice(); + auto& wData = frAlgo.fvWData[iThread]; + wData.HitKeyFlags().reset(frAlgo.fInputData.GetNhitKeys(), 0); for (int iS = 0; iS < frAlgo.GetParameters().GetNstationsActive(); ++iS) { - frAlgo.fvWData[iThread].TsHitIndices(iS).clear(); - frAlgo.fvWData[iThread].TsHitIndices(iS).reserve(frAlgo.fInputData.GetNhits()); + wData.TsHitIndices(iS).clear(); + wData.TsHitIndices(iS).reserve(frAlgo.fInputData.GetNhits()); } } - frAlgo.fvHitKeyFlags.reset(frAlgo.fInputData.GetNhitKeys(), 0); frAlgo.fHitTimeInfo.reset(frAlgo.fInputData.GetNhits()); @@ -94,6 +97,7 @@ void TrackFinder::FindTracks() const int nDataStreams = frAlgo.fInputData.GetNdataStreams(); + auto& wDataThread0 = frAlgo.fvWData[0]; // NOTE: Thread 0 must be always defined for (int iStream = 0; iStream < nDataStreams; ++iStream) { fscal maxTimeBeforeHit = std::numeric_limits<fscal>::lowest(); @@ -131,8 +135,11 @@ void TrackFinder::FindTracks() } if (info.fEventTimeMin > 500.e6 || info.fEventTimeMax < -500.) { // cut hits with bogus start time > 500 ms - frAlgo.fvHitKeyFlags[h.FrontKey()] = 1; - frAlgo.fvHitKeyFlags[h.BackKey()] = 1; + for (int iThread = 0; iThread < frAlgo.GetNofThreads(); ++iThread) { + auto& wData = frAlgo.fvWData[iThread]; + wData.IsHitKeyUsed(h.FrontKey()) = 1; + wData.IsHitKeyUsed(h.BackKey()) = 1; + } LOG(error) << "CATrackFinder: skip bogus hit " << h.ToString(); continue; } @@ -156,7 +163,7 @@ void TrackFinder::FindTracks() for (int ih = nStreamHits - 1; ih >= 0; --ih) { ca::HitIndex_t caHitId = frAlgo.fInputData.GetStreamStartIndex(iStream) + ih; const ca::Hit& h = frAlgo.fInputData.GetHit(caHitId); - if (frAlgo.fvHitKeyFlags[h.FrontKey()] || frAlgo.fvHitKeyFlags[h.BackKey()]) { + if (wDataThread0.IsHitKeyUsed(h.FrontKey()) || wDataThread0.IsHitKeyUsed(h.BackKey())) { continue; } // the hit is skipped CaHitTimeInfo& info = frAlgo.fHitTimeInfo[caHitId]; @@ -174,7 +181,7 @@ void TrackFinder::FindTracks() for (int ih = 0; (ih < nStreamHits) && (tmp < 10000); ++ih) { ca::HitIndex_t caHitId = frAlgo.fInputData.GetStreamStartIndex(iStream) + ih; const ca::Hit& h = frAlgo.fInputData.GetHit(caHitId); - if (frAlgo.fvHitKeyFlags[h.FrontKey()] || frAlgo.fvHitKeyFlags[h.BackKey()]) { + if (wDataThread0.IsHitKeyUsed(h.FrontKey()) || wDataThread0.IsHitKeyUsed(h.BackKey())) { continue; } // the hit is skipped CaHitTimeInfo& info = frAlgo.fHitTimeInfo[caHitId]; @@ -224,17 +231,19 @@ void TrackFinder::FindTracks() for (int iTh = 0; iTh < frAlgo.GetNofThreads(); ++iTh) { vThreadList.emplace_back(&TrackFinder::FindTracksThread, this, iTh); } - for (auto& th: vThreadList) { - if (th.joinable()) { th.join(); } + for (auto& th : vThreadList) { + if (th.joinable()) { + th.join(); + } } - auto Operation = [](size_t acc, const Vector<auto>& v) { return acc + v.size(); }; + auto Operation = [](size_t acc, const auto& v) { return acc + v.size(); }; int nRecoTracks = std::accumulate(fvRecoTracks.begin(), fvRecoTracks.end(), 0, Operation); int nRecoHits = std::accumulate(fvRecoHitIndices.begin(), fvRecoHitIndices.end(), 0, Operation); frAlgo.fRecoTracks.reserve(nRecoTracks); frAlgo.fRecoHits.reserve(nRecoHits); for (int iTh = 0; iTh < frAlgo.GetNofThreads(); ++iTh) { - frAlgo.fRecoTracks.insert(frAlgo.fRecoTracks.end(), fvRecoTracks[iTh].begin(), fvRecoTracks[iTh].end()); - frAlgo.fRecoHits.insert(frAlgo.fRecoHits.end(), fvRecoHitIndices[iTh].begin(), fvRecoHitIndices[iTh].end()); + frAlgo.fRecoTracks.insert(frAlgo.fRecoTracks.end(), fvRecoTracks[iTh].begin(), fvRecoTracks[iTh].end()); + frAlgo.fRecoHits.insert(frAlgo.fRecoHits.end(), fvRecoHitIndices[iTh].begin(), fvRecoHitIndices[iTh].end()); } } @@ -265,6 +274,7 @@ void TrackFinder::FindTracksThread(int iThread) const int nDataStreams = frAlgo.fInputData.GetNdataStreams(); bool areUntouchedDataLeft = true; // is the whole TS processed std::vector<HitIndex_t> sliceFirstHit(nDataStreams, 0); + auto& wData = frAlgo.fvWData[iThread]; // Define first hit, skip all the hits, which are before the first window for (int iStream = 0; iStream < nDataStreams; ++iStream) { @@ -274,7 +284,7 @@ void TrackFinder::FindTracksThread(int iThread) for (HitIndex_t caHitId = caFstId; caHitId < caLstId; ++caHitId) { CaHitTimeInfo& info = frAlgo.fHitTimeInfo[caHitId]; const ca::Hit& h = frAlgo.fInputData.GetHit(caHitId); - if (frAlgo.fvHitKeyFlags[h.FrontKey()] || frAlgo.fvHitKeyFlags[h.BackKey()]) { + if (wData.IsHitKeyUsed(h.FrontKey()) || wData.IsHitKeyUsed(h.BackKey())) { continue; } if (info.fMaxTimeBeforeHit < fvWindowStartThread[iThread]) { @@ -314,8 +324,8 @@ void TrackFinder::FindTracksThread(int iThread) nLoops1++; CaHitTimeInfo& info = frAlgo.fHitTimeInfo[caHitId]; const ca::Hit& h = frAlgo.fInputData.GetHit(caHitId); - if (frAlgo.fvHitKeyFlags[h.FrontKey()] - || frAlgo.fvHitKeyFlags[h.BackKey()]) { // the hit is already reconstructed + if (wData.IsHitKeyUsed(h.FrontKey()) || wData.IsHitKeyUsed(h.BackKey())) { + // the hit is already reconstructed continue; } if (info.fEventTimeMin @@ -409,20 +419,28 @@ void TrackFinder::FindTracksThread(int iThread) // release the track hits for (int i = 0; i < track.fNofHits; i++) { - int caHitId = frAlgo.fvWData[iThread].RecoHitIndex(trackFirstHit + i); - const auto& h = frAlgo.fInputData.GetHit(caHitId); - frAlgo.fvHitKeyFlags[h.FrontKey()] = 0; - frAlgo.fvHitKeyFlags[h.BackKey()] = 0; + int caHitId = frAlgo.fvWData[iThread].RecoHitIndex(trackFirstHit + i); + const auto& h = frAlgo.fInputData.GetHit(caHitId); + wData.IsHitKeyUsed(h.FrontKey()) = 0; + wData.IsHitKeyUsed(h.BackKey()) = 0; } } else { // save the track fvRecoTracks[iThread].push_back(track); + if (frAlgo.IfCollectTrackExtraInfo()) { + auto trkInfo = TrackExtraInfo(); + trkInfo.fTimeStart = (track.fParFirst.GetTime() - fStatTsStart) / (fStatTsEnd - fStatTsStart); + trkInfo.fTimeEnd = (track.fParLast.GetTime() - fStatTsStart) / (fStatTsEnd - fStatTsStart); + trkInfo.fThreadId = iThread; + fvRecoTrackExtraInfo[iThread].push_back(trkInfo); + } + // mark the track hits as used for (int i = 0; i < track.fNofHits; i++) { - int caHitId = frAlgo.fvWData[iThread].RecoHitIndex(trackFirstHit + i); - const auto& h = frAlgo.fInputData.GetHit(caHitId); - frAlgo.fvHitKeyFlags[h.FrontKey()] = 1; - frAlgo.fvHitKeyFlags[h.BackKey()] = 1; + int caHitId = frAlgo.fvWData[iThread].RecoHitIndex(trackFirstHit + i); + const auto& h = frAlgo.fInputData.GetHit(caHitId); + wData.IsHitKeyUsed(h.FrontKey()) = 1; + wData.IsHitKeyUsed(h.BackKey()) = 1; fvRecoHitIndices[iThread].push_back(caHitId); } } @@ -445,22 +463,30 @@ void TrackFinder::Init() assert(frAlgo.GetNofThreads() > 0); fvWindowStartThread.clear(); - fvWindowEndThread.clear(); - fvStatNwindows.clear(); - fvStatNhitsProcessed.clear(); - fvRecoTracks.clear(); - fvRecoHitIndices.clear(); - fvWindowStartThread.resize(frAlgo.GetNofThreads()); + fvWindowEndThread.clear(); fvWindowEndThread.resize(frAlgo.GetNofThreads()); + fvStatNwindows.clear(); fvStatNwindows.resize(frAlgo.GetNofThreads()); + fvStatNhitsProcessed.clear(); fvStatNhitsProcessed.resize(frAlgo.GetNofThreads()); + fvRecoTracks.clear(); fvRecoTracks.resize(frAlgo.GetNofThreads()); + fvRecoHitIndices.clear(); fvRecoHitIndices.resize(frAlgo.GetNofThreads()); + if (frAlgo.IfCollectTrackExtraInfo()) { + fvRecoTrackExtraInfo.clear(); + fvRecoTrackExtraInfo.resize(frAlgo.GetNofThreads()); + } + for (int iThread = 0; iThread < frAlgo.GetNofThreads(); ++iThread) { fvRecoTracks[iThread].SetName(std::string("TrackFinder::fvRecoTracks_") + std::to_string(iThread)); fvRecoHitIndices[iThread].SetName(std::string("TrackFinder::fvRecoHitIndices_") + std::to_string(iThread)); + if (frAlgo.IfCollectTrackExtraInfo()) { + fvRecoTrackExtraInfo[iThread].SetName(std::string("TrackFinder::fvRecoTrackExtraInfo_") + + std::to_string(iThread)); + } } fWindowLength = (ca::Framework::TrackingMode::kMcbm == frAlgo.fTrackingMode) ? 500 : 10000; diff --git a/algo/ca/core/tracking/CaTrackFinder.h b/algo/ca/core/tracking/CaTrackFinder.h index a8f33ab701376ca9f2841369dd6e97181e6cb0b5..107f81a3355a6322dc592a016720aa0cb7ef2e73 100644 --- a/algo/ca/core/tracking/CaTrackFinder.h +++ b/algo/ca/core/tracking/CaTrackFinder.h @@ -70,8 +70,9 @@ namespace cbm::algo::ca std::vector<int> fvStatNwindows; std::vector<int> fvStatNhitsProcessed; - std::vector<Vector<Track>> fvRecoTracks; ///< reconstructed tracks - std::vector<Vector<ca::HitIndex_t>> fvRecoHitIndices; ///< packed hits of reconstructed tracks + std::vector<Vector<Track>> fvRecoTracks; ///< reconstructed tracks + std::vector<Vector<TrackExtraInfo>> fvRecoTrackExtraInfo; ///< Extra info on tracks + std::vector<Vector<ca::HitIndex_t>> fvRecoHitIndices; ///< packed hits of reconstructed tracks float fWindowLength = 0.; float fWindowOverlap = 15.; // ns diff --git a/algo/ca/core/tracking/CaTrackFinderWindow.cxx b/algo/ca/core/tracking/CaTrackFinderWindow.cxx index 5e8f2676091102d6a8441263310984e6f093534f..646185bc5945c6950b87503ae124738a88d046d5 100644 --- a/algo/ca/core/tracking/CaTrackFinderWindow.cxx +++ b/algo/ca/core/tracking/CaTrackFinderWindow.cxx @@ -119,10 +119,7 @@ bool TrackFinderWindow::checkTripletMatch(const ca::Triplet& l, const ca::Triple return true; } -void TrackFinderWindow::InitTimeslice() -{ - fvHitKeyToTrack.reset(frAlgo.fInputData.GetNhitKeys(), -1); -} +void TrackFinderWindow::InitTimeslice() { fvHitKeyToTrack.reset(frAlgo.fInputData.GetNhitKeys(), -1); } // ************************************************************************************************** // * * @@ -250,7 +247,7 @@ void TrackFinderWindow::CaTrackFinderSlice() grid.StoreHits(frWData.Hits(), frWData.HitStartIndexOnStation(iS), frWData.NofHitsOnStation(iS), - frAlgo.fvHitKeyFlags); + frWData.HitKeyFlags()); /* clang-format on */ } @@ -267,7 +264,7 @@ void TrackFinderWindow::CaTrackFinderSlice() if (frAlgo.fCurrentIterationIndex > 0) { for (int ista = 0; ista < frAlgo.GetParameters().GetNstationsActive(); ++ista) { - frWData.Grid(ista).RemoveUsedHits(frWData.Hits(), frAlgo.fvHitKeyFlags); + frWData.Grid(ista).RemoveUsedHits(frWData.Hits(), frWData.HitKeyFlags()); } } @@ -447,7 +444,7 @@ void TrackFinderWindow::CaTrackFinderSlice() ca::Triplet& first_trip = (fvTriplets[istaF][itrip]); const auto& fstTripLHit = frWData.Hit(first_trip.GetLHit()); - if (frAlgo.fvHitKeyFlags[fstTripLHit.FrontKey()] || frAlgo.fvHitKeyFlags[fstTripLHit.BackKey()]) { + if (frWData.IsHitKeyUsed(fstTripLHit.FrontKey()) || frWData.IsHitKeyUsed(fstTripLHit.BackKey())) { continue; } @@ -639,12 +636,9 @@ void TrackFinderWindow::CaTrackFinderSlice() for (auto iHit : tr.Hits()) { const ca::Hit& hit = frAlgo.fInputData.GetHit(iHit); - /// used strips are marked - - frAlgo.fvHitKeyFlags[hit.FrontKey()] = 1; - frAlgo.fvHitKeyFlags[hit.BackKey()] = 1; - + frWData.IsHitKeyUsed(hit.FrontKey()) = 1; + frWData.IsHitKeyUsed(hit.BackKey()) = 1; frWData.RecoHitIndices().push_back(iHit); } Track t; @@ -669,8 +663,8 @@ void TrackFinderWindow::CaTrackFinderSlice() for (unsigned int ih = 0; ih < frWData.Hits().size(); ih++) { if (frWData.IsHitSuppressed(ih)) { const ca::Hit& hit = frWData.Hit(ih); - frAlgo.fvHitKeyFlags[hit.FrontKey()] = 1; - frAlgo.fvHitKeyFlags[hit.BackKey()] = 1; + frWData.IsHitKeyUsed(hit.FrontKey()) = 1; + frWData.IsHitKeyUsed(hit.BackKey()) = 1; } } } // ---- Loop over Track Finder iterations: END -----------------------------------------------------------// @@ -712,11 +706,11 @@ void TrackFinderWindow::CAFindTrack(int ista, ca::Branch& best_tr, const ca::Tri const auto& hitM = frWData.Hit(curr_trip->GetMHit()); const auto& hitR = frWData.Hit(curr_trip->GetRHit()); - if (!(frAlgo.fvHitKeyFlags[hitM.FrontKey()] || frAlgo.fvHitKeyFlags[hitM.BackKey()])) { + if (!(frWData.IsHitKeyUsed(hitM.FrontKey()) || frWData.IsHitKeyUsed(hitM.BackKey()))) { curr_tr.AddHit(hitM.Id()); } - if (!(frAlgo.fvHitKeyFlags[hitR.FrontKey()] || frAlgo.fvHitKeyFlags[hitR.BackKey()])) { + if (!(frWData.IsHitKeyUsed(hitR.FrontKey()) || frWData.IsHitKeyUsed(hitR.BackKey()))) { curr_tr.AddHit(hitR.Id()); } @@ -759,7 +753,7 @@ void TrackFinderWindow::CAFindTrack(int ista, ca::Branch& best_tr, const ca::Tri if (!checkTripletMatch(*curr_trip, new_trip, dchi2)) continue; const auto& hitL = frWData.Hit(new_trip.GetLHit()); // left hit of new triplet - if (frAlgo.fvHitKeyFlags[hitL.FrontKey()] || frAlgo.fvHitKeyFlags[hitL.BackKey()]) { //hits are used + if (frWData.IsHitKeyUsed(hitL.FrontKey()) || frWData.IsHitKeyUsed(hitL.BackKey())) { //hits are used // no used hits allowed -> compare and store track if ((curr_tr.NofHits() > best_tr.NofHits()) || ((curr_tr.NofHits() == best_tr.NofHits()) && (curr_tr.Chi2() < best_tr.Chi2()))) {