algo/ca/core/tracking/CaTrackFinder.h

+/* Copyright (C) 2022 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+   SPDX-License-Identifier: GPL-3.0-only
+   Authors: Sergei Zharko [committer], Maksym Zyzak */
+
+/// \file    CaTrackFinder.h
+
+#pragma once  // include this header only once per compilation unit
+
+#include "CaSimd.h"
+#include "CaTimesliceHeader.h"
+#include "CaTrackExtender.h"
+#include "CaTrackFinderWindow.h"
+#include "CaTrackFitter.h"
+#include "CaVector.h"
+#include "KfTrackParam.h"
+
+#include <vector>
+
+namespace cbm::algo::ca
+{
+
+  class Track;
+  class Triplet;
+  class Framework;
+
+  /// Class implements a clones merger algorithm for the CA track finder
+  ///
+  class TrackFinder {
+   public:
+    typedef std::pair<Vector<Track>, Vector<ca::HitIndex_t>> Output_t;
+
+    /// Default constructora
+    TrackFinder(const ca::Parameters<fvec>& pars, const fscal mass, const ca::TrackingMode& mode,
+                TrackingMonitorData& monitorData, int nThreads, double& recoTime);
+    /// Destructor
+    ~TrackFinder() = default;
+
+    /// Copy constructor
+    TrackFinder(const TrackFinder&) = delete;
+
+    /// Move constructor
+    TrackFinder(TrackFinder&&) = delete;
+
+    /// Copy assignment operator
+    TrackFinder& operator=(const TrackFinder&) = delete;
+
+    /// Move assignment operator
+    TrackFinder& operator=(TrackFinder&&) = delete;
+
+    Output_t FindTracks(const InputData& input, TimesliceHeader& tsHeader);
+
+    const auto& GetRecoTracksContainer(int iThread) const { return fvRecoTracks[iThread]; }
+    const auto& GetRecoHitIndicesContainer(int iThread) const { return fvRecoHitIndices[iThread]; }
+    TrackingMode GetTrackingMode() const { return fTrackingMode; }
+    const std::vector<ca::WindowData>& GetWData() const { return fvWData; }
+
+   private:
+    // -------------------------------
+    // Private methods
+    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;
+
+    // -------------------------------
+    // Data members
+
+    Vector<CaHitTimeInfo> fHitTimeInfo;
+
+    const Parameters<fvec>& fParameters;            ///< Object of Framework parameters class
+    fscal fDefaultMass{constants::phys::MuonMass};  ///< mass of the propagated particle [GeV/c2]
+    ca::TrackingMode fTrackingMode;
+
+    TrackingMonitorData& fMonitorData;                     ///< Tracking monitor data (statistics per call)
+    std::vector<TrackingMonitorData> fvMonitorDataThread;  ///< Tracking monitor data per thread
+
+    std::vector<ca::WindowData> fvWData;  ///< Intrnal data processed in a time-window
+
+    int fNofThreads;      ///< Number of threads to execute the track-finder
+    double& fCaRecoTime;  // time of the track finder + fitter
+
+    std::vector<Vector<Track>> fvRecoTracks;           ///< reconstructed tracks
+    std::vector<Vector<HitIndex_t>> fvRecoHitIndices;  ///< packed hits of reconstructed tracks
+
+    fscal fWindowLength  = 0.;   ///< Time window length [ns]
+    fscal fStatTsStart  = 0.;
+    fscal fStatTsEnd    = 0.;
+    int fStatNhitsTotal = 0;
+  };
+
+}  // namespace cbm::algo::ca

algo/ca/core/tracking/CaTrackFinderWindow.cxx

+/* Copyright (C) 2009-2021 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+   SPDX-License-Identifier: GPL-3.0-only
+   Authors: Valentina Akishina, Ivan Kisel, Sergey Gorbunov [committer], Igor Kulakov, 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 "CaTrackFinderWindow.h"
+
+#include "AlgoFairloggerCompat.h"
+#include "CaBranch.h"
+#include "CaFramework.h"
+#include "CaGrid.h"
+#include "CaGridEntry.h"
+#include "CaTrack.h"
+#include "CaTripletConstructor.h"
+
+#include <algorithm>
+#include <cstdio>
+#include <fstream>
+#include <iostream>
+#include <list>
+#include <map>
+#include <sstream>
+
+// #include "CaToolsDebugger.h"
+
+namespace cbm::algo::ca
+{
+  // -------------------------------------------------------------------------------------------------------------------
+  TrackFinderWindow::TrackFinderWindow(const ca::Parameters<fvec>& pars, const fscal mass, const ca::TrackingMode& mode,
+                                       ca::TrackingMonitorData& monitorData)
+    : fParameters(pars)
+    , fDefaultMass(mass)
+    , fTrackingMode(mode)
+    , frMonitorData(monitorData)
+    , fTrackExtender(pars, mass)
+    , fCloneMerger(pars, mass)
+    , fTrackFitter(pars, mass, mode)
+  {
+  }
+
+  // -------------------------------------------------------------------------------------------------------------------
+  bool TrackFinderWindow::checkTripletMatch(const ca::Triplet& l, const ca::Triplet& r, fscal& dchi2,
+                                            WindowData& wData) const
+  {
+    dchi2 = 1.e20;
+
+    if (r.GetMHit() != l.GetRHit()) return false;
+    if (r.GetLHit() != l.GetMHit()) return false;
+
+    if (r.GetMSta() != l.GetRSta()) return false;
+    if (r.GetLSta() != l.GetMSta()) return false;
+
+    const fscal tripletLinkChi2 = wData.CurrentIteration()->GetTripletLinkChi2();
+    if (r.IsMomentumFitted()) {
+      assert(l.IsMomentumFitted());
+
+      fscal dqp = l.GetQp() - r.GetQp();
+      fscal Cqp = l.GetCqp() + r.GetCqp();
+
+      if (!std::isfinite(dqp)) return false;
+      if (!std::isfinite(Cqp)) return false;
+
+      if (dqp * dqp > tripletLinkChi2 * Cqp) {
+        return false;  // bad neighbour // CHECKME why do we need recheck it?? (it really change result)
+      }
+      dchi2 = dqp * dqp / Cqp;
+    }
+    else {
+      fscal dtx = l.GetTx() - r.GetTx();
+      fscal Ctx = l.GetCtx() + r.GetCtx();
+
+      fscal dty = l.GetTy() - r.GetTy();
+      fscal Cty = l.GetCty() + r.GetCty();
+
+      // it shouldn't happen, but happens sometimes
+
+      if (!std::isfinite(dtx)) return false;
+      if (!std::isfinite(dty)) return false;
+      if (!std::isfinite(Ctx)) return false;
+      if (!std::isfinite(Cty)) return false;
+
+      if (dty * dty > tripletLinkChi2 * Cty) return false;
+      if (dtx * dtx > tripletLinkChi2 * Ctx) return false;
+
+      //dchi2 = 0.5f * (dtx * dtx / Ctx + dty * dty / Cty);
+      dchi2 = 0.;
+    }
+
+    if (!std::isfinite(dchi2)) return false;
+
+    return true;
+  }
+
+
+  // **************************************************************************************************
+  // *                                                                                                *
+  // *                            ------ CATrackFinder procedure ------                               *
+  // *                                                                                                *
+  // **************************************************************************************************
+
+  // -------------------------------------------------------------------------------------------------------------------
+  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;
+    for (int iS = 0; iS < fParameters.GetNstationsActive(); iS++) {
+      wData.HitStartIndexOnStation(iS) = nHits;
+      wData.NofHitsOnStation(iS)       = wData.TsHitIndices(iS).size();
+      nHits += wData.NofHitsOnStation(iS);
+    }
+    wData.HitStartIndexOnStation(fParameters.GetNstationsActive()) = nHits;
+    wData.ResetHitData(nHits);
+
+    for (int iS = 0; iS < fParameters.GetNstationsActive(); iS++) {
+      int iFstHit = wData.HitStartIndexOnStation(iS);
+      for (ca::HitIndex_t ih = 0; ih < wData.TsHitIndices(iS).size(); ++ih) {
+        ca::Hit h = hits[wData.TsHitIndex(iS, ih)];  /// D.S. 29.7.24: There is a copy operation here. Can be avoided?
+        h.SetId(wData.TsHitIndex(iS, ih));
+        wData.Hit(iFstHit + ih) = h;
+      }
+    }
+
+    if constexpr (fDebug) {
+      LOG(info) << "===== Sliding Window hits: ";
+      for (int i = 0; i < nHits; ++i) {
+        LOG(info) << "    " << wData.Hit(i).ToString();
+      }
+      LOG(info) << "===== ";
+    }
+
+    wData.RecoTracks().clear();
+    wData.RecoTracks().reserve(2 * nHits / fParameters.GetNstationsActive());
+
+    wData.RecoHitIndices().clear();
+    wData.RecoHitIndices().reserve(2 * nHits);
+  }
+
+  // -------------------------------------------------------------------------------------------------------------------
+  void TrackFinderWindow::PrepareGrid(const Vector<Hit>& hits, WindowData& wData)
+  {
+    for (int iS = 0; iS < fParameters.GetNstationsActive(); ++iS) {
+
+      fscal lasttime  = std::numeric_limits<fscal>::infinity();
+      fscal starttime = -std::numeric_limits<fscal>::infinity();
+      fscal gridMinX  = -0.1;
+      fscal gridMaxX  = 0.1;
+      fscal gridMinY  = -0.1;
+      fscal gridMaxY  = 0.1;
+
+      for (ca::HitIndex_t ih = 0; ih < wData.TsHitIndices(iS).size(); ++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());
+
+        const fscal time = h.T();
+        assert(std::isfinite(time));
+        lasttime  = std::min(lasttime, time);
+        starttime = std::max(starttime, time);
+      }
+
+      // TODO: changing the grid also changes the result. Investigate why it happens.
+      // TODO: find the optimal grid size
+
+      const int nSliceHits = wData.TsHitIndices(iS).size();
+      const fscal sizeY    = gridMaxY - gridMinY;
+      const fscal sizeX    = gridMaxX - gridMinX;
+      const int nBins2D    = 1 + nSliceHits;
+
+      // TODO: SG: the coefficients should be removed
+      const fscal scale    = fParameters.GetStation(iS).GetZ<fscal>() - fParameters.GetTargetPositionZ()[0];
+      const fscal maxScale = 0.3 * scale;
+      const fscal minScale = 0.01 * scale;
+
+      fscal yStep = 0.3 * sizeY / sqrt(nBins2D);
+      fscal xStep = 0.8 * sizeX / sqrt(nBins2D);
+      yStep       = std::clamp(yStep, minScale, maxScale);
+      xStep       = std::clamp(xStep, minScale, maxScale);
+
+      auto& grid = wData.Grid(iS);
+      grid.BuildBins(gridMinX, gridMaxX, gridMinY, gridMaxY, xStep, yStep);
+      /* clang-format off */
+      grid.StoreHits(wData.Hits(), 
+                     wData.HitStartIndexOnStation(iS), 
+                     wData.NofHitsOnStation(iS),
+                     wData.HitKeyFlags());
+      /* clang-format on */
+    }
+  }
+
+  // -------------------------------------------------------------------------------------------------------------------
+  void TrackFinderWindow::PrepareCAIteration(const ca::Iteration& caIteration, WindowData& wData, const bool isFirst)
+  {
+    wData.SetCurrentIteration(&caIteration);
+
+    // 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?
+
+    // --- 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);
+  }
+
+  // -------------------------------------------------------------------------------------------------------------------
+  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, *, *}
+
+    ca::TripletConstructor constructor(fParameters, wData, fDefaultMass, fTrackingMode);
+
+    // 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));
+      }
+    }
+
+    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());
+        }
+        vHitFirstTriplet[ihitl] = PackTripletId(istal, oldSize);
+        vHitNofTriplets[ihitl]  = vTriplets[istal].size() - oldSize;
+      }
+    }  // istal
+  }
+
+  // -------------------------------------------------------------------------------------------------------------------
+  void TrackFinderWindow::SearchNeighbors(WindowData& wData)
+  {
+    auto& [vHitFirstTriplet, vHitNofTriplets, vTriplets] = fTripletData;
+
+    for (int istal = fParameters.GetNstationsActive() - 2; istal >= wData.CurrentIteration()->GetFirstStationIndex();
+         istal--) {
+      // start with downstream chambers
+
+      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);
+
+        if (nNeighbours > 0) {
+          assert((int) neighStation >= istal + 1
+                 && (int) neighStation <= istal + 1 + wData.CurrentIteration()->GetMaxStationGap());
+        }
+
+        unsigned char level = 0;
+
+        for (unsigned int iN = 0; iN < nNeighbours; ++iN, ++neighTriplet, ++neighLocation) {
+
+          ca::Triplet& neighbour = vTriplets[neighStation][neighTriplet];
+
+          fscal dchi2 = 0.;
+          if (!checkTripletMatch(tr, neighbour, dchi2, wData)) continue;
+
+          if (tr.GetFNeighbour() == 0) {
+            tr.SetFNeighbour(neighLocation);
+          }
+          tr.SetNNeighbours(neighLocation - tr.GetFNeighbour() + 1);
+
+          level = std::max(level, static_cast<unsigned char>(neighbour.GetLevel() + 1));
+        }
+        tr.SetLevel(level);
+      }
+      frMonitorData.IncrementCounter(ECounter::Triplet, vTriplets[istal].size());
+    }
+  }
+
+  // -------------------------------------------------------------------------------------------------------------------
+  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);
+    }
+  }
+
+  // -------------------------------------------------------------------------------------------------------------------
+  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;
+
+    const unsigned char min_best_l =
+      (firstTripletLevel > min_level) ? firstTripletLevel + 2 : min_level + 3;  // loose maximum
+
+    // Uses persistent field to save memory allocations.
+    // fNewTr is only used here!
+    ca::Branch(&new_tr)[constants::size::MaxNstations] = fNewTr;
+
+    fvTrackCandidates.clear();
+    fvTrackCandidates.reserve(wData.Hits().size() / 10);
+
+    for (const auto& h : wData.Hits()) {
+      fvHitKeyToTrack[h.FrontKey()] = -1;
+      fvHitKeyToTrack[h.BackKey()]  = -1;
+    }
+
+    //== Loop over triplets with the required level, find and store track candidates
+    for (int istaF = wData.CurrentIteration()->GetFirstStationIndex();
+         istaF <= fParameters.GetNstationsActive() - 3 - firstTripletLevel; ++istaF) {
+
+      if (--nlevel == 0) break;  //TODO: SG: this is not needed
+
+      for (ca::Triplet& first_trip : vTriplets[istaF]) {
+
+        const auto& fstTripLHit = wData.Hit(first_trip.GetLHit());
+        if (wData.IsHitKeyUsed(fstTripLHit.FrontKey()) || wData.IsHitKeyUsed(fstTripLHit.BackKey())) {
+          continue;
+        }
+
+        //  skip track candidates that are too short
+
+        int minNhits = wData.CurrentIteration()->GetMinNhits();
+
+        if (fstTripLHit.Station() == 0) {
+          minNhits = wData.CurrentIteration()->GetMinNhitsStation0();
+        }
+        if (wData.CurrentIteration()->GetTrackFromTripletsFlag()) {
+          minNhits = 0;
+        }
+
+        if (3 + first_trip.GetLevel() < minNhits) {
+          continue;
+        }
+
+        // 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;
+        }
+
+        ca::Branch curr_tr;
+        curr_tr.AddHit(first_trip.GetLHit());
+        curr_tr.SetChi2(first_trip.GetChi2());
+
+        ca::Branch best_tr = curr_tr;
+
+        /// 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);
+
+        if (best_tr.NofHits() < firstTripletLevel + 2) continue;  // loose maximum one hit
+
+        if (best_tr.NofHits() < min_level + 3) continue;  // should find all hits for min_level
+
+        if (best_tr.NofHits() < minNhits) {
+          continue;
+        }
+
+        int ndf = best_tr.NofHits() * 2 - 5;
+
+        // TODO: automatize the NDF calculation
+
+        if (ca::TrackingMode::kGlobal == fTrackingMode || ca::TrackingMode::kMcbm == fTrackingMode) {
+          ndf = best_tr.NofHits() * 2 - 4;
+        }
+
+        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;
+          }
+        }
+        // NOTE: Temporarely disable the warnings, one has to provide more specific coefficient in the capacity
+        //                   reservation
+        fvTrackCandidates.push_back_no_warning(best_tr);
+        ca::Branch& tr = fvTrackCandidates.back();
+        tr.SetStation(istaF);
+        tr.SetId(fvTrackCandidates.size() - 1);
+        // Mark the candidate as dead. To became alive it should first pass the competition in DoCompetitionLoop
+        tr.SetAlive(false);
+        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
+  }
+
+  // -------------------------------------------------------------------------------------------------------------------
+  void TrackFinderWindow::DoCompetitionLoop(const WindowData& wData)
+  {
+
+    // look at the dead track candidates in fvTrackCandidates pool; select the best ones and make them alive
+
+    for (int iComp = 0; (iComp < 100); ++iComp) {
+
+      bool 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)) {
+                strip = tr.Id();
+              }
+              else {
+                return false;
+              }
+            }
+            else {
+              strip = tr.Id();
+            }
+            return true;
+          };
+
+          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 (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 = wData.Hit(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 = wData.Hit(hitId);
+            if (fvHitKeyToTrack[h.FrontKey()] == tr.Id()) {
+              fvHitKeyToTrack[h.FrontKey()] = -1;
+            }
+            if (fvHitKeyToTrack[h.BackKey()] == tr.Id()) {
+              fvHitKeyToTrack[h.BackKey()] = -1;
+            }
+          }
+        }
+        else {
+          repeatCompetition = true;
+        }
+      }  // itrack
+
+      if (!repeatCompetition) break;
+    }  // competitions
+  }
+
+  // -------------------------------------------------------------------------------------------------------------------
+  void TrackFinderWindow::SelectTracks(WindowData& wData)
+  {
+    for (Tindex iCandidate = 0; iCandidate < (Tindex) fvTrackCandidates.size(); ++iCandidate) {
+      ca::Branch& tr = fvTrackCandidates[iCandidate];
+
+      if constexpr (fDebug) {
+        LOG(info) << "iter " << wData.CurrentIteration()->GetName() << ", track candidate " << iCandidate
+                  << ": alive = " << tr.IsAlive();
+      }
+      if (!tr.IsAlive()) continue;
+
+      if (wData.CurrentIteration()->GetExtendTracksFlag()) {
+        if (tr.NofHits() < fParameters.GetNstationsActive()) {
+          fTrackExtender.ExtendBranch(tr, wData);
+        }
+      }
+
+      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         *
+       *                                                              *
+       *     I. Kisel                                    06.03.05     *
+       *     I.Kulakov                                    2012        *
+       *                                                              *
+       ****************************************************************/
+
+  // -------------------------------------------------------------------------------------------------------------------
+  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,
+                                      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
+  {
+
+    if (curr_trip->GetLevel() == 0)  // the end of the track -> check and store
+    {
+      // -- finish with current track
+      // add rest of hits
+      const auto& hitM = wData.Hit(curr_trip->GetMHit());
+      const auto& hitR = wData.Hit(curr_trip->GetRHit());
+
+      if (!(wData.IsHitKeyUsed(hitM.FrontKey()) || wData.IsHitKeyUsed(hitM.BackKey()))) {
+        curr_tr.AddHit(curr_trip->GetMHit());
+      }
+      if (!(wData.IsHitKeyUsed(hitR.FrontKey()) || wData.IsHitKeyUsed(hitR.BackKey()))) {
+        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
+      // TODO: SZh 21.08.2023: Replace hard-coded value with a parameter
+      int ndf = curr_tr.NofHits() * 2 - 5;
+
+      if (ca::TrackingMode::kGlobal == fTrackingMode || ca::TrackingMode::kMcbm == fTrackingMode) {
+        ndf = curr_tr.NofHits() * 2 - 4;
+      }
+      if (curr_tr.Chi2() > wData.CurrentIteration()->GetTrackChi2Cut() * ndf) {
+        return;
+      }
+
+      // -- select the best
+      if ((curr_tr.NofHits() > best_tr.NofHits())
+          || ((curr_tr.NofHits() == best_tr.NofHits()) && (curr_tr.Chi2() < best_tr.Chi2()))) {
+        best_tr = curr_tr;
+      }
+      return;
+    }
+    else  //MEANS level ! = 0
+    {
+      int N_neighbour = (curr_trip->GetNNeighbours());
+
+      for (Tindex in = 0; in < N_neighbour; in++) {
+
+        unsigned int ID      = curr_trip->GetFNeighbour() + in;
+        unsigned int Station = TripletId2Station(ID);
+        unsigned int Triplet = TripletId2Triplet(ID);
+
+        const ca::Triplet& new_trip = vTriplets[Station][Triplet];
+
+        fscal dchi2 = 0.;
+        if (!checkTripletMatch(*curr_trip, new_trip, dchi2, wData)) continue;
+
+        const auto& hitL = wData.Hit(new_trip.GetLHit());                                 // left hit of new triplet
+        if (wData.IsHitKeyUsed(hitL.FrontKey()) || wData.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()))) {
+            best_tr = curr_tr;
+          }
+        }
+        else {  // hit is not used: add the left hit from the new triplet to the current track
+
+          unsigned char new_L = curr_tr.NofHits() + 1;
+          fscal new_chi2      = curr_tr.Chi2() + dchi2;
+
+          if constexpr (0) {  //SGtrd2d debug!!
+            int mc01 = ca::Framework::GetMcTrackIdForWindowHit(curr_trip->GetLHit());
+            int mc02 = ca::Framework::GetMcTrackIdForWindowHit(curr_trip->GetMHit());
+            int mc03 = ca::Framework::GetMcTrackIdForWindowHit(curr_trip->GetRHit());
+            int mc11 = ca::Framework::GetMcTrackIdForWindowHit(new_trip.GetLHit());
+            int mc12 = ca::Framework::GetMcTrackIdForWindowHit(new_trip.GetMHit());
+            int mc13 = ca::Framework::GetMcTrackIdForWindowHit(new_trip.GetRHit());
+
+            if ((mc01 == mc02) && (mc02 == mc03)) {
+              LOG(info) << " sta " << ista << " mc0 " << mc01 << " " << mc02 << " " << mc03 << " mc1 " << mc11 << " "
+                        << mc12 << " " << mc13 << " chi2 " << curr_tr.Chi2() / (2 * (curr_tr.NofHits() + 2) - 4)
+                        << " new " << new_chi2 / (2 * (new_L + 2) - 4);
+              LOG(info) << "   hits " << curr_trip->GetLHit() << " " << curr_trip->GetMHit() << " "
+                        << curr_trip->GetRHit() << " " << new_trip.GetLHit();
+            }
+          }
+
+          int ndf = 2 * (new_L + 2) - 5;
+
+          if (ca::TrackingMode::kGlobal == fTrackingMode) {  //SGtrd2d!!!
+            ndf = 2 * (new_L + 2) - 4;
+          }
+          else if (ca::TrackingMode::kMcbm == fTrackingMode) {
+            ndf = 2 * (new_L + 2) - 4;
+          }
+          else {
+            ndf = new_L;  // TODO: SG:  2 * (new_L + 2) - 5
+          }
+
+          if (new_chi2 > wData.CurrentIteration()->GetTrackChi2Cut() * ndf) continue;
+
+          // add new hit
+          new_tr[ista] = curr_tr;
+          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, vTriplets);
+        }  // add triplet to track
+      }    // for neighbours
+    }      // level = 0
+  }
+
+}  // namespace cbm::algo::ca

