CaTrackFinder.cxx

Go to the documentation of this file.

/* Copyright (C) 2009-2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
   SPDX-License-Identifier: GPL-3.0-only
   Authors: Valentina Akishina, Ivan Kisel, Sergey Gorbunov [committer], Igor Kulakov, Sergei Zharko, Maksym Zyzak */
 
/*
 *=====================================================
 *
 *  CBM Level 1 4D Reconstruction
 *
 *  Authors: V.Akishina, I.Kisel,  S.Gorbunov, I.Kulakov, M.Zyzak
 *  Documentation: V.Akishina
 *
 *  e-mail : v.akishina@gsi.de
 *
 *=====================================================
 *
 *  Finds tracks using the Cellular Automaton algorithm
 *
 */
 
#include "CaTrackFinder.h"
 
#include "CaTimer.h"
#include "CaTrack.h"
#include "CaTriplet.h"
 
#include <chrono>
#include <fstream>
#include <sstream>
#include <thread>
 
 
namespace cbm::algo::ca
{
  using constants::phys::ProtonMassD;
  using constants::phys::SpeedOfLightInv;
  using constants::phys::SpeedOfLightInvD;
 
  // -------------------------------------------------------------------------------------------------------------------
  //
  TrackFinder::TrackFinder(const ca::Parameters<fvec>& pars, const fscal mass, const ca::TrackingMode& mode,
                           TrackingMonitorData& monitorData, int nThreads, double& recoTime)
    : fParameters(pars)
    , fDefaultMass(mass)
    , fTrackingMode(mode)
    , fMonitorData(monitorData)
    , fvMonitorDataThread(nThreads)
    , fvWData(nThreads)
    , fNofThreads(nThreads)
    , fCaRecoTime(recoTime)
    , fvRecoTracks(nThreads)
    , fvRecoHitIndices(nThreads)
    , fWindowLength((ca::TrackingMode::kMcbm == mode) ? 500 : 10000)
  {
    assert(fNofThreads > 0);
 
    for (int iThread = 0; iThread < fNofThreads; ++iThread) {
      fvRecoTracks[iThread].SetName(std::string("TrackFinder::fvRecoTracks_") + std::to_string(iThread));
      fvRecoHitIndices[iThread].SetName(std::string("TrackFinder::fvRecoHitIndices_") + std::to_string(iThread));
    }
  }
 
