diff --git a/DataFormats/Detectors/Upgrades/ALICE3/CMakeLists.txt b/DataFormats/Detectors/Upgrades/ALICE3/CMakeLists.txt index b3944c2e502d8..360b50d442d7d 100644 --- a/DataFormats/Detectors/Upgrades/ALICE3/CMakeLists.txt +++ b/DataFormats/Detectors/Upgrades/ALICE3/CMakeLists.txt @@ -10,3 +10,4 @@ # or submit itself to any jurisdiction. add_subdirectory(FD3) +add_subdirectory(TRK) diff --git a/DataFormats/Detectors/Upgrades/ALICE3/TRK/CMakeLists.txt b/DataFormats/Detectors/Upgrades/ALICE3/TRK/CMakeLists.txt new file mode 100644 index 0000000000000..c239a2a36845d --- /dev/null +++ b/DataFormats/Detectors/Upgrades/ALICE3/TRK/CMakeLists.txt @@ -0,0 +1,24 @@ +# Copyright 2019-2026 CERN and copyright holders of ALICE O2. +# See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +# All rights not expressly granted are reserved. +# +# This software is distributed under the terms of the GNU General Public +# License v3 (GPL Version 3), copied verbatim in the file "COPYING". +# +# In applying this license CERN does not waive the privileges and immunities +# granted to it by virtue of its status as an Intergovernmental Organization +# or submit itself to any jurisdiction. + +o2_add_library(DataFormatsTRK + SOURCES src/Cluster.cxx + src/ROFRecord.cxx + PUBLIC_LINK_LIBRARIES O2::CommonDataFormat + O2::DataFormatsITSMFT + O2::SimulationDataFormat +) + +o2_target_root_dictionary(DataFormatsTRK + HEADERS include/DataFormatsTRK/Cluster.h + include/DataFormatsTRK/ROFRecord.h + LINKDEF src/DataFormatsTRKLinkDef.h +) diff --git a/DataFormats/Detectors/Upgrades/ALICE3/TRK/include/DataFormatsTRK/Cluster.h b/DataFormats/Detectors/Upgrades/ALICE3/TRK/include/DataFormatsTRK/Cluster.h new file mode 100644 index 0000000000000..ec68191b3c43f --- /dev/null +++ b/DataFormats/Detectors/Upgrades/ALICE3/TRK/include/DataFormatsTRK/Cluster.h @@ -0,0 +1,38 @@ +// Copyright 2019-2026 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +#ifndef ALICEO2_DATAFORMATSTRK_CLUSTER_H +#define ALICEO2_DATAFORMATSTRK_CLUSTER_H + +#include +#include +#include + +namespace o2::trk +{ + +struct Cluster { + uint16_t chipID = 0; + uint16_t row = 0; + uint16_t col = 0; + uint16_t size = 1; + int16_t subDetID = -1; + int16_t layer = -1; + int16_t disk = -1; + + std::string asString() const; + + ClassDefNV(Cluster, 1); +}; + +} // namespace o2::trk + +#endif diff --git a/DataFormats/Detectors/Upgrades/ALICE3/TRK/include/DataFormatsTRK/ROFRecord.h b/DataFormats/Detectors/Upgrades/ALICE3/TRK/include/DataFormatsTRK/ROFRecord.h new file mode 100644 index 0000000000000..86ee31389fd5f --- /dev/null +++ b/DataFormats/Detectors/Upgrades/ALICE3/TRK/include/DataFormatsTRK/ROFRecord.h @@ -0,0 +1,75 @@ +// Copyright 2019-2026 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +#ifndef ALICEO2_DATAFORMATSTRK_ROFRECORD_H +#define ALICEO2_DATAFORMATSTRK_ROFRECORD_H + +#include "CommonDataFormat/InteractionRecord.h" +#include "CommonDataFormat/RangeReference.h" +#include +#include +#include + +namespace o2::trk +{ + +class ROFRecord +{ + public: + using EvIdx = o2::dataformats::RangeReference; + using BCData = o2::InteractionRecord; + using ROFtype = unsigned int; + + ROFRecord() = default; + ROFRecord(const BCData& bc, ROFtype rof, int idx, int n) + : mBCData(bc), mROFEntry(idx, n), mROFrame(rof) {} + + void setBCData(const BCData& bc) { mBCData = bc; } + void setROFrame(ROFtype rof) { mROFrame = rof; } + void setEntry(EvIdx entry) { mROFEntry = entry; } + void setFirstEntry(int idx) { mROFEntry.setFirstEntry(idx); } + void setNEntries(int n) { mROFEntry.setEntries(n); } + + const BCData& getBCData() const { return mBCData; } + BCData& getBCData() { return mBCData; } + EvIdx getEntry() const { return mROFEntry; } + EvIdx& getEntry() { return mROFEntry; } + int getNEntries() const { return mROFEntry.getEntries(); } + int getFirstEntry() const { return mROFEntry.getFirstEntry(); } + ROFtype getROFrame() const { return mROFrame; } + + std::string asString() const; + + private: + o2::InteractionRecord mBCData; + EvIdx mROFEntry; + ROFtype mROFrame = 0; + + ClassDefNV(ROFRecord, 1); +}; + +struct MC2ROFRecord { + using ROFtype = unsigned int; + + int eventRecordID = -1; + int rofRecordID = 0; + ROFtype minROF = 0; + ROFtype maxROF = 0; + + MC2ROFRecord() = default; + MC2ROFRecord(int evID, int rofRecID, ROFtype mnrof, ROFtype mxrof) : eventRecordID(evID), rofRecordID(rofRecID), minROF(mnrof), maxROF(mxrof) {} + + ClassDefNV(MC2ROFRecord, 1); +}; + +} // namespace o2::trk + +#endif diff --git a/DataFormats/Detectors/Upgrades/ALICE3/TRK/src/Cluster.cxx b/DataFormats/Detectors/Upgrades/ALICE3/TRK/src/Cluster.cxx new file mode 100644 index 0000000000000..6c96692ea5a9e --- /dev/null +++ b/DataFormats/Detectors/Upgrades/ALICE3/TRK/src/Cluster.cxx @@ -0,0 +1,28 @@ +// Copyright 2019-2026 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +#include "DataFormatsTRK/Cluster.h" +#include + +ClassImp(o2::trk::Cluster); + +namespace o2::trk +{ + +std::string Cluster::asString() const +{ + std::ostringstream stream; + stream << "chip=" << chipID << " row=" << row << " col=" << col << " size=" << size + << " subDet=" << subDetID << " layer=" << layer << " disk=" << disk; + return stream.str(); +} + +} // namespace o2::trk diff --git a/DataFormats/Detectors/Upgrades/ALICE3/TRK/src/DataFormatsTRKLinkDef.h b/DataFormats/Detectors/Upgrades/ALICE3/TRK/src/DataFormatsTRKLinkDef.h new file mode 100644 index 0000000000000..36528d9dd2c46 --- /dev/null +++ b/DataFormats/Detectors/Upgrades/ALICE3/TRK/src/DataFormatsTRKLinkDef.h @@ -0,0 +1,25 @@ +// Copyright 2019-2026 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +#ifdef __CLING__ + +#pragma link off all globals; +#pragma link off all classes; +#pragma link off all functions; + +#pragma link C++ class o2::trk::Cluster + ; +#pragma link C++ class std::vector < o2::trk::Cluster> + ; +#pragma link C++ class o2::trk::ROFRecord + ; +#pragma link C++ class std::vector < o2::trk::ROFRecord> + ; +#pragma link C++ class o2::trk::MC2ROFRecord + ; +#pragma link C++ class std::vector < o2::trk::MC2ROFRecord> + ; + +#endif diff --git a/DataFormats/Detectors/Upgrades/ALICE3/TRK/src/ROFRecord.cxx b/DataFormats/Detectors/Upgrades/ALICE3/TRK/src/ROFRecord.cxx new file mode 100644 index 0000000000000..79745f9854eb7 --- /dev/null +++ b/DataFormats/Detectors/Upgrades/ALICE3/TRK/src/ROFRecord.cxx @@ -0,0 +1,29 @@ +// Copyright 2019-2026 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +#include "DataFormatsTRK/ROFRecord.h" +#include + +ClassImp(o2::trk::ROFRecord); +ClassImp(o2::trk::MC2ROFRecord); + +namespace o2::trk +{ + +std::string ROFRecord::asString() const +{ + std::ostringstream stream; + stream << "IR=" << mBCData.asString() << " ROFrame=" << mROFrame + << " first=" << mROFEntry.getFirstEntry() << " n=" << mROFEntry.getEntries(); + return stream.str(); +} + +} // namespace o2::trk diff --git a/Detectors/Upgrades/ALICE3/TRK/macros/test/CMakeLists.txt b/Detectors/Upgrades/ALICE3/TRK/macros/test/CMakeLists.txt index d9908bbfeb1e5..edd9c785d89ce 100644 --- a/Detectors/Upgrades/ALICE3/TRK/macros/test/CMakeLists.txt +++ b/Detectors/Upgrades/ALICE3/TRK/macros/test/CMakeLists.txt @@ -28,4 +28,12 @@ o2_add_test_root_macro(CheckTracksCA.C O2::TRKBase O2::TRKSimulation O2::Steer - LABELS trk COMPILE_ONLY) \ No newline at end of file + LABELS trk COMPILE_ONLY) + +o2_add_test_root_macro(CheckClusters.C + PUBLIC_LINK_LIBRARIES O2::DataFormatsTRK + O2::SimulationDataFormat + O2::Framework + O2::TRKBase + O2::TRKSimulation + LABELS trk COMPILE_ONLY) diff --git a/Detectors/Upgrades/ALICE3/TRK/macros/test/CheckClusters.C b/Detectors/Upgrades/ALICE3/TRK/macros/test/CheckClusters.C new file mode 100644 index 0000000000000..327577102d86e --- /dev/null +++ b/Detectors/Upgrades/ALICE3/TRK/macros/test/CheckClusters.C @@ -0,0 +1,417 @@ +// Copyright 2019-2026 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +/// \file CheckClusters.C +/// \brief Macro to check TRK clusters and compare cluster positions to MC hit positions + +#if !defined(__CLING__) || defined(__ROOTCLING__) +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "DataFormatsTRK/Cluster.h" +#include "DataFormatsTRK/ROFRecord.h" +#include "TRKBase/GeometryTGeo.h" +#include "TRKBase/SegmentationChip.h" +#include "TRKSimulation/Hit.h" +#include "ITSMFTSimulation/AlpideSimResponse.h" +#include "CCDB/BasicCCDBManager.h" +#include "MathUtils/Cartesian.h" +#include "SimulationDataFormat/MCCompLabel.h" +#include "SimulationDataFormat/MCTruthContainer.h" +#include "DetectorsBase/GeometryManager.h" +#include "Framework/Logger.h" +#endif + +void CheckClusters(const std::string& clusfile = "o2clus_trk.root", + const std::string& hitfile = "o2sim_HitsTRK.root", + const std::string& inputGeom = "o2sim_geometry.root", + const std::string& ccdbUrl = "http://alice-ccdb.cern.ch", + long ccdbTimestamp = -1, + bool batch = false) +{ + gROOT->SetBatch(batch); + + using o2::MCCompLabel; + using ROFRec = o2::trk::ROFRecord; + using MC2ROF = o2::trk::MC2ROFRecord; + using HitVec = std::vector; + using MC2HITS_map = std::unordered_map; // maps (trackID << 32) + chipID -> hit index + + // ── Chip response (for hit-segment propagation to charge-collection plane) ── + // Fetches the same AlpideSimResponse from CCDB as the digitizer (IT3/Calib/APTSResponse) + // and computes Y-intersection planes with the same formulas from Digitizer::init() + auto& ccdbMgr = o2::ccdb::BasicCCDBManager::instance(); + ccdbMgr.setURL(ccdbUrl); + if (ccdbTimestamp > 0) { + ccdbMgr.setTimestamp(ccdbTimestamp); + } + auto* alpResp = ccdbMgr.get("IT3/Calib/APTSResponse"); + if (!alpResp) { + LOGP(fatal, "Cannot retrieve AlpideSimResponse from CCDB at {}", ccdbUrl); + return; + } + const float depthMax = alpResp->getDepthMax(); + + // ── Y-plane shifts: why VD and ML/OT need different values ──────────────── + // + // The APTS pixel response (AlpideSimResponse) uses an internal Y axis where: + // + // y = depthMax ── beam-entry (top) surface + // y = 0 ── charge-collection plane ← where clusters form + // y < 0 ── substrate (no response) + // + // The digitizer (Digitizer::init()) brings hit Y coordinates into this frame + // by adding a per-sub-detector shift before querying the response: + // + // y_APTS = y_local + shift [Digitizer.cxx ::processHit] + // + // The collection plane (y_APTS = 0) is therefore at y_local = −shift + // in the detector local frame. That is the Y value used here when + // propagating the MC hit segment to a single representative point. + // + // ── VD (vertex detector – curved sensors) ───────────────────────────────── + // After SegmentationChip::curvedToFlat() (convention: yFlat = dist − R): + // outer face (beam-entry): yFlat = +halfThickVD = +10 µm + // inner face (exit): yFlat = −halfThickVD = −10 µm + // The digitizer uses: + // + // mSimRespVDShift = depthMax − halfThickVD + // + // so the collection plane (y_APTS = 0) corresponds to: + // + // yPlaneVD = alice3resp::responseYShift = +5 µm + // + // i.e. 5 µm inside from the outer (entry) face. ✓ + // + // ── ML/OT (middle/outer tracker – flat sensors) ──────────────────────────── + // The local Y origin is at the GEOMETRIC CENTRE of the sensor volume. + // The outer (entry) surface is at y_local = +SiliconThicknessMLOT/2. + // The digitizer uses: + // + // mSimRespMLOTShift = depthMax − SiliconThicknessMLOT / 2 + // + // so the collection plane (y_APTS = 0) is at: + // + // yPlaneMLOT = SiliconThicknessMLOT/2 − depthMax + // + // ────────────────────────────────────────────────────────────────────────── + const float halfThicknessMLOT = o2::trk::SegmentationChip::SiliconThicknessMLOT / 2.f; + const float yPlaneVD = (float)o2::trk::constants::alice3resp::responseYShift; // VD: collection plane 5 µm inside outer (entry) face in flat local frame + const float yPlaneMLOT = halfThicknessMLOT - depthMax; // MLOT: entry @ +halfThick, collection depthMax below entry + LOGP(info, "Response depthMax = {:.4f} cm | VD Y-plane = {:.4f} cm | ML/OT Y-plane = {:.4f} cm", + depthMax, yPlaneVD, yPlaneMLOT); + + // ── Geometry ─────────────────────────────────────────────────────────────── + o2::base::GeometryManager::loadGeometry(inputGeom); + auto gman = o2::trk::GeometryTGeo::Instance(); + gman->fillMatrixCache(o2::math_utils::bit2Mask(o2::math_utils::TransformType::L2G)); + + // ── Hits ─────────────────────────────────────────────────────────────────── + TFile fileH(hitfile.data()); + auto* hitTree = dynamic_cast(fileH.Get("o2sim")); + if (!hitTree) { + LOGP(error, "Cannot find o2sim tree in {}", hitfile); + return; + } + std::vector mc2hitVec; + std::vector hitVecPool; + mc2hitVec.resize(hitTree->GetEntries()); + hitVecPool.resize(hitTree->GetEntries(), nullptr); + + // ── Clusters ─────────────────────────────────────────────────────────────── + TFile fileC(clusfile.data()); + auto* clusTree = dynamic_cast(fileC.Get("o2sim")); + if (!clusTree) { + LOGP(error, "Cannot find o2sim tree in {}", clusfile); + return; + } + + std::vector* clusArr = nullptr; + std::vector* rofRecVecP = nullptr; + std::vector* patternsPtr = nullptr; + clusTree->SetBranchAddress("TRKClusterComp", &clusArr); + clusTree->SetBranchAddress("TRKClustersROF", &rofRecVecP); + if (clusTree->GetBranch("TRKClusterPatt") != nullptr) { + clusTree->SetBranchAddress("TRKClusterPatt", &patternsPtr); + } + + o2::dataformats::MCTruthContainer* clusLabArr = nullptr; + std::vector mc2rofVec, *mc2rofVecP = &mc2rofVec; + bool hasMC = (clusTree->GetBranch("TRKClusterMCTruth") != nullptr); + if (hasMC) { + clusTree->SetBranchAddress("TRKClusterMCTruth", &clusLabArr); + clusTree->SetBranchAddress("TRKClustersMC2ROF", &mc2rofVecP); + } + + clusTree->GetEntry(0); + const unsigned int nROFRec = rofRecVecP ? (unsigned int)rofRecVecP->size() : 0u; + LOGP(info, "Number of ROF records: {}", nROFRec); + auto pattIt = patternsPtr ? patternsPtr->cbegin() : std::vector::const_iterator{}; + + // ── Build per-ROF MC event range ─────────────────────────────────────────── + std::vector mcEvMin(nROFRec, (int)hitTree->GetEntries()); + std::vector mcEvMax(nROFRec, -1); + if (hasMC) { + for (int imc = (int)mc2rofVec.size(); imc--;) { + const auto& mc2rof = mc2rofVec[imc]; + if (mc2rof.rofRecordID < 0) { + continue; + } + for (unsigned int irfd = mc2rof.maxROF - mc2rof.minROF + 1; irfd--;) { + unsigned int irof = mc2rof.rofRecordID + irfd; + if (irof >= nROFRec) { + continue; + } + if (mcEvMin[irof] > imc) { + mcEvMin[irof] = imc; + } + if (mcEvMax[irof] < imc) { + mcEvMax[irof] = imc; + } + } + } + } + + // ── Output ───────────────────────────────────────────────────────────────── + TFile fout("CheckClusters.root", "recreate"); + // columns: event, MC track label, + // local hit x/z (flat frame), global hit x/y/z (midpoint), + // global cluster x/y/z, local cluster x/z, + // residuals dx/dz (local, cluster - hit), + // ROF frame, cluster size, chipID, layer, disk, subDetID, row, col, pt + TNtuple nt("ntc", "TRK cluster ntuple", + "event:mcTrackID:hitLocX:hitLocZ:hitGlobX:hitGlobY:hitGlobZ:clusGlobX:clusGlobY:clusGlobZ:clusLocX:clusLocZ:rofFrame:clusSize:chipID:layer:disk:subdet:row:col:pt"); + + // ── Counters ─────────────────────────────────────────────────────────────── + long nTot{0}, nInvalidLabel{0}, nNoMCHit{0}, nValid{0}; + + // ── Main loop ────────────────────────────────────────────────────────────── + for (unsigned int irof = 0; irof < nROFRec; irof++) { + const auto& rofRec = (*rofRecVecP)[irof]; + + // Cache MC hit events for this ROF + if (hasMC) { + for (int im = mcEvMin[irof]; im <= mcEvMax[irof]; im++) { + if (hitVecPool[im] == nullptr) { + hitTree->SetBranchAddress("TRKHit", &hitVecPool[im]); + hitTree->GetEntry(im); + auto& mc2hit = mc2hitVec[im]; + const auto* hv = hitVecPool[im]; + for (int ih = (int)hv->size(); ih--;) { + const auto& hit = (*hv)[ih]; + uint64_t key = (uint64_t(hit.GetTrackID()) << 32) + hit.GetDetectorID(); + mc2hit.emplace(key, ih); + } + } + } + } + + for (int icl = 0; icl < rofRec.getNEntries(); icl++) { + const int clEntry = rofRec.getFirstEntry() + icl; + const auto& cluster = (*clusArr)[clEntry]; + nTot++; + + // ── Parse pattern → center-of-gravity within bounding box ────────── + // The cluster stores the bounding-box top-left pixel (row, col). + // The pattern stream encodes [rowSpan, colSpan, bitmap...] for each cluster. + // We accumulate pixel row/col offsets to obtain a sub-pixel CoG correction. + float cogDr{0.f}, cogDc{0.f}; // mean offsets from bbox origin (pixels) + if (patternsPtr) { + const uint8_t rowSpan = *pattIt++; + const uint8_t colSpan = *pattIt++; + const int nBytes = (rowSpan * colSpan + 7) / 8; + int nPix{0}, pixIdx{0}; + for (int ib = 0; ib < nBytes; ib++) { + const uint8_t byte = *pattIt++; + for (int bit = 7; bit >= 0 && pixIdx < rowSpan * colSpan; bit--, pixIdx++) { + if (byte & (1 << bit)) { + cogDr += pixIdx / colSpan; + cogDc += pixIdx % colSpan; + nPix++; + } + } + } + if (nPix > 1) { + cogDr /= nPix; + cogDc /= nPix; + } + } + + // ── Cluster local → global (CoG position) ───────────────────────────── + // Get local coords of the bounding-box corner pixel, then apply the + // fractional CoG displacement using the pixel pitch. + // Formula from detectorToLocalUnchecked: + // VD : xRow = 0.5*(width[lay]-pitchRow) - row*pitchRow → row↑ xRow↓ + // zCol = col*pitchCol + 0.5*(pitchCol-length) → col↑ zCol↑ + // MLOT: same structure with MLOT pitches + float clLocX{0.f}, clLocZ{0.f}; + o2::trk::SegmentationChip::detectorToLocalUnchecked( + cluster.row, cluster.col, clLocX, clLocZ, + cluster.subDetID, cluster.layer, cluster.disk); + const float pitchRow = (cluster.subDetID == 0) + ? o2::trk::SegmentationChip::PitchRowVD + : o2::trk::SegmentationChip::PitchRowMLOT; + const float pitchCol = (cluster.subDetID == 0) + ? o2::trk::SegmentationChip::PitchColVD + : o2::trk::SegmentationChip::PitchColMLOT; + clLocX -= cogDr * pitchRow; // increasing row → decreasing xRow + clLocZ += cogDc * pitchCol; // increasing col → increasing zCol + const float yResponse = (cluster.subDetID == 0) ? yPlaneVD : yPlaneMLOT; + // For VD the L2G matrix is built in the *curved* local frame (quasi-Cartesian, + // origin at the beam axis). Convert flat (clLocX, 0) → curved (xC, yC) first. + // For MLOT (flat sensors) the local frame is already Cartesian: pass directly. + // clLocX is already in the flat frame from detectorToLocalUnchecked + CoG and + // does NOT need any further transformation for the residual comparison. + o2::math_utils::Point3D locC; + if (cluster.subDetID == 0) { + auto cv = o2::trk::SegmentationChip::flatToCurved(cluster.layer, clLocX, 0.f); + locC = {cv.X(), cv.Y(), clLocZ}; + } else { + locC = {clLocX, yResponse, clLocZ}; + } + auto gloC = gman->getMatrixL2G(cluster.chipID)(locC); + + if (!hasMC || clusLabArr == nullptr) { + // No MC info: just fill geometry columns, leave residuals as 0 + std::array data = { + -1.f, -1.f, + 0.f, 0.f, 0.f, 0.f, 0.f, + (float)gloC.X(), (float)gloC.Y(), (float)gloC.Z(), + clLocX, clLocZ, + (float)rofRec.getROFrame(), (float)cluster.size, (float)cluster.chipID, + (float)cluster.layer, (float)cluster.disk, (float)cluster.subDetID, + (float)cluster.row, (float)cluster.col, -1.f}; + nt.Fill(data.data()); + continue; + } + + // ── MC label ─────────────────────────────────────────────────────── + const auto& labels = clusLabArr->getLabels(clEntry); + if (labels.empty() || !labels[0].isValid()) { + nInvalidLabel++; + continue; + } + const auto& lab = labels[0]; + const int trID = lab.getTrackID(); + const int evID = lab.getEventID(); + + // ── Find matching MC hit ──────────────────────────────────────────── + const auto& mc2hit = mc2hitVec[evID]; + uint64_t key = (uint64_t(trID) << 32) + cluster.chipID; + auto hitEntry = mc2hit.find(key); + if (hitEntry == mc2hit.end()) { + nNoMCHit++; + continue; + } + const auto& hit = (*hitVecPool[evID])[hitEntry->second]; + const float pt = TMath::Hypot(hit.GetPx(), hit.GetPy()); + + // ── Hit global midpoint ──────────────────────────────────────────── + const auto& gloHend = hit.GetPos(); + const auto& gloHsta = hit.GetPosStart(); + o2::math_utils::Point3D gloHmid( + 0.5f * (gloHend.X() + gloHsta.X()), + 0.5f * (gloHend.Y() + gloHsta.Y()), + 0.5f * (gloHend.Z() + gloHsta.Z())); + + // ── Hit global → local ───────────────────────────── + o2::math_utils::Point3D locHsta = gman->getMatrixL2G(cluster.chipID) ^ (gloHsta); // inverse L2G + o2::math_utils::Point3D locHend = gman->getMatrixL2G(cluster.chipID) ^ (gloHend); // inverse L2G + + // ── Propagate hit segment to the sensor response surface ─────────────── + // Rather than the geometric midpoint, find where the track segment crosses + // the response plane (y = responseYShift in the flat local frame). + // For VD (curved): convert both endpoints to flat frame first. + // For ML/OT (flat): use local coordinates directly. + float hitLocX{0.f}, hitLocZ{0.f}; + if (cluster.subDetID == 0) { // VD – curved sensor + auto flatSta = o2::trk::SegmentationChip::curvedToFlat(cluster.layer, locHsta.X(), locHsta.Y()); + auto flatEnd = o2::trk::SegmentationChip::curvedToFlat(cluster.layer, locHend.X(), locHend.Y()); + float x0 = flatSta.X(), y0 = flatSta.Y(), z0 = locHsta.Z(); + float dltx = flatEnd.X() - x0, dlty = flatEnd.Y() - y0, dltz = locHend.Z() - z0; + float r = (std::abs(dlty) > 1e-9f) ? (yPlaneVD - y0) / dlty : 0.5f; + hitLocX = x0 + r * dltx; + hitLocZ = z0 + r * dltz; + } else { // ML/OT – flat sensor + float x0 = locHsta.X(), y0 = locHsta.Y(), z0 = locHsta.Z(); + float dltx = locHend.X() - x0, dlty = locHend.Y() - y0, dltz = locHend.Z() - z0; + float r = (std::abs(dlty) > 1e-9f) ? (yPlaneMLOT - y0) / dlty : 0.5f; + hitLocX = x0 + r * dltx; + hitLocZ = z0 + r * dltz; + } + + nValid++; + std::array data = { + (float)evID, (float)trID, + hitLocX, hitLocZ, + (float)gloHmid.X(), (float)gloHmid.Y(), (float)gloHmid.Z(), + (float)gloC.X(), (float)gloC.Y(), (float)gloC.Z(), + clLocX, clLocZ, + (float)rofRec.getROFrame(), (float)cluster.size, (float)cluster.chipID, + (float)cluster.layer, (float)cluster.disk, (float)cluster.subDetID, + (float)cluster.row, (float)cluster.col, pt}; + nt.Fill(data.data()); + } + } + + // ── Summary ──────────────────────────────────────────────────────────────── + LOGP(info, "=== TRK Cluster vs Hit Summary ==="); + LOGP(info, "Total clusters: {}", nTot); + LOGP(info, "Valid (hit matched): {}", nValid); + LOGP(info, "Invalid/noise MC labels: {}", nInvalidLabel); + LOGP(info, "MC hit not found: {}", nNoMCHit); + // ── Visualisation ────────────────────────────────────────────────────────── + auto canvGlobal = new TCanvas("canvGlobal", "Cluster global positions", 1600, 800); + canvGlobal->Divide(2, 1); + canvGlobal->cd(1); + nt.Draw("clusGlobY:clusGlobX>>h_yx(500,-50,50,500,-50,50)", "", "colz"); + canvGlobal->cd(2); + nt.Draw("clusGlobY:clusGlobZ>>h_yz(500,-100,100,500,-50,50)", "", "colz"); + canvGlobal->SaveAs("trk_clusters_global.png"); + + auto canvRes = new TCanvas("canvRes", "Residuals (cluster - hit) [cm]", 1600, 1200); + canvRes->Divide(2, 3); + canvRes->cd(1)->SetLogy(); + nt.Draw("hitLocX-clusLocX>>h_dx_VD(200,-0.02,0.02)", "subdet==0&&event>=0"); + canvRes->cd(2)->SetLogy(); + nt.Draw("hitLocZ-clusLocZ>>h_dz_VD(200,-0.02,0.02)", "subdet==0&&event>=0"); + canvRes->cd(3)->SetLogy(); + nt.Draw("hitLocX-clusLocX>>h_dx_MLOT(200,-0.02,0.02)", "subdet==1&&event>=0"); + canvRes->cd(4)->SetLogy(); + nt.Draw("hitLocZ-clusLocZ>>h_dz_MLOT(200,-0.02,0.02)", "subdet==1&&event>=0"); + canvRes->cd(5)->SetLogz(); + nt.Draw("hitLocX-clusLocX:hitLocZ-clusLocZ>>h_dxdz_VD(200,-0.02,0.02,200,-0.02,0.02)", "subdet==0&&event>=0", "colz"); + canvRes->cd(6); + nt.Draw("hitLocX-clusLocX:hitLocZ-clusLocZ>>h_dxdz_MLOT(200,-0.02,0.02,200,-0.02,0.02)", "subdet==1&&event>=0", "colz"); + canvRes->SaveAs("trk_residuals.png"); + + auto canvResVsLayer = new TCanvas("canvResVsLayer", "Residuals vs layer", 1600, 600); + canvResVsLayer->Divide(2, 1); + canvResVsLayer->cd(1); + nt.Draw("hitLocX-clusLocX:layer>>h_dx_vs_lay(20,0,20,200,-0.02,0.02)", "event>=0", "prof"); + canvResVsLayer->cd(2); + nt.Draw("hitLocZ-clusLocZ:layer>>h_dz_vs_lay(20,0,20,200,-0.02,0.02)", "event>=0", "prof"); + canvResVsLayer->SaveAs("trk_residuals_vs_layer.png"); + + fout.cd(); + nt.Write(); + fout.Close(); + + LOGP(info, "Output saved to CheckClusters.root and PNG files"); +} diff --git a/Detectors/Upgrades/ALICE3/TRK/reconstruction/CMakeLists.txt b/Detectors/Upgrades/ALICE3/TRK/reconstruction/CMakeLists.txt index 01ddc783d192b..b9866c7d6aa4d 100644 --- a/Detectors/Upgrades/ALICE3/TRK/reconstruction/CMakeLists.txt +++ b/Detectors/Upgrades/ALICE3/TRK/reconstruction/CMakeLists.txt @@ -12,12 +12,14 @@ o2_add_library(TRKReconstruction TARGETVARNAME targetName SOURCES src/TimeFrame.cxx + src/Clusterer.cxx PUBLIC_LINK_LIBRARIES O2::ITStracking O2::GPUCommon Microsoft.GSL::GSL O2::CommonConstants O2::DataFormatsITSMFT + O2::DataFormatsTRK O2::SimulationDataFormat O2::ITSBase O2::ITSReconstruction @@ -31,4 +33,5 @@ o2_add_library(TRKReconstruction o2_target_root_dictionary(TRKReconstruction HEADERS include/TRKReconstruction/TimeFrame.h + include/TRKReconstruction/Clusterer.h LINKDEF src/TRKReconstructionLinkDef.h) diff --git a/Detectors/Upgrades/ALICE3/TRK/reconstruction/include/TRKReconstruction/Clusterer.h b/Detectors/Upgrades/ALICE3/TRK/reconstruction/include/TRKReconstruction/Clusterer.h new file mode 100644 index 0000000000000..abddafa312fb9 --- /dev/null +++ b/Detectors/Upgrades/ALICE3/TRK/reconstruction/include/TRKReconstruction/Clusterer.h @@ -0,0 +1,182 @@ +// Copyright 2019-2026 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +/// \file Clusterer.h +/// \brief Definition of the TRK cluster finder + +#ifndef ALICEO2_TRK_CLUSTERER_H +#define ALICEO2_TRK_CLUSTERER_H + +// uncomment to allow diagonal clusters, e.g. |* | +// | *| +#define _ALLOW_DIAGONAL_TRK_CLUSTERS_ + +#include "DataFormatsITSMFT/Digit.h" +#include "DataFormatsITSMFT/ROFRecord.h" +#include "DataFormatsITSMFT/ClusterPattern.h" +#include "DataFormatsTRK/Cluster.h" +#include "DataFormatsTRK/ROFRecord.h" +#include "SimulationDataFormat/ConstMCTruthContainer.h" +#include "SimulationDataFormat/MCCompLabel.h" +#include "SimulationDataFormat/MCTruthContainer.h" +#include "TRKBase/Specs.h" +#include +#include +#include +#include +#include +#include + +namespace o2::trk +{ + +class GeometryTGeo; + +class Clusterer +{ + public: + static constexpr int MaxLabels = 10; + static constexpr int MaxHugeClusWarn = 5; + + using Digit = o2::itsmft::Digit; + using DigROFRecord = o2::itsmft::ROFRecord; + using DigMC2ROFRecord = o2::itsmft::MC2ROFRecord; + using ClusterTruth = o2::dataformats::MCTruthContainer; + using ConstDigitTruth = o2::dataformats::ConstMCTruthContainerView; + using Label = o2::MCCompLabel; + + //---------------------------------------------- + struct BBox { + uint16_t chipID = 0xffff; + uint16_t rowMin = 0xffff, colMin = 0xffff; + uint16_t rowMax = 0, colMax = 0; + explicit BBox(uint16_t c) : chipID(c) {} + bool isInside(uint16_t r, uint16_t c) const { return r >= rowMin && r <= rowMax && c >= colMin && c <= colMax; } + uint16_t rowSpan() const { return rowMax - rowMin + 1; } + uint16_t colSpan() const { return colMax - colMin + 1; } + bool isAcceptableSize() const + { + return rowSpan() <= o2::itsmft::ClusterPattern::MaxRowSpan && + colSpan() <= o2::itsmft::ClusterPattern::MaxColSpan; + } + void adjust(uint16_t r, uint16_t c) + { + if (r < rowMin) { + rowMin = r; + } + if (r > rowMax) { + rowMax = r; + } + if (c < colMin) { + colMin = c; + } + if (c > colMax) { + colMax = c; + } + } + }; + + //---------------------------------------------- + struct ClustererThread { + Clusterer* parent = nullptr; + // column buffers (pre-cluster state); extra sentinel entries at [0] and [size-1] + int* column1 = nullptr; + int* column2 = nullptr; + int* curr = nullptr; ///< current column pre-cluster indices + int* prev = nullptr; ///< previous column pre-cluster indices + int size = constants::moduleMLOT::chip::nRows + 2; ///< reallocated per chip in initChip + + // pixels[i] = {next_in_chain, global_digit_index} + std::vector> pixels; + std::vector preClusterHeads; + std::vector preClusterIndices; + uint16_t currCol = 0xffff; + bool noLeftCol = true; + + std::array labelsBuff; ///< MC label buffer for one cluster + std::vector> pixArrBuff; ///< (row,col) pixel buffer for pattern + + // per-thread output (accumulated, then merged back by caller) + std::vector clusters; + std::vector patterns; + ClusterTruth labels; + + ///< reset column buffer + void resetColumn(int* buff) const { std::memset(buff, -1, sizeof(int) * (size - 2)); } + ///< swap current and previous column buffers + void swapColumnBuffers() { std::swap(prev, curr); } + + ///< append pixel ip to the pre-cluster headed at preClusterIndex + void expandPreCluster(uint32_t ip, uint16_t row, int preClusterIndex) + { + auto& firstIndex = preClusterHeads[preClusterIndices[preClusterIndex]]; + pixels.emplace_back(firstIndex, ip); + firstIndex = pixels.size() - 1; + curr[row] = preClusterIndex; + } + + ///< start a new pre-cluster with pixel ip at given row + void addNewPreCluster(uint32_t ip, uint16_t row) + { + preClusterHeads.push_back(pixels.size()); + pixels.emplace_back(-1, ip); + int lastIndex = preClusterIndices.size(); + preClusterIndices.push_back(lastIndex); + curr[row] = lastIndex; + } + + void fetchMCLabels(uint32_t digID, const ConstDigitTruth* labelsDig, int& nfilled); + void initChip(gsl::span digits, uint32_t first, GeometryTGeo* geom); + void updateChip(gsl::span digits, uint32_t ip); + void finishChip(gsl::span digits, + const ConstDigitTruth* labelsDigPtr, ClusterTruth* labelsClusPtr, + GeometryTGeo* geom); + void finishChipSingleHitFast(gsl::span digits, uint32_t hit, + const ConstDigitTruth* labelsDigPtr, ClusterTruth* labelsClusPtr, + GeometryTGeo* geom); + void processChip(gsl::span digits, int chipFirst, int chipN, + std::vector* clustersOut, std::vector* patternsOut, + const ConstDigitTruth* labelsDigPtr, ClusterTruth* labelsClusPtr, + GeometryTGeo* geom); + void streamCluster(const BBox& bbox, const std::vector>& pixbuf, + uint32_t totalCharge, bool doLabels, int nlab, + uint16_t chipID, int subDetID, int layer, int disk); + + ~ClustererThread() + { + delete[] column1; + delete[] column2; + } + explicit ClustererThread(Clusterer* par = nullptr) : parent(par) {} + ClustererThread(const ClustererThread&) = delete; + ClustererThread& operator=(const ClustererThread&) = delete; + }; + //---------------------------------------------- + + void process(gsl::span digits, + gsl::span digitROFs, + std::vector& clusters, + std::vector& patterns, + std::vector& clusterROFs, + const ConstDigitTruth* digitLabels = nullptr, + ClusterTruth* clusterLabels = nullptr, + gsl::span digMC2ROFs = {}, + std::vector* clusterMC2ROFs = nullptr); + + private: + int mNHugeClus = 0; + std::unique_ptr mThread; + std::vector mSortIdx; ///< reusable per-ROF sort buffer +}; + +} // namespace o2::trk + +#endif diff --git a/Detectors/Upgrades/ALICE3/TRK/reconstruction/src/Clusterer.cxx b/Detectors/Upgrades/ALICE3/TRK/reconstruction/src/Clusterer.cxx new file mode 100644 index 0000000000000..bdaa76319c1f2 --- /dev/null +++ b/Detectors/Upgrades/ALICE3/TRK/reconstruction/src/Clusterer.cxx @@ -0,0 +1,419 @@ +// Copyright 2019-2026 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +/// \file Clusterer.cxx +/// \brief Implementation of the TRK cluster finder + +#include "TRKReconstruction/Clusterer.h" +#include "TRKBase/GeometryTGeo.h" + +#include +#include + +namespace o2::trk +{ + +//__________________________________________________ +void Clusterer::process(gsl::span digits, + gsl::span digitROFs, + std::vector& clusters, + std::vector& patterns, + std::vector& clusterROFs, + const ConstDigitTruth* digitLabels, + ClusterTruth* clusterLabels, + gsl::span digMC2ROFs, + std::vector* clusterMC2ROFs) +{ + if (!mThread) { + mThread = std::make_unique(this); + } + + auto* geom = o2::trk::GeometryTGeo::Instance(); + + for (size_t iROF = 0; iROF < digitROFs.size(); ++iROF) { + const auto& inROF = digitROFs[iROF]; + const auto outFirst = static_cast(clusters.size()); + const int first = inROF.getFirstEntry(); + const int nEntries = inROF.getNEntries(); + + if (nEntries == 0) { + clusterROFs.emplace_back(inROF.getBCData(), inROF.getROFrame(), outFirst, 0); + continue; + } + + // Sort digit indices within this ROF by (chipID, col, row) so we can process + // chip by chip, column by column -- the same ordering the ALPIDE scanner expects. + mSortIdx.resize(nEntries); + std::iota(mSortIdx.begin(), mSortIdx.end(), first); + std::sort(mSortIdx.begin(), mSortIdx.end(), [&digits](int a, int b) { + const auto& da = digits[a]; + const auto& db = digits[b]; + if (da.getChipIndex() != db.getChipIndex()) { + return da.getChipIndex() < db.getChipIndex(); + } + if (da.getColumn() != db.getColumn()) { + return da.getColumn() < db.getColumn(); + } + return da.getRow() < db.getRow(); + }); + + // Process one chip at a time + int sliceStart = 0; + while (sliceStart < nEntries) { + const int chipFirst = sliceStart; + const uint16_t chipID = digits[mSortIdx[sliceStart]].getChipIndex(); + while (sliceStart < nEntries && digits[mSortIdx[sliceStart]].getChipIndex() == chipID) { + ++sliceStart; + } + const int chipN = sliceStart - chipFirst; + + mThread->processChip(digits, chipFirst, chipN, &clusters, &patterns, digitLabels, clusterLabels, geom); + } + + clusterROFs.emplace_back(inROF.getBCData(), inROF.getROFrame(), + outFirst, static_cast(clusters.size()) - outFirst); + } + + if (clusterMC2ROFs && !digMC2ROFs.empty()) { + clusterMC2ROFs->reserve(clusterMC2ROFs->size() + digMC2ROFs.size()); + for (const auto& in : digMC2ROFs) { + clusterMC2ROFs->emplace_back(in.eventRecordID, in.rofRecordID, in.minROF, in.maxROF); + } + } +} + +//__________________________________________________ +void Clusterer::ClustererThread::processChip(gsl::span digits, + int chipFirst, int chipN, + std::vector* clustersOut, + std::vector* patternsOut, + const ConstDigitTruth* labelsDigPtr, + ClusterTruth* labelsClusPtr, + GeometryTGeo* geom) +{ + // chipFirst and chipN are relative to mSortIdx (i.e. mSortIdx[chipFirst..chipFirst+chipN-1] + // are the global digit indices for this chip, already sorted by col then row). + // We use parent->mSortIdx to resolve the global index of each pixel. + const auto& sortIdx = parent->mSortIdx; + + if (chipN == 1) { + finishChipSingleHitFast(digits, sortIdx[chipFirst], labelsDigPtr, labelsClusPtr, geom); + } else { + initChip(digits, sortIdx[chipFirst], geom); + for (int i = chipFirst + 1; i < chipFirst + chipN; ++i) { + updateChip(digits, sortIdx[i]); + } + finishChip(digits, labelsDigPtr, labelsClusPtr, geom); + } + + // Flush per-thread output into the caller's containers + if (!clusters.empty()) { + clustersOut->insert(clustersOut->end(), clusters.begin(), clusters.end()); + clusters.clear(); + } + if (!patterns.empty()) { + patternsOut->insert(patternsOut->end(), patterns.begin(), patterns.end()); + patterns.clear(); + } + if (labelsClusPtr && labels.getNElements()) { + labelsClusPtr->mergeAtBack(labels); + labels.clear(); + } +} + +//__________________________________________________ +void Clusterer::ClustererThread::initChip(gsl::span digits, uint32_t first, GeometryTGeo* geom) +{ + const uint16_t chipID = digits[first].getChipIndex(); + + // Determine the number of rows for this chip's sensor type + size = constants::moduleMLOT::chip::nRows + 2; // default for ML/OT + if (geom) { + if (geom->getSubDetID(chipID) == 0) { // VD + const int layer = geom->getLayer(chipID); + size = constants::VD::petal::layer::nRows[layer] + 2; + } + } + + delete[] column1; + delete[] column2; + column1 = new int[size]; + column2 = new int[size]; + column1[0] = column1[size - 1] = -1; + column2[0] = column2[size - 1] = -1; + prev = column1 + 1; + curr = column2 + 1; + resetColumn(curr); + + pixels.clear(); + preClusterHeads.clear(); + preClusterIndices.clear(); + + const auto& pix = digits[first]; + currCol = pix.getColumn(); + curr[pix.getRow()] = 0; + preClusterHeads.push_back(0); + preClusterIndices.push_back(0); + pixels.emplace_back(-1, first); + noLeftCol = true; +} + +//__________________________________________________ +void Clusterer::ClustererThread::updateChip(gsl::span digits, uint32_t ip) +{ + const auto& pix = digits[ip]; + uint16_t row = pix.getRow(); + + if (currCol != pix.getColumn()) { + swapColumnBuffers(); + resetColumn(curr); + noLeftCol = false; + if (pix.getColumn() > currCol + 1) { + // gap: no connection with previous column + currCol = pix.getColumn(); + addNewPreCluster(ip, row); + noLeftCol = true; + return; + } + currCol = pix.getColumn(); + } + + bool orphan = true; + + if (noLeftCol) { + if (curr[row - 1] >= 0) { + expandPreCluster(ip, row, curr[row - 1]); + return; + } + } else { +#ifdef _ALLOW_DIAGONAL_TRK_CLUSTERS_ + int neighbours[]{curr[row - 1], prev[row], prev[row + 1], prev[row - 1]}; +#else + int neighbours[]{curr[row - 1], prev[row]}; +#endif + for (auto pci : neighbours) { + if (pci < 0) { + continue; + } + if (orphan) { + expandPreCluster(ip, row, pci); + orphan = false; + continue; + } + // merge two pre-clusters: assign the smaller index to both + if (preClusterIndices[pci] < preClusterIndices[curr[row]]) { + preClusterIndices[curr[row]] = preClusterIndices[pci]; + } else { + preClusterIndices[pci] = preClusterIndices[curr[row]]; + } + } + } + if (orphan) { + addNewPreCluster(ip, row); + } +} + +//__________________________________________________ +void Clusterer::ClustererThread::finishChip(gsl::span digits, + const ConstDigitTruth* labelsDigPtr, + ClusterTruth* labelsClusPtr, + GeometryTGeo* geom) +{ + const uint16_t chipID = digits[pixels[0].second].getChipIndex(); + + for (size_t i1 = 0; i1 < preClusterHeads.size(); ++i1) { + auto ci = preClusterIndices[i1]; + if (ci < 0) { + continue; + } + BBox bbox(chipID); + int nlab = 0; + uint32_t totalCharge = 0; + pixArrBuff.clear(); + + // Walk the linked list for this pre-cluster head + auto collectPixels = [&](int head) { + int next = head; + while (next >= 0) { + const auto& pixEntry = pixels[next]; + const auto& d = digits[pixEntry.second]; + uint16_t r = d.getRow(), c = d.getColumn(); + pixArrBuff.emplace_back(r, c); + bbox.adjust(r, c); + totalCharge += d.getCharge(); + if (labelsClusPtr) { + fetchMCLabels(pixEntry.second, labelsDigPtr, nlab); + } + next = pixEntry.first; + } + }; + + collectPixels(preClusterHeads[i1]); + preClusterIndices[i1] = -1; + + for (size_t i2 = i1 + 1; i2 < preClusterHeads.size(); ++i2) { + if (preClusterIndices[i2] != ci) { + continue; + } + collectPixels(preClusterHeads[i2]); + preClusterIndices[i2] = -1; + } + + // Determine geometry info + int subDetID = -1, layer = -1, disk = -1; + if (geom) { + subDetID = geom->getSubDetID(chipID); + layer = geom->getLayer(chipID); + disk = geom->getDisk(chipID); + } + + const bool doLabels = (labelsClusPtr != nullptr); + if (bbox.isAcceptableSize()) { + streamCluster(bbox, pixArrBuff, totalCharge, doLabels, nlab, chipID, subDetID, layer, disk); + } else { + // Huge cluster: split into MaxRowSpan x MaxColSpan tiles (same as ITS3) + auto warnLeft = MaxHugeClusWarn - parent->mNHugeClus; + if (warnLeft > 0) { + LOGP(warn, "Splitting huge TRK cluster: chipID {}, rows {}:{} cols {}:{}{}", + chipID, bbox.rowMin, bbox.rowMax, bbox.colMin, bbox.colMax, + warnLeft == 1 ? " (further warnings muted)" : ""); + parent->mNHugeClus++; + } + BBox bboxT(chipID); + bboxT.colMin = bbox.colMin; + do { + bboxT.rowMin = bbox.rowMin; + bboxT.colMax = std::min(bbox.colMax, uint16_t(bboxT.colMin + o2::itsmft::ClusterPattern::MaxColSpan - 1)); + do { + bboxT.rowMax = std::min(bbox.rowMax, uint16_t(bboxT.rowMin + o2::itsmft::ClusterPattern::MaxRowSpan - 1)); + std::vector> subPix; + uint32_t subCharge = 0; + for (const auto& [r, c] : pixArrBuff) { + if (bboxT.isInside(r, c)) { + subPix.emplace_back(r, c); + subCharge += 1; + } + } + if (!subPix.empty()) { + streamCluster(bboxT, subPix, subCharge, doLabels, nlab, chipID, subDetID, layer, disk); + } + bboxT.rowMin = bboxT.rowMax + 1; + } while (bboxT.rowMin <= bbox.rowMax); + bboxT.colMin = bboxT.colMax + 1; + } while (bboxT.colMin <= bbox.colMax); + } + } + // flush per-thread output to the caller via processChip +} + +//__________________________________________________ +void Clusterer::ClustererThread::finishChipSingleHitFast(gsl::span digits, uint32_t hit, + const ConstDigitTruth* labelsDigPtr, + ClusterTruth* labelsClusPtr, + GeometryTGeo* geom) +{ + const auto& d = digits[hit]; + const uint16_t chipID = d.getChipIndex(); + const uint16_t row = d.getRow(); + const uint16_t col = d.getColumn(); + + if (labelsClusPtr) { + int nlab = 0; + fetchMCLabels(hit, labelsDigPtr, nlab); + const auto cnt = static_cast(clusters.size()); + for (int i = nlab; i--;) { + labels.addElement(cnt, labelsBuff[i]); + } + } + + // 1×1 pattern: rowSpan=1, colSpan=1, one byte = 0x80 + patterns.emplace_back(1); + patterns.emplace_back(1); + patterns.emplace_back(0x80); + + Cluster cluster; + cluster.chipID = chipID; + cluster.row = row; + cluster.col = col; + cluster.size = 1; + if (geom) { + cluster.subDetID = geom->getSubDetID(chipID); + cluster.layer = geom->getLayer(chipID); + cluster.disk = geom->getDisk(chipID); + } + clusters.emplace_back(cluster); +} + +//__________________________________________________ +void Clusterer::ClustererThread::streamCluster(const BBox& bbox, + const std::vector>& pixbuf, + uint32_t totalCharge, + bool doLabels, int nlab, + uint16_t chipID, int subDetID, int layer, int disk) +{ + if (doLabels) { + const auto cnt = static_cast(clusters.size()); + for (int i = nlab; i--;) { + labels.addElement(cnt, labelsBuff[i]); // accumulate in thread-local buffer + } + } + + const uint16_t rowSpanW = bbox.rowSpan(); + const uint16_t colSpanW = bbox.colSpan(); + + // Encode the pixel pattern bitmap (rowSpan, colSpan, bytes...) + std::array patt{}; + for (const auto& [r, c] : pixbuf) { + uint32_t ir = r - bbox.rowMin, ic = c - bbox.colMin; + int nbit = ir * colSpanW + ic; + patt[nbit >> 3] |= (0x1 << (7 - (nbit % 8))); + } + patterns.emplace_back(static_cast(rowSpanW)); + patterns.emplace_back(static_cast(colSpanW)); + int nBytes = (rowSpanW * colSpanW + 7) / 8; + patterns.insert(patterns.end(), patt.begin(), patt.begin() + nBytes); + + Cluster cluster; + cluster.chipID = chipID; + cluster.row = bbox.rowMin; + cluster.col = bbox.colMin; + cluster.size = static_cast(pixbuf.size()); + cluster.subDetID = static_cast(subDetID); + cluster.layer = static_cast(layer); + cluster.disk = static_cast(disk); + clusters.emplace_back(cluster); +} + +//__________________________________________________ +void Clusterer::ClustererThread::fetchMCLabels(uint32_t digID, const ConstDigitTruth* labelsDig, int& nfilled) +{ + if (nfilled >= MaxLabels) { + return; + } + if (!labelsDig || digID >= labelsDig->getIndexedSize()) { + return; + } + const auto& lbls = labelsDig->getLabels(digID); + for (int i = lbls.size(); i--;) { + int ic = nfilled; + for (; ic--;) { + if (labelsBuff[ic] == lbls[i]) { + return; // already present + } + } + labelsBuff[nfilled++] = lbls[i]; + if (nfilled >= MaxLabels) { + break; + } + } +} + +} // namespace o2::trk diff --git a/Detectors/Upgrades/ALICE3/TRK/reconstruction/src/TRKReconstructionLinkDef.h b/Detectors/Upgrades/ALICE3/TRK/reconstruction/src/TRKReconstructionLinkDef.h index 09ab598ec626c..4eda22e350852 100644 --- a/Detectors/Upgrades/ALICE3/TRK/reconstruction/src/TRKReconstructionLinkDef.h +++ b/Detectors/Upgrades/ALICE3/TRK/reconstruction/src/TRKReconstructionLinkDef.h @@ -16,5 +16,6 @@ #pragma link off all functions; #pragma link C++ class o2::trk::TimeFrame < 11> + ; +#pragma link C++ class o2::trk::Clusterer + ; #endif diff --git a/Detectors/Upgrades/ALICE3/TRK/workflow/CMakeLists.txt b/Detectors/Upgrades/ALICE3/TRK/workflow/CMakeLists.txt index d6c8ea85c2bbd..42402fe6b62dc 100644 --- a/Detectors/Upgrades/ALICE3/TRK/workflow/CMakeLists.txt +++ b/Detectors/Upgrades/ALICE3/TRK/workflow/CMakeLists.txt @@ -13,6 +13,8 @@ o2_add_library(TRKWorkflow TARGETVARNAME targetName SOURCES src/DigitReaderSpec.cxx src/DigitWriterSpec.cxx + src/ClustererSpec.cxx + src/ClusterWriterSpec.cxx src/TrackerSpec.cxx src/TrackWriterSpec.cxx src/RecoWorkflow.cxx @@ -20,6 +22,7 @@ o2_add_library(TRKWorkflow O2::GPUWorkflow O2::SimConfig O2::DataFormatsITSMFT + O2::DataFormatsTRK O2::SimulationDataFormat O2::DPLUtils O2::TRKBase diff --git a/Detectors/Upgrades/ALICE3/TRK/workflow/include/TRKWorkflow/ClusterWriterSpec.h b/Detectors/Upgrades/ALICE3/TRK/workflow/include/TRKWorkflow/ClusterWriterSpec.h new file mode 100644 index 0000000000000..50d823b497bb9 --- /dev/null +++ b/Detectors/Upgrades/ALICE3/TRK/workflow/include/TRKWorkflow/ClusterWriterSpec.h @@ -0,0 +1,24 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +#ifndef O2_TRK_CLUSTERWRITER +#define O2_TRK_CLUSTERWRITER + +#include "Framework/DataProcessorSpec.h" + +namespace o2::trk +{ + +framework::DataProcessorSpec getClusterWriterSpec(bool useMC); + +} // namespace o2::trk + +#endif diff --git a/Detectors/Upgrades/ALICE3/TRK/workflow/include/TRKWorkflow/ClustererSpec.h b/Detectors/Upgrades/ALICE3/TRK/workflow/include/TRKWorkflow/ClustererSpec.h new file mode 100644 index 0000000000000..bacc1057c7b07 --- /dev/null +++ b/Detectors/Upgrades/ALICE3/TRK/workflow/include/TRKWorkflow/ClustererSpec.h @@ -0,0 +1,39 @@ +// Copyright 2019-2026 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +#ifndef O2_TRK_CLUSTERERDPL +#define O2_TRK_CLUSTERERDPL + +#include "Framework/DataProcessorSpec.h" +#include "Framework/Task.h" +#include "TRKReconstruction/Clusterer.h" + +namespace o2::trk +{ + +class ClustererDPL : public o2::framework::Task +{ + public: + ClustererDPL(bool useMC) : mUseMC(useMC) {} + void init(o2::framework::InitContext& ic) final; + void run(o2::framework::ProcessingContext& pc) final; + + private: + bool mUseMC = true; + int mNThreads = 1; + o2::trk::Clusterer mClusterer; +}; + +o2::framework::DataProcessorSpec getClustererSpec(bool useMC); + +} // namespace o2::trk + +#endif diff --git a/Detectors/Upgrades/ALICE3/TRK/workflow/src/ClusterWriterSpec.cxx b/Detectors/Upgrades/ALICE3/TRK/workflow/src/ClusterWriterSpec.cxx new file mode 100644 index 0000000000000..bc3a75c646198 --- /dev/null +++ b/Detectors/Upgrades/ALICE3/TRK/workflow/src/ClusterWriterSpec.cxx @@ -0,0 +1,65 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +/// @file ClusterWriterSpec.cxx + +#include + +#include "TRKWorkflow/ClusterWriterSpec.h" +#include "DPLUtils/MakeRootTreeWriterSpec.h" +#include "DataFormatsTRK/Cluster.h" +#include "DataFormatsTRK/ROFRecord.h" +#include "SimulationDataFormat/MCCompLabel.h" +#include "SimulationDataFormat/MCTruthContainer.h" + +using namespace o2::framework; + +namespace o2::trk +{ + +template +using BranchDefinition = MakeRootTreeWriterSpec::BranchDefinition; +using ClustersType = std::vector; +using PatternsType = std::vector; +using ROFrameType = std::vector; +using LabelsType = o2::dataformats::MCTruthContainer; +using ROFRecLblType = std::vector; + +DataProcessorSpec getClusterWriterSpec(bool useMC) +{ + auto clustersSize = std::make_shared(0); + auto clustersSizeGetter = [clustersSize](ClustersType const& clusters) { + *clustersSize = clusters.size(); + }; + auto logger = [clustersSize](ROFrameType const& rofs) { + LOG(info) << "TRKClusterWriter pulled " << *clustersSize << " clusters, in " << rofs.size() << " RO frames"; + }; + + return MakeRootTreeWriterSpec("trk-cluster-writer", + "o2clus_trk.root", + MakeRootTreeWriterSpec::TreeAttributes{"o2sim", "Tree with TRK clusters"}, + BranchDefinition{InputSpec{"compclus", "TRK", "COMPCLUSTERS", 0}, + "TRKClusterComp", + clustersSizeGetter}, + BranchDefinition{InputSpec{"patterns", "TRK", "PATTERNS", 0}, + "TRKClusterPatt"}, + BranchDefinition{InputSpec{"ROframes", "TRK", "CLUSTERSROF", 0}, + "TRKClustersROF", + logger}, + BranchDefinition{InputSpec{"labels", "TRK", "CLUSTERSMCTR", 0}, + "TRKClusterMCTruth", + (useMC ? 1 : 0)}, + BranchDefinition{InputSpec{"MC2ROframes", "TRK", "CLUSTERSMC2ROF", 0}, + "TRKClustersMC2ROF", + (useMC ? 1 : 0)})(); +} + +} // namespace o2::trk diff --git a/Detectors/Upgrades/ALICE3/TRK/workflow/src/ClustererSpec.cxx b/Detectors/Upgrades/ALICE3/TRK/workflow/src/ClustererSpec.cxx new file mode 100644 index 0000000000000..8aec63d69206b --- /dev/null +++ b/Detectors/Upgrades/ALICE3/TRK/workflow/src/ClustererSpec.cxx @@ -0,0 +1,99 @@ +// Copyright 2019-2026 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +#include "TRKWorkflow/ClustererSpec.h" +#include "DetectorsBase/GeometryManager.h" +#include "DataFormatsTRK/Cluster.h" +#include "DataFormatsTRK/ROFRecord.h" +#include "Framework/ConfigParamRegistry.h" +#include "Framework/Logger.h" +#include "SimulationDataFormat/ConstMCTruthContainer.h" + +namespace o2::trk +{ + +void ClustererDPL::init(o2::framework::InitContext& ic) +{ + mNThreads = std::max(1, ic.options().get("nthreads")); +} + +void ClustererDPL::run(o2::framework::ProcessingContext& pc) +{ + auto digits = pc.inputs().get>("digits"); + auto rofs = pc.inputs().get>("ROframes"); + + gsl::span mc2rofs; + gsl::span labelbuffer; + if (mUseMC) { + labelbuffer = pc.inputs().get>("labels"); + mc2rofs = pc.inputs().get>("MC2ROframes"); + } + o2::dataformats::ConstMCTruthContainerView labels(labelbuffer); + + std::vector clusters; + std::vector patterns; + std::vector clusterROFs; + std::unique_ptr> clusterLabels; + std::vector clusterMC2ROFs; + if (mUseMC) { + clusterLabels = std::make_unique>(); + } + o2::base::GeometryManager::loadGeometry("o2sim_geometry.root", false, true); + + mClusterer.process(digits, + rofs, + clusters, + patterns, + clusterROFs, + mUseMC ? &labels : nullptr, + clusterLabels.get(), + mc2rofs, + mUseMC ? &clusterMC2ROFs : nullptr); + + pc.outputs().snapshot(o2::framework::Output{"TRK", "COMPCLUSTERS", 0}, clusters); + pc.outputs().snapshot(o2::framework::Output{"TRK", "PATTERNS", 0}, patterns); + pc.outputs().snapshot(o2::framework::Output{"TRK", "CLUSTERSROF", 0}, clusterROFs); + + if (mUseMC) { + pc.outputs().snapshot(o2::framework::Output{"TRK", "CLUSTERSMCTR", 0}, *clusterLabels); + pc.outputs().snapshot(o2::framework::Output{"TRK", "CLUSTERSMC2ROF", 0}, clusterMC2ROFs); + } + + LOGP(info, "TRKClusterer pushed {} clusters in {} ROFs", clusters.size(), clusterROFs.size()); +} + +o2::framework::DataProcessorSpec getClustererSpec(bool useMC) +{ + std::vector inputs; + inputs.emplace_back("digits", "TRK", "DIGITS", 0, o2::framework::Lifetime::Timeframe); + inputs.emplace_back("ROframes", "TRK", "DIGITSROF", 0, o2::framework::Lifetime::Timeframe); + + std::vector outputs; + outputs.emplace_back("TRK", "COMPCLUSTERS", 0, o2::framework::Lifetime::Timeframe); + outputs.emplace_back("TRK", "PATTERNS", 0, o2::framework::Lifetime::Timeframe); + outputs.emplace_back("TRK", "CLUSTERSROF", 0, o2::framework::Lifetime::Timeframe); + + if (useMC) { + inputs.emplace_back("labels", "TRK", "DIGITSMCTR", 0, o2::framework::Lifetime::Timeframe); + inputs.emplace_back("MC2ROframes", "TRK", "DIGITSMC2ROF", 0, o2::framework::Lifetime::Timeframe); + outputs.emplace_back("TRK", "CLUSTERSMCTR", 0, o2::framework::Lifetime::Timeframe); + outputs.emplace_back("TRK", "CLUSTERSMC2ROF", 0, o2::framework::Lifetime::Timeframe); + } + + return o2::framework::DataProcessorSpec{ + "trk-clusterer", + inputs, + outputs, + o2::framework::AlgorithmSpec{o2::framework::adaptFromTask(useMC)}, + o2::framework::Options{{"nthreads", o2::framework::VariantType::Int, 1, {"Number of clustering threads"}}}}; +} + +} // namespace o2::trk diff --git a/Detectors/Upgrades/ALICE3/TRK/workflow/src/RecoWorkflow.cxx b/Detectors/Upgrades/ALICE3/TRK/workflow/src/RecoWorkflow.cxx index 5f6cbe2f96b04..d10feb4214f38 100644 --- a/Detectors/Upgrades/ALICE3/TRK/workflow/src/RecoWorkflow.cxx +++ b/Detectors/Upgrades/ALICE3/TRK/workflow/src/RecoWorkflow.cxx @@ -10,6 +10,9 @@ // or submit itself to any jurisdiction. #include "TRKWorkflow/RecoWorkflow.h" +#include "TRKWorkflow/ClustererSpec.h" +#include "TRKWorkflow/ClusterWriterSpec.h" +#include "TRKWorkflow/DigitReaderSpec.h" #include "TRKWorkflow/TrackerSpec.h" #include "TRKWorkflow/TrackWriterSpec.h" #include "Framework/CCDBParamSpec.h" @@ -28,10 +31,24 @@ framework::WorkflowSpec getWorkflow(bool useMC, o2::gpu::gpudatatypes::DeviceType dtype) { framework::WorkflowSpec specs; - specs.emplace_back(o2::trk::getTrackerSpec(useMC, hitRecoConfig, dtype)); + + if (!(upstreamDigits || upstreamClusters)) { + specs.emplace_back(o2::trk::getTRKDigitReaderSpec(useMC, false, "trkdigits.root")); + } + if (!upstreamClusters) { + specs.emplace_back(o2::trk::getClustererSpec(useMC)); + } if (!disableRootOutput) { - specs.emplace_back(o2::trk::getTrackWriterSpec(useMC)); + specs.emplace_back(o2::trk::getClusterWriterSpec(useMC)); + } + + if (!hitRecoConfig.empty()) { + LOGP(info, "Using hit reco config from file {}", hitRecoConfig); + specs.emplace_back(o2::trk::getTrackerSpec(useMC, hitRecoConfig, dtype)); + if (!disableRootOutput) { + specs.emplace_back(o2::trk::getTrackWriterSpec(useMC)); + } } return specs; diff --git a/Detectors/Upgrades/ALICE3/TRK/workflow/src/TrackerSpec.cxx b/Detectors/Upgrades/ALICE3/TRK/workflow/src/TrackerSpec.cxx index 8fc67f0fa5567..20bd45557dac5 100644 --- a/Detectors/Upgrades/ALICE3/TRK/workflow/src/TrackerSpec.cxx +++ b/Detectors/Upgrades/ALICE3/TRK/workflow/src/TrackerSpec.cxx @@ -396,9 +396,13 @@ DataProcessorSpec getTrackerSpec(bool useMC, const std::string& hitRecoConfig, o inputs.emplace_back("dummy", "TRK", "DUMMY", 0, Lifetime::Timeframe); - // inputs.emplace_back("compClusters", "TRK", "COMPCLUSTERS", 0, Lifetime::Timeframe); - // inputs.emplace_back("patterns", "TRK", "PATTERNS", 0, Lifetime::Timeframe); - // inputs.emplace_back("ROframes", "TRK", "CLUSTERSROF", 0, Lifetime::Timeframe); + constexpr bool expectClusterInputs = false; + if (expectClusterInputs) { + inputs.pop_back(); + inputs.emplace_back("compClusters", "TRK", "COMPCLUSTERS", 0, Lifetime::Timeframe); + inputs.emplace_back("patterns", "TRK", "PATTERNS", 0, Lifetime::Timeframe); + inputs.emplace_back("ROframes", "TRK", "CLUSTERSROF", 0, Lifetime::Timeframe); + } // inputs.emplace_back("itscldict", "TRK", "CLUSDICT", 0, Lifetime::Condition, ccdbParamSpec("ITS/Calib/ClusterDictionary")); // inputs.emplace_back("itsalppar", "TRK", "ALPIDEPARAM", 0, Lifetime::Condition, ccdbParamSpec("ITS/Config/AlpideParam"));