algo/ca/core/tracking/CaTrackFinderWindow.h

+/* Copyright (C) 2022 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+   SPDX-License-Identifier: GPL-3.0-only
+   Authors: Sergei Zharko [committer], Maksym Zyzak */
+
+/// \file    CaTrackFinderWindow.h
+/// \authors S.Zharko <s.zharko@gsi.de> (interface), M.Zyzak (original algorithm)
+/// \brief   A class wrapper over clones merger algorithm for the CA track finder (declaration)
+/// \since   22.07.2022
+
+#pragma once  // include this header only once per compilation unit
+
+#include "CaBranch.h"
+#include "CaCloneMerger.h"
+#include "CaGrid.h"
+#include "CaHit.h"
+#include "CaSimd.h"
+#include "CaTrackExtender.h"
+#include "CaTrackFitter.h"
+#include "CaTrackingMonitor.h"
+#include "CaVector.h"
+#include "CaWindowData.h"
+
+namespace cbm::algo::ca
+{
+
+  class Track;
+  class Triplet;
+  class Framework;
+
+  /// Class implements a clones merger algorithm for the CA track finder
+  ///
+  class TrackFinderWindow {
+   public:
+    /// Default constructor
+    TrackFinderWindow(const ca::Parameters<fvec>& pars, const fscal mass, const ca::TrackingMode& mode,
+                      ca::TrackingMonitorData& monitorData);
+    /// Destructor
+    ~TrackFinderWindow() = default;
+
+    /// Copy constructor
+    TrackFinderWindow(const TrackFinderWindow&) = default;
+
+    /// Move constructor
+    TrackFinderWindow(TrackFinderWindow&&) = default;
+
+    /// Copy assignment operator
+    TrackFinderWindow& operator=(const TrackFinderWindow&) = delete;
+
+    /// Move assignment operator
+    TrackFinderWindow& operator=(TrackFinderWindow&&) = delete;
+
+    void CaTrackFinderSlice(const ca::InputData& input, WindowData& wData);
+
+    /// \note The function initializes global arrays for a given thread
+    void InitTimeslice(size_t nHitKeys) { fvHitKeyToTrack.reset(nHitKeys, -1); }
+
+   private:
+    ///-------------------------------
+    /// 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, TripletArray_t& vTriplets);
+
+    bool checkTripletMatch(const ca::Triplet& l, const ca::Triplet& r, fscal& dchi2, WindowData& wData) const;
+
+    // ** Functions, which pack and unpack indexes of station and triplet **
+
+    // TODO: move to ca::Triplet class (S.Zharko)
+    /// Packs station and triplet indices to an unique triplet ID
+    /// \param  iStation  Index of station in the active stations array
+    /// \param  iTriplet  Index of triplet
+    static unsigned int PackTripletId(unsigned int iStation, unsigned int iTriplet);
+
+    /// Unpacks the triplet ID to its station index
+    /// \param  id  Unique triplet ID
+    static unsigned int TripletId2Station(unsigned int id);
+
+    /// Unpacks the triplet ID to its triplet index
+    /// \param  id  Unique triplet ID
+    static unsigned int TripletId2Triplet(unsigned int id);
+
+
+   private:
+    ///-------------------------------
+    /// Data members
+
+    static constexpr bool fDebug = false;  // print debug info
+
+    const Parameters<fvec>& fParameters;            ///< Object of Framework parameters class
+    fscal fDefaultMass{constants::phys::MuonMass};  ///< mass of the propagated particle [GeV/c2]
+    ca::TrackingMode fTrackingMode;
+
+    TrackingMonitorData& frMonitorData;  ///< Reference to monitor data
+    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;
+  };
+
+  // ********************************************
+  // ** Inline member functions implementation **
+  // ********************************************
+
+  // ---------------------------------------------------------------------------------------------------------------------
+  //
+  [[gnu::always_inline]] inline unsigned int TrackFinderWindow::PackTripletId(unsigned int iStation,
+                                                                              unsigned int iTriplet)
+  {
+    assert(iStation < constants::size::MaxNstations);
+    assert(iTriplet < constants::size::MaxNtriplets);
+    constexpr unsigned int kMoveStation = constants::size::TripletBits;
+    return (iStation << kMoveStation) + iTriplet;
+  }
+
+  // ---------------------------------------------------------------------------------------------------------------------
+  //
+  [[gnu::always_inline]] inline unsigned int TrackFinderWindow::TripletId2Station(unsigned int id)
+  {
+    constexpr unsigned int kMoveStation = constants::size::TripletBits;
+    return id >> kMoveStation;
+  }
+
+  // ---------------------------------------------------------------------------------------------------------------------
+  //
+  [[gnu::always_inline]] inline unsigned int TrackFinderWindow::TripletId2Triplet(unsigned int id)
+  {
+    constexpr unsigned int kTripletMask = (1u << constants::size::TripletBits) - 1u;
+    return id & kTripletMask;
+  }
+
+
+}  // namespace cbm::algo::ca

