diff --git a/Detectors/Upgrades/ALICE3/IOTOF/DataFormatsIOTOF/CMakeLists.txt b/Detectors/Upgrades/ALICE3/IOTOF/DataFormatsIOTOF/CMakeLists.txt index 9e075aabb2cc0..acdc927a6612b 100644 --- a/Detectors/Upgrades/ALICE3/IOTOF/DataFormatsIOTOF/CMakeLists.txt +++ b/Detectors/Upgrades/ALICE3/IOTOF/DataFormatsIOTOF/CMakeLists.txt @@ -13,7 +13,9 @@ o2_add_library(DataFormatsIOTOF SOURCES src/Digit.cxx # SOURCES src/MCLabel.cxx SOURCES src/Cluster.cxx - PUBLIC_LINK_LIBRARIES O2::DataFormatsITSMFT) + PUBLIC_LINK_LIBRARIES O2::DataFormatsITSMFT + O2::IOTOFBase + O2::FrameworkLogger) o2_target_root_dictionary(DataFormatsIOTOF HEADERS include/DataFormatsIOTOF/Digit.h diff --git a/Detectors/Upgrades/ALICE3/IOTOF/DataFormatsIOTOF/include/DataFormatsIOTOF/Cluster.h b/Detectors/Upgrades/ALICE3/IOTOF/DataFormatsIOTOF/include/DataFormatsIOTOF/Cluster.h index ad789c649c785..21028d21c9cde 100644 --- a/Detectors/Upgrades/ALICE3/IOTOF/DataFormatsIOTOF/include/DataFormatsIOTOF/Cluster.h +++ b/Detectors/Upgrades/ALICE3/IOTOF/DataFormatsIOTOF/include/DataFormatsIOTOF/Cluster.h @@ -9,28 +9,172 @@ // granted to it by virtue of its status as an Intergovernmental Organization // or submit itself to any jurisdiction. +/// \file Cluster.h +/// \brief Definition of the IOTOF cluster #ifndef ALICEO2_DATAFORMATSIOTOF_CLUSTER_H #define ALICEO2_DATAFORMATSIOTOF_CLUSTER_H -#include #include #include +#include +#include + +#include "Framework/Logger.h" + +namespace o2 +{ +namespace iotof +{ + +/// Compact encoding for ALICE3 IOTOF cluster parameters inside a single 64-bit word. +struct ClusterInfo { + // Bit widths (Total: 52 bits out of 64) + static constexpr int NBitsRow = 9; + static constexpr int NBitsCol = 8; + static constexpr int NBitsRowSpan = 4; + static constexpr int NBitsColSpan = 4; + static constexpr int NBitsPattern = 16; + static constexpr int NBitsTopology = 11; + + // Bit offsets (ordered logically from LSB to MSB) + static constexpr int ShiftRow = 0; + static constexpr int ShiftCol = ShiftRow + NBitsRow; // 9 + static constexpr int ShiftRowSpan = ShiftCol + NBitsCol; // 17 + static constexpr int ShiftColSpan = ShiftRowSpan + NBitsRowSpan; // 21 + static constexpr int ShiftPattern = ShiftColSpan + NBitsColSpan; // 25 + static constexpr int ShiftTopology = ShiftPattern + NBitsPattern; // 41 + + // Bit masks + static constexpr uint64_t MaskRow = (1ULL << NBitsRow) - 1; + static constexpr uint64_t MaskCol = (1ULL << NBitsCol) - 1; + static constexpr uint64_t MaskRowSpan = (1ULL << NBitsRowSpan) - 1; + static constexpr uint64_t MaskColSpan = (1ULL << NBitsColSpan) - 1; + static constexpr uint64_t MaskPattern = (1ULL << NBitsPattern) - 1; + static constexpr uint64_t MaskTopology = (1ULL << NBitsTopology) - 1; + + uint64_t data{0}; + + // Constructors + constexpr ClusterInfo() = default; + constexpr ClusterInfo(uint64_t d) : data(d) {} + + // Static packer + static constexpr uint64_t pack(uint32_t row, uint32_t col, uint32_t rowSpan, + uint32_t colSpan, uint32_t pattern, uint32_t topology) { + return ((static_cast(row) & MaskRow) << ShiftRow) | + ((static_cast(col) & MaskCol) << ShiftCol) | + ((static_cast(rowSpan) & MaskRowSpan) << ShiftRowSpan) | + ((static_cast(colSpan) & MaskColSpan) << ShiftColSpan) | + ((static_cast(pattern) & MaskPattern) << ShiftPattern) | + ((static_cast(topology) & MaskTopology) << ShiftTopology); + } -namespace o2::iotof + // Getters + constexpr uint32_t getRow() const { return (data >> ShiftRow) & MaskRow; } + constexpr uint32_t getCol() const { return (data >> ShiftCol) & MaskCol; } + constexpr uint32_t getRowSpan() const { return (data >> ShiftRowSpan) & MaskRowSpan; } + constexpr uint32_t getColSpan() const { return (data >> ShiftColSpan) & MaskColSpan; } + constexpr uint32_t getPattern() const { return (data >> ShiftPattern) & MaskPattern; } + constexpr uint32_t getTopology() const { return (data >> ShiftTopology) & MaskTopology; } + + // Setters + constexpr void setRow(uint32_t r) { + data = (data & ~(MaskRow << ShiftRow)) | ((static_cast(r) & MaskRow) << ShiftRow); + } + constexpr void setCol(uint32_t c) { + data = (data & ~(MaskCol << ShiftCol)) | ((static_cast(c) & MaskCol) << ShiftCol); + } + constexpr void setRowSpan(uint32_t rs) { + data = (data & ~(MaskRowSpan << ShiftRowSpan)) | ((static_cast(rs) & MaskRowSpan) << ShiftRowSpan); + } + constexpr void setColSpan(uint32_t cs) { + data = (data & ~(MaskColSpan << ShiftColSpan)) | ((static_cast(cs) & MaskColSpan) << ShiftColSpan); + } + constexpr void setPattern(uint32_t p) { + data = (data & ~(MaskPattern << ShiftPattern)) | ((static_cast(p) & MaskPattern) << ShiftPattern); + } + constexpr void setTopology(uint32_t t) { + data = (data & ~(MaskTopology << ShiftTopology)) | ((static_cast(t) & MaskTopology) << ShiftTopology); + } + + ClassDefNV(ClusterInfo, 1); +}; + +class Cluster { + public: + static constexpr uint16_t InvalidPatternID = static_cast(ClusterInfo::MaskPattern); + + Cluster() = default; + Cluster(UShort_t row, UShort_t col, UShort_t rowSpan, UShort_t colSpan, UShort_t patt, UShort_t topo, UShort_t chipID = 0, time_t time = 0.0f) + : mChipID(chipID), mTime(time) + { + mClusterInfo.data = ClusterInfo::pack(row, col, rowSpan, colSpan, patt, topo); + } + + void set(UShort_t row, UShort_t col, UShort_t rowSpan, UShort_t colSpan, UShort_t patt, UShort_t topo, UShort_t chipID, time_t time) + { + mClusterInfo.data = ClusterInfo::pack(row, col, rowSpan, colSpan, patt, topo); + mChipID = chipID; + mTime = time; + } -struct Cluster { - uint16_t chipID = 0; - uint16_t row = 0; - uint16_t col = 0; - uint16_t size = 1; - double time = 0.0; + // Unpack Getters + uint32_t getRow() const { return mClusterInfo.getRow(); } + uint32_t getCol() const { return mClusterInfo.getCol(); } + uint32_t getRowSpan() const { return mClusterInfo.getRowSpan(); } + uint32_t getColSpan() const { return mClusterInfo.getColSpan(); } + uint32_t getPattern() const { return mClusterInfo.getPattern(); } + uint32_t getTopology() const { return mClusterInfo.getTopology(); } + int getSize() const { + // Count the number of set bits in the pattern to determine the size of the cluster + uint32_t pattern = getPattern(); + int size = 0; + while (pattern) { + size += pattern & 1; + pattern >>= 1; + } + return size; + } + // BaseCluster / Interface Compatibility Getters + uint32_t getChipID() const { return mChipID; } + uint32_t getSensorID() const { return mChipID; } + time_t getTime() const { return mTime; } + uint64_t getPackedData() const { return mClusterInfo.data; } + + // Setters + void setRow(UShort_t r) { mClusterInfo.setRow(r); } + void setCol(UShort_t c) { mClusterInfo.setCol(c); } + void setRowSpan(UShort_t rs) { mClusterInfo.setRowSpan(rs); } + void setColSpan(UShort_t cs) { mClusterInfo.setColSpan(cs); } + void setPatternID(UShort_t p) { mClusterInfo.setPattern(p); } + void setTopology(UShort_t t) { mClusterInfo.setTopology(t); } + void setChipID(UShort_t c) { mChipID = c; } + void setTime(time_t t) { mTime = t; } + + // Operators & Debugging + bool operator==(const Cluster& cl) const + { + return mClusterInfo.data == cl.mClusterInfo.data && mChipID == cl.mChipID && mTime == cl.mTime; + } + + void print() const; std::string asString() const; - ClassDefNV(Cluster, 1); + private: + ClusterInfo mClusterInfo{}; ///< 64-bit packed structure containing geometry/topology + UShort_t mChipID{0}; ///< Chip / Sensor ID + float mTime{0.0f}; ///< Hit timing information + + void sanityCheck(); + + ClassDefNV(Cluster, 2); }; -} // namespace o2::iotof +} // namespace iotof +} // namespace o2 + +std::ostream& operator<<(std::ostream& stream, const o2::iotof::Cluster& cl); -#endif +#endif /* ALICEO2_DATAFORMATSIOTOF_CLUSTER_H */ diff --git a/Detectors/Upgrades/ALICE3/IOTOF/DataFormatsIOTOF/src/Cluster.cxx b/Detectors/Upgrades/ALICE3/IOTOF/DataFormatsIOTOF/src/Cluster.cxx index 6b5a4948900e7..22735a9225c19 100644 --- a/Detectors/Upgrades/ALICE3/IOTOF/DataFormatsIOTOF/src/Cluster.cxx +++ b/Detectors/Upgrades/ALICE3/IOTOF/DataFormatsIOTOF/src/Cluster.cxx @@ -9,19 +9,65 @@ // granted to it by virtue of its status as an Intergovernmental Organization // or submit itself to any jurisdiction. +/// \file Cluster.cxx +/// \brief Implementation of the IOTOF cluster + #include "DataFormatsIOTOF/Cluster.h" -#include +#include "Framework/Logger.h" +#include +#include +#include +// Root ClassImp macros for serialization metadata +ClassImp(o2::iotof::ClusterInfo); ClassImp(o2::iotof::Cluster); -namespace o2::iotof +namespace o2 +{ +namespace iotof { std::string Cluster::asString() const { - std::ostringstream stream; - stream << "chip=" << chipID << " row=" << row << " col=" << col << " size=" << size; - return stream.str(); + LOG(debug) << "[Cluster::asString] Converting Cluster to string"; + return std::format( + "chip: {:5d} | row: {:3d} col: {:3d} | span: {:2d}x{:2d} | pattern: {:5d} topology: {:4d}", + getChipID(), + getRow(), + getCol(), + getRowSpan(), + getColSpan(), + getPattern(), + getTopology() + ); +} + +//______________________________________________________________________________ +void Cluster::print() const +{ + std::cout << *this << "\n"; +} + +//______________________________________________________________________________ +void Cluster::sanityCheck() +{ + LOG(debug) << "[Cluster::sanityCheck] Performing sanity check on Cluster fields"; + + // Ensure extracted values fit within allowed bit masks + assert(getRow() <= ClusterInfo::MaskRow); + assert(getCol() <= ClusterInfo::MaskCol); + assert(getRowSpan() <= ClusterInfo::MaskRowSpan); + assert(getColSpan() <= ClusterInfo::MaskColSpan); + assert(getPattern() <= ClusterInfo::MaskPattern); + assert(getTopology() <= ClusterInfo::MaskTopology); } -} // namespace o2::iotof +} // namespace iotof +} // namespace o2 + +// Stream operator implementation +std::ostream& operator<<(std::ostream& stream, const o2::iotof::Cluster& cl) +{ + stream << cl.asString(); + return stream; +} diff --git a/Detectors/Upgrades/ALICE3/IOTOF/DataFormatsIOTOF/src/DataFormatsIOTOFLinkDef.h b/Detectors/Upgrades/ALICE3/IOTOF/DataFormatsIOTOF/src/DataFormatsIOTOFLinkDef.h index 7e121273d3fab..e639584ebfa75 100644 --- a/Detectors/Upgrades/ALICE3/IOTOF/DataFormatsIOTOF/src/DataFormatsIOTOFLinkDef.h +++ b/Detectors/Upgrades/ALICE3/IOTOF/DataFormatsIOTOF/src/DataFormatsIOTOFLinkDef.h @@ -18,6 +18,7 @@ #pragma link C++ class o2::iotof::Digit + ; #pragma link C++ class std::vector < o2::iotof::Digit> + ; +#pragma link C++ class o2::iotof::ClusterInfo + ; #pragma link C++ class o2::iotof::Cluster + ; #pragma link C++ class std::vector < o2::iotof::Cluster> + ; diff --git a/Detectors/Upgrades/ALICE3/IOTOF/base/CMakeLists.txt b/Detectors/Upgrades/ALICE3/IOTOF/base/CMakeLists.txt index 3b47b9451916d..c5c2b1c36bcab 100644 --- a/Detectors/Upgrades/ALICE3/IOTOF/base/CMakeLists.txt +++ b/Detectors/Upgrades/ALICE3/IOTOF/base/CMakeLists.txt @@ -11,10 +11,12 @@ o2_add_library(IOTOFBase SOURCES src/GeometryTGeo.cxx + src/Segmentation.cxx src/IOTOFBaseParam.cxx PUBLIC_LINK_LIBRARIES O2::DetectorsBase O2::MathUtils) o2_target_root_dictionary(IOTOFBase HEADERS include/IOTOFBase/GeometryTGeo.h + include/IOTOFBase/Segmentation.h include/IOTOFBase/IOTOFBaseParam.h) diff --git a/Detectors/Upgrades/ALICE3/IOTOF/simulation/include/IOTOFSimulation/Segmentation.h b/Detectors/Upgrades/ALICE3/IOTOF/base/include/IOTOFBase/Segmentation.h similarity index 93% rename from Detectors/Upgrades/ALICE3/IOTOF/simulation/include/IOTOFSimulation/Segmentation.h rename to Detectors/Upgrades/ALICE3/IOTOF/base/include/IOTOFBase/Segmentation.h index ddde28cf7dd7a..c726998fcd4bc 100644 --- a/Detectors/Upgrades/ALICE3/IOTOF/simulation/include/IOTOFSimulation/Segmentation.h +++ b/Detectors/Upgrades/ALICE3/IOTOF/base/include/IOTOFBase/Segmentation.h @@ -50,11 +50,11 @@ class Segmentation /// the center of the sensitive volulme. /// \param int iRow Detector x cell coordinate. Has the range 0 <= iRow < mNumberOfRows /// \param int iCol Detector z cell coordinate. Has the range 0 <= iCol < mNumberOfColumns - bool localToDetector(float x, float z, int& iRow, int& iCol, const int subDetectorID); + bool localToDetector(float x, float z, int& iRow, int& iCol, const int subDetectorID) const; /// same but w/o check for row/column range - void localToDetectorUnchecked(float xRow, float zCol, int& iRow, int& iCol, const int subDetectorID); + void localToDetectorUnchecked(float xRow, float zCol, int& iRow, int& iCol, const int subDetectorID) const; - /// Transformation from Detector cell coordiantes to Geant detector centered + /// Transformation from Detector cell coordinates to Geant detector centered /// local coordinates (cm) /// \param int iRow Detector x cell coordinate. Has the range 0 <= iRow < mNumberOfRows /// \param int iCol Detector z cell coordinate. Has the range 0 <= iCol < mNumberOfColumns @@ -67,7 +67,7 @@ class Segmentation // w/o check for row/col range template - void detectorToLocalUnchecked(L row, L col, T& xRow, T& zCol, const int subDetectorID) + void detectorToLocalUnchecked(L row, L col, T& xRow, T& zCol, const int subDetectorID) const { if (subDetectorID != 0 && subDetectorID != 1) { row = col = -1; @@ -78,7 +78,7 @@ class Segmentation zCol = col * specsConfig.PitchCol + getFirstColCoordinate(subDetectorID); } template - void detectorToLocalUnchecked(L row, L col, math_utils::Point3D& loc, const int subDetectorID) + void detectorToLocalUnchecked(L row, L col, math_utils::Point3D& loc, const int subDetectorID) const { if (subDetectorID != 0 && subDetectorID != 1) { row = col = -1; @@ -88,7 +88,7 @@ class Segmentation loc.SetCoordinates(getFirstRowCoordinate(subDetectorID) - row * specsConfig.PitchRow, T(0.), col * specsConfig.PitchCol + getFirstColCoordinate(subDetectorID)); } template - void detectorToLocalUnchecked(L row, L col, std::array& loc, const int subDetectorID) + void detectorToLocalUnchecked(L row, L col, std::array& loc, const int subDetectorID) const { if (subDetectorID != 0 && subDetectorID != 1) { row = col = -1; @@ -103,7 +103,7 @@ class Segmentation // same but with check for row/col range template - bool detectorToLocal(L row, L col, T& xRow, T& zCol, const int subDetectorID) + bool detectorToLocal(L row, L col, T& xRow, T& zCol, const int subDetectorID) const { if (subDetectorID != 0 && subDetectorID != 1) { row = col = -1; @@ -118,7 +118,7 @@ class Segmentation } template - bool detectorToLocal(L row, L col, math_utils::Point3D& loc, const int subDetectorID) + bool detectorToLocal(L row, L col, math_utils::Point3D& loc, const int subDetectorID) const { if (subDetectorID != 0 && subDetectorID != 1) { row = col = -1; @@ -132,7 +132,7 @@ class Segmentation return true; } template - bool detectorToLocal(L row, L col, std::array& loc, const int subDetectorID) + bool detectorToLocal(L row, L col, std::array& loc, const int subDetectorID) const { if (subDetectorID != 0 && subDetectorID != 1) { row = col = -1; @@ -146,12 +146,12 @@ class Segmentation return true; } - float getFirstRowCoordinate(const int subDetectorID) + float getFirstRowCoordinate(const int subDetectorID) const { const auto& specsConfig = ChipSpecificsParam::Instance(); return 0.5 * ((specsConfig.ActiveMatrixSizeRows() - specsConfig.PassiveEdgeTop + specsConfig.PassiveEdgeReadOut) - specsConfig.PitchRow); } - float getFirstColCoordinate(const int subDetectorID) + float getFirstColCoordinate(const int subDetectorID) const { const auto& specsConfig = ChipSpecificsParam::Instance(); return 0.5 * (specsConfig.PitchCol - specsConfig.ActiveMatrixSizeCols()); @@ -161,7 +161,7 @@ class Segmentation }; //_________________________________________________________________________________________________ -inline void Segmentation::localToDetectorUnchecked(float xRow, float zCol, int& iRow, int& iCol, const int subDetectorID) +inline void Segmentation::localToDetectorUnchecked(float xRow, float zCol, int& iRow, int& iCol, const int subDetectorID) const { // convert to row/col w/o over/underflow check if (subDetectorID != 0 && subDetectorID != 1) { @@ -187,7 +187,7 @@ inline void Segmentation::localToDetectorUnchecked(float xRow, float zCol, int& } //_________________________________________________________________________________________________ -inline bool Segmentation::localToDetector(float xRow, float zCol, int& iRow, int& iCol, const int subDetectorID) +inline bool Segmentation::localToDetector(float xRow, float zCol, int& iRow, int& iCol, const int subDetectorID) const { // convert to row/col if (subDetectorID != 0 && subDetectorID != 1) { diff --git a/Detectors/Upgrades/ALICE3/IOTOF/base/src/GeometryTGeo.cxx b/Detectors/Upgrades/ALICE3/IOTOF/base/src/GeometryTGeo.cxx index 8c8a36877eca8..09c73d21fdd67 100644 --- a/Detectors/Upgrades/ALICE3/IOTOF/base/src/GeometryTGeo.cxx +++ b/Detectors/Upgrades/ALICE3/IOTOF/base/src/GeometryTGeo.cxx @@ -313,13 +313,35 @@ void GeometryTGeo::Build(int loadTrans) } LOG(info) << "TF3 geometry: numberOfChipsITOF = " << mNumberOfChipsIOTOF[0] << ", numberOfChipsOTOF = " - << mNumberOfChipsIOTOF[1] << ", numberOfChips = " << numberOfChips << ", mNumberOfChipsPerStaveITOF" + << mNumberOfChipsIOTOF[1] << ", numberOfChips = " << numberOfChips << ", mNumberOfChipsPerStaveITOF = " << mNumberOfChipsPerStaveIOTOF[0]; setSize(numberOfChips); defineSensors(); fillTrackingFramesCache(); fillMatrixCache(loadTrans); + for (int j{0}; j < 2; ++j) { + LOG(info) << "Geometry for layer " << j << ": numberOfStaves = " << mNumberOfStavesIOTOF[j] << ", numberOfSubStaves = " << mNumberOfSubStavesIOTOF[j] + << ", numberOfModules = " << mNumberOfModulesIOTOF[j] << ", numberOfChipsPerModule = " << mNumberOfChipsPerModuleIOTOF[j] + << ", numberOfChipsPerSubStave = " << mNumberOfChipsPerSubStaveIOTOF[j] << ", numberOfChipsPerStave = " + << mNumberOfChipsPerStaveIOTOF[j] << ", numberOfChips = " << mNumberOfChipsIOTOF[j]; + } + const auto& specs = ChipSpecificsParam::Instance(); + LOG(info) << "specs.NCols = " << specs.NCols; + LOG(info) << "specs.NRows = " << specs.NRows; + LOG(info) << "specs.PitchCol = " << specs.PitchCol; + LOG(info) << "specs.PitchRow = " << specs.PitchRow; + LOG(info) << "specs.PassiveEdgeReadOut = " << specs.PassiveEdgeReadOut; + LOG(info) << "specs.PassiveEdgeTop = " << specs.PassiveEdgeTop; + LOG(info) << "specs.PassiveEdgeSide = " << specs.PassiveEdgeSide; + LOG(info) << "specs.PixelPassiveEdgeX = " << specs.PixelPassiveEdgeX; + LOG(info) << "specs.PixelPassiveEdgeZ = " << specs.PixelPassiveEdgeZ; + LOG(info) << "specs.SensorLayerThicknessEff = " << specs.SensorLayerThicknessEff; + LOG(info) << "specs.SensorLayerThickness = " << specs.SensorLayerThickness; + LOG(info) << "specs.NPixels = " << specs.NPixels(); + LOG(info) << "SensorSizeCols = " << specs.SensorSizeCols(); + LOG(info) << "SensorSizeRows = " << specs.SensorSizeRows(); + // fillMatrixCache(o2::math_utils::bit2Mask(o2::math_utils::TransformType::L2G)); } diff --git a/Detectors/Upgrades/ALICE3/IOTOF/base/src/IOTOFBaseLinkDef.h b/Detectors/Upgrades/ALICE3/IOTOF/base/src/IOTOFBaseLinkDef.h index cb5b047e72077..ba9457a4b96c9 100644 --- a/Detectors/Upgrades/ALICE3/IOTOF/base/src/IOTOFBaseLinkDef.h +++ b/Detectors/Upgrades/ALICE3/IOTOF/base/src/IOTOFBaseLinkDef.h @@ -16,6 +16,7 @@ #pragma link off all functions; #pragma link C++ class o2::iotof::GeometryTGeo + ; +#pragma link C++ class o2::iotof::Segmentation + ; #pragma link C++ class o2::iotof::IOTOFBaseParam + ; #pragma link C++ class o2::conf::ConfigurableParamHelper < o2::iotof::IOTOFBaseParam> + ; diff --git a/Detectors/Upgrades/ALICE3/IOTOF/simulation/src/Segmentation.cxx b/Detectors/Upgrades/ALICE3/IOTOF/base/src/Segmentation.cxx similarity index 96% rename from Detectors/Upgrades/ALICE3/IOTOF/simulation/src/Segmentation.cxx rename to Detectors/Upgrades/ALICE3/IOTOF/base/src/Segmentation.cxx index a7ec0d708c3b8..aa77bf50d069d 100644 --- a/Detectors/Upgrades/ALICE3/IOTOF/simulation/src/Segmentation.cxx +++ b/Detectors/Upgrades/ALICE3/IOTOF/base/src/Segmentation.cxx @@ -12,7 +12,7 @@ /// \file Segmentation.cxx /// \brief Implementation of the Segmentation class -#include "IOTOFSimulation/Segmentation.h" +#include "IOTOFBase/Segmentation.h" #include "IOTOFBase/IOTOFBaseParam.h" #include diff --git a/Detectors/Upgrades/ALICE3/IOTOF/macros/CheckClustersIOTOF.C b/Detectors/Upgrades/ALICE3/IOTOF/macros/CheckClustersIOTOF.C index 107e5a4d02bf8..cd8d7e31ad227 100644 --- a/Detectors/Upgrades/ALICE3/IOTOF/macros/CheckClustersIOTOF.C +++ b/Detectors/Upgrades/ALICE3/IOTOF/macros/CheckClustersIOTOF.C @@ -13,235 +13,1087 @@ /// \brief Simple macro to create clusters from TF3 digits #if !defined(__CLING__) || defined(__ROOTCLING__) + +#include + #include #include #include +#include #include #include #include -#include "IOTOFSimulation/Segmentation.h" #include "IOTOFBase/IOTOFBaseParam.h" #include "IOTOFBase/GeometryTGeo.h" +#include "IOTOFBase/Segmentation.h" +#include "IOTOFSimulation/Chip.h" +#include "IOTOFReconstruction/TopologyClassifier.h" +#include "ITSMFTSimulation/Hit.h" #include "DataFormatsIOTOF/Digit.h" #include "DataFormatsIOTOF/Cluster.h" #include "MathUtils/Utils.h" #include "SimulationDataFormat/ConstMCTruthContainer.h" #include "SimulationDataFormat/IOMCTruthContainerView.h" #include "SimulationDataFormat/MCCompLabel.h" +#include "SimulationDataFormat/MCTrack.h" +#include "SimulationDataFormat/TrackReference.h" +#include "SimulationDataFormat/MCEventHeader.h" #include "DetectorsBase/GeometryManager.h" #include "DataFormatsITSMFT/ROFRecord.h" #endif -#define ENABLE_UPGRADES +using namespace o2::base; +using namespace o2::iotof; +using o2::iotof::Digit; +using o2::iotof::Cluster; + +struct ClusterProperties { + int clsIdx = -1; + int eventID = -1; + int trackID = -1; + int chipID = -1; + int layer = -1; + uint16_t pattern = 0; + int rowStart = 0; + uint8_t rowSpan = 0; + int colStart = 0; + uint8_t colSpan = 0; + int size = 0; + bool isPrimary = false; + bool isFake = false; + bool isFakeDiffHits = false; + bool isFakeDiffTrks = false; + bool isFakeDiffEvts = false; + int hitIdx = -1; + Topologies topology = kOther; + uint32_t topoKey = 0; +}; + +struct HitData { + int hitIdx = -1; // In the hitsPerEvent[iEvt] array + std::vector assocClsIdxs{}; + std::vector assocDigitIdxs{}; +}; + +struct TrackData { + std::unordered_map> hitsByDetector; +}; + +void GetHitAvgPositionGlobal(const o2::itsmft::Hit& hit, o2::math_utils::Point3D& avgPos) { + + o2::math_utils::Point3D startPos = hit.GetPosStart(); + o2::math_utils::Point3D endPos = hit.GetPos(); + + avgPos = o2::math_utils::Point3D((startPos.X() + endPos.X()) / 2, (startPos.Y() + endPos.Y()) / 2, (startPos.Z() + endPos.Z()) / 2); +} + + + +void GetHitAvgPositionLocal(const o2::itsmft::Hit& hit, o2::iotof::GeometryTGeo* geom, o2::math_utils::Point3D& avgPos) { + + const int chipID = hit.GetDetectorID(); + + o2::math_utils::Point3D startPos = hit.GetPosStart(); + auto startPosLocal = geom->getMatrixL2G(chipID) ^ (startPos); + o2::math_utils::Point3D endPos = hit.GetPos(); + auto endPosLocal = geom->getMatrixL2G(chipID) ^ (endPos); + + avgPos = o2::math_utils::Point3D((startPosLocal.X() + endPosLocal.X()) / 2, (startPosLocal.Y() + endPosLocal.Y()) / 2, (startPosLocal.Z() + endPosLocal.Z()) / 2); +} + + +void GetDigitGlobalPos(const Digit& digit, + o2::math_utils::Point3D& globalPos, + o2::iotof::GeometryTGeo* geom, + o2::iotof::Segmentation* segm) { + const int chipID = digit.getChipIndex(); + const int layer = geom->getIOTOFLayer(chipID); + + float x = 0.f; + float z = 0.f; + if (layer >= 0) + segm->detectorToLocal(digit.getRow(), digit.getColumn(), x, z, layer); + + globalPos = geom->getMatrixL2G(chipID)(o2::math_utils::Point3D{x, 0.f, z}); +} + + +void PrintMcTrack(bool verbose, const o2::MCTrack& mcTrack) { + if (!verbose) { + return; + } + std::cout << "MCTrack: pdgCode = " << mcTrack.GetPdgCode() << ", isPrimary = " << mcTrack.isPrimary() << ", process: " << mcTrack.getProcess() << ", pt = " << mcTrack.GetPt() << ", eta = " << mcTrack.GetEta() << ", phi = " << mcTrack.GetPhi() << std::endl; +} + + +void PrintHit(bool verbose, o2::itsmft::Hit hit, o2::iotof::GeometryTGeo* iotofGeom) { + if (!verbose) { + return; + } + int layer{-1}, stave{-1}, subStave{-1}, module{-1}, chip{-1}; + iotofGeom->getIOTOFChipId(hit.GetDetectorID(), layer, stave, subStave, module, chip); + std::cout << "Hit: detectorID = " << hit.GetDetectorID() << ", layer = " << layer << ", stave = " << stave << ", subStave = " << subStave << ", module = " << module << ", chip = " << chip << ", trackID = " << hit.GetTrackID() << ", X = " << hit.GetX() << ", Y = " << hit.GetY() << ", Z = " << hit.GetZ() << ", time = " << hit.GetTime() << std::endl; +} + + +void PrintDigit(bool verbose, const o2::iotof::Digit& digit, auto& labels, o2::iotof::GeometryTGeo* iotofGeom, o2::iotof::Segmentation* segmInfo) { + if (!verbose) { + return; + } + + if (labels.empty()) { + std::cout << "Digit: no MCCompLabel associated, chipID = " << digit.getChipIndex() << ", row = " << digit.getRow() << ", col = " << digit.getColumn() << ", charge = " << digit.getCharge() << ", time = " << digit.getTime() << std::endl; + return; + } + const auto& evtTrackLabel = labels[0]; + if (!evtTrackLabel.isValid()) { + std::cout << "Digit: invalid MCCompLabel, chipID = " << digit.getChipIndex() << ", row = " << digit.getRow() << ", col = " << digit.getColumn() << ", charge = " << digit.getCharge() << ", time = " << digit.getTime() << std::endl; + return; + } + + const int eventID = evtTrackLabel.getEventID(); + const int trackID = evtTrackLabel.getTrackID(); + + int layer{-1}, stave{-1}, subStave{-1}, module{-1}, chip{-1}; + iotofGeom->getIOTOFChipId(digit.getChipIndex(), layer, stave, subStave, module, chip); + o2::math_utils::Point3D digitPos; + GetDigitGlobalPos(digit, digitPos, iotofGeom, segmInfo); + std::cout << "Digit: trackID = " << trackID << ", eventID = " << eventID << ", chipID = " + << digit.getChipIndex() << ", layer = " << layer << ", stave = " << stave + << ", subStave = " << subStave << ", module = " << module << ", chip = " << chip + << ", row = " << digit.getRow() << ", col = " << digit.getColumn() << ", charge = " << digit.getCharge() + << ", time = " << digit.getTime() << ", global position = (" << digitPos.X() << ", " << digitPos.Y() + << ", " << digitPos.Z() << ")" << std::endl; +} + + +void PrintCluster(bool verbose, + const o2::iotof::Cluster& cluster, + auto clsLabel, + o2::iotof::GeometryTGeo* iotofGeom, + o2::iotof::Segmentation* segmInfo) { + if (!verbose) { + return; + } + + if (clsLabel.empty()) + return; + + std::cout << "Cluster: " << clsLabel.size() << " MCCompLabels, chipID=" << cluster.getChipID() << ", row=" << cluster.getRow() << ", col=" << cluster.getCol() << ", rowSpan=" << cluster.getRowSpan() << ", colSpan=" << cluster.getColSpan() << ", topology=" << cluster.getTopology() << std::endl; + for (int iLabel = 0; iLabel < clsLabel.size(); ++iLabel) { + const auto& evtTrackLabel = clsLabel[iLabel]; + if (!evtTrackLabel.isValid()) + continue; + + const int eventID = evtTrackLabel.getEventID(); + const int trackID = evtTrackLabel.getTrackID(); + std::cout << " Label " << iLabel << ", eventID=" << eventID << ", trackID=" << trackID << std::endl; + } +} + + +template +void Print(bool verbose, Args&&... args) { + if (!verbose) { + return; + } + + (std::cout << ... << std::forward(args)) << std::endl; +} + + +void GetClusterGlobalPos(const o2::iotof::Cluster& cluster, + TopologyInfo topoInfo, + o2::math_utils::Point3D& globalPos, + o2::iotof::GeometryTGeo* iotofGeom, + o2::iotof::Segmentation* segmInfo){ + + float x = 0.f; + float y = 0.f; + float z = 0.f; + int rowCOG = cluster.getRow() + topoInfo.mOffsetXToCOG; + int colCOG = cluster.getCol() + topoInfo.mOffsetZToCOG; + segmInfo->detectorToLocal(rowCOG, colCOG, x, z, cluster.getChipID()); + globalPos = iotofGeom->getMatrixL2G(cluster.getChipID())(o2::math_utils::Point3D{x, 0.f, z}); +} -void CheckClustersIOTOF(std::string digiFilePath = "tf3digits.root", std::string clsFilePath = "tf3clusters.root", std::string inputGeomPath = "o2sim_geometry.root") -{ - gStyle->SetPalette(55); - using namespace o2::base; - using namespace o2::iotof; +int FindBestMatchingHit(const o2::iotof::Cluster& cluster, + TopologyInfo topoInfo, + std::vector& chipHitsIdxs, + std::vector* evtChipHits, + const std::vector* digitsArray, + o2::iotof::GeometryTGeo* iotofGeom, + o2::iotof::Segmentation* segmInfo){ + int bestHitIdx = -1; + float minDistanceSq = std::numeric_limits::max(); + o2::math_utils::Point3D clsPos; + GetClusterGlobalPos(cluster, topoInfo, clsPos, iotofGeom, segmInfo); - using o2::iotof::Cluster; - using o2::iotof::Digit; + for (int i = 0; i < chipHitsIdxs.size(); ++i) { + const auto& hit = (*evtChipHits)[chipHitsIdxs[i].hitIdx]; + + float dx = clsPos.X() - hit.GetX(); + float dy = clsPos.Y() - hit.GetY(); + float dz = clsPos.Z() - hit.GetZ(); + float distSq = dx*dx + dy*dy + dz*dz; + + if (distSq < minDistanceSq) { + minDistanceSq = distSq; + bestHitIdx = i; + } + } + + return bestHitIdx; // Returns -1 if no hit is within maxToleranceCm (true fake cluster) +} + + +void CheckClustersIOTOF(std::string kinefile = "o2sim_Kine.root", + std::string hitfile = "o2sim_HitsTF3.root", + std::string digiFilePath = "tf3digits.root", + std::string clsFilePath = "tf3clusters.root", + std::string clsFileTopoPath = "TF3ClustersTopologies.root", + std::string inputGeomPath = "o2sim_geometry.root", + bool verbose = false) +{ + Print(verbose, "CheckClustersTopologiesIOTOF: kinefile = ", kinefile, ", hitfile = ", hitfile, ", digiFilePath = ", digiFilePath, ", clsFilePath = ", clsFilePath, ", inputGeomPath = ", inputGeomPath); + gStyle->SetPalette(55); o2::conf::ConfigurableParam::updateFromString("IOTOFBase.segmentedInnerTOF=true;IOTOFBase.segmentedOuterTOF=true;IOTOFBase.enableForwardTOF=false;IOTOFBase.enableBackwardTOF=false"); - auto segGeom = o2::iotof::Segmentation::Instance(); + auto segmInfo = o2::iotof::Segmentation::Instance(); // Geometry o2::base::GeometryManager::loadGeometry(inputGeomPath); - auto* tofGeo = o2::iotof::GeometryTGeo::Instance(); - tofGeo->fillMatrixCache(o2::math_utils::bit2Mask(o2::math_utils::TransformType::L2G)); - - // Digits - TFile* digiFile = TFile::Open(digiFilePath.data()); - TTree* digiTree = (TTree*)digiFile->Get("o2sim"); - std::vector* digitsArray{nullptr}; - digiTree->SetBranchAddress("TF3Digit", &digitsArray); - std::vector* digiRofRecordsArr{nullptr}; - digiTree->SetBranchAddress("TF3DigitROF", &digiRofRecordsArr); - auto& digiRofArr = *digiRofRecordsArr; - o2::dataformats::IOMCTruthContainerView* digiLabelsArr{nullptr}; - digiTree->SetBranchAddress("TF3DigitMCTruth", &digiLabelsArr); - digiTree->GetEntry(0); - o2::dataformats::ConstMCTruthContainer digiLabels; - digiLabelsArr->copyandflatten(digiLabels); - - // Clusters - TFile* clsFile = TFile::Open(clsFilePath.data()); - TTree* clsTree = (TTree*)clsFile->Get("o2sim"); - std::vector* clsArray{nullptr}; - clsTree->SetBranchAddress("TF3ClusterComp", &clsArray); - std::vector* clsRofRecordsArr{nullptr}; - clsTree->SetBranchAddress("TF3ClusterROF", &clsRofRecordsArr); - auto& clsRofArr = *clsRofRecordsArr; - o2::dataformats::MCTruthContainer* clsLabels{nullptr}; - clsTree->SetBranchAddress("TF3ClusterMCTruth", &clsLabels); - clsTree->GetEntry(0); - - // Summary of entries in all branches + auto* iotofGeom = o2::iotof::GeometryTGeo::Instance(); + iotofGeom->fillMatrixCache(o2::math_utils::bit2Mask(o2::math_utils::TransformType::L2G)); + + // Cluster topologies dictionary + TFile* clsTopoFile = TFile::Open(clsFileTopoPath.data(), "READ"); + auto* clsTopoMapPtr = clsTopoFile->Get>("TF3ClusterTopologies"); + if (clsTopoMapPtr) { + std::cout << "Loaded " << clsTopoMapPtr->size() << " entries from TF3ClusterTopologies.root" << std::endl; + } else { + std::cerr << "Failed to load TF3ClusterTopologies from file!" << std::endl; + } + clsTopoFile->Close(); + std::cout << std::endl; - std::cout << "---> Number of digits: " << digitsArray->size() << std::endl; - std::cout << "---> Number of digit ROFs: " << digiRofArr.size() << std::endl; - std::cout << "---> Number of clusters: " << clsArray->size() << std::endl; - std::cout << "---> Number of cluster ROFs: " << clsRofArr.size() << std::endl; - std::cout << "---> Number of digits with MC label: " << digiLabels.getNElements() << std::endl; - std::cout << "---> Number of digits with MC label: " << digiLabels.getIndexedSize() << std::endl; - std::cout << "---> Number of clusters with MC label: " << clsLabels->getNElements() << std::endl; - std::cout << "---> Number of clusters with MC label: " << clsLabels->getIndexedSize() << std::endl; + std::cout << "Topologies summary: " << std::endl; + TopologyClassifier topoClassifier(*clsTopoMapPtr); + topoClassifier.print(); std::cout << std::endl; - auto clsTuple = new TNtuple("clsTuple", "clsTuple", "chip_id:x:y:z:row:col:time"); - clsTuple->SetDirectory(nullptr); - - TH1F* histXCoordCls = new TH1F("histXCoordCls", "histXCoordCls", 8000, -100, 100); - TH1F* histYCoordCls = new TH1F("histYCoordCls", "histYCoordCls", 8000, -100, 100); - TH1F* histZCoordCls = new TH1F("histZCoordCls", "histZCoordCls", 28000, -400, 400); - TH1F* histXCoordDigit = new TH1F("histXCoordDigit", "histXCoordDigit", 8000, -100, 100); - TH1F* histYCoordDigit = new TH1F("histYCoordDigit", "histYCoordDigit", 8000, -100, 100); - TH1F* histZCoordDigit = new TH1F("histZCoordDigit", "histZCoordDigit", 28000, -400, 400); - TH1F* histXCoordRes = new TH1F("histXCoordRes", "histXCoordRes", 100, -0.05, 0.05); - TH1F* histYCoordRes = new TH1F("histYCoordRes", "histYCoordRes", 100, -0.05, 0.05); - TH1F* histZCoordRes = new TH1F("histZCoordRes", "histZCoordRes", 100, -0.05, 0.05); - TH1F* histTimeRes = new TH1F("histTimeRes", "histTimeRes", 100, -0.05, 0.05); - - // Load all digits upfront and build a lookup map - int nDigits = digiTree->GetEntries(); - std::unordered_map digitsLabels; - for (int iDigit = 0; iDigit < digitsArray->size(); ++iDigit) { - auto label = digiLabels.getLabels(iDigit)[0]; - if (!label.isValid()) { + // Generated MC tracks and TrackRefs information + TFile* kineFile = TFile::Open(kinefile.data()); + TTree* kineTree = (TTree*)kineFile->Get("o2sim"); + const int nEvts = kineTree->GetEntries(); + std::vector*> mcTracksPerEvent(nEvts, nullptr); + std::vector*> mcTracksRefsPerEvent(nEvts, nullptr); + + // Hits information + TFile* hitFile = TFile::Open(hitfile.data()); + TTree* hitTree = (TTree*)hitFile->Get("o2sim"); + std::vector*> hitsPerEvent(nEvts, nullptr); + + // Digits information + TFile* digFile = TFile::Open(digiFilePath.data()); + TTree* digitsTree = (TTree*)digFile->Get("o2sim"); + std::vector* digitsArray = nullptr; + o2::dataformats::IOMCTruthContainerView* digitsLabelsArr = nullptr; + + digitsTree->SetBranchAddress("TF3Digit", &digitsArray); + digitsTree->SetBranchAddress("TF3DigitMCTruth", &digitsLabelsArr); + + // Clusters information + TFile* clsFile = TFile::Open(clsFilePath.data()); + TTree* clustersTree = (TTree*)clsFile->Get("o2sim"); + std::vector* clustersArray = nullptr; + std::vector* clustersPatternsArray = nullptr; + o2::dataformats::MCTruthContainer* clustersLabelsArr = nullptr; + + clustersTree->SetBranchAddress("TF3Cluster", &clustersArray); + clustersTree->SetBranchAddress("TF3ClusterPatt", &clustersPatternsArray); + clustersTree->SetBranchAddress("TF3ClusterMCTruth", &clustersLabelsArr); + + // Load hits and MC track refs, stored per-event + hitTree->SetBranchAddress("TF3Hit", &hitsPerEvent[0]); + kineTree->SetBranchAddress("MCTrack", &mcTracksPerEvent[0]); + kineTree->SetBranchAddress("TrackRefs", &mcTracksRefsPerEvent[0]); + for (int iEvt = 0; iEvt < nEvts; ++iEvt) { + hitTree->SetBranchAddress("TF3Hit", &hitsPerEvent[iEvt]); + hitTree->GetEntry(iEvt); + kineTree->SetBranchAddress("MCTrack", &mcTracksPerEvent[iEvt]); + kineTree->SetBranchAddress("TrackRefs", &mcTracksRefsPerEvent[iEvt]); + kineTree->GetEntry(iEvt); + Print(verbose, "Loaded hit event ", iEvt, " with ", hitsPerEvent[iEvt]->size(), " hits"); + } + + // Digits: TTree entries are not separated per-event, but all digits are stored in a single entry + digitsTree->GetEntry(0); + o2::dataformats::ConstMCTruthContainer digitsLabels; + digitsLabelsArr->copyandflatten(digitsLabels); + + // Clusters: TTree entries are not separated per-event, but all clusters are stored in a single entry + clustersTree->GetEntry(0); + o2::dataformats::ConstMCTruthContainer clustersLabels; + + // Store hit, digit and cluster properties for all tracks in all events + std::vector> allEvtsTrackData(nEvts); + TH2F* hEtaPhiHitsPrmTrkLayer0 = new TH2F("hEtaPhiHitsPrmTrkLayer0", "hEtaPhiHitsPrmTrkLayer0;#phi;#eta", 300, 0, 6.28319, 40, -2, 2); + TH2F* hEtaPhiHitsSecTrkLayer0 = new TH2F("hEtaPhiHitsSecTrkLayer0", "hEtaPhiHitsSecTrkLayer0;#phi;#eta", 300, 0, 6.28319, 40, -2, 2); + TH2F* hEtaPhiHitsPrmTrkLayer1 = new TH2F("hEtaPhiHitsPrmTrkLayer1", "hEtaPhiHitsPrmTrkLayer1;#phi;#eta", 300, 0, 6.28319, 40, -2, 2); + TH2F* hEtaPhiHitsSecTrkLayer1 = new TH2F("hEtaPhiHitsSecTrkLayer1", "hEtaPhiHitsSecTrkLayer1;#phi;#eta", 300, 0, 6.28319, 40, -2, 2); + // Load Hits, which are stored per-event + int nHits{0}, nHitsFromPrimaryTracks{0}, nHitsFromSecondaryTracks{0}; + Print(verbose, "\n\n----> Starting hits printouts ... "); + for (int iEvt = 0; iEvt < nEvts; ++iEvt) { + + Print(verbose, "Event ", iEvt, ": ", hitsPerEvent[iEvt]->size(), " hits"); + for (int iHit = 0; iHit < hitsPerEvent[iEvt]->size(); ++iHit) { + + const auto& hit = (*hitsPerEvent[iEvt])[iHit]; + const int trackID = hit.GetTrackID(); + const int chipIndex = hit.GetDetectorID(); + allEvtsTrackData[iEvt][trackID].hitsByDetector[chipIndex].push_back({iHit, {}, {}}); + nHits++; + + // Fill histograms + int hitLayer = iotofGeom->getIOTOFLayer(hit.GetDetectorID()); + auto& mcTrack = (*mcTracksPerEvent[iEvt])[trackID]; + bool isPrimary = mcTrack.isPrimary(); + if (isPrimary) nHitsFromPrimaryTracks++; + else nHitsFromSecondaryTracks++; + float genEta = mcTrack.GetEta(); + float genPhi = mcTrack.GetPhi(); + + if (hitLayer == 0 && isPrimary) { hEtaPhiHitsPrmTrkLayer0->Fill(genPhi, genEta); } + else if (hitLayer == 0 && !isPrimary) { hEtaPhiHitsSecTrkLayer0->Fill(genPhi, genEta); } + else if (hitLayer == 1 && isPrimary) { hEtaPhiHitsPrmTrkLayer1->Fill(genPhi, genEta); } + else { hEtaPhiHitsSecTrkLayer1->Fill(genPhi, genEta); } + + // PrintHit(verbose, hit, iotofGeom); + } + } + + // Debug prints for digits, use MCCompLabel to get event ID (getEventID()), track ID (getTrackID()) + Print(verbose, "\n\n----> Starting digits printouts ... "); + for (int iDigit = 0; iDigit < (int)digitsArray->size(); ++iDigit) { + + auto labels = digitsLabels.getLabels(iDigit); + if (labels.empty()) + continue; + const auto& evtTrackLabel = labels[0]; + if (!evtTrackLabel.isValid()) + continue; + + const int eventID = evtTrackLabel.getEventID(); + const int trackID = evtTrackLabel.getTrackID(); + + if (eventID < 0 || eventID >= nEvts) { + std::cerr << "WARNING: digit " << iDigit << " has invalid eventID=" << eventID << "\n"; + continue; + } + + const auto& digit = (*digitsArray)[iDigit]; + const auto& digitLabels = digitsLabels.getLabels(iDigit); + // PrintDigit(verbose, digit, digitLabels, iotofGeom, segmInfo); + auto& hitList = allEvtsTrackData[eventID][trackID].hitsByDetector[digit.getChipIndex()]; + for (auto& hit : hitList) { + hit.assocDigitIdxs.push_back(iDigit); + } + } + + // Debug prints for clusters, use MCCompLabel to get event ID (getEventID()), track ID (getTrackID()) + Print(verbose, "\n\n----> Starting clusters printouts ... "); + for (int iCls = 0; iCls < (int)clustersArray->size(); ++iCls) { + + const auto& cls = (*clustersArray)[iCls]; + const auto& clsLabels = clustersLabelsArr->getLabels(iCls); + + if (clsLabels.empty()) continue; + const auto& evtTrackLabel = clsLabels[0]; + if (!evtTrackLabel.isValid()) + continue; + + const int eventID = evtTrackLabel.getEventID(); + const int trackID = evtTrackLabel.getTrackID(); + + if (eventID < 0 || eventID >= nEvts) { + std::cerr << "WARNING: cluster " << iCls << " has invalid eventID=" << eventID << "\n"; + continue; + } + + // PrintCluster(verbose, cls, clsLabels, iotofGeom, segmInfo); + auto& hitList = allEvtsTrackData[eventID][trackID].hitsByDetector[cls.getChipID()]; + for (auto& hit : hitList) { + hit.assocClsIdxs.push_back(iCls); } - digitsLabels.emplace(label, iDigit); } - // LOOP on : ROFRecord array - for (unsigned int iROF = 0; iROF < clsRofArr.size(); ++iROF) { + // Debug print of allEvtsTrackData structure + for (int iEvt = 0; iEvt < nEvts; ++iEvt) { + for (const auto& [trackID, trackData] : allEvtsTrackData[iEvt]) { + Print(verbose, "\n\n\nEvent ", iEvt, ", Track ", trackID, ":"); + for (const auto& [chipID, hitsInfos] : trackData.hitsByDetector) { + Print(verbose, "-----------\n", "Chip ", chipID, ": ", hitsInfos.size(), " hits"); + for (const auto& hitInfo : hitsInfos) { + Print(verbose, "\nHit ", hitInfo.hitIdx, ": ", hitInfo.assocDigitIdxs.size(), " digits, ", hitInfo.assocClsIdxs.size(), " clusters"); + for (int iDigit=0; iDigitgetLabels(hitInfo.assocClsIdxs[iCls]); + PrintCluster(verbose, cls, clsLabels, iotofGeom, segmInfo); + } + } + } + } + } + + // Debug prints + std::cout << "\n***********************************" << std::endl; + Print(true, "Number of events: ", nEvts); + Print(true, "Number of hits: ", nHits); + Print(true, "-> from primary tracks: ", nHitsFromPrimaryTracks); + Print(true, "-> from secondary tracks: ", nHitsFromSecondaryTracks); + Print(true, "Number of digits: ", digitsArray->size()); + Print(true, "Number of digit labels: ", digitsLabels.getNElements()); + Print(true, "Number of entries in digit tree: ", digitsTree->GetEntries()); + Print(true, "Number of clusters: ", clustersArray->size()); + Print(true, "Number of clusters labels: ", clustersLabelsArr->getNElements()); + Print(true, "Number of entries in cluster tree: ", clustersTree->GetEntries()); + std::cout << "***********************************\n" << std::endl; + + // Create vectors of digits with same chip index, cluster candidates + TH2F* hCountHitMatchingType = new TH2F("hCountHitMatchingType", "hCountHitMatchingType;Hit matching type;#it{p}_{T}", 4, -0.5, 3.5, 50, 0, 10); + hCountHitMatchingType->GetXaxis()->SetBinLabel(1, "Primary, 1 to 1"); + hCountHitMatchingType->GetXaxis()->SetBinLabel(2, "Secondary, 1 to 1"); + hCountHitMatchingType->GetXaxis()->SetBinLabel(3, "Primary, min distance"); + hCountHitMatchingType->GetXaxis()->SetBinLabel(4, "Secondary, min distance"); + + std::vector clustersProperties; + clustersProperties.reserve(clustersArray->size()); // Pre-allocate memory - const unsigned int rofIndex = clsRofArr[iROF].getFirstEntry(); - const unsigned int rofNEntries = clsRofArr[iROF].getNEntries(); + for (int iCls = 0; iCls < (int)clustersArray->size(); ++iCls) { - // LOOP on : digits array - std::cout << "\n\n ----> Starting loop on digits for ROF " << iROF << " with index " << rofIndex << " and nEntries " << rofNEntries << std::endl; - for (unsigned int iDigit = rofIndex; iDigit < rofIndex + rofNEntries; iDigit++) { - if (iDigit % 10000 == 0) { - std::cout << "Reading digit " << iDigit << " / " << digitsArray->size() << std::endl; + const auto& cluster = (*clustersArray)[iCls]; + + // Cluster labels + const auto& clsLabels = clustersLabelsArr->getLabels(iCls); + std::cout << "Processing cluster " << iCls << " with " << clsLabels.size() << " MCCompLabels associated." << std::endl; + if (clsLabels.empty()) { + std::cout << "---> Empty cls label" << std::endl; + continue; + } + + const auto& firstEvtTrackLabel = clsLabels[0]; + if (!firstEvtTrackLabel.isValid()) { + std::cout << "---> Invalid first evt-track label" << std::endl; + continue; + } + const int eventID = firstEvtTrackLabel.getEventID(); + const int trackID = firstEvtTrackLabel.getTrackID(); + + if (eventID < 0 || eventID >= nEvts) { + std::cerr << "WARNING: cluster " << iCls << " has invalid eventID=" << eventID << "\n"; + continue; + } + + ClusterProperties clsProps; + clsProps.clsIdx = iCls; + + // Cluster geometric properties + clsProps.chipID = cluster.getChipID(); + clsProps.layer = iotofGeom->getIOTOFLayer(cluster.getChipID()); + clsProps.rowStart = cluster.getRow(); + clsProps.rowSpan = cluster.getRowSpan(); + clsProps.colStart = cluster.getCol(); + clsProps.colSpan = cluster.getColSpan(); + clsProps.pattern = cluster.getPattern(); + clsProps.size = cluster.getSize(); + clsProps.topology = static_cast(cluster.getTopology()); + uint32_t clsTopoKey = (static_cast(clsProps.rowSpan) << 24) | + (static_cast(clsProps.colSpan) << 16) | + static_cast(clsProps.pattern); + clsProps.topoKey = clsTopoKey; + TopologyInfo clsTopoInfo = topoClassifier.getTopologyFeatures(clsProps.topoKey); + + // Cluster association properties + clsProps.eventID = eventID; + clsProps.trackID = trackID; + clsProps.isPrimary = false; + clsProps.isFake = false; + clsProps.isFakeDiffHits = false; + clsProps.isFakeDiffTrks = false; + clsProps.isFakeDiffEvts = false; + clsProps.hitIdx = -1; + + // 1 to 1 hit-cluster correspondence, set eventID and trackID for the cluster + if (clsLabels.size() > 1) { + // Multiple hits associated with the cluster, + // check consistency of track and event IDs across + // all digits in the cluster to flag fake clusters + for (int iLabel = 1; iLabel < clsLabels.size(); ++iLabel) { + const auto& evtTrackLabel = clsLabels[iLabel]; + + if (!evtTrackLabel.isValid()) { + continue; + } + + const int eventID = firstEvtTrackLabel.getEventID(); + const int trackID = firstEvtTrackLabel.getTrackID(); + + if (eventID < 0 || eventID >= nEvts) { + std::cerr << "WARNING: cluster " << iCls << " has invalid eventID=" << eventID << "\n"; + continue; + } + + if (evtTrackLabel.getEventID() != eventID) { + std::cout << "Cluster " << iCls << " has inconsistent event IDs across labels: " << evtTrackLabel.getEventID() << " != " << eventID << std::endl; + clsProps.isFake = true; + clsProps.isFakeDiffEvts = true; + } + if (evtTrackLabel.getTrackID() != trackID) { + std::cout << "Cluster " << iCls << " has inconsistent track IDs across labels: " << evtTrackLabel.getTrackID() << " != " << trackID << std::endl; + clsProps.isFake = true; + clsProps.isFakeDiffTrks = true; + } + } + } + + // Cluster-hit matching + if (!clsProps.isFake) { + + const auto& mcTrack = (*mcTracksPerEvent[clsProps.eventID])[clsProps.trackID]; + clsProps.isPrimary = mcTrack.isPrimary(); + + auto& chipHitsIdxs = allEvtsTrackData[clsProps.eventID][clsProps.trackID].hitsByDetector[clsProps.chipID]; + if (chipHitsIdxs.empty()) { + clsProps.hitIdx = -1; + } else if (chipHitsIdxs.size() == 1) { + clsProps.hitIdx = 0; + hCountHitMatchingType->Fill(clsProps.isPrimary ? 0 : 2, mcTrack.GetPt()); + } else { + // Perform spatial matching for multi-hit candidates + clsProps.hitIdx = FindBestMatchingHit(cluster, clsTopoInfo, chipHitsIdxs, hitsPerEvent[clsProps.eventID], digitsArray, iotofGeom, segmInfo); + hCountHitMatchingType->Fill(clsProps.isPrimary ? 1 : 3, mcTrack.GetPt()); } - Int_t iRow = (*digitsArray)[iDigit].getRow(); - Int_t iCol = (*digitsArray)[iDigit].getColumn(); - Int_t iDetID = (*digitsArray)[iDigit].getChipIndex(); - Int_t chipID = (*digitsArray)[iDigit].getChipIndex(); - Int_t subDetID = tofGeo->getIOTOFLayer(iDetID); + if (clsProps.hitIdx != -1) { + chipHitsIdxs[clsProps.hitIdx].assocClsIdxs.push_back(clustersProperties.size()); + } else { + clsProps.isFake = true; + clsProps.isFakeDiffHits = true; + std::cout << "Cluster " << iCls << " has no matching hit, marked as fake." << std::endl; + } + } + + // PrintCluster(verbose, cluster, digitsArray, digitsLabels, hitsPerEvent, iotofGeom, segmInfo); + clustersProperties.push_back(clsProps); + } + Print(true, "----> Total number of clusters: ", clustersProperties.size()); + + // QA printouts and histograms + Print(true, "\n\n----> Starting QA logging ... "); + const char* trackName[2] = {"Prm", "Sec"}; + + // Count fake clusters + TH1F* hCountFakeClusters[2][2]; + for (int layer = 0; layer < 2; ++layer) { + for (int type = 0; type < 2; ++type) { + hCountFakeClusters[layer][type] = new TH1F(Form("hCountFakeClusters%sTrkLayer%d", trackName[type], layer), Form("Fake Cluster Counter %s Trk Layer %d", trackName[type], layer), 6, -0.5, 5.5); + hCountFakeClusters[layer][type]->GetXaxis()->SetBinLabel(1, "Total"); + hCountFakeClusters[layer][type]->GetXaxis()->SetBinLabel(2, "Real"); + hCountFakeClusters[layer][type]->GetXaxis()->SetBinLabel(3, "Fake"); + hCountFakeClusters[layer][type]->GetXaxis()->SetBinLabel(4, "Fake NoHit"); + hCountFakeClusters[layer][type]->GetXaxis()->SetBinLabel(5, "Fake DiffTrks"); + hCountFakeClusters[layer][type]->GetXaxis()->SetBinLabel(6, "Fake DiffEvts"); + } + } - Float_t x{0.f}, y{0.f}, z{0.f}; - if (subDetID >= 0) { - segGeom->detectorToLocal(iRow, iCol, x, z, subDetID); + // Loop over clusters and fill histograms + for (const auto& cluster : clustersProperties) { + int layer = cluster.layer; + int type = cluster.isPrimary ? 0 : 1; + hCountFakeClusters[layer][type]->Fill(0.f, 1); // Total clusters + if (cluster.isFake) { + hCountFakeClusters[layer][type]->Fill(2.f, 1); // Fake clusters + if (cluster.isFakeDiffHits) { + hCountFakeClusters[layer][type]->Fill(3.f, 1); // Fake NoHit + } + if (cluster.isFakeDiffTrks) { + hCountFakeClusters[layer][type]->Fill(4.f, 1); // Fake DiffTrks } + if (cluster.isFakeDiffEvts) { + hCountFakeClusters[layer][type]->Fill(5.f, 1); // Fake DiffEvts + } + } else { + hCountFakeClusters[layer][type]->Fill(1.f, 1); // Real clusters + } + } + + Print(true, "----> hCountFakeClusters filled"); + // Topology names + const std::array topologyNames = { + "kSingleDigit", "kLineOnRow", "kLineOnCol", "kDiagonal", "kSquare", + "kUpperTriangleLeft", "kUpperTriangleRight", "kLowerTriangleLeft", + "kLowerTriangleRight", "kSnake", "kSnakeRot90", "kSnakeRefl", + "kSnakeRot90Refl", "kHuge", "kOther"}; + + // Count topologies from frequency values in + // topologies dictionary and fill the summary histograms + TH1F* hTopoSummaryDictionary = new TH1F("hTopoSummaryDictionary", "Cluster Topology Count Summary;;Counts", kNTopologies, 0, kNTopologies); + for (const auto& [topoKey, topology] : topoClassifier.getTopologyMap()) { + hTopoSummaryDictionary->Fill(topology.mTopology, topology.mFrequency); + } - o2::math_utils::Point3D localDigitCoord(x, y, z); // local Digit + TH2F *hTrueClsSizeVsEta[2][2], *hTrueClsSizeVsPhi[2][2], *hFakeClsSizeVsEta[2][2], *hFakeClsSizeVsPhi[2][2], + *hClustersEtaPhi[2][2], *hTopoVsEta[2][2], *hClsSizeVsTopo[2][2], *hXRes[2][2], *hYRes[2][2], *hZRes[2][2], + *hTrackHitsXY[2][2], *hTrackDoubleHitsXY[2][2], *hTrackDoubleHitsPhiPt[2][2], *hTopoVsEtaPt[2][2][kNTopologies]; + TH1F *hNClustersFromHit[2][2], *hMeanTrueClsSizeVsEta[2][2], *hMeanTrueClsSizeVsPhi[2][2], *hMeanFakeClsSizeVsEta[2][2], + *hMeanFakeClsSizeVsPhi[2][2], *hRmsXRes[2][2], *hRmsYRes[2][2], *hRmsZRes[2][2], *hMeanXRes[2][2], *hMeanYRes[2][2], + *hMeanZRes[2][2]; + TH1F* hTopoSummaryTotal = new TH1F("hTopoSummaryTotal", "Cluster Topology Summary;;Counts", kNTopologies, 0, kNTopologies); + TH1F* hTopoSummaryReal = new TH1F("hTopoSummaryReal", "Cluster Topology Summary;;Counts", kNTopologies, 0, kNTopologies); + TH1F* hTopoSummaryFake = new TH1F("hTopoSummaryFake", "Cluster Topology Summary;;Counts", kNTopologies, 0, kNTopologies); - const auto globalDigitCoord = tofGeo->getMatrixL2G(chipID)(localDigitCoord); // convert to global - histXCoordDigit->Fill(globalDigitCoord.X()); - histYCoordDigit->Fill(globalDigitCoord.Y()); - histZCoordDigit->Fill(globalDigitCoord.Z()); - } // end loop on digits array + Print(true, "----> Defining histograms"); + for (int layer = 0; layer < 2; ++layer) { + for (int type = 0; type < 2; ++type) { + hClustersEtaPhi[layer][type] = new TH2F(Form("hNClsVsEtaPhi%sTrkLayer%d", trackName[type], layer), "Cluster #eta vs #phi;#phi;#eta", 300, 0, 6.28319, 40, -2, 2); + hTrueClsSizeVsEta[layer][type] = new TH2F(Form("hTrueClsSizeVsEta%sTrkLayer%d", trackName[type], layer), "True Cluster Size vs #eta;#eta", 300, -2, 2, 20, 0.5, 20.5); + hTrueClsSizeVsPhi[layer][type] = new TH2F(Form("hTrueClsSizeVsPhi%sTrkLayer%d", trackName[type], layer), "True Cluster Size vs #phi;#phi", 300, 0, 6.28319, 20, 0.5, 20.5); + hFakeClsSizeVsEta[layer][type] = new TH2F(Form("hFakeClsSizeVsEta%sTrkLayer%d", trackName[type], layer), "Fake Cluster Size vs #eta;#eta", 300, -2, 2, 20, 0.5, 20.5); + hFakeClsSizeVsPhi[layer][type] = new TH2F(Form("hFakeClsSizeVsPhi%sTrkLayer%d", trackName[type], layer), "Fake Cluster Size vs #phi;#phi", 300, 0, 6.28319, 20, 0.5, 20.5); + hNClustersFromHit[layer][type] = new TH1F(Form("hNClsPerHit%sTrkLayer%d", trackName[type], layer), ";N Cluster per Hit;Counts", 21, -0.5, 20.5); + hMeanTrueClsSizeVsEta[layer][type] = new TH1F(Form("hMeanTrueClsSizeVsEta%sTrkLayer%d", trackName[type], layer), "Mean True Cluster Size vs #eta;#eta", 300, -2, 2); + hMeanTrueClsSizeVsPhi[layer][type] = new TH1F(Form("hMeanTrueClsSizeVsPhi%sTrkLayer%d", trackName[type], layer), "Mean True Cluster Size vs #phi;#phi", 300, 0, 6.28319); + hMeanFakeClsSizeVsEta[layer][type] = new TH1F(Form("hMeanFakeClsSizeVsEta%sTrkLayer%d", trackName[type], layer), "Mean Fake Cluster Size vs #eta;#eta", 300, -2, 2); + hMeanFakeClsSizeVsPhi[layer][type] = new TH1F(Form("hMeanFakeClsSizeVsPhi%sTrkLayer%d", trackName[type], layer), "Mean Fake Cluster Size vs #phi;#phi", 300, 0, 6.28319); + hTopoVsEta[layer][type] = new TH2F(Form("hClsSizeVsEtaTopo%sTrkLayer%d", trackName[type], layer), "Cluster Topology vs #eta;;#eta", kNTopologies, 0, kNTopologies, 20, -2, 2); + hClsSizeVsTopo[layer][type] = new TH2F(Form("hClsSizeVsTopo%sTrkLayer%d", trackName[type], layer), "Cluster Topology vs N Digits;;N Digits", kNTopologies, 0, kNTopologies, 20, 0.5, 20.5); + hXRes[layer][type] = new TH2F(Form("hDeltaXClsHit%sTrkLayer%d", trackName[type], layer), ";#Delta X;#eta", 1000, -0.2, 0.2, 20, -2, 2); + hYRes[layer][type] = new TH2F(Form("hDeltaYClsHit%sTrkLayer%d", trackName[type], layer), ";#Delta Y;#eta", 1000, -0.2, 0.2, 20, -2, 2); + hZRes[layer][type] = new TH2F(Form("hDeltaZClsHit%sTrkLayer%d", trackName[type], layer), ";#Delta Z;#eta", 1000, -0.2, 0.2, 20, -2, 2); + hRmsXRes[layer][type] = new TH1F(Form("hRmsDeltaXClsHit%sTrkLayer%d", trackName[type], layer), ";#eta;RMS #Delta X", 20, -2, 2); + hRmsYRes[layer][type] = new TH1F(Form("hRmsDeltaYClsHit%sTrkLayer%d", trackName[type], layer), ";#eta;RMS #Delta Y", 20, -2, 2); + hRmsZRes[layer][type] = new TH1F(Form("hRmsDeltaZClsHit%sTrkLayer%d", trackName[type], layer), ";#eta;RMS #Delta Z", 20, -2, 2); + hMeanXRes[layer][type] = new TH1F(Form("hMeanDeltaXClsHit%sTrkLayer%d", trackName[type], layer), ";#eta;Mean #Delta X", 20, -2, 2); + hMeanYRes[layer][type] = new TH1F(Form("hMeanDeltaYClsHit%sTrkLayer%d", trackName[type], layer), ";#eta;Mean #Delta Y", 20, -2, 2); + hMeanZRes[layer][type] = new TH1F(Form("hMeanDeltaZClsHit%sTrkLayer%d", trackName[type], layer), ";#eta;Mean #Delta Z", 20, -2, 2); - // LOOP on : clusters array - std::cout << "\n\n ----> Starting loop on clusters for ROF " << iROF << " with index " << rofIndex << " and nEntries " << rofNEntries << std::endl; - for (unsigned int iCls = rofIndex; iCls < rofIndex + rofNEntries; iCls++) { - if (iCls % 10000 == 0) { - std::cout << "Reading cluster " << iCls << " / " << clsArray->size() << std::endl; + if (layer == 0) { + hTrackHitsXY[layer][type] = new TH2F(Form("hTrackHitsXY%sTrkLayer%d", trackName[type], layer), ";Hit X;Hit Y", 5000, -30, 30, 5000, -30, 30); + hTrackDoubleHitsXY[layer][type] = new TH2F(Form("hTrackDoubleHitsXY%sTrkLayer%d", trackName[type], layer), ";Hit X;Hit Y", 5000, -30, 30, 5000, -30, 30); + } else { + hTrackHitsXY[layer][type] = new TH2F(Form("hTrackHitsXY%sTrkLayer%d", trackName[type], layer), ";Hit X;Hit Y", 10000, -100, 100, 10000, -100, 100); + hTrackDoubleHitsXY[layer][type] = new TH2F(Form("hTrackDoubleHitsXY%sTrkLayer%d", trackName[type], layer), ";Hit X;Hit Y", 10000, -100, 100, 10000, -100, 100); } + hTrackDoubleHitsPhiPt[layer][type] = new TH2F(Form("hTrackDoubleHitsPhiPt%sTrkLayer%d", trackName[type], layer), ";#phi;p_{T}", 3000, 0, 6.28319, 50, 0, 10); - Int_t iRow = (*clsArray)[iCls].row; - Int_t iCol = (*clsArray)[iCls].col; - Int_t chipID = (*clsArray)[iCls].chipID; - Int_t subDetID = tofGeo->getIOTOFLayer(chipID); - Float_t time = (*clsArray)[iCls].time; + for (int topo = 0; topo < kNTopologies; ++topo) { + hTopoSummaryReal->GetXaxis()->SetBinLabel(topo + 1, topologyNames[topo].c_str()); + hTopoSummaryFake->GetXaxis()->SetBinLabel(topo + 1, topologyNames[topo].c_str()); + hTopoSummaryTotal->GetXaxis()->SetBinLabel(topo + 1, topologyNames[topo].c_str()); + hTopoSummaryDictionary->GetXaxis()->SetBinLabel(topo + 1, topologyNames[topo].c_str()); + hTopoVsEta[layer][type]->GetXaxis()->SetBinLabel(topo + 1, topologyNames[topo].c_str()); + hClsSizeVsTopo[layer][type]->GetXaxis()->SetBinLabel(topo + 1, topologyNames[topo].c_str()); + hTopoVsEtaPt[layer][type][topo] = new TH2F(Form("h%sVsEtaPt_%sTrk_TrkLayer%d", topologyNames[topo].c_str(), trackName[type], layer), Form("Cluster Topology %s vs Eta and Pt;#eta;p_{T}", topologyNames[topo].c_str()), 100, -2, 2, 20, 0, 10); + } + } + } + + // Loop over clusters + Print(true, "----> Looping over clusters and filling histograms"); + for (const auto& cls : clustersProperties) { + + const int layer = cls.layer; + const int topo = static_cast(cls.topology); + + const int chipID = cls.chipID; + const int eventID = cls.eventID; + const int trackID = cls.trackID; + + const auto& mcTrack = (*mcTracksPerEvent[eventID])[trackID]; + const float eta = mcTrack.GetEta(); + const float phi = mcTrack.GetPhi(); + const float pt = mcTrack.GetPt(); + const int type = cls.isPrimary ? 0 : 1; + const int size = cls.size; + + hTopoVsEtaPt[layer][type][topo]->Fill(eta, pt); + hTopoVsEta[layer][type]->Fill(topo, eta); + + hClsSizeVsTopo[layer][type]->Fill(topo, size); + + hTopoSummaryTotal->Fill(topo); + if (cls.isFake) { + hTopoSummaryFake->Fill(topo); + hFakeClsSizeVsEta[layer][type]->Fill(eta, size); + hFakeClsSizeVsPhi[layer][type]->Fill(phi, size); + } else { + hTopoSummaryReal->Fill(topo); + hTrueClsSizeVsEta[layer][type]->Fill(eta, size); + hTrueClsSizeVsPhi[layer][type]->Fill(phi, size); + } - Float_t x = 0.f, y = 0.f, z = 0.f; - if (subDetID >= 0) { - segGeom->detectorToLocal(iRow, iCol, x, z, subDetID); + if (cls.hitIdx < 0) { + continue; // Skip clusters without a matching hit + } + const auto& hitData = allEvtsTrackData[cls.eventID][cls.trackID].hitsByDetector[cls.chipID][cls.hitIdx]; + auto& hit = (*hitsPerEvent[cls.eventID])[hitData.hitIdx]; + hNClustersFromHit[layer][type]->Fill(hitData.assocClsIdxs.size()); + if (hitData.assocClsIdxs.size() > 0) + hClustersEtaPhi[layer][type]->Fill(phi, eta); + + o2::math_utils::Point3D clusterPos; + TopologyInfo clsTopoInfo = topoClassifier.getTopologyFeatures(cls.topoKey); + auto clsFull = clustersArray->at(cls.clsIdx); + GetClusterGlobalPos(clsFull, clsTopoInfo, clusterPos, iotofGeom, segmInfo); + o2::math_utils::Point3D avgPos; + GetHitAvgPositionGlobal(hit, avgPos); + hXRes[layer][type]->Fill(clusterPos.X() - avgPos.X(), eta); + hYRes[layer][type]->Fill(clusterPos.Y() - avgPos.Y(), eta); + hZRes[layer][type]->Fill(clusterPos.Z() - avgPos.Z(), eta); + } + + // Fill means and RMS of cluster size and residuals + for (int layer = 0; layer < 2; ++layer) { + for (int type = 0; type < 2; ++type) { + for (int etaBin = 1; etaBin <= hTrueClsSizeVsEta[layer][type]->GetNbinsX(); ++etaBin) { + // Project 1D histogram to get mean cluster size for this eta bin + TH1D* hClsSizeProj = hTrueClsSizeVsEta[layer][type]->ProjectionY(Form("hClsSizeProj_etaBin%d", etaBin), etaBin, etaBin); + hMeanTrueClsSizeVsEta[layer][type]->SetBinContent(etaBin, hClsSizeProj->GetMean()); + hMeanTrueClsSizeVsEta[layer][type]->SetBinError(etaBin, hClsSizeProj->GetMeanError()); + } + for (int phiBin = 1; phiBin <= hTrueClsSizeVsPhi[layer][type]->GetNbinsX(); ++phiBin) { + // Project 1D histogram to get mean cluster size for this eta bin + TH1D* hClsSizeProj = hTrueClsSizeVsPhi[layer][type]->ProjectionY(Form("hClsSizeProj_phiBin%d", phiBin), phiBin, phiBin); + hMeanTrueClsSizeVsPhi[layer][type]->SetBinContent(phiBin, hClsSizeProj->GetMean()); + hMeanTrueClsSizeVsPhi[layer][type]->SetBinError(phiBin, hClsSizeProj->GetMeanError()); + } + for (int etaBin = 1; etaBin <= hFakeClsSizeVsEta[layer][type]->GetNbinsX(); ++etaBin) { + // Project 1D histogram to get mean cluster size for this eta bin + TH1D* hClsSizeProj = hFakeClsSizeVsEta[layer][type]->ProjectionY(Form("hClsSizeProj_etaBin%d", etaBin), etaBin, etaBin); + hMeanFakeClsSizeVsEta[layer][type]->SetBinContent(etaBin, hClsSizeProj->GetMean()); + hMeanFakeClsSizeVsEta[layer][type]->SetBinError(etaBin, hClsSizeProj->GetMeanError()); + } + for (int phiBin = 1; phiBin <= hFakeClsSizeVsPhi[layer][type]->GetNbinsX(); ++phiBin) { + // Project 1D histogram to get mean cluster size for this eta bin + TH1D* hClsSizeProj = hFakeClsSizeVsPhi[layer][type]->ProjectionY(Form("hClsSizeProj_phiBin%d", phiBin), phiBin, phiBin); + hMeanFakeClsSizeVsPhi[layer][type]->SetBinContent(phiBin, hClsSizeProj->GetMean()); + hMeanFakeClsSizeVsPhi[layer][type]->SetBinError(phiBin, hClsSizeProj->GetMeanError()); } + for (int etaBin = 1; etaBin <= hXRes[layer][type]->GetNbinsY(); ++etaBin) { + TH1D* hXResProj = hXRes[layer][type]->ProjectionX(Form("hXResProj_etaBin%d", etaBin), etaBin, etaBin); + TH1D* hYResProj = hYRes[layer][type]->ProjectionX(Form("hYResProj_etaBin%d", etaBin), etaBin, etaBin); + TH1D* hZResProj = hZRes[layer][type]->ProjectionX(Form("hZResProj_etaBin%d", etaBin), etaBin, etaBin); + hRmsXRes[layer][type]->SetBinContent(etaBin, hXResProj->GetRMS()); + hRmsYRes[layer][type]->SetBinContent(etaBin, hYResProj->GetRMS()); + hRmsZRes[layer][type]->SetBinContent(etaBin, hZResProj->GetRMS()); + hRmsXRes[layer][type]->SetBinError(etaBin, hXResProj->GetRMSError()); + hRmsYRes[layer][type]->SetBinError(etaBin, hYResProj->GetRMSError()); + hRmsZRes[layer][type]->SetBinError(etaBin, hZResProj->GetRMSError()); + hMeanXRes[layer][type]->SetBinContent(etaBin, hXResProj->GetMean()); + hMeanYRes[layer][type]->SetBinContent(etaBin, hYResProj->GetMean()); + hMeanZRes[layer][type]->SetBinContent(etaBin, hZResProj->GetMean()); + hMeanXRes[layer][type]->SetBinError(etaBin, hXResProj->GetMeanError()); + hMeanYRes[layer][type]->SetBinError(etaBin, hYResProj->GetMeanError()); + hMeanZRes[layer][type]->SetBinError(etaBin, hZResProj->GetMeanError()); + } + } + } + + Print(true, "----> Looping over generated particles"); - o2::math_utils::Point3D localClsCoords(x, y, z); // local Digit - const auto globalClsCoords = tofGeo->getMatrixL2G(chipID)(localClsCoords); // convert to global - clsTuple->Fill((*clsArray)[iCls].chipID, - globalClsCoords.x(), - globalClsCoords.y(), - globalClsCoords.z(), - (*clsArray)[iCls].row, - (*clsArray)[iCls].col, - (*clsArray)[iCls].time); - histXCoordCls->Fill(globalClsCoords.x()); - histYCoordCls->Fill(globalClsCoords.y()); - histZCoordCls->Fill(globalClsCoords.z()); + // Generated particles + TH2F* hGenEtaPt[2] = {new TH2F("hGenEtaPtPrm", "Generated primary tracks;#eta;p_{T}", 100, -2, 2, 100, 0, 10), + new TH2F("hGenEtaPtSec", "Generated secondary tracks;#eta;p_{T}", 100, -2, 2, 100, 0, 10)}; + + for (int iEvt = 0; iEvt < nEvts; ++iEvt) { + for (const auto& mcTrack : *mcTracksPerEvent[iEvt]) { + const int type = mcTrack.isPrimary() ? 0 : 1; + hGenEtaPt[type]->Fill(mcTrack.GetEta(), mcTrack.GetPt()); + } + } - // Match to digit - auto digitLabelFromCls = (clsLabels->getLabels(iCls))[0]; - auto digitEntry = digitsLabels.find(digitLabelFromCls); + // Check eta and phi of tracks producing multiple hits, should reflect + // overlaps between staves and validate the geometry implementation + Print(true, "----> Looping over tracks producing multiple hits"); + for (int iEvt = 0; iEvt < nEvts; ++iEvt) { + for (const auto& [trackID, trackData] : allEvtsTrackData[iEvt]) { - if (digitEntry == digitsLabels.end()) { - LOG(error) << "No matching digit for cluster " << iCls << " with label " << digitLabelFromCls.getRawValue(); + const auto& mcTrack = (*mcTracksPerEvent[iEvt])[trackID]; + if (!mcTrack.isPrimary() || trackData.hitsByDetector.size() <= 1) { continue; } - int iDigit = digitEntry->second; - Int_t iRowFromDigit = (*digitsArray)[iDigit].getRow(); - Int_t iColFromDigit = (*digitsArray)[iDigit].getColumn(); - Int_t iChipIDFromDigit = (*digitsArray)[iDigit].getChipIndex(); - Int_t iSubDetIDFromDigit = tofGeo->getIOTOFLayer(iChipIDFromDigit); - Float_t timeFromDigit = (*digitsArray)[iDigit].getTime(); - - float xFromDigit = 0.f, yFromDigit = 0.f, zFromDigit = 0.f; - if (iSubDetIDFromDigit >= 0) { - segGeom->detectorToLocal(iRowFromDigit, iColFromDigit, xFromDigit, zFromDigit, iSubDetIDFromDigit); - } - - o2::math_utils::Point3D localDigitCoordFromDigit(xFromDigit, yFromDigit, zFromDigit); // local Digit - const auto globalDigitCoordFromDigit = tofGeo->getMatrixL2G(iChipIDFromDigit)(localDigitCoordFromDigit); // convert to global - histXCoordRes->Fill(globalClsCoords.x() - globalDigitCoordFromDigit.X()); - histYCoordRes->Fill(globalClsCoords.y() - globalDigitCoordFromDigit.Y()); - histZCoordRes->Fill(globalClsCoords.z() - globalDigitCoordFromDigit.Z()); - histTimeRes->Fill(time - timeFromDigit); - } // end loop on clusters array - } // end loop on ROFRecords - - std::cout << "Cluster array size: " << clsTuple->GetEntries() << std::endl; - - // cluster maps in the xy and yz planes - auto canvXY = new TCanvas("canvXY", "", 1600, 800); - canvXY->Divide(2, 1); - canvXY->cd(1); - clsTuple->Draw("y:x>>h_y_vs_x_IOTOF(1000, -100, 100, 1000, -100, 100)", "", "colz"); - canvXY->cd(2); - clsTuple->Draw("y:z>>h_y_vs_z_IOTOF(1000, -400, 400, 1000, -100, 100)", "", "colz"); - canvXY->SaveAs("clusters_digits_y_vs_x_vs_z.pdf"); - - // z distributions - auto canvZ = new TCanvas("canvZ", "", 800, 800); - canvZ->cd(); - clsTuple->Draw("z>>h_z_IOTOF(500, -70, 70)", ""); - canvZ->SaveAs("clusters_digits_z.pdf"); + // Index 0 -> Layer 0, Index 1 -> Layer 1 + std::vector distinctChips[2]; + for (const auto& [chipIdx, hitsVec] : trackData.hitsByDetector) { + + int layer{-1}, stave{-1}, subStave{-1}, module{-1}, chip{-1}; + iotofGeom->getIOTOFChipId(chipIdx, layer, stave, subStave, module, chip); + + // Check if current chip is a neighbor to any already accepted chip in this layer + // Required because the same track can produce multiple hits in adjacent chips, + // belonging to the same module/substave, therefore the double hit is not related + // to the detector geometry + const bool isNeighborToExisting = std::any_of( + distinctChips[layer].begin(), + distinctChips[layer].end(), + [&](int existingChipIdx) { + int layerA{-1}, staveA{-1}, subStaveA{-1}, moduleA{-1}, chipA{-1}; + iotofGeom->getIOTOFChipId(existingChipIdx, layerA, staveA, subStaveA, moduleA, chipA); + + // Reject adjacent modules in the same stave, substave + if (layer == layerA && stave == staveA && subStave == subStaveA && std::abs(module - moduleA) <= 1) { + return true; + } + // Reject adjacent chips with same stave, subStave, module but different chip index + if (layer == layerA && stave == staveA && subStave == subStaveA && module == moduleA &&std::abs(chip - chipA) <= 1) { + return true; + } + return false; + } + ); + + // Keep chip ONLY IF it is not an immediate neighbor to an existing one + if (!isNeighborToExisting) { + distinctChips[layer].push_back(chipIdx); + } + } + + // Fill histograms with properties of tracks producing multiple hits + for (int layer = 0; layer < 2; ++layer) { + for (const auto& [chipIdx, hitsVec] : trackData.hitsByDetector) { + if (iotofGeom->getIOTOFLayer(chipIdx) != layer) { + continue; + } + + for (const auto& hitData : hitsVec) { + if (hitData.hitIdx < 0) { + continue; // Skip if no matching hit + } + const auto& hit = (*hitsPerEvent[iEvt])[hitData.hitIdx]; + PrintHit(verbose, hit, iotofGeom); + + const int type = mcTrack.isPrimary() ? 0 : 1; + hTrackHitsXY[layer][type]->Fill(hit.GetX(), hit.GetY()); + } + } + } + + // Fill histograms with properties of tracks producing multiple hits + for (int layer = 0; layer < 2; ++layer) { + if (distinctChips[layer].size() <= 1) { + continue; + } + + for (const auto& [chipIdx, hitsVec] : trackData.hitsByDetector) { + if (iotofGeom->getIOTOFLayer(chipIdx) != layer) { + continue; + } + + for (const auto& hitData : hitsVec) { + if (hitData.hitIdx < 0) { + continue; // Skip if no matching hit + } + const auto& hit = (*hitsPerEvent[iEvt])[hitData.hitIdx]; + PrintHit(verbose, hit, iotofGeom); + + const int type = mcTrack.isPrimary() ? 0 : 1; + if (mcTrack.GetPt() > 5.0f) { + hTrackDoubleHitsXY[layer][type]->Fill(hit.GetX(), hit.GetY()); + } + hTrackDoubleHitsPhiPt[layer][type]->Fill(mcTrack.GetPhi(), mcTrack.GetPt()); + } + } + } + } + } + + Print(true, "----> Writing histograms"); + // Output TFile* outFile = new TFile("CheckClusters.root", "RECREATE"); - // Save all columns of the tuple as hists - clsTuple->Write(); - histXCoordCls->Write(); - histYCoordCls->Write(); - histZCoordCls->Write(); - histXCoordDigit->Write(); - histYCoordDigit->Write(); - histZCoordDigit->Write(); - histXCoordRes->Write(); - histYCoordRes->Write(); - histZCoordRes->Write(); - histTimeRes->Write(); - outFile->Write(); + for (int type = 0; type < 2; ++type) { + hGenEtaPt[type]->Write(); + } + + hEtaPhiHitsPrmTrkLayer0->Write(); + hEtaPhiHitsSecTrkLayer0->Write(); + hEtaPhiHitsPrmTrkLayer1->Write(); + hEtaPhiHitsSecTrkLayer1->Write(); + hTopoSummaryReal->Write(); + hTopoSummaryFake->Write(); + hTopoSummaryTotal->Write(); + hTopoSummaryDictionary->Write(); + hCountHitMatchingType->Write(); + + for (int layer = 0; layer < 2; ++layer) { + + for (int type = 0; type < 2; ++type) { + outFile->mkdir(Form("%sTrkLayer%d", trackName[type], layer)); + outFile->mkdir(Form("%sTrkLayer%d/Topologies", trackName[type], layer)); + outFile->cd(Form("%sTrkLayer%d", trackName[type], layer)); + + hCountFakeClusters[layer][type]->Write("hCountFakeClusters"); + + hClustersEtaPhi[layer][type]->Write("hClustersEtaPhi"); + hTrueClsSizeVsEta[layer][type]->Write("hTrueClsSizeVsEta"); + hTrueClsSizeVsPhi[layer][type]->Write("hTrueClsSizeVsPhi"); + hFakeClsSizeVsEta[layer][type]->Write("hFakeClsSizeVsEta"); + hFakeClsSizeVsPhi[layer][type]->Write("hFakeClsSizeVsPhi"); + + TH2F* hEfficiency = static_cast(hClustersEtaPhi[layer][type]->Clone(Form("hClusterEfficiencyVsEtaPhi%sTrkLayer%d", trackName[type], layer))); + TH2F* hHits = layer == 0 ? (type == 0 ? hEtaPhiHitsPrmTrkLayer0 : hEtaPhiHitsSecTrkLayer0) + : (type == 0 ? hEtaPhiHitsPrmTrkLayer1 : hEtaPhiHitsSecTrkLayer1); + hEfficiency->Divide(hHits); + hEfficiency->Write("hClsEfficiency"); + delete hEfficiency; + + hNClustersFromHit[layer][type]->Write("hNClustersFromHit"); + hClsSizeVsTopo[layer][type]->Write("hClsSizeVsTopo"); + hMeanTrueClsSizeVsEta[layer][type]->Write("hMeanTrueClsSizeVsEta"); + hMeanTrueClsSizeVsPhi[layer][type]->Write("hMeanTrueClsSizeVsPhi"); + hMeanFakeClsSizeVsEta[layer][type]->Write("hMeanFakeClsSizeVsEta"); + hMeanFakeClsSizeVsPhi[layer][type]->Write("hMeanFakeClsSizeVsPhi"); + hTopoVsEta[layer][type]->Write("hTopoVsEta"); + hXRes[layer][type]->Write("hXRes"); + hYRes[layer][type]->Write("hYRes"); + hZRes[layer][type]->Write("hZRes"); + hRmsXRes[layer][type]->Write("hRmsXRes"); + hRmsYRes[layer][type]->Write("hRmsYRes"); + hRmsZRes[layer][type]->Write("hRmsZRes"); + hMeanXRes[layer][type]->Write("hMeanXRes"); + hMeanYRes[layer][type]->Write("hMeanYRes"); + hMeanZRes[layer][type]->Write("hMeanZRes"); + + if (type == 0) { + hTrackHitsXY[layer][type]->Write("hTrackHitsXY"); + hTrackDoubleHitsXY[layer][type]->Write("hTrackDoubleHitsXY"); + hTrackDoubleHitsPhiPt[layer][type]->Write("hTrackDoubleHitsPhiPt"); + } + + outFile->cd(Form("%sTrkLayer%d/Topologies", trackName[type], layer)); + for (int topo = 0; topo < kNTopologies; ++topo) hTopoVsEtaPt[layer][type][topo]->Write(Form("%sVsEtaPt", topologyNames[topo].c_str())); + } + } + + // Create canvas overlapping hTrackHitsXY and hTrackDoubleHitsXY with + // different colors in a restricted range to visualize the double hits + + TCanvas* cTrackHitsXY[2][2]; + TCanvas* cTrackHitsXYZoom[2][2]; + for (int layer = 0; layer < 2; ++layer) { + for (int type = 0; type < 2; ++type) { + cTrackHitsXY[layer][type] = new TCanvas( + Form("cTrackHitsXY%sTrkLayer%d", trackName[type], layer), + Form("Track Hits XY %s Track Layer %d", trackName[type], layer), + 800, 600 + ); + + // Constrain in a box (xMin, xMax, yMin, yMax) to visualize the double hits + if (layer == 0) { + hTrackHitsXY[layer][type]->GetXaxis()->SetRangeUser(-22, 0); + hTrackHitsXY[layer][type]->GetYaxis()->SetRangeUser(-22, 0); + hTrackDoubleHitsXY[layer][type]->GetXaxis()->SetRangeUser(-22, 0); + hTrackDoubleHitsXY[layer][type]->GetYaxis()->SetRangeUser(-22, 0); + } else { + hTrackHitsXY[layer][type]->GetXaxis()->SetRangeUser(-50, -20); + hTrackHitsXY[layer][type]->GetYaxis()->SetRangeUser(-95, -75); + hTrackDoubleHitsXY[layer][type]->GetXaxis()->SetRangeUser(-50, -20); + hTrackDoubleHitsXY[layer][type]->GetYaxis()->SetRangeUser(-95, -75); + } + + // First histogram: normal track hits + hTrackHitsXY[layer][type]->SetLineColor(kBlue); + hTrackHitsXY[layer][type]->SetLineWidth(2); + hTrackHitsXY[layer][type]->SetFillStyle(0); + + // Draw only the histogram contours. + hTrackHitsXY[layer][type]->Draw("CONT3"); + + // Second histogram: double hits + hTrackDoubleHitsXY[layer][type]->SetLineColor(kRed); + hTrackDoubleHitsXY[layer][type]->SetLineWidth(2); + hTrackDoubleHitsXY[layer][type]->SetFillStyle(0); + + // Overlay the double-hit contours. + hTrackDoubleHitsXY[layer][type]->Draw("CONT3 SAME"); + + // Don't save stats panel + gStyle->SetOptStat(0); + + // Save + cTrackHitsXY[layer][type]->Write(); + cTrackHitsXY[layer][type]->SaveAs(Form("cTrackHitsXY%sTrkLayer%d.pdf", trackName[type], layer)); + } + } + + + // Check digit efficiency across pixel by print the local coordinates + // of hits without any cluster and digit associated to them + Print(true, "----> Checking digit efficiency across pixel"); + TH2F* hNotRecoHits[2][2]; + for (int layer = 0; layer < 2; ++layer) { + for (int type = 0; type < 2; ++type) { + hNotRecoHits[layer][type] = new TH2F(Form("hNotRecoHits%sTrkLayer%d", trackName[type], layer), "Hits with no clusters or digits", 6000, -3, 3, 600, 3, 3); + } + } + for (int iEvt = 0; iEvt < nEvts; ++iEvt) { + for (const auto& [trackID, trackData] : allEvtsTrackData[iEvt]) { + const auto& mcTrack = (*mcTracksPerEvent[iEvt])[trackID]; + const int type = mcTrack.isPrimary() ? 0 : 1; + + for (const auto& [chipIdx, hitsVec] : trackData.hitsByDetector) { + int layer{-1}, stave{-1}, subStave{-1}, module{-1}, chip{-1}; + iotofGeom->getIOTOFChipId(chipIdx, layer, stave, subStave, module, chip); + + for (const auto& hitData : hitsVec) { + if (hitData.hitIdx < 0) { + continue; // Skip if no matching hit + } + const auto& hit = (*hitsPerEvent[iEvt])[hitData.hitIdx]; + if (hitData.assocClsIdxs.empty() && hitData.assocDigitIdxs.empty()) { + Print(verbose, "Hit with no associated clusters or digits:"); + o2::math_utils::Point3D avgPos; + GetHitAvgPositionLocal(hit, iotofGeom, avgPos); + Print(verbose, Form("Local position: x = %.5f, y = %.5f, z = %.5f", avgPos.X(), avgPos.Y(), avgPos.Z())); + hNotRecoHits[layer][type]->Fill(avgPos.X(), avgPos.Y()); + } + } + } + } + } + // Write digit efficiency histograms + for (int layer = 0; layer < 2; ++layer) { + for (int type = 0; type < 2; ++type) { + outFile->cd(Form("%sTrkLayer%d", trackName[type], layer)); + hNotRecoHits[layer][type]->Write(); + } + } + outFile->Close(); + delete outFile; + + + // // Print all properties of fake clusters + // for (const auto& cluster : clusters) { + // if (cluster.isFakeDiffHits || cluster.isFakeDiffTrks || cluster.isFakeDiffEvts) { + // std::cout << "\n\n\nFake cluster properties: " << std::endl; + // PrintCluster(true, cluster, digitsArray, digitsLabels, hitsPerEvent, iotofGeom, segmInfo); + // } + // } + } diff --git a/Detectors/Upgrades/ALICE3/IOTOF/macros/CheckDigitsIOTOF.C b/Detectors/Upgrades/ALICE3/IOTOF/macros/CheckDigitsIOTOF.C index 26ffd08697d56..d6713a8ee0eee 100644 --- a/Detectors/Upgrades/ALICE3/IOTOF/macros/CheckDigitsIOTOF.C +++ b/Detectors/Upgrades/ALICE3/IOTOF/macros/CheckDigitsIOTOF.C @@ -22,7 +22,7 @@ #include #include -#include "IOTOFSimulation/Segmentation.h" +#include "IOTOFBase/Segmentation.h" #include "IOTOFBase/IOTOFBaseParam.h" #include "IOTOFBase/GeometryTGeo.h" #include "DataFormatsIOTOF/Digit.h" @@ -75,7 +75,10 @@ void addTLines(float pitch) gPad->Update(); } -void CheckDigitsIOTOF(std::string digifile = "tf3digits.root", std::string hitfile = "o2sim_HitsTF3.root", std::string inputGeom = "o2sim_geometry.root") +void CheckDigitsIOTOF(std::string digifile = "tf3digits.root", + std::string hitfile = "o2sim_HitsTF3.root", + std::string inputGeom = "o2sim_geometry.root", + std::string geomCfgStr = "IOTOFBase.segmentedInnerTOF=true;IOTOFBase.segmentedOuterTOF=true;IOTOFBase.enableForwardTOF=false;IOTOFBase.enableBackwardTOF=false;") { gStyle->SetPalette(55); @@ -85,7 +88,7 @@ void CheckDigitsIOTOF(std::string digifile = "tf3digits.root", std::string hitfi using o2::iotof::Digit; using o2::itsmft::Hit; - o2::conf::ConfigurableParam::updateFromString("IOTOFBase.segmentedInnerTOF=true;IOTOFBase.segmentedOuterTOF=true;IOTOFBase.enableForwardTOF=false;IOTOFBase.enableBackwardTOF=false"); + o2::conf::ConfigurableParam::updateFromString(geomCfgStr); auto seg = o2::iotof::Segmentation::Instance(); @@ -100,6 +103,8 @@ void CheckDigitsIOTOF(std::string digifile = "tf3digits.root", std::string hitfi o2::base::GeometryManager::loadGeometry(inputGeom); auto* gman = o2::iotof::GeometryTGeo::Instance(); gman->fillMatrixCache(o2::math_utils::bit2Mask(o2::math_utils::TransformType::L2G)); + std::cout << "Number of chips in ITOF: " << gman->getITOFNumberOfChips() << std::endl; + std::cout << "Number of chips in OTOF: " << gman->getOTOFNumberOfChips() << std::endl; // Hits TFile* hitFile = TFile::Open(hitfile.data()); @@ -250,13 +255,13 @@ void CheckDigitsIOTOF(std::string digifile = "tf3digits.root", std::string hitfi auto canvdXdZ = new TCanvas("canvdXdZ", "", 1600, 800); canvdXdZ->Divide(2, 1); canvdXdZ->cd(1); - nt->Draw("dx:dz>>h_dx_vs_dz_ITOF(600, -0.03, 0.03, 600, -0.03, 0.03)", "id >= 0 && id < 1920", "colz"); + nt->Draw("dx:dz>>h_dx_vs_dz_ITOF(1000, -0.05, 0.05, 1000, -0.05, 0.05)", "id >= 0 && id < 1920", "colz"); addTLines(0.01); auto h = (TH2F*)gPad->GetPrimitive("h_dx_vs_dz_ITOF"); Info("ITOF", "RMS(dx)=%.1f mu", h->GetRMS(2) * 1e4); Info("ITOF", "RMS(dz)=%.1f mu", h->GetRMS(1) * 1e4); canvdXdZ->cd(2); - nt->Draw("dx:dz>>h_dx_vs_dz_OTOF(600, -0.03, 0.03, 600, -0.03, 0.03)", "id >= 1920 && id < 53568", "colz"); + nt->Draw("dx:dz>>h_dx_vs_dz_OTOF(1000, -0.05, 0.05, 1000, -0.05, 0.05)", "id >= 1920 && id < 53568", "colz"); addTLines(0.01); h = (TH2F*)gPad->GetPrimitive("h_dx_vs_dz_OTOF"); Info("OTOF", "RMS(dx)=%.1f mu", h->GetRMS(2) * 1e4); diff --git a/Detectors/Upgrades/ALICE3/IOTOF/reconstruction/CMakeLists.txt b/Detectors/Upgrades/ALICE3/IOTOF/reconstruction/CMakeLists.txt index 9a887bff8127c..96979eab3b2f1 100644 --- a/Detectors/Upgrades/ALICE3/IOTOF/reconstruction/CMakeLists.txt +++ b/Detectors/Upgrades/ALICE3/IOTOF/reconstruction/CMakeLists.txt @@ -12,9 +12,19 @@ o2_add_library(IOTOFReconstruction TARGETVARNAME targetName SOURCES src/Clusterer.cxx + src/ClustererParam.cxx + src/TopologyClassifier.cxx PUBLIC_LINK_LIBRARIES Microsoft.GSL::GSL O2::DataFormatsIOTOF O2::IOTOFBase O2::IOTOFSimulation + O2::FrameworkLogger ) + +o2_target_root_dictionary( + IOTOFReconstruction + HEADERS include/IOTOFReconstruction/Clusterer.h + include/IOTOFReconstruction/ClustererParam.h + include/IOTOFReconstruction/TopologyClassifier.h + ) diff --git a/Detectors/Upgrades/ALICE3/IOTOF/reconstruction/include/IOTOFReconstruction/Clusterer.h b/Detectors/Upgrades/ALICE3/IOTOF/reconstruction/include/IOTOFReconstruction/Clusterer.h index 252ecf8917377..4b595145506fd 100644 --- a/Detectors/Upgrades/ALICE3/IOTOF/reconstruction/include/IOTOFReconstruction/Clusterer.h +++ b/Detectors/Upgrades/ALICE3/IOTOF/reconstruction/include/IOTOFReconstruction/Clusterer.h @@ -18,6 +18,9 @@ #include "DataFormatsIOTOF/Digit.h" #include "DataFormatsITSMFT/ROFRecord.h" #include "DataFormatsIOTOF/Cluster.h" +#include "IOTOFSimulation/DPLDigitizerParam.h" +#include "IOTOFReconstruction/ClustererParam.h" +#include "IOTOFReconstruction/TopologyClassifier.h" #include "SimulationDataFormat/ConstMCTruthContainer.h" #include "SimulationDataFormat/MCCompLabel.h" #include "SimulationDataFormat/MCTruthContainer.h" @@ -47,28 +50,32 @@ class Clusterer //---------------------------------------------- struct ClustererThread { - Clusterer* parent = nullptr; + Clusterer* mParent = nullptr; // Column buffers data members in TRK, for now not needed in TF3 // Further struct members in TRK, for now not needed in TF3 - std::array labelsBuff; ///< MC label buffer for one cluster + std::array mLabelsBuff; ///< MC label buffer for one cluster // per-thread output (accumulated, then merged back by caller) - std::vector clusters; - std::vector patterns; - ClusterTruth labels; + std::vector mClusters; + std::vector mPatterns; + ClusterTruth mLabels; // Further reset column buffer in TRK, not included for now in TF3 + TopologyClassifier mClsTopoClassifier; //! Convert the cluster topology to the corresponding entry in the dictionary. void fetchMCLabels(uint32_t digID, const ConstDigitTruth* labelsDig, int& nfilled); - void finishChipSingleHitFast(gsl::span digits, uint32_t digitIdx, - const ConstDigitTruth* labelsDigPtr, ClusterTruth* labelsClusPtr); + void findClustersSingleHit(gsl::span digits, uint32_t digitIdx, + const ConstDigitTruth* labelsDigPtr, ClusterTruth* labelsClusPtr); + void findClustersMultipleHits(gsl::span digits, gsl::span digitIdxs, + const ConstDigitTruth* labelsDigPtr, ClusterTruth* labelsClusPtr); void processChip(gsl::span digits, int chipFirst, int chipN, std::vector* clustersOut, std::vector* patternsOut, const ConstDigitTruth* labelsDigPtr, ClusterTruth* labelsClusPtr); + void writeTopologiesToFile(const char* filename); - explicit ClustererThread(Clusterer* par = nullptr) : parent(par) {} + explicit ClustererThread(Clusterer* par = nullptr) : mParent(par) {} ClustererThread(const ClustererThread&) = delete; ClustererThread& operator=(const ClustererThread&) = delete; }; @@ -84,6 +91,12 @@ class Clusterer gsl::span digMC2ROFs = {}, std::vector* clusterMC2ROFs = nullptr); + // ///< load the dictionary of cluster topologies + // void loadDictionary(const std::string& fileName) { mPattIdConverter.loadDictionary(fileName); } + // void setDictionary(const TopologyDictionary* dict) { mPattIdConverter.setDictionary(dict); } + // const TopologyDictionary& getDictionary() const { return mPattIdConverter.getDictionary(); } + // auto& getPattIdConverter() const { return mPattIdConverter; } + protected: std::unique_ptr mThread; std::vector mSortIdx; ///< reusable per-ROF sort buffer diff --git a/Detectors/Upgrades/ALICE3/IOTOF/reconstruction/include/IOTOFReconstruction/ClustererParam.h b/Detectors/Upgrades/ALICE3/IOTOF/reconstruction/include/IOTOFReconstruction/ClustererParam.h new file mode 100644 index 0000000000000..038cf639ba674 --- /dev/null +++ b/Detectors/Upgrades/ALICE3/IOTOF/reconstruction/include/IOTOFReconstruction/ClustererParam.h @@ -0,0 +1,43 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +/// \file ClustererParam.h +/// \brief Definition of the IOTOF clusterer settings + +#ifndef ALICEO2_IOTOFCLUSTERERPARAM_H_ +#define ALICEO2_IOTOFCLUSTERERPARAM_H_ + +#include "DetectorsCommonDataFormats/DetID.h" +#include "CommonUtils/ConfigurableParam.h" +#include "CommonUtils/ConfigurableParamHelper.h" +#include +#include + +// TO BE REMOVED BEFORE PUSH +#include "Framework/Logger.h" + +namespace o2 +{ +namespace iotof +{ +struct ClustererParam : public o2::conf::ConfigurableParamHelper { + + int maxTimeDiffNSigma = 3; ///< maximum time difference in nsigma for clustering + int maxFiredDigitsForCls = 16; ///< maximum time difference in nsigma for clustering + + // boilerplate stuff + make principal key + O2ParamDef(ClustererParam, "TF3ClustererParam"); +}; + +} // namespace iotof +} // namespace o2 + +#endif diff --git a/Detectors/Upgrades/ALICE3/IOTOF/reconstruction/include/IOTOFReconstruction/TopologyClassifier.h b/Detectors/Upgrades/ALICE3/IOTOF/reconstruction/include/IOTOFReconstruction/TopologyClassifier.h new file mode 100644 index 0000000000000..4d3c0162b9982 --- /dev/null +++ b/Detectors/Upgrades/ALICE3/IOTOF/reconstruction/include/IOTOFReconstruction/TopologyClassifier.h @@ -0,0 +1,106 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +/// \file TopologyClassifier.h +/// \brief Definition of the TopologyClassifier class. +/// +/// Short TopologyClassifier descritpion +/// +/// This class is for the association of the cluster +/// topology with the corresponding entry in the dictionary +/// + +#ifndef ALICEO2_IOTOF_TOPOLOGYCLASSIFIER_H +#define ALICEO2_IOTOF_TOPOLOGYCLASSIFIER_H + +#include +#include +#include + +#include + +// TO BE REMOVED BEFORE PUSH +#include "Framework/Logger.h" + +namespace o2 +{ +namespace iotof +{ + +enum Topologies : uint8_t { + kSingleDigit, + kLineOnRow, + kLineOnCol, + kSquare, + kRectangle, + kDiagonal, + kLowerTriangleLeft, + kLowerTriangleRight, + kUpperTriangleLeft, + kUpperTriangleRight, + kSnake, + kSnakeRefl, + kSnakeRot90, + kSnakeRot90Refl, + kHuge, + kOther, + kNTopologies +}; + +struct TopologyInfo { + int mSizeX = 0; + int mSizeZ = 0; + int mOffsetXToCOG = 0; + int mOffsetZToCOG = 0; + float mXMean = 0.f; + float mZMean = 0.f; + float mXSigma2 = 0.f; + float mZSigma2 = 0.f; + int mNPixels = 0; + int mFrequency = 0; + Topologies mTopology = Topologies::kNTopologies; + uint16_t mPattern; ///< Bitmask of fired pixels +}; + +class TopologyClassifier { + public: + // Define limits for domain validation + static constexpr uint8_t MaxRowSpan = 255; + static constexpr uint8_t MaxColSpan = 255; + static constexpr uint16_t MaxBitmask = 65535; + + TopologyClassifier() = default; + TopologyClassifier(std::unordered_map map) : mTopologyCache(std::move(map)) {} + + const std::unordered_map& getTopologyMap() const { return mTopologyCache; }; + void getTopology(uint16_t bitmask, uint16_t minRow, uint8_t spanRow, uint16_t minCol, uint8_t spanCol, uint8_t& topology); + TopologyInfo getTopologyFeatures(uint32_t key); + void accountTopology(uint16_t bitmask, uint16_t minRow, uint8_t spanRow, uint16_t minCol, uint8_t spanCol, uint8_t& topology); + void computeCOG(uint16_t bitmask, uint16_t minRow, uint8_t spanRow, uint16_t minCol, uint8_t spanCol, TopologyInfo& topoInfo); + + void saveCacheToFile(const char* filename); + void print(); + + private: + /// Packs: [ spanRow (8b) ][ spanCol (8b) ][ bitmask (16b) ] -> 32 bits total + [[nodiscard]] static constexpr uint32_t packKey(uint8_t spanRow, uint8_t spanCol, uint16_t bitmask) noexcept { + return (static_cast(spanRow) << 24) | + (static_cast(spanCol) << 16) | + static_cast(bitmask); + } + + std::unordered_map mTopologyCache; +}; + +} // namespace iotof +} // namespace o2 + +#endif // ALICEO2_IOTOF_TOPOLOGYCLASSIFIER_H diff --git a/Detectors/Upgrades/ALICE3/IOTOF/reconstruction/src/Clusterer.cxx b/Detectors/Upgrades/ALICE3/IOTOF/reconstruction/src/Clusterer.cxx index edb9f71ac7f04..7f1c93672bcac 100644 --- a/Detectors/Upgrades/ALICE3/IOTOF/reconstruction/src/Clusterer.cxx +++ b/Detectors/Upgrades/ALICE3/IOTOF/reconstruction/src/Clusterer.cxx @@ -40,21 +40,21 @@ void Clusterer::process(gsl::span digits, } for (size_t iROF = 0; iROF < digitROFs.size(); ++iROF) { - LOG(debug) << "Processing digit ROF " << iROF << "/" << digitROFs.size(); - const auto& inROF = digitROFs[iROF]; - const auto outFirst = static_cast(clusters.size()); - const int first = inROF.getFirstEntry(); - const int nEntries = inROF.getNEntries(); - - if (nEntries == 0) { - LOG(debug) << "Digit ROF " << iROF << " has no entries, skipping"; - clusterROFs.emplace_back(inROF.getBCData(), inROF.getROFrame(), outFirst, 0); + LOG(info) << "[Clusterer] Processing digit ROF " << iROF << "/" << digitROFs.size(); + const auto& digitsThisROF = digitROFs[iROF]; + const auto nStoredCls = static_cast(clusters.size()); + const int first = digitsThisROF.getFirstEntry(); + const int nDigits = digitsThisROF.getNEntries(); + + if (nDigits == 0) { + LOG(info) << "[Clusterer] Digit ROF " << iROF << " has no entries, skipping"; + clusterROFs.emplace_back(digitsThisROF.getBCData(), digitsThisROF.getROFrame(), nStoredCls, 0); continue; } - // Sort digit indices within this ROF by (chipID, col, row) - // chip by chip, column by column (taken from TRK). - mSortIdx.resize(nEntries); + // Sort digit indices within this ROF by (chipID, row, col, time) + // extended with time information from TRK. + mSortIdx.resize(nDigits); std::iota(mSortIdx.begin(), mSortIdx.end(), first); std::sort(mSortIdx.begin(), mSortIdx.end(), [&digits](int a, int b) { const auto& da = digits[a]; @@ -62,30 +62,35 @@ void Clusterer::process(gsl::span digits, if (da.getChipIndex() != db.getChipIndex()) { return da.getChipIndex() < db.getChipIndex(); } + if (da.getRow() != db.getRow()) { + return da.getRow() < db.getRow(); + } if (da.getColumn() != db.getColumn()) { return da.getColumn() < db.getColumn(); } - return da.getRow() < db.getRow(); + return da.getTime() < db.getTime(); }); - LOG(debug) << "Found " << nEntries << " digits for ROF " << iROF; - - // Process blocks of chips with the same chipID - int sliceStart = 0; - while (sliceStart < nEntries) { - const int chipFirst = sliceStart; - const uint16_t chipID = digits[mSortIdx[sliceStart]].getChipIndex(); - while (sliceStart < nEntries && digits[mSortIdx[sliceStart]].getChipIndex() == chipID) { - ++sliceStart; + LOG(debug) << "Found " << nDigits << " digits for ROF " << iROF; + + // Process blocks of digits within the same chip (marked by chipID) + int iDigit = 0; + while (iDigit < nDigits) { + const int firstDigit = iDigit; + const uint16_t chipID = digits[mSortIdx[iDigit]].getChipIndex(); + + // Define the span of digits featuring the same chipID + while (iDigit < nDigits && digits[mSortIdx[iDigit]].getChipIndex() == chipID) { + ++iDigit; } - const int chipN = sliceStart - chipFirst; + const int nDigitsThisChip = iDigit - firstDigit; - LOG(debug) << "Processing chip " << chipID << " with " << chipN << " digits, next chip start from index " << sliceStart; - mThread->processChip(digits, chipFirst, chipN, &clusters, &patterns, digitLabels, clusterLabels); + LOG(debug) << "Processing chip " << chipID << " with " << nDigitsThisChip << " digits, next digit starts from index " << iDigit; + mThread->processChip(digits, firstDigit, nDigitsThisChip, &clusters, &patterns, digitLabels, clusterLabels); } - LOG(debug) << "Finished processing digit ROF " << iROF << ", produced " << (clusters.size() - outFirst) << " clusters"; - clusterROFs.emplace_back(inROF.getBCData(), inROF.getROFrame(), - outFirst, static_cast(clusters.size()) - outFirst); + LOG(debug) << "Finished processing digit ROF " << iROF << ", produced " << (clusters.size() - nStoredCls) << " clusters"; + clusterROFs.emplace_back(digitsThisROF.getBCData(), digitsThisROF.getROFrame(), + nStoredCls, static_cast(clusters.size()) - nStoredCls); } LOG(info) << "Finished processing all digit ROFs, total clusters produced: " << clusters.size(); @@ -95,112 +100,277 @@ void Clusterer::process(gsl::span digits, clusterMC2ROFs->emplace_back(in.eventRecordID, in.rofRecordID, in.minROF, in.maxROF); } } + + LOG(info) << "Writing cluster topology map to file TF3ClusterTopologies.root"; + mThread->writeTopologiesToFile("TF3ClusterTopologies.root"); } //__________________________________________________ void Clusterer::ClustererThread::processChip(gsl::span digits, - int chipFirst, int chipN, + int firstDigitIdx, int nDigits, std::vector* clustersOut, std::vector* patternsOut, const ConstDigitTruth* labelsDigPtr, ClusterTruth* labelsClusPtr) { - // chipFirst and chipN are relative to mSortIdx (i.e. mSortIdx[chipFirst..chipFirst+chipN-1] - // are the global digit indices for this chip, already sorted by col then row). + // firstDigitIdx and nDigits are relative to mSortIdx (i.e. mSortIdx[firstDigitIdx..firstDigitIdx+nDigits-1] + // are the global digit indices for this chip, already sorted by time, col then row). // We use parent->mSortIdx to resolve the global index of each pixel. - const auto& sortIdx = parent->mSortIdx; + const auto& sortIdx = mParent->mSortIdx; + LOG(info) << ""; + LOG(info) << "----------------- NEW CHIP -----------------"; - // TRK has per-ROF readout, so multiple hits belonging to the same chip, i.e. chipN > 1, - // are handled with a preclusterer. TF3 still does not have per-ROF readout, so we - // use finishChipSingleHitFast on all hits for now. - for (auto i = 0; i < chipN; ++i) { - finishChipSingleHitFast(digits, sortIdx[chipFirst + i], labelsDigPtr, labelsClusPtr); + if (nDigits == 1) { + LOG(info) << "[Clusterer] Processing single hit chip"; + findClustersSingleHit(digits, sortIdx[firstDigitIdx], labelsDigPtr, labelsClusPtr); + } else { + LOG(info) << "[Clusterer] Processing multi-hit chip with " << nDigits << " hits"; + std::vector digitIdxs(nDigits); + std::iota(digitIdxs.begin(), digitIdxs.end(), firstDigitIdx); + findClustersMultipleHits(digits, gsl::span(digitIdxs), labelsDigPtr, labelsClusPtr); } - // // TRK logic for per-ROF readout, not used for TF3 yet. - // if (chipN == 1) { - // LOG(debug) << "Processing single hit chip"; - // finishChipSingleHitFast(digits, sortIdx[chipFirst], labelsDigPtr, labelsClusPtr); - // } else { - // LOG(debug) << "Processing multi-hit chip with " << chipN << " hits"; - // // Call to initChip() - // // Call to updateChip() - // // Call to finishChip() - // // Code for preclusters needed - // } - // Flush per-thread output into the caller's containers - if (!clusters.empty()) { - clustersOut->insert(clustersOut->end(), clusters.begin(), clusters.end()); - clusters.clear(); + if (!mClusters.empty()) { + clustersOut->insert(clustersOut->end(), mClusters.begin(), mClusters.end()); + mClusters.clear(); } - if (!patterns.empty()) { - patternsOut->insert(patternsOut->end(), patterns.begin(), patterns.end()); - patterns.clear(); + if (!mPatterns.empty()) { + patternsOut->insert(patternsOut->end(), mPatterns.begin(), mPatterns.end()); + mPatterns.clear(); } - if (labelsClusPtr && labels.getNElements()) { - labelsClusPtr->mergeAtBack(labels); - labels.clear(); + if (labelsClusPtr && mLabels.getNElements()) { + labelsClusPtr->mergeAtBack(mLabels); + mLabels.clear(); } } //__________________________________________________ -void Clusterer::ClustererThread::finishChipSingleHitFast(gsl::span digits, - uint32_t digitIdx, - const ConstDigitTruth* labelsDigPtr, - ClusterTruth* labelsClusPtr) +void Clusterer::ClustererThread::findClustersSingleHit(gsl::span digits, + uint32_t digitIdx, + const ConstDigitTruth* labelsDigPtr, + ClusterTruth* labelsClusPtr) { const auto& digit = digits[digitIdx]; const uint16_t chipID = digit.getChipIndex(); const uint16_t row = digit.getRow(); const uint16_t col = digit.getColumn(); - const double time = digit.getTime(); + const time_t time = digit.getTime(); if (labelsClusPtr) { - int nlab = 0; - fetchMCLabels(digitIdx, labelsDigPtr, nlab); - const auto cnt = static_cast(clusters.size()); - for (int i = nlab; i--;) { - labels.addElement(cnt, labelsBuff[i]); + int nMcLabels = 0; + fetchMCLabels(digitIdx, labelsDigPtr, nMcLabels); + const auto nStoredCls = static_cast(mClusters.size()); + for (int i = nMcLabels; i--;) { + mLabels.addElement(nStoredCls, mLabelsBuff[i]); } } - // 1×1 pattern: rowSpan=1, colSpan=1, one byte = 0x80 - patterns.emplace_back(1); - patterns.emplace_back(1); - patterns.emplace_back(0x80); - - Cluster cluster; - cluster.chipID = chipID; - cluster.row = row; - cluster.col = col; - cluster.size = 1; - cluster.time = time; - clusters.emplace_back(cluster); + const uint16_t minRow = row; + const uint16_t minCol = col; + uint8_t rowSpan{1}, colSpan{1}, clsTopology{0}; + constexpr uint16_t firedDigitsMask = (1U << 0); // 0x0001 (1) + mClsTopoClassifier.getTopology(firedDigitsMask, minRow, rowSpan, minCol, colSpan, clsTopology); + // Bit 0 corresponds to (rowOffset=0, colOffset=0) in row-major order + Cluster cluster(row, col, rowSpan, colSpan, firedDigitsMask, clsTopology, chipID, time); + + LOG(info) << "Pushing back cluster with row: " << row << ", col: " << col << ", rowSpan: " << rowSpan + << ", colSpan: " << colSpan << ", pattern: " << firedDigitsMask + << ", topology: " << clsTopology << ", chipID: " << chipID + << ", time: " << time; + + mClusters.emplace_back(cluster); +} + +//__________________________________________________ +void Clusterer::ClustererThread::findClustersMultipleHits(gsl::span digits, + gsl::span digitIdxs, + const ConstDigitTruth* labelsDigPtr, + ClusterTruth* labelsClusPtr) +{ + + // Constraints on time resolution + const auto& digitizerParams = o2::iotof::DPLDigitizerParam::Instance(); + float timeResolution = digitizerParams.timeResolution; // in ns + const auto& clustererParams = o2::iotof::ClustererParam::Instance(); + int maxTimeDiffNSigma = clustererParams.maxTimeDiffNSigma; // in nsigma + int maxFiredDigitsForCls = clustererParams.maxFiredDigitsForCls; // max fired digits in a cluster + + // Digits are ordered by (chipID, row, col, time) within the same chip, + // so we can group them into preclusters based on adjacency in row and column. + std::vector> preclusters; + int chipID = digits[digitIdxs[0]].getChipIndex(); + for (const auto& idx : digitIdxs) { + const auto& digit = digits[idx]; + const uint16_t row = digit.getRow(); + const uint16_t col = digit.getColumn(); + + bool addedToPrecluster = false; + for (auto& precluster : preclusters) { + const auto& lastDigitIdx = precluster.back(); + const auto& lastDigit = digits[lastDigitIdx]; + if (std::abs(static_cast(lastDigit.getRow()) - static_cast(row)) <= 1 && + std::abs(static_cast(lastDigit.getColumn()) - static_cast(col)) <= 1 && + std::abs(lastDigit.getTime() - digit.getTime()) <= maxTimeDiffNSigma*timeResolution) { + precluster.push_back(idx); + addedToPrecluster = true; + break; + } + } + if (!addedToPrecluster) { + preclusters.emplace_back(std::vector{idx}); + } + } + + // Debug preclusters + LOG(info) << "[Clusterer] Found " << preclusters.size() << " preclusters in chip " << chipID; + for (size_t i = 0; i < preclusters.size(); ++i) { + LOG(info) << "Precluster " << i << " has " << preclusters[i].size() << " digits"; + } + LOG(info) << ""; + + for (const auto& precluster : preclusters) { + LOG(info) << "[Clusterer] Processing precluster with " << precluster.size() << " digits"; + + const auto nStoredCls = static_cast(mClusters.size()); + + // Single-digit cluster in chip with multiple fired digits + if (precluster.size() == 1) { + LOG(info) << "[Clusterer] Processing single-digit precluster in multi-hit chip"; + const auto& digit = digits[precluster[0]]; + const uint16_t chipID = digit.getChipIndex(); + const uint16_t row = digit.getRow(); + const uint16_t col = digit.getColumn(); + const time_t time = digit.getTime(); + + if (labelsClusPtr) { + int nMcLabels = 0; + fetchMCLabels(precluster[0], labelsDigPtr, nMcLabels); + for (int i = nMcLabels; i--;) { + mLabels.addElement(nStoredCls, mLabelsBuff[i]); + } + } + + const uint16_t minRow = row; + const uint16_t minCol = col; + uint8_t rowSpan{1}, colSpan{1}, clsTopology{0}; + // Bit 0 corresponds to (rowOffset=0, colOffset=0) in row-major order + constexpr uint16_t firedDigitsMask = (1U << 0); // 0x0001 (1) + mClsTopoClassifier.getTopology(firedDigitsMask, minRow, rowSpan, minCol, colSpan, clsTopology); + // Bit 0 corresponds to (rowOffset=0, colOffset=0) in row-major order + Cluster cluster(row, col, rowSpan, colSpan, firedDigitsMask, clsTopology, chipID, time); + + LOG(info) << "Pushing back cluster with row: " << row << ", col: " << col << ", rowSpan: " << rowSpan + << ", colSpan: " << colSpan << ", pattern: " << firedDigitsMask + << ", topology: " << clsTopology << ", chipID: " << chipID + << ", time: " << time; + + mClusters.emplace_back(cluster); + } else { + LOG(info) << "[Clusterer] Processing multi-digit precluster with " << precluster.size() << " digits"; + // Retrieve min row, min col of the precluster + uint16_t minRow = std::numeric_limits::max(); + uint16_t maxRow = std::numeric_limits::min(); + uint16_t minCol = std::numeric_limits::max(); + uint16_t maxCol = std::numeric_limits::min(); + + int nMcLabels = 0; + + // Compute average time for digits in the precluster + time_t clsTime = 0.0; + for (const auto& idx : precluster) { + const auto& digit = digits[idx]; + minRow = std::min(minRow, digit.getRow()); + minCol = std::min(minCol, digit.getColumn()); + maxRow = std::max(maxRow, digit.getRow()); + maxCol = std::max(maxCol, digit.getColumn()); + clsTime += digit.getTime(); + fetchMCLabels(idx, labelsDigPtr, nMcLabels); + } + clsTime /= precluster.size(); + const uint8_t rowSpan = maxRow - minRow + 1; + const uint8_t colSpan = maxCol - minCol + 1; + + // Fired digits bitmask packed into a single 16-bit pattern variable + uint16_t firedDigitsMask = 0; + + if (rowSpan * colSpan > maxFiredDigitsForCls) { + LOG(warn) << "Adding huge precluster with rowSpan=" << rowSpan << ", colSpan=" << colSpan; + // Overflow precluster: pass InvalidPatternID (or 0) and kHuge topology flag + Cluster cluster(minRow, minCol, rowSpan, colSpan, Cluster::InvalidPatternID, Topologies::kHuge, chipID, clsTime); + mClusters.emplace_back(cluster); + continue; + } + + // Fill firedDigitsMask in Row-Major order (bit 0 = (minRow, minCol)) + for (const auto& idx : precluster) { + const auto& digit = digits[idx]; + const uint16_t rowOffset = digit.getRow() - minRow; + const uint16_t colOffset = digit.getColumn() - minCol; + + // Single bit position calculation + const uint16_t bitIndex = rowOffset * colSpan + colOffset; + + // Set bit in LSB-to-MSB order + if (bitIndex < ClusterInfo::NBitsPattern) { + firedDigitsMask |= (1U << bitIndex); + } + } + + uint8_t clsTopology{0}; + mClsTopoClassifier.getTopology(firedDigitsMask, minRow, rowSpan, minCol, colSpan, clsTopology); + + // Construct and add cluster using scalar pattern mask + // LOG(info) << "Number of MC labels for this cluster: " << nMcLabels; + for (int i = nMcLabels; i--;) { + // LOG(info) << "[Clusterer::findClustersMultipleHits] Adding MC label " << mLabelsBuff[i] << " to cluster at index " << nStoredCls; + mLabels.addElement(nStoredCls, mLabelsBuff[i]); + } + Cluster cluster(minRow, minCol, rowSpan, colSpan, firedDigitsMask, clsTopology, chipID, clsTime); + LOG(info) << "Pushing back cluster with row: " << minRow << ", col: " << minCol << ", rowSpan: " << rowSpan + << ", colSpan: " << colSpan << ", pattern: " << firedDigitsMask + << ", topology: " << Topologies::kSingleDigit << ", chipID: " << chipID + << ", time: " << clsTime; + mClusters.emplace_back(cluster); + } + } } //__________________________________________________ void Clusterer::ClustererThread::fetchMCLabels(uint32_t digID, const ConstDigitTruth* labelsDig, int& nfilled) { + // LOG(info) << "[Clusterer::ClustererThread::fetchMCLabels] Fetching MC labels for digit ID: " << digID; if (nfilled >= MaxLabels) { + // LOG(info) << "[Clusterer::ClustererThread::fetchMCLabels] Maximum number of labels (" << MaxLabels << ") already filled, skipping further labels."; return; } if (!labelsDig || digID >= labelsDig->getIndexedSize()) { + // LOG(info) << "[Clusterer::ClustererThread::fetchMCLabels] No labels found for digit ID: " << digID; return; } const auto& lbls = labelsDig->getLabels(digID); + // LOG(info) << "[Clusterer::ClustererThread::fetchMCLabels] Digit ID: " << digID << " has " << lbls.size() << " labels"; for (int i = lbls.size(); i--;) { int ic = nfilled; for (; ic--;) { - if (labelsBuff[ic] == lbls[i]) { + if (mLabelsBuff[ic] == lbls[i]) { + // LOG(info) << "[Clusterer::ClustererThread::fetchMCLabels] Label " << lbls[i] << " already present in buffer, skipping."; return; // already present } } - labelsBuff[nfilled++] = lbls[i]; + mLabelsBuff[nfilled++] = lbls[i]; if (nfilled >= MaxLabels) { + // LOG(info) << "[Clusterer::ClustererThread::fetchMCLabels] Reached maximum number of labels (" << MaxLabels << "), stopping further label fetching."; break; } } } +//__________________________________________________ +void Clusterer::ClustererThread::writeTopologiesToFile(const char* filename) +{ + mClsTopoClassifier.saveCacheToFile("TF3ClusterTopologies.root"); +} + + } // namespace o2::iotof diff --git a/Detectors/Upgrades/ALICE3/IOTOF/reconstruction/src/ClustererParam.cxx b/Detectors/Upgrades/ALICE3/IOTOF/reconstruction/src/ClustererParam.cxx new file mode 100644 index 0000000000000..88195400528ac --- /dev/null +++ b/Detectors/Upgrades/ALICE3/IOTOF/reconstruction/src/ClustererParam.cxx @@ -0,0 +1,24 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +#include "IOTOFReconstruction/ClustererParam.h" + +O2ParamImpl(o2::iotof::ClustererParam); + +namespace o2 +{ +namespace iotof +{ +// this makes sure that the constructor of the parameters is statically +// called so that these params are part of the parameter database +static auto& sClustererParamIOTOF = o2::iotof::ClustererParam::Instance(); +} // namespace iotof +} // namespace o2 diff --git a/Detectors/Upgrades/ALICE3/IOTOF/reconstruction/src/IOTOFReconstructionLinkDef.h b/Detectors/Upgrades/ALICE3/IOTOF/reconstruction/src/IOTOFReconstructionLinkDef.h new file mode 100644 index 0000000000000..46b6d93506d59 --- /dev/null +++ b/Detectors/Upgrades/ALICE3/IOTOF/reconstruction/src/IOTOFReconstructionLinkDef.h @@ -0,0 +1,27 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +#ifdef __CLING__ + +#pragma link off all globals; +#pragma link off all classes; +#pragma link off all functions; + +#pragma link C++ class o2::iotof::Clusterer + ; + +#pragma link C++ class o2::iotof::ClustererParam + ; + +#pragma link C++ class o2::iotof::TopologyClassifier + ; + +#pragma link C++ class o2::iotof::TopologyInfo+; +#pragma link C++ class std::unordered_map+; + +#endif diff --git a/Detectors/Upgrades/ALICE3/IOTOF/reconstruction/src/TopologyClassifier.cxx b/Detectors/Upgrades/ALICE3/IOTOF/reconstruction/src/TopologyClassifier.cxx new file mode 100644 index 0000000000000..f20287501ea3a --- /dev/null +++ b/Detectors/Upgrades/ALICE3/IOTOF/reconstruction/src/TopologyClassifier.cxx @@ -0,0 +1,291 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +/// \file TopologyClassifier.cxx +/// \brief Implementation of the TopologyClassifier class. + +#include "IOTOFReconstruction/TopologyClassifier.h" +#include "DataFormatsIOTOF/Cluster.h" + +// Include for bitset +#include + +ClassImp(o2::iotof::TopologyClassifier); + +using std::array; + +namespace o2 +{ +namespace iotof +{ + +void TopologyClassifier::getTopology(uint16_t bitmask, uint16_t minRow, uint8_t spanRow, uint16_t minCol, uint8_t spanCol, uint8_t& topology) +{ + + // 1. Guard against spans exceeding 8-bit representation for + // row, col span and 16-bit bitmasks + if (spanRow > MaxRowSpan || spanCol > MaxColSpan || bitmask > MaxBitmask) { + topology = Topologies::kHuge; + return; + } + + const uint32_t clsTopoKey = packKey(spanRow, spanCol, bitmask); + // Print the 16 bits of the bitmask for debugging + LOG(info) << "[TopologyClassifier::getTopology] Bitmask: " << std::bitset<16>(bitmask) << ", minRow: " << static_cast(minRow) << ", spanRow: " << static_cast(spanRow) + << ", minCol: " << static_cast(minCol) << ", spanCol: " << static_cast(spanCol); + LOG(info) << "[TopologyClassifier::getTopology] Packed key: " << clsTopoKey; + + // Check if the topology is already cached + auto it = mTopologyCache.find(clsTopoKey); + if (it != mTopologyCache.end()) { + topology = it->second.mTopology; + it->second.mFrequency++; + LOG(info) << "[TopologyClassifier::getTopology] Found cached topology: " << static_cast(topology); + return; + } + + // Classify the new topology and cache the result + accountTopology(bitmask, minRow, spanRow, minCol, spanCol, topology); +} + + +TopologyInfo TopologyClassifier::getTopologyFeatures(uint32_t key) +{ + auto it = mTopologyCache.find(key); + if (it != mTopologyCache.end()) { + return it->second; + } else { + LOG(info) << "[TopologyClassifier::getTopologyFeatures] No cached features found for key: " << key; + return TopologyInfo(); // Return default-constructed TopologyInfo if not found + } +} + +void TopologyClassifier::accountTopology(uint16_t bitmask, uint16_t minRow, uint8_t spanRow, uint16_t minCol, uint8_t spanCol, uint8_t& topology) +{ + LOG(info) << "[TopologyClassifier::accountTopology] Classifying topology for bitmask: " << std::bitset<16>(bitmask) << ", minRow: " << static_cast(minRow) << ", spanRow: " << static_cast(spanRow) + << ", minCol: " << static_cast(minCol) << ", spanCol: " << static_cast(spanCol); + + // New cluster topology features + TopologyInfo newTopo; + newTopo.mFrequency = 1; + newTopo.mPattern = bitmask; + newTopo.mSizeX = spanRow; + newTopo.mSizeZ = spanCol; + float xCOG{0.f}, zCOG{0.f}, mXMean{0.f}, mZMean{0.f}, mXSigma2{0.f}, mZSigma2{0.f}; + computeCOG(bitmask, minRow, spanRow, minCol, spanCol, newTopo); + + const int maxRow = minRow + spanRow - 1; + const int maxCol = minCol + spanCol - 1; + + const auto hasDigit = [bitmask, minRow, minCol, spanCol](int row, int col) -> bool { + const int bitIndex = (row - minRow) * spanCol + (col - minCol); + return (bitmask & (1U << bitIndex)) != 0; + }; + + // Basic shapes + if (spanRow == 1 && spanCol == 1) { + newTopo.mTopology = Topologies::kSingleDigit; + mTopologyCache[packKey(spanRow, spanCol, bitmask)] = newTopo; + return; + } + if (spanCol == 1) { + newTopo.mTopology = Topologies::kLineOnRow; + mTopologyCache[packKey(spanRow, spanCol, bitmask)] = newTopo; + return; + } + if (spanRow == 1) { + newTopo.mTopology = Topologies::kLineOnCol; + mTopologyCache[packKey(spanRow, spanCol, bitmask)] = newTopo; + return; + } + + // Calculate total active digits in the cluster mask + int firedDigits = 0; + for (int r = minRow; r <= maxRow; ++r) { + for (int c = minCol; c <= maxCol; ++c) { + if (hasDigit(r, c)) firedDigits++; + } + } + + // Square and rectangles: all pixels fired + if (firedDigits == spanRow * spanCol && spanRow == spanCol) { + newTopo.mTopology = Topologies::kSquare; + mTopologyCache[packKey(spanRow, spanCol, bitmask)] = newTopo; + return; + } + if (firedDigits == spanRow * spanCol && spanRow != spanCol) { + newTopo.mTopology = Topologies::kRectangle; + mTopologyCache[packKey(spanRow, spanCol, bitmask)] = newTopo; + return; + } + + // Corner occupancy + const bool hasBottomLeft = hasDigit(minRow, minCol); + const bool hasBottomRight = hasDigit(minRow, maxCol); + const bool hasTopLeft = hasDigit(maxRow, minCol); + const bool hasTopRight = hasDigit(maxRow, maxCol); + + // Diagonal and triangles + if (spanRow == spanCol) { + + // Triangles + const int nCorners = hasTopLeft + hasTopRight + hasBottomLeft + hasBottomRight; + if (nCorners == 3) { + const int missing = !hasTopLeft ? 0 : !hasTopRight ? 1 : !hasBottomLeft ? 2 : 3; + + switch (missing) { + case 0: newTopo.mTopology = Topologies::kLowerTriangleLeft; break; + case 1: newTopo.mTopology = Topologies::kLowerTriangleRight; break; + case 2: newTopo.mTopology = Topologies::kUpperTriangleLeft; break; + case 3: newTopo.mTopology = Topologies::kUpperTriangleRight; break; + } + mTopologyCache[packKey(spanRow, spanCol, bitmask)] = newTopo; + return; + } + + if ((firedDigits == spanRow && hasTopLeft && hasBottomRight && !hasTopRight && !hasBottomLeft) || + (firedDigits == spanRow && hasTopRight && hasBottomLeft && !hasTopLeft && !hasBottomRight)) { + newTopo.mTopology = Topologies::kDiagonal; + mTopologyCache[packKey(spanRow, spanCol, bitmask)] = newTopo; + return; + } + } + + // Snake: 3 x 2 + if (spanRow == 3 && spanCol == 2) { + const bool hasMiddleMin = hasDigit(minRow, minCol + 1); + const bool hasMiddleMax = hasDigit(minRow, maxCol + 1); + + if (hasMiddleMin && hasMiddleMax) { + if (!hasTopLeft && !hasBottomRight && hasTopRight && hasBottomLeft) { + newTopo.mTopology = Topologies::kSnake; + mTopologyCache[packKey(spanRow, spanCol, bitmask)] = newTopo; + return; + } + + if (hasTopLeft && hasBottomRight && !hasTopRight && !hasBottomLeft) { + newTopo.mTopology = Topologies::kSnakeRefl; + mTopologyCache[packKey(spanRow, spanCol, bitmask)] = newTopo; + return; + } + } + } + + // Snake rotated by 90 degrees: 2 x 3 + if (spanRow == 2 && spanCol == 3) { + const bool hasMiddleLeft = hasDigit(minRow + 1, minCol); + const bool hasMiddleRight = hasDigit(maxRow + 1, minCol); + + if (hasMiddleLeft && hasMiddleRight) { + if (!hasTopLeft && !hasBottomRight && hasTopRight && hasBottomLeft) { + newTopo.mTopology = Topologies::kSnakeRot90; + mTopologyCache[packKey(spanRow, spanCol, bitmask)] = newTopo; + return; + } + + if (hasTopLeft && hasBottomRight && !hasTopRight && !hasBottomLeft) { + newTopo.mTopology = Topologies::kSnakeRot90Refl; + mTopologyCache[packKey(spanRow, spanCol, bitmask)] = newTopo; + return; + } + } + } + + if (newTopo.mTopology == Topologies::kNTopologies) { + newTopo.mTopology = Topologies::kOther; + mTopologyCache[packKey(spanRow, spanCol, bitmask)] = newTopo; + return; + } +} + + +void TopologyClassifier::computeCOG(uint16_t bitmask, uint16_t minRow, uint8_t spanRow, uint16_t minCol, uint8_t spanCol, TopologyInfo& topoInfo) +{ + LOG(info) << "\n\nComputing COG"; + int xOffsetCOG = 0; + int zOffsetCOG = 0; + int firedPixels = 0; + + // Ensure nBits does not exceed the bitmask capacity (16 bits) + const int nBits = std::min(static_cast(spanRow * spanCol), 16); + + for (int iBit = 0; iBit < nBits; ++iBit) { + // Check if the pixel bit is set + if (bitmask & (1U << iBit)) { + int iRow = iBit / spanCol; + int iCol = iBit % spanCol; + + xOffsetCOG += minRow + iRow; + zOffsetCOG += minCol + iCol; + LOG(info) << "Fired pixel at (row, col): (" << (minRow + iRow) << ", " << (minCol + iCol) << ")"; + LOG(info) << "Current offsets: xOffsetCOG = " << xOffsetCOG << ", zOffsetCOG = " << zOffsetCOG; + ++firedPixels; + } + } + + topoInfo.mOffsetXToCOG = static_cast((static_cast(xOffsetCOG) / firedPixels) - static_cast(minRow)); + topoInfo.mOffsetZToCOG = static_cast((static_cast(zOffsetCOG) / firedPixels) - static_cast(minCol)); + LOG(info) << "Computed COG offsets: (" << topoInfo.mOffsetXToCOG << ", " << topoInfo.mOffsetZToCOG << ")"; + topoInfo.mNPixels = firedPixels; + + LOG(info) << "COG: (" << topoInfo.mOffsetXToCOG << ", " << topoInfo.mOffsetZToCOG << "), Fired Pixels: " << firedPixels; + + // TO BE IMPLEMENTED + topoInfo.mXMean = 0.f; + topoInfo.mZMean = 0.f; + topoInfo.mXSigma2 = 0.f; + topoInfo.mZSigma2 = 0.f; + + // const auto& chipSpecs = ChipSpecificsParam::Instance(); + // if (useDf) { + // topoInfo.mXmean = dX; + // topoInfo.mZmean = dZ; + // } else { // assign expected sigmas from the pixel X, Z sizes + // topoInfo.mXsigma2 = chipSpecs.PitchRow * chipSpecs.PitchRow / 12. / std::min(10, topoInfo.mSizeX); + // topoInfo.mZsigma2 = chipSpecs.PitchCol * chipSpecs.PitchCol / 12. / std::min(10, topoInfo.mSizeZ); + // } + +} + + +void TopologyClassifier::saveCacheToFile(const char* filename) { + TFile file(filename, "RECREATE"); + // Write directly using TObject::Write syntax with explicit class name handling + file.WriteObject(&mTopologyCache, "TF3ClusterTopologies"); + file.Close(); +} + + +void TopologyClassifier::print() { + LOG(info) << "Topology Cache Contents:"; + for (const auto& entry : mTopologyCache) { + const uint32_t key = entry.first; + const TopologyInfo& topoInfo = entry.second; + + uint8_t spanRow = (key >> 24) & 0xFF; + uint8_t spanCol = (key >> 16) & 0xFF; + uint16_t bitmask = key & 0xFFFF; + + LOG(info) << "Key: " << key + << ", SpanRow: " << static_cast(spanRow) + << ", SpanCol: " << static_cast(spanCol) + << ", Bitmask: " << std::bitset<16>(bitmask) + << ", Topology: " << static_cast(topoInfo.mTopology) + << ", COGx: " << topoInfo.mOffsetXToCOG + << ", COGz: " << topoInfo.mOffsetZToCOG + << ", NPixels: " << topoInfo.mNPixels + << ", Frequency: " << topoInfo.mFrequency; + } +} + + +} // namespace o2::iotof +} // namespace o2 diff --git a/Detectors/Upgrades/ALICE3/IOTOF/simulation/CMakeLists.txt b/Detectors/Upgrades/ALICE3/IOTOF/simulation/CMakeLists.txt index 3fbb27959a2a8..edf92ea533625 100644 --- a/Detectors/Upgrades/ALICE3/IOTOF/simulation/CMakeLists.txt +++ b/Detectors/Upgrades/ALICE3/IOTOF/simulation/CMakeLists.txt @@ -16,7 +16,6 @@ o2_add_library(IOTOFSimulation src/Digitizer.cxx src/DPLDigitizerParam.cxx #src/IOTOFServices.cxx - src/Segmentation.cxx PUBLIC_LINK_LIBRARIES O2::IOTOFBase O2::DataFormatsIOTOF O2::ITSMFTSimulation) @@ -28,4 +27,4 @@ o2_target_root_dictionary(IOTOFSimulation include/IOTOFSimulation/Digitizer.h include/IOTOFSimulation/DPLDigitizerParam.h #include/IOTOFSimulation/IOTOFServices.h - include/IOTOFSimulation/Segmentation.h) + ) diff --git a/Detectors/Upgrades/ALICE3/IOTOF/simulation/include/IOTOFSimulation/Digitizer.h b/Detectors/Upgrades/ALICE3/IOTOF/simulation/include/IOTOFSimulation/Digitizer.h index ae04346ea5de1..d5ede1547e0ed 100644 --- a/Detectors/Upgrades/ALICE3/IOTOF/simulation/include/IOTOFSimulation/Digitizer.h +++ b/Detectors/Upgrades/ALICE3/IOTOF/simulation/include/IOTOFSimulation/Digitizer.h @@ -34,7 +34,7 @@ #include "SimulationDataFormat/MCCompLabel.h" #include "SimulationDataFormat/MCTruthContainer.h" #include "IOTOFBase/GeometryTGeo.h" -#include "IOTOFSimulation/Segmentation.h" +#include "IOTOFBase/Segmentation.h" namespace o2::iotof { diff --git a/Detectors/Upgrades/ALICE3/IOTOF/simulation/src/IOTOFSimulationLinkDef.h b/Detectors/Upgrades/ALICE3/IOTOF/simulation/src/IOTOFSimulationLinkDef.h index a3cadccfc6d5a..651174de8db5c 100644 --- a/Detectors/Upgrades/ALICE3/IOTOF/simulation/src/IOTOFSimulationLinkDef.h +++ b/Detectors/Upgrades/ALICE3/IOTOF/simulation/src/IOTOFSimulationLinkDef.h @@ -23,7 +23,6 @@ #pragma link C++ class o2::base::DetImpl < o2::iotof::Detector> + ; #pragma link C++ class o2::iotof::Digitizer + ; -#pragma link C++ class o2::iotof::Segmentation + ; #pragma link C++ class o2::iotof::DPLDigitizerParam + ; #pragma link C++ class o2::conf::ConfigurableParamHelper < o2::iotof::DPLDigitizerParam> + ; diff --git a/Detectors/Upgrades/ALICE3/IOTOF/workflow/src/ClusterWriterSpec.cxx b/Detectors/Upgrades/ALICE3/IOTOF/workflow/src/ClusterWriterSpec.cxx index 4d63190be5d4c..8344ba70c0ac2 100644 --- a/Detectors/Upgrades/ALICE3/IOTOF/workflow/src/ClusterWriterSpec.cxx +++ b/Detectors/Upgrades/ALICE3/IOTOF/workflow/src/ClusterWriterSpec.cxx @@ -54,8 +54,8 @@ DataProcessorSpec getClusterWriterSpec(bool mctruth, bool dec, o2::header::DataO return MakeRootTreeWriterSpec((detStr + "ClusterWriter" + (dec ? "_dec" : "")).c_str(), (detStrL + "clusters.root").c_str(), MakeRootTreeWriterSpec::TreeAttributes{.name = "o2sim", .title = "Tree with TF3 clusters"}, - BranchDefinition{InputSpec{"tf3_compclus", detOrig, "COMPCLUSTERS", 0}, - (detStr + "ClusterComp").c_str(), + BranchDefinition{InputSpec{"tf3_clus", detOrig, "CLUSTERS", 0}, + (detStr + "Cluster").c_str(), logger}, BranchDefinition{InputSpec{"tf3_patterns", detOrig, "PATTERNS", 0}, (detStr + "ClusterPatt").c_str()}, diff --git a/Detectors/Upgrades/ALICE3/IOTOF/workflow/src/ClustererSpec.cxx b/Detectors/Upgrades/ALICE3/IOTOF/workflow/src/ClustererSpec.cxx index 79d823914727a..2b60219cff684 100644 --- a/Detectors/Upgrades/ALICE3/IOTOF/workflow/src/ClustererSpec.cxx +++ b/Detectors/Upgrades/ALICE3/IOTOF/workflow/src/ClustererSpec.cxx @@ -68,7 +68,7 @@ void ClustererDPL::run(o2::framework::ProcessingContext& pc) clusterLabels.get()); LOG(info) << "Clusterization produced " << clusters.size() << " clusters for layer " << iLayer; const auto subspec = static_cast(iLayer); - pc.outputs().snapshot(o2::framework::Output{"TF3", "COMPCLUSTERS", subspec}, clusters); + pc.outputs().snapshot(o2::framework::Output{"TF3", "CLUSTERS", subspec}, clusters); pc.outputs().snapshot(o2::framework::Output{"TF3", "PATTERNS", subspec}, patterns); pc.outputs().snapshot(o2::framework::Output{"TF3", "CLUSTERSROF", subspec}, clusterROFs); if (mUseMC) { @@ -92,7 +92,7 @@ o2::framework::DataProcessorSpec getClustererSpec(bool useMC) } std::vector outputs; - outputs.emplace_back("TF3", "COMPCLUSTERS", iLayer, o2::framework::Lifetime::Timeframe); + outputs.emplace_back("TF3", "CLUSTERS", iLayer, o2::framework::Lifetime::Timeframe); outputs.emplace_back("TF3", "PATTERNS", iLayer, o2::framework::Lifetime::Timeframe); outputs.emplace_back("TF3", "CLUSTERSROF", iLayer, o2::framework::Lifetime::Timeframe); if (useMC) { diff --git a/Framework/Core/src/CommonServices.cxx b/Framework/Core/src/CommonServices.cxx index c36a102bde80d..2b6d6023ac7d5 100644 --- a/Framework/Core/src/CommonServices.cxx +++ b/Framework/Core/src/CommonServices.cxx @@ -143,7 +143,7 @@ o2::framework::ServiceSpec CommonServices::monitoringSpec() // covers devices that quit themselves via readyToQuit(). .stop = [](ServiceRegistryRef, void* service) { auto* monitoring = reinterpret_cast(service); - monitoring->finalizeProcessMonitoring(); }, + monitoring->enableProcessMonitoring(); }, .exit = [](ServiceRegistryRef registry, void* service) { auto* monitoring = reinterpret_cast(service); monitoring->flushBuffer();