From e978883844b8bec03005ecbd036f752300028ab7 Mon Sep 17 00:00:00 2001 From: Felix Schlepper Date: Mon, 10 Feb 2025 11:31:58 +0100 Subject: [PATCH 01/10] ITS3: make segmentation fully constexpr Signed-off-by: Felix Schlepper --- .../include/ITS3Base/SegmentationMosaix.h | 60 +++++++++---------- 1 file changed, 28 insertions(+), 32 deletions(-) diff --git a/Detectors/Upgrades/ITS3/base/include/ITS3Base/SegmentationMosaix.h b/Detectors/Upgrades/ITS3/base/include/ITS3Base/SegmentationMosaix.h index 37bb6156685cd..cecec27c95db4 100644 --- a/Detectors/Upgrades/ITS3/base/include/ITS3Base/SegmentationMosaix.h +++ b/Detectors/Upgrades/ITS3/base/include/ITS3Base/SegmentationMosaix.h @@ -52,18 +52,20 @@ class SegmentationMosaix // | | | // x----------------------x public: - ~SegmentationMosaix() = default; - SegmentationMosaix(const SegmentationMosaix&) = default; - SegmentationMosaix(SegmentationMosaix&&) = delete; - SegmentationMosaix& operator=(const SegmentationMosaix&) = default; - SegmentationMosaix& operator=(SegmentationMosaix&&) = delete; - constexpr SegmentationMosaix(int layer) : mLayer{layer} {} + constexpr SegmentationMosaix(int layer) : mRadius(static_cast(constants::radiiMiddle[layer])) {} + constexpr ~SegmentationMosaix() = default; + constexpr SegmentationMosaix(const SegmentationMosaix&) = default; + constexpr SegmentationMosaix(SegmentationMosaix&&) = delete; + constexpr SegmentationMosaix& operator=(const SegmentationMosaix&) = default; + constexpr SegmentationMosaix& operator=(SegmentationMosaix&&) = delete; static constexpr int NCols{constants::pixelarray::nCols}; static constexpr int NRows{constants::pixelarray::nRows}; static constexpr int NPixels{NCols * NRows}; static constexpr float Length{constants::pixelarray::length}; + static constexpr float LengthH{Length / 2.f}; static constexpr float Width{constants::pixelarray::width}; + static constexpr float WidthH{Width / 2.f}; static constexpr float PitchCol{constants::pixelarray::length / static_cast(NCols)}; static constexpr float PitchRow{constants::pixelarray::width / static_cast(NRows)}; static constexpr float SensorLayerThickness{constants::totalThickness}; @@ -77,16 +79,16 @@ class SegmentationMosaix /// the center of the sensitive volume. /// \param yFlat Detector local flat coordinate y in cm with respect to /// the center of the sensitive volume. - void curvedToFlat(const float xCurved, const float yCurved, float& xFlat, float& yFlat) const noexcept + constexpr void curvedToFlat(const float xCurved, const float yCurved, float& xFlat, float& yFlat) const noexcept { // MUST align the flat surface with the curved surface with the original pixel array is on and account for metal // stack float dist = std::hypot(xCurved, yCurved); float phi = std::atan2(yCurved, xCurved); - xFlat = (getRadius() * phi) - Width / 2.f; + xFlat = (mRadius * phi) - WidthH; // the y position is in the silicon volume however we need the chip volume (silicon+metalstack) // this is accounted by a y shift - yFlat = dist - getRadius() + static_cast(constants::nominalYShift); + yFlat = dist + (static_cast(constants::nominalYShift) - mRadius); } /// Transformation from the flat surface to a curved surface @@ -99,15 +101,15 @@ class SegmentationMosaix /// the center of the sensitive volume. /// \param yCurved Detector local curved coordinate y in cm with respect to /// the center of the sensitive volume. - void flatToCurved(float xFlat, float yFlat, float& xCurved, float& yCurved) const noexcept + constexpr void flatToCurved(float xFlat, float yFlat, float& xCurved, float& yCurved) const noexcept { // MUST align the flat surface with the curved surface with the original pixel array is on and account for metal // stack // the y position is in the chip volume however we need the silicon volume // this is accounted by a -y shift - float dist = yFlat + getRadius() - static_cast(constants::nominalYShift); - xCurved = dist * std::cos((xFlat + Width / 2.f) / getRadius()); - yCurved = dist * std::sin((xFlat + Width / 2.f) / getRadius()); + float dist = yFlat + (mRadius - static_cast(constants::nominalYShift)); + xCurved = dist * std::cos((xFlat + WidthH) / mRadius); + yCurved = dist * std::sin((xFlat + WidthH) / mRadius); } /// Transformation from Geant detector centered local coordinates (cm) to @@ -121,7 +123,7 @@ class SegmentationMosaix /// the center of the sensitive volume. /// \param int iRow Detector x cell coordinate. /// \param int iCol Detector z cell coordinate. - bool localToDetector(float const xRow, float const zCol, int& iRow, int& iCol) const noexcept + constexpr bool localToDetector(float const xRow, float const zCol, int& iRow, int& iCol) const noexcept { localToDetectorUnchecked(xRow, zCol, iRow, iCol); if (!isValid(iRow, iCol)) { @@ -132,10 +134,10 @@ class SegmentationMosaix } // Same as localToDetector w.o. checks. - void localToDetectorUnchecked(float const xRow, float const zCol, int& iRow, int& iCol) const noexcept + constexpr void localToDetectorUnchecked(float const xRow, float const zCol, int& iRow, int& iCol) const noexcept { - iRow = std::floor((Width / 2. - xRow) / PitchRow); - iCol = std::floor((zCol + Length / 2.) / PitchCol); + iRow = static_cast(std::floor((WidthH - xRow) / PitchRow)); + iCol = static_cast(std::floor((zCol + LengthH) / PitchCol)); } /// Transformation from Detector cell coordinates to Geant detector centered @@ -148,7 +150,7 @@ class SegmentationMosaix /// center of the sensitive volume. /// If iRow and or iCol is outside of the segmentation range a value of -0.5*Dx() /// or -0.5*Dz() is returned. - bool detectorToLocal(int const iRow, int const iCol, float& xRow, float& zCol) const noexcept + constexpr bool detectorToLocal(int const iRow, int const iCol, float& xRow, float& zCol) const noexcept { if (!isValid(iRow, iCol)) { return false; @@ -159,13 +161,13 @@ class SegmentationMosaix // Same as detectorToLocal w.o. checks. // We position ourself in the middle of the pixel. - void detectorToLocalUnchecked(int const iRow, int const iCol, float& xRow, float& zCol) const noexcept + constexpr void detectorToLocalUnchecked(int const iRow, int const iCol, float& xRow, float& zCol) const noexcept { - xRow = -(static_cast(iRow) + 0.5f) * PitchRow + Width / 2.f; - zCol = (static_cast(iCol) + 0.5f) * PitchCol - Length / 2.f; + xRow = -(static_cast(iRow) + 0.5f) * PitchRow + WidthH; + zCol = (static_cast(iCol) + 0.5f) * PitchCol - LengthH; } - bool detectorToLocal(int const row, int const col, math_utils::Point3D& loc) const noexcept + bool detectorToLocal(int const row, int const col, math_utils::Point3D& loc) const noexcept { float xRow{0.}, zCol{0.}; if (!detectorToLocal(row, col, xRow, zCol)) { @@ -175,7 +177,7 @@ class SegmentationMosaix return true; } - void detectorToLocalUnchecked(int const row, int const col, math_utils::Point3D& loc) const noexcept + void detectorToLocalUnchecked(int const row, int const col, math_utils::Point3D& loc) const noexcept { float xRow{0.}, zCol{0.}; detectorToLocalUnchecked(row, col, xRow, zCol); @@ -184,22 +186,16 @@ class SegmentationMosaix private: template - [[nodiscard]] bool isValid(T const row, T const col) const noexcept + [[nodiscard]] constexpr bool isValid(T const row, T const col) const noexcept { if constexpr (std::is_floating_point_v) { // compares in local coord. - namespace cp = constants::pixelarray; - return !static_cast(row <= -cp::width / 2. || cp::width / 2. <= row || col <= -cp::length / 2. || cp::length / 2. <= col); + return (-WidthH < row && row < WidthH && -LengthH < col && col < LengthH); } else { // compares in rows/cols return !static_cast(row < 0 || row >= static_cast(NRows) || col < 0 || col >= static_cast(NCols)); } } - float getRadius() const noexcept - { - return static_cast(constants::radiiMiddle[mLayer]); - } - - int mLayer{0}; ///< chip layer + float mRadius; }; } // namespace o2::its3 From 9b1c1e7a44648313cf0b1be9ebda0ec8c3eefea0 Mon Sep 17 00:00:00 2001 From: Felix Schlepper Date: Mon, 10 Feb 2025 11:32:37 +0100 Subject: [PATCH 02/10] ITS3: minor changes in Digitizer for float Signed-off-by: Felix Schlepper --- .../include/ITS3Simulation/Digitizer.h | 10 ++++----- .../ITS3/simulation/src/Digitizer.cxx | 22 +++++++++---------- 2 files changed, 15 insertions(+), 17 deletions(-) diff --git a/Detectors/Upgrades/ITS3/simulation/include/ITS3Simulation/Digitizer.h b/Detectors/Upgrades/ITS3/simulation/include/ITS3Simulation/Digitizer.h index 508261f03ce38..efdf4f8096f82 100644 --- a/Detectors/Upgrades/ITS3/simulation/include/ITS3Simulation/Digitizer.h +++ b/Detectors/Upgrades/ITS3/simulation/include/ITS3Simulation/Digitizer.h @@ -106,14 +106,14 @@ class Digitizer : public TObject uint32_t mEventROFrameMin = 0xffffffff; ///< lowest RO frame for processed events (w/o automatic noise ROFs) uint32_t mEventROFrameMax = 0; ///< highest RO frame forfor processed events (w/o automatic noise ROFs) - const std::array mIBSegmentations{0, 1, 2}; + static constexpr std::array mIBSegmentations{0, 1, 2}; o2::itsmft::AlpideSimResponse* mSimRespIB = nullptr; // simulated response for IB o2::itsmft::AlpideSimResponse* mSimRespOB = nullptr; // simulated response for OB - float mSimRespIBShift{0.}; // adjusting the Y-shift in the IB response function to match sensor local coord. - float mSimRespIBScaleX{1.}; // scale x-local coordinate to response function x-coordinate - float mSimRespIBScaleZ{1.}; // scale z-local coordinate to response function z-coordinate - float mSimRespOBShift{0.}; // adjusting the Y-shift in the OB response function to match sensor local coord. + float mSimRespIBShift{0.f}; // adjusting the Y-shift in the IB response function to match sensor local coord. + float mSimRespIBScaleX{1.f}; // scale x-local coordinate to response function x-coordinate + float mSimRespIBScaleZ{1.f}; // scale z-local coordinate to response function z-coordinate + float mSimRespOBShift{0.f}; // adjusting the Y-shift in the OB response function to match sensor local coord. const o2::its::GeometryTGeo* mGeometry = nullptr; ///< ITS3 geometry diff --git a/Detectors/Upgrades/ITS3/simulation/src/Digitizer.cxx b/Detectors/Upgrades/ITS3/simulation/src/Digitizer.cxx index 74c7f31dad44a..53fee11421f0a 100644 --- a/Detectors/Upgrades/ITS3/simulation/src/Digitizer.cxx +++ b/Detectors/Upgrades/ITS3/simulation/src/Digitizer.cxx @@ -65,13 +65,13 @@ void Digitizer::init() if (const auto& func = ITS3Params::Instance().chipResponseFunction; func == "Alpide") { constexpr const char* responseFile = "$(O2_ROOT)/share/Detectors/ITSMFT/data/AlpideResponseData/AlpideResponseData.root"; loadSetResponseFunc("Alpide", responseFile, "response0", responseFile, "response1"); - mSimRespIBShift = mSimRespIB->getDepthMax() - SegmentationMosaix::SensorLayerThickness / 2.f + 10.e-4f; // TODO why this offset? + mSimRespIBShift = mSimRespIB->getDepthMax() - SegmentationMosaix::SensorLayerThickness / 2.f + 10.e-4f; mSimRespOBShift = mSimRespOB->getDepthMax() - SegmentationAlpide::SensorLayerThickness / 2.f; } else if (func == "APTS") { constexpr const char* responseFileIB = "$(O2_ROOT)/share/Detectors/Upgrades/ITS3/data/ITS3ChipResponseData/APTSResponseData.root"; constexpr const char* responseFileOB = "$(O2_ROOT)/share/Detectors/ITSMFT/data/AlpideResponseData/AlpideResponseData.root"; loadSetResponseFunc("APTS", responseFileIB, "response1", responseFileOB, "response1"); - mSimRespIBShift = mSimRespIB->getDepthMax() - 10.e-4f; + mSimRespIBShift = mSimRespIB->getDepthMax() - 6.5e-4f; mSimRespOBShift = mSimRespOB->getDepthMax() - SegmentationAlpide::SensorLayerThickness / 2.f; mSimRespIBScaleX = 0.5 * constants::pixelarray::pixels::apts::pitchX / SegmentationMosaix::PitchRow; mSimRespIBScaleZ = 0.5 * constants::pixelarray::pixels::apts::pitchZ / SegmentationMosaix::PitchCol; @@ -328,15 +328,13 @@ void Digitizer::processHit(const o2::itsmft::Hit& hit, uint32_t& maxFr, int evID const int maxNrows{innerBarrel ? SegmentationMosaix::NRows : SegmentationAlpide::NRows}; const int maxNcols{innerBarrel ? SegmentationMosaix::NCols : SegmentationAlpide::NCols}; - if (rowE >= maxNrows) { - rowE = maxNrows - 1; - } + + rowE = std::min(rowE, maxNrows - 1); colS -= AlpideRespSimMat::NPix / 2; colE += AlpideRespSimMat::NPix / 2; colS = std::max(colS, 0); - if (colE >= maxNcols) { - colE = maxNcols - 1; - } + colE = std::min(colE, maxNcols - 1); + int rowSpan = rowE - rowS + 1, colSpan = colE - colS + 1; // size of plaquet where some response is expected float respMatrix[rowSpan][colSpan]; // response accumulated here std::fill(&respMatrix[0][0], &respMatrix[0][0] + rowSpan * colSpan, 0.f); @@ -379,12 +377,12 @@ void Digitizer::processHit(const o2::itsmft::Hit& hit, uint32_t& maxFr, int evID float rowMax{}, colMax{}; const AlpideRespSimMat* rspmat{nullptr}; if (innerBarrel) { - rowMax = 0.5 * SegmentationMosaix::PitchRow; - colMax = 0.5 * SegmentationMosaix::PitchCol; + rowMax = 0.5f * SegmentationMosaix::PitchRow; + colMax = 0.5f * SegmentationMosaix::PitchCol; rspmat = mSimRespIB->getResponse(mSimRespIBScaleX * (xyzLocS.X() - cRowPix), mSimRespIBScaleZ * (xyzLocS.Z() - cColPix), xyzLocS.Y(), flipRow, flipCol, rowMax, colMax); } else { - rowMax = 0.5 * SegmentationAlpide::PitchRow; - colMax = 0.5 * SegmentationAlpide::PitchCol; + rowMax = 0.5f * SegmentationAlpide::PitchRow; + colMax = 0.5f * SegmentationAlpide::PitchCol; rspmat = mSimRespOB->getResponse(xyzLocS.X() - cRowPix, xyzLocS.Z() - cColPix, xyzLocS.Y(), flipRow, flipCol, rowMax, colMax); } From d5e54950d5bad22abeda0dc1419d4716cc94a975 Mon Sep 17 00:00:00 2001 From: Felix Schlepper Date: Tue, 11 Feb 2025 10:57:57 +0100 Subject: [PATCH 03/10] ITS3: move pixel specs down Signed-off-by: Felix Schlepper --- .../ITS3/base/include/ITS3Base/SpecsV2.h | 72 +++++++++++-------- 1 file changed, 44 insertions(+), 28 deletions(-) diff --git a/Detectors/Upgrades/ITS3/base/include/ITS3Base/SpecsV2.h b/Detectors/Upgrades/ITS3/base/include/ITS3Base/SpecsV2.h index 2b460d5c029d0..35c3db41ac77a 100644 --- a/Detectors/Upgrades/ITS3/base/include/ITS3Base/SpecsV2.h +++ b/Detectors/Upgrades/ITS3/base/include/ITS3Base/SpecsV2.h @@ -21,6 +21,11 @@ #include +// This files defines the design specifications of the chip. +// Each TGeoShape has the following properties +// length: dimension in z-axis +// width: dimension in xy-axes +// color: for visulisation namespace o2::its3::constants { constexpr double cm{1e+2}; // This is the default unit of TGeo so we use this as scale @@ -35,27 +40,6 @@ constexpr int nRows{442}; constexpr int nPixels{nRows * nCols}; constexpr EColor color{kGreen}; constexpr double area{width * length}; -namespace pixels -{ -namespace apts -{ -constexpr double pitchX{15.0 * mu}; -constexpr double pitchZ{15.0 * mu}; -} // namespace apts -namespace moss -{ -namespace top -{ -constexpr double pitchX{22.5 * mu}; -constexpr double pitchZ{22.5 * mu}; -} // namespace top -namespace bot -{ -constexpr double pitchX{18.0 * mu}; -constexpr double pitchZ{18.0 * mu}; -} // namespace bot -} // namespace moss -} // namespace pixels } // namespace pixelarray namespace tile { @@ -112,7 +96,7 @@ constexpr EColor color{kCyan}; } // namespace rec constexpr unsigned int nRSUs{12}; constexpr unsigned int nTilesPerSegment{nRSUs * rsu::nTiles}; -constexpr double length{nRSUs * rsu::length + lec::length + rec::length}; +constexpr double length{(nRSUs * rsu::length) + lec::length + rec::length}; constexpr double lengthSensitive{nRSUs * rsu::length}; } // namespace segment namespace carbonfoam @@ -136,9 +120,9 @@ constexpr EColor color{kBlack}; } // namespace metalstack namespace silicon { -constexpr double thickness{45 * mu}; // thickness of silicion -constexpr double radiusInner{(thickness + metalstack::thickness) / 2.}; // inner silicion radius correction -constexpr double radiusOuter{(thickness - metalstack::thickness) / 2.}; // outer silicion radius correction +constexpr double thickness{45 * mu}; // thickness of silicon +constexpr double thicknessIn{(thickness + metalstack::thickness) / 2.}; // inner silicon thickness +constexpr double thicknessOut{(thickness - metalstack::thickness) / 2.}; // outer silicon thickness } // namespace silicon constexpr unsigned int nLayers{3}; constexpr unsigned int nTotLayers{7}; @@ -147,10 +131,42 @@ constexpr double equatorialGap{1 * mm}; constexpr std::array nSegments{3, 4, 5}; constexpr double totalThickness{silicon::thickness + metalstack::thickness}; // total chip thickness constexpr std::array radii{19.0006 * mm, 25.228 * mm, 31.4554 * mm}; // nominal radius -constexpr std::array radiiInner{radii[0] - silicon::radiusInner, radii[1] - silicon::radiusInner, radii[2] - silicon::radiusInner}; // inner silicon radius -constexpr std::array radiiOuter{radii[0] + silicon::radiusOuter, radii[1] + silicon::radiusOuter, radii[2] + silicon::radiusOuter}; // outer silicon radius +constexpr std::array radiiInner{radii[0] - silicon::thicknessIn, radii[1] - silicon::thicknessIn, radii[2] - silicon::thicknessIn}; // inner silicon radius +constexpr std::array radiiOuter{radii[0] + silicon::thicknessOut, radii[1] + silicon::thicknessOut, radii[2] + silicon::thicknessOut}; // outer silicon radius constexpr std::array radiiMiddle{(radiiInner[0] + radiiOuter[0]) / 2., (radiiInner[1] + radiiOuter[1]) / 2., (radiiInner[2] + radiiOuter[2]) / 2.}; // middle silicon radius -constexpr double nominalYShift{-metalstack::thickness / 2.}; // shift to position in silicion volume to the chip volume (silicon+metalstack) +/*constexpr double nominalYShift{-metalstack::thickness / 2.}; // shift to position in silicion volume to the chip volume (silicon+metalstack)*/ +constexpr double nominalYShift{0}; // shift to position in silicion volume to the chip volume (silicon+metalstack) + +// extra information of pixels and their response functions +namespace pixelarray::pixels +{ +namespace mosaix +{ +constexpr double pitchX{width / static_cast(nRows)}; +constexpr double pitchZ{length / static_cast(nCols)}; +} // namespace mosaix +namespace apts +{ +constexpr double pitchX{15.0 * mu}; +constexpr double pitchZ{15.0 * mu}; +constexpr double responseUpperLimit{10 * mu}; +constexpr double responseYShift{responseUpperLimit - silicon::thicknessOut}; +} // namespace apts +namespace moss +{ +namespace top +{ +constexpr double pitchX{22.5 * mu}; +constexpr double pitchZ{22.5 * mu}; +} // namespace top +namespace bot +{ +constexpr double pitchX{18.0 * mu}; +constexpr double pitchZ{18.0 * mu}; +} // namespace bot +} // namespace moss +} // namespace pixelarray::pixels + namespace detID { constexpr unsigned int mDetIDs{2 * 12 * 12 * 12}; //< 2 Hemispheres * (3,4,5=12 segments in a layer) * 12 RSUs in a segment * 12 Tiles in a RSU From d11432befcba4f6b4280f5da9ca792a3efae51e4 Mon Sep 17 00:00:00 2001 From: Felix Schlepper Date: Tue, 11 Feb 2025 10:59:12 +0100 Subject: [PATCH 04/10] ITS3: account in segmentation trans for nominal shift Signed-off-by: Felix Schlepper --- .../ITS3/base/include/ITS3Base/SegmentationMosaix.h | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/Detectors/Upgrades/ITS3/base/include/ITS3Base/SegmentationMosaix.h b/Detectors/Upgrades/ITS3/base/include/ITS3Base/SegmentationMosaix.h index cecec27c95db4..f4c61db5fc25f 100644 --- a/Detectors/Upgrades/ITS3/base/include/ITS3Base/SegmentationMosaix.h +++ b/Detectors/Upgrades/ITS3/base/include/ITS3Base/SegmentationMosaix.h @@ -66,9 +66,10 @@ class SegmentationMosaix static constexpr float LengthH{Length / 2.f}; static constexpr float Width{constants::pixelarray::width}; static constexpr float WidthH{Width / 2.f}; - static constexpr float PitchCol{constants::pixelarray::length / static_cast(NCols)}; - static constexpr float PitchRow{constants::pixelarray::width / static_cast(NRows)}; + static constexpr float PitchCol{constants::pixelarray::pixels::mosaix::pitchZ}; + static constexpr float PitchRow{constants::pixelarray::pixels::mosaix::pitchX}; static constexpr float SensorLayerThickness{constants::totalThickness}; + static constexpr float NominalYShift{constants::nominalYShift}; /// Transformation from the curved surface to a flat surface /// \param xCurved Detector local curved coordinate x in cm with respect to @@ -88,7 +89,7 @@ class SegmentationMosaix xFlat = (mRadius * phi) - WidthH; // the y position is in the silicon volume however we need the chip volume (silicon+metalstack) // this is accounted by a y shift - yFlat = dist + (static_cast(constants::nominalYShift) - mRadius); + yFlat = dist - mRadius + NominalYShift; } /// Transformation from the flat surface to a curved surface @@ -107,7 +108,7 @@ class SegmentationMosaix // stack // the y position is in the chip volume however we need the silicon volume // this is accounted by a -y shift - float dist = yFlat + (mRadius - static_cast(constants::nominalYShift)); + float dist = yFlat - NominalYShift + mRadius; xCurved = dist * std::cos((xFlat + WidthH) / mRadius); yCurved = dist * std::sin((xFlat + WidthH) / mRadius); } @@ -173,7 +174,7 @@ class SegmentationMosaix if (!detectorToLocal(row, col, xRow, zCol)) { return false; } - loc.SetCoordinates(xRow, 0., zCol); + loc.SetCoordinates(xRow, NominalYShift, zCol); return true; } @@ -181,7 +182,7 @@ class SegmentationMosaix { float xRow{0.}, zCol{0.}; detectorToLocalUnchecked(row, col, xRow, zCol); - loc.SetCoordinates(xRow, 0., zCol); + loc.SetCoordinates(xRow, NominalYShift, zCol); } private: From d6e65db70856ef460bb386930bbe15e2286edfd3 Mon Sep 17 00:00:00 2001 From: Felix Schlepper Date: Tue, 11 Feb 2025 10:59:35 +0100 Subject: [PATCH 05/10] ITS3: Digitizer flip row for IB for APTS Signed-off-by: Felix Schlepper --- .../ITS3/simulation/include/ITS3Simulation/Digitizer.h | 1 + Detectors/Upgrades/ITS3/simulation/src/Digitizer.cxx | 9 +++++---- 2 files changed, 6 insertions(+), 4 deletions(-) diff --git a/Detectors/Upgrades/ITS3/simulation/include/ITS3Simulation/Digitizer.h b/Detectors/Upgrades/ITS3/simulation/include/ITS3Simulation/Digitizer.h index efdf4f8096f82..8d0f06a27343b 100644 --- a/Detectors/Upgrades/ITS3/simulation/include/ITS3Simulation/Digitizer.h +++ b/Detectors/Upgrades/ITS3/simulation/include/ITS3Simulation/Digitizer.h @@ -110,6 +110,7 @@ class Digitizer : public TObject o2::itsmft::AlpideSimResponse* mSimRespIB = nullptr; // simulated response for IB o2::itsmft::AlpideSimResponse* mSimRespOB = nullptr; // simulated response for OB + bool mSimRespIBOrientation{false}; // wether the orientation in the IB response function is flipped float mSimRespIBShift{0.f}; // adjusting the Y-shift in the IB response function to match sensor local coord. float mSimRespIBScaleX{1.f}; // scale x-local coordinate to response function x-coordinate float mSimRespIBScaleZ{1.f}; // scale z-local coordinate to response function z-coordinate diff --git a/Detectors/Upgrades/ITS3/simulation/src/Digitizer.cxx b/Detectors/Upgrades/ITS3/simulation/src/Digitizer.cxx index 53fee11421f0a..3c75bf3e8f680 100644 --- a/Detectors/Upgrades/ITS3/simulation/src/Digitizer.cxx +++ b/Detectors/Upgrades/ITS3/simulation/src/Digitizer.cxx @@ -71,10 +71,11 @@ void Digitizer::init() constexpr const char* responseFileIB = "$(O2_ROOT)/share/Detectors/Upgrades/ITS3/data/ITS3ChipResponseData/APTSResponseData.root"; constexpr const char* responseFileOB = "$(O2_ROOT)/share/Detectors/ITSMFT/data/AlpideResponseData/AlpideResponseData.root"; loadSetResponseFunc("APTS", responseFileIB, "response1", responseFileOB, "response1"); - mSimRespIBShift = mSimRespIB->getDepthMax() - 6.5e-4f; + mSimRespIBShift = mSimRespIB->getDepthMax() + (float)constants::pixelarray::pixels::apts::responseYShift; mSimRespOBShift = mSimRespOB->getDepthMax() - SegmentationAlpide::SensorLayerThickness / 2.f; - mSimRespIBScaleX = 0.5 * constants::pixelarray::pixels::apts::pitchX / SegmentationMosaix::PitchRow; - mSimRespIBScaleZ = 0.5 * constants::pixelarray::pixels::apts::pitchZ / SegmentationMosaix::PitchCol; + mSimRespIBScaleX = 0.5f * constants::pixelarray::pixels::apts::pitchX / SegmentationMosaix::PitchRow; + mSimRespIBScaleZ = 0.5f * constants::pixelarray::pixels::apts::pitchZ / SegmentationMosaix::PitchCol; + mSimRespIBOrientation = true; } else { LOGP(fatal, "ResponseFunction '{}' not implemented!", func); } @@ -401,7 +402,7 @@ void Digitizer::processHit(const o2::itsmft::Hit& hit, uint32_t& maxFr, int evID if (colDest < 0 || colDest >= colSpan) { continue; } - respMatrix[rowDest][colDest] += rspmat->getValue(irow, icol, flipRow, flipCol); + respMatrix[rowDest][colDest] += rspmat->getValue(irow, icol, ((innerBarrel && mSimRespIBOrientation) ? !flipRow : flipRow), flipCol); } } } From 546bb55b8c719f4b68c3164d78f762d618b6bda3 Mon Sep 17 00:00:00 2001 From: Felix Schlepper Date: Fri, 14 Feb 2025 10:03:52 +0100 Subject: [PATCH 06/10] ITS3: Add explanation to SegmentationMosaix Signed-off-by: Felix Schlepper --- .../include/ITS3Base/SegmentationMosaix.h | 20 +++++++++++++++++-- .../ITS3/base/include/ITS3Base/SpecsV2.h | 4 ++-- 2 files changed, 20 insertions(+), 4 deletions(-) diff --git a/Detectors/Upgrades/ITS3/base/include/ITS3Base/SegmentationMosaix.h b/Detectors/Upgrades/ITS3/base/include/ITS3Base/SegmentationMosaix.h index f4c61db5fc25f..2e3da358b4efe 100644 --- a/Detectors/Upgrades/ITS3/base/include/ITS3Base/SegmentationMosaix.h +++ b/Detectors/Upgrades/ITS3/base/include/ITS3Base/SegmentationMosaix.h @@ -32,6 +32,16 @@ class SegmentationMosaix // the more natural row,col layout: Also all the transformation between these // two. The class provides the transformation from the tile to TGeo // coordinates. + // In fact there exist three coordinate systems and one is transient. + // 1. The curved coordinate system. The chip's local coordinate system is + // defined with its center at the the mid-point of the tube. + // 2. The flat coordinate system. This is the tube segment projected onto a flat + // surface. In the projection we implicitly assume that the inner and outer + // stretch does not depend on the radius. + // Additionally, there is a difference between the flat geometrical center + // and the phyiscal center defined by the metal layer. + // 3. The detector coordinate system. Defined by the row and column segmentation + // defined at the upper edge in the flat coord. // row,col=0 // | @@ -71,7 +81,13 @@ class SegmentationMosaix static constexpr float SensorLayerThickness{constants::totalThickness}; static constexpr float NominalYShift{constants::nominalYShift}; - /// Transformation from the curved surface to a flat surface + /// Transformation from the curved surface to a flat surface. + /// Additionally a shift in the flat coordinates must be applied because + /// the center of the TGeoShap when projected will be higher than the + /// physical thickness of the chip (we add an additional hull to account for + /// the copper metal interconnection which is in reality part of the chip but in our + /// simulation the silicon and metal layer are separated). Thus we shift the projected center + /// down by this difference to align the coordinate systems. /// \param xCurved Detector local curved coordinate x in cm with respect to /// the center of the sensitive volume. /// \param yCurved Detector local curved coordinate y in cm with respect to @@ -93,7 +109,7 @@ class SegmentationMosaix } /// Transformation from the flat surface to a curved surface - /// It works only if the detector is not rototraslated + /// It works only if the detector is not rototraslated. /// \param xFlat Detector local flat coordinate x in cm with respect to /// the center of the sensitive volume. /// \param yFlat Detector local flat coordinate y in cm with respect to diff --git a/Detectors/Upgrades/ITS3/base/include/ITS3Base/SpecsV2.h b/Detectors/Upgrades/ITS3/base/include/ITS3Base/SpecsV2.h index 35c3db41ac77a..b604b8094363f 100644 --- a/Detectors/Upgrades/ITS3/base/include/ITS3Base/SpecsV2.h +++ b/Detectors/Upgrades/ITS3/base/include/ITS3Base/SpecsV2.h @@ -134,8 +134,8 @@ constexpr std::array radii{19.0006 * mm, 25.228 * mm, 31.4554 * constexpr std::array radiiInner{radii[0] - silicon::thicknessIn, radii[1] - silicon::thicknessIn, radii[2] - silicon::thicknessIn}; // inner silicon radius constexpr std::array radiiOuter{radii[0] + silicon::thicknessOut, radii[1] + silicon::thicknessOut, radii[2] + silicon::thicknessOut}; // outer silicon radius constexpr std::array radiiMiddle{(radiiInner[0] + radiiOuter[0]) / 2., (radiiInner[1] + radiiOuter[1]) / 2., (radiiInner[2] + radiiOuter[2]) / 2.}; // middle silicon radius -/*constexpr double nominalYShift{-metalstack::thickness / 2.}; // shift to position in silicion volume to the chip volume (silicon+metalstack)*/ -constexpr double nominalYShift{0}; // shift to position in silicion volume to the chip volume (silicon+metalstack) +constexpr double nominalYShift{-metalstack::thickness / 2.}; // shift to position in silicion volume to the chip volume (silicon+metalstack) +/*constexpr double nominalYShift{0}; // shift to position in silicion volume to the chip volume (silicon+metalstack)*/ // extra information of pixels and their response functions namespace pixelarray::pixels From 49a9c12925457b52984e911ed5f71ca0a3e578d5 Mon Sep 17 00:00:00 2001 From: Felix Schlepper Date: Thu, 27 Feb 2025 12:46:51 +0100 Subject: [PATCH 07/10] ITS3: split topo dict into IB and OB Signed-off-by: Felix Schlepper --- .../ITS3/base/include/ITS3Base/SpecsV2.h | 1 - .../ITS3/macros/test/CheckClusterSize.C | 5 +- .../ITS3/macros/test/CheckClustersITS3.C | 30 ++- .../macros/test/CompareClustersAndDigits.C | 4 +- .../ITS3/macros/test/CreateDictionariesITS3.C | 107 +++++--- .../BuildTopologyDictionary.h | 48 ++-- .../include/ITS3Reconstruction/Clusterer.h | 10 +- .../include/ITS3Reconstruction/IOUtils.h | 16 +- .../include/ITS3Reconstruction/LookUp.h | 20 +- .../ITS3Reconstruction/TopologyDictionary.h | 114 +++++--- .../src/BuildTopologyDictionary.cxx | 253 +++++++++++------- .../ITS3/reconstruction/src/Clusterer.cxx | 11 +- .../src/ITS3ReconstructionLinkDef.h | 1 + .../ITS3/reconstruction/src/LookUp.cxx | 21 +- .../reconstruction/src/TopologyDictionary.cxx | 164 ++++++++---- 15 files changed, 498 insertions(+), 307 deletions(-) diff --git a/Detectors/Upgrades/ITS3/base/include/ITS3Base/SpecsV2.h b/Detectors/Upgrades/ITS3/base/include/ITS3Base/SpecsV2.h index b604b8094363f..fedaad9182cce 100644 --- a/Detectors/Upgrades/ITS3/base/include/ITS3Base/SpecsV2.h +++ b/Detectors/Upgrades/ITS3/base/include/ITS3Base/SpecsV2.h @@ -135,7 +135,6 @@ constexpr std::array radiiInner{radii[0] - silicon::thicknessIn constexpr std::array radiiOuter{radii[0] + silicon::thicknessOut, radii[1] + silicon::thicknessOut, radii[2] + silicon::thicknessOut}; // outer silicon radius constexpr std::array radiiMiddle{(radiiInner[0] + radiiOuter[0]) / 2., (radiiInner[1] + radiiOuter[1]) / 2., (radiiInner[2] + radiiOuter[2]) / 2.}; // middle silicon radius constexpr double nominalYShift{-metalstack::thickness / 2.}; // shift to position in silicion volume to the chip volume (silicon+metalstack) -/*constexpr double nominalYShift{0}; // shift to position in silicion volume to the chip volume (silicon+metalstack)*/ // extra information of pixels and their response functions namespace pixelarray::pixels diff --git a/Detectors/Upgrades/ITS3/macros/test/CheckClusterSize.C b/Detectors/Upgrades/ITS3/macros/test/CheckClusterSize.C index 5fffed4fcc580..564b20350b883 100755 --- a/Detectors/Upgrades/ITS3/macros/test/CheckClusterSize.C +++ b/Detectors/Upgrades/ITS3/macros/test/CheckClusterSize.C @@ -255,13 +255,14 @@ void CheckClusterSize(std::string clusFileName = "o2clus_its.root", auto pattId = cluster.getPatternID(); auto id = cluster.getSensorID(); + auto ib = o2::its3::constants::detID::isDetITS3(id); int clusterSize{-1}; - if (pattId == o2::itsmft::CompCluster::InvalidPatternID || dict.isGroup(pattId)) { + if (pattId == o2::itsmft::CompCluster::InvalidPatternID || dict.isGroup(pattId, ib)) { o2::itsmft::ClusterPattern patt(pattIt); clusterSize = patt.getNPixels(); continue; } else { - clusterSize = dict.getNpixels(pattId); + clusterSize = dict.getNpixels(pattId, ib); } const auto& label = (clusLabArr->getLabels(clEntry))[0]; diff --git a/Detectors/Upgrades/ITS3/macros/test/CheckClustersITS3.C b/Detectors/Upgrades/ITS3/macros/test/CheckClustersITS3.C index a0f7dffc6ed1b..006271a1ea7bd 100644 --- a/Detectors/Upgrades/ITS3/macros/test/CheckClustersITS3.C +++ b/Detectors/Upgrades/ITS3/macros/test/CheckClustersITS3.C @@ -42,7 +42,7 @@ void CheckClustersITS3(const std::string& clusfile = "o2clus_its.root", const std::string& hitfile = "o2sim_HitsIT3.root", const std::string& inputGeom = "", - std::string dictfile = "../ccdb/IT3/Calib/ClusterDictionary/snapshot.root", + std::string dictfile = "./ccdb/IT3/Calib/ClusterDictionary/snapshot.root", bool batch = false) { gROOT->SetBatch(batch); @@ -63,7 +63,8 @@ void CheckClustersITS3(const std::string& clusfile = "o2clus_its.root", std::vector hitVecPool; std::vector mc2hitVec; - ULong_t cPattValid{0}, cPattInvalid{0}, cLabelInvalid{0}, cNoMC{0}; + ULong_t cPattValidIB{0}, cPattInvalidIB{0}, cLabelInvalidIB{0}, cNoMCIB{0}; + ULong_t cPattValidOB{0}, cPattInvalidOB{0}, cLabelInvalidOB{0}, cNoMCOB{0}; TFile fout("CheckClusters.root", "recreate"); TNtuple nt("ntc", "cluster ntuple", "ev:lab:hlx:hlz:hgx:hgz:tx:tz:cgx:cgy:cgz:clx:cly:clz:dx:dy:dz:ex:ez:patid:rof:npx:id:eta:row:col:lay"); @@ -103,6 +104,7 @@ void CheckClustersITS3(const std::string& clusfile = "o2clus_its.root", } else { LOG(info) << "Running without dictionary !"; } + dict.print(); // ROFrecords std::vector rofRecVec, *rofRecVecP = &rofRecVec; @@ -175,20 +177,18 @@ void CheckClustersITS3(const std::string& clusfile = "o2clus_its.root", auto isIB = o2::its3::constants::detID::isDetITS3(chipID); auto layer = o2::its3::constants::detID::getDetID2Layer(chipID); auto clusterSize{-1}; - if (pattID == o2::itsmft::CompCluster::InvalidPatternID || dict.isGroup(pattID)) { + if (pattID == o2::itsmft::CompCluster::InvalidPatternID || dict.isGroup(pattID, isIB)) { o2::itsmft::ClusterPattern patt(pattIt); locC = dict.getClusterCoordinates(cluster, patt, false); LOGP(debug, "I am invalid and I am on chip {}", chipID); - ++cPattInvalid; + (isIB) ? ++cPattInvalidIB : ++cPattInvalidOB; continue; } else { locC = dict.getClusterCoordinates(cluster); - errX = dict.getErrX(pattID); - errZ = dict.getErrZ(pattID); - errX *= (isIB) ? MosaixSegmentation::PitchRow : Segmentation::PitchRow; - errZ *= (isIB) ? MosaixSegmentation::PitchCol : Segmentation::PitchCol; - npix = dict.getNpixels(pattID); - ++cPattValid; + errX = dict.getErrX(pattID, isIB); + errZ = dict.getErrZ(pattID, isIB); + npix = dict.getNpixels(pattID, isIB); + (isIB) ? ++cPattValidIB : ++cPattValidOB; } // Transformation to the local --> global @@ -196,7 +196,7 @@ void CheckClustersITS3(const std::string& clusfile = "o2clus_its.root", const auto& lab = (clusLabArr->getLabels(clEntry))[0]; if (!lab.isValid()) { - ++cLabelInvalid; + (isIB) ? ++cLabelInvalidIB : ++cLabelInvalidOB; continue; } @@ -208,7 +208,7 @@ void CheckClustersITS3(const std::string& clusfile = "o2clus_its.root", auto hitEntry = mc2hit.find(key); if (hitEntry == mc2hit.end()) { LOG(debug) << "Failed to find MC hit entry for Tr" << trID << " chipID" << chipID; - ++cNoMC; + (isIB) ? ++cNoMCIB : ++cNoMCOB; continue; } const auto& hit = (*hitArray)[hitEntry->second]; @@ -263,8 +263,10 @@ void CheckClustersITS3(const std::string& clusfile = "o2clus_its.root", } } - LOGP(info, "There were {} valid PatternIDs and {} ({:.1f}%) invalid ones", cPattValid, cPattInvalid, ((float)cPattInvalid / (float)(cPattInvalid + cPattValid)) * 100); - LOGP(info, "There were {} invalid Labels and {} with No MC Hit information ", cLabelInvalid, cNoMC); + LOGP(info, "IB {} valid PatternIDs and {} ({:.1f}%) invalid ones", cPattValidIB, cPattInvalidIB, ((float)cPattInvalidIB / (float)(cPattInvalidIB + cPattValidIB)) * 100); + LOGP(info, "IB {} invalid Labels and {} with No MC Hit information ", cLabelInvalidIB, cNoMCIB); + LOGP(info, "OB {} valid PatternIDs and {} ({:.1f}%) invalid ones", cPattValidOB, cPattInvalidOB, ((float)cPattInvalidOB / (float)(cPattInvalidOB + cPattValidOB)) * 100); + LOGP(info, "OB {} invalid Labels and {} with No MC Hit information ", cLabelInvalidOB, cNoMCOB); auto canvCgXCgY = new TCanvas("canvCgXCgY", "", 1600, 1600); canvCgXCgY->Divide(2, 2); diff --git a/Detectors/Upgrades/ITS3/macros/test/CompareClustersAndDigits.C b/Detectors/Upgrades/ITS3/macros/test/CompareClustersAndDigits.C index 7770ede749f0c..c124481cc6f76 100644 --- a/Detectors/Upgrades/ITS3/macros/test/CompareClustersAndDigits.C +++ b/Detectors/Upgrades/ITS3/macros/test/CompareClustersAndDigits.C @@ -260,7 +260,7 @@ void CompareClustersAndDigits(std::string clusfile = "o2clus_it3.root", const auto pattID = cluster.getPatternID(); const auto isIB = o2::its3::constants::detID::isDetITS3(chipID); const auto layer = gman->getLayer(chipID); - if (pattID == o2::itsmft::CompCluster::InvalidPatternID || dict.isGroup(pattID)) { + if (pattID == o2::itsmft::CompCluster::InvalidPatternID || dict.isGroup(pattID, isIB)) { continue; } const auto& lab = (clusLabArr->getLabels(clEntry))[0]; @@ -316,7 +316,7 @@ void CompareClustersAndDigits(std::string clusfile = "o2clus_it3.root", data[chipID].cog->AddPoint(colC, rowC); constexpr float delta = 1e-2; - const auto& patt = dict.getPattern(cluster.getPatternID()); + const auto& patt = dict.getPattern(cluster.getPatternID(), isIB); auto box = new TBox( cluster.getCol() - delta - 0.5, cluster.getRow() - delta - 0.5, diff --git a/Detectors/Upgrades/ITS3/macros/test/CreateDictionariesITS3.C b/Detectors/Upgrades/ITS3/macros/test/CreateDictionariesITS3.C index 0116659940de1..cc241afb3357a 100644 --- a/Detectors/Upgrades/ITS3/macros/test/CreateDictionariesITS3.C +++ b/Detectors/Upgrades/ITS3/macros/test/CreateDictionariesITS3.C @@ -60,7 +60,7 @@ void CreateDictionariesITS3(bool saveDeltas = false, std::string collContextfile = "collisioncontext.root", std::string inputGeom = "", float checkOutliers = 2., // reject outliers (MC dX or dZ exceeds row/col span by a factor above the threshold) - float minPtMC = 0.01) // account only MC hits with pT above threshold + float minPtMC = 0.1) // account only MC hits with pT above threshold { const int QEDSourceID = 99; // Clusters from this MC source correspond to QED electrons @@ -84,10 +84,11 @@ void CreateDictionariesITS3(bool saveDeltas = false, std::array mMosaixSegmentations{0, 1, 2}; if (!clusDictFile.empty()) { clusDictOld.readFromFile(clusDictFile); - LOGP(info, "Loaded external cluster dictionary with {} entries from {}", clusDictOld.getSize(), clusDictFile); + LOGP(info, "Loaded external cluster dictionary with {} IB/{} OBentries from {}", clusDictOld.getSize(true), clusDictOld.getSize(false), clusDictFile); } - ULong_t cOk{0}, cOutliers{0}, cFailedMC{0}; + ULong_t cOkIB{0}, cOutliersIB{0}, cFailedMCIB{0}; + ULong_t cOkOB{0}, cOutliersOB{0}, cFailedMCOB{0}; TFile* fout = nullptr; TNtuple* nt = nullptr; @@ -233,17 +234,18 @@ void CreateDictionariesITS3(bool saveDeltas = false, const auto& cluster = (*clusArr)[clEntry]; o2::itsmft::ClusterPattern pattern; + bool ib = o2::its3::constants::detID::isDetITS3(cluster.getChipID()); if (cluster.getPatternID() != CompCluster::InvalidPatternID) { - if (clusDictOld.getSize() == 0) { + if (clusDictOld.getSize(ib) == 0) { LOG(error) << "Encountered patternID = " << cluster.getPatternID() << " != " << CompCluster::InvalidPatternID; LOG(error) << "Clusters have already been generated with a dictionary which was not provided"; return; } - if (clusDictOld.isGroup(cluster.getPatternID())) { + if (clusDictOld.isGroup(cluster.getPatternID(), ib)) { pattern.acquirePattern(pattIdx); } else { - pattern = clusDictOld.getPattern(cluster.getPatternID()); + pattern = clusDictOld.getPattern(cluster.getPatternID(), ib); } } else { pattern.acquirePattern(pattIdx); @@ -270,9 +272,8 @@ void CreateDictionariesITS3(bool saveDeltas = false, o2::math_utils::Vector3D xyzLocM; xyzLocM.SetCoordinates(0.5f * (xyzLocE.X() + xyzLocS.X()), 0.5f * (xyzLocE.Y() + xyzLocS.Y()), 0.5f * (xyzLocE.Z() + xyzLocS.Z())); auto locC = o2::its3::TopologyDictionary::getClusterCoordinates(cluster, pattern, false); - bool isIB = o2::its3::constants::detID::isDetITS3(chipID); int layer = gman->getLayer(chipID); - if (isIB) { + if (ib) { float xFlat{0.}, yFlat{0.}; mMosaixSegmentations[layer].curvedToFlat(xyzLocM.X(), xyzLocM.Y(), xFlat, yFlat); xyzLocM.SetCoordinates(xFlat, yFlat, xyzLocM.Z()); @@ -281,33 +282,33 @@ void CreateDictionariesITS3(bool saveDeltas = false, } dX = xyzLocM.X() - locC.X(); dZ = xyzLocM.Z() - locC.Z(); - dX /= (isIB) ? o2::its3::SegmentationMosaix::PitchRow : o2::itsmft::SegmentationAlpide::PitchRow; - dZ /= (isIB) ? o2::its3::SegmentationMosaix::PitchCol : o2::itsmft::SegmentationAlpide::PitchCol; + dX /= (ib) ? o2::its3::SegmentationMosaix::PitchRow : o2::itsmft::SegmentationAlpide::PitchRow; + dZ /= (ib) ? o2::its3::SegmentationMosaix::PitchCol : o2::itsmft::SegmentationAlpide::PitchCol; if (saveDeltas) { nt->Fill(topology.getHash(), dX, dZ); } if (checkOutliers > 0.) { if (bool bX = std::abs(dX) > topology.getRowSpan() * checkOutliers, bZ = std::abs(dZ) > topology.getColumnSpan() * checkOutliers; bX || bZ) { // ignore outlier - ++cOutliers; + (ib) ? ++cOutliersIB : ++cOutliersOB; LOGP(debug, "Ignored Value dX={} > {} * {} -> {}", dX, topology.getRowSpan(), checkOutliers, bX); LOGP(debug, "Ignored Value dZ={} > {} * {} -> {}", dZ, topology.getColumnSpan(), checkOutliers, bZ); dX = dZ = BuildTopologyDictionary::IgnoreVal; } else { - ++cOk; + (ib) ? ++cOkIB : ++cOkOB; } } } } else { /* LOGP(info, " Failed to find MC hit entry for Tr: {} chipID: {}", trID, chipID); */ /* lab.print(); */ - ++cFailedMC; + (ib) ? ++cFailedMCIB : ++cFailedMCOB; } - signalDictionary.accountTopology(topology, dX, dZ); + signalDictionary.accountTopology(topology, ib, dX, dZ); } else { - noiseDictionary.accountTopology(topology, dX, dZ); + noiseDictionary.accountTopology(topology, ib, dX, dZ); } } - completeDictionary.accountTopology(topology, dX, dZ); + completeDictionary.accountTopology(topology, ib, dX, dZ); } // clean MC cache for events which are not needed anymore @@ -323,12 +324,14 @@ void CreateDictionariesITS3(bool saveDeltas = false, } } - LOGP(info, "Clusters: {} okay (failed MCHit2Clus {}); outliers {}", cOk, cFailedMC, cOutliers); + LOGP(info, "IB Clusters: {} okay (failed MCHit2Clus {}); outliers {}", cOkIB, cFailedMCIB, cOutliersIB); + LOGP(info, "OB Clusters: {} okay (failed MCHit2Clus {}); outliers {}", cOkOB, cFailedMCOB, cOutliersOB); auto dID = o2::detectors::DetID::IT3; LOGP(info, "Complete Dictionary:"); - completeDictionary.setThreshold(probThreshold); + completeDictionary.setThreshold(probThreshold, true); + completeDictionary.setThreshold(probThreshold, false); completeDictionary.groupRareTopologies(); completeDictionary.printDictionaryBinary(o2::base::DetectorNameConf::getAlpideClusterDictionaryFileName(dID, "")); completeDictionary.printDictionary(o2::base::DetectorNameConf::getAlpideClusterDictionaryFileName(dID, "", "txt")); @@ -336,24 +339,34 @@ void CreateDictionariesITS3(bool saveDeltas = false, TFile histogramOutput("histograms.root", "recreate"); TCanvas* cComplete = new TCanvas("cComplete", "Distribution of all the topologies"); - cComplete->cd(); - cComplete->SetLogy(); - TH1F* hComplete = completeDictionary.getDictionary().getTopologyDistribution("hComplete"); - hComplete->SetDirectory(nullptr); - hComplete->Draw("hist"); - hComplete->Write(); + cComplete->Divide(2, 1); + cComplete->cd(1); + TH1F* hCompleteIB = completeDictionary.getDictionary().getTopologyDistribution("hCompleteInnerBarrel", true); + hCompleteIB->SetDirectory(nullptr); + hCompleteIB->Draw("hist"); + gPad->SetLogy(); + cComplete->cd(2); + TH1F* hCompleteOB = completeDictionary.getDictionary().getTopologyDistribution("hCompleteOuterBarrel", false); + hCompleteOB->SetDirectory(nullptr); + hCompleteOB->Draw("hist"); + gPad->SetLogy(); + histogramOutput.cd(); + hCompleteIB->Write(); + hCompleteOB->Write(); cComplete->Write(); if (clusLabArr) { LOGP(info, "Noise Dictionary:"); - noiseDictionary.setThreshold(0.0001); + noiseDictionary.setThreshold(0.0001, true); + noiseDictionary.setThreshold(0.0001, false); noiseDictionary.groupRareTopologies(); noiseDictionary.printDictionaryBinary(o2::base::DetectorNameConf::getAlpideClusterDictionaryFileName(dID, "noiseClusTopo")); noiseDictionary.printDictionary(o2::base::DetectorNameConf::getAlpideClusterDictionaryFileName(dID, "noiseClusTopo", "txt")); noiseDictionary.saveDictionaryRoot(o2::base::DetectorNameConf::getAlpideClusterDictionaryFileName(dID, "noiseClusTopo", "root")); LOGP(info, "Signal Dictionary:"); - signalDictionary.setThreshold(0.0001); + signalDictionary.setThreshold(0.0001, true); + signalDictionary.setThreshold(0.0001, false); signalDictionary.groupRareTopologies(); signalDictionary.printDictionaryBinary(o2::base::DetectorNameConf::getAlpideClusterDictionaryFileName(dID, "signal")); signalDictionary.printDictionary(o2::base::DetectorNameConf::getAlpideClusterDictionaryFileName(dID, "signal", "txt")); @@ -361,26 +374,42 @@ void CreateDictionariesITS3(bool saveDeltas = false, LOGP(info, "Plotting Channels"); auto cNoise = new TCanvas("cNoise", "Distribution of noise topologies"); - cNoise->cd(); - cNoise->SetLogy(); - auto hNoise = noiseDictionary.getDictionary().getTopologyDistribution("hNoise"); - hNoise->SetDirectory(nullptr); - hNoise->Draw("hist"); + cNoise->Divide(2, 1); + cNoise->cd(1); + auto hNoiseIB = noiseDictionary.getDictionary().getTopologyDistribution("hNoiseInnerBarrel", true); + hNoiseIB->SetDirectory(nullptr); + hNoiseIB->Draw("hist"); + gPad->SetLogy(); + cNoise->cd(2); + auto hNoiseOB = noiseDictionary.getDictionary().getTopologyDistribution("hNoiseOuterBarrel", false); + hNoiseOB->SetDirectory(nullptr); + hNoiseOB->Draw("hist"); + gPad->SetLogy(); histogramOutput.cd(); - hNoise->Write(); + hNoiseIB->Write(); + hNoiseOB->Write(); cNoise->Write(); + auto cSignal = new TCanvas("cSignal", "cSignal"); - cSignal->cd(); + cSignal->Divide(2, 1); + cSignal->cd(1); + auto hSignalIB = signalDictionary.getDictionary().getTopologyDistribution("hSignalInnerBarrel", true); + hSignalIB->SetDirectory(nullptr); + hSignalIB->Draw("hist"); + gPad->SetLogy(); + cSignal->cd(2); cSignal->SetLogy(); - auto hSignal = signalDictionary.getDictionary().getTopologyDistribution("hSignal"); - hSignal->SetDirectory(nullptr); - hSignal->Draw("hist"); + auto hSignalOB = signalDictionary.getDictionary().getTopologyDistribution("hSignalOuterBarrel", false); + hSignalOB->SetDirectory(nullptr); + hSignalOB->Draw("hist"); + gPad->SetLogy(); histogramOutput.cd(); - hSignal->Write(); + hSignalIB->Write(); + hSignalOB->Write(); cSignal->Write(); - sw.Stop(); - sw.Print(); } + sw.Stop(); + sw.Print(); if (saveDeltas) { fout->cd(); nt->Write(); diff --git a/Detectors/Upgrades/ITS3/reconstruction/include/ITS3Reconstruction/BuildTopologyDictionary.h b/Detectors/Upgrades/ITS3/reconstruction/include/ITS3Reconstruction/BuildTopologyDictionary.h index 7df603bb29fb2..662c58aeb2cd8 100644 --- a/Detectors/Upgrades/ITS3/reconstruction/include/ITS3Reconstruction/BuildTopologyDictionary.h +++ b/Detectors/Upgrades/ITS3/reconstruction/include/ITS3Reconstruction/BuildTopologyDictionary.h @@ -24,31 +24,47 @@ namespace o2::its3 class BuildTopologyDictionary { + using TopoInfo = std::unordered_map; + using TopoStat = std::map; + using TopoFreq = std::vector>; + public: static constexpr float IgnoreVal = 999.; - void accountTopology(const itsmft::ClusterTopology& cluster, float dX = IgnoreVal, float dZ = IgnoreVal); - void setNCommon(unsigned int nCommon); // set number of common topologies - void setThreshold(double thr); - void setThresholdCumulative(double cumulative); // Considering the integral + void accountTopology(const itsmft::ClusterTopology& cluster, bool IB, float dX = IgnoreVal, float dZ = IgnoreVal); + void setNCommon(unsigned int nCommon, bool IB); // set number of common topologies + void setThreshold(double thr, bool IB); + void setThresholdCumulative(double cumulative, bool IB); // Considering the integral void groupRareTopologies(); - friend std::ostream& operator<<(std::ostream& os, const BuildTopologyDictionary& BD); void printDictionary(const std::string& fname); void printDictionaryBinary(const std::string& fname); void saveDictionaryRoot(const std::string& fname); - unsigned int getTotClusters() const { return mTotClusters; } - unsigned int getNotInGroups() const { return mNCommonTopologies; } - TopologyDictionary getDictionary() const { return mDictionary; } + [[nodiscard]] unsigned int getTotClusters(bool IB) const { return (IB) ? mTotClustersIB : mTotClustersOB; } + [[nodiscard]] unsigned int getNotInGroups(bool IB) const { return (IB) ? mNCommonTopologiesIB : mNCommonTopologiesOB; } + [[nodiscard]] const TopologyDictionary& getDictionary() const { return mDictionary; } + + friend std::ostream& operator<<(std::ostream& os, const BuildTopologyDictionary& BD); private: - TopologyDictionary mDictionary; ///< Dictionary of topologies - std::map mTopologyMap; //! Temporary map of type - std::vector> mTopologyFrequency; //! , needed to define threshold - unsigned int mTotClusters{0}; - unsigned int mNCommonTopologies{0}; - double mFrequencyThreshold{0.}; - - std::unordered_map mMapInfo; + void accountTopologyImpl(const itsmft::ClusterTopology& cluster, TopoInfo& tinfo, TopoStat& tstat, unsigned int& ntot, float sigmaX, float sigmaZ, float dX, float dZ); + void setNCommonImpl(unsigned int ncom, TopoFreq& tfreq, TopoStat& tstat, unsigned int& ncommon, unsigned int ntot); + void setThresholdImpl(double thr, TopoFreq& tfreq, TopoInfo& tinfo, TopoStat& tstat, unsigned int& ncommon, double& freqthres, unsigned int ntot); + void setThresholdCumulativeImpl(double cumulative, TopoFreq& tfreq, unsigned int& ncommon, double& freqthres, unsigned int ntot); + void groupRareTopologiesImpl(TopoFreq& tfreq, TopoInfo& tinfo, TopoStat& tstat, unsigned int& ncommon, double& freqthres, TopologyDictionaryData& data, unsigned int ntot); + + TopologyDictionary mDictionary; ///< Dictionary of topologies + unsigned int mTotClustersIB{0}; + unsigned int mTotClustersOB{0}; + unsigned int mNCommonTopologiesIB{0}; + unsigned int mNCommonTopologiesOB{0}; + double mFrequencyThresholdIB{0.}; + double mFrequencyThresholdOB{0.}; + TopoInfo mMapInfoIB; + TopoInfo mMapInfoOB; + TopoStat mTopologyMapIB; //! IB Temporary map of type + TopoStat mTopologyMapOB; //! OB Temporary map of type + TopoFreq mTopologyFrequencyIB; //! IB , needed to define threshold + TopoFreq mTopologyFrequencyOB; //! OB , needed to define threshold ClassDefNV(BuildTopologyDictionary, 3); }; diff --git a/Detectors/Upgrades/ITS3/reconstruction/include/ITS3Reconstruction/Clusterer.h b/Detectors/Upgrades/ITS3/reconstruction/include/ITS3Reconstruction/Clusterer.h index 20acf07d4f547..a81db09217e9b 100644 --- a/Detectors/Upgrades/ITS3/reconstruction/include/ITS3Reconstruction/Clusterer.h +++ b/Detectors/Upgrades/ITS3/reconstruction/include/ITS3Reconstruction/Clusterer.h @@ -207,7 +207,7 @@ class Clusterer template static void streamCluster(const std::vector& pixbuf, const std::array* lblBuff, const BBox& bbox, const its3::LookUp& pattIdConverter, - VCLUS* compClusPtr, VPAT* patternsPtr, MCTruth* labelsClusPtr, int nlab, bool isHuge = false); + VCLUS* compClusPtr, VPAT* patternsPtr, MCTruth* labelsClusPtr, int nlab, bool isIB, bool isHuge = false); bool isContinuousReadOut() const { return mContinuousReadout; } void setContinuousReadOut(bool v) { mContinuousReadout = v; } @@ -230,7 +230,7 @@ class Clusterer ///< load the dictionary of cluster topologies void setDictionary(const its3::TopologyDictionary* dict) { - LOGP(info, "Setting TopologyDictionary with size={}", dict->getSize()); + LOGP(info, "Setting TopologyDictionary with IB size={} & OB size={}", dict->getSize(true), dict->getSize(false)); mPattIdConverter.setDictionary(dict); // dict->print(); } @@ -274,7 +274,7 @@ class Clusterer template void Clusterer::streamCluster(const std::vector& pixbuf, const std::array* lblBuff, const Clusterer::BBox& bbox, const its3::LookUp& pattIdConverter, - VCLUS* compClusPtr, VPAT* patternsPtr, MCTruth* labelsClusPtr, int nlab, bool isHuge) + VCLUS* compClusPtr, VPAT* patternsPtr, MCTruth* labelsClusPtr, int nlab, bool isIB, bool isHuge) { if (labelsClusPtr && lblBuff) { // MC labels were requested auto cnt = compClusPtr->size(); @@ -291,10 +291,10 @@ void Clusterer::streamCluster(const std::vector& pixbuf, const std::a int nbits = ir * colSpanW + ic; patt[nbits >> 3] |= (0x1 << (7 - (nbits % 8))); } - uint16_t pattID = (isHuge || pattIdConverter.size() == 0) ? CompCluster::InvalidPatternID : pattIdConverter.findGroupID(rowSpanW, colSpanW, patt.data()); + uint16_t pattID = (isHuge || pattIdConverter.size(isIB) == 0) ? CompCluster::InvalidPatternID : pattIdConverter.findGroupID(rowSpanW, colSpanW, isIB, patt.data()); uint16_t row = bbox.rowMin, col = bbox.colMin; LOGP(debug, "PattID: findGroupID({},{},{})={}", row, col, patt[0], pattID); - if (pattID == CompCluster::InvalidPatternID || pattIdConverter.isGroup(pattID)) { + if (pattID == CompCluster::InvalidPatternID || pattIdConverter.isGroup(pattID, isIB)) { if (pattID != CompCluster::InvalidPatternID) { // For groupped topologies, the reference pixel is the COG pixel float xCOG = 0., zCOG = 0.; diff --git a/Detectors/Upgrades/ITS3/reconstruction/include/ITS3Reconstruction/IOUtils.h b/Detectors/Upgrades/ITS3/reconstruction/include/ITS3Reconstruction/IOUtils.h index 4f865d11c90a9..b9e7fd0f6ec39 100644 --- a/Detectors/Upgrades/ITS3/reconstruction/include/ITS3Reconstruction/IOUtils.h +++ b/Detectors/Upgrades/ITS3/reconstruction/include/ITS3Reconstruction/IOUtils.h @@ -30,13 +30,14 @@ template o2::math_utils::Point3D extractClusterData(const itsmft::CompClusterExt& c, iterator& iter, const its3::TopologyDictionary* dict, T& sig2y, T& sig2z) { auto pattID = c.getPatternID(); + auto ib = constants::detID::isDetITS3(c.getSensorID()); // Dummy COG errors (about half pixel size) - sig2y = (constants::detID::isDetITS3(c.getSensorID())) ? DefClusError2Row : o2::its::ioutils::DefClusError2Row; - sig2z = (constants::detID::isDetITS3(c.getSensorID())) ? DefClusError2Col : o2::its::ioutils::DefClusError2Col; + sig2y = (ib) ? DefClusError2Row : o2::its::ioutils::DefClusError2Row; + sig2z = (ib) ? DefClusError2Col : o2::its::ioutils::DefClusError2Col; if (pattID != itsmft::CompCluster::InvalidPatternID) { - sig2y = dict->getErr2X(pattID) * sig2y; // Error is given in detector coordinates - sig2z = dict->getErr2Z(pattID) * sig2z; - if (!dict->isGroup(pattID)) { + sig2y = dict->getErr2X(pattID, ib); + sig2z = dict->getErr2Z(pattID, ib); + if (!dict->isGroup(pattID, ib)) { return dict->getClusterCoordinates(c); } else { o2::itsmft::ClusterPattern patt(iter); @@ -52,13 +53,14 @@ template o2::math_utils::Point3D extractClusterData(const itsmft::CompClusterExt& c, iterator& iter, const its3::TopologyDictionary* dict, T& sig2y, T& sig2z, uint8_t& cls) { auto pattID = c.getPatternID(); + auto ib = constants::detID::isDetITS3(c.getSensorID()); auto iterC = iter; unsigned int clusterSize{999}; - if (pattID == itsmft::CompCluster::InvalidPatternID || dict->isGroup(pattID)) { + if (pattID == itsmft::CompCluster::InvalidPatternID || dict->isGroup(pattID, ib)) { o2::itsmft::ClusterPattern patt(iterC); clusterSize = patt.getNPixels(); } else { - clusterSize = dict->getNpixels(pattID); + clusterSize = dict->getNpixels(pattID, ib); } cls = static_cast(std::clamp(clusterSize, static_cast(std::numeric_limits::min()), static_cast(std::numeric_limits::max()))); return extractClusterData(c, iter, dict, sig2y, sig2z); diff --git a/Detectors/Upgrades/ITS3/reconstruction/include/ITS3Reconstruction/LookUp.h b/Detectors/Upgrades/ITS3/reconstruction/include/ITS3Reconstruction/LookUp.h index 0fbecb41393ff..809a129a0debf 100644 --- a/Detectors/Upgrades/ITS3/reconstruction/include/ITS3Reconstruction/LookUp.h +++ b/Detectors/Upgrades/ITS3/reconstruction/include/ITS3Reconstruction/LookUp.h @@ -21,7 +21,6 @@ #ifndef ALICEO2_ITS3_LOOKUP_H #define ALICEO2_ITS3_LOOKUP_H -#include "DataFormatsITSMFT/ClusterTopology.h" #include "ITS3Reconstruction/TopologyDictionary.h" namespace o2::its3 @@ -32,20 +31,21 @@ class LookUp LookUp() = default; LookUp(std::string fileName); static int groupFinder(int nRow, int nCol); - int findGroupID(int nRow, int nCol, const unsigned char patt[itsmft::ClusterPattern::MaxPatternBytes]) const; - int getTopologiesOverThreshold() const { return mTopologiesOverThreshold; } + int findGroupID(int nRow, int nCol, bool IB, const unsigned char patt[itsmft::ClusterPattern::MaxPatternBytes]) const; + int getTopologiesOverThreshold(bool IB) const { return (IB) ? mTopologiesOverThresholdIB : mTopologiesOverThresholdOB; } void loadDictionary(std::string fileName); void setDictionary(const TopologyDictionary* dict); - bool isGroup(int id) const { return mDictionary.isGroup(id); } - int size() const { return mDictionary.getSize(); } - auto getPattern(int id) const { return mDictionary.getPattern(id); } - auto getDictionaty() const { return mDictionary; } + auto getDictionary() const { return mDictionary; } + bool isGroup(int id, bool IB) const { return mDictionary.isGroup(id, IB); } + int size(bool IB) const { return mDictionary.getSize(IB); } + auto getPattern(int id, bool IB) const { return mDictionary.getPattern(id, IB); } private: - TopologyDictionary mDictionary{}; - int mTopologiesOverThreshold{0}; + TopologyDictionary mDictionary; + int mTopologiesOverThresholdIB{0}; + int mTopologiesOverThresholdOB{0}; - ClassDefNV(LookUp, 2); + ClassDefNV(LookUp, 3); }; } // namespace o2::its3 diff --git a/Detectors/Upgrades/ITS3/reconstruction/include/ITS3Reconstruction/TopologyDictionary.h b/Detectors/Upgrades/ITS3/reconstruction/include/ITS3Reconstruction/TopologyDictionary.h index 1a9b4048fd06f..d5f5721170aa7 100644 --- a/Detectors/Upgrades/ITS3/reconstruction/include/ITS3Reconstruction/TopologyDictionary.h +++ b/Detectors/Upgrades/ITS3/reconstruction/include/ITS3Reconstruction/TopologyDictionary.h @@ -24,6 +24,18 @@ namespace o2::its3 class BuildTopologyDictionary; class LookUp; +struct TopologyDictionaryData { + static constexpr int STopoSize{(8 * 255) + 1}; + std::array mSmallTopologiesLUT{}; ///< Look-Up Table for the topologies with 1-byte linearised matrix + std::vector mVectorOfIDs; ///< Vector of topologies and groups + std::unordered_map mCommonMap; ///< Map of pair + std::unordered_map mGroupMap; ///< Map of pair + + void print() const noexcept; + + ClassDefNV(TopologyDictionaryData, 1); +}; + class TopologyDictionary { public: @@ -37,86 +49,103 @@ class TopologyDictionary static constexpr int MaxNumberOfRowClasses = 1 + ((itsmft::ClusterPattern::MaxRowSpan - 1) / RowClassSpan); ///< Maximum number of row classes for the groups of rare topologies static constexpr int MaxNumberOfColClasses = 1 + ((itsmft::ClusterPattern::MaxColSpan - 1) / ColClassSpan); ///< Maximum number of col classes for the groups of rare topologies static constexpr int NumberOfRareGroups = MaxNumberOfRowClasses * MaxNumberOfColClasses; ///< Number of entries corresponding to groups of rare topologies (those whos matrix exceed the max number of bytes are empty). + /// Resets internal structures + void reset() noexcept; + void resetMaps(bool IB = true) noexcept; /// Prints the dictionary friend std::ostream& operator<<(std::ostream& os, const its3::TopologyDictionary& dictionary); /// Prints the dictionary in a binary file void writeBinaryFile(const std::string& outputFile); /// Reads the dictionary from a binary file - int readBinaryFile(const std::string& fileName); - - int readFromFile(const std::string& fileName); + void readBinaryFile(const std::string& fileName); + void readFromFile(const std::string& fileName); + void print() const noexcept; /// Returns the x position of the COG for the n_th element - inline float getXCOG(int n) const + [[nodiscard]] float getXCOG(int n, bool IB = true) const { - assert(n >= 0 || n < (int)mVectorOfIDs.size()); - return mVectorOfIDs[n].mXCOG; + const auto& data = (IB) ? mDataIB : mDataOB; + assert(n >= 0 || n < (int)data.mVectorOfIDs.size()); + return data.mVectorOfIDs[n].mXCOG; } /// Returns the error on the x position of the COG for the n_th element - inline float getErrX(int n) const + [[nodiscard]] float getErrX(int n, bool IB = true) const { - assert(n >= 0 || n < (int)mVectorOfIDs.size()); - return mVectorOfIDs[n].mErrX; + const auto& data = (IB) ? mDataIB : mDataOB; + assert(n >= 0 || n < (int)data.mVectorOfIDs.size()); + return data.mVectorOfIDs[n].mErrX; } /// Returns the z position of the COG for the n_th element - inline float getZCOG(int n) const + [[nodiscard]] float getZCOG(int n, bool IB = true) const { - assert(n >= 0 || n < (int)mVectorOfIDs.size()); - return mVectorOfIDs[n].mZCOG; + const auto& data = (IB) ? mDataIB : mDataOB; + assert(n >= 0 || n < (int)data.mVectorOfIDs.size()); + return data.mVectorOfIDs[n].mZCOG; } /// Returns the error on the z position of the COG for the n_th element - inline float getErrZ(int n) const + [[nodiscard]] float getErrZ(int n, bool IB = true) const { - assert(n >= 0 || n < (int)mVectorOfIDs.size()); - return mVectorOfIDs[n].mErrZ; + const auto& data = (IB) ? mDataIB : mDataOB; + assert(n >= 0 || n < (int)data.mVectorOfIDs.size()); + return data.mVectorOfIDs[n].mErrZ; } /// Returns the error^2 on the x position of the COG for the n_th element - inline float getErr2X(int n) const + [[nodiscard]] float getErr2X(int n, bool IB = true) const { - assert(n >= 0 || n < (int)mVectorOfIDs.size()); - return mVectorOfIDs[n].mErr2X; + const auto& data = (IB) ? mDataIB : mDataOB; + assert(n >= 0 || n < (int)data.mVectorOfIDs.size()); + return data.mVectorOfIDs[n].mErr2X; } /// Returns the error^2 on the z position of the COG for the n_th element - inline float getErr2Z(int n) const + [[nodiscard]] float getErr2Z(int n, bool IB = true) const { - assert(n >= 0 || n < (int)mVectorOfIDs.size()); - return mVectorOfIDs[n].mErr2Z; + const auto& data = (IB) ? mDataIB : mDataOB; + assert(n >= 0 || n < (int)data.mVectorOfIDs.size()); + return data.mVectorOfIDs[n].mErr2Z; } /// Returns the hash of the n_th element - inline unsigned long getHash(int n) const + [[nodiscard]] unsigned long getHash(int n, bool IB = true) const { - assert(n >= 0 || n < (int)mVectorOfIDs.size()); - return mVectorOfIDs[n].mHash; + const auto& data = (IB) ? mDataIB : mDataOB; + assert(n >= 0 || n < (int)data.mVectorOfIDs.size()); + return data.mVectorOfIDs[n].mHash; } /// Returns the number of fired pixels of the n_th element - inline int getNpixels(int n) const + [[nodiscard]] int getNpixels(int n, bool IB = true) const { - assert(n >= 0 || n < (int)mVectorOfIDs.size()); - return mVectorOfIDs[n].mNpixels; + const auto& data = (IB) ? mDataIB : mDataOB; + assert(n >= 0 || n < (int)data.mVectorOfIDs.size()); + return data.mVectorOfIDs[n].mNpixels; } /// Returns the frequency of the n_th element; - inline double getFrequency(int n) const + [[nodiscard]] double getFrequency(int n, bool IB = true) const { - assert(n >= 0 || n < (int)mVectorOfIDs.size()); - return mVectorOfIDs[n].mFrequency; + const auto& data = (IB) ? mDataIB : mDataOB; + assert(n >= 0 || n < (int)data.mVectorOfIDs.size()); + return data.mVectorOfIDs[n].mFrequency; } /// Returns true if the element corresponds to a group of rare topologies - inline bool isGroup(int n) const + [[nodiscard]] bool isGroup(int n, bool IB = true) const { - assert(n >= 0 || n < (int)mVectorOfIDs.size()); - return mVectorOfIDs[n].mIsGroup; + const auto& data = (IB) ? mDataIB : mDataOB; + assert(n >= 0 || n < (int)data.mVectorOfIDs.size()); + return data.mVectorOfIDs[n].mIsGroup; } /// Returns the pattern of the topology - inline const itsmft::ClusterPattern& getPattern(int n) const + [[nodiscard]] const itsmft::ClusterPattern& getPattern(int n, bool IB = true) const { - assert(n >= 0 || n < (int)mVectorOfIDs.size()); - return mVectorOfIDs[n].mPattern; + const auto& data = (IB) ? mDataIB : mDataOB; + assert(n >= 0 || n < (int)data.mVectorOfIDs.size()); + return data.mVectorOfIDs[n].mPattern; } /// Fills a hostogram with the distribution of the IDs - TH1F* getTopologyDistribution(const std::string_view hname = "h_topo_dist") const; + [[nodiscard]] TH1F* getTopologyDistribution(const std::string_view hname, bool IB = true) const; /// Returns the number of elements in the dicionary; - int getSize() const { return (int)mVectorOfIDs.size(); } + [[nodiscard]] int getSize(bool IB) const + { + return static_cast((IB) ? mDataIB.mVectorOfIDs.size() : mDataOB.mVectorOfIDs.size()); + } /// Returns the local position of a compact cluster /// Returns the local position of a compact cluster @@ -133,13 +162,10 @@ class TopologyDictionary friend its3::LookUp; private: - static constexpr int STopoSize{8 * 255 + 1}; - std::unordered_map mCommonMap{}; ///< Map of pair - std::unordered_map mGroupMap{}; ///< Map of pair - int mSmallTopologiesLUT[STopoSize]{}; ///< Look-Up Table for the topologies with 1-byte linearised matrix - std::vector mVectorOfIDs{}; ///< Vector of topologies and groups + TopologyDictionaryData mDataIB; + TopologyDictionaryData mDataOB; - ClassDefNV(TopologyDictionary, 3); + ClassDefNV(TopologyDictionary, 4); }; } // namespace o2::its3 diff --git a/Detectors/Upgrades/ITS3/reconstruction/src/BuildTopologyDictionary.cxx b/Detectors/Upgrades/ITS3/reconstruction/src/BuildTopologyDictionary.cxx index 8667066f91fac..f7eec52f9434a 100644 --- a/Detectors/Upgrades/ITS3/reconstruction/src/BuildTopologyDictionary.cxx +++ b/Detectors/Upgrades/ITS3/reconstruction/src/BuildTopologyDictionary.cxx @@ -15,20 +15,34 @@ #include "ITS3Reconstruction/LookUp.h" #include "DataFormatsITSMFT/CompCluster.h" +#include "ITSMFTBase/SegmentationAlpide.h" +#include "ITS3Base/SegmentationMosaix.h" + #include "TFile.h" ClassImp(o2::its3::BuildTopologyDictionary); namespace o2::its3 { -void BuildTopologyDictionary::accountTopology(const itsmft::ClusterTopology& cluster, float dX, float dZ) +void BuildTopologyDictionary::accountTopology(const itsmft::ClusterTopology& cluster, bool IB, float dX, float dZ) { - mTotClusters++; + accountTopologyImpl(cluster, + ((IB) ? mMapInfoIB : mMapInfoOB), + ((IB) ? mTopologyMapIB : mTopologyMapOB), + ((IB) ? mTotClustersIB : mTotClustersOB), + ((IB) ? SegmentationMosaix::PitchRow : itsmft::SegmentationAlpide::PitchRow), + ((IB) ? SegmentationMosaix::PitchCol : itsmft::SegmentationAlpide::PitchCol), + dX, dZ); +} + +void BuildTopologyDictionary::accountTopologyImpl(const itsmft::ClusterTopology& cluster, TopoInfo& tinfo, TopoStat& tstat, unsigned int& tot, float sigmaX, float sigmaZ, float dX, float dZ) +{ + ++tot; bool useDf = dX < IgnoreVal / 2; // we may need to account the frequency but to not update the centroid // std::pair::iterator,bool> ret; // auto ret = mTopologyMap.insert(std::make_pair(cluster.getHash(), std::make_pair(cluster, 1))); - auto& topoStat = mTopologyMap[cluster.getHash()]; + auto& topoStat = tstat[cluster.getHash()]; topoStat.countsTotal++; if (topoStat.countsTotal == 1) { // a new topology is inserted topoStat.topology = cluster; @@ -44,14 +58,14 @@ void BuildTopologyDictionary::accountTopology(const itsmft::ClusterTopology& clu topInf.mZmean = dZ; topoStat.countsWithBias = 1; } else { // assign expected sigmas from the pixel X, Z sizes - topInf.mXsigma2 = 1.f / 12.f / (float)std::min(10, topInf.mSizeX); - topInf.mZsigma2 = 1.f / 12.f / (float)std::min(10, topInf.mSizeZ); + topInf.mXsigma2 = sigmaX * sigmaX / 12.f / (float)std::min(10, topInf.mSizeX); + topInf.mZsigma2 = sigmaZ * sigmaZ / (float)std::min(10, topInf.mSizeZ); } - mMapInfo.emplace(cluster.getHash(), topInf); + tinfo.emplace(cluster.getHash(), topInf); } else { if (useDf) { auto num = topoStat.countsWithBias++; - auto ind = mMapInfo.find(cluster.getHash()); + auto ind = tinfo.find(cluster.getHash()); float tmpxMean = ind->second.mXmean; float newxMean = ind->second.mXmean = ((tmpxMean)*num + dX) / (num + 1); float tmpxSigma2 = ind->second.mXsigma2; @@ -64,101 +78,135 @@ void BuildTopologyDictionary::accountTopology(const itsmft::ClusterTopology& clu } } -void BuildTopologyDictionary::setThreshold(double thr) +void BuildTopologyDictionary::setNCommon(unsigned int nCommon, bool IB) +{ + mDictionary.resetMaps(IB); + + auto& freqTopo = ((IB) ? mTopologyFrequencyIB : mTopologyFrequencyOB); + auto& freqThres = ((IB) ? mFrequencyThresholdIB : mFrequencyThresholdOB); + auto& comTopo = ((IB) ? mNCommonTopologiesIB : mNCommonTopologiesOB); + auto ntot = ((IB) ? mTotClustersIB : mTotClustersOB); + + setNCommonImpl(nCommon, + freqTopo, + ((IB) ? mTopologyMapIB : mTopologyMapOB), + comTopo, + ntot); + // Recaculate also the threshold + freqThres = ((double)freqTopo[comTopo - 1].first) / ntot; +} + +void BuildTopologyDictionary::setNCommonImpl(unsigned int ncom, TopoFreq& tfreq, TopoStat& tstat, unsigned int& ncommon, unsigned int ntot) { - mTopologyFrequency.clear(); - for (auto&& p : mTopologyMap) { // p is pair - mTopologyFrequency.emplace_back(p.second.countsTotal, p.first); + if (ncom >= itsmft::CompCluster::InvalidPatternID) { + LOGP(warning, "Redefining nCommon from {} to {} to be below InvalidPatternID", ncom, itsmft::CompCluster::InvalidPatternID - 1); + ncom = itsmft::CompCluster::InvalidPatternID - 1; + } + tfreq.clear(); + for (auto&& p : tstat) { // p os pair + tfreq.emplace_back(p.second.countsTotal, p.first); } - std::sort(mTopologyFrequency.begin(), mTopologyFrequency.end(), + std::sort(tfreq.begin(), tfreq.end(), [](const std::pair& couple1, const std::pair& couple2) { return (couple1.first > couple2.first); }); - mNCommonTopologies = 0; - mDictionary.mCommonMap.clear(); - mDictionary.mGroupMap.clear(); - mFrequencyThreshold = thr; - for (auto& q : mTopologyFrequency) { - if (((double)q.first) / mTotClusters > thr) { - mNCommonTopologies++; + ncommon = ncom; +} + +void BuildTopologyDictionary::setThreshold(double thr, bool IB) +{ + mDictionary.resetMaps(IB); + setThresholdImpl(thr, + ((IB) ? mTopologyFrequencyIB : mTopologyFrequencyOB), + ((IB) ? mMapInfoIB : mMapInfoOB), + ((IB) ? mTopologyMapIB : mTopologyMapOB), + ((IB) ? mNCommonTopologiesIB : mNCommonTopologiesOB), + ((IB) ? mFrequencyThresholdIB : mFrequencyThresholdOB), + ((IB) ? mTotClustersIB : mTotClustersOB)); +} + +void BuildTopologyDictionary::setThresholdImpl(double thr, TopoFreq& tfreq, TopoInfo& tinfo, TopoStat& tstat, unsigned int& ncommon, double& freqthres, unsigned int ntot) +{ + setNCommonImpl(0, tfreq, tstat, ncommon, ntot); + freqthres = thr; + for (auto& q : tfreq) { + if (((double)q.first) / ntot > thr) { + ++ncommon; } else { break; } } - if (mNCommonTopologies >= itsmft::CompCluster::InvalidPatternID) { - mFrequencyThreshold = ((double)mTopologyFrequency[itsmft::CompCluster::InvalidPatternID - 1].first) / mTotClusters; - LOGP(warning, "Redefining prob. threshould from {} to {} to be below InvalidPatternID (was {})", thr, mFrequencyThreshold, mNCommonTopologies); - mNCommonTopologies = itsmft::CompCluster::InvalidPatternID - 1; + if (ncommon >= itsmft::CompCluster::InvalidPatternID) { + freqthres = ((double)tfreq[itsmft::CompCluster::InvalidPatternID - 1].first) / ntot; + LOGP(warning, "Redefining prob. threshold from {} to {} to be below InvalidPatternID (was {})", thr, freqthres, ntot); + ncommon = itsmft::CompCluster::InvalidPatternID - 1; } } -void BuildTopologyDictionary::setNCommon(unsigned int nCommon) +void BuildTopologyDictionary::setThresholdCumulative(double cumulative, bool IB) { - if (nCommon >= itsmft::CompCluster::InvalidPatternID) { - LOGP(warning, "Redefining nCommon from {} to {} to be below InvalidPatternID", nCommon, itsmft::CompCluster::InvalidPatternID - 1); - nCommon = itsmft::CompCluster::InvalidPatternID - 1; - } - mTopologyFrequency.clear(); - for (auto&& p : mTopologyMap) { // p os pair - mTopologyFrequency.emplace_back(p.second.countsTotal, p.first); + if (cumulative <= 0. || cumulative >= 1.) { + cumulative = 0.99; } - std::sort(mTopologyFrequency.begin(), mTopologyFrequency.end(), - [](const std::pair& couple1, - const std::pair& couple2) { return (couple1.first > couple2.first); }); - mNCommonTopologies = nCommon; - mDictionary.mCommonMap.clear(); - mDictionary.mGroupMap.clear(); - mFrequencyThreshold = ((double)mTopologyFrequency[mNCommonTopologies - 1].first) / mTotClusters; + + auto& freqTopo = ((IB) ? mTopologyFrequencyIB : mTopologyFrequencyOB); + auto& freqThres = ((IB) ? mFrequencyThresholdIB : mFrequencyThresholdOB); + auto& statTopo = ((IB) ? mTopologyMapIB : mTopologyMapOB); + auto& comTopo = ((IB) ? mNCommonTopologiesIB : mNCommonTopologiesOB); + auto ntot = ((IB) ? mTotClustersIB : mTotClustersOB); + + mDictionary.resetMaps(IB); + setNCommonImpl(0, freqTopo, statTopo, comTopo, ntot); + setThresholdCumulativeImpl(cumulative, freqTopo, comTopo, freqThres, ntot); } -void BuildTopologyDictionary::setThresholdCumulative(double cumulative) +void BuildTopologyDictionary::setThresholdCumulativeImpl(double cumulative, TopoFreq& tfreq, unsigned int& ncommon, double& freqthres, unsigned int ntot) { - mTopologyFrequency.clear(); - if (cumulative <= 0. || cumulative >= 1.) { - cumulative = 0.99; - } double totFreq = 0.; - for (auto&& p : mTopologyMap) { // p os pair - mTopologyFrequency.emplace_back(p.second.countsTotal, p.first); - } - std::sort(mTopologyFrequency.begin(), mTopologyFrequency.end(), - [](const std::pair& couple1, - const std::pair& couple2) { return (couple1.first > couple2.first); }); - mNCommonTopologies = 0; - mDictionary.mCommonMap.clear(); - mDictionary.mGroupMap.clear(); - for (auto& q : mTopologyFrequency) { - totFreq += ((double)(q.first)) / mTotClusters; + for (auto& q : tfreq) { + totFreq += ((double)(q.first)) / ntot; if (totFreq < cumulative) { - mNCommonTopologies++; - if (mNCommonTopologies >= itsmft::CompCluster::InvalidPatternID) { - totFreq -= ((double)(q.first)) / mTotClusters; - mNCommonTopologies--; + ++ncommon; + if (ncommon >= itsmft::CompCluster::InvalidPatternID) { + totFreq -= ((double)(q.first)) / ntot; + --ncommon; LOGP(warning, "Redefining cumulative threshould from {} to {} to be below InvalidPatternID)", cumulative, totFreq); } } else { break; } } - mFrequencyThreshold = ((double)(mTopologyFrequency[--mNCommonTopologies].first)) / mTotClusters; - while (std::fabs(((double)mTopologyFrequency[mNCommonTopologies].first) / mTotClusters - mFrequencyThreshold) < 1.e-15) { - mNCommonTopologies--; + freqthres = ((double)(tfreq[--ncommon].first)) / ntot; + while (std::fabs(((double)tfreq[ncommon--].first) / ntot - freqthres) < 1.e-15) { } - mFrequencyThreshold = ((double)mTopologyFrequency[mNCommonTopologies++].first) / mTotClusters; + freqthres = ((double)tfreq[ncommon++].first) / ntot; } void BuildTopologyDictionary::groupRareTopologies() { LOG(info) << "Dictionary finalisation"; - LOG(info) << "Number of clusters: " << mTotClusters; + LOG(info) << "Number of IB clusters: " << mTotClustersIB; + LOG(info) << "Number of OB clusters: " << mTotClustersOB; + + groupRareTopologiesImpl(mTopologyFrequencyIB, mMapInfoIB, mTopologyMapIB, mNCommonTopologiesIB, mFrequencyThresholdIB, mDictionary.mDataIB, mNCommonTopologiesIB); + groupRareTopologiesImpl(mTopologyFrequencyOB, mMapInfoOB, mTopologyMapOB, mNCommonTopologiesOB, mFrequencyThresholdOB, mDictionary.mDataOB, mNCommonTopologiesOB); + + LOG(info) << "Dictionay finalised"; + LOG(info) << "IB:"; + mDictionary.mDataIB.print(); + LOG(info) << "OB:"; + mDictionary.mDataOB.print(); +} +void BuildTopologyDictionary::groupRareTopologiesImpl(TopoFreq& tfreq, TopoInfo& tinfo, TopoStat& tstat, unsigned int& ncommon, double& freqthres, TopologyDictionaryData& data, unsigned int ntot) +{ double totFreq = 0.; - for (unsigned int j = 0; j < mNCommonTopologies; j++) { + for (unsigned int j = 0; j < ncommon; j++) { itsmft::GroupStruct gr; - gr.mHash = mTopologyFrequency[j].second; - gr.mFrequency = ((double)(mTopologyFrequency[j].first)) / mTotClusters; + gr.mHash = tfreq[j].second; + gr.mFrequency = ((double)(tfreq[j].first)) / ntot; totFreq += gr.mFrequency; // rough estimation for the error considering a8 uniform distribution - const auto& topo = mMapInfo.find(gr.mHash)->second; + const auto& topo = tinfo.find(gr.mHash)->second; gr.mErrX = std::sqrt(topo.mXsigma2); gr.mErrZ = std::sqrt(topo.mZsigma2); gr.mErr2X = topo.mXsigma2; @@ -168,11 +216,11 @@ void BuildTopologyDictionary::groupRareTopologies() gr.mNpixels = topo.mNpixels; gr.mPattern = topo.mPattern; gr.mIsGroup = false; - mDictionary.mVectorOfIDs.push_back(gr); + data.mVectorOfIDs.push_back(gr); if (j == int(itsmft::CompCluster::InvalidPatternID - 1)) { LOGP(warning, "Limiting N unique topologies to {}, threshold freq. to {}, cumulative freq. to {} to be below InvalidPatternID", j, gr.mFrequency, totFreq); - mNCommonTopologies = j; - mFrequencyThreshold = gr.mFrequency; + ncommon = j; + freqthres = gr.mFrequency; break; } } @@ -192,8 +240,8 @@ void BuildTopologyDictionary::groupRareTopologies() // Create a structure for a group of rare topologies itsmft::GroupStruct gr; gr.mHash = (((unsigned long)(grNum)) << 32) & 0xffffffff00000000; - gr.mErrX = its3::TopologyDictionary::RowClassSpan / std::sqrt(12 * std::min(10, rowBinEdge)); - gr.mErrZ = its3::TopologyDictionary::ColClassSpan / std::sqrt(12 * std::min(10, colBinEdge)); + gr.mErrX = its3::TopologyDictionary::RowClassSpan / std::sqrt(12.f * (float)std::min(10, rowBinEdge)); + gr.mErrZ = its3::TopologyDictionary::ColClassSpan / std::sqrt(12.f * (float)std::min(10, colBinEdge)); gr.mErr2X = gr.mErrX * gr.mErrX; gr.mErr2Z = gr.mErrZ * gr.mErrZ; gr.mXCOG = 0; @@ -227,58 +275,65 @@ void BuildTopologyDictionary::groupRareTopologies() int rs{}, cs{}, index{}; // Updating the counts for the groups of rare topologies - for (auto j{mNCommonTopologies}; j < mTopologyFrequency.size(); j++) { - unsigned long hash1 = mTopologyFrequency[j].second; - rs = mTopologyMap.find(hash1)->second.topology.getRowSpan(); - cs = mTopologyMap.find(hash1)->second.topology.getColumnSpan(); + for (auto j{ncommon}; j < tfreq.size(); j++) { + unsigned long hash1 = tfreq[j].second; + rs = tstat.find(hash1)->second.topology.getRowSpan(); + cs = tstat.find(hash1)->second.topology.getColumnSpan(); index = its3::LookUp::groupFinder(rs, cs); - tmp_GroupMap[index].second += mTopologyFrequency[j].first; + tmp_GroupMap[index].second += tfreq[j].first; } for (auto&& p : tmp_GroupMap) { itsmft::GroupStruct& group = p.second.first; - group.mFrequency = ((double)p.second.second) / mTotClusters; - mDictionary.mVectorOfIDs.push_back(group); + group.mFrequency = ((double)p.second.second) / ntot; + data.mVectorOfIDs.push_back(group); } // Sorting the dictionary preserving all unique topologies - std::sort(mDictionary.mVectorOfIDs.begin(), mDictionary.mVectorOfIDs.end(), [](const itsmft::GroupStruct& a, const itsmft::GroupStruct& b) { + std::sort(data.mVectorOfIDs.begin(), data.mVectorOfIDs.end(), [](const itsmft::GroupStruct& a, const itsmft::GroupStruct& b) { return (!a.mIsGroup) && b.mIsGroup ? true : (a.mIsGroup && (!b.mIsGroup) ? false : (a.mFrequency > b.mFrequency)); }); - if (mDictionary.mVectorOfIDs.size() >= itsmft::CompCluster::InvalidPatternID - 1) { + if (data.mVectorOfIDs.size() >= itsmft::CompCluster::InvalidPatternID - 1) { LOGP(warning, "Max allowed {} patterns is reached, stopping", itsmft::CompCluster::InvalidPatternID - 1); - mDictionary.mVectorOfIDs.resize(itsmft::CompCluster::InvalidPatternID - 1); + data.mVectorOfIDs.resize(itsmft::CompCluster::InvalidPatternID - 1); } // Sorting the dictionary to final form - std::sort(mDictionary.mVectorOfIDs.begin(), mDictionary.mVectorOfIDs.end(), [](const itsmft::GroupStruct& a, const itsmft::GroupStruct& b) { return a.mFrequency > b.mFrequency; }); + std::sort(data.mVectorOfIDs.begin(), data.mVectorOfIDs.end(), [](const itsmft::GroupStruct& a, const itsmft::GroupStruct& b) { return a.mFrequency > b.mFrequency; }); // Creating the map for common topologies - for (int iKey = 0; iKey < mDictionary.getSize(); iKey++) { - itsmft::GroupStruct& gr = mDictionary.mVectorOfIDs[iKey]; + for (int iKey = 0; iKey < data.mVectorOfIDs.size(); iKey++) { + itsmft::GroupStruct& gr = data.mVectorOfIDs[iKey]; if (!gr.mIsGroup) { - mDictionary.mCommonMap.emplace(gr.mHash, iKey); + data.mCommonMap.emplace(gr.mHash, iKey); if (gr.mPattern.getUsedBytes() == 1) { - mDictionary.mSmallTopologiesLUT[(gr.mPattern.getColumnSpan() - 1) * 255 + (int)gr.mPattern.getByte(2)] = iKey; + data.mSmallTopologiesLUT[(gr.mPattern.getColumnSpan() - 1) * 255 + (int)gr.mPattern.getByte(2)] = iKey; } } else { - mDictionary.mGroupMap.emplace((int)(gr.mHash >> 32) & 0x00000000ffffffff, iKey); + data.mGroupMap.emplace((int)(gr.mHash >> 32) & 0x00000000ffffffff, iKey); } } - LOG(info) << "Dictionay finalised"; - LOG(info) << "Number of keys: " << mDictionary.getSize(); - LOG(info) << "Number of common topologies: " << mDictionary.mCommonMap.size(); - LOG(info) << "Number of groups of rare topologies: " << mDictionary.mGroupMap.size(); } std::ostream& operator<<(std::ostream& os, const BuildTopologyDictionary& DB) { - for (unsigned int i = 0; i < DB.mNCommonTopologies; i++) { - const unsigned long& hash = DB.mTopologyFrequency[i].second; + os << "--- InnerBarrel\n"; + for (unsigned int i = 0; i < DB.mNCommonTopologiesIB; i++) { + const unsigned long& hash = DB.mTopologyFrequencyIB[i].second; + os << "Hash: " << hash << '\n'; + os << "counts: " << DB.mTopologyMapIB.find(hash)->second.countsTotal; + os << " (with bias provided: " << DB.mTopologyMapIB.find(hash)->second.countsWithBias << ")" << '\n'; + os << "sigmaX: " << std::sqrt(DB.mMapInfoIB.find(hash)->second.mXsigma2) << '\n'; + os << "sigmaZ: " << std::sqrt(DB.mMapInfoIB.find(hash)->second.mZsigma2) << '\n'; + os << DB.mTopologyMapIB.find(hash)->second.topology; + } + os << "--- OuterBarrel\n"; + for (unsigned int i = 0; i < DB.mNCommonTopologiesOB; i++) { + const unsigned long& hash = DB.mTopologyFrequencyOB[i].second; os << "Hash: " << hash << '\n'; - os << "counts: " << DB.mTopologyMap.find(hash)->second.countsTotal; - os << " (with bias provided: " << DB.mTopologyMap.find(hash)->second.countsWithBias << ")" << '\n'; - os << "sigmaX: " << std::sqrt(DB.mMapInfo.find(hash)->second.mXsigma2) << '\n'; - os << "sigmaZ: " << std::sqrt(DB.mMapInfo.find(hash)->second.mZsigma2) << '\n'; - os << DB.mTopologyMap.find(hash)->second.topology; + os << "counts: " << DB.mTopologyMapOB.find(hash)->second.countsTotal; + os << " (with bias provided: " << DB.mTopologyMapOB.find(hash)->second.countsWithBias << ")" << '\n'; + os << "sigmaX: " << std::sqrt(DB.mMapInfoOB.find(hash)->second.mXsigma2) << '\n'; + os << "sigmaZ: " << std::sqrt(DB.mMapInfoOB.find(hash)->second.mZsigma2) << '\n'; + os << DB.mTopologyMapOB.find(hash)->second.topology; } return os; } diff --git a/Detectors/Upgrades/ITS3/reconstruction/src/Clusterer.cxx b/Detectors/Upgrades/ITS3/reconstruction/src/Clusterer.cxx index d59a4859b9990..bce17b3759340 100644 --- a/Detectors/Upgrades/ITS3/reconstruction/src/Clusterer.cxx +++ b/Detectors/Upgrades/ITS3/reconstruction/src/Clusterer.cxx @@ -14,7 +14,6 @@ #include -#include "Framework/Logger.h" #include "ITS3Reconstruction/Clusterer.h" #include "ITS3Base/SegmentationMosaix.h" #include "SimulationDataFormat/MCTruthContainer.h" @@ -252,7 +251,7 @@ void Clusterer::ClustererThread::finishChip(ChipPixelData* curChipData, CompClus preClusterIndices[i2] = -1; } if (bbox.isAcceptableSize()) { - parent->streamCluster(pixArrBuff, &labelsBuff, bbox, parent->mPattIdConverter, compClusPtr, patternsPtr, labelsClusPtr, nlab); + parent->streamCluster(pixArrBuff, &labelsBuff, bbox, parent->mPattIdConverter, compClusPtr, patternsPtr, labelsClusPtr, nlab, constants::detID::isDetITS3(curChipData->getChipID())); } else { auto warnLeft = MaxHugeClusWarn - parent->mNHugeClus; if (warnLeft > 0) { @@ -278,7 +277,7 @@ void Clusterer::ClustererThread::finishChip(ChipPixelData* curChipData, CompClus } } if (!pixbuf.empty()) { // Stream a piece of cluster only if the reduced bounding box is not empty - parent->streamCluster(pixbuf, &labelsBuff, bboxT, parent->mPattIdConverter, compClusPtr, patternsPtr, labelsClusPtr, nlab, true); + parent->streamCluster(pixbuf, &labelsBuff, bboxT, parent->mPattIdConverter, compClusPtr, patternsPtr, labelsClusPtr, nlab, constants::detID::isDetITS3(curChipData->getChipID()), true); pixbuf.clear(); } bboxT.rowMin = bboxT.rowMax + 1; @@ -305,10 +304,12 @@ void Clusterer::ClustererThread::finishChipSingleHitFast(uint32_t hit, ChipPixel } } + auto ib = constants::detID::isDetITS3(curChipData->getChipID()); + // add to compact clusters, which must be always filled unsigned char patt[ClusterPattern::MaxPatternBytes]{0x1 << (7 - (0 % 8))}; // unrolled 1 hit version of full loop in finishChip - uint16_t pattID = (parent->mPattIdConverter.size() == 0) ? CompCluster::InvalidPatternID : parent->mPattIdConverter.findGroupID(1, 1, patt); - if ((pattID == CompCluster::InvalidPatternID || parent->mPattIdConverter.isGroup(pattID)) && patternsPtr) { + uint16_t pattID = (parent->mPattIdConverter.size(ib) == 0) ? CompCluster::InvalidPatternID : parent->mPattIdConverter.findGroupID(1, 1, ib, patt); + if ((pattID == CompCluster::InvalidPatternID || parent->mPattIdConverter.isGroup(pattID, ib)) && patternsPtr) { patternsPtr->emplace_back(1); // rowspan patternsPtr->emplace_back(1); // colspan patternsPtr->insert(patternsPtr->end(), std::begin(patt), std::begin(patt) + 1); diff --git a/Detectors/Upgrades/ITS3/reconstruction/src/ITS3ReconstructionLinkDef.h b/Detectors/Upgrades/ITS3/reconstruction/src/ITS3ReconstructionLinkDef.h index f19a7fcaba9ca..2ebd89970d9a1 100644 --- a/Detectors/Upgrades/ITS3/reconstruction/src/ITS3ReconstructionLinkDef.h +++ b/Detectors/Upgrades/ITS3/reconstruction/src/ITS3ReconstructionLinkDef.h @@ -16,6 +16,7 @@ #pragma link off all functions; #pragma link C++ class o2::its3::Clusterer + ; +#pragma link C++ class o2::its3::TopologyDictionaryData + ; #pragma link C++ class o2::its3::TopologyDictionary + ; #pragma link C++ class o2::its3::BuildTopologyDictionary + ; #pragma link C++ class o2::its3::LookUp + ; diff --git a/Detectors/Upgrades/ITS3/reconstruction/src/LookUp.cxx b/Detectors/Upgrades/ITS3/reconstruction/src/LookUp.cxx index caabfa6f2decb..e137e091dc631 100644 --- a/Detectors/Upgrades/ITS3/reconstruction/src/LookUp.cxx +++ b/Detectors/Upgrades/ITS3/reconstruction/src/LookUp.cxx @@ -31,7 +31,8 @@ LookUp::LookUp(std::string fileName) void LookUp::loadDictionary(std::string fileName) { mDictionary.readFromFile(fileName); - mTopologiesOverThreshold = mDictionary.mCommonMap.size(); + mTopologiesOverThresholdIB = mDictionary.mDataIB.mCommonMap.size(); + mTopologiesOverThresholdOB = mDictionary.mDataOB.mCommonMap.size(); } void LookUp::setDictionary(const its3::TopologyDictionary* dict) @@ -39,7 +40,8 @@ void LookUp::setDictionary(const its3::TopologyDictionary* dict) if (dict != nullptr) { mDictionary = *dict; } - mTopologiesOverThreshold = mDictionary.mCommonMap.size(); + mTopologiesOverThresholdIB = mDictionary.mDataIB.mCommonMap.size(); + mTopologiesOverThresholdOB = mDictionary.mDataOB.mCommonMap.size(); } int LookUp::groupFinder(int nRow, int nCol) @@ -61,25 +63,26 @@ int LookUp::groupFinder(int nRow, int nCol) return grNum; } -int LookUp::findGroupID(int nRow, int nCol, const unsigned char patt[itsmft::ClusterPattern::MaxPatternBytes]) const +int LookUp::findGroupID(int nRow, int nCol, bool IB, const unsigned char patt[itsmft::ClusterPattern::MaxPatternBytes]) const { + const auto& data = (IB) ? mDictionary.mDataIB : mDictionary.mDataOB; int nBits = nRow * nCol; if (nBits < 9) { // Small unique topology - int ID = mDictionary.mSmallTopologiesLUT[(nCol - 1) * 255 + (int)patt[0]]; + int ID = data.mSmallTopologiesLUT[(nCol - 1) * 255 + (int)patt[0]]; if (ID >= 0) { return ID; } } else { // Big unique topology unsigned long hash = itsmft::ClusterTopology::getCompleteHash(nRow, nCol, patt); - auto ret = mDictionary.mCommonMap.find(hash); - if (ret != mDictionary.mCommonMap.end()) { + auto ret = data.mCommonMap.find(hash); + if (ret != data.mCommonMap.end()) { return ret->second; } } - if (!mDictionary.mGroupMap.empty()) { // rare valid topology group + if (!data.mGroupMap.empty()) { // rare valid topology group int index = groupFinder(nRow, nCol); - auto res = mDictionary.mGroupMap.find(index); - return res == mDictionary.mGroupMap.end() ? itsmft::CompCluster::InvalidPatternID : res->second; + auto res = data.mGroupMap.find(index); + return res == data.mGroupMap.end() ? itsmft::CompCluster::InvalidPatternID : res->second; } return itsmft::CompCluster::InvalidPatternID; } diff --git a/Detectors/Upgrades/ITS3/reconstruction/src/TopologyDictionary.cxx b/Detectors/Upgrades/ITS3/reconstruction/src/TopologyDictionary.cxx index e64f63071c837..61ab051ffb565 100644 --- a/Detectors/Upgrades/ITS3/reconstruction/src/TopologyDictionary.cxx +++ b/Detectors/Upgrades/ITS3/reconstruction/src/TopologyDictionary.cxx @@ -23,9 +23,16 @@ ClassImp(o2::its3::TopologyDictionary); namespace o2::its3 { +void TopologyDictionaryData::print() const noexcept +{ + LOG(info) << "Number of keys: " << mVectorOfIDs.size(); + LOG(info) << "Number of common topologies: " << mCommonMap.size(); + LOG(info) << "Number of groups of rare topologies: " << mGroupMap.size(); +} + TopologyDictionary::TopologyDictionary() { - memset(mSmallTopologiesLUT, -1, STopoSize * sizeof(int)); + reset(); } TopologyDictionary::TopologyDictionary(const std::string& fileName) @@ -33,10 +40,43 @@ TopologyDictionary::TopologyDictionary(const std::string& fileName) readFromFile(fileName); } +void TopologyDictionary::print() const noexcept +{ + LOG(info) << "ITS3 TopologyDictionary"; + LOG(info) << "InnerBarrel"; + mDataIB.print(); + LOG(info) << "OuterBarrel"; + mDataOB.print(); +} + +void TopologyDictionary::reset() noexcept +{ + mDataIB.mSmallTopologiesLUT.fill(-1); + mDataOB.mSmallTopologiesLUT.fill(-1); + mDataIB.mVectorOfIDs.clear(); + mDataOB.mVectorOfIDs.clear(); +} + +void TopologyDictionary::resetMaps(bool IB) noexcept +{ + auto& data = (IB) ? mDataIB : mDataOB; + data.mCommonMap.clear(); + data.mGroupMap.clear(); +} + std::ostream& operator<<(std::ostream& os, const its3::TopologyDictionary& dict) { int ID = 0; - for (auto& p : dict.mVectorOfIDs) { + os << "--- InnerBarrel:\n"; + for (auto& p : dict.mDataIB.mVectorOfIDs) { + os << "ID: " << ID++ << " Hash: " << p.mHash << " ErrX: " << p.mErrX << " ErrZ : " << p.mErrZ << " xCOG: " << p.mXCOG << " zCOG: " << p.mZCOG << " Npixles: " << p.mNpixels << " Frequency: " << p.mFrequency << " isGroup : " << std::boolalpha << p.mIsGroup << '\n' + << p.mPattern << '\n' + << "*********************************************************" << '\n' + << '\n'; + } + ID = 0; + os << "--- OuterBarrel:\n"; + for (auto& p : dict.mDataOB.mVectorOfIDs) { os << "ID: " << ID++ << " Hash: " << p.mHash << " ErrX: " << p.mErrX << " ErrZ : " << p.mErrZ << " xCOG: " << p.mXCOG << " zCOG: " << p.mZCOG << " Npixles: " << p.mNpixels << " Frequency: " << p.mFrequency << " isGroup : " << std::boolalpha << p.mIsGroup << '\n' << p.mPattern << '\n' << "*********************************************************" << '\n' @@ -48,24 +88,36 @@ std::ostream& operator<<(std::ostream& os, const its3::TopologyDictionary& dict) void TopologyDictionary::writeBinaryFile(const std::string& outputfile) { std::ofstream file_output(outputfile, std::ios::out | std::ios::binary); - for (auto& p : mVectorOfIDs) { - file_output.write(reinterpret_cast(&p.mHash), sizeof(unsigned long)); - file_output.write(reinterpret_cast(&p.mErrX), sizeof(float)); - file_output.write(reinterpret_cast(&p.mErrZ), sizeof(float)); - file_output.write(reinterpret_cast(&p.mErr2X), sizeof(float)); - file_output.write(reinterpret_cast(&p.mErr2Z), sizeof(float)); - file_output.write(reinterpret_cast(&p.mXCOG), sizeof(float)); - file_output.write(reinterpret_cast(&p.mZCOG), sizeof(float)); - file_output.write(reinterpret_cast(&p.mNpixels), sizeof(int)); - file_output.write(reinterpret_cast(&p.mFrequency), sizeof(double)); - file_output.write(reinterpret_cast(&p.mIsGroup), sizeof(bool)); - file_output.write(const_cast(reinterpret_cast(&p.mPattern.getPattern())), - sizeof(unsigned char) * (itsmft::ClusterPattern::kExtendedPatternBytes)); + if (!file_output) { + throw std::runtime_error(fmt::format("Cannot open output file %s!", outputfile)); } + + auto writeData = [](auto& file_output, auto& data) { + auto size = data.mVectorOfIDs.size(); + file_output.write(reinterpret_cast(&size), sizeof(size)); + for (auto& p : data.mVectorOfIDs) { + file_output.write(reinterpret_cast(&p.mHash), sizeof(unsigned long)); + file_output.write(reinterpret_cast(&p.mErrX), sizeof(float)); + file_output.write(reinterpret_cast(&p.mErrZ), sizeof(float)); + file_output.write(reinterpret_cast(&p.mErr2X), sizeof(float)); + file_output.write(reinterpret_cast(&p.mErr2Z), sizeof(float)); + file_output.write(reinterpret_cast(&p.mXCOG), sizeof(float)); + file_output.write(reinterpret_cast(&p.mZCOG), sizeof(float)); + file_output.write(reinterpret_cast(&p.mNpixels), sizeof(int)); + file_output.write(reinterpret_cast(&p.mFrequency), sizeof(double)); + file_output.write(reinterpret_cast(&p.mIsGroup), sizeof(bool)); + file_output.write(const_cast(reinterpret_cast(&p.mPattern.getPattern())), + sizeof(unsigned char) * (itsmft::ClusterPattern::kExtendedPatternBytes)); + } + }; + + writeData(file_output, mDataIB); + writeData(file_output, mDataOB); + file_output.close(); } -int TopologyDictionary::readFromFile(const std::string& fname) +void TopologyDictionary::readFromFile(const std::string& fname) { LOGP(info, "Reading TopologyDictionary from File '{}'", fname); if (o2::utils::Str::endsWith(fname, ".root")) { @@ -76,59 +128,63 @@ int TopologyDictionary::readFromFile(const std::string& fname) } else { throw std::runtime_error(fmt::format("Unrecognized format {}", fname)); } - return 0; } -int TopologyDictionary::readBinaryFile(const std::string& fname) +void TopologyDictionary::readBinaryFile(const std::string& fname) { - mVectorOfIDs.clear(); - mCommonMap.clear(); - for (auto& p : mSmallTopologiesLUT) { - p = -1; - } + reset(); + std::ifstream in(fname.data(), std::ios::in | std::ios::binary); - itsmft::GroupStruct gr; - int groupID = 0; if (!in.is_open()) { LOG(error) << "The file " << fname << " coud not be opened"; throw std::runtime_error("The file coud not be opened"); } else { - while (in.read(reinterpret_cast(&gr.mHash), sizeof(unsigned long))) { - in.read(reinterpret_cast(&gr.mErrX), sizeof(float)); - in.read(reinterpret_cast(&gr.mErrZ), sizeof(float)); - in.read(reinterpret_cast(&gr.mErr2X), sizeof(float)); - in.read(reinterpret_cast(&gr.mErr2Z), sizeof(float)); - in.read(reinterpret_cast(&gr.mXCOG), sizeof(float)); - in.read(reinterpret_cast(&gr.mZCOG), sizeof(float)); - in.read(reinterpret_cast(&gr.mNpixels), sizeof(int)); - in.read(reinterpret_cast(&gr.mFrequency), sizeof(double)); - in.read(reinterpret_cast(&gr.mIsGroup), sizeof(bool)); - in.read(const_cast(reinterpret_cast(&gr.mPattern.getPattern())), sizeof(unsigned char) * (itsmft::ClusterPattern::kExtendedPatternBytes)); - mVectorOfIDs.push_back(gr); - if (!gr.mIsGroup) { - mCommonMap.insert(std::make_pair(gr.mHash, groupID)); - if (gr.mPattern.getUsedBytes() == 1) { - mSmallTopologiesLUT[(gr.mPattern.getColumnSpan() - 1) * 255 + (int)gr.mPattern.getByte(2)] = groupID; + + auto readData = [](auto& in, auto& data) { + int groupID = 0; + std::size_t size{}, cur{}; + itsmft::GroupStruct gr; + in.read(reinterpret_cast(&size), sizeof(std::size_t)); + while (cur++ != size) { + in.read(reinterpret_cast(&gr.mHash), sizeof(unsigned long)); + in.read(reinterpret_cast(&gr.mErrX), sizeof(float)); + in.read(reinterpret_cast(&gr.mErrZ), sizeof(float)); + in.read(reinterpret_cast(&gr.mErr2X), sizeof(float)); + in.read(reinterpret_cast(&gr.mErr2Z), sizeof(float)); + in.read(reinterpret_cast(&gr.mXCOG), sizeof(float)); + in.read(reinterpret_cast(&gr.mZCOG), sizeof(float)); + in.read(reinterpret_cast(&gr.mNpixels), sizeof(int)); + in.read(reinterpret_cast(&gr.mFrequency), sizeof(double)); + in.read(reinterpret_cast(&gr.mIsGroup), sizeof(bool)); + in.read(const_cast(reinterpret_cast(&gr.mPattern.getPattern())), sizeof(unsigned char) * (itsmft::ClusterPattern::kExtendedPatternBytes)); + data.mVectorOfIDs.push_back(gr); + if (!gr.mIsGroup) { + data.mCommonMap.insert(std::make_pair(gr.mHash, groupID)); + if (gr.mPattern.getUsedBytes() == 1) { + data.mSmallTopologiesLUT[(gr.mPattern.getColumnSpan() - 1) * 255 + (int)gr.mPattern.getByte(2)] = groupID; + } + } else { + data.mGroupMap.insert(std::make_pair((int)(gr.mHash >> 32) & 0x00000000ffffffff, groupID)); } - } else { - mGroupMap.insert(std::make_pair((int)(gr.mHash >> 32) & 0x00000000ffffffff, groupID)); + groupID++; } - groupID++; - } + }; + + readData(in, mDataIB); + readData(in, mDataOB); } in.close(); - return 0; } -TH1F* TopologyDictionary::getTopologyDistribution(const std::string_view hname) const +TH1F* TopologyDictionary::getTopologyDistribution(const std::string_view hname, bool IB) const { - int dictSize = getSize(); - auto* histo = new TH1F(hname.data(), ";Topology ID;Frequency", dictSize, -0.5, dictSize - 0.5); + int dictSize = getSize(IB); + auto* histo = new TH1F(hname.data(), Form("%s;Topology ID;Frequency", (IB) ? "InnerBarrel" : "OuterBarrel"), dictSize, -0.5, dictSize - 0.5); histo->SetFillColor(kRed); histo->SetFillStyle(3005); histo->SetDrawOption("histo"); for (int i = 0; i < dictSize; i++) { - histo->Fill(i, getFrequency(i)); + histo->Fill(i, getFrequency(i, IB)); } return histo; } @@ -140,13 +196,13 @@ math_utils::Point3D TopologyDictionary::getClusterCoordinates(const itsmft::C math_utils::Point3D locCl; if (!its3::constants::detID::isDetITS3(cl.getSensorID())) { o2::itsmft::SegmentationAlpide::detectorToLocalUnchecked(cl.getRow(), cl.getCol(), locCl); - locCl.SetX(locCl.X() + this->getXCOG(cl.getPatternID()) * itsmft::SegmentationAlpide::PitchRow); - locCl.SetZ(locCl.Z() + this->getZCOG(cl.getPatternID()) * itsmft::SegmentationAlpide::PitchCol); + locCl.SetX(locCl.X() + this->getXCOG(cl.getPatternID(), false) * itsmft::SegmentationAlpide::PitchRow); + locCl.SetZ(locCl.Z() + this->getZCOG(cl.getPatternID(), false) * itsmft::SegmentationAlpide::PitchCol); } else { auto layer = its3::constants::detID::getDetID2Layer(cl.getSensorID()); mIBSegmentations[layer].detectorToLocalUnchecked(cl.getRow(), cl.getCol(), locCl); - locCl.SetX(locCl.X() + this->getXCOG(cl.getPatternID()) * its3::SegmentationMosaix::PitchRow); - locCl.SetZ(locCl.Z() + this->getZCOG(cl.getPatternID()) * its3::SegmentationMosaix::PitchCol); + locCl.SetX(locCl.X() + this->getXCOG(cl.getPatternID(), true) * its3::SegmentationMosaix::PitchRow); + locCl.SetZ(locCl.Z() + this->getZCOG(cl.getPatternID(), true) * its3::SegmentationMosaix::PitchCol); float xCurved{0.f}, yCurved{0.f}; mIBSegmentations[layer].flatToCurved(locCl.X(), locCl.Y(), xCurved, yCurved); locCl.SetXYZ(xCurved, yCurved, locCl.Z()); From 146c47a8631a2dd08e952fc2ef416c079e27581f Mon Sep 17 00:00:00 2001 From: Felix Schlepper Date: Wed, 2 Apr 2025 14:28:44 +0200 Subject: [PATCH 08/10] ITS3: add opt. sanitizer --- Detectors/Upgrades/ITS3/CMakeLists.txt | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/Detectors/Upgrades/ITS3/CMakeLists.txt b/Detectors/Upgrades/ITS3/CMakeLists.txt index 7cbb5cffe4243..73ad4b9d53e37 100644 --- a/Detectors/Upgrades/ITS3/CMakeLists.txt +++ b/Detectors/Upgrades/ITS3/CMakeLists.txt @@ -9,7 +9,8 @@ # granted to it by virtue of its status as an Intergovernmental Organization # or submit itself to any jurisdiction. -#add_compile_options(-O0 -g -fPIC) +#add_compile_options(-O0 -g -fPIC -fsanitize=address) +#add_link_options(-fsanitize=address) add_subdirectory(data) add_subdirectory(simulation) From 3b12d72aebe54b3ed7163acf08bb379387e86602 Mon Sep 17 00:00:00 2001 From: Felix Schlepper Date: Thu, 10 Apr 2025 08:59:15 +0200 Subject: [PATCH 09/10] From d3b480494d31e8402f86bf24706212df98c0664a Mon Sep 17 00:00:00 2001 From: ALICE Action Bot Date: Thu, 10 Apr 2025 10:01:35 +0000 Subject: [PATCH 10/10] Please consider the following formatting changes --- .../Upgrades/ITS3/base/include/ITS3Base/SegmentationMosaix.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Detectors/Upgrades/ITS3/base/include/ITS3Base/SegmentationMosaix.h b/Detectors/Upgrades/ITS3/base/include/ITS3Base/SegmentationMosaix.h index 2e3da358b4efe..f8d4a784120a0 100644 --- a/Detectors/Upgrades/ITS3/base/include/ITS3Base/SegmentationMosaix.h +++ b/Detectors/Upgrades/ITS3/base/include/ITS3Base/SegmentationMosaix.h @@ -184,7 +184,7 @@ class SegmentationMosaix zCol = (static_cast(iCol) + 0.5f) * PitchCol - LengthH; } - bool detectorToLocal(int const row, int const col, math_utils::Point3D& loc) const noexcept + bool detectorToLocal(int const row, int const col, math_utils::Point3D& loc) const noexcept { float xRow{0.}, zCol{0.}; if (!detectorToLocal(row, col, xRow, zCol)) { @@ -194,7 +194,7 @@ class SegmentationMosaix return true; } - void detectorToLocalUnchecked(int const row, int const col, math_utils::Point3D& loc) const noexcept + void detectorToLocalUnchecked(int const row, int const col, math_utils::Point3D& loc) const noexcept { float xRow{0.}, zCol{0.}; detectorToLocalUnchecked(row, col, xRow, zCol);