algo/ca/core/tracking/CaTrackFitter.cxx

+/* Copyright (C) 2010-2020 Frankfurt Institute for Advanced Studies, Goethe-Universitaet Frankfurt, Frankfurt
+   SPDX-License-Identifier: GPL-3.0-only
+   Authors: Ivan Kisel, Sergey Gorbunov, Igor Kulakov [committer], Valentina Akishina, Maksym Zyzak */
+
+#include "CaTrackFitter.h"
+
+#include "CaFramework.h"
+#include "KfTrackKalmanFilter.h"
+#include "KfTrackParam.h"
+
+#include <iostream>
+#include <vector>
+
+
+namespace cbm::algo::ca
+{
+  // -------------------------------------------------------------------------------------------------------------------
+  //
+  TrackFitter::TrackFitter(const ca::Parameters<fvec>& pars, const fscal mass, const ca::TrackingMode& mode)
+    : fParameters(pars)
+    , fSetup(fParameters.GetActiveSetup())
+    , fDefaultMass(mass)
+    , fTrackingMode(mode)
+  {
+  }
+
+  // -------------------------------------------------------------------------------------------------------------------
+  //
+  TrackFitter::~TrackFitter() {}
+
+  // -------------------------------------------------------------------------------------------------------------------
+  //
+  void TrackFitter::FitCaTracks(const ca::InputData& input, WindowData& wData)
+  {
+    //  LOG(info) << " Start CA Track Fitter ";
+    int start_hit = 0;  // for interation in wData.RecoHitIndices()
+
+    //  static kf::FieldValue fldB0, fldB1, fldB2 _fvecalignment;
+    //  static kf::FieldRegion fld _fvecalignment;
+    kf::FieldValue<fvec> fldB0, fldB1, fldB2 _fvecalignment;
+    kf::FieldRegion<fvec> fld _fvecalignment;
+
+
+    kf::FieldValue<fvec> fldB01, fldB11, fldB21 _fvecalignment;
+    kf::FieldRegion<fvec> fld1 _fvecalignment;
+
+    const int nStations = fParameters.GetNstationsActive();
+    int nTracks_SIMD    = fvec::size();
+
+    kf::TrackKalmanFilter<fvec> fit;  // fit parameters coresponding to the current track
+    TrackParamV& tr = fit.Tr();
+    fit.SetParticleMass(fDefaultMass);
+    fit.SetDoFitVelocity(true);
+
+    Track* t[fvec::size()]{nullptr};
+
+    const ca::Station<fvec>* sta = fParameters.GetStations().begin();
+
+    // Spatial-time position of a hit vs. station and track in the portion
+
+    fvec x[constants::size::MaxNstations];                       // Hit position along the x-axis [cm]
+    fvec y[constants::size::MaxNstations];                       // Hit position along the y-axis [cm]
+    kf::MeasurementXy<fvec> mxy[constants::size::MaxNstations];  // Covariance matrix for x,y
+
+    fvec z[constants::size::MaxNstations];  // Hit position along the z-axis (precised) [cm]
+
+    fvec time[constants::size::MaxNstations];  // Hit time [ns]
+    fvec dt2[constants::size::MaxNstations];   // Hit time uncertainty [ns] squared
+
+    fvec x_first;
+    fvec y_first;
+    kf::MeasurementXy<fvec> mxy_first;
+
+    fvec time_first;
+    fvec wtime_first;
+    fvec dt2_first;
+
+    fvec x_last;
+    fvec y_last;
+    kf::MeasurementXy<fvec> mxy_last;
+
+    fvec time_last;
+    fvec wtime_last;
+    fvec dt2_last;
+
+    fvec By[constants::size::MaxNstations];
+    fmask w[constants::size::MaxNstations];
+    fmask w_time[constants::size::MaxNstations];  // !!!
+
+    fvec y_temp;
+    fvec x_temp;
+    fvec fldZ0;
+    fvec fldZ1;
+    fvec fldZ2;
+    fvec z_start;
+    fvec z_end;
+
+    kf::FieldValue<fvec> fB[constants::size::MaxNstations], fB_temp _fvecalignment;
+
+
+    fvec ZSta[constants::size::MaxNstations];
+    for (int ista = 0; ista < nStations; ista++) {
+      ZSta[ista] = sta[ista].fZ;
+      mxy[ista].SetCov(1., 0., 1.);
+    }
+
+    unsigned short N_vTracks = wData.RecoTracks().size();  // number of tracks processed per one SSE register
+
+    for (unsigned short itrack = 0; itrack < N_vTracks; itrack += fvec::size()) {
+
+      if (N_vTracks - itrack < static_cast<unsigned short>(fvec::size())) nTracks_SIMD = N_vTracks - itrack;
+
+      for (int i = 0; i < nTracks_SIMD; i++) {
+        t[i] = &wData.RecoTrack(itrack + i);
+      }
+
+      // fill the rest of the SIMD vectors with something reasonable
+      for (uint i = nTracks_SIMD; i < fvec::size(); i++) {
+        t[i] = &wData.RecoTrack(itrack);
+      }
+
+      // get hits of current track
+      for (int ista = 0; ista < nStations; ista++) {
+        w[ista]      = fmask::Zero();
+        w_time[ista] = fmask::Zero();
+        z[ista]      = ZSta[ista];
+      }
+
+      //fmask isFieldPresent = fmask::Zero();
+
+      for (int iVec = 0; iVec < nTracks_SIMD; iVec++) {
+
+        int nHitsTrack = t[iVec]->fNofHits;
+        int iSta[constants::size::MaxNstations];
+
+        for (int ih = 0; ih < nHitsTrack; ih++) {
+
+          const ca::Hit& hit = input.GetHit(wData.RecoHitIndex(start_hit++));
+          const int ista     = hit.Station();
+          auto [detSystemId, iStLocal] = fParameters.GetActiveSetup().GetIndexMap().GlobalToLocal<EDetectorID>(ista);
+
+          //if (sta[ista].fieldStatus) { isFieldPresent[iVec] = true; }
+
+          iSta[ih]      = ista;
+          w[ista][iVec] = true;
+          if (sta[ista].timeInfo) {
+            w_time[ista][iVec] = true;
+          }
+          // subtract misalignment tolerances to get the original hit errors
+          float dX2Orig = hit.dX2() - fParameters.GetMisalignmentXsq(detSystemId);
+          float dY2Orig = hit.dY2() - fParameters.GetMisalignmentYsq(detSystemId);
+          float dXYOrig = hit.dXY();
+          if (dX2Orig < 0. || dY2Orig < 0. || fabs(dXYOrig / sqrt(dX2Orig * dY2Orig)) > 1.) {
+            dX2Orig = hit.dX2();
+            dY2Orig = hit.dY2();
+          }
+          float dT2Orig = hit.dT2() - fParameters.GetMisalignmentTsq(detSystemId);
+          if (dT2Orig < 0.) {
+            dT2Orig = hit.dT2();
+          }
+
+          x[ista][iVec]    = hit.X();  //x_temp[iVec];
+          y[ista][iVec]    = hit.Y();  //y_temp[iVec];
+          time[ista][iVec] = hit.T();
+          dt2[ista][iVec]  = dT2Orig;
+          if (!sta[ista].timeInfo) {
+            dt2[ista][iVec] = 1.e4;
+          }
+          z[ista][iVec]          = hit.Z();
+          fB_temp                = sta[ista].fieldSlice.GetFieldValue(x[ista], y[ista]);
+          mxy[ista].X()[iVec]    = hit.X();
+          mxy[ista].Y()[iVec]    = hit.Y();
+          mxy[ista].Dx2()[iVec]  = dX2Orig;
+          mxy[ista].Dy2()[iVec]  = dY2Orig;
+          mxy[ista].Dxy()[iVec]  = dXYOrig;
+          mxy[ista].NdfX()[iVec] = 1.;
+          mxy[ista].NdfY()[iVec] = 1.;
+
+          fB[ista].SetSimdEntry(fB_temp.GetBx()[iVec], fB_temp.GetBy()[iVec], fB_temp.GetBz()[iVec], iVec);
+
+          if (ih == 0) {
+            z_start[iVec]          = z[ista][iVec];
+            x_first[iVec]          = x[ista][iVec];
+            y_first[iVec]          = y[ista][iVec];
+            time_first[iVec]       = time[ista][iVec];
+            wtime_first[iVec]      = sta[ista].timeInfo ? 1. : 0.;
+            dt2_first[iVec]        = dt2[ista][iVec];
+            mxy_first.X()[iVec]    = mxy[ista].X()[iVec];
+            mxy_first.Y()[iVec]    = mxy[ista].Y()[iVec];
+            mxy_first.Dx2()[iVec]  = mxy[ista].Dx2()[iVec];
+            mxy_first.Dy2()[iVec]  = mxy[ista].Dy2()[iVec];
+            mxy_first.Dxy()[iVec]  = mxy[ista].Dxy()[iVec];
+            mxy_first.NdfX()[iVec] = mxy[ista].NdfX()[iVec];
+            mxy_first.NdfY()[iVec] = mxy[ista].NdfY()[iVec];
+          }
+          else if (ih == nHitsTrack - 1) {
+            z_end[iVec]           = z[ista][iVec];
+            x_last[iVec]          = x[ista][iVec];
+            y_last[iVec]          = y[ista][iVec];
+            mxy_last.X()[iVec]    = mxy[ista].X()[iVec];
+            mxy_last.Y()[iVec]    = mxy[ista].Y()[iVec];
+            mxy_last.Dx2()[iVec]  = mxy[ista].Dx2()[iVec];
+            mxy_last.Dy2()[iVec]  = mxy[ista].Dy2()[iVec];
+            mxy_last.Dxy()[iVec]  = mxy[ista].Dxy()[iVec];
+            mxy_last.NdfX()[iVec] = mxy[ista].NdfX()[iVec];
+            mxy_last.NdfY()[iVec] = mxy[ista].NdfY()[iVec];
+            time_last[iVec]       = time[ista][iVec];
+            dt2_last[iVec]        = dt2[ista][iVec];
+            wtime_last[iVec]      = sta[ista].timeInfo ? 1. : 0.;
+          }
+        }
+
+        for (int ih = nHitsTrack - 1; ih >= 0; ih--) {
+          const int ista              = iSta[ih];
+          const ca::Station<fvec>& st = fParameters.GetStation(ista);
+          By[ista][iVec]              = st.fieldSlice.GetFieldValue(0., 0.).GetBy()[0];
+        }
+      }
+
+      fit.GuessTrack(z_end, x, y, z, time, By, w, w_time, nStations);
+
+      if (ca::TrackingMode::kGlobal == fTrackingMode || ca::TrackingMode::kMcbm == fTrackingMode) {
+        tr.Qp() = fvec(1. / 1.1);
+      }
+
+      // tr.Qp() = iif(isFieldPresent, tr.Qp(), fvec(1. / 0.25));
+
+      for (int iter = 0; iter < 2; iter++) {  // 1.5 iterations
+
+        fit.SetQp0(tr.Qp());
+
+        // fit backward
+
+        int ista = nStations - 1;
+
+        time_last = iif(w_time[ista], time_last, fvec::Zero());
+        dt2_last  = iif(w_time[ista], dt2_last, fvec(1.e6));
+
+        tr.ResetErrors(mxy_last.Dx2(), mxy_last.Dy2(), 0.1, 0.1, 1.0, dt2_last, 1.e-2);
+        tr.C10()  = mxy_last.Dxy();
+        tr.X()    = mxy_last.X();
+        tr.Y()    = mxy_last.Y();
+        tr.Time() = time_last;
+        tr.Vi()   = constants::phys::SpeedOfLightInv;
+        tr.InitVelocityRange(0.5);
+        tr.Ndf()     = fvec(-5.) + fvec(2.);
+        tr.NdfTime() = fvec(-2.) + wtime_last;
+
+        fldZ1 = z[ista];
+
+        fldB1 = sta[ista].fieldSlice.GetFieldValue(tr.X(), tr.Y());
+
+        fldB1.SetSimdEntries(fB[ista], w[ista]);
+
+        fldZ2   = z[ista - 2];
+        fvec dz = fldZ2 - fldZ1;
+        fldB2   = sta[ista].fieldSlice.GetFieldValue(tr.X() + tr.Tx() * dz, tr.Y() + tr.Ty() * dz);
+        fldB2.SetSimdEntries(fB[ista - 2], w[ista - 2]);
+        fld.Set(fldB2, fldZ2, fldB1, fldZ1, fldB0, fldZ0);
+
+        for (--ista; ista >= 0; ista--) {
+
+          fldZ0 = z[ista];
+          dz    = (fldZ1 - fldZ0);
+          fldB0 = sta[ista].fieldSlice.GetFieldValue(tr.X() - tr.Tx() * dz, tr.Y() - tr.Ty() * dz);
+          fldB0.SetSimdEntries(fB[ista], w[ista]);
+          fld.Set(fldB0, fldZ0, fldB1, fldZ1, fldB2, fldZ2);
+
+          fmask initialised = (z[ista] < z_end) & (z_start <= z[ista]);
+
+          fld1 = fld;
+
+          fit.SetMask(initialised);
+          fit.Extrapolate(z[ista], fld1);
+          auto radThick = fSetup.GetMaterial(ista).GetThicknessX0(tr.X(), tr.Y());
+          fit.MultipleScattering(radThick);
+          fit.EnergyLossCorrection(radThick, kf::FitDirection::kUpstream);
+
+          fit.SetMask(initialised && w[ista]);
+          fit.FilterXY(mxy[ista]);
+          fit.FilterTime(time[ista], dt2[ista], fmask(sta[ista].timeInfo));
+
+
+          fldB2 = fldB1;
+          fldZ2 = fldZ1;
+          fldB1 = fldB0;
+          fldZ1 = fldZ0;
+        }
+
+        // extrapolate to the PV region
+
+        kf::TrackKalmanFilter fitpv = fit;
+        {
+          fitpv.SetMask(fmask::One());
+
+          kf::MeasurementXy<fvec> vtxInfo = wData.TargetMeasurement();
+          vtxInfo.SetDx2(1.e-8);
+          vtxInfo.SetDxy(0.);
+          vtxInfo.SetDy2(1.e-8);
+
+          if (ca::TrackingMode::kGlobal == fTrackingMode) {
+            kf::FieldRegion<fvec> fldFull(kf::GlobalField::fgOriginalFieldType, kf::GlobalField::fgOriginalField);
+            fitpv.SetMaxExtrapolationStep(1.);
+            for (int vtxIter = 0; vtxIter < 2; vtxIter++) {
+              fitpv.SetQp0(fitpv.Tr().Qp());
+              fitpv.Tr()      = fit.Tr();
+              fitpv.Tr().Qp() = fitpv.Qp0();
+              fitpv.Extrapolate(fParameters.GetTargetPositionZ(), fldFull);
+              fitpv.FilterXY(vtxInfo);
+            }
+          }
+          else {
+            fitpv.SetQp0(fitpv.Tr().Qp());
+            fitpv.Extrapolate(fParameters.GetTargetPositionZ(), fld);
+          }
+        }
+
+        //fit.SetQp0(tr.Qp());
+        //fit.SetMask(fmask::One());
+        //fit.MeasureVelocityWithQp();
+        //fit.FilterVi(TrackParamV::kClightNsInv);
+
+        TrackParamV Tf = fit.Tr();
+        if (ca::TrackingMode::kGlobal == fTrackingMode) {
+          Tf.Qp() = fitpv.Tr().Qp();
+        }
+
+        for (int iVec = 0; iVec < nTracks_SIMD; iVec++) {
+          t[iVec]->fParFirst.Set(Tf, iVec);
+        }
+
+        for (int iVec = 0; iVec < nTracks_SIMD; iVec++) {
+          t[iVec]->fParPV.Set(fitpv.Tr(), iVec);
+        }
+
+        if (iter == 1) {
+          break;
+        }  // only 1.5 iterations
+
+        // fit forward
+
+        ista = 0;
+
+        tr.ResetErrors(mxy_first.Dx2(), mxy_first.Dy2(), 0.1, 0.1, 1., dt2_first, 1.e-2);
+        tr.C10() = mxy_first.Dxy();
+
+        tr.X()    = mxy_first.X();
+        tr.Y()    = mxy_first.Y();
+        tr.Time() = time_first;
+        tr.Vi()   = constants::phys::SpeedOfLightInv;
+        tr.InitVelocityRange(0.5);
+
+        tr.Ndf()     = fvec(-5. + 2.);
+        tr.NdfTime() = fvec(-2.) + wtime_first;
+
+        fit.SetQp0(tr.Qp());
+
+        fldZ1 = z[ista];
+        fldB1 = sta[ista].fieldSlice.GetFieldValue(tr.X(), tr.Y());
+        fldB1.SetSimdEntries(fB[ista], w[ista]);
+
+        fldZ2 = z[ista + 2];
+        dz    = fldZ2 - fldZ1;
+        fldB2 = sta[ista].fieldSlice.GetFieldValue(tr.X() + tr.Tx() * dz, tr.Y() + tr.Ty() * dz);
+        fldB2.SetSimdEntries(fB[ista + 2], w[ista + 2]);
+        fld.Set(fldB2, fldZ2, fldB1, fldZ1, fldB0, fldZ0);
+
+        for (++ista; ista < nStations; ista++) {
+          fldZ0 = z[ista];
+          dz    = (fldZ1 - fldZ0);
+          fldB0 = sta[ista].fieldSlice.GetFieldValue(tr.X() - tr.Tx() * dz, tr.Y() - tr.Ty() * dz);
+          fldB0.SetSimdEntries(fB[ista], w[ista]);
+          fld.Set(fldB0, fldZ0, fldB1, fldZ1, fldB2, fldZ2);
+
+          fmask initialised = (z[ista] <= z_end) & (z_start < z[ista]);
+
+          fit.SetMask(initialised);
+          fit.Extrapolate(z[ista], fld);
+          auto radThick = fSetup.GetMaterial(ista).GetThicknessX0(tr.X(), tr.Y());
+          fit.MultipleScattering(radThick);
+          fit.EnergyLossCorrection(radThick, kf::FitDirection::kDownstream);
+          fit.SetMask(initialised && w[ista]);
+          fit.FilterXY(mxy[ista]);
+          fit.FilterTime(time[ista], dt2[ista], fmask(sta[ista].timeInfo));
+
+          fldB2 = fldB1;
+          fldZ2 = fldZ1;
+          fldB1 = fldB0;
+          fldZ1 = fldZ0;
+        }
+
+        //fit.SetQp0(tr.Qp());
+        //fit.SetMask(fmask::One());
+        //fit.MeasureVelocityWithQp();
+        //fit.FilterVi(TrackParamV::kClightNsInv);
+
+        TrackParamV Tl = fit.Tr();
+        if (ca::TrackingMode::kGlobal == fTrackingMode) {
+          Tl.Qp() = fitpv.Tr().Qp();
+        }
+
+        for (int iVec = 0; iVec < nTracks_SIMD; iVec++) {
+          t[iVec]->fParLast.Set(Tl, iVec);
+        }
+
+      }  // iter
+    }
+  }
+}  // namespace cbm::algo::ca

