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..b16a6e8bf2f39 100644 --- a/Detectors/Upgrades/ALICE3/IOTOF/DataFormatsIOTOF/include/DataFormatsIOTOF/Cluster.h +++ b/Detectors/Upgrades/ALICE3/IOTOF/DataFormatsIOTOF/include/DataFormatsIOTOF/Cluster.h @@ -9,28 +9,180 @@ // granted to it by virtue of its status as an Intergovernmental Organization // or submit itself to any jurisdiction. +/// \file Cluster.h +/// \brief Definition of the IOTOF cluster #ifndef ALICEO2_DATAFORMATSIOTOF_CLUSTER_H #define ALICEO2_DATAFORMATSIOTOF_CLUSTER_H -#include #include #include +#include +#include + +#include "Framework/Logger.h" + +namespace o2 +{ +namespace iotof +{ + +/// Compact encoding for ALICE3 IOTOF cluster parameters inside a single 64-bit word. +struct ClusterInfo { + // Bit widths (Total: 52 bits out of 64) + static constexpr int NBitsRow = 9; + static constexpr int NBitsCol = 8; + static constexpr int NBitsRowSpan = 4; + static constexpr int NBitsColSpan = 4; + static constexpr int NBitsPattern = 16; + static constexpr int NBitsTopology = 11; + + // Bit offsets (ordered logically from LSB to MSB) + static constexpr int ShiftRow = 0; + static constexpr int ShiftCol = ShiftRow + NBitsRow; // 9 + static constexpr int ShiftRowSpan = ShiftCol + NBitsCol; // 17 + static constexpr int ShiftColSpan = ShiftRowSpan + NBitsRowSpan; // 21 + static constexpr int ShiftPattern = ShiftColSpan + NBitsColSpan; // 25 + static constexpr int ShiftTopology = ShiftPattern + NBitsPattern; // 41 + + // Bit masks + static constexpr uint64_t MaskRow = (1ULL << NBitsRow) - 1; + static constexpr uint64_t MaskCol = (1ULL << NBitsCol) - 1; + static constexpr uint64_t MaskRowSpan = (1ULL << NBitsRowSpan) - 1; + static constexpr uint64_t MaskColSpan = (1ULL << NBitsColSpan) - 1; + static constexpr uint64_t MaskPattern = (1ULL << NBitsPattern) - 1; + static constexpr uint64_t MaskTopology = (1ULL << NBitsTopology) - 1; + + uint64_t data{0}; + + // Constructors + constexpr ClusterInfo() = default; + constexpr ClusterInfo(uint64_t d) : data(d) {} + + // Static packer + static constexpr uint64_t pack(uint32_t row, uint32_t col, uint8_t rowSpan, + uint8_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 uint8_t getRowSpan() const { return (data >> ShiftRowSpan) & MaskRowSpan; } + constexpr uint8_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(uint8_t rs) + { + data = (data & ~(MaskRowSpan << ShiftRowSpan)) | ((static_cast(rs) & MaskRowSpan) << ShiftRowSpan); + } + constexpr void setColSpan(uint8_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(); } + uint8_t getRowSpan() const { return mClusterInfo.getRowSpan(); } + uint8_t getColSpan() const { return mClusterInfo.getColSpan(); } + uint32_t getPattern() const { return mClusterInfo.getPattern(); } + uint32_t getTopology() const { return mClusterInfo.getTopology(); } + int getSize() const + { + // Count the number of set bits in the pattern to determine the size of the cluster + uint32_t pattern = getPattern(); + int size = 0; + while (pattern) { + size += pattern & 1; + pattern >>= 1; + } + return size; + } + // BaseCluster / Interface Compatibility Getters + uint32_t getChipID() const { return mChipID; } + uint32_t getSensorID() const { return mChipID; } + time_t getTime() const { return mTime; } + uint64_t getPackedData() const { return mClusterInfo.data; } + + // Setters + void setRow(UShort_t r) { mClusterInfo.setRow(r); } + void setCol(UShort_t c) { mClusterInfo.setCol(c); } + void setRowSpan(UShort_t rs) { mClusterInfo.setRowSpan(rs); } + void setColSpan(UShort_t cs) { mClusterInfo.setColSpan(cs); } + void setPatternID(UShort_t p) { mClusterInfo.setPattern(p); } + void setTopology(UShort_t t) { mClusterInfo.setTopology(t); } + void setChipID(UShort_t c) { mChipID = c; } + void setTime(time_t t) { mTime = t; } + + // Operators & Debugging + bool operator==(const Cluster& cl) const + { + return mClusterInfo.data == cl.mClusterInfo.data && mChipID == cl.mChipID && mTime == cl.mTime; + } + + void print() const; std::string asString() const; - ClassDefNV(Cluster, 1); + private: + ClusterInfo mClusterInfo{}; ///< 64-bit packed structure containing geometry/topology + UShort_t mChipID{0}; ///< Chip / Sensor ID + float mTime{0.0f}; ///< Hit timing information + + void sanityCheck(); + + ClassDefNV(Cluster, 2); }; -} // namespace o2::iotof +} // namespace iotof +} // namespace o2 + +std::ostream& operator<<(std::ostream& stream, const o2::iotof::Cluster& cl); -#endif +#endif /* ALICEO2_DATAFORMATSIOTOF_CLUSTER_H */ diff --git a/Detectors/Upgrades/ALICE3/IOTOF/DataFormatsIOTOF/src/Cluster.cxx b/Detectors/Upgrades/ALICE3/IOTOF/DataFormatsIOTOF/src/Cluster.cxx index 6b5a4948900e7..b0ea13477909f 100644 --- a/Detectors/Upgrades/ALICE3/IOTOF/DataFormatsIOTOF/src/Cluster.cxx +++ b/Detectors/Upgrades/ALICE3/IOTOF/DataFormatsIOTOF/src/Cluster.cxx @@ -9,19 +9,64 @@ // granted to it by virtue of its status as an Intergovernmental Organization // or submit itself to any jurisdiction. +/// \file Cluster.cxx +/// \brief Implementation of the IOTOF cluster + #include "DataFormatsIOTOF/Cluster.h" -#include +#include "Framework/Logger.h" +#include +#include +#include +// Root ClassImp macros for serialization metadata +ClassImp(o2::iotof::ClusterInfo); ClassImp(o2::iotof::Cluster); -namespace o2::iotof +namespace o2 +{ +namespace iotof { std::string Cluster::asString() const { - std::ostringstream stream; - stream << "chip=" << chipID << " row=" << row << " col=" << col << " size=" << size; - return stream.str(); + LOG(debug) << "[Cluster::asString] Converting Cluster to string"; + return std::format( + "chip: {:5d} | row: {:3d} col: {:3d} | span: {:2d}x{:2d} | pattern: {:5d} topology: {:4d}", + getChipID(), + getRow(), + getCol(), + getRowSpan(), + getColSpan(), + getPattern(), + getTopology()); +} + +//______________________________________________________________________________ +void Cluster::print() const +{ + std::cout << *this << "\n"; +} + +//______________________________________________________________________________ +void Cluster::sanityCheck() +{ + LOG(debug) << "[Cluster::sanityCheck] Performing sanity check on Cluster fields"; + + // Ensure extracted values fit within allowed bit masks + assert(getRow() <= ClusterInfo::MaskRow); + assert(getCol() <= ClusterInfo::MaskCol); + assert(getRowSpan() <= ClusterInfo::MaskRowSpan); + assert(getColSpan() <= ClusterInfo::MaskColSpan); + assert(getPattern() <= ClusterInfo::MaskPattern); + assert(getTopology() <= ClusterInfo::MaskTopology); } -} // namespace o2::iotof +} // namespace iotof +} // namespace o2 + +// Stream operator implementation +std::ostream& operator<<(std::ostream& stream, const o2::iotof::Cluster& cl) +{ + stream << cl.asString(); + return stream; +} diff --git a/Detectors/Upgrades/ALICE3/IOTOF/DataFormatsIOTOF/src/DataFormatsIOTOFLinkDef.h b/Detectors/Upgrades/ALICE3/IOTOF/DataFormatsIOTOF/src/DataFormatsIOTOFLinkDef.h index 7e121273d3fab..e639584ebfa75 100644 --- a/Detectors/Upgrades/ALICE3/IOTOF/DataFormatsIOTOF/src/DataFormatsIOTOFLinkDef.h +++ b/Detectors/Upgrades/ALICE3/IOTOF/DataFormatsIOTOF/src/DataFormatsIOTOFLinkDef.h @@ -18,6 +18,7 @@ #pragma link C++ class o2::iotof::Digit + ; #pragma link C++ class std::vector < o2::iotof::Digit> + ; +#pragma link C++ class o2::iotof::ClusterInfo + ; #pragma link C++ class o2::iotof::Cluster + ; #pragma link C++ class std::vector < o2::iotof::Cluster> + ; diff --git a/Detectors/Upgrades/ALICE3/IOTOF/base/CMakeLists.txt b/Detectors/Upgrades/ALICE3/IOTOF/base/CMakeLists.txt index 3b47b9451916d..c5c2b1c36bcab 100644 --- a/Detectors/Upgrades/ALICE3/IOTOF/base/CMakeLists.txt +++ b/Detectors/Upgrades/ALICE3/IOTOF/base/CMakeLists.txt @@ -11,10 +11,12 @@ o2_add_library(IOTOFBase SOURCES src/GeometryTGeo.cxx + src/Segmentation.cxx src/IOTOFBaseParam.cxx PUBLIC_LINK_LIBRARIES O2::DetectorsBase O2::MathUtils) o2_target_root_dictionary(IOTOFBase HEADERS include/IOTOFBase/GeometryTGeo.h + include/IOTOFBase/Segmentation.h include/IOTOFBase/IOTOFBaseParam.h) diff --git a/Detectors/Upgrades/ALICE3/IOTOF/simulation/include/IOTOFSimulation/Segmentation.h b/Detectors/Upgrades/ALICE3/IOTOF/base/include/IOTOFBase/Segmentation.h similarity index 93% rename from Detectors/Upgrades/ALICE3/IOTOF/simulation/include/IOTOFSimulation/Segmentation.h rename to Detectors/Upgrades/ALICE3/IOTOF/base/include/IOTOFBase/Segmentation.h index ddde28cf7dd7a..c726998fcd4bc 100644 --- a/Detectors/Upgrades/ALICE3/IOTOF/simulation/include/IOTOFSimulation/Segmentation.h +++ b/Detectors/Upgrades/ALICE3/IOTOF/base/include/IOTOFBase/Segmentation.h @@ -50,11 +50,11 @@ class Segmentation /// the center of the sensitive volulme. /// \param int iRow Detector x cell coordinate. Has the range 0 <= iRow < mNumberOfRows /// \param int iCol Detector z cell coordinate. Has the range 0 <= iCol < mNumberOfColumns - bool localToDetector(float x, float z, int& iRow, int& iCol, const int subDetectorID); + bool localToDetector(float x, float z, int& iRow, int& iCol, const int subDetectorID) const; /// same but w/o check for row/column range - void localToDetectorUnchecked(float xRow, float zCol, int& iRow, int& iCol, const int subDetectorID); + void localToDetectorUnchecked(float xRow, float zCol, int& iRow, int& iCol, const int subDetectorID) const; - /// Transformation from Detector cell coordiantes to Geant detector centered + /// Transformation from Detector cell coordinates to Geant detector centered /// local coordinates (cm) /// \param int iRow Detector x cell coordinate. Has the range 0 <= iRow < mNumberOfRows /// \param int iCol Detector z cell coordinate. Has the range 0 <= iCol < mNumberOfColumns @@ -67,7 +67,7 @@ class Segmentation // w/o check for row/col range template - void detectorToLocalUnchecked(L row, L col, T& xRow, T& zCol, const int subDetectorID) + void detectorToLocalUnchecked(L row, L col, T& xRow, T& zCol, const int subDetectorID) const { if (subDetectorID != 0 && subDetectorID != 1) { row = col = -1; @@ -78,7 +78,7 @@ class Segmentation zCol = col * specsConfig.PitchCol + getFirstColCoordinate(subDetectorID); } template - void detectorToLocalUnchecked(L row, L col, math_utils::Point3D& loc, const int subDetectorID) + void detectorToLocalUnchecked(L row, L col, math_utils::Point3D& loc, const int subDetectorID) const { if (subDetectorID != 0 && subDetectorID != 1) { row = col = -1; @@ -88,7 +88,7 @@ class Segmentation loc.SetCoordinates(getFirstRowCoordinate(subDetectorID) - row * specsConfig.PitchRow, T(0.), col * specsConfig.PitchCol + getFirstColCoordinate(subDetectorID)); } template - void detectorToLocalUnchecked(L row, L col, std::array& loc, const int subDetectorID) + void detectorToLocalUnchecked(L row, L col, std::array& loc, const int subDetectorID) const { if (subDetectorID != 0 && subDetectorID != 1) { row = col = -1; @@ -103,7 +103,7 @@ class Segmentation // same but with check for row/col range template - bool detectorToLocal(L row, L col, T& xRow, T& zCol, const int subDetectorID) + bool detectorToLocal(L row, L col, T& xRow, T& zCol, const int subDetectorID) const { if (subDetectorID != 0 && subDetectorID != 1) { row = col = -1; @@ -118,7 +118,7 @@ class Segmentation } template - bool detectorToLocal(L row, L col, math_utils::Point3D& loc, const int subDetectorID) + bool detectorToLocal(L row, L col, math_utils::Point3D& loc, const int subDetectorID) const { if (subDetectorID != 0 && subDetectorID != 1) { row = col = -1; @@ -132,7 +132,7 @@ class Segmentation return true; } template - bool detectorToLocal(L row, L col, std::array& loc, const int subDetectorID) + bool detectorToLocal(L row, L col, std::array& loc, const int subDetectorID) const { if (subDetectorID != 0 && subDetectorID != 1) { row = col = -1; @@ -146,12 +146,12 @@ class Segmentation return true; } - float getFirstRowCoordinate(const int subDetectorID) + float getFirstRowCoordinate(const int subDetectorID) const { const auto& specsConfig = ChipSpecificsParam::Instance(); return 0.5 * ((specsConfig.ActiveMatrixSizeRows() - specsConfig.PassiveEdgeTop + specsConfig.PassiveEdgeReadOut) - specsConfig.PitchRow); } - float getFirstColCoordinate(const int subDetectorID) + float getFirstColCoordinate(const int subDetectorID) const { const auto& specsConfig = ChipSpecificsParam::Instance(); return 0.5 * (specsConfig.PitchCol - specsConfig.ActiveMatrixSizeCols()); @@ -161,7 +161,7 @@ class Segmentation }; //_________________________________________________________________________________________________ -inline void Segmentation::localToDetectorUnchecked(float xRow, float zCol, int& iRow, int& iCol, const int subDetectorID) +inline void Segmentation::localToDetectorUnchecked(float xRow, float zCol, int& iRow, int& iCol, const int subDetectorID) const { // convert to row/col w/o over/underflow check if (subDetectorID != 0 && subDetectorID != 1) { @@ -187,7 +187,7 @@ inline void Segmentation::localToDetectorUnchecked(float xRow, float zCol, int& } //_________________________________________________________________________________________________ -inline bool Segmentation::localToDetector(float xRow, float zCol, int& iRow, int& iCol, const int subDetectorID) +inline bool Segmentation::localToDetector(float xRow, float zCol, int& iRow, int& iCol, const int subDetectorID) const { // convert to row/col if (subDetectorID != 0 && subDetectorID != 1) { diff --git a/Detectors/Upgrades/ALICE3/IOTOF/base/src/GeometryTGeo.cxx b/Detectors/Upgrades/ALICE3/IOTOF/base/src/GeometryTGeo.cxx index 8c8a36877eca8..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/base/src/IOTOFBaseLinkDef.h b/Detectors/Upgrades/ALICE3/IOTOF/base/src/IOTOFBaseLinkDef.h index cb5b047e72077..ba9457a4b96c9 100644 --- a/Detectors/Upgrades/ALICE3/IOTOF/base/src/IOTOFBaseLinkDef.h +++ b/Detectors/Upgrades/ALICE3/IOTOF/base/src/IOTOFBaseLinkDef.h @@ -16,6 +16,7 @@ #pragma link off all functions; #pragma link C++ class o2::iotof::GeometryTGeo + ; +#pragma link C++ class o2::iotof::Segmentation + ; #pragma link C++ class o2::iotof::IOTOFBaseParam + ; #pragma link C++ class o2::conf::ConfigurableParamHelper < o2::iotof::IOTOFBaseParam> + ; diff --git a/Detectors/Upgrades/ALICE3/IOTOF/simulation/src/Segmentation.cxx b/Detectors/Upgrades/ALICE3/IOTOF/base/src/Segmentation.cxx similarity index 96% rename from Detectors/Upgrades/ALICE3/IOTOF/simulation/src/Segmentation.cxx rename to Detectors/Upgrades/ALICE3/IOTOF/base/src/Segmentation.cxx index a7ec0d708c3b8..aa77bf50d069d 100644 --- a/Detectors/Upgrades/ALICE3/IOTOF/simulation/src/Segmentation.cxx +++ b/Detectors/Upgrades/ALICE3/IOTOF/base/src/Segmentation.cxx @@ -12,7 +12,7 @@ /// \file Segmentation.cxx /// \brief Implementation of the Segmentation class -#include "IOTOFSimulation/Segmentation.h" +#include "IOTOFBase/Segmentation.h" #include "IOTOFBase/IOTOFBaseParam.h" #include diff --git a/Detectors/Upgrades/ALICE3/IOTOF/macros/CheckClustersIOTOF.C b/Detectors/Upgrades/ALICE3/IOTOF/macros/CheckClustersIOTOF.C index 107e5a4d02bf8..01a069b59232e 100644 --- a/Detectors/Upgrades/ALICE3/IOTOF/macros/CheckClustersIOTOF.C +++ b/Detectors/Upgrades/ALICE3/IOTOF/macros/CheckClustersIOTOF.C @@ -9,239 +9,297 @@ // granted to it by virtue of its status as an Intergovernmental Organization // or submit itself to any jurisdiction. -/// \file CheckClustersIOTOF.C -/// \brief Simple macro to create clusters from TF3 digits +/// \file CheckClusters.C +/// \brief Simple macro to check TF3 clusters -#if !defined(__CLING__) || defined(__ROOTCLING__) #include #include -#include +#include #include +#include #include +#include #include -#include "IOTOFSimulation/Segmentation.h" -#include "IOTOFBase/IOTOFBaseParam.h" +#include "IOTOFBase/Segmentation.h" #include "IOTOFBase/GeometryTGeo.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 "IOTOFReconstruction/TopologyClassifier.h" +#include "ITSMFTSimulation/Hit.h" #include "DetectorsBase/GeometryManager.h" +#if !defined(__CLING__) || defined(__ROOTCLING__) +#include +#include +#include +#include +#include +#include +#include "IOTOFBase/IOTOFBaseParam.h" +#include "IOTOFBase/GeometryTGeo.h" +#include "DataFormatsIOTOF/Cluster.h" +#include "IOTOFReconstruction/TopologyClassifier.h" +#include "ITSMFTSimulation/Hit.h" #include "DataFormatsITSMFT/ROFRecord.h" - +#include "MathUtils/Cartesian.h" +#include "MathUtils/Utils.h" +#include "SimulationDataFormat/MCCompLabel.h" +#include "SimulationDataFormat/MCTruthContainer.h" +#include "DetectorsCommonDataFormats/DetectorNameConf.h" +#include "CCDB/BasicCCDBManager.h" #endif #define ENABLE_UPGRADES -void CheckClustersIOTOF(std::string digiFilePath = "tf3digits.root", std::string clsFilePath = "tf3clusters.root", std::string inputGeomPath = "o2sim_geometry.root") +void addTLines(float pitch) { - gStyle->SetPalette(55); + // Add grid lines at multiples of pitch on the current pad + if (!gPad) + return; + + gPad->Update(); + + Double_t xmin = gPad->GetUxmin(); + Double_t xmax = gPad->GetUxmax(); + Double_t ymin = gPad->GetUymin(); + Double_t ymax = gPad->GetUymax(); + + // Calculate the first vertical line position (multiple of pitch) + int nLinesX = 0; + for (float x = xmin; x <= xmax && nLinesX < 1000; x += pitch, nLinesX++) { + TLine* line = new TLine(x, ymin, x, ymax); + line->SetLineStyle(2); + line->SetLineColor(kGray); + line->Draw("same"); + } + + // Calculate the first horizontal line position (multiple of pitch) + int nLinesY = 0; + for (float y = ymin; y <= ymax && nLinesY < 1000; y += pitch, nLinesY++) { + TLine* line = new TLine(xmin, y, xmax, y); + line->SetLineStyle(2); + line->SetLineColor(kGray); + line->Draw("same"); + } + + gPad->Modified(); + gPad->Update(); +} + +void CheckClustersIOTOF(std::string clusfile = "tf3clusters.root", + std::string hitfile = "o2sim_HitsTF3.root", + std::string topodictfile = "TF3ClusterTopologies.root", + std::string inputGeom = "", + std::string cfgStr = "IOTOFBase.segmentedInnerTOF=true;IOTOFBase.segmentedOuterTOF=true;IOTOFBase.enableForwardTOF=false;IOTOFBase.enableBackwardTOF=false;") +{ + std::cout << "CheckClustersIOTOF: clusfile=" << clusfile << ", hitfile=" << hitfile << ", inputGeom=" << inputGeom << std::endl; + const int QEDSourceID = 99; // Clusters from this MC source correspond to QED electrons using namespace o2::base; using namespace o2::iotof; using o2::iotof::Cluster; - using o2::iotof::Digit; + using o2::itsmft::Hit; + + o2::conf::ConfigurableParam::updateFromString(cfgStr); + const auto& chipInfo = o2::iotof::ChipSpecificsParam::Instance(); + auto seg = o2::iotof::Segmentation::Instance(); - o2::conf::ConfigurableParam::updateFromString("IOTOFBase.segmentedInnerTOF=true;IOTOFBase.segmentedOuterTOF=true;IOTOFBase.enableForwardTOF=false;IOTOFBase.enableBackwardTOF=false"); + using ROFRec = o2::itsmft::ROFRecord; + using MC2ROF = o2::itsmft::MC2ROFRecord; + using HitVec = std::vector; + using MC2HITS_map = std::unordered_map; // maps (track_ID<<16 + chip_ID) to entry in the hit vector - auto segGeom = o2::iotof::Segmentation::Instance(); + std::vector hitVecPool; + std::vector mc2hitVec; + + TFile fout("CheckClusters.root", "recreate"); + TNtuple nt("ntc", "cluster ntuple", "chip:ev:lab:hlx:hlz:cgx:cgy:cgz:dx:dz"); // 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); + o2::base::GeometryManager::loadGeometry(inputGeom); + auto* gman = o2::iotof::GeometryTGeo::Instance(); + gman->fillMatrixCache(o2::math_utils::bit2Mask(o2::math_utils::TransformType::L2G)); + + // Cluster topologies dictionary + TFile* clsTopoFile = TFile::Open(topodictfile.data(), "READ"); + auto* clsTopoMapPtr = clsTopoFile->Get>("TF3ClusterTopologies"); + if (clsTopoMapPtr) { + std::cout << "Loaded " << clsTopoMapPtr->size() << " entries from " << topodictfile << std::endl; + } else { + std::cerr << "Failed to load TF3ClusterTopologies from " << topodictfile << std::endl; + } + // Construct map directly from the vector pairs + std::unordered_map topoMap(clsTopoMapPtr->begin(), clsTopoMapPtr->end()); + TopologyClassifier topoClassifier(std::move(topoMap)); + topoClassifier.setGeometry(gman); + topoClassifier.print(); + clsTopoFile->Close(); + + // Hits + TFile fileH(hitfile.data()); + TTree* hitTree = (TTree*)fileH.Get("o2sim"); + std::vector* hitArray = nullptr; + hitTree->SetBranchAddress("TF3Hit", &hitArray); + mc2hitVec.resize(hitTree->GetEntries()); + hitVecPool.resize(hitTree->GetEntries(), nullptr); + int nEvts = hitTree->GetEntries(); + std::cout << "CheckClustersIOTOF: hitTree has " << hitTree->GetEntries() << " entries" << std::endl; // 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 - 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 << 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()) { - continue; - } - digitsLabels.emplace(label, iDigit); + TFile fileC(clusfile.data()); + TTree* clusTree = (TTree*)fileC.Get("o2sim"); + clusTree->ls(); + std::vector* clusArr = nullptr; + clusTree->SetBranchAddress("TF3Cluster", &clusArr); + std::vector* patternsPtr = nullptr; + auto pattBranch = clusTree->GetBranch("TF3ClusterPatt"); + if (pattBranch) { + pattBranch->SetAddress(&patternsPtr); } + std::cout << "CheckClustersIOTOF: clusTree has " << clusTree->GetEntries() << " entries" << std::endl; - // LOOP on : ROFRecord array - for (unsigned int iROF = 0; iROF < clsRofArr.size(); ++iROF) { - - const unsigned int rofIndex = clsRofArr[iROF].getFirstEntry(); - const unsigned int rofNEntries = clsRofArr[iROF].getNEntries(); - - // 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; - } + // ROFrecords + std::vector rofRecVec, *rofRecVecP = &rofRecVec; + clusTree->SetBranchAddress("TF3ClusterROF", &rofRecVecP); + std::cout << "CheckClustersIOTOF: rofRecVec has " << rofRecVec.size() << " entries" << std::endl; - 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); - - Float_t x{0.f}, y{0.f}, z{0.f}; - if (subDetID >= 0) { - segGeom->detectorToLocal(iRow, iCol, x, z, subDetID); - } - - o2::math_utils::Point3D localDigitCoord(x, y, z); // local Digit - - 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 + // Cluster MC labels + o2::dataformats::MCTruthContainer* clusLabArr = nullptr; + if (hitTree && clusTree->GetBranch("TF3ClusterMCTruth")) { + clusTree->SetBranchAddress("TF3ClusterMCTruth", &clusLabArr); + } - // 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; + clusTree->GetEntry(0); + std::cout << "Number of clusters: " << clusArr->size() << std::endl; + std::cout << "Number of pattern bytes: " << (patternsPtr ? patternsPtr->size() : 0) << std::endl; + std::cout << "Number of label indices: " << (clusLabArr ? clusLabArr->getIndexedSize() : 0) << std::endl; + // return; + int nROFRec = (int)rofRecVec.size(); + + // << build min and max MC events used by each ROF + auto pattIt = patternsPtr->cbegin(); + int invalidPattIDCounter{0}; + for (int irof = 0; irof < nROFRec; irof++) { + const auto& rofRec = rofRecVec[irof]; + rofRec.print(); + + + // >> read and map MC events contributing to this ROF + for (int im = 0; im <= nEvts; im++) { + if (!hitVecPool[im]) { + hitTree->SetBranchAddress("TF3Hit", &hitVecPool[im]); + hitTree->GetEntry(im); + auto& mc2hit = mc2hitVec[im]; + const auto* hitArray = hitVecPool[im]; + for (int ih = hitArray->size(); ih--;) { + const auto& hit = (*hitArray)[ih]; + uint64_t key = (uint64_t(hit.GetTrackID()) << 32) + hit.GetDetectorID(); + mc2hit.emplace(key, ih); + } } + } - 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; - - Float_t x = 0.f, y = 0.f, z = 0.f; - if (subDetID >= 0) { - segGeom->detectorToLocal(iRow, iCol, x, z, subDetID); + // << cache MC events contributing to this ROF + for (int icl = 0; icl < rofRec.getNEntries(); icl++) { + int clEntry = icl; // entry of icl-th cluster of this ROF in the vector of clusters + std::cout << "Processing cluster " << icl << "/" << rofRec.getNEntries() << std::endl; + const auto& cluster = (*clusArr)[clEntry]; + + float errX{0.f}; + float errZ{0.f}; + int npix = 0; + uint16_t pattID = cluster.getPattern(); + uint8_t spanRow = cluster.getRowSpan(); + uint8_t spanCol = cluster.getColSpan(); + o2::math_utils::Point3D locC; + // std::cout << "CIAO1" << std::endl; + if (pattID == o2::iotof::Cluster::InvalidPatternID) { + invalidPattIDCounter++; + continue; } - - 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()); - - // Match to digit - auto digitLabelFromCls = (clsLabels->getLabels(iCls))[0]; - auto digitEntry = digitsLabels.find(digitLabelFromCls); - - if (digitEntry == digitsLabels.end()) { - LOG(error) << "No matching digit for cluster " << iCls << " with label " << digitLabelFromCls.getRawValue(); + // std::cout << "CIAO2" << std::endl; + + uint32_t topoKey = TopologyClassifier::makeKey(spanRow, spanCol, pattID); + errX = topoClassifier.getErrX(topoKey); + errZ = topoClassifier.getErrZ(topoKey); + npix = topoClassifier.getNPixels(topoKey); + auto chipID = cluster.getSensorID(); + // std::cout << "CIAO3" << std::endl; + + // Transformation to the local --> global + locC = topoClassifier.getClusterCoordinates(cluster); + // std::cout << "CIAO31" << std::endl; + auto gloC = gman->getMatrixL2G(chipID) * locC; + // std::cout << "CIAO32" << std::endl; + + // Check how many labels are there + if (clusLabArr->getLabels(clEntry).empty()) { 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); + const auto& lab = (clusLabArr->getLabels(clEntry))[0]; + // std::cout << "CIAO33" << std::endl; + + // std::cout << "CIAO4" << std::endl; + if (!lab.isValid() || lab.getSourceID() == QEDSourceID) + continue; + // std::cout << "CIAO5" << std::endl; + + // get MC info + int trID = lab.getTrackID(); + const auto& mc2hit = mc2hitVec[lab.getEventID()]; + const auto* hitArray = hitVecPool[lab.getEventID()]; + uint64_t key = (uint64_t(trID) << 32) + chipID; + auto hitEntry = mc2hit.find(key); + if (hitEntry == mc2hit.end()) { + LOG(error) << "Failed to find MC hit entry for Tr" << trID << " chipID" << chipID; + continue; } - - 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"); - - 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(); - outFile->Close(); + // std::cout << "CIAO6" << std::endl; + const auto& hit = (*hitArray)[hitEntry->second]; + // + float dx = 0, dz = 0; + int ievH = lab.getEventID(); + o2::math_utils::Point3D locH, locHsta; + + // mean local position of the hit + locH = gman->getMatrixL2G(chipID) ^ (hit.GetPos()); // inverse conversion from global to local + locHsta = gman->getMatrixL2G(chipID) ^ (hit.GetPosStart()); + // std::cout << "CIAO7" << std::endl; + auto x0 = locHsta.X(), dltx = locH.X() - x0; + auto y0 = locHsta.Y(), dlty = locH.Y() - y0; + auto z0 = locHsta.Z(), dltz = locH.Z() - z0; + auto r = (0.5 * (chipInfo.SensorLayerThickness - chipInfo.SensorLayerThicknessEff) - y0) / dlty; + locH.SetXYZ(x0 + r * dltx, y0 + r * dlty, z0 + r * dltz); + // locH.SetXYZ(0.5 * (locH.X() + locHsta.X()), 0.5 * (locH.Y() + locHsta.Y()), 0.5 * (locH.Z() + locHsta.Z())); + std::array data = {(float)chipID, (float)lab.getEventID(), (float)trID, + locH.X(), locH.Z(), + gloC.X(), gloC.Y(), gloC.Z(), + locC.X() - locH.X(), locC.Z() - locH.Z()}; + // std::cout << "CIAO8" << std::endl; + nt.Fill(data.data()); + } + } + std::cout << "CheckClustersIOTOF: Found " << invalidPattIDCounter << " clusters with invalid pattern ID" << std::endl; + + // distributions of differences between local positions of digits and hits in x and z + 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)", "chip >= 0 && chip < 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)", "chip >= 1920 && chip < 55488", "colz"); + addTLines(0.01); + h = (TH2F*)gPad->GetPrimitive("h_dx_vs_dz_OTOF"); + Info("OTOF", "RMS(dx)=%.1f mu", h->GetRMS(2) * 1e4); + Info("OTOF", "RMS(dz)=%.1f mu", h->GetRMS(1) * 1e4); + canvdXdZ->SaveAs("tf3clusters_dx_vs_dz.pdf"); + canvdXdZ->SaveAs("tf3clusters_dx_vs_dz.root"); + + fout.cd(); + nt.Write(); } diff --git a/Detectors/Upgrades/ALICE3/IOTOF/macros/CheckDigitsIOTOF.C b/Detectors/Upgrades/ALICE3/IOTOF/macros/CheckDigitsIOTOF.C index 6caf2eb471b50..4a34bff8a0a73 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" diff --git a/Detectors/Upgrades/ALICE3/IOTOF/macros/CheckTopologiesIOTOF.C b/Detectors/Upgrades/ALICE3/IOTOF/macros/CheckTopologiesIOTOF.C new file mode 100644 index 0000000000000..71d5a706ed722 --- /dev/null +++ b/Detectors/Upgrades/ALICE3/IOTOF/macros/CheckTopologiesIOTOF.C @@ -0,0 +1,146 @@ +#include +#include +#include +#include +#include +#include + +#include "TFile.h" +#include "TH2F.h" +#include "TGraphErrors.h" + +#include "Framework/Logger.h" +#include "IOTOFBase/IOTOFBaseParam.h" +#include "IOTOFReconstruction/TopologyClassifier.h" + +using namespace o2::iotof; + +void CheckTopologiesIOTOF(const char* topoFileName = "TF3ClusterTopologies.root", + std::string chipCfgStr = "", + const char* outFileName = "CheckTopologies.root") +{ + + o2::conf::ConfigurableParam::updateFromString(chipCfgStr); + const auto& chipInfo = o2::iotof::ChipSpecificsParam::Instance(); + + // Cluster topologies dictionary + TFile* clsTopoFile = TFile::Open(topoFileName, "READ"); + auto* clsTopoMapPtr = clsTopoFile->Get>("TF3ClusterTopologies"); + if (clsTopoMapPtr) { + std::cout << "Loaded " << clsTopoMapPtr->size() << " entries from " << topoFileName << std::endl; + } else { + std::cerr << "Failed to load TF3ClusterTopologies from " << topoFileName << std::endl; + } + + // 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) << ""; + } + + // Topology names + const std::array topologyNames = { + "kSingleDigit", "kLineOnRow", "kLineOnCol", "kSquare", "kRectangle", "kDiagonal", + "kLowerTriangleLeft", "kLowerTriangleRight", "kUpperTriangleLeft", "kUpperTriangleRight", + "kSnake", "kSnakeRefl", "kSnakeRot90", "kSnakeRot90Refl", "kHuge", "kOther"}; + + // Create output ROOT file + auto* outFile = TFile::Open(outFileName, "RECREATE"); + TH1F* hTopoSummaryDictionary = new TH1F("hTopoSummaryDictionary", "Cluster Topology Count Summary;;Counts", kNTopologies, 0, kNTopologies); + for (const auto& topoName : topologyNames) { + hTopoSummaryDictionary->GetXaxis()->SetBinLabel(&topoName - &topologyNames[0] + 1, topoName.c_str()); + } + for (const auto& [topoKey, topology] : topoClassifier.getTopologyMap()) { + hTopoSummaryDictionary->Fill(topology.mTopology, topology.mFrequency); + hTopoSummaryDictionary->SetBinError(topology.mTopology + 1, 0); + } + hTopoSummaryDictionary->Write(); + // Create directory structures for all categories + for (const auto& topoName : topologyNames) { + outFile->mkdir(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; + int frequency = topology.mFrequency; + + float minRowCoord = -1.5 * chipInfo.PitchRow; + float maxRowCoord = chipInfo.PitchRow * (spanRow + 0.5); + float minColCoord = -1.5 * chipInfo.PitchCol; + float maxColCoord = chipInfo.PitchCol * (spanCol + 0.5); + TH2F* hTopoDisplay = new TH2F(Form("spanRow_%i_spanCol_%i_key_%i_all", spanRow, spanCol, topoKey), Form("Cluster Topology %s;Row;Column", topoName.c_str()), + spanRow + 2, minRowCoord, maxRowCoord, spanCol + 2, minColCoord, maxColCoord); + + // One-point TGraph for COG + TGraphErrors* gTopoCOG = new TGraphErrors(1); + gTopoCOG->SetName(Form("spanRow_%i_spanCol_%i_key_%i_COG", spanRow, spanCol, topoKey)); + gTopoCOG->SetTitle(Form("Cluster Topology %s COG", topoName.c_str())); + gTopoCOG->SetPoint(0, topology.mXMean, topology.mZMean); + gTopoCOG->SetPointError(0, std::sqrt(topology.mXSigma2), std::sqrt(topology.mZSigma2)); + gTopoCOG->SetMarkerStyle(20); + gTopoCOG->SetMarkerColor(kBlue); + + // 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)) { + hTopoDisplay->SetBinContent(row+2, col+2, frequency); + } + } + } + outFile->cd(topoName.c_str()); + hTopoDisplay->Write(); + gTopoCOG->Write(); + delete hTopoDisplay; + delete gTopoCOG; + } + + outFile->Close(); + LOG(info) << "Successfully wrote topology displays to " << outFileName; +} 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..931e08512b6a2 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 fetchMCLabels(uint32_t digID, const ConstDigitTruth* labelsDig, int& nFilled); + 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..388fec83143b9 --- /dev/null +++ b/Detectors/Upgrades/ALICE3/IOTOF/reconstruction/include/IOTOFReconstruction/ClustererParam.h @@ -0,0 +1,43 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +/// \file ClustererParam.h +/// \brief Definition of the IOTOF clusterer settings + +#ifndef ALICEO2_IOTOFCLUSTERERPARAM_H_ +#define ALICEO2_IOTOFCLUSTERERPARAM_H_ + +#include "DetectorsCommonDataFormats/DetID.h" +#include "CommonUtils/ConfigurableParam.h" +#include "CommonUtils/ConfigurableParamHelper.h" +#include +#include + +// TO BE REMOVED BEFORE PUSH +#include "Framework/Logger.h" + +namespace o2 +{ +namespace iotof +{ +struct ClustererParam : public o2::conf::ConfigurableParamHelper { + + int maxTimeDiffNSigma = 3; ///< maximum time difference in nsigma for clustering + int maxFiredDigitsForCls = 16; ///< maximum time difference in nsigma for clustering + + // boilerplate stuff + make principal key + O2ParamDef(ClustererParam, "TF3ClustererParam"); +}; + +} // namespace iotof +} // namespace o2 + +#endif diff --git a/Detectors/Upgrades/ALICE3/IOTOF/reconstruction/include/IOTOFReconstruction/TopologyClassifier.h b/Detectors/Upgrades/ALICE3/IOTOF/reconstruction/include/IOTOFReconstruction/TopologyClassifier.h new file mode 100644 index 0000000000000..fbe8d1c71d515 --- /dev/null +++ b/Detectors/Upgrades/ALICE3/IOTOF/reconstruction/include/IOTOFReconstruction/TopologyClassifier.h @@ -0,0 +1,147 @@ +// 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 + +#include "IOTOFBase/GeometryTGeo.h" +#include "IOTOFBase/IOTOFBaseParam.h" +#include "IOTOFBase/Segmentation.h" +#include "DataFormatsIOTOF/Cluster.h" + +// TO BE REMOVED BEFORE PUSH +#include "Framework/Logger.h" + +namespace o2 +{ +namespace iotof +{ + +enum Topologies : uint8_t { + kSingleDigit, + kLineOnRow, + kLineOnCol, + kSquare, + kRectangle, + kDiagonal, + kLowerTriangleLeft, + kLowerTriangleRight, + kUpperTriangleLeft, + kUpperTriangleRight, + kSnake, + kSnakeRefl, + kSnakeRot90, + kSnakeRot90Refl, + kHuge, + kOther, + kNTopologies +}; + +struct TopologyInfo { + int mSizeX = 0; + int mSizeZ = 0; + int mOffsetXToCOG = 0; + int mOffsetZToCOG = 0; + float mXMean = 0.f; + float mZMean = 0.f; + float mXSigma2 = 0.f; + float mZSigma2 = 0.f; + int mNPixels = 0; + int mFrequency = 0; + Topologies mTopology = Topologies::kNTopologies; + uint16_t mPattern; ///< Bitmask of fired pixels + + 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 +{ + 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() { + sSegmentation = o2::iotof::Segmentation::Instance(); + } + TopologyClassifier(std::unordered_map map) : mTopologyCache(std::move(map)) { + sSegmentation = o2::iotof::Segmentation::Instance(); + } + + 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, uint32_t& topology); + TopologyInfo getTopologyFeatures(uint32_t key); + void accountTopology(uint16_t bitmask, uint16_t minRow, uint8_t spanRow, uint16_t minCol, uint8_t spanCol); + void computeCOG(uint16_t bitmask, uint16_t minRow, uint8_t spanRow, uint16_t minCol, uint8_t spanCol, TopologyInfo& topoInfo); + + math_utils::Point3D getClusterCoordinates(const Cluster& cluster); + + void saveCacheToFile(const char* filename); + void print(); + + float getErrX(uint32_t pattID) {return std::sqrt(getTopologyFeatures(pattID).mXSigma2);}; + float getErrZ(uint32_t pattID) {return std::sqrt(getTopologyFeatures(pattID).mZSigma2);}; + float getNPixels(uint32_t pattID) {return getTopologyFeatures(pattID).mNPixels;}; + + // Provide the common iotof::GeometryTGeo to access matrices and segmentation + void setGeometry(const o2::iotof::GeometryTGeo* gm) { mGeometry = gm; } + + static uint32_t makeKey(uint8_t spanRow, uint8_t spanCol, uint16_t bitmask) { + return (static_cast(spanRow) << 24) | + (static_cast(spanCol) << 16) | + static_cast(bitmask); + } + + 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; + const o2::iotof::GeometryTGeo* mGeometry = nullptr; ///< IOTOF geometry + static o2::iotof::Segmentation* sSegmentation; ///< IOTOF segmentation instance (singleton) + +}; + +} // 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 bb80ae17d2f62..a9c1ea579f14d 100644 --- a/Detectors/Upgrades/ALICE3/IOTOF/reconstruction/src/Clusterer.cxx +++ b/Detectors/Upgrades/ALICE3/IOTOF/reconstruction/src/Clusterer.cxx @@ -33,28 +33,27 @@ void Clusterer::process(gsl::span digits, gsl::span digMC2ROFs, std::vector* clusterMC2ROFs) { - LOG(info) << "Running clusterizer on " << digitROFs.size() << " ROFs, total digits: " << digits.size(); + LOG(info) << "RUNNING CLUSTERIZER ON " << digitROFs.size() << " ROFs, TOTAL DIGITS: " << digits.size(); if (!mThread) { mThread = std::make_unique(this); } 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(debug) << "Processing 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) { + 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,150 +61,310 @@ 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(); + LOG(info) << "FINISHED PROCESSING ALL DIGIT ROFS, TOTAL CLUSTERS PRODUCED: " << clusters.size(); if (clusterMC2ROFs && !digMC2ROFs.empty()) { clusterMC2ROFs->reserve(clusterMC2ROFs->size() + digMC2ROFs.size()); for (const auto& in : digMC2ROFs) { 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; - // 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) { + findClustersSingleHit(digits, sortIdx[firstDigitIdx], labelsDigPtr, labelsClusPtr); + } else { + std::vector digitIdxs(nDigits); - // // 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 - // } + for (int i = 0; i < nDigits; ++i) { + digitIdxs[i] = sortIdx[firstDigitIdx + i]; + } + + findClustersMultipleHits( + digits, + gsl::span(digitIdxs), + labelsDigPtr, + labelsClusPtr); + } // 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 = 0; i < nlab; i++) { - labels.addElement(cnt, labelsBuff[i]); + int nStoredLabels = 0; + fetchMCLabels(digitIdx, labelsDigPtr, nStoredLabels); + const auto nCls = static_cast(mClusters.size()); + for (int i = 0; i < nStoredLabels; i++) { + mLabels.addElement(nCls, 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}; + uint32_t 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(debug) << "Pushing back cluster with row: " << row << ", col: " << col << ", rowSpan: " << rowSpan + << ", colSpan: " << colSpan << ", pattern: " << firedDigitsMask + << ", topology: " << clsTopology << ", chipID: " << chipID + << ", time: " << time; + + mClusters.emplace_back(cluster); + mPatterns.emplace_back(static_cast(firedDigitsMask)); } //__________________________________________________ -void Clusterer::ClustererThread::fetchMCLabels(uint32_t digID, const ConstDigitTruth* labelsDig, int& nfilled) +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}); + } + } + + for (const auto& precluster : preclusters) { + + const auto nStoredCls = static_cast(mClusters.size()); + + // Single-digit cluster in chip with multiple fired digits + if (precluster.size() == 1) { + 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}; + uint32_t 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(minRow, minCol, rowSpan, colSpan, firedDigitsMask, clsTopology, chipID, time); + + LOG(debug) << "Pushing back cluster with row: " << row << ", col: " << col << ", rowSpan: " << rowSpan + << ", colSpan: " << colSpan << ", pattern: " << firedDigitsMask + << ", topology: " << clsTopology << ", chipID: " << chipID + << ", time: " << time; + + mClusters.emplace_back(cluster); + mPatterns.emplace_back(static_cast(firedDigitsMask)); + } else { + // 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) { + // 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); + mPatterns.emplace_back(Cluster::InvalidPatternID); + 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); + } + } + + uint32_t clsTopology{0}; + mClsTopoClassifier.getTopology(firedDigitsMask, minRow, rowSpan, minCol, colSpan, clsTopology); + + // Construct and add cluster using scalar pattern mask + for (int i = 0; i < nMcLabels; i++) { + mLabels.addElement(nStoredCls, mLabelsBuff[i]); + } + Cluster cluster(minRow, minCol, rowSpan, colSpan, firedDigitsMask, clsTopology, chipID, clsTime); + LOG(debug) << "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); + mPatterns.emplace_back(static_cast(firedDigitsMask)); + } + } +} + +//__________________________________________________ +void Clusterer::ClustererThread::fetchMCLabels(uint32_t digID, const ConstDigitTruth* labelsDig, int& nFilled) { if (!labelsDig || digID >= labelsDig->getIndexedSize()) { return; } - auto sortBuffer = [this]() { std::sort(this->labelsBuff.begin(), this->labelsBuff.end(), [](Label const& a, Label const& b) { return a.getTrackID() < b.getTrackID(); }); }; + auto sortBuffer = [this]() { std::sort(this->mLabelsBuff.begin(), this->mLabelsBuff.end(), [](Label const& a, Label const& b) { return a.getTrackID() < b.getTrackID(); }); }; for (const auto& label : labelsDig->getLabels(digID)) { bool skip = false; - for (int ic = 0; ic < nfilled; ic++) { - if (labelsBuff[ic] == label) { + for (int ic = 0; ic < nFilled; ic++) { + if (mLabelsBuff[ic] == label) { skip = true; break; } } if (!skip) { - if (nfilled < MaxLabels) { - labelsBuff[nfilled++] = label; - if (nfilled == MaxLabels) { + if (nFilled < MaxLabels) { + mLabelsBuff[nFilled++] = label; + if (nFilled == MaxLabels) { sortBuffer(); } - } else if (labelsBuff.back().getTrackID() > label.getTrackID()) { - labelsBuff.back() = label; + } else if (mLabelsBuff.back().getTrackID() > label.getTrackID()) { + mLabelsBuff.back() = label; sortBuffer(); } } } } +//__________________________________________________ +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..c38b8b7f02d9a --- /dev/null +++ b/Detectors/Upgrades/ALICE3/IOTOF/reconstruction/src/IOTOFReconstructionLinkDef.h @@ -0,0 +1,27 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +#ifdef __CLING__ + +#pragma link off all globals; +#pragma link off all classes; +#pragma link off all functions; + +#pragma link C++ class o2::iotof::Clusterer + ; + +#pragma link C++ class o2::iotof::ClustererParam + ; + +#pragma link C++ class o2::iotof::TopologyClassifier + ; + +#pragma link C++ class o2::iotof::TopologyInfo + ; +#pragma link C++ class std::unordered_map < uint32_t, o2::iotof::TopologyInfo> + ; + +#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..67ff0a7dffbe9 --- /dev/null +++ b/Detectors/Upgrades/ALICE3/IOTOF/reconstruction/src/TopologyClassifier.cxx @@ -0,0 +1,309 @@ +// 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 for bitset +#include + +ClassImp(o2::iotof::TopologyClassifier); + +using std::array; + +namespace o2 +{ +namespace iotof +{ + +o2::iotof::Segmentation* TopologyClassifier::sSegmentation = nullptr; + +void TopologyClassifier::getTopology(uint16_t bitmask, uint16_t minRow, uint8_t spanRow, uint16_t minCol, uint8_t spanCol, uint32_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); + + // Check if the topology is already cached + auto it = mTopologyCache.find(clsTopoKey); + if (it != mTopologyCache.end()) { + topology = it->second.mTopology; + it->second.mFrequency++; + LOG(debug) << "Found cached topology: " << static_cast(topology); + return; + } + + // Classify the new topology and cache the result + accountTopology(bitmask, minRow, spanRow, minCol, spanCol); + topology = mTopologyCache[clsTopoKey].mTopology; + LOG(debug) << "Classified new topology: " << static_cast(topology); +} + +TopologyInfo TopologyClassifier::getTopologyFeatures(uint32_t key) +{ + auto it = mTopologyCache.find(key); + if (it != mTopologyCache.end()) { + return it->second; + } else { + LOG(debug) << "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) +{ + LOG(debug) << "Classifying topology for bitmask: " << std::bitset<16>(bitmask) << ", minRow: " + << static_cast(minRow) << ", spanRow: " << static_cast(spanRow) + << ", minCol: " << static_cast(minCol) << ", spanCol: " << static_cast(spanCol); + + // New cluster topology features + TopologyInfo newTopo; + newTopo.mFrequency = 1; + newTopo.mPattern = bitmask; + newTopo.mSizeX = spanRow; + newTopo.mSizeZ = spanCol; + float xCOG{0.f}, zCOG{0.f}, mXMean{0.f}, mZMean{0.f}, mXSigma2{0.f}, mZSigma2{0.f}; + computeCOG(bitmask, minRow, spanRow, minCol, spanCol, newTopo); + + const int maxRow = minRow + spanRow - 1; + const int maxCol = minCol + spanCol - 1; + + const auto hasDigit = [bitmask, minRow, minCol, spanCol](int row, int col) -> bool { + const int bitIndex = (row - minRow) * spanCol + (col - minCol); + return (bitmask & (1U << bitIndex)) != 0; + }; + + // Basic shapes + if (spanRow == 1 && spanCol == 1) { + newTopo.mTopology = Topologies::kSingleDigit; + mTopologyCache[packKey(spanRow, spanCol, bitmask)] = newTopo; + return; + } + if (spanCol == 1) { + newTopo.mTopology = Topologies::kLineOnRow; + mTopologyCache[packKey(spanRow, spanCol, bitmask)] = newTopo; + return; + } + if (spanRow == 1) { + newTopo.mTopology = Topologies::kLineOnCol; + mTopologyCache[packKey(spanRow, spanCol, bitmask)] = newTopo; + return; + } + + // Calculate total active digits in the cluster mask + int firedDigits = 0; + for (int r = minRow; r <= maxRow; ++r) { + for (int c = minCol; c <= maxCol; ++c) { + if (hasDigit(r, c)) + firedDigits++; + } + } + + // Square and rectangles: all pixels fired + if (firedDigits == spanRow * spanCol && spanRow == spanCol) { + newTopo.mTopology = Topologies::kSquare; + mTopologyCache[packKey(spanRow, spanCol, bitmask)] = newTopo; + return; + } + if (firedDigits == spanRow * spanCol && spanRow != spanCol) { + newTopo.mTopology = Topologies::kRectangle; + mTopologyCache[packKey(spanRow, spanCol, bitmask)] = newTopo; + return; + } + + // Corner occupancy + const bool hasBottomLeft = hasDigit(minRow, minCol); + const bool hasBottomRight = hasDigit(minRow, maxCol); + const bool hasTopLeft = hasDigit(maxRow, minCol); + const bool hasTopRight = hasDigit(maxRow, maxCol); + + // Diagonal and triangles + if (spanRow == spanCol) { + + // Triangles + const int nCorners = hasTopLeft + hasTopRight + hasBottomLeft + hasBottomRight; + if (nCorners == 3) { + const int missing = !hasTopLeft ? 0 : !hasTopRight ? 1 + : !hasBottomLeft ? 2 + : 3; + + switch (missing) { + case 0: + newTopo.mTopology = Topologies::kLowerTriangleLeft; + break; + case 1: + newTopo.mTopology = Topologies::kLowerTriangleRight; + break; + case 2: + newTopo.mTopology = Topologies::kUpperTriangleLeft; + break; + case 3: + newTopo.mTopology = Topologies::kUpperTriangleRight; + break; + } + mTopologyCache[packKey(spanRow, spanCol, bitmask)] = newTopo; + return; + } + + if ((firedDigits == spanRow && hasTopLeft && hasBottomRight && !hasTopRight && !hasBottomLeft) || + (firedDigits == spanRow && hasTopRight && hasBottomLeft && !hasTopLeft && !hasBottomRight)) { + newTopo.mTopology = Topologies::kDiagonal; + mTopologyCache[packKey(spanRow, spanCol, bitmask)] = newTopo; + return; + } + } + + // Snake: 3 x 2 + if (spanRow == 3 && spanCol == 2) { + const bool hasMiddleMin = hasDigit(minRow, minCol + 1); + const bool hasMiddleMax = hasDigit(minRow, maxCol + 1); + + if (hasMiddleMin && hasMiddleMax) { + if (!hasTopLeft && !hasBottomRight && hasTopRight && hasBottomLeft) { + newTopo.mTopology = Topologies::kSnake; + mTopologyCache[packKey(spanRow, spanCol, bitmask)] = newTopo; + return; + } + + if (hasTopLeft && hasBottomRight && !hasTopRight && !hasBottomLeft) { + newTopo.mTopology = Topologies::kSnakeRefl; + mTopologyCache[packKey(spanRow, spanCol, bitmask)] = newTopo; + return; + } + } + } + + // Snake rotated by 90 degrees: 2 x 3 + if (spanRow == 2 && spanCol == 3) { + const bool hasMiddleLeft = hasDigit(minRow + 1, minCol); + const bool hasMiddleRight = hasDigit(maxRow + 1, minCol); + + if (hasMiddleLeft && hasMiddleRight) { + if (!hasTopLeft && !hasBottomRight && hasTopRight && hasBottomLeft) { + newTopo.mTopology = Topologies::kSnakeRot90; + mTopologyCache[packKey(spanRow, spanCol, bitmask)] = newTopo; + return; + } + + if (hasTopLeft && hasBottomRight && !hasTopRight && !hasBottomLeft) { + newTopo.mTopology = Topologies::kSnakeRot90Refl; + mTopologyCache[packKey(spanRow, spanCol, bitmask)] = newTopo; + return; + } + } + } + + if (newTopo.mTopology == Topologies::kNTopologies) { + newTopo.mTopology = Topologies::kOther; + mTopologyCache[packKey(spanRow, spanCol, bitmask)] = newTopo; + return; + } +} + +void TopologyClassifier::computeCOG(uint16_t bitmask, uint16_t minRow, uint8_t spanRow, uint16_t minCol, uint8_t spanCol, TopologyInfo& topoInfo) +{ + 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; + ++firedPixels; + } + } + + topoInfo.mOffsetXToCOG = static_cast((static_cast(xOffsetCOG) / firedPixels) - static_cast(minRow)); + topoInfo.mOffsetZToCOG = static_cast((static_cast(zOffsetCOG) / firedPixels) - static_cast(minCol)); + topoInfo.mNPixels = firedPixels; + + const auto& chipSpecs = ChipSpecificsParam::Instance(); + topoInfo.mXMean = (static_cast(xOffsetCOG) / firedPixels - minRow) * chipSpecs.PitchRow; + topoInfo.mZMean = (static_cast(zOffsetCOG) / firedPixels - minCol) * chipSpecs.PitchCol; + topoInfo.mXSigma2 = chipSpecs.PitchRow * chipSpecs.PitchRow / 12. / topoInfo.mSizeX; + topoInfo.mZSigma2 = chipSpecs.PitchCol * chipSpecs.PitchCol / 12. / topoInfo.mSizeZ; + + LOG(debug) << "Computed topology features"; + LOG(debug) << "COG offsets: (" << topoInfo.mOffsetXToCOG << ", " << topoInfo.mOffsetZToCOG << ")"; + LOG(debug) << "Shifts to mean: (" << topoInfo.mXMean << ", " << topoInfo.mZMean << ")"; + LOG(debug) << "Sigmas: (" << topoInfo.mXSigma2 << ", " << topoInfo.mZSigma2 << ")"; + LOG(debug) << "Fired Pixels: " << firedPixels; +} + +math_utils::Point3D TopologyClassifier::getClusterCoordinates(const Cluster& cluster) +{ + if (!mGeometry) { + LOG(fatal) << "Geometry not set in TopologyClassifier, cannot execute getClusterCoordinates!"; + return math_utils::Point3D{0.f, 0.f, 0.f}; + } + if (!sSegmentation) { + LOG(fatal) << "Segmentation not set in TopologyClassifier, cannot execute getClusterCoordinates!"; + return math_utils::Point3D{0.f, 0.f, 0.f}; + } + auto refRow = cluster.getRow(); + auto refCol = cluster.getCol(); + float x{0.f}; + float z{0.f}; + int layer = mGeometry->getIOTOFLayer(cluster.getChipID()); + sSegmentation->detectorToLocal(cluster.getRow(), cluster.getCol(), x, z, layer); + + uint32_t topoKey = cluster.getTopology(); + x += this->getTopologyFeatures(topoKey).mXMean; + z += this->getTopologyFeatures(topoKey).mZMean; + math_utils::Point3D locCl{x, 0.f, z}; + + return locCl; +} + +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() +{ + 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); + topoInfo.print(); + } +} + +} // namespace 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 6d55e33d5461b..9dccbe67652c9 100644 --- a/Detectors/Upgrades/ALICE3/IOTOF/simulation/include/IOTOFSimulation/Digitizer.h +++ b/Detectors/Upgrades/ALICE3/IOTOF/simulation/include/IOTOFSimulation/Digitizer.h @@ -35,7 +35,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/Digitizer.cxx b/Detectors/Upgrades/ALICE3/IOTOF/simulation/src/Digitizer.cxx index 3070851431fee..913e7bcdf0865 100644 --- a/Detectors/Upgrades/ALICE3/IOTOF/simulation/src/Digitizer.cxx +++ b/Detectors/Upgrades/ALICE3/IOTOF/simulation/src/Digitizer.cxx @@ -25,6 +25,8 @@ #include #include + +#include #include #include #include @@ -159,7 +161,6 @@ void Digitizer::processHit(const o2::itsmft::Hit& hit, int evID, int srcID) // Apply efficiency cut based on the hit segment mean position relative to the pixel center sSegmentation->detectorToLocal(rowIS, colIS, xPixelCenter, zPixelCenter, subdetectorID); if (!isEfficient(avgHitLocalX[irow][icol] - xPixelCenter, avgHitLocalZ[irow][icol] - zPixelCenter)) { - LOG(debug) << "Hit rejected by efficiency cut at pixel (" << rowIS << ", " << colIS << ") in chip " << chipID; continue; } @@ -183,19 +184,22 @@ void Digitizer::processHit(const o2::itsmft::Hit& hit, int evID, int srcID) void Digitizer::stepping(const o2::itsmft::Hit& hit, float**& respMatrix, float**& avgHitLocalX, float**& avgHitLocalZ, int& rowStart, int& colStart, int& rowSpan, int& colSpan) { - LOG(debug) << "\n\nPerforming stepping"; + LOG(debug) << "Stepping through hit for detector ID: " << hit.GetDetectorID(); const int chipID = hit.GetDetectorID(); const auto& matrix = mGeometry->getMatrixL2G(chipID); const int subdetectorID = mGeometry->getIOTOFLayer(chipID); + LOG(debug) << "Transforming hit positions to sensor frame"; auto xyzPositionStart(matrix ^ (hit.GetPosStart())); // start position in sensor frame auto xyzPositionEnd(matrix ^ (hit.GetPos())); // end position in sensor frame + LOG(debug) << "Hit start position in sensor frame: (" << xyzPositionStart.X() << ", " << xyzPositionStart.Y() << ", " << xyzPositionStart.Z() << ")"; const auto& digitizerParams = o2::iotof::DPLDigitizerParam::Instance(); const auto stepVector = (xyzPositionEnd - xyzPositionStart) / digitizerParams.nSimSteps; xyzPositionStart = xyzPositionStart + stepVector * 0.5f; // center the start position in the middle of the step xyzPositionEnd = xyzPositionEnd - stepVector * 0.5f; // center the end position in the middle of the step + LOG(debug) << "Stepping vector: (" << stepVector.X() << ", " << stepVector.Y() << ", " << stepVector.Z() << ")"; rowStart = -1; colStart = -1; int rowEnd = -1, colEnd = -1, nSkip = 0, nSteps = digitizerParams.nSimSteps; @@ -206,6 +210,7 @@ void Digitizer::stepping(const o2::itsmft::Hit& hit, float**& respMatrix, float* } xyzPositionStart += stepVector; } + LOG(debug) << "Hit start position in sensor frame after adjustment: (" << xyzPositionStart.X() << ", " << xyzPositionStart.Y() << ", " << xyzPositionStart.Z() << ")"; while (!sSegmentation->localToDetector(xyzPositionEnd.X(), xyzPositionEnd.Z(), rowEnd, colEnd, mGeometry->getIOTOFLayer(chipID))) { if (++nSkip > digitizerParams.nSimSteps) { // additional check to add: should we exclude something? @@ -214,6 +219,25 @@ void Digitizer::stepping(const o2::itsmft::Hit& hit, float**& respMatrix, float* } xyzPositionEnd -= stepVector; } + LOG(debug) << "Hit end position in sensor frame after adjustment: (" << xyzPositionEnd.X() << ", " << xyzPositionEnd.Y() << ", " << xyzPositionEnd.Z() << ")"; + + std::set crossedRows, crossedCols; + for (int iStep = nSteps; iStep--;) { + auto pixelCurrentPosLocal = xyzPositionStart + stepVector * iStep; + int row, col; + if (sSegmentation->localToDetector(pixelCurrentPosLocal.X(), pixelCurrentPosLocal.Z(), row, col, subdetectorID)) { + crossedRows.insert(row); + crossedCols.insert(col); + } + } + LOG(debug) << "Crossed rows: "; + for (const auto& row : crossedRows) { + LOG(debug) << row; + } + LOG(debug) << "Crossed cols: "; + for (const auto& col : crossedCols) { + LOG(debug) << col; + } if (rowStart > rowEnd) { std::swap(rowStart, rowEnd); @@ -227,12 +251,17 @@ void Digitizer::stepping(const o2::itsmft::Hit& hit, float**& respMatrix, float* rowEnd += digitizerParams.responseMatrixSize / 2; rowStart = std::max(rowStart, 0); colStart = std::max(colStart, 0); + LOG(debug) << "Row range: [" << rowStart << ", " << rowEnd << "], Col range: [" << colStart << ", " << colEnd << "]"; const auto& specsConfig = ChipSpecificsParam::Instance(); rowEnd = std::min(rowEnd, (specsConfig.NRows) - 1); colEnd = std::min(colEnd, (specsConfig.NCols) - 1); rowSpan = rowEnd - rowStart + 1; colSpan = colEnd - colStart + 1; + if (rowSpan <= 0 || colSpan <= 0) { + return; + } + LOG(debug) << "Final row range: [" << rowStart << ", " << rowEnd << "], Col range: [" << colStart << ", " << colEnd << "]"; respMatrix = new float*[rowSpan]; avgHitLocalX = new float*[rowSpan]; @@ -242,13 +271,16 @@ void Digitizer::stepping(const o2::itsmft::Hit& hit, float**& respMatrix, float* avgHitLocalX[i] = new float[colSpan](); avgHitLocalZ[i] = new float[colSpan](); } + LOG(debug) << "Allocated response matrix and average hit position arrays with size (" << rowSpan << ", " << colSpan << ")"; - if (!respMatrix || !avgHitLocalX || !avgHitLocalZ || rowSpan <= 0 || colSpan <= 0) { + if (!respMatrix || !avgHitLocalX || !avgHitLocalZ) { return; } + LOG(debug) << "Starting stepping through the hit with " << nSteps << " steps"; if (nSkip) { nSteps -= nSkip; } + LOG(debug) << "Adjusted number of steps after skipping: " << nSteps; int rowPrev = -1, colPrev = -1, row = 0, col = 0; auto pixelCurrentPosLocal = xyzPositionStart; @@ -264,6 +296,7 @@ void Digitizer::stepping(const o2::itsmft::Hit& hit, float**& respMatrix, float* // The step has reached another pixel, compute mean hit segment positions // for pixel efficiency evaluation and reset the start position for the next pixel + LOG(debug) << "iStep: " << iStep << ", Current pixel: (row,col) = (" << row << ", " << col << "), Previous pixel: (rowPrev,colPrev) = (" << rowPrev << ", " << colPrev << ")"; if (row != rowPrev || col != colPrev) { // Finalize the previous pixel @@ -296,14 +329,22 @@ void Digitizer::stepping(const o2::itsmft::Hit& hit, float**& respMatrix, float* } } } + LOG(debug) << "Finished stepping through the hit for detector ID: " << chipID; + LOG(debug) << "rowPrev: " << rowPrev << ", colPrev: " << colPrev << ", rowStart: " << rowStart << ", colStart: " << colStart; // Finalize the last pixel if (rowPrev != -1 && colPrev != -1) { const int irow = rowPrev - rowStart; const int icol = colPrev - colStart; + // Sizes of avgHitLocalX, avgHitLocalZ + LOG(debug) << "avgHitLocalX dimensions: " << rowSpan << " x " << colSpan; + LOG(debug) << "avgHitLocalZ dimensions: " << rowSpan << " x " << colSpan; + LOG(debug) << "Finalizing last pixel at (row,col) = (" << rowPrev << ", " << colPrev << ") with indices (irow,icol) = (" << irow << ", " << icol << ")"; avgHitLocalX[irow][icol] = 0.5f * (pixelStartPosLocal.X() + pixelCurrentPosLocal.X() - stepVector.X()); avgHitLocalZ[irow][icol] = 0.5f * (pixelStartPosLocal.Z() + pixelCurrentPosLocal.Z() - stepVector.Z()); + LOG(debug) << "Finalized last pixel average positions: avgHitLocalX = " << avgHitLocalX[irow][icol] << ", avgHitLocalZ = " << avgHitLocalZ[irow][icol]; } + LOG(debug) << "Finalized last pixel for detector ID: " << chipID; } //_______________________________________________________________________ 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..87f82e8b86ff2 100644 --- a/Detectors/Upgrades/ALICE3/IOTOF/workflow/src/ClustererSpec.cxx +++ b/Detectors/Upgrades/ALICE3/IOTOF/workflow/src/ClustererSpec.cxx @@ -67,15 +67,17 @@ void ClustererDPL::run(o2::framework::ProcessingContext& pc) mUseMC ? &labels : nullptr, clusterLabels.get()); LOG(info) << "Clusterization produced " << clusters.size() << " clusters for layer " << iLayer; + LOG(info) << "Clusterization produced " << patterns.size() << " patterns for layer " << iLayer; + LOG(info) << "Clusterization produced " << clusterROFs.size() << " ROFs 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) { pc.outputs().snapshot(o2::framework::Output{"TF3", "CLUSTERSMCTR", subspec}, *clusterLabels); } totalClusters += clusters.size(); - LOGP(info, "Pushed {} clusters in {} ROFs for layer {}", clusters.size(), clusterROFs.size(), iLayer); + LOGP(info, "Pushed {} clusters, {} patterns, in {} ROFs for layer {}", clusters.size(), patterns.size(), clusterROFs.size(), iLayer); LOGP(info, "Pushed {} MC labels for layer {}", mUseMC ? clusterLabels->getNElements() : 0, iLayer); } @@ -92,7 +94,7 @@ o2::framework::DataProcessorSpec getClustererSpec(bool useMC) } std::vector outputs; - outputs.emplace_back("TF3", "COMPCLUSTERS", iLayer, o2::framework::Lifetime::Timeframe); + outputs.emplace_back("TF3", "CLUSTERS", iLayer, o2::framework::Lifetime::Timeframe); outputs.emplace_back("TF3", "PATTERNS", iLayer, o2::framework::Lifetime::Timeframe); outputs.emplace_back("TF3", "CLUSTERSROF", iLayer, o2::framework::Lifetime::Timeframe); if (useMC) { diff --git a/Framework/Core/src/CommonServices.cxx b/Framework/Core/src/CommonServices.cxx index c36a102bde80d..2b6d6023ac7d5 100644 --- a/Framework/Core/src/CommonServices.cxx +++ b/Framework/Core/src/CommonServices.cxx @@ -143,7 +143,7 @@ o2::framework::ServiceSpec CommonServices::monitoringSpec() // covers devices that quit themselves via readyToQuit(). .stop = [](ServiceRegistryRef, void* service) { auto* monitoring = reinterpret_cast(service); - monitoring->finalizeProcessMonitoring(); }, + monitoring->enableProcessMonitoring(); }, .exit = [](ServiceRegistryRef registry, void* service) { auto* monitoring = reinterpret_cast(service); monitoring->flushBuffer();