algo/CMakeLists.txt

@@ -3,6 +3,7 @@ add_subdirectory(test)
 
 set(DEVICE_SRCS
   base/gpu/DeviceImage.cxx
+  base/gpu/Params.cxx
   detectors/sts/UnpackStsXpu.cxx
   detectors/sts/StsHitfinder.cxx
 )
@@ -10,12 +11,14 @@ set(DEVICE_SRCS
 set(SRCS
   ${DEVICE_SRCS}
   base/Options.cxx
+  base/util/TimingsFormat.cxx
+  data/sts/LandauTable.cxx
   evbuild/EventBuilder.cxx
-  global/Reco.cxx
   trigger/TimeClusterTrigger.cxx
   evselector/DigiEventSelector.cxx
   detectors/sts/StsHitfinderChain.cxx
   detectors/sts/StsReadoutConfig.cxx
+  detectors/sts/StsUnpackChain.cxx
   detectors/sts/UnpackSts.cxx
   detectors/much/MuchReadoutConfig.cxx
   detectors/much/UnpackMuch.cxx
@@ -26,6 +29,7 @@ set(SRCS
   detectors/bmon/UnpackBmon.cxx
   detectors/trd/TrdReadoutConfig.cxx
   detectors/trd/UnpackTrd.cxx
+  global/Reco.cxx
 )
 
 set(BUILD_INFO_CXX ${CMAKE_CURRENT_BINARY_DIR}/base/BuildInfo.cxx)
@@ -49,6 +53,7 @@ add_library(Algo SHARED ${SRCS})
 
 target_include_directories(Algo
   PUBLIC ${CMAKE_CURRENT_SOURCE_DIR}
+         ${CMAKE_CURRENT_SOURCE_DIR}/data
          ${CMAKE_CURRENT_SOURCE_DIR}/base
          ${CMAKE_CURRENT_SOURCE_DIR}/evbuild
          ${CMAKE_CURRENT_SOURCE_DIR}/global
@@ -67,6 +72,7 @@ target_link_libraries(Algo
             GSL
             Boost::program_options
             xpu
+            external::yaml-cpp
   INTERFACE FairLogger::FairLogger
             external::fles_ipc
 )
@@ -74,9 +80,8 @@ target_compile_definitions(Algo PUBLIC NO_ROOT)
 xpu_attach(Algo ${DEVICE_SRCS})
 
 install(TARGETS Algo DESTINATION lib)
-install(DIRECTORY base/gpu TYPE INCLUDE
-  FILES_MATCHING PATTERN "*.h"
-)
+install(DIRECTORY base/gpu TYPE INCLUDE FILES_MATCHING PATTERN "*.h")
+install(DIRECTORY data/sts TYPE INCLUDE FILES_MATCHING PATTERN "*.h")
 
 install(
   FILES       ca/simd/CaSimd.h

algo/base/ChainContext.h

@@ -5,11 +5,13 @@
 #define CBM_ALGO_BASE_CHAINCONTEXT_H
 
 #include "Options.h"
+#include "RecoParams.h"
 
 namespace cbm::algo
 {
   struct ChainContext {
     Options opts;
+    RecoParams recoParams;
   };
 }  // namespace cbm::algo
 

algo/base/SubChain.h

@@ -16,6 +16,7 @@ namespace cbm::algo
     void SetContext(ChainContext* ctx) { fContext = ctx; }
 
     const Options& Opts() const { return gsl::make_not_null(fContext)->opts; }
+    const RecoParams& Params() const { return gsl::make_not_null(fContext)->recoParams; }
 
   private:
     ChainContext* fContext = nullptr;

algo/base/config/Yaml.h

@@ -15,7 +15,7 @@
 #include "Prelude.h"
 #include "Property.h"
 
-namespace CORE_INTERFACE_NS::config
+namespace cbm::algo::config
 {
 
   template<typename T, T... Values, typename Func>
@@ -144,6 +144,6 @@ namespace CORE_INTERFACE_NS::config
     }
   };
 
-}  // namespace CORE_INTERFACE_NS::config
+}  // namespace cbm::algo::config
 
 #endif

algo/base/gpu/DeviceImage.h

@@ -7,7 +7,7 @@
 
 namespace cbm::algo
 {
-  struct GPUReco {
+  struct GPUReco : xpu::device_image {
   };
 }  // namespace cbm::algo
 

algo/base/gpu/Params.cxx

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

algo/base/gpu/Params.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_GPU_CONFIG_H
+#define CBM_ALGO_GPU_CONFIG_H
+
+#include <xpu/device.h>
+
+#include "DeviceImage.h"
+#include "Prelude.h"
+#include "RecoParams.h"
+
+namespace cbm::algo
+{
+  struct Params : xpu::constant<GPUReco, RecoParams> {
+  };
+}  // namespace cbm::algo
+
+#endif

algo/base/gpu/xpu_legacy.h

+/* Copyright (C) 2022 FIAS Frankfurt Institute for Advanced Studies, Frankfurt / Main
+   SPDX-License-Identifier: GPL-3.0-only
+   Authors: Felix Weiglhofer [committer]*/
+#ifndef CORE_COMPAT_XPU_LEGACY_H
+#define CORE_COMPAT_XPU_LEGACY_H
+
+#include <xpu/host.h>
+
+namespace xpu
+{
+
+  inline constexpr auto host_to_device = xpu::h2d;
+  inline constexpr auto device_to_host = xpu::d2h;
+
+  template<typename T>
+  class hd_buffer {
+
+  public:
+    hd_buffer() = default;
+    hd_buffer(size_t size) : m_buffer(size, xpu::buf_io) {}
+
+    T* h() { return xpu::h_view(m_buffer).begin(); }
+    T* d() { return m_buffer.get(); }
+
+    xpu::buffer<T>& underlying() { return m_buffer; }
+
+  private:
+    xpu::buffer<T> m_buffer;
+  };
+
+  template<typename T>
+  class d_buffer {
+
+  public:
+    d_buffer() = default;
+    d_buffer(size_t size) : m_buffer(size, xpu::buf_device) {}
+
+    T* d() { return m_buffer.get(); }
+
+    xpu::buffer<T>& underlying() { return m_buffer; }
+
+  private:
+    xpu::buffer<T> m_buffer;
+  };
+
+
+  template<typename T>
+  void copy(hd_buffer<T>& buf, direction dir)
+  {
+    static xpu::queue _Q;
+    _Q.copy(buf.underlying(), dir);
+    _Q.wait();
+  }
+
+  enum class side
+  {
+    host,
+    device
+  };
+
+  template<typename T, side S>
+  struct cmem_io {
+    using type = T*;
+  };
+
+  template<typename T>
+  struct cmem_io<T, side::host> {
+    using type = hd_buffer<T>;
+  };
+
+  template<typename T, side S>
+  using cmem_io_t = typename cmem_io<T, S>::type;
+
+  template<typename T, side S>
+  struct cmem_device {
+    using type = T*;
+  };
+
+  template<typename T>
+  struct cmem_device<T, side::host> {
+    using type = d_buffer<T>;
+  };
+
+  template<typename T, side S>
+  using cmem_device_t = typename cmem_device<T, S>::type;
+
+}  // namespace xpu
+
+#define XPU_BLOCK_SIZE_1D(...)
+
+#define XPU_EXPORT_KERNEL(Image, Kernel, ...) XPU_EXPORT_KERNEL_II(Image, Kernel, xpu::no_smem, 64, ##__VA_ARGS__)
+
+#define XPU_EXPORT_KERNEL_II(Image, Kernel, SMEM, BlockSize, ...)                                                      \
+  struct Kernel : xpu::kernel<Image> {                                                                                 \
+    using block_size = xpu::block_size<BlockSize>;                                                                     \
+    using context    = xpu::kernel_context<SMEM>;                                                                      \
+    XPU_D void operator()(context& ctx, ##__VA_ARGS__);                                                                \
+  }
+
+#define XPU_KERNEL(Kernel, smemIgnored, ...)                                                                           \
+  XPU_EXPORT(Kernel);                                                                                                  \
+  XPU_D void Kernel::operator()(context& ctx, ##__VA_ARGS__)
+
+#endif

algo/base/util/TimingsFormat.cxx

+/* Copyright (C) 2022 FIAS Frankfurt Institute for Advanced Studies, Frankfurt / Main
+   SPDX-License-Identifier: GPL-3.0-only
+   Authors: Felix Weiglhofer [committer]*/
+#include "TimingsFormat.h"
+
+#include <iomanip>
+#include <sstream>
+
+#include <fmt/format.h>
+#include <xpu/host.h>
+
+namespace cbm::algo
+{
+
+  class TimingsFormat {
+
+  public:
+    void Begin(size_t align)
+    {
+      fAlign = align;
+      fSS    = {};
+    }
+
+    void Title(std::string_view title) { fSS << fmt::format("{:<{}}\n", title, fAlign); }
+
+    std::string Finalize() { return fSS.str(); }
+
+    void Fmt(const xpu::timings& t)
+    {
+      fIndent += 2;
+      Measurement("Memcpy(h2d)", t.copy(xpu::h2d), t.throughput_copy(xpu::h2d));
+      NewLine();
+      Measurement("Memcpy(d2h)", t.copy(xpu::d2h), t.throughput_copy(xpu::d2h));
+      NewLine();
+      Measurement("Memset", t.memset(), t.throughput_memset());
+      NewLine();
+
+      for (xpu::timings& st : t.children()) {
+        if (st.kernels().empty()) {
+          Measurement(st.name(), st.wall(), st.throughput());
+          NewLine();
+        }
+        else {
+          Title(st.name());
+          Fmt(st);
+          NewLine();
+        }
+      }
+
+      for (xpu::kernel_timings& kt : t.kernels()) {
+        Measurement(KernelName(kt), kt.total(), kt.throughput());
+        NewLine();
+      }
+
+      Measurement("Kernel time", t.kernel_time(), t.throughput_kernels());
+      NewLine();
+      Measurement("Wall time", t.wall(), t.throughput());
+      fIndent -= 2;
+    }
+
+    void FmtSubtimers(const xpu::timings& t)
+    {
+      const auto subtimes = t.children();
+      for (auto it = subtimes.begin(); it != subtimes.end(); ++it) {
+        Title(it->name());
+        Fmt(*it);
+        if (std::next(it) != subtimes.end()) { NewLine(); }
+      }
+    }
+
+    void FmtSummary(const xpu::timings& t)
+    {
+      fIndent += 2;
+      Measurement("Memcpy(h2d)", t.copy(xpu::h2d), t.throughput_copy(xpu::h2d));
+      NewLine();
+      Measurement("Memcpy(d2h)", t.copy(xpu::d2h), t.throughput_copy(xpu::d2h));
+      NewLine();
+      Measurement("Memset", t.memset(), t.throughput_memset());
+      NewLine();
+      Measurement("Kernel time", t.kernel_time(), t.throughput_kernels());
+      NewLine();
+      Measurement("Wall time", t.wall(), t.throughput());
+      fIndent -= 2;
+    }
+
+    void NewLine() { fSS << "\n"; }
+
+  private:
+    void Measurement(std::string_view name, f64 time, f64 throughput)
+    {
+      fSS << std::setw(fIndent) << std::setfill(' ') << std::right << "";
+      fSS << std::setw(fAlign) << std::setfill(' ') << std::left << fmt::format("{}:", name);
+      Real(time, 10, 3, "ms");
+      if (std::isnormal(throughput)) {
+        fSS << " (";
+        Real(throughput, 7, 3, "GB/s");
+        fSS << ")";
+      }
+    }
+
+    void Real(double x, int width, int precision, std::string_view unit)
+    {
+      fSS << std::setw(width) << std::setfill(' ') << std::right << std::fixed << std::setprecision(precision) << x
+          << " " << unit;
+    }
+
+    std::string_view KernelName(xpu::kernel_timings& kt)
+    {
+      std::string prefix = "cbm::algo::";
+      if (kt.name().compare(0, prefix.size(), prefix) == 0) { return kt.name().substr(11); }
+      else {
+        return kt.name();
+      }
+    }
+
+    size_t fAlign  = 0;
+    size_t fIndent = 0;
+    std::stringstream fSS;
+
+  };  // class TimingsFormat
+
+  std::string MakeReport(std::string_view title, const xpu::timings& t, size_t align)
+  {
+    TimingsFormat tf;
+    tf.Begin(align);
+    tf.Title(title);
+    tf.Fmt(t);
+    return tf.Finalize();
+  }
+
+  std::string MakeReportSubtimers(std::string_view title, const xpu::timings& t, size_t align)
+  {
+    TimingsFormat tf;
+    tf.Begin(align);
+    tf.Title(title);
+    tf.FmtSubtimers(t);
+    return tf.Finalize();
+  }
+
+  std::string MakeReportSummary(std::string_view title, const xpu::timings& t, size_t align)
+  {
+    TimingsFormat tf;
+    tf.Begin(align);
+    tf.Title(title);
+    tf.FmtSummary(t);
+    return tf.Finalize();
+  }
+
+}  // namespace cbm::algo

algo/base/util/TimingsFormat.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_BASE_UTIL_TIMINGSFORMAT_H
+#define CBM_ALGO_BASE_UTIL_TIMINGSFORMAT_H
+
+#include <string>
+#include <string_view>
+
+#include "Prelude.h"
+
+namespace xpu
+{
+  class timings;
+}
+
+namespace cbm::algo
+{
+  std::string MakeReport(std::string_view title, const xpu::timings& t, size_t align = 40);
+  std::string MakeReportSubtimers(std::string_view title, const xpu::timings& t, size_t align = 40);
+  std::string MakeReportSummary(std::string_view, const xpu::timings& t, size_t align = 40);
+}  // namespace cbm::algo
+
+#endif

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 "Prelude.h"
+
+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]
+  };
+
+}  // namespace cbm::algo::sts
+
+#endif  // CBM_ALGO_DATA_STS_CLUSTER_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 "Prelude.h"
+
+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 fDz;         ///< Z position error [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]
+  };
+
+}  // namespace cbm::algo::sts
+
+#endif  // CBM_ALGO_DATA_STS_HIT_H