algo/ca/core/tracking/CaTrackFitter.h

+/* Copyright (C) 2023-2024 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+   SPDX-License-Identifier: GPL-3.0-only
+   Authors: Sergei Zharko [committer], Maksym Zyzak */
+
+/// \file    TrackFitter.h
+/// \brief   A class wrapper over clones merger algorithm for the Ca track finder (declaration)
+/// \since   19.10.2023
+
+#pragma once  // include this header only once per compilation unit
+
+#include "CaBranch.h"
+#include "CaHit.h"
+#include "CaInputData.h"
+#include "CaParameters.h"
+#include "CaSimd.h"
+#include "CaVector.h"
+#include "CaWindowData.h"
+#include "KfTrackParam.h"
+
+
+namespace cbm::algo::ca
+{
+  class Track;
+
+  /// Class implements a track fit the CA track finder
+  ///
+  class TrackFitter {
+   public:
+    /// Default constructor
+    TrackFitter(const ca::Parameters<fvec>& pars, const fscal mass, const ca::TrackingMode& mode);
+
+    /// Destructor
+    ~TrackFitter();
+
+    /// Copy constructor
+    TrackFitter(const TrackFitter&) = default;
+
+    /// Move constructor
+    TrackFitter(TrackFitter&&) = default;
+
+    /// Copy assignment operator
+    TrackFitter& operator=(const TrackFitter&) = delete;
+
+    /// Move assignment operator
+    TrackFitter& operator=(TrackFitter&&) = delete;
+
+    /// Fit tracks, found by the CA tracker
+    void FitCaTracks(const ca::InputData& input, WindowData& wData);
+
+
+   private:
+    ///-------------------------------
+    /// Data members
+
+    const Parameters<fvec>& fParameters;            ///< Object of Framework parameters class
+    const cbm::algo::kf::Setup<fvec>& fSetup;       ///< Setup instance
+    fscal fDefaultMass{constants::phys::MuonMass};  ///< mass of the propagated particle [GeV/c2]
+    ca::TrackingMode fTrackingMode;
+  };
+
+}  // namespace cbm::algo::ca

algo/ca/core/tracking/CaTripletConstructor.cxx

