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

algo/ca/core/utils/CaTrackingMonitor.h

+/* Copyright (C) 2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+   SPDX-License-Identifier: GPL-3.0-only
+   Authors: Sergei Zharko [committer] */
+
+/// \file   CaTrackingMonitor.h
+/// \brief  Monitor specialization for the tracking algorithm
+/// \since  19.10.2023
+/// \author S.Zharko <s.zharko@gsi.de>
+
+// NOTE: SZh: #prama once does not work properly in ROOT macros, so to use there enums from the header one
+//            have to use old approach to protect from multiple includes
+//#pragma once
+
+#ifndef CaTrackingMonitor_h
+#define CaTrackingMonitor_h 1
+
+#include "CaMonitor.h"
+
+#include <boost/serialization/base_object.hpp>
+
+namespace cbm::algo::ca
+{
+  /// \enum  ECounter
+  /// \brief Counter keys for the CA algo monitor
+  enum class ECounter
+  {
+    TrackingCall,  ///< number of the routine calls
+    SubTS,         ///< number of sub time-slices
+    RecoTrack,     ///< number of reconstructed tracks
+    RecoHit,       ///< number of reconstructed hits
+    RecoHitUsed,   ///< number of used reconstructed hits
+    Triplet,       ///< number of triplets
+    // TODO: Provide counters vs. detector ID
+    RecoMvdHit,        ///< number of MVD hits in tracks
+    RecoStsHit,        ///< number of STS hits in tracks
+    RecoMuchHit,       ///< number of MUCH hits in tracks
+    RecoTrdHit,        ///< number of TRD hits in tracks
+    RecoTofHit,        ///< number of TOF hits in tracks
+    UndefinedMvdHit,   ///< number of undefined MVD hits
+    UndefinedStsHit,   ///< number of undefined STS hits
+    UndefinedMuchHit,  ///< number of undefined MuCh hits
+    UndefinedTrdHit,   ///< number of undefined TRD hits
+    UndefinedTofHit,   ///< number of undefined TOF hits
+    END
+  };
+
+  /// \enum  ETimer
+  /// \brief Timer keys for the CA algo monitor
+  /* clang-format off */
+  // NOTE: SZh 21.03.2024: Disabling clang-format to indicate the scope of timers using the indent
+  enum class ETimer
+  {
+    TrackingChain,
+      PrepareInputData,
+        PrepareStsHits,
+        PrepareTrdHits,
+        PrepareTofHits,
+        InputDataTransmission,
+        CaHitCreation,
+      Tracking,
+        PrepareTimeslice,
+        TrackingThread,
+          PrepareThread,
+          PrepareWindow,
+          TrackingWindow,
+            InitWindow,
+            PrepareGrid,
+            FindTracks, /// (iterations loop)
+              PrepareIteration,
+              ConstructTriplets,
+              SearchNeighbours,
+              CreateTracks,
+              SuppressHitKeys,
+            FitTracks,
+            MergeClones,
+          StoreTracksWindow,
+        StoreTracksFinal,
+      Qa,
+    END
+  };
+  /* clang-format on */
+
+  using TrackingMonitorData = MonitorData<ECounter, ETimer>;
+
+  /// \class cbm::algo::ca::TrackingMonitor
+  /// \brief A monitor specialization for cbm::algo::ca::Framework class
+  class TrackingMonitor : public Monitor<ECounter, ETimer> {
+   public:
+    /// \brief Default constructor
+    TrackingMonitor() : Monitor<ECounter, ETimer>("CA Framework Monitor")
+    {
+      SetCounterName(ECounter::TrackingCall, "full routine calls");
+      SetCounterName(ECounter::RecoTrack, "reco tracks");
+      SetCounterName(ECounter::RecoHit, "reco hits");
+      SetCounterName(ECounter::Triplet, "triplets");
+      SetCounterName(ECounter::RecoHitUsed, "used reco hits");
+      SetCounterName(ECounter::SubTS, "sub-timeslices");
+      SetCounterName(ECounter::RecoMvdHit, "MVD hits in tracks");
+      SetCounterName(ECounter::RecoStsHit, "STS hits in tracks");
+      SetCounterName(ECounter::RecoMuchHit, "MUCH hits in tracks");
+      SetCounterName(ECounter::RecoTrdHit, "TRD hits in tracks");
+      SetCounterName(ECounter::RecoTofHit, "TOF hits in tracks");
+      SetCounterName(ECounter::UndefinedMvdHit, "undefined MVD hits");
+      SetCounterName(ECounter::UndefinedStsHit, "undefined STS hits");
+      SetCounterName(ECounter::UndefinedMuchHit, "undefined MuCh hits");
+      SetCounterName(ECounter::UndefinedTrdHit, "undefined TRD hits");
+      SetCounterName(ECounter::UndefinedTofHit, "undefined TOF hits");
+
+      SetTimerName(ETimer::TrackingChain, "tracking chain");
+      SetTimerName(ETimer::PrepareInputData, "input data preparation");
+      SetTimerName(ETimer::PrepareStsHits, "STS hits preparation");
+      SetTimerName(ETimer::PrepareTrdHits, "TRD hits preparation");
+      SetTimerName(ETimer::PrepareTofHits, "TOF hits preparation");
+      SetTimerName(ETimer::InputDataTransmission, "input data transmission");
+      SetTimerName(ETimer::CaHitCreation, "CA hit creation");
+      SetTimerName(ETimer::Tracking, "algorithm execution");
+      SetTimerName(ETimer::PrepareTimeslice, "timeslice preparation");
+      SetTimerName(ETimer::TrackingThread, "tracking on one thread");
+      SetTimerName(ETimer::PrepareThread, "thread preparation");
+      SetTimerName(ETimer::PrepareWindow, "window preparation");
+      SetTimerName(ETimer::TrackingWindow, "tracking in one window");
+      SetTimerName(ETimer::InitWindow, "window initialization");
+      SetTimerName(ETimer::PrepareGrid, "grid preparation");
+      SetTimerName(ETimer::FindTracks, "track finding");
+      SetTimerName(ETimer::PrepareIteration, "iteration preparation");
+      SetTimerName(ETimer::ConstructTriplets, "triplet constrcution");
+      SetTimerName(ETimer::SearchNeighbours, "neighbors search");
+      SetTimerName(ETimer::CreateTracks, "track creation");
+      SetTimerName(ETimer::SuppressHitKeys, "used hit key suppression");
+      SetTimerName(ETimer::FitTracks, "track fitter");
+      SetTimerName(ETimer::MergeClones, "clone merger");
+      SetTimerName(ETimer::StoreTracksWindow, "track storing in window");
+      SetTimerName(ETimer::StoreTracksFinal, "final track storing");
+      SetTimerName(ETimer::Qa, "QA");
+      SetRatioKeys({ECounter::TrackingCall, ECounter::SubTS, ECounter::RecoTrack});
+    }
+
+   private:
+    friend class boost::serialization::access;
+    template<typename Archive>
+    void serialize(Archive& ar, const unsigned int /*version*/)
+    {
+      ar& boost::serialization::base_object<Monitor<ECounter, ETimer>>(*this);
+    }
+  };
+
+}  // namespace cbm::algo::ca
+
+#endif

algo/ca/core/utils/CaUtils.cxx

+/* Copyright (C) 2022-2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+   SPDX-License-Identifier: GPL-3.0-only
+   Authors: Sergey Gorbunov, Sergei Zharko [committer] */
+
+#include "CaUtils.h"
+
+/// Namespace contains compile-time constants definition for the CA tracking algorithm
+///
+namespace cbm::algo::ca::utils
+{
+}  // namespace cbm::algo::ca::utils

algo/ca/core/utils/CaUtils.h