algo/data/sts/HitfinderPars.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_HITFINDERPARS_H
+#define CBM_ALGO_DATA_STS_HITFINDERPARS_H
+
+#include "LandauTable.h"
+#include "Prelude.h"
+#include "config/Property.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);
+      }
+
+      static constexpr auto Properties = std::make_tuple(
+        config::Property(&Asic::nAdc, "nAdc"), config::Property(&Asic::dynamicRange, "dynamicRange"),
+        config::Property(&Asic::threshold, "threshold"), config::Property(&Asic::timeResolution, "timeResolution"),
+        config::Property(&Asic::deadTime, "deadTime"), config::Property(&Asic::noise, "noise"),
+        config::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;
+
+      static constexpr auto Properties = std::make_tuple(
+        config::Property(&ModuleTransform::rotation, "rotation", "Rotation matrix", YAML::Flow),
+        config::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;
+
+      static constexpr auto Properties = std::make_tuple(
+        config::Property(&Module::address, "address", "Hardware Address", YAML::Hex),
+        config::Property(&Module::dY, "dY"), config::Property(&Module::pitch, "pitch"),
+        config::Property(&Module::stereoF, "stereoF"), config::Property(&Module::stereoB, "stereoB"),
+        config::Property(&Module::lorentzF, "lorentzF"), config::Property(&Module::lorentzB, "lorentzB"),
+        config::Property(&Module::localToGlobal, "localToGlobal"));
+    };
+
+    Asic asic;
+    int nChannels;
+    std::vector<Module> modules;
+    LandauTable landauTable;  // Landau table for hitfinder, read from a seperate file
+
+    static constexpr auto Properties = std::make_tuple(
+      config::Property(&HitfinderPars::asic, "asic",
+                       "Asic definitions. Currently assumes same parameters for all asics."),
+      config::Property(&HitfinderPars::nChannels, "nChannels",
+                       "Total number of channels per module. Hitfinder assumes nChannels / 2 channels per side."),
+      config::Property(&HitfinderPars::modules, "modules"));
+  };
+}  // namespace cbm::algo::sts
+
+#endif

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 <fairlogger/Logger.h>
+
+#include <fstream>
+
+using namespace cbm::algo;
+
+sts::LandauTable sts::LandauTable::FromFile(std::filesystem::path path)
+{
+  sts::LandauTable table;
+
+  std::vector<f32> charge;
+  std::vector<f32> prob;
+  std::ifstream file(path);
+
+  while (!file.eof()) {
+
+    f32 q, p;
+    file >> q >> p;
+    charge.push_back(q);
+    prob.push_back(p);
+
+    LOG(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 <filesystem>
+#include <vector>
+
+#include "Prelude.h"
+
+namespace cbm::algo::sts
+{
+
+  struct LandauTable {
+
+    static LandauTable FromFile(std::filesystem::path path);
+
+    std::vector<f32> values;
+    f32 stepSize = 0;
+  };
+
+}  // namespace cbm::algo::sts
+
+#endif  // CBM_ALGO_DATA_STS_LANDAUTABLE_H

algo/detectors/sts/StsHitfinderChain.h

@@ -14,46 +14,16 @@
 #include <optional>
 #include <vector>
 
-#include <xpu/host.h>
-
 #include "StsHitfinder.h"
+#include "SubChain.h"
+#include "gpu/xpu_legacy.h"
+#include "sts/HitfinderPars.h"
+#include "sts/LandauTable.h"
 
-namespace cbm::algo
+namespace cbm::algo::sts
 {
 
-  struct StsModuleTransformationMatrix {
-    // 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.
-    std::array<float, 9> rotation;  // 3x3 matrix
-    std::array<float, 3> translation;
-  };
-
-  struct StsModulePar {
-    int32_t address;
-    float dY;
-    float pitch;
-    float stereoF;
-    float stereoB;
-    float lorentzF;
-    float lorentzB;
-    StsModuleTransformationMatrix localToGlobal;
-  };
-
-  struct StsLandauTable {
-    std::vector<float> values;
-    float stepSize;
-  };
-
-  struct StsHitfinderPar {
-    StsAsicPar asic;  // Asic definitions. Currently assumes same parameters for all asics.
-    int nChannels;    // Total number of channels per module. Hitfinder assumes nChannels / 2 channels per side.
-    std::vector<StsModulePar> modules;
-    StsLandauTable landauTable;
-  };
-
-  struct StsHitfinderTimes {
+  struct HitfinderTimes {
     double timeIO          = 0;
     double timeSortDigi    = 0;
     double timeCluster     = 0;
@@ -67,22 +37,14 @@ namespace cbm::algo
    * handles memory transfer for input / output data and conversion to the
    * regular CBM STS types.
    */
-  class StsHitfinderChain {
+  class HitfinderChain : public SubChain {
 
   public:
-    void SetParameters(const StsHitfinderPar& parameters);
-    const StsHitfinderPar& GetParameters() const { return *fPars; }
+    void SetParameters(const sts::HitfinderPars& parameters);
+    const sts::HitfinderPars& GetParameters() const { return *fPars; }
 
     void operator()(gsl::span<const CbmStsDigi>);
 
-    StsHitfinderTimes GetHitfinderTimes() const { return fHitfinderTimes; }
-
-    /**
-       * Returns access to (host copies of) the raw buffers used by the gpu hitfinder.
-       * This is currently the only way to access the produced hits and cluster.
-       */
-    const StsHitfinderHost& GetHitfinderBuffers() const { return fHitfinderHost; }
-
   private:
     // Shorthands to map module-address onto digis in that module
     using DigiMapSide = std::map<int, std::vector<CbmStsDigi>>;
@@ -96,8 +58,8 @@ namespace cbm::algo
        *  TODO: These values might be wildly off. Look for better estimations.
        *  TODO: Should be configurable at runtime.
        */
-    static constexpr float kClustersPerDigiFactor = 0.5f;
-    static constexpr float kHitsPerClustersFactor = 4.f;
+    static constexpr float kClustersPerDigiFactor = 1.f;
+    static constexpr float kHitsPerClustersFactor = 1.5f;
 
     /**
        * Ensure parameters were set. Raises log(fatal) otherwise.
@@ -127,21 +89,24 @@ namespace cbm::algo
     void FlattenDigis(DigiMap& digiMap);
     void FlattenDigisSide(DigiMapSide& digis, bool isFront);
 
-    void CollectTimingInfo();
-
     /**
        * Transfer Hits / Clusters back to host and convert to common CBM types.
        */
-    void AppendClustersModule(int module, bool isFront, std::vector<StsGpuCluster>&);
-    void AppendHitsModule(int module, std::vector<StsGpuHit>&);
+    void AppendClustersModule(int module, bool isFront, std::vector<sts::Cluster>&);
+    void AppendHitsModule(int module, std::vector<sts::Hit>&);
+
+    // Debug functions, ensure reco produces sane results
+    void EnsureDigiOffsets(DigiMap&);
+    void EnsureDigisSorted();
+    void EnsureChannelOffsets(span<u32>);
+    void EnsureClustersSane(span<GpuClusterIdx>, span<int>);
+    void EnsureClustersSorted();
 
-    std::optional<const StsHitfinderPar> fPars;
-    StsHitfinderTimes fHitfinderTimes;
+    std::optional<const sts::HitfinderPars> fPars;
 
-    StsHitfinderHost fHitfinderHost;
-    StsHitfinderGpu fHitfinderGpu;
+    Hitfinder fHitfinder;
   };
 
-}  // namespace cbm::algo
+}  // namespace cbm::algo::sts
 
 #endif  // #ifndef CBM_ALGO_STS_HITFINDER_CHAIN_H

algo/detectors/sts/StsUnpackChain.cxx

+/* Copyright (C) 2022 FIAS Frankfurt Institute for Advanced Studies, Frankfurt / Main
+   SPDX-License-Identifier: GPL-3.0-only
+   Authors: Felix Weiglhofer [committer]*/
+#include "StsUnpackChain.h"
+
+#include "CbmStsDigi.h"
+
+#include <Timeslice.hpp>
+
+#include <fairlogger/Logger.h>
+
+#include <xpu/host.h>
+
+using namespace cbm::algo;
+
+void sts::UnpackChain::Init(StsReadoutConfig config)
+{
+  fConfig = config;
+
+  // Copied from CbmTaskUnpack in CBMRoot.
+
+  // --- Common parameters for all components for STS
+  uint32_t numChansPerAsicSts   = 128;  // R/O channels per ASIC for STS
+  uint32_t numAsicsPerModuleSts = 16;   // Number of ASICs per module for STS
+
+  // Create one algorithm per component for STS and configure it with parameters
+  auto equipIdsSts = fConfig->GetEquipmentIds();
+  for (auto equip : equipIdsSts) {
+    std::unique_ptr<UnpackStsPar> par(new UnpackStsPar());
+    par->fNumChansPerAsic   = numChansPerAsicSts;
+    par->fNumAsicsPerModule = numAsicsPerModuleSts;
+    const size_t numElinks  = fConfig->GetNumElinks(equip);
+    for (size_t elink = 0; elink < numElinks; elink++) {
+      UnpackStsElinkPar elinkPar;
+      auto mapEntry        = fConfig->Map(equip, elink);
+      elinkPar.fAddress    = mapEntry.first;   // Module address for this elink
+      elinkPar.fAsicNr     = mapEntry.second;  // ASIC number within module
+      elinkPar.fTimeOffset = 0.;
+      elinkPar.fAdcOffset  = 1.;
+      elinkPar.fAdcGain    = 1.;
+      // elinkPar.fIsPulser   = mapEntry.pulser;
+      // TODO: Add parameters for time and ADC calibration
+      par->fElinkParams.push_back(elinkPar);
+    }
+    fAlgoSts[equip].SetParams(std::move(par));
+    LOG(debug) << "--- Configured STS equipment " << equip << " with " << numElinks << " elinks";
+  }  //# equipments
+}
+
+std::vector<CbmStsDigi> sts::UnpackChain::Run(const fles::Timeslice& timeslice)
+{
+  assert(fConfig.has_value() && "UnpackChain: Configuration not initialized (call Init first))");
+
+  xpu::scoped_timer t_("STS Unpacker");
+
+  std::vector<CbmStsDigi> digis;
+  size_t numPulsers = 0;
+
+  for (uint64_t comp = 0; comp < timeslice.num_components(); comp++) {
+
+    // --- Component log
+    size_t numDigisInComp = 0;
+    uint64_t numMsInComp  = 0;
+
+    auto systemId = static_cast<fles::SubsystemIdentifier>(timeslice.descriptor(comp, 0).sys_id);
+
+    if (systemId == fles::SubsystemIdentifier::STS) {
+      const uint16_t equipmentId = timeslice.descriptor(comp, 0).eq_id;
+      // std::cout << "STS Unpacker: Component " << comp << ", equipment ID " << equipmentId << std::endl;
+      // for (auto e : fAlgoSts) {
+      //   std::cout << "STS Unpacker: Equipment ID " << e.first << std::endl;
+      // }
+      const auto algoIt = fAlgoSts.find(equipmentId);
+      assert(algoIt != fAlgoSts.end());
+
+      // The current algorithm works for the STS data format version 0x20 used in 2021.
+      // Other versions are not yet supported.
+      // In the future, different data formats will be supported by instantiating different
+      // algorithms depending on the version.
+      assert(timeslice.descriptor(comp, 0).sys_ver == 0x20);
+
+      // --- Microslice loop
+      numMsInComp = timeslice.num_microslices(comp);
+      for (uint64_t mslice = 0; mslice < numMsInComp; mslice++) {
+        const auto msDescriptor = timeslice.descriptor(comp, mslice);
+        const auto msContent    = timeslice.content(comp, mslice);
+
+        xpu::t_add_bytes(msDescriptor.size);
+
+        auto result = (algoIt->second)(msContent, msDescriptor, timeslice.start_time());
+        LOG(debug1) << "STS Unpacker: Component " << comp << ", microslice " << mslice << ", digis "
+                    << result.first.size() << ", errors " << result.second.fNumNonHitOrTsbMessage << " | "
+                    << result.second.fNumErrElinkOutOfRange << " | " << result.second.fNumErrInvalidFirstMessage
+                    << " | " << result.second.fNumErrInvalidMsSize << " | " << result.second.fNumErrTimestampOverflow
+                    << " | ";
+
+        // numPulsers += result.second.fNumPulserHits;
+        numDigisInComp += result.first.size();
+
+        digis.insert(digis.end(), result.first.begin(), result.first.end());
+      }  //# microslice
+    }    // system STS
+
+  }  // component
+
+  LOG(info) << "Timeslice contains " << digis.size() << " STS digis (discarded " << numPulsers << " pulser hits)";
+  return digis;
+}