+/* Copyright (C) 2010-2021 Frankfurt Institute for Advanced Studies, Goethe-Universitaet Frankfurt, Frankfurt
+   SPDX-License-Identifier: GPL-3.0-only
+   Authors: Igor Kulakov [committer], Maksym Zyzak, Valentina Akishina */
+
+#include "CaTripletConstructor.h"
+
+#include "CaDefines.h"
+#include "CaFramework.h"
+#include "CaGridArea.h"
+
+#include <algorithm>
+#include <iostream>
+// #include "CaToolsDebugger.h"
+#include "AlgoFairloggerCompat.h"
+#include "KfTrackKalmanFilter.h"
+#include "KfTrackParam.h"
+
+// using cbm::ca::tools::Debugger;
+
+namespace cbm::algo::ca
+{
+  TripletConstructor::TripletConstructor(const ca::Parameters<fvec>& pars, WindowData& wData, const fscal mass,
+                                         const ca::TrackingMode& mode)
+    : fParameters(pars)
+    , fSetup(fParameters.GetActiveSetup())
+    , frWData(wData)
+    , fDefaultMass(mass)
+    , fTrackingMode(mode)
+  {
+    // FIXME: SZh 24.08.2022
+    //        This approach is suitable only for a case, when all the stations inside a magnetic field come right before
+    //        all the stations outside of the field!
+    fNfieldStations = std::lower_bound(fParameters.GetStations().cbegin(),
+                                       fParameters.GetStations().cbegin() + fParameters.GetNstationsActive(),
+                                       0,  // we are looking for the first zero element
+                                       [](const ca::Station<fvec>& s, int edge) { return bool(s.fieldStatus) > edge; })
+                      - fParameters.GetStations().cbegin();
+
+    fIsTargetField = !frWData.TargB().IsZero();
+  }
+
+
+  bool TripletConstructor::InitStations(int istal, int istam, int istar)
+  {
+    fIstaL = istal;
+    fIstaM = istam;
+    fIstaR = istar;
+
+    if (fIstaM >= fParameters.GetNstationsActive()) {
+      return false;
+    }
+    fStaL = &fParameters.GetStation(fIstaL);
+    fStaM = &fParameters.GetStation(fIstaM);
+    fStaR = &fParameters.GetStation(fIstaR);
+
+    {  // two stations for approximating the field between the target and the left hit
+      const int sta1 = (fNfieldStations <= 1) ? 1 : std::clamp(fIstaL, 1, fNfieldStations - 1);
+      const int sta0 = sta1 / 2;  // station between fIstaL and the target
+
+      assert(0 <= sta0 && sta0 < sta1 && sta1 < fParameters.GetNstationsActive());
+      fFld0Sta[0] = &fParameters.GetStation(sta0);
+      fFld0Sta[1] = &fParameters.GetStation(sta1);
+    }
+
+    {  // three stations for approximating the field between the left and the right hit
+
+      int sta0 = fIstaL;
+      int sta1 = fIstaM;
+      int sta2 = fIstaM + 1;
+      if (sta2 >= fParameters.GetNstationsActive()) {
+        sta2 = fIstaM;
+        sta1 = fIstaM - 1;
+        sta0 = (sta1 <= 0) ? 0 : std::clamp(sta0, 0, sta1 - 1);
+      }
+      if (fParameters.GetNstationsActive() >= 3) {
+        assert(0 <= sta0 && sta0 < sta1 && sta1 < sta2 && sta2 < fParameters.GetNstationsActive());
+      }
+
+      fFld1Sta[0] = &fParameters.GetStation(sta0);
+      fFld1Sta[1] = &fParameters.GetStation(sta1);
+      fFld1Sta[2] = &fParameters.GetStation(sta2);
+    }
+    return true;
+  }
+
+
+  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;
+    }
+
+    kf::TrackKalmanFilter<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.
+      //kf::FieldValue<fvec> lB, mB, cB, rB _fvecalignment;           currently not used
+      //kf::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)
+      {
+        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
+
+      // 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
+      kf::FieldRegion<fvec> fld0;
+      {
+        kf::FieldValue<fvec> B0 = fFld0Sta[0]->fieldSlice.GetFieldValueForLine(T);
+        kf::FieldValue<fvec> B1 = fFld0Sta[1]->fieldSlice.GetFieldValueForLine(T);
+        fld0.Set(frWData.TargB(), fParameters.GetTargetPositionZ(), B0, fFld0Sta[0]->fZ, B1, fFld0Sta[1]->fZ);
+      }
+
+      {  // field, made by the left hit, the middle station and the right station
+        // Will be used for extrapolation to the right hit
+        kf::FieldValue<fvec> B0 = fFld1Sta[0]->fieldSlice.GetFieldValueForLine(T);
+        kf::FieldValue<fvec> B1 = fFld1Sta[1]->fieldSlice.GetFieldValueForLine(T);
+        kf::FieldValue<fvec> B2 = fFld1Sta[2]->fieldSlice.GetFieldValueForLine(T);
+        fFldL.Set(B0, fFld1Sta[0]->fZ, B1, fFld1Sta[1]->fZ, B2, fFld1Sta[2]->fZ);
+      }
+
+      // add the target constraint
+      fit.FilterWithTargetAtLine(fParameters.GetTargetPositionZ(), frWData.TargetMeasurement(), fld0);
+      fit.MultipleScattering(fSetup.GetMaterial(fIstaL).GetThicknessX0(T.GetX(), T.GetY()));
+
+      // extrapolate to the middle hit
+      fit.ExtrapolateLine(fStaM->fZ, fFldL);
+    }
+
+    /// Find the doublets. Reformat data into portions of doublets.
+    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;
+    }
+
+    FindTripletHits();
+    FindTriplets();
+    SelectTriplets(tripletsOut);
+  }
+
+
+  void TripletConstructor::FindDoublets(kf::TrackKalmanFilter<fvec>& fit)
+  {
+    // ---- Add the middle hits to parameters estimation ----
+    Vector<TrackParamV>& tracks   = fDoubletData.first;
+    Vector<ca::HitIndex_t>& hitsM = fDoubletData.second;
+
+    tracks.clear();
+    tracks.reserve(hitsM.size());
+
+    const bool isMomentumFitted = (fIsTargetField || (fStaL->fieldStatus != 0) || (fStaM->fieldStatus != 0));
+
+    const TrackParamV Tr = fit.Tr();  // copy contents of fit
+
+    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 = fit.Tr();
+      T2              = Tr;
+      fit.SetQp0(fvec(0.f));
+
+      fvec z_2 = hitm.Z();
+      kf::MeasurementXy<fvec> m_2(hitm.X(), hitm.Y(), hitm.dX2(), hitm.dY2(), hitm.dXY(), fvec::One(), fvec::One());
+      fvec t_2   = hitm.T();
+      fvec dt2_2 = hitm.dT2();
+
+      // add the middle hit
+      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(fSetup.GetMaterial(fIstaM).GetThicknessX0(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) {
+        // TODO: SG: adjust cuts, put them to parameters
+
+        const fscal tx = T2.Tx()[0];
+        const fscal ty = T2.Ty()[0];
+        const fscal tt = T2.Vi()[0] * sqrt(1. + tx * tx + ty * ty);  // dt/dl * dl/dz
+
+        for (auto itClone = it + 1; itClone != hitsM.end(); itClone++) {
+
+          const int indClone      = *itClone;
+          const ca::Hit& hitClone = frWData.Hit(indClone);
+
+          const fscal dz = hitClone.Z() - T2.Z()[0];
+
+          if ((fStaM->timeInfo) && (T2.NdfTime()[0] >= 0)) {
+            const fscal dt = T2.Time()[0] + tt * dz - hitClone.T();
+            if (!(fabs(dt) <= 3.5 * sqrt(T2.C55()[0]) + hitClone.RangeT())) {
+              continue;
+            }
+          }
+
+          const fscal dx = T2.GetX()[0] + tx * dz - hitClone.X();
+          if (!(fabs(dx) <= 3.5 * sqrt(T2.C00()[0]) + hitClone.RangeX())) {
+            continue;
+          }
+
+          const fscal dy = T2.Y()[0] + ty * dz - hitClone.Y();
+          if (!(fabs(dy) <= 3.5 * sqrt(T2.C11()[0]) + hitClone.RangeY())) {
+            continue;
+          }
+
+          if (fParameters.DevIsSuppressOverlapHitsViaMc()) {
+            const int iMC = ca::Framework::GetMcTrackIdForWindowHit(fIhitL);
+            if ((iMC != ca::Framework::GetMcTrackIdForWindowHit(indM))
+                || (iMC != ca::Framework::GetMcTrackIdForWindowHit(indClone))) {
+              continue;
+            }
+          }
+
+          frWData.SuppressHit(indClone);
+        }
+      }
+
+      tracks.push_back(T2);
+      *it2 = indM;
+      it2++;
+    }  // it
+    hitsM.shrink(std::distance(hitsM.begin(), it2));
+  }
+
+
+  void TripletConstructor::FindTripletHits()
+  {
+    //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;
+
+    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.
+    kf::TrackKalmanFilter<fvec> fit(fmask::One(), true);
+    fit.SetParticleMass(frWData.CurrentIteration()->GetElectronFlag() ? constants::phys::ElectronMass : fDefaultMass);
+    fit.SetQp0(fvec(0.));
+
+    {
+      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 < hitsM_2.size(); i2++) {
+
+      fit.SetTrack(tracks_2[i2]);
+      TrackParamV& T2 = fit.Tr();
+
+      // extrapolate to the right hit station
+      fit.Extrapolate(fStaR->fZ, fFldL);
+
+      if constexpr (fDebugDublets) {
+        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 = frWData.Hit(iwhit[0]);
+        const ca::Hit& h1 = frWData.Hit(iwhit[1]);
+
+        LOG(info) << "\n======  Extrapolated Doublet : "
+                  << "  iter " << frWData.CurrentIteration()->GetName() << " hits: {" << fIstaL << "/" << ih0 << " "
+                  << fIstaM << "/" << ih1 << " " << fIstaR << "/?} xyz: {" << h0.X() << " " << h0.Y() << " " << h0.Z()
+                  << "}, {" << h1.X() << " " << h1.Y() << " " << h1.Z() << "} chi2 " << T2.GetChiSq()[0] << " ndf "
+                  << T2.Ndf()[0] << " chi2time " << T2.ChiSqTime()[0] << " ndfTime " << T2.NdfTime()[0];
+        LOG(info) << "  extr. track: " << T2.ToString(0);
+      }
+
+      // ---- Find the triplets(right hit). Reformat data in the portion of triplets. ----
+      int iMC = -1;
+
+      auto rejectDoublet = [&]() -> bool {
+        if (ca::TrackingMode::kSts == fTrackingMode && (T2.C44()[0] < 0)) {
+          return true;
+        }
+        if (T2.C00()[0] < 0 || T2.C11()[0] < 0 || T2.C22()[0] < 0 || T2.C33()[0] < 0 || T2.C55()[0] < 0) {
+          return true;
+        }
+        if (fabs(T2.Tx()[0]) > maxSlope) {
+          return true;
+        }
+        if (fabs(T2.Ty()[0]) > maxSlope) {
+          return true;
+        }
+        if (fParameters.DevIsMatchTripletsViaMc()) {
+          int indM = hitsM_2[i2];
+          iMC      = ca::Framework::GetMcTrackIdForWindowHit(fIhitL);
+          if (iMC < 0 || iMC != ca::Framework::GetMcTrackIdForWindowHit(indM)) {
+            return true;
+          }
+        }
+        return false;
+      };
+      const bool isDoubletGood = !rejectDoublet();
+
+      if constexpr (fDebugDublets) {
+        if (isDoubletGood) {
+          LOG(info) << "  extrapolated doublet accepted";
+        }
+        else {
+          LOG(info) << "  extrapolated doublet rejected";
+          LOG(info) << "======== end of extrapolated doublet ==== \n";
+        }
+      }
+
+      if (!isDoubletGood) {
+        continue;
+      }
+
+      Vector<ca::HitIndex_t> collectedHits;
+      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
+        // of our tests and finally results in crashes on run4
+        //      LOG(debug) << "Ca: GetMaxTripletPerDoublets==" << fParameters.GetMaxTripletPerDoublets()
+        //                 << " reached in findTripletsStep0()";
+        // reject already created triplets for this doublet
+        collectedHits.clear();
+      }
+      if constexpr (fDebugDublets) {
+        LOG(info) << " collected " << collectedHits.size() << " hits on the right station ";
+      }
+      for (ca::HitIndex_t& irh : collectedHits) {
+        if constexpr (fDebugDublets) {
+          const ca::Hit& h    = frWData.Hit(irh);
+          ca::HitIndex_t ihit = h.Id();
+          LOG(info) << "  hit " << ihit << " " << h.ToString();
+        }
+        if (frWData.IsHitSuppressed(irh)) {
+          if constexpr (fDebugDublets) {
+            LOG(info) << "  the hit is suppressed";
+          }
+          continue;
+        }
+        // pack the triplet
+        hitsM_3.push_back(hitsM_2[i2]);
+        hitsR_3.push_back(irh);
+      }  // search area
+      if constexpr (fDebugDublets) {
+        LOG(info) << "======== end of extrapolated doublet ==== \n";
+      }
+    }  // i2
+  }
+
+  void TripletConstructor::FindTriplets()
+  {
+    constexpr int nIterations = 2;
+
+    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());
+
+    tracks.clear();
+    tracks.reserve(hitsM.size());
+
+    /// Refit Triplets
+    if constexpr (fDebugTriplets) {
+      //cbm::ca::tools::Debugger::Instance().AddNtuple(
+      //  "triplets", "ev:iter:i0:x0:y0:z0:i1:x1:y1:z1:i2:x2:y2:z2:mc:sta:p:vx:vy:vz:chi2:ndf:chi2time:ndfTime");
+    }
+
+    kf::TrackKalmanFilter<fvec> fit;
+    fit.SetMask(fmask::One());
+    fit.SetParticleMass(frWData.CurrentIteration()->GetElectronFlag() ? constants::phys::ElectronMass : fDefaultMass);
+
+    // prepare data
+    const int NHits       = 3;  // triplets
+    const int ista[NHits] = {fIstaL, fIstaM, fIstaR};
+
+    const ca::Station<fvec> sta[NHits] = {fParameters.GetStation(ista[0]), fParameters.GetStation(ista[1]),
+                                          fParameters.GetStation(ista[2])};
+
+    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
+
+    for (size_t i3 = 0; i3 < hitsM.size(); ++i3) {
+
+      // prepare data
+      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()};
+
+      if (fParameters.DevIsMatchTripletsViaMc()) {
+        int mc1 = ca::Framework::GetMcTrackIdForCaHit(ihit[0]);
+        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;
+        }
+      }
+
+      fscal x[NHits], y[NHits], z[NHits], t[NHits];
+      fscal dt2[NHits];
+      kf::MeasurementXy<fvec> mxy[NHits];
+
+      for (int ih = 0; ih < NHits; ++ih) {
+        const ca::Hit& hit = frWData.Hit(iwhit[ih]);
+        mxy[ih] = kf::MeasurementXy<fvec>(hit.X(), hit.Y(), hit.dX2(), hit.dY2(), hit.dXY(), fvec::One(), fvec::One());
+
+        x[ih]   = hit.X();
+        y[ih]   = hit.Y();
+        z[ih]   = hit.Z();
+        t[ih]   = hit.T();
+        dt2[ih] = hit.dT2();
+      };
+
+      // find the field along the track
+
+      kf::FieldValue<fvec> B[3] _fvecalignment;
+      kf::FieldRegion<fvec> fld _fvecalignment;
+      kf::FieldRegion<fvec> fldTarget _fvecalignment;
+
+      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);
+      };
+
+      fld.Set(B[0], sta[0].fZ, B[1], sta[1].fZ, B[2], sta[2].fZ);
+      fldTarget.Set(frWData.TargB(), fParameters.GetTargetPositionZ(), B[0], sta[0].fZ, B[1], sta[1].fZ);
+
+      TrackParamV& T = fit.Tr();
+
+      // initial parameters
+      {
+        fvec dz01 = 1. / (z[1] - z[0]);
+        T.Tx()    = (x[1] - x[0]) * dz01;
+        T.Ty()    = (y[1] - y[0]) * dz01;
+        T.Qp()    = 0.;
+        T.Vi()    = 0.;
+      }
+
+      // repeat several times in order to improve the precision
+      for (int iiter = 0; iiter < nIterations; ++iiter) {
+
+        auto fitTrack = [&](int startIdx, int endIdx, int step, kf::FitDirection direction) {
+          const fvec maxQp = frWData.CurrentIteration()->GetMaxQp();
+          fit.SetQp0(T.Qp());
+          fit.Qp0()(fit.Qp0() > maxQp)  = maxQp;
+          fit.Qp0()(fit.Qp0() < -maxQp) = -maxQp;
+
+          int ih0  = startIdx;
+          T.X()    = x[ih0];
+          T.Y()    = y[ih0];
+          T.Z()    = z[ih0];
+          T.Time() = t[ih0];
+          T.Qp()   = 0.;
+          T.Vi()   = 0.;
+
+          T.ResetErrors(1., 1., 1., 1., 100., (sta[ih0].timeInfo ? dt2[ih0] : 1.e6), 1.e2);
+          T.C00() = mxy[ih0].Dx2();
+          T.C10() = mxy[ih0].Dxy();
+          T.C11() = mxy[ih0].Dy2();
+
+          T.Ndf()     = -ndfTrackModel + 2;
+          T.NdfTime() = sta[ih0].timeInfo ? 0 : -1;
+
+          if (startIdx == 0) {  //only for the forward fit
+            fit.FilterWithTargetAtLine(fParameters.GetTargetPositionZ(), frWData.TargetMeasurement(), fldTarget);
+          }
+
+          for (int ih = startIdx + step; ih != endIdx; ih += step) {
+            fit.Extrapolate(z[ih], fld);
+            auto radThick = fSetup.GetMaterial(ista[ih]).GetThicknessX0(T.X(), T.Y());
+            fit.MultipleScattering(radThick);
+            fit.EnergyLossCorrection(radThick, direction);
+            fit.FilterXY(mxy[ih]);
+            fit.FilterTime(t[ih], dt2[ih], fmask(sta[ih].timeInfo));
+          }
+        };
+
+        // Fit downstream
+        fitTrack(0, NHits, 1, kf::FitDirection::kDownstream);
+
+        if (iiter == nIterations - 1) break;
+
+        // Fit upstream
+        fitTrack(NHits - 1, -1, -1, kf::FitDirection::kUpstream);
+      }  // for iiter
+
+      tracks.push_back(T);
+
+      if constexpr (fDebugTriplets) {
+        int ih0 = ihit[0];
+        int ih1 = ihit[1];
+        int ih2 = ihit[2];
+        int mc1 = ca::Framework::GetMcTrackIdForCaHit(ih0);
+        int mc2 = ca::Framework::GetMcTrackIdForCaHit(ih1);
+        int mc3 = ca::Framework::GetMcTrackIdForCaHit(ih2);
+
+        if (1 || (mc1 >= 0) && (mc1 == mc2) && (mc1 == mc3)) {
+          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 << " "
+                    << fIstaM << "/" << ih1 << " " << fIstaR << "/" << ih2 << "} xyz: {" << h0.X() << " " << h0.Y()
+                    << " " << h0.Z() << "}, {" << h1.X() << " " << h1.Y() << " " << h1.Z() << "}, {" << h2.X() << ", "
+                    << h2.Y() << ", " << h2.Z() << "} chi2 " << T.GetChiSq()[0] << " ndf " << T.Ndf()[0] << " chi2time "
+                    << T.ChiSqTime()[0] << " ndfTime " << T.NdfTime()[0];
+          /*
+          cbm::ca::tools::Debugger::Instance().FillNtuple(
+            "triplets", mctr.iEvent, frAlgo.fCurrentIterationIndex, ih0, h0.X(), h0.Y(), h0.Z(), ih1, h1.X(), h1.Y(),
+            h1.Z(), ih2, h2.X(), h2.Y(), h2.Z(), mc1, fIstaL, mctr.p, mctr.x, mctr.y, mctr.z, (fscal) T.GetChiSq()[0],
+            (fscal) T.Ndf()[0], (fscal) T.ChiSqTime()[0], (fscal) T.NdfTime()[0]);
+          */
+        }
+      }
+    }  //i3
+  }    // FindTriplets
+
+
+  void TripletConstructor::SelectTriplets(Vector<ca::Triplet>& tripletsOut)
+  {
+    /// Selects good triplets and saves them into tripletsOut.
+    /// Finds neighbouring triplets at the next station.
+
+    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));
+
+    tripletsOut.clear();
+    tripletsOut.reserve(hitsM.size());
+
+    for (size_t i3 = 0; i3 < hitsM.size(); ++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 = 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 > frWData.CurrentIteration()->GetTripletFinalChi2Cut()) {
+          continue;
+        }
+      }
+      // assert(std::isfinite(chi2) && chi2 > 0);
+
+      if (!std::isfinite(chi2) || chi2 < 0) {
+        continue;
+      }
+      //if (!T3.IsEntryConsistent(true, 0)) { continue; }
+
+      fscal qp  = T3.Qp()[0];
+      fscal Cqp = T3.C44()[0];
+
+      // TODO: SG: a magic correction that comes from the legacy code
+      // removing it leads to a higher ghost ratio
+      Cqp += 0.001;
+
+      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(Vector<ca::HitIndex_t>& collectedHits, kf::TrackKalmanFilter<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);
+
+    TrackParamV& T = fit.Tr();
+    //LOG(info) << T.chi2[0] ;
+    T.ChiSq() = 0.;
+
+    // if make it bigger the found hits will be rejected later because of the chi2 cut.
+    const fvec Pick_m22    = (fvec(chi2Cut) - T.GetChiSq());
+    const fscal timeError2 = T.C55()[0];
+    const fscal time       = T.Time()[0];
+
+    const auto& grid = frWData.Grid(iSta);
+    const fvec maxDZ = frWData.CurrentIteration()->GetMaxDZ();
+    ca::GridArea area(grid, T.X()[0], T.Y()[0],
+                      (sqrt(Pick_m22 * T.C00()) + grid.GetMaxRangeX() + maxDZ * kf::utils::fabs(T.Tx()))[0],
+                      (sqrt(Pick_m22 * T.C11()) + grid.GetMaxRangeY() + maxDZ * kf::utils::fabs(T.Ty()))[0]);
+    if constexpr (fDebugCollectHits) {
+      LOG(info) << "search area: " << T.X()[0] << " " << T.Y()[0] << " "
+                << (sqrt(Pick_m22 * T.C00()) + grid.GetMaxRangeX() + maxDZ * kf::utils::fabs(T.Tx()))[0] << " "
+                << (sqrt(Pick_m22 * T.C11()) + grid.GetMaxRangeY() + maxDZ * kf::utils::fabs(T.Ty()))[0];
+    }
+    if (fParameters.DevIsIgnoreHitSearchAreas()) {
+      area.DoLoopOverEntireGrid();
+    }
+
+    // loop over station hits (index incremented in condition)
+    for (ca::HitIndex_t ih = 0; area.GetNextObjectId(ih) && ((int) collectedHits.size() < maxNhits);) {
+
+      if (frWData.IsHitSuppressed(ih)) {
+        continue;
+      }
+
+      const ca::Hit& hit = frWData.Hit(ih);
+      if constexpr (fDebugCollectHits) {
+        LOG(info) << "hit in the search area: " << hit.ToString();
+        LOG(info) << "   check the hit.. ";
+      }
+
+      if (iMC >= 0) {  // match via MC
+        if (iMC != ca::Framework::GetMcTrackIdForWindowHit(ih)) {
+          if constexpr (fDebugCollectHits) {
+            LOG(info) << "   hit mc does not match";
+          }
+          continue;
+        }
+      }
+
+      // check time-boundaries
+
+      //TODO: move hardcoded cuts to parameters
+      if ((sta.timeInfo) && (T.NdfTime()[0] >= 0)) {
+        if (fabs(time - hit.T()) > 1.4 * (3.5 * sqrt(timeError2) + hit.RangeT())) {
+          if constexpr (fDebugCollectHits) {
+            LOG(info) << "   hit time does not match";
+          }
+          continue;
+        }
+        // if (fabs(time - hit.T()) > 30) continue;
+      }
+
+      // - check whether hit belong to the window ( track position +\- errors ) -
+      const fscal z = hit.Z();
+
+      // check y-boundaries
+      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];
+
+      if (fabs(dY) > dy_est) {
+        if constexpr (fDebugCollectHits) {
+          LOG(info) << "   hit Y does not match";
+        }
+        continue;
+      }
+
+      // check x-boundaries
+      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];
+
+      if (fabs(dX) > dx_est) {
+        if constexpr (fDebugCollectHits) {
+          LOG(info) << "   hit X does not match";
+        }
+        continue;
+      }
+
+      // check chi2
+      kf::MeasurementXy<fvec> mxy(hit.X(), hit.Y(), hit.dX2(), hit.dY2(), hit.dXY(), fvec::One(), fvec::One());
+
+      const fvec C10            = fit.ExtrapolateLineDxy(z);
+      const auto [chi2x, chi2u] = kf::TrackKalmanFilter<fvec>::GetChi2XChi2U(mxy, x, y, C00, C10, C11);
+
+      // TODO: adjust the cut, cut on chi2x & chi2u separately
+      if (!frWData.CurrentIteration()->GetTrackFromTripletsFlag()) {
+        if (chi2x[0] > chi2Cut) {
+          if constexpr (fDebugCollectHits) {
+            LOG(info) << "   hit chi2X is too large";
+          }
+          continue;
+        }
+        if (chi2u[0] > chi2Cut) {
+          if constexpr (fDebugCollectHits) {
+            LOG(info) << "   hit chi2U is too large";
+          }
+          continue;
+        }
+      }
+      if constexpr (fDebugCollectHits) {
+        LOG(info) << "   hit passed all the checks";
+      }
+      collectedHits.push_back(ih);
+
+    }  // loop over the hits in the area
+  }
+}  // namespace cbm::algo::ca

