From 60d114eb22fd85e33d875368033e1f0e2ade32b9 Mon Sep 17 00:00:00 2001 From: Slavomir Hnatic <hnatics@jinr.ru> Date: Wed, 26 Mar 2025 21:50:09 +0300 Subject: [PATCH 1/2] Fast clusterhitfinder: v2.0.8 (author: V.Krylov) --- detectors/tpc/clusterHitFinder/CMakeLists.txt | 29 +- .../clusterHitFinder/fast/Tpc2dClusterFast.h | 2 +- .../tpc/clusterHitFinder/fast/TpcAdcHit.h | 68 -- .../tpc/clusterHitFinder/fast/TpcCluster.h | 520 ------------- .../tpc/clusterHitFinder/fast/TpcClusterHit.h | 52 -- .../fast/TpcClusterHitFinderFast.cxx | 96 ++- .../fast/TpcClusterHitFinderFast.h | 5 +- .../clusterHitFinder/fast/TpcClusterUnits.h | 683 ------------------ .../clusterHitFinder/fast/TpcEventClusters.h | 102 --- .../tpc/clusterHitFinder/fast/TpcEventObj.cxx | 339 +++++++++ .../tpc/clusterHitFinder/fast/TpcEventObj.h | 89 +++ .../tpc/clusterHitFinder/fast/TpcPadCluster.h | 370 ---------- .../tpc/clusterHitFinder/fast/TpcPadHit.h | 52 -- .../clusterHitFinder/fast/TpcRowClusters.h | 276 ------- .../clusterHitFinder/fast/TpcSectorClusters.h | 142 ---- .../clusterHitFinder/fast/TpcSectorObj.cxx | 401 ++++++++++ .../tpc/clusterHitFinder/fast/TpcSectorObj.h | 77 ++ .../fast/clustering/TpcAdcGauss.h | 317 ++++++++ .../fast/clustering/TpcAdcHit.h | 77 ++ .../fast/clustering/TpcAdcStat.h | 52 ++ .../fast/clustering/TpcCluster.h | 313 ++++++++ .../fast/clustering/TpcClusters.cxx | 16 + .../fast/clustering/TpcClusters.h | 190 +++++ .../fast/clustering/TpcEventPayload.h | 144 ++++ .../clusterHitFinder/fast/clustering/TpcHit.h | 108 +++ .../fast/clustering/TpcPadCluster.h | 308 ++++++++ .../fast/clustering/TpcPadHit.h | 79 ++ .../fast/clustering/TpcRowClusters.h | 406 +++++++++++ .../fast/clustering/TpcSampaPack.h | 173 +++++ .../fast/clustering/TpcSectorGeo.h | 62 ++ .../fast/clustering/TpcTrack.h | 209 ++++++ .../tpc/clusterHitFinder/fast/utility.cxx | 92 +++ detectors/tpc/clusterHitFinder/fast/utility.h | 37 + 33 files changed, 3579 insertions(+), 2307 deletions(-) delete mode 100644 detectors/tpc/clusterHitFinder/fast/TpcAdcHit.h delete mode 100644 detectors/tpc/clusterHitFinder/fast/TpcCluster.h delete mode 100644 detectors/tpc/clusterHitFinder/fast/TpcClusterHit.h delete mode 100644 detectors/tpc/clusterHitFinder/fast/TpcClusterUnits.h delete mode 100644 detectors/tpc/clusterHitFinder/fast/TpcEventClusters.h create mode 100644 detectors/tpc/clusterHitFinder/fast/TpcEventObj.cxx create mode 100644 detectors/tpc/clusterHitFinder/fast/TpcEventObj.h delete mode 100644 detectors/tpc/clusterHitFinder/fast/TpcPadCluster.h delete mode 100644 detectors/tpc/clusterHitFinder/fast/TpcPadHit.h delete mode 100644 detectors/tpc/clusterHitFinder/fast/TpcRowClusters.h delete mode 100644 detectors/tpc/clusterHitFinder/fast/TpcSectorClusters.h create mode 100644 detectors/tpc/clusterHitFinder/fast/TpcSectorObj.cxx create mode 100644 detectors/tpc/clusterHitFinder/fast/TpcSectorObj.h create mode 100644 detectors/tpc/clusterHitFinder/fast/clustering/TpcAdcGauss.h create mode 100644 detectors/tpc/clusterHitFinder/fast/clustering/TpcAdcHit.h create mode 100644 detectors/tpc/clusterHitFinder/fast/clustering/TpcAdcStat.h create mode 100644 detectors/tpc/clusterHitFinder/fast/clustering/TpcCluster.h create mode 100644 detectors/tpc/clusterHitFinder/fast/clustering/TpcClusters.cxx create mode 100644 detectors/tpc/clusterHitFinder/fast/clustering/TpcClusters.h create mode 100644 detectors/tpc/clusterHitFinder/fast/clustering/TpcEventPayload.h create mode 100644 detectors/tpc/clusterHitFinder/fast/clustering/TpcHit.h create mode 100644 detectors/tpc/clusterHitFinder/fast/clustering/TpcPadCluster.h create mode 100644 detectors/tpc/clusterHitFinder/fast/clustering/TpcPadHit.h create mode 100644 detectors/tpc/clusterHitFinder/fast/clustering/TpcRowClusters.h create mode 100644 detectors/tpc/clusterHitFinder/fast/clustering/TpcSampaPack.h create mode 100644 detectors/tpc/clusterHitFinder/fast/clustering/TpcSectorGeo.h create mode 100644 detectors/tpc/clusterHitFinder/fast/clustering/TpcTrack.h create mode 100644 detectors/tpc/clusterHitFinder/fast/utility.cxx create mode 100644 detectors/tpc/clusterHitFinder/fast/utility.h diff --git a/detectors/tpc/clusterHitFinder/CMakeLists.txt b/detectors/tpc/clusterHitFinder/CMakeLists.txt index aa751b94..c75fd768 100644 --- a/detectors/tpc/clusterHitFinder/CMakeLists.txt +++ b/detectors/tpc/clusterHitFinder/CMakeLists.txt @@ -31,6 +31,10 @@ link_directories(${LINK_DIRECTORIES}) # List of source files set(SRCS fast/TpcClusterHitFinderFast.cxx + fast/TpcEventObj.cxx + fast/TpcSectorObj.cxx + fast/utility.cxx + fast/clustering/TpcClusters.cxx hitProducer/MpdTpcHitProducer.cxx mlem/TpcClusterHitFinderMlem.cxx mlem/MpdTpc2dCluster.cxx @@ -40,15 +44,22 @@ set(SRCS Set(HEADERS fast/TpcClusterHitFinderFast.h fast/Tpc2dClusterFast.h - fast/TpcAdcHit.h - fast/TpcCluster.h - fast/TpcClusterHit.h - fast/TpcClusterUnits.h - fast/TpcEventClusters.h - fast/TpcPadCluster.h - fast/TpcPadHit.h - fast/TpcRowClusters.h - fast/TpcSectorClusters.h + fast/TpcEventObj.h + fast/TpcSectorObj.h + fast/utility.h + fast/clustering/TpcAdcGauss.h + fast/clustering/TpcAdcHit.h + fast/clustering/TpcAdcStat.h + fast/clustering/TpcCluster.h + fast/clustering/TpcClusters.h + fast/clustering/TpcEventPayload.h + fast/clustering/TpcHit.h + fast/clustering/TpcPadCluster.h + fast/clustering/TpcPadHit.h + fast/clustering/TpcRowClusters.h + fast/clustering/TpcSampaPack.h + fast/clustering/TpcSectorGeo.h + fast/clustering/TpcTrack.h hitProducer/MpdTpcHitProducer.h mlem/TpcClusterHitFinderMlem.h mlem/MpdTpc2dCluster.h diff --git a/detectors/tpc/clusterHitFinder/fast/Tpc2dClusterFast.h b/detectors/tpc/clusterHitFinder/fast/Tpc2dClusterFast.h index b3cf1865..f742b1e9 100644 --- a/detectors/tpc/clusterHitFinder/fast/Tpc2dClusterFast.h +++ b/detectors/tpc/clusterHitFinder/fast/Tpc2dClusterFast.h @@ -18,7 +18,7 @@ #define TPC2DCLUSTERFAST_HH #include "AbstractTpc2dCluster.h" -#include "TpcEventClusters.h" +#include "TpcEventObj.h" class Tpc2dClusterFast : public AbstractTpc2dCluster { public: diff --git a/detectors/tpc/clusterHitFinder/fast/TpcAdcHit.h b/detectors/tpc/clusterHitFinder/fast/TpcAdcHit.h deleted file mode 100644 index 1f445ccd..00000000 --- a/detectors/tpc/clusterHitFinder/fast/TpcAdcHit.h +++ /dev/null @@ -1,68 +0,0 @@ -// -// Joint Institute for Nuclear Research (JINR), Dubna, Russia, Jun-2024 -// MPD TPC Clustering objects for NICA project -// author: Viktor A.Krylov JINR/LNP -// email: kryman@jinr.ru -// -#ifndef ADCHIT_H -#define ADCHIT_H - -#include <stdlib.h> -#include <iomanip> -#include <iostream> -#include <sstream> -#include <ostream> -#include <cassert> - -namespace tpcClustering { -class AdcHit { - uint _nTimeBin; // timebin number [0..310) - float _fAdc; - int _nTrackId; // Track number for TPC hit -public: - AdcHit() : _fAdc(0), _nTimeBin(0) {} - AdcHit(int n, float f) : _nTimeBin(n), _fAdc(f) {} - AdcHit(int n, float f, int nTrackId) : _nTimeBin(n), _fAdc(f), _nTrackId(nTrackId) {} - AdcHit(const AdcHit &r) : _nTimeBin(r._nTimeBin), _fAdc(r._fAdc), _nTrackId(r._nTrackId) {} - ~AdcHit() {} - const AdcHit &operator=(const AdcHit &r) - { - if (this != &r) { - _nTimeBin = r._nTimeBin; - _fAdc = r._fAdc; - _nTrackId = r._nTrackId; //_nTrackId - } - return *this; - } - inline int getTrackID() const { return _nTrackId; } - inline float getAdc() const { return _fAdc; } - inline void setAdc(float f) { _fAdc = f; } - inline uint getTimeBin() const { return _nTimeBin; } - inline int getHitWay(float fDistance, uint nBins) const - { - return _nTimeBin * fDistance / nBins; - } // distance: m, nBins: 310 - inline int getHitTime(float fDistance, uint nBins, float fRate) const - { - return getHitWay(fDistance, nBins) / fRate; - } // rate: m/c - std::ostringstream &Print(std::ostringstream &oss) const - { - oss << _nTimeBin << ":" << _fAdc; - return oss; - } - friend std::ostringstream &operator<<(std::ostringstream &, const AdcHit &); - friend std::ostringstream &operator<<(std::ostringstream &, const AdcHit *); -}; -inline std::ostringstream &operator<<(std::ostringstream &oss, const AdcHit &r) -{ - oss << r.getTimeBin() << ":" << r.getAdc(); - return oss; -} -inline std::ostringstream &operator<<(std::ostringstream &oss, const AdcHit *p) -{ - oss << p->getTimeBin() << ":" << p->getAdc(); - return oss; -} -} // namespace tpcClustering -#endif \ No newline at end of file diff --git a/detectors/tpc/clusterHitFinder/fast/TpcCluster.h b/detectors/tpc/clusterHitFinder/fast/TpcCluster.h deleted file mode 100644 index 94941354..00000000 --- a/detectors/tpc/clusterHitFinder/fast/TpcCluster.h +++ /dev/null @@ -1,520 +0,0 @@ -// -// Joint Institute for Nuclear Research (JINR), Dubna, Russia, Jun-2024 -// MPD TPC Clustering objects for NICA project -// author: Viktor A.Krylov JINR/LNP -// email: kryman@jinr.ru -// -#ifndef CLUSTER_H -#define CLUSTER_H - -#include "TpcClusterHit.h" -#include "TpcPadCluster.h" - -#define CLUSTER_PAD_MAX 5 - -namespace tpcClustering { -class Cluster : public ClusterHit { - friend class PadCluster; - uint _nTimeBinMin; // PadClusters timebin min range - uint _nTimeBinMax; // PadClusters timebin max range - uint _nPadMin; // PadClusters pad min range - uint _nPadMax; // PadClusters pad max range - std::list<const PadCluster *> _lstPadClusters; - HitRange<uint, float> _AdcMaxTimeRange; - HitRange<uint, float> _AdcMaxPadRange; - HitRange<uint, float> _AdcSumTimeRange; - HitRange<uint, float> _AdcSumPadRange; - -public: - Cluster() : ClusterHit(0), _nTimeBinMin(UINT_MAX), _nTimeBinMax(UINT_MAX), _nPadMin(UINT_MAX), _nPadMax(UINT_MAX) {} - Cluster(int nPad, const AdcHit &r) : ClusterHit() - { - PadCluster *pPadCluster = new PadCluster(nPad, r); - _nTimeBinMin = pPadCluster->getTimeBinMin(); - _nTimeBinMax = pPadCluster->getTimeBinMax(); - _nPadMin = _nPadMax = pPadCluster->getPad(); - _lstPadClusters.push_back(pPadCluster); - } - Cluster(uint nTimeBinMin, uint nTimeBinMax, uint nPadMin, uint nPadMax) : ClusterHit(0) - { - _nTimeBinMin = nTimeBinMin; - _nTimeBinMax = nTimeBinMax; - _nPadMin = nPadMin; - _nPadMax = nPadMax; - } - Cluster(uint nId, const PadCluster *pPadCluster) : ClusterHit(nId) - { - _nTimeBinMin = pPadCluster->getTimeBinMin(); - _nTimeBinMax = pPadCluster->getTimeBinMax(); - _nPadMin = _nPadMax = pPadCluster->getPad(); - _lstPadClusters.push_back(pPadCluster); - } - virtual ~Cluster() - { - while (!_lstPadClusters.empty()) { - delete _lstPadClusters.front(); - _lstPadClusters.pop_front(); - } - } - inline const std::list<const PadCluster *> &getPadClusters() const { return _lstPadClusters; } - bool isOut(uint nPad) const { return _nPadMin && nPad < _nPadMin - 1 || nPad > _nPadMax + 1; } - void resetStat() { _nTimeBinMin = _nTimeBinMax = _nPadMin = _nPadMax = UINT_MAX; } - void calcStat(const PadCluster *pPadCluster) - { - uint nTimeBinMin = pPadCluster->getTimeBinMin(); - uint nTimeBinMax = pPadCluster->getTimeBinMax(); - if (_nTimeBinMin == UINT_MAX && _nTimeBinMax == UINT_MAX) { - _nTimeBinMin = nTimeBinMin; - _nTimeBinMax = nTimeBinMax; - } else { - if (nTimeBinMin < _nTimeBinMin) _nTimeBinMin = nTimeBinMin; - if (nTimeBinMax > _nTimeBinMax) _nTimeBinMax = nTimeBinMax; - } - uint nPad = pPadCluster->getPad(); - if (_nPadMin == UINT_MAX && _nPadMax == UINT_MAX) - _nPadMin = _nPadMax = nPad; - else if (nPad < _nPadMin) - _nPadMin = nPad; - else if (nPad > _nPadMax) - _nPadMax = nPad; - } - void Add(const PadCluster *pPadCluster) - { - calcStat(pPadCluster); - _lstPadClusters.push_back(pPadCluster); - } - // void Add( const Cluster* pCluster ) { std::copy( pCluster->_lstPadClusters.begin(), - // pCluster->_lstPadClusters.end(), _lstPadClusters.end() ); } - void Add(const Cluster *pCluster) - { - std::for_each(pCluster->_lstPadClusters.begin(), pCluster->_lstPadClusters.end(), - [this](const PadCluster *pPadCluster) { - ((PadCluster *)pPadCluster)->setCluster(this); - this->Add(pPadCluster); - }); - } - inline void Clear() { _lstPadClusters.clear(); } - inline uint getTimeBinMin() const { return _nTimeBinMin; } - inline uint getTimeBinMax() const { return _nTimeBinMax; } - inline uint getPadMin() const { return _nPadMin; } - inline uint getPadMax() const { return _nPadMax; } - inline void loadAdcSum(const Cluster *pCluster) - { - _AdcSumTimeRange.SetRange(_nTimeBinMin, _nTimeBinMax); - _AdcSumPadRange.SetRange(_nPadMin, _nPadMax); - if (pCluster) { - const PadCluster *pPadCluster = pCluster->_lstPadClusters.back(); - _AdcSumTimeRange.AddTimeHitAdc(pPadCluster); - _AdcSumPadRange.AddPadHitAdcSum(pPadCluster); - } - for (std::list<const PadCluster *>::const_iterator i = _lstPadClusters.begin(); i != _lstPadClusters.end(); i++) { - _AdcSumTimeRange.AddTimeHitAdc(*i); - _AdcSumPadRange.AddPadHitAdcSum(*i); - } - } - inline float getAdcSum() const - { - float fAdcSum = 0; - for (std::list<const PadCluster *>::const_iterator i = _lstPadClusters.begin(); i != _lstPadClusters.end(); i++) - fAdcSum += (*i)->getAdcHitSum(); - return fAdcSum; - } - // http://www.pm298.ru/dekart.php - // x = (X - a) * cos + (Y - b) * sin - // y = (Y - b) * cos - (X - a) * sin - // X = x * cos - y * sin + a - // Y = y * cos + x * sin + b - - // x = X * cos - Y * sin - // y = X * sin + Y * cos - void TestRotate(int x, int y, float fCos, float fSin, const char *szName) - { - float xCos = x * fCos; - float xSin = x * fSin; - float yCos = y * fCos; - float ySin = y * fSin; - int nX = round(xCos + ySin); - int nY = round(yCos - xSin); - std::cout << szName << "xy:(" << x << ',' << y << ") XY:(" << nX << ',' << nY << ')'; - } - void TestRotate(int x1, int y1, int x2, int y2) - { - int x21 = x2 - x1; - int y21 = y2 - y1; - float f = sqrt(x21 * x21 + y21 * y21); - float fCos = x21 / f; - float fSin = y21 / f; - if (std::signbit(fCos)) { // change rotate direction to clockwise - fSin = -fSin; - fCos = -fCos; - } - if (fCos < 0.5f) // > 60 - return; - TestRotate(x1, y1, fCos, fSin, "\t1: "); - TestRotate(x2, y2, fCos, fSin, "\t2: "); - std::cout << "\t(" << fCos << ',' << fSin << ')' << std::flush << std::endl; - } - void Rotate() - { - if (_lstPadClusters.size() < CLUSTER_PAD_MAX) // short pad base: do not rotate of axes - return; - // for( std::list<const PadCluster*>::const_iterator i = _lstPadClusters.begin(); i != _lstPadClusters.end(); i++ - // ) - // calcStat(*i); - const PadCluster *pFront = _lstPadClusters.front(); - int nFrontTime = pFront->getAdcHitMax().getTimeBin(); - int nFrontPad = pFront->getPad(); - const PadCluster *pBack = _lstPadClusters.back(); - int nBackTime = pBack->getAdcHitMax().getTimeBin(); - int nBackPad = pBack->getPad(); - int x21 = nBackTime - nFrontTime; - int y21 = nBackPad - nFrontPad; - - float f = sqrt(x21 * x21 + y21 * y21); - float fCos = x21 / f; - float fSin = y21 / f; - if (std::signbit(fCos)) { // change rotate direction to clockwise - fCos = -fCos; - fSin = -fSin; - }; - if (fCos < 0.5f) // if the angle > 60: do not touch - return; - // int a = _nTimeBinMin; // axis x shift - // int b = std::signbit(fSin) ? _nPadMax : _nPadMin; // axis y shift - int nXMin = INT_MAX; - int nXMax = 0; - int nYMin = INT_MAX; - int nYMax = 0; - std::list<PadCluster *> lstPadClusters; - PadCluster *pXYPadCluster = NULL; - for (std::list<const PadCluster *>::const_iterator i = _lstPadClusters.begin(); i != _lstPadClusters.end(); i++) { - const PadCluster *pPadCluster = *i; - int nPad = pPadCluster->getPad(); - float fYCos = nPad * fCos; - float fYSin = nPad * fSin; - size_t nSize = pPadCluster->getSize(); - for (uint n = 0; n < nSize; n++) { - const AdcHit &hit = pPadCluster->operator[](n); - float fAdc = hit.getAdc(); - int nTimeBin = hit.getTimeBin(); - float fXCos = nTimeBin * fCos; - float fXSin = nTimeBin * fSin; - int nX = round(fXCos + fYSin); - int nY = round(fYCos - fXSin); - AdcHit adcHit(nX, fAdc); - if (pXYPadCluster && pXYPadCluster->getPad() == nY) - pXYPadCluster->Add(adcHit); - else { - pXYPadCluster = new PadCluster(nY, adcHit); - lstPadClusters.push_back(pXYPadCluster); - } - // pRowClusters->Load(nX, nY, hit.getAdc()); - // std::cout << nTimeBin << '>' << nX << ':' << nY << ' '; - std::cout << nX << ' ' << nY << ' ' << fAdc << std::endl; - // std::cout << hit.getTimeBin() << ' ' << pPadCluster->getPad() << ' ' << nX << ' ' << nY << ' ' << fAdc - // << std::endl; - if (nX < nXMin) - nXMin = nX; - else if (nX > nXMax) - nXMax = nX; - if (nY < nYMin) - nYMin = nY; - else if (nY > nYMax) - nYMax = nY; - } - std::cout << std::flush << std::endl; - } - std::cout << "front{" << nFrontTime << ',' << nFrontPad << "} back:{" << nBackTime << ',' << nBackPad - << "} shift:{" << fCos << ',' << fSin << "} x:{" << nXMin << ',' << nXMax << "} y:{" << nYMin << ',' - << nYMax << '}' << std::flush << std::endl; - } - Cluster *Split(const Cluster *pPreviousCluster) - { // split current Cluster on new ones depends on AdcMax values + new array of ADC max values - bool bAdcMin = false; - std::list<const PadCluster *>::const_iterator iAdcMin; - std::list<const PadCluster *>::const_iterator iAdcMax; - std::list<const PadCluster *>::iterator i = _lstPadClusters.begin(); - const PadCluster *pPadCluster = *i; - calcStat(pPadCluster); - _AdcMaxTimeRange.Push_back(pPadCluster->getAdcHitMax().getTimeBin(), pPadCluster->getAdcMax()); - _AdcMaxPadRange.Push_back(pPadCluster->getPad(), pPadCluster->getAdcMax()); - iAdcMin = iAdcMax = i++; - for (; i != _lstPadClusters.end(); i++) { - pPadCluster = *i; - float fAdc = pPadCluster->getAdcMax(); - if (fAdc > (*iAdcMin)->getAdcMax() && bAdcMin) { - Cluster *pCluster = new Cluster(_nTimeBinMin, _nTimeBinMax, _nPadMin, _nPadMax); - pCluster->_AdcMaxTimeRange = _AdcMaxTimeRange; - pCluster->_AdcMaxPadRange = _AdcMaxPadRange; - pCluster->_AdcSumTimeRange = _AdcSumTimeRange; - pCluster->_AdcSumPadRange = _AdcSumPadRange; - - float fAdcMin = (*iAdcMin)->getAdcMax(); - float fLeft = (*iAdcMax)->getAdcMax() - fAdcMin; - assert(fLeft >= 0); - float fRight = fAdc - fAdcMin; - assert(fRight >= 0); - float f = fLeft + fRight; - fLeft /= f; // left aspect ratio - pCluster->_AdcMaxTimeRange.CorrectBackValue(fLeft); - pCluster->_AdcMaxTimeGauss.operator=(pCluster->_AdcMaxTimeRange.CalcGaussRes()); - pCluster->_AdcMaxPadRange.CorrectBackValue(fLeft); - pCluster->_AdcMaxPadGauss.operator=(pCluster->_AdcMaxPadRange.CalcGaussRes()); - - pCluster->_AdcSumTimeRange.setBackAspectRatio(fLeft); - pCluster->_AdcSumPadRange.setBackAspectRatio(fLeft); - - resetStat(); - calcStat(*iAdcMin); - - fRight /= f; // right aspect ratio - _AdcMaxTimeRange.Reset(); - _AdcMaxTimeRange.Push_back((*iAdcMin)->getAdcHitMax().getTimeBin(), fAdcMin * fRight); - _AdcMaxPadRange.Reset(); - _AdcMaxPadRange.Push_back((*iAdcMin)->getPad(), fAdcMin * fRight); - - _AdcSumTimeRange.setFrontAspectRatio(fRight); - _AdcSumPadRange.setFrontAspectRatio(fRight); - - pCluster->_lstPadClusters.splice(pCluster->_lstPadClusters.end(), _lstPadClusters, _lstPadClusters.begin(), - i); - std::for_each(pCluster->_lstPadClusters.begin(), pCluster->_lstPadClusters.end(), - [pCluster](const PadCluster *p) { ((PadCluster *)p)->setCluster(pCluster); }); - - pCluster->loadAdcSum(pPreviousCluster); - pCluster->_AdcSumTimeRange - .CorrectEdgeValues(); // update ADC min value for the combinet cluster by scaling factor - pCluster->_AdcSumTimeGauss.operator=(pCluster->_AdcSumTimeRange.CalcGaussRes()); - pCluster->_AdcSumPadRange - .CorrectEdgeValues(); // update ADC min value for the combinet cluster by scaling factor - pCluster->_AdcSumPadGauss.operator=(pCluster->_AdcSumPadRange.CalcGaussRes()); - - // _AdcMaxRange.Push_back( (*i)->getPad(), fAdc ); // ?! - return pCluster; - } else if (fAdc < (*iAdcMin)->getAdcMax()) - bAdcMin = true; - calcStat(pPadCluster); - _AdcMaxTimeRange.Push_back(pPadCluster->getAdcHitMax().getTimeBin(), fAdc); - _AdcMaxPadRange.Push_back(pPadCluster->getPad(), fAdc); - iAdcMax = iAdcMin; - iAdcMin = i; - } - _AdcMaxTimeGauss.operator=(_AdcMaxTimeRange.CalcGaussRes()); - _AdcMaxPadGauss.operator=(_AdcMaxPadRange.CalcGaussRes()); - loadAdcSum(pPreviousCluster); - _AdcSumTimeGauss.operator=(_AdcSumTimeRange.CalcGaussRes()); - _AdcSumPadGauss.operator=(_AdcSumPadRange.CalcGaussRes()); - return NULL; - } -#ifdef GSL_WAVELET - void WaveletTransform() - { - uint nTimeBins = _nTimeBinMax - _nTimeBinMin + 1; - uint nPads = _nPadMax - _nPadMin + 1; - int k = 1; - bool bPrint = 0; - while ((nTimeBins >= pow(2, k)) or (nPads >= pow(2, k))) k++; - int nSize = pow(2, k); - std::vector<double> abscoeff(nSize * nSize); - std::vector<size_t> p(nSize * nSize); - // std::cout << "Time: " << nTimeBins << " Pads: " << nPads << " k: " << k << " tmin: " << _nTimeBinMin << " pmin - // " << _nPadMin << std::endl; - gsl_matrix *data = gsl_matrix_calloc(nSize, nSize); - for (std::list<const PadCluster *>::iterator it1 = _lstPadClusters.begin(); it1 != _lstPadClusters.end(); it1++) { - PadCluster *pPadCluster = (PadCluster *)(*it1); - for (std::vector<AdcHit>::iterator it2 = (pPadCluster->_vAdcHits).begin(); - it2 != (pPadCluster->_vAdcHits).end(); it2++) { - gsl_matrix_set(data, (pPadCluster->getPad() - _nPadMin), (it2->getTimeBin() - _nTimeBinMin), it2->getAdc()); - // std::cout << "PadCl: " << (pPadCluster->getPad() - _nPadMin) << " : " << pPadCluster->getPad() << ", - // Time: " << (it2->getTimeBin()-_nTimeBinMin) << " : " << it2->getTimeBin() << std::endl; - } - } - if (bPrint) { - std::cout << "MATRIX:" << std::endl; - for (int i = 0; i < nSize; i++) { - for (int j = 0; j < nSize; j++) std::cout << gsl_matrix_get(data, i, j) << "\t"; - std::cout << std::endl; - } - } - // std::cout << std::endl; - gsl_wavelet *w = gsl_wavelet_alloc(gsl_wavelet_daubechies_centered, 4); - // gsl_wavelet * w = gsl_wavelet_alloc(gsl_wavelet_bspline_centered, 204); - gsl_wavelet_workspace *work = gsl_wavelet_workspace_alloc(data->tda); - gsl_wavelet2d_transform_matrix_forward(w, data, work); - - for (int i = 0; i < nSize; i++) - for (int j = 0; j < nSize; j++) abscoeff[i * nSize + j] = fabs(gsl_matrix_get(data, i, j)); - gsl_sort_index(p.data(), abscoeff.data(), 1, nSize * nSize); - int nc = 10; // count of largest components of the wavelet transform - for (int ip = 0; (ip + nc) < nSize * nSize; ip++) { - int i = p[ip] % nSize; - int j = p[ip] / nSize; - gsl_matrix_set(data, i, j, 0.); - } - - gsl_wavelet2d_transform_matrix_inverse(w, data, work); - if (bPrint) { - std::cout << "MATRIX+:" << std::endl; - for (int i = 0; i < nSize; i++) { - for (int j = 0; j < nSize; j++) { - if (abs(gsl_matrix_get(data, i, j)) > 0.01) - std::cout << gsl_matrix_get(data, i, j) << "\t"; - else - std::cout << 0 << "\t"; - } - std::cout << std::endl; - } - } - for (std::list<const PadCluster *>::iterator it1 = _lstPadClusters.begin(); it1 != _lstPadClusters.end(); it1++) { - PadCluster *pPadCluster = (PadCluster *)(*it1); - for (std::vector<AdcHit>::iterator it2 = (pPadCluster->_vAdcHits).begin(); - it2 != (pPadCluster->_vAdcHits).end(); it2++) { - it2->setAdc(gsl_matrix_get(data, (pPadCluster->getPad() - _nPadMin), (it2->getTimeBin() - _nTimeBinMin))); - } - pPadCluster->ReCalcStats(); - } - // for( std::list<const PadCluster*>::iterator it1 = _lstPadClusters.begin(); it1 != _lstPadClusters.end(); it1++ - // ) { - // std::vector<AdcHit> vnAdcHits = (*it1)->_vAdcHits; - // for( std::vector<AdcHit>::iterator it2 = vnAdcHits.begin(); it2 != vnAdcHits.end(); it2++ ) - // std::cout << it2->getAdc() << " "; - // } - // std::cout << std::endl; - } -#endif - std::ostringstream &Draw(std::ostringstream &oss) const - { - ClusterHit::Print(oss); - int nTable = (nADC_DISCRECITY + 3) * 2 + _nTimeBinMax - _nTimeBinMin + 8 + 4; - oss << std::string(nTable, '=') << std::endl; - oss << "| pads |" << std::string(_nTimeBinMax - _nTimeBinMin + 3, ' ') << "| " << _RED_MSG_("adcMax:") - << std::string(nADC_DISCRECITY - 6, ' ') << "| " << "adcSum:" << std::string(nADC_DISCRECITY - 6, ' ') << '|' - << std::endl; - oss << std::string(nTable, '=') << std::endl; - ClusterStat stat; - stat.Calc(_lstPadClusters); - for (std::list<const PadCluster *>::const_iterator i = _lstPadClusters.begin(); i != _lstPadClusters.end(); i++) { - const PadCluster *pPadCluster = *i; - if (i != _lstPadClusters.begin()) oss << std::endl; - oss << '|' << std::setw(4) << pPadCluster->getPad() << ": | "; - - pPadCluster->Draw(oss, _nTimeBinMin, _nTimeBinMax, this); - oss << " |"; - - float fAdc = pPadCluster->getAdcHitMax().getAdc(); - long nAdc = 0; - if (stat._fAdcMax > stat._fAdcMin) { - float f = (stat._fAdcMax - stat._fAdcMin) / nADC_DISCRECITY; - nAdc = lroundf((fAdc - stat._fAdcMin) / f); - } else - nAdc = nADC_DISCRECITY / 2; - assert(nAdc <= nADC_DISCRECITY); - - if (nAdc) oss << std::string(nAdc, '-'); - oss << (fAdc == stat._fAdcMax ? '>' : '+'); - uint nSpace = nADC_DISCRECITY - nAdc + 1; - if (nSpace) oss << std::string(nSpace, ' '); - oss << '|'; - - float fAdcSum = pPadCluster->getAdcHitSum(); - long nAdcSum = 0; - if (stat._fSumMax > stat._fSumMin) { - float f = (stat._fSumMax - stat._fSumMin) / nADC_DISCRECITY; - nAdcSum = lroundf((fAdcSum - stat._fSumMin) / f); - } else - nAdcSum = nADC_DISCRECITY / 2; - assert(nAdcSum <= nADC_DISCRECITY); - - if (nAdcSum) oss << std::string(nAdcSum, '-'); - oss << (fAdcSum == stat._fSumMax ? '>' : '+'); - - nSpace = nADC_DISCRECITY - nAdcSum + 1; - if (nSpace) oss << std::string(nSpace, ' '); - oss << '|'; - } - oss << std::endl << std::string(nTable, '=') << std::endl; - ; - const uint nTAB = 4; - oss << std::string(nTAB, ' ') << "Cluster region:{" << " timebin:{ min:" << _nTimeBinMin - << " max:" << _nTimeBinMax << " } pad:{ min:" << _nPadMin << " max:" << _nPadMax << " } }" << std::endl; - - _AdcMaxPadRange.Draw(oss, " pad", "AdcMax"); - _AdcMaxPadRange.Print(oss, "pad", "AdcMax"); - - _AdcMaxTimeRange.Draw(oss, "time", "AdcMax"); - _AdcMaxTimeRange.Print(oss, "time", "AdcMax"); - - _AdcSumPadRange.Draw(oss, " pad", "AdcSum"); - _AdcSumPadRange.Print(oss, "pad", "AdcSum"); - - _AdcSumTimeRange.Draw(oss, "time", "AdcSum"); - _AdcSumTimeRange.Print(oss, "time", "AdcSum"); - return oss; - } - std::ostringstream &Print(std::ostringstream &oss, const std::list<const PadCluster *> &r) const - { - const uint nTAB = 4; - if (!r.empty()) { - oss << std::string(nTAB, ' ') << "======== "; - for (std::list<const PadCluster *>::const_iterator i = r.begin(); i != r.end(); i++) { - if (i != r.begin()) oss << std::string(nTAB, ' ') << "-------- "; - if (&r == &_lstPadClusters) - (*i)->PrintPad(oss); - else - (*i)->PrintTime(oss); - } - } - return oss; - } - std::ostringstream &Print(std::ostringstream &oss) const - { - ClusterHit::Print(oss); - const uint nTAB = 4; - oss << std::string(nTAB, ' ') << "Cluster region:{" << " timebin:{ min:" << _nTimeBinMin - << " max:" << _nTimeBinMax << " } pad:{ min:" << _nPadMin << " max:" << _nPadMax << " } }" << std::endl; - _AdcMaxPadRange.Print(oss, "pad", "AdcMax"); - _AdcMaxTimeRange.Print(oss, "time", "AdcMax"); - _AdcSumPadRange.Print(oss, "pad", "AdcSum"); - _AdcSumTimeRange.Print(oss, "time", "AdcSum"); - Print(oss, _lstPadClusters); - return oss; - } - void DrawPadClusters(std::ostringstream &oss) const - { - for (std::list<const PadCluster *>::const_iterator i = _lstPadClusters.begin(); i != _lstPadClusters.end(); i++) { - const PadCluster *pPadCluster = (*i); - oss << std::setw(3) << pPadCluster->getPad() << ":"; - uint nPos = 0; - uint nMin = pPadCluster->getTimeBinMin(); - if (nMin < nPos) { // shift back to the 'n' positions (for combined PadClusters!!!) - long int nShift = nPos - nMin; - oss.seekp(oss.tellp() - nShift); - nPos -= nShift; - } - uint nGap = nMin - nPos; - if (nGap) { - oss << std::string(nGap, '-'); - nPos += nGap; - } - char ch = SZ_ALPHABET[getId() % strlen(SZ_ALPHABET)]; - pPadCluster->Draw(oss, this, ch); - nPos += pPadCluster->getSize(); - oss << std::endl; - } - } -}; -uint PadCluster::getClusterId(const Cluster *pCluster) const -{ - return pCluster == NULL ? 0 : ((const ClusterHit *)pCluster)->getId(); -} -float PadCluster::getClusterTime(const Cluster *pCluster) const -{ - return pCluster == NULL ? 0 : ((const ClusterHit *)pCluster)->getTimeHitMean(); -} -float PadCluster::getClusterPad(const Cluster *pCluster) const -{ - return pCluster == NULL ? 0 : ((const ClusterHit *)pCluster)->getPadHitMean(); -} -} // namespace tpcClustering - -#endif diff --git a/detectors/tpc/clusterHitFinder/fast/TpcClusterHit.h b/detectors/tpc/clusterHitFinder/fast/TpcClusterHit.h deleted file mode 100644 index c26afff8..00000000 --- a/detectors/tpc/clusterHitFinder/fast/TpcClusterHit.h +++ /dev/null @@ -1,52 +0,0 @@ -// -// Joint Institute for Nuclear Research (JINR), Dubna, Russia, Jun-2024 -// MPD TPC Clustering objects for NICA project -// author: Viktor A.Krylov JINR/LNP -// email: kryman@jinr.ru -// -#ifndef CLUSTERHIT_H -#define CLUSTERHIT_H - -#include "TpcClusterUnits.h" - -namespace tpcClustering { -class ClusterHit { - uint _nId; // cluster id -public: - Gauss _AdcMaxTimeGauss; - Gauss _AdcMaxPadGauss; - Gauss _AdcSumTimeGauss; - Gauss _AdcSumPadGauss; - - ClusterHit() : _nId(0) {} - ClusterHit(uint nId) : _nId(nId) {} - ~ClusterHit() {} - inline void setId(uint nId) { _nId = nId; } - inline uint getId() const { return _nId; } - inline float getPadHitMean() const { return _AdcMaxPadGauss.getMean(); } - inline float getTimeHitMean() const { return _AdcMaxTimeGauss.getMean(); } - std::ostringstream &Print(std::ostringstream &oss) const - { - char ch = SZ_ALPHABET[getId() % strlen(SZ_ALPHABET)]; - oss << "Cluster:: id:{ " << _nId << " '" << ch << "' }" << std::endl; - uint nTAB = 4; - oss << std::string(nTAB, ' ') << "Cluster hit:" << std::endl; - - oss << std::string(nTAB * 2, ' ') << "Time plane gauss:"; - _AdcMaxTimeGauss.Print(oss); - oss << std::endl; - oss << std::string(nTAB * 2, ' ') << "Pad plane gauss:"; - _AdcMaxPadGauss.Print(oss); - oss << std::endl; - - oss << std::string(nTAB * 2, ' ') << "Time plane sum gauss:"; - _AdcSumTimeGauss.Print(oss); - oss << std::endl; - oss << std::string(nTAB * 2, ' ') << "Pad plane sum gauss:"; - _AdcSumPadGauss.Print(oss); - oss << std::endl; - return oss; - } -}; -} // namespace tpcClustering -#endif diff --git a/detectors/tpc/clusterHitFinder/fast/TpcClusterHitFinderFast.cxx b/detectors/tpc/clusterHitFinder/fast/TpcClusterHitFinderFast.cxx index 09f85977..4cd1d613 100644 --- a/detectors/tpc/clusterHitFinder/fast/TpcClusterHitFinderFast.cxx +++ b/detectors/tpc/clusterHitFinder/fast/TpcClusterHitFinderFast.cxx @@ -1,5 +1,7 @@ #include "TpcClusterHitFinderFast.h" #include "Tpc2dClusterFast.h" +#include "clustering/TpcClusters.h" +#include "TpcSectorObj.h" #include "MpdTpcDigit.h" #include "MpdTpcHit.h" @@ -9,14 +11,13 @@ #include "FairLogger.h" // C/C++ Headers ---------------------- -#include <iostream> #include <math.h> #include <vector> #define fV_DRIFT 0.0055 // cm/ns using namespace std; -using namespace tpcClustering; +using namespace TpcClustering; //__________________________________________________________________________ @@ -26,45 +27,80 @@ void TpcClusterHitFinderFast::FindHits() LOG(info) << "TpcClusterHitFinderFast::Exec started: Event " << nEvent; - /* - void EventClusters::Load(uint nSector, uint nRow, uint nPad, uint nTimeBin, float fAdc, int nTrackId) - - initial loading of digits - - partial sorting of digits to pre-clusters - - charge correction at the boundary between pre-clusters - - void EventClusters::DefineClusters(std::ostringstream& debugInfo) - - converting pre-clusters into clusters, where they are grouped with neighboring digits - until cluster closure - - on cluster closure the local coordinate of the Hit is evaluated - */ - std::ostringstream oss; - EventClusters *pEventClusters = new EventClusters(nEvent); - - int nDigits = digiArray->GetEntriesFast(); - for (int i = 0; i < nDigits; ++i) { - MpdTpcDigit *pdigit = (MpdTpcDigit *)digiArray->UncheckedAt(i); // Input digit - pEventClusters->Load(pdigit->GetSector(), pdigit->GetRow(), pdigit->GetPad(), pdigit->GetTimeBin(), - pdigit->GetAdc(), pdigit->GetOrigin()); - } + // Load digit collection separately for each sector + + TpcEventObj tpcEvent(nEvent); + int startSector = ((MpdTpcDigit *)digiArray->UncheckedAt(0))->GetSector(); + int nextSector = startSector; + int totalSectors = 24; + int startDigit = 0; + int nDigits = digiArray->GetEntriesFast(); + + for (int iSec = startSector; iSec < totalSectors; ++iSec) { + if (iSec != nextSector) continue; + + TpcSectorObj *pTpcSector = new TpcSectorObj(iSec); + pTpcSector->_pTpcClusters = new TpcClusters(iSec); + for (int i = startDigit; i < nDigits; ++i) { + MpdTpcDigit *pdigit = (MpdTpcDigit *)digiArray->UncheckedAt(i); // Input digit + if (pdigit->GetSector() != iSec) { + nextSector = pdigit->GetSector(); + startDigit = i; + break; + } + pTpcSector->_pTpcClusters->Load(pdigit->GetRow(), pdigit->GetPad(), pdigit->GetTimeBin(), pdigit->GetAdc()); + } - if (!pEventClusters->isEmpty()) { - pEventClusters->DefineClusters(oss); - TransformOutput(pEventClusters); + tpcEvent._vpTpcSectors.push_back(std::move(pTpcSector)); + delete pTpcSector; } - delete pEventClusters; + // Perform calculations (clusters & hits) + tpcEvent.Calc(); + + // Extract & transform hits to proper representation + TransformOutput(tpcEvent._vpTpcSectors); } -void TpcClusterHitFinderFast::TransformOutput(const EventClusters *pEventClusters) +void TpcClusterHitFinderFast::TransformOutput(const std::vector<TpcSectorObj *> &tpcSectors) { + int padID; + int nHits = 0; + for (int iSec = 0; iSec < tpcSectors.size(); ++iSec) { + int currentSector = tpcSectors[iSec]->_nSector; + for (int iHit = 0; iHit < tpcSectors[iSec]->_vTpcAdcHits.size(); ++iHit) { + const auto ¤tHit = tpcSectors[iSec]->_vTpcAdcHits[iHit]; + padID = (currentSector % secGeo->GetSectorCountHalf()) | (currentHit._nRow << 5); + TVector3 p3glob(currentHit._fX, currentHit._fY, currentHit._fZ); + + MpdTpcHit *hit = new ((*hitArray)[nHits]) MpdTpcHit(padID, p3glob, TVector3(0.05, 0.0, 0.2), nHits); + + hit->SetLayer(currentHit._nRow); + hit->SetSector(currentSector); + // hit->SetLocalPosition(p3loc); + hit->SetEnergyLoss(currentHit._fAdc); + hit->SetModular(1); + // hit->SetPad(int(fPad)); + // hit->SetBin(int(fTime)); + // hit->SetPadCoordinate(fPad); + // hit->SetTimeBinCoordinate(fTime); + // hit->SetNdigits(rlstPadClusters.size()); + // hit->SetStep(secGeo->GetPadHeight()[currentPadArea]); + + ++nHits; + } + } + /* - store clusters in clusArray - transform ClusterHit from local to global coordinate - store transformed hits in hitArray */ - int nClusters = 0; + // hits are tpcEvent._vpTpcSectors->_vTpcAdcHits + + /*int nClusters = 0; int nHits = 0; const list<SectorClusters *> &rlstEventClusters = pEventClusters->_lstSectorClusters; @@ -119,7 +155,7 @@ void TpcClusterHitFinderFast::TransformOutput(const EventClusters *pEventCluster double phSec = nSector * secGeo->GetSectorPhiRad(); p3glob.RotateZ(phSec); - /* write hit(s) */ + // write hit(s) int padID = (nSector % secGeo->GetSectorCountHalf()) | (nRow << 5); MpdTpcHit *hit = new ((*hitArray)[nHits]) MpdTpcHit(padID, p3glob, TVector3(0.05, 0.0, 0.2), nClusters - 1); hit->SetLayer(nRow); @@ -137,6 +173,6 @@ void TpcClusterHitFinderFast::TransformOutput(const EventClusters *pEventCluster ++nHits; } } - } + }*/ } ClassImp(TpcClusterHitFinderFast); diff --git a/detectors/tpc/clusterHitFinder/fast/TpcClusterHitFinderFast.h b/detectors/tpc/clusterHitFinder/fast/TpcClusterHitFinderFast.h index 9b09012e..12489eab 100644 --- a/detectors/tpc/clusterHitFinder/fast/TpcClusterHitFinderFast.h +++ b/detectors/tpc/clusterHitFinder/fast/TpcClusterHitFinderFast.h @@ -11,7 +11,8 @@ #include "AbstractTpcClusterHitFinder.h" // Fast ClusterHitFinder library -#include "TpcEventClusters.h" +#include "TpcEventObj.h" +using namespace TpcClustering; // Collaborating Class Declarations -- class TClonesArray; @@ -32,7 +33,7 @@ public: void FindHits(); private: - void TransformOutput(const tpcClustering::EventClusters *pEventClusters); + void TransformOutput(const std::vector<TpcSectorObj *> &tpcSectors); ClassDef(TpcClusterHitFinderFast, 1); }; diff --git a/detectors/tpc/clusterHitFinder/fast/TpcClusterUnits.h b/detectors/tpc/clusterHitFinder/fast/TpcClusterUnits.h deleted file mode 100644 index 94339108..00000000 --- a/detectors/tpc/clusterHitFinder/fast/TpcClusterUnits.h +++ /dev/null @@ -1,683 +0,0 @@ -// -// Joint Institute for Nuclear Research (JINR), Dubna, Russia, Nov-2021/May-2023 -// MPD TPC Clustering objects for NICA project -// author: Viktor A.Krylov JINR/LNP -// email: kryman@jinr.ru -// -#ifndef TPC_CLUSTERING_H -#define TPC_CLUSTERING_H - -#include <stdlib.h> -#include <iomanip> -#include <iostream> -#include <sstream> -#include <ostream> -#include <cassert> -#include <math.h> -#include <string> -#include <cstring> -#include <vector> -#include <array> -#include <list> -#include <algorithm> -#include <iterator> -#include <climits> -#include <float.h> -// #include <gsl/gsl_sort.h> -// #include <gsl/gsl_wavelet2d.h> -// #include <gsl/gsl_matrix.h> -#include <pthread.h> - -#include "TpcPadCluster.h" -#include "TpcAdcHit.h" -#include "TpcPadHit.h" - -#define HIT_FINDER "v0.2.0b" -#define AUXHIT_DEBUG -#define HIT_ADJUSTMENT -#define nTPC_SECTORS 24 -#define CLUSTER_PAD_MAX 5 -#undef DEBUG_OUTPUT -#undef GSL_WAVELET - -// const float fADC_DOORSTEP = 2.5; // adc sensitivity - -// BASH COLOR CODES: {black:30,red:31,green:32,yellow:33,blue:34,magenta:35,cyan:36,white:37}; background:+10; -// bright:+60 BASH TEXT DECORATION: {regular:0,bold:1,italic:3,underscore:4} BASH CODE STRUCTURE: <PREFIX>;<COLOR>;<TEXT -// DECORATION> -#define _RED_ "0;31" -#define _GREEN_ "0;32" -#define _YELLOW_ "0;33" -#define _BASH_CODE_(_CODE) "\e[" _CODE "m" -#define _RED_MSG_(_MSG) _BASH_CODE_(_RED_) _MSG _BASH_CODE_() -// -// A Fast, Accurate, and Separable Method for Fitting a Gaussian Function -// Ibrahim Al-Nahhal, Octavia A. Dobre, Ertugrul Basar, Cecilia Moloney, Salama Ikki -// DOI: 10.1109/MSP.2019.2927685 -// arXiv:1907.07241v1 -// -// GUO’S algorithm without iterative procedure -// Solution of AX = B -> {A = LU} -> LUX = B -> {UX = Y} -> LY = B -> Y -// UX = Y -> X -// | Sum(y^2) Sum(x*y^2) Sum(x^2*y^2) | -// MATRIX A = | Sum(x*y^2) Sum(x^2*y^2) Sum(x^3*y^2) | -// | Sum(x^2*y^2) Sum(x^3*y^2) Sum(x^4*y^2) | -// -// | Sum(y^2*Log(y)) | -// MATRIX B = | Sum(x*y^2*Log(y)) | -// | Sum(x^2*y^2*Log(y)) | -// -// | l11 0 0 | | l11 l21 l31 | -// L = | l21 l22 0 | U = | 0 l22 l32 | -// | l31 l32 l33 | | 0 0 l33 | -// -template <typename X, typename Y> -struct XY { - X x; - Y y; -}; -template <typename X, typename Y> -double calcMatrix(const std::vector<XY<X, Y>> &rvXY, std::array<double, 5> &radSum, std::array<double, 3> &radSumLog) -{ - double dSum0 = 0, dSum1 = 0, dSum2 = 0, dSum3 = 0, dSum4 = 0; - double dSumLog0 = 0, dSumLog1 = 0, dSumLog2 = 0; - double yMax = 0; - X xMax = 0, xMin = 0; - double dGravityAxis = 0, dGravitySum = 0; - size_t nSize = rvXY.size(); - for (uint i = 0; i < nSize; i++) { - X x = rvXY[i].x; - X x2 = x * x; - X x3 = x2 * x; - X x4 = x3 * x; - if (!i || x < xMin) xMin = x; - - double y = rvXY[i].y; - if (y > yMax) { - yMax = y; - xMax = x; - } - dGravitySum += y; - dGravityAxis += x * y; - double y2 = y * y; - - dSum0 += y2; - dSum1 += x * y2; - dSum2 += x2 * y2; - dSum3 += x3 * y2; - dSum4 += x4 * y2; - - double dYLog = y2 * log(y); - dSumLog0 += dYLog; - dSumLog1 += x * dYLog; - dSumLog2 += x2 * dYLog; - } - radSum = {dSum0, dSum1, dSum2, dSum3, dSum4}; - radSumLog = {dSumLog0, dSumLog1, dSumLog2}; - dGravityAxis = dGravitySum > 0.0 ? dGravityAxis / dGravitySum : 0.0; - return dGravityAxis <= xMin ? xMax : dGravityAxis; -} -const double dSQRT_2_PI = sqrt(M_PI_2); -inline double calcGaussArea(double x, double y, double d, double xMin, double xMax) -{ // Gauss area from xMin to xMax - double d2 = M_SQRT2 * d; // deviation - double fErfMin = std::erf((xMin - x) / d2); - double fErfMax = std::erf((xMax - x) / d2); - return dSQRT_2_PI * y * d * (fErfMax - fErfMin); -} -template <typename X, typename Y> -std::array<double, 4> calcGaussRes(const std::vector<XY<X, Y>> &rvXY) -{ - assert(rvXY.size() > 2); - - std::array<double, 5> adSum; - std::array<double, 3> adSumLog; - double xMax = calcMatrix<X, Y>(rvXY, adSum, adSumLog); - - double dL11 = sqrt(adSum[0]); - double dL21 = adSum[1] / dL11; - double dL22 = sqrt(adSum[2] - dL21 * dL21); - double dL31 = adSum[2] / dL11; - double dL32 = (adSum[3] - dL21 * dL31) / dL22; - double dL33 = sqrt(adSum[4] - dL31 * dL31 - dL32 * dL32); - - double dY1 = adSumLog[0] / dL11; - double dY2 = (adSumLog[1] - dL21 * dY1) / dL22; - double dY3 = (adSumLog[2] - dL31 * dY1 - dL32 * dY2) / dL33; - - double dX2 = dY3 / dL33; - double dX1 = (dY2 - dL32 * dX2) / dL22; - double dX0 = (dY1 - dL31 * dX2 - dL21 * dX1) / dL11; - - double dX25 = -0.5 / dX2; - double x = dX1 * dX25; - double y = exp(dX0 + 0.5 * dX1 * x); - double d = dX25 > 0 ? sqrt(dX25) : 0; // deviation - return {x, y, d, xMax}; -} -template <typename X, typename Y> -double calcRsquared(const std::vector<XY<X, Y>> &rvXY, double x, double y, double d, double dYSum) -{ // coefficient of determination - size_t nSize = rvXY.size(); - if (nSize > 2) { // deviation should be defined!!! - double dYMean = dYSum / nSize; - double d2Deviation2 = 2 * d * d; - double dSumTotal = 0; // Sum total = Sum[(y - ymean)^2] - double dSumRes = 0; // Sum of residuals = Sum[(y - f(x))^2] - for (uint i = 0; i < nSize; i++) { - double dX = rvXY[i].x - x; - double dGauss = y * exp(-dX * dX / d2Deviation2); - double dY = rvXY[i].y; - double dTot = dY - dYMean; - double dRes = dY - dGauss; - - dSumTotal += dTot * dTot; - dSumRes += dRes * dRes; - } - return 1 - dSumRes / dSumTotal; - } else - return 0; -} - -namespace tpcClustering { -struct Gauss { - std::array<double, 4> - _array; // array[0: mean value, 1: gaussian peak value (ADC), 2: standard deviation, 3: adjustment] - double _dArea; // gaussian area of the normal distribution (should be defined for reference only!) - double _dRsquared; // r-squared value of the normal distribution (should be defined for reference only!) -public: - Gauss() : _array{0, 0, 0, 0}, _dArea(0), _dRsquared(0) {} - Gauss(const Gauss &r) : _dArea(r._dArea), _dRsquared(r._dRsquared) { _array = r._array; } - Gauss &operator=(const Gauss &r) - { - if (this != &r) { - _array = r._array; - _dArea = r._dArea; - _dRsquared = r._dRsquared; - } - return *this; - } - void Set(const std::array<double, 4> &r) { _array = r; } - - void setMean(double d) { _array[0] = d; } - double getMean() const { return _array[0]; } - - void setMax(double d) { _array[1] = d; } - double getMax() const { return _array[1]; } - - void setDeviation(double d) { _array[2] = d; } - double getDeviation() const { return _array[2]; } - - double getAdjustment() const - { - return _array[3]; - } // adjustment value for unreasonable fit results to max ADC value!!! - - void setArea(double d) { _dArea = d; } - double getArea() const { return _dArea; } - - void setRsquared(double d) { _dRsquared = d; } - double getRsquared() const { return _dRsquared; } - - std::ostringstream &Print(std::ostringstream &oss) const - { - oss << "{ mean:" << getMean() << " max:" << getMax() << " deviation:" << getDeviation() - << " adjustment:" << getAdjustment() << " area:" << _dArea << " R-squared:" << _dRsquared << " }"; - return oss; - } -}; - -struct ClusterStat { - float _fAdcMin; - float _fAdcMax; - float _fSumMin; - float _fSumMax; - uint _nMin; - uint _nMax; - ClusterStat() : _fAdcMin(0), _fAdcMax(0), _fSumMin(0), _fSumMax(0), _nMin(UINT_MAX), _nMax(UINT_MAX) {} - void Reset() - { - _fAdcMin = _fAdcMax = _fSumMin = _fSumMax = 0; - _nMin = _nMax = UINT_MAX; - } - void Calc(const std::list<const PadCluster *> &rList) - { - Reset(); - for (std::list<const PadCluster *>::const_iterator i = rList.begin(); i != rList.end(); i++) { - const PadCluster *pPadCluster = *i; - float fAdcMin = pPadCluster->getAdcMin(); - float fAdcMax = pPadCluster->getAdcMax(); - if (_fAdcMin == 0 && _fAdcMax == 0) { - _fAdcMin = fAdcMin; - _fAdcMax = fAdcMax; - } else { - if (fAdcMin < _fAdcMin) _fAdcMin = fAdcMin; - if (fAdcMax > _fAdcMax) _fAdcMax = fAdcMax; - } - float fAdcSum = pPadCluster->getAdcHitSum(); - if (fAdcSum > 0) { - if (_fSumMin == 0 && _fSumMax == 0) - _fSumMin = _fSumMax = fAdcSum; - else if (fAdcSum < _fSumMin) - _fSumMin = fAdcSum; - else if (fAdcSum > _fSumMax) - _fSumMax = fAdcSum; - } - uint nTimeBinMin = pPadCluster->getTimeBinMin(); - uint nTimeBinMax = pPadCluster->getTimeBinMax(); - if (_nMin == UINT_MAX && _nMax == UINT_MAX) { - _nMin = nTimeBinMin; - _nMax = nTimeBinMax; - } else { - if (nTimeBinMin < _nMin) _nMin = nTimeBinMin; - if (nTimeBinMax > _nMax) _nMax = nTimeBinMax; - } - }; - } -}; -const uint nADC_DISCRECITY = 40; // should be more than 9 -struct ClusterBounds { - uint _nTimeBinMin; // timebin min range - uint _nTimeBinMax; // timebin max range - uint _nPadMin; // pad min range - uint _nPadMax; // pad max range - ClusterBounds() { Reset(); } - ClusterBounds &operator=(const ClusterBounds &r) - { - if (this != &r) { - _nTimeBinMin = r._nTimeBinMin; - _nTimeBinMax = r._nTimeBinMax; - _nPadMin = r._nPadMin; - _nPadMax = r._nPadMax; - } - return *this; - } - void Reset() { _nTimeBinMin = _nTimeBinMax = _nPadMin = _nPadMax = UINT_MAX; } - void SetTimeBin(uint nTimeBin) - { - if (_nTimeBinMin == UINT_MAX && _nTimeBinMax == UINT_MAX) - _nTimeBinMin = _nTimeBinMax = nTimeBin; - else if (nTimeBin < _nTimeBinMin) - _nTimeBinMin = nTimeBin; - else if (nTimeBin > _nTimeBinMax) - _nTimeBinMax = nTimeBin; - } - void SetPad(uint nPad) - { - if (_nPadMin == UINT_MAX && _nPadMax == UINT_MAX) - _nPadMin = _nPadMax = nPad; - else if (nPad < _nPadMin) - _nPadMin = nPad; - else if (nPad > _nPadMax) - _nPadMax = nPad; - } - void Set(const PadCluster *pPadCluster) - { - SetTimeBin(pPadCluster->getTimeBinMin()); - SetTimeBin(pPadCluster->getTimeBinMax()); - SetPad(pPadCluster->getPad()); - } - void Set(uint nTimeBin, uint nPad) - { - SetTimeBin(nTimeBin); - SetPad(nPad); - } -}; -template <typename X, typename Y> -struct HitRange { - X _xMin; // x min range - X _xMax; // x max range - Y _yMin; // y min range - Y _yMax; // y max range - Y _ySum; // y sum accumulator - float _fFront; // front aspect ratio - float _fBack; // back aspect ratio - std::vector<XY<X, Y>> _vXY; // full range hit values - HitRange() { Reset(); } - HitRange<X, Y> &operator=(const HitRange<X, Y> &r) - { - if (this != &r) { - _xMin = r._xMin; - _xMax = r._xMax; - _yMin = r._yMin; - _yMax = r._yMax; - _fFront = r._fFront; - _fBack = r._fBack; - _ySum = r._ySum; - _vXY.assign(r._vXY.begin(), r._vXY.end()); - } - return *this; - } - const XY<X, Y> &operator[](uint n) const - { - assert(n < _vXY.size()); - return _vXY[n]; - } - void Reset() - { - _xMin = _xMax = UINT_MAX; - _yMin = _yMax = _ySum = 0; - _fFront = _fBack = 1.0f; - _vXY.clear(); - } - void SetRange(X xMin, X xMax) - { - assert(xMax >= xMin); - _xMin = xMin; - _xMax = xMax; - _vXY.resize(xMax - xMin + 1); - } - void setFrontAspectRatio(float f) - { - assert(f <= 1.0f); - _fFront = f; - } - void setBackAspectRatio(float f) - { - assert(f <= 1.0f); - _fBack = f; - } - void CorrectFrontValue() - { - if (_fFront != 1.0f) { - Y y = _vXY.front().y * _fFront; - if (y < _yMin) _yMin = y; - _ySum -= _vXY.front().y - y; - _vXY.front().y = y; - } - } - void CorrectFrontValue(float f) - { - Y y = _vXY.front().y * f; - if (y < _yMin) _yMin = y; - _ySum -= _vXY.front().y - y; - _vXY.front().y = y; - } - void CorrectBackValue() - { - if (_fBack != 1.0f) { - Y y = _vXY.back().y * _fBack; - if (y < _yMin) _yMin = y; - _ySum -= _vXY.back().y - y; - _vXY.back().y = y; - } - } - void CorrectBackValue(float f) - { - Y y = _vXY.back().y * f; - if (y < _yMin) _yMin = y; - _ySum -= _vXY.back().y - y; - _vXY.back().y = y; - } - void CorrectEdgeValues() - { - CorrectFrontValue(); - CorrectBackValue(); - } - void Add(X x, Y y) - { - assert(x >= _xMin && x <= _xMax); - X n = x - _xMin; - assert(n < _vXY.size()); - - XY<X, Y> &r = _vXY[n]; - if (r.x == 0 && r.y == 0) { // define 1st hit - r.x = x; - r.y = y; - } else { - assert(r.x == x); - r.y += y; - } - if (r.y > _yMax) _yMax = r.y; - } - void AddTimeHitAdc(const PadCluster *pPadCluster) - { - assert(pPadCluster); - size_t nSize = pPadCluster->getSize(); - for (uint i = 0; i < nSize; i++) { - const AdcHit &rAdcHit = pPadCluster->operator[](i); - Add(rAdcHit.getTimeBin(), rAdcHit.getAdc()); - } - } - void AddPadHitAdcSum(const PadCluster *pPadCluster) - { - assert(pPadCluster); - Add(pPadCluster->getPad(), pPadCluster->getAdcHitSum()); - } - void Push_back(X x, Y y) - { - if (_xMin == UINT_MAX && _xMax == UINT_MAX) - _xMin = _xMax = x; - else if (x < _xMin) - _xMin = x; - else if (x > _xMax) - _xMax = x; - if (_yMin == 0 && _yMax == 0) - _yMin = _yMax = y; - else if (y < _yMin) - _yMin = y; - else if (y > _yMax) - _yMax = y; - _ySum += y; - - if (_vXY.size() && _vXY.back().x == x) { // is the ADC value belong to previous channel - if (y > _vXY.back().y) // take max value - _vXY.back().y = y; - } else - _vXY.push_back({.x = x, .y = y}); - } - // void CorrectBackMinValue( Y y ) { - // if( y < _yMin ) - // _yMin = y; - // _ySum -= _vXY.back().y - y; - // _vXY.back().y = y; - // } - // float CalcTimeBinMean( float fPadMean ) { - // float fMean = (fPadMean - _nPadMin) / (_nPadMax - _nPadMin); - // float f = (_nTimeBinMax - _nTimeBinMin) * fMean; - // return _n1stTimeBin > _nLastTimeBin ? _n1stTimeBin - f : _n1stTimeBin + f; - // } - Gauss CalcGaussRes() - { - Gauss gauss; - size_t nSize = _vXY.size(); - if (nSize > 2) { - if (nSize == _xMax - _xMin + 1) - gauss.Set(::calcGaussRes<X, Y>(_vXY)); - else { - std::vector<XY<X, Y>> vXY; - for (uint i = 0; i < nSize; i++) { - const XY<X, Y> &r = _vXY[i]; - if (r.x != 0 && r.y != 0) vXY.push_back(r); - } - if (vXY.size()) gauss.Set(::calcGaussRes<X, Y>(vXY)); - } -#ifdef AUXHIT_DEBUG - if (gauss.getDeviation() != 0) { - gauss.setArea(::calcGaussArea(gauss.getMean(), gauss.getMax(), gauss.getDeviation(), 0, nSize)); - gauss.setRsquared( - calcRsquared<X, Y>(_vXY, gauss.getMean(), gauss.getMax(), gauss.getDeviation(), double(_ySum))); - } else { // it cannot fit the Gauss - gauss.setArea(0); - gauss.setRsquared(0); - } -#endif -#ifdef HIT_ADJUSTMENT - double dPadMean = gauss.getMean(); - if (dPadMean < _xMin || dPadMean > _xMax || gauss.getDeviation() == 0) gauss.setMean(gauss.getAdjustment()); -#endif - } else if (nSize == 2) { - const XY<X, Y> &r0 = _vXY[0]; - const XY<X, Y> &r1 = _vXY[1]; - float fMean = (r0.x * r0.y + r1.x * r1.y) / (r0.y + r1.y); - gauss.setMean(fMean); // shift to real value - gauss.setMax(r0.y > r1.y ? r0.y : r1.y); - gauss.setDeviation(0); - gauss.setArea(0); - gauss.setRsquared(0); - } else if (nSize == 1) { - const XY<X, Y> &r = _vXY[0]; - gauss.setMean(r.x); - gauss.setMax(r.y); - gauss.setDeviation(0); - gauss.setArea(0); - gauss.setRsquared(0); - } else - assert(false); - return gauss; - } - // Gauss CalcGaussRes() { - // Gauss gauss; - // size_t nSize = _vXY.size(); - // if( nSize > 2 ) { - // gauss.Set(::calcGaussRes<uint,float>(_vXnYadc)); - // #ifdef AUXHIT_DEBUG - // if( gauss.getDeviation() != 0 ) { - // gauss.setArea(::calcGaussArea( gauss.getMean(), gauss.getMax(), gauss.getDeviation(), 0, nSize )); - // gauss.setRsquared( calcRsquared<uint,float>( _vXnYadc, gauss.getMean(), gauss.getMax(), - // gauss.getDeviation(), double(_fAdcSum) ) ); - // } else { // it cannot fit the Gauss - // gauss.setArea(0); - // gauss.setRsquared(0); - // } - // #endif - // #ifdef HIT_ADJUSTMENT - // double dPadMean = gauss.getMean(); - // if( dPadMean < _nMin || dPadMean > _nMax || gauss.getDeviation() == 0 ) - // gauss.setMean( gauss.getAdjustment() ); - // #endif - // } else if( nSize == 2 ) { - // XY<uint,float> xy0 = _vXnYadc[0]; - // XY<uint,float> xy1 = _vXnYadc[1]; - // float fMean = (xy0.x * xy0.y + xy1.x * xy1.y)/(xy0.y + xy1.y); - // gauss.setMean( fMean ); // shift to real value - // gauss.setMax( xy0.y > xy1.y ? xy0.y : xy1.y ); - // gauss.setDeviation(0); - // gauss.setArea(0); - // gauss.setRsquared(0); - // } else if( nSize == 1 ) { - // gauss.setMean( 0 ); - // gauss.setMax( _vXnYadc[0].y ); - // gauss.setDeviation(0); - // gauss.setArea(0); - // gauss.setRsquared(0); - // } else - // assert( false ); - // return gauss; - // } - std::ostringstream &Draw(std::ostringstream &oss, const char *szName, const char *szTag) const - { - uint nSize = _vXY.size(); - if (!nSize) return oss; - uint nADC_DISCRECITY2 = nADC_DISCRECITY * 2; - int nTable = nADC_DISCRECITY2 + 3 + 8; - oss << std::string(nTable, '=') << std::endl; - oss << "| " << szName << " | " << szTag << std::string(nADC_DISCRECITY2 - 6, ' ') << '|' << std::endl; - oss << std::string(nTable, '=') << std::endl; - for (uint i = 0; i < nSize; i++) { - if (i) oss << std::endl; - oss << '|' << std::setw(4) << _vXY[i].x << ": |"; - Y y = _vXY[i].y; - X x = 0; - if (_yMax > _yMin) { - float f = (_yMax - _yMin) / nADC_DISCRECITY2; - x = lroundf((y - _yMin) / f); - } else - x = nADC_DISCRECITY2 / 2; - assert(x <= nADC_DISCRECITY2); - - if (x) oss << std::string(x, '-'); - oss << (y == _yMax ? '>' : '+'); - uint nSpace = nADC_DISCRECITY2 - x + 1; - if (nSpace) oss << std::string(nSpace, ' '); - oss << '|'; - } - oss << std::endl << std::string(nTable, '=') << std::endl; - return oss; - } - std::ostringstream &PrintAdc(std::ostringstream &oss, const char *szName, const char *szTag) const - { - size_t nSize = _vXY.size(); - if (nSize) { - uint nTAB = 4; - oss << std::string(nTAB, ' ') << szTag << " " << szName << " hits:{ "; - for (uint i = 0; i < nSize; i++) oss << _vXY[i].x << ':' << _vXY[i].y << ' '; - oss << '}' << std::endl; - } - return oss; - } - std::ostringstream &PrintReg(std::ostringstream &oss, const char *szName, const char *szTag) const - { - if (_vXY.size()) { - uint nTAB = 4; - oss << std::string(nTAB, ' ') << szTag << " " << szName << " region:{ " << szName << ":{ min:" << _xMin - << " max:" << _xMax << " } adc:{ min:" << _yMin << " max:" << _yMax << " sum:" << _ySum - << " } aspect ratio:{ front:" << _fFront << " back:" << _fBack << " } }" << std::endl; - } - return oss; - } - std::ostringstream &Print(std::ostringstream &oss, const char *szName, const char *szTag) const - { - PrintReg(oss, szName, szTag); - PrintAdc(oss, szName, szTag); - return oss; - } -}; - -static const char *szDECADE = "0123456789"; -const uint nTIMEBIN_SCALE_SIZE = 30; -inline void DrawTimeBinTopScale(std::ostringstream &oss, uint nSize) -{ - oss << std::setw(4) << ":"; - for (uint i = 0; i < nSize; i++) oss << std::left << std::setw(10) << i * 10; - oss << std::endl; - oss << std::internal << std::setw(4) << ":"; - for (uint i = 0; i < nSize; i++) oss << szDECADE; - oss << std::internal << std::endl; -} -inline void DrawTimeBinBottomScale(std::ostringstream &oss, uint nSize) -{ - oss << std::setw(4) << ":"; - for (uint i = 0; i < nSize; i++) oss << szDECADE; - oss << std::endl << std::internal << std::setw(4) << ":"; - for (uint i = 0; i < nSize; i++) oss << std::left << std::setw(10) << i * 10; - oss << std::internal << std::endl; -} -struct RowStat { - uint _nPadMin; // input row statistic - uint _nPadMax; - uint _nTimeBinMin; - uint _nTimeBinMax; - float _fAdcMin; - float _fAdcMax; - void Init() - { - _nPadMin = _nTimeBinMin = UINT_MAX; - _nPadMax = _nTimeBinMax = 0; - _fAdcMin = FLT_MAX; - _fAdcMax = 0; - } - RowStat() { Init(); } - RowStat(uint nPad, uint nTimeBin, float fAdc) - { - _nPadMin = _nPadMax = nPad; - _nTimeBinMin = _nTimeBinMax = nTimeBin; - _fAdcMin = _fAdcMax = fAdc; - } - void PickOut(uint nPad, uint nTimeBin, float fAdc) - { - if (nPad > _nPadMax) _nPadMax = nPad; - if (nPad < _nPadMin) _nPadMin = nPad; - if (nTimeBin > _nTimeBinMax) _nTimeBinMax = nTimeBin; - if (nTimeBin < _nTimeBinMin) _nTimeBinMin = nTimeBin; - if (fAdc > _fAdcMax) _fAdcMax = fAdc; - if (fAdc < _fAdcMin) _fAdcMin = fAdc; - } - void Print(std::ostringstream &oss) const - { - oss << "RowStat: pad:{ min:" << _nPadMin << " max:" << _nPadMax << " } timebin:{ min:" << _nTimeBinMin - << " max:" << _nTimeBinMax << " } adc:{ min:" << _fAdcMin << " max:" << _fAdcMax << " }"; - } -}; -} // namespace tpcClustering -#endif \ No newline at end of file diff --git a/detectors/tpc/clusterHitFinder/fast/TpcEventClusters.h b/detectors/tpc/clusterHitFinder/fast/TpcEventClusters.h deleted file mode 100644 index ef03c7cd..00000000 --- a/detectors/tpc/clusterHitFinder/fast/TpcEventClusters.h +++ /dev/null @@ -1,102 +0,0 @@ -// -// Joint Institute for Nuclear Research (JINR), Dubna, Russia, Jun-2024 -// MPD TPC Clustering objects for NICA project -// author: Viktor A.Krylov JINR/LNP -// email: kryman@jinr.ru -// -#ifndef EVENT_CLUSTERS_H -#define EVENT_CLUSTERS_H - -#include "TpcSectorClusters.h" - -namespace tpcClustering { -class EventClusters { - int _nEvent; - -public: - std::list<SectorClusters *> _lstSectorClusters; - EventClusters() : _nEvent(0) {} - EventClusters(int nEvent) : _nEvent(nEvent) {} - virtual ~EventClusters() { Clear(); } - inline void Clear() - { - while (!_lstSectorClusters.empty()) { - delete _lstSectorClusters.front(); - _lstSectorClusters.pop_front(); - } - } - void DefineClusters(std::ostringstream &oss) - { -#ifdef PARALLEL - std::size_t nSize = _lstSectorClusters.size(); - pthread_t *pThreads = new pthread_t[nSize]; - uint i = 0; - for (std::list<SectorClusters *>::const_iterator iSectorCLusters = _lstSectorClusters.begin(); - iSectorCLusters != _lstSectorClusters.end(); iSectorCLusters++, i++) - int res = pthread_create(&pThreads[i], NULL, ThreadFunc, (void *)*iSectorCLusters); - for (i = 0; i < nSize; i++) pthread_join(pThreads[i], NULL); -#else - for (std::list<SectorClusters *>::const_iterator iSectorCLusters = _lstSectorClusters.begin(); - iSectorCLusters != _lstSectorClusters.end(); iSectorCLusters++) - (*iSectorCLusters)->DefineClusters(oss); -#endif - } - inline int getEvent() const { return _nEvent; } - inline void setEvent(int nEvent) { _nEvent = nEvent; } - inline bool isEmpty() const { return _lstSectorClusters.empty(); } - inline std::size_t getSize() const { return _lstSectorClusters.size(); } - inline SectorClusters *getFront() { return _lstSectorClusters.empty() ? NULL : _lstSectorClusters.front(); } - inline SectorClusters *getBack() { return _lstSectorClusters.empty() ? NULL : _lstSectorClusters.back(); } - inline void Add(SectorClusters *pSectorClusters) { _lstSectorClusters.push_back(pSectorClusters); } - inline void Load(uint nSector, uint nRow, uint nPad, uint nTimeBin, float fAdc, int nTrackId) - { - if (_lstSectorClusters.empty()) - _lstSectorClusters.push_back(new SectorClusters(nSector, nRow, nPad, nTimeBin, fAdc, nTrackId)); - else { - SectorClusters *pSectorClusters = _lstSectorClusters.back(); - if (pSectorClusters->getSector() != nSector) - _lstSectorClusters.push_back(new SectorClusters(nSector, nRow, nPad, nTimeBin, fAdc, nTrackId)); - else - pSectorClusters->Load(nRow, nPad, nTimeBin, fAdc, nTrackId); - } - } - void Type(std::ostringstream &oss) const - { - oss << "Event: " << getEvent() << std::endl; - for (std::list<SectorClusters *>::const_iterator iSectorClusters = _lstSectorClusters.begin(); - iSectorClusters != _lstSectorClusters.end(); iSectorClusters++) { - oss << " Sector: " << (*iSectorClusters)->getSector() << std::endl; - (*iSectorClusters)->Type(oss); - } - } - void Print(std::ostringstream &oss) const - { - oss << "Print:EventClusters:: event: " << getEvent() << std::endl; - for (std::list<SectorClusters *>::const_iterator iSectorClusters = _lstSectorClusters.begin(); - iSectorClusters != _lstSectorClusters.end(); iSectorClusters++) - (*iSectorClusters)->Print(oss); - } - void PrintPadCluster(std::ostringstream &oss) const - { - oss << "PrintPadClusters:EventClusters:: event: " << getEvent() << std::endl; - for (std::list<SectorClusters *>::const_iterator iSectorClusters = _lstSectorClusters.begin(); - iSectorClusters != _lstSectorClusters.end(); iSectorClusters++) - (*iSectorClusters)->PrintPadClusters(oss); - } - void Draw(std::ostringstream &oss) const - { - oss << "Draw:EventClusters:: event: " << getEvent() << std::endl; - for (std::list<SectorClusters *>::const_iterator iSectorClusters = _lstSectorClusters.begin(); - iSectorClusters != _lstSectorClusters.end(); iSectorClusters++) - (*iSectorClusters)->Draw(oss); - } - void DrawPadClusters(std::ostringstream &oss) const - { - oss << "DrawPadClusters:EventClusters:: event: " << getEvent() << std::endl; - for (std::list<SectorClusters *>::const_iterator iSectorClusters = _lstSectorClusters.begin(); - iSectorClusters != _lstSectorClusters.end(); iSectorClusters++) - (*iSectorClusters)->DrawPadClusters(oss); - } -}; -} // namespace tpcClustering -#endif \ No newline at end of file diff --git a/detectors/tpc/clusterHitFinder/fast/TpcEventObj.cxx b/detectors/tpc/clusterHitFinder/fast/TpcEventObj.cxx new file mode 100644 index 00000000..035482c9 --- /dev/null +++ b/detectors/tpc/clusterHitFinder/fast/TpcEventObj.cxx @@ -0,0 +1,339 @@ +// +// Joint Institute for Nuclear Research (JINR), Dubna, Russia, Oct-2024 +// TPC Event structures for MPD/NICA project +// author: Viktor Krylov JINR/LHEP +// email: kryman@jinr.ru +// + +#include <chrono> +#include <thread> +#include <vector> + +#include <TDirectory.h> + +#include "utility.h" +#include "TpcEventObj.h" +#include "TpcSectorObj.h" + +#include "clustering/TpcClusters.h" + +using namespace TpcClustering; + +bool bVERBOSE = false; // print additional info + +const char *szTPC_EVENTS = "TpcEvents"; // */TpcEvents +const char *szTPC_RECO_FOLDER = "reco"; // */TpcEvents/reco +const char *szTPC_MC_FOLDER = "mc"; // */TpcEvents/mc +const char *szTPC_MC_VERTEX = "vertex"; // */TpcEvents/mc/vertex +const char *szTPC_MC_TRACKS = "tracks"; // */TpcEvents/mc/tracks +const char *szTPC_MC_TRACK = "track"; // */TpcEvents/mc/tracks/id/track +const char *szTPC_MC_POINTS = "points"; // */TpcEvents/mc/tracks/id/points +const char *szTPC_MC_SAMPA = "sampa"; // */TpcEvents/mc/sectors/0-23/sampa +const char *szTPC_SECTORS = "sectors"; // */TpcEvents/*/sectors +const char *szTPC_DIGITS = "digits"; // */TpcEvents/*/sectors/0-23/digits + +TpcEventObj::TpcEventObj(TDirectory *pDir) +{ + _nEvent = atoi(pDir->GetName()); + Load(pDir); +} +TpcEventObj::~TpcEventObj() +{ + for (TpcSectorObj *pTpcSector : _vpTpcSectors) delete pTpcSector; + _vpTpcSectors.clear(); + for (TpcTrack *pTpcTrack : _vpTpcTracks) delete pTpcTrack; + _vpTpcTracks.clear(); +} +const TpcSectorObj *TpcEventObj::Find(uint nSector) const +{ + std::vector<TpcSectorObj *>::const_iterator iTpcSector = + std::find_if(_vpTpcSectors.begin(), _vpTpcSectors.end(), + [nSector](const TpcSectorObj *pTpcSector) { return pTpcSector->_nSector == nSector; }); + return iTpcSector != _vpTpcSectors.end() ? *iTpcSector : NULL; +} + +void TpcEventObj::loadTrack(TDirectory *pDir) +{ + try { + uint nTrackID = static_cast<uint>(std::stoi(pDir->GetName())); + TKey *pKey = pDir->FindKey(::szTPC_MC_TRACK); + if (pKey) { + std::vector<uint8_t> *pv; + pDir->GetObject(pKey->GetName(), pv); + if (pv) { + TKey *pPointsKey = pDir->FindKey(::szTPC_MC_POINTS); + if (pPointsKey) { + std::vector<uint8_t> *pvTpcPoints; + pDir->GetObject(pPointsKey->GetName(), pvTpcPoints); + if (pvTpcPoints) { + TpcTrack *pTpcTrack = new TpcTrack(nTrackID, TpcTrackInfo(*(const TpcTrackInfo *)pv->data())); + pTpcTrack->_pvTpcPoints = (std::vector<TpcPoint> *)pvTpcPoints; + _vpTpcTracks.push_back(pTpcTrack); + } + } + pv->clear(); + delete pv; + } else + std::cout << "Error:: tpc digits: " << ::szTPC_MC_TRACK << " has wrong object: " << pKey->GetName() + << std::endl; + } + } catch (std::invalid_argument const &ex) { + std::cout << "Error:: invalid argument: " << ex.what() << '\n'; + } +} +void TpcEventObj::loadTracks(TDirectory *pDir) +{ + assert(pDir); + + TIter iKey(pDir->GetListOfKeys()); + while (TKey *pKey = (TKey *)iKey()) { + if (pKey->IsFolder()) { + TDirectory *pTrackDir = pDir->GetDirectory(pKey->GetName()); + if (pTrackDir) loadTrack(pTrackDir); + } else + std::cout << "Error:: the folder: " << szTPC_EVENTS << " has unexpected key type: " << pKey->GetName() + << std::endl; + } +} +void TpcEventObj::loadSectors(TDirectory *pDir, uint nSector) +{ + assert(pDir); + TIter iKey(pDir->GetListOfKeys()); + while (TKey *pKey = (TKey *)iKey()) { + if (pKey->IsFolder()) { + TDirectory *pSectorDir = pDir->GetDirectory(pKey->GetName()); + if (pSectorDir) { + try { + uint nSectorDir = static_cast<uint>(std::stoi(pSectorDir->GetName())); + if (nSector == UINT32_MAX || nSectorDir == nSector) { + TpcSectorObj *pTpcSector = new TpcSectorObj(nSectorDir); + if (pTpcSector->loadSampaBuf(pSectorDir)) { + pTpcSector->loadSampaBuf(); + pTpcSector->loadDigits(pSectorDir); // it is needed for statistical analysis + } else if (pTpcSector->loadDigits(pSectorDir)) { + pTpcSector->loadDigits(); + // pTpcSector->createSampaBuf(); // create SampaBuf to store it in a new root file + } + _vpTpcSectors.push_back(pTpcSector); + } + } catch (std::invalid_argument const &ex) { + std::cout << "Error:: loadSectors: invalid argument: " << ex.what() << std::endl; + } + } else + std::cout << "Error:: loadSectors: getDirectory!\n"; + } else + std::cout << "Error:: loadSectors: wrong key type!\n"; + } +} + +#define PTHREAD_FUNC +// #undef PTHREAD_FUNC +#ifdef PTHREAD_FUNC +#include <pthread.h> +void *ThreadFuncDefClusters(void *pArgs) +{ + TpcSectorObj *pTpcSector = (TpcSectorObj *)pArgs; + pTpcSector->_pTpcClusters->defClusters(); + pTpcSector->_pTpcClusters->defDigits(pTpcSector->_pvTpcDigits); + pTpcSector->defAdcHits(); + pthread_exit(NULL); +} +#endif + +void TpcEventObj::Calc() +{ +#ifdef PTHREAD_FUNC + std::size_t nSize = _vpTpcSectors.size(); + pthread_t *pThreads = new pthread_t[nSize]; + // std::cout << "start working threads:" << nSize << std::endl; + for (uint i = 0; i < nSize; i++) { + int res = pthread_create(&pThreads[i], NULL, ThreadFuncDefClusters, (void *)_vpTpcSectors[i]); + } + for (uint i = 0; i < nSize; i++) pthread_join(pThreads[i], NULL); + delete[] pThreads; +// std::cout << "finish working threads:" << nSize << std::endl; +#else + for (TpcSectorObj *pTpcSector : _vpTpcSectors) { + pTpcSector->_pTpcClusters->defClusters(); + pTpcSector->_pTpcClusters->defDigits(pTpcSector->_pvTpcDigits); + pTpcSector->defAdcHits(); + } +#endif +} +// void TpcEvent::Calc() { // for debug only! +// std::ostringstream oss; +// for( TpcSector* pTpcSector: _vpTpcSector ) { +// pTpcSector->_pTpcClusters->defClusters(oss); +// pTpcSector->_pTpcClusters->defDigits(pTpcSector->_pvMcDigits); +// } +// std::cout << oss.str(); +// oss.clear(); +// } +void TpcEventObj::Load(TDirectory *pDir, uint nSector /* = UINT32_MAX*/) +{ + TKey *pMcKey = pDir->FindKey(szTPC_MC_FOLDER); + if (pMcKey && pMcKey->IsFolder()) { + TDirectory *pMcDir = pDir->GetDirectory(pMcKey->GetName()); + TIter iKey(pMcDir->GetListOfKeys()); + while (TKey *pKey = (TKey *)iKey()) { + if (pKey->IsFolder()) { + if (!std::strcmp(pKey->GetName(), szTPC_MC_TRACKS)) { + TDirectory *pTrackDir = pMcDir->GetDirectory(pKey->GetName()); + if (pTrackDir) loadTracks(pTrackDir); + } else if (!std::strcmp(pKey->GetName(), szTPC_SECTORS)) { + TDirectory *pSectorDir = pMcDir->GetDirectory(pKey->GetName()); + if (pSectorDir) loadSectors(pSectorDir, nSector); + } else + std::cout << "Error:: loadEvent: has wrong folder name: " << pKey->GetName() << std::endl; + } else if (!std::strcmp(pKey->GetName(), szTPC_MC_VERTEX)) { + std::vector<uint8_t> *pv; + pMcDir->GetObject(pKey->GetName(), pv); + if (pv) { + _TpcVertex = *(TpcVertex *)pv->data(); + pv->clear(); + delete pv; + } else + std::cout << "Error:: McVertex: " << szTPC_MC_VERTEX << " has wrong object data!" << std::endl; + } else + std::cout << "Error:: the key: '" << pKey->GetName() << "' has wrong name!" << std::endl; + } + } else + std::cout << "Error:: the folder: '" << szTPC_MC_FOLDER << "' cannot be found!" << std::endl; +} +void TpcEventObj::McStore(TDirectory *pDir) const +{ + TDirectory *pMcDir = pDir->mkdir(szTPC_MC_FOLDER, "Tpc Monte Carlo", kTRUE); + if (pMcDir && pMcDir->cd()) { + const TpcVertex *pTpcVertex = &_TpcVertex; + std::vector<uint8_t> vVertex((uint8_t *)pTpcVertex, (uint8_t *)pTpcVertex + sizeof(TpcVertex)); + pMcDir->WriteObject(&vVertex, szTPC_MC_VERTEX); + + TDirectory *pTracksDir = pMcDir->mkdir(szTPC_MC_TRACKS, "Tpc Event McTracks", kTRUE); + if (pTracksDir && pTracksDir->cd()) { + for (const TpcTrack *pTpcTrack : _vpTpcTracks) pTpcTrack->Store(pTracksDir); + } else + std::cout << "Error:: TpcEvent::MC: the folder: '" << szTPC_MC_TRACKS << "' cannot be found or created!" + << std::endl; + + TDirectory *pSectorsDir = pMcDir->mkdir(szTPC_SECTORS, "Tpc Sectors of Digits", kTRUE); + if (pSectorsDir && pSectorsDir->cd()) { + for (const TpcSectorObj *pTpcSector : _vpTpcSectors) pTpcSector->McStore(pSectorsDir); + } else + std::cout << "Error:: TpcEvent::MC: the folder: '" << szTPC_SECTORS << "' cannot be found or created!" + << std::endl; + } +} +void TpcEventObj::Store(TDirectory *pDir) const +{ + TDirectory *pRecoDir = pDir->mkdir(szTPC_RECO_FOLDER, "Tpc Reconstruction", kTRUE); + if (pRecoDir && pRecoDir->cd()) { + TDirectory *pSectorsDir = pRecoDir->mkdir(szTPC_SECTORS, "Tpc Sector Clusters", kTRUE); + if (pSectorsDir && pSectorsDir->cd()) { + for (const TpcSectorObj *pTpcSector : _vpTpcSectors) pTpcSector->Store(pSectorsDir, _vpTpcTracks); + } else + std::cout << "Error:: TpcEvent::Store: the folder: '" << szTPC_SECTORS << "' cannot be found or created!" + << std::endl; + } else + std::cout << "Error:: TpcEvent::Store: the folder: '" << szTPC_RECO_FOLDER << "' cannot be found or created!" + << std::endl; +} +void readTpcEvent(TDirectory *pDir, TDirectory *pStoreDir, uint nEvent) +{ + TKey *pKey = pDir->FindKey(std::to_string(nEvent).c_str()); + if (pKey && pKey->IsFolder()) { + TDirectory *pEventDir = pDir->GetDirectory(pKey->GetName()); + if (pEventDir) { + std::cout << "TPC Event:" << nEvent << std::endl; + TpcEventObj tpcEvent(nEvent); + tpcEvent.Load(pEventDir); // data loader + tpcEvent.Calc(); // define clusters and their hits + if (pStoreDir) { + if (pStoreDir != pDir) { + TDirectory *pEventDir_ = pStoreDir->mkdir(pEventDir->GetName(), "TPC Event", kTRUE); + if (pEventDir_ && pEventDir_->cd()) { + tpcEvent.McStore(pEventDir_); // store mc tracks to the new file directory + tpcEvent.Store(pEventDir_); // store results to the new file directory + } + } else + tpcEvent.Store(pEventDir); // store results to the file directory + } + } + } else + std::cout << "Error:: the folder: " << szTPC_EVENTS << ": hasn't event number: " << nEvent + << " or has unexpected key type: " << pKey->GetName() << std::endl; +} +void readTpcEvents(TDirectory *pDir, TDirectory *pStoreDir) +{ + if (pDir->GetNkeys()) { + TIter iKey(pDir->GetListOfKeys()); + while (TKey *pKey = (TKey *)iKey()) { + if (pKey->IsFolder()) { + TDirectory *pEventDir = pDir->GetDirectory(pKey->GetName()); + if (pEventDir) { + try { + uint nEvent = static_cast<uint>(std::stoi(pEventDir->GetName())); + std::cout << "TPC Event:" << nEvent << std::endl; + TpcEventObj tpcEvent(nEvent); + CpuTimer ct; + ct.Start(); + tpcEvent.Load(pEventDir); // data loader + ct.Stop(); + ct.PrintDuration("Event load time,ms:"); + ct.Start(); + tpcEvent.Calc(); // define clusters and their hits + ct.Stop(); + ct.PrintDuration("Event calc time,ms:"); + if (pStoreDir) { + if (pStoreDir != pDir) { + TDirectory *pEventDir_ = pStoreDir->mkdir(pEventDir->GetName(), "TPC Event", kTRUE); + if (pEventDir_ && pEventDir_->cd()) { + tpcEvent.McStore(pEventDir_); // store mc tracks to the new file directory + tpcEvent.Store(pEventDir_); // store results to the new file directory + } + } else + tpcEvent.Store(pEventDir); // store adc hits to the file directory + } + } catch (std::invalid_argument const &ex) { + std::cout << "Error:: invalid argument: " << ex.what() << '\n'; + } + } + } else + std::cout << "Error:: the folder: " << szTPC_EVENTS << " has unexpected key type: " << pKey->GetName() + << std::endl; + } + } else + std::cout << "Error:: the folder: " << szTPC_EVENTS << ": has no one keys!" << std::endl; +} +void readTpcRootFile(TFile *pFile, const char *szAction, const char *szPostfix, uint nEvent) +{ + TKey *pKey = pFile->FindKey(szTPC_EVENTS); + if (pKey) { + TDirectory *pDir = pFile->GetDirectory(pKey->GetName()); + if (pDir && pDir->cd()) { + if (!std::strcmp(szAction, "NEW")) { + std::string strNewFile = pFile->GetName(); + size_t nPos = strNewFile.rfind('.'); + if (nPos != std::string::npos) strNewFile.insert(nPos, szPostfix); + TFile *pNewFile = new TFile(strNewFile.c_str(), "NEW"); + if (pNewFile && pNewFile->IsOpen()) { + TDirectory *pStoreDir = pNewFile->mkdir(szTPC_EVENTS, "TPC Events", kTRUE); + if (pStoreDir && pStoreDir->cd()) + (nEvent == UINT32_MAX) ? readTpcEvents(pDir, pStoreDir) : readTpcEvent(pDir, pStoreDir, nEvent); + else + std::cout << "Error:: TpcEvent:: the folder: '" << szTPC_EVENTS << "' cannot be created!" + << std::endl; + pNewFile->Close(); + delete pNewFile; + } else + std::cout << "Error:: the root file: '" << strNewFile << "' cannot be created!" << std::endl; + } else { + TDirectory *pStoreDir = !std::strcmp(szAction, "UPDATE") ? pDir : NULL; + (nEvent == UINT32_MAX) ? readTpcEvents(pDir, pStoreDir) : readTpcEvent(pDir, pStoreDir, nEvent); + } + } else + std::cout << "Error:: the TPC root file: '" << pFile->GetName() << "' do not have the folder: " << szTPC_EVENTS + << std::endl; + } else + std::cout << "Error:: the TPC root file: '" << pFile->GetName() << "' do not have the key: " << szTPC_EVENTS + << std::endl; +} diff --git a/detectors/tpc/clusterHitFinder/fast/TpcEventObj.h b/detectors/tpc/clusterHitFinder/fast/TpcEventObj.h new file mode 100644 index 00000000..a562f188 --- /dev/null +++ b/detectors/tpc/clusterHitFinder/fast/TpcEventObj.h @@ -0,0 +1,89 @@ +// +// Joint Institute for Nuclear Research (JINR), Dubna, Russia, Oct-2024 +// TPC Event structure for MPD/NICA project +// author: Viktor Krylov JINR/LHEP +// email: kryman@jinr.ru +// +#ifndef TPC_EVENT_H +#define TPC_EVENT_H + +#include <iomanip> +#include <iostream> +#include <fstream> +#include <stdlib.h> +#include <filesystem> +#include <cassert> +#include <cstring> +#include <set> +#include <vector> +#include <numeric> + +#include <TFile.h> +#include <TDirectory.h> +#include <TKey.h> +#include <TH2S.h> + +#include "clustering/TpcTrack.h" + +#define _RED_ "0;31" +#define _GREEN_ "0;32" +#define _YELLOW_ "0;33" +#define _BASH_CODE_(_CODE) "\e[" _CODE "m" +#define _RED_MSG_(_MSG) _BASH_CODE_(_RED_) _MSG _BASH_CODE_() +#define _MSG_1_(_0, _1) std::cout << _BASH_CODE_(_0) << _1 << _BASH_CODE_() +#define _MSG_2_(_0, _1, _2) std::cout << _BASH_CODE_(_0) << _1 << _BASH_CODE_() << _2 << std::endl +#define _MSG_3_(_0, _1, _2, _3) std::cout << _BASH_CODE_(_0) << _1 << _2 << _3 << _BASH_CODE_() + +class TFile; +class TDirectory; +class TpcSectorObj; + +struct TpcEventObj { + uint _nEvent; + TpcClustering::TpcVertex _TpcVertex; + std::vector<TpcClustering::TpcTrack *> _vpTpcTracks; + std::vector<TpcSectorObj *> _vpTpcSectors; + TpcEventObj() : _nEvent(0) {} + TpcEventObj(uint nEvent) : _nEvent(nEvent) {} + TpcEventObj(TDirectory *pDir); + // { + // _nEvent = atoi(pDir->GetName()); + // Load(pDir); + // } + // TpcEventObj( TDirectory* pDir ) { + // try { + // _nEvent = static_cast<uint>( std::stoi(pDir->GetName()) ); + // } catch (std::invalid_argument const& ex) { + // std::cout << "Error:: invalid argument: '" << pDir->GetName() << "' " << ex.what() << '\n'; + // } + // } + virtual ~TpcEventObj(); + // { + // for( TpcSector* pTpcSector : _vpTpcSectors ) + // delete pTpcSector; + // _vpTpcSectors.clear(); + // for( TpcTrack* pTpcTrack : _vpTpcTracks ) + // delete pTpcTrack; + // _vpTpcTracks.clear(); + // } + void loadTrack(TDirectory *pDir); + void loadTracks(TDirectory *pDir); + void loadSectors(TDirectory *pDir, uint nSector); + // void loadEvent( TDirectory* pDir, uint nSector = UINT32_MAX ); + void Calc(); + void McStore(TDirectory *pDir) const; + void Store(TDirectory *pDir) const; + void Load(TDirectory *pDir, uint nSector = UINT32_MAX); + // void Load( TDirectory* pDir ); + bool IsEmpty() const { return _vpTpcSectors.size() == 0; } + const TpcSectorObj *Find(uint nSector) const; + // { + // std::vector<TpcSector*>::const_iterator iTpcSector = std::find_if(_vpTpcSectors.begin(), _vpTpcSectors.end(), + // [nSector](const TpcSector* pTpcSector){ return pTpcSector->_nSector == nSector; } + // ); + // return iTpcSector != _vpTpcSectors.end() ? *iTpcSector : NULL; + // } +}; +void readTpcRootFile(TFile *pFile, const char *szAction, const char *szPostfix = NULL, uint nEvent = UINT32_MAX); + +#endif // TPC_EVENT_H diff --git a/detectors/tpc/clusterHitFinder/fast/TpcPadCluster.h b/detectors/tpc/clusterHitFinder/fast/TpcPadCluster.h deleted file mode 100644 index 1a2be1f0..00000000 --- a/detectors/tpc/clusterHitFinder/fast/TpcPadCluster.h +++ /dev/null @@ -1,370 +0,0 @@ -// -// Joint Institute for Nuclear Research (JINR), Dubna, Russia, Jun-2024 -// MPD TPC Clustering objects for NICA project -// author: Viktor A.Krylov JINR/LNP -// email: kryman@jinr.ru -// -#ifndef PAD_CLUSTER_H -#define PAD_CLUSTER_H - -#include <vector> -#include <cstring> -#include <algorithm> - -#include "TpcAdcHit.h" -#include "TpcPadHit.h" - -// BASH COLOR CODES: {black:30,red:31,green:32,yellow:33,blue:34,magenta:35,cyan:36,white:37}; background:+10; -// bright:+60 BASH TEXT DECORATION: {regular:0,bold:1,italic:3,underscore:4} BASH CODE STRUCTURE: <PREFIX>;<COLOR>;<TEXT -// DECORATION> -#define _RED_ "0;31" -#define _GREEN_ "0;32" -#define _YELLOW_ "0;33" -#define _BASH_CODE_(_CODE) "\e[" _CODE "m" -#define _RED_MSG_(_MSG) _BASH_CODE_(_RED_) _MSG _BASH_CODE_() - -const float fADC_DOORSTEP = 2.5; // adc sensitivity -static const char *SZ_ALPHABET = "ABCDEFJHIJKLMNOPQRSTUVWXYZ"; - -namespace tpcClustering { -class Cluster; -class PadCluster : public PadHit { - Cluster *_pCluster; // link to a Cluster where the PadCluster belongs - const PadCluster *_pPadClusterBefore; // link to prior PadCluster (!= null: combined cluster) - const PadCluster *_pPadClusterAfter; // link to next PadCluster (!= null: combined cluster) - - uint _nTimeBinMin; // timebin min range - uint _nTimeBinMax; // timebin max range - - float _fAdcMax; // max adc value in the AdcHit vector - float _fAdcMin; // min adc value in the AdcHit vector - - uint _nTimeBinLapMin; // overlapping timebin min range in neighbour pad - uint _nTimeBinLapMax; // overlapping timebin max range in neighbour pad - - uint _iAdcMax; // index of the max adc value in the AdcHit vector - uint _iAdcMin; // index of the last min adc value in the AdcHit vector - -public: - std::vector<AdcHit> _vAdcHits; // vector of ADC values - - PadCluster() - : PadHit(), _pCluster(NULL), _pPadClusterBefore(NULL), _pPadClusterAfter(NULL), _nTimeBinMin(UINT_MAX), - _nTimeBinMax(UINT_MAX), _fAdcMin(0), _fAdcMax(0), _nTimeBinLapMin(0), _nTimeBinLapMax(0), _iAdcMin(0), - _iAdcMax(0) - { - } - PadCluster(int nPad) - : PadHit(nPad), _pCluster(NULL), _pPadClusterBefore(NULL), _pPadClusterAfter(NULL), _nTimeBinMin(UINT_MAX), - _nTimeBinMax(UINT_MAX), _fAdcMin(0), _fAdcMax(0), _nTimeBinLapMin(0), _nTimeBinLapMax(0), _iAdcMin(0), - _iAdcMax(0) - { - } - PadCluster(int nPad, const AdcHit &r) - : PadHit(nPad, r), _pCluster(NULL), _pPadClusterBefore(NULL), _pPadClusterAfter(NULL), _nTimeBinLapMin(0), - _nTimeBinLapMax(0), _iAdcMin(0), _iAdcMax(0) - { - _nTimeBinMin = _nTimeBinMax = r.getTimeBin(); - _fAdcMin = _fAdcMax = r.getAdc(); - _vAdcHits.push_back(r); - } - PadCluster(int nPad, const AdcHit &r, Cluster *pCluster) - : PadHit(nPad, r), _pCluster(pCluster), _pPadClusterBefore(NULL), _pPadClusterAfter(NULL), _nTimeBinLapMin(0), - _nTimeBinLapMax(0), _iAdcMin(0), _iAdcMax(0) - { - _nTimeBinMin = _nTimeBinMax = r.getTimeBin(); - _fAdcMin = _fAdcMax = r.getAdc(); - _vAdcHits.push_back(r); - } - PadCluster(const PadCluster &r) - : PadHit(r), _pCluster(r._pCluster), _pPadClusterBefore(r._pPadClusterBefore), - _pPadClusterAfter(r._pPadClusterAfter), _nTimeBinMin(r._nTimeBinMin), _nTimeBinMax(r._nTimeBinMax), - _fAdcMin(r._fAdcMin), _fAdcMax(r._fAdcMax), _nTimeBinLapMin(r._nTimeBinLapMin), - _nTimeBinLapMax(r._nTimeBinLapMax), _iAdcMin(r._iAdcMin), _iAdcMax(r._iAdcMax) - { - _vAdcHits = r._vAdcHits; - } - virtual ~PadCluster() { _vAdcHits.clear(); } - AdcHit operator[](size_t n) const - { - assert(n < _vAdcHits.size()); - return _vAdcHits[n]; - } - const AdcHit &operator[](size_t n) - { - assert(n < _vAdcHits.size()); - return _vAdcHits[n]; - } - inline uint getClusterId(const Cluster *pCluster) const; - inline float getClusterTime(const Cluster *pCluster) const; - inline float getClusterPad(const Cluster *pCluster) const; - - inline const std::vector<AdcHit> &getAdcHits() const { return _vAdcHits; } - - inline size_t getSize() const { return _vAdcHits.size(); } - inline Cluster *getCluster() const { return _pCluster; } - inline void setCluster(Cluster *pCluster) { _pCluster = pCluster; } - inline AdcHit getAdcHitMax() const { return _vAdcHits.empty() ? AdcHit(0, 0.0) : _vAdcHits[_iAdcMax]; } - inline AdcHit getAdcHitMin() const { return _vAdcHits.empty() ? AdcHit(0, 0.0) : _vAdcHits[_iAdcMin]; } - inline uint getTimeBinMin() const { return _nTimeBinMin; } - inline uint getTimeBinMax() const { return _nTimeBinMax; } - inline float getAdcMin() const { return _fAdcMin; } - inline float getAdcMax() const { return _fAdcMax; } - inline void setAdcMax(float f) - { - _fAdcMax = f; - } // correct AdcMax PadCluster value for splitted cluster by Scaling factor - inline uint getTimeBinLapMin() const { return _nTimeBinLapMin; } - inline uint getTimeBinLapMax() const { return _nTimeBinLapMax; } - - inline const PadCluster *getPadClusterBefore() const { return _pPadClusterBefore; } - inline void setPadClusterBefore(const PadCluster *pPadCluster) { _pPadClusterBefore = pPadCluster; } - inline const PadCluster *getPadClusterAfter() const { return _pPadClusterAfter; } - inline void setPadClusterAfter(const PadCluster *pPadCluster) { _pPadClusterAfter = pPadCluster; } - - inline void setPadTimeOverlapping(const PadCluster *pPadCluster) - { - assert(!(isBefore(pPadCluster) || isAfter(pPadCluster))); - - _nTimeBinLapMin = _nTimeBinMin > pPadCluster->getTimeBinMin() ? _nTimeBinMin : pPadCluster->getTimeBinMin(); - _nTimeBinLapMax = _nTimeBinMax < pPadCluster->getTimeBinMax() ? _nTimeBinMax : pPadCluster->getTimeBinMax(); - } - inline bool isOut(const AdcHit &r) const { return r.getTimeBin() > _nTimeBinMax + 1; } - inline bool isCombined() const { return _pPadClusterBefore != NULL || _pPadClusterAfter != NULL; } - inline bool isBelongTo(const PadCluster *pPadCluster) const - { - if (isCombined() && pPadCluster->isCombined()) { // 1:1 - uint nTimeBinMin = _pPadClusterBefore == NULL ? _nTimeBinMin : _nTimeBinMin + 1; - uint nTimeBinMax = _pPadClusterAfter == NULL ? _nTimeBinMax : _nTimeBinMax - 1; - uint nTimeBinMinExt = pPadCluster->getPadClusterBefore() == NULL ? pPadCluster->getTimeBinMin() - : pPadCluster->getTimeBinMin() + 1; - uint nTimeBinMaxExt = - pPadCluster->getPadClusterAfter() == NULL ? pPadCluster->getTimeBinMax() : pPadCluster->getTimeBinMax() - 1; - uint nAdcHitMaxTimeBin = getAdcHitMax().getTimeBin(); - uint nAdcHitMaxTimeBinExt = pPadCluster->getAdcHitMax().getTimeBin(); - - bool bMin = nTimeBinMin < nTimeBinMinExt && pPadCluster->getPadClusterBefore() != NULL; - bool bMax = nTimeBinMax > nTimeBinMaxExt && pPadCluster->getPadClusterAfter() != NULL; - bool bMinMax = bMin || bMax; - bool bMinExt = nTimeBinMinExt < nTimeBinMin && _pPadClusterBefore != NULL; - bool bMaxExt = nTimeBinMaxExt > nTimeBinMax && _pPadClusterAfter != NULL; - bool bMinMaxExt = bMinExt || bMaxExt; - if (bMinMax && bMinMaxExt) - return nAdcHitMaxTimeBin >= nTimeBinMinExt && nAdcHitMaxTimeBin <= nTimeBinMaxExt && - nAdcHitMaxTimeBinExt >= nTimeBinMin && nAdcHitMaxTimeBinExt <= nTimeBinMax; - else if (bMinMax) - return nAdcHitMaxTimeBin >= nTimeBinMinExt && nAdcHitMaxTimeBin <= nTimeBinMaxExt; - else if (bMinMaxExt) - return nAdcHitMaxTimeBinExt >= nTimeBinMin && nAdcHitMaxTimeBinExt <= nTimeBinMax; - else - return true; - } else if (isCombined()) { // 1:0 - uint nTimeBinMin = _pPadClusterBefore == NULL ? _nTimeBinMin : _nTimeBinMin + 1; - uint nTimeBinMax = _pPadClusterAfter == NULL ? _nTimeBinMax : _nTimeBinMax - 1; - uint nAdcHitMaxTimeBinExt = pPadCluster->getAdcHitMax().getTimeBin(); - return (_pPadClusterBefore == NULL && pPadCluster->getTimeBinMax() <= _nTimeBinMax) || - (_pPadClusterAfter == NULL && pPadCluster->getTimeBinMin() >= _nTimeBinMin) || - (nAdcHitMaxTimeBinExt >= nTimeBinMin && nAdcHitMaxTimeBinExt <= nTimeBinMax); // max hit - } else if (pPadCluster->isCombined()) { // 0:1 - uint nTimeBinMinExt = pPadCluster->getPadClusterBefore() == NULL ? pPadCluster->getTimeBinMin() - : pPadCluster->getTimeBinMin() + 1; - uint nTimeBinMaxExt = - pPadCluster->getPadClusterAfter() == NULL ? pPadCluster->getTimeBinMax() : pPadCluster->getTimeBinMax() - 1; - uint nAdcHitMaxTimeBin = getAdcHitMax().getTimeBin(); - return (pPadCluster->getPadClusterBefore() == NULL && _nTimeBinMax <= nTimeBinMaxExt) || - (pPadCluster->getPadClusterAfter() == NULL && _nTimeBinMin >= nTimeBinMinExt) || - (nAdcHitMaxTimeBin >= nTimeBinMinExt && nAdcHitMaxTimeBin <= nTimeBinMaxExt); // max hit - } else // 0:0 - return true; - } - inline bool isAfter(const PadCluster *pPadCluster) const { return _nTimeBinMin > pPadCluster->getTimeBinMax(); } - inline bool isBefore(const PadCluster *pPadCluster) const { return _nTimeBinMax < pPadCluster->getTimeBinMin(); } - AdcHit splitAdcHit(const AdcHit &rAdcHit) - { // proportionally split AdcHit value between the two combined clusters; - std::size_t nSize = _vAdcHits.size(); - assert(nSize > 1); - - AdcHit adcHit = _vAdcHits.back(); - float adcMid = adcHit.getAdc(); - float adcLeft = _vAdcHits[nSize - 2].getAdc() - adcMid; - float adcRight = rAdcHit.getAdc() - adcMid; - assert(adcLeft > 0 && adcRight > 0); - - float fRatio = 1 + adcRight / adcLeft; // coefficient of proportionality (ratio: left to right) - _vAdcHits.back().setAdc(adcMid / fRatio); - adcHit.setAdc(adcMid - _vAdcHits.back().getAdc()); - PadHit::ReduceAdcSum(adcHit.getAdc()); - return adcHit; - } - inline bool isCombined(const AdcHit &r) const - { - return _iAdcMax < _iAdcMin && r.getAdc() > getAdcHitMin().getAdc() + fADC_DOORSTEP; - } - inline void Add(const AdcHit &r) - { - int nTimeBin = r.getTimeBin(); - if (_nTimeBinMin == UINT_MAX && _nTimeBinMax == UINT_MAX) - _nTimeBinMin = _nTimeBinMax = nTimeBin; - else if (nTimeBin < _nTimeBinMin) - _nTimeBinMin = nTimeBin; - else if (nTimeBin > _nTimeBinMax) - _nTimeBinMax = nTimeBin; - - float fAdc = r.getAdc(); - if (_fAdcMin == 0 && _fAdcMax == 0) - _fAdcMin = _fAdcMax = fAdc; - else if (fAdc < _fAdcMin) - _fAdcMin = fAdc; - else if (fAdc > _fAdcMax) - _fAdcMax = fAdc; - - size_t nSize = _vAdcHits.size(); - if (_iAdcMax > _iAdcMin) { - if (fAdc > getAdcHitMax().getAdc()) - _iAdcMax = nSize; - else - _iAdcMin = nSize; - } else { - if (fAdc < getAdcHitMin().getAdc()) - _iAdcMin = nSize; - else - _iAdcMax = nSize; - } - PadHit::AddAdc(fAdc); - _vAdcHits.push_back(r); - } - inline void ReCalcStats() - { - PadHit::clearAdcHitSum(); - _nTimeBinMin = _nTimeBinMax = UINT_MAX; - _fAdcMin = _fAdcMax = 0; - _iAdcMin = _iAdcMax = 0; - int i = 0; - for (std::vector<AdcHit>::const_iterator it = _vAdcHits.begin(); it != _vAdcHits.end(); it++, i++) { - int nTimeBin = it->getTimeBin(); - if (_nTimeBinMin == UINT_MAX && _nTimeBinMax == UINT_MAX) - _nTimeBinMin = _nTimeBinMax = nTimeBin; - else if (nTimeBin < _nTimeBinMin) - _nTimeBinMin = nTimeBin; - else if (nTimeBin > _nTimeBinMax) - _nTimeBinMax = nTimeBin; - - float fAdc = it->getAdc(); - if (_fAdcMin == 0 && _fAdcMax == 0) - _fAdcMin = _fAdcMax = fAdc; - else if (fAdc < _fAdcMin) { - _iAdcMin = i; - _fAdcMin = fAdc; - } else if (fAdc > _fAdcMax) { - _iAdcMax = i; - _fAdcMax = fAdc; - } - PadHit::AddAdc(fAdc); - } - } - std::ostringstream &Draw(std::ostringstream &oss, uint nTimeBinMin, uint nTimeBinMax, const Cluster *pCluster, - char ch = '+') const - { - uint nTimeBin = pCluster ? lroundf(getClusterTime(pCluster)) : UINT_MAX; - uint nPad = pCluster ? lroundf(getClusterPad(pCluster)) : UINT_MAX; - uint nMax = getAdcHitMax().getTimeBin(); - size_t nSize = _vAdcHits.size(); - for (uint i = nTimeBinMin, j = 0; i <= nTimeBinMax; i++) - if (j < nSize && i == _vAdcHits[j].getTimeBin()) { - if (i == nTimeBin && nPad == PadHit::getPad()) - oss << _RED_MSG_("*"); // red cluster hit star '*' - else if (i == nMax && nSize > 1) - oss << '^'; - else - oss << ch; - j++; - } else - oss << ' '; - return oss; - } - std::ostringstream &Draw2(std::ostringstream &oss, uint nPadMin, uint nPadMax, char ch = '+') const - { - uint nTimeBin = _pCluster ? lroundf(getClusterTime(_pCluster)) : UINT_MAX; - uint nPad = _pCluster ? lroundf(getClusterPad(_pCluster)) : UINT_MAX; - uint nMax = getAdcHitMax().getTimeBin(); - size_t nSize = _vAdcHits.size(); - for (uint i = nPadMin, j = 0; i <= nPadMax; i++) - if (j < nSize && i == _vAdcHits[j].getTimeBin()) { - if (i == nPad && nTimeBin == PadHit::getPad()) - oss << _RED_MSG_("*"); // red cluster hit star '*' - else if (i == nMax && nSize > 1) - oss << '^'; - else - oss << ch; - j++; - } else - oss << ' '; - return oss; - } - std::ostringstream &Draw(std::ostringstream &oss, char ch = '+') const - { - uint nTimeBin = _pCluster ? lroundf(getClusterTime(_pCluster)) : UINT_MAX; - uint nPad = _pCluster ? lroundf(getClusterPad(_pCluster)) : UINT_MAX; - uint nMax = getAdcHitMax().getTimeBin(); - size_t nSize = _vAdcHits.size(); - for (uint i = 0; i < nSize; i++) { - uint n = _vAdcHits[i].getTimeBin(); - if (n == nTimeBin && nPad == PadHit::getPad()) - oss << _RED_MSG_("*"); // red cluster hit star '*' - else if (n == nMax && nSize > 1) - oss << '^'; - else - oss << ch; - } - return oss; - } - std::ostringstream &Draw(std::ostringstream &oss, const Cluster *pCluster, char ch = '+') const - { - assert(pCluster); - uint nTimeBin = pCluster ? lroundf(getClusterTime(pCluster)) : UINT_MAX; - uint nPad = _pCluster ? lroundf(getClusterPad(pCluster)) : UINT_MAX; - uint nMax = getAdcHitMax().getTimeBin(); - size_t nSize = _vAdcHits.size(); - for (uint i = 0; i < nSize; i++) { - uint n = _vAdcHits[i].getTimeBin(); - if (n == nTimeBin && nPad == PadHit::getPad()) - oss << _RED_MSG_("*"); // red cluster hit star '*' - else if (n == nMax && nSize > 1) - oss << '^'; - else - oss << ch; - } - return oss; - } - std::ostringstream &PrintPad(std::ostringstream &oss) const - { - PadHit::PrintPad(oss); - uint nClusterId = getClusterId(_pCluster); - char chClusterId = _pCluster == NULL ? '-' : SZ_ALPHABET[nClusterId % std::strlen(SZ_ALPHABET)]; - const uint nTab = 4; - oss << std::string(nTab, ' ') << "PadCluster: timebins:[" << _nTimeBinMin << "," << _nTimeBinMax - << "], adcSpikes:{ min[" << _iAdcMin << "], max[" << _iAdcMax << "] }, combined:{" - << (_pPadClusterBefore != NULL ? '+' : '-') << (_pPadClusterAfter != NULL ? '+' : '-') << "}, overlapping:[" - << _nTimeBinLapMin << "," << _nTimeBinLapMax << ']' << std::endl; - oss << std::string(nTab, ' ') << "AdcHits: { "; - for (uint i = 0; i < _vAdcHits.size(); i++) oss << _vAdcHits[i] << ' '; - oss << '}' << std::endl; - return oss; - } - std::ostringstream &PrintTime(std::ostringstream &oss) const - { - PadHit::PrintTime(oss); - uint nClusterId = getClusterId(_pCluster); - char chClusterId = _pCluster == NULL ? '-' : SZ_ALPHABET[nClusterId % strlen(SZ_ALPHABET)]; - const uint nTab = 4; - oss << std::string(nTab, ' ') << "TimeCluster:: pads: [" << _nTimeBinMin << "," << _nTimeBinMax << "]" - << ", adcSpikes: { min[" << _iAdcMin << "], max[" << _iAdcMax << "] }" << ", cluster: " << nClusterId << " '" - << chClusterId << "'" << std::endl; - oss << std::string(nTab, ' ') << "AdcHits: { "; - for (uint i = 0; i < _vAdcHits.size(); i++) oss << _vAdcHits[i] << ' '; - oss << '}' << std::endl; - return oss; - } -}; -} // namespace tpcClustering -#endif diff --git a/detectors/tpc/clusterHitFinder/fast/TpcPadHit.h b/detectors/tpc/clusterHitFinder/fast/TpcPadHit.h deleted file mode 100644 index 8762bbd3..00000000 --- a/detectors/tpc/clusterHitFinder/fast/TpcPadHit.h +++ /dev/null @@ -1,52 +0,0 @@ -// -// Joint Institute for Nuclear Research (JINR), Dubna, Russia, Jun-2024 -// MPD TPC Clustering objects for NICA project -// author: Viktor A.Krylov JINR/LNP -// email: kryman@jinr.ru -// -#ifndef PADHIT_H -#define PADHIT_H - -namespace tpcClustering { -class PadHit { - uint _nPad; // Pad id number - float _fAdcHitSum; // total ADC sum -public: - PadHit() : _nPad(UINT_MAX), _fAdcHitSum(0) {} - PadHit(uint nPad) : _nPad(nPad), _fAdcHitSum(0) {} - PadHit(uint nPad, const AdcHit &r) : _nPad(nPad) { _fAdcHitSum = r.getAdc(); } - PadHit(const PadHit &r) : _nPad(r._nPad), _fAdcHitSum(r._fAdcHitSum) {} - ~PadHit() {} - inline void clearAdcHitSum() { _fAdcHitSum = 0; } - PadHit &operator=(const PadHit &r) - { - if (this != &r) { - _nPad = r._nPad; - _fAdcHitSum = r._fAdcHitSum; - } - return *this; - } - PadHit &operator=(const AdcHit &r) - { - ((AdcHit *)this)->operator=(r); - return *this; - } - inline uint getPad() const { return _nPad; } - inline uint getTimeBin() const { return _nPad; } - inline void setAdcHitSum(float f) { _fAdcHitSum = f; } - inline float getAdcHitSum() const { return _fAdcHitSum; } - inline void AddAdc(float f) { _fAdcHitSum += f; } - inline void ReduceAdcSum(float f) { _fAdcHitSum -= f; } - std::ostringstream &PrintPad(std::ostringstream &oss) const - { - oss << "pad: " << _nPad << ", AdcSum: " << _fAdcHitSum << std::endl; - return oss; - } - std::ostringstream &PrintTime(std::ostringstream &oss) const - { - oss << "timebin: " << _nPad << ", AdcSum: " << _fAdcHitSum << std::endl; - return oss; - } -}; -} // namespace tpcClustering -#endif diff --git a/detectors/tpc/clusterHitFinder/fast/TpcRowClusters.h b/detectors/tpc/clusterHitFinder/fast/TpcRowClusters.h deleted file mode 100644 index 0e7093e0..00000000 --- a/detectors/tpc/clusterHitFinder/fast/TpcRowClusters.h +++ /dev/null @@ -1,276 +0,0 @@ -// -// Joint Institute for Nuclear Research (JINR), Dubna, Russia, Jun-2024 -// MPD TPC Clustering objects for NICA project -// author: Viktor A.Krylov JINR/LNP -// email: kryman@jinr.ru -// -#ifndef ROW_CLUSTERS_H -#define ROW_CLUSTERS_H - -#include "TpcClusterUnits.h" -#include "TpcCluster.h" - -namespace tpcClustering { -class RowClusters { - uint _nRow; // sector row number - uint _nClusterId; // current cluster id for internal usage - RowStat _RowStat; // row statistic - - std::list<Cluster *> _lstClusters; - void Joint(std::vector<std::list<PadCluster *>::iterator> &rviPadClusters, uint nPadCluster, PadCluster *pPadCluster) - { - Cluster *pCluster = NULL; - uint nTimeBinMin = 0; - uint nTimeBinMax = 0; - uint nPadMin = 0; - uint nPadMax = 0; - uint nPad = pPadCluster->getPad(); - std::list<PadCluster *>::iterator iPadCluster1st = rviPadClusters[nPadCluster - 1]; - std::list<PadCluster *>::iterator iPadCluster2nd = rviPadClusters[nPadCluster]; - for (std::list<PadCluster *>::iterator iPadCluster = iPadCluster1st; iPadCluster != iPadCluster2nd; - iPadCluster++) { - PadCluster *pPadCluster0 = (*iPadCluster); - Cluster *pCluster0 = pPadCluster0->getCluster(); - assert(pCluster0); - - if (pCluster0->isOut(nPad) || pPadCluster->isBefore(pPadCluster0)) - break; - else if (pPadCluster->isAfter(pPadCluster0)) - continue; - else if (pPadCluster->isBelongTo(pPadCluster0)) { - pPadCluster->setPadTimeOverlapping(pPadCluster0); - pCluster = pPadCluster->getCluster(); - if (pCluster == NULL) { - pCluster = pPadCluster0->getCluster(); - assert(pCluster != NULL); - - nTimeBinMin = pCluster->getTimeBinMin(); - nTimeBinMax = pCluster->getTimeBinMax(); - nPadMin = pCluster->getPadMin(); - nPadMax = pCluster->getPadMax(); - - pPadCluster->setCluster(pCluster); - pCluster->Add(pPadCluster); - } else if (pCluster0->getTimeBinMax() >= nTimeBinMin && pCluster0->getTimeBinMin() <= nTimeBinMax && - pCluster0->getPadMax() == - pCluster0->getPadMin()) { // if it is truth then we should add whole PadClusters list to the - // Cluster and remove Cluster0 - pCluster->Add(pCluster0); - pCluster0->Clear(); - _lstClusters.remove_if([pCluster0](Cluster *p) { return p == pCluster0; }); - delete pCluster0; - } - } - } - if (pCluster == NULL) { - pCluster = new Cluster(_nClusterId++, pPadCluster); - pPadCluster->setCluster(pCluster); - _lstClusters.push_back(pCluster); - } - } - -public: - std::list<PadCluster *> _lstPadClusters; // list of fired ADC hits - std::vector<std::list<PadCluster *>::iterator> _viPadClusters; - RowClusters() : _nRow(0), _nClusterId(0), _RowStat() {} - RowClusters(uint nRow) : _nRow(nRow), _nClusterId(0), _RowStat() {} - RowClusters(uint nRow, uint nPad, uint nTimeBin, float fAdc, int nTrackId) - : _nRow(nRow), _nClusterId(0), _RowStat(nPad, nTimeBin, fAdc) - { - AdcHit adcHit(nTimeBin, fAdc, nTrackId); // new ADC hit - _lstPadClusters.push_back(new PadCluster(nPad, adcHit)); - _viPadClusters.push_back(_lstPadClusters.begin()); - } - virtual ~RowClusters() - { - _lstPadClusters.clear(); - _viPadClusters.clear(); - while (!_lstClusters.empty()) { - delete _lstClusters.front(); - _lstClusters.pop_front(); - } - } - inline uint getRow() const { return _nRow; } - inline const std::list<Cluster *> &getClusters() const { return _lstClusters; } - inline bool isEmpty() const { return _lstClusters.empty(); } - inline void Load(uint nPad, uint nTimeBin, float fAdc, int nTrackId) - { - _RowStat.PickOut(nPad, nTimeBin, fAdc); - AdcHit adcHit(nTimeBin, fAdc, nTrackId); // new ADC hit - - assert(!_lstPadClusters.empty()); - - PadCluster *pPadCluster = _lstPadClusters.back(); - if (nPad != pPadCluster->getPad()) { // new PadCluster - _lstPadClusters.push_back(new PadCluster(nPad, adcHit)); - _viPadClusters.push_back(std::prev(_lstPadClusters.end())); - } else if (pPadCluster->isOut(adcHit)) // new cluster in the pad - _lstPadClusters.push_back(new PadCluster(nPad, adcHit)); - else if (pPadCluster->isCombined(adcHit)) { // combined PadCluster - PadCluster *pCombinedPadCluster = new PadCluster(nPad, pPadCluster->splitAdcHit(adcHit)); - pCombinedPadCluster->setPadClusterBefore(pPadCluster); - pPadCluster->setPadClusterAfter(pCombinedPadCluster); - pCombinedPadCluster->Add(adcHit); - _lstPadClusters.push_back(pCombinedPadCluster); - } else // simple PadCluster - pPadCluster->Add(adcHit); - } - void Joint(std::vector<std::list<PadCluster *>::iterator> &rviPadClusters) - { - assert(rviPadClusters.size() > 1); - size_t nPadSize = rviPadClusters.size() - 1; - for (uint i = 0; i < nPadSize; i++) { - std::list<PadCluster *>::iterator iPadCluster1st = rviPadClusters[i]; - std::list<PadCluster *>::iterator iPadCluster2nd = rviPadClusters[i + 1]; - for (std::list<PadCluster *>::iterator iPadCluster = iPadCluster1st; iPadCluster != iPadCluster2nd; - iPadCluster++) { - PadCluster *pPadCluster = (*iPadCluster); - if (i == 0) { // init pad clusters seed - Cluster *pCluster = new Cluster(_nClusterId++, pPadCluster); - pPadCluster->setCluster(pCluster); - _lstClusters.push_back(pCluster); - } else - Joint(rviPadClusters, i, pPadCluster); - } - } - } - void Joint() - { - assert(_viPadClusters.size()); - _viPadClusters.push_back(_lstPadClusters.end()); // add the end of the PadClusters before processing - size_t nPadSize = _viPadClusters.size() - 1; - for (uint i = 0; i < nPadSize; i++) { - std::list<PadCluster *>::iterator iPadCluster1st = _viPadClusters[i]; - std::list<PadCluster *>::iterator iPadCluster2nd = _viPadClusters[i + 1]; - for (std::list<PadCluster *>::iterator iPadCluster = iPadCluster1st; iPadCluster != iPadCluster2nd; - iPadCluster++) { - PadCluster *pPadCluster = (*iPadCluster); - if (i == 0) { // init pad clusters seed - Cluster *pCluster = new Cluster(_nClusterId++, pPadCluster); - pPadCluster->setCluster(pCluster); - _lstClusters.push_back(pCluster); - } else - Joint(_viPadClusters, i, pPadCluster); - } - } - } - void Split() - { - for (std::list<Cluster *>::iterator i = _lstClusters.begin(); i != _lstClusters.end(); ++i) { - const Cluster *pCluster = NULL; - // (*i)->WaveletTransform(); - (*i)->resetStat(); - // (*i)->Rotate(); - while (Cluster *pNewCluster = (*i)->Split(pCluster)) { - pNewCluster->setId(_nClusterId++); - _lstClusters.insert(i, pNewCluster); - pCluster = pNewCluster; - } - } - } - void Draw(std::ostringstream &oss) const - { - oss << "Draw:RowClusters:: row: " << _nRow << ' '; - _RowStat.Print(oss); - oss << std::endl; - for (std::list<Cluster *>::const_iterator i = _lstClusters.begin(); i != _lstClusters.end(); ++i) { - (*i)->Draw(oss); - oss << std::endl; - } - } - void DrawPadClusters(std::ostringstream &oss) const - { - oss << "DrawPadClusters:RowClusters:: row: " << _nRow << ' '; - _RowStat.Print(oss); - oss << std::endl; - uint nTimeBinSize = _RowStat._nTimeBinMax / 10 + 1; - DrawTimeBinTopScale(oss, nTimeBinSize); - - assert(_viPadClusters.size() > 1); - size_t nPadSize = _viPadClusters.size() - 1; - for (uint i = 0; i < nPadSize; i++) { - std::list<PadCluster *>::iterator iPadCluster = _viPadClusters[i]; - oss << std::setw(3) << (*iPadCluster)->getPad() << ":"; - uint nPos = 0; - for (; iPadCluster != _viPadClusters[i + 1]; iPadCluster++) { - PadCluster *pPadCluster = *iPadCluster; - uint nMin = pPadCluster->getTimeBinMin(); - if (nMin < nPos) { // shift back to the 'n' positions (for combined PadClusters!!!) - long int nShift = nPos - nMin; - oss.seekp(oss.tellp() - nShift); - nPos -= nShift; - } - uint nGap = nMin - nPos; - if (nGap) { - oss << std::string(nGap, '-'); - nPos += nGap; - } - char ch = - pPadCluster->getCluster() ? SZ_ALPHABET[pPadCluster->getCluster()->getId() % strlen(SZ_ALPHABET)] : '+'; - pPadCluster->Draw(oss, ch); - nPos += pPadCluster->getSize(); - } - oss << std::endl; - } - DrawTimeBinBottomScale(oss, nTimeBinSize); - } - void DrawRowClusters(std::ostringstream &oss) const - { - oss << "Draw RowClusters:: row: " << _nRow << ' '; - _RowStat.Print(oss); - oss << std::endl; - uint nTimeBinSize = _RowStat._nTimeBinMax / 10 + 1; - DrawTimeBinTopScale(oss, nTimeBinSize); - for (std::list<Cluster *>::const_iterator i = _lstClusters.begin(); i != _lstClusters.end(); i++) - (*i)->DrawPadClusters(oss); - DrawTimeBinBottomScale(oss, nTimeBinSize); - } - void Print(std::ostringstream &oss, std::vector<std::list<PadCluster *>::iterator> &rviPadClusters) - { - oss << "RowClusters:: row: " << _nRow << ' '; - _RowStat.Print(oss); - oss << std::endl; - if (rviPadClusters.size() > 1) { - size_t nSize1 = rviPadClusters.size() - 1; - for (uint i = 0; i < nSize1; i++) { - std::list<PadCluster *>::iterator iPadCluster = rviPadClusters[i]; - std::list<PadCluster *>::iterator iPadClusterLast = rviPadClusters[i + 1]; - oss << " ---------------- id: " << i << std::endl; - for (; iPadCluster != iPadClusterLast; iPadCluster++) { - oss << "\t"; - (*iPadCluster)->PrintPad(oss); - } - } - } - oss << " ----------------" << std::endl; - } - void PrintPadClusters(std::ostringstream &oss) - { - oss << "PrintPadClusters:RowClusters:: row: " << _nRow << ' '; - _RowStat.Print(oss); - oss << std::endl; - if (_viPadClusters.size() > 1) { - size_t nSize1 = _viPadClusters.size() - 1; - for (uint i = 0; i < nSize1; i++) { - std::list<PadCluster *>::iterator iPadCluster = _viPadClusters[i]; - std::list<PadCluster *>::iterator iPadClusterLast = _viPadClusters[i + 1]; - oss << " ---------------- id: " << i << std::endl; - for (; iPadCluster != iPadClusterLast; iPadCluster++) { - oss << "\t"; - (*iPadCluster)->PrintPad(oss); - } - } - } - oss << " ----------------" << std::endl; - } - void Print(std::ostringstream &oss) const - { - oss << "Print:RowClusters:: row: " << getRow() << ' '; - _RowStat.Print(oss); - oss << std::endl; - for (std::list<Cluster *>::const_iterator i = _lstClusters.begin(); i != _lstClusters.end(); ++i) - (*i)->Print(oss); - } -}; -} // namespace tpcClustering -#endif diff --git a/detectors/tpc/clusterHitFinder/fast/TpcSectorClusters.h b/detectors/tpc/clusterHitFinder/fast/TpcSectorClusters.h deleted file mode 100644 index c79656cf..00000000 --- a/detectors/tpc/clusterHitFinder/fast/TpcSectorClusters.h +++ /dev/null @@ -1,142 +0,0 @@ -// -// Joint Institute for Nuclear Research (JINR), Dubna, Russia, Jun-2024 -// MPD TPC Clustering objects for NICA project -// author: Viktor A.Krylov JINR/LNP -// email: kryman@jinr.ru -// -#ifndef SECTOR_CLUSTERS_H -#define SECTOR_CLUSTERS_H - -#include "TpcRowClusters.h" - -#define nTPC_ROWS 53 - -namespace tpcClustering { - -class SectorClusters { - int _nSector; - -public: - std::list<RowClusters *> _lstRowClusters; - SectorClusters() : _nSector(0) {} - SectorClusters(int nSector) : _nSector(nSector) {} - SectorClusters(int nSector, uint nRow, uint nPad, uint nTimeBin, float fAdc, int nTrackId) : _nSector(nSector) - { - _lstRowClusters.push_back(new RowClusters(nRow, nPad, nTimeBin, fAdc, nTrackId)); - } - virtual ~SectorClusters() - { - while (!_lstRowClusters.empty()) { - delete _lstRowClusters.front(); - _lstRowClusters.pop_front(); - } - } - inline int getSector() const { return _nSector; } - inline void setSector(int nSector) { _nSector = nSector; } - inline bool isEmpty() const { return _lstRowClusters.empty(); } - inline std::size_t getSize() const { return _lstRowClusters.size(); } -#ifdef PARALLEL - void DefineClusters() - { - for (std::list<RowClusters *>::const_iterator iRowCLusters = _lstRowClusters.begin(); - iRowCLusters != _lstRowClusters.end(); iRowCLusters++) { - (*iRowCLusters)->Joint(); - (*iRowCLusters)->Split(); - } - } -#else - void DefineClusters(std::ostringstream &oss) - { - for (std::list<RowClusters *>::const_iterator iRowCLusters = _lstRowClusters.begin(); - iRowCLusters != _lstRowClusters.end(); iRowCLusters++) { - -#ifdef DEBUG_OUTPUT - oss << "Before joint" << std::endl; - // DrawPadClusters(oss); - // PrintPadClusters(oss); -#endif - (*iRowCLusters)->Joint(); -#ifdef DEBUG_OUTPUT - oss << "Before split" << std::endl; - DrawPadClusters(oss); - // Draw(oss); - // Print(oss); - // return; -#endif - (*iRowCLusters)->Split(); -#ifdef DEBUG_OUTPUT - oss << "After split" << std::endl; - DrawPadClusters(oss); - // DrawRowClusters(oss); - Draw(oss); - Print(oss); -#endif - } - } -#endif - inline void Add(RowClusters *pRowClusters) { _lstRowClusters.push_back(pRowClusters); } - inline void Load(uint nRow, uint nPad, uint nTimeBin, float fAdc, int nTrackId) - { - if (_lstRowClusters.empty()) - _lstRowClusters.push_back(new RowClusters(nRow, nPad, nTimeBin, fAdc, nTrackId)); - else { - RowClusters *pRowClusters = _lstRowClusters.back(); - if (pRowClusters->getRow() != nRow) - _lstRowClusters.push_back(new RowClusters(nRow, nPad, nTimeBin, fAdc, nTrackId)); - else - pRowClusters->Load(nPad, nTimeBin, fAdc, nTrackId); - } - } - void Type(std::ostringstream &oss) const - { - std::string strRows = std::string(nTPC_ROWS, '-'); - for (std::list<RowClusters *>::const_iterator i = _lstRowClusters.begin(); i != _lstRowClusters.end(); i++) { - uint nRow = (*i)->getRow(); - strRows[nRow] = szDECADE[nRow % 10]; - } - oss << " Rows: " << strRows << std::endl; - } - void Print(std::ostringstream &oss) const - { - oss << "SectorClusters:: sector: " << getSector() << std::endl; - for (std::list<RowClusters *>::const_iterator i = _lstRowClusters.begin(); i != _lstRowClusters.end(); i++) - (*i)->Print(oss); - } - void PrintPadClusters(std::ostringstream &oss) const - { - oss << "SectorClusters:: sector: " << getSector() << std::endl; - for (std::list<RowClusters *>::const_iterator i = _lstRowClusters.begin(); i != _lstRowClusters.end(); i++) - (*i)->PrintPadClusters(oss); - } - void Draw(std::ostringstream &oss) const - { - oss << "SectorClusters:: sector: " << getSector() << std::endl; - for (std::list<RowClusters *>::const_iterator i = _lstRowClusters.begin(); i != _lstRowClusters.end(); i++) - (*i)->Draw(oss); - } - void DrawPadClusters(std::ostringstream &oss) const - { - oss << "SectorClusters:: sector: " << getSector() << std::endl; - for (std::list<RowClusters *>::const_iterator i = _lstRowClusters.begin(); i != _lstRowClusters.end(); i++) - (*i)->DrawPadClusters(oss); - } - void DrawRowClusters(std::ostringstream &oss) const - { - oss << "SectorClusters:: sector: " << getSector() << std::endl; - for (std::list<RowClusters *>::const_iterator i = _lstRowClusters.begin(); i != _lstRowClusters.end(); i++) - (*i)->DrawRowClusters(oss); - } -}; - -#undef PARALLEL -#ifdef PARALLEL - -void *ThreadFunc(void *pArgs) -{ - ((SectorClusters *)pArgs)->DefineClusters(); - pthread_exit(NULL); -} - -#endif -} // namespace tpcClustering -#endif diff --git a/detectors/tpc/clusterHitFinder/fast/TpcSectorObj.cxx b/detectors/tpc/clusterHitFinder/fast/TpcSectorObj.cxx new file mode 100644 index 00000000..0b739818 --- /dev/null +++ b/detectors/tpc/clusterHitFinder/fast/TpcSectorObj.cxx @@ -0,0 +1,401 @@ +// +// Joint Institute for Nuclear Research (JINR), Dubna, Russia, Mar-2025 +// TPC Sector class for MPD/NICA project +// author: Viktor Krylov JINR/LHEP +// email: kryman@jinr.ru +// +#include <vector> + +#include <TKey.h> +#include <TH1F.h> +#include <TH2F.h> +#include <TH2I.h> +#include <TH3S.h> +#include <TGraph2D.h> +#include <TDirectoryFile.h> + +#include "TpcSectorObj.h" +#include "clustering/TpcClusters.h" +#include "clustering/TpcTrack.h" +#include "clustering/TpcSectorGeo.h" +#include "clustering/TpcSampaPack.h" + +extern const char *szTPC_MC_SAMPA; +extern const char *szTPC_DIGITS; + +bool bH1FXYZ_DIGIT_DEV = false; // histogram TH1F MC track / digit deviation +bool bH1FXYZ_HIT_DEV = false; // histogram TH1F MC track / hit deviation +bool bH2FXYZ_HIT_DEV = false; // histogram TH1F MC track / hit deviation +bool bH2I_H3S_DIGITS = false; // graph TGraph2D adc hit distribution +bool bG2D_HITS = false; // graph TGraph2D adc hit distribution + +using namespace TpcClustering; + +TpcSectorObj::~TpcSectorObj() +{ + delete[] _pSampaBuf; + _vTpcAdcHits.clear(); + _pvTpcDigits->clear(); + delete _pvTpcDigits; + delete _pTpcClusters; +} + +bool TpcSectorObj::loadSampaBuf(TDirectory *pDir) +{ + assert(pDir); + + TKey *pKey = pDir->FindKey(szTPC_MC_SAMPA); + if (pKey) { + const std::vector<uint8_t> *pv; + pDir->GetObject(pKey->GetName(), pv); + if (pv) { + _nSampaSize = pv->size(); + _pSampaBuf = pv->data(); + return true; + } else + std::cout << "Error:: sampa: " << szTPC_MC_SAMPA << " has wrong object: " << pKey->GetName() << std::endl; + } else + std::cout << "Warning:: sector:" << _nSector << " has not SAMPA buffer!" << std::endl; + return false; +} +void TpcSectorObj::loadSampaBuf() +{ + assert(!_pTpcClusters); + + _pTpcClusters = new TpcClusters(_nSector); + const std::uint8_t *pLast = _pSampaBuf + _nSampaSize; + const std::uint8_t *pi = _pSampaBuf; + while (pi < pLast) { + TpcSampaPack tsp; + tsp.HeaderDecode(pi); + pi = tsp.PayloadDecode(pi); + if (!tsp.isEmpty()) _pTpcClusters->Load(tsp._nRow, tsp._nPad, tsp._TpcAdcTimebins, tsp._nSize); + } + delete _pSampaBuf; + _pSampaBuf = nullptr; + _nSampaSize = 0; +}; +void TpcSectorObj::createSampaBuf() +{ + assert(_pvTpcDigits && !_pSampaBuf); + + _nSampaSize = nTPC_SECTOR_PADS * nSAMPA_BUFFER_SIZE; + // std::uint8_t* pSampaBuf = new std::uint8_t[_nSampaSize]{0}; + // assert(pSampaBuf); + + // memset(pSampaBuf, 0, _nSampaSize); + uint nRow = 0; + uint nPad = 0; + uint nTimeBin = 0; + std::uint16_t nAdc = 0; + // std::uint8_t* pPadBuf = (anTPC_ROW_PAD_SHIFT[nRow] + nPad) * nSAMPA_BUFFER_SIZE; + // TpcSampaPack tsp; + // std::uint8_t* pPadBuf = pSampaBuf; + // for( const TpcDigit& rTpcDigit : *_pvTpcDigits ) { + // if( rTpcDigit._nRow == nRow && rTpcDigit._nPad == nPad ) { + // assert(rTpcDigit._nAdc); + + // if( rTpcDigit._nTimeBin == nTimeBin ) + // nAdc += rTpcDigit._nAdc; + // else { + // if( pPadBuf == NULL ) + // pPadBuf = pSampaBuf + (anTPC_ROW_PAD_SHIFT[nRow] + nPad) * nSAMPA_BUFFER_SIZE; + // tsp.EncodeAdc(pPadBuf, nTimeBin, nAdc); + // nAdc = rTpcDigit._nAdc; + // } + // } else { + // if( pPadBuf != NULL ) { + // tsp.EncodeAdc(pPadBuf, nTimeBin, nAdc); + // } + + // } + // } + // if( pPadBuf == NULL ) + // pPadBuf = pSampaBuf + (anTPC_ROW_PAD_SHIFT[nRow] + nPad) * nPadSize; + // std::uint8_t* pPadBuf = anTPC_ROW_PAD_SHIFT[nRow] + nPad; + + // tsp.EncodeAdc(pSampaBuf, rTpcDigit._nTimeBin, rTpcDigit._nAdc); + // std::cout << "encode row:" << tpcDigit._nRow << " pad:" << rTpcDigit._nPad << " timebin:" << + // rTpcDigit._nTimeBin << " adc:" << rTpcDigit._nAdc << std::endl; + // } else { + // for( uint i = 0; i < nSAMPA_PAYLOAD_WORDS; i++ ) { + // uint16_t nAdc = tsp.DecodeAdc(pSampaBuf, i); + // if( nAdc ) + // std::cout << "decode timebin:" << i << " adc:" << nAdc << std::endl; + // } + // tpcDigit = rTpcDigit; + // memset(pSampaBuf, 0, nSampaSize); + // } + // } +} +bool TpcSectorObj::loadDigits(TDirectory *pDir) +{ + assert(pDir); + + TKey *pKey = pDir->FindKey(szTPC_DIGITS); + if (pKey) { + std::vector<uint8_t> *pv; + pDir->GetObject(pKey->GetName(), pv); + if (pv) { + _pvTpcDigits = (std::vector<TpcDigit> *)pv; + return true; + } else + std::cout << "Error:: tpc digits: " << szTPC_DIGITS << " has wrong object: " << pKey->GetName() << std::endl; + } else + std::cout << "Warning:: sector:" << _nSector << " has not 'digits' folder!" << std::endl; + return false; +} +void TpcSectorObj::loadDigits() +{ + assert(!_pTpcClusters && _pvTpcDigits); + + _pTpcClusters = new TpcClusters(_nSector); + for (const TpcDigit &rTpcDigit : *_pvTpcDigits) + if (rTpcDigit._nAdc != 0) + _pTpcClusters->Load(rTpcDigit._nRow, rTpcDigit._nPad, rTpcDigit._nTimeBin, rTpcDigit._nAdc); +} +void TpcSectorObj::defAdcHits() +{ + for (TpcRowClusters *pTpcRowClusters : _pTpcClusters->_lstTpcRowClusters) { + pTpcRowClusters->defAdcHits(_nSector, _vTpcAdcHits); + } +} +void TpcSectorObj::McStore(TDirectory *pDir) const +{ + TDirectory *pSectorDir = pDir->mkdir(std::to_string(_nSector).c_str(), "TPC Sector Digits", kTRUE); + if (pSectorDir && pSectorDir->cd() && _pvTpcDigits) { + const TpcDigit *pTpcDigit = _pvTpcDigits->data(); + std::vector<uint8_t> vDigits((uint8_t *)pTpcDigit, + (uint8_t *)pTpcDigit + (sizeof(TpcDigit) * _pvTpcDigits->size())); + pSectorDir->WriteObject(&vDigits, szTPC_DIGITS); + } else + std::cout << "Error:: TpcSector::Store: the folder: '" << _nSector << "' cannot be found or created!" + << std::endl; +} +void TpcSectorObj::StoreG2D(TDirectory *pDir) const +{ + if (_vTpcAdcHits.size()) { + std::string strSector = std::to_string(_nSector); + TGraph2D g2d(_vTpcAdcHits.size()); + g2d.SetTitle(("ADC Cluster Hits sector:" + strSector).c_str()); + g2d.GetXaxis()->SetTitle("X"); + g2d.GetYaxis()->SetTitle("Z"); + g2d.GetZaxis()->SetTitle("Y"); + uint nPoints = 0; + for (const TpcAdcHit &rTpcAdcHit : _vTpcAdcHits) + g2d.SetPoint(nPoints++, rTpcAdcHit._fX, rTpcAdcHit._fZ, rTpcAdcHit._fY); + if (nPoints) pDir->WriteObject(&g2d, ("G2D:" + strSector).c_str()); + } +} +void TpcSectorObj::StoreH2I_H3S(TDirectory *pDir) const +{ + if (_pvTpcDigits && _pvTpcDigits->size()) { + std::string strSector = std::to_string(_nSector); + TH2I h2i(("H2I:" + strSector).c_str(), ("TPC Sector:" + strSector).c_str(), nTPC_ROW_MAX_PADS, 0, + nTPC_ROW_MAX_PADS, nTPC_SECTOR_ROWS, 0, nTPC_SECTOR_ROWS); + h2i.GetXaxis()->SetTitle("Pads[0-123]"); + h2i.GetYaxis()->SetTitle("Rows[0-52]"); + h2i.GetZaxis()->SetTitle("ADC"); + + TH3S h3s(("H3S:" + strSector).c_str(), ("TPC Sector:" + strSector).c_str(), nTPC_ROW_MAX_PADS, 0, + nTPC_ROW_MAX_PADS, nTPC_TIMEBINS, 0, nTPC_TIMEBINS, nTPC_SECTOR_ROWS, 0, nTPC_SECTOR_ROWS); + h3s.GetXaxis()->SetTitle("Pads[0-123]"); + h3s.GetYaxis()->SetTitle("Timebins[0-310]"); + h3s.GetZaxis()->SetTitle("Rows[0-52]"); + for (const TpcDigit &rTpcDigit : *_pvTpcDigits) { + uint16_t nPad = rTpcDigit._nPad + 1; + uint16_t nRow = rTpcDigit._nRow + 1; + uint16_t nTimeBin = rTpcDigit._nTimeBin + 1; + h2i.SetBinContent(nPad, nRow, h2i.GetBinContent(nPad, nRow) + rTpcDigit._nAdc); + h3s.SetBinContent(nPad, nTimeBin, nRow, h3s.GetBinContent(nPad, nTimeBin, nRow) + rTpcDigit._nAdc); + } + if (h2i.GetEntries()) pDir->WriteObject(&h2i, ("H2I:" + strSector).c_str()); + if (h3s.GetEntries()) pDir->WriteObject(&h3s, ("H3S:" + strSector).c_str()); + } +} +void TpcSectorObj::StoreH1FxyzDigits(TDirectory *pDir, const std::vector<TpcTrack *> &rvpTpcTracks) const +{ + assert(_pvTpcDigits); + + TH1F h1fx("H1Fx:", "TPC MC Track 'X' Digits Dispersion:", 250, -5.0f, 5.0f); + h1fx.GetXaxis()->SetTitle("X cm"); + h1fx.GetYaxis()->SetTitle("Digits"); + h1fx.SetFillColor(45); + TH1F h1fy("H1Fy:", "TPC MC Track 'Y' Digits Dispersion:", 250, -5.0f, 5.0f); + h1fy.GetXaxis()->SetTitle("Y cm"); + h1fy.GetYaxis()->SetTitle("Digits"); + h1fy.SetFillColor(45); + TH1F h1fz("H1Fz:", "TPC MC Track 'Z' Digits Dispersion:", 250, -5.0f, 5.0f); + h1fz.GetXaxis()->SetTitle("Z cm"); + h1fz.GetYaxis()->SetTitle("Digits"); + h1fz.SetFillColor(45); + + for (const TpcDigit &rTpcDigit : *_pvTpcDigits) { // A + TpcPoint tpcLocal(rTpcDigit._nRow, rTpcDigit._nPad, rTpcDigit._nTimeBin); + TpcPoint tpcGlobal(tpcLocal.getGlobal(_nSector)); + uint nTrackId = rTpcDigit._nTrack; + std::vector<TpcTrack *>::const_iterator iTpcTrack = std::find_if( + rvpTpcTracks.begin(), rvpTpcTracks.end(), [nTrackId](const TpcTrack *p) { return p->_nTrackId == nTrackId; }); + if (iTpcTrack != rvpTpcTracks.end()) { + TpcTrack *pTpcTrack = *iTpcTrack; + size_t nSize = pTpcTrack->_pvTpcPoints->size(); + if (nSize > 1) { + float fDev = FLT_MAX; // deviation + TpcPoint tpcNormal; + uint n = 0; + for (uint i = 1; i < nSize; i++) { + const TpcPoint &rTpcPoint1 = pTpcTrack->_pvTpcPoints->operator[](i - 1); // B + const TpcPoint &rTpcPoint2 = pTpcTrack->_pvTpcPoints->operator[](i); // C + TpcPoint tpcN = tpcGlobal.getNormal(rTpcPoint1, rTpcPoint2); + float f = + tpcN.isDefined() ? sqrt(tpcN._fX * tpcN._fX + tpcN._fY * tpcN._fY + tpcN._fZ * tpcN._fZ) : FLT_MAX; + if (f != FLT_MAX) { + if (fDev > f) { + fDev = f; + tpcNormal = tpcN; + n = i; + } + } else { + std::cout << "Wrong track points! track id:" << pTpcTrack->_nTrackId << std::endl + << " point[" << i - 1 << "]:" << rTpcPoint1 << std::endl + << " point[" << i << "]:" << rTpcPoint2 << std::endl; + } + } + assert(n); + + h1fx.Fill(tpcNormal._fX); + h1fy.Fill(tpcNormal._fY); + h1fz.Fill(tpcNormal._fZ); + } + } + } + if (h1fx.GetEntries()) pDir->WriteObject(&h1fx, "H1FxD:"); + if (h1fy.GetEntries()) pDir->WriteObject(&h1fy, "H1FyD:"); + if (h1fz.GetEntries()) pDir->WriteObject(&h1fz, "H1FzD:"); +} + +void TpcSectorObj::StoreH1FxyzHits(TDirectory *pDir, const std::vector<TpcTrack *> &rvpTpcTracks) const +{ + assert(_pvTpcDigits); + + TH1F h1fx("H1Fx:", "TPC MC Track 'X' Hit Dispersion:", 250, -2.5f, 2.5f); + h1fx.GetXaxis()->SetTitle("X cm"); + h1fx.GetYaxis()->SetTitle("Hits"); + h1fx.SetFillColor(45); + TH1F h1fy("H1Fy:", "TPC MC Track 'Y' Hit Dispersion:", 250, -2.5f, 2.5f); + h1fy.GetXaxis()->SetTitle("Y cm"); + h1fy.GetYaxis()->SetTitle("Hits"); + h1fy.SetFillColor(45); + TH1F h1fz("H1Fz:", "TPC MC Track 'Z' Hit Dispersion:", 250, -2.5f, 2.5f); + h1fz.GetXaxis()->SetTitle("Z cm"); + h1fz.GetYaxis()->SetTitle("Hits"); + h1fz.SetFillColor(45); + + for (const TpcAdcHit &rTpcAdcHit : _vTpcAdcHits) { // A + uint nRow = rTpcAdcHit._nRow; + uint nClusterId = rTpcAdcHit._nClusterId; + TpcPoint tpcHitPoint = {rTpcAdcHit._fX, rTpcAdcHit._fY, rTpcAdcHit._fZ}; + const TpcCluster *pTpcCluster = _pTpcClusters->findTpcCluster(rTpcAdcHit._nRow, rTpcAdcHit._nClusterId); + if (pTpcCluster && _pvTpcDigits) { + for (uint nTpcDigitId : pTpcCluster->_vnTpcDigits) { + uint nTrackId = _pvTpcDigits->operator[](nTpcDigitId)._nTrack; + std::vector<TpcTrack *>::const_iterator iTpcTrack = + std::find_if(rvpTpcTracks.begin(), rvpTpcTracks.end(), + [nTrackId](const TpcTrack *p) { return p->_nTrackId == nTrackId; }); + if (iTpcTrack != rvpTpcTracks.end()) { + TpcTrack *pTpcTrack = *iTpcTrack; + size_t nSize = pTpcTrack->_pvTpcPoints->size(); + if (nSize > 1) { + float fDev = FLT_MAX; // min hit deviation + TpcPoint tpcNormal; + uint n = 0; + for (uint i = 1; i < nSize; i++) { + const TpcPoint &rTpcPoint1 = pTpcTrack->_pvTpcPoints->operator[](i - 1); // B + const TpcPoint &rTpcPoint2 = pTpcTrack->_pvTpcPoints->operator[](i); // C + TpcPoint tpcN = tpcHitPoint.getNormal(rTpcPoint1, rTpcPoint2); + float f = tpcN.isDefined() ? sqrt(tpcN._fX * tpcN._fX + tpcN._fY * tpcN._fY + tpcN._fZ * tpcN._fZ) + : FLT_MAX; + if (f != FLT_MAX) { + if (fDev > f) { + fDev = f; + tpcNormal = tpcN; + n = i; + } + } else { + std::cout << "Wrong track points! track id:" << pTpcTrack->_nTrackId << std::endl + << " point[" << i - 1 << "]:" << rTpcPoint1 << std::endl + << " point[" << i << "]:" << rTpcPoint2 << std::endl; + } + } + assert(n); + + h1fx.Fill(tpcNormal._fX); + h1fy.Fill(tpcNormal._fY); + h1fz.Fill(tpcNormal._fZ); + break; + } + } + } + } + } + if (h1fx.GetEntries()) pDir->WriteObject(&h1fx, "H1FxH:"); + if (h1fy.GetEntries()) pDir->WriteObject(&h1fy, "H1FyH:"); + if (h1fz.GetEntries()) pDir->WriteObject(&h1fz, "H1FzH:"); +} +void TpcSectorObj::StoreH2FxyzHits(TDirectory *pDir, const std::vector<TpcTrack *> &rvpTpcTracks) const +{ + assert(_pvTpcDigits); + + size_t nTracks = rvpTpcTracks.size(); + TH2F h2fx("H2Fx:", "TPC MC Track 'X' Hit Dispersion:", 250, -2.5f, 2.5f, nTracks, 0, nTracks); + h2fx.GetXaxis()->SetTitle("X cm"); + h2fx.GetYaxis()->SetTitle("Track ID"); + h2fx.GetZaxis()->SetTitle("Hits"); + h2fx.SetFillColor(45); + TH2F h2fy("H1Fy:", "TPC MC Track 'Y' Hit Dispersion:", 250, -2.5f, 2.5f, nTracks, 0, nTracks); + h2fy.GetXaxis()->SetTitle("Y cm"); + h2fy.GetYaxis()->SetTitle("Track ID"); + h2fy.GetZaxis()->SetTitle("Hits"); + h2fy.SetFillColor(45); + TH2F h2fz("H1Fz:", "TPC MC Track 'Z' Hit Dispersion:", 250, -2.5f, 2.5f, nTracks, 0, nTracks); + h2fz.GetXaxis()->SetTitle("Z cm"); + h2fz.GetYaxis()->SetTitle("Track ID"); + h2fz.GetZaxis()->SetTitle("Hits"); + h2fz.SetFillColor(45); + + for (const TpcTrack *pTpcTrack : rvpTpcTracks) { + uint32_t nTrackId = pTpcTrack->getTrackId(); + std::vector<TpcDigit>::const_iterator iTpcDigitBegin = _pvTpcDigits->begin(); + std::vector<TpcDigit>::const_iterator iTpcDigit = iTpcDigitBegin; + std::vector<TpcDigit>::const_iterator iTpcDigitEnd = _pvTpcDigits->end(); + while ((iTpcDigit = std::find_if(iTpcDigit, iTpcDigitEnd, [nTrackId](const TpcDigit &r) { + return r._nTrack == nTrackId; + })) != iTpcDigitEnd) { + uint nTpcDigit = std::distance(iTpcDigitBegin, iTpcDigit++); + const TpcDigit &rTpcDigit = _pvTpcDigits->operator[](nTpcDigit); + // const TpcCluster* pTpcCluster = _pTpcClusters->findTpcCluster(rTpcDigit._nRow, rTpcAdcHit._nClusterId); + } + } + if (h2fx.GetEntries()) pDir->WriteObject(&h2fx, "H2FxH:"); + if (h2fy.GetEntries()) pDir->WriteObject(&h2fy, "H2FyH:"); + if (h2fz.GetEntries()) pDir->WriteObject(&h2fz, "H2FzH:"); +} +void TpcSectorObj::Store(TDirectory *pDir, const std::vector<TpcTrack *> &rvpTpcTracks) const +{ + std::string strSector = std::to_string(_nSector); + TDirectory *pSectorDir = pDir->mkdir(strSector.c_str(), "TPC Sector Clusters", kTRUE); + if (pSectorDir && pSectorDir->cd()) { + const TpcAdcHit *pTpcAdcHit = _vTpcAdcHits.data(); + std::vector<uint8_t> vHits((uint8_t *)pTpcAdcHit, + (uint8_t *)pTpcAdcHit + (sizeof(pTpcAdcHit) * _vTpcAdcHits.size())); + pSectorDir->WriteObject(&vHits, szTPC_HITS); + if (_pTpcClusters) _pTpcClusters->Store(pSectorDir); + if (bH1FXYZ_DIGIT_DEV) StoreH1FxyzDigits(pSectorDir, rvpTpcTracks); + if (bH1FXYZ_HIT_DEV) StoreH1FxyzHits(pSectorDir, rvpTpcTracks); + if (bH2FXYZ_HIT_DEV) StoreH2FxyzHits(pSectorDir, rvpTpcTracks); + if (bH2I_H3S_DIGITS) StoreH2I_H3S(pSectorDir); + if (bG2D_HITS) StoreG2D(pSectorDir); // TGraph2D does not work in JSRoot File Viewer! + } else + std::cout << "Error:: TpcSector::Store: the folder: '" << _nSector << "' cannot be found or created!" + << std::endl; +} \ No newline at end of file diff --git a/detectors/tpc/clusterHitFinder/fast/TpcSectorObj.h b/detectors/tpc/clusterHitFinder/fast/TpcSectorObj.h new file mode 100644 index 00000000..cfe2309a --- /dev/null +++ b/detectors/tpc/clusterHitFinder/fast/TpcSectorObj.h @@ -0,0 +1,77 @@ +// +// Joint Institute for Nuclear Research (JINR), Dubna, Russia, Mar-2025 +// TPC Sector class for MPD/NICA project +// author: Viktor Krylov JINR/LHEP +// email: kryman@jinr.ru +// +#ifndef TPC_SECTOR_H +#define TPC_SECTOR_H + +#include "clustering/TpcEventPayload.h" +#include "clustering/TpcAdcHit.h" +#include "clustering/TpcSampaPack.h" +// #include "clustering/TpcTrack.h" +#include "clustering/TpcClusters.h" + +class TDirectory; + +struct TpcSectorObj { + uint _nSector; // sector number + std::streampos _nSampaSize; // _pSampa size: 1624240 + const std::uint8_t *_pSampaBuf; // size: 1624240 | pads: 4112 | pad size: 395 bytes == 3160 bits + std::vector<TpcClustering::TpcDigit> *_pvTpcDigits; // digit vector + TpcClustering::TpcClusters *_pTpcClusters; // row clusters + std::vector<TpcClustering::TpcAdcHit> _vTpcAdcHits; // cluster's adc hits + + TpcSectorObj() : _nSector(0), _nSampaSize(0), _pSampaBuf(nullptr), _pvTpcDigits(nullptr), _pTpcClusters(nullptr) {} + TpcSectorObj(uint nSector) + : _nSector(nSector), _nSampaSize(0), _pSampaBuf(nullptr), _pvTpcDigits(nullptr), _pTpcClusters(nullptr) + { + } + virtual ~TpcSectorObj(); + bool loadDigits(TDirectory *pDir); // load digits from root file to '_pvTpcDigits' + void loadDigits(); // load digits data to '_pTpcClusters' + bool loadSampaBuf(TDirectory *pDir); // load sampa buffer from root file to '_pSampa' + void loadSampaBuf(); // decode and load sampa data to '_pTpcClusters' (release _pSampa buffer!) + void createSampaBuf(); // create SAMPA buf for the sector to store it in a root file + // void defClusters(TDirectory* pDir); + // void defClusters(); // for debug only! + // void defDigits(); // for debug tracking only! + void defAdcHits(); + void McStore(TDirectory *pDir) const; // store mctracks to root file + void Store(TDirectory *pDir, const std::vector<TpcClustering::TpcTrack *> &rvpTpcTracks) const; + void StoreG2D(TDirectory *pDir) const; + void StoreH2I_H3S(TDirectory *pDir) const; + void StoreH1FxyzDigits(TDirectory *pDir, const std::vector<TpcClustering::TpcTrack *> &rvpTpcTracks) const; + void StoreH1FxyzHits(TDirectory *pDir, const std::vector<TpcClustering::TpcTrack *> &rvpTpcTracks) const; + void StoreH2FxyzHits(TDirectory *pDir, const std::vector<TpcClustering::TpcTrack *> &rvpTpcTracks) const; + // TpcClusters* getClusters(); + // void defClusters( TpcClusters* pClusters, uint nRow ); + // inline bool operator<( const TpcSector* pRight) const { return GetSector() < pRight->GetSector(); } + // std::uint16_t GetNum( char chLeft, char chRight ) const { + // assert( _strName.size() ); + + // std::size_t nLeft = _strName.find_last_of(chLeft); + // assert( nLeft != std::string::npos ); + // nLeft++; + + // std::size_t nRight = _strName.find_last_of(chRight); + // assert( nRight != std::string::npos && nRight > nLeft ); + + // return std::stoi(_strName.substr(nLeft, nRight - nLeft)); + // } + // inline std::uint16_t GetSector() const { return GetNum('_', '.'); } + void Decode(TpcClustering::TpcSectorPayload *pTpcSectorPayload) const + { + const std::uint8_t *pLast = _pSampaBuf + _nSampaSize; + const std::uint8_t *pi = _pSampaBuf; + while (pi < pLast) { + TpcClustering::TpcSampaPack tsp; + tsp.HeaderDecode(pi); + pi = tsp.PayloadDecode(pi); + if (!tsp.isEmpty()) pTpcSectorPayload->Push(tsp); + } + } +}; + +#endif // TPC_SECTOR_H diff --git a/detectors/tpc/clusterHitFinder/fast/clustering/TpcAdcGauss.h b/detectors/tpc/clusterHitFinder/fast/clustering/TpcAdcGauss.h new file mode 100644 index 00000000..7b8a3d09 --- /dev/null +++ b/detectors/tpc/clusterHitFinder/fast/clustering/TpcAdcGauss.h @@ -0,0 +1,317 @@ +// Joint Institute for Nuclear Research (JINR), Dubna, Russia, Jan-2024 +// TPC ADC templates and functions of gaussian normal distribution for MPD/NICA project +// author: Viktor A.Krylov JINR/LHEP +// email: kryman@jinr.ru +// +#ifndef TPC_ADC_GAUSS_H +#define TPC_ADC_GAUSS_H + +#include <iostream> +#include <array> +#include <vector> +#include <math.h> + +namespace TpcClustering { + +// +// A Fast, Accurate, and Separable Method for Fitting a Gaussian Function +// Ibrahim Al-Nahhal, Octavia A. Dobre, Ertugrul Basar, Cecilia Moloney, Salama Ikki +// DOI: 10.1109/MSP.2019.2927685 +// arXiv:1907.07241v1 +// +// GUO’S algorithm without iterative procedure +// Solution of AX = B -> {A = LU} -> LUX = B -> {UX = Y} -> LY = B -> Y +// UX = Y -> X +// | Sum(y^2) Sum(x*y^2) Sum(x^2*y^2) | +// MATRIX A = | Sum(x*y^2) Sum(x^2*y^2) Sum(x^3*y^2) | +// | Sum(x^2*y^2) Sum(x^3*y^2) Sum(x^4*y^2) | +// +// | Sum(y^2*Log(y)) | +// MATRIX B = | Sum(x*y^2*Log(y)) | +// | Sum(x^2*y^2*Log(y)) | +// +// | l11 0 0 | | l11 l21 l31 | +// L = | l21 l22 0 | U = | 0 l22 l32 | +// | l31 l32 l33 | | 0 0 l33 | +// +const double dSQRT_2_PI = sqrt(M_PI_2); +inline double calcGaussArea(double x, double y, double d, double xMin, double xMax) +{ // Gauss area from xMin to xMax + double d2 = M_SQRT2 * d; // deviation + double fErfMin = std::erf((xMin - x) / d2); + double fErfMax = std::erf((xMax - x) / d2); + return dSQRT_2_PI * y * d * (fErfMax - fErfMin); +}; + +template <typename Y> +double calcMatrix(const std::vector<Y> &rvy, std::array<double, 5> &radSum, std::array<double, 3> &radSumLog) +{ + double dSum0 = 0, dSum1 = 0, dSum2 = 0, dSum3 = 0, dSum4 = 0; + double dSumLog0 = 0, dSumLog1 = 0, dSumLog2 = 0; + uint iMax = 0; + Y yMax = 0; + double dGravityAxis = 0, dGravitySum = 0; + size_t nSize = rvy.size(); + assert(nSize < 256); + + for (uint i = 0; i < nSize; i++) { + uint i2 = i * i; + uint i3 = i2 * i; + uint i4 = i3 * i; + + Y y = rvy[i]; + if (y > yMax) { + yMax = y; + iMax = i; + } + dGravitySum += y; + dGravityAxis += i * y; + double d2 = y * y; + + dSum0 += d2; + dSum1 += i * d2; + dSum2 += i2 * d2; + dSum3 += i3 * d2; + dSum4 += i4 * d2; + + double dYLog = d2 * log(y); + dSumLog0 += dYLog; + dSumLog1 += i * dYLog; + dSumLog2 += i2 * dYLog; + } + radSum = {dSum0, dSum1, dSum2, dSum3, dSum4}; + radSumLog = {dSumLog0, dSumLog1, dSumLog2}; + dGravityAxis = dGravitySum > 0.0 ? dGravityAxis / dGravitySum : 0.0; + return dGravityAxis > 0 ? dGravityAxis : iMax; +} +template <typename Y> +std::array<double, 4> calcGaussRes(const std::vector<Y> &rvy) +{ + assert(rvy.size() > 2); + + std::array<double, 5> adSum; + std::array<double, 3> adSumLog; + double xMax = calcMatrix(rvy, adSum, adSumLog); + + double dL11 = sqrt(adSum[0]); + double dL21 = adSum[1] / dL11; + double dL22 = sqrt(adSum[2] - dL21 * dL21); + double dL31 = adSum[2] / dL11; + double dL32 = (adSum[3] - dL21 * dL31) / dL22; + double dL33sqrt = adSum[4] - dL31 * dL31 - dL32 * dL32; // should be > 0.0 (double precision error!) + // example: 2969669.53125 - 2969669.5312500002727183440017704 < 0.0! + double dL33 = dL33sqrt > 0.0 ? sqrt(dL33sqrt) : 0.0; + + double dY1 = adSumLog[0] / dL11; + double dY2 = (adSumLog[1] - dL21 * dY1) / dL22; + double dY3 = dL33 > 0.0 ? (adSumLog[2] - dL31 * dY1 - dL32 * dY2) / dL33 : 0.0; + + double dX2 = dL33 > 0.0 ? dY3 / dL33 : 0.0; + double dX1 = (dY2 - dL32 * dX2) / dL22; + double dX0 = (dY1 - dL31 * dX2 - dL21 * dX1) / dL11; + if (dX2 != 0.0) { // some times when adc has const values we should define 'x' in the middle of array size + double dX25 = -0.5 / dX2; + double x = dX1 * dX25; + double y = exp(dX0 + 0.5 * dX1 * x); + double d = dX25 > 0 ? sqrt(dX25) : 0; // deviation + return {x, y, d, xMax}; + } else + return {rvy.size() / 2.0, xMax, 0.0, xMax}; // in the middle of array +} +template <typename Y> +double calcRsquared(const std::vector<Y> &rvy, double x, double y, double d, double dYSum) +{ // coefficient of determination + size_t nSize = rvy.size(); + if (nSize < 3 || d == 0.0) return 0; + // deviation should be defined!!! + double dYMean = dYSum / nSize; + double d2Deviation2 = 2 * d * d; + double dSumTotal = 0; // Sum total = Sum[(y - ymean)^2] + double dSumRes = 0; // Sum of residuals = Sum[(y - f(x))^2] + for (uint i = 0; i < nSize; i++) { + double dX = rvy[i] - x; + double dGauss = y * exp(-dX * dX / d2Deviation2); + double dY = rvy[i]; + double dTot = dY - dYMean; + double dRes = dY - dGauss; + + dSumTotal += dTot * dTot; + dSumRes += dRes * dRes; + } + return 1 - dSumRes / dSumTotal; +}; +template <typename X, typename Y> +struct XY { + X x; + Y y; +}; +template <typename X, typename Y> +double calcMatrix(const std::vector<XY<X, Y>> &rvXY, std::array<double, 5> &radSum, std::array<double, 3> &radSumLog) +{ + double dSum0 = 0, dSum1 = 0, dSum2 = 0, dSum3 = 0, dSum4 = 0; + double dSumLog0 = 0, dSumLog1 = 0, dSumLog2 = 0; + double yMax = 0; + X xMax = 0, xMin = 0; + double dGravityAxis = 0, dGravitySum = 0; + size_t nSize = rvXY.size(); + for (uint i = 0; i < nSize; i++) { + X x = rvXY[i].x; + X x2 = x * x; + X x3 = x2 * x; + X x4 = x3 * x; + if (!i || x < xMin) xMin = x; + + double y = rvXY[i].y; + if (y > yMax) { + yMax = y; + xMax = x; + } + dGravitySum += y; + dGravityAxis += x * y; + double y2 = y * y; + + dSum0 += y2; + dSum1 += x * y2; + dSum2 += x2 * y2; + dSum3 += x3 * y2; + dSum4 += x4 * y2; + + double dYLog = y2 * log(y); + dSumLog0 += dYLog; + dSumLog1 += x * dYLog; + dSumLog2 += x2 * dYLog; + } + radSum = {dSum0, dSum1, dSum2, dSum3, dSum4}; + radSumLog = {dSumLog0, dSumLog1, dSumLog2}; + dGravityAxis = dGravitySum > 0.0 ? dGravityAxis / dGravitySum : 0.0; + return dGravityAxis <= xMin ? xMax : dGravityAxis; +} +template <typename X, typename Y> +std::array<double, 4> calcGaussRes(const std::vector<XY<X, Y>> &rvXY) +{ + assert(rvXY.size() > 2); + + std::array<double, 5> adSum; + std::array<double, 3> adSumLog; + double xMax = calcMatrix<X, Y>(rvXY, adSum, adSumLog); + + double dL11 = sqrt(adSum[0]); + double dL21 = adSum[1] / dL11; + double dL22 = sqrt(adSum[2] - dL21 * dL21); + double dL31 = adSum[2] / dL11; + double dL32 = (adSum[3] - dL21 * dL31) / dL22; + double dL33 = sqrt(adSum[4] - dL31 * dL31 - dL32 * dL32); + + double dY1 = adSumLog[0] / dL11; + double dY2 = (adSumLog[1] - dL21 * dY1) / dL22; + double dY3 = (adSumLog[2] - dL31 * dY1 - dL32 * dY2) / dL33; + + double dX2 = dY3 / dL33; + double dX1 = (dY2 - dL32 * dX2) / dL22; + double dX0 = (dY1 - dL31 * dX2 - dL21 * dX1) / dL11; + + double dX25 = -0.5 / dX2; + double x = dX1 * dX25; + double y = exp(dX0 + 0.5 * dX1 * x); + double d = dX25 > 0 ? sqrt(dX25) : 0; // deviation + return {x, y, d, xMax}; +} +template <typename X, typename Y> +double calcRsquared(const std::vector<XY<X, Y>> &rvXY, double x, double y, double d, double dYSum) +{ // coefficient of determination + size_t nSize = rvXY.size(); + if (nSize > 2) { // deviation should be defined!!! + double dYMean = dYSum / nSize; + double d2Deviation2 = 2 * d * d; + double dSumTotal = 0; // Sum total = Sum[(y - ymean)^2] + double dSumRes = 0; // Sum of residuals = Sum[(y - f(x))^2] + for (uint i = 0; i < nSize; i++) { + double dX = rvXY[i].x - x; + double dGauss = y * exp(-dX * dX / d2Deviation2); + double dY = rvXY[i].y; + double dTot = dY - dYMean; + double dRes = dY - dGauss; + + dSumTotal += dTot * dTot; + dSumRes += dRes * dRes; + } + return 1 - dSumRes / dSumTotal; + } else + return 0; +}; +class TpcAdcGauss { + std::array<double, 4> + _arr; // array[0: mean value, 1: gaussian peak value (ADC), 2: standard deviation, 3: adjustment] + double _dArea; // gaussian area of the normal distribution (should be defined for reference only!) + double _dRsquared; // r-squared value of the normal distribution (should be defined for reference only!) +public: + TpcAdcGauss() : _arr{0, 0, 0, 0}, _dArea(0), _dRsquared(0) {} + TpcAdcGauss(const std::array<double, 4> &r) : _arr(r), _dArea(0), _dRsquared(0) {} + TpcAdcGauss(const TpcAdcGauss &r) + : _arr(r._arr), _dArea(r._dArea), _dRsquared(r._dRsquared) { /*_array = r._array;*/ } + TpcAdcGauss &operator=(const TpcAdcGauss &r) + { + if (this != &r) { + _arr = r._arr; + _dArea = r._dArea; + _dRsquared = r._dRsquared; + } + return *this; + } + void Set(const std::array<double, 4> &r) { _arr = r; } + + void setMean(double d) { _arr[0] = d; } + double getMean() const { return _arr[0]; } + + void setMax(double d) { _arr[1] = d; } + double getMax() const { return _arr[1]; } + + void setDeviation(double d) { _arr[2] = d; } + double getDeviation() const { return _arr[2]; } + + double getAdjustment() const { return _arr[3]; } // adjustment value for unreasonable fit results to max ADC value!!! + void calc(const std::vector<float> &rv) + { + size_t nSize = rv.size(); + if (nSize == 1) + _arr = {0, double(rv.front()), 0, 0}; + else if (nSize == 2) { + float f1st = rv.front(); + float f2nd = rv.back(); + double dMax = f2nd > f1st ? f2nd : f1st; + double dMean = f2nd / (f1st + f2nd); + _arr = {dMean, dMax, 0, 0}; + } else { + _arr = calcGaussRes<float>(rv); +#ifdef AUXHIT_DEBUG + if (_AdcGauss.getDeviation() != 0) { + _AdcGauss.setArea(::calcGaussArea(gauss.getMean(), gauss.getMax(), gauss.getDeviation(), _nTimeBin, nSize)); + _AdcGauss.setRsquared( + calcRsquared<X, Y>(_vXY, gauss.getMean(), gauss.getMax(), gauss.getDeviation(), double(_ySum))); + } else { // it cannot fit the Gauss + _AdcGauss.setArea(0); + _AdcGauss.setRsquared(0); + } +#endif + } + } + void setArea(double d) { _dArea = d; } + double getArea() const { return _dArea; } + void calcArea(uint nSize) { _dArea = calcGaussArea(_arr[0], _arr[1], _arr[2], 0, nSize - 1); } + + void setRsquared(double d) { _dRsquared = d; } + double getRsquared() const { return _dRsquared; } + // void calcRsquared() { + // _dRsquared = ::calcRsquared<float>( _vXY, getMean(), getMax(), getDeviation(), double(_ySum) ) ); + // } + std::ostringstream &Print(std::ostringstream &oss) const + { + oss << " {mean:" << getMean() << " max:" << getMax() << " deviation:" << getDeviation() + << " adjustment:" << getAdjustment() << " area:" << _dArea << " R-squared:" << _dRsquared << "}"; + return oss; + } +}; + +} // namespace TpcClustering + +#endif // TPC_ADC_GAUSS_H \ No newline at end of file diff --git a/detectors/tpc/clusterHitFinder/fast/clustering/TpcAdcHit.h b/detectors/tpc/clusterHitFinder/fast/clustering/TpcAdcHit.h new file mode 100644 index 00000000..bcff1c1b --- /dev/null +++ b/detectors/tpc/clusterHitFinder/fast/clustering/TpcAdcHit.h @@ -0,0 +1,77 @@ +// +// Joint Institute for Nuclear Research (JINR), Dubna, Russia, Jun-2024 +// TPC ADC hit structure for MPD/NICA project +// author: Viktor A.Krylov JINR/LNP +// email: kryman@jinr.ru +// +#ifndef ADC_HIT_H +#define ADC_HIT_H + +#include <iostream> +#include <climits> + +#include "TpcSectorGeo.h" +#include "TpcTrack.h" + +namespace TpcClustering { + +struct TpcAdcHit { + float _fX; // global coords:{ x, y, z } + float _fY; + float _fZ; + uint _nRow; // sector row number + uint _nClusterId; // cluster id + float _fAdc; // cluster adc value + + TpcAdcHit() : _fX(0), _fY(0), _fZ(0), _fAdc(0), _nRow(0), _nClusterId(0) {} + TpcAdcHit(uint nRow, uint nClusterId, float fAdc) + : _fX(0), _fY(0), _fZ(0), _nRow(nRow), _nClusterId(nClusterId), _fAdc(fAdc) + { + } + TpcAdcHit(float fX, float fY, float fZ, uint nRow, uint nClusterId, float fAdc) + : _fX(fX), _fY(fY), _fZ(fZ), _nRow(nRow), _nClusterId(nClusterId), _fAdc(fAdc) + { + } + TpcAdcHit &operator=(const TpcAdcHit &r) + { + return this == &r ? *this : *(TpcAdcHit *)memcpy(this, (const void *)&r, sizeof(TpcAdcHit)); + } + void operator=(const TpcPoint &r) + { + _fX = r._fX; + _fY = r._fY; + _fZ = r._fZ; + } + void Set(uint nRow, uint nClusterId, float fAdc) + { + _nRow = nRow; + _nClusterId = nClusterId; + _fAdc = fAdc; + } + void Set(float fX, float fY, float fZ) + { + _fX = fX; + _fY = fY; + _fZ = fZ; + } + void setGlobal(uint nSector, uint nRow, float fPad, float fTimeBin) + { + this->operator=(TpcPoint(nRow, fPad, fTimeBin).getGlobal(nSector)); + } + void shiftSampaDelay(uint nSector, float fDelay) { _fZ += nSector < 12 ? fDelay : -fDelay; } + friend std::ostream &operator<<(std::ostream &os, const TpcAdcHit &r) + { + os << "AdcHit:: x:" << r._fX << "; y:" << r._fY << "; z:" << r._fZ << "; row:" << r._nRow + << "; cluster:" << r._nClusterId << "; adc:" << r._fAdc << std::endl; + return os; + } + void Print() const + { + std::cout << "AdcHit:: x:" << _fX << "; y:" << _fY << "; z:" << _fZ << "; row:" << _nRow + << "; cluster:" << _nClusterId << "; adc:" << _fAdc << std::endl; + } +}; + +} // namespace TpcClustering + +#endif // ADC_HIT_H diff --git a/detectors/tpc/clusterHitFinder/fast/clustering/TpcAdcStat.h b/detectors/tpc/clusterHitFinder/fast/clustering/TpcAdcStat.h new file mode 100644 index 00000000..133e8ff5 --- /dev/null +++ b/detectors/tpc/clusterHitFinder/fast/clustering/TpcAdcStat.h @@ -0,0 +1,52 @@ +#ifndef TPC_ADC_STAT_H +#define TPC_ADC_STAT_H + +#include <iostream> +#include <climits> +#include <float.h> + +namespace TpcClustering { + +struct TpcAdcStat { + uint _nPadMin; + uint _nPadMax; + uint _nTimeBinMin; + uint _nTimeBinMax; + float _fAdcMin; + float _fAdcMax; + TpcAdcStat() + : _nPadMin(UINT_MAX), _nPadMax(0), _nTimeBinMin(UINT_MAX), _nTimeBinMax(0), _fAdcMin(FLT_MAX), _fAdcMax(0) + { + } + TpcAdcStat(uint nPad, uint nTimeBin, float fAdc) + { + _nPadMin = _nPadMax = nPad; + _nTimeBinMin = _nTimeBinMax = nTimeBin; + _fAdcMin = _fAdcMax = fAdc; + } + void Reset() + { + _nPadMin = _nTimeBinMin = UINT_MAX; + _nPadMax = _nTimeBinMax = 0; + _fAdcMin = FLT_MAX; + _fAdcMax = 0; + } + void PickOut(uint nPad, uint nTimeBin, float fAdc) + { + if (nPad > _nPadMax) _nPadMax = nPad; + if (nPad < _nPadMin) _nPadMin = nPad; + if (nTimeBin > _nTimeBinMax) _nTimeBinMax = nTimeBin; + if (nTimeBin < _nTimeBinMin) _nTimeBinMin = nTimeBin; + if (fAdc > _fAdcMax) _fAdcMax = fAdc; + if (fAdc < _fAdcMin) _fAdcMin = fAdc; + } + void Print(std::ostringstream &oss) const + { + oss << "TpcAdcStat:{ pad:[" << _nPadMin << "-" << _nPadMax << "], timebin:[" << _nTimeBinMin << "-" + << _nTimeBinMax << "], adc:[" << _fAdcMin << "-" << _fAdcMax << "] }"; + } +}; + +} // namespace TpcClustering + +#endif // TPC_ADC_STAT_H \ No newline at end of file diff --git a/detectors/tpc/clusterHitFinder/fast/clustering/TpcCluster.h b/detectors/tpc/clusterHitFinder/fast/clustering/TpcCluster.h new file mode 100644 index 00000000..1e22099a --- /dev/null +++ b/detectors/tpc/clusterHitFinder/fast/clustering/TpcCluster.h @@ -0,0 +1,313 @@ +// +// Joint Institute for Nuclear Research (JINR), Dubna, Russia, Jun-2024 +// TPCCluster object for MPD/NICA project +// author: Viktor A.Krylov JINR/LNP +// email: kryman@jinr.ru +// +#ifndef TPC_CLUSTER_H_ +#define TPC_CLUSTER_H_ + +#include <list> +#include <algorithm> +#include <iterator> +#include <iomanip> + +#include "TpcHit.h" +#include "TpcPadCluster.h" +#include "TpcSectorGeo.h" +// #include "../TpcEventObj.h" + +#define CLUSTER_PAD_MAX 5 + +namespace TpcClustering { + +class TpcCluster : public TpcHit { + friend class PadCluster; + uint _nTimeBinMin; // PadClusters timebin min range + uint _nTimeBinMax; // PadClusters timebin max range + uint _nPadMin; // PadClusters pad min range + uint _nPadMax; // PadClusters pad max range + std::list<const TpcPadCluster *> _lstTpcPadClusters; + +public: + std::vector<uint> _vnTpcDigits; + TpcCluster() : TpcHit(), _nTimeBinMin(UINT_MAX), _nTimeBinMax(UINT_MAX), _nPadMin(UINT_MAX), _nPadMax(UINT_MAX) {} + TpcCluster(uint nTimeBinMin, uint nTimeBinMax, uint nPadMin, uint nPadMax) : TpcHit() + { + _nTimeBinMin = nTimeBinMin; + _nTimeBinMax = nTimeBinMax; + _nPadMin = nPadMin; + _nPadMax = nPadMax; + } + TpcCluster(uint nId, const TpcPadCluster *pTpcPadCluster) : TpcHit(nId, pTpcPadCluster->getAdcSum()) + { + _nTimeBinMin = pTpcPadCluster->getTimeBinMin(); + _nTimeBinMax = pTpcPadCluster->getTimeBinMax(); + _nPadMin = _nPadMax = pTpcPadCluster->getPad(); + _lstTpcPadClusters.push_back(pTpcPadCluster); + ((TpcPadCluster *)pTpcPadCluster)->setCluster(this); + if (pTpcPadCluster->isUndefined()) TpcHit::setUndefined(); // the tpc hit is undefined + } + virtual ~TpcCluster() + { + for (const TpcPadCluster *pTpcPadCluster : _lstTpcPadClusters) delete pTpcPadCluster; + _lstTpcPadClusters.clear(); + _vnTpcDigits.clear(); + } + const std::list<const TpcPadCluster *> &getPadClusters() const { return _lstTpcPadClusters; } + inline uint getTimeBinMin() const { return _nTimeBinMin; } + inline uint getTimeBinMax() const { return _nTimeBinMax; } + inline uint getPadMin() const { return _nPadMin; } + inline uint getPadMax() const { return _nPadMax; } + bool isOut(uint nPad) const { return _nPadMin && nPad < _nPadMin - 1 || nPad > _nPadMax + 1; } + + void defDigits(const std::vector<TpcDigit> *pvTpcDigits, uint nRow) + { + for (const TpcPadCluster *pTpcPadCluster : _lstTpcPadClusters) { + uint nPad = pTpcPadCluster->getPad(); + uint nTimeBinMin = pTpcPadCluster->getTimeBinMin(); + uint nTimeBinMax = pTpcPadCluster->getTimeBinMax(); + std::vector<TpcDigit>::const_iterator it = pvTpcDigits->begin(); + while ((it = std::find_if(it, pvTpcDigits->end(), [&](const auto &r) { + return r._nRow == nRow && r._nPad == nPad && r._nTimeBin >= nTimeBinMin && + r._nTimeBin <= nTimeBinMax; + })) != pvTpcDigits->end()) + _vnTpcDigits.push_back(it++ - pvTpcDigits->begin()); + } + } + void setBinContent(TH2S &rH2S) const + { +#ifdef TPC_CLUSTER_H2S + for (const TpcPadCluster *pTpcPadCluster : _lstTpcPadClusters) pTpcPadCluster->setBinContent(rH2S); +#endif + } + void setBinContent(TH2F &rH2F) const + { +#ifdef TPC_CLUSTER_H2F + for (const TpcPadCluster *pTpcPadCluster : _lstTpcPadClusters) pTpcPadCluster->setBinContent(rH2F); +#endif + } + void resetStat() { _nTimeBinMin = _nTimeBinMax = _nPadMin = _nPadMax = UINT_MAX; } + void calcStat() + { + resetStat(); + float fAdcSum = 0; + for (const auto &pTpcPadCluster : _lstTpcPadClusters) { + calcStat(pTpcPadCluster); + ((TpcPadCluster *)pTpcPadCluster)->setCluster(this); + fAdcSum += pTpcPadCluster->getAdcSum(); + } + TpcHit::setAdcSum(fAdcSum); + } + void calcStat(const TpcPadCluster *pTpcPadCluster) + { + if (_nTimeBinMin == UINT_MAX && _nTimeBinMax == UINT_MAX) { + _nTimeBinMin = pTpcPadCluster->getTimeBinMin(); + _nTimeBinMax = pTpcPadCluster->getTimeBinMax(); + } else { + uint nTimeBinMin = pTpcPadCluster->getTimeBinMin(); + if (nTimeBinMin < _nTimeBinMin) _nTimeBinMin = nTimeBinMin; + uint nTimeBinMax = pTpcPadCluster->getTimeBinMax(); + if (nTimeBinMax > _nTimeBinMax) _nTimeBinMax = nTimeBinMax; + } + if (_nPadMin == UINT_MAX && _nPadMax == UINT_MAX) + _nPadMin = _nPadMax = pTpcPadCluster->getPad(); + else { + uint nPad = pTpcPadCluster->getPad(); + if (nPad < _nPadMin) + _nPadMin = nPad; + else if (nPad > _nPadMax) + _nPadMax = nPad; + } + } + void Calc(const std::vector<float> &rvfAdc, uint nRow, const TpcPadCluster *pTpcPadCluster) + { + assert(rvfAdc.size()); + + if (pTpcPadCluster->isUndefined()) TpcHit::setUndefined(); + + const TpcPadCluster *pTpcPadClusterFront = _lstTpcPadClusters.front(); + assert(pTpcPadClusterFront); + + TpcHit::Calc(rvfAdc, nRow, pTpcPadClusterFront->getPad(), pTpcPadClusterFront->getTimeBin(), + pTpcPadCluster->getTimeBin()); + } + void Add(const TpcPadCluster *pTpcPadCluster) + { + calcStat(pTpcPadCluster); + TpcHit::addAdcSum(pTpcPadCluster->getAdcSum()); + _lstTpcPadClusters.push_back(pTpcPadCluster); + ((TpcPadCluster *)pTpcPadCluster)->setCluster(this); + } + // inline void Clear() { // should be used with attention to except memory leaks! + // setAdcSum(0); + // _lstPadClusters.clear(); + // } + // Rotate coordinate axes + // http://www.pm298.ru/dekart.php + // x = (X - a) * cos + (Y - b) * sin + // y = (Y - b) * cos - (X - a) * sin + // X = x * cos - y * sin + a + // Y = y * cos + x * sin + b + + // x = X * cos - Y * sin + // y = X * sin + Y * cos + TpcCluster *Split(uint nRow) + { // split current Cluster on new ones depends on AdcMax values + new array of ADC max values + assert(_lstTpcPadClusters.size()); + + if (isUndefined()) // the tpc cluster is well undefined + return NULL; + + std::vector<float> vfAdc; // pad plane {pad:adc} + uint iMin; // index of the last adc min value in the vector; + std::list<const TpcPadCluster *>::const_iterator iAdcMin; + for (std::list<const TpcPadCluster *>::iterator i = _lstTpcPadClusters.begin(); i != _lstTpcPadClusters.end(); + i++) { + const TpcPadCluster *pTpcPadCluster = *i; + assert(!pTpcPadCluster->isUndefined()); + + float fTimeBin = pTpcPadCluster->getTimeBin(); + float fAdc = pTpcPadCluster->getAdcMax(); + assert(fAdc > 0.0f && fAdc < 65535.0f && !isnan(fAdc) && !isnan(fTimeBin)); + + if (vfAdc.size() == 0) { // pad cluster adc max vector seed + iMin = 0; + iAdcMin = i; // iterator of minimal pad cluster + vfAdc.push_back(fAdc); + continue; + } else if (fAdc < vfAdc.back()) { + iMin = vfAdc.size(); + iAdcMin = i; // iterator of minimal pad cluster + vfAdc.push_back(fAdc); + continue; + } else if (fAdc > vfAdc.back() && iMin != 0) { // new tpc cluster + assert(vfAdc.size() > 1); + + uint nMin = vfAdc.size() - iMin; + if (nMin == 1) { // simple minimum /\./ + const TpcPadCluster *pTpcPadClusterMin = *iAdcMin; + TpcPadCluster *pTpcPadClusterNew = new TpcPadCluster(*pTpcPadClusterMin); + float fAdcLeft = vfAdc[iMin - 1]; + float fScale = fAdcLeft / (fAdcLeft + fAdc); // ^.^ + vfAdc.back() *= fScale; // + pTpcPadClusterNew->Split(fScale); + TpcCluster *pTpcClusterNew = new TpcCluster(); + pTpcClusterNew->_lstTpcPadClusters.splice(pTpcClusterNew->_lstTpcPadClusters.begin(), _lstTpcPadClusters, + _lstTpcPadClusters.begin(), iAdcMin); + pTpcClusterNew->_lstTpcPadClusters.push_back(pTpcPadClusterNew); + pTpcClusterNew->calcStat(); + pTpcClusterNew->Calc(vfAdc, nRow, pTpcPadClusterNew); + ((TpcPadCluster *)pTpcPadClusterMin)->Split(1.0f - fScale); + calcStat(); + return pTpcClusterNew; + } else { // undefined minimum /\.../ + uint nMin2 = nMin >> 1; + vfAdc.resize(vfAdc.size() - nMin2); + const TpcPadCluster *pTpcPadClusterLast = nullptr; + for (uint i = 0; i < nMin2; i++) pTpcPadClusterLast = *(iAdcMin++); + assert(pTpcPadClusterLast); + + TpcCluster *pTpcClusterNew = new TpcCluster(); + pTpcClusterNew->_lstTpcPadClusters.splice(pTpcClusterNew->_lstTpcPadClusters.begin(), _lstTpcPadClusters, + _lstTpcPadClusters.begin(), iAdcMin); + if (nMin & 0x1) { // odd: 3 and more + const TpcPadCluster *pTpcPadClusterMin = *iAdcMin; + ((TpcPadCluster *)pTpcPadClusterMin)->Split(0.5f); // devide all adc values between two clusters + TpcPadCluster *pTpcPadClusterNew = new TpcPadCluster(*pTpcPadClusterMin); + pTpcClusterNew->_lstTpcPadClusters.push_back(pTpcPadClusterNew); + pTpcClusterNew->calcStat(); + vfAdc.back() *= 0.5f; // devide max adc value between two clusters + pTpcClusterNew->Calc(vfAdc, nRow, pTpcPadClusterNew); + } else { // even: 2 and more + pTpcClusterNew->calcStat(); + pTpcClusterNew->Calc(vfAdc, nRow, pTpcPadClusterLast); + } + calcStat(); + return pTpcClusterNew; + } + } else + vfAdc.push_back(fAdc); + } + assert(!vfAdc.empty()); + + Calc(vfAdc, nRow, *iAdcMin); + return NULL; + // std::ranges::sort( vxy, [](XY<float, float> xy1st, XY<float, float> xy2nd) {return xy1st.x < xy2nd.x;} ); + // std::sort(vxy.begin(), vxy.end(), [](XY<float, float> xy1st, XY<float, float> xy2nd) {return xy1st.x < + // xy2nd.x;} ); + } + std::ostringstream &Print(std::ostringstream &oss, const std::list<const TpcPadCluster *> &r) const + { + for (std::list<const TpcPadCluster *>::const_iterator i = r.begin(); i != r.end(); i++) { + oss << std::string(4, ' '); + (*i)->Print(oss); + } + return oss; + } + std::ostringstream &Print(std::ostringstream &oss) const + { + TpcHit::Print(oss); + const uint nTAB = 2; + oss << std::string(nTAB, ' ') << "Cluster region: " + << "timebin:{ min:" << _nTimeBinMin << " max:" << _nTimeBinMax << " }, pad:{ min:" << _nPadMin + << " max:" << _nPadMax << " }" << std::endl; + Print(oss, _lstTpcPadClusters); + return oss; + } + void Draw(std::ostringstream &oss) const + { + for (std::list<const TpcPadCluster *>::const_iterator i = _lstTpcPadClusters.begin(); + i != _lstTpcPadClusters.end(); i++) { + const TpcPadCluster *pTpcPadCluster = (*i); + oss << std::setw(3) << pTpcPadCluster->getPad() << ":"; + uint nPos = 0; + uint nMin = pTpcPadCluster->getTimeBinMin(); + if (nMin < nPos) { // shift back to the 'n' positions (for combined PadClusters!!!) + long int nShift = nPos - nMin; + oss.seekp(oss.tellp() - nShift); + nPos -= nShift; + } + uint nGap = nMin - nPos; + if (nGap) { + oss << std::string(nGap, '-'); + nPos += nGap; + } + char ch = SZ_ALPHABET[getId() % strlen(SZ_ALPHABET)]; + // pPadCluster->Draw(oss, this, ch); + pTpcPadCluster->Draw(oss, ch); + oss << std::endl; + nPos += pTpcPadCluster->getSize(); + } + } + void Draw(std::ostringstream &oss, bool bAdc = false) const + { + for (std::list<const TpcPadCluster *>::const_iterator i = _lstTpcPadClusters.begin(); + i != _lstTpcPadClusters.end(); i++) { + const TpcPadCluster *pTpcPadCluster = (*i); + uint nMin = pTpcPadCluster->getTimeBinMin(); + if (bAdc) pTpcPadCluster->Draw(oss, nMin + 4); + oss << std::setw(3) << pTpcPadCluster->getPad() << ":"; + if (nMin) oss << std::string(nMin, '-'); + char ch = SZ_ALPHABET[getId() % strlen(SZ_ALPHABET)]; + pTpcPadCluster->Draw(oss, ch); + oss << std::endl; + } + } +}; +uint TpcPadCluster::getClusterId(const TpcCluster *pTpcCluster) const +{ + return pTpcCluster == NULL ? 0 : ((const TpcHit *)pTpcCluster)->getId(); +} +float TpcPadCluster::getClusterTime(const TpcCluster *pTpcCluster) const +{ + return pTpcCluster == NULL ? 0 : ((const TpcHit *)pTpcCluster)->getTimeBin(); +} +float TpcPadCluster::getClusterPad(const TpcCluster *pTpcCluster) const +{ + return pTpcCluster == NULL ? 0 : ((const TpcHit *)pTpcCluster)->getPad(); +} + +} // namespace TpcClustering + +#endif // TPC_CLUSTER_H \ No newline at end of file diff --git a/detectors/tpc/clusterHitFinder/fast/clustering/TpcClusters.cxx b/detectors/tpc/clusterHitFinder/fast/clustering/TpcClusters.cxx new file mode 100644 index 00000000..3ba9182d --- /dev/null +++ b/detectors/tpc/clusterHitFinder/fast/clustering/TpcClusters.cxx @@ -0,0 +1,16 @@ +// +// Joint Institute for Nuclear Research (JINR), Dubna, Russia, Jun-2024 +// TPCClusters object for MPD/NICA project +// author: Viktor A.Krylov JINR/LNP +// email: kryman@jinr.ru +// +#include "TpcClusters.h" + +namespace TpcClustering { + +TpcClusters::TpcClusters(int nSector, uint nRow, uint nPad, uint nTimeBin, float fAdc) : _nSector(nSector) +{ + _lstTpcRowClusters.push_back(new TpcRowClusters(nRow, nPad, nTimeBin, fAdc)); +} + +} // namespace TpcClustering \ No newline at end of file diff --git a/detectors/tpc/clusterHitFinder/fast/clustering/TpcClusters.h b/detectors/tpc/clusterHitFinder/fast/clustering/TpcClusters.h new file mode 100644 index 00000000..ba7d3293 --- /dev/null +++ b/detectors/tpc/clusterHitFinder/fast/clustering/TpcClusters.h @@ -0,0 +1,190 @@ +// +// Joint Institute for Nuclear Research (JINR), Dubna, Russia, Jun-2024 +// TPCClusters object for MPD/NICA project +// author: Viktor A.Krylov JINR/LNP +// email: kryman@jinr.ru +// +#ifndef TPC_CLUSTERS_H +#define TPC_CLUSTERS_H + +#include <algorithm> + +#include <TDirectory.h> +#include <TH2I.h> +#include <TH3S.h> + +#include "TpcSectorGeo.h" +#include "TpcAdcHit.h" +#include "TpcRowClusters.h" +#include "TpcSampaPack.h" + +#define DEBUG_OUTPUT + +#define szTPC_CLUSTERS "clusters" +#define szTPC_HITS "hits" + +namespace TpcClustering { + +class TpcClusters { + int _nSector; + +public: + std::list<TpcRowClusters *> _lstTpcRowClusters; + + TpcClusters() : _nSector(0) {} + TpcClusters(int nSector) : _nSector(nSector) {} + TpcClusters(int nSector, uint nRow, uint nPad, uint nTimeBin, float fAdc); + // { + // _lstTpcRowClusters.push_back( new TpcRowClusters( nRow, nPad, nTimeBin, fAdc ) ); + // } + virtual ~TpcClusters() + { + for (const TpcRowClusters *pTpcRowClusters : _lstTpcRowClusters) delete pTpcRowClusters; + _lstTpcRowClusters.clear(); + } + inline int getSector() const { return _nSector; } + inline void setSector(int nSector) { _nSector = nSector; } + inline bool isEmpty() const { return _lstTpcRowClusters.empty(); } + inline std::size_t getSize() const { return _lstTpcRowClusters.size(); } + const TpcCluster *findTpcCluster(uint nRow, uint nClusterId) const + { + std::list<TpcRowClusters *>::const_iterator iTpcRowClusters = + std::find_if(_lstTpcRowClusters.begin(), _lstTpcRowClusters.end(), + [nRow](const TpcRowClusters *p) { return p->getRow() == nRow; }); + if (iTpcRowClusters != _lstTpcRowClusters.end()) { + const TpcRowClusters *pTpcRowClusters = *iTpcRowClusters; + std::list<TpcCluster *>::const_iterator iTpcCluster = + std::find_if(pTpcRowClusters->getClusters().begin(), pTpcRowClusters->getClusters().end(), + [nClusterId](const TpcCluster *p) { return p->getId() == nClusterId; }); + if (iTpcCluster != pTpcRowClusters->getClusters().end()) return *iTpcCluster; + } + return NULL; + } + void defDigits(const std::vector<TpcDigit> *pvTpcDigits) + { + for (std::list<TpcRowClusters *>::const_iterator i = _lstTpcRowClusters.begin(); i != _lstTpcRowClusters.end(); + i++) + (*i)->defDigits(pvTpcDigits); + } + void Store(TDirectory *pDir) + { + TDirectory *pClustersDir = pDir->mkdir(szTPC_CLUSTERS, "Tpc Sector Clusters", kTRUE); + if (pClustersDir && pClustersDir->cd()) + for (TpcRowClusters *pTpcRowClusters : _lstTpcRowClusters) pTpcRowClusters->Store(pClustersDir); + } + void defClusters() + { // for release + for (TpcRowClusters *pTpcRowClusters : _lstTpcRowClusters) { + pTpcRowClusters->calcPadClusters(); // define pad cluster hits in the row + pTpcRowClusters->Joint(); + pTpcRowClusters->Split(); + } + } + void defClusters(std::ostringstream &oss) + { // for debug only! + for (std::list<TpcRowClusters *>::const_iterator i = _lstTpcRowClusters.begin(); i != _lstTpcRowClusters.end(); + i++) { + TpcRowClusters *pTpcRowClusters = *i; + pTpcRowClusters->calcPadClusters(); // define pad cluster hit in the row +#ifdef DEBUG_OUTPUT + oss << "Before joint" << std::endl; + // DrawPadClusters(oss); + pTpcRowClusters->DrawPadClusters(oss); + // PrintPadClusters(oss); + // pRowClusters->PrintPadClusterMatrix(oss); + // pRowClusters->Draw(oss); + // DrawRowClusters(oss); + // PrintPadClusters(oss); + std::cout << oss.str(); + oss.str(""); // clear output string buffer + oss.clear(); +#endif + pTpcRowClusters->Joint(); +#ifdef DEBUG_OUTPUT + oss << "Before split" << std::endl; + // pRowClusters->DrawPadClusters(oss); + // pRowClusters->DrawRowTpcClusters(oss); + pTpcRowClusters->Draw(oss); + // Draw(oss); + // Print(oss); + // return; + std::cout << oss.str(); + oss.str(""); // clear string buffer + oss.clear(); +#endif + pTpcRowClusters->Split(); +#ifdef DEBUG_OUTPUT + oss << "After split" << std::endl; + // pRowClusters->DrawPadClusters(oss); + // pRowClusters->DrawRowTpcClusters(oss); + pTpcRowClusters->Draw(oss); + pTpcRowClusters->Draw(oss, true); + + Print(oss); + std::cout << oss.str(); + oss.str(""); // clear string buffer + oss.clear(); +#endif + // break; // just one row + } + } + inline void Add(TpcRowClusters *pTpcRowClusters) { _lstTpcRowClusters.push_back(pTpcRowClusters); } + inline void Load(uint nRow, uint nPad, uint nTimeBin, float fAdc) + { + if (_lstTpcRowClusters.empty() || _lstTpcRowClusters.back()->getRow() != nRow) + _lstTpcRowClusters.push_back(new TpcRowClusters(nRow, nPad, nTimeBin, fAdc)); + else + _lstTpcRowClusters.back()->Load(nPad, nTimeBin, fAdc); + } + inline void Load(uint nRow, uint nPad, TpcAdcTimebin *pTimeBins, uint nSize) + { + if (_lstTpcRowClusters.empty() || _lstTpcRowClusters.back()->getRow() != nRow) + _lstTpcRowClusters.push_back(new TpcRowClusters(nRow)); + for (uint i = 0; i < nSize; i++) { + TpcAdcTimebin &r = pTimeBins[i]; + _lstTpcRowClusters.back()->Load(nPad, r._nTimebin, r._nAdc); + } + } + void Type(std::ostringstream &oss) const + { + std::string strRows = std::string(nTPC_SECTOR_ROWS, '-'); + for (std::list<TpcRowClusters *>::const_iterator i = _lstTpcRowClusters.begin(); i != _lstTpcRowClusters.end(); + i++) { + uint nRow = (*i)->getRow(); + strRows[nRow] = szDECADE[nRow % 10]; + } + oss << " Rows: " << strRows << std::endl; + } + void Print(std::ostringstream &oss) const + { + oss << "TpcClusters:: sector: " << getSector() << std::endl; + for (std::list<TpcRowClusters *>::const_iterator i = _lstTpcRowClusters.begin(); i != _lstTpcRowClusters.end(); + i++) + (*i)->Print(oss); + } + void PrintPadClusters(std::ostringstream &oss) const + { + oss << "TpcClusters:: sector: " << getSector() << std::endl; + for (std::list<TpcRowClusters *>::const_iterator i = _lstTpcRowClusters.begin(); i != _lstTpcRowClusters.end(); + i++) + (*i)->PrintPadClusters(oss); + } + void DrawPadClusters(std::ostringstream &oss) const + { + oss << "TpcClusters:: sector: " << getSector() << std::endl; + for (std::list<TpcRowClusters *>::const_iterator i = _lstTpcRowClusters.begin(); i != _lstTpcRowClusters.end(); + i++) + (*i)->DrawPadClusters(oss); + } + void Draw(std::ostringstream &oss) const + { + oss << "TpcClusters:: sector: " << getSector() << std::endl; + for (std::list<TpcRowClusters *>::const_iterator i = _lstTpcRowClusters.begin(); i != _lstTpcRowClusters.end(); + i++) + (*i)->Draw(oss); + } +}; + +} // namespace TpcClustering + +#endif // TPC_CLUSTERS_H \ No newline at end of file diff --git a/detectors/tpc/clusterHitFinder/fast/clustering/TpcEventPayload.h b/detectors/tpc/clusterHitFinder/fast/clustering/TpcEventPayload.h new file mode 100644 index 00000000..7c5358cf --- /dev/null +++ b/detectors/tpc/clusterHitFinder/fast/clustering/TpcEventPayload.h @@ -0,0 +1,144 @@ +// +// Joint Institute for Nuclear Research (JINR), Dubna, Russia, June-2024 +// TPC Event Payload structures for MPD/NICA project +// author: Viktor Krylov JINR/LHEP +// email: kryman@jinr.ru +// +#ifndef TPC_EVENT_PAYLOAD_H +#define TPC_EVENT_PAYLOAD_H + +#include <iostream> +#include <iomanip> +#include <fstream> +#include <vector> + +#include "TpcSampaPack.h" + +namespace TpcClustering { + +struct TpcPadPayload { + uint _nPad; // row pad number [0-123] + std::vector<TpcAdcTimebin> _vAdcTimebins; + TpcPadPayload() : _nPad(0) {} + TpcPadPayload(const TpcSampaPack &rSP) + : _nPad(rSP._nPad), _vAdcTimebins(rSP._TpcAdcTimebins, rSP._TpcAdcTimebins + rSP._nSize) + { + } + ~TpcPadPayload() { _vAdcTimebins.clear(); } + void Print() const + { + std::printf("%4d: ", _nPad); + for (uint i = 0; i < nSAMPA_PAYLOAD_WORDS; i++) { + auto it = std::find_if(_vAdcTimebins.begin(), _vAdcTimebins.end(), + [i](const TpcAdcTimebin &r) { return r._nTimebin == i; }); + std::cout << (it != _vAdcTimebins.end() ? '+' : '-'); + } + std::printf("\n"); + } +}; +struct TpcRowPayload { + uint _nRow; // sector row number: [0-52] + std::vector<TpcPadPayload> _vPads; + TpcRowPayload() : _nRow(0) {} + TpcRowPayload(const TpcSampaPack &r) : _nRow(r._nRow), _vPads{TpcPadPayload(r)} {} + ~TpcRowPayload() { _vPads.clear(); } + void Push(const TpcSampaPack &r) + { + assert(_vPads.back()._nPad != r._nPad); + + if (!r.isEmpty()) _vPads.push_back(TpcPadPayload(r)); + } + void Print(uint nSize) const + { + for (uint i = 0; i < nSize; i++) { + auto it = std::find_if(_vPads.begin(), _vPads.end(), [i](const TpcPadPayload &r) { return r._nPad == i; }); + std::cout << (it != _vPads.end() ? '+' : '-'); + } + std::cout << std::endl; + } + void Print() const + { + uint nSize = anTPC_ROW_PADS[_nRow]; + for (uint i = 0; i < nSize; i++) { + auto it = std::find_if(_vPads.begin(), _vPads.end(), [i](const TpcPadPayload &r) { return r._nPad == i; }); + if (it != _vPads.end()) (*it).Print(); + // else + // std::cout << std::setw(4) << i << ": " << std::string(nSAMPA_PAYLOAD_WORDS, '-') << std::endl; + } + } +}; +struct TpcSectorPayload { + uint _nSector; // sector number: [0-23] + std::vector<TpcRowPayload> _vRows; + TpcSectorPayload() : _nSector(0) {} + TpcSectorPayload(uint nSector) : _nSector(nSector) {} + ~TpcSectorPayload() { _vRows.clear(); } + void Push(const TpcSampaPack &r) + { + if (_vRows.empty() || _vRows.back()._nRow != r._nRow) + _vRows.push_back(TpcRowPayload(r)); + else + _vRows.back().Push(r); + } + void Print(uint nRow) const + { + std::cout << " Sector: " << _nSector; + if (nRow < nTPC_SECTOR_ROWS) { + std::cout << " Row: " << nRow; + auto it = + std::find_if(_vRows.begin(), _vRows.end(), [nRow](const TpcRowPayload &r) { return r._nRow == nRow; }); + if (it != _vRows.end()) { + std::cout << std::endl; + (*it).Print(); + } else + std::cout << " <all pads are empty!>" << std::endl; + } else { + std::cout << std::endl; + for (uint i = 0; i < nTPC_SECTOR_ROWS; i++) { + uint nSize = anTPC_ROW_PADS[i]; + uint nBlank = (nTPC_ROW_MAX_PADS - nSize) / 2; + std::cout << std::setw(4) << i << ": " << std::string(nBlank, ' '); + + auto it = std::find_if(_vRows.begin(), _vRows.end(), [i](const TpcRowPayload &r) { return r._nRow == i; }); + if (it != _vRows.end()) + (*it).Print(nSize); + else + std::cout << std::string(nSize, '-') << std::endl; + } + } + } +}; +struct TpcEventPayload { + uint _nEvent; // event number + std::vector<const TpcSectorPayload *> _vpSectors; + TpcEventPayload() : _nEvent(0) {} + TpcEventPayload(uint nEvent) : _nEvent(nEvent) {} + ~TpcEventPayload() { Clear(); } + void Clear() + { + for (const TpcSectorPayload *pTpcSectorPayload : _vpSectors) delete pTpcSectorPayload; + _vpSectors.clear(); + } + const TpcSectorPayload *Find(uint nSector) const + { + std::vector<const TpcSectorPayload *>::const_iterator iBuf = + std::find_if(_vpSectors.begin(), _vpSectors.end(), + [nSector](const TpcSectorPayload *pBuf) { return pBuf->_nSector == nSector; }); + return iBuf != _vpSectors.end() ? *iBuf : NULL; + } + void Print(uint nSector, uint nRow) const + { + assert(nSector < nTPC_SECTORS); + + auto it = std::find_if(_vpSectors.begin(), _vpSectors.end(), + [nSector](const TpcSectorPayload *p) { return p->_nSector == nSector; }); + if (it != _vpSectors.end()) + (*it)->Print(nRow); + else + std::cout << " Sector: " << nSector << " <completely empty!>" << std::endl; + } +}; + +} // namespace TpcClustering + +#endif // TPC_EVENT_PAYLOAD_H \ No newline at end of file diff --git a/detectors/tpc/clusterHitFinder/fast/clustering/TpcHit.h b/detectors/tpc/clusterHitFinder/fast/clustering/TpcHit.h new file mode 100644 index 00000000..d6bd402d --- /dev/null +++ b/detectors/tpc/clusterHitFinder/fast/clustering/TpcHit.h @@ -0,0 +1,108 @@ +// +// Joint Institute for Nuclear Research (JINR), Dubna, Russia, Jun-2024 +// TPCHit object for MPD/NICA project +// author: Viktor A.Krylov JINR/LNP +// email: kryman@jinr.ru +// +#ifndef TPC_HIT_H +#define TPC_HIT_H + +#include <cassert> +#include <sstream> +#include <cmath> + +#include "TpcAdcGauss.h" + +namespace TpcClustering { + +// const uint nADC_DISCRECITY = 40; // should be more than 9 +// #define HEX(ch) std::setw(2) << std::setfill('0') << std::hex << int(ch) +static const char *SZ_ALPHABET = "0123456789ABCDEFJHIJKLMNOPQRSTUVWXYZ"; + +class TpcHit { + uint _nId; // cluster id + float _fTimeBin; // cluster timebin hit + float _fTimeBinDev; // deviation of the timebin + float _fPad; // cluster pad hit + float _fPadDev; // deviation of the pad + float _fAdcMax; // adc max value + float _fAdcSum; // total ADC sum of all pad clusters + bool _bUndefined; // is not well defined or it has ultimate position in a tpc sector +public: + TpcHit() : _nId(0), _fTimeBin(0), _fTimeBinDev(0), _fPad(0), _fPadDev(0), _fAdcSum(0), _bUndefined(false) {} + TpcHit(uint nId, float fAdc) + : _nId(nId), _fTimeBin(0), _fTimeBinDev(0), _fPad(0), _fPadDev(0), _fAdcSum(fAdc), _bUndefined(false) + { + } + virtual ~TpcHit() {} + inline void setId(uint nId) { _nId = nId; } + inline uint getId() const { return _nId; } + inline void setTimeBin(float fTimeBin) { _fTimeBin = fTimeBin; } + inline float getTimeBin() const { return _fTimeBin; } + inline void setTimeBinDev(float fDev) { _fTimeBinDev = fDev; } + inline float getTimeBinDev() const { return _fTimeBinDev; } + inline void setPad(float fPad) { _fPad = fPad; } + inline float getPad() const { return _fPad; } + inline void setPadDev(float fDev) { _fPadDev = fDev; } + inline float getPadDev() const { return _fPadDev; } + inline void setAdcMax(float fAdc) { _fAdcMax = fAdc; } + inline float getAdcMax() const { return _fAdcMax; } + + inline void setAdcSum(float fAdc) { _fAdcSum = fAdc; } + inline float getAdcSum() const { return _fAdcSum; } + inline void addAdcSum(float fAdc) { _fAdcSum += fAdc; } + inline void reduceAdcSum(float fAdc) + { + assert(fAdc < _fAdcSum); + _fAdcSum -= fAdc; + } + + inline bool isUndefined() const { return _bUndefined; } + inline void setUndefined() { _bUndefined = true; } + void Calc(const std::vector<float> &rvfAdc, uint nRow, uint nPad, float fTimeBin, float fTimeBinLast) + { + size_t nSize = rvfAdc.size(); + assert(nSize); + + TpcAdcGauss gauss; + gauss.calc(rvfAdc); + float fMean = gauss.getMean(); + if (fMean < 0.0f || fMean > nSize - 1) { + fMean = gauss.getAdjustment(); // gauss adjustment for fitting problems + if (fMean == 0.0f) + TpcHit::setUndefined(); + else if (fMean < 0.0f) { + fMean = 0.0f; + _fAdcMax = rvfAdc[0]; + } else if (nSize > 1 && fMean > nSize - 1) { + fMean = nSize - 1; + _fAdcMax = rvfAdc[nSize - 1]; + } + } else + _fAdcMax = gauss.getMax(); + _fPadDev = gauss.getDeviation(); + _fPad = nPad + fMean; + float fScale = nSize > 1 ? fMean / (nSize - 1) : 1.0f; + float fWidth = fTimeBinLast - fTimeBin; // timebin width +/- + _fTimeBin = fTimeBin + fScale * fWidth; + _fTimeBinDev = fabs(fWidth / nSize) * _fPadDev; + + assert(_fPad < anTPC_ROW_PADS[nRow] && _fTimeBin < nTPC_TIMEBINS); + } + std::ostringstream &Print(std::ostringstream &oss) const + { + std::streamsize ss = std::cout.precision(); + oss << std::setprecision(2); + char ch = SZ_ALPHABET[getId() % strlen(SZ_ALPHABET)]; + oss << "Cluster:: id:" << _nId << "/'" << ch << "', adcSum:" << _fAdcSum << ", adcMax:" << _fAdcMax + << ", timebin:" << _fTimeBin << "(+/-)" << _fTimeBinDev << ", pad:" << _fPad << "(+/-)" << _fPadDev; + if (isUndefined()) oss << ", undefined"; + oss << std::endl; + oss << std::setprecision(ss); + return oss; + } +}; + +} // namespace TpcClustering + +#endif // TPC_HIT_H \ No newline at end of file diff --git a/detectors/tpc/clusterHitFinder/fast/clustering/TpcPadCluster.h b/detectors/tpc/clusterHitFinder/fast/clustering/TpcPadCluster.h new file mode 100644 index 00000000..b65669fc --- /dev/null +++ b/detectors/tpc/clusterHitFinder/fast/clustering/TpcPadCluster.h @@ -0,0 +1,308 @@ +// +// Joint Institute for Nuclear Research (JINR), Dubna, Russia, Jun-2024 +// TPC PadCluster class for MPD/NICA project +// author: Viktor A.Krylov JINR/LNP +// email: kryman@jinr.ru +// +#ifndef TPC_PAD_CLUSTER_H_ +#define TPC_PAD_CLUSTER_H_ + +#include <vector> +#include <cstring> +#include <algorithm> +#include <cassert> + +#include "TpcPadHit.h" + +#define _RED_ "0;31" +#define _GREEN_ "0;32" +#define _YELLOW_ "0;33" +#define _BASH_CODE_(_CODE) "\e[" _CODE "m" +#define _RED_MSG_(_MSG) _BASH_CODE_(_RED_) _MSG _BASH_CODE_() + +namespace TpcClustering { + +class TpcCluster; +// const float fADC_DOORSTEP = 2.5; // adc sensitivity + +class TpcPadCluster : public TpcPadHit { + TpcCluster *_pTpcCluster; // link to a TPC Cluster which the PadCluster belongs + + uint _nTimeBin; // timebin value correspond to _vfAdc[0] + uint _iAdcMax; // index of the max adc value in the Adc vector + uint _iAdcMin; // index of the min adc value in the Adc vector +public: + std::vector<float> _vfAdc; // vector of ADC hits + TpcPadCluster() : TpcPadHit(), _pTpcCluster(NULL), _iAdcMax(0), _nTimeBin(0) {} + TpcPadCluster(uint nPad, uint nTimeBin, float fAdc) + : TpcPadHit(nPad, fAdc), _nTimeBin(nTimeBin), _pTpcCluster(NULL), _iAdcMin(0), _iAdcMax(0) + { + assert(fAdc > 0.0f); + + _vfAdc.push_back(fAdc); + } + TpcPadCluster(const TpcPadCluster &rTpcPadCluster) : TpcPadHit(rTpcPadCluster) + { + _nTimeBin = rTpcPadCluster.getTimeBinMin(); + _iAdcMax = rTpcPadCluster.getAdcMaxPos(); + _iAdcMin = rTpcPadCluster.getAdcMinPos(); + _vfAdc = rTpcPadCluster._vfAdc; + for (float fAdc : _vfAdc) assert(fAdc > 0.0f && fAdc < 65535.0f); + _pTpcCluster = rTpcPadCluster.getCluster(); + } + virtual ~TpcPadCluster() { _vfAdc.clear(); } + float operator[](size_t i) const + { + assert(i < _vfAdc.size()); + return _vfAdc[i]; + } + inline uint getClusterId(const TpcCluster *pTpcCluster) const; + inline float getClusterTime(const TpcCluster *pTpcCluster) const; + inline float getClusterPad(const TpcCluster *pTpcCluster) const; + inline TpcCluster *getCluster() const { return _pTpcCluster; } + inline void setCluster(TpcCluster *pCluster) { _pTpcCluster = pCluster; } + + inline size_t getSize() const { return _vfAdc.size(); } + inline uint getTimeBinAdcMax() const { return _nTimeBin + _iAdcMax; } + + inline uint getTimeBinMin() const { return _nTimeBin; } + inline uint getTimeBinMax() const { return _nTimeBin + _vfAdc.size() - 1; } + inline uint getAdcMinPos() const { return _iAdcMin; } + inline uint getAdcMaxPos() const { return _iAdcMax; } + + inline void Split(float fScale) + { + assert(fScale > 0); + + TpcPadHit::Split(fScale); + for (uint i = 0; i < _vfAdc.size(); i++) _vfAdc[i] *= fScale; + } + // for joint clusters + inline bool isOut(const TpcPadCluster *pTpcPadCluster) const { return pTpcPadCluster->getPad() > getPad() + 1; } + inline bool isBefore(const TpcPadCluster *pTpcPadCluster) const + { + return pTpcPadCluster->getTimeBinAdcMax() > getTimeBinAdcMax() + 1; + } + inline bool isAfter(const TpcPadCluster *pTpcPadCluster) const + { + return pTpcPadCluster->getTimeBinAdcMax() + 1 < getTimeBinAdcMax(); + } + inline bool isIn(const TpcPadCluster *pTpcPadCluster) const + { + return pTpcPadCluster->getTimeBinAdcMax() + 1 >= getTimeBinAdcMax() && + pTpcPadCluster->getTimeBinAdcMax() <= getTimeBinAdcMax() + 1; + } + // for cluster split + inline bool isOut(uint nTimeBin) const { return nTimeBin > getTimeBinMax() + 1; } + inline bool isIn(float fAdc) const { return fAdc < _vfAdc.back() || _iAdcMax + 1 == _vfAdc.size(); } + inline void addAdcHit(float fAdc) + { + assert(fAdc > 0.0f); + + float fAdcBack = _vfAdc.back(); + if (fAdc < fAdcBack) + _iAdcMin = _vfAdc.size(); + else if (fAdc >= fAdcBack) // =! + _iAdcMax = _vfAdc.size(); + addAdcSum(fAdc); + _vfAdc.push_back(fAdc); + } + TpcPadCluster *splitAdcHit(uint nPad, uint nTimeBin, float fAdc) + { // proportionally split AdcHit value between two combined clusters; + assert(_vfAdc.size() > 1 && _iAdcMin > 0 && fAdc > 0.0f); + + float fAdcMin = _vfAdc[_iAdcMin]; + size_t nSize = _vfAdc.size(); + uint nMin = nSize - _iAdcMin; // position of the 1st adc min hit downstream + if (nMin == 1) { // simple minimum /\./ + float fAdcLeft = _vfAdc[_iAdcMin - 1]; // a prior Adc hit value + float fAdcRight = fAdcMin * fAdc / (fAdcLeft + fAdc); + _vfAdc.back() -= fAdcRight; + reduceAdcSum(fAdcRight); + TpcPadCluster *pTpcPadCluster = new TpcPadCluster(nPad, nTimeBin - 1, fAdcRight); + pTpcPadCluster->addAdcHit(fAdc); + return pTpcPadCluster; + } else { // udefined minimum /\.../ + assert(nMin); + + uint nMin2 = nMin >> 1; + _vfAdc.resize(nSize - nMin2); + reduceAdcSum(fAdcMin * nMin2); + float fAdcMean = fAdcMin; // middle hit adc value + if (nMin & 0x1) { // odd: 3 hits and more + fAdcMean /= 2; + _vfAdc.back() = fAdcMean; + reduceAdcSum(fAdcMean); + nMin2++; + } + TpcPadCluster *pTpcPadCluster = new TpcPadCluster(nPad, nTimeBin - nMin2, fAdcMean); + while (--nMin2) pTpcPadCluster->addAdcHit(fAdcMin); + pTpcPadCluster->addAdcHit(fAdc); + return pTpcPadCluster; + } + } + void calcGauss() + { + size_t nSize = _vfAdc.size(); + assert(nSize); + + TpcAdcGauss gauss; + gauss.calc(_vfAdc); + float fMean = gauss.getMean(); + if (fMean < 0.0f || fMean > nSize - 1) { + fMean = gauss.getAdjustment(); // gauss adjustment for fitting problems + if (fMean == 0.0f) + TpcPadHit::setUndefined(); + else if (fMean < 0.0f) + fMean = 0.0f; + else if (fMean > nSize - 1) + fMean = nSize - 1; + TpcPadHit::setAdcMax(_vfAdc[_iAdcMax]); // it is not correct in normal case but sometimes fAdcMax is far more + // then _vfAdc[_iAdcMax] + } else + TpcPadHit::setAdcMax(gauss.getMax()); + TpcPadHit::setTimeBin(_nTimeBin + fMean); + float fTimeBin = TpcPadHit::getTimeBin(); + if (fTimeBin <= 0.0f || fTimeBin >= nTPC_TIMEBINS - 1) TpcPadHit::setUndefined(); + TpcPadHit::setDev(gauss.getDeviation()); + } + void setBinContent(TH2F &rH2F) const + { + size_t nSize = _vfAdc.size(); + for (uint i = 0; i < nSize; i++) { + uint nPad = getPad() + 1; + uint nTimeBin = _nTimeBin + i + 1; + rH2F.SetBinContent(nPad, nTimeBin, rH2F.GetBinContent(nPad, nTimeBin) + _vfAdc[i]); + // std::cout << " pad: " << getPad() << " timebin: " << _nTimeBin + i << " adc: " << _vfAdc[i] << std::endl; + } + } + void setBinContent(TH2S &rH2S) const + { + size_t nSize = _vfAdc.size(); + for (uint i = 0; i < nSize; i++) { + uint nPad = getPad() + 1; + uint nTimeBin = _nTimeBin + i + 1; + rH2S.SetBinContent(nPad, nTimeBin, rH2S.GetBinContent(nPad, nTimeBin) + _vfAdc[i]); + // std::cout << " pad: " << getPad() << " timebin: " << _nTimeBin + i << " adc: " << _vfAdc[i] << std::endl; + } + } + std::ostringstream &Draw(std::ostringstream &oss, char ch = '+') const + { + uint nTimeBin = _pTpcCluster ? lroundf(getClusterTime(_pTpcCluster)) : UINT_MAX; + uint nPad = _pTpcCluster ? lroundf(getClusterPad(_pTpcCluster)) : UINT_MAX; + uint nMax = getTimeBinAdcMax(); + size_t nSize = _vfAdc.size(); + for (uint i = 0; i < nSize; i++) { + uint n = _nTimeBin + i; + if (n == nTimeBin && nPad == TpcPadHit::getPad()) + // oss << _BASH_CODE_(_RED_) << ch << _BASH_CODE_(); // red cluster hit + oss << _RED_MSG_("*"); // red cluster hit star '*' + else if (n == nMax && nSize > 1) + // oss << _BASH_CODE_(_GREEN_) << ch << _BASH_CODE_(); // yellow pad cluster hit + oss << '^'; + else + oss << ch; + } + // oss << std::endl; + return oss; + } + std::ostringstream &Draw(std::ostringstream &oss, uint nPos) const + { + assert(_vfAdc.size()); + + uint nTimeBin = _pTpcCluster ? lroundf(getClusterTime(_pTpcCluster)) : UINT_MAX; + uint nPad = _pTpcCluster ? lroundf(getClusterPad(_pTpcCluster)) : UINT_MAX; + uint nMax = getTimeBinAdcMax(); + float fAdcMax = _vfAdc[_iAdcMax]; + float fAdcMin = _vfAdc[_iAdcMin]; + if (fAdcMin > _vfAdc[0]) fAdcMin = _vfAdc[0]; + uint iLeft = _iAdcMax; + uint iRight = _iAdcMax; + bool bLeft = true; // the end of the left arm if false + bool bRight = true; // the end of the right arm if false + size_t nSize = _vfAdc.size(); + float fBin = nSize > 1 ? (fAdcMax - fAdcMin) / (nSize - 1) : 0; + uint i = nSize; + while (i-- && (bLeft || bRight)) { + std::string str(nSize, ' '); + float fAdcLeft = FLT_MAX; + float fAdcRight = FLT_MAX; + if (iLeft == iRight) { + str[_iAdcMax] = '^'; + iLeft ? iLeft-- : bLeft = false; + iRight + 1 < nSize ? iRight++ : bRight = false; + } + float fAdc = fAdcMin + fBin * i; + if (bLeft) { + bool b = false; + while (fBin == 0.0f || _vfAdc[iLeft] >= fAdc) { + str[iLeft] = b ? '_' : '/'; + fAdcLeft = b ? 0 : _vfAdc[iLeft]; + b = true; + if (!iLeft) { + bLeft = false; + break; + } else + iLeft--; + } + } + if (bRight) { + bool b = false; + while (fBin == 0.0f || _vfAdc[iRight] >= fAdc) { + str[iRight] = b ? '_' : '\\'; + fAdcRight = b ? 0 : _vfAdc[iRight]; + b = true; + if (iRight + 1 == nSize) { + bRight = false; + break; + } else + iRight++; + } + } + if (str.find_first_not_of(' ') != std::string::npos) { + if (fAdcLeft == FLT_MAX && fAdcRight == FLT_MAX) + fAdc = fAdcMax; + else if (fAdcLeft != FLT_MAX && fAdcLeft != 0 && fAdcRight == FLT_MAX) + fAdc = fAdcLeft; + else if (fAdcRight != FLT_MAX && fAdcRight != 0 && fAdcLeft == FLT_MAX) + fAdc = fAdcRight; + std::streamsize ss = std::cout.precision(); + oss << std::string(nPos, ' ') << str << std::setw(7) << std::fixed << std::setprecision(2) << fAdc + << std::endl; + oss << std::setprecision(ss); + } + } + return oss; + } + std::ostringstream &Print(std::ostringstream &oss) const + { + uint nClusterId = getClusterId(_pTpcCluster); + oss << "PadCluster:: pad:" << getPad() << " timebins:[" << _nTimeBin << "-" << _nTimeBin + _vfAdc.size() - 1 + << "], adcSpikes:{ max[" << _iAdcMax << "],min[" << _iAdcMin << "] }, adcHits:{ "; + std::streamsize ss = std::cout.precision(); + oss << std::setprecision(0); + for (uint i = 0; i < _vfAdc.size(); i++) oss << _vfAdc[i] << ' '; + oss << '}' << std::endl; + oss << std::setprecision(ss); + oss << std::string(6, ' '); + TpcPadHit::Print(oss); + return oss; + } + void Print() const + { + uint nClusterId = getClusterId(_pTpcCluster); + std::cout << "PadCluster:: pad:" << getPad() << " timebins:[" << _nTimeBin << "-" << _nTimeBin + _vfAdc.size() - 1 + << "], adcSpikes:{ max[" << _iAdcMax << "],min[" << _iAdcMin << "] }, adcHits:{ "; + std::streamsize ss = std::cout.precision(); + std::cout << std::setprecision(0) << std::fixed; + for (uint i = 0; i < _vfAdc.size(); i++) std::cout << _vfAdc[i] << ' '; + std::cout << std::setprecision(ss); + std::cout << '}' << std::endl; + std::cout << std::string(6, ' '); + TpcPadHit::Print(); + } +}; + +} // namespace TpcClustering + +#endif // TPC_PAD_CLUSTER_H \ No newline at end of file diff --git a/detectors/tpc/clusterHitFinder/fast/clustering/TpcPadHit.h b/detectors/tpc/clusterHitFinder/fast/clustering/TpcPadHit.h new file mode 100644 index 00000000..a83ba6bc --- /dev/null +++ b/detectors/tpc/clusterHitFinder/fast/clustering/TpcPadHit.h @@ -0,0 +1,79 @@ +// +// Joint Institute for Nuclear Research (JINR), Dubna, Russia, Jun-2024 +// TPC PadHit object for MPD/NICA project +// author: Viktor A.Krylov JINR/LNP +// email: kryman@jinr.ru +// +#ifndef TPC_PAD_HIT_H +#define TPC_PAD_HIT_H + +#include <iosfwd> + +#include "TpcAdcGauss.h" + +namespace TpcClustering { + +class TpcPadHit { + uint _nPad; // Pad id number + float _fTimeBin; // pad cluster mean gauss timebin + float _fDev; // deviation of the timebin + float _fAdcMax; // gauss max value + float _fAdcSum; // total ADC sum of all hits + bool _bUndefined; // the hit is not well defined or it has ultimate position in the sector +public: + TpcPadHit() : _nPad(UINT_MAX), _fTimeBin(0), _fDev(0), _fAdcMax(0), _fAdcSum(0), _bUndefined(false) {} + TpcPadHit(uint nPad, float fAdc) : _nPad(nPad), _fTimeBin(0), _fAdcSum(fAdc), _bUndefined(false) {} + TpcPadHit(const TpcPadHit &rPadHit) + : _nPad(rPadHit.getPad()), _fTimeBin(rPadHit.getTimeBin()), _fDev(rPadHit.getDev()), + _fAdcMax(rPadHit.getAdcMax()), _fAdcSum(rPadHit.getAdcSum()), _bUndefined(rPadHit.isUndefined()) + { + } + ~TpcPadHit() {} + inline uint getPad() const { return _nPad; } + inline void setTimeBin(float f) { _fTimeBin = f; } + inline float getTimeBin() const { return _fTimeBin; } + inline void setDev(float f) { _fDev = f; } + inline float getDev() const { return _fDev; } + inline void setAdcMax(float f) + { + assert(f > 0.0f && f < 65535.0f); + _fAdcMax = f; + } + inline float getAdcMax() const { return _fAdcMax; } + inline void setAdcSum(float f) { _fAdcSum = f; } + inline float getAdcSum() const { return _fAdcSum; } + inline void addAdcSum(float f) { _fAdcSum += f; } + inline void reduceAdcSum(float f) { _fAdcSum -= f; } + inline bool isUndefined() const { return _bUndefined; } + inline void setUndefined() { _bUndefined = true; } + inline void Split(float f) + { + _fAdcMax *= f; + _fAdcSum *= f; + } + std::ostringstream &Print(std::ostringstream &oss) const + { + std::streamsize ss = std::cout.precision(); + oss << std::setprecision(2); + oss << "TpcPadHit:: timebin:" << _fTimeBin << "(+/-)" << _fDev << ", adcMax:" << _fAdcMax + << ", adcSum:" << _fAdcSum; + if (isUndefined()) oss << ", undefined"; + oss << std::endl; + oss << std::setprecision(ss); + return oss; + } + void Print() const + { + std::streamsize ss = std::cout.precision(); + std::cout << std::setprecision(2); + std::cout << "TpcPadHit:: timebin:" << _fTimeBin << "(+/-)" << _fDev << ", adcMax:" << _fAdcMax + << ", adcSum:" << _fAdcSum; + std::cout << std::setprecision(ss); + if (isUndefined()) std::cout << ", undefined"; + std::cout << std::endl; + } +}; + +} // namespace TpcClustering + +#endif // TPC_PAD_HIT_H \ No newline at end of file diff --git a/detectors/tpc/clusterHitFinder/fast/clustering/TpcRowClusters.h b/detectors/tpc/clusterHitFinder/fast/clustering/TpcRowClusters.h new file mode 100644 index 00000000..30def841 --- /dev/null +++ b/detectors/tpc/clusterHitFinder/fast/clustering/TpcRowClusters.h @@ -0,0 +1,406 @@ +// +// Joint Institute for Nuclear Research (JINR), Dubna, Russia, Jun-2024 +// TPC RowClusters object for MPD/NICA project +// author: Viktor A.Krylov JINR/LNP +// email: kryman@jinr.ru +// +#ifndef TPC_ROW_CLUSTERS_H_ +#define TPC_ROW_CLUSTERS_H_ + +#include <iomanip> +#include <ios> +#include <sstream> +#include <list> + +#include <TDirectory.h> + +#include "TpcSectorGeo.h" +#include "TpcAdcStat.h" +#include "TpcCluster.h" + +namespace TpcClustering { + +static const char *szDECADE = "0123456789"; +inline void DrawTimeBinTopScale(std::ostringstream &oss, uint nSize) +{ + oss << std::setw(4) << ":"; + for (uint i = 0; i < nSize; i++) oss << std::left << std::setw(10) << i * 10; + oss << std::endl; + oss << std::internal << std::setw(4) << ":"; + for (uint i = 0; i < nSize; i++) oss << szDECADE; + oss << std::internal << std::endl; +}; +inline void DrawTimeBinBottomScale(std::ostringstream &oss, uint nSize) +{ + oss << std::setw(4) << ":"; + for (uint i = 0; i < nSize; i++) oss << szDECADE; + oss << std::endl << std::internal << std::setw(4) << ":"; + for (uint i = 0; i < nSize; i++) oss << std::left << std::setw(10) << i * 10; + oss << std::internal << std::endl; +}; + +class TpcRowClusters { + uint _nRow; // sector row number + uint _nClusterId; // current cluster id for internal usage + TpcAdcStat _tpcAdcStat; // row statistic (only for output, nothing more!) + std::list<TpcCluster *> _lstTpcClusters; + +public: + std::list<TpcPadCluster *> _lstTpcPadClusters; // list of fired ADC hits + std::vector<std::list<TpcPadCluster *>::iterator> _viTpcPadClusters; + TpcRowClusters() : _nRow(0), _nClusterId(0), _tpcAdcStat() {} + TpcRowClusters(uint nRow) : _nRow(nRow), _nClusterId(0), _tpcAdcStat() {} + TpcRowClusters(uint nRow, uint nPad, uint nTimeBin, float fAdc) + : _nRow(nRow), _nClusterId(0), _tpcAdcStat(nPad, nTimeBin, fAdc) + { + _lstTpcPadClusters.push_back(new TpcPadCluster(nPad, nTimeBin, fAdc)); + _viTpcPadClusters.push_back(_lstTpcPadClusters.begin()); + } + virtual ~TpcRowClusters() + { + std::for_each(_lstTpcPadClusters.begin(), _lstTpcPadClusters.end(), [this](const TpcPadCluster *pTpcPadCluster) { + if (pTpcPadCluster->getCluster() == NULL) delete pTpcPadCluster; + // delete pPadCluster->getCluster(); + }); + _lstTpcPadClusters.clear(); + _viTpcPadClusters.clear(); + std::for_each(_lstTpcClusters.begin(), _lstTpcClusters.end(), + [this](const TpcCluster *pTpcCluster) { delete pTpcCluster; }); + _lstTpcClusters.clear(); + } + // inline bool operator==(uint nRow) const {return _nRow == nRow; } + inline uint getRow() const { return _nRow; } + inline uint getClusterId() const { return _nClusterId; } + inline const std::list<TpcCluster *> &getClusters() const { return _lstTpcClusters; } + bool isAllPadClustersBelongTo() + { + for (const TpcPadCluster *pTpcPadCluster : _lstTpcPadClusters) + if (pTpcPadCluster->getCluster() == NULL) return true; + return false; + } + inline bool isEmpty() const { return _lstTpcClusters.empty(); } + void defDigits(const std::vector<TpcDigit> *pvTpcDigits) + { + for (std::list<TpcCluster *>::const_iterator i = _lstTpcClusters.begin(); i != _lstTpcClusters.end(); i++) + (*i)->defDigits(pvTpcDigits, _nRow); + } + void defAdcHits(uint nSector, std::vector<TpcAdcHit> &rvTpcAdcHits) + { + for (TpcCluster *pTpcCluster : _lstTpcClusters) { + if (pTpcCluster->isUndefined()) // remove all ill-defined tpc clusters! + continue; + TpcAdcHit tpcAdcHit(_nRow, pTpcCluster->getId(), pTpcCluster->getAdcSum()); + tpcAdcHit = TpcPoint(_nRow, pTpcCluster->getPad(), pTpcCluster->getTimeBin()) + .getGlobal(nSector); // convert to global coords + tpcAdcHit.shiftSampaDelay(nSector, fTAU_SAMPA_DELAYns * fDRIFT_VELOCITYcmns); + rvTpcAdcHits.push_back(tpcAdcHit); + } + } + inline void Load(uint nPad, uint nTimeBin, float fAdc) + { + _tpcAdcStat.PickOut(nPad, nTimeBin, fAdc); + if (_lstTpcPadClusters.empty()) { + _lstTpcPadClusters.push_back(new TpcPadCluster(nPad, nTimeBin, fAdc)); + _viTpcPadClusters.push_back(_lstTpcPadClusters.begin()); + } else { + TpcPadCluster *pTpcPadCluster = _lstTpcPadClusters.back(); + if (nPad != pTpcPadCluster->getPad()) { // new PadCluster in the row + _lstTpcPadClusters.push_back(new TpcPadCluster(nPad, nTimeBin, fAdc)); + _viTpcPadClusters.push_back(--_lstTpcPadClusters.end()); + } else if (pTpcPadCluster->isOut(nTimeBin)) // new PadCluster in the pad + _lstTpcPadClusters.push_back(new TpcPadCluster(nPad, nTimeBin, fAdc)); + else if (pTpcPadCluster->isIn(fAdc)) // add new adc hit + pTpcPadCluster->addAdcHit(fAdc); + else // split adc hit + _lstTpcPadClusters.push_back(pTpcPadCluster->splitAdcHit(nPad, nTimeBin, fAdc)); + } + } + void calcPadClusters() + { + assert(_lstTpcPadClusters.size()); + + for (TpcPadCluster *pTpcPadCluster : _lstTpcPadClusters) pTpcPadCluster->calcGauss(); + } + void Joint() + { + assert(_viTpcPadClusters.size()); + + size_t nSize = _viTpcPadClusters.size(); + for (uint i = 0; i < nSize; i++) { + std::list<TpcPadCluster *>::iterator iTpcPadCluster = _viTpcPadClusters[i]; + uint i1 = i + 1; + std::list<TpcPadCluster *>::iterator iTpcPadCluster1st = + i1 < nSize ? _viTpcPadClusters[i1] : _lstTpcPadClusters.end(); + if (iTpcPadCluster1st == _lstTpcPadClusters.end() || + (*iTpcPadCluster)->isOut(*iTpcPadCluster1st)) // gap between pads or it is a final pad in the row + do { + const TpcPadCluster *pTpcPadCluster = *iTpcPadCluster; + if (pTpcPadCluster->getCluster() == nullptr) // is it not belong to a cluster yet? + _lstTpcClusters.push_back(new TpcCluster(_nClusterId++, pTpcPadCluster)); + } while (++iTpcPadCluster != iTpcPadCluster1st); + else { + uint i2 = i + 2; + std::list<TpcPadCluster *>::iterator iTpcPadCluster2nd = + i2 < nSize ? _viTpcPadClusters[i2] : _lstTpcPadClusters.end(); + std::list<TpcPadCluster *>::iterator iTpcPadClusterNext = iTpcPadCluster1st; // next pad + while (iTpcPadCluster != iTpcPadCluster1st) { + TpcPadCluster *pTpcPadCluster = *iTpcPadCluster; + TpcCluster *pTpcCluster = pTpcPadCluster->getCluster(); + if (pTpcCluster == nullptr) // is it not defined yet? + _lstTpcClusters.push_back(pTpcCluster = new TpcCluster(_nClusterId++, pTpcPadCluster)); + if (pTpcPadCluster->isUndefined()) { // ill-defined pad cluster should be self-contained in the undefined + // tpc cluster! + iTpcPadCluster++; + continue; + } + while (iTpcPadClusterNext != iTpcPadCluster2nd) { // try to joint it with next pad cluster + const TpcPadCluster *pTpcPadClusterNext = *iTpcPadClusterNext; + if (pTpcPadClusterNext->isUndefined() || pTpcPadCluster->isAfter(pTpcPadClusterNext)) { + iTpcPadClusterNext++; // take next pad cluster in the pad + continue; + } else if (pTpcPadCluster->isBefore(pTpcPadClusterNext)) + break; // try to check new pad cluster + else { // joint to the tpc cluster + pTpcCluster->Add(pTpcPadClusterNext); + iTpcPadClusterNext++; + break; + } + } + iTpcPadCluster++; + } + } + } + assert(!isAllPadClustersBelongTo()); + + _lstTpcPadClusters + .clear(); // clear all pad clusters in the list (all pad clusters in the list should be belong to TpcClusters!) + _viTpcPadClusters.clear(); // clear all row pad iterators in the vector + } + void Split() + { + for (std::list<TpcCluster *>::iterator i = _lstTpcClusters.begin(); i != _lstTpcClusters.end(); i++) { + while (TpcCluster *pTpcCluster = (*i)->Split(_nRow)) { + pTpcCluster->setId(_nClusterId++); + _lstTpcClusters.insert(i, pTpcCluster); + } + } + } + void Store(TDirectory *pDir) + { + std::string strRow = std::to_string(_nRow); + TDirectory *pRowDir = pDir->mkdir(strRow.c_str(), ("row:" + strRow).c_str(), kTRUE); + if (pRowDir && pRowDir->cd()) { + for (const TpcCluster *pTpcCluster : _lstTpcClusters) { + std::string strId = std::to_string(pTpcCluster->getId()); + TDirectory *pClusterDir = pRowDir->mkdir(strId.c_str(), ("cluster:" + strId).c_str(), kTRUE); + if (pClusterDir && pClusterDir->cd()) { + pClusterDir->WriteObject(&pTpcCluster->_vnTpcDigits, "digits"); + // TpcHit* pHit = (TpcHit*)pTpcCluster; + TpcHit *pHit = (TpcHit *)pTpcCluster; + std::vector<uint8_t> vHit((uint8_t *)pHit, (uint8_t *)pHit + sizeof(TpcHit)); + pClusterDir->WriteObject(&vHit, "hit"); + } + } + } + } + // void Store( TDirectory* pDir ) { // for debug in order to save histograms in root file + // std::string strRow = std::to_string(_nRow); + // TDirectory* pRowDir = pDir->mkdir( strRow.c_str(), ("row:" + strRow).c_str(), kTRUE ); + // uint nPadMin = UINT_MAX; + // uint nPadMax = UINT_MAX; + // uint nTimeBinMin = UINT_MAX; + // uint nTimeBinMax = UINT_MAX; + // if( pRowDir && pRowDir->cd() ) { + // for( const TpcCluster* pTpcCluster : _lstTpcClusters ) { + // std::string strId = std::to_string(pTpcCluster->getId()); + // TDirectory* pClusterDir = pRowDir->mkdir( strId.c_str(), strId.c_str(), kTRUE ); + // if( pClusterDir && pClusterDir->cd() ) { + // pClusterDir->WriteObject( &pTpcCluster->_vnTpcDigits, "digits" ); + // TpcHit* pHit = (TpcHit*)pTpcCluster; + // std::vector<uint8_t> vHit((uint8_t*)pHit, (uint8_t*)pHit + sizeof(TpcHit)); + // pClusterDir->WriteObject( &vHit, "hit" ); + // } + // if( nTimeBinMin == UINT_MAX && nTimeBinMax == UINT_MAX ) { + // nTimeBinMin = pTpcCluster->getTimeBinMin(); + // nTimeBinMax = pTpcCluster->getTimeBinMax(); + // } else { + // if( nTimeBinMin > pTpcCluster->getTimeBinMin() ) + // nTimeBinMin = pTpcCluster->getTimeBinMin(); + // if( nTimeBinMax < pTpcCluster->getTimeBinMax() ) + // nTimeBinMax = pTpcCluster->getTimeBinMax(); + // } + // if( nPadMin == UINT_MAX && nPadMax == UINT_MAX ) { + // nPadMin = pTpcCluster->getPadMin(); + // nPadMax = pTpcCluster->getPadMax(); + // } else { + // if( nPadMin > pTpcCluster->getPadMin() ) + // nPadMin = pTpcCluster->getPadMin(); + // if( nPadMax < pTpcCluster->getPadMax() ) + // nPadMax = pTpcCluster->getPadMax(); + // } + // } + // #ifdef TPC_CLUSTER_H2X + // if( nPadMin != UINT_MAX && nPadMax != UINT_MAX && nTimeBinMin != UINT_MAX && nTimeBinMax != UINT_MAX ) + // { + // uint nPads = anTPC_ROW_PADS[_nRow]; + // std::string strRow = std::to_string(_nRow); + // TH2F h2f( "H2F", ("TPC Clusters row:" + strRow).c_str(), nPads, 0, nPads, nTPC_TIMEBINS, 0, + // nTPC_TIMEBINS ); h2f.GetXaxis()->SetTitle("Pads"); + // // h2f.GetXaxis()->SetRange(nPadMin, nPadMax + 1); + // h2f.GetYaxis()->SetTitle("TimeBins"); + // // h2f.GetYaxis()->SetRange(nTimeBinMin, nTimeBinMax + 1); + // h2f.GetZaxis()->SetTitle("ADC"); + // // h2f.GetZaxis()->SetRange(0, 1023); + // TH2S h2s( "H2S", ("TPC Undefined Clusters row:" + strRow).c_str(), nPads, 0, nPads, nTPC_TIMEBINS, + // 0, nTPC_TIMEBINS ); h2s.GetXaxis()->SetTitle("Pads"); h2s.GetYaxis()->SetTitle("TimeBins"); + // h2s.GetZaxis()->SetTitle("ADC"); + // for( const TpcCluster* pTpcCluster : _lstTpcClusters ) + // if( pTpcCluster->isUndefined() ) + // pTpcCluster->setBinContent(h2s); + // else + // pTpcCluster->setBinContent(h2f); + // if( h2s.GetEntries() ) + // pRowDir->WriteObject( &h2s, ("H2S row:" + strRow).c_str() ); + // if( h2f.GetEntries() ) + // pRowDir->WriteObject( &h2f, ("H2F row:" + strRow).c_str() ); + // } + // #endif + // } + // } + void PrintPadClusterMatrix(std::ostringstream &oss) + { + oss << "DrawPadClusters:RowClusters:: row: " << _nRow << ' '; + _tpcAdcStat.Print(oss); + oss << std::endl; + for (TpcPadCluster *pTpcPadCluster : _lstTpcPadClusters) { + size_t nSize = pTpcPadCluster->getSize(); + for (uint i = 0; i < nSize; i++) { + float fAdc = pTpcPadCluster->operator[](i); + oss << pTpcPadCluster->getPad() << ' ' << pTpcPadCluster->getTimeBinMin() + i << ' ' + << pTpcPadCluster->operator[](i) << std::endl; + } + } + } + void DrawPadClusters(std::ostringstream &oss) const + { + oss << "DrawPadClusters:RowClusters:: row: " << _nRow << ' '; + _tpcAdcStat.Print(oss); + oss << std::endl; + uint nTimeBinSize = _tpcAdcStat._nTimeBinMax / 10 + 1; + DrawTimeBinTopScale(oss, nTimeBinSize); + + size_t nPadSize = _viTpcPadClusters.size(); + for (uint i = 0; i < nPadSize; i++) { + std::list<TpcPadCluster *>::iterator iTpcPadCluster = _viTpcPadClusters[i]; + oss << std::setw(3) << (*iTpcPadCluster)->getPad() << ":"; + uint nPos = 0; + for (iTpcPadCluster; iTpcPadCluster != _viTpcPadClusters[i + 1] && iTpcPadCluster != _lstTpcPadClusters.end(); + iTpcPadCluster++) { + TpcPadCluster *pTpcPadCluster = *iTpcPadCluster; + uint nMin = pTpcPadCluster->getTimeBinMin(); + if (nMin < nPos) { // shift back to the 'n' positions (for combined PadClusters!!!) + long int nShift = nPos - nMin; + oss.seekp(oss.tellp() - nShift); + nPos -= nShift; + } + uint nGap = nMin - nPos; + if (nGap) { + oss << std::string(nGap, '-'); + nPos += nGap; + } + char ch = pTpcPadCluster->getCluster() + ? SZ_ALPHABET[pTpcPadCluster->getCluster()->getId() % strlen(SZ_ALPHABET)] + : '+'; + pTpcPadCluster->Draw(oss, ch); + nPos += pTpcPadCluster->getSize(); + } + oss << std::endl; + } + DrawTimeBinBottomScale(oss, nTimeBinSize); + } + // void Draw( std::ostringstream& oss ) const { + // oss << "Draw RowClusters:: row: " << _nRow << ' '; + // _AdcStat.Print(oss); + // oss << std::endl; + // uint nTimeBinSize = _AdcStat._nTimeBinMax / 10 + 1; + // DrawTimeBinTopScale(oss, nTimeBinSize); + // uint nPadMax = 0; + // for( std::list<TpcCluster*>::const_iterator i = _lstTpcClusters.begin(); i != _lstTpcClusters.end(); i++ ) { + // const TpcCluster* pTpcCluster = *i; + // pTpcCluster->Draw(oss); + // // oss << std::endl; + // if( pTpcCluster->getPadMax() > nPadMax ) + // nPadMax == pTpcCluster->getPadMax(); + // } + // if( nPadMax != anMPD_TPC_ROW_PADS[_nRow] ) // print empty last pad in the row + // oss << std::setw(3) << anMPD_TPC_ROW_PADS[_nRow] << ":---- last pad in the row ----" << std::endl; + // DrawTimeBinBottomScale(oss, nTimeBinSize); + // } + void Draw(std::ostringstream &oss, bool bAdc = false) const + { // bAdc: draw adc pad cluster histogram + oss << "Draw RowClusters:: row: " << _nRow << ' '; + _tpcAdcStat.Print(oss); + oss << std::endl; + uint nTimeBinSize = _tpcAdcStat._nTimeBinMax / 10 + 1; + DrawTimeBinTopScale(oss, nTimeBinSize); + uint nPadMax = 0; + for (std::list<TpcCluster *>::const_iterator i = _lstTpcClusters.begin(); i != _lstTpcClusters.end(); i++) { + const TpcCluster *pTpcCluster = *i; + pTpcCluster->Draw(oss, bAdc); + if (pTpcCluster->getPadMax() > nPadMax) nPadMax == pTpcCluster->getPadMax(); + } + if (nPadMax != anTPC_ROW_PADS[_nRow]) // print empty last pad in the row + oss << std::setw(3) << anTPC_ROW_PADS[_nRow] << ":---- last pad in the row ----" << std::endl; + if (!bAdc) DrawTimeBinBottomScale(oss, nTimeBinSize); + } + void Print(std::ostringstream &oss, std::vector<std::list<TpcPadCluster *>::iterator> &rviTpcPadClusters) + { + oss << "RowClusters:: row: " << _nRow << ' '; + _tpcAdcStat.Print(oss); + oss << std::endl; + if (rviTpcPadClusters.size() > 1) { + size_t nSize1 = rviTpcPadClusters.size() - 1; + for (uint i = 0; i < nSize1; i++) { + std::list<TpcPadCluster *>::iterator iTpcPadCluster = rviTpcPadClusters[i]; + std::list<TpcPadCluster *>::iterator iTpcPadClusterLast = rviTpcPadClusters[i + 1]; + oss << " ---------------- id: " << i << std::endl; + for (iTpcPadCluster; iTpcPadCluster != iTpcPadClusterLast; iTpcPadCluster++) { + oss << "\t"; + (*iTpcPadCluster)->Print(oss); + } + } + } + oss << " ----------------" << std::endl; + } + void PrintPadClusters(std::ostringstream &oss) + { + oss << "PrintPadClusters:RowClusters:: row: " << _nRow << ' '; + _tpcAdcStat.Print(oss); + oss << std::endl; + if (_viTpcPadClusters.size() > 1) { + size_t nSize1 = _viTpcPadClusters.size() - 1; + for (uint i = 0; i < nSize1; i++) { + std::list<TpcPadCluster *>::iterator iTpcPadCluster = _viTpcPadClusters[i]; + std::list<TpcPadCluster *>::iterator iTpcPadClusterLast = _viTpcPadClusters[i + 1]; + oss << " ---------------- id: " << i << std::endl; + for (iTpcPadCluster; iTpcPadCluster != iTpcPadClusterLast; iTpcPadCluster++) { + oss << "\t"; + (*iTpcPadCluster)->Print(oss); + } + } + } + oss << " ----------------" << std::endl; + } + void Print(std::ostringstream &oss) const + { + oss << "Print:RowClusters:: row: " << getRow() << ' '; + _tpcAdcStat.Print(oss); + oss << std::endl; + for (std::list<TpcCluster *>::const_iterator i = _lstTpcClusters.begin(); i != _lstTpcClusters.end(); i++) + (*i)->Print(oss); + } +}; + +} // namespace TpcClustering + +#endif // TPC_ROW_CLUSTERS_H_ \ No newline at end of file diff --git a/detectors/tpc/clusterHitFinder/fast/clustering/TpcSampaPack.h b/detectors/tpc/clusterHitFinder/fast/clustering/TpcSampaPack.h new file mode 100644 index 00000000..3b90b108 --- /dev/null +++ b/detectors/tpc/clusterHitFinder/fast/clustering/TpcSampaPack.h @@ -0,0 +1,173 @@ +// +// Joint Institute for Nuclear Research (JINR), Dubna, Russia, June-2024 +// TPC SAMPA Pack structure for MPD/NICA project +// author: Viktor Krylov JINR/LHEP +// email: kryman@jinr.ru +// +#ifndef TPC_SAMPA_PACK_H +#define TPC_SAMPA_PACK_H + +#include <cassert> +#include <arpa/inet.h> + +#include "TpcSectorGeo.h" + +#define nSAMPA_HEADER_BITS 50 +#define nSAMPA_PAYLOAD_BITS 10 +#define nSAMPA_PAYLOAD_WORDS 311 +#define nSAMPA_BUFFER_SIZE ((nSAMPA_HEADER_BITS + nSAMPA_PAYLOAD_WORDS * nSAMPA_PAYLOAD_BITS) >> 3) + +// SAMPA Header +// ------ 6 Hamming code +// - 1 Header parity bit +// - -- 3 Packet type +// | ------ ---- 10 Payload 10 bit words +// | | ---- 4 Hardware address id +// | | ----- 5 Channel address id +// | | | --- -------- -------- - Bunch-crossing counter +// | | | | | | - Data payload parity bit +// --------|--------|--------|--------|--------|--------|-- 50 bits +// 0 1 2 3 4 5 6 +// 01234567 89ABCDEF 01234567 89ABCDEF 01234567 89ABCDEF 01 +// +// Simple emulate of the SAMPA Header +// 00000000 00000000 00000000 000rrrrr rppppppp tttttttt tt +// t: Timebins +// p: Pad number +// r: Row number +// 0: ? +// 0 1 2 3 4 5 6 7 +// 0123456 7 8 9 A B C D E F +// 0808080123 4567 89AB CDEF 0123 4567 89AB CDEF 0123 4567 89AB CDEF 0123 4567 89AB CDEF 0123 4567 0 +// --oo oooo oooo ---- | | | | | | +// [header] ---- oooo oooo oo-- | | | | | +// ---- --oo oooo oooo | | | | | +// oooo oooo oo-- ---- | | | +// --oo oooo oooo ---- | | +// ---- oooo oooo oo-- | +// ---- --oo oooo oooo | + +namespace TpcClustering { + +struct TpcAdcTimebin { + std::uint16_t _nTimebin; // timebin number: [0-310] + std::uint16_t _nAdc; // ADC value > 0 + TpcAdcTimebin() : _nTimebin(0), _nAdc(0) {} + TpcAdcTimebin(std::uint16_t nTimebin, std::uint16_t nAdc) : _nTimebin(nTimebin), _nAdc(nAdc) {} +}; +struct TpcSampaPack { // header: 50 bits| + std::uint32_t _nTBD; // | 27 bits + std::uint8_t _nRow; // | 6 bits + std::uint8_t _nPad; // | 7 bits + std::uint16_t _nTimebins; // | 10 bits + uint _nSize; // number of significant ADC values in TpcAdcTimebin array + TpcAdcTimebin _TpcAdcTimebins[nSAMPA_PAYLOAD_WORDS]; + + TpcSampaPack() = default; + void HeaderDecode(const std::uint8_t *pBuf) + { + assert(pBuf); + + // _nTBD = (*pBuf << 24) + (*(pBuf + 1) << 16) + (*(pBuf + 2) << 8) + *(pBuf + 3); + // _nTBD >>= 5; + + _nRow = (*(pBuf + 3) & 0x1F) << 1 | (*(pBuf + 4) & 0x80) >> 7; + assert(_nRow < nTPC_SECTOR_ROWS); + + _nPad = *(pBuf + 4) & 0x7F; + _nTimebins = (*(pBuf + 5) << 8) + *(pBuf + 6); + _nTimebins >>= 6; + assert(_nTimebins == nSAMPA_PAYLOAD_WORDS); + } + void HeaderEncode(const std::uint8_t *pBuf) + { + assert(pBuf); + + // _nTBD = (*pBuf << 24) + (*(pBuf + 1) << 16) + (*(pBuf + 2) << 8) + *(pBuf + 3); + // _nTBD >>= 5; + + _nRow = (*(pBuf + 3) & 0x1F) << 1 | (*(pBuf + 4) & 0x80) >> 7; + assert(_nRow < nTPC_SECTOR_ROWS); + + _nPad = *(pBuf + 4) & 0x7F; + _nTimebins = (*(pBuf + 5) << 8) + *(pBuf + 6); + _nTimebins >>= 6; + assert(_nTimebins == nSAMPA_PAYLOAD_WORDS); + } + std::uint16_t DecodeADC(const std::uint8_t *pBuf, uint nTimebin) const + { + std::uint16_t nAdc = (*pBuf << 8) + *(pBuf + 1); + switch (nTimebin & 0x3) { + case 0: // --oo oooo oooo ---- + return (nAdc & 0x3FFF) >> 4; + case 1: // ---- oooo oooo oo-- + return (nAdc & 0xFFF) >> 2; + case 2: // ---- --oo oooo oooo + return nAdc & 0x3FF; + case 3: // oooo oooo oo-- ---- + return nAdc >>= 6; + } + return nAdc; + } + std::uint16_t DecodeAdc(const std::uint8_t *pBuf, uint nTimeBin) const + { + assert(pBuf); + + std::uint16_t i = nTimeBin + (nTimeBin >> 2); + const std::uint8_t *pAdcBuf = pBuf + (nSAMPA_HEADER_BITS >> 3) + i; + switch (nTimeBin & 0x3) { + case 0: // --oo oooo oooo ---- + return (ntohs(*(const uint16_t *)pAdcBuf) & 0x3FFF) >> 4; + case 1: // ---- oooo oooo oo-- + return (ntohs(*(const uint16_t *)pAdcBuf) & 0xFFF) >> 2; + case 2: // ---- --oo oooo oooo + return ntohs(*(const uint16_t *)pAdcBuf) & 0x3FF; + case 3: // oooo oooo oo-- ---- + return ntohs(*(const uint16_t *)(pAdcBuf + 1)) >> 6; + } + return 0; + } + void EncodeAdc(std::uint8_t *pBuf, std::uint16_t nTimeBin, std::uint16_t nAdc) const + { + assert(pBuf); + + std::uint16_t nAdc10 = nAdc & 0x3FF; + std::uint16_t i = nTimeBin + (nTimeBin >> 2); + std::uint8_t *pAdcBuf = pBuf + (nSAMPA_HEADER_BITS >> 3) + i; + switch (nTimeBin & 0x3) { + case 0: // --oo oooo oooo ---- + *(uint16_t *)pAdcBuf |= htons(nAdc10 << 4); + break; + case 1: // ---- oooo oooo oo-- + *(uint16_t *)pAdcBuf |= htons(nAdc10 << 2); + break; + case 2: // ---- --oo oooo oooo + *(uint16_t *)pAdcBuf |= htons(nAdc10); + break; + case 3: // oooo oooo oo-- ---- + *(uint16_t *)(pAdcBuf + 1) |= htons(nAdc10 << 6); + break; + } + } + const std::uint8_t *PayloadDecode(const std::uint8_t *pBuf) + { + const std::uint8_t *pi = pBuf + 6; // exclude SAMPA header 50 bits + _nSize = 0; + for (uint i = 0; i < _nTimebins; i++) { + std::uint16_t nAdc = DecodeADC(pi, i); + if (nAdc) // ADC > 0 + _TpcAdcTimebins[_nSize++] = TpcAdcTimebin(i, nAdc); + pi += (i & 0x3) == 2 ? 2 : 1; + } + return pi; + } + inline bool isEmpty() const { return _nSize == 0; } + void Print() const + { + std::printf("\tHeader: TBD: %d; nROW: %d; PAD: %d; TBS: %d\n", _nTBD, _nRow, _nPad, _nTimebins); + } +}; + +} // namespace TpcClustering + +#endif // SAMPA_PACK_H \ No newline at end of file diff --git a/detectors/tpc/clusterHitFinder/fast/clustering/TpcSectorGeo.h b/detectors/tpc/clusterHitFinder/fast/clustering/TpcSectorGeo.h new file mode 100644 index 00000000..5f020a00 --- /dev/null +++ b/detectors/tpc/clusterHitFinder/fast/clustering/TpcSectorGeo.h @@ -0,0 +1,62 @@ +// +// Joint Institute for Nuclear Research (JINR), Dubna, Russia, Jan-2025 +// TPC Sector Geometry definitions for MPD/NICA project +// author: Viktor Krylov JINR/LHEP +// email: kryman@jinr.ru +// +#ifndef TPC_SECTOR_GEO_H +#define TPC_SECTOR_GEO_H + +#include <cassert> +#include <algorithm> + +#define TPC_HITS_H3F +#define TPC_DIGITS_H3S +#define TPC_HITDEV_H2F +#define TPC_ROWPAD_H2I +#define TPC_CLUSTER_H2X +#define TPC_CLUSTER_H2F +#define TPC_CLUSTER_H2S + +#define nTPC_SECTOR_ROWS 53 + +#define nTPC_TIMEBINS 311 +#define nTPC_SECTOR_PADS 4112 +#define nTPC_SECTOR_SAMPA_SIZE 1624240 +#define nTPC_SECTORS 24 +#define nTPC_SECTORS2 12 +#define nTPC_ROW_MAX_PADS 124 +#define nTPC_INNER_ROWS 27 +#define nTPC_OUTER_ROWS 26 + +#define nTIMEBIN_COUNT 310 +#define fTIMEBIN_LENGTHns 100.0f + +#define fINNER_PAD_WIDTHcm 0.5f +#define fINNER_PAD_HEIGHTcm 1.2f +#define fOUTER_PAD_WIDTHcm 0.5f +#define fOUTER_PAD_HEIGHTcm 1.8f + +#define fPAD_PLANE_OFFSET 40.3f +#define fPAD_AREA_OFFSET 0.4f +#define fDRIFT_LENGTHcm 163.879f +// ADC transfer function (TF) tau param, ns, needed to compensate TF offset in drift time! +#define fTAU_SAMPA_DELAYns 160.0f +#define fDRIFT_VELOCITYcmns 0.0055f + +const uint anTPC_ROW_PADS[nTPC_SECTOR_ROWS] = { + // sectors: 24; rows: 53, pads: 1530 + 2582 = 4112 * 24 = 98,688 pads + 40, 42, 42, 44, 46, 46, 48, 48, 50, 52, 52, 54, 56, 56, + 58, 60, 60, 62, 64, 64, 66, 66, 68, 70, 70, 72, 74, // inner short pad rows: 0-26 + 76, 78, 80, 82, 82, 84, 86, 88, 90, 92, 94, 96, 98, 100, + 102, 104, 106, 108, 110, 112, 114, 116, 118, 120, 122, 124 // outer long pad rows: 27-52 +}; +const uint anTPC_ROW_PAD_SHIFT[nTPC_SECTOR_ROWS] = { + // sectors: 24; rows: 53, pads: 1530 + 2582 = 4112 * 24 = 98,688 pads + 0, 40, 82, 124, 168, 214, 260, 308, 356, 406, 458, 510, 564, 620, + 676, 734, 794, 854, 916, 980, 1044, 1110, 1176, 1244, 1314, 1384, 1456, // inner short pad rows: 0-26 + 1530, 1606, 1684, 1764, 1846, 1928, 2012, 2098, 2186, 2276, 2368, 2462, 2558, 2656, + 2756, 2858, 2962, 3068, 3176, 3286, 3398, 3512, 3628, 3746, 3866, 3988 // outer long pad rows: 27-52 +}; + +#endif // TPC_SECTOR_GEO_H \ No newline at end of file diff --git a/detectors/tpc/clusterHitFinder/fast/clustering/TpcTrack.h b/detectors/tpc/clusterHitFinder/fast/clustering/TpcTrack.h new file mode 100644 index 00000000..6c07c526 --- /dev/null +++ b/detectors/tpc/clusterHitFinder/fast/clustering/TpcTrack.h @@ -0,0 +1,209 @@ +// +// Joint Institute for Nuclear Research (JINR), Dubna, Russia, Feb-2025 +// TPC MC Track objects to read track info from root file for MPD/NICA project +// author: Viktor Krylov JINR/LHEP +// email: kryman@jinr.ru +// +#ifndef TPC_MCTRACK_H +#define TPC_MCTRACK_H + +#include <stdlib.h> +#include <iostream> +#include <stdint.h> +#include <float.h> +#include <vector> + +#include "TpcSectorGeo.h" + +extern const char *szTPC_MC_POINTS; +extern const char *szTPC_MC_TRACK; + +namespace TpcClustering { + +struct TpcDigit { + uint16_t _nRow; // row number + uint16_t _nPad; // pad number + uint16_t _nTimeBin; // timebin value + uint16_t _nAdc; // adc value + uint32_t _nTrack; // track id + TpcDigit() : _nRow(0), _nPad(0), _nTimeBin(0), _nAdc(0), _nTrack(0) {} + void Print() const + { + std::cout << "TpcDigit:: row:" << _nRow << "; pad:" << _nPad << "; timebin:" << _nTimeBin << "; adc:" << _nAdc + << "; track:" << _nTrack << std::endl; + } + friend std::ostream &operator<<(std::ostream &sout, const TpcDigit &r) + { + return sout << "{ row:" << r._nRow << "; pad:" << r._nPad << "; timebin:" << r._nTimeBin << "; adc:" << r._nAdc + << "; track:" << r._nTrack << " }"; + } +}; +struct TpcPoint { + float _fX; + float _fY; + float _fZ; + TpcPoint() : _fX(0), _fY(0), _fZ(0) {} + TpcPoint(float fX, float fY, float fZ) : _fX(fX), _fY(fY), _fZ(fZ) {} + TpcPoint(uint nRow, float fPad, float fTimeBin) { setLocal(nRow, fPad, fTimeBin); } // convert to local coords + TpcPoint(uint16_t nRow, uint16_t nPad, uint16_t nTimeBin) + { + setLocal(nRow, nPad, nTimeBin); + } // convert to local coords + TpcPoint(const TpcPoint &r) { *this = r; } + + bool isDefined() const { return _fX != 0 || _fY != 0 || _fZ != 0; } + const TpcPoint &operator=(const TpcPoint &r) + { + if (this != &r) std::memcpy(this, &r, sizeof(TpcPoint)); + return *this; + } + void setLocal(uint nRow, float fPad, float fTimeBin) + { + assert(nRow < nTPC_SECTOR_ROWS && fPad < anTPC_ROW_PADS[nRow]); + + float fX = fPad - (anTPC_ROW_PADS[nRow] >> 1); + float fY = nRow + 0.5f; // + 0.5 in the middle of the row + if (nRow > nTPC_INNER_ROWS) { + _fX = fX * fOUTER_PAD_WIDTHcm; + _fY = (fY - nTPC_INNER_ROWS) * fOUTER_PAD_HEIGHTcm + nTPC_INNER_ROWS * fINNER_PAD_HEIGHTcm; + } else { + _fX = fX * fINNER_PAD_WIDTHcm; + _fY = fY * fINNER_PAD_HEIGHTcm; + } + _fZ = fTimeBin * fTIMEBIN_LENGTHns * fDRIFT_VELOCITYcmns; + } + TpcPoint getGlobal(uint nSector) const + { + bool bSector = nSector < nTPC_SECTORS2; + // rotate 180 around Y + float fX = bSector ? -_fX : _fX; + float fY = _fY + fPAD_PLANE_OFFSET + fPAD_AREA_OFFSET; + float fZ = bSector ? fDRIFT_LENGTHcm - _fZ : _fZ - fDRIFT_LENGTHcm; + // rotate around Z + float fRad = M_PI / 180; + float fSectorPhiRad = -30 * fRad * nSector; + float fSin = sinf(fSectorPhiRad); + float fCos = cosf(fSectorPhiRad); + + return {fCos * fX - fSin * fY, fSin * fX + fCos * fY, fZ}; + } + TpcPoint getNormal(const TpcPoint &rTpcPoint1, const TpcPoint &rTpcPoint2) const + { + float fX = _fX - rTpcPoint1._fX; // BA + float fY = _fY - rTpcPoint1._fY; + float fZ = _fZ - rTpcPoint1._fZ; + + float fSx = rTpcPoint2._fX - rTpcPoint1._fX; // BC + float fSy = rTpcPoint2._fY - rTpcPoint1._fY; + float fSz = rTpcPoint2._fZ - rTpcPoint1._fZ; + float fS = fSx * fSx + fSy * fSy + fSz * fSz; + if (fS > 0.0f) { + float fT = (fSx * fX + fSy * fY + fSz * fZ) / fS; + + float fHx = fSx * fT + rTpcPoint1._fX; // perpendicular to deviation point + float fHy = fSy * fT + rTpcPoint1._fY; + float fHz = fSz * fT + rTpcPoint1._fZ; + + float fNx = _fX - fHx; // normal to hit point + float fNy = _fY - fHy; + float fNz = _fZ - fHz; + // float f = fNx * fSx + fNy * fSy + fNz * fSz; // check orthogonality + return {fNx, fNy, fNz}; + } else + return {0.0f, 0.0f, 0.0f}; + } + float getDeviation(const TpcPoint &rTpcPoint1, const TpcPoint &rTpcPoint2) const + { + TpcPoint tpcNormal = getNormal(rTpcPoint1, rTpcPoint2); + return isDefined() + ? sqrt(tpcNormal._fX * tpcNormal._fX + tpcNormal._fY * tpcNormal._fY + tpcNormal._fZ * tpcNormal._fZ) + : FLT_MAX; + } + void Print() const { std::cout << "TPC_Point:: x:" << _fX << "; y:" << _fY << "; z:" << _fZ << std::endl; } + friend std::ostream &operator<<(std::ostream &sout, const TpcPoint &r) + { + return sout << "{ x:" << r._fX << "; y:" << r._fY << "; z:" << r._fZ << " }"; + } +}; +// std::ostream& operator<<(std::ostream& sout, const McPoint& r) { +// return sout << "{ x:" << r._fX << "; y:" << r._fY << "; z:" << r._fZ << " }"; +// } +struct TpcVertex : TpcPoint { + TpcVertex() : TpcPoint() {} + TpcVertex(float fX, float fY, float fZ) : TpcPoint(fX, fY, fZ) {} + TpcVertex(const TpcPoint &r) : TpcPoint(r) {} + void Print() const { std::cout << "TpcVertex:: x:" << _fX << "; y:" << _fY << "; z:" << _fZ << std::endl; } +}; +// struct RootTrackInfo { int id;int pdg;float px;float py;float pz; }; // record format in root file +struct TpcTrackInfo { + int _nMotherId; // -1: primary track + int _nPDG; + float _fPx; + float _fPy; + float _fPz; + + TpcTrackInfo() : _nMotherId(0), _nPDG(0), _fPx(0), _fPy(0), _fPz(0) {} + TpcTrackInfo(int nMotherId, int nPDG, float fPx, float fPy, float fPz) + : _nMotherId(nMotherId), _nPDG(nPDG), _fPx(fPx), _fPy(fPy), _fPz(fPz) + { + } + TpcTrackInfo(const TpcTrackInfo &r) + : _nMotherId(r._nMotherId), _nPDG(r._nPDG), _fPx(r._fPx), _fPy(r._fPy), _fPz(r._fPz) + { + } + friend std::ostream &operator<<(std::ostream &sout, const TpcTrackInfo &r) + { + return sout << "TpcTrackInfo:: mother_id:" << r._nMotherId << "; pdg:" << r._nPDG << "; { px:" << r._fPx + << "; py:" << r._fPy << "; pz:" << r._fPz << " }" << std::endl; + } + void Print() const + { + std::cout << *this; + // std::cout << "TpcTrackInfo:: mother_id:" << _nMotherId << "; pdg:" << _nPDG << "; { px:" << _fPx << "; py:" << + // _fPy << "; pz:" << _fPz << " }" << std::endl; + } +}; +struct TpcTrack : TpcTrackInfo { + int _nTrackId; + std::vector<TpcPoint> *_pvTpcPoints; + TpcTrack() : TpcTrackInfo(), _nTrackId(0) {} + TpcTrack(int nTrackId, const TpcTrackInfo &r) : TpcTrackInfo(r), _nTrackId(nTrackId) {} + virtual ~TpcTrack() + { + _pvTpcPoints->clear(); + delete _pvTpcPoints; + } + friend std::ostream &operator<<(std::ostream &sout, const TpcTrack &r) + { + sout << "TpcTrack:: id:" << r._nTrackId << ' ' << static_cast<const TpcTrackInfo &>(r); + sout << std::string(2, ' ') << "points:" << std::endl; + for (const TpcPoint &r : *r._pvTpcPoints) sout << std::string(4, ' ') << r << std::endl; + return sout; + } + int getTrackId() const { return _nTrackId; } + void Store(TDirectory *pDir) const + { + TDirectory *pTrackDir = pDir->mkdir(std::to_string(_nTrackId).c_str(), "Tpc Event McTrack", kTRUE); + if (pTrackDir && pTrackDir->cd() && _pvTpcPoints) { + const TpcPoint *pTpcPoint = _pvTpcPoints->data(); + std::vector<uint8_t> vPoints((uint8_t *)pTpcPoint, + (uint8_t *)pTpcPoint + (sizeof(TpcPoint) * _pvTpcPoints->size())); + pTrackDir->WriteObject(&vPoints, ::szTPC_MC_POINTS); + + const TpcTrackInfo *pTrackInfo = this; + std::vector<uint8_t> vTrack((uint8_t *)pTrackInfo, (uint8_t *)pTrackInfo + sizeof(TpcTrackInfo)); + pTrackDir->WriteObject(&vTrack, ::szTPC_MC_TRACK); + } + } + void Print() const + { + std::cout << "TpcTrack:: id:" << _nTrackId << ' '; + TpcTrackInfo::Print(); + std::cout << std::string(2, ' ') << "points:" << std::endl; + for (const TpcPoint &r : *_pvTpcPoints) std::cout << std::string(4, ' ') << r << std::endl; + } +}; + +} // namespace TpcClustering + +#endif // TPC_MCTRACK_H \ No newline at end of file diff --git a/detectors/tpc/clusterHitFinder/fast/utility.cxx b/detectors/tpc/clusterHitFinder/fast/utility.cxx new file mode 100644 index 00000000..0b3a2c18 --- /dev/null +++ b/detectors/tpc/clusterHitFinder/fast/utility.cxx @@ -0,0 +1,92 @@ +// +// Joint Institute for Nuclear Research (JINR), Dubna, Russia, Nov-2024 +// Useful utilities for MPD/NICA project +// author: Viktor Krylov JINR/LHEP +// email: kryman@jinr.ru +// +#include <cstring> +#include <iostream> +#include <iomanip> +#include <unistd.h> +#include <string> +#include <termios.h> + +#include "utility.h" + +int getKey() +{ + // set stdin to unbuffered in order to read directly from the keyboard. + struct termios term; + tcgetattr(STDIN_FILENO, &term); + // ECHO: Enable echo + // ICANON: Canonical input (erase and kill processing) + term.c_lflag &= ~(ECHO | ICANON); + tcsetattr(STDIN_FILENO, TCSANOW, &term); + int key = std::cin.get(); // read keyboard + // restore c_lflag + term.c_lflag |= ECHO | ICANON; + tcsetattr(STDIN_FILENO, TCSANOW, &term); + return key; +} +// std::string szTrim( const char* sz, char ch ) { +// if( !strlen(sz) ) return std::string(sz, 0); +// const char* p = sz; +// while( *p == ch ) p++; +// size_t n = strlen(p); +// if( n ) while( *(p + n - 1) == ch ) n--; +// return std::string(p, n); +// } +std::string szTrim(char const *sz, char ch) +{ + if (!strlen(sz)) return std::string(sz, 0); + const char *p = sz; + while (*p == ch) p++; + size_t n = strlen(p); + if (n) + while (*(p + n - 1) == ch) n--; + return std::string(p, n); +} +std::string GetCCver() +{ + switch (__cplusplus) { + case 202101L: + case 202302L: return "'C++23'"; + case 202002L: return "'C++20'"; + case 201703L: return "'C++17'"; + case 201402L: return "'C++14'"; + case 201103L: return "'C++11'"; + case 199711L: return "'C++98'"; + default: return '\'' + __cplusplus + "'"; + } +} +// https://favtutor.com/blogs/wildcard-pattern-matching +// bool isMatch(const std::string& pattern, const std::string& target) { +// if (pattern.empty() && target.empty()) +// return true; +// if (pattern.empty()) +// return false; +// if (target.empty()) { +// if (pattern[0] == '*') +// return isMatch(pattern.substr(1), target); +// return false; +// } + +// if (pattern[0] == target[0] || pattern[0] == '?') +// return isMatch(pattern.substr(1), target.substr(1)); +// if (pattern[0] == '*') +// return isMatch(pattern.substr(1), target) || isMatch(pattern, target.substr(1)); + +// return false; +// } +bool isMatch(const char *szPattern, const char *szTarget) +{ + size_t nPattern = strlen(szPattern); + if (!nPattern) // all filenames are matched + return true; + size_t nTarget = strlen(szTarget); + if (!nTarget) // empty target cannot be matched + return false; + if (*szPattern == *szTarget || *szPattern == '?') return isMatch(szPattern + 1, szTarget + 1); + if (*szPattern == '*') return isMatch(szPattern + 1, szTarget) || isMatch(szPattern, szTarget + 1); + return false; +} diff --git a/detectors/tpc/clusterHitFinder/fast/utility.h b/detectors/tpc/clusterHitFinder/fast/utility.h new file mode 100644 index 00000000..e2eddd4c --- /dev/null +++ b/detectors/tpc/clusterHitFinder/fast/utility.h @@ -0,0 +1,37 @@ +// +// Joint Institute for Nuclear Research (JINR), Dubna, Russia, Nov-2024 +// Useful utilities for MPD/NICA project +// author: Viktor Krylov JINR/LHEP +// email: kryman@jinr.ru +// +#ifndef SAMPA_UTILITY_H +#define SAMPA_UTILITY_H + +#include <sstream> +#include <chrono> +#include <math.h> +#include <string> + +int getKey(); +// std::string szTrim( const char* sz, char ch ); +std::string szTrim(char const *sz, char ch); +std::string GetCCver(); +bool isMatch(const char *szPattern, const char *szTarget); +std::string GetCCver(); + +class CpuTimer { + std::chrono::high_resolution_clock::time_point _tpStart; + std::chrono::high_resolution_clock::time_point _tpStop; + +public: + CpuTimer() {} + ~CpuTimer() {} + void Start() { _tpStart = std::chrono::high_resolution_clock::now(); } + void Stop() { _tpStop = std::chrono::high_resolution_clock::now(); } + const std::chrono::duration<float, std::milli> GetDuration() const + { + return std::chrono::duration_cast<std::chrono::duration<float, std::milli>>(_tpStop - _tpStart); + } + void PrintDuration(const char *szMsg) const { std::printf("%s%f\n", szMsg, GetDuration().count()); } +}; +#endif // SAMPA_UTILITY_H -- GitLab From 137ecc2acc0cddb34e89f21e12ed180c555087ff Mon Sep 17 00:00:00 2001 From: Slavomir Hnatic <hnatics@jinr.ru> Date: Tue, 1 Apr 2025 22:08:26 +0300 Subject: [PATCH 2/2] Fast clusterhitfinder: rename TpcPoint to TpcPointObj to resolve ambiguity --- .../tpc/clusterHitFinder/fast/TpcEventObj.cxx | 2 +- .../clusterHitFinder/fast/TpcSectorObj.cxx | 32 ++++++------- .../fast/clustering/TpcAdcHit.h | 4 +- .../fast/clustering/TpcRowClusters.h | 2 +- .../fast/clustering/TpcTrack.h | 48 +++++++++---------- 5 files changed, 44 insertions(+), 44 deletions(-) diff --git a/detectors/tpc/clusterHitFinder/fast/TpcEventObj.cxx b/detectors/tpc/clusterHitFinder/fast/TpcEventObj.cxx index 035482c9..8842fefa 100644 --- a/detectors/tpc/clusterHitFinder/fast/TpcEventObj.cxx +++ b/detectors/tpc/clusterHitFinder/fast/TpcEventObj.cxx @@ -67,7 +67,7 @@ void TpcEventObj::loadTrack(TDirectory *pDir) pDir->GetObject(pPointsKey->GetName(), pvTpcPoints); if (pvTpcPoints) { TpcTrack *pTpcTrack = new TpcTrack(nTrackID, TpcTrackInfo(*(const TpcTrackInfo *)pv->data())); - pTpcTrack->_pvTpcPoints = (std::vector<TpcPoint> *)pvTpcPoints; + pTpcTrack->_pvTpcPoints = (std::vector<TpcPointObj> *)pvTpcPoints; _vpTpcTracks.push_back(pTpcTrack); } } diff --git a/detectors/tpc/clusterHitFinder/fast/TpcSectorObj.cxx b/detectors/tpc/clusterHitFinder/fast/TpcSectorObj.cxx index 0b739818..4537b298 100644 --- a/detectors/tpc/clusterHitFinder/fast/TpcSectorObj.cxx +++ b/detectors/tpc/clusterHitFinder/fast/TpcSectorObj.cxx @@ -231,8 +231,8 @@ void TpcSectorObj::StoreH1FxyzDigits(TDirectory *pDir, const std::vector<TpcTrac h1fz.SetFillColor(45); for (const TpcDigit &rTpcDigit : *_pvTpcDigits) { // A - TpcPoint tpcLocal(rTpcDigit._nRow, rTpcDigit._nPad, rTpcDigit._nTimeBin); - TpcPoint tpcGlobal(tpcLocal.getGlobal(_nSector)); + TpcPointObj tpcLocal(rTpcDigit._nRow, rTpcDigit._nPad, rTpcDigit._nTimeBin); + TpcPointObj tpcGlobal(tpcLocal.getGlobal(_nSector)); uint nTrackId = rTpcDigit._nTrack; std::vector<TpcTrack *>::const_iterator iTpcTrack = std::find_if( rvpTpcTracks.begin(), rvpTpcTracks.end(), [nTrackId](const TpcTrack *p) { return p->_nTrackId == nTrackId; }); @@ -240,14 +240,14 @@ void TpcSectorObj::StoreH1FxyzDigits(TDirectory *pDir, const std::vector<TpcTrac TpcTrack *pTpcTrack = *iTpcTrack; size_t nSize = pTpcTrack->_pvTpcPoints->size(); if (nSize > 1) { - float fDev = FLT_MAX; // deviation - TpcPoint tpcNormal; - uint n = 0; + float fDev = FLT_MAX; // deviation + TpcPointObj tpcNormal; + uint n = 0; for (uint i = 1; i < nSize; i++) { - const TpcPoint &rTpcPoint1 = pTpcTrack->_pvTpcPoints->operator[](i - 1); // B - const TpcPoint &rTpcPoint2 = pTpcTrack->_pvTpcPoints->operator[](i); // C - TpcPoint tpcN = tpcGlobal.getNormal(rTpcPoint1, rTpcPoint2); - float f = + const TpcPointObj &rTpcPoint1 = pTpcTrack->_pvTpcPoints->operator[](i - 1); // B + const TpcPointObj &rTpcPoint2 = pTpcTrack->_pvTpcPoints->operator[](i); // C + TpcPointObj tpcN = tpcGlobal.getNormal(rTpcPoint1, rTpcPoint2); + float f = tpcN.isDefined() ? sqrt(tpcN._fX * tpcN._fX + tpcN._fY * tpcN._fY + tpcN._fZ * tpcN._fZ) : FLT_MAX; if (f != FLT_MAX) { if (fDev > f) { @@ -294,7 +294,7 @@ void TpcSectorObj::StoreH1FxyzHits(TDirectory *pDir, const std::vector<TpcTrack for (const TpcAdcHit &rTpcAdcHit : _vTpcAdcHits) { // A uint nRow = rTpcAdcHit._nRow; uint nClusterId = rTpcAdcHit._nClusterId; - TpcPoint tpcHitPoint = {rTpcAdcHit._fX, rTpcAdcHit._fY, rTpcAdcHit._fZ}; + TpcPointObj tpcHitPoint = {rTpcAdcHit._fX, rTpcAdcHit._fY, rTpcAdcHit._fZ}; const TpcCluster *pTpcCluster = _pTpcClusters->findTpcCluster(rTpcAdcHit._nRow, rTpcAdcHit._nClusterId); if (pTpcCluster && _pvTpcDigits) { for (uint nTpcDigitId : pTpcCluster->_vnTpcDigits) { @@ -306,13 +306,13 @@ void TpcSectorObj::StoreH1FxyzHits(TDirectory *pDir, const std::vector<TpcTrack TpcTrack *pTpcTrack = *iTpcTrack; size_t nSize = pTpcTrack->_pvTpcPoints->size(); if (nSize > 1) { - float fDev = FLT_MAX; // min hit deviation - TpcPoint tpcNormal; - uint n = 0; + float fDev = FLT_MAX; // min hit deviation + TpcPointObj tpcNormal; + uint n = 0; for (uint i = 1; i < nSize; i++) { - const TpcPoint &rTpcPoint1 = pTpcTrack->_pvTpcPoints->operator[](i - 1); // B - const TpcPoint &rTpcPoint2 = pTpcTrack->_pvTpcPoints->operator[](i); // C - TpcPoint tpcN = tpcHitPoint.getNormal(rTpcPoint1, rTpcPoint2); + const TpcPointObj &rTpcPoint1 = pTpcTrack->_pvTpcPoints->operator[](i - 1); // B + const TpcPointObj &rTpcPoint2 = pTpcTrack->_pvTpcPoints->operator[](i); // C + TpcPointObj tpcN = tpcHitPoint.getNormal(rTpcPoint1, rTpcPoint2); float f = tpcN.isDefined() ? sqrt(tpcN._fX * tpcN._fX + tpcN._fY * tpcN._fY + tpcN._fZ * tpcN._fZ) : FLT_MAX; if (f != FLT_MAX) { diff --git a/detectors/tpc/clusterHitFinder/fast/clustering/TpcAdcHit.h b/detectors/tpc/clusterHitFinder/fast/clustering/TpcAdcHit.h index bcff1c1b..a74b6d37 100644 --- a/detectors/tpc/clusterHitFinder/fast/clustering/TpcAdcHit.h +++ b/detectors/tpc/clusterHitFinder/fast/clustering/TpcAdcHit.h @@ -36,7 +36,7 @@ struct TpcAdcHit { { return this == &r ? *this : *(TpcAdcHit *)memcpy(this, (const void *)&r, sizeof(TpcAdcHit)); } - void operator=(const TpcPoint &r) + void operator=(const TpcPointObj &r) { _fX = r._fX; _fY = r._fY; @@ -56,7 +56,7 @@ struct TpcAdcHit { } void setGlobal(uint nSector, uint nRow, float fPad, float fTimeBin) { - this->operator=(TpcPoint(nRow, fPad, fTimeBin).getGlobal(nSector)); + this->operator=(TpcPointObj(nRow, fPad, fTimeBin).getGlobal(nSector)); } void shiftSampaDelay(uint nSector, float fDelay) { _fZ += nSector < 12 ? fDelay : -fDelay; } friend std::ostream &operator<<(std::ostream &os, const TpcAdcHit &r) diff --git a/detectors/tpc/clusterHitFinder/fast/clustering/TpcRowClusters.h b/detectors/tpc/clusterHitFinder/fast/clustering/TpcRowClusters.h index 30def841..997dfb5b 100644 --- a/detectors/tpc/clusterHitFinder/fast/clustering/TpcRowClusters.h +++ b/detectors/tpc/clusterHitFinder/fast/clustering/TpcRowClusters.h @@ -90,7 +90,7 @@ public: if (pTpcCluster->isUndefined()) // remove all ill-defined tpc clusters! continue; TpcAdcHit tpcAdcHit(_nRow, pTpcCluster->getId(), pTpcCluster->getAdcSum()); - tpcAdcHit = TpcPoint(_nRow, pTpcCluster->getPad(), pTpcCluster->getTimeBin()) + tpcAdcHit = TpcPointObj(_nRow, pTpcCluster->getPad(), pTpcCluster->getTimeBin()) .getGlobal(nSector); // convert to global coords tpcAdcHit.shiftSampaDelay(nSector, fTAU_SAMPA_DELAYns * fDRIFT_VELOCITYcmns); rvTpcAdcHits.push_back(tpcAdcHit); diff --git a/detectors/tpc/clusterHitFinder/fast/clustering/TpcTrack.h b/detectors/tpc/clusterHitFinder/fast/clustering/TpcTrack.h index 6c07c526..dbd7702b 100644 --- a/detectors/tpc/clusterHitFinder/fast/clustering/TpcTrack.h +++ b/detectors/tpc/clusterHitFinder/fast/clustering/TpcTrack.h @@ -38,23 +38,23 @@ struct TpcDigit { << "; track:" << r._nTrack << " }"; } }; -struct TpcPoint { +struct TpcPointObj { float _fX; float _fY; float _fZ; - TpcPoint() : _fX(0), _fY(0), _fZ(0) {} - TpcPoint(float fX, float fY, float fZ) : _fX(fX), _fY(fY), _fZ(fZ) {} - TpcPoint(uint nRow, float fPad, float fTimeBin) { setLocal(nRow, fPad, fTimeBin); } // convert to local coords - TpcPoint(uint16_t nRow, uint16_t nPad, uint16_t nTimeBin) + TpcPointObj() : _fX(0), _fY(0), _fZ(0) {} + TpcPointObj(float fX, float fY, float fZ) : _fX(fX), _fY(fY), _fZ(fZ) {} + TpcPointObj(uint nRow, float fPad, float fTimeBin) { setLocal(nRow, fPad, fTimeBin); } // convert to local coords + TpcPointObj(uint16_t nRow, uint16_t nPad, uint16_t nTimeBin) { setLocal(nRow, nPad, nTimeBin); } // convert to local coords - TpcPoint(const TpcPoint &r) { *this = r; } + TpcPointObj(const TpcPointObj &r) { *this = r; } - bool isDefined() const { return _fX != 0 || _fY != 0 || _fZ != 0; } - const TpcPoint &operator=(const TpcPoint &r) + bool isDefined() const { return _fX != 0 || _fY != 0 || _fZ != 0; } + const TpcPointObj &operator=(const TpcPointObj &r) { - if (this != &r) std::memcpy(this, &r, sizeof(TpcPoint)); + if (this != &r) std::memcpy(this, &r, sizeof(TpcPointObj)); return *this; } void setLocal(uint nRow, float fPad, float fTimeBin) @@ -72,7 +72,7 @@ struct TpcPoint { } _fZ = fTimeBin * fTIMEBIN_LENGTHns * fDRIFT_VELOCITYcmns; } - TpcPoint getGlobal(uint nSector) const + TpcPointObj getGlobal(uint nSector) const { bool bSector = nSector < nTPC_SECTORS2; // rotate 180 around Y @@ -87,7 +87,7 @@ struct TpcPoint { return {fCos * fX - fSin * fY, fSin * fX + fCos * fY, fZ}; } - TpcPoint getNormal(const TpcPoint &rTpcPoint1, const TpcPoint &rTpcPoint2) const + TpcPointObj getNormal(const TpcPointObj &rTpcPoint1, const TpcPointObj &rTpcPoint2) const { float fX = _fX - rTpcPoint1._fX; // BA float fY = _fY - rTpcPoint1._fY; @@ -112,15 +112,15 @@ struct TpcPoint { } else return {0.0f, 0.0f, 0.0f}; } - float getDeviation(const TpcPoint &rTpcPoint1, const TpcPoint &rTpcPoint2) const + float getDeviation(const TpcPointObj &rTpcPoint1, const TpcPointObj &rTpcPoint2) const { - TpcPoint tpcNormal = getNormal(rTpcPoint1, rTpcPoint2); + TpcPointObj tpcNormal = getNormal(rTpcPoint1, rTpcPoint2); return isDefined() ? sqrt(tpcNormal._fX * tpcNormal._fX + tpcNormal._fY * tpcNormal._fY + tpcNormal._fZ * tpcNormal._fZ) : FLT_MAX; } void Print() const { std::cout << "TPC_Point:: x:" << _fX << "; y:" << _fY << "; z:" << _fZ << std::endl; } - friend std::ostream &operator<<(std::ostream &sout, const TpcPoint &r) + friend std::ostream &operator<<(std::ostream &sout, const TpcPointObj &r) { return sout << "{ x:" << r._fX << "; y:" << r._fY << "; z:" << r._fZ << " }"; } @@ -128,10 +128,10 @@ struct TpcPoint { // std::ostream& operator<<(std::ostream& sout, const McPoint& r) { // return sout << "{ x:" << r._fX << "; y:" << r._fY << "; z:" << r._fZ << " }"; // } -struct TpcVertex : TpcPoint { - TpcVertex() : TpcPoint() {} - TpcVertex(float fX, float fY, float fZ) : TpcPoint(fX, fY, fZ) {} - TpcVertex(const TpcPoint &r) : TpcPoint(r) {} +struct TpcVertex : TpcPointObj { + TpcVertex() : TpcPointObj() {} + TpcVertex(float fX, float fY, float fZ) : TpcPointObj(fX, fY, fZ) {} + TpcVertex(const TpcPointObj &r) : TpcPointObj(r) {} void Print() const { std::cout << "TpcVertex:: x:" << _fX << "; y:" << _fY << "; z:" << _fZ << std::endl; } }; // struct RootTrackInfo { int id;int pdg;float px;float py;float pz; }; // record format in root file @@ -164,8 +164,8 @@ struct TpcTrackInfo { } }; struct TpcTrack : TpcTrackInfo { - int _nTrackId; - std::vector<TpcPoint> *_pvTpcPoints; + int _nTrackId; + std::vector<TpcPointObj> *_pvTpcPoints; TpcTrack() : TpcTrackInfo(), _nTrackId(0) {} TpcTrack(int nTrackId, const TpcTrackInfo &r) : TpcTrackInfo(r), _nTrackId(nTrackId) {} virtual ~TpcTrack() @@ -177,7 +177,7 @@ struct TpcTrack : TpcTrackInfo { { sout << "TpcTrack:: id:" << r._nTrackId << ' ' << static_cast<const TpcTrackInfo &>(r); sout << std::string(2, ' ') << "points:" << std::endl; - for (const TpcPoint &r : *r._pvTpcPoints) sout << std::string(4, ' ') << r << std::endl; + for (const TpcPointObj &r : *r._pvTpcPoints) sout << std::string(4, ' ') << r << std::endl; return sout; } int getTrackId() const { return _nTrackId; } @@ -185,9 +185,9 @@ struct TpcTrack : TpcTrackInfo { { TDirectory *pTrackDir = pDir->mkdir(std::to_string(_nTrackId).c_str(), "Tpc Event McTrack", kTRUE); if (pTrackDir && pTrackDir->cd() && _pvTpcPoints) { - const TpcPoint *pTpcPoint = _pvTpcPoints->data(); + const TpcPointObj *pTpcPoint = _pvTpcPoints->data(); std::vector<uint8_t> vPoints((uint8_t *)pTpcPoint, - (uint8_t *)pTpcPoint + (sizeof(TpcPoint) * _pvTpcPoints->size())); + (uint8_t *)pTpcPoint + (sizeof(TpcPointObj) * _pvTpcPoints->size())); pTrackDir->WriteObject(&vPoints, ::szTPC_MC_POINTS); const TpcTrackInfo *pTrackInfo = this; @@ -200,7 +200,7 @@ struct TpcTrack : TpcTrackInfo { std::cout << "TpcTrack:: id:" << _nTrackId << ' '; TpcTrackInfo::Print(); std::cout << std::string(2, ' ') << "points:" << std::endl; - for (const TpcPoint &r : *_pvTpcPoints) std::cout << std::string(4, ' ') << r << std::endl; + for (const TpcPointObj &r : *_pvTpcPoints) std::cout << std::string(4, ' ') << r << std::endl; } }; -- GitLab