From 56819b96292a94732629aeffa6156e992a3e4d13 Mon Sep 17 00:00:00 2001
From: Dominik Smith <smith@th.physik.uni-frankfurt.de>
Date: Thu, 18 Apr 2024 12:08:08 +0200
Subject: [PATCH] - Added online capable TRD hitfinders (1D and 2D) in
cbm::algo. - Added ROOT-free implementations of CbmTrdCluster and CbmTrdHit.
- Added class CbmTaskTrdHitFinder as wrapper for new hitfinders. - Added
class CbmTaskTrdHitFinderParWrite to produce YAML input files. - Added
example macro trd_hitfinder_run.C in beamtime/mcbm2022 folder to run new
hitfinder. - Added a required getter function to CbmTrdParModGeo.h. - Updated
CMakeLists and LinkDef files to enable compilation.
---
algo/CMakeLists.txt | 9 +
algo/detectors/trd/Cluster.cxx | 24 +
algo/detectors/trd/Cluster.h | 85 ++
algo/detectors/trd/Cluster2D.cxx | 149 ++++
algo/detectors/trd/Cluster2D.h | 174 ++++
algo/detectors/trd/Clusterizer.cxx | 203 +++++
algo/detectors/trd/Clusterizer.h | 57 ++
algo/detectors/trd/Clusterizer2D.cxx | 119 +++
algo/detectors/trd/Clusterizer2D.h | 52 ++
algo/detectors/trd/DigiRec.cxx | 123 +++
algo/detectors/trd/DigiRec.h | 79 ++
algo/detectors/trd/Hit.cxx | 83 ++
algo/detectors/trd/Hit.h | 172 ++++
algo/detectors/trd/HitFinder.cxx | 233 +++++
algo/detectors/trd/HitFinder.h | 68 ++
algo/detectors/trd/HitFinder2D.cxx | 934 ++++++++++++++++++++
algo/detectors/trd/HitFinder2D.h | 242 +++++
algo/detectors/trd/HitFinder2DPars.h | 33 +
algo/detectors/trd/HitFinderPars.h | 33 +
algo/detectors/trd/Hitfind.cxx | 138 +++
algo/detectors/trd/Hitfind.h | 84 ++
algo/detectors/trd/Hitfind2DSetup.h | 64 ++
algo/detectors/trd/HitfindSetup.h | 65 ++
core/detectors/trd/CbmTrdParModGeo.h | 6 +-
macro/beamtime/mcbm2022/trd_hitfinder_run.C | 216 +++++
reco/tasks/CMakeLists.txt | 2 +
reco/tasks/CbmRecoTasksLinkDef.h | 2 +
reco/tasks/CbmTaskTrdHitFinder.cxx | 153 ++++
reco/tasks/CbmTaskTrdHitFinder.h | 93 ++
reco/tasks/CbmTaskTrdHitFinderParWrite.cxx | 213 +++++
reco/tasks/CbmTaskTrdHitFinderParWrite.h | 63 ++
31 files changed, 3969 insertions(+), 2 deletions(-)
create mode 100644 algo/detectors/trd/Cluster.cxx
create mode 100644 algo/detectors/trd/Cluster.h
create mode 100644 algo/detectors/trd/Cluster2D.cxx
create mode 100644 algo/detectors/trd/Cluster2D.h
create mode 100644 algo/detectors/trd/Clusterizer.cxx
create mode 100644 algo/detectors/trd/Clusterizer.h
create mode 100644 algo/detectors/trd/Clusterizer2D.cxx
create mode 100644 algo/detectors/trd/Clusterizer2D.h
create mode 100644 algo/detectors/trd/DigiRec.cxx
create mode 100644 algo/detectors/trd/DigiRec.h
create mode 100644 algo/detectors/trd/Hit.cxx
create mode 100644 algo/detectors/trd/Hit.h
create mode 100644 algo/detectors/trd/HitFinder.cxx
create mode 100644 algo/detectors/trd/HitFinder.h
create mode 100644 algo/detectors/trd/HitFinder2D.cxx
create mode 100644 algo/detectors/trd/HitFinder2D.h
create mode 100644 algo/detectors/trd/HitFinder2DPars.h
create mode 100644 algo/detectors/trd/HitFinderPars.h
create mode 100644 algo/detectors/trd/Hitfind.cxx
create mode 100644 algo/detectors/trd/Hitfind.h
create mode 100644 algo/detectors/trd/Hitfind2DSetup.h
create mode 100644 algo/detectors/trd/HitfindSetup.h
create mode 100644 macro/beamtime/mcbm2022/trd_hitfinder_run.C
create mode 100644 reco/tasks/CbmTaskTrdHitFinder.cxx
create mode 100644 reco/tasks/CbmTaskTrdHitFinder.h
create mode 100644 reco/tasks/CbmTaskTrdHitFinderParWrite.cxx
create mode 100644 reco/tasks/CbmTaskTrdHitFinderParWrite.h
diff --git a/algo/CMakeLists.txt b/algo/CMakeLists.txt
index 7cbf4a011a..aaadf755c8 100644
--- a/algo/CMakeLists.txt
+++ b/algo/CMakeLists.txt
@@ -117,6 +117,15 @@ set(SRCS
detectors/bmon/ReadoutConfig.cxx
detectors/bmon/Unpack.cxx
detectors/bmon/UnpackMS.cxx
+ detectors/trd/Cluster.cxx
+ detectors/trd/Clusterizer.cxx
+ detectors/trd/Cluster2D.cxx
+ detectors/trd/Clusterizer2D.cxx
+ detectors/trd/DigiRec.cxx
+ detectors/trd/Hit.cxx
+ detectors/trd/Hitfind.cxx
+ detectors/trd/HitFinder.cxx
+ detectors/trd/HitFinder2D.cxx
detectors/trd/ReadoutConfig.cxx
detectors/trd/Unpack.cxx
detectors/trd/UnpackMS.cxx
diff --git a/algo/detectors/trd/Cluster.cxx b/algo/detectors/trd/Cluster.cxx
new file mode 100644
index 0000000000..867cfe1ef9
--- /dev/null
+++ b/algo/detectors/trd/Cluster.cxx
@@ -0,0 +1,24 @@
+/* Copyright (C) 2024 Facility for Antiproton and Ion Research in Europe, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Dominik Smith [committer], Alexandru Bercuci */
+
+#include "Cluster.h"
+
+#include "CbmTrdDigi.h"
+
+namespace cbm::algo::trd
+{
+ //____________________________________________________________________
+ Cluster::Cluster(const std::vector<int32_t>& indices, const std::vector<const CbmTrdDigi*>& digis, int32_t address,
+ uint16_t ncols, uint16_t nrows)
+ : fDigiInd()
+ , fDigis()
+ , fAddress(address)
+ , fNCols(ncols)
+ {
+ fDigiInd.assign(indices.begin(), indices.end());
+ fDigis.assign(digis.begin(), digis.end());
+ SetNRows(nrows);
+ }
+
+} // namespace cbm::algo::trd
diff --git a/algo/detectors/trd/Cluster.h b/algo/detectors/trd/Cluster.h
new file mode 100644
index 0000000000..ea92a66a8f
--- /dev/null
+++ b/algo/detectors/trd/Cluster.h
@@ -0,0 +1,85 @@
+/* Copyright (C) 2024 Facility for Antiproton and Ion Research in Europe, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Dominik Smith [committer], Alexandru Bercuci */
+
+#pragma once
+
+#include <cstddef>
+#include <cstdint>
+#include <vector>
+
+class CbmTrdDigi;
+
+namespace cbm::algo::trd
+{
+ /**
+ * \class Cluster
+ * \author Dominik Smith <d.smith@gsi.de>
+ * \brief Data container for TRD clusters.
+ */
+ class Cluster {
+
+ public:
+ /**
+ * \brief Default constructor.
+ */
+ Cluster() = delete;
+ Cluster(const std::vector<int32_t>& indices, const std::vector<const CbmTrdDigi*>& digis, int32_t address,
+ uint16_t ncols, uint16_t nrows);
+
+ /**
+ * \brief Destructor.
+ */
+ virtual ~Cluster(){};
+
+ /**
+ * \brief Number of digis in cluster.
+ * \return Number of digis in cluster.
+ */
+ size_t GetNofDigis() const { return fDigiInd.size(); }
+
+ /**
+ * \brief Get digi at position index.
+ * \param[in] index Position of digi in array.
+ * \return Digi index in TClonesArray.
+ */
+ int32_t GetDigi(int32_t index) const { return fDigiInd[index]; }
+
+ /**
+ * \brief Get array of digi pointers.
+ * \return Array of digi pointers.
+ */
+ const std::vector<const CbmTrdDigi*>& GetDigis() const { return fDigis; }
+
+ /**
+ * \brief Get array of digi indices.
+ * \return Array of digi indices in TClonesArray.
+ */
+ const std::vector<int32_t>& GetDigiIndices() const { return fDigiInd; }
+
+ /** Accessors **/
+ int32_t GetAddress() const { return fAddress; }
+ uint16_t GetStartCh() const { return fStartCh; }
+ uint16_t GetNCols() const { return fNCols; }
+ uint16_t GetNRows() const { return fNRows & 0x1f; }
+ uint16_t GetNRowsRaw() const { return fNRows; }
+ uint32_t GetStartTime() const { return fStartTime; }
+ bool HasFaspDigis() const { return false; }
+
+ private:
+ std::vector<int32_t> fDigiInd; ///< Array of digi indices
+ std::vector<const CbmTrdDigi*> fDigis; ///< Array of digi pointers
+
+ int32_t fAddress = 0; ///< Unique detector ID
+ uint8_t fNCols = 0; //< number of columns with charge above threshold
+ uint8_t fNRows = 0x1f; //< cluster row info plus extra meta data. Use dedicated getters for the correct value
+ uint16_t fStartCh = 0xffff; //< channel address of first channel
+ uint32_t fStartTime = 0xffffffff; //! start time of cluster in clk units wrt buffer start
+
+ void SetNRows(uint16_t nrows)
+ {
+ fNRows &= (7 << 5);
+ fNRows |= (nrows & 0x1f);
+ }
+ };
+} // namespace cbm::algo::trd
diff --git a/algo/detectors/trd/Cluster2D.cxx b/algo/detectors/trd/Cluster2D.cxx
new file mode 100644
index 0000000000..35253c2bfa
--- /dev/null
+++ b/algo/detectors/trd/Cluster2D.cxx
@@ -0,0 +1,149 @@
+/* Copyright (C) 2024 Facility for Antiproton and Ion Research in Europe, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Dominik Smith [committer], Alexandru Bercuci */
+
+#include "Cluster2D.h"
+
+#include "CbmTrdDigi.h"
+
+#include <cmath>
+
+using std::vector;
+
+namespace cbm::algo::trd
+{
+ //____________________________________________________________________
+ Cluster2D::Cluster2D(const Cluster2D& ref)
+ : fDigis(ref.fDigis)
+ , fDigiIndices(ref.fDigiIndices)
+ , fAddress(ref.fAddress)
+ , fNCols(ref.fNCols)
+ , fNRows(ref.fNRows)
+ , fStartCh(ref.fStartCh)
+ , fStartTime(ref.fStartTime)
+ {
+ }
+
+ //____________________________________________________________________
+ Cluster2D::Cluster2D(int32_t address, int32_t idx, const CbmTrdDigi* digi, uint16_t chT, uint16_t chR, int32_t row,
+ int32_t time)
+ : fAddress(address)
+ , fNCols(0)
+ , fStartCh(0xffff)
+ , fStartTime(time)
+ {
+ SetNRows(row);
+ SetStart(false);
+ SetStop(false);
+ AddDigi(idx, digi, chT, chR);
+ }
+
+ //____________________________________________________________________
+ bool Cluster2D::AddDigi(int32_t idx, const CbmTrdDigi* digi, uint16_t chT, uint16_t chR, int32_t dt)
+ {
+ /** Extend basic functionality of Cluster2D::AddDigi() for the case of 2D.
+ * If chT < 0 use the basic functionality [default].
+ *
+ * For the 2D the parameters are intergpreted as follows
+ * chT : tilted paired channel [default 0x0fffffff]
+ * chR : rectangular paired channel
+ * dt : offset in clks of the prompt signal
+ *
+ * if chT and chR positive the (chT, chR) are interpreted as the 2 channels
+ * of the digi specific to the 2D version. The following specific cases
+ * can be distinguished :
+ * - ch == 0 : no data, cluster signal sequence terminator
+ * - ch == -ch : no data, channel masked in HW
+ */
+
+ if (chT == 0xffff) { // basic functionality for rectangular pads
+ AddDigiIdxPair(idx, digi);
+ return true;
+ }
+
+ uint16_t chMin = (chT != 0 ? chT : chR), chMax = (chR != 0 ? chR : chT);
+
+ // assume triangular pads only
+ if (!fNCols) { // first digi
+ fStartCh = chMin;
+ AddDigiIdxPair(idx, digi);
+ }
+ else if (chMin > GetEndCh()) { // digi @ end
+ //if (HasStop()) return false;
+ AddDigiIdxPair(idx, digi);
+ }
+ else if (chMax < fStartCh) { // digi @ beginning
+ //if (HasStart()) return false;
+ fStartCh = chMin;
+ vector<int32_t> idx_vec = GetDigiIndices();
+ vector<const CbmTrdDigi*> digi_vec = GetDigis();
+ ClearDigis();
+ AddDigiIdxPair(idx, digi);
+ AddDigiIdxPairs(idx_vec, digi_vec);
+ }
+ int nch(0);
+ if (chT == 0)
+ SetStart();
+ else
+ nch++;
+ if (chR == 0)
+ SetStop();
+ else
+ nch++;
+
+ fNCols += nch;
+ if (dt > 0) fStartTime -= dt;
+
+ return true;
+ }
+
+ //____________________________________________________________________
+ int32_t Cluster2D::IsChannelInRange(uint16_t chT, uint16_t chR) const
+ {
+ if (!fNCols) return -2;
+ // if(IsTerminatedLeft() && fAddressCh[0]>ch) return -1;
+ // if(IsTerminatedRight() && fAddressCh[clSize-1]<ch) return 1;
+
+ uint16_t chMin = (chT != 0 ? chT : chR), chMax = (chR != 0 ? chR : chT);
+ if (fStartCh > chMax + 1) return -1;
+ if (fStartCh + fNCols < chMin) return 1;
+ return 0;
+ }
+
+ //____________________________________________________________________
+ bool Cluster2D::Merge(Cluster2D* second)
+ {
+ if (GetRow() != second->GetRow()) return false;
+ // time difference condition
+ if (fNCols == 1 || second->fNCols == 1) {
+ if (abs(int32_t(second->fStartTime - fStartTime)) > 50) return false;
+ }
+ else if (abs(int32_t(second->fStartTime - fStartTime)) > 20)
+ return false;
+ // look before current
+ if (second->fStartCh + second->fNCols == fStartCh) {
+ fStartCh = second->fStartCh;
+ fNCols += second->fNCols;
+ fStartTime = std::min(fStartTime, second->fStartTime);
+
+ vector<int32_t> idx_vec = GetDigiIndices();
+ vector<const CbmTrdDigi*> digi_vec = GetDigis();
+ ClearDigis();
+ AddDigiIdxPairs(second->GetDigiIndices(), second->GetDigis());
+ AddDigiIdxPairs(idx_vec, digi_vec);
+ if (second->HasStart()) SetStart();
+ return true;
+ }
+
+ // look after current
+ if (fStartCh + fNCols == second->fStartCh) {
+ fNCols += second->fNCols;
+ fStartTime = std::min(fStartTime, second->fStartTime);
+ AddDigiIdxPairs(second->GetDigiIndices(), second->GetDigis());
+ if (second->HasStop()) SetStop();
+ return true;
+ }
+ return false;
+ }
+
+} // namespace cbm::algo::trd
diff --git a/algo/detectors/trd/Cluster2D.h b/algo/detectors/trd/Cluster2D.h
new file mode 100644
index 0000000000..b2b362e350
--- /dev/null
+++ b/algo/detectors/trd/Cluster2D.h
@@ -0,0 +1,174 @@
+/* Copyright (C) 2024 Facility for Antiproton and Ion Research in Europe, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Dominik Smith [committer], Alexandru Bercuci */
+
+#pragma once
+
+#include <cstdint>
+#include <vector> // for vector
+
+#define BIT(n) (1ULL << (n))
+#define SETBIT(n, i) ((n) |= BIT(i))
+#define CLRBIT(n, i) ((n) &= ~BIT(i))
+#define TESTBIT(n, i) ((bool) (((n) &BIT(i)) != 0))
+
+class CbmTrdDigi;
+
+namespace cbm::algo::trd
+{
+ /**
+ * \class Cluster2D
+ * \author Dominik Smith <d.smith@gsi.de>
+ * \brief Data Container for TRD clusters.
+ */
+ class Cluster2D {
+
+ public:
+ enum eClusterDef
+ {
+ kFasp = 5 ///< set type of FEE digis contained
+ ,
+ kStart ///< only for triangular if no T in first col
+ ,
+ kStop ///< only for triangular if no R in last col
+ };
+
+ /**
+ * \brief Default constructor.
+ */
+ Cluster2D() = delete;
+ Cluster2D(const Cluster2D& ref);
+ /**
+ * \brief Constructor starting from first digit (FASP specific).
+ * \param[in] address global module address
+ * \param[in] idx global digi index in the TClonesArray
+ * \param[in] digi pointer to digi
+ * \param[in] chT RO channel address within the module for tilt pairing
+ * \param[in] chR RO channel address within the module for rect pairing
+ * \param[in] r module row for the RO channel
+ * \param[in] time relative buffer DAQ time
+ */
+ Cluster2D(int32_t address, int32_t idx, const CbmTrdDigi* digi, uint16_t chT, uint16_t chR, int32_t r,
+ int32_t time);
+ /**
+ * \brief Destructor.
+ */
+ virtual ~Cluster2D(){};
+
+ /**
+ * \brief Add digi to cluster.
+ * \param[in] index Digi index in TClonesArray.
+ * \param[in] digi pointer to digi
+ */
+ void AddDigiIdxPair(int32_t index, const CbmTrdDigi* digi)
+ {
+ fDigiIndices.push_back(index);
+ fDigis.push_back(digi);
+ }
+
+ /**
+ * \brief Add array of digi to cluster.
+ * \param[in] indices Array of digi indices in TClonesArray.
+ */
+ void AddDigiIdxPairs(const std::vector<int32_t>& indices, const std::vector<const CbmTrdDigi*> digis)
+ {
+ fDigiIndices.insert(fDigiIndices.end(), indices.begin(), indices.end());
+ fDigis.insert(fDigis.end(), digis.begin(), digis.end());
+ }
+
+ /**
+ * \brief Number of digis in cluster.
+ * \return Number of digis in cluster.
+ */
+ int32_t GetNofDigis() const { return fDigiIndices.size(); }
+
+ /**
+ * \brief Get digi at position index.
+ * \param[in] index Position of digi in array.
+ * \return Digi index in TClonesArray.
+ */
+ int32_t GetDigi(int32_t index) const { return fDigiIndices[index]; }
+
+ /**
+ * \brief Get array of digi indices.
+ * \return Array of digi indices in TClonesArray.
+ */
+ const std::vector<int32_t>& GetDigiIndices() const { return fDigiIndices; }
+
+ /**
+ * \brief Get array of digi pointers.
+ * \return Array of digi pointers
+ */
+ std::vector<const CbmTrdDigi*>& GetDigis() { return fDigis; }
+
+
+ /**
+ * \brief Remove all digis.
+ */
+ void ClearDigis()
+ {
+ fDigiIndices.clear();
+ fDigis.clear();
+ }
+
+ /** Accessors **/
+ int32_t GetAddress() const { return fAddress; }
+
+ /** \brief Append digi to cluster
+ * \param[in] idx index of digi in TClonesArray
+ * \param[in] digi pointer to digi
+ * \param[in] chT RO channel for digi (tilt pairing for FASP) default 0xffff (SPADIC)
+ * \param[in] chR RO channel for rect pairing (only for FASP)
+ * \param[in] dt update start time of cluster if current digi is prompt
+ * \return true if successful
+ */
+ bool AddDigi(int32_t idx, const CbmTrdDigi* digi, uint16_t chT = 0xffff, uint16_t chR = 0, int32_t dt = 0);
+
+ /** Accessors **/
+ uint16_t GetNCols() const { return fNCols; }
+ uint16_t GetNRows() const { return fNRows & 0x1f; }
+ uint16_t GetNRowsRaw() const { return fNRows; }
+ uint16_t GetEndCh() const { return fStartCh + fNCols - 1; }
+ uint16_t GetRow() const { return GetNRows(); }
+ uint16_t GetSize() const { return GetNCols(); }
+ uint16_t GetStartCh() const { return fStartCh; }
+ uint32_t GetStartTime() const { return fStartTime; }
+ bool HasFaspDigis() const { return true; }
+ bool HasStart() const { return TESTBIT(fNRows, kStart); }
+ bool HasStop() const { return TESTBIT(fNRows, kStop); }
+
+ /** \brief Query on RO channels list
+ * \param[in] chT tilted RO channel for digi
+ * \param[in] chR rectangular RO channel for digi
+ * \return -1 before range; 0 in range; 1 after range; -2 cluster empty of digits
+ */
+ int32_t IsChannelInRange(uint16_t chT, uint16_t chR) const;
+
+ /** \brief Merge current cluster with info from second
+ * \param[in] second cluster to be added
+ * \param[in] typ the type of pad-plane of the source chamber; true if Trd2d
+ * \return success or fail
+ */
+ bool Merge(Cluster2D* second);
+
+ /** Setters **/
+ void SetNCols(uint16_t ncols) { fNCols = ncols; }
+ void SetNRows(uint16_t nrows)
+ {
+ fNRows &= (7 << 5);
+ fNRows |= (nrows & 0x1f);
+ }
+ void SetStart(bool set = true) { set ? SETBIT(fNRows, kStart) : CLRBIT(fNRows, kStart); }
+ void SetStop(bool set = true) { set ? SETBIT(fNRows, kStop) : CLRBIT(fNRows, kStop); }
+
+ private:
+ std::vector<const CbmTrdDigi*> fDigis; ///< Array of digi pointers
+ std::vector<int32_t> fDigiIndices; ///< Array of digi indices
+
+ int32_t fAddress = 0; ///< Unique detector ID
+ uint8_t fNCols = 0; //< number of columns with charge above threshold
+ uint8_t fNRows = 0x1f; //< cluster row info plus extra meta data. Use dedicated getters for the correct value
+ uint16_t fStartCh = 0xffff; //< channel address of first channel
+ uint32_t fStartTime = 0xffffffff; //! start time of cluster in clk units wrt buffer start
+ };
+} // namespace cbm::algo::trd
diff --git a/algo/detectors/trd/Clusterizer.cxx b/algo/detectors/trd/Clusterizer.cxx
new file mode 100644
index 0000000000..3332d6d221
--- /dev/null
+++ b/algo/detectors/trd/Clusterizer.cxx
@@ -0,0 +1,203 @@
+/* Copyright (C) 2024 Facility for Antiproton and Ion Research in Europe, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Dominik Smith [committer], Etienne Bechtel, Florian Uhlig */
+
+#include "Clusterizer.h"
+
+#include <algorithm>
+#include <unordered_map>
+
+namespace cbm::algo::trd
+{
+
+ //_______________________________________________________________________________
+ std::vector<Cluster> Clusterizer::operator()(const std::vector<std::pair<CbmTrdDigi, int32_t>>& inVec) const
+ {
+ std::vector<inputType> inputData; //digis storage;
+ std::vector<Cluster> clustersOut;
+ std::vector<inputType*> digiPtrVec; //digis pointer storage (kSelf only);
+ std::vector<std::vector<inputType*>> digiMapSelf; //channel-wise sorted input
+ std::vector<std::vector<inputType*>> digiMapNeig; //channel-wise sorted input
+
+ const size_t numCols = fParams.rowPar[0].padPar.size();
+ const size_t numRows = fParams.rowPar.size();
+ const size_t numChan = numCols * numRows;
+
+ for (auto& input : inVec) {
+ inputData.push_back(std::make_tuple(input.second, &input.first, input.first.GetTime()));
+ }
+ ///// DRAW PICTURE HERE (LOGIC OF MAPS)
+
+ // Initialize the digi maps
+ digiMapSelf.resize(numChan);
+ digiMapNeig.resize(numChan);
+
+ for (auto& digidata : inputData) {
+ const CbmTrdDigi* digi = (const CbmTrdDigi*) std::get<1>(digidata);
+ if (!digi) continue; // Skip invalid digis
+ const CbmTrdDigi::eTriggerType type = static_cast<CbmTrdDigi::eTriggerType>(digi->GetTriggerType());
+ const int channel = digi->GetAddressChannel();
+ if (type == CbmTrdDigi::eTriggerType::kSelf) {
+ digiMapSelf[channel].push_back(&digidata);
+ digiPtrVec.push_back(&digidata);
+ }
+ if (type == CbmTrdDigi::eTriggerType::kNeighbor) {
+ digiMapNeig[channel].push_back(&digidata);
+ }
+ }
+
+ //Store last position in input channels to avoid unnecessary checks.
+ std::vector<std::vector<inputType*>::iterator> lastChanPosSelf, lastChanPosNeig;
+ for (size_t chan = 0; chan < numChan; chan++) {
+ lastChanPosSelf.push_back(digiMapSelf[chan].begin());
+ lastChanPosNeig.push_back(digiMapNeig[chan].begin());
+ }
+
+ // search for an unprocessed main triggered digi and then start a subloop to
+ // directly construct the cluster (search for main-trigger then add the neighbors)
+ for (auto mainit = digiPtrVec.begin(); mainit != digiPtrVec.end(); mainit++) {
+
+ // Skip if digi already processed
+ const CbmTrdDigi* digi = (const CbmTrdDigi*) std::get<1>(**mainit);
+ if (digi == nullptr) continue;
+
+ // get digi time
+ const double time = std::get<2>(**mainit);
+
+ // variety of neccessary address information; uses the "combiId" for the
+ // comparison of digi positions
+ const int digiId = std::get<0>(**mainit);
+ const int channel = digi->GetAddressChannel();
+ const int ncols = fParams.rowPar[0].padPar.size();
+
+ // some logic information which is used to process and find the clusters
+ int lowcol = channel;
+ int highcol = channel;
+
+ // the "seal" bools register when the logical end of the cluster was found
+ bool sealtopcol = false;
+ bool sealbotcol = false;
+
+ // //vector which contains the actual cluster
+ std::vector<std::pair<int, const CbmTrdDigi*>> cluster;
+ cluster.emplace_back(digiId, digi);
+ std::get<1>(**mainit) = nullptr;
+
+ // loop to find the other pads corresponding to the main trigger
+ // is exited either if the implemented trigger logic is fullfilled
+ // or if there are no more adjacend pads due to edges,etc.
+ while (true) {
+
+ // counter which is used to easily break clusters which are at the edge and
+ // therefore do not fullfill the classical look
+ const size_t oldSize = cluster.size();
+
+ const int col = channel % ncols;
+ if (col == ncols) sealtopcol = true; ///// Axel comment: change to ncols - 1, also for main!!
+ if (col == 0) sealbotcol = true;
+
+ if (!sealbotcol && 0 <= lowcol - 1) {
+ if (TryAddDigi(&digiMapSelf, lowcol - 1, &lastChanPosSelf, &cluster, time)) {
+ lowcol--;
+ }
+ }
+ if (!sealtopcol && highcol + 1 < numChan) {
+ if (TryAddDigi(&digiMapSelf, highcol + 1, &lastChanPosSelf, &cluster, time)) {
+ highcol++;
+ }
+ }
+ if (!sealbotcol && 0 <= lowcol - 1) {
+ if (TryAddDigi(&digiMapNeig, lowcol - 1, &lastChanPosNeig, &cluster, time)) {
+ sealbotcol = true;
+ }
+ }
+ if (!sealtopcol && highcol + 1 < numChan) {
+ if (TryAddDigi(&digiMapNeig, highcol + 1, &lastChanPosNeig, &cluster, time)) {
+ sealtopcol = true;
+ }
+ }
+
+ ////// Row merging was formerly here. To do: Implement as post-processing
+
+ // some finish criteria
+ if (cluster.size() - oldSize == 0) break;
+ if (sealbotcol && sealtopcol) break;
+ } //! while (true)
+
+ addClusters(cluster, &clustersOut);
+ } //! for (auto mainit = inputData.begin(); mainit != inputData.end(); mainit++)
+
+ /// Row merging could be here in principle.
+
+ return clustersOut;
+ }
+
+ bool Clusterizer::TryAddDigi(std::vector<std::vector<inputType*>>* digimap, int chan,
+ std::vector<std::vector<inputType*>::iterator>* lastPos,
+ std::vector<std::pair<int, const CbmTrdDigi*>>* cluster, const double digiTime) const
+ {
+ const double interval = CbmTrdDigi::Clk(CbmTrdDigi::eCbmTrdAsicType::kSPADIC);
+
+ // find the FN digis of main trigger or adjacent main triggers
+ for (auto FNit = (*lastPos)[chan]; FNit != (*digimap)[chan].end(); FNit++) {
+
+ // update starting time
+ const double newtime = std::get<2>(**FNit);
+ if (newtime < digiTime - interval) {
+ (*lastPos)[chan]++;
+ continue;
+ }
+
+ // break if out of range
+ if (newtime > digiTime + interval) break;
+
+ // Skip already processed digis
+ const CbmTrdDigi* d = (const CbmTrdDigi*) std::get<1>(**FNit);
+ if (d == nullptr) continue;
+
+ // add digi to cluster
+ const int digiid = std::get<0>(**FNit);
+ cluster->emplace_back(digiid, d);
+ std::get<1>(**FNit) = nullptr;
+ return true;
+ }
+ return false;
+ }
+
+ //_____________________________________________________________________
+ void Clusterizer::addClusters(std::vector<std::pair<int, const CbmTrdDigi*>> cluster,
+ std::vector<Cluster>* clustersOut) const
+ {
+ // create vector for indice matching
+ std::vector<int> digiIndices;
+ digiIndices.reserve(cluster.size());
+
+ // add digi ids to vector
+ std::transform(cluster.begin(), cluster.end(), std::back_inserter(digiIndices),
+ [](const auto& pair) { return pair.first; });
+
+ // create vector for pointers
+ std::vector<const CbmTrdDigi*> digis;
+ digis.reserve(cluster.size());
+
+ // add digi ids to vector
+ std::transform(cluster.begin(), cluster.end(), std::back_inserter(digis),
+ [](const auto& pair) { return pair.second; });
+
+ // Count rows and columns (TO DO: needs to be re-computed / moved if row merging is implemented)
+ std::unordered_map<int, bool> cols, rows;
+ for (const auto& pair : cluster) {
+ const int ncols = fParams.rowPar[0].padPar.size();
+ int digiChannel = pair.second->GetAddressChannel();
+ int colId = digiChannel % ncols;
+ int globalRow = digiChannel / ncols;
+ int combiId = globalRow * ncols + colId;
+ cols[combiId] = true;
+ rows[globalRow] = true;
+ }
+
+ // add the cluster to the Array
+ clustersOut->emplace_back(digiIndices, digis, fParams.address, cols.size(), rows.size());
+ }
+
+} // namespace cbm::algo::trd
diff --git a/algo/detectors/trd/Clusterizer.h b/algo/detectors/trd/Clusterizer.h
new file mode 100644
index 0000000000..35a203b740
--- /dev/null
+++ b/algo/detectors/trd/Clusterizer.h
@@ -0,0 +1,57 @@
+/* Copyright (C) 2024 Facility for Antiproton and Ion Research in Europe, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Dominik Smith [committer], Etienne Bechtel, Florian Uhlig */
+
+#pragma once
+
+#include "CbmTrdDigi.h"
+#include "Cluster.h"
+#include "HitFinderPars.h"
+
+#include <tuple>
+#include <vector>
+
+namespace cbm::algo::trd
+{
+
+ //// TO DO: Add monitoring data here!!
+
+ /** @class Clusterizer
+ ** @brief Algo class for TRD cluster building
+ ** @author Dominik Smith <d.smith@gsi.de>
+ ** @since 09.04.2024
+ **
+ **/
+
+ class Clusterizer {
+ public:
+ typedef std::tuple<int, const CbmTrdDigi*, double>
+ inputType; //digi index, pointer to digi (null if processed), digi time
+
+ /** \brief Default constructor.*/
+ Clusterizer() = default;
+
+ /** \brief Default constructor.*/
+ Clusterizer(HitFinderModPar par) : fParams(par){};
+
+ /** @brief Destructor **/
+ virtual ~Clusterizer(){};
+
+ /** @brief Execution
+ ** @param inVec Digi data for one module
+ ** @return Vector of constructed clusters
+ **/
+ std::vector<Cluster> operator()(const std::vector<std::pair<CbmTrdDigi, int32_t>>& inVec) const;
+
+ protected:
+ private:
+ HitFinderModPar fParams; ///< Parameter container
+
+ bool TryAddDigi(std::vector<std::vector<inputType*>>* digimap, int chan,
+ std::vector<std::vector<inputType*>::iterator>* lastPos,
+ std::vector<std::pair<int, const CbmTrdDigi*>>* cluster, const double digiTime) const;
+
+ void addClusters(std::vector<std::pair<int, const CbmTrdDigi*>> cluster, std::vector<Cluster>* clustersOut) const;
+ };
+
+} // namespace cbm::algo::trd
diff --git a/algo/detectors/trd/Clusterizer2D.cxx b/algo/detectors/trd/Clusterizer2D.cxx
new file mode 100644
index 0000000000..3782ee8681
--- /dev/null
+++ b/algo/detectors/trd/Clusterizer2D.cxx
@@ -0,0 +1,119 @@
+/* Copyright (C) 2024 Facility for Antiproton and Ion Research in Europe, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Dominik Smith [committer], Alexandru Bercuci */
+
+#include "Clusterizer2D.h"
+
+#include <algorithm>
+
+namespace cbm::algo::trd
+{
+
+ //_______________________________________________________________________________
+ std::vector<Cluster2D> Clusterizer2D::operator()(const std::vector<std::pair<CbmTrdDigi, int32_t>>& inVec,
+ const uint64_t t0) const
+ {
+ const size_t numRows = fParams.rowPar.size();
+ const size_t numCols = fParams.rowPar[0].padPar.size();
+
+ std::vector<std::vector<Cluster2D*>> buffer(numRows); //row-wise organized clusters
+ std::vector<size_t> bufferPos(numRows); //storage of starting positon in buffers
+ std::vector<Cluster2D> clustersOut; //clusters to be returned
+ std::vector<inputType> inputData; //formatted input
+
+ // Apply input formatting
+ for (auto& input : inVec) {
+ const CbmTrdDigi* d = &input.first;
+ const size_t index = input.second;
+ int pad = d->GetAddressChannel();
+ uint16_t chT = pad << 1;
+ uint16_t chR = chT + 1;
+ int dtime;
+ double t;
+ double r = d->GetCharge(t, dtime); //subtract timeslice start time. maybe not needed anymore?
+ int tm = d->GetTimeDAQ() - t0;
+ const int row = pad / numCols;
+ const int col = pad % numCols;
+ if (dtime < 0) tm += dtime; // correct for the difference between tilt and rect
+ if (r < 1 && !fParams.rowPar[row].padPar[col].chRMasked) chR = 0;
+ if (t < 1 && !fParams.rowPar[row].padPar[col].chTMasked) chT = 0;
+ inputData.emplace_back(chT, chR, tm, row, index, d);
+ }
+
+ // Sort input by time. TO DO: Use insert-sort like in TOF hitfinder
+ std::sort(inputData.begin(), inputData.end(),
+ [](const auto& a, const auto& b) -> bool { return std::get<2>(a) < std::get<2>(b); });
+
+ // Fill row buffers
+ for (auto& in : inputData) {
+ auto& [chT, chR, tm, row, id, digi] = in;
+
+ // check for saved
+ if (buffer[row].empty()) {
+ buffer[row].push_back(new Cluster2D(fParams.address, id, digi, chT, chR, row, tm));
+ bufferPos[row] = 0;
+ continue;
+ }
+ /////// To do: Check whether position storing truly helps, after Cluster2D::Merge() is optimized
+
+ // Get stored starting position
+ auto itc = (buffer[row].begin() + bufferPos[row]);
+
+ //Update cluster
+ for (; itc != buffer[row].end(); itc++) {
+ if ((int64_t)(*itc)->GetStartTime() < tm - 5) {
+ bufferPos[row]++;
+ continue;
+ }
+ if ((*itc)->IsChannelInRange(chT, chR) == 0) {
+ const uint tc = (*itc)->GetStartTime();
+ (*itc)->AddDigi(id, digi, chT, chR, tc - tm);
+ break;
+ }
+ }
+ if (itc == buffer[row].end()) {
+ buffer[row].push_back(new Cluster2D(fParams.address, id, digi, chT, chR, row, tm));
+ }
+ }
+
+ // Build clusters
+ for (auto& clRow : buffer) {
+
+ // TO DO: Process the following remark from original code:
+ // TODO look left and right for masked channels. If they exists add them to cluster.
+ // if (AddClusterEdges(*itc0) && CWRITE(0)) cout << " edge cl[0] : " << (*itc0)->ToString();
+
+ // Phase 0 : try merge clusters if more than one on a row
+ for (auto itc0 = clRow.begin(); itc0 != clRow.end(); itc0++) {
+ if (*itc0 == nullptr) {
+ continue;
+ }
+ for (auto itc1 = std::next(itc0); itc1 != clRow.end(); itc1++) {
+ if (*itc1 == nullptr) {
+ continue;
+ }
+ if ((*itc1)->GetStartTime() - (*itc0)->GetStartTime() > 50) {
+ break;
+ }
+ if ((*itc0)->Merge((*itc1))) {
+ delete (*itc1);
+ *itc1 = nullptr;
+ itc1 = itc0;
+ }
+ }
+ }
+
+ // Phase 1 : copy older clusters from the buffer to the module wise storage
+ for (auto& pcluster : clRow) {
+ if (pcluster == nullptr) {
+ continue;
+ }
+ clustersOut.emplace_back(*pcluster);
+ delete pcluster;
+ pcluster = nullptr;
+ }
+ }
+ return clustersOut;
+ }
+
+} // namespace cbm::algo::trd
diff --git a/algo/detectors/trd/Clusterizer2D.h b/algo/detectors/trd/Clusterizer2D.h
new file mode 100644
index 0000000000..1fcf3134ed
--- /dev/null
+++ b/algo/detectors/trd/Clusterizer2D.h
@@ -0,0 +1,52 @@
+/* Copyright (C) 2024 Facility for Antiproton and Ion Research in Europe, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Dominik Smith [committer], Alexandru Bercuci */
+
+#pragma once
+
+#include "CbmTrdDigi.h"
+#include "Cluster2D.h"
+#include "HitFinder2DPars.h"
+
+#include <tuple>
+#include <vector>
+
+namespace cbm::algo::trd
+{
+
+ //// TO DO: Add monitoring data here!!
+
+ /** @class Clusterizer2D
+ ** @brief Algo class for TRD2D cluster building
+ ** @author Dominik Smith <d.smith@gsi.de>
+ ** @since 05.04.2024
+ **
+ **/
+
+ class Clusterizer2D {
+ public:
+ typedef std::tuple<uint16_t, uint16_t, int, int, size_t, const CbmTrdDigi*>
+ inputType; //Tuple: chT, chR, tm, row, id, digi
+
+ /** \brief Default constructor.*/
+ Clusterizer2D() = default;
+
+ /** \brief Default constructor.*/
+ Clusterizer2D(HitFinder2DModPar par) : fParams(par){};
+
+ /** @brief Destructor **/
+ virtual ~Clusterizer2D(){};
+
+ /** @brief Execution
+ ** @param inVec Digi data for one module
+ ** @param t0 Start time of timeslice
+ ** @return Vector of constructed clusters
+ **/
+ std::vector<Cluster2D> operator()(const std::vector<std::pair<CbmTrdDigi, int32_t>>& inVec, uint64_t t0) const;
+
+ protected:
+ private:
+ HitFinder2DModPar fParams; ///< Parameter container
+ };
+
+} // namespace cbm::algo::trd
diff --git a/algo/detectors/trd/DigiRec.cxx b/algo/detectors/trd/DigiRec.cxx
new file mode 100644
index 0000000000..f74ccbce31
--- /dev/null
+++ b/algo/detectors/trd/DigiRec.cxx
@@ -0,0 +1,123 @@
+/* Copyright (C) 2018-2020 Horia Hulubei National Institute of Physics and Nuclear Engineering, Bucharest
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Alexandru Bercuci[committer] */
+
+#include "DigiRec.h"
+
+#include <cstring>
+
+namespace cbm::algo::trd
+{
+
+ //// Taken from CbmTrdFASP. TO DO: Remove these magic numbers!!
+ float DigiRec::fgBaseline = 0.25; //[V] AB : match HW on mCBM22
+ float DigiRec::fgOutGain = 2.025; // [V/ 4095 ADC chs]
+
+ //_____________________________________________________________________________
+ DigiRec::DigiRec() : CbmTrdDigi(), fStatus(0)
+ {
+ fG[0] = 1.;
+ fG[1] = 1.;
+ memset(fT, 0, 3 * sizeof(double));
+ }
+
+ //_____________________________________________________________________________
+ DigiRec::DigiRec(const CbmTrdDigi& d, double* G, double* T) : CbmTrdDigi(d), fStatus(0)
+ {
+ if (G) {
+ fG[0] = G[0];
+ fG[1] = G[1];
+ }
+ else {
+ fG[0] = 1.;
+ fG[1] = 1.;
+ }
+ if (T)
+ memcpy(fT, T, 3 * sizeof(double));
+ else
+ memset(fT, 0, 3 * sizeof(double));
+
+ int dt;
+ double t, r = CbmTrdDigi::GetCharge(t, dt);
+ if (t > 4094) SETBIT(fStatus, 0);
+ if (r > 4094) SETBIT(fStatus, 1);
+ }
+
+ //_____________________________________________________________________________
+ DigiRec::DigiRec(const CbmTrdDigi& d, const CbmTrdDigi& dr, double* G, double* T) : CbmTrdDigi(d), fStatus(0)
+ {
+ /** Constructor and RAW digi merger
+*/
+ if (G) {
+ fG[0] = G[0];
+ fG[1] = G[1];
+ }
+ else {
+ fG[0] = 1.;
+ fG[1] = 1.;
+ }
+ if (T)
+ memcpy(fT, T, 3 * sizeof(double));
+ else
+ memset(fT, 0, 3 * sizeof(double));
+
+ int dt;
+ double t, r = dr.GetCharge(t, dt);
+ if (r > 4094) SETBIT(fStatus, 1);
+ CbmTrdDigi::GetCharge(t, dt);
+ if (t > 4094) SETBIT(fStatus, 0);
+ dt = dr.GetTimeDAQ() - GetTimeDAQ();
+ SetCharge(t, r, dt);
+ int rtrg(dr.GetTriggerType() & 2), ttrg(GetTriggerType() & 1);
+ SetTriggerType(rtrg | ttrg); //merge the triggers
+ }
+
+ //_____________________________________________________________________________
+ double DigiRec::GetCharge(int typ, bool& on) const
+ {
+ int dT;
+ double T, R = CbmTrdDigi::GetCharge(T, dT);
+ on = true;
+ if (typ) { // rectangular pairing
+ if (R > 0) {
+ R -= GetBaselineCorr();
+ return fG[1] * R;
+ }
+ else
+ on = false;
+ }
+ else { // tilt pairing
+ if (T > 0) {
+ T -= GetBaselineCorr();
+ return fG[0] * T;
+ }
+ else
+ on = false;
+ }
+ return 0.;
+ }
+
+ //_____________________________________________________________________________
+ double DigiRec::GetTime(int typ) const
+ {
+ int dT;
+ double T, R = CbmTrdDigi::GetCharge(T, dT);
+ if (typ) { // rectangular pairing
+ double R0 = R - fT[2];
+ return GetTimeDAQ() + dT + fT[0] * R0 * R0 + fT[1] * R0;
+ }
+ else { // tilt pairing
+ double T0 = T - fT[2];
+ return GetTimeDAQ() + fT[0] * T0 * T0 + fT[1] * T0;
+ }
+ }
+
+ //_____________________________________________________________________________
+ void DigiRec::Init(double g[2], double t[3])
+ {
+ fG[0] = g[0];
+ fG[1] = g[1];
+ memcpy(fT, t, 3 * sizeof(double));
+ }
+
+} // namespace cbm::algo::trd
diff --git a/algo/detectors/trd/DigiRec.h b/algo/detectors/trd/DigiRec.h
new file mode 100644
index 0000000000..86e7136920
--- /dev/null
+++ b/algo/detectors/trd/DigiRec.h
@@ -0,0 +1,79 @@
+/* Copyright (C) 2018-2020 Horia Hulubei National Institute of Physics and Nuclear Engineering, Bucharest
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Alexandru Bercuci[committer] */
+
+#pragma once
+
+#include "CbmTrdDigi.h"
+
+#define BIT(n) (1ULL << (n))
+#define SETBIT(n, i) ((n) |= BIT(i))
+#define CLRBIT(n, i) ((n) &= ~BIT(i))
+#define TESTBIT(n, i) ((bool) (((n) &BIT(i)) != 0))
+
+/** @class DigiRec
+ ** @brief Extend the TRD(2D) digi class to incorporate FEE calibration.
+ ** @author Alexandru Bercucic <abercuci@niham.nipne.ro>
+ ** @since 01.10.2021
+ ** @date 01.10.2021
+ **
+ ** The digi class contains the information as it is produced by the FEE (ASIC/GETS)
+ ** The variation from channel to channel is captured by running the pulser on anode wires
+ ** using various signal values, frequencies, etc. The calibrated baselines, gains, jitter,
+ ** etc. are transported via the parameter files and are applied to the data within the digRec
+ ** class which is in the end used to calculate the TRD hit parameters.
+ **/
+
+namespace cbm::algo::trd
+{
+
+ class DigiRec : public CbmTrdDigi {
+ friend class HitFinder2D;
+
+ public:
+ /** \brief Default constructor*/
+ DigiRec();
+ /** \brief Wrap CbmTrdDigi constructor*/
+ DigiRec(const CbmTrdDigi& d, double* g = NULL, double* t = NULL);
+ virtual ~DigiRec() { ; }
+
+ /** \brief Return calibrated tilt signal
+ * \param[out] on flag signal exists */
+ double GetTiltCharge(bool& on) const { return GetCharge(0, on); }
+ /** \brief Return calibrated tilt time [ns]*/
+ double GetTiltTime() const { return GetTime(0); }
+ /** \brief Return calibrated rect signal
+ * \param[out] on flag signal exists */
+ double GetRectCharge(bool& on) const { return GetCharge(1, on); }
+ /** \brief Return calibrated rect time [ns]*/
+ double GetRectTime() const { return GetTime(1); }
+ /** \brief Return calibrated signal
+ * \param[in] typ tilt [0], rect [1]
+ * \param[out] on flag signal exists
+ */
+ double GetCharge(int typ, bool& on) const;
+ /** \brief Return calibrated time
+ * \param[in] typ tilt [0], rect [1]
+ */
+ double GetTime(int typ) const;
+ bool HasRectOvf() const { return TESTBIT(fStatus, 1); }
+ bool HasTiltOvf() const { return TESTBIT(fStatus, 0); }
+ /** \brief Init FEE gain and time walk corrections */
+ void Init(double g[2], double t[3]);
+
+ static float GetBaselineCorr() { return 4095. * fgBaseline / fgOutGain; }
+
+ protected:
+ /** \brief Constructor and merger*/
+ DigiRec(const CbmTrdDigi& dt, const CbmTrdDigi& dr, double* g = NULL, double* t = NULL);
+
+ private:
+ unsigned char fStatus; //< bit map to store calibration flags
+ double fG[2]; //< FEE gain correction for channel T & R
+ double fT[3]; //< FEE time walk correction as function of charge
+
+ static float fgBaseline; ///< FASP baseline [V]
+ static float fgOutGain; ///< FASP -> ADC gain [V/4095 ADC]
+ };
+
+} // namespace cbm::algo::trd
diff --git a/algo/detectors/trd/Hit.cxx b/algo/detectors/trd/Hit.cxx
new file mode 100644
index 0000000000..363bf40f8b
--- /dev/null
+++ b/algo/detectors/trd/Hit.cxx
@@ -0,0 +1,83 @@
+/* Copyright (C) 2024 Facility for Antiproton and Ion Research in Europe, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Dominik Smith [committer], Matus Kalisky, Florian Uhlig, Andrey Lebedev */
+
+/**
+ * \file Hit.cxx
+ * \author Dominik Smith <d.smith@gsi.de>
+ * \date 2024
+ **/
+#include "Hit.h"
+
+namespace cbm::algo::trd
+{
+
+ Hit::Hit() : Hit(-1, 0., 0., 0., 0., 0., 0., 0., -1, -1., -1.) {}
+
+ Hit::Hit(int32_t address, const ROOT::Math::XYZVector& pos, const ROOT::Math::XYZVector& dpos, double dxy,
+ int32_t refId, double eLoss, double time, double timeError)
+ : Hit(address, pos.X(), pos.Y(), pos.Z(), dpos.X(), dpos.Y(), dpos.Z(), dxy, refId, time, timeError)
+ {
+ fELoss = eLoss;
+ }
+
+ Hit::Hit(int32_t address, double x, double y, double z, double dx, double dy, double dz, double dxy, int32_t refId,
+ double time, double timeError)
+ : fX(x)
+ , fY(y)
+ , fZ(z)
+ , fDx(dx)
+ , fDy(dy)
+ , fDz(dz)
+ , fDxy(dxy)
+ , fRefId(refId)
+ , fAddress(address)
+ , fTime(time)
+ , fTimeError(timeError)
+ , fDefine(0)
+ , fNeighborId(-1)
+ , fELoss(-1.)
+ {
+ }
+
+ Hit& Hit::operator=(const Hit& rhs)
+ {
+ if (this != &rhs) {
+ fX = rhs.fX;
+ fY = rhs.fY;
+ fZ = rhs.fZ;
+ fDx = rhs.fDx;
+ fDy = rhs.fDy;
+ fDz = rhs.fDz;
+ fDxy = rhs.fDxy;
+ fRefId = rhs.fRefId;
+ fAddress = rhs.fAddress;
+ fTime = rhs.fTime;
+ fTimeError = rhs.fTimeError;
+ fDefine = rhs.fDefine;
+ fNeighborId = rhs.fNeighborId;
+ fELoss = rhs.fELoss;
+ }
+ return *this;
+ }
+
+ void Hit::Position(ROOT::Math::XYZVector& pos) const { pos.SetXYZ(GetX(), GetY(), GetZ()); }
+
+ void Hit::PositionError(ROOT::Math::XYZVector& dpos) const { dpos.SetXYZ(GetDx(), GetDy(), GetDz()); }
+
+ void Hit::SetPosition(const ROOT::Math::XYZVector& pos)
+ {
+ SetX(pos.X());
+ SetY(pos.Y());
+ SetZ(pos.Z());
+ }
+
+ void Hit::SetPositionError(const ROOT::Math::XYZVector& dpos)
+ {
+ SetDx(dpos.X());
+ SetDy(dpos.Y());
+ SetDz(dpos.Z());
+ }
+
+
+} // namespace cbm::algo::trd
diff --git a/algo/detectors/trd/Hit.h b/algo/detectors/trd/Hit.h
new file mode 100644
index 0000000000..2f99d8e6fe
--- /dev/null
+++ b/algo/detectors/trd/Hit.h
@@ -0,0 +1,172 @@
+/* Copyright (C) 2024 Facility for Antiproton and Ion Research in Europe, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Dominik Smith [committer], Matus Kalisky, Florian Uhlig, Andrey Lebedev */
+
+/**
+ * \file Hit.h
+ * \author Dominik Smith <d.smith@gsi.de>
+ * \date 2024
+ * \brief A light-weight TRD hit class for online reconstruction, based on CbmTrdHit. .
+ **/
+
+#pragma once
+
+#include "CbmTrdAddress.h" // for CbmTrdAddress
+#include "Math/Rotation3D.h"
+#include "Math/Vector3Dfwd.h"
+
+#include <Rtypes.h> // for CLRBIT, SETBIT, TESTBIT, ClassDef
+#include <RtypesCore.h> // for Double32_t
+
+#include <cstdint>
+#include <string> // for string
+
+namespace cbm::algo::trd
+{
+
+ /** \class Hit
+ * \brief A light-weight TRD hit class for online reconstruction, based on CbmTrdHit. .
+ */
+
+ class Hit {
+ public:
+ enum HitDef
+ {
+ kType = 0 ///< set type of pad layout
+ ,
+ kMaxType ///< set type of pad on which the maximum charge is found
+ ,
+ kRowCross ///< mark hit defined by 2 clusters
+ ,
+ kOvfl ///< mark over-flow in the data
+ };
+
+
+ /**
+ * \brief Default constructor.
+ */
+ Hit();
+ // Hit(const Hit&);
+ Hit& operator=(const Hit&);
+
+ /**
+ * \brief Standard constructor.
+ *\param address Unique detector ID
+ *\param pos Position in global c.s. [cm]
+ *\param dpos Errors of position in global c.s. [cm]
+ *\param dxy XY correlation of the hit
+ *\param refId Index of digi or cluster
+ *\param eLoss TR + dEdx
+ **/
+ Hit(int32_t address, const ROOT::Math::XYZVector& pos, const ROOT::Math::XYZVector& dpos, double dxy, int32_t refId,
+ double eLoss, double time = 0., double timeError = 0.);
+
+ /**
+ * \brief Standard constructor.
+ * \param[in] address Detector unique identifier.
+ * \param[in] x X position of the hit [cm].
+ * \param[in] y Y position of the hit [cm].
+ * \param[in] z Z position of the hit [cm].
+ * \param[in] dx X position error of the hit [cm].
+ * \param[in] dy Y position error of the hit [cm].
+ * \param[in] dz Z position error of the hit [cm].
+ * \param[in] dxy X-Y covariance of the hit [cm**2].
+ * \param[in] refId Some reference ID.
+ * \param[in] time Hit time [ns].
+ * \param[in] timeError Error of hit time [ns].
+ **/
+ Hit(int32_t address, double x, double y, double z, double dx, double dy, double dz, double dxy, int32_t refId,
+ double time = -1., double timeError = -1.);
+
+
+ /** \brief Destructor. */
+ ~Hit(){};
+
+ /* Accessors */
+ int32_t GetPlaneId() const { return CbmTrdAddress::GetLayerId(GetAddress()); }
+ double GetELoss() const { return fELoss; }
+ bool GetClassType() const { return TESTBIT(fDefine, kType); }
+ bool GetMaxType() const { return TESTBIT(fDefine, kMaxType); }
+ bool HasOverFlow() const { return TESTBIT(fDefine, kOvfl); }
+ bool IsRowCross() const { return TESTBIT(fDefine, kRowCross); }
+ bool IsUsed() const { return (GetRefId() < 0); }
+ double GetX() const { return fX; }
+ double GetY() const { return fY; }
+ double GetZ() const { return fZ; }
+ double GetDx() const { return fDx; }
+ double GetDy() const { return fDy; }
+ double GetDz() const { return fDz; }
+ double GetDxy() const { return fDxy; }
+ int32_t GetRefId() const { return fRefId; }
+ int32_t GetAddress() const { return fAddress; }
+ double GetTime() const { return fTime; }
+ double GetTimeError() const { return fTimeError; }
+
+ /**
+ * \brief Copies hit position to pos.
+ * \param pos Output hit position.
+ */
+ void Position(ROOT::Math::XYZVector& pos) const;
+
+ /**
+ * \brief Copies hit position error to pos.
+ * \param pos Output hit position error.
+ */
+ void PositionError(ROOT::Math::XYZVector& dpos) const;
+
+ /* Setters */
+ void SetELoss(double loss) { fELoss = loss; }
+ /** \brief Mark overflow in one or more digits which define the hit.*/
+ void SetOverFlow(bool set = true) { set ? SETBIT(fDefine, kOvfl) : CLRBIT(fDefine, kOvfl); }
+ /** \brief Mark hit reconstructed between pad rows.*/
+ void SetRowCross(bool set = true) { set ? SETBIT(fDefine, kRowCross) : CLRBIT(fDefine, kRowCross); }
+ /** \brief Type of pad layout used in reconstruction R[0], T[1]*/
+ void SetClassType(bool set = true) { set ? SETBIT(fDefine, kType) : CLRBIT(fDefine, kType); }
+ /** \brief Extra bool definition for the hit (e.g. the type of maximum for triangular pads).*/
+ void SetMaxType(bool set = true) { set ? SETBIT(fDefine, kMaxType) : CLRBIT(fDefine, kMaxType); }
+ void SetX(double x) { fX = x; }
+ void SetY(double y) { fY = y; }
+ void SetZ(double z) { fZ = z; }
+ void SetDx(double dx) { fDx = dx; }
+ void SetDy(double dy) { fDy = dy; }
+ void SetDz(double dz) { fDz = dz; }
+ void SetDxy(double dxy) { fDxy = dxy; }
+ void SetRefId(int32_t refId) { fRefId = refId; }
+ void SetAddress(int32_t address) { fAddress = address; }
+ void SetTime(double time) { fTime = time; }
+ void SetTime(double time, double error)
+ {
+ fTime = time;
+ fTimeError = error;
+ }
+ void SetTimeError(double error) { fTimeError = error; }
+
+ /**
+ * \brief Sets position of the hit.
+ * \param pos new hit position.
+ **/
+ void SetPosition(const ROOT::Math::XYZVector& pos);
+
+ /**
+ * \breif Sets position error of the hit.
+ * \param dpos new hit position error
+ **/
+ void SetPositionError(const ROOT::Math::XYZVector& dpos);
+
+
+ private:
+ double fX, fY, fZ; ///< X, Y, Z positions of hit [cm]
+ double fDx, fDy, fDz; ///< X, Y, Z errors [cm]
+ double fDxy; ///< X-Y covariance of the hit [cm**2]
+ int32_t fRefId; ///< some reference id (usually to cluster, digi or MC point)
+ int32_t fAddress; ///< detector unique identifier
+ double fTime; ///< Hit time [ns]
+ double fTimeError; ///< Error of hit time [ns]
+
+ uint8_t fDefine; // hit extra info
+ int32_t fNeighborId; // refId in case of row cross clusters
+ Double32_t fELoss; // Energy deposit due to TR + dEdx
+ };
+
+
+} // namespace cbm::algo::trd
diff --git a/algo/detectors/trd/HitFinder.cxx b/algo/detectors/trd/HitFinder.cxx
new file mode 100644
index 0000000000..4faeaf608e
--- /dev/null
+++ b/algo/detectors/trd/HitFinder.cxx
@@ -0,0 +1,233 @@
+/* Copyright (C) 2024 Facility for Antiproton and Ion Research in Europe, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Dominik Smith [committer], Etienne Bechtel, Florian Uhlig */
+
+#include "HitFinder.h"
+
+#include "CbmTrdDigi.h"
+
+namespace cbm::algo::trd
+{
+ constexpr double HitFinder::kxVar_Value[2][5];
+ constexpr double HitFinder::kyVar_Value[2][5];
+
+ //_______________________________________________________________________________
+ HitFinder::HitFinder(HitFinderModPar params) : fParams(params) {}
+
+ //_______________________________________________________________________________
+ std::vector<Hit> HitFinder::operator()(std::vector<Cluster>* clusters)
+ {
+ const int nclusters = clusters->size();
+ std::vector<Hit> hits;
+
+ for (int icluster = 0; icluster < nclusters; icluster++) {
+ Cluster* cluster = &clusters->at(icluster);
+ auto hit = MakeHit(icluster, cluster, &cluster->GetDigis(), nclusters);
+ if (hit.GetAddress() == -1) continue;
+ hits.push_back(hit);
+ }
+ return hits;
+ }
+
+ //_______________________________________________________________________________
+ Hit HitFinder::MakeHit(int clusterId, const Cluster* cluster, const std::vector<const CbmTrdDigi*>* digis, size_t)
+ {
+ ROOT::Math::XYZVector hit_posV(0.0, 0.0, 0.0);
+
+ double xVar = 0;
+ double yVar = 0;
+ double totalCharge = 0;
+ double time = 0.;
+ int errorclass = 0.;
+ bool EB = false;
+ bool EBP = false;
+
+ size_t skipped = 0;
+ for (auto& digi : *digis) {
+ if (!digi) {
+ skipped++;
+ continue;
+ std::cout << " no digi " << std::endl;
+ }
+
+ const double digiCharge = digi->GetCharge();
+
+ // D.Smith 15.4.24: These get overwritten for each digi. Should be OR?
+ errorclass = digi->GetErrorClass();
+ EB = digi->IsFlagged(0);
+ EBP = digi->IsFlagged(1);
+
+ if (digiCharge <= 0.05) {
+ skipped++;
+ continue;
+ }
+
+ time += digi->GetTime();
+ totalCharge += digiCharge;
+
+ int digiCol;
+ int digiRow = GetPadRowCol(digi->GetAddressChannel(), digiCol);
+
+ const ROOT::Math::XYZVector& local_pad_posV = fParams.rowPar[digiRow].padPar[digiCol].pos;
+ const ROOT::Math::XYZVector& local_pad_dposV = fParams.rowPar[digiRow].padPar[digiCol].pos;
+
+ double xMin = local_pad_posV.X() - local_pad_dposV.X();
+ double xMax = local_pad_posV.X() + local_pad_dposV.X();
+ xVar += (xMax * xMax + xMax * xMin + xMin * xMin) * digiCharge;
+
+ double yMin = local_pad_posV.Y() - local_pad_dposV.Y();
+ double yMax = local_pad_posV.Y() + local_pad_dposV.Y();
+ yVar += (yMax * yMax + yMax * yMin + yMin * yMin) * digiCharge;
+
+ hit_posV.SetX(hit_posV.X() + local_pad_posV.X() * digiCharge);
+ hit_posV.SetY(hit_posV.Y() + local_pad_posV.Y() * digiCharge);
+ hit_posV.SetZ(hit_posV.Z() + local_pad_posV.Z() * digiCharge);
+ }
+ time /= digis->size() - skipped;
+
+ if (totalCharge <= 0) return Hit();
+
+ hit_posV.SetX(hit_posV.X() / totalCharge);
+ hit_posV.SetY(hit_posV.Y() / totalCharge);
+ hit_posV.SetZ(hit_posV.Z() / totalCharge);
+
+ if (EB) {
+ xVar = kxVar_Value[0][errorclass];
+ yVar = kyVar_Value[0][errorclass];
+ }
+ else {
+ if (EBP) time -= 46; //due to the event time of 0 in the EB mode and the ULong in the the digi time
+ //TODO: move to parameter file
+ xVar = kxVar_Value[1][errorclass];
+ yVar = kyVar_Value[1][errorclass];
+ }
+
+ ROOT::Math::XYZVector global = fParams.translation + fParams.rotation(hit_posV);
+
+ // TO DO: REMOVE MAGIC NUMBERS!
+ if (!EB) { // preliminary correction for angle dependence in the position
+ // reconsutrction
+ global.SetX(global.X() + (0.00214788 + global.X() * 0.000195394));
+ global.SetY(global.Y() + (0.00370566 + global.Y() * 0.000213235));
+ }
+
+ ROOT::Math::XYZVector cluster_pad_dposV(xVar, yVar, 0);
+ TransformHitError(cluster_pad_dposV); //Gets overwritten below?
+
+ // TODO: get momentum for more exact spacial error
+ if ((fParams.orientation == 1) || (fParams.orientation == 3)) {
+ cluster_pad_dposV.SetX(sqrt(fParams.padSizeErrY)); // Original version uses padSizeErrY here for some reason
+ //cluster_pad_dposV.SetX( sqrt(fParams.padSizeErrX ) ); // x-component correct?
+ }
+ else {
+ cluster_pad_dposV.SetY(sqrt(fParams.padSizeErrY));
+ }
+
+ // Set charge of incomplete clusters (missing NTs) to -1 (not deleting them because they are still relevant for tracking)
+ if (!IsClusterComplete(cluster)) totalCharge = -1.0;
+
+ return Hit(fParams.address, global, cluster_pad_dposV, 0, clusterId, totalCharge / 1e6, time,
+ double(8.5)); // TODO: move to parameter file
+ }
+
+ void HitFinder::TransformHitError(ROOT::Math::XYZVector& hitErr) const
+ {
+ double x, y;
+ x = hitErr.X();
+ y = hitErr.Y();
+
+ if ((fParams.orientation == 1) || (fParams.orientation == 3)) { // for orientations 1 or 3
+ hitErr.SetX(y); // swap errors
+ hitErr.SetY(x); // swap errors
+ }
+ }
+
+
+ double HitFinder::GetSpaceResolution(double val)
+ {
+
+ std::pair<double, double> res[12] = {
+ std::make_pair(0.5, 0.4), std::make_pair(1, 0.35), std::make_pair(2, 0.3), std::make_pair(2.5, 0.3),
+ std::make_pair(3.5, 0.28), std::make_pair(4.5, 0.26), std::make_pair(5.5, 0.26), std::make_pair(6.5, 0.26),
+ std::make_pair(7.5, 0.26), std::make_pair(8.5, 0.26), std::make_pair(8.5, 0.26), std::make_pair(9.5, 0.26)};
+
+ double selval = 0.;
+
+ for (int n = 0; n < 12; n++) {
+ if (val < res[0].first) selval = res[0].second;
+ if (n == 11) {
+ selval = res[11].second;
+ break;
+ }
+ if (val >= res[n].first && val <= res[n + 1].first) {
+ double dx = res[n + 1].first - res[n].first;
+ double dy = res[n + 1].second - res[n].second;
+ double slope = dy / dx;
+ selval = (val - res[n].first) * slope + res[n].second;
+ break;
+ }
+ }
+
+ return selval;
+ }
+
+ bool HitFinder::IsClusterComplete(const Cluster* cluster)
+ {
+ int colMin = fParams.rowPar[0].padPar.size(); //numCols
+ int rowMin = fParams.rowPar.size(); //numRows
+
+ for (int i = 0; i < cluster->GetNofDigis(); ++i) {
+ // const CbmTrdDigi* digi = CbmDigiManager::Instance()->Get<CbmTrdDigi>(cluster->GetDigi(i));
+ const CbmTrdDigi* digi = (cluster->GetDigis())[i];
+
+ int digiCol;
+ int digiRow = GetPadRowCol(digi->GetAddressChannel(), digiCol);
+
+ if (digiCol < colMin) colMin = digiCol;
+ if (digiRow < rowMin) rowMin = digiRow;
+ }
+
+ const uint16_t nCols = cluster->GetNCols();
+ const uint16_t nRows = cluster->GetNRows();
+
+ CbmTrdDigi* digiMap[nRows][nCols]; //create array on stack for optimal performance
+ memset(digiMap, 0, sizeof(CbmTrdDigi*) * nCols * nRows); //init with nullpointers
+
+ for (int i = 0; i < cluster->GetNofDigis(); ++i) {
+ // const CbmTrdDigi* digi = CbmDigiManager::Instance()->Get<CbmTrdDigi>(cluster->GetDigi(i));
+ const CbmTrdDigi* digi = (cluster->GetDigis())[i];
+
+ int digiCol;
+ int digiRow = GetPadRowCol(digi->GetAddressChannel(), digiCol);
+
+ if (digiMap[digiRow - rowMin][digiCol - colMin])
+ return false; // To be investigated why this sometimes happens (Redmin Issue 2914)
+
+ digiMap[digiRow - rowMin][digiCol - colMin] = const_cast<CbmTrdDigi*>(digi);
+ }
+
+ // check if each row of the cluster starts and ends with a kNeighbor digi
+ for (int iRow = 0; iRow < nRows; ++iRow) {
+ int colStart = 0;
+ while (digiMap[iRow][colStart] == nullptr)
+ ++colStart;
+ if (digiMap[iRow][colStart]->GetTriggerType() != static_cast<int>(CbmTrdDigi::eTriggerType::kNeighbor))
+ return false;
+
+ int colStop = nCols - 1;
+ while (digiMap[iRow][colStop] == nullptr)
+ --colStop;
+ if (digiMap[iRow][colStop]->GetTriggerType() != static_cast<int>(CbmTrdDigi::eTriggerType::kNeighbor))
+ return false;
+ }
+
+ return true;
+ }
+
+ int HitFinder::GetPadRowCol(int address, int& c)
+ {
+ c = address % fParams.rowPar[0].padPar.size();
+ return address / fParams.rowPar[0].padPar.size();
+ }
+
+} // namespace cbm::algo::trd
diff --git a/algo/detectors/trd/HitFinder.h b/algo/detectors/trd/HitFinder.h
new file mode 100644
index 0000000000..3f4475fcc9
--- /dev/null
+++ b/algo/detectors/trd/HitFinder.h
@@ -0,0 +1,68 @@
+/* Copyright (C) 2024 Facility for Antiproton and Ion Research in Europe, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Dominik Smith [committer], Etienne Bechtel, Florian Uhlig */
+
+#pragma once
+
+#include "Cluster.h"
+#include "Hit.h"
+#include "HitFinderPars.h"
+#include "Math/Rotation3D.h"
+#include "Math/Vector3Dfwd.h"
+
+#include <tuple>
+#include <vector>
+
+class CbmTrdDigi;
+
+namespace cbm::algo::trd
+{
+ /**
+ * \brief Rectangular pad module; Cluster finding and hit reconstruction algorithms
+ **/
+ class HitFinder {
+ public:
+ HitFinder(){};
+ /**
+ * \brief Constructor with placement
+ **/
+ HitFinder(HitFinderModPar params);
+ virtual ~HitFinder(){};
+
+ /* \brief Steering routine for building hits */
+ std::vector<Hit> operator()(std::vector<Cluster>* clusters);
+
+ double GetSpaceResolution(double val = 3.0);
+ bool IsClusterComplete(const Cluster* cluster);
+
+
+ protected:
+ private:
+ typedef std::tuple<int, const CbmTrdDigi*, double>
+ inputType; //digi index, pointer to digi (null if processed), digi time
+
+ HitFinder(const HitFinder& ref);
+ const HitFinder& operator=(const HitFinder& ref);
+
+ Hit MakeHit(int cId, const Cluster* c, const std::vector<const CbmTrdDigi*>* digis, size_t);
+
+ /** \brief Addressing ASIC on module based on id
+ * \param[in] id module wise ASIC identifier
+ * \return ASIC address within experiment
+ */
+ inline int GetPadRowCol(int address, int& c);
+
+ HitFinderModPar fParams; ///< Parameter container
+
+ void TransformHitError(ROOT::Math::XYZVector& hitErr) const;
+
+ // different error classes for the position resolution based on the simulation results
+ // the error classes are defined for the different module types
+ // TODO: move to parameter file
+ static constexpr double kxVar_Value[2][5] = {{0.0258725, 0.0267693, 0.0344325, 0.0260322, 0.040115},
+ {0.0426313, 0.0426206, 0.0636962, 0.038981, 0.0723851}};
+ static constexpr double kyVar_Value[2][5] = {{0.024549, 0.025957, 0.0250713, 0.0302682, 0.0291146},
+ {0.0401438, 0.0407502, 0.0397242, 0.0519485, 0.0504586}};
+ };
+
+} // namespace cbm::algo::trd
diff --git a/algo/detectors/trd/HitFinder2D.cxx b/algo/detectors/trd/HitFinder2D.cxx
new file mode 100644
index 0000000000..d3b8a0628f
--- /dev/null
+++ b/algo/detectors/trd/HitFinder2D.cxx
@@ -0,0 +1,934 @@
+/* Copyright (C) 2024 Facility for Antiproton and Ion Research in Europe, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Dominik Smith [committer], Alexandru Bercuci */
+
+#include "HitFinder2D.h"
+
+#include "log.hpp"
+
+#include <numeric>
+
+using std::cout;
+using std::endl;
+using std::vector;
+
+namespace cbm::algo::trd
+{
+
+ //////////// TO DO: fdEdep is not used!!! Will be needed at some point.
+ double HitFinder2D::fgDT[] = {4.181e-6, 1586, 24};
+
+ //_______________________________________________________________________________
+ HitFinder2D::HitFinder2D(HitFinder2DModPar params) : fParams(params), fHit(params.rowPar[0].padPar.size()) {}
+
+ //_______________________________________________________________________________
+ std::vector<Hit> HitFinder2D::operator()(std::vector<Cluster2D>* clusters)
+ {
+ const int nclusters = clusters->size();
+ std::vector<Hit> hits;
+
+ for (int icluster = 0; icluster < nclusters; icluster++) {
+ Cluster2D* cluster = &clusters->at(icluster);
+ auto hit = MakeHit(icluster, cluster, &cluster->GetDigis(), nclusters);
+ if (hit.GetAddress() == -1) continue;
+ hits.push_back(hit);
+ }
+
+ int a0, nhits = hits.size();
+ float Dx(2 * fParams.padSizeX), Dy(2 * fParams.padSizeY);
+ for (int ih(0); ih < nhits; ih++) {
+ Hit* h0 = &hits[ih];
+ if (h0->IsUsed()) continue;
+
+ for (int jh(ih + 1); jh < nhits; jh++) {
+ Hit* h1 = &hits[jh];
+ if (h1->IsUsed()) continue;
+
+ // basic check on Time
+ if (h1->GetTime() < 4000 - h0->GetTime()) continue; // TODO check with preoper time calibration
+ // skip next hits for being too far (10us) in the future
+ if (h1->GetTime() > 10000 + h0->GetTime()) break;
+
+ // basic check on Row
+ if (std::abs(h1->GetY() - h0->GetY()) > Dy) continue;
+
+ // basic check on Col
+ if (std::abs(h1->GetX() - h0->GetX()) > Dx) continue;
+
+ // go to topologic checks
+ if (!(a0 = CheckMerge(h0->GetRefId(), h1->GetRefId()))) continue;
+
+ ProjectDigis(a0 < 0 ? h0->GetRefId() : h1->GetRefId(), a0 < 0 ? h1->GetRefId() : h0->GetRefId());
+
+ // call the working algorithm
+ if (MergeHits(h0, a0)) {
+ h0->SetRowCross((bool) h1);
+ h1->SetRefId(-1);
+ }
+ }
+ }
+
+ const auto ret = std::remove_if(hits.begin(), hits.end(), [](const auto& obj) { return obj.IsUsed(); });
+ hits.erase(ret, hits.end());
+
+ // clear all calibrated digis
+ for (auto& vec : fDigis) {
+ for (auto& digiPtr : vec)
+ delete digiPtr;
+ vec.clear();
+ }
+ fDigis.clear();
+
+ return hits;
+ }
+
+ //_______________________________________________________________________________
+ int HitFinder2D::CheckMerge(int cid, int cjd)
+ {
+ /** Check topologic conditions if the 2 clusters (in digi representation) can be merged.
+ * The first pair is always from the bottom row
+ * return the anode candidate wrt boundary 1, 2, 3 for the first 3 anodes of the upper row; -1, -2, -3 for the bottom row (r0) or 0 if the check fails
+ */
+
+ bool on(kFALSE);
+ int vc[2] = {-1, -1}, vm[2] = {0};
+ double M[2] = {-1., -1.}, S[2] = {0.};
+ vector<DigiRec*>::iterator jp[2];
+
+ for (int rowId(0); rowId < 2; rowId++) {
+ double rtMax = 0.;
+ double mdMax = 0.;
+ const int vcid[2] = {cid, cjd};
+ for (auto id = fDigis[vcid[rowId]].begin(); id != fDigis[vcid[rowId]].end(); id++) {
+ int col;
+ GetPadRowCol((*id)->GetAddressChannel(), col);
+
+ // mark max position and type
+ const double t = (*id)->GetTiltCharge(on);
+ if (on && t > rtMax) {
+ vc[rowId] = col;
+ vm[rowId] = 0;
+ rtMax = t;
+ }
+ const double r = (*id)->GetRectCharge(on);
+ if (on && r > rtMax) {
+ vc[rowId] = col;
+ vm[rowId] = 1;
+ rtMax = r;
+ }
+
+ double m = 0.;
+ double d = 0.;
+ if (!rowId) { // compute TR pairs on the bottom row
+ m = 0.5 * (t + r);
+ d = r - t;
+ }
+ else { // compute RT pairs on the upper row
+ auto jd = std::next(id);
+ double T = 0;
+ if (jd != fDigis[vcid[rowId]].end()) T = (*jd)->GetTiltCharge(on);
+ m = 0.5 * (r + T);
+ d = r - T;
+ }
+ if (std::abs(m) > 0.) d = 1.e2 * d / m;
+ if (m > mdMax) {
+ mdMax = m;
+ M[rowId] = m;
+ S[rowId] = d;
+ jp[rowId] = id;
+ }
+ }
+ }
+ int rowMax = M[0] > M[1] ? 0 : 1;
+
+ // basic check on col of the max signal
+ const int dc = vc[1] - vc[0];
+ if (dc < 0 || dc > 1) return 0;
+
+ // special care for both tilt maxima : the TT case
+ // recalculate values on the min row on neighbor column
+ if (!vm[0] && !vm[1]) {
+ if (rowMax == 0) { // modify r=1
+ double r = 0.;
+ double T = 0.;
+ if (M[1] >= 0) {
+ if (jp[1] != fDigis[cjd].end()) jp[1]++;
+ if (jp[1] != fDigis[cjd].end()) {
+ r = (*jp[1])->GetRectCharge(on);
+ jp[1]++;
+ if (jp[1] != fDigis[cjd].end()) T = (*jp[1])->GetTiltCharge(on);
+ }
+ }
+ M[1] = 0.5 * (r + T);
+ S[1] = r - T;
+ }
+ else { // modify r=0
+ double r = 0.;
+ double t = 0.;
+ if (M[0] >= 0) {
+ if (jp[0] != fDigis[cid].begin()) jp[0]--;
+ if (jp[0] != fDigis[cid].begin()) {
+ r = (*jp[0])->GetRectCharge(on);
+ t = (*jp[0])->GetTiltCharge(on);
+ }
+ }
+ M[0] = 0.5 * (t + r);
+ S[0] = r - t;
+ }
+ }
+ rowMax = M[0] > M[1] ? 0 : 1;
+
+ // Build the ratio of the charge
+ const float mM = M[rowMax ? 0 : 1] / M[rowMax];
+ const float mS = std::abs(S[rowMax]), mM_l[3] = {0.15, 0.5, 1}, mM_r[3] = {0, 0.28, 1}, mS_r[3] = {43, 27, 20};
+ float dSdM[2], S0[2];
+
+ for (int i(0); i < 2; i++) {
+ dSdM[i] = (mS_r[i + 1] - mS_r[i]) / (mM_r[i + 1] - mM_r[i]);
+ S0[i] = mS_r[i] - dSdM[i] * mM_r[i];
+ }
+ const int irange = mM < mM_r[1] ? 0 : 1;
+ if (mS > S0[irange] + dSdM[irange] * mM) return 0;
+
+ for (int ia(0); ia < 3; ia++) {
+ if (mM < mM_l[ia]) return (rowMax ? 1 : -1) * (3 - ia);
+ }
+ return 0;
+ }
+
+
+ //_______________________________________________________________________________
+ bool HitFinder2D::MergeHits(Hit* h, int a0)
+ {
+ int n0(fHit.fSignal.size() - 2);
+ auto [dx, dy] = fHit.GetDxDy(n0);
+
+ int typ(fHit.GetHitClass());
+ // get position correction [pw]
+ double xcorr = fHit.GetXcorr(dx, typ, 1) / fParams.padSizeX, xcorrBias(xcorr);
+ if (fHit.IsBiasX()) {
+ typ = fHit.GetHitRcClass(a0);
+ int xmap = fHit.vyM & 0xff;
+ switch (n0) {
+ case 4:
+ if (dx < 0)
+ xcorrBias += (fHit.IsBiasXleft() ? -0.12 : 0.176);
+ else
+ xcorrBias += (xmap == 53 || xmap == 80 || xmap == 113 || xmap == 117 ? -0.176 : 0.12);
+ break;
+ case 5:
+ case 7:
+ if (typ < 0) break;
+ if (xmap == 50 || xmap == 58 || xmap == 146 || xmap == 154) {
+ if (typ == 2)
+ xcorr += 0.0813;
+ else if (typ == 3) {
+ xcorr -= 0.0813;
+ typ = 2;
+ }
+ dx -= xcorr;
+ fHit.RecenterXoffset(dx);
+ xcorrBias = fHit.GetXcorr(dx, typ, 2) / fParams.padSizeX;
+ }
+ else {
+ dx -= xcorr;
+ fHit.RecenterXoffset(dx);
+ if (typ < 2)
+ xcorrBias = fHit.GetXcorr(dx, typ, 2) / fParams.padSizeX;
+ else if (typ == 2)
+ xcorrBias = 0.12;
+ else if (typ == 3)
+ xcorrBias = -0.12;
+ }
+ break;
+ default:
+ if (typ < 0)
+ break;
+ else if (typ == 2)
+ xcorr += 0.0813;
+ else if (typ == 3) {
+ xcorr -= 0.0813;
+ typ = 2;
+ }
+ dx -= xcorr;
+ fHit.RecenterXoffset(dx);
+ xcorrBias = fHit.GetXcorr(dx, typ, 2) / fParams.padSizeX;
+ break;
+ }
+ }
+ else {
+ if (typ) xcorrBias += (dx < 0 ? 1 : -1) * 0.0293;
+ }
+ std::tie(dx, dy) = CorrectPosition(dx, dy, xcorrBias);
+
+ double edx(1), edy(1), edt(60), time(-21), tdrift(100), e(200);
+ CalibrateHit(h, dx, dy, edx, edy, edt, time, tdrift, e);
+
+ return kTRUE;
+ }
+
+ //_______________________________________________________________________________
+ bool HitFinder2D::BuildHit(Hit* h)
+ {
+ const int n0 = fHit.fSignal.size() - 2;
+ auto [dx, dy] = fHit.GetDxDy(n0);
+
+ // get position correction
+ const double xcorr = fHit.GetXcorr(dx, fHit.GetHitClass()) / fParams.padSizeX;
+ std::tie(dx, dy) = CorrectPosition(dx, dy, xcorr);
+
+ // get anode wire offset
+ for (int ia = 0; ia <= NANODE; ia++) {
+ const float ya = -1.35 + ia * 0.3; //anode position in local pad coordinates
+ if (dy > ya + 0.15) continue;
+ break;
+ }
+
+ // Error parametrization X : parabolic model on cl size
+ const double parX[] = {0.713724, -0.318667, 0.0366036};
+ const double parY[] = {0.0886413, 0., 0.0435141};
+ const int nex = std::min(n0, 7);
+ const double edx = (n0 < 3) ? 1. : parX[0] + parX[1] * nex + parX[2] * nex * nex;
+ const double edy = (n0 < 3) ? 1. : parY[0] + parY[2] * dy * dy;
+ const double edt = (n0 < 3) ? 60. : 26.33; // should be parametrized as function of da TODO
+
+ // COMPUTE TIME
+ double t_avg = 0.;
+ for (int idx = 1; idx <= n0; idx++) {
+ HitFactory2D::signal& sig = fHit.fSignal[idx];
+ const double dtFEE = fgDT[0] * (sig.s - fgDT[1]) * (sig.s - fgDT[1]) / fClk;
+ t_avg += (sig.t - dtFEE);
+ if (sig.xe > 0) sig.x += dy / fParams.padSizeY;
+ }
+
+ const double time = (n0 > 1) ? t_avg / n0 : -21.;
+ const double tdrift = 100.; // should depend on Ua
+
+ // COMPUTE ENERGY (conversion to GeV)
+ const double e_avg = 1.e-9
+ * std::accumulate(fHit.fSignal.begin(), fHit.fSignal.end(), 0,
+ [](int sum, const auto& sig) { return sum + sig.s; });
+
+ CalibrateHit(h, dx, dy, edx, edy, edt, time, tdrift, e_avg);
+ return kTRUE;
+ }
+
+ //_______________________________________________________________________________
+ std::pair<double, double> HitFinder2D::CorrectPosition(double dx, double dy, const double xcorr)
+ {
+ dx -= xcorr;
+ fHit.RecenterXoffset(dx);
+ dy = dx - dy;
+ fHit.RecenterYoffset(dy);
+
+ const double ycorr = fHit.GetYcorr(dy) / fParams.padSizeY;
+ dy += ycorr;
+ fHit.RecenterYoffset(dy);
+ dx *= fParams.padSizeX;
+ dy *= fParams.padSizeY;
+ return std::make_pair(dx, dy);
+ }
+
+ //_______________________________________________________________________________
+ void HitFinder2D::CalibrateHit(Hit* h, const double dx, const double dy, const double edx, const double edy,
+ const double edt, const double time, const double tdrift, const double eloss)
+ {
+ const ROOT::Math::XYZVector& local_pad = fParams.rowPar[fHit.vrM].padPar[fHit.vcM].pos;
+ const ROOT::Math::XYZVector local(local_pad.X() + dx, local_pad.Y() + dy, local_pad.Z());
+ const ROOT::Math::XYZVector global = fParams.translation + fParams.rotation(local);
+ h->SetX(global.X());
+ h->SetY(global.Y());
+ h->SetZ(global.Z());
+
+ h->SetDx(edx);
+ h->SetDy(edy);
+ h->SetDz(0.);
+ h->SetDxy(0.);
+ h->SetTime(fClk * (fHit.vt0 + time) - tdrift + 30.29 + fHitTimeOff, edt);
+ h->SetELoss(eloss);
+ h->SetClassType();
+ h->SetMaxType(fHit.IsMaxTilt());
+ h->SetOverFlow(fHit.HasOvf());
+ }
+
+ //_______________________________________________________________________________
+ Hit HitFinder2D::MakeHit(int ic, const Cluster2D* cl, std::vector<const CbmTrdDigi*>* digis, size_t nclusters)
+ {
+ if (fDigis.empty()) {
+ fDigis.resize(nclusters);
+ }
+
+ Hit hit;
+ if (!LoadDigis(digis, ic)) return hit;
+ if (!ProjectDigis(ic)) return hit;
+
+ hit.SetAddress(fParams.address);
+ hit.SetRefId(ic);
+ BuildHit(&hit);
+ return hit;
+ }
+
+ //_______________________________________________________________________________
+ int HitFinder2D::LoadDigis(vector<const CbmTrdDigi*>* din, int cid)
+ {
+ /** Load RAW digis info in calibration aware strucuture CbmTrdDigiReco
+ * Do basic sanity checks; also incomplete adjacent digi and if found merge them.
+ */
+
+ int last_col(-1);
+ for (auto i = din->begin(); i != din->end(); i++) {
+ if ((*i) == nullptr) continue;
+ int colT(-1), colR(-1);
+ const CbmTrdDigi* dgT = (*i);
+ GetPadRowCol(dgT->GetAddressChannel(), colT);
+
+ // check column order for cluster /// TO DO: we can probably drop this check
+ if (last_col >= 0 && colT != last_col + 1) {
+ L_(error) << "TrdModuleRec2D::LoadDigis : digis in cluster " << cid << " not in increasing order !";
+ return 0;
+ }
+
+ last_col = colT;
+ const CbmTrdDigi* dgR = nullptr;
+
+ auto j = std::next(i);
+ if (j != din->end()) {
+ dgR = (*j);
+ GetPadRowCol(dgR->GetAddressChannel(), colR);
+ }
+ if (colR == colT && dgR != nullptr) {
+ fDigis[cid].emplace_back(new DigiRec(*dgT, *dgR));
+ (*j) = nullptr;
+ }
+ else
+ fDigis[cid].emplace_back(new DigiRec(*dgT));
+ }
+ return fDigis[cid].size();
+ }
+
+ //_______________________________________________________________________________
+ int HitFinder2D::ProjectDigis(int cid, int cjd)
+ {
+ /** Load digis information in working vectors Digis are represented in the normal coordinate system of
+ * (pad width [pw], DAQ time [clk], signal [ADC chs]) with rectangular signals at integer
+ * positions.
+ */
+
+ if (fDigis[cid].empty()) {
+ L_(debug) << "TrdModuleRec2D::ProjectDigis : Request cl id " << cid << " not found.";
+ return 0;
+ }
+
+ std::vector<HitFactory2D::signal>& hitSig = fHit.fSignal;
+ fHit.reset();
+
+ bool on(0); // flag signal transmition
+ int n0(0), nsr(0), // no of signals in the cluster (all/rect)
+ NR(0), // no of rect signals/channel in case of RC
+ NT(0), // no of tilt signals/channel in case of RC
+ ovf(1); // over-flow flag for at least one of the digis
+ Char_t dt0(0); // cluster time offset wrt arbitrary t0
+ double err, // final error parametrization for signal
+ xc(-0.5), // running signal-pad position
+ max(0.); // max signal
+ const double re = 100.; // rect signal
+ const double te = 100.; // tilt signal
+
+ int j(0), col(-1), col0(0), col1(0), step(0), row1;
+
+ // link second row if needed
+ vector<DigiRec*>::iterator i1;
+ if (cjd >= 0) {
+ if (fDigis[cjd].empty()) {
+ L_(debug) << "TrdModuleRec2D::ProjectDigis : Request cl id " << cjd << " not found. Skip second row.";
+ cjd = -1;
+ }
+ else
+ i1 = fDigis[cjd].begin();
+ }
+
+ for (auto i = fDigis[cid].begin(); i != fDigis[cid].end(); i++, j++) {
+ const DigiRec* dg = (*i);
+ const DigiRec* dg0 = nullptr;
+ const DigiRec* dg1 = nullptr;
+
+ // initialize
+ if (col < 0) {
+ fHit.vrM = GetPadRowCol(dg->GetAddressChannel(), col);
+ fHit.vt0 = dg->GetTimeDAQ(); // set arbitrary t0 to avoid double digis loop
+ }
+ GetPadRowCol(dg->GetAddressChannel(), col0);
+ int nt = 0; // no of tilt signals/channel in case of RC
+ int nr = 0; // no of rect signals/channel in case of RC
+
+ // read calibrated signals
+ double t = dg->GetTiltCharge(on);
+ if (on) nt = 1;
+ double r = dg->GetRectCharge(on);
+ if (on) nr = 1;
+ // look for matching neighbor digis in case of pad row cross
+ if (cjd >= 0) {
+ if ((dg0 = (i1 != fDigis[cjd].end()) ? (*i1) : nullptr)) {
+ row1 = GetPadRowCol(dg0->GetAddressChannel(), col1);
+ if (!step) step = fHit.vrM - row1;
+ if (col1 == col0) {
+ r += dg0->GetRectCharge(on);
+ if (on) nr++;
+ }
+ else
+ dg0 = nullptr;
+ }
+ if (step == 1 && i1 != fDigis[cjd].begin()) {
+ dg1 = (*(i1 - 1));
+ GetPadRowCol(dg1->GetAddressChannel(), col1);
+ if (col1 == col0 - 1) {
+ t += dg1->GetTiltCharge(on);
+ if (on) nt++;
+ }
+ else
+ dg1 = nullptr;
+ }
+ if (step == -1 && i1 != fDigis[cjd].end() && i1 + 1 != fDigis[cjd].end()) {
+ dg1 = (*(i1 + 1));
+ GetPadRowCol(dg1->GetAddressChannel(), col1);
+ if (col1 == col0 + 1) {
+ t += dg1->GetTiltCharge(on);
+ if (on) nt++;
+ }
+ else
+ dg1 = nullptr;
+ }
+ if (dg0) i1++;
+ }
+
+ //TO DO: two blocks below might be mergable
+ // process tilt signal/time
+ char ddt = dg->GetTiltTime() - fHit.vt0; // signal time offset wrt prompt
+ if (ddt < dt0) dt0 = ddt;
+ if (abs(t) > 0) {
+ if (nt > 1) t *= 0.5;
+ err = te * (nt > 1 ? 0.707 : 1);
+ if (dg->HasTiltOvf()) {
+ ovf = -1;
+ err = 150.;
+ }
+ if (t > max) {
+ max = t;
+ fHit.vcM = j;
+ fHit.SetMaxTilt(1);
+ fHit.viM = hitSig.size();
+ }
+ }
+ else
+ err = 300.;
+
+ hitSig.emplace_back(t, err, ddt, xc, 0.035);
+ xc += 0.5;
+
+ // process rect signal/time
+ ddt = dg->GetRectTime() - fHit.vt0; // signal time offset wrt prompt
+ if (ddt < dt0) dt0 = ddt;
+ if (abs(r) > 0) {
+ nsr++;
+ if (nr > 1) r *= 0.5;
+ err = re * (nr > 1 ? 0.707 : 1);
+ if (dg->HasRectOvf()) {
+ ovf = -1;
+ err = 150.;
+ }
+ if (r > max) {
+ max = r;
+ fHit.vcM = j;
+ fHit.SetMaxTilt(0);
+ fHit.viM = hitSig.size();
+ }
+ }
+ else
+ err = 300.;
+
+ hitSig.emplace_back(r, err, ddt, xc, 0.);
+ xc += 0.5;
+ NR += nr;
+ NT += nt;
+ }
+
+ //TO DO: two blocks below might be mergable
+ // add front and back anchor points if needed
+ if (std::abs(hitSig[0].s) > 1.e-3) {
+ xc = hitSig[0].x;
+ char ddt = hitSig[0].t;
+ hitSig.emplace(hitSig.begin(), 0., 300., ddt, xc - 0.5, 0.);
+ fHit.viM++;
+ }
+ int n(hitSig.size() - 1);
+ if (std::abs(hitSig[n].s) > 1.e-3) {
+ xc = hitSig[n].x + 0.5;
+ char ddt = hitSig[n].t;
+ hitSig.emplace_back(0., 300., ddt, xc, 0.035);
+ }
+
+ n0 = hitSig.size() - 2;
+ // compute cluster asymmetry
+ int nR = n0 + 1 - fHit.viM;
+ if (nR == fHit.viM) {
+ fHit.SetSymmHit(1);
+ if (hitSig.size() % 2) {
+ double LS(0.), RS(0.);
+ for (UChar_t idx(0); idx < fHit.viM; idx++)
+ LS += hitSig[idx].s;
+ for (uint idx(fHit.viM + 1); idx < hitSig.size(); idx++)
+ RS += hitSig[idx].s;
+ fHit.SetLeftSgn(LS < RS ? 0 : 1);
+ }
+ }
+ else {
+ fHit.SetSymmHit(0);
+ if (fHit.viM < nR)
+ fHit.SetLeftHit(0);
+ else if (fHit.viM > nR)
+ fHit.SetLeftHit(1);
+ }
+ // recenter time and space profile
+ fHit.vt0 += dt0;
+ for (auto& sig : hitSig) {
+ sig.t -= dt0;
+ sig.x -= fHit.vcM;
+ }
+ fHit.vcM += col;
+
+ // check if all signals have same significance
+ int nmiss = 2 * nsr - NR;
+ if (cjd >= 0 && nmiss) {
+ fHit.SetBiasX(1);
+ for (UChar_t idx(1); idx < fHit.viM; idx++) {
+ if (hitSig[idx].xe > 0.) continue; // select rect signals
+ if (hitSig[idx].se > re * 0.8) fHit.SetBiasXleft(1);
+ }
+ if (hitSig[fHit.viM].xe <= 0. && hitSig[fHit.viM].se > re * 0.8) fHit.SetBiasXmid(1);
+ for (UChar_t idx(fHit.viM + 1); idx < hitSig.size() - 1; idx++) {
+ if (hitSig[idx].xe > 0.) continue; // select rect signals
+ if (hitSig[idx].se > re * 0.8) fHit.SetBiasXright(1);
+ }
+ }
+ else
+ fHit.SetBiasX(0);
+ nmiss = 2 * n0 - 2 * nsr - NT;
+ if (cjd >= 0 && nmiss) {
+ fHit.SetBiasY();
+ for (UChar_t idx(1); idx < fHit.viM; idx++) {
+ if (hitSig[idx].xe > 0. && hitSig[idx].se > te * 0.8) fHit.SetBiasYleft(1); // select tilt signals
+ }
+ if (hitSig[fHit.viM].xe > 0. && hitSig[fHit.viM].se > te * 0.8) fHit.SetBiasYmid(1);
+ for (UChar_t idx(fHit.viM + 1); idx < hitSig.size() - 1; idx++) {
+ if (hitSig[idx].xe > 0. && hitSig[idx].se > te * 0.8) fHit.SetBiasYright(1); // select tilt signals
+ }
+ }
+ else
+ fHit.SetBiasY(0);
+
+ if (ovf < 0) fHit.SetOvf();
+ return ovf * (hitSig.size() - 2);
+ }
+
+ int HitFinder2D::GetPadRowCol(int address, int& c)
+ {
+ c = address % fParams.rowPar[0].padPar.size();
+ return address / fParams.rowPar[0].padPar.size();
+ }
+
+ ///////////////////////////////////////////
+
+ //_______________________________________________________________________________
+ std::pair<double, double> HitFactory2D::GetDxDy(const int n0)
+ {
+ double dx, dy;
+ switch (n0) {
+ case 1:
+ if (IsMaxTilt()) { // T
+ dx = -0.5;
+ dy = 0;
+ }
+ else { // R
+ dx = 0.5;
+ dy = 0;
+ }
+ break;
+ case 2:
+ if (IsOpenLeft() && IsOpenRight()) { // RT
+ dx = viM == 1 ? 0. : -1;
+ dy = -0.5;
+ }
+ else { // TR
+ dx = 0.;
+ dy = 0.5;
+ }
+ break;
+ case 3:
+ if (IsMaxTilt() && !IsSymmHit()) { // TRT asymm
+ dx = viM == 1 ? 0. : -1;
+ dy = GetYoffset();
+ }
+ else if (!IsMaxTilt() && IsSymmHit()) { // TRT symm
+ dx = 0.;
+ dy = GetYoffset();
+ }
+ else if (IsMaxTilt() && IsSymmHit()) { // RTR symm
+ dx = GetXoffset();
+ dy = 0.;
+ }
+ else if (!IsMaxTilt() && !IsSymmHit()) { // RTR asymm
+ dx = GetXoffset();
+ dy = viM == 1 ? -0.5 : 0.5;
+ }
+ break;
+ default:
+ dx = GetXoffset();
+ dy = GetYoffset();
+ break;
+ }
+ RecenterXoffset(dx);
+ return std::make_pair(dx, dy);
+ }
+
+ //_______________________________________________________________________________
+ double HitFactory2D::GetXoffset(int n0) const
+ {
+ double dx(0.), R(0.);
+ int n(n0 ? n0 : fSignal.size());
+ for (int ir(0); ir < n; ir++) {
+ if (fSignal[ir].xe > 0) continue; // select rectangular coupling
+ R += fSignal[ir].s;
+ dx += fSignal[ir].s * fSignal[ir].x;
+ }
+ if (std::abs(R) > 0) return dx / R;
+ //L_(debug) << "HitFactory2D::GetXoffset : Null total charge for hit size " << n;
+ return 0.;
+ }
+
+ //_______________________________________________________________________________
+ double HitFactory2D::GetYoffset(int n0) const
+ {
+ double dy(0.), T(0.);
+ int n(n0 ? n0 : fSignal.size());
+ for (int it(0); it < n; it++) {
+ if (fSignal[it].xe > 0) { // select tilted coupling
+ T += fSignal[it].s;
+ dy += fSignal[it].s * fSignal[it].x;
+ }
+ }
+ if (std::abs(T) > 0) return dy / T;
+ // L_(debug) << "HitFactory2D::GetYoffset : Null total charge for hit size " << n;
+ //if (CWRITE(1))
+ return 0.;
+ }
+
+ //_______________________________________________________________________________
+ void HitFactory2D::RecenterXoffset(double& dx)
+ {
+ /** Shift graph representation to fit dx[pw] in [-0.5, 0.5]
+ */
+
+ if (dx >= -0.5 && dx < 0.5) return;
+ int ishift = int(dx - 0.5) + (dx > 0.5 ? 1 : 0);
+ if (vcM + ishift < 0)
+ ishift = -vcM;
+ else if (vcM + ishift >= nCols)
+ ishift = nCols - vcM - 1;
+
+ dx -= ishift;
+ vcM += ishift;
+ for (uint idx(0); idx < fSignal.size(); idx++)
+ fSignal[idx].x -= ishift;
+ }
+
+ //_______________________________________________________________________________
+ void HitFactory2D::RecenterYoffset(double& dy)
+ {
+ /** Shift graph representation to fit dy[ph] in [-0.5, 0.5]
+ */
+
+ if (dy >= -0.5 && dy < 0.5) return;
+ int ishift = int(dy - 0.5) + (dy > 0.5 ? 1 : 0);
+ dy -= ishift;
+ }
+
+ //_______________________________________________________________________________
+ int HitFactory2D::GetHitClass() const
+ {
+ /** Incapsulate hit classification criteria based on signal topology
+ * [0] : center hit type
+ * [1] : side hit type
+ */
+
+ int n0(fSignal.size() - 2);
+ if ((n0 == 5 && ((IsMaxTilt() && IsSymmHit()) || (!IsMaxTilt() && !IsSymmHit()))) || // TRTRT symm/asymm
+ n0 == 4 || (n0 == 3 && ((IsMaxTilt() && IsSymmHit()) || (!IsMaxTilt() && !IsSymmHit()))))
+ return 1; // RTR symm/asymm
+ else if (n0 > 5 && HasOvf())
+ return 2;
+ return 0;
+ }
+
+ //_______________________________________________________________________________
+ int HitFactory2D::GetHitRcClass(int a0) const
+ {
+ int a0m = std::abs(a0);
+ uint8_t xmap = vyM & 0xff;
+ if (a0m == 2 && IsBiasXleft() && IsBiasXright() && !IsBiasXmid())
+ return 0;
+ else if (a0m == 3 && ((IsBiasXleft() && IsBiasXright()) || xmap == 116 || xmap == 149 || xmap == 208))
+ return 1;
+ else if (!IsBiasXleft()
+ && (a0m == 2
+ || (a0m == 3 && ((!IsBiasXright() && IsBiasXmid()) || xmap == 209 || xmap == 212 || xmap == 145))))
+ return 2;
+ else if (!IsBiasXright()
+ && (a0m == 2 || (a0m == 3 && ((!IsBiasXleft() && IsBiasXmid()) || xmap == 112 || xmap == 117))))
+ return 3;
+ else
+ return -1;
+ }
+
+ //_______________________________________________________________________________
+ double HitFactory2D::GetXcorr(double dxIn, int typ, int cls) const
+ {
+ /** Give the linear interpolation of SYS correction for current position offset "dx" based
+ * on LUT calculated wrt MC EbyE data. The position offset is expresed in [pw] units
+ * while the output is in [cm]
+ */
+
+ if (typ < 0 || typ > 2) {
+ //L_(error)<< GetName() << "::GetXcorr : type in-param "<<typ<<" out of range.";
+ return 0;
+ }
+ double dx = std::abs(dxIn);
+ int ii = std::max(0, Nint(dx / fgCorrXdx) - 1), i0; // i1;
+
+ if (ii < 0 || ii > NBINSCORRX) {
+ // L_(debug) << GetName() << "::GetXcorr : LUT idx " << ii << " out of range for dx=" << dxIn;
+ return 0;
+ }
+ if (dx < fgCorrXdx * ii) {
+ i0 = std::max(0, ii - 1);
+ /*i1=ii;*/
+ }
+ else {
+ i0 = ii;
+ /*i1=TMath::Min(NBINSCORRX-1,ii+1);*/
+ }
+
+ float* xval = &fgCorrXval[typ][i0];
+ if (cls == 1)
+ xval = &fgCorrRcXval[typ][i0];
+ else if (cls == 2)
+ xval = &fgCorrRcXbiasXval[typ][i0];
+ double DDx = (xval[1] - xval[0]), a = DDx / fgCorrXdx, b = xval[0] - DDx * (i0 + 0.5);
+ return (dxIn > 0 ? 1 : -1) * b + a * dxIn;
+ }
+
+ //_______________________________________________________________________________
+ double HitFactory2D::GetYcorr(double dy, int /* cls*/) const
+ {
+ /** Process y offset. Apply systematic correction for y (MC derived).
+ * The position offset is expresed in [pw] units while the output is in [cm]
+ */
+ float fdy(1.), yoff(0.);
+ int n0(fSignal.size() - 2);
+ switch (n0) {
+ case 3:
+ fdy = fgCorrYval[0][0];
+ yoff = fgCorrYval[0][1];
+ if (IsMaxTilt() && IsSymmHit()) {
+ fdy = 0.;
+ yoff = (dy > 0 ? -1 : 1) * 1.56;
+ }
+ else if (!IsMaxTilt() && !IsSymmHit()) {
+ fdy = 0.;
+ yoff = (dy > 0 ? -1 : 1) * 1.06;
+ }
+ else if (!IsMaxTilt() && IsSymmHit()) {
+ fdy = 2.114532;
+ yoff = -0.263;
+ }
+ else /*if(IsMaxTilt()&&!IsSymmHit())*/ {
+ fdy = 2.8016010;
+ yoff = -1.38391;
+ }
+ break;
+ case 4:
+ fdy = fgCorrYval[1][0];
+ yoff = fgCorrYval[1][1];
+ if ((!IsMaxTilt() && IsLeftHit()) || (IsMaxTilt() && !IsLeftHit())) yoff *= -1;
+ break;
+ case 5:
+ case 7:
+ case 9:
+ case 11:
+ fdy = fgCorrYval[2][0];
+ yoff = fgCorrYval[2][1];
+ break;
+ case 6:
+ case 8:
+ case 10:
+ fdy = fgCorrYval[3][0];
+ yoff = fgCorrYval[3][1];
+ if ((!IsMaxTilt() && IsLeftHit()) || (IsMaxTilt() && !IsLeftHit())) yoff *= -1;
+ break;
+ }
+ return dy * fdy + yoff;
+ }
+
+ //_______________________________________________________________________________
+ bool HitFactory2D::IsOpenRight() const
+ {
+ int nR = fSignal.size() - 1 - viM;
+ return (nR % 2 && IsMaxTilt()) || (!(nR % 2) && !IsMaxTilt());
+ }
+
+ float HitFactory2D::fgCorrXdx = 0.01;
+ float HitFactory2D::fgCorrXval[3][NBINSCORRX] = {
+ {-0.001, -0.001, -0.002, -0.002, -0.003, -0.003, -0.003, -0.004, -0.004, -0.006, -0.006, -0.006, -0.007,
+ -0.007, -0.008, -0.008, -0.008, -0.009, -0.009, -0.011, -0.011, -0.011, -0.012, -0.012, -0.012, -0.012,
+ -0.013, -0.013, -0.013, -0.013, -0.014, -0.014, -0.014, -0.014, -0.014, -0.016, -0.016, -0.016, -0.016,
+ -0.017, -0.017, -0.017, -0.018, -0.018, -0.018, -0.018, -0.018, 0.000, 0.000, 0.000},
+ {0.467, 0.430, 0.396, 0.364, 0.335, 0.312, 0.291, 0.256, 0.234, 0.219, 0.207, 0.191, 0.172,
+ 0.154, 0.147, 0.134, 0.123, 0.119, 0.109, 0.122, 0.113, 0.104, 0.093, 0.087, 0.079, 0.073,
+ 0.067, 0.063, 0.058, 0.053, 0.049, 0.046, 0.042, 0.038, 0.036, 0.032, 0.029, 0.027, 0.024,
+ 0.022, 0.019, 0.017, 0.014, 0.013, 0.011, 0.009, 0.007, 0.004, 0.003, 0.001},
+ {0.001, 0.001, 0.001, 0.001, 0.002, 0.002, 0.001, 0.002, 0.004, 0.003, 0.002, 0.002, 0.002,
+ 0.002, 0.002, 0.002, 0.003, 0.004, 0.003, 0.004, 0.004, 0.007, 0.003, 0.004, 0.002, 0.002,
+ -0.011, -0.011, -0.012, -0.012, -0.012, -0.013, -0.013, -0.013, -0.014, -0.014, -0.014, -0.016, -0.016,
+ -0.016, -0.017, -0.017, -0.017, -0.018, -0.018, -0.018, -0.019, 0.029, 0.018, 0.001}};
+ float HitFactory2D::fgCorrYval[NBINSCORRY][2] = {{2.421729, 0.},
+ {0.629389, -0.215285},
+ {0.23958, 0.},
+ {0.151913, 0.054404}};
+ float HitFactory2D::fgCorrRcXval[2][NBINSCORRX] = {
+ {-0.00050, -0.00050, -0.00150, -0.00250, -0.00250, -0.00350, -0.00450, -0.00450, -0.00550, -0.00650,
+ -0.00650, -0.00750, -0.00850, -0.00850, -0.00850, -0.00950, -0.00950, -0.00950, -0.01050, -0.01150,
+ -0.01150, -0.01150, -0.01250, -0.01250, -0.01250, -0.01250, -0.01350, -0.01350, -0.01350, -0.01350,
+ -0.01450, -0.01450, -0.01450, -0.01550, -0.01550, -0.01550, -0.01550, -0.01650, -0.01650, -0.01550,
+ -0.01650, -0.01614, -0.01620, -0.01624, -0.01626, -0.01627, -0.01626, -0.01624, -0.01620, -0.01615},
+ {0.36412, 0.34567, 0.32815, 0.31152, 0.29574, 0.28075, 0.26652, 0.25302, 0.24020, 0.22803,
+ 0.21647, 0.21400, 0.19400, 0.18520, 0.17582, 0.16600, 0.14600, 0.13800, 0.14280, 0.14200,
+ 0.13400, 0.12600, 0.12200, 0.11000, 0.10200, 0.09400, 0.09000, 0.08600, 0.08200, 0.07400,
+ 0.07000, 0.06600, 0.06600, 0.06200, 0.05800, 0.05400, 0.05400, 0.05000, 0.04600, 0.04600,
+ 0.04200, 0.03800, 0.03800, 0.03400, 0.03400, 0.03000, 0.03000, 0.02600, 0.02200, 0.02200}};
+ float HitFactory2D::fgCorrRcXbiasXval[3][NBINSCORRX] = {
+ {0.00100, 0.00260, 0.00540, 0.00740, 0.00900, 0.01060, 0.01300, 0.01460, 0.01660, 0.01900,
+ 0.02060, 0.02260, 0.02420, 0.02700, 0.02860, 0.02980, 0.03220, 0.03340, 0.03540, 0.03620,
+ 0.03820, 0.04020, 0.04180, 0.04340, 0.04460, 0.04620, 0.04740, 0.04941, 0.05088, 0.05233,
+ 0.05375, 0.05515, 0.05653, 0.05788, 0.05921, 0.06052, 0.06180, 0.06306, 0.06430, 0.06551,
+ 0.06670, 0.06786, 0.06901, 0.07012, 0.07122, 0.07229, 0.07334, 0.07436, 0.07536, 0.07634},
+ {0.00100, 0.00380, 0.00780, 0.00900, 0.01220, 0.01460, 0.01860, 0.01940, 0.02260, 0.02540,
+ 0.02820, 0.03060, 0.03220, 0.03660, 0.03980, 0.04094, 0.04420, 0.04620, 0.04824, 0.04980,
+ 0.05298, 0.05532, 0.05740, 0.05991, 0.06217, 0.06500, 0.06540, 0.06900, 0.07096, 0.07310,
+ 0.07380, 0.07729, 0.07935, 0.08139, 0.08340, 0.08538, 0.08734, 0.08928, 0.08900, 0.09307,
+ 0.09493, 0.09340, 0.09858, 0.09620, 0.09740, 0.10386, 0.09980, 0.10726, 0.10892, 0.11056},
+ {0.00011, 0.00140, 0.00340, 0.00420, 0.00500, 0.00620, 0.00820, 0.00860, 0.01060, 0.01100,
+ 0.01220, 0.01340, 0.01500, 0.01540, 0.01700, 0.01820, 0.01900, 0.02060, 0.02180, 0.02260,
+ 0.02340, 0.02420, 0.02500, 0.02500, 0.02660, 0.02740, 0.02820, 0.02900, 0.03020, 0.03180,
+ 0.03300, 0.03260, 0.03380, 0.03460, 0.03500, 0.03580, 0.03780, 0.03820, 0.03860, 0.03900,
+ 0.04100, 0.04180, 0.04060, 0.04300, 0.04340, 0.04340, 0.04380, 0.04460, 0.04580, 0.04540}};
+
+} // namespace cbm::algo::trd
diff --git a/algo/detectors/trd/HitFinder2D.h b/algo/detectors/trd/HitFinder2D.h
new file mode 100644
index 0000000000..5398b58d8c
--- /dev/null
+++ b/algo/detectors/trd/HitFinder2D.h
@@ -0,0 +1,242 @@
+/* Copyright (C) 2024 Facility for Antiproton and Ion Research in Europe, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Dominik Smith [committer], Alexandru Bercuci */
+
+#pragma once
+
+#include "CbmTrdDigi.h"
+#include "Cluster2D.h"
+#include "DigiRec.h"
+#include "Hit.h"
+#include "HitFinder2DPars.h"
+#include "Math/Rotation3D.h"
+#include "Math/Vector3Dfwd.h"
+
+#include <cstdint>
+#include <tuple>
+#include <unordered_map>
+#include <vector>
+
+#define NBINSCORRX 50 //! no of bins in the discretized correction LUT
+#define NBINSCORRY 4 //! no of bins in the parametrization correction
+
+#define NANODE 9
+
+using std::vector;
+class DigiRec;
+class CbmTrdParFaspChannel;
+
+
+namespace cbm::algo::trd
+{
+
+ // working representation of a hit on which the reconstruction is performed
+ class HitFactory2D {
+
+ public:
+ struct signal {
+ double s; //! working copy of signals from cluster
+ double se; //! working copy of signal errors from cluster
+ char t; //! working copy of signal relative timing
+ double x; //! working copy of signal relative positions
+ double xe; //! working copy of signal relative position errors
+ signal(double _s, double _se, char _t, double _x, double _xe) : s(_s), se(_se), t(_t), x(_x), xe(_xe) {}
+ signal() : s(0), se(0), t(0), x(0), xe(0) {}
+ };
+
+ std::vector<signal> fSignal;
+ int nCols;
+
+ uint64_t vt0 = 0; //! start time of current hit [clk]
+ uint8_t vcM = 0; //! maximum col
+ uint8_t vrM = 0; //! maximum row
+ uint8_t viM = 0; //! index of maximum signal in the projection
+ uint16_t vyM = 0; //! bit map for cluster topology classification
+
+ HitFactory2D(int ncols) : nCols(ncols), vt0(0), vcM(0), vrM(0), viM(0), vyM(0) {}
+
+ void reset()
+ {
+ fSignal.clear();
+ vt0 = 0;
+ vcM = 0;
+ vrM = 0;
+ viM = 0;
+ vyM = 0;
+ }
+
+ bool HasLeftSgn() const { return TESTBIT(vyM, 3); }
+ bool HasOvf() const { return TESTBIT(vyM, 12); }
+ bool IsBiasX() const { return TESTBIT(vyM, 4); }
+ bool IsBiasXleft() const { return TESTBIT(vyM, 5); }
+ bool IsBiasXmid() const { return TESTBIT(vyM, 6); }
+ bool IsBiasXright() const { return TESTBIT(vyM, 7); }
+ bool IsBiasY() const { return TESTBIT(vyM, 8); }
+ bool IsBiasYleft() const { return TESTBIT(vyM, 9); }
+ bool IsLeftHit() const { return TESTBIT(vyM, 2); }
+ bool IsMaxTilt() const { return TESTBIT(vyM, 0); }
+ bool IsOpenLeft() const { return (viM % 2 && !IsMaxTilt()) || (!(viM % 2) && IsMaxTilt()); }
+ inline bool IsOpenRight() const;
+ bool IsSymmHit() const { return TESTBIT(vyM, 1); }
+ void SetBiasX(bool set = 1) { set ? SETBIT(vyM, 4) : CLRBIT(vyM, 4); }
+ void SetBiasXleft(bool set = 1) { set ? SETBIT(vyM, 5) : CLRBIT(vyM, 5); }
+ void SetBiasXmid(bool set = 1) { set ? SETBIT(vyM, 6) : CLRBIT(vyM, 6); }
+ void SetBiasXright(bool set = 1) { set ? SETBIT(vyM, 7) : CLRBIT(vyM, 7); }
+ void SetBiasY(bool set = 1) { set ? SETBIT(vyM, 8) : CLRBIT(vyM, 8); }
+ void SetBiasYleft(bool set = 1) { set ? SETBIT(vyM, 9) : CLRBIT(vyM, 9); }
+ void SetBiasYmid(bool set = 1) { set ? SETBIT(vyM, 10) : CLRBIT(vyM, 10); }
+ void SetBiasYright(bool set = 1) { set ? SETBIT(vyM, 11) : CLRBIT(vyM, 11); }
+ void SetLeftSgn(bool set = 1) { set ? SETBIT(vyM, 3) : CLRBIT(vyM, 3); }
+ void SetLeftHit(bool set = 1) { set ? SETBIT(vyM, 2) : CLRBIT(vyM, 2); }
+ void SetSymmHit(bool set = 1) { set ? SETBIT(vyM, 1) : CLRBIT(vyM, 1); }
+ void SetMaxTilt(bool set = 1) { set ? SETBIT(vyM, 0) : CLRBIT(vyM, 0); }
+ void SetOvf(bool set = 1) { set ? SETBIT(vyM, 12) : CLRBIT(vyM, 12); }
+ uint16_t GetHitMap() const { return vyM; }
+
+ std::pair<double, double> GetDxDy(const int n0);
+ double GetXoffset(int n0 = 0) const;
+ double GetYoffset(int n0 = 0) const;
+
+ void RecenterXoffset(double& dx);
+ /** \brief Shift graph representation to [-0.5, 0.5]
+ * \param[in] dy offset wrt center pad [ph]
+ */
+ void RecenterYoffset(double& dy);
+ /** \brief Hit classification wrt center pad
+ * \return hit type : center hit [0]; side hit [1]
+ */
+
+ int GetHitClass() const;
+ /** \brief Hit classification wrt signal bias
+ * \return hit type : center hit [0]; side hit [1]
+ */
+ int GetHitRcClass(int a0) const;
+
+ double GetXcorr(double dx, int typ, int cls = 0) const;
+ /** \brief y position correction based on LUT
+ * \param[in] dy offset computed on charge sharing expressed in [ph]
+ * \param[in] cls correction class
+ * \return correction expresed in [cm]
+ */
+ double GetYcorr(double dy, int cls = 0) const;
+ /** \brief Shift graph representation to [-0.5, 0.5]
+ * \param[in] dx offset wrt center pad [pw]
+ */
+
+ static float fgCorrXdx; //! step of the discretized correction LUT
+ static float fgCorrXval[3][NBINSCORRX]; //! discretized correction LUT
+ static float fgCorrYval[NBINSCORRY][2]; //! discretized correction params
+ static float fgCorrRcXval[2][NBINSCORRX]; //! discretized correction LUT
+ static float fgCorrRcXbiasXval[3][NBINSCORRX]; //! discretized correction LUT
+
+ // Nint function copied from TMath until better option is available
+ template<typename T>
+ inline int Nint(T x) const
+ {
+ int i;
+ if (x >= 0) {
+ i = int(x + 0.5);
+ if (i & 1 && x + 0.5 == T(i)) i--;
+ }
+ else {
+ i = int(x - 0.5);
+ if (i & 1 && x - 0.5 == T(i)) i++;
+ }
+ return i;
+ }
+ };
+
+
+ /** @class HitFinder2D
+ ** @brief Cluster finding and hit reconstruction algorithms for the TRD(2D) module.
+ ** @author Alexandru Bercucic <abercuci@niham.nipne.ro>
+ ** @since 01.02.2019
+ ** @date 01.10.2021
+ **
+ ** Extend the TRD module reconstructor for the particular case of inner TRD-2D modules.
+ ** The class is a collection of algorithms to :
+ ** - identify time and spatially correlated digis and build clusters
+ ** - identify the type of clusters and apply further merging/deconvolution
+ ** - apply FEE (channel gain, baseline) and detector (energy gain, maps, etc) calibration
+ ** - steer the calculation of hit 4D parameters (x, y, t, E)
+ **/
+
+ class HitFinder2D {
+ public:
+ typedef std::tuple<uint16_t, uint16_t, int, int, size_t> inputType; //Tuple: chT, chR, tm, row, id
+
+ /** \brief Default constructor.*/
+ HitFinder2D() : fHit(0){};
+ /** \brief Constructor with placement */
+ HitFinder2D(HitFinder2DModPar params);
+
+ virtual ~HitFinder2D(){};
+
+ /** \brief Steering routine for building hits **/
+ std::vector<Hit> operator()(std::vector<Cluster2D>* clusters);
+
+ /** \brief Load RAW digis into working array of RECO digis
+ * \param[in] din list of RAW digis in increasing order of column no
+ * \param[in] cid cluster index in the cluster array
+ * \return no of digis loaded
+ */
+
+ int LoadDigis(vector<const CbmTrdDigi*>* din, int cid);
+ int ProjectDigis(int cid, int cjd = -1);
+ /** \brief Implement topologic cuts for hit merging
+ * \return index of anode wire wrt to boundary or 0 if check fails
+ */
+ int CheckMerge(int cid, int cjd);
+ /** \brief Algorithm for hit merging
+ * \param[in] h hit to be modified by addition of hp.
+ * \param[in] a0 anode hypothesis around boundary (see CheckMerge function).
+ * \return TRUE if succesful.
+ */
+ bool MergeHits(Hit* h, int a0);
+ bool BuildHit(Hit* h);
+ /** \brief x position correction based on LUT
+ * \param[in] dx offset computed on charge sharing expressed in [pw]
+ * \param[in] typ hit type central pad [0] or edge [1]
+ * \param[in] cls correction class x wrt center [0] row-cross [1] row-cross biassed x [2]
+ * \return correction expresed in [cm]
+ */
+
+ /** \brief Time offset to synchronize TRD2D hits to the rest of detectors
+ * \param dt offset in [ns]
+ */
+ void SetHitTimeOffset(int dt) { fHitTimeOff = dt; }
+
+ protected:
+ private:
+ HitFinder2D(const HitFinder2D& ref);
+ const HitFinder2D& operator=(const HitFinder2D& ref);
+
+ Hit MakeHit(int cId, const Cluster2D* c, std::vector<const CbmTrdDigi*>* digis, size_t nclusters);
+
+ const float fClk = CbmTrdDigi::Clk(CbmTrdDigi::eCbmTrdAsicType::kFASP);
+
+ HitFinder2DModPar fParams; ///< Parameter container
+
+ /** \brief Addressing ASIC on module based on id
+ * \param[in] id module wise ASIC identifier
+ * \return ASIC address within experiment
+ */
+ inline int GetPadRowCol(int address, int& c);
+
+ std::pair<double, double> CorrectPosition(double dx, double dy, const double xcorr);
+
+ void CalibrateHit(Hit* h, const double dx, const double dy, const double edx, const double edy, const double edt,
+ const double time, const double tdrift, const double eloss);
+
+ std::vector<vector<DigiRec*>> fDigis; //!cluster-wise organized calibrated digi
+
+ int fHitTimeOff = 0; //! hit time offset for synchronization
+
+ // working representation of a hit on which the reconstruction is performed
+ HitFactory2D fHit;
+
+ static Double_t fgDT[3]; //! FASP delay wrt signal
+ };
+
+
+} // namespace cbm::algo::trd
diff --git a/algo/detectors/trd/HitFinder2DPars.h b/algo/detectors/trd/HitFinder2DPars.h
new file mode 100644
index 0000000000..142c7bb3b3
--- /dev/null
+++ b/algo/detectors/trd/HitFinder2DPars.h
@@ -0,0 +1,33 @@
+/* Copyright (C) 2022 Facility for Antiproton and Ion Research in Europe, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Dominik Smith [committer], Pierre-Alain Loizeau, Felix Weiglhofer */
+#pragma once
+
+#include "Math/Rotation3D.h"
+#include "Math/Vector3Dfwd.h"
+
+#include <vector>
+
+namespace cbm::algo::trd
+{
+
+ struct HitFinder2DPadPar {
+ ROOT::Math::XYZVector pos;
+ ROOT::Math::XYZVector posErr; /// TO DO: probably not needed
+ bool chRMasked = false;
+ bool chTMasked = false;
+ };
+
+ struct HitFinder2DRowPar {
+ std::vector<HitFinder2DPadPar> padPar;
+ };
+
+ struct HitFinder2DModPar {
+ double padSizeX = 0;
+ double padSizeY = 0;
+ uint16_t address = 0;
+ std::vector<HitFinder2DRowPar> rowPar;
+ ROOT::Math::XYZVector translation;
+ ROOT::Math::Rotation3D rotation;
+ };
+} // namespace cbm::algo::trd
diff --git a/algo/detectors/trd/HitFinderPars.h b/algo/detectors/trd/HitFinderPars.h
new file mode 100644
index 0000000000..404d6e5f4c
--- /dev/null
+++ b/algo/detectors/trd/HitFinderPars.h
@@ -0,0 +1,33 @@
+/* Copyright (C) 2022 Facility for Antiproton and Ion Research in Europe, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Dominik Smith [committer], Pierre-Alain Loizeau, Felix Weiglhofer */
+#pragma once
+
+#include "Math/Rotation3D.h"
+#include "Math/Vector3Dfwd.h"
+
+#include <vector>
+
+namespace cbm::algo::trd
+{
+ struct HitFinderPadPar {
+ ROOT::Math::XYZVector pos;
+ ROOT::Math::XYZVector posErr; /// TO DO: probably not needed
+ };
+
+ struct HitFinderRowPar {
+ std::vector<HitFinderPadPar> padPar;
+ };
+
+ struct HitFinderModPar {
+ double padSizeX = 0;
+ double padSizeY = 0;
+ double padSizeErrX = 0;
+ double padSizeErrY = 0;
+ uint16_t address = 0;
+ int orientation = 0;
+ std::vector<HitFinderRowPar> rowPar;
+ ROOT::Math::XYZVector translation;
+ ROOT::Math::Rotation3D rotation;
+ };
+} // namespace cbm::algo::trd
diff --git a/algo/detectors/trd/Hitfind.cxx b/algo/detectors/trd/Hitfind.cxx
new file mode 100644
index 0000000000..f1328a68ec
--- /dev/null
+++ b/algo/detectors/trd/Hitfind.cxx
@@ -0,0 +1,138 @@
+/* Copyright (C) 2023 Facility for Antiproton and Ion Research in Europe, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Dominik Smith [committer] */
+
+#include "Hitfind.h"
+
+#include "compat/OpenMP.h"
+#include "log.hpp"
+#include "util/TimingsFormat.h"
+
+#include <chrono>
+
+using namespace std;
+using fles::Subsystem;
+
+namespace cbm::algo::trd
+{
+ // ----- Constructor ------------------------------------------------------
+ Hitfind::Hitfind(trd::HitfindSetup setup, trd::Hitfind2DSetup setup2D) // : fStorDigi(fNbSm.size())
+ {
+
+ // Create one algorithm per module for TRD and configure it with parameters
+ for (size_t mod = 0; mod < setup.modules.size(); mod++) {
+
+ cbm::algo::trd::HitfindSetup::Mod& module = setup.modules[mod];
+ cbm::algo::trd::HitFinderModPar params;
+ params.rowPar.resize(module.rowPar.size());
+
+ const double* tra_ptr = module.translation.data();
+ const double* rot_ptr = module.rotation.data();
+ params.translation = ROOT::Math::XYZVector(tra_ptr[0], tra_ptr[1], tra_ptr[2]);
+ params.rotation = ROOT::Math::Rotation3D(&rot_ptr[0], &rot_ptr[9]);
+ params.address = module.address;
+ params.padSizeX = module.padSizeX;
+ params.padSizeY = module.padSizeY;
+ params.padSizeErrX = module.padSizeErrX;
+ params.padSizeErrY = module.padSizeErrY;
+ params.orientation = module.orientation;
+
+ for (size_t row = 0; row < module.rowPar.size(); row++) {
+ cbm::algo::trd::HitfindSetup::Row& rowPar = module.rowPar[row];
+ params.rowPar[row].padPar.resize(rowPar.padPar.size());
+
+ for (size_t col = 0; col < rowPar.padPar.size(); col++) {
+ cbm::algo::trd::HitfindSetup::Pad& pad = rowPar.padPar[col];
+ cbm::algo::trd::HitFinderPadPar& padPar = params.rowPar[row].padPar[col];
+
+ const double* pos_ptr = pad.position.data();
+ const double* posErr_ptr = pad.positionError.data();
+ padPar.pos = ROOT::Math::XYZVector(pos_ptr[0], pos_ptr[1], pos_ptr[2]);
+ padPar.posErr = ROOT::Math::XYZVector(posErr_ptr[0], posErr_ptr[1], posErr_ptr[2]);
+ }
+ }
+ fHitFind[module.address] = std::make_unique<cbm::algo::trd::HitFinder>(params);
+ fClusterBuild[module.address] = std::make_unique<cbm::algo::trd::Clusterizer>(params);
+ }
+
+ // Create one algorithm per module for TRD2D and configure it with parameters
+ for (size_t mod = 0; mod < setup2D.modules.size(); mod++) {
+
+ cbm::algo::trd::Hitfind2DSetup::Mod& module = setup2D.modules[mod];
+ cbm::algo::trd::HitFinder2DModPar params;
+ params.rowPar.resize(module.rowPar.size());
+
+ const double* tra_ptr = module.translation.data();
+ const double* rot_ptr = module.rotation.data();
+ params.translation = ROOT::Math::XYZVector(tra_ptr[0], tra_ptr[1], tra_ptr[2]);
+ params.rotation = ROOT::Math::Rotation3D(&rot_ptr[0], &rot_ptr[9]);
+ params.address = module.address;
+ params.padSizeX = module.padSizeX;
+ params.padSizeY = module.padSizeY;
+
+ for (size_t row = 0; row < module.rowPar.size(); row++) {
+ cbm::algo::trd::Hitfind2DSetup::Row& rowPar = module.rowPar[row];
+ params.rowPar[row].padPar.resize(rowPar.padPar.size());
+
+ for (size_t col = 0; col < rowPar.padPar.size(); col++) {
+ cbm::algo::trd::Hitfind2DSetup::Pad& pad = rowPar.padPar[col];
+ cbm::algo::trd::HitFinder2DPadPar& padPar = params.rowPar[row].padPar[col];
+
+ const double* pos_ptr = pad.position.data();
+ const double* posErr_ptr = pad.positionError.data();
+ padPar.pos = ROOT::Math::XYZVector(pos_ptr[0], pos_ptr[1], pos_ptr[2]);
+ padPar.posErr = ROOT::Math::XYZVector(posErr_ptr[0], posErr_ptr[1], posErr_ptr[2]);
+ padPar.chRMasked = pad.chRMasked;
+ padPar.chTMasked = pad.chTMasked;
+ }
+ }
+ fHitFind2d[module.address] = std::make_unique<cbm::algo::trd::HitFinder2D>(params);
+ fClusterBuild2d[module.address] = std::make_unique<cbm::algo::trd::Clusterizer2D>(params);
+ }
+
+ L_(info) << "--- Configured hitfinder algorithms for TRD.";
+ L_(info) << "==================================================";
+ }
+ // ----------------------------------------------------------------------------
+
+
+ // ----- Execution -------------------------------------------------------
+ Hitfind::resultType Hitfind::operator()(gsl::span<CbmTrdDigi> digiIn)
+ {
+ // --- Output data
+ resultType result = {};
+ std::vector<Hit>& hitsOut = result.first;
+
+ // Intermediate storage variables (digi, index) per module
+ std::unordered_map<int, std::vector<std::pair<CbmTrdDigi, int32_t>>> digiBuffer; //[modAddress]
+
+ const size_t maxdigis = 100000000001;
+
+ // Sort digis by module address
+ for (size_t idigi = 0; idigi < digiIn.size(); idigi++) {
+ const CbmTrdDigi* digi = &digiIn[idigi];
+ digiBuffer[digi->GetAddressModule()].emplace_back(*digi, idigi);
+ if (idigi == maxdigis) {
+ break;
+ }
+ }
+
+ // Process 1D digis
+ for (auto imod = fHitFind.begin(); imod != fHitFind.end(); imod++) {
+ auto clusters = (*fClusterBuild[imod->first])(digiBuffer[imod->first]);
+ auto hits = (*fHitFind[imod->first])(&clusters);
+ hitsOut.insert(hitsOut.end(), hits.begin(), hits.end());
+ }
+
+ // Process 2D digis
+ for (auto imod = fHitFind2d.begin(); imod != fHitFind2d.end(); imod++) {
+ auto clusters = (*fClusterBuild2d[imod->first])(digiBuffer[imod->first], 0.); // Number is TS start time (T0)
+ auto hits = (*fHitFind2d[imod->first])(&clusters);
+ hitsOut.insert(hitsOut.end(), hits.begin(), hits.end());
+ }
+
+ return result;
+ }
+ // ----------------------------------------------------------------------------
+
+} // namespace cbm::algo::trd
diff --git a/algo/detectors/trd/Hitfind.h b/algo/detectors/trd/Hitfind.h
new file mode 100644
index 0000000000..c3fcfa2745
--- /dev/null
+++ b/algo/detectors/trd/Hitfind.h
@@ -0,0 +1,84 @@
+/* Copyright (C) 2024 Facility for Antiproton and Ion Research in Europe, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Dominik Smith [committer] */
+
+#pragma once
+
+#include "CbmTrdDigi.h"
+#include "PODVector.h"
+#include "PartitionedVector.h"
+#include "trd/Clusterizer.h"
+#include "trd/Clusterizer2D.h"
+#include "trd/HitFinder.h"
+#include "trd/HitFinder2D.h"
+#include "trd/Hitfind2DSetup.h"
+#include "trd/HitfindSetup.h"
+
+#include <gsl/span>
+#include <optional>
+#include <sstream>
+#include <vector>
+
+#include <xpu/host.h>
+
+
+namespace cbm::algo::trd
+{
+
+ /** @struct HitfindMonitorData
+ ** @author Dominik Smith <d.smith@gsi.de>
+ ** @since 17 Apr 2024
+ ** @brief Monitoring data for hitfinding
+ **/
+ struct HitfindMonitorData {
+ xpu::timings fTime;
+ xpu::timings fSortTime;
+ size_t fNumDigis = 0;
+ size_t fNumHits = 0;
+
+ std::string print() const
+ {
+ std::stringstream ss;
+ ss << "Hitfind stats: num digis " << fNumDigis << ", time " << fTime.wall() << " ms ( " << fTime.throughput()
+ << " GB/s ), sort time " << fSortTime.wall() << " ms, num hits " << fNumHits << std::endl;
+ return ss.str();
+ }
+ };
+
+ /** @class Hitfind
+ ** @brief Algo class for hitfinding
+ ** @author Dominik Smith <d.smith@gsi.de>
+ ** @since 17.04.2024
+ **
+ **/
+ class Hitfind {
+
+ public:
+ //typedef std::tuple<PartitionedVector<Hit>, HitfindMonitorData, PODVector<i32>> resultType;
+ typedef std::pair<std::vector<Hit>, HitfindMonitorData> resultType;
+
+
+ /** @brief Algorithm execution
+ ** @param fles timeslice to hitfind
+ ** @return pair: digi timeslice, monitoring data
+ **
+ ** @note Modifies input digis for time calibration
+ **/
+ resultType operator()(gsl::span<CbmTrdDigi> digiIn);
+
+ /** @brief Constructor **/
+ explicit Hitfind(trd::HitfindSetup, trd::Hitfind2DSetup);
+
+ private: // members
+ /** @brief Cluster building algorithms per module */
+ std::unordered_map<int, std::unique_ptr<cbm::algo::trd::Clusterizer2D>> fClusterBuild2d;
+ std::unordered_map<int, std::unique_ptr<cbm::algo::trd::Clusterizer>> fClusterBuild;
+
+ /** @brief Hit finding algorithms per module */
+ std::unordered_map<int, std::unique_ptr<cbm::algo::trd::HitFinder2D>> fHitFind2d;
+ std::unordered_map<int, std::unique_ptr<cbm::algo::trd::HitFinder>> fHitFind;
+
+ /** @brief Intermediate storage variables (digi, index) **/
+ std::vector<std::vector<std::vector<std::pair<CbmTrdDigi, int32_t>>>> fStorDigi; //[nbType][nbSm*nbRpc][nDigis]
+ };
+} // namespace cbm::algo::trd
diff --git a/algo/detectors/trd/Hitfind2DSetup.h b/algo/detectors/trd/Hitfind2DSetup.h
new file mode 100644
index 0000000000..a045abbc75
--- /dev/null
+++ b/algo/detectors/trd/Hitfind2DSetup.h
@@ -0,0 +1,64 @@
+/* Copyright (C) 2024 Facility for Antiproton and Ion Research in Europe, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Dominik Smith [committer] */
+#pragma once
+
+#include "Definitions.h"
+#include "yaml/Property.h"
+
+#include <array>
+#include <map>
+#include <string>
+#include <vector>
+
+namespace cbm::algo::trd
+{
+
+ /**
+ * @brief Hitfind setup / Hardware cabling for TRD2D
+ * Used to create the hardware mapping for the TRD2D hitfinder.
+ */
+ struct Hitfind2DSetup {
+
+ struct Pad {
+ std::array<double, 3> position;
+ std::array<double, 3> positionError;
+ bool chRMasked;
+ bool chTMasked;
+
+ CBM_YAML_PROPERTIES(yaml::Property(&Pad::position, "position", "Local position", YAML::Flow),
+ yaml::Property(&Pad::positionError, "positionError", "Local position error", YAML::Flow),
+ yaml::Property(&Pad::chRMasked, "chRMasked", "Is channel R component masked"),
+ yaml::Property(&Pad::chTMasked, "chTMasked", "Is channel T component masked"));
+ };
+
+ struct Row {
+ std::vector<Pad> padPar;
+
+ CBM_YAML_PROPERTIES(yaml::Property(&Row::padPar, "padPar", "pad parameters"));
+ };
+
+
+ struct Mod {
+ double padSizeX;
+ double padSizeY;
+ u16 address;
+ std::vector<Row> rowPar;
+ std::array<double, 3> translation;
+ std::array<double, 9> rotation;
+
+ CBM_YAML_PROPERTIES(yaml::Property(&Mod::padSizeX, "padSizeX", "X size of pads"),
+ yaml::Property(&Mod::padSizeY, "padSizeY", "Y size of pads"),
+ yaml::Property(&Mod::address, "address", "module address"),
+ yaml::Property(&Mod::translation, "translation", "Module position", YAML::Flow),
+ yaml::Property(&Mod::rotation, "rotation", "Module rotation", YAML::Flow),
+ yaml::Property(&Mod::rowPar, "rowPar", "row parameters"));
+ };
+
+ std::vector<Mod> modules;
+
+ CBM_YAML_PROPERTIES(yaml::Property(&Hitfind2DSetup::modules, "modules", "Parameters of modules"));
+ };
+
+
+} // namespace cbm::algo::trd
diff --git a/algo/detectors/trd/HitfindSetup.h b/algo/detectors/trd/HitfindSetup.h
new file mode 100644
index 0000000000..68214432c7
--- /dev/null
+++ b/algo/detectors/trd/HitfindSetup.h
@@ -0,0 +1,65 @@
+/* Copyright (C) 2024 Facility for Antiproton and Ion Research in Europe, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Dominik Smith [committer] */
+#pragma once
+
+#include "Definitions.h"
+#include "yaml/Property.h"
+
+#include <array>
+#include <map>
+#include <string>
+#include <vector>
+
+namespace cbm::algo::trd
+{
+
+ /**
+ * @brief Hitfind setup / Hardware cabling for TRD
+ * Used to create the hardware mapping for the TRD hitfinder.
+ */
+ struct HitfindSetup {
+
+ struct Pad {
+ std::array<double, 3> position;
+ std::array<double, 3> positionError;
+
+ CBM_YAML_PROPERTIES(yaml::Property(&Pad::position, "position", "Local position", YAML::Flow),
+ yaml::Property(&Pad::positionError, "positionError", "Local position error", YAML::Flow));
+ };
+
+ struct Row {
+ std::vector<Pad> padPar;
+
+ CBM_YAML_PROPERTIES(yaml::Property(&Row::padPar, "padPar", "pad parameters"));
+ };
+
+
+ struct Mod {
+ double padSizeX;
+ double padSizeY;
+ double padSizeErrX;
+ double padSizeErrY;
+ u16 address;
+ int orientation;
+ std::vector<Row> rowPar;
+ std::array<double, 3> translation;
+ std::array<double, 9> rotation;
+
+ CBM_YAML_PROPERTIES(yaml::Property(&Mod::padSizeX, "padSizeX", "X size of pads"),
+ yaml::Property(&Mod::padSizeY, "padSizeY", "Y size of pads"),
+ yaml::Property(&Mod::padSizeErrX, "padSizeErrX", "error of X size of pads"),
+ yaml::Property(&Mod::padSizeErrY, "padSizeErrY", "error of Y size of pads"),
+ yaml::Property(&Mod::address, "address", "module address"),
+ yaml::Property(&Mod::orientation, "orientation", "module orientation"),
+ yaml::Property(&Mod::translation, "translation", "Module position", YAML::Flow),
+ yaml::Property(&Mod::rotation, "rotation", "Module rotation", YAML::Flow),
+ yaml::Property(&Mod::rowPar, "rowPar", "row parameters"));
+ };
+
+ std::vector<Mod> modules;
+
+ CBM_YAML_PROPERTIES(yaml::Property(&HitfindSetup::modules, "modules", "Parameters of modules"));
+ };
+
+} // namespace cbm::algo::trd
diff --git a/core/detectors/trd/CbmTrdParModGeo.h b/core/detectors/trd/CbmTrdParModGeo.h
index 42a4da3d66..49559f66ff 100644
--- a/core/detectors/trd/CbmTrdParModGeo.h
+++ b/core/detectors/trd/CbmTrdParModGeo.h
@@ -14,7 +14,7 @@ class TGeoPhysicalNode;
/** \brief Definition of geometry for one TRD module **/
class CbmTrdParModGeo : public CbmTrdParMod {
-public:
+ public:
CbmTrdParModGeo(const char* name = "CbmTrdParModGeo", const char* title = "TRD module geometry");
virtual ~CbmTrdParModGeo();
virtual Double_t GetDX() const;
@@ -43,7 +43,9 @@ public:
virtual void LocalToMaster(Double_t in[3], Double_t out[3]) const;
bool SetNode();
-private:
+ TGeoPhysicalNode* GetNode() { return fNode; }
+
+ private:
CbmTrdParModGeo(const CbmTrdParModGeo&);
const CbmTrdParModGeo& operator=(const CbmTrdParModGeo&);
diff --git a/macro/beamtime/mcbm2022/trd_hitfinder_run.C b/macro/beamtime/mcbm2022/trd_hitfinder_run.C
new file mode 100644
index 0000000000..0b7cb3fb6c
--- /dev/null
+++ b/macro/beamtime/mcbm2022/trd_hitfinder_run.C
@@ -0,0 +1,216 @@
+/* Copyright (C) 2020 Facility for Antiproton and Ion Research in Europe, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Pierre-Alain Loizeau, Adrian Weber [committer], Dominik Smith */
+
+
+#include <math.h>
+#include <stdio.h>
+#include <string.h>
+
+/// FIXME: Disable clang formatting to keep easy parameters overview
+/* clang-format off */
+Bool_t trd_hitfinder_run(UInt_t uRunId = 2457,
+ Int_t nTimeslices = 5,
+ // Int_t nTimeslices = 1,
+ TString sInpDir = "/home/dsmith/cbmroot/macro/beamtime/mcbm2022/data/",
+ TString sOutDir = "rec/",
+ Int_t iUnpFileIndex = -1)
+{
+ /// FIXME: Re-enable clang formatting after parameters initial values setting
+ /* clang-format on */
+
+ // --- Logger settings ----------------------------------------------------
+ TString logLevel = "INFO";
+ TString logVerbosity = "LOW";
+ // ------------------------------------------------------------------------
+
+
+ // ----- Environment --------------------------------------------------
+ TString myName = "trd_hitfinder_run"; // this macro's name for screen output
+ TString srcDir = gSystem->Getenv("VMCWORKDIR"); // top source directory
+ // ------------------------------------------------------------------------
+
+
+ // ----- In- and output file names ------------------------------------
+ /// Standardized RUN ID
+ TString sRunId = TString::Format("%04u", uRunId);
+
+ /// Initial pattern
+ //TString inFile = sInpDir + "/" + sRunId + "_faspMixed_10ts.digi";
+ TString inFile = sInpDir + "/" + sRunId + ".digi";
+
+ TString parFileOut = sOutDir + "/reco_event_mcbm_params_" + sRunId;
+ TString outFile = sOutDir + "/reco_event_mcbm_" + sRunId;
+
+ // Your folder with the Tof Calibration files;
+ TString TofFileFolder = "/home/dsmith/cbmroot/macro/run/data/";
+
+ /// Add index of splitting at unpacking level if needed
+ if (0 <= iUnpFileIndex) {
+ inFile += TString::Format("_%02u", iUnpFileIndex);
+ // the input par file is not split during unpacking!
+ parFileOut += TString::Format("_%02u", iUnpFileIndex);
+ outFile += TString::Format("_%02u", iUnpFileIndex);
+ } // if ( 0 <= uUnpFileIndex )
+ /// Add ROOT file suffix
+ inFile += ".root";
+ parFileOut += ".root";
+ outFile += ".root";
+ // ---------------------------------------------
+
+
+ // --- Load the geometry setup ----
+ // This is currently only required by the TRD (parameters)
+ cbm::mcbm::ToForceLibLoad dummy; /// Needed to trigger loading of the library as no fct dict in ROOT6 and CLING
+ TString geoSetupTag = "";
+ try {
+ geoSetupTag = cbm::mcbm::GetSetupFromRunId(uRunId);
+ }
+ catch (const std::invalid_argument& e) {
+ std::cout << "Error in mapping from runID to setup name: " << e.what() << std::endl;
+ return kFALSE;
+ }
+
+ TString geoFile = srcDir + "/macro/mcbm/data/" + geoSetupTag + ".geo.root";
+ CbmSetup* geoSetup = CbmSetup::Instance();
+ geoSetup->LoadSetup(geoSetupTag);
+
+ // You can modify the pre-defined setup by using
+ geoSetup->SetActive(ECbmModuleId::kTof, kTRUE);
+
+ //----- Load Parameters --------------------------------------------------
+ TList* parFileList = new TList();
+
+ // ----- TRD digitisation parameters -------------------------------------
+ TString geoTagTrd;
+ if (geoSetup->IsActive(ECbmModuleId::kTrd)) {
+ if (geoSetup->GetGeoTag(ECbmModuleId::kTrd, geoTagTrd)) {
+ TString paramFilesTrd(Form("%s/parameters/trd/trd_%s", srcDir.Data(), geoTagTrd.Data()));
+ std::vector<TString> paramFilesVecTrd = {"asic", "digi", "gas", "gain"};
+ for (auto parIt : paramFilesVecTrd) {
+ parFileList->Add(new TObjString(Form("%s.%s.par", paramFilesTrd.Data(), parIt.Data())));
+ }
+ }
+ for (auto parFileVecIt : *parFileList) {
+ std::cout << Form("TrdParams - %s - added to parameter file list", parFileVecIt->GetName()) << std::endl;
+ }
+ }
+
+ // ----- Timer --------------------------------------------------------
+ TStopwatch timer;
+ timer.Start();
+ // ------------------------------------------------------------------------
+
+
+ // ----- FairRunAna ---------------------------------------------------
+ FairRunAna* run = new FairRunAna();
+ FairFileSource* inputSource = new FairFileSource(inFile);
+ run->SetSource(inputSource);
+
+ FairRootFileSink* outputSink = new FairRootFileSink(outFile);
+ run->SetSink(outputSink);
+ run->SetGeomFile(geoFile);
+
+ // Define output file for FairMonitor histograms
+ TString monitorFile{outFile};
+ monitorFile.ReplaceAll("reco", "reco.monitor");
+ FairMonitor::GetMonitor()->EnableMonitor(kTRUE, monitorFile);
+ // ------------------------------------------------------------------------
+
+
+ // ----- Logger settings ----------------------------------------------
+ FairLogger::GetLogger()->SetLogScreenLevel(logLevel.Data());
+ FairLogger::GetLogger()->SetLogVerbosityLevel(logVerbosity.Data());
+
+
+ // =========================================================================
+ // === local TRD Reconstruction ===
+ // =========================================================================
+
+ bool bTRD = true;
+ bool bTRD2d = true;
+ if ((bTRD && geoSetup->IsActive(ECbmModuleId::kTrd))
+ || (bTRD2d && (geoSetup->IsActive(ECbmModuleId::kTrd2d) || geoSetup->IsActive(ECbmModuleId::kTrd)))) {
+ //Double_t triggerThreshold = 0.5e-6; // SIS100
+
+ CbmTaskTrdHitFinderParWrite* trdHitfinderPar = new CbmTaskTrdHitFinderParWrite();
+ run->AddTask(trdHitfinderPar);
+
+ CbmTaskTrdHitFinder* trdHitfinder = new CbmTaskTrdHitFinder();
+ run->AddTask(trdHitfinder);
+
+ // CbmTrdClusterFinder* trdCluster;
+ // trdCluster = new CbmTrdClusterFinder();
+ // trdCluster->SetNeighbourEnable(true, false);
+ // trdCluster->SetMinimumChargeTH(triggerThreshold);
+ // trdCluster->SetRowMerger(true);
+ // run->AddTask(trdCluster);
+ // std::cout << "-I- : Added task " << trdCluster->GetName() << std::endl;
+
+ // CbmTrdHitProducer* trdHit = new CbmTrdHitProducer();
+ // run->AddTask(trdHit);
+ // std::cout << "-I- : Added task " << trdHit->GetName() << std::endl;
+ }
+
+ // ----- Parameter database --------------------------------------------
+ std::cout << std::endl << std::endl;
+ std::cout << "-I- " << myName << ": Set runtime DB" << std::endl;
+ FairRuntimeDb* rtdb = run->GetRuntimeDb();
+ FairParRootFileIo* parIo1 = new FairParRootFileIo();
+ FairParAsciiFileIo* parIo2 = new FairParAsciiFileIo();
+ FairParRootFileIo* parIo3 = new FairParRootFileIo();
+ parIo2->open(parFileList, "in");
+ rtdb->setSecondInput(parIo2);
+ parIo3->open(parFileOut.Data(), "RECREATE");
+ // ------------------------------------------------------------------------
+ rtdb->setOutput(parIo3);
+ rtdb->saveOutput();
+ rtdb->print();
+
+ // ----- Run initialisation -------------------------------------------
+ std::cout << std::endl;
+ std::cout << "-I- " << myName << ": Initialise run" << std::endl;
+ run->Init();
+
+
+ // ----- Start run ----------------------------------------------------
+ std::cout << std::endl << std::endl;
+ std::cout << "-I- " << myName << ": Starting run" << std::endl;
+ run->Run(0, nTimeslices);
+ // ------------------------------------------------------------------------
+
+ // ----- Finish -------------------------------------------------------
+ timer.Stop();
+ FairMonitor::GetMonitor()->Print();
+ Double_t rtime = timer.RealTime();
+ Double_t ctime = timer.CpuTime();
+ std::cout << std::endl << std::endl;
+ std::cout << "Macro finished successfully." << std::endl;
+ std::cout << "Output file is " << outFile << std::endl;
+ std::cout << "Parameter file is " << parFileOut << std::endl;
+ std::cout << "Real time " << rtime << " s, CPU time " << ctime << " s" << std::endl;
+ std::cout << std::endl;
+ // ------------------------------------------------------------------------
+
+
+ // ----- Resource monitoring ------------------------------------------
+ // Extract the maximal used memory an add is as Dart measurement
+ // This line is filtered by CTest and the value send to CDash
+ FairSystemInfo sysInfo;
+ Float_t maxMemory = sysInfo.GetMaxMemory();
+ std::cout << "<DartMeasurement name=\"MaxMemory\" type=\"numeric/double\">";
+ std::cout << maxMemory;
+ std::cout << "</DartMeasurement>" << std::endl;
+
+ Float_t cpuUsage = ctime / rtime;
+ std::cout << "<DartMeasurement name=\"CpuLoad\" type=\"numeric/double\">";
+ std::cout << cpuUsage;
+ std::cout << "</DartMeasurement>" << std::endl;
+ // ------------------------------------------------------------------------
+
+ /// --- Screen output for automatic tests
+ std::cout << " Test passed" << std::endl;
+ std::cout << " All ok " << std::endl;
+
+ return kTRUE;
+}
diff --git a/reco/tasks/CMakeLists.txt b/reco/tasks/CMakeLists.txt
index 533e37764f..b07e9fea5d 100644
--- a/reco/tasks/CMakeLists.txt
+++ b/reco/tasks/CMakeLists.txt
@@ -18,6 +18,8 @@ set(SRCS
CbmTaskTofHitFinder.cxx
CbmTaskTofClusterizer.cxx
CbmTaskTofClusterizerParWrite.cxx
+ CbmTaskTrdHitFinder.cxx
+ CbmTaskTrdHitFinderParWrite.cxx
CbmTaskUnpack.cxx
)
diff --git a/reco/tasks/CbmRecoTasksLinkDef.h b/reco/tasks/CbmRecoTasksLinkDef.h
index fbff7c46c3..7242e33ab9 100644
--- a/reco/tasks/CbmRecoTasksLinkDef.h
+++ b/reco/tasks/CbmRecoTasksLinkDef.h
@@ -21,6 +21,8 @@
#pragma link C++ class CbmTaskTofHitFinder + ;
#pragma link C++ class CbmTaskTofClusterizer + ;
#pragma link C++ class CbmTaskTofClusterizerParWrite + ;
+#pragma link C++ class CbmTaskTrdHitFinderParWrite + ;
+#pragma link C++ class CbmTaskTrdHitFinder + ;
#pragma link C++ class CbmTaskTriggerDigi + ;
#pragma link C++ class CbmTaskUnpack + ;
//#pragma link C++ class cbm::algo::MainConfig + ;
diff --git a/reco/tasks/CbmTaskTrdHitFinder.cxx b/reco/tasks/CbmTaskTrdHitFinder.cxx
new file mode 100644
index 0000000000..9e505e9e3b
--- /dev/null
+++ b/reco/tasks/CbmTaskTrdHitFinder.cxx
@@ -0,0 +1,153 @@
+/* Copyright (C) 2010-2021 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Florian Uhlig [committer], Pascal Raisig, Alexandru Bercuci */
+
+#include "CbmTaskTrdHitFinder.h"
+
+#include "CbmDigiManager.h"
+
+#include <FairRootManager.h>
+#include <FairRunAna.h>
+#include <FairRuntimeDb.h>
+#include <Logger.h>
+
+// C++ Classes and includes
+#include "compat/Filesystem.h"
+#include "yaml/Yaml.h"
+
+#include <TStopwatch.h>
+#include <TVector3.h>
+
+#include <cmath>
+#include <iomanip>
+#include <iostream>
+
+using std::fixed;
+using std::left;
+using std::right;
+using std::setprecision;
+using std::setw;
+using std::stringstream;
+
+//_____________________________________________________________________
+CbmTaskTrdHitFinder::CbmTaskTrdHitFinder() : FairTask("TrdClusterFinder", 1) {}
+
+// ---- Destructor ----------------------------------------------------
+CbmTaskTrdHitFinder::~CbmTaskTrdHitFinder()
+{
+ if (fClusters) {
+ fClusters->clear();
+ delete fClusters;
+ }
+ if (fHits) {
+ fHits->clear();
+ delete fHits;
+ }
+}
+
+//_____________________________________________________________________
+void CbmTaskTrdHitFinder::SetParContainers()
+{
+ // FairRuntimeDb* rtdb = FairRunAna::Instance()->GetRuntimeDb();
+}
+
+//_____________________________________________________________________
+InitStatus CbmTaskTrdHitFinder::Init()
+{
+ CbmDigiManager::Instance()->Init();
+ if (!CbmDigiManager::Instance()->IsPresent(ECbmModuleId::kTrd)) LOG(fatal) << GetName() << "Missing Trd digi branch.";
+
+ fClusters = new std::vector<CbmTrdCluster>();
+ FairRootManager* ioman = FairRootManager::Instance();
+ ioman->RegisterAny("TrdCluster", fClusters, true);
+
+ fHits = new std::vector<CbmTrdHit>();
+ ioman->RegisterAny("TrdHit", fHits, true);
+
+ // Read hitfinder parameters and initialize algo
+ fAlgo = std::make_unique<cbm::algo::trd::Hitfind>(
+ cbm::algo::yaml::ReadFromFile<cbm::algo::trd::HitfindSetup>("TrdHitfinderPar.yaml"),
+ cbm::algo::yaml::ReadFromFile<cbm::algo::trd::Hitfind2DSetup>("TrdHitfinder2DPar.yaml"));
+
+ return kSUCCESS;
+}
+
+//_____________________________________________________________________
+void CbmTaskTrdHitFinder::Exec(Option_t* /*option*/)
+{
+
+ size_t numHitsCall = 0;
+ fClusters->clear();
+ fHits->clear();
+
+ TStopwatch timerTs;
+ timerTs.Start();
+
+ std::vector<CbmTrdDigi> digiVec;
+
+ for (int32_t iDigi = 0; iDigi < CbmDigiManager::Instance()->GetNofDigis(ECbmModuleId::kTrd); iDigi++) {
+ const CbmTrdDigi* tDigi = CbmDigiManager::Instance()->Get<CbmTrdDigi>(iDigi);
+ digiVec.push_back(CbmTrdDigi(*tDigi));
+ }
+
+ auto [hits, monitor] = (*fAlgo)(digiVec);
+ AddHits(&hits);
+
+ timerTs.Stop();
+ fProcessTime += timerTs.RealTime();
+
+ stringstream logOut;
+ logOut << setw(20) << left << GetName() << " [";
+ logOut << fixed << setw(8) << setprecision(1) << right << timerTs.RealTime() * 1000. << " ms] ";
+ logOut << "TS " << fNrTs;
+ LOG(info) << logOut.str();
+}
+
+
+//_____________________________________________________________________
+void CbmTaskTrdHitFinder::AddHits(std::vector<cbm::algo::trd::Hit>* hits)
+{
+ for (auto& hit : *hits) {
+ TVector3 pos(hit.GetX(), hit.GetY(), hit.GetZ());
+ TVector3 dpos(hit.GetDx(), hit.GetDy(), hit.GetDz());
+
+ CbmTrdHit& hitSave =
+ fHits->emplace_back(hit.GetAddress(), pos, dpos, hit.GetDxy(), hit.GetRefId(), hit.GetELoss(), hit.GetTime(),
+ hit.GetTimeError()); // TODO implement copy constructor
+ hitSave.SetOverFlow(hit.HasOverFlow());
+ hitSave.SetRowCross(hit.IsRowCross());
+ hitSave.SetClassType(hit.GetClassType());
+ hitSave.SetMaxType(hit.GetMaxType());
+ }
+}
+
+
+//_____________________________________________________________________
+void CbmTaskTrdHitFinder::Finish()
+{
+ std::cout << std::endl;
+ LOG(info) << "=====================================";
+ LOG(info) << GetName() << ": Finish run";
+ LOG(info) << GetName() << ": Run summary ";
+ LOG(info) << GetName() << ": Processing time : " << std::fixed << std::setprecision(3) << fProcessTime;
+ LOG(info) << GetName() << ": Nr of events : " << fNrEvents;
+ LOG(info) << GetName() << ": Nr of input digis : " << fNrDigis;
+ LOG(info) << GetName() << ": Nr of produced clusters : " << fNrClusters;
+ LOG(info) << GetName() << ": Nr of clusters / event : " << std::fixed << std::setprecision(2)
+ << (fNrEvents > 0 ? fNrClusters / (Double_t) fNrEvents : 0);
+ LOG(info) << GetName() << ": Nr of digis / cluster : " << std::fixed << std::setprecision(2)
+ << (fNrClusters > 0 ? fNrDigis / (Double_t) fNrClusters : 0);
+ LOG(info) << "=====================================";
+
+ LOG(info) << GetName() << ": Nr of events : " << fNrEvents;
+ LOG(info) << GetName() << ": Nr of input clusters : " << fNrClusters;
+ LOG(info) << GetName() << ": Nr of produced hits : " << fNrHits;
+ LOG(info) << GetName() << ": Nr of hits / event : " << std::fixed << std::setprecision(2)
+ << (fNrEvents > 0 ? fNrHits / (Double_t) fNrEvents : 0);
+ LOG(info) << GetName() << ": Nr of hits / clusters: " << std::fixed << std::setprecision(2)
+ << (fNrClusters > 0 ? fNrHits / (Double_t) fNrClusters : 0);
+ LOG(info) << "=====================================";
+ std::cout << std::endl;
+}
+
+ClassImp(CbmTaskTrdHitFinder)
diff --git a/reco/tasks/CbmTaskTrdHitFinder.h b/reco/tasks/CbmTaskTrdHitFinder.h
new file mode 100644
index 0000000000..15ca0995f4
--- /dev/null
+++ b/reco/tasks/CbmTaskTrdHitFinder.h
@@ -0,0 +1,93 @@
+/* Copyright (C) 2023 Facility for Antiproton and Ion Research in Europe, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Dominik Smith [committer], Alexandru Bercuci */
+
+#ifndef CBMTASKTRDHITFINDER_H
+#define CBMTASKTRDHITFINDER_H
+
+#include "CbmTrdCluster.h"
+#include "CbmTrdDigi.h"
+#include "CbmTrdHit.h"
+#include "FairTask.h"
+#include "trd/Hit.h"
+#include "trd/Hitfind.h"
+
+#include <map>
+#include <set>
+#include <vector>
+
+/** CbmTaskTrdHitFinder.h
+ *@author Dominik Smith <d.smith@gsi.de>
+ **
+ ** Task to build TRD and TRD2D hits from yaml configuratin.
+ **/
+class CbmTaskTrdHitFinder : public FairTask {
+
+ public:
+ /**
+ * \brief Default constructor.
+ */
+ CbmTaskTrdHitFinder();
+
+ /**
+ * \brief Default destructor.
+ */
+ ~CbmTaskTrdHitFinder();
+
+ /** Initialisation **/
+ virtual InitStatus Init();
+ virtual void SetParContainers();
+
+ /** \brief Executed task **/
+ virtual void Exec(Option_t* option);
+
+ /** Finish task **/
+ virtual void Finish();
+
+ private:
+ CbmTaskTrdHitFinder(const CbmTaskTrdHitFinder&);
+ CbmTaskTrdHitFinder& operator=(const CbmTaskTrdHitFinder&);
+
+ /**
+ * @brief Build hits from clusters for a given module
+ */
+ template<class TModule, class TCluster>
+ void BuildHits(TModule* mod, std::vector<TCluster>* clusters);
+
+ template<class TCluster>
+ void AddClusters(std::vector<TCluster>* clusters);
+ void AddHits(std::vector<cbm::algo::trd::Hit>* hits);
+
+ /** @brief Create one algo object for each RPC **/
+ bool InitAlgos();
+
+ /** Output array of CbmTrdCluster **/
+ std::vector<CbmTrdCluster>* fClusters = nullptr;
+
+ /** @brief Output array of CbmTrdHit */
+ std::vector<CbmTrdHit>* fHits = nullptr;
+
+ /** @brief Hit finding algorithm */
+ std::unique_ptr<cbm::algo::trd::Hitfind> fAlgo;
+
+ /** @brief Number of processed time slices */
+ UInt_t fNrTs = 0;
+
+ /** @brief Number of processed events (without CbmEvent corresponds to nr of exec calls) */
+ UInt_t fNrEvents = 0;
+
+ /** @brief Number of digis as input for the hit production. */
+ UInt_t fNrDigis = 0;
+
+ /** @brief Number of produced clusters. */
+ UInt_t fNrClusters = 0;
+
+ /** @brief Number of produced hits. */
+ UInt_t fNrHits = 0;
+
+ /** @brief Total processing time [RealTime]. */
+ Float_t fProcessTime = 0;
+
+ ClassDef(CbmTaskTrdHitFinder, 1);
+};
+#endif
diff --git a/reco/tasks/CbmTaskTrdHitFinderParWrite.cxx b/reco/tasks/CbmTaskTrdHitFinderParWrite.cxx
new file mode 100644
index 0000000000..9d23f6aa99
--- /dev/null
+++ b/reco/tasks/CbmTaskTrdHitFinderParWrite.cxx
@@ -0,0 +1,213 @@
+/* Copyright (C) 2010-2021 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Florian Uhlig [committer], Pascal Raisig, Alexandru Bercuci */
+
+#include "CbmTaskTrdHitFinderParWrite.h"
+
+#include "CbmTrdParAsic.h"
+#include "CbmTrdParFasp.h"
+#include "CbmTrdParModAsic.h"
+#include "CbmTrdParModDigi.h"
+#include "CbmTrdParModGeo.h"
+#include "CbmTrdParSetAsic.h"
+#include "CbmTrdParSetDigi.h"
+#include "CbmTrdParSetGeo.h"
+#include "TGeoPhysicalNode.h"
+#include "TVector3.h"
+#include "trd/Hitfind2DSetup.h"
+#include "trd/HitfindSetup.h"
+#include "yaml/Yaml.h"
+
+#include <FairRootManager.h>
+#include <FairRunAna.h>
+#include <FairRuntimeDb.h>
+#include <Logger.h>
+
+#include <iomanip>
+#include <iostream>
+#include <vector>
+
+using std::fixed;
+using std::left;
+using std::right;
+using std::setprecision;
+using std::setw;
+using std::stringstream;
+
+//_____________________________________________________________________
+CbmTaskTrdHitFinderParWrite::CbmTaskTrdHitFinderParWrite() : FairTask("TrdClusterFinder", 1) {}
+
+//_____________________________________________________________________
+void CbmTaskTrdHitFinderParWrite::SetParContainers()
+{
+ FairRuntimeDb* rtdb = FairRunAna::Instance()->GetRuntimeDb();
+ fAsicPar = static_cast<CbmTrdParSetAsic*>(rtdb->getContainer("CbmTrdParSetAsic"));
+ fDigiPar = static_cast<CbmTrdParSetDigi*>(rtdb->getContainer("CbmTrdParSetDigi"));
+ fGeoPar = static_cast<CbmTrdParSetGeo*>(rtdb->getContainer("CbmTrdParSetGeo"));
+}
+
+//_____________________________________________________________________
+InitStatus CbmTaskTrdHitFinderParWrite::Init()
+{
+ // Get the full geometry information of the detector gas layers and store
+ // them with the CbmTrdModuleRec. This information can then be used for
+ // transformation calculations
+ if (!fGeoPar->LoadAlignVolumes()) {
+ LOG(error) << GetName() << "::Init: GEO info for modules unavailable !";
+ return kFATAL;
+ }
+ if (fDigiPar->GetNrOfModules() != fGeoPar->GetNrOfModules()) {
+ LOG(fatal) << "CbmTaskTrdHitFinderParWrite::Init() - Geometry and parameter files"
+ << " have different number of modules.";
+ }
+
+ // Sets of module IDs
+ std::vector<uint16_t> trdModules;
+ std::vector<uint16_t> trd2DModules;
+
+ // Store IDs of 1D and 2D modules
+ for (auto entry : fDigiPar->GetModuleMap()) {
+ const auto moduleId = entry.first;
+
+ // Get ASIC parameters for this module
+ const CbmTrdParModAsic* asicPar = static_cast<const CbmTrdParModAsic*>(fAsicPar->GetModulePar(moduleId));
+ if (!asicPar) continue;
+
+ if (asicPar->GetAsicType() == CbmTrdDigi::eCbmTrdAsicType::kFASP) {
+ trd2DModules.push_back(entry.first);
+ }
+ else {
+ trdModules.push_back(entry.first);
+ }
+ }
+
+ // Create setup files for each type
+ cbm::algo::trd::HitfindSetup setup;
+ cbm::algo::trd::Hitfind2DSetup setup2D;
+ setup.modules.resize(trdModules.size());
+ setup2D.modules.resize(trd2DModules.size());
+
+ // Fill 2D setup files
+ for (size_t mod; mod < trd2DModules.size(); mod++) {
+ const uint16_t moduleId = trd2DModules[mod];
+
+ // Get Geometry parameters for module
+ CbmTrdParModGeo* pGeo = static_cast<CbmTrdParModGeo*>(fGeoPar->GetModulePar(moduleId));
+
+ // Get ASIC parameters for this module
+ const CbmTrdParModAsic* asicPar = static_cast<const CbmTrdParModAsic*>(fAsicPar->GetModulePar(moduleId));
+
+ if (!asicPar) continue;
+ const CbmTrdParModDigi* digiPar = static_cast<const CbmTrdParModDigi*>(fDigiPar->GetModulePar(moduleId));
+
+ cbm::algo::trd::Hitfind2DSetup::Mod& module = setup2D.modules[mod];
+ module.padSizeX = digiPar->GetPadSizeX(0);
+ module.padSizeY = digiPar->GetPadSizeY(0);
+ module.address = moduleId;
+ module.rowPar.resize(digiPar->GetNofRows());
+
+ const double* tra_ptr = pGeo->GetNode()->GetMatrix()->GetTranslation();
+ const double* rot_ptr = pGeo->GetNode()->GetMatrix()->GetRotationMatrix();
+ memcpy(module.translation.data(), tra_ptr, 3 * sizeof(double));
+ memcpy(module.rotation.data(), rot_ptr, 9 * sizeof(double));
+
+ for (int row = 0; row < digiPar->GetNofRows(); row++) {
+
+ cbm::algo::trd::Hitfind2DSetup::Row& rowPar = module.rowPar[row];
+ rowPar.padPar.resize(digiPar->GetNofColumns());
+
+ int rowInSector = 0;
+ const int sector = digiPar->GetSectorRow(row, rowInSector); // Second is ref. TO DO: Return std::pair instead.
+
+ for (int col = 0; col < digiPar->GetNofColumns(); col++) {
+ cbm::algo::trd::Hitfind2DSetup::Pad& padPar = rowPar.padPar[col];
+ TVector3 pos, posErr;
+ digiPar->GetPadPosition(sector, col, rowInSector, pos, posErr);
+ const std::array<double, 3> pos_ptr = {pos[0], pos[1], pos[2]};
+ const std::array<double, 3> posErr_ptr = {posErr[0], posErr[1], posErr[2]};
+ memcpy(padPar.position.data(), pos_ptr.data(), 3 * sizeof(double));
+ memcpy(padPar.positionError.data(), posErr_ptr.data(), 3 * sizeof(double));
+
+ const size_t chan = row * digiPar->GetNofColumns() + col;
+ const CbmTrdParFaspChannel *daqFaspChT(nullptr), *daqFaspChR(nullptr);
+ asicPar->GetFaspChannelPar(chan, daqFaspChT, daqFaspChR);
+ if (!daqFaspChR) {
+ padPar.chRMasked = false; // TODO implement case for not installed
+ }
+ else {
+ padPar.chRMasked = daqFaspChR->IsMasked();
+ }
+ if (!daqFaspChT) {
+ padPar.chTMasked = false; // TODO implement case for not installed
+ }
+ else {
+ padPar.chTMasked = daqFaspChT->IsMasked();
+ }
+ }
+ }
+ } // for (size_t mod; mod < trd2DModules.size(); mod++) {
+
+
+ // Fill 1D setup files
+ for (size_t mod; mod < trdModules.size(); mod++) {
+ const uint16_t moduleId = trdModules[mod];
+
+ // Get Geometry parameters for module
+ CbmTrdParModGeo* pGeo = static_cast<CbmTrdParModGeo*>(fGeoPar->GetModulePar(moduleId));
+
+ // Get ASIC parameters for this module
+ const CbmTrdParModAsic* asicPar = static_cast<const CbmTrdParModAsic*>(fAsicPar->GetModulePar(moduleId));
+
+ if (!asicPar) continue;
+ const CbmTrdParModDigi* digiPar = static_cast<const CbmTrdParModDigi*>(fDigiPar->GetModulePar(moduleId));
+
+ cbm::algo::trd::HitfindSetup::Mod& module = setup.modules[mod];
+ module.padSizeX = digiPar->GetPadSizeX(0);
+ module.padSizeY = digiPar->GetPadSizeY(0);
+ module.padSizeErrX = digiPar->GetPadSizeX(1);
+ module.padSizeErrY = digiPar->GetPadSizeY(1);
+ module.orientation = digiPar->GetOrientation();
+ module.address = moduleId;
+ module.rowPar.resize(digiPar->GetNofRows());
+
+ const double* tra_ptr = pGeo->GetNode()->GetMatrix()->GetTranslation();
+ const double* rot_ptr = pGeo->GetNode()->GetMatrix()->GetRotationMatrix();
+ memcpy(module.translation.data(), tra_ptr, 3 * sizeof(double));
+ memcpy(module.rotation.data(), rot_ptr, 9 * sizeof(double));
+
+ for (int row = 0; row < digiPar->GetNofRows(); row++) {
+
+ cbm::algo::trd::HitfindSetup::Row& rowPar = module.rowPar[row];
+ rowPar.padPar.resize(digiPar->GetNofColumns());
+
+ int rowInSector = 0;
+ const int sector = digiPar->GetSectorRow(row, rowInSector); // Second is ref. TO DO: Return std::pair instead.
+
+ for (int col = 0; col < digiPar->GetNofColumns(); col++) {
+ cbm::algo::trd::HitfindSetup::Pad& padPar = rowPar.padPar[col];
+ TVector3 pos, posErr;
+ digiPar->GetPadPosition(sector, col, rowInSector, pos, posErr);
+ const std::array<double, 3> pos_ptr = {pos[0], pos[1], pos[2]};
+ const std::array<double, 3> posErr_ptr = {posErr[0], posErr[1], posErr[2]};
+ memcpy(padPar.position.data(), pos_ptr.data(), 3 * sizeof(double));
+ memcpy(padPar.positionError.data(), posErr_ptr.data(), 3 * sizeof(double));
+ }
+ }
+ } // for (size_t mod; mod < trdModules.size(); mod++) {
+
+
+ /* Write Yaml files */
+
+ cbm::algo::yaml::Dump dump;
+ std::ofstream fout("TrdHitfinderPar.yaml");
+ fout << dump(setup);
+ fout.close();
+
+ std::ofstream f2Dout("TrdHitfinder2DPar.yaml");
+ f2Dout << dump(setup2D);
+ f2Dout.close();
+
+ return kSUCCESS;
+}
+
+ClassImp(CbmTaskTrdHitFinderParWrite)
diff --git a/reco/tasks/CbmTaskTrdHitFinderParWrite.h b/reco/tasks/CbmTaskTrdHitFinderParWrite.h
new file mode 100644
index 0000000000..24ff1fbce5
--- /dev/null
+++ b/reco/tasks/CbmTaskTrdHitFinderParWrite.h
@@ -0,0 +1,63 @@
+/* Copyright (C) 2024 Facility for Antiproton and Ion Research in Europe, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Dominik Smith [committer] */
+
+#ifndef CBMTRDHITFINDERPARWRITE_H
+#define CBMTRDHITFINDERPARWRITE_H
+
+#include "FairTask.h"
+#include "trd/Clusterizer.h"
+#include "trd/Clusterizer2D.h"
+#include "trd/Hit.h"
+#include "trd/HitFinder.h"
+#include "trd/HitFinder2D.h"
+
+#include <map>
+#include <set>
+#include <vector>
+
+class CbmTrdParSetAsic;
+class CbmTrdParSetGeo;
+class CbmTrdParSetDigi;
+
+/** CbmTaskTrdHitFinderParWrite.h
+ *@author Dominik Smith <d.smith@gsi.de>
+ **
+ ** Task to create YAML files for TRD hitfinders
+ **
+ **/
+class CbmTaskTrdHitFinderParWrite : public FairTask {
+
+ public:
+ /**
+ * \brief Default constructor.
+ */
+ CbmTaskTrdHitFinderParWrite();
+
+ /**
+ * \brief Default destructor.
+ */
+ ~CbmTaskTrdHitFinderParWrite(){};
+
+ /** Initialisation **/
+ virtual InitStatus Init();
+ virtual void SetParContainers();
+
+ /** \brief Executed task **/
+ virtual void Exec(Option_t* option){};
+
+ /** Finish task **/
+ virtual void Finish(){};
+
+ private:
+ CbmTaskTrdHitFinderParWrite(const CbmTaskTrdHitFinderParWrite&);
+ CbmTaskTrdHitFinderParWrite& operator=(const CbmTaskTrdHitFinderParWrite&);
+
+ //==================================================================
+ CbmTrdParSetAsic* fAsicPar = nullptr; ///< parameter list for ASIC characterization
+ CbmTrdParSetDigi* fDigiPar = nullptr; ///< parameter list for read-out geometry
+ CbmTrdParSetGeo* fGeoPar = nullptr; ///< parameter list for modules geometry
+
+ ClassDef(CbmTaskTrdHitFinderParWrite, 1);
+};
+#endif
--
GitLab