diff --git a/algo/ca/core/tracking/CaTrackFinder.cxx b/algo/ca/core/tracking/CaTrackFinder.cxx index 9e82e6c46d8ba4d1f647603899a802765135e583..e7679e3efe809a3d45eb55a6fc2b1fb24c35b413 100644 --- a/algo/ca/core/tracking/CaTrackFinder.cxx +++ b/algo/ca/core/tracking/CaTrackFinder.cxx @@ -19,6 +19,7 @@ */ #include "CaFramework.h" +#include "CaTimer.h" #include "CaTrack.h" #include <chrono> @@ -82,11 +83,9 @@ void TrackFinder::FindTracks() fscal targY = frAlgo.GetParameters().GetTargetPositionY()[0]; fscal targZ = frAlgo.GetParameters().GetTargetPositionZ()[0]; - fscal tsStart = std::numeric_limits<fscal>::max(); // starting time of sub-TS - - fscal statTsStart = 0.; // end time of the TS - fscal statTsEnd = 0.; // end time of the TS - int statNhitsTotal = 0; + fStatTsStart = std::numeric_limits<fscal>::max(); // end time of the TS + fStatTsEnd = 0.; // end time of the TS + fStatNhitsTotal = 0; // calculate possible event time for the hits (frAlgo.fHitTimeInfo array) @@ -97,7 +96,7 @@ void TrackFinder::FindTracks() fscal maxTimeBeforeHit = std::numeric_limits<fscal>::lowest(); int nStreamHits = frAlgo.fInputData.GetStreamNhits(iStream); - statNhitsTotal += nStreamHits; + fStatNhitsTotal += nStreamHits; for (int ih = 0; ih < nStreamHits; ++ih) { @@ -140,11 +139,11 @@ void TrackFinder::FindTracks() } info.fMaxTimeBeforeHit = maxTimeBeforeHit; - if (tsStart > info.fEventTimeMax - 5) { - tsStart = info.fEventTimeMax - 5; // 5 ns margin + if (fStatTsStart > info.fEventTimeMax - 5) { + fStatTsStart = info.fEventTimeMax - 5; // 5 ns margin (FIXME: define in constants!!!) } - if (statTsEnd < info.fEventTimeMin) { - statTsEnd = info.fEventTimeMin; + if (fStatTsEnd < info.fEventTimeMin) { + fStatTsEnd = info.fEventTimeMin; } } @@ -187,202 +186,247 @@ void TrackFinder::FindTracks() } } - statTsStart = tsStart; - if (statTsEnd < statTsStart) { // all hits belong to one sub-timeslice - // statTsEnd = statTsStart; + if (fStatTsEnd < fStatTsStart) { // all hits belong to one sub-timeslice + // statTsEnd = statTsStart; } - if (statTsEnd > statTsStart + 500.e6) { // 500 ms maximal length of the TS - statTsEnd = statTsStart + 500.e6; + if (fStatTsEnd > fStatTsStart + 500.e6) { // 500 ms maximal length of the TS + fStatTsEnd = fStatTsStart + 500.e6; } - LOG(info) << "CA tracker process time slice " << statTsStart / 1.e6 << " ms ... " << statTsEnd / 1.e6 << " ms with " - << statNhitsTotal << " hits"; + LOG(info) << "CA tracker process time slice " << fStatTsStart / 1.e6 << " ms ... " << fStatTsEnd / 1.e6 << " ms with " + << fStatNhitsTotal << " hits"; // cut data into sub-timeslices and process them one by one - - std::vector<int> statNwindows(frAlgo.GetNofThreads(), 0); - std::vector<int> statNhitsProcessed(frAlgo.GetNofThreads(), 0); + float threadDuration = (fStatTsEnd - fStatTsStart) / frAlgo.GetNofThreads(); + for (int iThread = 0; iThread < frAlgo.GetNofThreads(); ++iThread) { + fvWindowStartThread[iThread] = fStatTsStart + iThread * threadDuration; + fvWindowEndThread[iThread] = fvWindowStartThread[iThread] + threadDuration; + LOG(info) << "Thread: " << iThread << " from " << fvWindowStartThread[iThread] / 1.e6 << " ms to " + << fvWindowEndThread[iThread] / 1.e6 << " ms"; + } for (int iThread = 0; iThread < frAlgo.GetNofThreads(); ++iThread) { - bool areUntouchedDataLeft = true; // is the whole TS processed - ca::HitIndex_t sliceFirstHit[nDataStreams]; + this->FindTracksThread(iThread); + } - for (int iStream = 0; iStream < nDataStreams; ++iStream) { - sliceFirstHit[iStream] = frAlgo.fInputData.GetStreamStartIndex(iStream); - } + frAlgo.fMonitorData.StopTimer(ETimer::FindTracks); + frAlgo.fMonitorData.IncrementCounter(ECounter::RecoTrack, frAlgo.fRecoTracks.size()); + frAlgo.fMonitorData.IncrementCounter(ECounter::RecoHitUsed, frAlgo.fRecoHits.size()); - statNwindows[iThread] = 0; - statNhitsProcessed[iThread] = 0; + auto timerEnd = std::chrono::high_resolution_clock::now(); + frAlgo.fCaRecoTime = (double) (std::chrono::duration<double>(timerEnd - timerStart).count()); + + int statNhitsProcessedTotal = 0; + int statNwindowsTotal = 0; + for (int iThread = 0; iThread < frAlgo.GetNofThreads(); ++iThread) { + statNhitsProcessedTotal += fvStatNhitsProcessed[iThread]; + statNwindowsTotal += fvStatNwindows[iThread]; + } + + LOG(info) << "CA tracker: time slice finished. Reconstructed " << frAlgo.fRecoTracks.size() << " tracks with " + << frAlgo.fRecoHits.size() << " hits. Processed " << statNhitsProcessedTotal << " hits in " + << statNwindowsTotal << " time windows. Reco time " << frAlgo.fCaRecoTime / 1.e9 << " s"; +} - int statLastLogTimeChunk = -1; +// --------------------------------------------------------------------------------------------------------------------- +// +void TrackFinder::FindTracksThread(int iThread) +{ + LOG(info) << "---- CA: searching for tracks on thread " << iThread; + const int nDataStreams = frAlgo.fInputData.GetNdataStreams(); + bool areUntouchedDataLeft = true; // is the whole TS processed + std::vector<HitIndex_t> sliceFirstHit(nDataStreams, 0); - while (areUntouchedDataLeft) { - frAlgo.fMonitorData.IncrementCounter(ECounter::SubTS); - // select the sub-slice hits - for (int iS = 0; iS < frAlgo.GetParameters().GetNstationsActive(); ++iS) { - frAlgo.fvWData[iThread].TsHitIndices(iS).clear(); + // Define first hit, skip all the hits, which are before the first window + for (int iStream = 0; iStream < nDataStreams; ++iStream) { + sliceFirstHit[iStream] = frAlgo.fInputData.GetStreamStartIndex(iStream); + HitIndex_t caLstId = frAlgo.fInputData.GetStreamStopIndex(iStream); + HitIndex_t caFstId = sliceFirstHit[iStream]; + 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()]) { + continue; + } + if (info.fMaxTimeBeforeHit < fvWindowStartThread[iThread]) { + // this hit and all hits before are before the overlap + sliceFirstHit[iStream] = caHitId + 1; } + } + } - areUntouchedDataLeft = false; + fvStatNwindows[iThread] = 0; + fvStatNhitsProcessed[iThread] = 0; - // TODO: SG: skip empty regions and start the subslice with the earliest hit + int statLastLogTimeChunk = -1; - statNwindows[iThread]++; - int statNwindowHits = 0; + while (areUntouchedDataLeft) { + frAlgo.fMonitorData.IncrementCounter(ECounter::SubTS); + // select the sub-slice hits + for (int iS = 0; iS < frAlgo.GetParameters().GetNstationsActive(); ++iS) { + frAlgo.fvWData[iThread].TsHitIndices(iS).clear(); + } - for (int iStream = 0; iStream < nDataStreams; ++iStream) { + areUntouchedDataLeft = false; - for (ca::HitIndex_t caHitId = sliceFirstHit[iStream]; caHitId < frAlgo.fInputData.GetStreamStopIndex(iStream); - ++caHitId) { + // TODO: SG: skip empty regions and start the subslice with the earliest hit - CaHitTimeInfo& info = frAlgo.fHitTimeInfo[caHitId]; - const ca::Hit& h = frAlgo.fInputData.GetHit(caHitId); + fvStatNwindows[iThread]++; - if (frAlgo.fvHitKeyFlags[h.FrontKey()] - || frAlgo.fvHitKeyFlags[h.BackKey()]) { // the hit is already reconstructed - continue; + int statNwindowHits = 0; + + Timer tHitInit; + tHitInit.Start(); + int nLoops1 = 0; + int nLoops2 = 0; + for (int iStream = 0; iStream < nDataStreams; ++iStream) { + for (ca::HitIndex_t caHitId = sliceFirstHit[iStream]; caHitId < frAlgo.fInputData.GetStreamStopIndex(iStream); + ++caHitId) { + 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 + continue; + } + if (info.fEventTimeMin + > fvWindowStartThread[iThread] + fWindowLength) { // the hit is too late for the sub slice + areUntouchedDataLeft = true; + if (info.fMinTimeAfterHit > fvWindowStartThread[iThread] + fWindowLength) { + // this hit and all later hits are out of the sub-slice + break; } - if (info.fEventTimeMin > tsStart + fWindowLength) { // the hit is too late for the sub slice - areUntouchedDataLeft = true; - if (info.fMinTimeAfterHit > tsStart + fWindowLength) { - // this hit and all later hits are out of the sub-slice - break; + } + else { + if (fvWindowStartThread[iThread] <= info.fEventTimeMax) { // the hit belongs to the sub-slice + frAlgo.fvWData[iThread].TsHitIndices(h.Station()).push_back(caHitId); + statNwindowHits++; + if (info.fMaxTimeBeforeHit < fvWindowStartThread[iThread] + fWindowLength - fWindowOverlap) { + // this hit and all hits before are before the overlap + sliceFirstHit[iStream] = caHitId + 1; } + //LOG(warning) << "event time " << ((float) info.fEventTimeMax )<< " " << h.ToString(); } - else { - if (tsStart <= info.fEventTimeMax) { // the hit belongs to the sub-slice - frAlgo.fvWData[iThread].TsHitIndices(h.Station()).push_back(caHitId); - statNwindowHits++; - if (info.fMaxTimeBeforeHit < tsStart + fWindowLength - fWindowOverlap) { - // this hit and all hits before are before the overlap - sliceFirstHit[iStream] = caHitId + 1; - } - //LOG(warning) << "event time " << ((float) info.fEventTimeMax )<< " " << h.ToString(); - } - } // else the hit has been alread processed in previous sub-slices - } + } // else the hit has been alread processed in previous sub-slices + nLoops2++; } - statNhitsProcessed[iThread] += statNwindowHits; - - // print the LOG for every 10 ms of data processed - { - int currentChunk = (int) ((tsStart - statTsStart) / 10.e6); - if (!areUntouchedDataLeft || currentChunk > statLastLogTimeChunk) { - statLastLogTimeChunk = currentChunk; - double dataRead = 100. * (tsStart + fWindowLength - statTsStart) / (statTsEnd - statTsStart); - if (dataRead > 100.) { - dataRead = 100.; - } - LOG(info) << "CA tracker process sliding window N " << statNwindows[iThread] << ": time " << tsStart / 1.e6 - << " ms + " << fWindowLength / 1.e3 << " us) with " << statNwindowHits << " hits. " - << " Processing " << dataRead << " % of the TS time and " - << 100. * statNhitsProcessed[iThread] / statNhitsTotal << " % of TS hits." - << " Already reconstructed " << frAlgo.fRecoTracks.size() << " tracks "; + } + fvStatNhitsProcessed[iThread] += statNwindowHits; + tHitInit.Stop(); + + // print the LOG for every 10 ms of data processed + { + int currentChunk = (int) ((fvWindowStartThread[iThread] - fStatTsStart) / 10.e6); + if (!areUntouchedDataLeft || currentChunk > statLastLogTimeChunk) { + statLastLogTimeChunk = currentChunk; + double dataRead = + 100. * (fvWindowStartThread[iThread] + fWindowLength - fStatTsStart) / (fStatTsEnd - fStatTsStart); + if (dataRead > 100.) { + dataRead = 100.; } + LOG(info) << "CA tracker process sliding window N " << fvStatNwindows[iThread] << ": time " + << fvWindowStartThread[iThread] / 1.e6 << " ms + " << fWindowLength / 1.e3 << " us) with " + << statNwindowHits << " hits. " + << " Processing " << dataRead << " % of the TS time and " + << 100. * fvStatNhitsProcessed[iThread] / fStatNhitsTotal << " % of TS hits." + << " Already reconstructed " << frAlgo.fRecoTracks.size() << " tracks "; } + } - 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(); - } + 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(); } } + } - frAlgo.fMonitorData.StartTimer(ETimer::TrackFinder); - frAlgo.fvTrackFinderWindow[iThread].CaTrackFinderSlice(); - frAlgo.fMonitorData.StopTimer(ETimer::TrackFinder); - - - // save reconstructed tracks with no hits in the overlap region - - tsStart += 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 - - int trackFirstHit = 0; - for (const auto& track : frAlgo.fvWData[iThread].RecoTracks()) { - bool isTrackCompletelyInOverlap = true; - for (int ih = 0; ih < track.fNofHits; ih++) { - int caHitId = frAlgo.fvWData[iThread].RecoHitIndex(trackFirstHit + ih); - CaHitTimeInfo& info = frAlgo.fHitTimeInfo[caHitId]; - if (info.fEventTimeMax < tsStart) { // this hit is before the overlap - isTrackCompletelyInOverlap = false; - break; - } + Timer tTracking; + tTracking.Start(); + frAlgo.fMonitorData.StartTimer(ETimer::TrackFinder); + frAlgo.fvTrackFinderWindow[iThread].CaTrackFinderSlice(); + frAlgo.fMonitorData.StopTimer(ETimer::TrackFinder); + tTracking.Stop(); + + // save reconstructed tracks with no hits in the overlap region + //if (fvWindowStartThread[iThread] > 13.23e6 && fvWindowStartThread[iThread] < 13.26e6) { + fvWindowStartThread[iThread] += 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 + + int trackFirstHit = 0; + for (const auto& track : frAlgo.fvWData[iThread].RecoTracks()) { + bool isTrackCompletelyInOverlap = true; + for (int ih = 0; ih < track.fNofHits; ih++) { + int caHitId = frAlgo.fvWData[iThread].RecoHitIndex(trackFirstHit + ih); + CaHitTimeInfo& info = frAlgo.fHitTimeInfo[caHitId]; + if (info.fEventTimeMax < fvWindowStartThread[iThread]) { // this hit is before the overlap + isTrackCompletelyInOverlap = false; + break; } + } - if (!areUntouchedDataLeft) { // don't reject tracks in the overlap when no more data are left - isTrackCompletelyInOverlap = 0; - } + if (!areUntouchedDataLeft) { // don't reject tracks in the overlap when no more data are left + isTrackCompletelyInOverlap = 0; + } - if (isTrackCompletelyInOverlap) { - // - // Don't save tracks from the overlap region, since they might have additional hits in the next subslice - // - - // 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; - } + if (isTrackCompletelyInOverlap) { + // + // Don't save tracks from the overlap region, since they might have additional hits in the next subslice + // + + // 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; } - else { // save the track - frAlgo.fRecoTracks.push_back(track); - // 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; - frAlgo.fRecoHits.push_back(caHitId); - } + } + else { // save the track + frAlgo.fRecoTracks.push_back(track); + // 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; + frAlgo.fRecoHits.push_back(caHitId); } - trackFirstHit += track.fNofHits; - } // sub-timeslice tracks + } + trackFirstHit += track.fNofHits; + } // sub-timeslice tracks - tsStart -= 5; // do 5 ns margin + if (fvWindowStartThread[iThread] > fvWindowEndThread[iThread]) { + break; + } + else { + fvWindowStartThread[iThread] -= 5; // do 5 ns margin } } - - frAlgo.fMonitorData.StopTimer(ETimer::FindTracks); - frAlgo.fMonitorData.IncrementCounter(ECounter::RecoTrack, frAlgo.fRecoTracks.size()); - frAlgo.fMonitorData.IncrementCounter(ECounter::RecoHitUsed, frAlgo.fRecoHits.size()); - - auto timerEnd = std::chrono::high_resolution_clock::now(); - frAlgo.fCaRecoTime = (double) (std::chrono::duration<double>(timerEnd - timerStart).count()); - - int statNhitsProcessedTotal = 0; - int statNwindowsTotal = 0; - for (int iThread = 0; iThread < frAlgo.GetNofThreads(); ++iThread) { - statNhitsProcessedTotal += statNhitsProcessed[iThread]; - statNwindowsTotal += statNwindows[iThread]; - } - - LOG(info) << "CA tracker: time slice finished. Reconstructed " << frAlgo.fRecoTracks.size() << " tracks with " - << frAlgo.fRecoHits.size() << " hits. Processed " << statNhitsProcessedTotal << " hits in " - << statNwindowsTotal << " time windows. Reco time " << frAlgo.fCaRecoTime / 1.e9 << " s"; } // --------------------------------------------------------------------------------------------------------------------- // void TrackFinder::Init() { - fvTsStartThread.clear(); - fvTsEndThread.clear(); + assert(frAlgo.GetNofThreads() > 0); + fvWindowStartThread.clear(); + fvWindowEndThread.clear(); + fvStatNwindows.clear(); + fvStatNhitsProcessed.clear(); - fvTsStartThread.resize(frAlgo.GetNofThreads()); - fvTsEndThread.resize(frAlgo.GetNofThreads()); fvWindowStartThread.resize(frAlgo.GetNofThreads()); + fvWindowEndThread.resize(frAlgo.GetNofThreads()); + fvStatNwindows.resize(frAlgo.GetNofThreads()); + fvStatNhitsProcessed.resize(frAlgo.GetNofThreads()); + 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 9bfe0f9bcb77cc3c8f3b935d1bb9e66d92c20cef..ba5eca97e06eb0d213c63a92419b4a8914577988 100644 --- a/algo/ca/core/tracking/CaTrackFinder.h +++ b/algo/ca/core/tracking/CaTrackFinder.h @@ -48,6 +48,7 @@ namespace cbm::algo::ca TrackFinder& operator=(TrackFinder&&) = delete; void FindTracks(); + void FindTracksThread(int iThread); void Init(); private: @@ -63,12 +64,18 @@ namespace cbm::algo::ca ca::Framework& frAlgo; ///< Reference to the main track finder algorithm class - std::vector<float> fvTsStartThread; - std::vector<float> fvTsEndThread; std::vector<float> fvWindowStartThread; + std::vector<float> fvWindowEndThread; + + std::vector<int> fvStatNwindows; + std::vector<int> fvStatNhitsProcessed; float fWindowLength = 0.; float fWindowOverlap = 15.; // ns + + float fStatTsStart = 0.; + float fStatTsEnd = 0.; + int fStatNhitsTotal = 0; }; } // namespace cbm::algo::ca