+/* Copyright (C) 2022-2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+   SPDX-License-Identifier: GPL-3.0-only
+   Authors: Sergey Gorbunov, Sergei Zharko [committer] */
+
+/// \file   CaUtils.h
+/// \brief  Compile-time constants definition for the CA tracking algorithm
+/// \since  02.06.2022
+/// \author S.Zharko <s.zharko@gsi.de>
+
+#pragma once  // include this header only once per compilation unit
+
+#include "CaHit.h"
+#include "KfDefs.h"
+#include "KfMeasurementTime.h"
+#include "KfMeasurementXy.h"
+#include "KfSimd.h"
+#include "KfTrackKalmanFilter.h"
+
+
+namespace cbm::algo::ca::utils
+{
+
+  template<typename T>
+  inline void FilterHit(kf::TrackKalmanFilter<T>& fit, const ca::Hit& hit, const kf::utils::masktype<T>& timeInfo)
+  {
+    kf::MeasurementXy<T> m;
+    m.SetDx2(hit.dX2());
+    m.SetDy2(hit.dY2());
+    m.SetDxy(hit.dXY());
+    m.SetX(hit.X());
+    m.SetY(hit.Y());
+    m.SetNdfX(T(1.));
+    m.SetNdfY(T(1.));
+    fit.FilterXY(m);
+    fit.FilterTime(hit.T(), hit.dT2(), timeInfo);
+  }
+
+}  // namespace cbm::algo::ca::utils

algo/ca/core/utils/CaVector.h

+/* Copyright (C) 2021-2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+   SPDX-License-Identifier: GPL-3.0-only
+   Authors: Sergey Gorbunov [committer], Sergei Zharko */
+
+/// \file CaVector.h
+/// \author Sergey Gorbunov
+/// \date 2021-06-16
+
+//#pragma once  // include this header only once per compilation unit (not robust)
+
+#ifndef CA_CORE_CaVector_h
+#define CA_CORE_CaVector_h 1
+
+#include "AlgoFairloggerCompat.h"
+
+#include <boost/serialization/access.hpp>
+#include <boost/serialization/base_object.hpp>
+#include <boost/serialization/string.hpp>
+#include <boost/serialization/vector.hpp>
+
+#include <memory>
+#include <sstream>
+
+namespace cbm::algo::ca
+{
+  /// \class Vector
+  ///
+  /// ca::Vector class is a wrapper around std::vector.
+  /// It does the following:
+  /// 1. gives names to vectors for better debugging
+  /// 2. controls the out-of-range access in debug mode
+  /// 3. supresses methods that are currently not controlled
+  /// 4. warns when slow memory operations are called,
+  ///    i.e. when the preallocated capacity is reached and the entire vector should be copied to a new place
+  /// 5. blocks usage of boolean vectors, as they have a special
+  ///    space-optimized but slow implementation in std. (Use ca::Vector<char> instead).
+  ///
+  template<class T>
+  class Vector : private std::vector<T> {
+    friend class boost::serialization::access;
+
+   public:
+    using Tbase                  = std::vector<T>;
+    using value_type             = T;
+    using allocator_type         = typename Tbase::allocator_type;
+    using pointer                = typename std::allocator_traits<allocator_type>::pointer;
+    using const_pointer          = typename std::allocator_traits<allocator_type>::const_pointer;
+    using reference              = value_type&;
+    using const_reference        = const value_type&;
+    using size_type              = typename Tbase::size_type;
+    using difference_type        = typename Tbase::difference_type;
+    using iterator               = typename Tbase::iterator;
+    using const_iterator         = typename Tbase::const_iterator;
+    using reverse_iterator       = std::reverse_iterator<iterator>;
+    using const_reverse_iterator = std::reverse_iterator<const_iterator>;
+
+
+    /// \brief Generic constructor from vairadic parameter list
+    template<typename... Tinput>
+    Vector(Tinput... value) : Tbase(value...)
+    {
+    }
+
+    /// \brief Generic constructor from vairadic parameter list including the name of the vector
+    template<typename... Tinput>
+    Vector(const char* name, Tinput... value) : Tbase(value...)
+                                              , fName(name)
+    {
+    }
+
+    /// \brief Constructor to make initializations like ca::Vector<int> myVector {"MyVector", {1, 2, 3}}
+    Vector(const std::string& name, std::initializer_list<T> init) : Tbase(init), fName(name) {}
+
+    /// \brief Copy constructor
+    Vector(const Vector& v) : Tbase() { *this = v; }
+
+    /// \brief Move constructor
+    Vector(Vector&& v) noexcept : Tbase(std::move(v)), fName(std::move(v.fName)) {}
+
+    /// \brief Copy assignment operator
+    Vector& operator=(const Vector& v)
+    {
+      if (this != &v) {
+        fName = v.fName;
+        Tbase::reserve(v.capacity());  // make sure that the capacity is transmitted
+        Tbase::assign(v.begin(), v.end());
+      }
+      return *this;
+    }
+
+    /// \brief Move assignment operator
+    Vector& operator=(Vector&& v) noexcept
+    {
+      if (this != &v) {
+        std::swap(fName, v.fName);
+        Tbase::swap(v);
+      }
+      return *this;
+    }
+
+    /// \brief Swap operator
+    void swap(Vector& v) noexcept
+    {
+      if (this != &v) {
+        std::swap(fName, v.fName);
+        Tbase::swap(v);
+      }
+    }
+
+    /// \brief Sets the name of the vector
+    /// \param s  Name of the vector
+    void SetName(const std::string& s) { fName = s; }
+
+    /// \brief Sets the name of the vector
+    /// \param s  Name of the vector (string stream)
+    void SetName(const std::basic_ostream<char>& s)
+    {
+      // helps to set a composed name in a single line via:
+      // SetName(std::stringstream()<<"my name "<<..whatever..);
+      fName = dynamic_cast<const std::stringstream&>(s).str();
+    }
+
+    /// \brief  Gets name of the vector
+    std::string GetName() const
+    {
+      std::string s = "ca::Vector<";
+      s += fName + "> ";
+      return s;
+    }
+
+    /// \brief Clears vector and resizes it to the selected size with selected values
+    /// \param count  New size of the vector
+    /// \param value  Variadic list of the parameters to pass to the base std::vector::resize function
+    template<typename... Tinput>
+    void reset(std::size_t count, Tinput... value)
+    {
+      // does the same as Tbase::assign(), but works with the default T constructor too
+      // (no second parameter)
+      Tbase::clear();
+      Tbase::resize(count, value...);
+    }
+
+    /// \brief Enlarges the vector to the new size
+    /// \param count  New size of the vector
+    /// \param value  Value of the new elements
+    template<typename... Tinput>
+    void enlarge(std::size_t count, Tinput... value)
+    {
+      if (count < Tbase::size()) {
+        LOG(fatal) << "ca::Vector \"" << fName << "\"::enlarge(" << count
+                   << "): the new size is smaller than the current one " << Tbase::size() << ", something goes wrong.";
+        assert(count >= Tbase::size());
+      }
+      if ((!Tbase::empty()) && (count > Tbase::capacity())) {
+        LOG(warning) << "ca::Vector \"" << fName << "\"::enlarge(" << count << "): allocated capacity of "
+                     << Tbase::capacity() << " is reached, the vector of size " << Tbase::size()
+                     << " will be copied to the new place.";
+      }
+      Tbase::resize(count, value...);
+    }
+
+    /// \brief Reduces the vector to a given size
+    /// \param count  Size of the new vector
+    void shrink(std::size_t count)
+    {
+      if (count > Tbase::size()) {
+        LOG(fatal) << "ca::Vector \"" << fName << "\"::shrink(" << count
+                   << "): the new size is bigger than the current one " << Tbase::size() << ", something goes wrong.";
+        assert(count < Tbase::size());
+      }
+      Tbase::resize(count);
+    }
+
+    /// \brief Reserves a new size for the vector
+    /// \param count  New size of the vector
+    void reserve(std::size_t count)
+    {
+      if (!Tbase::empty()) {
+        LOG(fatal) << "ca::Vector \"" << fName << "\"::reserve(" << count << "): the vector is not empty; "
+                   << " it will be copied to the new place.";
+        assert(Tbase::empty());
+      }
+      Tbase::reserve(count);
+    }
+
+    /// \brief Pushes back a value to the vector
+    /// \param value  New value
+    /// \note  Raises a warning, if the vector re-alocates memory
+    template<typename Tinput>
+    void push_back(Tinput value)
+    {
+      if (Tbase::size() >= Tbase::capacity()) {
+        LOG(warning) << "ca::Vector \"" << fName << "\"::push_back(): allocated capacity of " << Tbase::capacity()
+                     << " is reached, re-allocate and copy.";
+      }
+      Tbase::push_back(value);
+    }
+
+    /// \brief Pushes back a value to the vector without testing for the memory re-alocation
+    /// \param value  New value
+    template<typename Tinput>
+    void push_back_no_warning(Tinput value)
+    {
+      Tbase::push_back(value);
+    }
+
+    /// \brief Creates a parameter in the end of the vector
+    /// \param value  Variadic list of the parameters, which are passed to the element constructor
+    template<typename... Tinput>
+    void emplace_back(Tinput&&... value)
+    {
+      if (Tbase::size() >= Tbase::capacity()) {
+        LOG(warning) << "ca::Vector \"" << fName << "\"::emplace_back(): allocated capacity of " << Tbase::capacity()
+                     << " is reached, re-allocate and copy.";
+      }
+      Tbase::emplace_back(value...);
+    }
+
+    /// \brief Mutable access to the element by its index
+    /// \param pos  Index of the element
+    T& operator[](std::size_t pos)
+    {
+      if (pos >= Tbase::size()) {
+        LOG(fatal) << "ca::Vector \"" << fName << "\": trying to access element " << pos
+                   << " outside of the vector of the size of " << Tbase::size();
+        assert(pos < Tbase::size());
+      }
+      return Tbase::operator[](pos);
+    }
+
+    /// \brief Constant access to the element by its index
+    /// \param pos  Index of the element
+    const T& operator[](std::size_t pos) const
+    {
+      if (pos >= Tbase::size()) {
+        LOG(fatal) << "ca::Vector \"" << fName << "\": trying to access element " << pos
+                   << " outside of the vector of the size of " << Tbase::size();
+        assert(pos < Tbase::size());
+      }
+      return Tbase::operator[](pos);
+    }
+
+    /// \brief Mutable access to the last element of the vector
+    T& back()
+    {
+      if (Tbase::size() == 0) {
+        LOG(fatal) << "ca::Vector \"" << fName << "\": trying to access element of an empty vector";
+        assert(Tbase::size() > 0);
+      }
+      return Tbase::back();
+    }
+
+    /// \brief Constant access to the last element of the vector
+    const T& back() const
+    {
+      if (Tbase::size() == 0) {
+        LOG(fatal) << "ca::Vector \"" << fName << "\": trying to access element of an empty vector";
+        assert(Tbase::size() > 0);
+      }
+      return Tbase::back();
+    }
+
+    using Tbase::begin;
+    using Tbase::capacity;
+    using Tbase::cbegin;
+    using Tbase::cend;
+    using Tbase::clear;
+    using Tbase::end;
+    using Tbase::insert;  //TODO:: make it private
+    using Tbase::pop_back;
+    using Tbase::rbegin;
+    using Tbase::reserve;
+    using Tbase::shrink_to_fit;
+    using Tbase::size;
+
+
+   private:
+    std::string fName{"no name"};  ///< Name of the vector
+    using Tbase::assign;           // use reset() instead
+    using Tbase::at;
+    using Tbase::resize;
+
+    /// \brief Serialization function for the vector
+    template<class Archive>
+    void serialize(Archive& ar, const unsigned int /*version*/)
+    {
+      ar& boost::serialization::base_object<Tbase>(*this);
+      ar& fName;
+    }
+  };
+
+  ///
+  /// std::vector<bool> has a special implementation that is space-optimized
+  /// and therefore slow and not thread-safe.
+  /// That is why one should use ca::Vector<char> instead.
+  ///
+  template<>
+  class Vector<bool> {
+  };
+}  // namespace cbm::algo::ca
+
+#endif  // CA_CORE_CaVector_h