  // -------------------------------------------------------------------------------------------------------------------
  //CBM Level 1 4D Reconstruction
  //Finds tracks using the Cellular Automaton algorithm
  //
  TrackFinder::Output_t TrackFinder::FindTracks(const InputData& input, TimesliceHeader& tsHeader)
  {
    Output_t output;
    Vector<Track>& recoTracks        = output.first;   //reconstructed tracks
    Vector<ca::HitIndex_t>& recoHits = output.second;  //packed hits of reconstructed tracks
 
    if (input.GetNhits() < 1) {
      LOG(warn) << "No hits were passed to the ca::TrackFinder. Stopping the routine";
      return output;
    }
 
    //
    // The main CaTrackFinder routine
    // It splits the input data into sub-timeslices
    // and runs the track finder over the sub-slices
    //
    fMonitorData.StartTimer(ETimer::Tracking);
    fMonitorData.StartTimer(ETimer::PrepareTimeslice);
    fMonitorData.IncrementCounter(ECounter::TrackingCall);
    fMonitorData.IncrementCounter(ECounter::RecoHit, input.GetNhits());
 
    auto timerStart = std::chrono::high_resolution_clock::now();
 
    auto& wDataThread0 = fvWData[0];  // NOTE: Thread 0 must be always defined
 
    // ----- Reset data arrays -----------------------------------------------------------------------------------------
 
    wDataThread0.HitKeyFlags().reset(input.GetNhitKeys(), 0);
 
    fHitTimeInfo.reset(input.GetNhits());
 
    // TODO: move these values to Parameters namespace (S.Zharko)
 
    // length of sub-TS
    const fscal minProtonMomentum = 0.1;
    const fscal preFactor         = sqrt(1. + ProtonMassD * ProtonMassD / (minProtonMomentum * minProtonMomentum));
    const fscal targX             = fParameters.GetTargetPositionX()[0];
    const fscal targY             = fParameters.GetTargetPositionY()[0];
    const fscal targZ             = fParameters.GetTargetPositionZ()[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 (fHitTimeInfo array)
    for (int iStream = 0; iStream < input.GetNdataStreams(); ++iStream) {
 
      fscal maxTimeBeforeHit = std::numeric_limits<fscal>::lowest();
      const int nStreamHits  = input.GetStreamNhits(iStream);
      fStatNhitsTotal += nStreamHits;
 
      for (int ih = 0; ih < nStreamHits; ++ih) {
 
        ca::HitIndex_t caHitId      = input.GetStreamStartIndex(iStream) + ih;
        const ca::Hit& h            = input.GetHit(caHitId);
        const ca::Station<fvec>& st = fParameters.GetStation(h.Station());
        const fscal dx              = h.X() - targX;
        const fscal dy              = h.Y() - targY;
        const fscal dz              = h.Z() - targZ;
        const fscal l               = sqrt(dx * dx + dy * dy + dz * dz);
        const fscal timeOfFlightMin = l * SpeedOfLightInv;
        const fscal timeOfFlightMax = 1.5 * l * preFactor * SpeedOfLightInvD;
        const fscal dt              = h.RangeT();
        // TODO: Is it possible, that the proton mass selection affects the search of heavier particles?
 
        CaHitTimeInfo& info = fHitTimeInfo[caHitId];
        info.fEventTimeMin  = st.timeInfo ? (h.T() - dt - timeOfFlightMax) : -1.e10;
        info.fEventTimeMax  = st.timeInfo ? (h.T() + dt - timeOfFlightMin) : 1.e10;
 
        // NOTE: if not a MT part, use wDataThread0.IsHitKeyUsed, it will be later copied to other threads
        if (info.fEventTimeMin > 500.e6 || info.fEventTimeMax < -500.) {  // cut hits with bogus start time > 500 ms
          wDataThread0.IsHitKeyUsed(h.FrontKey()) = 1;
          wDataThread0.IsHitKeyUsed(h.BackKey())  = 1;
          LOG(error) << "CATrackFinder: skip bogus hit " << h.ToString();
          continue;
        }
        maxTimeBeforeHit       = std::max(maxTimeBeforeHit, info.fEventTimeMax);
        info.fMaxTimeBeforeHit = maxTimeBeforeHit;
        fStatTsStart           = std::min(fStatTsStart, info.fEventTimeMax);
        fStatTsEnd             = std::max(fStatTsEnd, info.fEventTimeMin);
      }
 
      fscal minTimeAfterHit = std::numeric_limits<fscal>::max();
      // loop in the reverse order to fill CaHitTimeInfo::fMinTimeAfterHit fields
 
      for (int ih = nStreamHits - 1; ih >= 0; --ih) {
        ca::HitIndex_t caHitId = input.GetStreamStartIndex(iStream) + ih;
        const ca::Hit& h       = input.GetHit(caHitId);
        if (wDataThread0.IsHitKeyUsed(h.FrontKey()) || wDataThread0.IsHitKeyUsed(h.BackKey())) {
          continue;
        }  // the hit is skipped
        CaHitTimeInfo& info   = fHitTimeInfo[caHitId];
        minTimeAfterHit       = std::min(minTimeAfterHit, info.fEventTimeMin);
        info.fMinTimeAfterHit = minTimeAfterHit;
      }
 
      if (0) {
        static int tmp = 0;
        if (tmp < 10000) {
          tmp++;
          LOG(warning) << "\n\n stream " << iStream << " hits " << nStreamHits << "\n\n";
          for (int ih = 0; (ih < nStreamHits) && (tmp < 10000); ++ih) {
            ca::HitIndex_t caHitId = input.GetStreamStartIndex(iStream) + ih;
            const ca::Hit& h       = input.GetHit(caHitId);
            if (wDataThread0.IsHitKeyUsed(h.FrontKey()) || wDataThread0.IsHitKeyUsed(h.BackKey())) {
              continue;
            }  // the hit is skipped
            CaHitTimeInfo& info = fHitTimeInfo[caHitId];
            if (h.Station() < 4) {
              tmp++;
              LOG(warning) << " hit sta " << h.Station() << " stream " << iStream << " time " << h.T() << " event time "
                           << info.fEventTimeMin << " .. " << info.fEventTimeMax << " max time before hit "
                           << info.fMaxTimeBeforeHit << " min time after hit " << info.fMinTimeAfterHit;
            }
          }
        }
      }
    }
    // all hits belong to one sub-timeslice; 1 s is the maximal length of the TS
    fStatTsEnd = std::clamp(fStatTsEnd, fStatTsStart, fStatTsStart + 1.e9f);
 
    LOG(debug) << "CA tracker process time slice " << fStatTsStart * 1.e-6 << " -- " << fStatTsEnd * 1.e-6
               << " [ms] with " << fStatNhitsTotal << " hits";
 
    int nWindows = static_cast<int>((fStatTsEnd - fStatTsStart) / fWindowLength) + 1;
    if (nWindows < 1) {  // Situation, when fStatTsEnd == fStatTsStart
      nWindows = 1;
    }
 
    // int nWindowsThread = nWindows / fNofThreads;
    // LOG(info) << "CA: estimated number of time windows: " << nWindows;
 
    std::vector<std::pair<fscal, fscal>> vWindowRangeThread(fNofThreads);
    {  // Estimation of number of hits in time windows
      //Timer time;
      //time.Start();
      const HitIndex_t nHitsTot = input.GetNhits();
      const int nSt             = fParameters.GetNstationsActive();
 
      // Count number of hits per window and station
      std::vector<HitIndex_t> nHitsWindowSta(nWindows * nSt, 0);
 
      for (HitIndex_t iHit = 0; iHit < nHitsTot; ++iHit) {
        const auto& hit  = input.GetHit(iHit);
        const auto& info = fHitTimeInfo[iHit];
        int iWindow      = static_cast<int>((info.fEventTimeMin - fStatTsStart) / fWindowLength);
        if (iWindow < 0) {
          iWindow = 0;
        }
        if (iWindow >= nWindows) {
          LOG(error) << "ca: Hit out of range: iHit = " << iHit << ", min. event time = " << info.fEventTimeMin * 1.e-6
                     << " ms, window = " << iWindow;
          continue;
        }
        ++nHitsWindowSta[hit.Station() + iWindow * nSt];
      }
 
      // Remove hits from the "monster events"
      if (ca::TrackingMode::kMcbm == fTrackingMode) {
        const 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);
      HitIndex_t nHitsCollected = 0;
      for (auto it = nHitsWindowSta.begin(); it != nHitsWindowSta.end(); it += nSt) {
        nHitsCollected = nHitsWindow.emplace_back(std::accumulate(it, it + nSt, nHitsCollected));
      }
 
      // Get time range for threads
      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);
        const size_t iWbegin               = std::distance(nHitsWindow.begin(), windowIt) + 1;
        vWindowRangeThread[iTh].first      = fStatTsStart + iWbegin * fWindowLength;
        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 * fWindowLength;
    //  vWindowEndThread[iThread] =
    //    vWindowStartThread[iThread] + nWindowsThread * fWindowLength;
    //}
 
    for (int iThread = 0; iThread < fNofThreads; ++iThread) {
      auto& entry  = vWindowRangeThread[iThread];
      double start = entry.first * 1.e-6;
      double end   = entry.second * 1.e-6;
      LOG(debug) << "Thread: " << iThread << " from " << start << " ms  to " << end << " ms (delta = " << end - start
                 << " ms)";
    }
 
    // Statistics for monitoring
    std::vector<int> vStatNwindows(fNofThreads), vStatNhitsProcessed(fNofThreads);
 
    fMonitorData.StopTimer(ETimer::PrepareTimeslice);
    // Save tracks
    if (fNofThreads == 1) {
      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]);
      fMonitorData.StopTimer(ETimer::StoreTracksFinal);
    }
    else {
      std::vector<std::thread> vThreadList;
      vThreadList.reserve(fNofThreads);
      for (int iTh = 0; iTh < fNofThreads; ++iTh) {
        vThreadList.emplace_back(&TrackFinder::FindTracksThread, this, std::ref(input), iTh,
                                 std::ref(vWindowRangeThread[iTh]), std::ref(vStatNwindows[iTh]),
                                 std::ref(vStatNhitsProcessed[iTh]));
      }
      for (auto& th : vThreadList) {
        if (th.joinable()) {
          th.join();
        }
      }
      fMonitorData.StartTimer(ETimer::StoreTracksFinal);
      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);
      recoTracks.reserve(nRecoTracks);
      recoHits.reserve(nRecoHits);
      for (int iTh = 0; iTh < fNofThreads; ++iTh) {
        recoTracks.insert(recoTracks.end(), fvRecoTracks[iTh].begin(), fvRecoTracks[iTh].end());
        recoHits.insert(recoHits.end(), fvRecoHitIndices[iTh].begin(), fvRecoHitIndices[iTh].end());
      }
      fMonitorData.StopTimer(ETimer::StoreTracksFinal);
    }
 
    fMonitorData.IncrementCounter(ECounter::RecoTrack, recoTracks.size());
    fMonitorData.IncrementCounter(ECounter::RecoHitUsed, recoHits.size());
 
    auto timerEnd = std::chrono::high_resolution_clock::now();
    fCaRecoTime   = (double) (std::chrono::duration<double>(timerEnd - timerStart).count());
 
    // Add thread monitors to the main monitor
    for (auto& monitor : fvMonitorDataThread) {
      fMonitorData.AddMonitorData(monitor, true);
      //fMonitorData.AddMonitorData(monitor);
      monitor.Reset();
    }
 
    const int statNhitsProcessedTotal = std::accumulate(vStatNhitsProcessed.begin(), vStatNhitsProcessed.end(), 0);
    const int statNwindowsTotal       = std::accumulate(vStatNwindows.begin(), vStatNwindows.end(), 0);
 
    // Filling TS headear
    tsHeader.Start() = fStatTsStart;
    tsHeader.End()   = fStatTsEnd;
 
    fMonitorData.StopTimer(ETimer::Tracking);
 
    LOG(debug) << "CA tracker: time slice finished. Reconstructed " << recoTracks.size() << " tracks with "
               << recoHits.size() << " hits. Processed " << statNhitsProcessedTotal << " hits in " << statNwindowsTotal
               << " time windows. Reco time " << fCaRecoTime / 1.e9 << " s";
    return output;
  }
 
  // -------------------------------------------------------------------------------------------------------------------
  //
  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());
    Timer timer;
    timer.Start();
 
    auto& monitor = fvMonitorDataThread[iThread];
    monitor.StartTimer(ETimer::TrackingThread);
    monitor.StartTimer(ETimer::PrepareThread);
 
    // Init vectors
    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;
      tracks.clear();
      tracks.reserve(nTracksExpected / fNofThreads);
      hitIndices.clear();
      hitIndices.reserve(nHitsExpected / fNofThreads);
      if (iThread != 0) {
        wData.HitKeyFlags() = fvWData[0].HitKeyFlags();
      }
      for (int iS = 0; iS < nStations; ++iS) {
        wData.TsHitIndices(iS).clear();
        wData.TsHitIndices(iS).reserve(nHitsTot);
      }
    }
 
    // 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 (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 < windowRange.first) {
          range.first = caHitId + 1;
        }
        if (info.fMinTimeAfterHit > windowRange.second) {
          range.second = caHitId;
        }
      }
    }
 
    int statLastLogTimeChunk = -1;
 
    // Track finder algorithm for the time window
    ca::TrackFinderWindow trackFinderWindow(fParameters, fDefaultMass, fTrackingMode, monitor);
    trackFinderWindow.InitTimeslice(input.GetNhitKeys());
 
    monitor.StopTimer(ETimer::PrepareThread);
 
    while (true) {
      monitor.IncrementCounter(ECounter::SubTS);
      // select the sub-slice hits
      for (int iS = 0; iS < fParameters.GetNstationsActive(); ++iS) {
        wData.TsHitIndices(iS).clear();
      }
      bool areUntouchedDataLeft = false;  // is the whole TS processed
 
      // TODO: SG: skip empty regions and start the subslice with the earliest hit
 
      statNwindows++;
      //out << statNwindows << ' ';
 
      monitor.StartTimer(ETimer::PrepareWindow);
 
      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;
          }
          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;
            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) {
            range.first = caHitId + 1;  // this hit and all hits before are before the overlap
          }
        }
      }
 
      //out << statNwindowHits << ' ';
      //if (statNwindowHits == 0) {  // Empty window
      //  monitor.StopTimer(ETimer::PrepareWindow);
      //  out << 0 << ' ' << 0 << ' ' << 0 << '\n';
      //  continue;
      //}
 
      if (ca::TrackingMode::kMcbm == fTrackingMode) {
        // cut at 50 hits per station per 1 us.
        int maxStationHits = (int) (50 * fWindowLength / 1.e3);
        for (int ista = 0; ista < fParameters.GetNstationsActive(); ++ista) {
          int nHitsSta = static_cast<int>(wData.TsHitIndices(ista).size());
          if (nHitsSta > maxStationHits) {
            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) ((windowRange.first - fStatTsStart) / 10.e6);
        if (!areUntouchedDataLeft || currentChunk > statLastLogTimeChunk) {
          statLastLogTimeChunk = currentChunk;
          double dataRead = 100. * (windowRange.first + fWindowLength - fStatTsStart) / (fStatTsEnd - fStatTsStart);
          if (dataRead > 100.) {
            dataRead = 100.;
          }
          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 " << tracks.size() << " tracks on thread #" << iThread;
        }
      }
 
      //out << statNwindowHits << ' ';
      monitor.StopTimer(ETimer::PrepareWindow);
 
      //Timer trackingInWindow;  //DBG
      //trackingInWindow.Start();
      monitor.StartTimer(ETimer::TrackingWindow);
      trackFinderWindow.CaTrackFinderSlice(input, wData);
      monitor.StopTimer(ETimer::TrackingWindow);
      //trackingInWindow.Stop();
      //out << trackingInWindow.GetTotalMs() << ' ';
 
      // save reconstructed tracks with no hits in the overlap region
      //if (windowRange.first > 13.23e6 && windowRange.first < 13.26e6) {
      windowRange.first += fWindowLength;
      // 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';
 
      monitor.StartTimer(ETimer::StoreTracksWindow);
      auto trackFirstHit = wData.RecoHitIndices().begin();
 
      for (const auto& track : wData.RecoTracks()) {
 
        const bool isTrackCompletelyInOverlap =
          std::all_of(trackFirstHit, trackFirstHit + track.fNofHits, [&](int caHitId) {
            CaHitTimeInfo& info = fHitTimeInfo[caHitId];
            return info.fEventTimeMax >= windowRange.first;
          });
 
        // Don't save tracks from the overlap region, since they might have additional hits in the next subslice.
        // Don't reject tracks in the overlap when no more data are left
        const bool useFlag = !isTrackCompletelyInOverlap || !areUntouchedDataLeft;
 
        for (int i = 0; i < track.fNofHits; i++) {
          const int caHitId                = *(trackFirstHit + i);
          const auto& h                    = input.GetHit(caHitId);
          wData.IsHitKeyUsed(h.FrontKey()) = static_cast<int>(useFlag);
          wData.IsHitKeyUsed(h.BackKey())  = static_cast<int>(useFlag);
 
          if (useFlag) {
            hitIndices.push_back(caHitId);
          }
        }
        if (useFlag) {
          tracks.push_back(track);
        }
 
        trackFirstHit += track.fNofHits;
      }  // sub-timeslice tracks
      monitor.StopTimer(ETimer::StoreTracksWindow);
 
      if (windowRange.first > windowRange.second) {
        break;
      }
      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: " << hitIndices.size() << ')';
  }
}  // namespace cbm::algo::ca