diff --git a/algo/CMakeLists.txt b/algo/CMakeLists.txt
index b84c480fa2928b21e0c7c67ebe5879cd43dee29f..0390d9e536445b0aaeb632fa76fef9ee100866ef 100644
--- a/algo/CMakeLists.txt
+++ b/algo/CMakeLists.txt
@@ -131,8 +131,11 @@ set(SRCS
detectors/trd/DigiRec.cxx
detectors/trd/Hit.cxx
detectors/trd/Hitfind.cxx
+ detectors/trd/HitFactory2D.cxx
detectors/trd/HitFinder.cxx
detectors/trd/HitFinder2D.cxx
+ detectors/trd/HitMerger.cxx
+ detectors/trd/HitMerger2D.cxx
detectors/trd/ReadoutConfig.cxx
detectors/trd/TrackingInterface.cxx
detectors/trd/Unpack.cxx
diff --git a/algo/detectors/trd/Cluster2D.cxx b/algo/detectors/trd/Cluster2D.cxx
index 35253c2bfabb95a62fc7630d29d5377cd01e1e29..b5e409c0041b20328231f00fffb0515df3dee573 100644
--- a/algo/detectors/trd/Cluster2D.cxx
+++ b/algo/detectors/trd/Cluster2D.cxx
@@ -5,6 +5,7 @@
#include "Cluster2D.h"
#include "CbmTrdDigi.h"
+#include "log.hpp"
#include <cmath>
@@ -14,7 +15,8 @@ namespace cbm::algo::trd
{
//____________________________________________________________________
Cluster2D::Cluster2D(const Cluster2D& ref)
- : fDigis(ref.fDigis)
+ : fRecDigis(ref.fRecDigis)
+ , fDigis(ref.fDigis)
, fDigiIndices(ref.fDigiIndices)
, fAddress(ref.fAddress)
, fNCols(ref.fNCols)
@@ -146,4 +148,39 @@ namespace cbm::algo::trd
return false;
}
+ bool Cluster2D::Finalize(size_t numCols)
+ {
+ /** 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 = fDigis.begin(); i != fDigis.end(); i++) {
+ if ((*i) == nullptr) continue;
+ int colR(-1);
+ const CbmTrdDigi* dgT = (*i);
+ int colT = dgT->GetAddressChannel() % numCols;
+
+ // 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 not in increasing order !";
+ return false;
+ }
+ last_col = colT;
+ const CbmTrdDigi* dgR = nullptr;
+
+ auto j = std::next(i);
+ if (j != fDigis.end()) {
+ dgR = (*j);
+ colR = dgR->GetAddressChannel() % numCols;
+ }
+ if (colR == colT && dgR != nullptr) {
+ fRecDigis.emplace_back(*dgT, *dgR);
+ i++;
+ }
+ else {
+ fRecDigis.emplace_back(*dgT);
+ }
+ }
+ return true;
+ }
} // namespace cbm::algo::trd
diff --git a/algo/detectors/trd/Cluster2D.h b/algo/detectors/trd/Cluster2D.h
index bd3f4cf29331d6003d39e20049a92cb556ddd82e..40cf184475fcd5447d783710872fbd8eaafedd4c 100644
--- a/algo/detectors/trd/Cluster2D.h
+++ b/algo/detectors/trd/Cluster2D.h
@@ -4,6 +4,7 @@
#pragma once
+#include "DigiRec.h"
#include "compat/RTypes.h"
#include <cstdint>
@@ -98,6 +99,11 @@ namespace cbm::algo::trd
*/
std::vector<const CbmTrdDigi*>& GetDigis() { return fDigis; }
+ /**
+ * \brief Get array of calibrated digis
+ * \return Array of calibrated digis
+ */
+ const std::vector<DigiRec>& GetRecDigis() const { return fRecDigis; }
/**
* \brief Remove all digis.
@@ -106,6 +112,7 @@ namespace cbm::algo::trd
{
fDigiIndices.clear();
fDigis.clear();
+ fRecDigis.clear();
}
/** Accessors **/
@@ -148,6 +155,12 @@ namespace cbm::algo::trd
*/
bool Merge(Cluster2D* second);
+ /** \brief Fill array of calibrated digis
+ * \param[in] number of columns
+ * \return success or fail
+ */
+ bool Finalize(const size_t numCols);
+
/** Setters **/
void SetNCols(uint16_t ncols) { fNCols = ncols; }
void SetNRows(uint16_t nrows)
@@ -159,6 +172,7 @@ namespace cbm::algo::trd
void SetStop(bool set = true) { set ? SETBIT(fNRows, kStop) : CLRBIT(fNRows, kStop); }
private:
+ std::vector<DigiRec> fRecDigis; ///< Array of calibrated digis
std::vector<const CbmTrdDigi*> fDigis; ///< Array of digi pointers
std::vector<int32_t> fDigiIndices; ///< Array of digi indices
diff --git a/algo/detectors/trd/Clusterizer.cxx b/algo/detectors/trd/Clusterizer.cxx
index 962ed2fdd19ec5d2a705ec112f76bd80285fdb8d..dff15aa23eab38e215402fcf86449d1e28acd1de 100644
--- a/algo/detectors/trd/Clusterizer.cxx
+++ b/algo/detectors/trd/Clusterizer.cxx
@@ -135,8 +135,6 @@ namespace cbm::algo::trd
addClusters(cluster, &clustersOut);
} //! for (auto mainit = inputData.begin(); mainit != inputData.end(); mainit++)
- /// Row merging could be here in principle.
-
return clustersOut;
}
diff --git a/algo/detectors/trd/Clusterizer2D.cxx b/algo/detectors/trd/Clusterizer2D.cxx
index 3782ee8681000463471bdc12d479b677478c027d..f0028509a2a87859ea68cf99b9da8a14edd8ef74 100644
--- a/algo/detectors/trd/Clusterizer2D.cxx
+++ b/algo/detectors/trd/Clusterizer2D.cxx
@@ -109,10 +109,12 @@ namespace cbm::algo::trd
continue;
}
clustersOut.emplace_back(*pcluster);
+ clustersOut.back().Finalize(fParams.rowPar[0].padPar.size());
delete pcluster;
pcluster = nullptr;
}
}
+
return clustersOut;
}
diff --git a/algo/detectors/trd/DigiRec.h b/algo/detectors/trd/DigiRec.h
index a6ef0cfe5f4217a7bc0a72f224a8528d2c50a245..135171e26d8d84627597c05f053e25cf8567a6ef 100644
--- a/algo/detectors/trd/DigiRec.h
+++ b/algo/detectors/trd/DigiRec.h
@@ -22,8 +22,8 @@
namespace cbm::algo::trd
{
-
class DigiRec : public CbmTrdDigi {
+ friend class Cluster2D;
friend class HitFinder2D;
public:
@@ -31,6 +31,9 @@ namespace cbm::algo::trd
DigiRec();
/** \brief Wrap CbmTrdDigi constructor*/
DigiRec(const CbmTrdDigi& d, double* g = NULL, double* t = NULL);
+ /** \brief Constructor and merger*/
+ DigiRec(const CbmTrdDigi& dt, const CbmTrdDigi& dr, double* g = NULL, double* t = NULL);
+
virtual ~DigiRec() { ; }
/** \brief Return calibrated tilt signal
@@ -59,10 +62,6 @@ namespace cbm::algo::trd
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
diff --git a/algo/detectors/trd/Hit.cxx b/algo/detectors/trd/Hit.cxx
index c82d5f8b1dbcd8da91614b1526ec97b7ea39953d..d452f365784c616ba37dc1dfad158175b1b7e6d3 100644
--- a/algo/detectors/trd/Hit.cxx
+++ b/algo/detectors/trd/Hit.cxx
@@ -40,27 +40,6 @@ namespace cbm::algo::trd
{
}
- 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(X(), Y(), Z()); }
void Hit::PositionError(ROOT::Math::XYZVector& dpos) const { dpos.SetXYZ(Dx(), Dy(), Dz()); }
diff --git a/algo/detectors/trd/Hit.h b/algo/detectors/trd/Hit.h
index 844dfde9fd5dda4ff8317681d74e98bd2ea48673..c4d78b06e8eea0e10fc04d885eef45ade60a2c59 100644
--- a/algo/detectors/trd/Hit.h
+++ b/algo/detectors/trd/Hit.h
@@ -48,8 +48,6 @@ namespace cbm::algo::trd
* \brief Default constructor.
*/
Hit();
- // Hit(const Hit&);
- Hit& operator=(const Hit&);
/**
* \brief Standard constructor.
diff --git a/algo/detectors/trd/HitFactory2D.cxx b/algo/detectors/trd/HitFactory2D.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..1d94a05df0703092c964ff2db24683819b4e4856
--- /dev/null
+++ b/algo/detectors/trd/HitFactory2D.cxx
@@ -0,0 +1,328 @@
+/* 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
+{
+
+ //_______________________________________________________________________________
+ std::pair<double, double> HitFactory2D::CorrectPosition(double dx, double dy, const double xcorr,
+ const double padSizeX, const double padSizeY)
+ {
+ dx -= xcorr;
+ RecenterXoffset(dx);
+ dy = dx - dy;
+ RecenterYoffset(dy);
+
+ const double ycorr = GetYcorr(dy) / padSizeY;
+ dy += ycorr;
+ RecenterYoffset(dy);
+ dx *= padSizeX;
+ dy *= padSizeY;
+ return std::make_pair(dx, dy);
+ }
+
+
+ //_______________________________________________________________________________
+ 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/HitFactory2D.h b/algo/detectors/trd/HitFactory2D.h
new file mode 100644
index 0000000000000000000000000000000000000000..c2fa53a7dce9b67c343580516f6d5b2b15449d12
--- /dev/null
+++ b/algo/detectors/trd/HitFactory2D.h
@@ -0,0 +1,137 @@
+/* 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>
+
+#define NBINSCORRX 50 //! no of bins in the discretized correction LUT
+#define NBINSCORRY 4 //! no of bins in the parametrization correction
+
+using std::vector;
+
+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]
+ */
+
+ std::pair<double, double> CorrectPosition(double dx, double dy, const double xcorr, const double padSizeX,
+ const double padSizeY);
+
+ 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;
+ }
+ };
+
+} // namespace cbm::algo::trd
diff --git a/algo/detectors/trd/HitFinder.cxx b/algo/detectors/trd/HitFinder.cxx
index 84a57e6287866e0fd16768ad52bc40665fd6dbc4..358cfc68604af832ef9d7c4728ec5cf72735a475 100644
--- a/algo/detectors/trd/HitFinder.cxx
+++ b/algo/detectors/trd/HitFinder.cxx
@@ -15,18 +15,18 @@ namespace cbm::algo::trd
HitFinder::HitFinder(HitFinderModPar params) : fParams(params) {}
//_______________________________________________________________________________
- std::vector<Hit> HitFinder::operator()(std::vector<Cluster>* clusters)
+ std::vector<trd::HitFinder::resultType> HitFinder::operator()(std::vector<Cluster>* clusters)
{
const int nclusters = clusters->size();
- std::vector<Hit> hits;
+ std::vector<resultType> hitData;
for (int icluster = 0; icluster < nclusters; icluster++) {
Cluster* cluster = &clusters->at(icluster);
auto hit = MakeHit(icluster, cluster, &cluster->GetDigis(), nclusters);
if (hit.Address() == -1) continue;
- hits.push_back(hit);
+ hitData.emplace_back(hit, std::vector<DigiRec>()); //DigiRec currently not used in 1D but in 2D
}
- return hits;
+ return hitData;
}
//_______________________________________________________________________________
@@ -178,10 +178,7 @@ namespace cbm::algo::trd
int colMax = 0;
int rowMax = 0;
- for (int i = 0; i < cluster->GetNofDigis(); ++i) {
- // const CbmTrdDigi* digi = CbmDigiManager::Instance()->Get<CbmTrdDigi>(cluster->GetDigi(i));
- const CbmTrdDigi* digi = (cluster->GetDigis())[i];
-
+ for (auto& digi : cluster->GetDigis()) {
int digiCol;
int digiRow = GetPadRowCol(digi->GetAddressChannel(), digiCol);
@@ -197,10 +194,7 @@ namespace cbm::algo::trd
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];
-
+ for (auto& digi : cluster->GetDigis()) {
int digiCol;
int digiRow = GetPadRowCol(digi->GetAddressChannel(), digiCol);
diff --git a/algo/detectors/trd/HitFinder.h b/algo/detectors/trd/HitFinder.h
index 3f4475fcc92901c5e9ff8dfa9ae84c2616640b7c..db5503c2db01fa879f662501ec705dff76bcd43b 100644
--- a/algo/detectors/trd/HitFinder.h
+++ b/algo/detectors/trd/HitFinder.h
@@ -5,6 +5,7 @@
#pragma once
#include "Cluster.h"
+#include "DigiRec.h"
#include "Hit.h"
#include "HitFinderPars.h"
#include "Math/Rotation3D.h"
@@ -22,6 +23,8 @@ namespace cbm::algo::trd
**/
class HitFinder {
public:
+ typedef std::pair<Hit, std::vector<DigiRec>> resultType;
+
HitFinder(){};
/**
* \brief Constructor with placement
@@ -30,7 +33,7 @@ namespace cbm::algo::trd
virtual ~HitFinder(){};
/* \brief Steering routine for building hits */
- std::vector<Hit> operator()(std::vector<Cluster>* clusters);
+ std::vector<resultType> operator()(std::vector<Cluster>* clusters);
double GetSpaceResolution(double val = 3.0);
bool IsClusterComplete(const Cluster* cluster);
diff --git a/algo/detectors/trd/HitFinder2D.cxx b/algo/detectors/trd/HitFinder2D.cxx
index 0adb74bb54cff675218a961e35d673d2c7e58dd6..37fa9bafe371a88e2b2ea26c3ffc04982ec234d0 100644
--- a/algo/detectors/trd/HitFinder2D.cxx
+++ b/algo/detectors/trd/HitFinder2D.cxx
@@ -19,263 +19,48 @@ namespace cbm::algo::trd
double HitFinder2D::fgDT[] = {4.181e-6, 1586, 24};
//_______________________________________________________________________________
- HitFinder2D::HitFinder2D(HitFinder2DModPar params) : fParams(params), fHit(params.rowPar[0].padPar.size()) {}
+ HitFinder2D::HitFinder2D(HitFinder2DModPar params) : fParams(params) {}
//_______________________________________________________________________________
- std::vector<Hit> HitFinder2D::operator()(std::vector<Cluster2D>* clusters)
+ std::vector<trd::HitFinder2D::resultType> HitFinder2D::operator()(std::vector<Cluster2D>* clusters)
{
const int nclusters = clusters->size();
- std::vector<Hit> hits;
+ std::vector<resultType> hitData;
for (int icluster = 0; icluster < nclusters; icluster++) {
- Cluster2D* cluster = &clusters->at(icluster);
- auto hit = MakeHit(icluster, cluster, &cluster->GetDigis(), nclusters);
+ auto& cluster = clusters->at(icluster);
+ auto hit = MakeHit(icluster, &cluster);
if (hit.Address() == -1) continue;
- hits.push_back(hit);
+ hitData.emplace_back(hit, cluster.GetRecDigis());
}
-
- 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->Time() < 4000 - h0->Time()) continue; // TODO check with preoper time calibration
- // skip next hits for being too far (10us) in the future
- if (h1->Time() > 10000 + h0->Time()) break;
-
- // basic check on Row
- if (std::abs(h1->Y() - h0->Y()) > Dy) continue;
-
- // basic check on Col
- if (std::abs(h1->X() - h0->X()) > 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;
+ return hitData;
}
//_______________________________________________________________________________
- int HitFinder2D::CheckMerge(int cid, int cjd)
+ Hit HitFinder2D::MakeHit(int ic, const Cluster2D* cluster)
{
- /** 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];
+ Hit hit;
+ if (cluster->GetRecDigis().size() == 0) return hit;
- 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;
+ auto [hitFlag, hitF] = ProjectDigis(cluster);
+ if (!hitFlag) return hit;
- for (int ia(0); ia < 3; ia++) {
- if (mM < mM_l[ia]) return (rowMax ? 1 : -1) * (3 - ia);
- }
- return 0;
+ hit.SetAddress(fParams.address);
+ hit.SetRefId(ic);
+ BuildHit(&hit, hitF);
+ return hit;
}
//_______________________________________________________________________________
- 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)
+ bool HitFinder2D::BuildHit(Hit* h, HitFactory2D& hitF)
{
- const int n0 = fHit.fSignal.size() - 2;
- auto [dx, dy] = fHit.GetDxDy(n0);
+ const int n0 = hitF.fSignal.size() - 2;
+ auto [dx, dy] = hitF.GetDxDy(n0);
// get position correction
- const double xcorr = fHit.GetXcorr(dx, fHit.GetHitClass()) / fParams.padSizeX;
- std::tie(dx, dy) = CorrectPosition(dx, dy, xcorr);
+ const double xcorr = hitF.GetXcorr(dx, hitF.GetHitClass()) / fParams.padSizeX;
+ std::tie(dx, dy) = hitF.CorrectPosition(dx, dy, xcorr, fParams.padSizeX, fParams.padSizeY);
// get anode wire offset
for (int ia = 0; ia <= NANODE; ia++) {
@@ -295,7 +80,7 @@ namespace cbm::algo::trd
// COMPUTE TIME
double t_avg = 0.;
for (int idx = 1; idx <= n0; idx++) {
- HitFactory2D::signal& sig = fHit.fSignal[idx];
+ HitFactory2D::signal& sig = hitF.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;
@@ -306,121 +91,51 @@ namespace cbm::algo::trd
// COMPUTE ENERGY (conversion to GeV)
const double e_avg = 1.e-9
- * std::accumulate(fHit.fSignal.begin(), fHit.fSignal.end(), 0,
+ * std::accumulate(hitF.fSignal.begin(), hitF.fSignal.end(), 0,
[](int sum, const auto& sig) { return sum + sig.s; });
- CalibrateHit(h, dx, dy, edx, edy, edt, time, tdrift, e_avg);
+ CalibrateHit(h, dx, dy, edx, edy, edt, time, tdrift, e_avg, hitF);
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 double edt, const double time, const double tdrift, const double eloss,
+ const HitFactory2D& hitF)
{
- const ROOT::Math::XYZVector& local_pad = fParams.rowPar[fHit.vrM].padPar[fHit.vcM].pos;
+ const ROOT::Math::XYZVector& local_pad = fParams.rowPar[hitF.vrM].padPar[hitF.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->SetTime(fClk * (hitF.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();
+ h->SetMaxType(hitF.IsMaxTilt());
+ h->SetOverFlow(hitF.HasOvf());
}
//_______________________________________________________________________________
- int HitFinder2D::ProjectDigis(int cid, int cjd)
+ std::pair<int, HitFactory2D> HitFinder2D::ProjectDigis(const Cluster2D* cluster)
{
/** 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;
+ if (cluster->GetRecDigis().empty()) {
+ L_(debug) << "TrdModuleRec2D::ProjectDigis : Requested cluster not found.";
+ return std::make_pair(0, HitFactory2D(0));
}
+ HitFactory2D hitF(fParams.rowPar[0].padPar.size());
- std::vector<HitFactory2D::signal>& hitSig = fHit.fSignal;
- fHit.reset();
+ std::vector<HitFactory2D::signal>& hitSig = hitF.fSignal;
+ hitF.reset();
bool on(0); // flag signal transmition
int n0(0), nsr(0), // no of signals in the cluster (all/rect)
@@ -434,28 +149,15 @@ namespace cbm::algo::trd
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();
- }
+ int j(0), col(-1), col0(0);
- for (auto i = fDigis[cid].begin(); i != fDigis[cid].end(); i++, j++) {
- const DigiRec* dg = (*i);
- const DigiRec* dg0 = nullptr;
- const DigiRec* dg1 = nullptr;
+ for (auto i = cluster->GetRecDigis().begin(); i != cluster->GetRecDigis().end(); i++, j++) {
+ const DigiRec* dg = &(*i);
// initialize
if (col < 0) {
- fHit.vrM = GetPadRowCol(dg->GetAddressChannel(), col);
- fHit.vt0 = dg->GetTimeDAQ(); // set arbitrary t0 to avoid double digis loop
+ hitF.vrM = GetPadRowCol(dg->GetAddressChannel(), col);
+ hitF.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
@@ -466,44 +168,10 @@ namespace cbm::algo::trd
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
+ char ddt = dg->GetTiltTime() - hitF.vt0; // signal time offset wrt prompt
if (ddt < dt0) dt0 = ddt;
if (abs(t) > 0) {
if (nt > 1) t *= 0.5;
@@ -514,9 +182,9 @@ namespace cbm::algo::trd
}
if (t > max) {
max = t;
- fHit.vcM = j;
- fHit.SetMaxTilt(1);
- fHit.viM = hitSig.size();
+ hitF.vcM = j;
+ hitF.SetMaxTilt(1);
+ hitF.viM = hitSig.size();
}
}
else
@@ -526,7 +194,7 @@ namespace cbm::algo::trd
xc += 0.5;
// process rect signal/time
- ddt = dg->GetRectTime() - fHit.vt0; // signal time offset wrt prompt
+ ddt = dg->GetRectTime() - hitF.vt0; // signal time offset wrt prompt
if (ddt < dt0) dt0 = ddt;
if (abs(r) > 0) {
nsr++;
@@ -538,9 +206,9 @@ namespace cbm::algo::trd
}
if (r > max) {
max = r;
- fHit.vcM = j;
- fHit.SetMaxTilt(0);
- fHit.viM = hitSig.size();
+ hitF.vcM = j;
+ hitF.SetMaxTilt(0);
+ hitF.viM = hitSig.size();
}
}
else
@@ -558,7 +226,7 @@ namespace cbm::algo::trd
xc = hitSig[0].x;
char ddt = hitSig[0].t;
hitSig.emplace(hitSig.begin(), 0., 300., ddt, xc - 0.5, 0.);
- fHit.viM++;
+ hitF.viM++;
}
int n(hitSig.size() - 1);
if (std::abs(hitSig[n].s) > 1.e-3) {
@@ -569,366 +237,45 @@ namespace cbm::algo::trd
n0 = hitSig.size() - 2;
// compute cluster asymmetry
- int nR = n0 + 1 - fHit.viM;
- if (nR == fHit.viM) {
- fHit.SetSymmHit(1);
+ int nR = n0 + 1 - hitF.viM;
+ if (nR == hitF.viM) {
+ hitF.SetSymmHit(1);
if (hitSig.size() % 2) {
double LS(0.), RS(0.);
- for (UChar_t idx(0); idx < fHit.viM; idx++)
+ for (UChar_t idx(0); idx < hitF.viM; idx++)
LS += hitSig[idx].s;
- for (uint idx(fHit.viM + 1); idx < hitSig.size(); idx++)
+ for (uint idx(hitF.viM + 1); idx < hitSig.size(); idx++)
RS += hitSig[idx].s;
- fHit.SetLeftSgn(LS < RS ? 0 : 1);
+ hitF.SetLeftSgn(LS < RS ? 0 : 1);
}
}
else {
- fHit.SetSymmHit(0);
- if (fHit.viM < nR)
- fHit.SetLeftHit(0);
- else if (fHit.viM > nR)
- fHit.SetLeftHit(1);
+ hitF.SetSymmHit(0);
+ if (hitF.viM < nR)
+ hitF.SetLeftHit(0);
+ else if (hitF.viM > nR)
+ hitF.SetLeftHit(1);
}
// recenter time and space profile
- fHit.vt0 += dt0;
+ hitF.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
- }
+ sig.x -= hitF.vcM;
}
- else
- fHit.SetBiasY(0);
+ hitF.vcM += col;
+
+ hitF.SetBiasX(0);
+ hitF.SetBiasY(0);
- if (ovf < 0) fHit.SetOvf();
- return ovf * (hitSig.size() - 2);
+ if (ovf < 0) hitF.SetOvf();
+ return std::make_pair(ovf * (hitSig.size() - 2), hitF);
}
+
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
index 5398b58d8c17e3c8c48d199532212e5ee4c0b974..973f2f8ed97fa5bbdd5d886ec75ad600e22459c1 100644
--- a/algo/detectors/trd/HitFinder2D.h
+++ b/algo/detectors/trd/HitFinder2D.h
@@ -8,6 +8,7 @@
#include "Cluster2D.h"
#include "DigiRec.h"
#include "Hit.h"
+#include "HitFactory2D.h"
#include "HitFinder2DPars.h"
#include "Math/Rotation3D.h"
#include "Math/Vector3Dfwd.h"
@@ -17,136 +18,13 @@
#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>
@@ -163,55 +41,35 @@ namespace cbm::algo::trd
class HitFinder2D {
public:
- typedef std::tuple<uint16_t, uint16_t, int, int, size_t> inputType; //Tuple: chT, chR, tm, row, id
+ typedef std::pair<Hit, std::vector<DigiRec>> resultType;
/** \brief Default constructor.*/
- HitFinder2D() : fHit(0){};
+ HitFinder2D(){};
+
/** \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
- */
+ std::vector<resultType> operator()(std::vector<Cluster2D>* clusters);
- int LoadDigis(vector<const CbmTrdDigi*>* din, int cid);
- int ProjectDigis(int cid, int cjd = -1);
+ std::pair<int, HitFactory2D> ProjectDigis(const Cluster2D* cluster);
/** \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]
- */
+ bool BuildHit(Hit* h, HitFactory2D& hitF);
/** \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);
+ Hit MakeHit(int cId, const Cluster2D* cluster);
const float fClk = CbmTrdDigi::Clk(CbmTrdDigi::eCbmTrdAsicType::kFASP);
@@ -223,18 +81,11 @@ namespace cbm::algo::trd
*/
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
+ const double time, const double tdrift, const double eloss, const HitFactory2D& hitF);
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
};
diff --git a/algo/detectors/trd/HitMerger.cxx b/algo/detectors/trd/HitMerger.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..9adddadc589d476fe0a93eb855f5cbd04ede75ec
--- /dev/null
+++ b/algo/detectors/trd/HitMerger.cxx
@@ -0,0 +1,21 @@
+/* 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 "HitMerger.h"
+
+namespace cbm::algo::trd
+{
+
+ //_______________________________________________________________________________
+ HitMerger::HitMerger(HitFinderModPar params) : fParams(params) {}
+
+ HitMerger::outputType HitMerger::operator()(std::vector<inputType>& hitsRow1, std::vector<inputType>& hitsRow2)
+ {
+ //// TO DO: Implement something here!
+
+ return std::make_pair(std::move(hitsRow1), std::move(hitsRow2));
+ }
+
+
+} // namespace cbm::algo::trd
diff --git a/algo/detectors/trd/HitMerger.h b/algo/detectors/trd/HitMerger.h
new file mode 100644
index 0000000000000000000000000000000000000000..4d01c20e325eb4ede4c26b43d02f3dfeca4df097
--- /dev/null
+++ b/algo/detectors/trd/HitMerger.h
@@ -0,0 +1,43 @@
+/* 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 "DigiRec.h"
+#include "Hit.h"
+#include "HitFinderPars.h"
+
+#include <tuple>
+#include <vector>
+
+namespace cbm::algo::trd
+{
+ /**
+ * \brief Rectangular pad module; Hit merging
+ **/
+ class HitMerger {
+ public:
+ typedef std::pair<Hit, std::vector<DigiRec>> inputType;
+ typedef std::pair<std::vector<inputType>, std::vector<inputType>> outputType;
+
+ HitMerger(){};
+ /**
+ * \brief Constructor with placement
+ **/
+ HitMerger(HitFinderModPar params);
+ virtual ~HitMerger(){};
+
+ /* \brief Steering routine for building hits */
+ outputType operator()(std::vector<inputType>& hitsRow1, std::vector<inputType>& hitsRow2);
+
+
+ protected:
+ private:
+ HitMerger(const HitMerger& ref);
+ const HitMerger& operator=(const HitMerger& ref);
+
+ HitFinderModPar fParams; ///< Parameter container
+ };
+
+} // namespace cbm::algo::trd
diff --git a/algo/detectors/trd/HitMerger2D.cxx b/algo/detectors/trd/HitMerger2D.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..b30a2d9082f57327004d2194f963a477beb8eea2
--- /dev/null
+++ b/algo/detectors/trd/HitMerger2D.cxx
@@ -0,0 +1,520 @@
+/* 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 "HitMerger2D.h"
+
+#include "log.hpp"
+
+#include <functional>
+#include <numeric>
+
+
+using std::cout;
+using std::endl;
+using std::vector;
+
+namespace cbm::algo::trd
+{
+
+ //_______________________________________________________________________________
+ HitMerger2D::HitMerger2D(HitFinder2DModPar params) : fParams(params) {}
+
+ //_______________________________________________________________________________
+ HitMerger2D::outputType HitMerger2D::operator()(std::vector<inputType>& hitsRow1, std::vector<inputType>& hitsRow2)
+ {
+ std::vector<inputType*> hitData;
+
+ for (auto& elem : hitsRow1) {
+ hitData.push_back(&elem);
+ }
+ for (auto& elem : hitsRow2) {
+ hitData.push_back(&elem);
+ }
+
+ // Sort step below is in principle needed, but wasn't present in original implementation
+ // TO DO: Clarify
+ //std::sort(hitData.begin(), hitData.end(), [](const auto& h0, const auto& h1) { return h0->first.Time() < h1->first.Time(); });
+
+ for (size_t ih = 0; ih < hitData.size(); ih++) {
+ Hit* h0 = &hitData[ih]->first;
+ std::vector<DigiRec>* h0digis = &hitData[ih]->second;
+
+ if (h0->IsUsed()) continue;
+
+ for (size_t jh = ih + 1; jh < hitData.size(); jh++) {
+ Hit* h1 = &hitData[jh]->first;
+ std::vector<DigiRec>* h1digis = &hitData[jh]->second;
+
+ if (h1->IsUsed()) continue;
+
+ // basic check on Time
+ if (h1->Time() < 4000 - h0->Time()) continue; // TODO check with preoper time calibration
+ // skip next hits for being too far (10us) in the future
+ if (h1->Time() > 10000 + h0->Time()) break;
+
+ // basic check on Col
+ if (std::abs(h1->X() - h0->X()) > 2 * fParams.padSizeX) continue;
+
+ // basic check on Row
+ if (std::abs(h1->Y() - h0->Y()) > 2 * fParams.padSizeY) continue;
+
+ // go to topologic checks
+ const int a0 = CheckMerge(h0digis, h1digis);
+ if (a0 == 0) continue;
+
+ auto [hitFlag, hitF] = ProjectDigis(a0 < 0 ? h0digis : h1digis, a0 < 0 ? h1digis : h0digis);
+ // call the working algorithm
+ if (MergeHits(h0, a0, hitF)) {
+ h0->SetRowCross((bool) h1);
+ h1->SetRefId(-1);
+ }
+ }
+ }
+
+ const auto ret1 =
+ std::remove_if(hitsRow1.begin(), hitsRow1.end(), [](const auto& obj) { return obj.first.IsUsed(); });
+ hitsRow1.erase(ret1, hitsRow1.end());
+
+ const auto ret2 =
+ std::remove_if(hitsRow2.begin(), hitsRow2.end(), [](const auto& obj) { return obj.first.IsUsed(); });
+ hitsRow2.erase(ret2, hitsRow2.end());
+
+ return std::make_pair(std::move(hitsRow1), std::move(hitsRow2));
+ }
+
+
+ //_______________________________________________________________________________
+ int HitMerger2D::CheckMerge(std::vector<DigiRec>* cid, std::vector<DigiRec>* 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>::const_iterator jp[2];
+
+ for (int rowId(0); rowId < 2; rowId++) {
+ double rtMax = 0.;
+ double mdMax = 0.;
+ const std::vector<DigiRec>* vcid[2] = {cid, cjd};
+ for (auto id = vcid[rowId]->begin(); id != 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 != 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] != cjd->end()) jp[1]++;
+ if (jp[1] != cjd->end()) {
+ r = (*jp[1]).GetRectCharge(on);
+ jp[1]++;
+ if (jp[1] != 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] != cid->begin()) jp[0]--;
+ if (jp[0] != 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 HitMerger2D::MergeHits(Hit* h, int a0, HitFactory2D& hitF)
+ {
+ int n0(hitF.fSignal.size() - 2);
+ auto [dx, dy] = hitF.GetDxDy(n0);
+
+ int typ(hitF.GetHitClass());
+ // get position correction [pw]
+ double xcorr = hitF.GetXcorr(dx, typ, 1) / fParams.padSizeX, xcorrBias(xcorr);
+ if (hitF.IsBiasX()) {
+ typ = hitF.GetHitRcClass(a0);
+ int xmap = hitF.vyM & 0xff;
+ switch (n0) {
+ case 4:
+ if (dx < 0)
+ xcorrBias += (hitF.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;
+ hitF.RecenterXoffset(dx);
+ xcorrBias = hitF.GetXcorr(dx, typ, 2) / fParams.padSizeX;
+ }
+ else {
+ dx -= xcorr;
+ hitF.RecenterXoffset(dx);
+ if (typ < 2)
+ xcorrBias = hitF.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;
+ hitF.RecenterXoffset(dx);
+ xcorrBias = hitF.GetXcorr(dx, typ, 2) / fParams.padSizeX;
+ break;
+ }
+ }
+ else {
+ if (typ) xcorrBias += (dx < 0 ? 1 : -1) * 0.0293;
+ }
+ std::tie(dx, dy) = hitF.CorrectPosition(dx, dy, xcorrBias, fParams.padSizeX, fParams.padSizeY);
+
+ double edx(1), edy(1), edt(60), time(-21), tdrift(100), e(200);
+ CalibrateHit(h, dx, dy, edx, edy, edt, time, tdrift, e, hitF);
+
+ return kTRUE;
+ }
+
+ //_______________________________________________________________________________
+ void HitMerger2D::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 HitFactory2D& hitF)
+ {
+ const ROOT::Math::XYZVector& local_pad = fParams.rowPar[hitF.vrM].padPar[hitF.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 * (hitF.vt0 + time) - tdrift + 30.29 + fHitTimeOff, edt);
+ h->SetELoss(eloss);
+ h->SetClassType();
+ h->SetMaxType(hitF.IsMaxTilt());
+ h->SetOverFlow(hitF.HasOvf());
+ }
+
+
+ //_______________________________________________________________________________
+ std::pair<int, HitFactory2D> HitMerger2D::ProjectDigis(std::vector<DigiRec>* cid, std::vector<DigiRec>* 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 (cid->empty()) {
+ L_(debug) << "TrdModuleRec2D::ProjectDigis : Request cl id " << cid << " not found.";
+ return std::make_pair(0, HitFactory2D(0));
+ }
+ HitFactory2D hitF(fParams.rowPar[0].padPar.size());
+
+ std::vector<HitFactory2D::signal>& hitSig = hitF.fSignal;
+ hitF.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;
+
+ // check integrity of input
+ if (cid == nullptr || cjd == nullptr || cid->empty() || cjd->empty()) {
+ L_(debug) << "TrdModuleRec2D::ProjectDigis : Requested cluster not found.";
+ return std::make_pair(0, HitFactory2D(0));
+ }
+
+ vector<DigiRec>::const_iterator i1 = cjd->begin();
+
+ for (auto i = cid->begin(); i != cid->end(); i++, j++) {
+ const DigiRec* dg = &(*i);
+ const DigiRec* dg0 = nullptr;
+
+ // initialize
+ if (col < 0) {
+ hitF.vrM = GetPadRowCol(dg->GetAddressChannel(), col);
+ hitF.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 ((dg0 = (i1 != cjd->end()) ? &(*i1) : nullptr)) {
+ row1 = GetPadRowCol(dg0->GetAddressChannel(), col1);
+ if (!step) step = hitF.vrM - row1;
+ if (col1 == col0) {
+ r += dg0->GetRectCharge(on);
+ if (on) nr++;
+ }
+ else
+ dg0 = nullptr;
+ }
+ if (step == 1 && i1 != cjd->begin()) {
+ const auto dg1 = &(*(i1 - 1));
+ GetPadRowCol(dg1->GetAddressChannel(), col1);
+ if (col1 == col0 - 1) {
+ t += dg1->GetTiltCharge(on);
+ if (on) nt++;
+ }
+ }
+ if (step == -1 && i1 != cjd->end() && i1 + 1 != cjd->end()) {
+ const auto dg1 = &(*(i1 + 1));
+ GetPadRowCol(dg1->GetAddressChannel(), col1);
+ if (col1 == col0 + 1) {
+ t += dg1->GetTiltCharge(on);
+ if (on) nt++;
+ }
+ }
+ if (dg0) i1++;
+
+ //TO DO: two blocks below might be mergable
+ // process tilt signal/time
+ char ddt = dg->GetTiltTime() - hitF.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;
+ hitF.vcM = j;
+ hitF.SetMaxTilt(1);
+ hitF.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() - hitF.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;
+ hitF.vcM = j;
+ hitF.SetMaxTilt(0);
+ hitF.viM = hitSig.size();
+ }
+ }
+ else
+ err = 300.;
+
+ hitSig.emplace_back(r, err, ddt, xc, 0.);
+ xc += 0.5;
+ NR += nr;
+ NT += nt;
+ } // for (auto i = cid->begin(); i != cid->end(); i++, j++)
+
+
+ //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.);
+ hitF.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 - hitF.viM;
+ if (nR == hitF.viM) {
+ hitF.SetSymmHit(1);
+ if (hitSig.size() % 2) {
+ double LS(0.), RS(0.);
+ for (UChar_t idx(0); idx < hitF.viM; idx++)
+ LS += hitSig[idx].s;
+ for (uint idx(hitF.viM + 1); idx < hitSig.size(); idx++)
+ RS += hitSig[idx].s;
+ hitF.SetLeftSgn(LS < RS ? 0 : 1);
+ }
+ }
+ else {
+ hitF.SetSymmHit(0);
+ if (hitF.viM < nR)
+ hitF.SetLeftHit(0);
+ else if (hitF.viM > nR)
+ hitF.SetLeftHit(1);
+ }
+ // recenter time and space profile
+ hitF.vt0 += dt0;
+ for (auto& sig : hitSig) {
+ sig.t -= dt0;
+ sig.x -= hitF.vcM;
+ }
+ hitF.vcM += col;
+
+ // check if all signals have same significance
+ const int nmissX = 2 * nsr - NR;
+ if (nmissX) {
+ hitF.SetBiasX(1);
+ for (UChar_t idx(1); idx < hitF.viM; idx++) {
+ if (hitSig[idx].xe > 0.) continue; // select rect signals
+ if (hitSig[idx].se > re * 0.8) hitF.SetBiasXleft(1);
+ }
+ if (hitSig[hitF.viM].xe <= 0. && hitSig[hitF.viM].se > re * 0.8) hitF.SetBiasXmid(1);
+ for (UChar_t idx(hitF.viM + 1); idx < hitSig.size() - 1; idx++) {
+ if (hitSig[idx].xe > 0.) continue; // select rect signals
+ if (hitSig[idx].se > re * 0.8) hitF.SetBiasXright(1);
+ }
+ }
+ else {
+ hitF.SetBiasX(0);
+ }
+
+ const int nmissY = 2 * n0 - 2 * nsr - NT;
+ if (nmissY) {
+ hitF.SetBiasY();
+ for (UChar_t idx(1); idx < hitF.viM; idx++) {
+ if (hitSig[idx].xe > 0. && hitSig[idx].se > te * 0.8) hitF.SetBiasYleft(1); // select tilt signals
+ }
+ if (hitSig[hitF.viM].xe > 0. && hitSig[hitF.viM].se > te * 0.8) hitF.SetBiasYmid(1);
+ for (UChar_t idx(hitF.viM + 1); idx < hitSig.size() - 1; idx++) {
+ if (hitSig[idx].xe > 0. && hitSig[idx].se > te * 0.8) hitF.SetBiasYright(1); // select tilt signals
+ }
+ }
+ else {
+ hitF.SetBiasY(0);
+ }
+
+ if (ovf < 0) hitF.SetOvf();
+ return std::make_pair(ovf * (hitSig.size() - 2), hitF);
+ }
+
+ int HitMerger2D::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/HitMerger2D.h b/algo/detectors/trd/HitMerger2D.h
new file mode 100644
index 0000000000000000000000000000000000000000..113f029e0d7e2150ae6e4fbcf81fb0ef3dd5f7b3
--- /dev/null
+++ b/algo/detectors/trd/HitMerger2D.h
@@ -0,0 +1,97 @@
+/* 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 "DigiRec.h"
+#include "Hit.h"
+#include "HitFactory2D.h"
+#include "HitFinder2DPars.h"
+#include "Math/Rotation3D.h"
+#include "Math/Vector3Dfwd.h"
+
+#include <cstdint>
+#include <tuple>
+#include <unordered_map>
+#include <vector>
+
+using std::vector;
+class DigiRec;
+
+namespace cbm::algo::trd
+{
+
+ /** @class HitMerger2D
+ ** @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 HitMerger2D {
+ public:
+ typedef std::pair<Hit, std::vector<DigiRec>> inputType;
+ typedef std::pair<std::vector<inputType>, std::vector<inputType>> outputType;
+
+ /** \brief Default constructor.*/
+ HitMerger2D(){};
+
+ /** \brief Constructor with placement */
+ HitMerger2D(HitFinder2DModPar params);
+
+ virtual ~HitMerger2D(){};
+
+ /** \brief Steering routine for building hits **/
+ outputType operator()(std::vector<inputType>& hitsRow1, std::vector<inputType>& hitsRow2);
+
+ std::pair<int, HitFactory2D> ProjectDigis(std::vector<DigiRec>* cid, std::vector<DigiRec>* cjd);
+ /** \brief Implement topologic cuts for hit merging
+ * \return index of anode wire wrt to boundary or 0 if check fails
+ */
+ int CheckMerge(std::vector<DigiRec>* cid, std::vector<DigiRec>* 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, HitFactory2D& hitF);
+
+
+ /** \brief Time offset to synchronize TRD2D hits to the rest of detectors
+ * \param dt offset in [ns]
+ */
+ void SetHitTimeOffset(int dt) { fHitTimeOff = dt; }
+
+ private:
+ HitMerger2D(const HitMerger2D& ref);
+ const HitMerger2D& operator=(const HitMerger2D& ref);
+
+ Hit MakeHit(int cId);
+
+ const float fClk = CbmTrdDigi::Clk(CbmTrdDigi::eCbmTrdAsicType::kFASP);
+
+ HitFinder2DModPar fParams; ///< Parameter container
+
+ 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, const HitFactory2D& hitF);
+
+ /** \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);
+
+ int fHitTimeOff = 0; //! hit time offset for synchronization
+ };
+
+
+} // namespace cbm::algo::trd
diff --git a/algo/detectors/trd/Hitfind.cxx b/algo/detectors/trd/Hitfind.cxx
index e0f0cb9043e8d959cd8b7dc2ba9d78c5e1f69c28..f1e3efc0e870d8fc3c89aa12d569caa0566f86ab 100644
--- a/algo/detectors/trd/Hitfind.cxx
+++ b/algo/detectors/trd/Hitfind.cxx
@@ -53,6 +53,7 @@ namespace cbm::algo::trd
fRowList.emplace_back(module.address, false, row);
}
fHitFind[module.address] = std::make_unique<cbm::algo::trd::HitFinder>(params);
+ fHitMerge[module.address] = std::make_unique<cbm::algo::trd::HitMerger>(params);
fClusterBuild[module.address] = std::make_unique<cbm::algo::trd::Clusterizer>(params);
fModId[module.address] = fModList.size();
fModList.emplace_back(module.address, false, module.rowPar.size(), module.rowPar[0].padPar.size());
@@ -91,6 +92,7 @@ namespace cbm::algo::trd
fRowList.emplace_back(module.address, true, row);
}
fHitFind2d[module.address] = std::make_unique<cbm::algo::trd::HitFinder2D>(params);
+ fHitMerge2d[module.address] = std::make_unique<cbm::algo::trd::HitMerger2D>(params);
fClusterBuild2d[module.address] = std::make_unique<cbm::algo::trd::Clusterizer2D>(params);
fModId[module.address] = fModList.size();
fModList.emplace_back(module.address, true, module.rowPar.size(), module.rowPar[0].padPar.size());
@@ -104,7 +106,7 @@ namespace cbm::algo::trd
// ----- Execution -------------------------------------------------------
Hitfind::resultType Hitfind::operator()(gsl::span<CbmTrdDigi> digiIn)
{
- constexpr bool DebugCheckInput = false;
+ constexpr bool DebugCheckInput = true;
// --- Output data
resultType result = {};
@@ -112,26 +114,17 @@ namespace cbm::algo::trd
auto& monitor = std::get<1>(result);
// Intermediate digi storage variables (digi, index) per module and row
- std::unordered_map<int, std::vector<std::vector<std::pair<CbmTrdDigi, int32_t>>>> digiBuffer; //[modAddress][row]
+ std::vector<std::vector<std::vector<std::pair<CbmTrdDigi, int32_t>>>> digiBuffer; //[modAddress][row]
+ digiBuffer.resize(fModList.size());
- // Intermediate 1D clusters per module and row
- std::unordered_map<int, std::vector<std::vector<Cluster>>> clusterBuffer; //[modAddress][row]
-
- // Intermediate 2D clusters per module and row
- std::unordered_map<int, std::vector<std::vector<Cluster2D>>> clusterBuffer2d; //[modAddress][row]
+ // Intermediate hits per module and row
+ std::vector<std::vector<hitDataType>> hitBuffer; //[row]
+ hitBuffer.resize(fRowList.size());
// Initialize storage buffers
for (size_t mod = 0; mod < fModList.size(); mod++) {
- const int address = std::get<0>(fModList[mod]);
- const bool is2D = std::get<1>(fModList[mod]);
const size_t numRows = std::get<2>(fModList[mod]);
- digiBuffer[address].resize(numRows);
- if (is2D) {
- clusterBuffer2d[address].resize(numRows);
- }
- else {
- clusterBuffer[address].resize(numRows);
- }
+ digiBuffer[mod].resize(numRows);
}
// Loop over the digis array and store the digis in separate vectors for
@@ -157,42 +150,23 @@ namespace cbm::algo::trd
const size_t modId = fModId[address];
const size_t numCols = std::get<3>(fModList[modId]);
const int row = digi->GetAddressChannel() / numCols;
- digiBuffer[address][row].emplace_back(*digi, idigi);
+ digiBuffer[modId][row].emplace_back(*digi, idigi);
}
monitor.sortTime = xpu::pop_timer();
- xpu::push_timer("BuildClusters");
- xpu::t_add_bytes(digiIn.size_bytes());
-
- // Cluster building
- CBM_PARALLEL_FOR(schedule(dynamic))
- for (size_t row = 0; row < fRowList.size(); row++) {
-
- const int address = std::get<0>(fRowList[row]);
- const bool is2D = std::get<1>(fRowList[row]);
- const size_t rowInMod = std::get<2>(fRowList[row]);
- const auto& digiInput = digiBuffer[address][rowInMod];
- if (is2D) {
- clusterBuffer2d[address][rowInMod] =
- (*fClusterBuild2d[address])(digiInput, 0.); // Number is TS start time (T0)
- }
- else {
- clusterBuffer[address][rowInMod] = (*fClusterBuild[address])(digiInput);
- }
- }
- monitor.timeClusterize = xpu::pop_timer();
-
- // Hit finding
+ // Hit finding results
PODVector<Hit> hitsFlat; // hit storage
- PODVector<size_t> modSizes; // nHits per modules
- PODVector<uint> modAddresses; // address of modules
+ PODVector<size_t> rowSizes; // nHits per row
+ PODVector<uint> rowAddresses; // address of row
// Prefix array for parallelization
std::vector<size_t> hitsPrefix;
std::vector<size_t> sizePrefix;
std::vector<size_t> addrPrefix;
- xpu::push_timer("FindHits");
+ xpu::push_timer("BuildHits");
+ xpu::t_add_bytes(digiIn.size_bytes());
+
CBM_PARALLEL()
{
const int ithread = openmp::GetThreadNum();
@@ -205,44 +179,84 @@ namespace cbm::algo::trd
addrPrefix.resize(nthreads + 1);
}
- std::vector<Hit> local_hits;
- std::vector<size_t> local_sizes;
- std::vector<uint> local_addresses;
+ // Cluster and hit building
+ CBM_OMP(for schedule(dynamic))
+ for (size_t row = 0; row < fRowList.size(); row++) {
+ const int address = std::get<0>(fRowList[row]);
+ const bool is2D = std::get<1>(fRowList[row]);
+ const size_t rowInMod = std::get<2>(fRowList[row]);
+ const size_t modId = fModId[address];
+ const auto& digiInput = digiBuffer[modId][rowInMod];
+ if (is2D) {
+ auto clusters = (*fClusterBuild2d[address])(digiInput, 0.); // Number is TS start time (T0)
+ hitBuffer[row] = (*fHitFind2d[address])(&clusters);
+ }
+ else {
+ auto clusters = (*fClusterBuild[address])(digiInput);
+ hitBuffer[row] = (*fHitFind[address])(&clusters);
+ }
+ }
- CBM_OMP(for schedule(dynamic) nowait)
- for (size_t mod = 0; mod < fModList.size(); mod++) {
- const int address = std::get<0>(fModList[mod]);
- const bool is2D = std::get<1>(fModList[mod]);
-
- // Lambda expression for vector concatenation
- auto concatVec = [](auto& acc, const auto& innerVec) {
- acc.insert(acc.end(), innerVec.begin(), innerVec.end());
- return std::move(acc);
- };
-
- // Flatten the input cluster vector of vectors and build hits
- std::vector<Hit> mod_hits;
+ // Row-merging for even rows
+ CBM_OMP(for schedule(dynamic))
+ for (size_t row = 0; row < fRowList.size() / 2; row++) {
+ const size_t row1 = 2 * row;
+ const size_t row2 = 2 * row + 1;
+ const int address = std::get<0>(fRowList[row1]);
+ const bool is2D = std::get<1>(fRowList[row1]);
+ if (row2 >= fRowList.size() || std::get<0>(fRowList[row2]) != address) {
+ continue;
+ }
if (is2D) {
- auto& buffer = clusterBuffer2d[address];
- auto clusters = std::accumulate(buffer.begin(), buffer.end(), std::vector<Cluster2D>(), concatVec);
- //sort(clusters.begin(), clusters.end(), [](const auto& a, const auto& b) { return a.GetStartTime() < b.GetStartTime(); });
- mod_hits = (*fHitFind2d[address])(&clusters);
+ std::tie(hitBuffer[row1], hitBuffer[row2]) = (*fHitMerge2d[address])(hitBuffer[row1], hitBuffer[row2]);
}
else {
- auto& buffer = clusterBuffer[address];
- auto clusters = std::accumulate(buffer.begin(), buffer.end(), std::vector<Cluster>(), concatVec);
- //sort(clusters.begin(), clusters.end(), [](const auto& a, const auto& b) { return a.GetStartTime() < b.GetStartTime(); });
- mod_hits = (*fHitFind[address])(&clusters);
+ std::tie(hitBuffer[row1], hitBuffer[row2]) = (*fHitMerge[address])(hitBuffer[row1], hitBuffer[row2]);
+ }
+ }
+
+ // Row-merging for odd rows
+ CBM_OMP(for schedule(dynamic))
+ for (size_t row = 0; row < fRowList.size() / 2; row++) {
+ const size_t row1 = 2 * row + 1;
+ const size_t row2 = 2 * row + 2;
+ if (row2 >= fRowList.size()) {
+ continue;
}
+ const int address = std::get<0>(fRowList[row1]);
+ const bool is2D = std::get<1>(fRowList[row1]);
+ if (std::get<0>(fRowList[row2]) != address) {
+ continue;
+ }
+ if (is2D) {
+ std::tie(hitBuffer[row1], hitBuffer[row2]) = (*fHitMerge2d[address])(hitBuffer[row1], hitBuffer[row2]);
+ }
+ else {
+ std::tie(hitBuffer[row1], hitBuffer[row2]) = (*fHitMerge[address])(hitBuffer[row1], hitBuffer[row2]);
+ }
+ }
+
+ std::vector<Hit> local_hits;
+ std::vector<size_t> local_sizes;
+ std::vector<uint> local_addresses;
+
+ // Buffer merging
+ CBM_OMP(for schedule(dynamic) nowait)
+ for (size_t row = 0; row < fRowList.size(); row++) {
+ const int address = std::get<0>(fRowList[row]);
+ auto& hits = hitBuffer[row];
+ std::vector<Hit> row_hits;
+ std::transform(hits.begin(), hits.end(), std::back_inserter(row_hits), [](const auto& p) { return p.first; });
+
// store partition size
- local_sizes.push_back(mod_hits.size());
+ local_sizes.push_back(row_hits.size());
// store hw address of partition
local_addresses.push_back(address);
// Append clusters to output
- local_hits.insert(local_hits.end(), std::make_move_iterator(mod_hits.begin()),
- std::make_move_iterator(mod_hits.end()));
+ local_hits.insert(local_hits.end(), std::make_move_iterator(row_hits.begin()),
+ std::make_move_iterator(row_hits.end()));
}
hitsPrefix[ithread + 1] = local_hits.size();
@@ -257,20 +271,21 @@ namespace cbm::algo::trd
addrPrefix[i] += addrPrefix[i - 1];
}
hitsFlat.resize(hitsPrefix[nthreads]);
- modSizes.resize(sizePrefix[nthreads]);
- modAddresses.resize(addrPrefix[nthreads]);
+ rowSizes.resize(sizePrefix[nthreads]);
+ rowAddresses.resize(addrPrefix[nthreads]);
}
std::move(local_hits.begin(), local_hits.end(), hitsFlat.begin() + hitsPrefix[ithread]);
- std::move(local_sizes.begin(), local_sizes.end(), modSizes.begin() + sizePrefix[ithread]);
- std::move(local_addresses.begin(), local_addresses.end(), modAddresses.begin() + addrPrefix[ithread]);
+ std::move(local_sizes.begin(), local_sizes.end(), rowSizes.begin() + sizePrefix[ithread]);
+ std::move(local_addresses.begin(), local_addresses.end(), rowAddresses.begin() + addrPrefix[ithread]);
}
+
// Monitoring
monitor.timeHitfind = xpu::pop_timer();
monitor.numDigis = digiIn.size();
monitor.numHits = hitsFlat.size();
// Create ouput vector
- hitsOut = PartitionedVector(std::move(hitsFlat), modSizes, modAddresses);
+ hitsOut = PartitionedVector(std::move(hitsFlat), rowSizes, rowAddresses);
// Ensure hits are time sorted
CBM_PARALLEL_FOR(schedule(dynamic))
diff --git a/algo/detectors/trd/Hitfind.h b/algo/detectors/trd/Hitfind.h
index 93dc592f647c4899f26eb3fc0cefc7b30e4294ec..ee2566317fc1f1fed5839b403170cb438e62cf63 100644
--- a/algo/detectors/trd/Hitfind.h
+++ b/algo/detectors/trd/Hitfind.h
@@ -5,12 +5,15 @@
#pragma once
#include "CbmTrdDigi.h"
+#include "DigiRec.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/HitMerger.h"
+#include "trd/HitMerger2D.h"
#include "trd/Hitfind2DSetup.h"
#include "trd/HitfindSetup.h"
@@ -55,6 +58,7 @@ namespace cbm::algo::trd
public:
typedef std::tuple<PartitionedVector<Hit>, HitfindMonitorData> resultType;
+ typedef std::pair<Hit, std::vector<DigiRec>> hitDataType;
/** @brief Algorithm execution
** @param fles timeslice to hitfind
@@ -76,6 +80,10 @@ namespace cbm::algo::trd
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 Hit merging algorithms per module */
+ std::unordered_map<int, std::unique_ptr<cbm::algo::trd::HitMerger2D>> fHitMerge2d;
+ std::unordered_map<int, std::unique_ptr<cbm::algo::trd::HitMerger>> fHitMerge;
+
/** @brief List of modules (address, type flag (true = 2D), numRows, numCols) */
std::vector<std::tuple<int, bool, size_t, size_t>> fModList;
diff --git a/reco/tasks/CbmTaskTrdHitFinder.cxx b/reco/tasks/CbmTaskTrdHitFinder.cxx
index 3902bbd450419c2c2e31809a20369feb9ffd280d..b9785e20b15aa968f7a3682ef998ed2928e9ba87 100644
--- a/reco/tasks/CbmTaskTrdHitFinder.cxx
+++ b/reco/tasks/CbmTaskTrdHitFinder.cxx
@@ -75,8 +75,6 @@ InitStatus CbmTaskTrdHitFinder::Init()
//_____________________________________________________________________
void CbmTaskTrdHitFinder::Exec(Option_t* /*option*/)
{
-
- size_t numHitsCall = 0;
fClusters->clear();
fHits->clear();
diff --git a/reco/tasks/CbmTaskTrdHitFinderParWrite.cxx b/reco/tasks/CbmTaskTrdHitFinderParWrite.cxx
index 75ba2be3de7f32c2faf5f68e81266dc55288328f..247f73e738888386a51e3440416ce8ccb1c62ef3 100644
--- a/reco/tasks/CbmTaskTrdHitFinderParWrite.cxx
+++ b/reco/tasks/CbmTaskTrdHitFinderParWrite.cxx
@@ -90,7 +90,7 @@ InitStatus CbmTaskTrdHitFinderParWrite::Init()
setup2D.modules.resize(trd2DModules.size());
// Fill 2D setup files
- for (size_t mod; mod < trd2DModules.size(); mod++) {
+ for (size_t mod = 0; mod < trd2DModules.size(); mod++) {
const uint16_t moduleId = trd2DModules[mod];
// Get Geometry parameters for module
@@ -151,7 +151,7 @@ InitStatus CbmTaskTrdHitFinderParWrite::Init()
// Fill 1D setup files
- for (size_t mod; mod < trdModules.size(); mod++) {
+ for (size_t mod = 0; mod < trdModules.size(); mod++) {
const uint16_t moduleId = trdModules[mod];
// Get Geometry parameters for module
diff --git a/reco/tasks/CbmTaskTrdHitFinderParWrite.h b/reco/tasks/CbmTaskTrdHitFinderParWrite.h
index 24ff1fbce5dc9ab5df72f5efe3c0af0e7219ead2..6543c220669067a4eff512cdd3005fa581241b6e 100644
--- a/reco/tasks/CbmTaskTrdHitFinderParWrite.h
+++ b/reco/tasks/CbmTaskTrdHitFinderParWrite.h
@@ -44,7 +44,7 @@ class CbmTaskTrdHitFinderParWrite : public FairTask {
virtual void SetParContainers();
/** \brief Executed task **/
- virtual void Exec(Option_t* option){};
+ virtual void Exec(Option_t*){};
/** Finish task **/
virtual void Finish(){};