From 3b631a5f2b94d146e24344b74a8bb5da40821876 Mon Sep 17 00:00:00 2001 From: Marcello Di Costanzo Date: Mon, 31 Aug 2026 10:37:47 +0200 Subject: [PATCH 1/3] Implement cluster finder with spatial and timing constraints --- .../IOTOF/DataFormatsIOTOF/CMakeLists.txt | 4 +- .../include/DataFormatsIOTOF/Cluster.h | 172 ++++++++- .../IOTOF/DataFormatsIOTOF/src/Cluster.cxx | 60 +++- .../src/DataFormatsIOTOFLinkDef.h | 1 + .../Upgrades/ALICE3/IOTOF/base/CMakeLists.txt | 2 + .../include/IOTOFBase}/Segmentation.h | 24 +- .../ALICE3/IOTOF/base/src/IOTOFBaseLinkDef.h | 1 + .../{simulation => base}/src/Segmentation.cxx | 4 +- .../ALICE3/IOTOF/macros/CheckClustersIOTOF.C | 2 +- .../ALICE3/IOTOF/macros/CheckDigitsIOTOF.C | 3 +- .../IOTOF/reconstruction/CMakeLists.txt | 10 + .../include/IOTOFReconstruction/Clusterer.h | 29 +- .../IOTOFReconstruction/ClustererParam.h | 49 +++ .../IOTOFReconstruction/TopologyClassifier.h | 105 ++++++ .../IOTOF/reconstruction/src/Clusterer.cxx | 334 +++++++++++++----- .../reconstruction/src/ClustererParam.cxx | 24 ++ .../src/IOTOFReconstructionLinkDef.h | 25 ++ .../reconstruction/src/TopologyClassifier.cxx | 278 +++++++++++++++ .../ALICE3/IOTOF/simulation/CMakeLists.txt | 3 +- .../include/IOTOFSimulation/Digitizer.h | 2 +- .../simulation/src/IOTOFSimulationLinkDef.h | 1 - .../IOTOF/workflow/src/ClusterWriterSpec.cxx | 4 +- .../IOTOF/workflow/src/ClustererSpec.cxx | 4 +- 23 files changed, 1005 insertions(+), 136 deletions(-) rename Detectors/Upgrades/ALICE3/IOTOF/{simulation/include/IOTOFSimulation => base/include/IOTOFBase}/Segmentation.h (94%) rename Detectors/Upgrades/ALICE3/IOTOF/{simulation => base}/src/Segmentation.cxx (93%) create mode 100644 Detectors/Upgrades/ALICE3/IOTOF/reconstruction/include/IOTOFReconstruction/ClustererParam.h create mode 100644 Detectors/Upgrades/ALICE3/IOTOF/reconstruction/include/IOTOFReconstruction/TopologyClassifier.h create mode 100644 Detectors/Upgrades/ALICE3/IOTOF/reconstruction/src/ClustererParam.cxx create mode 100644 Detectors/Upgrades/ALICE3/IOTOF/reconstruction/src/IOTOFReconstructionLinkDef.h create mode 100644 Detectors/Upgrades/ALICE3/IOTOF/reconstruction/src/TopologyClassifier.cxx 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..8eb1cf3b03944 100644 --- a/Detectors/Upgrades/ALICE3/IOTOF/DataFormatsIOTOF/include/DataFormatsIOTOF/Cluster.h +++ b/Detectors/Upgrades/ALICE3/IOTOF/DataFormatsIOTOF/include/DataFormatsIOTOF/Cluster.h @@ -1,4 +1,4 @@ -// Copyright 2019-2026 CERN and copyright holders of ALICE O2. +// 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. // @@ -9,28 +9,172 @@ // granted to it by virtue of its status as an Intergovernmental Organization // or submit itself to any jurisdiction. -#ifndef ALICEO2_DATAFORMATSIOTOF_CLUSTER_H -#define ALICEO2_DATAFORMATSIOTOF_CLUSTER_H +/// \file Cluster.h +/// \brief Definition of the IOTOF cluster +#ifndef ALICEO2_IOTOF_CLUSTER_H +#define ALICEO2_IOTOF_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_IOTOF_CLUSTER_H */ \ No newline at end of file diff --git a/Detectors/Upgrades/ALICE3/IOTOF/DataFormatsIOTOF/src/Cluster.cxx b/Detectors/Upgrades/ALICE3/IOTOF/DataFormatsIOTOF/src/Cluster.cxx index 6b5a4948900e7..6d1ff66bf9a3a 100644 --- a/Detectors/Upgrades/ALICE3/IOTOF/DataFormatsIOTOF/src/Cluster.cxx +++ b/Detectors/Upgrades/ALICE3/IOTOF/DataFormatsIOTOF/src/Cluster.cxx @@ -1,4 +1,4 @@ -// Copyright 2019-2026 CERN and copyright holders of ALICE O2. +// 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. // @@ -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"; } -} // namespace o2::iotof +//______________________________________________________________________________ +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 iotof +} // namespace o2 + +// Stream operator implementation +std::ostream& operator<<(std::ostream& stream, const o2::iotof::Cluster& cl) +{ + stream << cl.asString(); + return stream; +} \ No newline at end of file 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 94% rename from Detectors/Upgrades/ALICE3/IOTOF/simulation/include/IOTOFSimulation/Segmentation.h rename to Detectors/Upgrades/ALICE3/IOTOF/base/include/IOTOFBase/Segmentation.h index ddde28cf7dd7a..504b050486fc2 100644 --- a/Detectors/Upgrades/ALICE3/IOTOF/simulation/include/IOTOFSimulation/Segmentation.h +++ b/Detectors/Upgrades/ALICE3/IOTOF/base/include/IOTOFBase/Segmentation.h @@ -50,9 +50,9 @@ 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 /// local coordinates (cm) @@ -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/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 93% rename from Detectors/Upgrades/ALICE3/IOTOF/simulation/src/Segmentation.cxx rename to Detectors/Upgrades/ALICE3/IOTOF/base/src/Segmentation.cxx index a7ec0d708c3b8..31517139e2279 100644 --- a/Detectors/Upgrades/ALICE3/IOTOF/simulation/src/Segmentation.cxx +++ b/Detectors/Upgrades/ALICE3/IOTOF/base/src/Segmentation.cxx @@ -12,8 +12,8 @@ /// \file Segmentation.cxx /// \brief Implementation of the Segmentation class -#include "IOTOFSimulation/Segmentation.h" -#include "IOTOFBase/IOTOFBaseParam.h" +#include +#include #include namespace o2 diff --git a/Detectors/Upgrades/ALICE3/IOTOF/macros/CheckClustersIOTOF.C b/Detectors/Upgrades/ALICE3/IOTOF/macros/CheckClustersIOTOF.C index 107e5a4d02bf8..d85797a38a75c 100644 --- a/Detectors/Upgrades/ALICE3/IOTOF/macros/CheckClustersIOTOF.C +++ b/Detectors/Upgrades/ALICE3/IOTOF/macros/CheckClustersIOTOF.C @@ -20,7 +20,7 @@ #include #include -#include "IOTOFSimulation/Segmentation.h" +#include "IOTOFBase/Segmentation.h" #include "IOTOFBase/IOTOFBaseParam.h" #include "IOTOFBase/GeometryTGeo.h" #include "DataFormatsIOTOF/Digit.h" diff --git a/Detectors/Upgrades/ALICE3/IOTOF/macros/CheckDigitsIOTOF.C b/Detectors/Upgrades/ALICE3/IOTOF/macros/CheckDigitsIOTOF.C index 26ffd08697d56..af4e59de827f8 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" @@ -77,6 +77,7 @@ void addTLines(float pitch) void CheckDigitsIOTOF(std::string digifile = "tf3digits.root", std::string hitfile = "o2sim_HitsTF3.root", std::string inputGeom = "o2sim_geometry.root") { + std::cout << "\ndigifile=" << digifile << "\nhitfile=" << hitfile << "\ninputGeom=" << inputGeom << std::endl; gStyle->SetPalette(55); using namespace o2::base; 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..af007ada4c530 --- /dev/null +++ b/Detectors/Upgrades/ALICE3/IOTOF/reconstruction/include/IOTOFReconstruction/ClustererParam.h @@ -0,0 +1,49 @@ +// 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"); + + private: + static constexpr float DEFNoisePerPixel() + { + return 1e-8; // ITS/MFT values here!! + } +}; + +} // 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..d837dae3948d2 --- /dev/null +++ b/Detectors/Upgrades/ALICE3/IOTOF/reconstruction/include/IOTOFReconstruction/TopologyClassifier.h @@ -0,0 +1,105 @@ +// 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, + 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..8d6b3e3e1fa14 --- /dev/null +++ b/Detectors/Upgrades/ALICE3/IOTOF/reconstruction/src/IOTOFReconstructionLinkDef.h @@ -0,0 +1,25 @@ +// 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::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..c0cf00c8be9c4 --- /dev/null +++ b/Detectors/Upgrades/ALICE3/IOTOF/reconstruction/src/TopologyClassifier.cxx @@ -0,0 +1,278 @@ +// 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; + } + + // Corner occupancy + const bool hasTopLeft = hasDigit(minRow, minCol); + const bool hasTopRight = hasDigit(minRow, maxCol); + const bool hasBottomLeft = hasDigit(maxRow, minCol); + const bool hasBottomRight = hasDigit(maxRow, maxCol); + + // Diagonal and square + if (spanRow == spanCol) { + + if ((hasTopLeft && hasBottomRight && !hasTopRight && !hasBottomLeft) || + (!hasTopLeft && !hasBottomRight && hasTopRight && hasBottomLeft)) { + newTopo.mTopology = Topologies::kDiagonal; + mTopologyCache[packKey(spanRow, spanCol, bitmask)] = newTopo; + return; + } + + if (hasTopLeft && hasTopRight && hasBottomLeft && hasBottomRight) { + newTopo.mTopology = Topologies::kSquare; + mTopologyCache[packKey(spanRow, spanCol, bitmask)] = newTopo; + return; + } + + // Triangles (exactly one missing corner) + 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; + } + } + + // Snake: 3 x 2 + if (spanRow == 3 && spanCol == 2) { + const bool hasMiddleMin = hasDigit(minRow + 1, minCol); + const bool hasMiddleMax = hasDigit(minRow + 1, maxCol); + + if (hasMiddleMin && hasMiddleMax) { + if (hasTopLeft && hasBottomRight) { + newTopo.mTopology = Topologies::kSnake; + mTopologyCache[packKey(spanRow, spanCol, bitmask)] = newTopo; + return; + } + + if (!hasTopLeft && !hasBottomRight) { + 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, minCol + 1); + const bool hasMiddleRight = hasDigit(maxRow, minCol + 1); + + if (hasMiddleLeft && hasMiddleRight) { + if (hasTopLeft && hasBottomRight) { + newTopo.mTopology = Topologies::kSnakeRot90; + mTopologyCache[packKey(spanRow, spanCol, bitmask)] = newTopo; + return; + } + + if (!hasTopLeft && !hasBottomRight) { + 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; + } + // Insert in map +} + + +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) { From c0bab8f7a50314d53f0536a6ab116fb03ff3c626 Mon Sep 17 00:00:00 2001 From: Marcello Di Costanzo Date: Tue, 1 Sep 2026 19:33:42 +0200 Subject: [PATCH 2/3] Update QA macro --- .../include/DataFormatsIOTOF/Cluster.h | 8 +- .../IOTOF/DataFormatsIOTOF/src/Cluster.cxx | 4 +- .../ALICE3/IOTOF/base/src/Segmentation.cxx | 4 +- .../ALICE3/IOTOF/macros/CheckClustersIOTOF.C | 1196 ++++++++++++++--- .../ALICE3/IOTOF/macros/CheckDigitsIOTOF.C | 1 - 5 files changed, 1032 insertions(+), 181 deletions(-) diff --git a/Detectors/Upgrades/ALICE3/IOTOF/DataFormatsIOTOF/include/DataFormatsIOTOF/Cluster.h b/Detectors/Upgrades/ALICE3/IOTOF/DataFormatsIOTOF/include/DataFormatsIOTOF/Cluster.h index 8eb1cf3b03944..21028d21c9cde 100644 --- a/Detectors/Upgrades/ALICE3/IOTOF/DataFormatsIOTOF/include/DataFormatsIOTOF/Cluster.h +++ b/Detectors/Upgrades/ALICE3/IOTOF/DataFormatsIOTOF/include/DataFormatsIOTOF/Cluster.h @@ -1,4 +1,4 @@ -// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// Copyright 2019-2026 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. // @@ -11,8 +11,8 @@ /// \file Cluster.h /// \brief Definition of the IOTOF cluster -#ifndef ALICEO2_IOTOF_CLUSTER_H -#define ALICEO2_IOTOF_CLUSTER_H +#ifndef ALICEO2_DATAFORMATSIOTOF_CLUSTER_H +#define ALICEO2_DATAFORMATSIOTOF_CLUSTER_H #include #include @@ -177,4 +177,4 @@ class Cluster std::ostream& operator<<(std::ostream& stream, const o2::iotof::Cluster& cl); -#endif /* ALICEO2_IOTOF_CLUSTER_H */ \ No newline at end of file +#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 6d1ff66bf9a3a..22735a9225c19 100644 --- a/Detectors/Upgrades/ALICE3/IOTOF/DataFormatsIOTOF/src/Cluster.cxx +++ b/Detectors/Upgrades/ALICE3/IOTOF/DataFormatsIOTOF/src/Cluster.cxx @@ -1,4 +1,4 @@ -// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// Copyright 2019-2026 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. // @@ -70,4 +70,4 @@ std::ostream& operator<<(std::ostream& stream, const o2::iotof::Cluster& cl) { stream << cl.asString(); return stream; -} \ No newline at end of file +} diff --git a/Detectors/Upgrades/ALICE3/IOTOF/base/src/Segmentation.cxx b/Detectors/Upgrades/ALICE3/IOTOF/base/src/Segmentation.cxx index 31517139e2279..aa77bf50d069d 100644 --- a/Detectors/Upgrades/ALICE3/IOTOF/base/src/Segmentation.cxx +++ b/Detectors/Upgrades/ALICE3/IOTOF/base/src/Segmentation.cxx @@ -12,8 +12,8 @@ /// \file Segmentation.cxx /// \brief Implementation of the Segmentation class -#include -#include +#include "IOTOFBase/Segmentation.h" +#include "IOTOFBase/IOTOFBaseParam.h" #include namespace o2 diff --git a/Detectors/Upgrades/ALICE3/IOTOF/macros/CheckClustersIOTOF.C b/Detectors/Upgrades/ALICE3/IOTOF/macros/CheckClustersIOTOF.C index d85797a38a75c..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 "IOTOFBase/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; -void CheckClustersIOTOF(std::string digiFilePath = "tf3digits.root", std::string clsFilePath = "tf3clusters.root", std::string inputGeomPath = "o2sim_geometry.root") -{ - gStyle->SetPalette(55); +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}); +} + + +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); + + for (int i = 0; i < chipHitsIdxs.size(); ++i) { + const auto& hit = (*evtChipHits)[chipHitsIdxs[i].hitIdx]; - using namespace o2::base; - using namespace o2::iotof; + 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) +} - using o2::iotof::Cluster; - using o2::iotof::Digit; + +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); } - digitsLabels.emplace(label, iDigit); } - // LOOP on : ROFRecord array - for (unsigned int iROF = 0; iROF < clsRofArr.size(); ++iROF) { + // 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); + } + } + + // 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"); - const unsigned int rofIndex = clsRofArr[iROF].getFirstEntry(); - const unsigned int rofNEntries = clsRofArr[iROF].getNEntries(); + std::vector clustersProperties; + clustersProperties.reserve(clustersArray->size()); // Pre-allocate memory - // 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; + for (int iCls = 0; iCls < (int)clustersArray->size(); ++iCls) { + + 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()); + } + + 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"}; - 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); + // 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 + } + } - o2::math_utils::Point3D localDigitCoord(x, y, z); // local Digit + 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); + } - 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 + 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); - // 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; + 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); + + 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); + } + } + } - Float_t x = 0.f, y = 0.f, z = 0.f; - if (subDetID >= 0) { - segGeom->detectorToLocal(iRow, iCol, x, z, subDetID); + // 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); + } + + 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()); + } + } + } - 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()); + Print(true, "----> Looping over generated particles"); + + // 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 af4e59de827f8..581c93a236c98 100644 --- a/Detectors/Upgrades/ALICE3/IOTOF/macros/CheckDigitsIOTOF.C +++ b/Detectors/Upgrades/ALICE3/IOTOF/macros/CheckDigitsIOTOF.C @@ -77,7 +77,6 @@ void addTLines(float pitch) void CheckDigitsIOTOF(std::string digifile = "tf3digits.root", std::string hitfile = "o2sim_HitsTF3.root", std::string inputGeom = "o2sim_geometry.root") { - std::cout << "\ndigifile=" << digifile << "\nhitfile=" << hitfile << "\ninputGeom=" << inputGeom << std::endl; gStyle->SetPalette(55); using namespace o2::base; From e5e7678d4add8d4c7244eecc236e4acad0f74196 Mon Sep 17 00:00:00 2001 From: Marcello Di Costanzo Date: Thu, 3 Sep 2026 10:21:14 +0200 Subject: [PATCH 3/3] Before ITS topology implementation --- .../base/include/IOTOFBase/Segmentation.h | 2 +- .../ALICE3/IOTOF/base/src/GeometryTGeo.cxx | 2 +- .../ALICE3/IOTOF/macros/CheckClustersIOTOF.C | 1423 ++++++++++------- .../ALICE3/IOTOF/macros/CheckDigitsIOTOF.C | 158 +- .../IOTOFReconstruction/ClustererParam.h | 6 - .../IOTOFReconstruction/TopologyClassifier.h | 12 + .../IOTOF/reconstruction/src/Clusterer.cxx | 32 +- .../src/IOTOFReconstructionLinkDef.h | 2 + .../reconstruction/src/TopologyClassifier.cxx | 77 +- .../ALICE3/IOTOF/simulation/src/Digitizer.cxx | 2 +- Framework/Core/src/CommonServices.cxx | 2 +- 11 files changed, 1075 insertions(+), 643 deletions(-) diff --git a/Detectors/Upgrades/ALICE3/IOTOF/base/include/IOTOFBase/Segmentation.h b/Detectors/Upgrades/ALICE3/IOTOF/base/include/IOTOFBase/Segmentation.h index 504b050486fc2..c726998fcd4bc 100644 --- a/Detectors/Upgrades/ALICE3/IOTOF/base/include/IOTOFBase/Segmentation.h +++ b/Detectors/Upgrades/ALICE3/IOTOF/base/include/IOTOFBase/Segmentation.h @@ -54,7 +54,7 @@ class Segmentation /// same but w/o check for row/column range 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 diff --git a/Detectors/Upgrades/ALICE3/IOTOF/base/src/GeometryTGeo.cxx b/Detectors/Upgrades/ALICE3/IOTOF/base/src/GeometryTGeo.cxx index 8c8a36877eca8..e54e21e07df56 100644 --- a/Detectors/Upgrades/ALICE3/IOTOF/base/src/GeometryTGeo.cxx +++ b/Detectors/Upgrades/ALICE3/IOTOF/base/src/GeometryTGeo.cxx @@ -313,7 +313,7 @@ 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); diff --git a/Detectors/Upgrades/ALICE3/IOTOF/macros/CheckClustersIOTOF.C b/Detectors/Upgrades/ALICE3/IOTOF/macros/CheckClustersIOTOF.C index cd8d7e31ad227..5a6a4b51134b5 100644 --- a/Detectors/Upgrades/ALICE3/IOTOF/macros/CheckClustersIOTOF.C +++ b/Detectors/Upgrades/ALICE3/IOTOF/macros/CheckClustersIOTOF.C @@ -10,11 +10,12 @@ // or submit itself to any jurisdiction. /// \file CheckClustersIOTOF.C -/// \brief Simple macro to create clusters from TF3 digits +/// \brief QA macro for TF3 clusters #if !defined(__CLING__) || defined(__ROOTCLING__) #include +#include #include #include @@ -50,37 +51,6 @@ 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) { @@ -91,7 +61,6 @@ void GetHitAvgPositionGlobal(const o2::itsmft::Hit& hit, o2::math_utils::Point3D } - void GetHitAvgPositionLocal(const o2::itsmft::Hit& hit, o2::iotof::GeometryTGeo* geom, o2::math_utils::Point3D& avgPos) { const int chipID = hit.GetDetectorID(); @@ -105,22 +74,6 @@ void GetHitAvgPositionLocal(const o2::itsmft::Hit& hit, o2::iotof::GeometryTGeo* } -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; @@ -135,44 +88,41 @@ void PrintHit(bool verbose, o2::itsmft::Hit hit, o2::iotof::GeometryTGeo* iotofG } 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; + o2::math_utils::Point3D avgPos; + GetHitAvgPositionGlobal(hit, avgPos); + std::cout << "Hit: detectorID = " << hit.GetDetectorID() << ", avgPos = (" << avgPos.X() << ", " + << avgPos.Y() << ", " << avgPos.Z() << ")" << ", 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(); +// Fare residuo in-chip +void GetClusterGlobalPos(const o2::iotof::Cluster& cluster, + TopologyInfo topoInfo, + o2::math_utils::Point3D& globalPos, + o2::iotof::GeometryTGeo* iotofGeom, + o2::iotof::Segmentation* segmInfo){ - 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; + std::cout << "Computing cluster global position for cluster with bottom left corner at (row = " << cluster.getRow() << ", col = " << cluster.getCol() << "), chipID = " << cluster.getChipID() << std::endl; + float x = 0.f; + float y = 0.f; + float z = 0.f; + int rowCOG = cluster.getRow() + topoInfo.mOffsetXToCOG; + int colCOG = cluster.getCol() + topoInfo.mOffsetZToCOG; + topoInfo.print(); + std::cout << "Cluster COG at (row = " << rowCOG << ", col = " << colCOG << ")" << std::endl; + segmInfo->detectorToLocal(rowCOG, colCOG, x, z, cluster.getChipID()); + globalPos = iotofGeom->getMatrixL2G(cluster.getChipID())(o2::math_utils::Point3D{x, 0.f, z}); } void PrintCluster(bool verbose, const o2::iotof::Cluster& cluster, auto clsLabel, + TopologyInfo topoInfo, o2::iotof::GeometryTGeo* iotofGeom, o2::iotof::Segmentation* segmInfo) { if (!verbose) { @@ -182,16 +132,28 @@ void PrintCluster(bool verbose, 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; + o2::math_utils::Point3D clsPos; + GetClusterGlobalPos(cluster, topoInfo, clsPos, iotofGeom, segmInfo); - const int eventID = evtTrackLabel.getEventID(); - const int trackID = evtTrackLabel.getTrackID(); - std::cout << " Label " << iLabel << ", eventID=" << eventID << ", trackID=" << trackID << std::endl; - } + int layer{-1}, stave{-1}, subStave{-1}, module{-1}, chip{-1}; + iotofGeom->getIOTOFChipId(cluster.getChipID(), layer, stave, subStave, module, chip); + + std::cout << "Cluster: chipID=" << cluster.getChipID() << ", row=" << cluster.getRow() << ", col=" << cluster.getCol() + << ", layer=" << layer << ", stave=" << stave << ", subStave=" << subStave + << ", module=" << module << ", chip=" << chip + << ", rowSpan=" << cluster.getRowSpan() << ", colSpan=" << cluster.getColSpan() + << ", size=" << cluster.getSize() << ", labels=" << clsLabel.size() << " MCCompLabels" + << ", topology=" << cluster.getTopology() << ", time=" << cluster.getTime() + << 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 << " Contributing track to cls, label " << iLabel << ", eventID=" << eventID << ", trackID=" << trackID << std::endl; + // } } @@ -205,27 +167,29 @@ void Print(bool verbose, Args&&... args) { } -void GetClusterGlobalPos(const o2::iotof::Cluster& cluster, - TopologyInfo topoInfo, - o2::math_utils::Point3D& globalPos, - o2::iotof::GeometryTGeo* iotofGeom, - o2::iotof::Segmentation* segmInfo){ +void GetClusterLocalPos(const o2::iotof::Cluster& cluster, + TopologyInfo topoInfo, + o2::math_utils::Point3D& localPos, + o2::iotof::GeometryTGeo* iotofGeom, + o2::iotof::Segmentation* segmInfo){ + std::cout << "Computing cluster global position for cluster with bottom left corner at (row = " << cluster.getRow() << ", col = " << cluster.getCol() << "), chipID = " << cluster.getChipID() << std::endl; float x = 0.f; float y = 0.f; float z = 0.f; int rowCOG = cluster.getRow() + topoInfo.mOffsetXToCOG; int colCOG = cluster.getCol() + topoInfo.mOffsetZToCOG; + topoInfo.print(); + std::cout << "Cluster COG at (row = " << rowCOG << ", col = " << colCOG << ")" << std::endl; segmInfo->detectorToLocal(rowCOG, colCOG, x, z, cluster.getChipID()); - globalPos = iotofGeom->getMatrixL2G(cluster.getChipID())(o2::math_utils::Point3D{x, 0.f, z}); + localPos = o2::math_utils::Point3D{x, 0.f, z}; } int FindBestMatchingHit(const o2::iotof::Cluster& cluster, - TopologyInfo topoInfo, - std::vector& chipHitsIdxs, - std::vector* evtChipHits, - const std::vector* digitsArray, + TopologyInfo topoInfo, + const std::vector& chipHitsIdxs, + const std::vector* evtChipHits, o2::iotof::GeometryTGeo* iotofGeom, o2::iotof::Segmentation* segmInfo){ int bestHitIdx = -1; @@ -234,7 +198,7 @@ int FindBestMatchingHit(const o2::iotof::Cluster& cluster, GetClusterGlobalPos(cluster, topoInfo, clsPos, iotofGeom, segmInfo); for (int i = 0; i < chipHitsIdxs.size(); ++i) { - const auto& hit = (*evtChipHits)[chipHitsIdxs[i].hitIdx]; + const auto& hit = (*evtChipHits)[chipHitsIdxs[i]]; float dx = clsPos.X() - hit.GetX(); float dy = clsPos.Y() - hit.GetY(); @@ -251,18 +215,47 @@ int FindBestMatchingHit(const o2::iotof::Cluster& 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; + int nAssocPrimaries = 0; // More than one primary MC particle from the same event is associated to the cluster + bool isShared = false; // More than one primary MC particle from the same event is associated to the cluster + bool isFake = false; // More than one primary MC particle from different events is associated to the cluster + int hitIdx = -1; + Topologies topology = kOther; + uint32_t topoKey = 0; +}; + +struct DetectorData { + std::vector hitIndicesL0; + std::vector hitIndicesL1; + std::vector clsIndicesL0; + std::vector clsIndicesL1; +}; + + 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", + std::string geomCfgStr = "IOTOFBase.segmentedInnerTOF=true;IOTOFBase.segmentedOuterTOF=true;IOTOFBase.enableForwardTOF=false;IOTOFBase.enableBackwardTOF=false;", bool verbose = false) { - Print(verbose, "CheckClustersTopologiesIOTOF: kinefile = ", kinefile, ", hitfile = ", hitfile, ", digiFilePath = ", digiFilePath, ", clsFilePath = ", clsFilePath, ", inputGeomPath = ", inputGeomPath); + Print(verbose, "CheckClustersTopologiesIOTOF: kinefile = ", kinefile, ", hitfile = ", hitfile, ", clsFilePath = ", clsFilePath, ", inputGeomPath = ", inputGeomPath); gStyle->SetPalette(55); - o2::conf::ConfigurableParam::updateFromString("IOTOFBase.segmentedInnerTOF=true;IOTOFBase.segmentedOuterTOF=true;IOTOFBase.enableForwardTOF=false;IOTOFBase.enableBackwardTOF=false"); + o2::conf::ConfigurableParam::updateFromString(geomCfgStr); auto segmInfo = o2::iotof::Segmentation::Instance(); @@ -275,17 +268,57 @@ void CheckClustersIOTOF(std::string kinefile = "o2sim_Kine.root", 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; + std::cout << "Loaded " << clsTopoMapPtr->size() << " entries from " << clsFileTopoPath << std::endl; } else { - std::cerr << "Failed to load TF3ClusterTopologies from file!" << std::endl; + std::cerr << "Failed to load TF3ClusterTopologies from " << clsFileTopoPath << std::endl; } - clsTopoFile->Close(); - std::cout << std::endl; - std::cout << "Topologies summary: " << std::endl; - TopologyClassifier topoClassifier(*clsTopoMapPtr); + // Construct map directly from the vector pairs + std::unordered_map topoMap(clsTopoMapPtr->begin(), clsTopoMapPtr->end()); + std::cout << "\nTopologies summary:" << std::endl; + TopologyClassifier topoClassifier(std::move(topoMap)); topoClassifier.print(); std::cout << std::endl; + clsTopoFile->Close(); + + // Sorted topology map by spanRow, spanCol, and then by bitmask for better organization in the output file + auto topologyMap = topoClassifier.getTopologyMap(); + std::vector> sortedTopoMap(topologyMap.begin(), topologyMap.end()); + std::sort(sortedTopoMap.begin(), sortedTopoMap.end(), [](const auto& a, const auto& b) { + int topoA = a.second.mTopology; + int topoB = b.second.mTopology; + uint8_t spanRowA = (a.first >> 24) & 0xFF; + uint8_t spanColA = (a.first >> 16) & 0xFF; + uint8_t spanRowB = (b.first >> 24) & 0xFF; + uint8_t spanColB = (b.first >> 16) & 0xFF; + int nPixelsA = a.second.mNPixels; + int nPixelsB = b.second.mNPixels; + int frequencyA = a.second.mFrequency; + int frequencyB = b.second.mFrequency; + if (topoA != topoB) return topoA < topoB; + if (frequencyA != frequencyB) return frequencyA > frequencyB; + if (spanRowA != spanRowB) return spanRowA < spanRowB; + if (spanColA != spanColB) return spanColA < spanColB; + if (nPixelsA != nPixelsB) return nPixelsA < nPixelsB; + return a.first < b.first; // Finally sort by bitmask if spans are equal + }); + + // Print the sorted topology map + for (const auto& entry : sortedTopoMap) { + 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); + topoInfo.print(); + LOG(info) << ""; + } // Generated MC tracks and TrackRefs information TFile* kineFile = TFile::Open(kinefile.data()); @@ -299,26 +332,6 @@ void CheckClustersIOTOF(std::string kinefile = "o2sim_Kine.root", 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]); @@ -332,146 +345,163 @@ void CheckClustersIOTOF(std::string kinefile = "o2sim_Kine.root", 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 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); - // 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 ... "); + // Store hit, cluster and MC particles properties for all tracks in all events + std::unordered_map tracksHitCls; + + // Load hits and MC tracks, which are stored per-event + // Generated particles + TH2F* hGenEtaPt[2] = {new TH2F("hGenEtaPtPrm", "Generated primary tracks;#eta;p_{T}", 40, -2, 2, 100, 0, 10), + new TH2F("hGenEtaPtSec", "Generated secondary tracks;#eta;p_{T}", 40, -2, 2, 100, 0, 10)}; for (int iEvt = 0; iEvt < nEvts; ++iEvt) { + for (int iTrack = 0; iTrack < mcTracksPerEvent[iEvt]->size(); ++iTrack) { + const auto& mcTrack = (*mcTracksPerEvent[iEvt])[iTrack]; + // if (!mcTrack.isPrimary()) + // continue; + const int type = mcTrack.isPrimary() ? 0 : 1; + hGenEtaPt[type]->Fill(mcTrack.GetEta(), mcTrack.GetPt()); + const uint64_t trackKey = (static_cast(iEvt) << 32) | static_cast(iTrack); + tracksHitCls[trackKey] = DetectorData(); + } + } - Print(verbose, "Event ", iEvt, ": ", hitsPerEvent[iEvt]->size(), " hits"); + int nHits{0}, nHitsFromPrimaryTracks{0}, nHitsFromSecondaryTracks{0}; + TH2F* hEtaPhiHitsPrmTrkLayer0 = new TH2F("hEtaPhiHitsPrmTrkLayer0", "hEtaPhiHitsPrmTrkLayer0;#phi;#eta", 64, 0, 6.28319, 40, -2, 2); + TH2F* hEtaPhiHitsSecTrkLayer0 = new TH2F("hEtaPhiHitsSecTrkLayer0", "hEtaPhiHitsSecTrkLayer0;#phi;#eta", 64, 0, 6.28319, 40, -2, 2); + TH2F* hEtaPhiHitsPrmTrkLayer1 = new TH2F("hEtaPhiHitsPrmTrkLayer1", "hEtaPhiHitsPrmTrkLayer1;#phi;#eta", 64, 0, 6.28319, 40, -2, 2); + TH2F* hEtaPhiHitsSecTrkLayer1 = new TH2F("hEtaPhiHitsSecTrkLayer1", "hEtaPhiHitsSecTrkLayer1;#phi;#eta", 64, 0, 6.28319, 40, -2, 2); + TH2F* hEtaPtHitsPrmTrkLayer0 = new TH2F("hEtaPtHitsPrmLayer0", "Generated primary tracks;#eta;p_{T}", 40, -2, 2, 50, 0, 10); + TH2F* hEtaPtHitsSecTrkLayer0 = new TH2F("hEtaPtHitsSecLayer0", "Generated secondary tracks;#eta;p_{T}", 40, -2, 2, 50, 0, 10); + TH2F* hEtaPtHitsPrmTrkLayer1 = new TH2F("hEtaPtHitsPrmLayer1", "Generated primary tracks;#eta;p_{T}", 40, -2, 2, 50, 0, 10); + TH2F* hEtaPtHitsSecTrkLayer1 = new TH2F("hEtaPtHitsSecLayer1", "Generated secondary tracks;#eta;p_{T}", 40, -2, 2, 50, 0, 10); + for (int iEvt = 0; iEvt < nEvts; ++iEvt) { for (int iHit = 0; iHit < hitsPerEvent[iEvt]->size(); ++iHit) { - + nHits++; 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++; + const uint64_t trackKey = (static_cast(iEvt) << 32) | static_cast(trackID); + if (tracksHitCls.find(trackKey) == tracksHitCls.end()) { + continue; + } + int layer = iotofGeom->getIOTOFLayer(hit.GetDetectorID()); + if (layer == 0) { + tracksHitCls[trackKey].hitIndicesL0.push_back(iHit); + } else if (layer == 1) { + tracksHitCls[trackKey].hitIndicesL1.push_back(iHit); + } - // Fill histograms - int hitLayer = iotofGeom->getIOTOFLayer(hit.GetDetectorID()); - auto& mcTrack = (*mcTracksPerEvent[iEvt])[trackID]; + auto &mcTrack = (*mcTracksPerEvent[iEvt])[trackID]; bool isPrimary = mcTrack.isPrimary(); if (isPrimary) nHitsFromPrimaryTracks++; else nHitsFromSecondaryTracks++; + float genEta = mcTrack.GetEta(); float genPhi = mcTrack.GetPhi(); + float genPt = mcTrack.GetPt(); - 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); + int hitLayer = iotofGeom->getIOTOFLayer(hit.GetDetectorID()); + if (hitLayer == 0 && isPrimary) { + hEtaPhiHitsPrmTrkLayer0->Fill(genPhi, genEta); + hEtaPtHitsPrmTrkLayer0->Fill(genEta, genPt); + } + else if (hitLayer == 0 && !isPrimary) { + hEtaPhiHitsSecTrkLayer0->Fill(genPhi, genEta); + hEtaPtHitsSecTrkLayer0->Fill(genEta, genPt); + } + else if (hitLayer == 1 && isPrimary) { + hEtaPhiHitsPrmTrkLayer1->Fill(genPhi, genEta); + hEtaPtHitsPrmTrkLayer1->Fill(genEta, genPt); + } + else { + hEtaPhiHitsSecTrkLayer1->Fill(genPhi, genEta); + hEtaPtHitsSecTrkLayer1->Fill(genEta, genPt); + } } } - // 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()) + TH2F* hEtaPhiClsPrmTrkLayer0 = new TH2F("hEtaPhiClsPrmTrkLayer0", "hEtaPhiClsPrmTrkLayer0;#phi;#eta", 64, 0, 6.28319, 40, -2, 2); + TH2F* hEtaPhiClsSecTrkLayer0 = new TH2F("hEtaPhiClsSecTrkLayer0", "hEtaPhiClsSecTrkLayer0;#phi;#eta", 64, 0, 6.28319, 40, -2, 2); + TH2F* hEtaPhiClsPrmTrkLayer1 = new TH2F("hEtaPhiClsPrmTrkLayer1", "hEtaPhiClsPrmTrkLayer1;#phi;#eta", 64, 0, 6.28319, 40, -2, 2); + TH2F* hEtaPhiClsSecTrkLayer1 = new TH2F("hEtaPhiClsSecTrkLayer1", "hEtaPhiClsSecTrkLayer1;#phi;#eta", 64, 0, 6.28319, 40, -2, 2); + for (int iCls = 0; iCls < clustersArray->size(); ++iCls) { + const auto& cls = (*clustersArray)[iCls]; + const auto& clsLabels = clustersLabelsArr->getLabels(iCls); + if (clsLabels.empty()) continue; - const int eventID = evtTrackLabel.getEventID(); - const int trackID = evtTrackLabel.getTrackID(); + // Link the cluster to the primary track + int clsEventID{-1}, clsTrackID{-1}; + bool hasValidLabels{false}; + int nAssoc{0}; + for (int iLabel=0; iLabel= nEvts) { - std::cerr << "WARNING: digit " << iDigit << " has invalid eventID=" << eventID << "\n"; - continue; - } + const auto& label = clsLabels[iLabel]; + if (!label.isValid()) + continue; + hasValidLabels = true; - 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); - } - } + int iEvtID = label.getEventID(); + int iTrkID = label.getTrackID(); - // 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) { + auto &iMcTrack = (*mcTracksPerEvent[iEvtID])[iTrkID]; - const auto& cls = (*clustersArray)[iCls]; - const auto& clsLabels = clustersLabelsArr->getLabels(iCls); + // Do not update the track label of the cluster + // if there are multiple labels and at least one + // of them is a primary track + if (!iMcTrack.isPrimary() && nAssoc > 0) { + continue; + } + clsEventID = iEvtID; + clsTrackID = iTrkID; + nAssoc++; + } - if (clsLabels.empty()) + if (!hasValidLabels) { + std::cerr << "WARNING: cluster " << iCls << " has no valid labels\n"; continue; - const auto& evtTrackLabel = clsLabels[0]; - if (!evtTrackLabel.isValid()) + } + + if (clsEventID < 0 || clsEventID >= nEvts) { + std::cerr << "WARNING: cluster " << iCls << " has invalid eventID=" << clsEventID << "\n"; continue; + } - const int eventID = evtTrackLabel.getEventID(); - const int trackID = evtTrackLabel.getTrackID(); + const auto& mcTrack = (*mcTracksPerEvent[clsEventID])[clsTrackID]; + float genEta = mcTrack.GetEta(); + float genPhi = mcTrack.GetPhi(); + bool isPrimary = mcTrack.isPrimary(); - if (eventID < 0 || eventID >= nEvts) { - std::cerr << "WARNING: cluster " << iCls << " has invalid eventID=" << eventID << "\n"; + const uint64_t trackKey = (static_cast(clsEventID) << 32) | static_cast(clsTrackID); + if (tracksHitCls.find(trackKey) == tracksHitCls.end()) { continue; } - - // PrintCluster(verbose, cls, clsLabels, iotofGeom, segmInfo); - auto& hitList = allEvtsTrackData[eventID][trackID].hitsByDetector[cls.getChipID()]; - for (auto& hit : hitList) { - hit.assocClsIdxs.push_back(iCls); + int clsLayer = iotofGeom->getIOTOFLayer(cls.getChipID()); + if (clsLayer == 0) { + tracksHitCls[trackKey].clsIndicesL0.push_back(iCls); + } else if (clsLayer == 1) { + tracksHitCls[trackKey].clsIndicesL1.push_back(iCls); } - } - // 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); - } - } - } - } - } + if (clsLayer == 0 && isPrimary) { hEtaPhiClsPrmTrkLayer0->Fill(genPhi, genEta); } + else if (clsLayer == 0 && !isPrimary) { hEtaPhiClsSecTrkLayer0->Fill(genPhi, genEta); } + else if (clsLayer == 1 && isPrimary) { hEtaPhiClsPrmTrkLayer1->Fill(genPhi, genEta); } + else { hEtaPhiClsSecTrkLayer1->Fill(genPhi, genEta); } - // 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); @@ -482,45 +512,72 @@ void CheckClustersIOTOF(std::string kinefile = "o2sim_Kine.root", std::vector clustersProperties; clustersProperties.reserve(clustersArray->size()); // Pre-allocate memory + Print(verbose, "\n\n----> Starting clusters printouts ... "); for (int iCls = 0; iCls < (int)clustersArray->size(); ++iCls) { - const auto& cluster = (*clustersArray)[iCls]; - - // Cluster labels + const auto& cls = (*clustersArray)[iCls]; 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; + if (clsLabels.empty()) continue; + + // Update with primary track if multiple labels are present + int nAssocPrimaries{0}, nAssocSecondaries{0}; + int clsEventID{-1}, clsTrackID{-1}; + bool hasValidLabels{false}; + std::set uniqueEventIDs; + + for (int iLabel=0; iLabel 0) { + continue; + } + clsEventID = iEvtID; + clsTrackID = iTrkID; } - 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"; + if (clsEventID < 0 || clsEventID >= nEvts) { + std::cerr << "WARNING: cluster " << iCls << " has invalid eventID=" << clsEventID << "\n"; continue; } + const auto& mcTrack = (*mcTracksPerEvent[clsEventID])[clsTrackID]; + 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()); + clsProps.chipID = cls.getChipID(); + clsProps.layer = iotofGeom->getIOTOFLayer(cls.getChipID()); + clsProps.rowStart = cls.getRow(); + clsProps.rowSpan = cls.getRowSpan(); + clsProps.colStart = cls.getCol(); + clsProps.colSpan = cls.getColSpan(); + clsProps.pattern = cls.getPattern(); + clsProps.size = cls.getSize(); + clsProps.topology = static_cast(cls.getTopology()); uint32_t clsTopoKey = (static_cast(clsProps.rowSpan) << 24) | (static_cast(clsProps.colSpan) << 16) | static_cast(clsProps.pattern); @@ -528,95 +585,125 @@ void CheckClustersIOTOF(std::string kinefile = "o2sim_Kine.root", 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.eventID = clsEventID; + clsProps.trackID = clsTrackID; + clsProps.nAssocPrimaries = nAssocPrimaries; + clsProps.isShared = (nAssocPrimaries > 1); + clsProps.isFake = (uniqueEventIDs.size() > 1); clsProps.hitIdx = -1; + clsProps.isPrimary = mcTrack.isPrimary(); + float genPt = mcTrack.GetPt(); - // 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; - } + uint64_t trackKey = (static_cast(clsEventID) << 32) | static_cast(clsTrackID); - 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; - } - } + // Get hits in the chip associated to the cluster's track and event + std::vector chipHitsIdxs; + if (clsProps.layer == 0) { + chipHitsIdxs = tracksHitCls[trackKey].hitIndicesL0; + } else if (clsProps.layer == 1) { + chipHitsIdxs = tracksHitCls[trackKey].hitIndicesL1; } - - // 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()); - } - - 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; - } + if (chipHitsIdxs.size() == 0) { + std::cout << "No hits by this track and event in this chip: " << clsProps.chipID << std::endl; + continue; + } if (chipHitsIdxs.size() == 1) { + clsProps.hitIdx = 0; + hCountHitMatchingType->Fill(clsProps.isPrimary ? 0 : 2, genPt); + } else { + // Perform spatial matching for chips with multiple hits + clsProps.hitIdx = FindBestMatchingHit(cls, clsTopoInfo, chipHitsIdxs, hitsPerEvent[clsProps.eventID], iotofGeom, segmInfo); + hCountHitMatchingType->Fill(clsProps.isPrimary ? 1 : 3, genPt); } + // // Print cluster information // PrintCluster(verbose, cluster, digitsArray, digitsLabels, hitsPerEvent, iotofGeom, segmInfo); clustersProperties.push_back(clsProps); + } Print(true, "----> Total number of clusters: ", clustersProperties.size()); + // Print features of all clusters and hits of the particles in the events + for (const auto& trackProperties : tracksHitCls) { + const uint64_t trackKey = trackProperties.first; + const auto eventID = static_cast(trackKey >> 32); + const auto trackID = static_cast(trackKey & 0xFFFFFFFF); + Print(verbose, "\n\n"); + PrintMcTrack(verbose, (*mcTracksPerEvent[eventID])[trackID]); + + Print(verbose, "Layer 0"); + for (const auto& hitIdx : trackProperties.second.hitIndicesL0) { + const auto& hit = (*hitsPerEvent[eventID])[hitIdx]; + PrintHit(verbose, hit, iotofGeom); + } + for (const auto& clsIdx : trackProperties.second.clsIndicesL0) { + const auto& cls = (*clustersArray)[clsIdx]; + const auto& clsLabels = clustersLabelsArr->getLabels(clsIdx); + uint32_t clsTopoKey = (static_cast(cls.getRowSpan()) << 24) | + (static_cast(cls.getColSpan()) << 16) | + static_cast(cls.getPattern()); + TopologyInfo clsTopoInfo = topoClassifier.getTopologyFeatures(clsTopoKey); + PrintCluster(verbose, cls, clsLabels, clsTopoInfo, iotofGeom, segmInfo); + } + Print(verbose, "Layer 1"); + for (const auto& hitIdx : trackProperties.second.hitIndicesL1) { + const auto& hit = (*hitsPerEvent[eventID])[hitIdx]; + PrintHit(verbose, hit, iotofGeom); + } + for (const auto& clsIdx : trackProperties.second.clsIndicesL1) { + const auto& cls = (*clustersArray)[clsIdx]; + const auto& clsLabels = clustersLabelsArr->getLabels(clsIdx); + uint32_t clsTopoKey = (static_cast(cls.getRowSpan()) << 24) | + (static_cast(cls.getColSpan()) << 16) | + static_cast(cls.getPattern()); + TopologyInfo clsTopoInfo = topoClassifier.getTopologyFeatures(clsTopoKey); + PrintCluster(verbose, cls, clsLabels, clsTopoInfo, 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 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; + // QA printouts and histograms Print(true, "\n\n----> Starting QA logging ... "); const char* trackName[2] = {"Prm", "Sec"}; + // Output + TFile* outFile = new TFile("CheckClusters.root", "RECREATE"); + for (int type = 0; type < 2; ++type) { + hGenEtaPt[type]->Write(); + } + + hEtaPhiHitsPrmTrkLayer0->Write(); + hEtaPhiHitsSecTrkLayer0->Write(); + hEtaPhiHitsPrmTrkLayer1->Write(); + hEtaPhiHitsSecTrkLayer1->Write(); + hEtaPtHitsPrmTrkLayer0->Write(); + hEtaPtHitsSecTrkLayer0->Write(); + hEtaPtHitsPrmTrkLayer1->Write(); + hEtaPtHitsSecTrkLayer1->Write(); + hEtaPhiClsPrmTrkLayer0->Write(); + hEtaPhiClsSecTrkLayer0->Write(); + hEtaPhiClsPrmTrkLayer1->Write(); + hEtaPhiClsSecTrkLayer1->Write(); + + // Count fake clusters - TH1F* hCountFakeClusters[2][2]; + TH1F* hCountClsTypes[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"); + hCountClsTypes[layer][type] = new TH1F(Form("hCountClsTypes%sTrkLayer%d", trackName[type], layer), Form("Cluster Counter %s Trk Layer %d", trackName[type], layer), 4, -0.5, 3.5); + hCountClsTypes[layer][type]->GetXaxis()->SetBinLabel(1, "Total"); + hCountClsTypes[layer][type]->GetXaxis()->SetBinLabel(2, "Real"); + hCountClsTypes[layer][type]->GetXaxis()->SetBinLabel(3, "Shared"); + hCountClsTypes[layer][type]->GetXaxis()->SetBinLabel(4, "Fake"); } } @@ -624,30 +711,23 @@ void CheckClustersIOTOF(std::string kinefile = "o2sim_Kine.root", for (const auto& cluster : clustersProperties) { int layer = cluster.layer; int type = cluster.isPrimary ? 0 : 1; - hCountFakeClusters[layer][type]->Fill(0.f, 1); // Total clusters + hCountClsTypes[layer][type]->Fill(0.f, 1); // Total clusters + if (cluster.isShared) { + hCountClsTypes[layer][type]->Fill(2.f, 1); // Shared clusters + continue; + } 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 + hCountClsTypes[layer][type]->Fill(3.f, 1); // Fake clusters + continue; } + hCountClsTypes[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"}; + "kSingleDigit", "kLineOnRow", "kLineOnCol", "kSquare", "kRectangle", "kDiagonal", + "kLowerTriangleLeft", "kLowerTriangleRight", "kUpperTriangleLeft", "kUpperTriangleRight", + "kSnake", "kSnakeRefl", "kSnakeRot90", "kSnakeRot90Refl", "kHuge", "kOther"}; // Count topologies from frequency values in // topologies dictionary and fill the summary histograms @@ -656,40 +736,50 @@ void CheckClustersIOTOF(std::string kinefile = "o2sim_Kine.root", hTopoSummaryDictionary->Fill(topology.mTopology, topology.mFrequency); } - 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]; + TH2F *hTrueClsSizeVsEta[2][2], *hTrueClsSizeVsPhi[2][2], *hFakeClsSizeVsEta[2][2], *hFakeClsSizeVsPhi[2][2], + *hClustersEtaPhi[2][2], *hTopoVsEta[2][2], *hClsSizeVsTopo[2][2], + *hXResVsEta[2][2], *hYResVsEta[2][2], *hZResVsEta[2][2], *hXResVsTopo[2][2], *hYResVsTopo[2][2], *hZResVsTopo[2][2], + *hTrackHitsXY[2][2], *hTrackDoubleHitsXY[2][2], *hTrackDoubleHitsPhiPt[2][2], *hTopoVsEtaPt[2][2][kNTopologies], + *hRecoClsEtaPt[2][2]; + TH1F *hMeanTrueClsSizeVsEta[2][2], *hMeanTrueClsSizeVsPhi[2][2], *hMeanFakeClsSizeVsEta[2][2], *hMeanFakeClsSizeVsPhi[2][2], + *hRmsXResVsEta[2][2], *hRmsYResVsEta[2][2], *hRmsZResVsEta[2][2], *hMeanXResVsEta[2][2], *hMeanYResVsEta[2][2], *hMeanZResVsEta[2][2], + *hRmsXResVsTopo[2][2], *hRmsYResVsTopo[2][2], *hRmsZResVsTopo[2][2], *hMeanXResVsTopo[2][2], *hMeanYResVsTopo[2][2], *hMeanZResVsTopo[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); - 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); + hClustersEtaPhi[layer][type] = new TH2F(Form("hNClsVsEtaPhi%sTrkLayer%d", trackName[type], layer), "Cluster #eta vs #phi;#phi;#eta", 64, 0, 6.28319, 40, -2, 2); + hTrueClsSizeVsEta[layer][type] = new TH2F(Form("hTrueClsSizeVsEta%sTrkLayer%d", trackName[type], layer), "True Cluster Size vs #eta;#eta", 40, -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", 64, 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", 40, -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", 64, 0, 6.28319, 20, 0.5, 20.5); + hMeanTrueClsSizeVsEta[layer][type] = new TH1F(Form("hMeanTrueClsSizeVsEta%sTrkLayer%d", trackName[type], layer), "Mean True Cluster Size vs #eta;#eta", 40, -2, 2); + hMeanTrueClsSizeVsPhi[layer][type] = new TH1F(Form("hMeanTrueClsSizeVsPhi%sTrkLayer%d", trackName[type], layer), "Mean True Cluster Size vs #phi;#phi", 64, 0, 6.28319); + hMeanFakeClsSizeVsEta[layer][type] = new TH1F(Form("hMeanFakeClsSizeVsEta%sTrkLayer%d", trackName[type], layer), "Mean Fake Cluster Size vs #eta;#eta", 40, -2, 2); + hMeanFakeClsSizeVsPhi[layer][type] = new TH1F(Form("hMeanFakeClsSizeVsPhi%sTrkLayer%d", trackName[type], layer), "Mean Fake Cluster Size vs #phi;#phi", 64, 0, 6.28319); + hRecoClsEtaPt[layer][type] = new TH2F(Form("hRecoClsEtaPt%sTrkLayer%d", trackName[type], layer), "Reconstructed Cluster vs p_{T};#eta;p_{T}", 40, -2, 2, 50, 0, 10); + hTopoVsEta[layer][type] = new TH2F(Form("hClsSizeVsEtaTopo%sTrkLayer%d", trackName[type], layer), "Cluster Topology vs #eta;;#eta", kNTopologies, 0, kNTopologies, 40, -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); + hXResVsEta[layer][type] = new TH2F(Form("hDeltaXClsHit%sTrkLayer%d", trackName[type], layer), ";#Delta X;#eta", 1000, -2, 2, 40, -2, 2); + hYResVsEta[layer][type] = new TH2F(Form("hDeltaYClsHit%sTrkLayer%d", trackName[type], layer), ";#Delta Y;#eta", 1000, -2, 2, 40, -2, 2); + hZResVsEta[layer][type] = new TH2F(Form("hDeltaZClsHit%sTrkLayer%d", trackName[type], layer), ";#Delta Z;#eta", 1000, -2, 2, 40, -2, 2); + hRmsXResVsEta[layer][type] = new TH1F(Form("hRmsDeltaXClsHit%sTrkLayer%d", trackName[type], layer), ";#eta;RMS #Delta X", 40, -2, 2); + hRmsYResVsEta[layer][type] = new TH1F(Form("hRmsDeltaYClsHit%sTrkLayer%d", trackName[type], layer), ";#eta;RMS #Delta Y", 40, -2, 2); + hRmsZResVsEta[layer][type] = new TH1F(Form("hRmsDeltaZClsHit%sTrkLayer%d", trackName[type], layer), ";#eta;RMS #Delta Z", 40, -2, 2); + hMeanXResVsEta[layer][type] = new TH1F(Form("hMeanDeltaXClsHit%sTrkLayer%d", trackName[type], layer), ";#eta;Mean #Delta X", 40, -2, 2); + hMeanYResVsEta[layer][type] = new TH1F(Form("hMeanDeltaYClsHit%sTrkLayer%d", trackName[type], layer), ";#eta;Mean #Delta Y", 40, -2, 2); + hMeanZResVsEta[layer][type] = new TH1F(Form("hMeanDeltaZClsHit%sTrkLayer%d", trackName[type], layer), ";#eta;Mean #Delta Z", 40, -2, 2); + hXResVsTopo[layer][type] = new TH2F(Form("hDeltaXClsHit%sTrkLayer%d", trackName[type], layer), ";#Delta X;Cluster Topology", 1000, -2, 2, kNTopologies, 0, kNTopologies); + hYResVsTopo[layer][type] = new TH2F(Form("hDeltaYClsHit%sTrkLayer%d", trackName[type], layer), ";#Delta Y;Cluster Topology", 1000, -2, 2, kNTopologies, 0, kNTopologies); + hZResVsTopo[layer][type] = new TH2F(Form("hDeltaZClsHit%sTrkLayer%d", trackName[type], layer), ";#Delta Z;Cluster Topology", 1000, -2, 2, kNTopologies, 0, kNTopologies); + hRmsXResVsTopo[layer][type] = new TH1F(Form("hRmsDeltaXClsHit%sTrkLayer%d", trackName[type], layer), ";Cluster Topology;RMS #Delta X", kNTopologies, 0, kNTopologies); + hRmsYResVsTopo[layer][type] = new TH1F(Form("hRmsDeltaYClsHit%sTrkLayer%d", trackName[type], layer), ";Cluster Topology;RMS #Delta Y", kNTopologies, 0, kNTopologies); + hRmsZResVsTopo[layer][type] = new TH1F(Form("hRmsDeltaZClsHit%sTrkLayer%d", trackName[type], layer), ";Cluster Topology;RMS #Delta Z", kNTopologies, 0, kNTopologies); + hMeanXResVsTopo[layer][type] = new TH1F(Form("hMeanDeltaXClsHit%sTrkLayer%d", trackName[type], layer), ";Cluster Topology;Mean #Delta X", kNTopologies, 0, kNTopologies); + hMeanYResVsTopo[layer][type] = new TH1F(Form("hMeanDeltaYClsHit%sTrkLayer%d", trackName[type], layer), ";Cluster Topology;Mean #Delta Y", kNTopologies, 0, kNTopologies); + hMeanZResVsTopo[layer][type] = new TH1F(Form("hMeanDeltaZClsHit%sTrkLayer%d", trackName[type], layer), ";Cluster Topology;Mean #Delta Z", kNTopologies, 0, kNTopologies); 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); @@ -707,18 +797,140 @@ void CheckClustersIOTOF(std::string kinefile = "o2sim_Kine.root", 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); + 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()), 40, -2, 2, 20, 0, 10); } } } + // Check digit efficiency across pixel by print the local coordinates + // of hits without any cluster and digit associated to them + TH2F* hNotRecoHits[2][2]; + for (int l = 0; l < 2; ++l) { + for (int t = 0; t < 2; ++t) { + hNotRecoHits[l][t] = new TH2F( + Form("hNotRecoHits_layer%d_type%d", l, t), + Form("Local Position Unreconstructed Hits L%d Type%d;x (cm);y (cm)", l, t), + 1000, -2., 2., // Adjust binning/ranges to your sensor dimensions + 1000, -2., 2. + ); + } + } + + for (const auto& trackProperties : tracksHitCls) { + int eventID = static_cast(trackProperties.first >> 32); + int trackID = static_cast(trackProperties.first & 0xFFFFFFFF); + const auto& mcTrack = (*mcTracksPerEvent[eventID])[trackID]; + + if (!mcTrack.isPrimary()) { continue; } + const int type = 0; + + if (trackProperties.second.hitIndicesL0.empty() && trackProperties.second.hitIndicesL1.empty()) { + continue; // Skip tracks without hits in both layers + } + + // Helper lambda to check if two chips are adjacent + auto areChipsAdjacent = [&](int chipID1, int chipID2) -> bool { + int l1{-1}, s1{-1}, ss1{-1}, m1{-1}, c1{-1}; + int l2{-1}, s2{-1}, ss2{-1}, m2{-1}, c2{-1}; + + iotofGeom->getIOTOFChipId(chipID1, l1, s1, ss1, m1, c1); + iotofGeom->getIOTOFChipId(chipID2, l2, s2, ss2, m2, c2); + + bool sameLayer = (l1 == l2); + bool sameStave = (s1 == s2); + bool sameSubStave = (ss1 == ss2); + bool sameModule = (m1 == m2); + bool sameChip = (c1 == c2); + + // // Same chip check + // if (sameLayer && sameStave && sameSubStave && sameModule && sameChip) { + // return true; + // } + // Adjacent module check on same stave & substave + if (sameLayer && sameStave && sameSubStave && std::abs(m1 - m2) <= 1) { + return true; + } + // Adjacent chip check on same module + if (sameLayer && sameStave && sameSubStave && sameModule && std::abs(c1 - c2) <= 1) { + return true; + } + return false; + }; + + // Helper lambda to check a layer's cluster list + auto hasNonAdjacentDoubleClusters = [&](const std::vector& hitIndices) -> bool { + if (hitIndices.size() < 2) return false; + + // Collect unique chip IDs for this layer + std::vector chips; + for (int hitIdx : hitIndices) { + chips.push_back((*hitsPerEvent[eventID])[hitIdx].GetDetectorID()); + } + std::sort(chips.begin(), chips.end()); + chips.erase(std::unique(chips.begin(), chips.end()), chips.end()); + + if (chips.size() < 2) return false; // All clusters are on the exact same chip + + // Check if ANY pair of chips is non-adjacent + for (size_t i = 0; i < chips.size(); ++i) { + for (size_t j = i + 1; j < chips.size(); ++j) { + if (!areChipsAdjacent(chips[i], chips[j])) { + return true; // Found double clusters on non-adjacent chips! + } + } + } + return false; + }; + + // Evaluate for L0 and L1 directly using DetectorData + bool hasDoubleClustersL0 = hasNonAdjacentDoubleClusters(trackProperties.second.hitIndicesL0); + for (const auto& hitIdx : trackProperties.second.hitIndicesL0) { + const auto& hit = (*hitsPerEvent[eventID])[hitIdx]; + PrintHit(verbose, hit, iotofGeom); + hTrackHitsXY[0][type]->Fill(hit.GetX(), hit.GetY()); + + // Check for non-reconstructed hits + if (trackProperties.second.clsIndicesL0.empty()) { + o2::math_utils::Point3D avgPos; + GetHitAvgPositionLocal(hit, iotofGeom, avgPos); + hNotRecoHits[0][type]->Fill(avgPos.X(), avgPos.Z()); + } + if (hasDoubleClustersL0) { + hTrackDoubleHitsPhiPt[0][type]->Fill(mcTrack.GetPhi(), mcTrack.GetPt()); + if (mcTrack.GetPt() > 5.0f) { + hTrackDoubleHitsXY[0][type]->Fill(hit.GetX(), hit.GetY()); + } + } + } + bool hasDoubleClustersL1 = hasNonAdjacentDoubleClusters(trackProperties.second.hitIndicesL1); + for (const auto& hitIdx : trackProperties.second.hitIndicesL1) { + const auto& hit = (*hitsPerEvent[eventID])[hitIdx]; + PrintHit(verbose, hit, iotofGeom); + hTrackHitsXY[1][type]->Fill(hit.GetX(), hit.GetY()); + + // Check for non-reconstructed hits + if (trackProperties.second.clsIndicesL1.empty()) { + o2::math_utils::Point3D avgPos; + GetHitAvgPositionLocal(hit, iotofGeom, avgPos); + hNotRecoHits[1][type]->Fill(avgPos.X(), avgPos.Z()); + } + if (hasDoubleClustersL1) { + hTrackDoubleHitsPhiPt[1][type]->Fill(mcTrack.GetPhi(), mcTrack.GetPt()); + if (mcTrack.GetPt() > 5.0f) { + hTrackDoubleHitsXY[1][type]->Fill(hit.GetX(), hit.GetY()); + } + } + } + } // end event loop + // Loop over clusters Print(true, "----> Looping over clusters and filling histograms"); + std::cout << "\n\n\n\n\n\n\n" << std::endl; 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; @@ -729,10 +941,10 @@ void CheckClustersIOTOF(std::string kinefile = "o2sim_Kine.root", 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); @@ -740,30 +952,43 @@ void CheckClustersIOTOF(std::string kinefile = "o2sim_Kine.root", 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); + continue; // Skip clusters without a matching hit } - if (cls.hitIdx < 0) { - continue; // Skip clusters without a matching hit + hTopoSummaryReal->Fill(topo); + hTrueClsSizeVsEta[layer][type]->Fill(eta, size); + hTrueClsSizeVsPhi[layer][type]->Fill(phi, size); + float weight = cls.isShared ? cls.nAssocPrimaries : 1.0f; // Weight for shared clusters + hClustersEtaPhi[layer][type]->Fill(phi, eta, weight); + hRecoClsEtaPt[layer][type]->Fill(eta, pt, weight); + + int hitIdx{-1}; + uint64_t trackKey = (static_cast(eventID) << 32) | static_cast(trackID); + if (layer == 0) { + hitIdx = tracksHitCls[trackKey].hitIndicesL0[cls.hitIdx]; + } else if (layer == 1) { + hitIdx = tracksHitCls[trackKey].hitIndicesL1[cls.hitIdx]; } - 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); + auto& hit = (*hitsPerEvent[cls.eventID])[hitIdx]; o2::math_utils::Point3D clusterPos; TopologyInfo clsTopoInfo = topoClassifier.getTopologyFeatures(cls.topoKey); auto clsFull = clustersArray->at(cls.clsIdx); - GetClusterGlobalPos(clsFull, clsTopoInfo, clusterPos, iotofGeom, segmInfo); + GetClusterLocalPos(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); + GetHitAvgPositionLocal(hit, iotofGeom, avgPos); + // if (clusterPos.X() - avgPos.X() > 1) { + + // } + // PrintHit(true, hit, iotofGeom); + // std::cout << "Hit average position " << avgPos.X() << ", " << avgPos.Y() << ", " << avgPos.Z() << std::endl; + // std::cout << std::endl; + hXResVsEta[layer][type]->Fill(clusterPos.X() - avgPos.X(), eta); + hYResVsEta[layer][type]->Fill(clusterPos.Y() - avgPos.Y(), eta); + hZResVsEta[layer][type]->Fill(clusterPos.Z() - avgPos.Z(), eta); + hXResVsTopo[layer][type]->Fill(clusterPos.X() - avgPos.X(), topo); + hYResVsTopo[layer][type]->Fill(clusterPos.Y() - avgPos.Y(), topo); + hZResVsTopo[layer][type]->Fill(clusterPos.Z() - avgPos.Z(), topo); } // Fill means and RMS of cluster size and residuals @@ -793,147 +1018,43 @@ void CheckClustersIOTOF(std::string kinefile = "o2sim_Kine.root", 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()); + for (int etaBin = 1; etaBin <= hXResVsEta[layer][type]->GetNbinsY(); ++etaBin) { + TH1D* hXResProjVsEta = hXResVsEta[layer][type]->ProjectionX(Form("hXResProj_etaBin%d", etaBin), etaBin, etaBin); + TH1D* hYResProjVsEta = hYResVsEta[layer][type]->ProjectionX(Form("hYResProj_etaBin%d", etaBin), etaBin, etaBin); + TH1D* hZResProjVsEta = hZResVsEta[layer][type]->ProjectionX(Form("hZResProj_etaBin%d", etaBin), etaBin, etaBin); + hRmsXResVsEta[layer][type]->SetBinContent(etaBin, hXResProjVsEta->GetRMS()); + hRmsYResVsEta[layer][type]->SetBinContent(etaBin, hYResProjVsEta->GetRMS()); + hRmsZResVsEta[layer][type]->SetBinContent(etaBin, hZResProjVsEta->GetRMS()); + hRmsXResVsEta[layer][type]->SetBinError(etaBin, hXResProjVsEta->GetRMSError()); + hRmsYResVsEta[layer][type]->SetBinError(etaBin, hYResProjVsEta->GetRMSError()); + hRmsZResVsEta[layer][type]->SetBinError(etaBin, hZResProjVsEta->GetRMSError()); + hMeanXResVsEta[layer][type]->SetBinContent(etaBin, hXResProjVsEta->GetMean()); + hMeanYResVsEta[layer][type]->SetBinContent(etaBin, hYResProjVsEta->GetMean()); + hMeanZResVsEta[layer][type]->SetBinContent(etaBin, hZResProjVsEta->GetMean()); + hMeanXResVsEta[layer][type]->SetBinError(etaBin, hXResProjVsEta->GetMeanError()); + hMeanYResVsEta[layer][type]->SetBinError(etaBin, hYResProjVsEta->GetMeanError()); + hMeanZResVsEta[layer][type]->SetBinError(etaBin, hZResProjVsEta->GetMeanError()); } - } - } - - Print(true, "----> Looping over generated particles"); - - // 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()); - } - } - - // 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]) { - - const auto& mcTrack = (*mcTracksPerEvent[iEvt])[trackID]; - if (!mcTrack.isPrimary() || trackData.hitsByDetector.size() <= 1) { - continue; - } - - // 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()); - } - } + for (int topoBin = 1; topoBin <= hXResVsTopo[layer][type]->GetNbinsY(); ++topoBin) { + TH1D* hXResProjVsTopo = hXResVsTopo[layer][type]->ProjectionX(Form("hXResProj_topoBin%d", topoBin), topoBin, topoBin); + TH1D* hYResProjVsTopo = hYResVsTopo[layer][type]->ProjectionX(Form("hYResProj_topoBin%d", topoBin), topoBin, topoBin); + TH1D* hZResProjVsTopo = hZResVsTopo[layer][type]->ProjectionX(Form("hZResProj_topoBin%d", topoBin), topoBin, topoBin); + hRmsXResVsTopo[layer][type]->SetBinContent(topoBin, hXResProjVsTopo->GetRMS()); + hRmsYResVsTopo[layer][type]->SetBinContent(topoBin, hYResProjVsTopo->GetRMS()); + hRmsZResVsTopo[layer][type]->SetBinContent(topoBin, hZResProjVsTopo->GetRMS()); + hRmsXResVsTopo[layer][type]->SetBinError(topoBin, hXResProjVsTopo->GetRMSError()); + hRmsYResVsTopo[layer][type]->SetBinError(topoBin, hYResProjVsTopo->GetRMSError()); + hRmsZResVsTopo[layer][type]->SetBinError(topoBin, hZResProjVsTopo->GetRMSError()); + hMeanXResVsTopo[layer][type]->SetBinContent(topoBin, hXResProjVsTopo->GetMean()); + hMeanYResVsTopo[layer][type]->SetBinContent(topoBin, hYResProjVsTopo->GetMean()); + hMeanZResVsTopo[layer][type]->SetBinContent(topoBin, hZResProjVsTopo->GetMean()); + hMeanXResVsTopo[layer][type]->SetBinError(topoBin, hXResProjVsTopo->GetMeanError()); + hMeanYResVsTopo[layer][type]->SetBinError(topoBin, hYResProjVsTopo->GetMeanError()); + hMeanZResVsTopo[layer][type]->SetBinError(topoBin, hZResProjVsTopo->GetMeanError()); } } } - Print(true, "----> Writing histograms"); - // Output - TFile* outFile = new TFile("CheckClusters.root", "RECREATE"); - for (int type = 0; type < 2; ++type) { - hGenEtaPt[type]->Write(); - } - - hEtaPhiHitsPrmTrkLayer0->Write(); - hEtaPhiHitsSecTrkLayer0->Write(); - hEtaPhiHitsPrmTrkLayer1->Write(); - hEtaPhiHitsSecTrkLayer1->Write(); hTopoSummaryReal->Write(); hTopoSummaryFake->Write(); hTopoSummaryTotal->Write(); @@ -947,7 +1068,7 @@ void CheckClustersIOTOF(std::string kinefile = "o2sim_Kine.root", outFile->mkdir(Form("%sTrkLayer%d/Topologies", trackName[type], layer)); outFile->cd(Form("%sTrkLayer%d", trackName[type], layer)); - hCountFakeClusters[layer][type]->Write("hCountFakeClusters"); + hCountClsTypes[layer][type]->Write("hCountClsTypes"); hClustersEtaPhi[layer][type]->Write("hClustersEtaPhi"); hTrueClsSizeVsEta[layer][type]->Write("hTrueClsSizeVsEta"); @@ -955,29 +1076,110 @@ void CheckClustersIOTOF(std::string kinefile = "o2sim_Kine.root", 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; + TH2F* hEffEtaPhi = static_cast(hClustersEtaPhi[layer][type]->Clone(Form("hClsEffVsEtaPhi%sTrkLayer%d", trackName[type], layer))); + TH2F* hEtaPhiHits = layer == 0 ? (type == 0 ? hEtaPhiHitsPrmTrkLayer0 : hEtaPhiHitsSecTrkLayer0) + : (type == 0 ? hEtaPhiHitsPrmTrkLayer1 : hEtaPhiHitsSecTrkLayer1); + hEffEtaPhi->Divide(hEtaPhiHits); + hEffEtaPhi->Write("hClsEffEtaPhi"); + delete hEffEtaPhi; + + TH2F* hEffEtaPt = static_cast(hRecoClsEtaPt[layer][type]->Clone(Form("hClsEffVsEtaPhi%sTrkLayer%d", trackName[type], layer))); + TH2F* hEtaPtHits = layer == 0 ? (type == 0 ? hEtaPtHitsPrmTrkLayer0 : hEtaPtHitsSecTrkLayer0) + : (type == 0 ? hEtaPtHitsPrmTrkLayer1 : hEtaPtHitsSecTrkLayer1); + hEffEtaPt->Divide(hEtaPtHits); + hEffEtaPt->Write("hClsEffEtaPt"); + delete hEffEtaPt; + + // Compute Efficiency vs Eta + TH1D* hClustersEta = hClustersEtaPhi[layer][type]->ProjectionY(Form("hClsEta_%sTrkLayer%d", trackName[type], layer)); + TH1D* hHitsEta = hEtaPhiHits->ProjectionY(Form("hHitsEta_%sTrkLayer%d", trackName[type], layer)); + + TH1F* hEffVsEta = static_cast(hClustersEta->Clone(Form("hClsEffVsEta%sTrkLayer%d", trackName[type], layer))); + hEffVsEta->Divide(hHitsEta); + + // Compute proper binomial uncertainties + for (int bin = 1; bin <= hEffVsEta->GetNbinsX(); ++bin) { + double eff = hEffVsEta->GetBinContent(bin); + double nHits = hHitsEta->GetBinContent(bin); + + if (nHits > 0) { + // Clamp eff between 0 and 1 to prevent sqrt of negative numbers due to numerical precision + eff = std::clamp(eff, 0.0, 1.0); + double err = std::sqrt(eff * (1.0 - eff) / nHits); + hEffVsEta->SetBinError(bin, err); + } else { + hEffVsEta->SetBinError(bin, 0); + } + } + hEffVsEta->Write("hClsEffVsEta"); + + TH1D* hClustersPt = hRecoClsEtaPt[layer][type]->ProjectionY(Form("hClsPt_%sTrkLayer%d", trackName[type], layer)); + TH1D* hHitsPt = hEtaPtHits->ProjectionY(Form("hHitsPt_%sTrkLayer%d", trackName[type], layer)); + + TH1F* hEffVsPt = static_cast(hClustersPt->Clone(Form("hClsEffVsPt%sTrkLayer%d", trackName[type], layer))); + hEffVsPt->Divide(hHitsPt); + // Compute proper binomial uncertainties + for (int bin = 1; bin <= hEffVsPt->GetNbinsX(); ++bin) { + double eff = hEffVsPt->GetBinContent(bin); + double nHits = hHitsPt->GetBinContent(bin); + + if (nHits > 0) { + // Clamp eff between 0 and 1 to prevent sqrt of negative numbers due to numerical precision + eff = std::clamp(eff, 0.0, 1.0); + double err = std::sqrt(eff * (1.0 - eff) / nHits); + hEffVsPt->SetBinError(bin, err); + } else { + hEffVsPt->SetBinError(bin, 0); + } + } + hEffVsPt->Write("hClsEffVsPt"); + + // Compute Efficiency vs Phi + TH1D* hClustersPhi = hClustersEtaPhi[layer][type]->ProjectionX(Form("hClsPhi_%sTrkLayer%d", trackName[type], layer)); + TH1D* hHitsPhi = hEtaPhiHits->ProjectionX(Form("hHitsPhi_%sTrkLayer%d", trackName[type], layer)); + + TH1F* hEffVsPhi = static_cast(hClustersPhi->Clone(Form("hClsEffVsPhi%sTrkLayer%d", trackName[type], layer))); + hEffVsPhi->Divide(hHitsPhi); + + // Compute proper binomial uncertainties + for (int bin = 1; bin <= hEffVsPhi->GetNbinsX(); ++bin) { + double eff = hEffVsPhi->GetBinContent(bin); + double nHits = hHitsPhi->GetBinContent(bin); + + if (nHits > 0) { + eff = std::clamp(eff, 0.0, 1.0); + double err = std::sqrt(eff * (1.0 - eff) / nHits); + hEffVsPhi->SetBinError(bin, err); + } else { + hEffVsPhi->SetBinError(bin, 0); + } + } + hEffVsPhi->Write("hClsEffVsPhi"); - 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"); + hXResVsEta[layer][type]->Write("hXResVsEta"); + hYResVsEta[layer][type]->Write("hYResVsEta"); + hZResVsEta[layer][type]->Write("hZResVsEta"); + hRmsXResVsEta[layer][type]->Write("hRmsXResVsEta"); + hRmsYResVsEta[layer][type]->Write("hRmsYResVsEta"); + hRmsZResVsEta[layer][type]->Write("hRmsZResVsEta"); + hMeanXResVsEta[layer][type]->Write("hMeanXResVsEta"); + hMeanYResVsEta[layer][type]->Write("hMeanYResVsEta"); + hMeanZResVsEta[layer][type]->Write("hMeanZResVsEta"); + hXResVsTopo[layer][type]->Write("hXResVsTopo"); + hYResVsTopo[layer][type]->Write("hYResVsTopo"); + hZResVsTopo[layer][type]->Write("hZResVsTopo"); + hRmsXResVsTopo[layer][type]->Write("hRmsXResVsTopo"); + hRmsYResVsTopo[layer][type]->Write("hRmsYResVsTopo"); + hRmsZResVsTopo[layer][type]->Write("hRmsZResVsTopo"); + hMeanXResVsTopo[layer][type]->Write("hMeanXResVsTopo"); + hMeanYResVsTopo[layer][type]->Write("hMeanYResVsTopo"); + hMeanZResVsTopo[layer][type]->Write("hMeanZResVsTopo"); if (type == 0) { hTrackHitsXY[layer][type]->Write("hTrackHitsXY"); @@ -1041,59 +1243,118 @@ void CheckClustersIOTOF(std::string kinefile = "o2sim_Kine.root", } } - - // 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]; + // Write digit efficiency histograms 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); + outFile->cd(Form("%sTrkLayer%d", trackName[type], layer)); + hNotRecoHits[layer][type]->Write(); } } - 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 + // Map all found cluster topologies to histograms + // Create directories of all topologies + for (const auto& topoName : topologyNames) { + outFile->mkdir(Form("TopologyDictionary/All/%s", topoName.c_str())); + outFile->mkdir(Form("TopologyDictionary/Real/%s", topoName.c_str())); + outFile->mkdir(Form("TopologyDictionary/Fake/%s", topoName.c_str())); + } + for (int iMapEntry = 0; iMapEntry < sortedTopoMap.size(); ++iMapEntry) { + const auto& [topoKey, topology] = sortedTopoMap[iMapEntry]; + std::string topoName = topologyNames[topology.mTopology]; + int spanRow = topology.mSizeX; + int spanCol = topology.mSizeZ; + uint16_t bitmask = topology.mPattern; + + TH2F* hTopoDisplayAll = new TH2F(Form("spanRow_%i_spanCol_%i_key_%i_all", spanRow, spanCol, topoKey), Form("Cluster Topology %s;Row;Column", topoName.c_str()), + spanRow + 2, -1.5, spanRow + 0.5, spanCol + 2, -1.5, spanCol + 0.5); + TH2F* hTopoDisplayReal = new TH2F(Form("spanRow_%i_spanCol_%i_key_%i_real", spanRow, spanCol, topoKey), Form("Cluster Topology %s;Row;Column", topoName.c_str()), + spanRow + 2, -1.5, spanRow + 0.5, spanCol + 2, -1.5, spanCol + 0.5); + TH2F* hTopoDisplayFake = new TH2F(Form("spanRow_%i_spanCol_%i_key_%i_fake", spanRow, spanCol, topoKey), Form("Cluster Topology %s;Row;Column", topoName.c_str()), + spanRow + 2, -1.5, spanRow + 0.5, spanCol + 2, -1.5, spanCol + 0.5); + + TH2F* hTopoCOGAll = new TH2F(Form("spanRow_%i_spanCol_%i_key_%i_all_COG", spanRow, spanCol, topoKey), Form("Cluster Topology %s;Row;Column", topoName.c_str()), + spanRow + 2, -1.5, spanRow + 0.5, spanCol + 2, -1.5, spanCol + 0.5); + TH2F* hTopoCOGReal = new TH2F(Form("spanRow_%i_spanCol_%i_key_%i_real_COG", spanRow, spanCol, topoKey), Form("Cluster Topology %s;Row;Column", topoName.c_str()), + spanRow + 2, -1.5, spanRow + 0.5, spanCol + 2, -1.5, spanCol + 0.5); + TH2F* hTopoCOGFake = new TH2F(Form("spanRow_%i_spanCol_%i_key_%i_fake_COG", spanRow, spanCol, topoKey), Form("Cluster Topology %s;Row;Column", topoName.c_str()), + spanRow + 2, -1.5, spanRow + 0.5, spanCol + 2, -1.5, spanCol + 0.5); + + int frequency = topology.mFrequency; + hTopoCOGAll->Fill(topology.mOffsetXToCOG, topology.mOffsetZToCOG, frequency); + int countFakeThisTopo = std::count_if(clustersProperties.begin(), clustersProperties.end(), + [topoKey](const ClusterProperties& cls) + { return cls.topoKey == topoKey && cls.isFake; }); + hTopoCOGAll->Fill(topology.mOffsetXToCOG, topology.mOffsetZToCOG, countFakeThisTopo); + int countRealThisTopo = std::count_if(clustersProperties.begin(), clustersProperties.end(), + [topoKey](const ClusterProperties& cls) + { return cls.topoKey == topoKey && !cls.isFake; }); + hTopoCOGAll->Fill(topology.mOffsetXToCOG, topology.mOffsetZToCOG, countRealThisTopo); + + // Loop over the bits of bitmask and fill the histogram + for (int row = 0; row < spanRow; ++row) { + for (int col = 0; col < spanCol; ++col) { + int bitIndex = row * spanCol + col; + if (bitmask & (1 << bitIndex)) { + hTopoDisplayAll->Fill(row, col, frequency); + if (countRealThisTopo > 0) { + hTopoDisplayReal->Fill(row, col, countFakeThisTopo); } - 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()); + if (countFakeThisTopo > 0) { + hTopoDisplayFake->Fill(row, col, countRealThisTopo); } } } } - } - // 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->cd(Form("TopologyDictionary/All/%s", topoName.c_str())); + hTopoDisplayAll->Write(); + hTopoCOGAll->Write(); + if (countRealThisTopo > 0) { + outFile->cd(Form("TopologyDictionary/Real/%s", topoName.c_str())); + hTopoDisplayReal->Write(); + hTopoCOGReal->Write(); + } + if (countFakeThisTopo > 0) { + outFile->cd(Form("TopologyDictionary/Fake/%s", topoName.c_str())); + hTopoDisplayFake->Write(); + hTopoCOGFake->Write(); } + delete hTopoDisplayAll; + delete hTopoDisplayReal; + delete hTopoDisplayFake; + delete hTopoCOGAll; + delete hTopoCOGReal; + delete hTopoCOGFake; } outFile->Close(); delete outFile; + // Print all hits without any cluster associated to them + Print(verbose, "----> Printing all hits without any cluster associated to them"); + for (const auto& trackProperties : tracksHitCls) { + int eventID = static_cast(trackProperties.first >> 32); + int trackID = static_cast(trackProperties.first & 0xFFFFFFFF); + const auto& mcTrack = (*mcTracksPerEvent[eventID])[trackID]; + if (!mcTrack.isPrimary()) { continue; } - // // 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); - // } - // } + if (!(trackProperties.second.hitIndicesL0.size() > 0 && trackProperties.second.clsIndicesL0.empty()) || + !(trackProperties.second.hitIndicesL1.size() > 0 && trackProperties.second.clsIndicesL1.empty())) { + continue; // Skip tracks with reconstructed clusters + } + Print(verbose, "\nTrack has not reconstructed clusters for hits in both layers. Printing track and hit information:"); + PrintMcTrack(verbose, mcTrack); + if (trackProperties.second.clsIndicesL0.empty()) { + for (const auto& hitIdx : trackProperties.second.hitIndicesL0) { + const auto& hit = (*hitsPerEvent[eventID])[hitIdx]; + PrintHit(verbose, hit, iotofGeom); + } + } + if (trackProperties.second.clsIndicesL1.empty()) { + for (const auto& hitIdx : trackProperties.second.hitIndicesL1) { + const auto& hit = (*hitsPerEvent[eventID])[hitIdx]; + PrintHit(verbose, hit, iotofGeom); + } + } + } } diff --git a/Detectors/Upgrades/ALICE3/IOTOF/macros/CheckDigitsIOTOF.C b/Detectors/Upgrades/ALICE3/IOTOF/macros/CheckDigitsIOTOF.C index 581c93a236c98..65a921a814d8f 100644 --- a/Detectors/Upgrades/ALICE3/IOTOF/macros/CheckDigitsIOTOF.C +++ b/Detectors/Upgrades/ALICE3/IOTOF/macros/CheckDigitsIOTOF.C @@ -16,21 +16,27 @@ #include #include #include +#include #include #include #include #include #include +#include + #include "IOTOFBase/Segmentation.h" #include "IOTOFBase/IOTOFBaseParam.h" #include "IOTOFBase/GeometryTGeo.h" +#include "IOTOFSimulation/Digitizer.h" #include "DataFormatsIOTOF/Digit.h" #include "ITSMFTSimulation/Hit.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 "DetectorsBase/GeometryManager.h" #include "CCDB/BasicCCDBManager.h" @@ -75,7 +81,11 @@ 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 kinefile = "o2sim_Kine.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 +95,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(); @@ -123,15 +133,48 @@ void CheckDigitsIOTOF(std::string digifile = "tf3digits.root", std::string hitfi digTree->GetEntry(0); + // MC tracks + TFile* kineFile = TFile::Open(kinefile.data()); + TTree* kineTree = (TTree*)kineFile->Get("o2sim"); + std::vector*> mcTracksPerEvent(nevH, nullptr); + std::vector*> mcTracksRefsPerEvent(nevH, nullptr); + kineTree->SetBranchAddress("MCTrack", &mcTracksPerEvent[0]); + kineTree->SetBranchAddress("TrackRefs", &mcTracksRefsPerEvent[0]); + + TH1F* hGenHitsEta[2][2] = {{ + new TH1F("hGenHitsEtaPrmL0", "hGenHitsEtaPrmL0", 40, -2, 2), + new TH1F("hGenHitsEtaSecL0", "hGenHitsEtaSecL0", 40, -2, 2), + }, { + new TH1F("hGenHitsEtaPrmL1", "hGenHitsEtaPrmL1", 40, -2, 2), + new TH1F("hGenHitsEtaSecL1", "hGenHitsEtaSecL1", 40, -2, 2), + }}; + // Load all MC hit events upfront and build the hit lookup map. for (int im = 0; im < nevH; ++im) { hitTree->SetBranchAddress("TF3Hit", &hitArray[im]); hitTree->GetEntry(im); + kineTree->SetBranchAddress("MCTrack", &mcTracksPerEvent[im]); + kineTree->SetBranchAddress("TrackRefs", &mcTracksRefsPerEvent[im]); + kineTree->GetEntry(im); auto& mc2hit = mc2hitVec[im]; for (int ih = hitArray[im]->size(); ih--;) { const auto& hit = (*hitArray[im])[ih]; uint64_t key = (uint64_t(hit.GetTrackID()) << 32) + hit.GetDetectorID(); mc2hit.emplace(key, ih); + + auto &mcTrack = mcTracksPerEvent[im]->at(hit.GetTrackID()); + bool isPrimary = mcTrack.isPrimary(); + + int layer = gman->getIOTOFLayer(hit.GetDetectorID()); + if (layer == 0 && isPrimary) { + hGenHitsEta[0][0]->Fill(mcTrack.GetEta()); + } else if (layer == 0 && !isPrimary) { + hGenHitsEta[0][1]->Fill(mcTrack.GetEta()); + } else if (layer == 1 && isPrimary) { + hGenHitsEta[1][0]->Fill(mcTrack.GetEta()); + } else if (layer == 1 && !isPrimary) { + hGenHitsEta[1][1]->Fill(mcTrack.GetEta()); + } } } @@ -142,11 +185,21 @@ void CheckDigitsIOTOF(std::string digifile = "tf3digits.root", std::string hitfi plabelsArr->copyandflatten(labels); // LOOP on : ROFRecord array + TH1F* hRecoDigitEta[2][2] = {{ + new TH1F("hRecoDigitEtaPrmL0", "hRecoDigitEtaPrmL0", 40, -2, 2), + new TH1F("hRecoDigitEtaSecL0", "hRecoDigitEtaSecL0", 40, -2, 2), + }, { + new TH1F("hRecoDigitEtaPrmL1", "hRecoDigitEtaPrmL1", 40, -2, 2), + new TH1F("hRecoDigitEtaSecL1", "hRecoDigitEtaSecL1", 40, -2, 2), + }}; + + std::unordered_map> hitDigitMap; for (unsigned int iROF = 0; iROF < rofArr.size(); ++iROF) { const unsigned int rofIndex = rofArr[iROF].getFirstEntry(); const unsigned int rofNEntries = rofArr[iROF].getNEntries(); + std::unordered_map> tracksWithDigits; // LOOP on : digits array for (unsigned int iDigit = rofIndex; iDigit < rofIndex + rofNEntries; iDigit++) { if (iDigit % 1000 == 0) { @@ -176,10 +229,11 @@ void CheckDigitsIOTOF(std::string digifile = "tf3digits.root", std::string hitfi } int trID = lab.getTrackID(); + int evtID = lab.getEventID(); const auto gloD = gman->getMatrixL2G(chipID)(locD); // convert to global - std::unordered_map* mc2hit = &mc2hitVec[lab.getEventID()]; + std::unordered_map* mc2hit = &mc2hitVec[evtID]; // get MC info uint64_t key = (uint64_t(trID) << 32) + chipID; @@ -191,7 +245,7 @@ void CheckDigitsIOTOF(std::string digifile = "tf3digits.root", std::string hitfi } ////// HITS - Hit& hit = (*hitArray[lab.getEventID()])[hitEntry->second]; + Hit& hit = (*hitArray[evtID])[hitEntry->second]; auto xyzLocE = gman->getMatrixL2G(chipID) ^ (hit.GetPos()); // inverse conversion from global to local auto xyzLocS = gman->getMatrixL2G(chipID) ^ (hit.GetPosStart()); @@ -220,6 +274,21 @@ void CheckDigitsIOTOF(std::string digifile = "tf3digits.root", std::string hitfi locH.X() - locD.X(), locH.Z() - locD.Z()); /// difference in x and z between the hit and the digit in the local frame nt2->Fill(chipID, gloD.Z(), locHS.X() - locHE.X(), locHS.Z() - locHE.Z()); /// differences between local hit start and hit end positions + // Check if key is already in the set of tracks with digits, + // else we double count digits in efficiency calculation + // when using stepping + if (tracksWithDigits[evtID].find(key) == tracksWithDigits[evtID].end()) { + tracksWithDigits[evtID].insert(key); + int digitLayer = gman->getIOTOFLayer(chipID); + auto& mcTrack = mcTracksPerEvent[evtID]->at(trID); + bool isPrimary = mcTrack.isPrimary(); + hRecoDigitEta[digitLayer][isPrimary ? 0 : 1]->Fill(mcTrack.GetEta()); + } + + // Fill the hitDigitMap for later analysis + // Hit key from event ID and hit index + uint64_t hitKey = (uint64_t(evtID) << 32) + hitEntry->second; + hitDigitMap[hitKey].push_back(iDigit); } // end loop on digits array } // end loop on ROFRecords @@ -228,21 +297,21 @@ void CheckDigitsIOTOF(std::string digifile = "tf3digits.root", std::string hitfi auto canvXY = new TCanvas("canvXY", "", 1600, 800); canvXY->Divide(2, 1); canvXY->cd(1); - nt->Draw("y:x>>h_y_vs_x_IOTOF(1000, -100, 100, 1000, -100, 100)", "id >= 0 && id < 53568", "colz"); + nt->Draw("y:x>>h_y_vs_x_IOTOF(1000, -100, 100, 1000, -100, 100)", "id >= 0 && id < 55488", "colz"); canvXY->cd(2); - nt->Draw("y:z>>h_y_vs_z_IOTOF(1000, -400, 400, 1000, -100, 100)", "id >= 0 && id < 53568", "colz"); + nt->Draw("y:z>>h_y_vs_z_IOTOF(1000, -400, 400, 1000, -100, 100)", "id >= 0 && id < 55488", "colz"); canvXY->SaveAs("tf3digits_y_vs_x_vs_z.pdf"); // z distributions auto canvZ = new TCanvas("canvZ", "", 800, 800); canvZ->cd(); - nt->Draw("z>>h_z_IOTOF(500, -70, 70)", "id >= 0 && id < 53568 "); + nt->Draw("z>>h_z_IOTOF(500, -70, 70)", "id >= 0 && id < 55488 "); canvZ->SaveAs("tf3digits_z.pdf"); // dz distributions (difference between local position of digits and hits in x and z) auto canvdZ = new TCanvas("canvdZ", "", 800, 800); canvdZ->cd(); - nt->Draw("dz>>h_dz_ML(500, -0.05, 0.05)", "id >= 0 && id < 53568 "); + nt->Draw("dz>>h_dz_ML(500, -0.05, 0.05)", "id >= 0 && id < 55488 "); canvdZ->SaveAs("tf3digits_dz.pdf"); canvdZ->SaveAs("tf3digits_dz.root"); @@ -250,13 +319,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 < 55488", "colz"); addTLines(0.01); h = (TH2F*)gPad->GetPrimitive("h_dx_vs_dz_OTOF"); Info("OTOF", "RMS(dx)=%.1f mu", h->GetRMS(2) * 1e4); @@ -275,7 +344,7 @@ void CheckDigitsIOTOF(std::string digifile = "tf3digits.root", std::string hitfi Info("ITOF", "RMS(dxH)=%.1f mu", h->GetRMS(2) * 1e4); Info("ITOF", "RMS(dzH)=%.1f mu", h->GetRMS(1) * 1e4); canvdXdZHit->cd(2); - nt2->Draw("dxH:dzH>>h_dxH_vs_dzH_OTOF(300, -0.03, 0.03, 300, -0.03, 0.03)", "id >= 1920 && id < 53568", "colz"); + nt2->Draw("dxH:dzH>>h_dxH_vs_dzH_OTOF(300, -0.03, 0.03, 300, -0.03, 0.03)", "id >= 1920 && id < 55488", "colz"); addTLines(0.01); h = (TH2F*)gPad->GetPrimitive("h_dxH_vs_dzH_OTOF"); Info("OTOF", "RMS(dxH)=%.1f mu", h->GetRMS(2) * 1e4); @@ -283,5 +352,72 @@ void CheckDigitsIOTOF(std::string digifile = "tf3digits.root", std::string hitfi canvdXdZHit->SaveAs("trkdigits_dxH_vs_dzH.pdf"); f->Write(); + + std::string trackName[2] = {"Prm", "Sec"}; + f->mkdir("PrmTrkLayer0"); + f->mkdir("SecTrkLayer0"); + f->mkdir("PrmTrkLayer1"); + f->mkdir("SecTrkLayer1"); + for (int layer = 0; layer < 2; ++layer) { + for (int type = 0; type < 2; ++type) { + f->cd(Form("%sTrkLayer%d", trackName[type].c_str(), layer)); + hGenHitsEta[layer][type]->Write(); + hRecoDigitEta[layer][type]->Write(); + TH1F* hEffDigitEta = static_cast(hRecoDigitEta[layer][type]->Clone("hEffDigitEta")); + hEffDigitEta->Divide(hGenHitsEta[layer][type]); + // Set errors + for (int bin = 1; bin <= hEffDigitEta->GetNbinsX(); ++bin) { + double eff = hEffDigitEta->GetBinContent(bin); + double nGen = hGenHitsEta[layer][type]->GetBinContent(bin); + double err = 0.0; + if (nGen > 0) { + err = std::sqrt(eff * (1 - eff) / nGen); + } + hEffDigitEta->SetBinError(bin, err); + } + hEffDigitEta->SetTitle(";#eta;Digit Efficiency"); + hEffDigitEta->Write(); + delete hEffDigitEta; + } + } + + // Plot avg fraction of charge collected by digits for + // each hit vs eta, should reflect the digit efficiency + for (int layer = 0; layer < 2; ++layer) { + for (int type = 0; type < 2; ++type) { + + f->cd(Form("%sTrkLayer%d", trackName[type].c_str(), layer)); + TH2F* hFracCharge = new TH2F(Form("hFracCharge_Layer%d_Type%d", layer, type), ";Fraction of charge collected by digits;Entries", 40, -2, 2, 200, 0, 1); + + for (const auto& hitDigitPair : hitDigitMap) { + + uint64_t hitKey = hitDigitPair.first; + int evtID = static_cast(hitKey >> 32); + int hitIndex = static_cast(hitKey & 0xFFFFFFFF); + const auto& hit = (*hitArray[evtID])[hitIndex]; + + int hitLayer = gman->getIOTOFLayer(hit.GetDetectorID()); + if (hitLayer != layer) continue; + + float energyLoss = hit.GetEnergyLoss(); // in GeV + int charge = static_cast(energyLoss * 2.77778e+08); + + auto& mcTrack = mcTracksPerEvent[evtID]->at(hit.GetTrackID()); + bool isPrimary = mcTrack.isPrimary(); + if ((isPrimary ? 0 : 1) != type) continue; + + const auto& digitIndices = hitDigitPair.second; + float totalDigitCharge = 0.0f; + for (int digitIndex : digitIndices) { + totalDigitCharge += (*digArr)[digitIndex].getCharge(); + } + float fracCharge = totalDigitCharge / charge; + hFracCharge->Fill(mcTrack.GetEta(), fracCharge); + } + hFracCharge->Write(); + delete hFracCharge; + } + } + f->Close(); } diff --git a/Detectors/Upgrades/ALICE3/IOTOF/reconstruction/include/IOTOFReconstruction/ClustererParam.h b/Detectors/Upgrades/ALICE3/IOTOF/reconstruction/include/IOTOFReconstruction/ClustererParam.h index af007ada4c530..038cf639ba674 100644 --- a/Detectors/Upgrades/ALICE3/IOTOF/reconstruction/include/IOTOFReconstruction/ClustererParam.h +++ b/Detectors/Upgrades/ALICE3/IOTOF/reconstruction/include/IOTOFReconstruction/ClustererParam.h @@ -35,12 +35,6 @@ struct ClustererParam : public o2::conf::ConfigurableParamHelper // boilerplate stuff + make principal key O2ParamDef(ClustererParam, "TF3ClustererParam"); - - private: - static constexpr float DEFNoisePerPixel() - { - return 1e-8; // ITS/MFT values here!! - } }; } // namespace iotof diff --git a/Detectors/Upgrades/ALICE3/IOTOF/reconstruction/include/IOTOFReconstruction/TopologyClassifier.h b/Detectors/Upgrades/ALICE3/IOTOF/reconstruction/include/IOTOFReconstruction/TopologyClassifier.h index d837dae3948d2..5995515ea06a8 100644 --- a/Detectors/Upgrades/ALICE3/IOTOF/reconstruction/include/IOTOFReconstruction/TopologyClassifier.h +++ b/Detectors/Upgrades/ALICE3/IOTOF/reconstruction/include/IOTOFReconstruction/TopologyClassifier.h @@ -40,6 +40,7 @@ enum Topologies : uint8_t { kLineOnRow, kLineOnCol, kSquare, + kRectangle, kDiagonal, kLowerTriangleLeft, kLowerTriangleRight, @@ -67,6 +68,17 @@ struct TopologyInfo { int mFrequency = 0; Topologies mTopology = Topologies::kNTopologies; uint16_t mPattern; ///< Bitmask of fired pixels + + void print() const { + LOG(info) << "---> TopologyInfo: Topology = " << static_cast(mTopology) + << ", SizeX = " << mSizeX << ", SizeZ = " << mSizeZ + << ", OffsetXToCOG = " << mOffsetXToCOG << ", OffsetZToCOG = " << mOffsetZToCOG + << ", XMean = " << mXMean << ", ZMean = " << mZMean + << ", XSigma2 = " << mXSigma2 << ", ZSigma2 = " << mZSigma2 + << ", NPixels = " << mNPixels + << ", Frequency = " << mFrequency + << ", Pattern (bitmask) = 0x" << std::hex << mPattern; + } }; class TopologyClassifier { diff --git a/Detectors/Upgrades/ALICE3/IOTOF/reconstruction/src/Clusterer.cxx b/Detectors/Upgrades/ALICE3/IOTOF/reconstruction/src/Clusterer.cxx index 7f1c93672bcac..973a8911f68f7 100644 --- a/Detectors/Upgrades/ALICE3/IOTOF/reconstruction/src/Clusterer.cxx +++ b/Detectors/Upgrades/ALICE3/IOTOF/reconstruction/src/Clusterer.cxx @@ -117,18 +117,36 @@ void Clusterer::ClustererThread::processChip(gsl::span digits, // 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 = mParent->mSortIdx; - LOG(info) << ""; - LOG(info) << "----------------- NEW CHIP -----------------"; - + // LOG(info) << ""; + // LOG(info) << "----------------- NEW CHIP -----------------"; + // for (int i = 0; i < nDigits; ++i) { + // const auto& digit = digits[sortIdx[firstDigitIdx + i]]; + // LOG(info) << "[Clusterer] Digit " << i << "/" << nDigits << ": chipID=" << digit.getChipIndex() + // << ", row=" << digit.getRow() << ", col=" << digit.getColumn() + // << ", charge=" << digit.getCharge() << ", time=" << digit.getTime(); + // } 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); + + for (int i = 0; i < nDigits; ++i) { + digitIdxs[i] = sortIdx[firstDigitIdx + i]; + } + + findClustersMultipleHits( + digits, + gsl::span(digitIdxs), + 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); + // } // Flush per-thread output into the caller's containers if (!mClusters.empty()) { @@ -258,7 +276,7 @@ void Clusterer::ClustererThread::findClustersMultipleHits(gsl::span 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); + Cluster cluster(minRow, minCol, rowSpan, colSpan, firedDigitsMask, clsTopology, chipID, time); LOG(info) << "Pushing back cluster with row: " << row << ", col: " << col << ", rowSpan: " << rowSpan << ", colSpan: " << colSpan << ", pattern: " << firedDigitsMask diff --git a/Detectors/Upgrades/ALICE3/IOTOF/reconstruction/src/IOTOFReconstructionLinkDef.h b/Detectors/Upgrades/ALICE3/IOTOF/reconstruction/src/IOTOFReconstructionLinkDef.h index 8d6b3e3e1fa14..46b6d93506d59 100644 --- a/Detectors/Upgrades/ALICE3/IOTOF/reconstruction/src/IOTOFReconstructionLinkDef.h +++ b/Detectors/Upgrades/ALICE3/IOTOF/reconstruction/src/IOTOFReconstructionLinkDef.h @@ -17,6 +17,8 @@ #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+; diff --git a/Detectors/Upgrades/ALICE3/IOTOF/reconstruction/src/TopologyClassifier.cxx b/Detectors/Upgrades/ALICE3/IOTOF/reconstruction/src/TopologyClassifier.cxx index c0cf00c8be9c4..2daf6357b63b9 100644 --- a/Detectors/Upgrades/ALICE3/IOTOF/reconstruction/src/TopologyClassifier.cxx +++ b/Detectors/Upgrades/ALICE3/IOTOF/reconstruction/src/TopologyClassifier.cxx @@ -107,29 +107,36 @@ void TopologyClassifier::accountTopology(uint16_t bitmask, uint16_t minRow, uint 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 hasTopLeft = hasDigit(minRow, minCol); - const bool hasTopRight = hasDigit(minRow, maxCol); - const bool hasBottomLeft = hasDigit(maxRow, minCol); - const bool hasBottomRight = hasDigit(maxRow, maxCol); + 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 square + // Diagonal and triangles if (spanRow == spanCol) { - if ((hasTopLeft && hasBottomRight && !hasTopRight && !hasBottomLeft) || - (!hasTopLeft && !hasBottomRight && hasTopRight && hasBottomLeft)) { - newTopo.mTopology = Topologies::kDiagonal; - mTopologyCache[packKey(spanRow, spanCol, bitmask)] = newTopo; - return; - } - - if (hasTopLeft && hasTopRight && hasBottomLeft && hasBottomRight) { - newTopo.mTopology = Topologies::kSquare; - mTopologyCache[packKey(spanRow, spanCol, bitmask)] = newTopo; - return; - } - - // Triangles (exactly one missing corner) + // Triangles const int nCorners = hasTopLeft + hasTopRight + hasBottomLeft + hasBottomRight; if (nCorners == 3) { const int missing = !hasTopLeft ? 0 : !hasTopRight ? 1 : !hasBottomLeft ? 2 : 3; @@ -143,21 +150,28 @@ void TopologyClassifier::accountTopology(uint16_t bitmask, uint16_t minRow, uint 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 + 1, minCol); - const bool hasMiddleMax = hasDigit(minRow + 1, maxCol); + const bool hasMiddleMin = hasDigit(minRow, minCol + 1); + const bool hasMiddleMax = hasDigit(minRow, maxCol + 1); if (hasMiddleMin && hasMiddleMax) { - if (hasTopLeft && hasBottomRight) { + if (!hasTopLeft && !hasBottomRight && hasTopRight && hasBottomLeft) { newTopo.mTopology = Topologies::kSnake; mTopologyCache[packKey(spanRow, spanCol, bitmask)] = newTopo; return; } - if (!hasTopLeft && !hasBottomRight) { + if (hasTopLeft && hasBottomRight && !hasTopRight && !hasBottomLeft) { newTopo.mTopology = Topologies::kSnakeRefl; mTopologyCache[packKey(spanRow, spanCol, bitmask)] = newTopo; return; @@ -167,17 +181,17 @@ void TopologyClassifier::accountTopology(uint16_t bitmask, uint16_t minRow, uint // Snake rotated by 90 degrees: 2 x 3 if (spanRow == 2 && spanCol == 3) { - const bool hasMiddleLeft = hasDigit(minRow, minCol + 1); - const bool hasMiddleRight = hasDigit(maxRow, minCol + 1); + const bool hasMiddleLeft = hasDigit(minRow + 1, minCol); + const bool hasMiddleRight = hasDigit(maxRow + 1, minCol); if (hasMiddleLeft && hasMiddleRight) { - if (hasTopLeft && hasBottomRight) { + if (!hasTopLeft && !hasBottomRight && hasTopRight && hasBottomLeft) { newTopo.mTopology = Topologies::kSnakeRot90; mTopologyCache[packKey(spanRow, spanCol, bitmask)] = newTopo; return; } - if (!hasTopLeft && !hasBottomRight) { + if (hasTopLeft && hasBottomRight && !hasTopRight && !hasBottomLeft) { newTopo.mTopology = Topologies::kSnakeRot90Refl; mTopologyCache[packKey(spanRow, spanCol, bitmask)] = newTopo; return; @@ -190,7 +204,6 @@ void TopologyClassifier::accountTopology(uint16_t bitmask, uint16_t minRow, uint mTopologyCache[packKey(spanRow, spanCol, bitmask)] = newTopo; return; } - // Insert in map } @@ -264,12 +277,8 @@ void TopologyClassifier::print() { 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; + << ", Bitmask: " << std::bitset<16>(bitmask); + topoInfo.print(); } } diff --git a/Detectors/Upgrades/ALICE3/IOTOF/simulation/src/Digitizer.cxx b/Detectors/Upgrades/ALICE3/IOTOF/simulation/src/Digitizer.cxx index 5de01639312c4..2cac7e471df19 100644 --- a/Detectors/Upgrades/ALICE3/IOTOF/simulation/src/Digitizer.cxx +++ b/Detectors/Upgrades/ALICE3/IOTOF/simulation/src/Digitizer.cxx @@ -194,7 +194,7 @@ void Digitizer::stepping(const o2::itsmft::Hit& hit, float**& respMatrix, int& r LOG(debug) << "Hit position out of bounds for detector ID " << chipID; return; // hit is outside the active area } - xyzPositionEnd += stepVector; + xyzPositionEnd -= stepVector; } if (rowStart > rowEnd) { 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();