diff --git a/algo/ca/core/tracking/CaTrackFinder.cxx b/algo/ca/core/tracking/CaTrackFinder.cxx index df883016cedf441f30ef633dbb2e8dece609adfa..695dabeaefd6ed6777c86132883580841ed1db79 100644 --- a/algo/ca/core/tracking/CaTrackFinder.cxx +++ b/algo/ca/core/tracking/CaTrackFinder.cxx @@ -23,6 +23,8 @@ #include "CaTrack.h" #include <chrono> +#include <fstream> +#include <sstream> #include <thread> using namespace cbm::algo::ca; @@ -193,19 +195,85 @@ void TrackFinder::FindTracks() int nWindows = static_cast<int>((fStatTsEnd - fStatTsStart - fWindowOverlap - fWindowMargin) / (fWindowLength - fWindowOverlap - fWindowMargin)) + 1; - LOG(info) << "CA: estimated number of time windows: " << nWindows; - int nWindowsThread = nWindows / frAlgo.GetNofThreads(); float windowDelta = fWindowLength - fWindowOverlap - fWindowMargin; + //LOG(info) << "CA: estimated number of time windows: " << nWindows; + + + { // Estimation of number of hits in time windows + //Timer time; + //time.Start(); + int nSt = frAlgo.GetParameters().GetNstationsActive(); + + // Count number of hits per window and station + double a = fStatTsStart + fWindowOverlap + fWindowMargin; + double b = fWindowLength - fWindowOverlap - fWindowMargin; + HitIndex_t nHitsTot = frAlgo.fInputData.GetNhits(); + std::vector<HitIndex_t> nHitsWindowSta(nWindows * nSt, 0); + for (HitIndex_t iHit = 0; iHit < nHitsTot; ++iHit) { + const auto& hit = frAlgo.fInputData.GetHit(iHit); + //const auto& timeInfo = frAlgo.fHitTimeInfo[iHit]; + 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::Framework::TrackingMode::kMcbm == frAlgo.fTrackingMode) { + auto maxNofHitsSta = static_cast<HitIndex_t>(50 * fWindowLength / 1.e3); + for (auto& content : nHitsWindowSta) { + if (content > maxNofHitsSta) { + content = 0; + } + } + } + + // Integrate number of hits per window + std::vector<HitIndex_t> nHitsWindow; + nHitsWindow.reserve(nWindows); + auto windowStaIt = nHitsWindowSta.begin(); + HitIndex_t nHitsCollected = 0; + while (windowStaIt != nHitsWindowSta.end()) { + nHitsCollected = nHitsWindow.emplace_back(std::accumulate(windowStaIt, windowStaIt + nSt, nHitsCollected)); + windowStaIt = windowStaIt + nSt; + } + + // Get time range for threads + HitIndex_t nHitsPerThread = nHitsCollected / frAlgo.GetNofThreads(); + auto windowIt = nHitsWindow.begin(); + size_t iWbegin = 0; + fvWindowStartThread[0] = fStatTsStart; + for (int iTh = 1; iTh < frAlgo.GetNofThreads(); ++iTh) { + windowIt = std::lower_bound(windowIt, nHitsWindow.end(), iTh * nHitsPerThread); + iWbegin = std::distance(nHitsWindow.begin(), windowIt) + 1; + fvWindowStartThread[iTh] = fStatTsStart + iWbegin * windowDelta; + fvWindowEndThread[iTh - 1] = fvWindowStartThread[iTh]; + } + //time.Stop(); + //LOG(info) << "Thread boarders estimation time: " << time.GetTotalMs() << " ms"; + LOG(info) << "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 < frAlgo.GetNofThreads(); ++iThread) { + // fvWindowStartThread[iThread] = fStatTsStart + iThread * nWindowsThread * windowDelta; + // fvWindowEndThread[iThread] = + // fvWindowStartThread[iThread] + nWindowsThread * windowDelta + fWindowOverlap + fWindowMargin; + //} + fvWindowEndThread[frAlgo.GetNofThreads() - 1] = fStatTsEnd; + for (int iThread = 0; iThread < frAlgo.GetNofThreads(); ++iThread) { - fvWindowStartThread[iThread] = fStatTsStart + iThread * nWindowsThread * windowDelta; - fvWindowEndThread[iThread] = - fvWindowStartThread[iThread] + nWindowsThread * windowDelta + fWindowOverlap + fWindowMargin; - LOG(info) << "Thread: " << iThread << " from " << fvWindowStartThread[iThread] / 1.e6 << " ms to " - << fvWindowEndThread[iThread] / 1.e6 << " ms"; + double start = fvWindowStartThread[iThread] * 1.e-6; + double end = fvWindowEndThread[iThread] * 1.e-6; + LOG(info) << "Thread: " << iThread << " from " << start << " ms to " << end << " ms (delta = " << end - start + << " ms)"; } - fvWindowEndThread[frAlgo.GetNofThreads() - 1] = fStatTsEnd; frAlgo.fMonitorData.StopTimer(ETimer::PrepareTimeslice); // Save tracks @@ -251,6 +319,7 @@ void TrackFinder::FindTracks() // Add thread monitors to the main monitor for (auto& monitor : frAlgo.fvMonitorDataThread) { frAlgo.fMonitorData.AddMonitorData(monitor, true); + //frAlgo.fMonitorData.AddMonitorData(monitor); monitor.Reset(); } @@ -277,7 +346,12 @@ void TrackFinder::FindTracks() // void TrackFinder::FindTracksThread(int iThread) { + //std::stringstream filename; + //filename << "./dbg_caTrackFinder::FindTracksThread_" << iThread << ".txt"; + //std::ofstream out(filename.str()); auto& monitor = frAlgo.fvMonitorDataThread[iThread]; + Timer timer; + timer.Start(); monitor.StartTimer(ETimer::TrackingThread); monitor.StartTimer(ETimer::PrepareThread); @@ -353,6 +427,7 @@ void TrackFinder::FindTracksThread(int iThread) // TODO: SG: skip empty regions and start the subslice with the earliest hit fvStatNwindows[iThread]++; + //out << fvStatNwindows[iThread] << ' '; int statNwindowHits = 0; @@ -385,10 +460,30 @@ void TrackFinder::FindTracksThread(int iThread) } // else the hit has been alread processed in previous sub-slices } } + + //out << statNwindowHits << ' '; + //if (statNwindowHits == 0) { // Empty window + // frAlgo.fvMonitorDataThread[iThread].StopTimer(ETimer::PrepareWindow); + // out << 0 << ' ' << 0 << ' ' << 0 << '\n'; + // continue; + //} + + if (ca::Framework::TrackingMode::kMcbm == frAlgo.fTrackingMode) { + // cut at 50 hits per station per 1 us. + int maxStationHits = (int) (50 * fWindowLength / 1.e3); + for (int ista = 0; ista < frAlgo.GetParameters().GetNstationsActive(); ++ista) { + int nHitsSta = static_cast<int>(frAlgo.fvWData[iThread].TsHitIndices(ista).size()); + if (nHitsSta > maxStationHits) { + statNwindowHits -= nHitsSta; + frAlgo.fvWData[iThread].TsHitIndices(ista).clear(); + } + } + } + fvStatNhitsProcessed[iThread] += statNwindowHits; // print the LOG for every 10 ms of data processed - { + if constexpr (0) { int currentChunk = (int) ((fvWindowStartThread[iThread] - fStatTsStart) / 10.e6); if (!areUntouchedDataLeft || currentChunk > statLastLogTimeChunk) { statLastLogTimeChunk = currentChunk; @@ -406,21 +501,18 @@ void TrackFinder::FindTracksThread(int iThread) } } - if (ca::Framework::TrackingMode::kMcbm == frAlgo.fTrackingMode) { - // cut at 50 hits per station per 1 us. - int maxStationHits = (int) (50 * fWindowLength / 1.e3); - for (int ista = 0; ista < frAlgo.GetParameters().GetNstationsActive(); ++ista) { - if ((int) frAlgo.fvWData[iThread].TsHitIndices(ista).size() > maxStationHits) { - frAlgo.fvWData[iThread].TsHitIndices(ista).clear(); - } - } - } + + //out << statNwindowHits << ' '; frAlgo.fvMonitorDataThread[iThread].StopTimer(ETimer::PrepareWindow); + //Timer trackingInWindow; //DBG + //trackingInWindow.Start(); frAlgo.fvMonitorDataThread[iThread].StartTimer(ETimer::TrackingWindow); frAlgo.fvTrackFinderWindow[iThread].CaTrackFinderSlice(); frAlgo.fvMonitorDataThread[iThread].StopTimer(ETimer::TrackingWindow); + //trackingInWindow.Stop(); + //out << trackingInWindow.GetTotalMs() << ' '; // save reconstructed tracks with no hits in the overlap region //if (fvWindowStartThread[iThread] > 13.23e6 && fvWindowStartThread[iThread] < 13.26e6) { @@ -428,6 +520,8 @@ void TrackFinder::FindTracksThread(int iThread) // 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 << frAlgo.fvWData[iThread].RecoHitIndices().size() << ' '; + //out << frAlgo.fvWData[iThread].RecoTracks().size() << '\n'; frAlgo.fvMonitorDataThread[iThread].StartTimer(ETimer::StoreTracksWindow); int trackFirstHit = 0; @@ -482,6 +576,10 @@ void TrackFinder::FindTracksThread(int iThread) } } frAlgo.fvMonitorDataThread[iThread].StopTimer(ETimer::TrackingThread); + //timer.Stop(); + //LOG(info) << "CA: finishing tracking on thread " << iThread << " (time: " << timer.GetTotalMs() << " ms, " + // << "hits processed: " << fvStatNhitsProcessed[iThread] << ", " + // << "hits used: " << fvRecoHitIndices[iThread].size() << ')'; } // ---------------------------------------------------------------------------------------------------------------------