algo/ca/core/tracking/CaTripletConstructor.h

+/* Copyright (C) 2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+   SPDX-License-Identifier: GPL-3.0-only
+   Authors: Sergey Gorbunov [committer] */
+
+/// \file CaTripletConstructor.h
+/// \brief Triplet constructor for the CA tracker
+/// \author Sergey Gorbunov
+
+#pragma once  // include this header only once per compilation unit
+
+#include "CaFramework.h"
+#include "CaGridEntry.h"
+#include "CaStation.h"
+#include "CaTriplet.h"
+#include "CaVector.h"
+#include "CaWindowData.h"
+#include "KfFieldRegion.h"
+#include "KfSetup.h"
+#include "KfSimd.h"
+#include "KfTrackKalmanFilter.h"
+#include "KfTrackParam.h"
+
+namespace cbm::algo::ca
+{
+
+  /// Construction of triplets for the CA tracker
+  ///
+  class alignas(kf::VcMemAlign) TripletConstructor {
+
+   public:
+    ///                             ------  Constructors and destructor ------
+
+    /// Constructor
+    /// \param nThreads  Number of threads for multi-threaded mode
+    TripletConstructor(const ca::Parameters<fvec>& pars, WindowData& wData, const fscal mass,
+                       const ca::TrackingMode& mode);
+
+    /// Copy constructor
+    TripletConstructor(const TripletConstructor&) = delete;
+
+    /// Move constructor
+    TripletConstructor(TripletConstructor&&) = delete;
+
+    /// Copy assignment operator
+    TripletConstructor& operator=(const TripletConstructor&) = delete;
+
+    /// Move assignment operator
+    TripletConstructor& operator=(TripletConstructor&&) = delete;
+
+    /// Destructor
+    ~TripletConstructor() = default;
+
+    ///                             ------  FUNCTIONAL PART ------
+    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(kf::TrackKalmanFilter<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 FindTripletHits();
+
+    /// Find triplets on station
+    void FindTriplets();
+
+    /// Select good triplets.
+    void SelectTriplets(Vector<ca::Triplet>& tripletsOut);
+
+    void CollectHits(Vector<ca::HitIndex_t>& collectedHits, kf::TrackKalmanFilter<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 cbm::algo::kf::Setup<fvec>& fSetup;  ///< Reference to the setup
+    WindowData& frWData;
+    fscal fDefaultMass{constants::phys::MuonMass};  ///< mass of the propagated particle [GeV/c2]
+    ca::TrackingMode fTrackingMode;
+
+    bool fIsTargetField{false};  ///< is the magnetic field present at the target
+
+    // Number of stations situated in field region (MVD + STS in CBM)
+    int fNfieldStations{0};
+
+    int fIstaL{-1};  ///< left station index
+    int fIstaM{-1};  ///< middle station index
+    int fIstaR{-1};  ///< right station index
+
+    const ca::Station<fvec>* fStaL{nullptr};  ///< left station
+    const ca::Station<fvec>* fStaM{nullptr};  ///< mid station
+    const ca::Station<fvec>* fStaR{nullptr};  ///< right station
+
+    const ca::Station<fvec>*
+      fFld0Sta[2];  // two stations for approximating the field between the target and the left hit
+    const ca::Station<fvec>*
+      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
+    kf::FieldRegion<fvec> fFldL;
+
+    // Persistent storage for triplet tracks and hits
+    Triplet_t fTripletData;
+
+    // Persistent storage for doublet tracks and hits
+    Doublet_t fDoubletData;
+
+   private:
+    static constexpr bool fDebugDublets     = false;  // print debug info for dublets
+    static constexpr bool fDebugTriplets    = false;  // print debug info for triplets
+    static constexpr bool fDebugCollectHits = false;  // print debug info for CollectHits
+  };
+
+}  // namespace cbm::algo::ca

algo/ca/core/utils/CaDefines.h

+/* Copyright (C) 2010-2021 Frankfurt Institute for Advanced Studies, Goethe-Universitaet Frankfurt, Frankfurt
+   SPDX-License-Identifier: GPL-3.0-only
+   Authors: Maksym Zyzak, Igor Kulakov [committer], Sergey Gorbunov */
+
+/// \file CaDefines.h
+/// \brief Macros for the CA tracking algorithm
+/// \details Try to minimize the amount of macros. Try to call this header only from cxx files.
+
+#pragma once  // include this header only once per compilation unit
+
+#include "AlgoFairloggerCompat.h"
+
+#include <cassert>
+#include <sstream>
+
+// #define CBMCA_DEBUG_MODE
+
+#if defined(CBMCA_DEBUG_MODE)
+
+#define CBMCA_DEBUG_ASSERT(v)                                                                                          \
+  if (!(v)) {                                                                                                          \
+    LOG(error) << __FILE__ << ":" << __LINE__ << " assertion failed: " << #v << " = " << (v);                          \
+    assert(v);                                                                                                         \
+  }
+
+#define CBMCA_DEBUG_SHOW(expr)                                                                                         \
+  LOG(info) << __FILE__ << ":" << __LINE__ << ": \033[01;38;5;208m" << (#expr) << "\033[0m = " << (expr);
+
+#define CBMCA_DEBUG_SHOWF(msg)                                                                                         \
+  LOG(info) << "(!) " << __FILE__ << ":" << __LINE__ << ": \033[01;38;5;208m" << (#msg) << "\033[0m";
+
+#define CBMCA_DEBUG_SHOWCONTAINER(cont)                                                                                \
+  std::stringstream ss;                                                                                                \
+  ss << __FILE__ << ":" << __LINE__ << ": \033[01;38;5;208m" << (#cont) << "\033[0m: ";                                \
+  std::for_each(cont.cbegin(), cont.cend(), [](const auto& el) { ss << el << " "; });                                  \
+  LOG(info) << ss.str();
+
+#else  // not CBMCA_DEBUG_MODE
+
+#define CBMCA_DEBUG_ASSERT(v)
+
+#define CBMCA_DEBUG_SHOW(expr)
+
+#define CBMCA_DEBUG_SHOWF(msg)
+
+#define CBMCA_DEBUG_SHOWCONTAINER(cont)
+
+#endif  // CBMCA_DEBUG_MODE

algo/ca/core/utils/CaEnumArray.h

+/* Copyright (C) 2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+   SPDX-License-Identifier: GPL-3.0-only
+   Authors: Sergei Zharko [committer] */
+
+/// @file   CaEnumArray.h
+/// @brief  Implementation of cbm::algo::ca::EnumArray class
+/// @since  02.05.2023
+/// @author Sergei Zharko <s.zharko@gsi.de>
+
+#pragma once  // include this header only once per compilation unit
+
+#include <boost/serialization/array.hpp>
+#include <boost/serialization/base_object.hpp>
+
+#include <array>
+
+namespace cbm::algo::ca
+{
+  /// \enum  EDummy
+  /// \brief Dummy enum, representing an array with zero elements
+  ///
+  /// This enumeration is to be used as the default template parameter
+  enum class EDummy
+  {
+    END
+  };
+
+  /// \class cbm::algo::ca::EnumArray
+  /// \brief Class of arrays, which can be accessed by an enum class entry as an index
+  /// \note  The enum-array must contain an entry END, which represents the number of the enumeration entries and
+  ///        is used, as an array size
+  /// \tparam  E  The enum class type
+  /// \tparam  T  Type of data in the underlying array
+  template<class E, class T>
+  class EnumArray : public std::array<T, static_cast<std::size_t>(E::END)> {
+    using U   = typename std::underlying_type<E>::type;  ///< Underlying type of enumeration
+    using Arr = std::array<T, static_cast<std::size_t>(E::END)>;
+
+   public:
+    /// \brief Mutable access operator, indexed by enum entry
+    T& operator[](const E entry)
+    {
+      return std::array<T, static_cast<std::size_t>(E::END)>::operator[](static_cast<U>(entry));
+    }
+
+    /// \brief Mutable access operator, indexed by underlying index type
+    T& operator[](U index) { return std::array<T, static_cast<std::size_t>(E::END)>::operator[](index); }
+
+    /// \brief Constant access operator, indexed by enum entry
+    const T& operator[](const E& entry) const
+    {
+      return std::array<T, static_cast<std::size_t>(E::END)>::operator[](static_cast<U>(entry));
+    }
+
+    /// \brief Constant access operator, indexed by underlying index type
+    /// \param index  Index of the element
+    const T& operator[](U index) const { return std::array<T, static_cast<std::size_t>(E::END)>::operator[](index); }
+
+    /// \brief Convertion operator to the base array class
+    operator Arr&() const { return *this; }
+
+   private:
+    friend class boost::serialization::access;
+    template<typename Archive>
+    void serialize(Archive& ar, const unsigned int /*version*/)
+    {
+      ar& boost::serialization::base_object<Arr>(*this);
+    }
+  };
+}  // namespace cbm::algo::ca

algo/ca/core/utils/CaMonitor.h

+/* Copyright (C) 2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+   SPDX-License-Identifier: GPL-3.0-only
+   Authors: Sergei Zharko [committer] */
+
+/// \file   CaMonitor.h
+/// \brief  CA Tracking monitor class
+/// \since  25.08.2023
+/// \author S.Zharko <s.zharko@gsi.de>
+
+#pragma once  // include this header only once per compilation unit
+
+#include "CaEnumArray.h"
+#include "CaMonitorData.h"
+
+#include <boost/serialization/access.hpp>
+#include <boost/serialization/string.hpp>
+#include <boost/serialization/vector.hpp>
+
+#include <algorithm>
+#include <iomanip>
+#include <sstream>
+#include <string>
+#include <string_view>
+#include <vector>
+
+
+namespace cbm::algo::ca
+{
+  /// \class  Monitor
+  /// \brief  Monitor class for the CA tracking
+  /// \tparam ECounterKey  A enum class, containing keys for monitorables
+  /// \tparam ETimerKey    A enum class, containing keys for timers
+  ///
+  template<class ECounterKey, class ETimerKey = EDummy>
+  class Monitor {
+   public:
+    /// \brief Alias to array, indexed by the monitorable enumeration
+    template<typename T>
+    using CounterArray = EnumArray<ECounterKey, T>;
+
+    /// \brief Alias to array, indexed by the timer enumeration
+    template<typename T>
+    using TimerArray = EnumArray<ETimerKey, T>;
+
+    /// \brief Constructor
+    /// \param name Name of monitor
+    Monitor(std::string_view name) : fsName(name) {}
+
+    /// \brief Default constructor
+    Monitor() = default;
+
+    /// \brief Destructor
+    ~Monitor() = default;
+
+    /// \brief Copy constructor
+    Monitor(const Monitor&) = delete;
+
+    /// \brief Move constructor
+    Monitor(Monitor&&) = delete;
+
+    /// \brief Copy assignment operator
+    Monitor& operator=(const Monitor&) = delete;
+
+    /// \brief Move assignment operator
+    Monitor& operator=(Monitor&&) = delete;
+
+    /// \brief Adds the other monitor data to this
+    /// \param data     Reference to the other MonitorData object to add
+    /// \param parallel True, if the monitor data were obtained in parallel (See Timer::AddMonitorData)
+    void AddMonitorData(const MonitorData<ECounterKey, ETimerKey>& data, bool parallel = false)
+    {
+      fMonitorData.AddMonitorData(data, parallel);
+    }
+
+    /// \brief Gets counter name
+    /// \param key Counter key
+    const std::string& GetCounterName(ECounterKey key) const { return faCounterNames[key]; }
+
+    /// \brief Gets counter value
+    /// \param key
+    int GetCounterValue(ECounterKey key) const { return fMonitorData.GetCounterValue(key); }
+
+    /// \brief Gets monitor data
+    const MonitorData<ECounterKey, ETimerKey>& GetMonitorData() const { return fMonitorData; }
+
+    /// \brief Gets name of the monitor
+    const std::string& GetName() const { return fsName; }
+
+    /// \brief Gets timer
+    /// \param key Timer key
+    const Timer& GetTimer(ETimerKey key) const { return fMonitorData.GetTimer(key); }
+
+    /// \brief Gets timer name
+    /// \param key Timer key
+    const std::string& GetTimerName(ETimerKey key) const { return faTimerNames[key]; }
+
+    /// \brief Increments key counter by 1
+    /// \param key Counter key
+    void IncrementCounter(ECounterKey key) { fMonitorData.IncrementCounter(key); };
+
+    /// \brief Increments key counter by a number
+    /// \param key  Counter key
+    /// \param num  Number to add
+    void IncrementCounter(ECounterKey key, int num) { fMonitorData.IncrementCounter(key, num); }
+
+    /// \brief Resets the counters
+    void Reset() { fMonitorData.Reset(); }
+
+    /// \brief Sets keys of counters, which are used as denominators for ratios
+    /// \param vKeys Vector of keys
+    void SetRatioKeys(const std::vector<ECounterKey>& vKeys) { fvCounterRatioKeys = vKeys; }
+
+    /// \brief Sets name of counter
+    /// \param name  Name of counter
+    void SetCounterName(ECounterKey key, std::string_view name) { faCounterNames[key] = name; }
+
+    /// \brief Sets monitor data
+    /// \param data  The MonitorData class objet
+    void SetMonitorData(const MonitorData<ECounterKey, ETimerKey>& data) { fMonitorData = data; }
+
+    /// \brief Sets name of the monitor
+    /// \param name  Name of the monitor
+    void SetName(std::string_view name) { fsName = name; }
+
+    /// \brief Sets name of the timer
+    /// \param name  Name of the timer
+    void SetTimerName(ETimerKey key, std::string_view name) { faTimerNames[key] = name; }
+
+    /// \brief Starts timer
+    /// \param key  Timer key
+    void StartTimer(ETimerKey key) { fMonitorData.StartTimer(key); }
+
+    /// \brief Stops timer
+    /// \param key  Timer key
+    void StopTimer(ETimerKey key) { fMonitorData.StopTimer(key); }
+
+    /// \brief Prints counters summary to string
+    std::string ToString() const;
+
+    friend class boost::serialization::access;
+    template<typename Archive>
+    void serialize(Archive& ar, const unsigned int /*version*/)
+    {
+      ar& fMonitorData;
+      ar& faTimerNames;
+      ar& faCounterNames;
+      ar& fvCounterRatioKeys;
+      ar& fsName;
+    }
+
+   private:
+    MonitorData<ECounterKey, ETimerKey> fMonitorData;  ///< Monitor data
+    TimerArray<std::string> faTimerNames{};            ///< Array of timer names
+    CounterArray<std::string> faCounterNames{};        ///< Array of counter names
+    std::vector<ECounterKey> fvCounterRatioKeys{};     ///< List of keys, which are used as denominators in ratios
+    std::string fsName;                                ///< Name of the monitor
+  };
+
+
+  // *****************************************
+  // **  Template function implementations  **
+  // *****************************************
+
+  // ---------------------------------------------------------------------------------------------------------------------
+  //
+  template<class ECounterKey, class ETimerKey>
+  std::string Monitor<ECounterKey, ETimerKey>::ToString() const
+  {
+    using std::left;
+    using std::right;
+    using std::setfill;
+    using std::setw;
+    std::stringstream msg;
+    constexpr size_t width    = 24;
+    auto Cmp                  = [](const std::string& l, const std::string& r) { return l.size() < r.size(); };
+    msg << "\n===== Monitor: " << fsName << "\n";
+    if constexpr (!std::is_same_v<ETimerKey, EDummy>) {
+      size_t widthKeyTimer{std::max_element(faTimerNames.begin(), faTimerNames.end(), Cmp)->size() + 1};
+      msg << "\n----- Timers:\n";
+      msg << setw(widthKeyTimer) << left << "Key" << ' ';
+      msg << setw(width) << left << "N Calls" << ' ';
+      msg << setw(width) << left << "Average [s]" << ' ';
+      msg << setw(width) << left << "Min [s]" << ' ';
+      msg << setw(width) << left << "Max [s]" << ' ';
+      msg << setw(width) << left << "Total [s]" << '\n';
+      msg << right;
+      for (int iKey = 0; iKey < fMonitorData.GetNofTimers(); ++iKey) {
+        const Timer& timer = fMonitorData.GetTimer(static_cast<ETimerKey>(iKey));
+        msg << setw(widthKeyTimer) << faTimerNames[iKey] << ' ';
+        msg << setw(width) << timer.GetNofCalls() << ' ';
+        msg << setw(width) << timer.GetAverage() << ' ';
+        msg << setw(width) << timer.GetMin() << ' ';
+        msg << setw(width) << timer.GetMax() << ' ';
+        msg << setw(width) << timer.GetTotal() << '\n';
+      }
+    }
+
+    msg << "\n----- Counters:\n";
+    size_t widthKeyCounter{std::max_element(faCounterNames.begin(), faCounterNames.end(), Cmp)->size() + 1};
+    msg << setw(widthKeyCounter) << left << "Key" << ' ';
+    msg << setw(width) << left << "Total" << ' ';
+    for (auto key : fvCounterRatioKeys) {
+      msg << setw(width) << left << std::string("per ") + faCounterNames[key] << ' ';
+    }
+    msg << '\n';
+    for (int iKey = 0; iKey < fMonitorData.GetNofCounters(); ++iKey) {
+      auto counterValue = fMonitorData.GetCounterValue(static_cast<ECounterKey>(iKey));
+      msg << setw(widthKeyCounter) << left << faCounterNames[iKey] << ' ';
+      msg << setw(width) << right << counterValue << ' ';
+      for (auto keyDen : fvCounterRatioKeys) {
+        auto ratio = static_cast<double>(counterValue) / fMonitorData.GetCounterValue(keyDen);
+        msg << setw(width) << right << ratio << ' ';
+      }
+      msg << '\n';
+    }
+    return msg.str();
+  }
+}  // namespace cbm::algo::ca

algo/ca/core/utils/CaMonitorData.h

+/* Copyright (C) 2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+   SPDX-License-Identifier: GPL-3.0-only
+   Authors: Sergei Zharko [committer] */
+
+/// \file   CaMonitorData.h
+/// \brief  A block of data for cbm::algo::ca::Monitor
+/// \since  19.10.2023
+/// \author S.Zharko <s.zharko@gsi.de>
+
+#pragma once
+
+#include "CaEnumArray.h"
+#include "CaTimer.h"
+
+#include <boost/serialization/access.hpp>
+
+#include <algorithm>
+
+namespace cbm::algo::ca
+{
+  /// \class MonitorData
+  /// \brief Monitor data block
+  /// \tparam ECounterKey  A enum class, containing keys for monitorables
+  /// \tparam ETimerKey    A enum class, containing keys for timers
+  ///
+  template<class ECounterKey, class ETimerKey>
+  class MonitorData {
+   public:
+    /// \brief Alias to array, indexed by the counter enumeration
+    template<typename T>
+    using CounterArray = EnumArray<ECounterKey, T>;
+
+    /// \brief Alias to array, indexed by the timer enumeration
+    template<typename T>
+    using TimerArray = EnumArray<ETimerKey, T>;
+
+    /// \brief Default constructor
+    MonitorData() = default;
+
+    /// \brief Copy constructor
+    MonitorData(const MonitorData&) = default;
+
+    /// \brief Move constructor
+    MonitorData(MonitorData&&) = default;
+
+    /// \brief Destructor
+    ~MonitorData() = default;
+
+    /// \brief Copy assignment operator
+    MonitorData& operator=(const MonitorData&) = default;
+
+    /// \brief Move assignment operator
+    MonitorData& operator=(MonitorData&&) = default;
+
+    /// \brief Adds the other monitor data to this
+    /// \param other    Reference to the other MonitorData object to add
+    /// \param parallel If the monitors were filled in parallel (See CaTimer::AddTimer)
+    void AddMonitorData(const MonitorData& other, bool parallel = false);
+
+    /// \brief Gets counter value
+    /// \param key
+    int GetCounterValue(ECounterKey key) const { return faCounters[key]; }
+
+    /// \brief Gets number of counters
+    int GetNofCounters() const { return static_cast<int>(ECounterKey::END); }
+
+    /// \brief Gets number of timers
+    int GetNofTimers() const { return static_cast<int>(ETimerKey::END); }
+
+    /// \brief Gets timer
+    const Timer& GetTimer(ETimerKey key) const { return faTimers[key]; }
+
+    /// \brief Increments key counter by 1
+    /// \param key Counter key
+    void IncrementCounter(ECounterKey key) { ++faCounters[key]; };
+
+    /// \brief Increments key counter by a number
+    /// \param key  Counter key
+    /// \param num  Number to add
+    void IncrementCounter(ECounterKey key, int num) { faCounters[key] += num; }
+
+    /// \brief Resets all the counters and timers
+    void Reset();
+
+    /// \brief Resets a particular counter
+    /// \param key  Counter key
+    void ResetCounter(ECounterKey key) { faCounters[key] = 0; }
+
+    /// \brief Starts timer
+    /// \param key  Timer key
+    void StartTimer(ETimerKey key) { faTimers[key].Start(); }
+
+    /// \brief Stops timer
+    /// \param key  Timer key
+    void StopTimer(ETimerKey key) { faTimers[key].Stop(); }
+
+   private:
+    friend class boost::serialization::access;
+    template<typename Archive>
+    void serialize(Archive& ar, const unsigned int /*version*/)
+    {
+      ar& faTimers;
+      ar& faCounters;
+    }
+
+    TimerArray<Timer> faTimers{};    ///< Array of timers
+    CounterArray<int> faCounters{};  ///< Array of counters
+  };
+
+
+  // *****************************************
+  // **  Template function implementations  **
+  // *****************************************
+
+  // ---------------------------------------------------------------------------------------------------------------------
+  //
+  template<class ECounterKey, class ETimerKey>
+  inline void MonitorData<ECounterKey, ETimerKey>::AddMonitorData(const MonitorData<ECounterKey, ETimerKey>& other,
+                                                                  bool parallel)
+  {
+    for (size_t iCounter = 0; iCounter < faCounters.size(); ++iCounter) {
+      faCounters[iCounter] += other.faCounters[iCounter];
+    }
+    for (size_t iTimer = 0; iTimer < faTimers.size(); ++iTimer) {
+      faTimers[iTimer].AddTimer(other.faTimers[iTimer], parallel);
+    }
+  }
+
+  // ---------------------------------------------------------------------------------------------------------------------
+  //
+  template<class ECounterKey, class ETimerKey>
+  inline void MonitorData<ECounterKey, ETimerKey>::Reset()
+  {
+    faCounters.fill(0);
+    std::for_each(faTimers.begin(), faTimers.end(), [](auto& timer) { timer.Reset(); });
+  }
+}  // namespace cbm::algo::ca

algo/ca/core/utils/CaObjectInitController.h

+/* Copyright (C) 2022 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+   SPDX-License-Identifier: GPL-3.0-only
+   Authors: Sergey Gorbunov, Sergei Zharko [committer] */
+
+#pragma once  // include this header only once per compilation unit
+
+/// @file   CaObjectInitController.h
+/// @author Sergei Zharko
+/// @date   22.02.2022
+
+#include "AlgoFairloggerCompat.h"
+
+#include <bitset>
+#include <cassert>
+#include <iomanip>
+#include <iostream>
+#include <sstream>
+#include <string>
+
+namespace cbm::algo::ca
+{
+  /// ObjectInitController is a class, which provides flags system
+  /// and functionality needed for L1 algorithm objects initialization
+  ///
+  /// ObjectInitController is a class, which provides flags system
+  /// and functionality needed for L1 algorithm objects initialization
+  ///
+  // TODO: Possible improvement: introduce another template parameter, which represents a local enum class...
+  template<int NumberOfFlags, class InitKeyEnum>
+  class ObjectInitController {
+   public:
+    /// Gets an initialization status of the flag placed at bitIndex
+    /// \param bitIndex    index of bit
+    bool GetFlag(InitKeyEnum bitKey) const;
+    /// Checks, if the object is finalized, i.e. all its fields were set up
+    bool IsFinalized() const { return fInitFlags.size() == fInitFlags.count(); }
+    /// Sets an initialization status of the flag placed at bitIndex
+    /// \param bitIndex    index of bit
+    /// \param newStatus   flag value (true is default)
+    void SetFlag(InitKeyEnum bitKey, bool newStatus = true);
+    /// String representation of initialization flags contents
+    /// \param indentLevel   number of indent charachets int output
+    std::string ToString(int indentLevel = 0) const;
+
+   private:
+    std::bitset<NumberOfFlags> fInitFlags{};  ///< object of flags sets
+  };
+
+  //
+  //------------------------------------------------------------------------------------------------------------------------
+  //
+  template<int NumberOfFlags, class InitKeyEnum>
+  bool ObjectInitController<NumberOfFlags, InitKeyEnum>::GetFlag(InitKeyEnum bitKey) const
+  {
+    int bitIndex = static_cast<int>(bitKey);
+    if (bitIndex >= NumberOfFlags || bitIndex < 0) {
+      LOG(fatal) << "L1OnjectInitController::GetFlagStatus: attempt of flag access with index = " << bitIndex
+                 << " outside the range [0, " << NumberOfFlags << ']';
+      //assert((!(bitIndex >= NumberOfFlags || bitIndex < 0)));
+    }
+    return fInitFlags[static_cast<int>(bitKey)];
+  }
+  //
+  //------------------------------------------------------------------------------------------------------------------------
+  //
+  template<int NumberOfFlags, class InitKeyEnum>
+  void ObjectInitController<NumberOfFlags, InitKeyEnum>::SetFlag(InitKeyEnum bitKey, bool newStatus)
+  {
+    int bitIndex = static_cast<int>(bitKey);
+    if (bitIndex >= NumberOfFlags || bitIndex < 0) {
+      LOG(fatal) << "L1OnjectInitController::GetFlagStatus: attempt of flag access with index = " << bitIndex
+                 << " outside the range [0, " << NumberOfFlags << ']';
+      //assert((!(bitIndex >= NumberOfFlags || bitIndex < 0)));
+    }
+    fInitFlags[static_cast<int>(bitKey)] = newStatus;
+  }
+  //
+  //------------------------------------------------------------------------------------------------------------------------
+  //
+  template<int NumberOfFlags, class InitKeyEnum>
+  std::string ObjectInitController<NumberOfFlags, InitKeyEnum>::ToString(int indentLevel) const
+  {
+    std::stringstream aStream{};
+    constexpr char indentChar = '\t';
+    std::string indent(indentLevel, indentChar);
+    aStream << indent << "CaObjectInitController: flag values";
+    aStream << '\n' << indent << "index: ";
+    for (int idx = 0; idx < NumberOfFlags; ++idx) {
+      aStream << std::setw(3) << std::setfill(' ') << idx << ' ';
+    }
+    aStream << '\n' << indent << "value: ";
+    for (int idx = 0; idx < NumberOfFlags; ++idx) {
+      aStream << "  " << static_cast<int>(fInitFlags[idx]) << ' ';
+    }
+    return aStream.str();
+  }
+
+}  // namespace cbm::algo::ca

algo/ca/core/utils/CaSimd.h

+/* Copyright (C) 2010-2014 Frankfurt Institute for Advanced Studies, Goethe-Universitaet Frankfurt, Frankfurt
+   SPDX-License-Identifier: GPL-3.0-only
+   Authors: Sergey Gorbunov [committer]*/
+
+#pragma once  // include this header only once per compilation unit
+
+#include "KfSimd.h"
+#include "KfUtils.h"
+
+namespace cbm::algo::ca
+{
+  // Backward compatebility to ca::fvec etc.
+  using fvec  = kf::fvec;
+  using fscal = kf::fscal;
+  using fmask = kf::fmask;
+
+  // Utils namespace
+  namespace kfutils = cbm::algo::kf::utils;
+  namespace kfdefs  = cbm::algo::kf::defs;
+}  // namespace cbm::algo::ca

algo/ca/core/utils/CaTimer.h

+/* Copyright (C) 2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+   SPDX-License-Identifier: GPL-3.0-only
+   Authors: Sergei Zharko [committer] */
+
+/// \file   CaTimer.h
+/// \brief  Timer class for CA tracking (header)
+/// \since  18.10.2023
+/// \author S.Zharko <s.zharko@gsi.de>
+
+#pragma once
+
+#include <boost/serialization/serialization.hpp>
+
+#include <chrono>
+#include <cstdint>
+#include <iostream>
+#include <limits>
+#include <string>
+
+namespace cbm::algo::ca
+{
+  /// \class  Timer
+  /// \brief  A timer class for the monitor
+  ///
+  class Timer {
+   public:
+    using Clock         = std::chrono::high_resolution_clock;
+    using Duration      = std::chrono::nanoseconds;
+    using DurationCount = Duration::rep;
+    using TimePoint     = std::chrono::time_point<Clock, Duration>;
+
+    /// \brief Default constructor
+    Timer() = default;
+
+    /// \brief Destructor
+    ~Timer() = default;
+
+    /// \brief Copy constructor
+    Timer(const Timer&) = default;
+
+    /// \brief Move constructor
+    Timer(Timer&&) = default;
+
+    /// \brief Copy assignment operator
+    Timer& operator=(const Timer&) = default;
+
+    /// \brief Move assignment operator
+    Timer& operator=(Timer&&) = default;
+
+    /// \brief Adds another timer
+    /// \param other    Reference to the other Timer object to add
+    /// \param parallel Bool: if the timers were executed in parallel
+    ///
+    /// If the parallel flag is true then the resulting fTotal time is taken as a maximum of each total time of
+    /// the appended timers. If the parallel flag is false, the resulting fTotal is a sum of all timers.
+    void AddTimer(const Timer& other, bool parallel);
+
+    /// \brief Gets average time [s]
+    double GetAverage() const { return static_cast<double>(fTotal) / fNofCalls * 1.e-9; }
+
+    /// \brief Gets time of the longest call [s]
+    double GetMax() const { return static_cast<double>(fMax) * 1.e-9; }
+
+    /// \brief Gets time of the shortest call [s]
+    double GetMin() const { return static_cast<double>(fMin) * 1.e-9; }
+
+    /// \brief Gets number of calls
+    int GetNofCalls() const { return fNofCalls; }
+
+    /// \brief Gets index of the longest call
+    int GetMaxCallIndex() const { return fMaxCallIndex; }
+
+    /// \brief Gets index of the longest call
+    int GetMinCallIndex() const { return fMinCallIndex; }
+
+    /// \brief Gets total time [s]
+    double GetTotal() const { return static_cast<double>(fTotal) * 1.e-9; }
+
+    /// \brief Gets total time [ms]
+    double GetTotalMs() const { return GetTotal() * 1.e3; }
+
+    /// \brief Resets the timer
+    void Reset();
+
+    /// \brief Starts the timer
+    void Start() { fStart = Clock::now(); }
+
+    /// \brief Stops the timer
+    void Stop();
+
+   private:
+    friend class boost::serialization::access;
+    template<typename Archive>
+    void serialize(Archive& ar, const unsigned int /*version*/)
+    {
+      ar& fMin;
+      ar& fMax;
+      ar& fTotal;
+      ar& fNofCalls;
+      ar& fMinCallIndex;
+      ar& fMaxCallIndex;
+    }
+
+    TimePoint fStart     = TimePoint();
+    DurationCount fMin   = std::numeric_limits<DurationCount>::max();  ///< Minimal time period
+    DurationCount fMax   = std::numeric_limits<DurationCount>::min();  ///< Maximal time period
+    DurationCount fTotal = DurationCount(0);                           ///< Total measured time period [ns]
+    int fNofCalls        = 0;                                          ///< Number of timer calls [ns]
+    int fMinCallIndex    = -1;                                         ///< Index of the shortest call [ns]
+    int fMaxCallIndex    = -1;                                         ///< Index of the longest call [ns]
+  };                                                                   // class Timer
+
+
+  // ************************************
+  // **   Inline function definition   **
+  // ************************************
+
+  // -------------------------------------------------------------------------------------------------------------------
+  //
+  inline void Timer::AddTimer(const Timer& other, bool parallel)
+  {
+    if (other.fMin < fMin) {
+      fMin          = other.fMin;
+      fMinCallIndex = other.fMinCallIndex + fNofCalls;
+    }
+    if (other.fMax > fMax) {
+      fMax          = other.fMax;
+      fMaxCallIndex = other.fMinCallIndex + fNofCalls;
+    }
+    if (parallel) {
+      if (fTotal < other.fTotal) {
+        fTotal = other.fTotal;
+      }
+    }
+    else {
+      fTotal += other.fTotal;
+    }
+    fNofCalls += other.fNofCalls;
+  }
+
+  // -------------------------------------------------------------------------------------------------------------------
+  //
+  inline void Timer::Reset()
+  {
+    fStart        = TimePoint();
+    fMin          = std::numeric_limits<DurationCount>::max();
+    fMax          = std::numeric_limits<DurationCount>::min();
+    fMinCallIndex = -1;
+    fMaxCallIndex = -1;
+    fTotal        = DurationCount(0);
+    fNofCalls     = 0;
+  }
+
+  // -------------------------------------------------------------------------------------------------------------------
+  //
+  inline void Timer::Stop()
+  {
+    auto stop = Clock::now();
+    auto time = std::chrono::duration_cast<Duration>(stop - fStart).count();
+    if (fMin > time) {
+      fMin          = time;
+      fMinCallIndex = fNofCalls;
+    }
+    if (fMax < time) {
+      fMax          = time;
+      fMaxCallIndex = fNofCalls;
+    }
+    fTotal += time;
+    ++fNofCalls;
+  }
+}  // namespace cbm::algo::ca