algo/ca/qa/CaQa.cxx

+/* Copyright (C) 2023-2025 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+   SPDX-License-Identifier: GPL-3.0-only
+   Authors: Sergei Zharko [committer] */
+
+/// \file   CaQa.cxx
+/// \date   20.11.2023
+/// \brief  A QA module for CA tracking (implementation)
+/// \author S.Zharko <s.zharko@gsi.de>
+
+#include "CaQa.h"
+
+#include "CaDefs.h"
+#include "CaInputData.h"
+#include "CaParameters.h"
+#include "CaTrack.h"
+#include "TrackingDefs.h"
+
+#include <algorithm>
+#include <fstream>
+#include <limits>
+
+#include <fmt/format.h>
+
+using cbm::algo::ca::Hit;
+using cbm::algo::ca::Qa;
+using cbm::algo::ca::constants::math::Pi;
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+bool Qa::CheckInit() const
+{
+  bool res = true;
+  if (!fpParameters) {
+    L_(error) << "cbm::algo::ca::OutputQa::Build(): parameters object is undefined";
+    res = false;
+  }
+  if (!fpInputData) {
+    L_(error) << "cbm::algo::ca::OutputQa::Build(): input data object is undefined";
+    res = false;
+  }
+  if (!fpvTracks) {
+    L_(error) << "cbm::algo::ca::OutputQa::Build(): track vector is undefined";
+    res = false;
+  }
+  if (!fpvRecoHits) {
+    L_(error) << "cbm::algo::ca::OutputQa::Build(): used hit index vector is undefined";
+    res = false;
+  }
+  return res;
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+void Qa::Init()
+{
+  using cbm::algo::qa::CanvasConfig;
+  using cbm::algo::qa::Data;
+  using cbm::algo::qa::H1D;
+  using cbm::algo::qa::H2D;
+  using cbm::algo::qa::PadConfig;
+  using cbm::algo::qa::Prof2D;
+  using fmt::format;
+
+  if (!IsActive()) {
+    return;
+  }
+
+  // ********************
+  // ** Hit distributions
+  int nSt = fpParameters->GetNstationsActive();  // Number of active tracking stations
+  {
+    // Occupancy distributions
+    fvphHitOccupXY.resize(nSt);
+    fvphHitOccupZX.resize(nSt);
+    if (kDebug) {
+      fvphHitOccupZY.resize(nSt);
+      fvphHitUsageXY.resize(nSt);
+    }
+
+    // Station sizes in transverse plane
+    std::vector<double> vMinX(nSt);
+    std::vector<double> vMaxX(nSt);
+    std::vector<double> vMinY(nSt);
+    std::vector<double> vMaxY(nSt);
+
+    int nBinsXY = 200;
+    int nBinsZ  = 400;
+
+    for (int iSt = 0; iSt < nSt; ++iSt) {
+      const auto& station = fpParameters->GetStation(iSt);
+      vMaxX[iSt]          = station.GetXmax<double>();
+      vMinX[iSt]          = station.GetXmin<double>();
+      vMaxY[iSt]          = station.GetYmax<double>();
+      vMinY[iSt]          = station.GetYmin<double>();
+      double dy           = (vMaxY[iSt] - vMinY[iSt]) * kXYZMargin;
+      double dx           = (vMaxX[iSt] - vMinX[iSt]) * kXYZMargin;
+      vMaxX[iSt] += dx;
+      vMinX[iSt] -= dx;
+      vMaxY[iSt] += dy;
+      vMinY[iSt] -= dy;
+    }
+    // Station max
+    double xMinA = *std::min_element(vMinX.begin(), vMinX.end());
+    double xMaxA = *std::max_element(vMaxX.begin(), vMaxX.end());
+    double yMinA = *std::min_element(vMinY.begin(), vMinY.end());
+    double yMaxA = *std::max_element(vMaxY.begin(), vMaxY.end());
+    double zMinA = fpParameters->GetStation(0).GetZ<double>();
+    double zMaxA = fpParameters->GetStation(nSt - 1).GetZ<double>();
+    {
+      double dz = (zMaxA - zMinA) * kXYZMargin;
+      zMinA -= dz;
+      zMaxA += dz;
+    }
+    for (int iSt = 0; iSt < nSt; ++iSt) {
+      int iStGeo           = fpParameters->GetActiveToGeoMap()[iSt];
+      auto [detID, iStLoc] = fpParameters->GetGeoToLocalIdMap()[iStGeo];
+      for (auto hitSet : kHitSets) {
+        auto setNm = EHitSet::Input == hitSet ? "input" : "used";
+        auto setTl = EHitSet::Input == hitSet ? "Input" : "Used";
+        {
+          auto name = format("hit_{}_occup_xy_sta_{}", setNm, iSt);
+          auto titl = format("{} hit occupancy in XY plane for station {} ({}{});x [cm];y [cm]", setTl, iSt,
+                             kDetName[detID], iStLoc);
+          fvphHitOccupXY[iSt][hitSet] =
+            MakeObj<H2D>(name, titl, nBinsXY, vMinX[iSt], vMaxX[iSt], nBinsXY, vMinY[iSt], vMaxY[iSt]);
+        }
+        {
+          auto name = format("hit_{}_occup_zx_sta_{}", setNm, iSt);
+          auto titl = format("{} hit occupancy in ZX plane for station {} ({}{});z [cm];x [cm]", setTl, iSt,
+                             kDetName[detID], iStLoc);
+          fvphHitOccupZX[iSt][hitSet] = MakeObj<H2D>(name, titl, nBinsZ, zMinA, zMaxA, nBinsXY, xMinA, xMaxA);
+        }
+        if (kDebug) {
+          auto name = format("hit_{}_occup_zy_sta_{}", setNm, iSt);
+          auto titl = format("{} hit occupancy in ZY plane for station {} ({}{});z [cm];y [cm]", setTl, iSt,
+                             kDetName[detID], iStLoc);
+          fvphHitOccupZY[iSt][hitSet] = MakeObj<H2D>(name, titl, nBinsZ, zMinA, zMaxA, nBinsXY, yMinA, yMaxA);
+        }
+      }
+      if (kDebug) {
+        auto name = format("hit_usage_xy_sta_{}", iSt);
+        auto titl = format("Hit usage in XY plane for station {} ({}{});x [cm];y [cm]", iSt, kDetName[detID], iStLoc);
+        fvphHitUsageXY[iSt] =
+          MakeObj<Prof2D>(name, titl, nBinsXY, vMinX[iSt], vMaxX[iSt], nBinsXY, vMinY[iSt], vMaxY[iSt], 0., 1.);
+      }
+    }
+    if (kDebug) {
+      for (auto hitSet : kHitSets) {
+        constexpr int NBins = 1000000;
+        auto setNm          = EHitSet::Input == hitSet ? "input" : "used";
+        auto setTl          = EHitSet::Input == hitSet ? "Input" : "Used";
+        {
+          auto name                    = format("hit_{}_front_key_index", setNm);
+          auto titl                    = format("{} hit front key index;ID_{{key}}/N_{{keys}};Count", setTl);
+          fvphHitFrontKeyIndex[hitSet] = MakeObj<H1D>(name, titl, NBins, 0., 1.);
+        }
+        {
+          auto name                   = format("hit_{}_back_key_index", setNm);
+          auto titl                   = format("{} hit back key index;ID_{{key}}/N_{{keys}};Count", setTl);
+          fvphHitBackKeyIndex[hitSet] = MakeObj<H1D>(name, titl, NBins, 0., 1.);
+        }
+      }
+    }
+
+    // Hit time distributions
+    if (kDebug) {
+      fvphHitTime.resize(nSt + 1);  // [nSt] - total over stations
+      for (int iSt = 0; iSt < nSt + 1; ++iSt) {
+        auto staNm = iSt == nSt ? "" : format("_sta_{}", iSt);
+        auto staTl = iSt == nSt ? "" : format(" in station {}", iSt);
+        for (auto hitSet : kHitSets) {
+          auto setNm               = EHitSet::Input == hitSet ? "input" : "used";
+          auto setTl               = EHitSet::Input == hitSet ? "Input" : "Used";
+          auto name                = format("hit_{}_rel_time{}", setNm, staNm);
+          auto titl                = format("{} hit relative time{}; #delta t_{{hit}};Count", setTl, staTl);
+          fvphHitTime[iSt][hitSet] = MakeObj<H1D>(name, titl, 10000, -0.1, 1.1);
+        }
+      }
+    }
+  }
+
+  // **********************
+  // ** Track distributions
+  std::vector<std::string> vsPointName = {"first", "last"};
+  for (int i = 0; i < knTrkParPoints; ++i) {
+    if (!kDebug && i > 0) {
+      continue;
+    }
+    {
+      auto sName      = format("{}_track_{}_theta", GetTaskName(), vsPointName[i]);
+      auto sTitl      = format("#theta at {} hit; #theta", vsPointName[i]);
+      fvphTrkTheta[i] = MakeObj<H1D>(sName, sTitl, 62, 0., 90.);
+    }
+    {
+      auto sName    = format("{}_track_{}_phi", GetTaskName(), vsPointName[i]);
+      auto sTitl    = format("#phi at {} hit; #phi", vsPointName[i]);
+      fvphTrkPhi[i] = MakeObj<H1D>(sName, sTitl, 62, -180., 180.);
+    }
+    {
+      auto sName         = format("{}_track_{}_thata_phi", GetTaskName(), vsPointName[i]);
+      auto sTitl         = format("#theta vs #phi at {} hit; #phi; #theta", vsPointName[i]);
+      fvphTrkPhiTheta[i] = MakeObj<H2D>(sName, sTitl, 62, -180., 180., 62, 0., 90.);
+    }
+    {
+      auto sName        = format("{}_track_{}_chi2_ndf", GetTaskName(), vsPointName[i]);
+      auto sTitl        = format("#chi^{{2}}/NDF at {} hit; #chi^{{2}}/NDF", vsPointName[i]);
+      fvphTrkChi2Ndf[i] = MakeObj<H1D>(sName, sTitl, 100, 0., 20.);
+    }
+  }
+  {
+    int nBins   = knStaMax;
+    double xMin = -0.5;
+    double xMax = double(knStaMax) - 0.5;
+    {
+      auto sName      = format("{}_track_fst_lst_sta", GetTaskName());
+      auto sTitl      = "First vs. last station index;ID^{last}_{station};ID^{first}_{station}";
+      fphTrkFstLstSta = MakeObj<H2D>(sName, sTitl, nBins, xMin, xMax, nBins, xMin, xMax);
+    }
+    {
+      auto sName     = format("{}_track_origin", GetTaskName());
+      auto sTitl     = "Track origin;x [cm];y [cm]";
+      fphTrkOriginXY = MakeObj<H2D>(sName, sTitl, kOriginB, kOriginL, kOriginU, kOriginB, kOriginL, kOriginU);
+    }
+    fphTrkNofHits = MakeObj<H1D>("n_hits", "Number of hits;N_{hit}", nBins, xMin, xMax);
+  }
+
+  // ---- Init canvases
+  {
+    // Input hits canvas
+    {
+      for (auto hitSet : kHitSets) {
+        auto setNm = EHitSet::Input == hitSet ? "input" : "used";
+        auto setTl = EHitSet::Input == hitSet ? "Input" : "Used";
+        {  // XY
+          auto name = format("{}_ca_hit_{}_occupancy_xy", GetTaskName(), setNm);
+          auto titl = format("{} hit occupancy in different stations in XY plane", setNm);
+          auto canv = CanvasConfig(name, titl);
+          for (int iSt = 0; iSt < nSt; ++iSt) {
+            auto pad = PadConfig(false, false, false, false, true);
+            pad.RegisterHistogram(fvphHitOccupXY[iSt][hitSet], "colz");
+            canv.AddPadConfig(pad);
+          }
+          AddCanvasConfig(canv);
+        }
+        {  // ZX and ZY
+          auto name = format("{}_ca_hit_{}_occupancy_zx_zy", GetTaskName(), setNm);
+          auto titl = format("{} hit occupancy in different stations in ZX and ZY planes", setTl);
+          auto canv = CanvasConfig(name, titl);
+          {  // ZX
+            auto pad = PadConfig();
+            for (int iSt = 0; iSt < nSt; ++iSt) {
+              pad.RegisterHistogram(fvphHitOccupZX[iSt][hitSet], (iSt == 0 ? "colz" : "cols same"));
+            }
+            canv.AddPadConfig(pad);
+          }
+          if (kDebug) {  // ZY
+            auto pad = PadConfig();
+            for (int iSt = 0; iSt < nSt; ++iSt) {
+              pad.RegisterHistogram(fvphHitOccupZY[iSt][hitSet], (iSt == 0 ? "colz" : "cols same"));
+            }
+            canv.AddPadConfig(pad);
+          }
+          AddCanvasConfig(canv);
+        }
+      }
+      if (kDebug) {
+        auto name = format("{}_ca_hit_usage_xy", GetTaskName());
+        auto titl = format("Hit usage in different stations in XY plane");
+        auto canv = CanvasConfig(name, titl);
+        for (int iSt = 0; iSt < nSt; ++iSt) {
+          auto pad = PadConfig();
+          pad.RegisterHistogram(fvphHitUsageXY[iSt], "colz");
+          canv.AddPadConfig(pad);
+        }
+        AddCanvasConfig(canv);
+      }
+    }
+
+    // Tracks canvas
+    {
+      auto canv = CanvasConfig(format("{}_ca_tracks", GetTaskName()), "Tracking output QA");
+      {
+        auto pad = PadConfig(true, true, false, false, true);
+        pad.RegisterHistogram(fvphTrkPhiTheta[0], "colz");
+        canv.AddPadConfig(pad);
+      }
+      {
+        auto pad = PadConfig(true, true, false, false, false);
+        pad.RegisterHistogram(fvphTrkChi2Ndf[0], "hist");
+        canv.AddPadConfig(pad);
+      }
+      {
+        auto pad = PadConfig(true, true, false, true, false);
+        pad.RegisterHistogram(fphTrkNofHits, "hist");
+        canv.AddPadConfig(pad);
+      }
+      {
+        auto pad = PadConfig(true, true, false, false, true);
+        pad.RegisterHistogram(fphTrkFstLstSta, "colz");
+        canv.AddPadConfig(pad);
+      }
+      {
+        auto pad = PadConfig(true, true, false, false, false);
+        pad.RegisterHistogram(fphTrkOriginXY, "colz");
+        canv.AddPadConfig(pad);
+      }
+      AddCanvasConfig(canv);
+    }
+  }
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+void Qa::Exec()
+{
+  if (!IsActive()) {
+    return;
+  }
+
+  if (!CheckInit()) {
+    L_(fatal) << "ca::Qa: instance is not initialized";
+    assert(false);
+  }
+
+  int nHitsInput = fpInputData->GetNhits();
+  // Map: if hit used in tracking
+  std::vector<unsigned char> vbHitUsed(nHitsInput, false);
+  for (int iH : (*fpvRecoHits)) {
+    vbHitUsed[iH] = true;
+  }
+
+  // Calculate max and min hit time
+  if constexpr (kDebug) {
+    const auto& hits = fpInputData->GetHits();
+    auto [minTimeIt, maxTimeIt] =
+      std::minmax_element(hits.cbegin(), hits.cend(), [](const auto& h1, const auto& h2) { return h1.T() < h2.T(); });
+    fMinHitTime = minTimeIt->T();
+    fMaxHitTime = maxTimeIt->T();
+  }
+
+  // Fill input hit histograms
+  {
+    for (int iH = 0; iH < nHitsInput; ++iH) {
+      const auto& hit = fpInputData->GetHit(iH);
+      FillHitDistributionsForHitSet(EHitSet::Input, hit);
+      if (vbHitUsed[iH]) {
+        FillHitDistributionsForHitSet(EHitSet::Used, hit);
+      }
+      if constexpr (kDebug) {
+        int iSt  = hit.Station();
+        double x = hit.X();
+        double y = hit.Y();
+        fvphHitUsageXY[iSt]->Fill(x, y, static_cast<double>(vbHitUsed[iH]));
+      }
+    }
+  }
+
+  // Fill track histograms
+  {
+    int trkFirstHit = 0;  // Index of hit in fpvRecoHits
+    for (auto trkIt = fpvTracks->begin(); trkIt != fpvTracks->end(); ++trkIt) {
+      auto& track = *trkIt;
+      int nHits   = track.fNofHits;
+      // Indices of hits in fpInputData->GetHits()
+      int iFstHit = (*fpvRecoHits)[trkFirstHit];
+      int iLstHit = (*fpvRecoHits)[trkFirstHit + nHits - 1];
+
+      // Distributions in different track points
+      for (int ip = 0; ip < knTrkParPoints; ++ip) {
+        if (!kDebug && ip > 0) {
+          continue;
+        }
+        //int iHit           = (ip == 0 ? iFstHit : iLstHit);
+        //const auto& hit    = fpInputData->GetHit(iHit);
+        const auto& trkPar = (ip == 0 ? track.fParFirst : track.fParLast);
+        fvphTrkTheta[ip]->Fill(trkPar.GetTheta() * 180. / Pi);
+        fvphTrkPhi[ip]->Fill(trkPar.GetPhi() * 180. / Pi);
+        fvphTrkPhiTheta[ip]->Fill(trkPar.GetPhi() * 180. / Pi, trkPar.GetTheta() * 180. / Pi);
+        fvphTrkChi2Ndf[ip]->Fill(trkPar.GetChiSq() / trkPar.GetNdf());
+      }
+      // Other distributions
+      fphTrkFstLstSta->Fill(fpInputData->GetHit(iLstHit).Station(), fpInputData->GetHit(iFstHit).Station());
+      fphTrkNofHits->Fill(nHits);
+      fphTrkOriginXY->Fill(track.fParPV.X(), track.fParPV.Y());
+      trkFirstHit += nHits;
+    }
+  }
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+void Qa::FillHitDistributionsForHitSet(Qa::EHitSet hitSet, const Hit& hit)
+{
+  int nSt  = fpParameters->GetNstationsActive();
+  int iSt  = hit.Station();
+  double x = hit.X();
+  double y = hit.Y();
+  double z = hit.Z();
+  fvphHitOccupXY[iSt][hitSet]->Fill(x, y);
+  fvphHitOccupZX[iSt][hitSet]->Fill(z, x);
+  if constexpr (kDebug) {
+    fvphHitOccupZY[iSt][hitSet]->Fill(z, y);
+    auto nKeys = static_cast<double>(fpInputData->GetNhitKeys());
+    fvphHitFrontKeyIndex[hitSet]->Fill(hit.FrontKey() / nKeys);
+    fvphHitBackKeyIndex[hitSet]->Fill(hit.BackKey() / nKeys);
+    double relTime = (hit.T() - fMinHitTime) / (fMaxHitTime - fMinHitTime);
+    fvphHitTime[iSt][hitSet]->Fill(relTime);
+    fvphHitTime[nSt][hitSet]->Fill(relTime);
+  }
+}

algo/ca/qa/CaQa.h

+/* Copyright (C) 2023-2025 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+   SPDX-License-Identifier: GPL-3.0-only
+   Authors: Sergei Zharko [committer] */
+
+/// \file   CaQa.h
+/// \date   20.11.2023
+/// \brief  A QA module for CA tracking (header)
+/// \author S.Zharko <s.zharko@gsi.de>
+
+#pragma once
+
+#include "CaEnumArray.h"
+#include "CaHit.h"  // for HitIndex_t
+#include "CaTimesliceHeader.h"
+#include "CaVector.h"
+#include "qa/QaTaskHeader.h"
+
+namespace cbm::algo
+{
+  namespace qa
+  {
+    class H1D;
+    class H2D;
+    class Manager;
+  }  // namespace qa
+  namespace ca
+  {
+    template<typename DataT>
+    class Parameters;
+    class InputData;
+    class Track;
+  }  // namespace ca
+}  // namespace cbm::algo
+
+namespace cbm::algo::ca
+{
+  /// \class cbm::algo::ca::qa::Qa
+  /// \brief Qa class for the CA tracking QA (header)
+  ///
+  class Qa : public qa::TaskHeader {
+    /// \brief Hit set entries
+    enum class EHitSet
+    {
+      Input,  ///< Input hits
+      Used,   ///< Hits used in tracks
+      END
+    };
+
+    /// \brief Definition of enum array over EHitSet entries
+    template<typename T>
+    using HitSetArray_t = EnumArray<EHitSet, T>;
+
+    /// \brief Array of EHitSet entries for iteration
+    static constexpr HitSetArray_t<EHitSet> kHitSets = {EHitSet::Input, EHitSet::Used};
+
+   public:
+    /// \brief Constructor
+    /// \param pManager  Pointer to the QA manager
+    /// \param name      Name of the QA (directory)
+    Qa(const std::unique_ptr<qa::Manager>& pManager, std::string_view name) : qa::TaskHeader(pManager, name) {}
+
+    /// \brief Constructor from the configuration object
+    /// \param config  QA configuration object
+    Qa() = default;
+
+    /// \brief Copy constructor
+    Qa(const Qa&) = delete;
+
+    /// \brief Move constructor
+    Qa(Qa&&) = delete;
+
+    /// \brief Destructor
+    ~Qa() = default;
+
+    /// \brief Copy assignment operator
+    Qa& operator=(const Qa&) = delete;
+
+    /// \brief Move assignment operator
+    Qa& operator=(Qa&&) = delete;
+
+    /// \brief QA execution function
+    void Exec();
+
+    /// \brief Initializes the QA
+    void Init();
+
+    /// \brief Check initialization
+    /// \return true   All variables are initialized
+    /// \return false  Some of are not initialized
+    bool CheckInit() const;
+
+    /// \brief Registers tracking input data object
+    /// \note  Call per TS
+    void RegisterInputData(const InputData* pInputData) { fpInputData = pInputData; }
+
+    /// \brief Registers track vector
+    /// \note  Call per TS
+    void RegisterTracks(const Vector<Track>* pvTracks) { fpvTracks = pvTracks; }
+
+    /// \brief Registers reco hits indices vector
+    /// \note  Call per TS
+    void RegisterRecoHitIndices(const Vector<HitIndex_t>* pvRecoHits) { fpvRecoHits = pvRecoHits; }
+
+    /// \brief Registers tracking parameters object
+    /// \note  Call per run
+    void RegisterParameters(const Parameters<fvec>* pParameters) { fpParameters = pParameters; }
+
+   private:
+    /// \brief Fills hit distributions
+    /// \param  hitSet  Hit set enum entry
+    /// \param  hit     Reference to hit
+    void FillHitDistributionsForHitSet(EHitSet hitSet, const ca::Hit& hit);
+
+    // parameters
+    static constexpr double kXYZMargin  = 0.05;   ///< Margin for occupancy distributions in XY plane
+    static constexpr int knHitSets      = 2;      ///< Number of hit sets: input/used
+    static constexpr int knTrkParPoints = 2;      ///< Number of track points to build par distributions
+    static constexpr int knStaMax       = 16;     ///< Max number of stations (histogram binning)
+    static constexpr bool kDebug        = false;  ///< Additional histograms
+    static constexpr int kOriginB       = 400;    ///< Track X(Y) at origin: n bins
+    static constexpr double kOriginL    = -1.;    ///< Track X(Y) at origin: lower bound [cm]
+    static constexpr double kOriginU    = +1.;    ///< Track X(Y) at origin: upper bound [cm]
+
+    double fMinHitTime = std::numeric_limits<double>::max();
+    double fMaxHitTime = std::numeric_limits<double>::lowest();
+
+    const Parameters<fvec>* fpParameters      = nullptr;  ///< Pointer to tracking parameters
+    const InputData* fpInputData              = nullptr;  ///< Pointer to input data
+    const Vector<Track>* fpvTracks            = nullptr;  ///< Pointer to tracks vector
+    const Vector<HitIndex_t>* fpvRecoHits     = nullptr;  ///< Pointer to reco hit indices
+
+    // Hit distributions
+    using OccupHistContainer_t = std::vector<HitSetArray_t<qa::H2D*>>;
+    OccupHistContainer_t fvphHitOccupXY;  ///< hist: Hit occupancy in different stations in XY plane
+    OccupHistContainer_t fvphHitOccupZX;  ///< hist: Hit occupancy in different stations in ZX plane
+    OccupHistContainer_t fvphHitOccupZY;  ///< hist: Hit occupancy in different stations in ZY plane
+
+    std::vector<qa::Prof2D*> fvphHitUsageXY;  ///< prof: Hit usage in different stations in XY plane
+
+    HitSetArray_t<qa::H1D*> fvphHitFrontKeyIndex = {nullptr, nullptr};  ///< Indices of front hit keys
+    HitSetArray_t<qa::H1D*> fvphHitBackKeyIndex  = {nullptr, nullptr};  ///< Indices of back hit keys
+
+    std::vector<HitSetArray_t<qa::H1D*>> fvphHitTime;  ///< Time distribution of hits
+
+    // Track distributions
+    std::array<qa::H1D*, knTrkParPoints> fvphTrkTheta    = {{0}};  ///< hist: theta at first/last hit
+    std::array<qa::H1D*, knTrkParPoints> fvphTrkPhi      = {{0}};  ///< hist: phi at first/last hit
+    std::array<qa::H1D*, knTrkParPoints> fvphTrkChi2Ndf  = {{0}};  ///< hist: chi2/NDF at first/last hit
+    std::array<qa::H2D*, knTrkParPoints> fvphTrkPhiTheta = {{0}};  ///< hist: theta vs. phi at first/last hit
+
+    qa::H2D* fphTrkOriginXY  = nullptr;  ///< hist: origin of tracks in x-y plane [cm]
+    qa::H2D* fphTrkFstLstSta = nullptr;  ///< hist: fst vs lst station index
+    qa::H1D* fphTrkNofHits   = nullptr;  ///< hist: number of hits in track
+  };
+}  // namespace cbm::algo::ca

algo/data/sts/Cluster.h

+/* Copyright (C) 2022 FIAS Frankfurt Institute for Advanced Studies, Frankfurt / Main
+   SPDX-License-Identifier: GPL-3.0-only
+   Authors: Felix Weiglhofer [committer]*/
+#ifndef CBM_ALGO_DATA_STS_CLUSTER_H
+#define CBM_ALGO_DATA_STS_CLUSTER_H
+
+#include "Definitions.h"
+
+#include <boost/serialization/access.hpp>
+
+namespace cbm::algo::sts
+{
+
+  struct Cluster {
+    real fCharge;         ///< Total charge
+    i32 fSize;            ///< Difference between first and last channel
+    real fPosition;       ///< Cluster centre in channel number units
+    real fPositionError;  ///< Cluster centre error (r.m.s.) in channel number units
+    u32 fTime;            ///< cluster time [ns]
+    real fTimeError;      ///< Error of cluster time [ns]
+
+   private:  // serialization
+    friend class boost::serialization::access;
+
+    template<class Archive>
+    void serialize(Archive& ar, unsigned int /*version*/)
+    {
+      ar& fCharge;
+      ar& fSize;
+      ar& fPosition;
+      ar& fPositionError;
+      ar& fTime;
+      ar& fTimeError;
+    }
+  };
+
+}  // namespace cbm::algo::sts
+
+#endif  // CBM_ALGO_DATA_STS_CLUSTER_H

algo/data/sts/Digi.h

+/* Copyright (C) 2022 FIAS Frankfurt Institute for Advanced Studies, Frankfurt / Main
+   SPDX-License-Identifier: GPL-3.0-only
+   Authors: Felix Weiglhofer [committer]*/
+#ifndef CBM_ALGO_DATA_STS_DIGI_H
+#define CBM_ALGO_DATA_STS_DIGI_H
+
+#include "CbmStsDigi.h"
+
+namespace cbm::algo::sts
+{
+
+  using Digi = CbmStsDigi;
+
+}  // namespace cbm::algo::sts
+
+#endif  // CBM_ALGO_DATA_STS_DIGI_H

algo/data/sts/Hit.h

+/* Copyright (C) 2023 FIAS Frankfurt Institute for Advanced Studies, Frankfurt / Main
+   SPDX-License-Identifier: GPL-3.0-only
+   Authors: Felix Weiglhofer [committer] */
+
+#ifndef CBM_ALGO_DATA_STS_HIT_H
+#define CBM_ALGO_DATA_STS_HIT_H
+
+#include "Definitions.h"
+
+#include <boost/serialization/access.hpp>
+
+namespace cbm::algo::sts
+{
+
+  struct Hit {
+    real fX, fY;      ///< X, Y positions of hit [cm]
+    real fZ;          ///< Z position of hit [cm]
+    u32 fTime;        ///< Hit time [ns]
+    real fDx, fDy;    ///< X, Y errors [cm]
+    real fDxy;        ///< XY correlation
+    real fTimeError;  ///< Error of hit time [ns]
+    real fDu;         ///< Error of coordinate across front-side strips [cm]
+    real fDv;         ///< Error of coordinate across back-side strips [cm]
+
+    u32 fFrontClusterId;  ///< Index of front-side cluster, used by tracking to reduce combinatorics
+    u32 fBackClusterId;   ///< Index of back-side cluster, used by tracking to reduce combinatorics
+
+    // Interface for tracking
+    double X() const { return fX; }
+    double Y() const { return fY; }
+    double Z() const { return fZ; }
+    double Time() const { return fTime; }
+    double Dx() const { return fDx; }
+    double Dy() const { return fDy; }
+    double TimeError() const { return fTimeError; }
+
+   private:  // serialization
+    friend class boost::serialization::access;
+
+    template<class Archive>
+    void serialize(Archive& ar, unsigned int /*version*/)
+    {
+      ar& fX;
+      ar& fY;
+      ar& fZ;
+      ar& fTime;
+      ar& fDx;
+      ar& fDy;
+      ar& fDxy;
+      ar& fTimeError;
+      ar& fDu;
+      ar& fDv;
+
+      ar& fFrontClusterId;
+      ar& fBackClusterId;
+    }
+  };
+
+}  // namespace cbm::algo::sts
+
+#endif  // CBM_ALGO_DATA_STS_HIT_H

algo/data/sts/HitfinderPars.cxx

+/* Copyright (C) 2024 FIAS Frankfurt Institute for Advanced Studies, Frankfurt / Main
+   SPDX-License-Identifier: GPL-3.0-only
+   Authors: Felix Weiglhofer [committer] */
+
+#include "HitfinderPars.h"
+
+using namespace cbm::algo;
+CBM_YAML_INSTANTIATE(sts::HitfinderPars);

algo/data/sts/HitfinderPars.h

+/* Copyright (C) 2023-2025 FIAS Frankfurt Institute for Advanced Studies, Frankfurt / Main
+   SPDX-License-Identifier: GPL-3.0-only
+   Authors: Felix Weiglhofer [committer] */
+#pragma once
+
+#include "CbmYaml.h"
+#include "Definitions.h"
+#include "LandauTable.h"
+
+#include <xpu/defines.h>
+
+namespace cbm::algo::sts
+{
+  struct HitfinderPars {
+    struct Asic {
+      int nAdc;
+      float dynamicRange;
+      float threshold;
+      float timeResolution;
+      float deadTime;
+      float noise;
+      float zeroNoiseRate;
+
+      XPU_D float AdcToCharge(unsigned short adc) const
+      {
+        return threshold + dynamicRange / float(nAdc) * (float(adc) + 0.5f);
+      }
+
+      CBM_YAML_PROPERTIES(
+        yaml::Property(&Asic::nAdc, "nAdc"),
+        yaml::Property(&Asic::dynamicRange, "dynamicRange"),
+        yaml::Property(&Asic::threshold, "threshold"),
+        yaml::Property(&Asic::timeResolution, "timeResolution"),
+        yaml::Property(&Asic::deadTime, "deadTime"),
+        yaml::Property(&Asic::noise, "noise"),
+        yaml::Property(&Asic::zeroNoiseRate, "zeroNoiseRate")
+      );
+    };
+
+    struct ModuleTransform {
+      // Rotation + translation matrix to transform
+      // local module coordinates into global coordinate system.
+      // No need for fancy math types here. These values are just copied
+      // and moved to the GPU.
+      // TODO: thats a lie, should use glm::mat3x4
+      std::array<float, 9> rotation;  // 3x3 matrix
+      std::array<float, 3> translation;
+
+      CBM_YAML_PROPERTIES(yaml::Property(&ModuleTransform::rotation, "rotation", "Rotation matrix", YAML::Flow),
+                          yaml::Property(&ModuleTransform::translation, "translation", "Translation vector",
+                                         YAML::Flow));
+    };
+
+    struct Module {
+      int32_t address;
+      float dY;
+      float pitch;
+      float stereoF;
+      float stereoB;
+      float lorentzF;
+      float lorentzB;
+      ModuleTransform localToGlobal;
+
+      CBM_YAML_PROPERTIES(yaml::Property(&Module::address, "address", "Hardware Address", YAML::Hex),
+                          yaml::Property(&Module::dY, "dY"), yaml::Property(&Module::pitch, "pitch"),
+                          yaml::Property(&Module::stereoF, "stereoF"), yaml::Property(&Module::stereoB, "stereoB"),
+                          yaml::Property(&Module::lorentzF, "lorentzF"), yaml::Property(&Module::lorentzB, "lorentzB"),
+                          yaml::Property(&Module::localToGlobal, "localToGlobal"));
+    };
+
+    Asic asic;
+    int nChannels;
+    std::vector<Module> modules;
+    LandauTable landauTable;  // Landau table for hitfinder, read from a seperate file
+
+    CBM_YAML_PROPERTIES(
+      yaml::Property(&HitfinderPars::asic, "asic",
+                     "Asic definitions. Currently assumes same parameters for all asics."),
+      yaml::Property(&HitfinderPars::nChannels, "nChannels",
+                     "Total number of channels per module. Hitfinder assumes nChannels / 2 channels per side."),
+      yaml::Property(&HitfinderPars::modules, "modules"));
+  };
+}  // namespace cbm::algo::sts
+
+CBM_YAML_EXTERN_DECL(cbm::algo::sts::HitfinderPars);

algo/data/sts/LandauTable.cxx

+/* Copyright (C) 2023 FIAS Frankfurt Institute for Advanced Studies, Frankfurt / Main
+   SPDX-License-Identifier: GPL-3.0-only
+   Authors: Felix Weiglhofer [committer] */
+#include "LandauTable.h"
+
+#include "AlgoFairloggerCompat.h"
+
+#include <fstream>
+
+using namespace cbm::algo;
+
+sts::LandauTable sts::LandauTable::FromFile(fs::path path)
+{
+  sts::LandauTable table;
+
+  std::vector<f32> charge;
+  std::vector<f32> prob;
+  std::ifstream file(path.string());
+
+  while (!file.eof()) {
+
+    f32 q, p;
+    file >> q >> p;
+    charge.push_back(q);
+    prob.push_back(p);
+
+    L_(trace) << "Reading Landau table " << path << " q=" << q << " p=" << p;
+  }
+
+  // TODO: check if charge is monotonically increasing, also more than 2 entries
+
+  table.stepSize = charge[1] - charge[0];
+  table.values   = std::move(prob);
+
+  return table;
+}

algo/data/sts/LandauTable.h

+/* Copyright (C) 2023 FIAS Frankfurt Institute for Advanced Studies, Frankfurt / Main
+   SPDX-License-Identifier: GPL-3.0-only
+   Authors: Felix Weiglhofer [committer] */
+#ifndef CBM_ALGO_DATA_STS_LANDAUTABLE_H
+#define CBM_ALGO_DATA_STS_LANDAUTABLE_H
+
+#include "Definitions.h"
+#include "compat/Filesystem.h"
+
+#include <vector>
+
+namespace cbm::algo::sts
+{
+
+  struct LandauTable {
+
+    static LandauTable FromFile(fs::path path);
+
+    std::vector<f32> values;
+    f32 stepSize = 0;
+  };
+
+}  // namespace cbm::algo::sts
+
+#endif  // CBM_ALGO_DATA_STS_LANDAUTABLE_H

algo/detectors/bmon/Calibrate.cxx

+/* Copyright (C) 2025 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+   SPDX-License-Identifier: GPL-3.0-only
+   Authors: Sergei Zharko [committer] */
+
+/// \file   Calibrate.h
+/// \brief  Calibrator for the BMON digis (implementation)
+/// \since  04.02.2025
+/// \author Sergei Zharko <s.zharko@gsi.de>
+
+#include "Calibrate.h"
+
+#include "AlgoFairloggerCompat.h"
+#include "CbmTofAddress.h"
+#include "util/TimingsFormat.h"
+
+#include <bitset>
+#include <chrono>
+
+using cbm::algo::bmon::Calibrate;
+using cbm::algo::bmon::CalibrateSetup;
+using fles::Subsystem;
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+Calibrate::Calibrate(CalibrateSetup setup) : fSetup(setup), fSelectionBitsOffset(0)
+{
+  // Initialize offset arrays for channel deadtime check
+  int nDiamondsInSetup = fSetup.diamonds.size();
+  if (nDiamondsInSetup == 0) {
+    throw std::runtime_error("No diamonds found in the BMON calibration config");
+  }
+  if (!(fSetup.selectionMask) != (nDiamondsInSetup == 1)) {
+    throw std::runtime_error("Wrong diamond selection mask: for a single diamond it must be zero, and for multiple"
+                             " diamonds it must be non-zero");
+  }
+  if (nDiamondsInSetup > 1) {
+    // Define the selection bit offset
+    while (!((fSetup.selectionMask >> fSelectionBitsOffset) % 2)) {
+      ++fSelectionBitsOffset;
+    }
+
+    // Sort the diamonds in the setup by their SM or Side or other distinguishing index
+    std::sort(fSetup.diamonds.begin(), fSetup.diamonds.end(), [&](const auto& lhs, const auto& rhs) {
+      return GetDiamondIndex(lhs.refAddress) < GetDiamondIndex(rhs.refAddress);
+    });
+  }
+
+  // Remove the channel information from the reference address
+  for (auto& diamond : fSetup.diamonds) {
+    diamond.refAddress &= ~CbmTofAddress::GetChannelIdBitmask();
+  }
+
+  // Calculate the channel offsets, needed for the dead time
+  fChannelOffset = std::vector<size_t>(nDiamondsInSetup + 1, 0);  // the last element -- total number of channels
+  for (int32_t iD = 0; iD < nDiamondsInSetup; ++iD) {
+    fChannelOffset[iD + 1] = fChannelOffset[iD] + fSetup.diamonds[iD].chanPar.size();
+  }
+  fChannelDeadTime = std::vector<double>(fChannelOffset.back(), std::numeric_limits<double>::quiet_NaN());
+
+  // **** DEBUG: END
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+Calibrate::resultType Calibrate::operator()(gsl::span<const CbmBmonDigi> digiIn)
+{
+  xpu::push_timer("BmonCalibrate");
+  xpu::t_add_bytes(digiIn.size_bytes());
+
+  // --- Output data
+  resultType result = {};
+
+  auto& calDigiOut = std::get<0>(result);
+  auto& monitor    = std::get<1>(result);
+  calDigiOut.reserve(digiIn.size());
+
+  // Reset the channel dead time
+  std::fill(fChannelDeadTime.begin(), fChannelDeadTime.end(), std::numeric_limits<double>::quiet_NaN());
+
+  const auto& diamonds = fSetup.diamonds;
+  for (const auto& digi : digiIn) {
+    uint32_t address = static_cast<uint32_t>(digi.GetAddress());
+    int32_t iChannel = CbmTofAddress::GetChannelId(address);
+    size_t iDiamond  = GetDiamondIndex(address);
+    if (iDiamond >= diamonds.size()
+        || (address & ~CbmTofAddress::GetChannelIdBitmask()) != diamonds[iDiamond].refAddress) {
+      monitor.fDigiCalibUnknownRPC++;
+      L_(error) << "Unknown BMON digi address: " << CbmTofAddress::ToString(address) << ", iDiamond = " << iDiamond;
+      continue;
+    }
+
+    const auto& diamondPar = diamonds[iDiamond];
+    const auto& channelPar = diamondPar.chanPar[iChannel];
+
+    // Check dead time
+    const size_t iChannelGlb = fChannelOffset[iDiamond] + iChannel;
+    const double deadTime    = fChannelDeadTime[iChannelGlb];
+
+    if (!std::isnan(deadTime) && digi.GetTime() <= deadTime) {
+      fChannelDeadTime[iChannelGlb] = digi.GetTime() + diamondPar.channelDeadtime;
+      monitor.fDigiDeadTimeCount++;
+      continue;
+    }
+    fChannelDeadTime[iChannelGlb] = digi.GetTime() + diamondPar.channelDeadtime;
+
+    // Create calibrated digi
+    CbmBmonDigi& pCalDigi = calDigiOut.emplace_back(digi);
+    pCalDigi.SetAddress(address);
+
+    // calibrate Digi Time
+    pCalDigi.SetTime(pCalDigi.GetTime() - channelPar.vCPTOff);
+
+    // subtract Offset
+    const double dTot = std::max(pCalDigi.GetCharge() - channelPar.vCPTotOff, 0.001);
+
+    // calibrate Digi charge (corresponds to TOF ToT)
+    pCalDigi.SetCharge(dTot * channelPar.vCPTotGain);
+
+    // walk correction
+    const std::vector<double>& walk = channelPar.vCPWalk;
+    const double dTotBinSize        = (diamondPar.TOTMax - diamondPar.TOTMin) / diamondPar.numClWalkBinX;
+    int32_t iWx                     = std::max((int32_t)((pCalDigi.GetCharge() - diamondPar.TOTMin) / dTotBinSize), 0);
+    iWx                             = std::min(iWx, diamondPar.numClWalkBinX - 1);
+
+    const double dDTot = (pCalDigi.GetCharge() - diamondPar.TOTMin) / dTotBinSize - (double) iWx - 0.5;
+    double dWT         = walk[iWx];
+
+    // linear interpolation to next bin
+    if (dDTot > 0) {
+      if (iWx < diamondPar.numClWalkBinX - 1) {
+        dWT += dDTot * (walk[iWx + 1] - walk[iWx]);
+      }
+    }
+    else {
+      if (0 < iWx) {
+        dWT -= dDTot * (walk[iWx - 1] - walk[iWx]);
+      }
+    }
+    pCalDigi.SetTime(pCalDigi.GetTime() - dWT);  // calibrate Digi Time
+  }
+
+  /// Sort the buffers of hits due to the time offsets applied
+  // (insert-sort faster than std::sort due to pre-sorting)
+  for (std::size_t i = 1; i < calDigiOut.size(); ++i) {
+    CbmTofDigi digi = calDigiOut[i];
+    std::size_t j   = i;
+    while (j > 0 && calDigiOut[j - 1].GetTime() > digi.GetTime()) {
+      calDigiOut[j] = calDigiOut[j - 1];
+      --j;
+    }
+    calDigiOut[j] = digi;
+  }
+
+  //  Kept for possible unsorted input
+  //  std::sort(calDigiOut.begin(), calDigiOut.end(),
+  //  [](const CbmTofDigi& a, const CbmTofDigi& b) -> bool { return a.GetTime() < b.GetTime(); });
+
+  monitor.fTime     = xpu::pop_timer();
+  monitor.fNumDigis = digiIn.size();
+  return result;
+}