From 03745403d875d69cbed35336dde865899844657d Mon Sep 17 00:00:00 2001 From: Andrea Sofia Triolo Date: Sat, 7 Nov 2026 03:56:00 +0100 Subject: [PATCH 01/11] [ALICE3] VD: fixed digitization for the kIRIS4a and kIRISFullCyl geometries (#16063) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-authored-by: Andrea Sofia Triolo --- .../base/include/TRKBase/SegmentationChip.h | 34 +++++++++++++------ .../TRKFT3/TRK/base/include/TRKBase/Specs.h | 24 ++++++++----- .../TRKFT3/TRK/base/src/GeometryTGeo.cxx | 3 ++ .../common/macros/CheckTopologyDictionary.C | 3 +- .../common/reconstruction/src/Clusterer.cxx | 2 +- .../TRKFT3Simulation/ChipDigitsContainer.h | 2 +- .../include/TRKFT3Simulation/Digitizer.h | 2 +- 7 files changed, 47 insertions(+), 23 deletions(-) diff --git a/Detectors/Upgrades/ALICE3/TRKFT3/TRK/base/include/TRKBase/SegmentationChip.h b/Detectors/Upgrades/ALICE3/TRKFT3/TRK/base/include/TRKBase/SegmentationChip.h index fac5198966fd6..117a2a2c2c590 100644 --- a/Detectors/Upgrades/ALICE3/TRKFT3/TRK/base/include/TRKBase/SegmentationChip.h +++ b/Detectors/Upgrades/ALICE3/TRKFT3/TRK/base/include/TRKBase/SegmentationChip.h @@ -72,6 +72,20 @@ class SegmentationChip static constexpr std::array radiiVD = constants::VD::petal::layer::radii; + static inline int sNPetalsVD = -1; // not set + static void setNumberOfPetalsVD(int n) { sNPetalsVD = n; } + static float getWidthVD(int layer) // get the width of the petal according to the simulated geometry (supported: kIRIS4 (4 petals), kIRIS4a (3 petals), kIRISFullCyl) + { + if (sNPetalsVD < 1) { + LOGP(fatal, "SegmentationChip: VD number of petals not set. GeometryTGeo::Build() must be called before using the VD segmentation"); + } + if (sNPetalsVD == 1) { // FullCyl: 360°, no gap + return 2 * M_PI * constants::VD::petal::layer::radii[layer]; + } + return constants::VD::petal::layer::radii[layer] * 2 * M_PI / sNPetalsVD - constants::VD::petal::layer::gaps[layer]; + } + static int getNRowsVD(int layer) { return static_cast(getWidthVD(layer) / PitchRowVD); } + /// Transformation from Geant detector centered local coordinates (cm) to /// Pixel cell numbers iRow and iCol. /// Returns kTRUE if point x,z is inside sensitive volume, kFALSE otherwise. @@ -113,7 +127,7 @@ class SegmentationChip if (subDetID == 0) { pitchRow = PitchRowVD; pitchCol = PitchColVD; - maxWidth = constants::VD::petal::layer::width[layer]; + maxWidth = getWidthVD(layer); maxLength = constants::VD::petal::layer::length; // TODO: change this to use the layer and disk } else if (subDetID == 1 || subDetID == 2) { @@ -128,11 +142,11 @@ class SegmentationChip }; // Check local coordinates (cm) validity. - static constexpr bool isValidLoc(float x, float z, int subDetID, int layer) noexcept + static bool isValidLoc(float x, float z, int subDetID, int layer) noexcept { float maxWidth(0), maxLength(0); if (subDetID == 0) { - maxWidth = constants::VD::petal::layer::width[layer]; + maxWidth = getWidthVD(layer); maxLength = constants::VD::petal::layer::length; // TODO: change this to use the layer and disk } else if (subDetID == 1 || subDetID == 2) { // ML/OT @@ -143,12 +157,12 @@ class SegmentationChip } // Check detector coordinates validity. - static constexpr bool isValidDet(int row, int col, int subDetID, int layer) noexcept + static bool isValidDet(int row, int col, int subDetID, int layer) noexcept { // Check if the row and column are within the valid range int nRows(0), nCols(0); if (subDetID == 0) { - nRows = constants::VD::petal::layer::nRows[layer]; + nRows = getNRowsVD(layer); nCols = constants::VD::petal::layer::nCols; // TODO: change this to use the layer and disk } else if (subDetID == 1 || subDetID == 2) { @@ -171,7 +185,7 @@ class SegmentationChip /// \param int subDetID Sub-detector ID (0 for VD, 1 for ML/OT) /// \param int layer Layer number (0 to 2 for VD, 0 to 7 for ML/OT) /// \param int disk Disk number (0 to 5 for VD) - static constexpr bool detectorToLocal(int iRow, int iCol, float& xRow, float& zCol, int subDetID, int layer, int disk) noexcept + static bool detectorToLocal(int iRow, int iCol, float& xRow, float& zCol, int subDetID, int layer, int disk) noexcept { if (!isValidDet(iRow, iCol, subDetID, layer)) { LOGP(debug, "Detector coordinates not valid: iRow = {}, iCol = {}", iRow, iCol); @@ -194,7 +208,7 @@ class SegmentationChip /// xRow = half chip width - iRow(center) * pitch /// zCol = iCol * pitch - half chip lenght if (subDetID == 0) { - xRow = 0.5 * (constants::VD::petal::layer::width[layer] - PitchRowVD) - (row * PitchRowVD); + xRow = 0.5 * (getWidthVD(layer) - PitchRowVD) - (row * PitchRowVD); zCol = col * PitchColVD + 0.5 * (PitchColVD - constants::VD::petal::layer::length); } else if (subDetID == 1 || subDetID == 2) { // ML/OT xRow = 0.5 * (constants::moduleMLOT::chip::width - constants::moduleMLOT::chip::passiveEdgeReadOut - PitchRowMLOT) - (row * PitchRowMLOT); @@ -251,9 +265,9 @@ class SegmentationChip /// Print segmentation info static void Print() noexcept { - LOG(info) << "Number of rows:\nVD L0: " << constants::VD::petal::layer::nRows[0] - << "\nVD L1: " << constants::VD::petal::layer::nRows[1] - << "\nVD L2: " << constants::VD::petal::layer::nRows[2] + LOG(info) << "Number of rows:\nVD L0: " << getNRowsVD(0) + << "\nVD L1: " << getNRowsVD(1) + << "\nVD L2: " << getNRowsVD(2) << "\nML/OT chip: " << constants::moduleMLOT::chip::nRows; LOG(info) << "Number of cols:\nVD: " << constants::VD::petal::layer::nCols diff --git a/Detectors/Upgrades/ALICE3/TRKFT3/TRK/base/include/TRKBase/Specs.h b/Detectors/Upgrades/ALICE3/TRKFT3/TRK/base/include/TRKBase/Specs.h index 9ac1f789d43ba..866178f820231 100644 --- a/Detectors/Upgrades/ALICE3/TRKFT3/TRK/base/include/TRKBase/Specs.h +++ b/Detectors/Upgrades/ALICE3/TRKFT3/TRK/base/include/TRKBase/Specs.h @@ -48,15 +48,13 @@ constexpr int nLayers{3}; // number of layers in each VD petal constexpr int nDisks{6}; // number of disks in each VD petal namespace layer { -constexpr double pitchX{10 * mu}; // pitch of the row -constexpr double pitchZ{10 * mu}; // pitch of the column -constexpr double totalThickness{silicon::thickness + metalstack::thickness}; // total thickness of the chip -constexpr std::array gaps{1.63 * mm, 1.2 * mm, 1.2 * mm}; // gaps between two consecutive petals -constexpr std::array radii{0.5 * cm, 1.2 * cm, 2.5 * cm}; // radius of layer in cm -constexpr std::array width{radii[0] * 2 * M_PI / 4 - gaps[0], radii[1] * 2 * M_PI / 4 - gaps[1], radii[2] * 2 * M_PI / 4 - gaps[2]}; // width of the quarter of layer in cm -constexpr double length{50 * cm}; // length of the layer -constexpr int nCols{static_cast(length / pitchZ)}; // number of columns in the chip -constexpr std::array nRows{static_cast(width[0] / pitchX), static_cast(width[1] / pitchX), static_cast(width[2] / pitchX)}; // number of rows in the chip. For the moment is different for each layer since a siner segmentation in repetitive units is stil to be implemented +constexpr double pitchX{10 * mu}; // pitch of the row +constexpr double pitchZ{10 * mu}; // pitch of the column +constexpr double totalThickness{silicon::thickness + metalstack::thickness}; // total thickness of the chip +constexpr std::array gaps{1.63 * mm, 1.2 * mm, 1.2 * mm}; // gaps between two consecutive petals +constexpr std::array radii{0.5 * cm, 1.2 * cm, 2.5 * cm}; // radius of layer in cm +constexpr double length{50 * cm}; // length of the layer +constexpr int nCols{static_cast(length / pitchZ)}; // number of columns in the chip } // namespace layer namespace disk @@ -223,6 +221,14 @@ constexpr double responseYShift{5 * mu}; /// center of the epitaxial layer constexpr double thickness{20 * mu}; } // namespace alice3resp +namespace alpideresp /// parameters for the alpide chip response +{ +constexpr double pitchX{26.88 * mu}; +constexpr double pitchZ{29.24 * mu}; +constexpr double responseYShift{5 * mu}; /// center of the epitaxial layer +constexpr double thickness{45 * mu}; +} // namespace alpideresp + namespace MLOTDisks { constexpr int nLayers{12}; // number of disks in the ML and OT diff --git a/Detectors/Upgrades/ALICE3/TRKFT3/TRK/base/src/GeometryTGeo.cxx b/Detectors/Upgrades/ALICE3/TRKFT3/TRK/base/src/GeometryTGeo.cxx index ddfc844cc964d..699b65dd31df6 100644 --- a/Detectors/Upgrades/ALICE3/TRKFT3/TRK/base/src/GeometryTGeo.cxx +++ b/Detectors/Upgrades/ALICE3/TRKFT3/TRK/base/src/GeometryTGeo.cxx @@ -84,6 +84,9 @@ void GeometryTGeo::Build(int loadTrans) mNumberOfLayersMLOT = extractNumberOfLayersMLOT(); mNumberOfPetalsVD = extractNumberOfPetalsVD(); + SegmentationChip::setNumberOfPetalsVD(mNumberOfPetalsVD); + LOGP(info, "VD segmentation set for {} petals", mNumberOfPetalsVD); + mNumberOfActivePartsVD = extractNumberOfActivePartsVD(); mNumberOfLayersVD = extractNumberOfLayersVD(); mNumberOfDisksVD = extractNumberOfDisksVD(); diff --git a/Detectors/Upgrades/ALICE3/TRKFT3/common/macros/CheckTopologyDictionary.C b/Detectors/Upgrades/ALICE3/TRKFT3/common/macros/CheckTopologyDictionary.C index f97a586373308..04e4a5ab686da 100644 --- a/Detectors/Upgrades/ALICE3/TRKFT3/common/macros/CheckTopologyDictionary.C +++ b/Detectors/Upgrades/ALICE3/TRKFT3/common/macros/CheckTopologyDictionary.C @@ -39,6 +39,7 @@ #include "SimulationDataFormat/MCCompLabel.h" #include "SimulationDataFormat/MCTruthContainer.h" #include "TRKBase/GeometryTGeo.h" +#include "TRKBase/SegmentationChip.h" #include "TRKBase/Specs.h" #include "TRKReconstruction/BuildTopologyDictionary.h" #include @@ -190,7 +191,7 @@ void forEachCluster(const std::string& directory, F&& f) const auto* pattern = candidate.pattern; // Left out: at the matrix's edge, or one of several pieces of a crossing. const bool vd = cluster.subDetID == 0; - const int nRows = vd ? o2::trk::constants::VD::petal::layer::nRows[cluster.layer] : o2::trk::constants::moduleMLOT::chip::nRows; + const int nRows = vd ? o2::trk::SegmentationChip::getNRowsVD(cluster.layer) : o2::trk::constants::moduleMLOT::chip::nRows; const int nCols = vd ? o2::trk::constants::VD::petal::layer::nCols : o2::trk::constants::moduleMLOT::chip::nCols; if (cluster.row == 0 || cluster.col == 0 || cluster.row + pattern[0] >= nRows || cluster.col + pattern[1] >= nCols || nClusters[candidate.key] > static_cast(candidate.segments.size())) { diff --git a/Detectors/Upgrades/ALICE3/TRKFT3/common/reconstruction/src/Clusterer.cxx b/Detectors/Upgrades/ALICE3/TRKFT3/common/reconstruction/src/Clusterer.cxx index 03e7e5285a31f..8e37df7bfe10f 100644 --- a/Detectors/Upgrades/ALICE3/TRKFT3/common/reconstruction/src/Clusterer.cxx +++ b/Detectors/Upgrades/ALICE3/TRKFT3/common/reconstruction/src/Clusterer.cxx @@ -181,7 +181,7 @@ void Clusterer::ClustererThread::initChip(gsl::span digits, if (geom) { if (geom->getSubDetID(chipID) == 0) { // VD const int layer = geom->getLayer(chipID); - size = constants::VD::petal::layer::nRows[layer] + 2; + size = SegmentationChip::getNRowsVD(layer) + 2; } } diff --git a/Detectors/Upgrades/ALICE3/TRKFT3/common/simulation/include/TRKFT3Simulation/ChipDigitsContainer.h b/Detectors/Upgrades/ALICE3/TRKFT3/common/simulation/include/TRKFT3Simulation/ChipDigitsContainer.h index 0fcbfd2591a0e..5176e8cbbcfe7 100644 --- a/Detectors/Upgrades/ALICE3/TRKFT3/common/simulation/include/TRKFT3Simulation/ChipDigitsContainer.h +++ b/Detectors/Upgrades/ALICE3/TRKFT3/common/simulation/include/TRKFT3Simulation/ChipDigitsContainer.h @@ -96,7 +96,7 @@ void o2::trkft3::ChipDigitsContainer::addNoise(UInt_t rofMin, UInt_t rofMax, con int maxCols = 0; if (subDetID == 0) { - maxRows = o2::trk::constants::VD::petal::layer::nRows[layer]; + maxRows = o2::trk::SegmentationChip::getNRowsVD(layer); maxCols = o2::trk::constants::VD::petal::layer::nCols; } else { maxRows = o2::trk::constants::moduleMLOT::chip::nRows; diff --git a/Detectors/Upgrades/ALICE3/TRKFT3/common/simulation/include/TRKFT3Simulation/Digitizer.h b/Detectors/Upgrades/ALICE3/TRKFT3/common/simulation/include/TRKFT3Simulation/Digitizer.h index d0cf085d971a0..65629cf7de6fc 100644 --- a/Detectors/Upgrades/ALICE3/TRKFT3/common/simulation/include/TRKFT3Simulation/Digitizer.h +++ b/Detectors/Upgrades/ALICE3/TRKFT3/common/simulation/include/TRKFT3Simulation/Digitizer.h @@ -130,7 +130,7 @@ class Digitizer int getNRows(int subDetID, int layer) { if (subDetID == 0) { // VD - return o2::trk::constants::VD::petal::layer::nRows[layer]; + return o2::trk::SegmentationChip::getNRowsVD(layer); } else if (subDetID == 1 || subDetID == 2) { // ML/OT return o2::trk::constants::moduleMLOT::chip::nRows; } From 9d6f193590a61a8e870acc7a16cd6e5d922bc4ef Mon Sep 17 00:00:00 2001 From: Maximiliano Puccio Date: Tue, 10 Nov 2026 14:42:01 +0100 Subject: [PATCH 02/11] ITSMFT: fix MC labels for single pixel clusters in the labels --- Detectors/ITSMFT/common/reconstruction/src/Clusterer.cxx | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/Detectors/ITSMFT/common/reconstruction/src/Clusterer.cxx b/Detectors/ITSMFT/common/reconstruction/src/Clusterer.cxx index 2168c4b8d3308..8f872361d35d8 100644 --- a/Detectors/ITSMFT/common/reconstruction/src/Clusterer.cxx +++ b/Detectors/ITSMFT/common/reconstruction/src/Clusterer.cxx @@ -315,7 +315,10 @@ void Clusterer::ClustererThread::finishChipSingleHitFast(uint32_t hit, ChipPixel if (labelsClusPtr) { // MC labels were requested int nlab = 0; - fetchMCLabels(curChipData->getStartID() + hit, labelsDigPtr, nlab); + // Squashing can remove earlier digits on the chip, so the surviving + // pixel's source index need not be contiguous with the chip's start. + const auto digitID = parent->mSquashingDepth ? curChipData->getOrderedPixId(hit) : curChipData->getStartID() + hit; + fetchMCLabels(digitID, labelsDigPtr, nlab); auto cnt = compClusPtr->size(); for (int i = 0; i < nlab; i++) { labelsClusPtr->addElement(cnt, labelsBuff[i]); From b9e9eddab17b0d4e4c5441dbf75b3e9a90c9eb09 Mon Sep 17 00:00:00 2001 From: Maximiliano Puccio Date: Mon, 9 Nov 2026 14:01:45 +0100 Subject: [PATCH 03/11] GenericTracking: Fix compilation for VD --- Detectors/GenericTracking/tools/src/importGeometry.cxx | 2 +- Detectors/GenericTracking/workflow/src/TRK.cxx | 5 ++--- 2 files changed, 3 insertions(+), 4 deletions(-) diff --git a/Detectors/GenericTracking/tools/src/importGeometry.cxx b/Detectors/GenericTracking/tools/src/importGeometry.cxx index e7323dcbda0dc..a357acedbd0c4 100644 --- a/Detectors/GenericTracking/tools/src/importGeometry.cxx +++ b/Detectors/GenericTracking/tools/src/importGeometry.cxx @@ -95,7 +95,7 @@ double checkTRKPixels(std::span sensors) const int subDet = disc ? ft3.getSubDetID(id) : geometry.getSubDetID(id); const int layer = disc ? ft3.getLayer(id) : geometry.getLayer(id); // within the sub-detector const auto& spec = Geometry::layerSpec(sensors[chip].layer); - const int nRows = subDet == 0 ? o2::trk::constants::VD::petal::layer::nRows[static_cast(layer)] : o2::trk::constants::moduleMLOT::chip::nRows; + const int nRows = subDet == 0 ? Segmentation::getNRowsVD(layer) : o2::trk::constants::moduleMLOT::chip::nRows; const int nCols = subDet == 0 ? o2::trk::constants::VD::petal::layer::nCols : o2::trk::constants::moduleMLOT::chip::nCols; for (const auto [row, col] : std::array, 5>{{{0, 0}, {nRows - 1, 0}, {0, nCols - 1}, {nRows - 1, nCols - 1}, {nRows / 2, nCols / 2}}}) { float xRow = 0.F; diff --git a/Detectors/GenericTracking/workflow/src/TRK.cxx b/Detectors/GenericTracking/workflow/src/TRK.cxx index a6aeb78294e95..d87d91c31d052 100644 --- a/Detectors/GenericTracking/workflow/src/TRK.cxx +++ b/Detectors/GenericTracking/workflow/src/TRK.cxx @@ -69,11 +69,10 @@ namespace // TRKFT3.h restates O2's TRK sensors; here they are held to O2's, to the bit. namespace specs = o2::trk::constants; +// The four-petal VD widths are derived from these radii and gaps. O2 now +// computes widths at run time from the loaded geometry's petal count. static_assert(detectors::trkft3::kVDRadius == specs::VD::petal::layer::radii); static_assert(detectors::trkft3::kVDPetalGap == specs::VD::petal::layer::gaps); -static_assert(exactlyEqual(detectors::trkft3::vdWidth(0), specs::VD::petal::layer::width[0])); -static_assert(exactlyEqual(detectors::trkft3::vdWidth(1), specs::VD::petal::layer::width[1])); -static_assert(exactlyEqual(detectors::trkft3::vdWidth(2), specs::VD::petal::layer::width[2])); static_assert(exactlyEqual(detectors::trkft3::kVDLength, specs::VD::petal::layer::length)); static_assert(exactlyEqual(detectors::trkft3::kVDThickness, specs::VD::silicon::thickness)); static_assert(exactlyEqual(detectors::trkft3::kVDPitch, specs::VD::petal::layer::pitchX) && exactlyEqual(detectors::trkft3::kVDPitch, specs::VD::petal::layer::pitchZ)); From 6401e3994262f55b37bbf5582b00f7d83aaa5e9d Mon Sep 17 00:00:00 2001 From: Maximiliano Puccio Date: Tue, 10 Nov 2026 15:12:19 +0100 Subject: [PATCH 04/11] GenericTracking: support staggered ITS and MFT cluster streams Read per-layer clusters, labels and ROFs in a shared dense ordering while preserving composed O2 cluster references in the output adapters. Propagate staggering options through reconstruction and study workflows, use per-layer frame timing and select the fastest output clock. Decode native per-layer cluster references for common-truth QA. Initialize ALICE3 converter geometry before using the curved VD segmentation. Validation: build ITS/MFT reconstruction and study workflows and TRK study; retain the completed uniform/staggered MC device-study comparison evidence. --- .../study/src/MeasurementSinkSpec.cxx | 31 ++- .../study/tools/StudyWorkflow.h | 4 +- .../study/tools/gtrk-its-study-workflow.cxx | 2 +- .../study/tools/gtrk-mft-study-workflow.cxx | 2 +- .../GenericTracking/Workflow/AdaptorSpec.h | 4 +- .../Workflow/ClusterConverterSpec.h | 5 +- .../GenericTracking/Workflow/Conditions.h | 22 +- .../GenericTracking/Workflow/ITSMFTOptions.h | 39 ++++ .../workflow/src/AdaptorDevice.h | 213 ++++++------------ .../workflow/src/AdaptorHelpers.h | 199 ++++++++++++++++ .../workflow/src/AdaptorITS.cxx | 2 +- .../workflow/src/AdaptorMFT.cxx | 2 +- .../workflow/src/ClusterConverterITSMFT.h | 58 +++-- .../GenericTracking/workflow/src/ITS.cxx | 4 +- .../workflow/src/ITSMFTInput.h | 181 +++++++++++++++ .../GenericTracking/workflow/src/MFT.cxx | 4 +- .../workflow/src/RecoWorkflow.h | 4 +- .../GenericTracking/workflow/src/TRK.cxx | 24 +- .../workflow/src/gtrk-its-reco-workflow.cxx | 5 +- .../workflow/src/gtrk-mft-reco-workflow.cxx | 5 +- 20 files changed, 605 insertions(+), 205 deletions(-) create mode 100644 Detectors/GenericTracking/workflow/include/GenericTracking/Workflow/ITSMFTOptions.h create mode 100644 Detectors/GenericTracking/workflow/src/AdaptorHelpers.h create mode 100644 Detectors/GenericTracking/workflow/src/ITSMFTInput.h diff --git a/Detectors/GenericTracking/study/src/MeasurementSinkSpec.cxx b/Detectors/GenericTracking/study/src/MeasurementSinkSpec.cxx index 48088f1403cf7..fffb992e06006 100644 --- a/Detectors/GenericTracking/study/src/MeasurementSinkSpec.cxx +++ b/Detectors/GenericTracking/study/src/MeasurementSinkSpec.cxx @@ -9,6 +9,7 @@ // granted to it by virtue of its status as an Intergovernmental Organization // or submit itself to any jurisdiction. +#include "DataFormatsITSMFT/ClusterID.h" #include "GenericTracking/Study/MCVertexReaderSpec.h" #include "GenericTracking/Study/Digest.h" #include "GenericTracking/Study/MeasurementSinkSpec.h" @@ -3638,7 +3639,7 @@ class MeasurementSinkTask : public Task tracksOfRow[view.clusters[track].index[i]].push_back(static_cast(track)); } } - const auto theirs = mO2Tracks->read(static_cast(mTimeframe)); + const auto theirs = o2Clusters(frame); constexpr auto kMaxClusters = static_cast(Detector::nLayers) + 1U; std::array o2Total{}; std::array o2Found{}; @@ -3842,12 +3843,36 @@ class MeasurementSinkTask : public Task return rowOfCluster; } + /// Decode O2's per-layer references to the converter's dense cluster order. + /// Monolithic references (stream zero) already are dense indices. + [[nodiscard]] TrackClusters o2Clusters(const Frame& frame) const + { + auto tracks = mO2Tracks->read(static_cast(mTimeframe)); + std::array::nLayers + 1U> offsets{}; + for (Frame::Index row = 0U; row < frame.size(); ++row) { + ++offsets[mInterface.geometry().sensor(frame.measurement(row).surface).layer + 1U]; + } + std::partial_sum(offsets.begin(), offsets.end(), offsets.begin()); + for (auto& track : tracks) { + for (auto& cluster : track) { + const auto layer = static_cast(o2::itsmft::clusID2Layer(static_cast(cluster))); + if (layer == 0U) { + continue; + } + const auto index = static_cast(o2::itsmft::clusID2Index(static_cast(cluster))); + cluster = layer + 1U < offsets.size() && index < offsets[layer + 1U] - offsets[layer] + ? offsets[layer] + index : std::numeric_limits::max(); + } + } + return tracks; + } + /// O2's tracks of this timeframe as our measurement rows: O2's tracks name - /// O2's cluster indices, which are our measurements' clusterIndex. + /// cluster references, decoded to our measurements' clusterIndex. [[nodiscard]] TrackRows o2TrackRows(const Frame& frame) const { const auto rowOfCluster = rowsOfClusters(frame); - const auto clusters = mO2Tracks->read(static_cast(mTimeframe)); + const auto clusters = o2Clusters(frame); TrackRows theirs; theirs.reserve(clusters.size()); for (const auto& track : clusters) { diff --git a/Detectors/GenericTracking/study/tools/StudyWorkflow.h b/Detectors/GenericTracking/study/tools/StudyWorkflow.h index 3d870b820f4c5..e5584c393dba6 100644 --- a/Detectors/GenericTracking/study/tools/StudyWorkflow.h +++ b/Detectors/GenericTracking/study/tools/StudyWorkflow.h @@ -37,6 +37,7 @@ #include "GenericTracking/Study/MCVertexReaderSpec.h" #include "GenericTracking/Study/MeasurementSinkSpec.h" #include "GenericTracking/Workflow/ClusterConverterSpec.h" +#include "GenericTracking/Workflow/ITSMFTOptions.h" #include "GenericTracking/Workflow/MeasurementProducerSpec.h" #include "GenericTracking/Workflow/TrackerSpec.h" @@ -54,6 +55,7 @@ namespace o2::gtrk::study /// Every study workflow's options. inline void addStudyOptions(std::vector& options) { + o2::itsmft::DPLAlpideParamInitializer::addConfigOption(options); using o2::framework::ConfigParamSpec; using o2::framework::VariantType; options.push_back(ConfigParamSpec{"disable-mc", VariantType::Bool, false, {"do not propagate MC labels"}}); @@ -86,7 +88,7 @@ template if (!context.options().get("clusters-from-upstream")) { specs.insert(specs.end(), readers.begin(), readers.end()); } - specs.push_back(gw::getClusterConverterSpec(useMC)); + specs.push_back(gw::getClusterConverterSpec(useMC, gw::staggeredInput(context))); specs.push_back(gw::getMeasurementProducerSpec()); specs.push_back(getMeasurementSinkSpec(useMC, mcVertices, products)); if (products != TrackerProducts::kNone) { diff --git a/Detectors/GenericTracking/study/tools/gtrk-its-study-workflow.cxx b/Detectors/GenericTracking/study/tools/gtrk-its-study-workflow.cxx index b249b9d84feb0..3593854de694a 100644 --- a/Detectors/GenericTracking/study/tools/gtrk-its-study-workflow.cxx +++ b/Detectors/GenericTracking/study/tools/gtrk-its-study-workflow.cxx @@ -40,7 +40,7 @@ WorkflowSpec defineDataProcessing(const ConfigContext& context) o2::conf::ConfigurableParam::updateFromString(context.options().get("configKeyValues")); const bool useMC = !context.options().get("disable-mc"); - auto specs = o2::gtrk::study::studyChain(context, {o2::itsmft::getITSClusterReaderSpec(useMC, false, true, false)}); + auto specs = o2::gtrk::study::studyChain(context, {o2::itsmft::getITSClusterReaderSpec(useMC, o2::gtrk::workflow::staggeredInput(context), true, false)}); const o2::raw::HBFUtilsInitializer hbfIni(context, specs); diff --git a/Detectors/GenericTracking/study/tools/gtrk-mft-study-workflow.cxx b/Detectors/GenericTracking/study/tools/gtrk-mft-study-workflow.cxx index 2a5e060e598d8..e1aebcc237470 100644 --- a/Detectors/GenericTracking/study/tools/gtrk-mft-study-workflow.cxx +++ b/Detectors/GenericTracking/study/tools/gtrk-mft-study-workflow.cxx @@ -40,7 +40,7 @@ WorkflowSpec defineDataProcessing(const ConfigContext& context) o2::conf::ConfigurableParam::updateFromString(context.options().get("configKeyValues")); const bool useMC = !context.options().get("disable-mc"); - auto specs = o2::gtrk::study::studyChain(context, {o2::itsmft::getMFTClusterReaderSpec(useMC, false, true, false)}); + auto specs = o2::gtrk::study::studyChain(context, {o2::itsmft::getMFTClusterReaderSpec(useMC, o2::gtrk::workflow::staggeredInput(context), true, false)}); const o2::raw::HBFUtilsInitializer hbfIni(context, specs); diff --git a/Detectors/GenericTracking/workflow/include/GenericTracking/Workflow/AdaptorSpec.h b/Detectors/GenericTracking/workflow/include/GenericTracking/Workflow/AdaptorSpec.h index e3a0602f8bc90..2340881c6f6d4 100644 --- a/Detectors/GenericTracking/workflow/include/GenericTracking/Workflow/AdaptorSpec.h +++ b/Detectors/GenericTracking/workflow/include/GenericTracking/Workflow/AdaptorSpec.h @@ -25,8 +25,10 @@ namespace o2::gtrk::workflow /// The adaptor of \a Det (the ITS or the MFT); with \a useMC the track labels /// from the clusters' (and, for the ITS, the vertices' outputs, without truth). +/// With staggered input, subscribe to each layer and publish composed O2 cluster +/// references; GTRK measurement indices remain dense in layer-stream order. template -o2::framework::DataProcessorSpec getAdaptorSpec(bool useMC); +o2::framework::DataProcessorSpec getAdaptorSpec(bool useMC, bool staggered = false); } // namespace o2::gtrk::workflow diff --git a/Detectors/GenericTracking/workflow/include/GenericTracking/Workflow/ClusterConverterSpec.h b/Detectors/GenericTracking/workflow/include/GenericTracking/Workflow/ClusterConverterSpec.h index e3d64259b0687..19b30553bca53 100644 --- a/Detectors/GenericTracking/workflow/include/GenericTracking/Workflow/ClusterConverterSpec.h +++ b/Detectors/GenericTracking/workflow/include/GenericTracking/Workflow/ClusterConverterSpec.h @@ -24,10 +24,9 @@ namespace o2::gtrk::workflow { /// Defined in the detector's own unit (ITS.cxx, MFT.cxx, TRK.cxx). The ITS's -/// and MFT's readout frames must not be staggered: only subspecification 0 is -/// read. +/// and MFT's staggered input has one stream per layer, concatenated in layer order. template -o2::framework::DataProcessorSpec getClusterConverterSpec(bool useMC); +o2::framework::DataProcessorSpec getClusterConverterSpec(bool useMC, bool staggered = false); } // namespace o2::gtrk::workflow diff --git a/Detectors/GenericTracking/workflow/include/GenericTracking/Workflow/Conditions.h b/Detectors/GenericTracking/workflow/include/GenericTracking/Workflow/Conditions.h index 244febffe34d3..667d32a25b5bf 100644 --- a/Detectors/GenericTracking/workflow/include/GenericTracking/Workflow/Conditions.h +++ b/Detectors/GenericTracking/workflow/include/GenericTracking/Workflow/Conditions.h @@ -157,24 +157,20 @@ class Conditions return false; } - /// The readout frames of the run's `DPLAlpideParam`, as O2's ITS tracker - /// takes them (`ITSTrackingInterface::updateTimeDependentParams`): frame `k` - /// at `k * length + delay + bias` BC, one per layer, every layer the same. - /// Throws std::runtime_error for staggered frames: the converter reads one - /// ROF stream. + /// Frame k starts at k * length + delay + bias BC, independently per layer. [[nodiscard]] std::vector frameTiming() const { const auto& param = AlpideParam::Instance(); - const FrameTiming result{.lengthBC = param.getROFLengthInBC(0), .firstStartBC = param.getROFDelayInBC(0) + param.getROFBiasInBC(0)}; - for (int layer = 1; layer < AlpideParam::getNLayers(); ++layer) { - if (param.getROFLengthInBC(layer) != result.lengthBC || param.getROFDelayInBC(layer) + param.getROFBiasInBC(layer) != result.firstStartBC) { - throw std::runtime_error(std::format("{}'s readout frames are staggered (layer {}): not read yet", Det::name, layer)); + std::vector result; + result.reserve(Detector::nLayers); + for (int layer = 0; layer < Detector::nLayers; ++layer) { + const auto length = param.getROFLengthInBC(layer); + if (length <= 0) { + throw std::runtime_error(std::format("{} layer {}: AlpideParam gives readout frames of {} BC", Det::name, layer, length)); } + result.push_back({.lengthBC = length, .firstStartBC = param.getROFDelayInBC(layer) + param.getROFBiasInBC(layer)}); } - if (result.lengthBC <= 0) { - throw std::runtime_error(std::format("{}'s AlpideParam gives readout frames of {} BC", Det::name, result.lengthBC)); - } - return std::vector(static_cast(Detector::nLayers), result); + return result; } /// The run's orbits per timeframe, its GRPECS's, as O2's ITS tracking diff --git a/Detectors/GenericTracking/workflow/include/GenericTracking/Workflow/ITSMFTOptions.h b/Detectors/GenericTracking/workflow/include/GenericTracking/Workflow/ITSMFTOptions.h new file mode 100644 index 0000000000000..4412af16c03cc --- /dev/null +++ b/Detectors/GenericTracking/workflow/include/GenericTracking/Workflow/ITSMFTOptions.h @@ -0,0 +1,39 @@ +// 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 GTRK_WORKFLOW_ITSMFTOPTIONS_H_ +#define GTRK_WORKFLOW_ITSMFTOPTIONS_H_ + +#include "GenericTracking/Detectors/ITS.h" +#include "GenericTracking/Detectors/MFT.h" +#include "DataFormatsITSMFT/DPLAlpideParamInitializer.h" +#include "Framework/ConfigContext.h" + +#include + +namespace o2::gtrk::workflow +{ + +template +[[nodiscard]] bool staggeredInput(const o2::framework::ConfigContext& context) +{ + if constexpr (std::same_as) { + return o2::itsmft::DPLAlpideParamInitializer::isITSStaggeringEnabled(context); + } else if constexpr (std::same_as) { + return o2::itsmft::DPLAlpideParamInitializer::isMFTStaggeringEnabled(context); + } else { + return false; + } +} + +} // namespace o2::gtrk::workflow + +#endif // GTRK_WORKFLOW_ITSMFTOPTIONS_H_ diff --git a/Detectors/GenericTracking/workflow/src/AdaptorDevice.h b/Detectors/GenericTracking/workflow/src/AdaptorDevice.h index b0452a5533e2a..832af96fa070c 100644 --- a/Detectors/GenericTracking/workflow/src/AdaptorDevice.h +++ b/Detectors/GenericTracking/workflow/src/AdaptorDevice.h @@ -23,6 +23,8 @@ #ifndef GTRK_WORKFLOW_ADAPTORDEVICE_H_ #define GTRK_WORKFLOW_ADAPTORDEVICE_H_ +#include "AdaptorHelpers.h" + #include "GenericTracking/Alice/Clusters.h" #include "GenericTracking/Detectors/ITS.h" #include "GenericTracking/Detectors/MFT.h" @@ -80,101 +82,6 @@ namespace o2::gtrk::workflow namespace adaptor { -/// A particle's label, the fake flag unless \a good; O2's noise label for no -/// particle. -[[nodiscard]] inline o2::MCCompLabel labelOf(std::uint64_t particle, bool good) { return particle == kNoTrackParticle ? o2::MCCompLabel(true) : labelOfParticle(particle, !good); } - -/// A timeframe's clusters as the adaptor needs them: each one's readout frame -/// (its record in CLUSTERSROF) and that frame's start in BC from the -/// timeframe's start, without the delay and bias. -struct ClusterFrames { - std::vector record; ///< per cluster - std::vector recordBC; ///< per record -}; - -[[nodiscard]] inline ClusterFrames clusterFrames(std::span rofs, std::size_t clusters, std::uint32_t tfFirstOrbit) -{ - ClusterFrames result; - result.record.assign(clusters, std::numeric_limits::max()); - result.recordBC.reserve(rofs.size()); - for (std::size_t r = 0; r < rofs.size(); ++r) { - const auto& ir = rofs[r].getBCData(); - result.recordBC.push_back((static_cast(ir.orbit - tfFirstOrbit) * o2::constants::lhc::LHCMaxBunches) + ir.bc); - for (int i = 0; i < rofs[r].getNEntries(); ++i) { - result.record[static_cast(rofs[r].getFirstEntry()) + static_cast(i)] = static_cast(r); - } - } - return result; -} - -/// A track's time as ITS takes it (`TrackerTraits`): the intersection of its -/// clusters' frames `[start + delay + bias, + length)` in BC from the -/// timeframe's start; where they do not intersect -- a track across frames, -/// which only a margin lets the tracker make -- the record's box, the frames -/// widened by the margin. [lower, upper) in BC. -struct BCRange { - std::int64_t lower = 0; - std::int64_t upper = 0; -}; - -template -[[nodiscard]] BCRange trackTime(const Record& record, std::span clusterIndex, const ClusterFrames& frames, const FrameTiming& timing) -{ - BCRange nominal{.lower = std::numeric_limits::min(), .upper = std::numeric_limits::max()}; - for (std::uint8_t i = 0; i < record.count; ++i) { - const auto start = frames.recordBC[frames.record[clusterIndex[record.clusters.index[i]]]] + timing.firstStartBC; - nominal.lower = std::max(nominal.lower, start); - nominal.upper = std::min(nominal.upper, start + timing.lengthBC); - } - if (nominal.lower < nominal.upper) { - return nominal; - } - return {.lower = bcOfPS(record.timePS), .upper = bcOfPS(static_cast(record.timePS) + record.widthPS)}; -} - -/// The earliest frame (its CLUSTERSROF record) among a track's clusters. -template -[[nodiscard]] std::uint32_t earliestRecord(const Record& record, std::span clusterIndex, const ClusterFrames& frames) -{ - std::uint32_t earliest = std::numeric_limits::max(); - for (std::uint8_t i = 0; i < record.count; ++i) { - earliest = std::min(earliest, frames.record[clusterIndex[record.clusters.index[i]]]); - } - return earliest; -} - -/// Each cluster's pixel count, the converter's decoding (`localClusters`). -[[nodiscard]] inline std::vector pixelCounts(std::span compact, std::span patterns, const o2::itsmft::TopologyDictionary& dictionary) -{ - const auto local = alice::localClusters(compact, patterns, dictionary); - std::vector result(local.size()); - std::ranges::transform(local, result.begin(), [](const auto& cluster) { return cluster.nPixels; }); - return result; -} - -/// The track's particle by the shared rule (TrackLabel.h), its rows' -/// clusters' labels, innermost first; \a carried gets, per row, whether it -/// carries that particle. -template -[[nodiscard]] o2::MCCompLabel trackLabel(const Record& record, std::span clusterIndex, const o2::dataformats::MCTruthContainer& labels, std::span carried) -{ - constexpr std::size_t kRows = std::extent_v>; - std::array rows{}; - std::array track{}; - for (std::uint8_t i = 0; i < record.count; ++i) { - const auto all = labels.getLabels(clusterIndex[record.clusters.index[i]]); - rows[i] = rowParticles(std::span{all.data(), all.size()}); - track[i] = i; - } - const std::span order{track.data(), record.count}; - const auto majority = majorityParticle(rows, order); - const bool good = majority.particle != kNoTrackParticle && majority.carried == record.count; - for (std::uint8_t i = 0; i < record.count; ++i) { - carried[i] = majority.particle != kNoTrackParticle && rows[i].carries(majority.particle); - } - return labelOf(majority.particle, good); -} - /// The description of O2's tracker metadata of \a Det, META/TRACKER. template [[nodiscard]] o2::header::DataDescription metaDescription() @@ -209,10 +116,10 @@ class AdaptorTask : public o2::framework::Task { public: using Record = TrackRecord; - using Labels = o2::dataformats::MCTruthContainer; static constexpr bool kITS = std::same_as; - AdaptorTask(bool useMC, std::shared_ptr> conditions) : mUseMC(useMC), mConditions(std::move(conditions)) {} + AdaptorTask(bool useMC, bool staggered, std::shared_ptr> conditions) + : mUseMC(useMC), mStaggered(staggered), mConditions(std::move(conditions)) {} void init(o2::framework::InitContext& ic) final { mConditions->init(ic); } @@ -237,18 +144,13 @@ class AdaptorTask : public o2::framework::Task const auto orbit = pc.services().get().firstTForbit; const auto records = pc.inputs().get>("tracks"); const auto rows = pc.inputs().get>("clusterIndex"); - const auto rofs = pc.inputs().get>("rofs"); - const auto compact = pc.inputs().get>("compclusters"); - const auto patterns = pc.inputs().get>("patterns"); - std::unique_ptr labels; - if (mUseMC) { - labels = pc.inputs().get("labels"); - } + const auto input = readITSMFTInput(pc, mUseMC, mStaggered ? Detector::nLayers : 1); const std::span clusterIndex{rows.data(), rows.size()}; - const std::span frameRecords{rofs.data(), rofs.size()}; - const auto frames = adaptor::clusterFrames(frameRecords, compact.size(), orbit); - const auto pixels = compact.empty() ? std::vector{} : adaptor::pixelCounts({compact.data(), compact.size()}, {patterns.data(), patterns.size()}, *mDictionary); - const auto timing = mConditions->frameTiming().front(); + const auto frames = adaptor::clusterFrames(input, orbit); + const auto pixels = adaptor::pixelCounts(input, *mDictionary); + const auto timings = mConditions->frameTiming(); + const auto clock = mStaggered ? adaptor::clockLayer(timings) : 0U; + const auto frameRecords = input.streams()[clock].rofs; const double bz = nominalBzKiloGauss(); // The fitted tracks, in the tracker's order. @@ -259,17 +161,17 @@ class AdaptorTask : public o2::framework::Task fitted.push_back(track); } } - const auto layerOf = [&](std::uint32_t row) { return mGeometry->surface(compact[clusterIndex[row]].getSensorID()).layer; }; + const auto layerOf = [&](std::uint32_t row) { return mGeometry->surface(input.cluster(clusterIndex[row]).getSensorID()).layer; }; if constexpr (kITS) { - runITS(pc, records, fitted, clusterIndex, frameRecords, frames, pixels, labels.get(), timing, bz, orbit, layerOf); + runITS(pc, records, fitted, clusterIndex, frameRecords, frames, pixels, input, timings, clock, bz, orbit, layerOf); } else { - runMFT(pc, records, fitted, clusterIndex, frameRecords, frames, pixels, labels.get(), bz, orbit, layerOf); + runMFT(pc, records, fitted, clusterIndex, frameRecords, frames, pixels, input, timings, clock, bz, orbit, layerOf); } } private: template - void runITS(o2::framework::ProcessingContext& pc, gsl::span records, std::span fitted, std::span clusterIndex, std::span rofs, const adaptor::ClusterFrames& frames, std::span pixels, const Labels* labels, const FrameTiming& timing, double bz, std::uint32_t orbit, const LayerOf& layerOf) + void runITS(o2::framework::ProcessingContext& pc, gsl::span records, std::span fitted, std::span clusterIndex, std::span rofs, const adaptor::ClusterFrames& frames, std::span pixels, const ITSMFTInput& input, std::span timings, std::size_t clock, double bz, std::uint32_t orbit, const LayerOf& layerOf) { using o2::framework::Output; using AlpideParam = o2::itsmft::DPLAlpideParam; @@ -281,14 +183,14 @@ class AdaptorTask : public o2::framework::Task auto& vertices = pc.outputs().make>(Output{origin, "VERTICES", 0}); auto& vertexROFs = pc.outputs().make>(Output{origin, "VERTICESROF", 0}); - const auto length = static_cast(timing.lengthBC); - const auto delay = static_cast(AlpideParam::Instance().getROFDelayInBC(0)); + const auto length = static_cast(timings[clock].lengthBC); + const auto delay = static_cast(AlpideParam::Instance().getROFDelayInBC(static_cast(clock))); using BC = o2::itsmft::tracking::ROFTimingLayer::BCType; // ITS's frame of a time: its lower edge's, 0 before the first. const o2::itsmft::tracking::ROFTimingLayer rofTiming{ .mROFLength = static_cast(length), .mROFDelay = static_cast(delay), - .mROFBias = static_cast(AlpideParam::Instance().getROFBiasInBC(0)), + .mROFBias = static_cast(AlpideParam::Instance().getROFBiasInBC(static_cast(clock))), }; // A vertex's frame as a track's: the lower edge of its time with the half // width at most half a frame. Its box's own lower edge lies a margin @@ -328,19 +230,19 @@ class AdaptorTask : public o2::framework::Task out.track.getParamOut() = toTrackParCov(record.outwardState, mGeometry->planeOf(record.outwardState), outerAlpha, bz); out.track.setChi2(record.chi2); out.track.setNumberOfClusters(record.count); - if (labels != nullptr) { - out.label = adaptor::trackLabel(record, clusterIndex, *labels, std::span{carried.data(), record.count}); + if (mUseMC) { + out.label = adaptor::trackLabel(record, clusterIndex, input, std::span{carried.data(), record.count}); } std::uint32_t pattern = 0U; for (std::uint8_t i = 0; i < record.count; ++i) { const auto layer = layerOf(record.clusters.index[i]); pattern |= 1U << layer; - pattern |= labels != nullptr && carried[i] ? 1U << (16U + layer) : 0U; + pattern |= mUseMC && carried[i] ? 1U << (16U + layer) : 0U; out.track.setClusterSize(layer, static_cast(pixels[clusterIndex[record.clusters.index[i]]])); } out.track.setPattern(pattern); out.track.setSharedClusters(innermost[clusterIndex[record.clusters.index[0]]] > 1U); - const auto time = adaptor::trackTime(record, clusterIndex, frames, timing); + const auto time = adaptor::trackTime(record, clusterIndex, frames, timings, layerOf); const float half = std::min(static_cast(time.upper - time.lower) / 2.F, static_cast(length) / 2.F); const float centre = static_cast(time.lower + time.upper) / 2.F; out.track.getTimeStamp() = o2::its::TimeStamp(centre, half); @@ -379,16 +281,8 @@ class AdaptorTask : public o2::framework::Task for (const auto& vertex : vertices) { frameCount = std::max(frameCount, vertexFrame(vertex) + 1); } - for (auto* frameRecords : {&trackROFs, &vertexROFs}) { - frameRecords->resize(static_cast(frameCount)); - for (std::int64_t i = 0; i < frameCount; ++i) { - auto& rof = (*frameRecords)[static_cast(i)]; - rof.setBCData(o2::InteractionRecord(static_cast(((i * length) % o2::constants::lhc::LHCMaxBunches) + delay), orbit + static_cast(i * length / o2::constants::lhc::LHCMaxBunches))); - rof.setROFrame(static_cast(i)); - rof.setNEntries(0); - rof.setFirstEntry(-1); - } - } + adaptor::fillClockRecords(trackROFs, static_cast(frameCount), static_cast(length), static_cast(delay), orbit); + vertexROFs = trackROFs; // As ITS: the standard vertexing mode on the input's frames, with the vertexer. if (recoParam().vertexer) { for (std::size_t i = 0; i < std::min(rofs.size(), vertexROFs.size()); ++i) { @@ -402,7 +296,7 @@ class AdaptorTask : public o2::framework::Task const auto& record = records[c.record]; c.track.setFirstClusterEntry(static_cast(clusterIds.size())); for (std::uint8_t i = record.count; i-- > 0;) { // outermost first - clusterIds.push_back(static_cast(clusterIndex[record.clusters.index[i]])); + clusterIds.push_back(input.externalIndex(clusterIndex[record.clusters.index[i]])); } tracks.push_back(c.track); ++entries[static_cast(c.frame)]; @@ -440,7 +334,7 @@ class AdaptorTask : public o2::framework::Task } template - void runMFT(o2::framework::ProcessingContext& pc, gsl::span records, std::span fitted, std::span clusterIndex, std::span rofs, const adaptor::ClusterFrames& frames, std::span pixels, const Labels* labels, double bz, std::uint32_t orbit, const LayerOf& layerOf) + void runMFT(o2::framework::ProcessingContext& pc, gsl::span records, std::span fitted, std::span clusterIndex, std::span rofs, const adaptor::ClusterFrames& frames, std::span pixels, const ITSMFTInput& input, std::span timings, std::size_t clock, double bz, std::uint32_t orbit, const LayerOf& layerOf) { using o2::framework::Output; const auto& origin = detectorOrigin(); @@ -448,6 +342,14 @@ class AdaptorTask : public o2::framework::Task auto& clusterIds = pc.outputs().make>(Output{origin, "TRACKCLSID", 0}); auto& tracks = pc.outputs().make>(Output{origin, "TRACKS", 0}); + using AlpideParam = o2::itsmft::DPLAlpideParam; + using BC = o2::itsmft::tracking::ROFTimingLayer::BCType; + const auto length = timings[clock].lengthBC; + const auto delay = AlpideParam::Instance().getROFDelayInBC(static_cast(clock)); + const o2::itsmft::tracking::ROFTimingLayer rofTiming{ + .mROFLength = static_cast(length), .mROFDelay = static_cast(delay), + .mROFBias = static_cast(AlpideParam::Instance().getROFBiasInBC(static_cast(clock)))}; + struct Converted { o2::mft::TrackMFT track; o2::MCCompLabel label; @@ -467,21 +369,37 @@ class AdaptorTask : public o2::framework::Task out.track.setOutParam(toTrackParCovFwd(record.outwardState, mGeometry->planeOf(record.outwardState), bz, record.outwardChi2)); out.track.setCA(true); out.track.setNumberOfPoints(record.count); - if (labels != nullptr) { - out.label = adaptor::trackLabel(record, clusterIndex, *labels, std::span{carried.data(), record.count}); + if (mUseMC) { + out.label = adaptor::trackLabel(record, clusterIndex, input, std::span{carried.data(), record.count}); } for (std::uint8_t i = 0; i < record.count; ++i) { out.track.setClusterSize(layerOf(record.clusters.index[i]), static_cast(pixels[clusterIndex[record.clusters.index[i]]])); } - // TrackMFT has no time: a track across frames is timed by the earlier. - out.frame = adaptor::earliestRecord(record, clusterIndex, frames); - for (std::uint8_t i = 0; i < record.count; ++i) { - if (frames.record[clusterIndex[record.clusters.index[i]]] != out.frame) { - ++crossing; - break; + // Preserve the monolithic MFT's earliest-record convention. Per-layer + // record numbers have different clocks: use the interval's lower edge + // on the fastest layer's clock instead. + if (mStaggered) { + const auto time = adaptor::trackTime(record, clusterIndex, frames, timings, layerOf); + const float half = std::min(static_cast(time.upper - time.lower) / 2.F, static_cast(length) / 2.F); + const float centre = static_cast(time.lower + time.upper) / 2.F; + out.frame = rofTiming.getROF(o2::its::TimeStamp(centre, half)); + } else { + out.frame = adaptor::earliestRecord(record, clusterIndex, frames); + for (std::uint8_t i = 0; i < record.count; ++i) { + if (frames.record[clusterIndex[record.clusters.index[i]]] != out.frame) { + ++crossing; + break; + } } } } + if (mStaggered) { + std::size_t count = static_cast(o2::constants::lhc::LHCMaxBunches / length) * mConditions->orbitsPerTimeframe(); + for (const auto& track : converted) { + count = std::max(count, static_cast(track.frame) + 1U); + } + adaptor::fillClockRecords(trackROFs, count, length, delay, orbit); + } std::ranges::stable_sort(converted, {}, &Converted::frame); std::vector entries(trackROFs.size() + 1U, 0); @@ -490,7 +408,7 @@ class AdaptorTask : public o2::framework::Task const auto& record = records[c.record]; c.track.setExternalClusterIndexOffset(static_cast(clusterIds.size())); for (std::uint8_t i = 0; i < record.count; ++i) { // innermost first - clusterIds.push_back(static_cast(clusterIndex[record.clusters.index[i]])); + clusterIds.push_back(input.externalIndex(clusterIndex[record.clusters.index[i]])); } tracks.push_back(c.track); ++entries[c.frame]; @@ -507,10 +425,12 @@ class AdaptorTask : public o2::framework::Task trackLabels.push_back(c.label); } } - LOGP(info, "{}: the timeframe at orbit {}: {} tracks ({} failed refits left out, {} across readout frames, timed by the earlier), {} readout frames", Det::name, orbit, tracks.size(), records.size() - fitted.size(), crossing, trackROFs.size()); + LOGP(info, "{}: the timeframe at orbit {}: {} tracks ({} failed refits left out), {} readout frames; {}", Det::name, orbit, tracks.size(), records.size() - fitted.size(), trackROFs.size(), + mStaggered ? std::format("clock layer {}", clock) : std::format("{} across readout frames, timed by the earlier", crossing)); } bool mUseMC = false; + bool mStaggered = false; bool mConfigurationPublished = false; std::shared_ptr> mConditions; ///< the run's geometry, field, frames and orbits per timeframe std::unique_ptr> mGeometry; ///< the planes the tracker's states are bound to @@ -520,12 +440,16 @@ class AdaptorTask : public o2::framework::Task } // namespace adaptor template -o2::framework::DataProcessorSpec getAdaptorSpec(bool useMC) +o2::framework::DataProcessorSpec getAdaptorSpec(bool useMC, bool staggered) { using namespace o2::framework; constexpr bool kITS = std::same_as; const auto origin = detectorOrigin(); - std::vector inputs{InputSpec{"tracks", gDataOrigin, description(kTracks), 0, Lifetime::Timeframe}, InputSpec{"clusterIndex", gDataOrigin, description(kMeasurements), static_cast(MeasurementColumn::clusterIndex), Lifetime::Timeframe}, InputSpec{"rofs", origin, "CLUSTERSROF", 0, Lifetime::Timeframe}, InputSpec{"compclusters", origin, "COMPCLUSTERS", 0, Lifetime::Timeframe}, InputSpec{"patterns", origin, "PATTERNS", 0, Lifetime::Timeframe}, InputSpec{"dictionary", origin, "CLUSDICT", 0, Lifetime::Condition, ccdbParamSpec(origin.template as() + "/Calib/ClusterDictionary")}}; + std::vector inputs{ + InputSpec{"tracks", gDataOrigin, description(kTracks), 0, Lifetime::Timeframe}, + InputSpec{"clusterIndex", gDataOrigin, description(kMeasurements), static_cast(MeasurementColumn::clusterIndex), Lifetime::Timeframe}, + InputSpec{"dictionary", origin, "CLUSDICT", 0, Lifetime::Condition, ccdbParamSpec(origin.template as() + "/Calib/ClusterDictionary")}}; + addITSMFTClusterInputs(inputs, origin, useMC, staggered ? Detector::nLayers : 1); std::vector outputs{{origin, "TRACKS", 0, Lifetime::Timeframe}, {origin, "TRACKCLSID", 0, Lifetime::Timeframe}, {"META", adaptor::metaDescription(), 0, Lifetime::Sporadic}}; if constexpr (kITS) { inputs.emplace_back("measurementGeometry", gDataOrigin, description(kMeasurementGeometry), 0, Lifetime::Timeframe); @@ -539,7 +463,6 @@ o2::framework::DataProcessorSpec getAdaptorSpec(bool useMC) outputs.emplace_back(origin, "MFTTrackROF", 0, Lifetime::Timeframe); } if (useMC) { - inputs.emplace_back("labels", origin, "CLUSTERSMCTR", 0, Lifetime::Timeframe); outputs.emplace_back(origin, "TRACKSMCTR", 0, Lifetime::Timeframe); if constexpr (kITS) { outputs.emplace_back(origin, "VERTICESMCTR", 0, Lifetime::Timeframe); @@ -551,7 +474,7 @@ o2::framework::DataProcessorSpec getAdaptorSpec(bool useMC) .name = deviceName("adaptor"), .inputs = inputs, .outputs = outputs, - .algorithm = AlgorithmSpec{adaptFromTask>(useMC, conditions)}, + .algorithm = AlgorithmSpec{adaptFromTask>(useMC, staggered, conditions)}, }; } diff --git a/Detectors/GenericTracking/workflow/src/AdaptorHelpers.h b/Detectors/GenericTracking/workflow/src/AdaptorHelpers.h new file mode 100644 index 0000000000000..136fb86eb1e81 --- /dev/null +++ b/Detectors/GenericTracking/workflow/src/AdaptorHelpers.h @@ -0,0 +1,199 @@ +// 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 GTRK_WORKFLOW_ADAPTORHELPERS_H_ +#define GTRK_WORKFLOW_ADAPTORHELPERS_H_ + +#include "ITSMFTInput.h" +#include "GenericTracking/Alice/Clusters.h" +#include "GenericTracking/Base/FrameTiming.h" +#include "GenericTracking/Workflow/TrackLabel.h" +#include "DataFormatsITSMFT/TopologyDictionary.h" + +#include +#include +#include +#include +#include +#include +#include + +namespace o2::gtrk::workflow::adaptor +{ + +/// A particle's label, the fake flag unless \a good; O2's noise label for no +/// particle. +[[nodiscard]] inline o2::MCCompLabel labelOf(std::uint64_t particle, bool good) { return particle == kNoTrackParticle ? o2::MCCompLabel(true) : labelOfParticle(particle, !good); } + +/// A timeframe's clusters as the adaptor needs them: each one's readout frame +/// (its record in CLUSTERSROF) and that frame's start in BC from the +/// timeframe's start, without the delay and bias. +struct ClusterFrames { + std::vector record; ///< per cluster + std::vector recordBC; ///< per record +}; + +[[nodiscard]] inline ClusterFrames clusterFrames(std::span rofs, std::size_t clusters, std::uint32_t tfFirstOrbit) +{ + ClusterFrames result; + result.record.assign(clusters, std::numeric_limits::max()); + result.recordBC.reserve(rofs.size()); + std::size_t next = 0U; + for (std::size_t r = 0; r < rofs.size(); ++r) { + if (rofs[r].getFirstEntry() < 0 || static_cast(rofs[r].getFirstEntry()) != next || + rofs[r].getNEntries() < 0 || static_cast(rofs[r].getNEntries()) > clusters - next) { + throw std::runtime_error("cluster ROFs do not cover the input in order"); + } + next += static_cast(rofs[r].getNEntries()); + const auto& ir = rofs[r].getBCData(); + result.recordBC.push_back(((static_cast(ir.orbit) - tfFirstOrbit) * o2::constants::lhc::LHCMaxBunches) + ir.bc); + for (int i = 0; i < rofs[r].getNEntries(); ++i) { + result.record[static_cast(rofs[r].getFirstEntry()) + static_cast(i)] = static_cast(r); + } + } + if (next != clusters) { + throw std::runtime_error("cluster ROFs do not cover the input"); + } + return result; +} + +/// A track's time as ITS takes it (`TrackerTraits`): the intersection of its +/// clusters' frames `[start + delay + bias, + length)` in BC from the +/// timeframe's start; where they do not intersect -- a track across frames, +/// which only a margin lets the tracker make -- the record's box, the frames +/// widened by the margin. [lower, upper) in BC. +struct BCRange { + std::int64_t lower = 0; + std::int64_t upper = 0; +}; + +template +[[nodiscard]] BCRange trackTime(const Record& record, std::span clusterIndex, const ClusterFrames& frames, std::span timings, const LayerOf& layerOf) +{ + BCRange nominal{.lower = std::numeric_limits::min(), .upper = std::numeric_limits::max()}; + for (std::uint8_t i = 0; i < record.count; ++i) { + const auto row = record.clusters.index[i]; + const auto& timing = timings[layerOf(row)]; + const auto start = frames.recordBC[frames.record[clusterIndex[row]]] + timing.firstStartBC; + nominal.lower = std::max(nominal.lower, start); + nominal.upper = std::min(nominal.upper, start + timing.lengthBC); + } + if (nominal.lower < nominal.upper) { + return nominal; + } + return {.lower = bcOfPS(record.timePS), .upper = bcOfPS(static_cast(record.timePS) + record.widthPS)}; +} + +/// The earliest frame (its CLUSTERSROF record) among a track's clusters. +template +[[nodiscard]] std::uint32_t earliestRecord(const Record& record, std::span clusterIndex, const ClusterFrames& frames) +{ + std::uint32_t earliest = std::numeric_limits::max(); + for (std::uint8_t i = 0; i < record.count; ++i) { + earliest = std::min(earliest, frames.record[clusterIndex[record.clusters.index[i]]]); + } + return earliest; +} + +/// Each cluster's pixel count, the converter's decoding (`localClusters`). +[[nodiscard]] inline std::vector pixelCounts(std::span compact, std::span patterns, const o2::itsmft::TopologyDictionary& dictionary) +{ + const auto local = alice::localClusters(compact, patterns, dictionary); + std::vector result(local.size()); + std::ranges::transform(local, result.begin(), [](const auto& cluster) { return cluster.nPixels; }); + return result; +} + +/// The track's particle by the shared rule (TrackLabel.h), its rows' +/// clusters' labels, innermost first; \a carried gets, per row, whether it +/// carries that particle. +template +[[nodiscard]] o2::MCCompLabel trackLabel(const Record& record, std::span clusterIndex, const ITSMFTInput& input, std::span carried) +{ + constexpr std::size_t kRows = std::extent_v>; + std::array rows{}; + std::array track{}; + for (std::uint8_t i = 0; i < record.count; ++i) { + const auto all = input.labels(clusterIndex[record.clusters.index[i]]); + rows[i] = rowParticles(std::span{all.data(), all.size()}); + track[i] = i; + } + const std::span order{track.data(), record.count}; + const auto majority = majorityParticle(rows, order); + const bool good = majority.particle != kNoTrackParticle && majority.carried == record.count; + for (std::uint8_t i = 0; i < record.count; ++i) { + carried[i] = majority.particle != kNoTrackParticle && rows[i].carries(majority.particle); + } + return labelOf(majority.particle, good); +} + +/// Flatten ROF identities in the same layer order as the converter's clusters. +[[nodiscard]] inline ClusterFrames clusterFrames(const ITSMFTInput& input, std::uint32_t orbit) +{ + ClusterFrames result; + result.record.reserve(input.size()); + for (const auto& stream : input.streams()) { + auto frames = clusterFrames(stream.rofs, stream.compact.size(), orbit); + const auto offset = static_cast(result.recordBC.size()); + for (const auto record : frames.record) { + result.record.push_back(offset + record); + } + result.recordBC.insert(result.recordBC.end(), frames.recordBC.begin(), frames.recordBC.end()); + } + return result; +} + +[[nodiscard]] inline std::vector pixelCounts(const ITSMFTInput& input, const o2::itsmft::TopologyDictionary& dictionary) +{ + std::vector result; + result.reserve(input.size()); + for (const auto& stream : input.streams()) { + if (!stream.compact.empty()) { + const auto pixels = pixelCounts(stream.compact, stream.patterns, dictionary); + result.insert(result.end(), pixels.begin(), pixels.end()); + } + } + return result; +} + +/// O2's output clock is the layer with the most frames per timeframe. Equal +/// lengths retain the first layer, independently of the input's occupancy. +[[nodiscard]] inline std::size_t clockLayer(std::span timings) +{ + if (timings.empty() || std::ranges::any_of(timings, [](const auto& timing) { return timing.lengthBC <= 0; })) { + throw std::invalid_argument("output clock needs positive per-layer frame lengths"); + } + return static_cast(std::ranges::min_element(timings, {}, &FrameTiming::lengthBC) - timings.begin()); +} + +/// Frame records carry the delay, whereas frame intervals also carry the bias. +/// Normalise a delayed record across an orbit boundary. +template +inline void fillClockRecords(Records& result, std::size_t count, int lengthBC, int delayBC, std::uint32_t orbit) +{ + if (lengthBC <= 0 || delayBC < 0) { + throw std::invalid_argument("invalid output clock record"); + } + result.resize(count); + for (std::size_t i = 0; i < count; ++i) { + const auto bc = static_cast(i) * lengthBC + delayBC; + auto& rof = result[i]; + rof.setBCData(o2::InteractionRecord(static_cast(bc % o2::constants::lhc::LHCMaxBunches), + orbit + static_cast(bc / o2::constants::lhc::LHCMaxBunches))); + rof.setROFrame(static_cast(i)); + rof.setNEntries(0); + rof.setFirstEntry(-1); + } +} + +} // namespace o2::gtrk::workflow::adaptor + +#endif // GTRK_WORKFLOW_ADAPTORHELPERS_H_ diff --git a/Detectors/GenericTracking/workflow/src/AdaptorITS.cxx b/Detectors/GenericTracking/workflow/src/AdaptorITS.cxx index 4dd2b5194eace..6db81098acd95 100644 --- a/Detectors/GenericTracking/workflow/src/AdaptorITS.cxx +++ b/Detectors/GenericTracking/workflow/src/AdaptorITS.cxx @@ -14,6 +14,6 @@ namespace o2::gtrk { -template o2::framework::DataProcessorSpec workflow::getAdaptorSpec(bool useMC); +template o2::framework::DataProcessorSpec workflow::getAdaptorSpec(bool useMC, bool staggered); } // namespace o2::gtrk diff --git a/Detectors/GenericTracking/workflow/src/AdaptorMFT.cxx b/Detectors/GenericTracking/workflow/src/AdaptorMFT.cxx index 658d1bd9ae11e..29bc019b55c55 100644 --- a/Detectors/GenericTracking/workflow/src/AdaptorMFT.cxx +++ b/Detectors/GenericTracking/workflow/src/AdaptorMFT.cxx @@ -14,6 +14,6 @@ namespace o2::gtrk { -template o2::framework::DataProcessorSpec workflow::getAdaptorSpec(bool useMC); +template o2::framework::DataProcessorSpec workflow::getAdaptorSpec(bool useMC, bool staggered); } // namespace o2::gtrk diff --git a/Detectors/GenericTracking/workflow/src/ClusterConverterITSMFT.h b/Detectors/GenericTracking/workflow/src/ClusterConverterITSMFT.h index a9604105ec6a9..2354bcee31aa0 100644 --- a/Detectors/GenericTracking/workflow/src/ClusterConverterITSMFT.h +++ b/Detectors/GenericTracking/workflow/src/ClusterConverterITSMFT.h @@ -16,6 +16,8 @@ #ifndef GTRK_WORKFLOW_CLUSTERCONVERTERITSMFT_H_ #define GTRK_WORKFLOW_CLUSTERCONVERTERITSMFT_H_ +#include "ITSMFTInput.h" + #include "GenericTracking/Alice/Clusters.h" #include "GenericTracking/Clustering/Cluster.h" #include "GenericTracking/Workflow/Conditions.h" @@ -46,15 +48,18 @@ namespace o2::gtrk::workflow /// O2's compact clusters of the ITS or the MFT, as the O2 cluster reader /// publishes them, turned into GTRK/CLUSTERS for the measurement producer /// and GTRK/CLUSMCTR for the MC truth. The order is O2's, so a -/// measurement's cluster index is O2's cluster index. +/// measurement's cluster index addresses streams concatenated in layer order. template class ClusterConverterTask : public o2::framework::Task { public: - explicit ClusterConverterTask(bool useMC) : mUseMC(useMC) {} + ClusterConverterTask(bool useMC, bool staggered) : mUseMC(useMC), mStreams(staggered ? Detector::nLayers : 1) {} void finaliseCCDB(o2::framework::ConcreteDataMatcher& matcher, void* obj) final { + if (matcher == o2::framework::ConcreteDataMatcher(Conditions::origin(), "ALPIDEPARAM", 0)) { + return; // the CCDB fetch updates the singleton + } if (matcher == o2::framework::ConcreteDataMatcher(Conditions::origin(), "CLUSDICT", 0)) { mDictionary = static_cast(obj); LOGP(info, "{} topology dictionary from the CCDB: {} patterns", Det::name, mDictionary->getSize()); @@ -64,52 +69,61 @@ class ClusterConverterTask : public o2::framework::Task void run(o2::framework::ProcessingContext& pc) final { pc.inputs().get("dictionary"); // delivers a new one through finaliseCCDB - const auto compact = pc.inputs().get>("compclusters"); - const auto patterns = pc.inputs().get>("patterns"); - const auto rofs = pc.inputs().get>("rofs"); - const auto& timing = pc.services().get(); - - std::vector clusters; - if (!rofs.empty()) { - const auto local = alice::localClusters(std::span{compact.data(), compact.size()}, std::span{patterns.data(), patterns.size()}, *mDictionary); - const std::span frames{rofs.data(), rofs.size()}; - clusters = alice::toClusters(local, frames, alice::readoutFrames(frames), timing.firstTForbit); + using AlpideParam = o2::itsmft::DPLAlpideParam::value>; + pc.inputs().get("alpideParam"); + const auto input = readITSMFTInput(pc, mUseMC, mStreams); + const auto& param = AlpideParam::Instance(); + const auto orbit = pc.services().get().firstTForbit; + + std::vector timings; + timings.reserve(Detector::nLayers); + for (int layer = 0; layer < Detector::nLayers; ++layer) { + timings.push_back({.lengthBC = param.getROFLengthInBC(layer), .firstStartBC = param.getROFDelayInBC(layer) + param.getROFBiasInBC(layer)}); + if (mStreams == 1 && timings.back() != timings.front()) { + throw std::runtime_error("staggered ITS/MFT timing needs per-layer cluster input (--enable-its-staggering or --enable-mft-staggering)"); + } } - LOGP(info, "{}: {} clusters in {} readout frames for the timeframe at orbit {}", Det::name, clusters.size(), rofs.size(), timing.firstTForbit); + const auto clusters = convertITSMFTInput(input, *mDictionary, timings, orbit); + std::size_t rofCount = 0U; + for (const auto& stream : input.streams()) { + rofCount += stream.rofs.size(); + } + LOGP(info, "{}: {} clusters in {} readout frames from {} streams for the timeframe at orbit {}", Det::name, clusters.size(), rofCount, mStreams, orbit); pc.outputs().snapshot(o2::framework::Output{gDataOrigin, description(kClusters), 0}, clusters); if (mUseMC) { - const auto labels = pc.inputs().get*>("labels"); - if (labels->getIndexedSize() != clusters.size() && !clusters.empty()) { - LOGP(fatal, "{}: {} labelled clusters for {} clusters", Det::name, labels->getIndexedSize(), clusters.size()); - } auto& flat = pc.outputs().make>(o2::framework::Output{gDataOrigin, description(kClustersMCTruth), 0}); - labels->flatten_to(flat); + concatenateLabels(input).flatten_to(flat); } } private: bool mUseMC; + int mStreams; const o2::itsmft::TopologyDictionary* mDictionary = nullptr; ///< the CCDB cache's, from the run's ITS/MFT Calib/ClusterDictionary }; /// The converter's device for the ITS or the MFT. template -o2::framework::DataProcessorSpec getITSMFTClusterConverterSpec(bool useMC) +o2::framework::DataProcessorSpec getITSMFTClusterConverterSpec(bool useMC, bool staggered) { using namespace o2::framework; const auto origin = Conditions::origin(); - std::vector inputs{InputSpec{"compclusters", origin, "COMPCLUSTERS", 0, Lifetime::Timeframe}, InputSpec{"patterns", origin, "PATTERNS", 0, Lifetime::Timeframe}, InputSpec{"rofs", origin, "CLUSTERSROF", 0, Lifetime::Timeframe}, InputSpec{"dictionary", origin, "CLUSDICT", 0, Lifetime::Condition, ccdbParamSpec(origin.template as() + "/Calib/ClusterDictionary")}}; + const int streams = staggered ? Detector::nLayers : 1; + std::vector inputs; + addITSMFTClusterInputs(inputs, origin, useMC, streams); + const auto detector = origin.template as(); + inputs.emplace_back("dictionary", origin, "CLUSDICT", 0, Lifetime::Condition, ccdbParamSpec(detector + "/Calib/ClusterDictionary")); + inputs.emplace_back("alpideParam", origin, "ALPIDEPARAM", 0, Lifetime::Condition, ccdbParamSpec(detector + "/Config/AlpideParam")); std::vector outputs{OutputSpec{gDataOrigin, description(kClusters), 0, Lifetime::Timeframe}}; if (useMC) { - inputs.emplace_back("labels", origin, "CLUSTERSMCTR", 0, Lifetime::Timeframe); outputs.emplace_back(gDataOrigin, description(kClustersMCTruth), 0, Lifetime::Timeframe); } return DataProcessorSpec{ .name = deviceName("cluster-converter"), .inputs = inputs, .outputs = outputs, - .algorithm = AlgorithmSpec{adaptFromTask>(useMC)}, + .algorithm = AlgorithmSpec{adaptFromTask>(useMC, staggered)}, .options = Options{}, }; } diff --git a/Detectors/GenericTracking/workflow/src/ITS.cxx b/Detectors/GenericTracking/workflow/src/ITS.cxx index 46aad4dc52184..9032eb7887b35 100644 --- a/Detectors/GenericTracking/workflow/src/ITS.cxx +++ b/Detectors/GenericTracking/workflow/src/ITS.cxx @@ -36,9 +36,9 @@ std::vector alice::placeChips() } template <> -o2::framework::DataProcessorSpec workflow::getClusterConverterSpec(bool useMC) +o2::framework::DataProcessorSpec workflow::getClusterConverterSpec(bool useMC, bool staggered) { - return workflow::getITSMFTClusterConverterSpec(useMC); + return workflow::getITSMFTClusterConverterSpec(useMC, staggered); } template o2::framework::DataProcessorSpec workflow::getMeasurementProducerSpec(); diff --git a/Detectors/GenericTracking/workflow/src/ITSMFTInput.h b/Detectors/GenericTracking/workflow/src/ITSMFTInput.h new file mode 100644 index 0000000000000..0e5db921d16b1 --- /dev/null +++ b/Detectors/GenericTracking/workflow/src/ITSMFTInput.h @@ -0,0 +1,181 @@ +// 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 GTRK_WORKFLOW_ITSMFTINPUT_H_ +#define GTRK_WORKFLOW_ITSMFTINPUT_H_ + +#include "GenericTracking/Alice/Clusters.h" +#include "GenericTracking/Base/FrameTiming.h" +#include "DataFormatsITSMFT/ClusterID.h" +#include "DataFormatsITSMFT/CompCluster.h" +#include "DataFormatsITSMFT/ROFRecord.h" +#include "Framework/InputSpec.h" +#include "Framework/InputRecord.h" +#include "Framework/ProcessingContext.h" +#include "SimulationDataFormat/MCCompLabel.h" +#include "SimulationDataFormat/MCTruthContainer.h" + +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace o2::gtrk::workflow +{ + +inline std::string clusterInputBinding(const char* name, int stream, int streams) +{ + return streams == 1 ? name : std::format("{}_{}", name, stream); +} + +/// The same subscription and ordering for the converter and the output adaptor. +inline void addITSMFTClusterInputs(std::vector& inputs, o2::header::DataOrigin origin, bool useMC, int streams) +{ + using o2::framework::Lifetime; + for (int stream = 0; stream < streams; ++stream) { + inputs.emplace_back(clusterInputBinding("compclusters", stream, streams), origin, "COMPCLUSTERS", stream, Lifetime::Timeframe); + inputs.emplace_back(clusterInputBinding("patterns", stream, streams), origin, "PATTERNS", stream, Lifetime::Timeframe); + inputs.emplace_back(clusterInputBinding("rofs", stream, streams), origin, "CLUSTERSROF", stream, Lifetime::Timeframe); + if (useMC) { + inputs.emplace_back(clusterInputBinding("labels", stream, streams), origin, "CLUSTERSMCTR", stream, Lifetime::Timeframe); + } + } +} + +/// GTRK indices address the streams concatenated in layer order. O2 references +/// retain the input stream and its local index, including when a stream is empty. +class ITSMFTInput +{ + public: + using Labels = o2::dataformats::MCTruthContainer; + struct Stream { + std::span compact; + std::span patterns; + std::span rofs; + std::unique_ptr labels; + }; + struct Source { + std::size_t stream; + std::uint32_t index; + }; + + explicit ITSMFTInput(std::vector streams) : mStreams(std::move(streams)) + { + if (mStreams.empty() || mStreams.size() > o2::itsmft::MaxClusLayers) { + throw std::invalid_argument("ITS/MFT cluster input: invalid stream count"); + } + for (const auto& stream : mStreams) { + if (stream.compact.size() > static_cast(o2::itsmft::ClusIndexMask) + 1U || + stream.compact.size() > std::numeric_limits::max() - mOffsets.back()) { + throw std::invalid_argument("ITS/MFT cluster input: too many clusters"); + } + if (stream.labels && stream.labels->getIndexedSize() != stream.compact.size() && !stream.compact.empty()) { + throw std::invalid_argument("ITS/MFT cluster input: labels do not cover the clusters"); + } + mOffsets.push_back(mOffsets.back() + static_cast(stream.compact.size())); + } + } + + [[nodiscard]] std::span streams() const { return mStreams; } + [[nodiscard]] std::uint32_t size() const { return mOffsets.back(); } + [[nodiscard]] std::uint32_t offset(std::size_t stream) const { return mOffsets.at(stream); } + + [[nodiscard]] Source source(std::uint32_t index) const + { + if (index >= size()) { + throw std::out_of_range("ITS/MFT cluster input: cluster index out of range"); + } + const auto stream = static_cast(std::upper_bound(mOffsets.begin(), mOffsets.end(), index) - mOffsets.begin() - 1); + return {stream, index - mOffsets[stream]}; + } + + [[nodiscard]] int externalIndex(std::uint32_t index) const + { + const auto from = source(index); + return o2::itsmft::composeClusID(static_cast(from.stream), static_cast(from.index)); + } + + [[nodiscard]] const o2::itsmft::CompClusterExt& cluster(std::uint32_t index) const + { + const auto from = source(index); + return mStreams[from.stream].compact[from.index]; + } + + [[nodiscard]] auto labels(std::uint32_t index) const + { + const auto from = source(index); + return mStreams[from.stream].labels->getLabels(from.index); + } + + private: + std::vector mStreams; + std::vector mOffsets{0U}; +}; + +/// Decode streams separately, then concatenate clusters and their labels in +/// the same order. The run's frame lengths also cover single-ROF inputs. +inline std::vector convertITSMFTInput(const ITSMFTInput& input, const o2::itsmft::TopologyDictionary& dictionary, + std::span timings, std::uint32_t orbit) +{ + if (timings.size() < input.streams().size()) { + throw std::invalid_argument("ITS/MFT input: missing per-layer timing"); + } + std::vector result; + result.reserve(input.size()); + for (std::size_t stream = 0; stream < input.streams().size(); ++stream) { + const auto& source = input.streams()[stream]; + if (source.rofs.size() > 1U && alice::readoutFrames(source.rofs).lengthBC != timings[stream].lengthBC) { + throw std::runtime_error(std::format("ITS/MFT stream {}: ROF spacing differs from the run's frame length {} BC", stream, timings[stream].lengthBC)); + } + const auto local = alice::localClusters(source.compact, source.patterns, dictionary); + const auto clusters = alice::toClusters(local, source.rofs, {.lengthBC = timings[stream].lengthBC}, orbit); + result.insert(result.end(), clusters.begin(), clusters.end()); + } + return result; +} + +inline ITSMFTInput::Labels concatenateLabels(const ITSMFTInput& input) +{ + ITSMFTInput::Labels labels; + for (const auto& stream : input.streams()) { + if (!stream.compact.empty()) { + labels.mergeAtBack(*stream.labels); + } + } + return labels; +} + +inline ITSMFTInput readITSMFTInput(o2::framework::ProcessingContext& pc, bool useMC, int streams) +{ + std::vector result; + result.reserve(streams); + for (int stream = 0; stream < streams; ++stream) { + const auto compact = pc.inputs().get>(clusterInputBinding("compclusters", stream, streams).c_str()); + const auto patterns = pc.inputs().get>(clusterInputBinding("patterns", stream, streams).c_str()); + const auto rofs = pc.inputs().get>(clusterInputBinding("rofs", stream, streams).c_str()); + result.push_back({.compact = {compact.data(), compact.size()}, .patterns = {patterns.data(), patterns.size()}, .rofs = {rofs.data(), rofs.size()}, + .labels = useMC ? pc.inputs().get(clusterInputBinding("labels", stream, streams).c_str()) : nullptr}); + } + return ITSMFTInput{std::move(result)}; +} + +} // namespace o2::gtrk::workflow + +#endif // GTRK_WORKFLOW_ITSMFTINPUT_H_ diff --git a/Detectors/GenericTracking/workflow/src/MFT.cxx b/Detectors/GenericTracking/workflow/src/MFT.cxx index 70176703632a4..6f191c05cca3d 100644 --- a/Detectors/GenericTracking/workflow/src/MFT.cxx +++ b/Detectors/GenericTracking/workflow/src/MFT.cxx @@ -36,9 +36,9 @@ std::vector alice::placeChips() } template <> -o2::framework::DataProcessorSpec workflow::getClusterConverterSpec(bool useMC) +o2::framework::DataProcessorSpec workflow::getClusterConverterSpec(bool useMC, bool staggered) { - return workflow::getITSMFTClusterConverterSpec(useMC); + return workflow::getITSMFTClusterConverterSpec(useMC, staggered); } template o2::framework::DataProcessorSpec workflow::getMeasurementProducerSpec(); diff --git a/Detectors/GenericTracking/workflow/src/RecoWorkflow.h b/Detectors/GenericTracking/workflow/src/RecoWorkflow.h index 78edc26f25a51..d9387a3638ce2 100644 --- a/Detectors/GenericTracking/workflow/src/RecoWorkflow.h +++ b/Detectors/GenericTracking/workflow/src/RecoWorkflow.h @@ -17,6 +17,7 @@ #define GTRK_WORKFLOW_RECOWORKFLOW_H_ #include "GenericTracking/Workflow/ClusterConverterSpec.h" +#include "GenericTracking/Workflow/ITSMFTOptions.h" #include "GenericTracking/Workflow/MeasurementProducerSpec.h" #include "GenericTracking/Workflow/TrackerSpec.h" @@ -44,6 +45,7 @@ inline void addRecoOptions(std::vector& options, options.push_back(ConfigParamSpec{"clusters-from-upstream", VariantType::Bool, false, {"the clusters from the workflow before this one in the pipe, not O2's readers"}}); options.push_back(ConfigParamSpec{"vertices-from-upstream", VariantType::Bool, false, {"the native vertex seeds (GTRK/VERTICES) from the workflow before this one in the pipe"}}); if (o2Format) { + o2::itsmft::DPLAlpideParamInitializer::addConfigOption(options); options.push_back(ConfigParamSpec{"disable-adaptor", VariantType::Bool, false, {"only the native tracks and vertices (GTRK/...), not O2's formats nor their files"}}); options.push_back(ConfigParamSpec{"disable-root-output", VariantType::Bool, false, {"O2's formats without O2's writers"}}); options.push_back(ConfigParamSpec{"tracking-in-gpu-reco", @@ -79,7 +81,7 @@ template if (!context.options().get("clusters-from-upstream")) { specs.insert(specs.end(), readers.begin(), readers.end()); } - specs.push_back(getClusterConverterSpec(useMC)); + specs.push_back(getClusterConverterSpec(useMC, staggeredInput(context))); specs.push_back(getMeasurementProducerSpec()); const bool inGPUReco = trackingInGPUReco(context); if (inGPUReco && context.options().get("vertices-from-upstream")) { diff --git a/Detectors/GenericTracking/workflow/src/TRK.cxx b/Detectors/GenericTracking/workflow/src/TRK.cxx index d87d91c31d052..fcebdf2f14e03 100644 --- a/Detectors/GenericTracking/workflow/src/TRK.cxx +++ b/Detectors/GenericTracking/workflow/src/TRK.cxx @@ -56,9 +56,11 @@ #include #include #include +#include #include #include #include +#include #include namespace o2::gtrk::alice @@ -246,12 +248,24 @@ class TrkClusterConverterTask : public Task [[nodiscard]] static o2::header::DataOrigin originOf(int layer) { return layer < kBarrelLayers ? o2::header::DataOrigin{"TRK"} : o2::header::DataOrigin{"FT3"}; } [[nodiscard]] static o2::header::DataHeader::SubSpecificationType subspecOf(int layer) { return static_cast(layer < kBarrelLayers ? layer : layer - kBarrelLayers); } - TrkClusterConverterTask(bool useMC, TrkDictionaryConditions dictionaries) : mUseMC(useMC), mDictionaries(dictionaries) {} + TrkClusterConverterTask(bool useMC, TrkDictionaryConditions dictionaries, std::shared_ptr geometryRequest) + : mUseMC(useMC), mDictionaries(dictionaries), mGeometryRequest(std::move(geometryRequest)) {} - void finaliseCCDB(ConcreteDataMatcher& matcher, void* obj) final { mDictionaries.finaliseCCDB(matcher, obj); } + void init(InitContext&) final { o2::base::GRPGeomHelper::instance().setRequest(mGeometryRequest); } + + void finaliseCCDB(ConcreteDataMatcher& matcher, void* obj) final + { + if (!o2::base::GRPGeomHelper::instance().finaliseCCDB(matcher, obj)) { + mDictionaries.finaliseCCDB(matcher, obj); + } + } void run(ProcessingContext& pc) final { + o2::base::GRPGeomHelper::instance().checkUpdates(pc); + // Dictionary coordinates use the VD segmentation's petal count, which + // GeometryTGeo derives from the loaded geometry in this device's process. + o2::trk::GeometryTGeo::Instance(); TrkDictionaryConditions::checkUpdates(pc); const auto& timing = pc.services().get(); std::vector clusters; @@ -299,12 +313,13 @@ class TrkClusterConverterTask : public Task private: bool mUseMC; TrkDictionaryConditions mDictionaries; + std::shared_ptr mGeometryRequest; }; } // namespace template <> -DataProcessorSpec getClusterConverterSpec(bool useMC) +DataProcessorSpec getClusterConverterSpec(bool useMC, bool /*staggered*/) { std::vector inputs; std::vector outputs{OutputSpec{gDataOrigin, description(kClusters), 0, Lifetime::Timeframe}}; @@ -322,11 +337,12 @@ DataProcessorSpec getClusterConverterSpec(bool useMC) outputs.emplace_back(gDataOrigin, description(kClustersMCTruth), 0, Lifetime::Timeframe); } const TrkDictionaryConditions dictionaries(inputs); + auto geometryRequest = std::make_shared(false, false, false, false, false, o2::base::GRPGeomRequest::Aligned, inputs); return DataProcessorSpec{ .name = deviceName("cluster-converter"), .inputs = inputs, .outputs = outputs, - .algorithm = AlgorithmSpec{adaptFromTask(useMC, dictionaries)}, + .algorithm = AlgorithmSpec{adaptFromTask(useMC, dictionaries, geometryRequest)}, .options = Options{}, }; } diff --git a/Detectors/GenericTracking/workflow/src/gtrk-its-reco-workflow.cxx b/Detectors/GenericTracking/workflow/src/gtrk-its-reco-workflow.cxx index 36e47dec18fa9..9e0c9e13786ff 100644 --- a/Detectors/GenericTracking/workflow/src/gtrk-its-reco-workflow.cxx +++ b/Detectors/GenericTracking/workflow/src/gtrk-its-reco-workflow.cxx @@ -65,10 +65,11 @@ WorkflowSpec defineDataProcessing(const ConfigContext& context) namespace gw = o2::gtrk::workflow; using Det = o2::gtrk::detectors::ITS; const bool useMC = !context.options().get("disable-mc"); - auto specs = gw::trackingChain(context, {o2::itsmft::getITSClusterReaderSpec(useMC, false, true, false)}); + const bool staggered = gw::staggeredInput(context); + auto specs = gw::trackingChain(context, {o2::itsmft::getITSClusterReaderSpec(useMC, staggered, true, false)}); // The adaptor reads the tracks: after the GPU reconstruction when the tracking runs inside it. if (!context.options().get("disable-adaptor") && !gw::trackingInGPUReco(context)) { - specs.push_back(gw::getAdaptorSpec(useMC)); + specs.push_back(gw::getAdaptorSpec(useMC, staggered)); if (!context.options().get("disable-root-output")) { specs.push_back(o2::its::getTrackWriterSpec(useMC)); specs.push_back(o2::globaltracking::getIRFrameWriterSpec("irfr:ITS/IRFRAMES/0", "o2_its_irframe.root", "irframe-writer-its")); diff --git a/Detectors/GenericTracking/workflow/src/gtrk-mft-reco-workflow.cxx b/Detectors/GenericTracking/workflow/src/gtrk-mft-reco-workflow.cxx index 105754b5fbc30..9b2bb5e1a02a4 100644 --- a/Detectors/GenericTracking/workflow/src/gtrk-mft-reco-workflow.cxx +++ b/Detectors/GenericTracking/workflow/src/gtrk-mft-reco-workflow.cxx @@ -63,10 +63,11 @@ WorkflowSpec defineDataProcessing(const ConfigContext& context) namespace gw = o2::gtrk::workflow; using Det = o2::gtrk::detectors::MFT; const bool useMC = !context.options().get("disable-mc"); - auto specs = gw::trackingChain(context, {o2::itsmft::getMFTClusterReaderSpec(useMC, false, true, false)}); + const bool staggered = gw::staggeredInput(context); + auto specs = gw::trackingChain(context, {o2::itsmft::getMFTClusterReaderSpec(useMC, staggered, true, false)}); // The adaptor reads the tracks: after the GPU reconstruction when the tracking runs inside it. if (!context.options().get("disable-adaptor") && !gw::trackingInGPUReco(context)) { - specs.push_back(gw::getAdaptorSpec(useMC)); + specs.push_back(gw::getAdaptorSpec(useMC, staggered)); if (!context.options().get("disable-root-output")) { specs.push_back(o2::mft::getTrackWriterSpec(useMC)); } From a7855698fda74fbe7a00466076764fb58048c6ef Mon Sep 17 00:00:00 2001 From: Maximiliano Puccio Date: Tue, 10 Nov 2026 18:13:07 +0100 Subject: [PATCH 05/11] GenericTracking: publish ITS vertex MC truth --- .../tracking/GPU/cuda/DeviceTracker.cuh | 13 ++- .../tracking/GPU/cuda/VertexGlue.cuh | 40 ++++++++- .../Tracking/GPU/DeviceTracker.h | 10 +-- .../test/device/testDeviceTracker.cxx | 17 +++- .../GenericTracking/Workflow/DataSpecs.h | 4 + .../Workflow/GPURecoTracking.h | 2 +- .../GenericTracking/Workflow/TrackerSpec.h | 4 +- .../Workflow/TrackingInterface.h | 9 +- .../GenericTracking/Workflow/VertexLabel.h | 84 +++++++++++++++++++ .../workflow/src/AdaptorDevice.h | 51 +++++++++-- .../workflow/src/GPURecoTracking.cxx | 6 +- .../workflow/src/RecoWorkflow.h | 2 +- .../workflow/src/TrackerDevice.h | 38 +++++---- .../workflow/src/TrackerITS.cxx | 4 +- .../workflow/src/TrackerMFT.cxx | 4 +- .../workflow/src/TrackerTRKFT3.cxx | 2 +- .../src/GPUWorkflowGenericTracking.cxx | 4 +- 17 files changed, 245 insertions(+), 49 deletions(-) create mode 100644 Detectors/GenericTracking/workflow/include/GenericTracking/Workflow/VertexLabel.h diff --git a/Detectors/GenericTracking/tracking/GPU/cuda/DeviceTracker.cuh b/Detectors/GenericTracking/tracking/GPU/cuda/DeviceTracker.cuh index 0af400d69b538..f816f226c8e81 100644 --- a/Detectors/GenericTracking/tracking/GPU/cuda/DeviceTracker.cuh +++ b/Detectors/GenericTracking/tracking/GPU/cuda/DeviceTracker.cuh @@ -641,7 +641,7 @@ class DeviceTrackerImpl final : public DeviceTracker mStream.synchronize(); } - [[nodiscard]] StageStatus findVertices(std::size_t pass, const VertexerSettings& settings, const VertexConstraint& beamSpot, VertexerResult& result) final + [[nodiscard]] StageStatus findVertices(std::size_t pass, const VertexerSettings& settings, const VertexConstraint& beamSpot, VertexerResult& result, std::vector* claims) final { requireNoIteration("findVertices"); if (pass >= mIterationCount) { @@ -649,6 +649,9 @@ class DeviceTrackerImpl final : public DeviceTracker } // Slice by slice, each its own frames' vertices; the lists in time order. result = {}; + if (claims != nullptr) { + claims->clear(); + } const auto& slices = mPlan->slices; StageStatus status = StageStatus::kComplete; for (std::size_t slice = 0; slice < slices.size() && status == StageStatus::kComplete; ++slice) { @@ -658,14 +661,20 @@ class DeviceTrackerImpl final : public DeviceTracker DeviceStageOutput triplets; status = tripletStages(pass, pass, slice, tracklets, permutation, triplets); VertexerResult found; + const auto claimBegin = claims == nullptr ? 0U : claims->size(); if (status == StageStatus::kComplete) { - status = vertex_glue::findVertices(*mRuntime, mArena, mStream, std::as_const(mColumns).view(), static_cast(mRows), constView(tracklets), constView(triplets), triplets.size(), mParameters.data() + pass, settings, beamSpot, ownedFrames(mHostFrame->grid().geometry(), slices[slice].slice.coreBlocks), found); + status = vertex_glue::findVertices(*mRuntime, mArena, mStream, std::as_const(mColumns).view(), static_cast(mRows), constView(tracklets), constView(triplets), triplets.size(), mParameters.data() + pass, settings, beamSpot, ownedFrames(mHostFrame->grid().geometry(), slices[slice].slice.coreBlocks), found, claims); } // A complete vertexer has synchronised after its last readback. if (status != StageStatus::kComplete) { mStream.synchronize(); } closeFrom(Scope::kSlice); + if (claims != nullptr) { + for (std::size_t c = claimBegin; c < claims->size(); ++c) { + (*claims)[c].vertex += static_cast(result.vertices.size()); + } + } appendVertices(result, found); } return status; diff --git a/Detectors/GenericTracking/tracking/GPU/cuda/VertexGlue.cuh b/Detectors/GenericTracking/tracking/GPU/cuda/VertexGlue.cuh index 42482ebf39af3..293ac87ebd241 100644 --- a/Detectors/GenericTracking/tracking/GPU/cuda/VertexGlue.cuh +++ b/Detectors/GenericTracking/tracking/GPU/cuda/VertexGlue.cuh @@ -680,14 +680,29 @@ struct RoundStorage { } // namespace } // namespace detail +/// Download only contributor identities; vertex finding stays on the device. +struct ClaimRows { + const Index* sortedTriplets; + const std::int32_t* claims; + TrackletConstView tracklets; + TripletConstView triplets; + TripletClaim* output; + + GPUdi() void operator()(Index i) const + { + output[i] = {.rows = tripletRows(tracklets, triplets, sortedTriplets[i]), .vertex = claims[i]}; + } +}; + /// The vertices of the pass's \a tripletCount triplets on \a measurements, /// \a rows of them, into \a result (`DeviceTracker::findVertices`); \a pass /// the pass's parameters on the device, for the candidates' errors; the /// vertices of \a owned's frames, a slice's, as `Vertexer::find` finds them. /// Its arrays on the base end, in the caller's scope; kOutOfMemory if the /// arena refuses. +/// With \a keptClaims, append the claimed triplets in the output vertex order. template -[[nodiscard]] StageStatus findVertices(DeviceRuntime& runtime, DeviceArena& arena, DeviceStream& stream, MeasurementView measurements, Index rows, TrackletConstView tracklets, TripletConstView triplets, Index tripletCount, const IterationParameters* pass, const VertexerSettings& settings, const VertexConstraint& beamSpot, const OwnedFrames& owned, VertexerResult& result) +[[nodiscard]] StageStatus findVertices(DeviceRuntime& runtime, DeviceArena& arena, DeviceStream& stream, MeasurementView measurements, Index rows, TrackletConstView tracklets, TripletConstView triplets, Index tripletCount, const IterationParameters* pass, const VertexerSettings& settings, const VertexConstraint& beamSpot, const OwnedFrames& owned, VertexerResult& result, std::vector* keptClaims) { using namespace detail; result = {}; @@ -1059,12 +1074,24 @@ template download(selectedScatterSum, scatterSum.data(), vertexSlots, stream); download(selectedCount, members.data(), vertexSlots, stream); } + std::vector downloadedClaims; + DeviceSpan claimRows; + if (keptClaims != nullptr) { + if (!claimRows.allocate(arena, count)) { + return StageStatus::kOutOfMemory; + } + algorithms::forEach(runtime, stream, count, ClaimRows{std::as_const(sortedTriplets).data(), std::as_const(claims).data(), tracklets, triplets, claimRows.data()}); + downloadedClaims.resize(count); + download(claimRows, downloadedClaims.data(), count, stream); + } stream.synchronize(); + std::vector vertexOfSlot(keptClaims == nullptr ? 0U : vertexSlots, VertexerResult::kNone); std::uint32_t mostInFrame = 0U; for (Index f = 0; f < frames; ++f) { struct Found { VertexSeed seed; Index contributors = 0U; + std::size_t slot = 0U; }; std::vector inFrame; for (std::size_t r = 0; r < found; ++r) { @@ -1080,15 +1107,24 @@ template .earliestPS = earliest[at], .latestPS = latest[at], }; - inFrame.push_back({.seed = vertexSeedOf(sums, beamSpot), .contributors = contributors[at]}); + inFrame.push_back({.seed = vertexSeedOf(sums, beamSpot), .contributors = contributors[at], .slot = at}); } mostInFrame = std::max(mostInFrame, static_cast(inFrame.size())); std::ranges::stable_sort(inFrame, {}, [](const Found& vertex) { return vertex.seed.z; }); for (const auto& vertex : inFrame) { + if (keptClaims != nullptr) { + vertexOfSlot[vertex.slot] = static_cast(result.vertices.size()); + } result.vertices.push_back(vertex.seed); result.contributors.push_back(vertex.contributors); } } + for (auto claim : downloadedClaims) { + if (claim.vertex >= 0) { + claim.vertex = vertexOfSlot[static_cast(claim.vertex)]; + keptClaims->push_back(claim); + } + } result.rounds = mostInFrame + 1U; return StageStatus::kComplete; } diff --git a/Detectors/GenericTracking/tracking/include/GenericTracking/Tracking/GPU/DeviceTracker.h b/Detectors/GenericTracking/tracking/include/GenericTracking/Tracking/GPU/DeviceTracker.h index d2d13a614fa31..f2c29de310dd5 100644 --- a/Detectors/GenericTracking/tracking/include/GenericTracking/Tracking/GPU/DeviceTracker.h +++ b/Detectors/GenericTracking/tracking/include/GenericTracking/Tracking/GPU/DeviceTracker.h @@ -235,10 +235,10 @@ class DeviceTracker /// (`OwnedFrames`), \a beamSpot giving the vertices' transverse spread. /// The vertices, sorted by /// time and each frame's by z, and their contributors into \a result, as - /// `VertexingPass::find` gives them; its counts too, but no per-triplet - /// claims and no beam estimate. Not with an + /// `VertexingPass::find` gives them; its counts too, and optionally their + /// claimed triplets in \a claims; no beam estimate. Not with an /// iteration open; the used mask and the tracks stay as they are. - [[nodiscard]] virtual StageStatus findVertices(std::size_t pass, const VertexerSettings& settings, const VertexConstraint& beamSpot, VertexerResult& result) = 0; + [[nodiscard]] virtual StageStatus findVertices(std::size_t pass, const VertexerSettings& settings, const VertexConstraint& beamSpot, VertexerResult& result, std::vector* claims = nullptr) = 0; /// Replace the open timeframe's plan with \a plan, of the same frame and /// the same run (std::invalid_argument before anything changes @@ -269,7 +269,7 @@ class DeviceTracker /// the pass, iteration \a pass, into \a result (findVertices). The /// timeframe stays open for trackTimeframe; \a frame and \a plan must /// outlive it. - [[nodiscard]] StageStatus findTimeframeVertices(const TimeFrame& frame, const TimeframePlan& plan, std::size_t pass, const VertexerSettings& settings, const VertexConstraint& beamSpot, VertexerResult& result, TimeframeReport& report) + [[nodiscard]] StageStatus findTimeframeVertices(const TimeFrame& frame, const TimeframePlan& plan, std::size_t pass, const VertexerSettings& settings, const VertexConstraint& beamSpot, VertexerResult& result, TimeframeReport& report, std::vector* claims = nullptr) { report = {}; auto start = Clock::now(); @@ -279,7 +279,7 @@ class DeviceTracker return status; } start = Clock::now(); - status = findVertices(pass, settings, beamSpot, result); + status = findVertices(pass, settings, beamSpot, result, claims); report.vertexSeconds = since(start); return status; } diff --git a/Detectors/GenericTracking/tracking/test/device/testDeviceTracker.cxx b/Detectors/GenericTracking/tracking/test/device/testDeviceTracker.cxx index f67c63d8544b1..e65c5913e89c2 100644 --- a/Detectors/GenericTracking/tracking/test/device/testDeviceTracker.cxx +++ b/Detectors/GenericTracking/tracking/test/device/testDeviceTracker.cxx @@ -764,7 +764,7 @@ BOOST_AUTO_TEST_CASE(TheVertexingPassGivesTheHostsListAndTracks) // The host: the pass's list, and the tracking from it. VertexingPass pass{passSettings(), VertexerSettings{}, Fixture::kBeamConstraint, scattering, pool, 2U, 2}; - const auto expected = pass.find(*fixture.frame, 13, 0); + const auto expected = pass.find(*fixture.frame, 13, 0, true); BOOST_REQUIRE(expected.status == StageStatus::kComplete); BOOST_REQUIRE(expected.found.vertices.size() == kCollisionsNS.size()); const auto host = makeHost({.source = VertexSource::kVertexList}); @@ -786,7 +786,8 @@ BOOST_AUTO_TEST_CASE(TheVertexingPassGivesTheHostsListAndTracks) const auto empty = planTimeframe(*fixture.frame, parameters, windows, host->beamLine(), {}, 13, 0); BOOST_REQUIRE((tracker->beginTimeframe(*fixture.frame, empty) == StageStatus::kComplete)); VertexerResult found; - BOOST_REQUIRE((tracker->findVertices(host->iterationCount(), VertexerSettings{}, Fixture::kBeamConstraint, found) == StageStatus::kComplete)); + std::vector claims; + BOOST_REQUIRE((tracker->findVertices(host->iterationCount(), VertexerSettings{}, Fixture::kBeamConstraint, found, &claims) == StageStatus::kComplete)); BOOST_REQUIRE(found.vertices.size() == expected.found.vertices.size()); BOOST_TEST(found.candidates == expected.found.candidates); for (std::size_t v = 0; v < found.vertices.size(); ++v) { @@ -799,6 +800,18 @@ BOOST_AUTO_TEST_CASE(TheVertexingPassGivesTheHostsListAndTracks) } } + // Contributor identities must survive both vertex sorting and slice offsets. + auto expectedClaims = expected.claims; + std::erase_if(expectedClaims, [](const TripletClaim& claim) { return claim.vertex < 0; }); + const auto order = [](const TripletClaim& claim) { return std::pair{claim.vertex, claim.rows}; }; + std::ranges::sort(expectedClaims, {}, order); + std::ranges::sort(claims, {}, order); + BOOST_REQUIRE(claims.size() == expectedClaims.size()); + for (std::size_t c = 0; c < claims.size(); ++c) { + BOOST_TEST(claims[c].vertex == expectedClaims[c].vertex); + BOOST_TEST(claims[c].rows == expectedClaims[c].rows, boost::test_tools::per_element()); + } + // The tracking's plan with the device's list replaces the empty one. const auto trackingPlan = planTimeframe(*fixture.frame, parameters, windows, host->beamLine(), found.vertices, 13, 0); BOOST_REQUIRE((tracker->replacePlan(trackingPlan) == StageStatus::kComplete)); diff --git a/Detectors/GenericTracking/workflow/include/GenericTracking/Workflow/DataSpecs.h b/Detectors/GenericTracking/workflow/include/GenericTracking/Workflow/DataSpecs.h index 271c8b4b5130b..1c10c13d366c4 100644 --- a/Detectors/GenericTracking/workflow/include/GenericTracking/Workflow/DataSpecs.h +++ b/Detectors/GenericTracking/workflow/include/GenericTracking/Workflow/DataSpecs.h @@ -102,6 +102,10 @@ inline constexpr std::string_view kVertexSeeds{"VTX"}; /// Each of those vertices' contributors, a std::uint32_t per vertex. inline constexpr std::string_view kVertexContributors{"VTXCONTR"}; +/// MC contributor identities: TripletClaim records with measurement rows and +/// indices into kVertexSeeds. Empty when MC propagation is disabled. +inline constexpr std::string_view kVertexClaims{"VTXCLAIMS"}; + /// Detector \a Det's name in lower case: "its", "mft", "trk". template [[nodiscard]] std::string detectorName() diff --git a/Detectors/GenericTracking/workflow/include/GenericTracking/Workflow/GPURecoTracking.h b/Detectors/GenericTracking/workflow/include/GenericTracking/Workflow/GPURecoTracking.h index 507e7d3ab1e2a..812b0d3d18ff2 100644 --- a/Detectors/GenericTracking/workflow/include/GenericTracking/Workflow/GPURecoTracking.h +++ b/Detectors/GenericTracking/workflow/include/GenericTracking/Workflow/GPURecoTracking.h @@ -73,7 +73,7 @@ class GPURecoTracking /// The tracking of \a detector, "its" or "mft"; std::invalid_argument for /// another. -[[nodiscard]] std::unique_ptr makeGPURecoTracking(std::string_view detector); +[[nodiscard]] std::unique_ptr makeGPURecoTracking(std::string_view detector, bool useMC = false); } // namespace o2::gtrk::workflow diff --git a/Detectors/GenericTracking/workflow/include/GenericTracking/Workflow/TrackerSpec.h b/Detectors/GenericTracking/workflow/include/GenericTracking/Workflow/TrackerSpec.h index dd6e3840b9b51..3458ab1a4f082 100644 --- a/Detectors/GenericTracking/workflow/include/GenericTracking/Workflow/TrackerSpec.h +++ b/Detectors/GenericTracking/workflow/include/GenericTracking/Workflow/TrackerSpec.h @@ -37,14 +37,14 @@ using TrackRecord = typename TrackRefit::Tracks::Record; /// The device; with \a withVertices the native vertex seeds (kVertices) are /// an input. template -o2::framework::DataProcessorSpec getTrackerSpec(bool withVertices); +o2::framework::DataProcessorSpec getTrackerSpec(bool withVertices, bool useMC = false); class GPURecoTracking; /// \a Det's tracking inside the GPU reconstruction (GPURecoTracking.h), the /// same timeframe as the device's. template -std::unique_ptr makeGPURecoTrackingOf(); +std::unique_ptr makeGPURecoTrackingOf(bool useMC); } // namespace o2::gtrk::workflow diff --git a/Detectors/GenericTracking/workflow/include/GenericTracking/Workflow/TrackingInterface.h b/Detectors/GenericTracking/workflow/include/GenericTracking/Workflow/TrackingInterface.h index fcb9c756acb4b..4bb17002c6815 100644 --- a/Detectors/GenericTracking/workflow/include/GenericTracking/Workflow/TrackingInterface.h +++ b/Detectors/GenericTracking/workflow/include/GenericTracking/Workflow/TrackingInterface.h @@ -289,6 +289,7 @@ class TrackingInterface std::span seeds; const VertexerResult* found = nullptr; ///< the own vertexer's list and contributors, host or device; none without it const typename Vertexing::Result* pass = nullptr; ///< the own vertexer's on the host, with its claims; none without it or on a GPU + std::span claims; ///< the own vertexer's claimed triplets when requested double seconds = 0.; double planSeconds = 0.; ///< the initial device plan, before upload and vertex finding bool dropped = false; ///< out of memory, dropped (dropTFUponFailure): no tracking @@ -308,19 +309,19 @@ class TrackingInterface if (!beginDeviceTimeframe()) { return {.seconds = std::chrono::duration(std::chrono::steady_clock::now() - start).count(), .planSeconds = planSeconds, .dropped = true}; } - const auto status = mDevice->findTimeframeVertices(frame, mDevicePlan, mDeviceIterations.size() - 1U, mVertexing->settings(), beamSpotConstraint(mBeamLine, param), mDeviceFound, mDeviceReport); + const auto status = mDevice->findTimeframeVertices(frame, mDevicePlan, mDeviceIterations.size() - 1U, mVertexing->settings(), beamSpotConstraint(mBeamLine, param), mDeviceFound, mDeviceReport, keepClaims ? &mDeviceClaims : nullptr); const double seconds = std::chrono::duration(std::chrono::steady_clock::now() - start).count(); if (status != StageStatus::kComplete) { return {.seconds = seconds, .planSeconds = planSeconds, .dropped = dropOrFail(status, "the vertexing pass")}; } - return {.seeds = mDeviceFound.vertices, .found = &mDeviceFound, .seconds = seconds, .planSeconds = planSeconds}; + return {.seeds = mDeviceFound.vertices, .found = &mDeviceFound, .claims = keepClaims ? std::span{mDeviceClaims} : std::span{}, .seconds = seconds, .planSeconds = planSeconds}; } mPass.emplace(mVertexing->find(frame, frame.grid().geometry().blocksPerSlice(param.sliceWidthPS), param.sliceHaloPS, keepClaims)); const double seconds = std::chrono::duration(std::chrono::steady_clock::now() - start).count(); if (mPass->status != StageStatus::kComplete) { return {.pass = &*mPass, .seconds = seconds, .dropped = dropOrFail(mPass->status, "the vertexing pass")}; } - return {.seeds = mPass->found.vertices, .found = &mPass->found, .pass = &*mPass, .seconds = seconds}; + return {.seeds = mPass->found.vertices, .found = &mPass->found, .pass = &*mPass, .claims = mPass->claims, .seconds = seconds}; } /// The tracker's run on \a frame from \a seeds. Valid until endTimeframe. @@ -363,6 +364,7 @@ class TrackingInterface mDevicePlan = {}; mDeviceTrackingPlan = {}; mDeviceFound = {}; + mDeviceClaims.clear(); mDeviceTimeframe = false; mVertexed = false; } @@ -553,6 +555,7 @@ class TrackingInterface /// tracks downloaded into the pool, the calls' report. TimeframePlan mDevicePlan; TimeframePlan mDeviceTrackingPlan; + std::vector mDeviceClaims; ///< optional MC contributor identities VertexerResult mDeviceFound; std::shared_ptr mDevicePool; typename TrackRefit::Tracks mDeviceTracks; diff --git a/Detectors/GenericTracking/workflow/include/GenericTracking/Workflow/VertexLabel.h b/Detectors/GenericTracking/workflow/include/GenericTracking/Workflow/VertexLabel.h new file mode 100644 index 0000000000000..6dd3fea663b70 --- /dev/null +++ b/Detectors/GenericTracking/workflow/include/GenericTracking/Workflow/VertexLabel.h @@ -0,0 +1,84 @@ +// 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 GTRK_WORKFLOW_VERTEXLABEL_H_ +#define GTRK_WORKFLOW_VERTEXLABEL_H_ + +#include "GenericTracking/Tracking/VertexStage.h" +#include "GenericTracking/Workflow/TrackLabel.h" + +#include +#include +#include +#include +#include +#include + +namespace o2::gtrk::workflow +{ + +struct VertexTruth { + o2::MCCompLabel label; + float purity = 0.F; +}; + +/// One vote per claimed, cluster-disjoint triplet, from the particle common +/// to its three measurements. Vote for its (source, event), ignoring track ID. +/// Mixed/noise/QED triplets cast no vote but remain in the purity denominator. +/// The largest vote wins (smallest source/event on ties); without a strict +/// majority its label is fake. No valid votes gives an unset label, purity 0. +/// RowLabels maps measurement rows to their original clusters' MC labels. +/// Claims must cover all contributors, before any adapter vertex reordering. +template +[[nodiscard]] std::vector vertexTruth(std::span claims, std::span contributors, const RowLabels& rowLabels) +{ + using Collision = std::pair; // source, event + std::vector> votes(contributors.size()); + std::vector counts(contributors.size(), 0U); + for (const auto& claim : claims) { + if (claim.vertex < 0) { + continue; + } + const auto v = static_cast(claim.vertex); + if (v >= contributors.size()) { + throw std::runtime_error("vertex MC claim refers to a missing vertex"); + } + ++counts[v]; + std::array rows; + for (std::size_t i = 0; i < rows.size(); ++i) { + rows[i] = rowParticles(rowLabels(claim.rows[i])); + } + constexpr std::array order{0U, 1U, 2U}; + if (const auto particle = goodParticle(rows, order); particle != kNoTrackParticle) { + const auto label = labelOfParticle(particle); + ++votes[v][{label.getSourceID(), label.getEventID()}]; + } + } + std::vector result(contributors.size()); + for (std::size_t v = 0; v < result.size(); ++v) { + if (counts[v] != contributors[v]) { + throw std::runtime_error("vertex MC claims do not match its contributor count"); + } + std::uint32_t best = 0U; + for (const auto& [collision, count] : votes[v]) { + if (count > best) { + best = count; + result[v].label = {o2::MCCompLabel::maxTrackID(), collision.second, collision.first, count <= contributors[v] / 2U}; + } + } + result[v].purity = contributors[v] == 0U ? 0.F : static_cast(best) / static_cast(contributors[v]); + } + return result; +} + +} // namespace o2::gtrk::workflow + +#endif // GTRK_WORKFLOW_VERTEXLABEL_H_ diff --git a/Detectors/GenericTracking/workflow/src/AdaptorDevice.h b/Detectors/GenericTracking/workflow/src/AdaptorDevice.h index 832af96fa070c..49766b4c902f0 100644 --- a/Detectors/GenericTracking/workflow/src/AdaptorDevice.h +++ b/Detectors/GenericTracking/workflow/src/AdaptorDevice.h @@ -24,6 +24,7 @@ #define GTRK_WORKFLOW_ADAPTORDEVICE_H_ #include "AdaptorHelpers.h" +#include "GenericTracking/Workflow/VertexLabel.h" #include "GenericTracking/Alice/Clusters.h" #include "GenericTracking/Detectors/ITS.h" @@ -251,19 +252,45 @@ class AdaptorTask : public o2::framework::Task // ITS's order: by the time's lower edge, the tracker's order within. std::ranges::stable_sort(converted, {}, [](const Converted& c) { return c.track.getTimeStamp().getTimeStamp() - c.track.getTimeStamp().getTimeStampError(); }); + // Keep each vertex's MC truth beside it through the output ordering. + struct ConvertedVertex { + o2::its::Vertex vertex; + VertexTruth truth; + }; + std::vector convertedVertices; // The own vertexer's vertices, at the beam line's x and y. const auto shipped = pc.inputs().get>("measurementGeometry"); const auto seeds = pc.inputs().get>("vertices"); const auto contributors = pc.inputs().get>("contributors"); + if (contributors.size() != seeds.size()) { + throw std::runtime_error("vertex seeds and contributor counts differ in size"); + } + std::vector truth; + if (mUseMC) { + const auto claims = pc.inputs().get>("vertexClaims"); + truth = vertexTruth({claims.data(), claims.size()}, {contributors.data(), contributors.size()}, [&](std::uint32_t row) { + if (row >= clusterIndex.size()) { + throw std::runtime_error("vertex MC claim refers to a missing measurement"); + } + const auto labels = input.labels(clusterIndex[row]); + return std::span{labels.data(), labels.size()}; + }); + } + convertedVertices.reserve(seeds.size()); for (std::size_t v = 0; v < seeds.size(); ++v) { const auto& seed = seeds[v]; const std::array cov{seed.varTransverse, 0.F, seed.varTransverse, 0.F, 0.F, seed.varZ}; - auto& vertex = vertices.emplace_back(o2::math_utils::Point3D(static_cast(shipped.beamLine.x), static_cast(shipped.beamLine.y), seed.z), cov, static_cast(std::min(contributors[v], std::numeric_limits::max())), 0.F); + auto& convertedVertex = convertedVertices.emplace_back(); + convertedVertex.truth = mUseMC ? truth[v] : VertexTruth{}; + auto& vertex = convertedVertex.vertex; + vertex = o2::its::Vertex(o2::math_utils::Point3D(static_cast(shipped.beamLine.x), static_cast(shipped.beamLine.y), seed.z), cov, static_cast(std::min(contributors[v], std::numeric_limits::max())), 0.F); const auto lower = bcOfPS(seed.timePS); vertex.setTimeStamp(o2::its::TimeEstBC(static_cast(lower), static_cast(bcOfPS(static_cast(seed.timePS) + seed.widthPS) - lower))); } // By frame, as the frame records count them, then by time and contributors. - std::ranges::stable_sort(vertices, [&](const o2::its::Vertex& a, const o2::its::Vertex& b) { + std::ranges::stable_sort(convertedVertices, [&](const ConvertedVertex& ca, const ConvertedVertex& cb) { + const auto& a = ca.vertex; + const auto& b = cb.vertex; const auto frameA = vertexFrame(a); const auto frameB = vertexFrame(b); if (frameA != frameB) { @@ -272,6 +299,11 @@ class AdaptorTask : public o2::framework::Task return a.getTimeStamp().lower() != b.getTimeStamp().lower() ? a.getTimeStamp().lower() < b.getTimeStamp().lower() : a.getNContributors() > b.getNContributors(); }); + vertices.reserve(convertedVertices.size()); + for (const auto& convertedVertex : convertedVertices) { + vertices.push_back(convertedVertex.vertex); + } + // The clock's frames of the timeframe, as ITS's (`setBCData`), and more // if a track or vertex lies beyond them. std::int64_t frameCount = (o2::constants::lhc::LHCMaxBunches / length) * mConditions->orbitsPerTimeframe(); @@ -326,9 +358,14 @@ class AdaptorTask : public o2::framework::Task for (const auto& c : converted) { trackLabels.push_back(c.label); } - // No vertex truth: the products do not say which clusters contributed. - pc.outputs().make>(Output{origin, "VERTICESMCTR", 0}, vertices.size()); - pc.outputs().make>(Output{origin, "VERTICESMCPUR", 0}, vertices.size(), 0.F); + auto& vertexLabels = pc.outputs().make>(Output{origin, "VERTICESMCTR", 0}); + auto& vertexPurities = pc.outputs().make>(Output{origin, "VERTICESMCPUR", 0}); + vertexLabels.reserve(convertedVertices.size()); + vertexPurities.reserve(convertedVertices.size()); + for (const auto& convertedVertex : convertedVertices) { + vertexLabels.push_back(convertedVertex.truth.label); + vertexPurities.push_back(convertedVertex.truth.purity); + } } LOGP(info, "{}: the timeframe at orbit {}: {} tracks ({} failed refits left out, {} in a frame turned by pi), {} vertices, {} readout frames", Det::name, orbit, tracks.size(), records.size() - fitted.size(), turned, vertices.size(), trackROFs.size()); } @@ -347,8 +384,7 @@ class AdaptorTask : public o2::framework::Task const auto length = timings[clock].lengthBC; const auto delay = AlpideParam::Instance().getROFDelayInBC(static_cast(clock)); const o2::itsmft::tracking::ROFTimingLayer rofTiming{ - .mROFLength = static_cast(length), .mROFDelay = static_cast(delay), - .mROFBias = static_cast(AlpideParam::Instance().getROFBiasInBC(static_cast(clock)))}; + .mROFLength = static_cast(length), .mROFDelay = static_cast(delay), .mROFBias = static_cast(AlpideParam::Instance().getROFBiasInBC(static_cast(clock)))}; struct Converted { o2::mft::TrackMFT track; @@ -465,6 +501,7 @@ o2::framework::DataProcessorSpec getAdaptorSpec(bool useMC, bool staggered) if (useMC) { outputs.emplace_back(origin, "TRACKSMCTR", 0, Lifetime::Timeframe); if constexpr (kITS) { + inputs.emplace_back("vertexClaims", gDataOrigin, description(kVertexClaims), 0, Lifetime::Timeframe); outputs.emplace_back(origin, "VERTICESMCTR", 0, Lifetime::Timeframe); outputs.emplace_back(origin, "VERTICESMCPUR", 0, Lifetime::Timeframe); } diff --git a/Detectors/GenericTracking/workflow/src/GPURecoTracking.cxx b/Detectors/GenericTracking/workflow/src/GPURecoTracking.cxx index d52520edb214e..a3a619f20e273 100644 --- a/Detectors/GenericTracking/workflow/src/GPURecoTracking.cxx +++ b/Detectors/GenericTracking/workflow/src/GPURecoTracking.cxx @@ -23,13 +23,13 @@ namespace o2::gtrk::workflow { -std::unique_ptr makeGPURecoTracking(std::string_view detector) +std::unique_ptr makeGPURecoTracking(std::string_view detector, bool useMC) { if (detector == "its") { - return makeGPURecoTrackingOf(); + return makeGPURecoTrackingOf(useMC); } if (detector == "mft") { - return makeGPURecoTrackingOf(); + return makeGPURecoTrackingOf(useMC); } throw std::invalid_argument("GenericTracking: no tracking inside the GPU reconstruction for " + std::string{detector}); } diff --git a/Detectors/GenericTracking/workflow/src/RecoWorkflow.h b/Detectors/GenericTracking/workflow/src/RecoWorkflow.h index d9387a3638ce2..48ccfbf5bbd17 100644 --- a/Detectors/GenericTracking/workflow/src/RecoWorkflow.h +++ b/Detectors/GenericTracking/workflow/src/RecoWorkflow.h @@ -88,7 +88,7 @@ template throw std::invalid_argument("--tracking-in-gpu-reco takes no native vertex seeds (--vertices-from-upstream)"); } if (!inGPUReco) { - specs.push_back(getTrackerSpec(context.options().get("vertices-from-upstream"))); + specs.push_back(getTrackerSpec(context.options().get("vertices-from-upstream"), useMC)); } return specs; } diff --git a/Detectors/GenericTracking/workflow/src/TrackerDevice.h b/Detectors/GenericTracking/workflow/src/TrackerDevice.h index 2594ee157889a..dab7e47f39231 100644 --- a/Detectors/GenericTracking/workflow/src/TrackerDevice.h +++ b/Detectors/GenericTracking/workflow/src/TrackerDevice.h @@ -60,7 +60,7 @@ namespace o2::gtrk::workflow { /// The tracking device's products made in \a pc's outputs, empty, published on -/// every way out: (tracks, vertices, contributors), DPL's vectors. +/// every way out: (tracks, vertices, contributors, claims), DPL's vectors. template [[nodiscard]] auto makeTrackerProducts(o2::framework::ProcessingContext& pc) { @@ -70,7 +70,9 @@ template auto& tracks = pc.outputs().make>>(Output{gDataOrigin, description(kTracks), 0}); auto& vertices = pc.outputs().make>(Output{gDataOrigin, description(kVertexSeeds), 0}); auto& contributors = pc.outputs().make>(Output{gDataOrigin, description(kVertexContributors), 0}); - return std::tie(tracks, vertices, contributors); + static_assert(o2::framework::is_messageable::value); + auto& claims = pc.outputs().make>(Output{gDataOrigin, description(kVertexClaims), 0}); + return std::tie(tracks, vertices, contributors, claims); } /// The products' output specs. @@ -78,7 +80,7 @@ template [[nodiscard]] std::vector trackerOutputs() { using o2::framework::Lifetime; - return {{gDataOrigin, description(kTracks), 0, Lifetime::Timeframe}, {gDataOrigin, description(kVertexSeeds), 0, Lifetime::Timeframe}, {gDataOrigin, description(kVertexContributors), 0, Lifetime::Timeframe}}; + return {{gDataOrigin, description(kTracks), 0, Lifetime::Timeframe}, {gDataOrigin, description(kVertexSeeds), 0, Lifetime::Timeframe}, {gDataOrigin, description(kVertexContributors), 0, Lifetime::Timeframe}, {gDataOrigin, description(kVertexClaims), 0, Lifetime::Timeframe}}; } /// One timeframe of \a Det's tracking, on \a interface's host or device @@ -86,11 +88,11 @@ template /// tracks, the products. The tracking device's run, and the GPU /// reconstruction's (GPURecoTracking.h). template -void trackTimeframe(o2::framework::ProcessingContext& pc, TrackingInterface& interface, RunConditions& conditions, bool withVertices) +void trackTimeframe(o2::framework::ProcessingContext& pc, TrackingInterface& interface, RunConditions& conditions, bool withVertices, bool useMC) { conditions.checkUpdates(pc); const auto orbit = pc.services().get().firstTForbit; - auto [tracks, vertices, contributors] = makeTrackerProducts(pc); + auto [tracks, vertices, contributors, claims] = makeTrackerProducts(pc); const auto ended = interface.endTimeframeOnExit(); const auto input = readTimeframe(pc, withVertices); @@ -105,7 +107,7 @@ void trackTimeframe(o2::framework::ProcessingContext& pc, TrackingInterface } const TimeFrame frame{measurementView(input.columns), interface.grid(), input.offsets}; - const auto found = interface.vertices(frame, input.seeds, false); + const auto found = interface.vertices(frame, input.seeds, useMC); if (found.dropped) { return; } @@ -117,6 +119,11 @@ void trackTimeframe(o2::framework::ProcessingContext& pc, TrackingInterface if (found.found != nullptr) { vertices.assign(found.found->vertices.begin(), found.found->vertices.end()); contributors.assign(found.found->contributors.begin(), found.found->contributors.end()); + for (const auto& claim : found.claims) { + if (claim.vertex >= 0) { + claims.push_back(claim); + } + } } const auto& result = *tracked.result.tracks; tracks.reserve(result.size()); @@ -153,8 +160,9 @@ class TrackerTask : public o2::framework::Task public: using RunConditions = Conditions; - TrackerTask(bool withVertices, std::shared_ptr conditions) + TrackerTask(bool withVertices, bool useMC, std::shared_ptr conditions) : mWithVertices(withVertices), + mUseMC(useMC), mInterface(withVertices), mConditions(std::move(conditions)) { @@ -176,18 +184,19 @@ class TrackerTask : public o2::framework::Task void finaliseCCDB(o2::framework::ConcreteDataMatcher& matcher, void* obj) final { mConditions->finaliseCCDB(matcher, obj); } - void run(o2::framework::ProcessingContext& pc) final { trackTimeframe(pc, mInterface, *mConditions, mWithVertices); } + void run(o2::framework::ProcessingContext& pc) final { trackTimeframe(pc, mInterface, *mConditions, mWithVertices, mUseMC); } void endOfStream(o2::framework::EndOfStreamContext& /*ec*/) final { reportDropped(mInterface); } private: bool mWithVertices = false; + bool mUseMC = false; TrackingInterface mInterface; std::shared_ptr mConditions; ///< the run's field, geometry, MeanVertex and material LUT (Conditions.h) }; template -o2::framework::DataProcessorSpec getTrackerSpec(bool withVertices) +o2::framework::DataProcessorSpec getTrackerSpec(bool withVertices, bool useMC) { using namespace o2::framework; auto inputs = measurementInputs(withVertices); @@ -196,7 +205,7 @@ o2::framework::DataProcessorSpec getTrackerSpec(bool withVertices) .name = deviceName("tracker"), .inputs = inputs, .outputs = trackerOutputs(), - .algorithm = AlgorithmSpec{adaptFromTask>(withVertices, conditions)}, + .algorithm = AlgorithmSpec{adaptFromTask>(withVertices, useMC, conditions)}, .options = Options{{"gtrk-gpu-backend", VariantType::String, "", {"the GPU backend to track on, cuda or hip; the host if empty"}}, {"gtrk-gpu-device", VariantType::Int, 0, {"the GPU backend's device index"}}, {"gtrk-gpu-arena-mb", VariantType::Int, 0, {"the device tracker's arena in O2's GPU pool [MB]; required with a GPU backend"}}}, }; } @@ -212,7 +221,7 @@ class GPURecoTrackingOf final : public GPURecoTracking /// No native seeds: the vertex list, where an iteration needs one, is the /// own vertexer's. - GPURecoTrackingOf() : mInterface(false) {} + explicit GPURecoTrackingOf(bool useMC) : mUseMC(useMC), mInterface(false) {} void appendInputs(std::vector& inputs) final { @@ -245,7 +254,7 @@ class GPURecoTrackingOf final : public GPURecoTracking // Given back after the timeframe's end, which detaches the tracker. const auto borrowed = mInterface.context().borrow(tag()); mInterface.useBlock(borrowed->block()); - trackTimeframe(pc, mInterface, *mConditions, false); + trackTimeframe(pc, mInterface, *mConditions, false, mUseMC); LOGP(info, "{}: the GPU reconstruction's pool lent {:.1f} MB, the arena's peak {:.1f} MB", Det::name, static_cast(borrowed->block().bytes) / 1.e6, static_cast(mInterface.lastArenaPeakBytes()) / 1.e6); } @@ -266,14 +275,15 @@ class GPURecoTrackingOf final : public GPURecoTracking return packed; } + bool mUseMC = false; TrackingInterface mInterface; std::unique_ptr mConditions; ///< the detector's own conditions; the run-wide ones are the GPU reconstruction's }; template -std::unique_ptr makeGPURecoTrackingOf() +std::unique_ptr makeGPURecoTrackingOf(bool useMC) { - return std::make_unique>(); + return std::make_unique>(useMC); } } // namespace o2::gtrk::workflow diff --git a/Detectors/GenericTracking/workflow/src/TrackerITS.cxx b/Detectors/GenericTracking/workflow/src/TrackerITS.cxx index ac0802b711b4d..033287c486157 100644 --- a/Detectors/GenericTracking/workflow/src/TrackerITS.cxx +++ b/Detectors/GenericTracking/workflow/src/TrackerITS.cxx @@ -21,7 +21,7 @@ namespace o2::gtrk { -template o2::framework::DataProcessorSpec workflow::getTrackerSpec(bool withVertices); -template std::unique_ptr workflow::makeGPURecoTrackingOf(); +template o2::framework::DataProcessorSpec workflow::getTrackerSpec(bool withVertices, bool useMC); +template std::unique_ptr workflow::makeGPURecoTrackingOf(bool useMC); } // namespace o2::gtrk diff --git a/Detectors/GenericTracking/workflow/src/TrackerMFT.cxx b/Detectors/GenericTracking/workflow/src/TrackerMFT.cxx index 6167ddfda6944..0a9df8abdb37e 100644 --- a/Detectors/GenericTracking/workflow/src/TrackerMFT.cxx +++ b/Detectors/GenericTracking/workflow/src/TrackerMFT.cxx @@ -21,7 +21,7 @@ namespace o2::gtrk { -template o2::framework::DataProcessorSpec workflow::getTrackerSpec(bool withVertices); -template std::unique_ptr workflow::makeGPURecoTrackingOf(); +template o2::framework::DataProcessorSpec workflow::getTrackerSpec(bool withVertices, bool useMC); +template std::unique_ptr workflow::makeGPURecoTrackingOf(bool useMC); } // namespace o2::gtrk diff --git a/Detectors/GenericTracking/workflow/src/TrackerTRKFT3.cxx b/Detectors/GenericTracking/workflow/src/TrackerTRKFT3.cxx index f91794c10d631..0d6946d216083 100644 --- a/Detectors/GenericTracking/workflow/src/TrackerTRKFT3.cxx +++ b/Detectors/GenericTracking/workflow/src/TrackerTRKFT3.cxx @@ -20,6 +20,6 @@ namespace o2::gtrk { -template o2::framework::DataProcessorSpec workflow::getTrackerSpec(bool withVertices); +template o2::framework::DataProcessorSpec workflow::getTrackerSpec(bool withVertices, bool useMC); } // namespace o2::gtrk diff --git a/GPU/Workflow/src/GPUWorkflowGenericTracking.cxx b/GPU/Workflow/src/GPUWorkflowGenericTracking.cxx index ec91e0ee2ea50..32ada64ae827c 100644 --- a/GPU/Workflow/src/GPUWorkflowGenericTracking.cxx +++ b/GPU/Workflow/src/GPUWorkflowGenericTracking.cxx @@ -22,10 +22,10 @@ namespace o2::gpu void GPURecoWorkflowSpec::createGenericTracking() { if (mSpecConfig.runGenericTrackingITS) { - mGenericTracking.push_back(o2::gtrk::workflow::makeGPURecoTracking("its")); + mGenericTracking.push_back(o2::gtrk::workflow::makeGPURecoTracking("its", mSpecConfig.processMC)); } if (mSpecConfig.runGenericTrackingMFT) { - mGenericTracking.push_back(o2::gtrk::workflow::makeGPURecoTracking("mft")); + mGenericTracking.push_back(o2::gtrk::workflow::makeGPURecoTracking("mft", mSpecConfig.processMC)); } } From ca5e27ae2a9cb75d02f8f8209d188f73a11d26bd Mon Sep 17 00:00:00 2001 From: Maximiliano Puccio Date: Tue, 10 Nov 2026 20:04:47 +0100 Subject: [PATCH 06/11] ITS: apply vertex multiplicity cuts to ROF masks --- .../ITS/tracking/include/ITStracking/FastMultEst.h | 13 +++++++++---- 1 file changed, 9 insertions(+), 4 deletions(-) diff --git a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/FastMultEst.h b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/FastMultEst.h index d283d6eed2d45..d79d86e98975c 100644 --- a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/FastMultEst.h +++ b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/FastMultEst.h @@ -24,6 +24,7 @@ #include "ITStracking/FastMultEstConfig.h" #include "ITSMFTTracking/ROFLookupTables.h" #include +#include #include namespace o2::its @@ -62,10 +63,14 @@ struct FastMultEst { if (!multEstConf.isPassingVtxMultCut(vertex.getNContributors())) { const auto& timestamp{vertex.getTimeStamp()}; for (int layer = 0; layer < NLayers; ++layer) { - uint32_t startROF = sel.getLayer(layer).getROF(timestamp.lower()); - uint32_t endROF = sel.getLayer(layer).getROF(timestamp.upper()); - for (uint32_t rof = startROF; rof <= endROF; ++rof) { - sel.setROFsEnabled(layer, rof, 0); + const auto& timing = sel.getLayer(layer); + if (timing.mNROFsTF == 0) { + continue; + } + const auto startROF = timing.getROF(timestamp.lower()); + const auto endROF = std::min(timing.getROF(timestamp.upper()), timing.mNROFsTF - 1U); + if (startROF <= endROF && startROF < timing.mNROFsTF) { + sel.setROFsEnabled(layer, startROF, endROF - startROF + 1U, 0); } } } From eec3fca8a42960babc44c14150625a2fc93b4195 Mon Sep 17 00:00:00 2001 From: Maximiliano Puccio Date: Tue, 10 Nov 2026 20:04:47 +0100 Subject: [PATCH 07/11] GenericTracking: apply ITS ROF and multiplicity selection --- .../GenericTracking/workflow/CMakeLists.txt | 4 +- .../GenericTracking/Workflow/DataSpecs.h | 4 + .../Workflow/ITSROFSelection.h | 68 +++++ .../GenericTracking/Workflow/RecoParam.h | 4 + .../Workflow/TrackingInterface.h | 4 +- .../workflow/src/AdaptorDevice.h | 12 +- .../workflow/src/ITSROFSelection.cxx | 247 ++++++++++++++++++ .../workflow/src/TrackerDevice.h | 95 +++++-- .../workflow/src/gtrk-its-reco-workflow.cxx | 18 +- 9 files changed, 431 insertions(+), 25 deletions(-) create mode 100644 Detectors/GenericTracking/workflow/include/GenericTracking/Workflow/ITSROFSelection.h create mode 100644 Detectors/GenericTracking/workflow/src/ITSROFSelection.cxx diff --git a/Detectors/GenericTracking/workflow/CMakeLists.txt b/Detectors/GenericTracking/workflow/CMakeLists.txt index 67af4f0273036..fcb9eb0b4d33c 100644 --- a/Detectors/GenericTracking/workflow/CMakeLists.txt +++ b/Detectors/GenericTracking/workflow/CMakeLists.txt @@ -15,11 +15,11 @@ # unit, ALICE 3's (TRK.cxx) with the upgrades, and the tracking device # (TrackerDevice.h), one unit per detector, each instantiating the tracker. set(sources src/Clusters.cxx src/Placement.cxx src/RecoParam.cxx src/ITS.cxx src/MFT.cxx src/TrackerITS.cxx src/TrackerMFT.cxx src/AdaptorITS.cxx - src/AdaptorMFT.cxx src/GPURecoTracking.cxx) + src/AdaptorMFT.cxx src/GPURecoTracking.cxx src/ITSROFSelection.cxx) set(trackerUnits src/TrackerITS.cxx src/TrackerMFT.cxx) set(publicLibraries O2::GTrkTracking O2::Framework O2::CommonUtils O2::DataFormatsCalibration O2::DetectorsBase O2::DataFormatsITSMFT O2::ITSBase O2::MFTBase O2::DetectorsRaw ROOT::Core) -set(privateLibraries O2::ITSMFTBase O2::ITSMFTTracking O2::DataFormatsITS O2::DataFormatsMFT O2::SimulationDataFormat O2::CommonConstants +set(privateLibraries O2::ITStracking O2::DataFormatsTRD O2::ITSMFTBase O2::ITSMFTTracking O2::DataFormatsITS O2::DataFormatsMFT O2::SimulationDataFormat O2::CommonConstants O2::MathUtils O2::ReconstructionDataFormats ROOT::Geom ROOT::RIO ROOT::Tree) if(ENABLE_UPGRADES) list(APPEND sources src/TRK.cxx src/TrackerTRKFT3.cxx) diff --git a/Detectors/GenericTracking/workflow/include/GenericTracking/Workflow/DataSpecs.h b/Detectors/GenericTracking/workflow/include/GenericTracking/Workflow/DataSpecs.h index 1c10c13d366c4..bf8c736a25fa3 100644 --- a/Detectors/GenericTracking/workflow/include/GenericTracking/Workflow/DataSpecs.h +++ b/Detectors/GenericTracking/workflow/include/GenericTracking/Workflow/DataSpecs.h @@ -106,6 +106,10 @@ inline constexpr std::string_view kVertexContributors{"VTXCONTR"}; /// indices into kVertexSeeds. Empty when MC propagation is disabled. inline constexpr std::string_view kVertexClaims{"VTXCLAIMS"}; +/// ITS clock-ROF acceptance after cluster/random and vertex multiplicity cuts. +/// One uint8_t per clock ROF; empty for a dropped/empty timeframe. +inline constexpr std::string_view kSelectedROFs{"ROFSEL"}; + /// Detector \a Det's name in lower case: "its", "mft", "trk". template [[nodiscard]] std::string detectorName() diff --git a/Detectors/GenericTracking/workflow/include/GenericTracking/Workflow/ITSROFSelection.h b/Detectors/GenericTracking/workflow/include/GenericTracking/Workflow/ITSROFSelection.h new file mode 100644 index 0000000000000..f706a9ca7f799 --- /dev/null +++ b/Detectors/GenericTracking/workflow/include/GenericTracking/Workflow/ITSROFSelection.h @@ -0,0 +1,68 @@ +// 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 GTRK_WORKFLOW_ITSROFSELECTION_H_ +#define GTRK_WORKFLOW_ITSROFSELECTION_H_ + +#include "GenericTracking/Base/FrameTiming.h" +#include "GenericTracking/Tracking/Measurement.h" +#include "GenericTracking/Tracking/PairCut.h" +#include "Framework/InputSpec.h" + +#include +#include +#include +#include + +namespace o2::framework +{ +class ProcessingContext; +} + +namespace o2::gtrk::workflow +{ + +/// ITS's native fastMultConfig selection: cluster/random selection before +/// vertexing, vertex-contributor selection before tracking. It consumes raw +/// cluster ROFs so that expanded measurement time margins cannot change a +/// cluster's ROF identity. Both host and device receive the selected rows. +class ITSROFSelection +{ + public: + static void addInputs(std::vector& inputs); + ITSROFSelection(o2::framework::ProcessingContext& pc, std::span timing, int orbitsPerTF); + ~ITSROFSelection(); + ITSROFSelection(const ITSROFSelection&) = delete; + ITSROFSelection& operator=(const ITSROFSelection&) = delete; + + /// Vertices remain in the output, as in native ITS; the cut disables their + /// overlapping ROFs for tracking. Returns whether a vertex cut is active. + bool selectWithVertices(std::span seeds, std::span contributors); + [[nodiscard]] std::vector clockMask() const; + + struct Measurements { + MeasurementStore columns; + std::vector offsets; + std::vector originalRows; ///< selected row -> original input measurement row + }; + + /// Preserve bucket ordering and original measurement identities. No copy + /// when every row passes. originalRows composes mappings after vertex cuts. + [[nodiscard]] std::optional filter(MeasurementView source, std::span offsets, std::span originalRows = {}) const; + + private: + struct State; + std::unique_ptr mState; +}; + +} // namespace o2::gtrk::workflow + +#endif // GTRK_WORKFLOW_ITSROFSELECTION_H_ diff --git a/Detectors/GenericTracking/workflow/include/GenericTracking/Workflow/RecoParam.h b/Detectors/GenericTracking/workflow/include/GenericTracking/Workflow/RecoParam.h index 560bc3cd8da32..497b523d82e44 100644 --- a/Detectors/GenericTracking/workflow/include/GenericTracking/Workflow/RecoParam.h +++ b/Detectors/GenericTracking/workflow/include/GenericTracking/Workflow/RecoParam.h @@ -269,6 +269,10 @@ struct RecoSettings { /// empty products and is counted, instead of failing the run. Off, the /// default, as ITS's DropTFUponFailure: the run fails. bool dropTFUponFailure = false; + + /// ITS ROF selection uses native fastMultConfig. Trigger input: 0 none, + /// 1 ITS physics triggers, 2 TRD tracklet triggers. Other detectors ignore it. + int selectWithTriggers = 0; }; /// Detector \a N's settings (`o2::detectors::DetID`): `GTRKReco`. diff --git a/Detectors/GenericTracking/workflow/include/GenericTracking/Workflow/TrackingInterface.h b/Detectors/GenericTracking/workflow/include/GenericTracking/Workflow/TrackingInterface.h index 4bb17002c6815..509e4d616e9f3 100644 --- a/Detectors/GenericTracking/workflow/include/GenericTracking/Workflow/TrackingInterface.h +++ b/Detectors/GenericTracking/workflow/include/GenericTracking/Workflow/TrackingInterface.h @@ -497,9 +497,9 @@ class TrackingInterface [[nodiscard]] Tracked trackOnDevice(const Frame& frame, std::span seeds, std::int32_t blocksPerSlice, std::chrono::steady_clock::time_point start) { const auto& param = recoParam(); - const bool replaced = mVertexing != nullptr; + const bool replaced = mVertexed; const std::size_t iterations = mTracker->iterationCount(); - mDeviceTrackingPlan = replaced ? planTimeframe(frame, mDeviceIterations, mDeviceWindows, mTracker->beamLine(), seeds, blocksPerSlice, param.sliceHaloPS) : mTracker->planTimeframe(frame, seeds, blocksPerSlice, param.sliceHaloPS); + mDeviceTrackingPlan = mVertexing != nullptr ? planTimeframe(frame, mDeviceIterations, mDeviceWindows, mTracker->beamLine(), seeds, blocksPerSlice, param.sliceHaloPS) : mTracker->planTimeframe(frame, seeds, blocksPerSlice, param.sliceHaloPS); Tracked tracked; tracked.planSeconds = std::chrono::duration(std::chrono::steady_clock::now() - start).count(); if (!beginDeviceTimeframe()) { diff --git a/Detectors/GenericTracking/workflow/src/AdaptorDevice.h b/Detectors/GenericTracking/workflow/src/AdaptorDevice.h index 49766b4c902f0..12ff2ce1d7f56 100644 --- a/Detectors/GenericTracking/workflow/src/AdaptorDevice.h +++ b/Detectors/GenericTracking/workflow/src/AdaptorDevice.h @@ -108,6 +108,9 @@ void publishConfiguration(o2::framework::ProcessingContext& pc, bool& published) TMap md; md.SetOwnerKeyValue(); md.Add(new TObjString(name.c_str()), new TObjString(o2::conf::ConfigurableParam::asJSON(name).c_str())); + if constexpr (std::same_as) { + md.Add(new TObjString("fastMultConfig"), new TObjString(o2::conf::ConfigurableParam::asJSON("fastMultConfig").c_str())); + } pc.outputs().snapshot(o2::framework::Output{"META", metaDescription(), 0}, md); } @@ -177,6 +180,7 @@ class AdaptorTask : public o2::framework::Task using o2::framework::Output; using AlpideParam = o2::itsmft::DPLAlpideParam; const auto& origin = detectorOrigin(); + const auto selectedROFs = pc.inputs().get>("selectedROFs"); auto& irFrames = pc.outputs().make>(Output{origin, "IRFRAMES", 0}); auto& clusterIds = pc.outputs().make>(Output{origin, "TRACKCLSID", 0}); auto& tracks = pc.outputs().make>(Output{origin, "TRACKS", 0}); @@ -338,9 +342,10 @@ class AdaptorTask : public o2::framework::Task trackROFs[i].setFirstEntry(entries[i]); trackROFs[i].setNEntries(entries[i + 1] - entries[i]); trackROFs[i].setFlags(vertexROFs[i].getFlags()); - // Every frame is processed: no multiplicity selection here. - auto& irFrame = irFrames.emplace_back(trackROFs[i].getBCData(), trackROFs[i].getBCData() + (length - 1)); - irFrame.info = static_cast(trackROFs[i].getNEntries()); + if (i < selectedROFs.size() && selectedROFs[i] != 0U) { + auto& irFrame = irFrames.emplace_back(trackROFs[i].getBCData(), trackROFs[i].getBCData() + (length - 1)); + irFrame.info = static_cast(trackROFs[i].getNEntries()); + } } std::ranges::fill(entries, 0); for (const auto& vertex : vertices) { @@ -490,6 +495,7 @@ o2::framework::DataProcessorSpec getAdaptorSpec(bool useMC, bool staggered) if constexpr (kITS) { inputs.emplace_back("measurementGeometry", gDataOrigin, description(kMeasurementGeometry), 0, Lifetime::Timeframe); inputs.emplace_back("vertices", gDataOrigin, description(kVertexSeeds), 0, Lifetime::Timeframe); + inputs.emplace_back("selectedROFs", gDataOrigin, description(kSelectedROFs), 0, Lifetime::Timeframe); inputs.emplace_back("contributors", gDataOrigin, description(kVertexContributors), 0, Lifetime::Timeframe); outputs.emplace_back(origin, "ITSTrackROF", 0, Lifetime::Timeframe); outputs.emplace_back(origin, "VERTICES", 0, Lifetime::Timeframe); diff --git a/Detectors/GenericTracking/workflow/src/ITSROFSelection.cxx b/Detectors/GenericTracking/workflow/src/ITSROFSelection.cxx new file mode 100644 index 0000000000000..6db18eb4f2f21 --- /dev/null +++ b/Detectors/GenericTracking/workflow/src/ITSROFSelection.cxx @@ -0,0 +1,247 @@ +// 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 "GenericTracking/Workflow/ITSROFSelection.h" +#include "GenericTracking/Workflow/RecoParam.h" +#include "DataFormatsITSMFT/DPLAlpideParam.h" +#include "DataFormatsITSMFT/PhysTrigger.h" +#include "DataFormatsTRD/TriggerRecord.h" +#include "Framework/DataRefUtils.h" +#include "Framework/InputRecord.h" +#include "Framework/InputRecordWalker.h" +#include "Framework/ProcessingContext.h" +#include "Framework/TimingInfo.h" +#include "ITStracking/FastMultEst.h" + +#include +#include +#include +#include + +namespace o2::gtrk::workflow +{ + +namespace +{ +using Estimator = o2::its::FastMultEst; +using Overlap = Estimator::ROFOverlapTableN; +using Mask = Estimator::ROFMaskTableN; +using Timing = o2::itsmft::tracking::ROFTimingLayer; +using Alpide = o2::itsmft::DPLAlpideParam; +constexpr int kLayers = Estimator::NLayers; + +int triggerType() +{ + const auto type = recoParam().selectWithTriggers; + if (type < 0 || type > 2) { + throw std::invalid_argument("GTRKRecoITS.selectWithTriggers must be 0 (none), 1 (ITS physics) or 2 (TRD)"); + } + return type; +} +} // namespace + +struct ITSROFSelection::State { + struct ClusterROF { + int layer; + int rof; + }; + Overlap overlap; + Mask mask; + std::vector clusterROFs; +}; + +void ITSROFSelection::addInputs(std::vector& inputs) +{ + using namespace o2::framework; + inputs.emplace_back("its_selection_clusters", ConcreteDataTypeMatcher{"ITS", "COMPCLUSTERS"}, Lifetime::Timeframe); + inputs.emplace_back("its_selection_rofs", ConcreteDataTypeMatcher{"ITS", "CLUSTERSROF"}, Lifetime::Timeframe); + if (triggerType() == 1) { + inputs.emplace_back("its_selection_triggers", "ITS", "PHYSTRIG", 0, Lifetime::Timeframe); + } else if (triggerType() == 2) { + inputs.emplace_back("its_selection_triggers", "TRD", "TRKTRGRD", 0, Lifetime::Timeframe); + } + if (o2::its::FastMultEstConfig::Instance().isVtxMultCutRequested() && !recoParam().vertexer) { + throw std::invalid_argument("fastMultConfig vertex cuts require GTRKRecoITS.vertexer=true (upstream VertexSeed inputs have no contributor counts)"); + } +} + +ITSROFSelection::ITSROFSelection(o2::framework::ProcessingContext& pc, std::span timing, int orbitsPerTF) + : mState(std::make_unique()) +{ + using namespace o2::framework; + std::array, kLayers> clusters{}; + std::array, kLayers> rofs{}; + std::array seen{}; + for (const auto& ref : InputRecordWalker{pc.inputs()}) { + const auto* header = DataRefUtils::getHeader(ref); + if (header->dataOrigin != o2::header::DataOrigin{"ITS"}) { + continue; + } + const auto layer = header->subSpecification; + if (header->dataDescription == o2::header::DataDescription{"COMPCLUSTERS"}) { + if (layer >= kLayers) { + throw std::runtime_error("ITS selection: invalid cluster stream"); + } + clusters[layer] = pc.inputs().get>(ref); + } else if (header->dataDescription == o2::header::DataDescription{"CLUSTERSROF"}) { + if (layer >= kLayers) { + throw std::runtime_error("ITS selection: invalid ROF stream"); + } + rofs[layer] = pc.inputs().get>(ref); + seen[layer] = true; + } + } + const auto streams = std::ranges::count(seen, true); + if (streams != 1 && streams != kLayers) { + throw std::runtime_error("ITS selection requires one monolithic stream or seven layer streams"); + } + const bool staggered = streams == kLayers; + if (timing.size() != kLayers || !seen[0]) { + throw std::runtime_error("ITS selection: missing layer timing or stream zero"); + } + auto& state = *mState; + const auto& alpide = Alpide::Instance(); + const auto& config = recoParam(); + for (int layer = 0; layer < kLayers; ++layer) { + const int nROFs = std::max(orbitsPerTF * (o2::constants::lhc::LHCMaxBunches / timing[layer].lengthBC), static_cast(rofs[staggered ? layer : 0].size())); + state.overlap.defineLayer(layer, {static_cast(nROFs), static_cast(timing[layer].lengthBC), static_cast(alpide.getROFDelayInBC(layer)), static_cast(alpide.getROFBiasInBC(layer)), static_cast(config.frameMarginBC.empty() ? 0 : config.frameMarginBC[layer])}); + } + state.overlap.init(); + state.mask = Mask{state.overlap}; + static const o2::itsmft::ChipMappingITS mapping; + for (int stream = 0; stream < streams; ++stream) { + const auto offset = state.clusterROFs.size(); + state.clusterROFs.resize(offset + clusters[stream].size(), {-1, -1}); + std::size_t next = 0; + for (std::size_t r = 0; r < rofs[stream].size(); ++r) { + const auto& record = rofs[stream][r]; + if (record.getFirstEntry() < 0 || static_cast(record.getFirstEntry()) != next || record.getNEntries() < 0 || static_cast(record.getNEntries()) > clusters[stream].size() - next) { + throw std::runtime_error("ITS selection: cluster ROFs do not cover the stream in order"); + } + for (int i = 0; i < record.getNEntries(); ++i, ++next) { + const int layer = staggered ? stream : mapping.getLayer(clusters[stream][next].getSensorID()); + state.clusterROFs[offset + next] = {layer, static_cast(r)}; + } + } + if (next != clusters[stream].size()) { + throw std::runtime_error("ITS selection: cluster ROFs do not cover the stream"); + } + } + const auto firstOrbit = pc.services().get().firstTForbit; + const o2::InteractionRecord firstIR{0, firstOrbit}; + std::vector triggers; + if (triggerType() == 1) { + const auto input = pc.inputs().get>("its_selection_triggers"); + triggers.assign(input.begin(), input.end()); + } else if (triggerType() == 2) { + auto lastIR = firstIR; + const auto input = pc.inputs().get>("its_selection_triggers"); + for (const auto& trigger : input) { + if (trigger.getBCData() >= lastIR && trigger.getNumberOfTracklets()) { + lastIR = trigger.getBCData(); + triggers.push_back({.ir = lastIR, .data = 0}); + } + } + } + // The native estimator indexes the selection array with the trigger ROF; + // discard triggers outside this input before handing it that array. + const auto& mult = o2::its::FastMultEstConfig::Instance(); + const int layer = mult.isMultCutRequested() ? std::clamp(mult.cutMultClusLayer, 0, kLayers - 1) : state.overlap.getView().getClock(); + const auto count = rofs[staggered ? layer : 0].size(); + std::erase_if(triggers, [&](const auto& trigger) { + const auto bc = trigger.ir.differenceInBC(firstIR); + return bc < 0 || bc > std::numeric_limits::max() || static_cast(state.overlap.getLayer(layer).getROF(static_cast(bc))) >= count; + }); + Estimator estimator; + const auto selected = estimator.selectROFs(rofs, clusters, triggers, firstOrbit, staggered, state.overlap.getView(), state.mask); + LOGP(info, "ITS ROF selection: {} of {} ROFs on selection layer {}", selected, count, layer); +} + +ITSROFSelection::~ITSROFSelection() = default; + +bool ITSROFSelection::selectWithVertices(std::span seeds, std::span contributors) +{ + if (!o2::its::FastMultEstConfig::Instance().isVtxMultCutRequested()) { + return false; + } + if (seeds.size() != contributors.size()) { + throw std::runtime_error("ITS vertex selection: contributor counts do not match the vertices"); + } + std::vector vertices; + vertices.reserve(seeds.size()); + for (std::size_t v = 0; v < seeds.size(); ++v) { + auto& vertex = vertices.emplace_back(); + vertex.setNContributors(static_cast(std::min(contributors[v], std::numeric_limits::max()))); + const auto start = bcOfPS(seeds[v].timePS); + const auto end = bcOfPS(static_cast(seeds[v].timePS) + seeds[v].widthPS); + vertex.setTimeStamp({static_cast(start), static_cast(end - start)}); + } + Estimator{}.selectROFsWithVertices(vertices, mState->overlap.getView(), mState->mask); + return true; +} + +std::vector ITSROFSelection::clockMask() const +{ + const auto clock = mState->overlap.getView().getClock(); + std::vector result(mState->overlap.getLayer(clock).mNROFsTF); + const auto mask = mState->mask.getView(); + for (std::size_t r = 0; r < result.size(); ++r) { + result[r] = mask.isROFEnabled(clock, static_cast(r)); + } + return result; +} + +std::optional ITSROFSelection::filter(MeasurementView source, std::span offsets, std::span originalRows) const +{ + const auto mask = mState->mask.getView(); + const auto accepted = [&](std::uint32_t row) { + const auto cluster = source.clusterIndex[row]; + if (cluster >= mState->clusterROFs.size()) { + throw std::runtime_error("ITS selection: measurement cluster index out of range"); + } + const auto [layer, rof] = mState->clusterROFs[cluster]; + return mask.isROFEnabled(layer, rof); + }; + std::uint32_t firstRejected = 0U; + while (firstRejected < offsets.back() && accepted(firstRejected)) { + ++firstRejected; + } + if (firstRejected == offsets.back()) { + return std::nullopt; + } + std::vector rows; + for (std::uint32_t r = 0; r < firstRejected; ++r) { + rows.push_back(r); + } + for (std::uint32_t r = firstRejected + 1U; r < offsets.back(); ++r) { + if (accepted(r)) { + rows.push_back(r); + } + } + Measurements selected; + resizeMeasurementColumns(selected.columns, rows.size()); + selected.originalRows.reserve(rows.size()); + for (std::size_t i = 0; i < rows.size(); ++i) { + Measurement measurement; + measurement.apply(source[rows[i]], [](auto& to, const auto& from) -> auto& { to = from; return to; }); + writeMeasurement(selected.columns, i, measurement); + selected.originalRows.push_back(originalRows.empty() ? rows[i] : originalRows[rows[i]]); + } + selected.offsets.reserve(offsets.size()); + auto at = rows.begin(); + for (const auto offset : offsets) { + at = std::lower_bound(at, rows.end(), offset); + selected.offsets.push_back(static_cast(at - rows.begin())); + } + return selected; +} + +} // namespace o2::gtrk::workflow diff --git a/Detectors/GenericTracking/workflow/src/TrackerDevice.h b/Detectors/GenericTracking/workflow/src/TrackerDevice.h index dab7e47f39231..9d334b7f45d9d 100644 --- a/Detectors/GenericTracking/workflow/src/TrackerDevice.h +++ b/Detectors/GenericTracking/workflow/src/TrackerDevice.h @@ -31,6 +31,7 @@ #include "GenericTracking/Workflow/DataSpecs.h" #include "GenericTracking/Workflow/GPURecoTracking.h" #include "GenericTracking/Workflow/GridConfig.h" +#include "GenericTracking/Workflow/ITSROFSelection.h" #include "GenericTracking/Workflow/RecoParam.h" #include "GenericTracking/Workflow/TrackerSpec.h" #include "GenericTracking/Workflow/TrackingInterface.h" @@ -60,7 +61,7 @@ namespace o2::gtrk::workflow { /// The tracking device's products made in \a pc's outputs, empty, published on -/// every way out: (tracks, vertices, contributors, claims), DPL's vectors. +/// every way out: (tracks, vertices, contributors, claims, selectedROFs), DPL's vectors. template [[nodiscard]] auto makeTrackerProducts(o2::framework::ProcessingContext& pc) { @@ -72,7 +73,8 @@ template auto& contributors = pc.outputs().make>(Output{gDataOrigin, description(kVertexContributors), 0}); static_assert(o2::framework::is_messageable::value); auto& claims = pc.outputs().make>(Output{gDataOrigin, description(kVertexClaims), 0}); - return std::tie(tracks, vertices, contributors, claims); + auto& selectedROFs = pc.outputs().make>(Output{gDataOrigin, description(kSelectedROFs), 0}); + return std::tie(tracks, vertices, contributors, claims, selectedROFs); } /// The products' output specs. @@ -80,7 +82,7 @@ template [[nodiscard]] std::vector trackerOutputs() { using o2::framework::Lifetime; - return {{gDataOrigin, description(kTracks), 0, Lifetime::Timeframe}, {gDataOrigin, description(kVertexSeeds), 0, Lifetime::Timeframe}, {gDataOrigin, description(kVertexContributors), 0, Lifetime::Timeframe}, {gDataOrigin, description(kVertexClaims), 0, Lifetime::Timeframe}}; + return {{gDataOrigin, description(kTracks), 0, Lifetime::Timeframe}, {gDataOrigin, description(kVertexSeeds), 0, Lifetime::Timeframe}, {gDataOrigin, description(kVertexContributors), 0, Lifetime::Timeframe}, {gDataOrigin, description(kVertexClaims), 0, Lifetime::Timeframe}, {gDataOrigin, description(kSelectedROFs), 0, Lifetime::Timeframe}}; } /// One timeframe of \a Det's tracking, on \a interface's host or device @@ -92,7 +94,11 @@ void trackTimeframe(o2::framework::ProcessingContext& pc, TrackingInterface { conditions.checkUpdates(pc); const auto orbit = pc.services().get().firstTForbit; - auto [tracks, vertices, contributors, claims] = makeTrackerProducts(pc); + auto [tracks, vertices, contributors, claims, selectedROFs] = makeTrackerProducts(pc); + // Keep copied inputs alive until device work is synchronized on every + // exit, including an allocation failure during upload or vertexing. + std::optional beforeVertices, afterVertices; + std::vector trackingSeeds; const auto ended = interface.endTimeframeOnExit(); const auto input = readTimeframe(pc, withVertices); @@ -106,30 +112,81 @@ void trackTimeframe(o2::framework::ProcessingContext& pc, TrackingInterface return; } - const TimeFrame frame{measurementView(input.columns), interface.grid(), input.offsets}; - const auto found = interface.vertices(frame, input.seeds, useMC); - if (found.dropped) { - return; + std::unique_ptr selection; + auto measurements = measurementView(input.columns); + auto offsets = input.offsets; + std::span originalRows; + if constexpr (std::same_as) { + selection = std::make_unique(pc, conditions.frameTiming(), conditions.orbitsPerTimeframe()); + beforeVertices = selection->filter(measurements, offsets); + if (beforeVertices) { + measurements = measurementView(beforeVertices->columns); + offsets = beforeVertices->offsets; + originalRows = beforeVertices->originalRows; + } + if (offsets.back() == 0U) { + const auto mask = selection->clockMask(); + selectedROFs.assign(mask.begin(), mask.end()); + return; + } } - const auto tracked = interface.track(frame, found.seeds); - if (tracked.dropped) { + const TimeFrame vertexFrame{measurements, interface.grid(), offsets}; + const auto found = interface.vertices(vertexFrame, input.seeds, useMC); + if (found.dropped) { return; } - // The vertices with the tracks only: a dropped timeframe publishes nothing. + // Snapshot before a vertex-multiplicity cut can end the vertexing pass + // and rebuild the tracking input. Publish nothing if tracking then fails. + std::vector foundVertices; + std::vector foundContributors; + std::vector foundClaims; + trackingSeeds.assign(found.seeds.begin(), found.seeds.end()); if (found.found != nullptr) { - vertices.assign(found.found->vertices.begin(), found.found->vertices.end()); - contributors.assign(found.found->contributors.begin(), found.found->contributors.end()); - for (const auto& claim : found.claims) { + foundVertices = found.found->vertices; + foundContributors = found.found->contributors; + for (auto claim : found.claims) { if (claim.vertex >= 0) { - claims.push_back(claim); + if (!originalRows.empty()) { + for (auto& row : claim.rows) { + row = originalRows[row]; + } + } + foundClaims.push_back(claim); } } } + if (selection && selection->selectWithVertices(foundVertices, foundContributors)) { + afterVertices = selection->filter(measurements, offsets, originalRows); + if (afterVertices) { + interface.endTimeframe(); + measurements = measurementView(afterVertices->columns); + offsets = afterVertices->offsets; + originalRows = afterVertices->originalRows; + } + } + const TimeFrame frame{measurements, interface.grid(), offsets}; + const auto tracked = interface.track(frame, trackingSeeds); + if (tracked.dropped) { + return; + } + vertices.assign(foundVertices.begin(), foundVertices.end()); + contributors.assign(foundContributors.begin(), foundContributors.end()); + claims.assign(foundClaims.begin(), foundClaims.end()); + if (selection) { + const auto mask = selection->clockMask(); + selectedROFs.assign(mask.begin(), mask.end()); + } const auto& result = *tracked.result.tracks; tracks.reserve(result.size()); std::size_t failed = 0U; for (std::size_t track = 0; track < result.size(); ++track) { - tracks.push_back(result.read(track)); + auto record = result.read(track); + if (!originalRows.empty()) { + for (std::uint8_t c = 0; c < record.count; ++c) { + record.clusters.index[c] = originalRows[record.clusters.index[c]]; + } + } + tracks.push_back(record); failed += tracks.back().status == PropagationStatus::kOk ? 0U : 1U; } LOGP(info, "{}: the timeframe at orbit {}: {} tracks ({} with a failed refit), {} vertices, {:.3f} s tracking, {:.3f} s vertexing", Det::name, orbit, tracks.size(), failed, vertices.size(), tracked.seconds, found.seconds); @@ -200,6 +257,9 @@ o2::framework::DataProcessorSpec getTrackerSpec(bool withVertices, bool useMC) { using namespace o2::framework; auto inputs = measurementInputs(withVertices); + if constexpr (std::same_as) { + ITSROFSelection::addInputs(inputs); + } auto conditions = std::make_shared>(inputs); return DataProcessorSpec{ .name = deviceName("tracker"), @@ -226,6 +286,9 @@ class GPURecoTrackingOf final : public GPURecoTracking void appendInputs(std::vector& inputs) final { auto own = measurementInputs(false); + if constexpr (std::same_as) { + ITSROFSelection::addInputs(own); + } mConditions = std::make_unique(own); inputs.insert(inputs.end(), own.begin(), own.end()); } diff --git a/Detectors/GenericTracking/workflow/src/gtrk-its-reco-workflow.cxx b/Detectors/GenericTracking/workflow/src/gtrk-its-reco-workflow.cxx index 9e0c9e13786ff..9c3b51cf3f6b5 100644 --- a/Detectors/GenericTracking/workflow/src/gtrk-its-reco-workflow.cxx +++ b/Detectors/GenericTracking/workflow/src/gtrk-its-reco-workflow.cxx @@ -45,15 +45,21 @@ void customize(std::vector& policies) } // NOLINTNEXTLINE(misc-use-internal-linkage) -void customize(std::vector& workflowOptions) { o2::gtrk::workflow::addRecoOptions(workflowOptions, true); } +void customize(std::vector& workflowOptions) +{ + o2::gtrk::workflow::addRecoOptions(workflowOptions, true); + workflowOptions.push_back({"select-with-triggers", VariantType::String, "none", {"use triggers to prescale processed ROFs: phys, trd, none"}}); +} #include "GenericTracking/Detectors/ITS.h" #include "GenericTracking/Workflow/AdaptorSpec.h" +#include "GenericTracking/Workflow/RecoParam.h" #include "GlobalTrackingWorkflowWriters/IRFrameWriterSpec.h" #include "ITSMFTWorkflow/ClusterReaderSpec.h" #include "ITSWorkflow/TrackWriterSpec.h" +#include #include #include "Framework/runDataProcessing.h" // IWYU pragma: keep @@ -62,11 +68,19 @@ WorkflowSpec defineDataProcessing(const ConfigContext& context) { o2::conf::ConfigurableParam::updateFromString(context.options().get("configKeyValues")); + const auto trigger = context.options().get("select-with-triggers"); + if (trigger != "none") { + if (trigger != "phys" && trigger != "trd") { + throw std::invalid_argument("--select-with-triggers expects phys, trd or none"); + } + o2::conf::ConfigurableParam::updateFromString("GTRKRecoITS.selectWithTriggers=" + std::string(trigger == "phys" ? "1" : "2")); + } + namespace gw = o2::gtrk::workflow; using Det = o2::gtrk::detectors::ITS; const bool useMC = !context.options().get("disable-mc"); const bool staggered = gw::staggeredInput(context); - auto specs = gw::trackingChain(context, {o2::itsmft::getITSClusterReaderSpec(useMC, staggered, true, false)}); + auto specs = gw::trackingChain(context, {o2::itsmft::getITSClusterReaderSpec(useMC, staggered, true, o2::gtrk::recoParam().selectWithTriggers == 1)}); // The adaptor reads the tracks: after the GPU reconstruction when the tracking runs inside it. if (!context.options().get("disable-adaptor") && !gw::trackingInGPUReco(context)) { specs.push_back(gw::getAdaptorSpec(useMC, staggered)); From 8fd95bd28fda64c5eb6aef96c9f4495a1cfd52ba Mon Sep 17 00:00:00 2001 From: Maximiliano Puccio Date: Wed, 11 Nov 2026 10:58:48 +0100 Subject: [PATCH 08/11] GenericTracking: record generated-momentum QA for matched tracks --- .../study/include/GenericTracking/Study/TrackingQA.h | 11 ++++++++++- .../study/src/MeasurementSinkSpec.cxx | 1 + Detectors/GenericTracking/study/src/TrackingQA.cxx | 12 +++++++++++- 3 files changed, 22 insertions(+), 2 deletions(-) diff --git a/Detectors/GenericTracking/study/include/GenericTracking/Study/TrackingQA.h b/Detectors/GenericTracking/study/include/GenericTracking/Study/TrackingQA.h index ff6cb3ea1e004..d99e868d44ab4 100644 --- a/Detectors/GenericTracking/study/include/GenericTracking/Study/TrackingQA.h +++ b/Detectors/GenericTracking/study/include/GenericTracking/Study/TrackingQA.h @@ -25,7 +25,9 @@ /// (StageFunnel.h); its PDG codes, decay geometry, layers, /// the funnel's verdict and the other tracker's /// (`QAParticle`) -/// tracks: the fit's pT, eta and phi, its iteration, fake and clone +/// tracks: the fit's pT, eta and phi, its iteration, fake and clone; +/// fit status, generated pT/eta, origin and relative inverse-pT +/// residual for matched tracks, without requiring hit files /// /// Every efficiency follows from the two; for iteration i, `bit` = 1 << i: /// @@ -95,6 +97,7 @@ struct TrackKinematics { std::uint8_t iteration = 0U; float chi2 = 0.F; std::int32_t ndf = 0; + bool fitValid = true; ///< the refit's propagation status is kOk }; /// Iterations a particle's `foundBy` holds. @@ -231,6 +234,12 @@ struct QATrack { bool clone = false; float chi2 = 0.F; ///< the fit's; TMath::Prob(chi2, ndf) is its probability std::int32_t ndf = 0; + bool fitValid = true; + std::uint64_t mcParticle = kNoTrackParticle; ///< source/event/track identity, local to this timeframe + float ptGen = std::numeric_limits::quiet_NaN(); ///< generated momentum at production [GeV] + float etaGen = std::numeric_limits::quiet_NaN(); + std::uint8_t origin = 0xFFU; ///< ParticleOrigin; unknown for unlabelled tracks + float invPtResidual = std::numeric_limits::quiet_NaN(); ///< (1/pT - 1/pTgen) * pTgen }; /// A good track against its particle's MC crossing of its innermost diff --git a/Detectors/GenericTracking/study/src/MeasurementSinkSpec.cxx b/Detectors/GenericTracking/study/src/MeasurementSinkSpec.cxx index fffb992e06006..839042a2230ec 100644 --- a/Detectors/GenericTracking/study/src/MeasurementSinkSpec.cxx +++ b/Detectors/GenericTracking/study/src/MeasurementSinkSpec.cxx @@ -3348,6 +3348,7 @@ class MeasurementSinkTask : public Task .iteration = view.iteration[track], .chi2 = view.chi2[track], .ndf = view.ndf[track], + .fitValid = view.status[track] == PropagationStatus::kOk, }; // A good track against its particle's crossing, with the hits. const auto particle = goodParticle(rows, ours[track]); diff --git a/Detectors/GenericTracking/study/src/TrackingQA.cxx b/Detectors/GenericTracking/study/src/TrackingQA.cxx index 1568cfb075124..0346e622423da 100644 --- a/Detectors/GenericTracking/study/src/TrackingQA.cxx +++ b/Detectors/GenericTracking/study/src/TrackingQA.cxx @@ -221,8 +221,18 @@ TrackingQA trackingQA(std::span rows, const ParticleTable& truth .fake = particle == kNoTrackParticle, .chi2 = k.chi2, .ndf = k.ndf, + .fitValid = k.fitValid, + .mcParticle = particle, }; if (!record.fake) { + if (const auto generated = truth.find(particle); generated != truth.end()) { + record.ptGen = generated->second.pt; + record.etaGen = generated->second.eta; + record.origin = static_cast(generated->second.origin); + if (record.fitValid && std::isfinite(k.pt) && k.pt > 0.F && std::isfinite(record.ptGen) && record.ptGen > 0.F) { + record.invPtResidual = static_cast((static_cast(record.ptGen) / static_cast(k.pt)) - 1.); + } + } auto& mask = foundBy[particle]; record.clone = mask != 0U; mask |= std::uint64_t{1} << k.iteration; @@ -280,7 +290,7 @@ void stream(o2::utils::TreeStreamRedirector& output, const TrackingQA& qa) output << "particles" << "pt=" << p.pt << "eta=" << p.eta << "phi=" << p.phi << "origin=" << p.origin << "foundBy=" << p.foundBy << "seededBy=" << p.seededBy << "pdg=" << p.pdg << "motherPdg=" << p.motherPdg << "vr=" << e.vr << "vz=" << p.vz << "d0=" << e.d0 << "z0=" << e.z0 << "openingAngle=" << e.openingAngle << "motherMiss=" << e.motherMiss << "layers=" << p.layers << "lostAt=" << e.lostAt << "pairLoss=" << e.pairLoss << "roadLoss=" << e.roadLoss << "fitLoss=" << e.fitLoss << "fitRoadLayers=" << e.fitRoadLayers << "roadLosses=" << e.roadLosses << "fitLosses=" << e.fitLosses << "raisedFit=" << e.raisedFit << "foundByO2=" << p.foundByO2 << "findable=" << p.findable << "outwardLayers=" << e.outwardLayers << "\n"; } for (const auto& t : qa.tracks) { - output << "tracks" << "pt=" << t.pt << "eta=" << t.eta << "phi=" << t.phi << "iteration=" << t.iteration << "fake=" << t.fake << "clone=" << t.clone << "chi2=" << t.chi2 << "ndf=" << t.ndf << "\n"; + output << "tracks" << "pt=" << t.pt << "eta=" << t.eta << "phi=" << t.phi << "iteration=" << t.iteration << "fake=" << t.fake << "clone=" << t.clone << "chi2=" << t.chi2 << "ndf=" << t.ndf << "fitValid=" << t.fitValid << "mcParticle=" << t.mcParticle << "ptGen=" << t.ptGen << "etaGen=" << t.etaGen << "origin=" << t.origin << "dInvPtRel=" << t.invPtResidual << "\n"; } } From 841ff3694ccab1956f9f1d13f9737abf42c17e01 Mon Sep 17 00:00:00 2001 From: Maximiliano Puccio Date: Wed, 11 Nov 2026 11:20:06 +0100 Subject: [PATCH 09/11] GenericTracking: compare native momentum QA and record full-length coverage --- .../GenericTracking/Study/Alice/Tracks.h | 8 +++- .../GenericTracking/Study/TrackingQA.h | 14 +++++-- .../study/src/MeasurementSinkSpec.cxx | 22 +++++++++- .../GenericTracking/study/src/TrackingQA.cxx | 17 ++++++-- .../GenericTracking/study/src/Tracks.cxx | 41 +++++++++++++++++++ 5 files changed, 94 insertions(+), 8 deletions(-) diff --git a/Detectors/GenericTracking/study/include/GenericTracking/Study/Alice/Tracks.h b/Detectors/GenericTracking/study/include/GenericTracking/Study/Alice/Tracks.h index 04a434bb60bc7..e5ecdc29ffa0b 100644 --- a/Detectors/GenericTracking/study/include/GenericTracking/Study/Alice/Tracks.h +++ b/Detectors/GenericTracking/study/include/GenericTracking/Study/Alice/Tracks.h @@ -11,7 +11,7 @@ /// \file Tracks.h /// \brief O2's ITS and MFT tracks, as the clusters they hold, their MC labels' -/// fake flags and, for the ITS, their states at the innermost cluster. +/// fake flags, fitted kinematics and ITS innermost-cluster states. /// /// All a comparison with our tracks needs, because both trackers read the /// same clusters in the same order (Clusters.h): a track is its list of O2 @@ -34,6 +34,8 @@ class TTree; namespace o2::gtrk::study { +struct TrackKinematics; + /// One timeframe's tracks, each as its cluster indices. using TrackClusters = std::vector>; @@ -56,6 +58,10 @@ class TrackFile /// Whether each track's own MC label is flagged fake, as O2's tracker set it. [[nodiscard]] std::vector fakeFlags(long long entry); + /// Fitted inward kinematics in the same order as read(). Native outputs + /// contain accepted tracks; fitValid checks finite, positive momentum. + [[nodiscard]] std::vector kinematics(long long entry); + /// Every ITS track's state at its innermost cluster. Throws /// std::runtime_error for the MFT, whose tracks are parameterised in z. [[nodiscard]] std::vector innerStates(long long entry); diff --git a/Detectors/GenericTracking/study/include/GenericTracking/Study/TrackingQA.h b/Detectors/GenericTracking/study/include/GenericTracking/Study/TrackingQA.h index d99e868d44ab4..57dcc235b28a2 100644 --- a/Detectors/GenericTracking/study/include/GenericTracking/Study/TrackingQA.h +++ b/Detectors/GenericTracking/study/include/GenericTracking/Study/TrackingQA.h @@ -150,7 +150,8 @@ struct ParticleExtras { /// own part, a FitLoss, kOther if the fit takes it; kQANotObserved /// otherwise. std::uint8_t raisedFit = kQANotObserved; - std::uint8_t outwardLayers = 0U; ///< the layers of its MC hits' outward leg (`outwardLayers`) + std::uint8_t outwardLayers = 0U; ///< the layers of its MC hits' outward leg (`outwardLayers`) + std::uint64_t expectedLayerMask = 0U; ///< full-length reference: detector layers, or ALICE3 outward MC traversal }; /// A findable particle of the `particles` tree. @@ -235,6 +236,8 @@ struct QATrack { float chi2 = 0.F; ///< the fit's; TMath::Prob(chi2, ndf) is its probability std::int32_t ndf = 0; bool fitValid = true; + std::uint64_t expectedLayerMask = 0U; ///< zero when a full-length reference is unavailable + std::uint64_t layerMask = 0U; ///< distinct reconstructed detector layers std::uint64_t mcParticle = kNoTrackParticle; ///< source/event/track identity, local to this timeframe float ptGen = std::numeric_limits::quiet_NaN(); ///< generated momentum at production [GeV] float etaGen = std::numeric_limits::quiet_NaN(); @@ -352,13 +355,18 @@ using SeededBy = std::unordered_map; /// \a qa's records into \a output's `particles` and `tracks` trees. void stream(o2::utils::TreeStreamRedirector& output, const TrackingQA& qa); +/// Write the same track schema for a comparison tracker, without duplicating +/// the generated-particle tree. Truth and clone rules are those of trackingQA(). +void streamTracks(o2::utils::TreeStreamRedirector& output, std::span tracks, const char* tree = "tracks"); + /// \a entries into \a output's `resolution` tree. void streamResolution(o2::utils::TreeStreamRedirector& output, std::span entries); /// The run's one entry of the `run` tree: what the macro's axes need, /// \a settings' findable rules (`minLayers`, `minPtGeV`, `minOutwardLayers`) -/// for its definitions, and the run's MC \a events. -void streamRun(o2::utils::TreeStreamRedirector& output, int iterations, EtaRange eta, const ComparisonSettings& settings, std::size_t events); +/// for its definitions, and the run's MC \a events. A nonzero fullLengthMask +/// requires all detector layers; zero denotes per-particle MC traversal masks. +void streamRun(o2::utils::TreeStreamRedirector& output, int iterations, EtaRange eta, const ComparisonSettings& settings, std::size_t events, std::uint64_t fullLengthMask = 0U); } // namespace o2::gtrk::study diff --git a/Detectors/GenericTracking/study/src/MeasurementSinkSpec.cxx b/Detectors/GenericTracking/study/src/MeasurementSinkSpec.cxx index 839042a2230ec..3449b52350cd4 100644 --- a/Detectors/GenericTracking/study/src/MeasurementSinkSpec.cxx +++ b/Detectors/GenericTracking/study/src/MeasurementSinkSpec.cxx @@ -3238,7 +3238,7 @@ class MeasurementSinkTask : public Task const std::string name = mQAOutput->GetFile()->GetName(); if (mInterface.configured()) { // The MC events: the signal source's, as the truth is. - streamRun(*mQAOutput, workflow::iterationCount(recoParam()), etaReach(PlacedEnvelopes::all()), comparisonSettings(), mKinematics ? mKinematics->getNEvents(0) : 0U); + streamRun(*mQAOutput, workflow::iterationCount(recoParam()), etaReach(PlacedEnvelopes::all()), comparisonSettings(), mKinematics ? mKinematics->getNEvents(0) : 0U, std::is_same_v ? 0U : (std::uint64_t{1} << Detector::nLayers) - 1U); } mQAOutput->Close(); mQAOutput.reset(); @@ -3399,6 +3399,7 @@ class MeasurementSinkTask : public Task // Each signal particle's outward leg in its MC hits (TrackingQA.h), its // event's segments ordered by track once. std::unordered_map outward; + std::unordered_map outwardMasks; std::unordered_map> legMomenta; // --species-oracle if (mHits) { std::unordered_map> byEvent; @@ -3429,6 +3430,13 @@ class MeasurementSinkTask : public Task const auto& truth = particles.at(id); const auto leg = outwardLeg(crossings, {static_cast(truth.vx), static_cast(truth.vy), static_cast(truth.vz)}); outward[id] = outwardLayers(leg); + if constexpr (std::is_same_v) { + if (outward[id] >= comparisonSettings().minLayers) { + for (const auto& crossing : leg) { + outwardMasks[id] |= std::uint64_t{1} << crossing.layer; + } + } + } if (mSpeciesOracle && (std::abs(truth.pdg) == 11 || (std::abs(truth.pdg) == 211 && id % 16U == 0U))) { legMomenta[id] = momentaAlongLeg(leg); } @@ -3443,6 +3451,13 @@ class MeasurementSinkTask : public Task } const auto& helix = label->second.helix; auto& extra = extras[id]; + if constexpr (std::is_same_v) { + if (const auto mask = outwardMasks.find(id); mask != outwardMasks.end()) { + extra.expectedLayerMask = mask->second; + } + } else { + extra.expectedLayerMask = (std::uint64_t{1} << Detector::nLayers) - 1U; + } if (const auto leg = outward.find(id); leg != outward.end()) { extra.outwardLayers = leg->second; } @@ -3481,6 +3496,11 @@ class MeasurementSinkTask : public Task } } stream(*mQAOutput, trackingQA(rows, particles, ours, kinematics, funnel != nullptr ? funnel->seededBy() : SeededBy{}, extras, theirs, workflow::iterationCount(recoParam()), comparisonSettings())); + if (theirs != nullptr) { + const auto nativeKinematics = mO2Tracks->kinematics(static_cast(mTimeframe)); + const auto nativeQA = trackingQA(rows, particles, *theirs, nativeKinematics, {}, extras, nullptr, 1, comparisonSettings()); + streamTracks(*mQAOutput, nativeQA.tracks, "tracksO2"); + } streamResolution(*mQAOutput, resolution); streamOracle(*mQAOutput, oracle); if (!legMomenta.empty()) { diff --git a/Detectors/GenericTracking/study/src/TrackingQA.cxx b/Detectors/GenericTracking/study/src/TrackingQA.cxx index 0346e622423da..518958c82a8fa 100644 --- a/Detectors/GenericTracking/study/src/TrackingQA.cxx +++ b/Detectors/GenericTracking/study/src/TrackingQA.cxx @@ -224,7 +224,13 @@ TrackingQA trackingQA(std::span rows, const ParticleTable& truth .fitValid = k.fitValid, .mcParticle = particle, }; + for (const auto row : tracks[track]) { + record.layerMask |= std::uint64_t{1} << rows[row].layer; + } if (!record.fake) { + if (const auto extra = extras.find(particle); extra != extras.end()) { + record.expectedLayerMask = extra->second.expectedLayerMask; + } if (const auto generated = truth.find(particle); generated != truth.end()) { record.ptGen = generated->second.pt; record.etaGen = generated->second.eta; @@ -289,8 +295,13 @@ void stream(o2::utils::TreeStreamRedirector& output, const TrackingQA& qa) const auto& e = p.extras; output << "particles" << "pt=" << p.pt << "eta=" << p.eta << "phi=" << p.phi << "origin=" << p.origin << "foundBy=" << p.foundBy << "seededBy=" << p.seededBy << "pdg=" << p.pdg << "motherPdg=" << p.motherPdg << "vr=" << e.vr << "vz=" << p.vz << "d0=" << e.d0 << "z0=" << e.z0 << "openingAngle=" << e.openingAngle << "motherMiss=" << e.motherMiss << "layers=" << p.layers << "lostAt=" << e.lostAt << "pairLoss=" << e.pairLoss << "roadLoss=" << e.roadLoss << "fitLoss=" << e.fitLoss << "fitRoadLayers=" << e.fitRoadLayers << "roadLosses=" << e.roadLosses << "fitLosses=" << e.fitLosses << "raisedFit=" << e.raisedFit << "foundByO2=" << p.foundByO2 << "findable=" << p.findable << "outwardLayers=" << e.outwardLayers << "\n"; } - for (const auto& t : qa.tracks) { - output << "tracks" << "pt=" << t.pt << "eta=" << t.eta << "phi=" << t.phi << "iteration=" << t.iteration << "fake=" << t.fake << "clone=" << t.clone << "chi2=" << t.chi2 << "ndf=" << t.ndf << "fitValid=" << t.fitValid << "mcParticle=" << t.mcParticle << "ptGen=" << t.ptGen << "etaGen=" << t.etaGen << "origin=" << t.origin << "dInvPtRel=" << t.invPtResidual << "\n"; + streamTracks(output, qa.tracks); +} + +void streamTracks(o2::utils::TreeStreamRedirector& output, std::span tracks, const char* tree) +{ + for (const auto& t : tracks) { + output << tree << "pt=" << t.pt << "eta=" << t.eta << "phi=" << t.phi << "iteration=" << t.iteration << "fake=" << t.fake << "clone=" << t.clone << "chi2=" << t.chi2 << "ndf=" << t.ndf << "expectedLayerMask=" << t.expectedLayerMask << "layerMask=" << t.layerMask << "fitValid=" << t.fitValid << "mcParticle=" << t.mcParticle << "ptGen=" << t.ptGen << "etaGen=" << t.etaGen << "origin=" << t.origin << "dInvPtRel=" << t.invPtResidual << "\n"; } } @@ -316,6 +327,6 @@ void streamMomenta(o2::utils::TreeStreamRedirector& output, std::span #include +#include #include #include +#include namespace o2::gtrk::study { @@ -63,6 +66,34 @@ TrackClusters readEntry(TTree& tree, long long entry, const char* trackBranch, c return result; } +template +std::vector readKinematics(TTree& tree, long long entry, const char* branch) +{ + std::vector tracks; + auto* pointer = &tracks; + tree.SetBranchAddress(branch, &pointer); + tree.GetEntry(entry); + tree.ResetBranchAddresses(); + std::vector result; + result.reserve(tracks.size()); + for (const auto& track : tracks) { + const auto pt = static_cast(track.getPt()); + const auto eta = static_cast(track.getEta()); + const auto phi = static_cast(track.getPhi()); + float chi2 = 0.F; + int ndf = 0; + if constexpr (std::is_same_v) { + chi2 = track.getChi2(); + ndf = 2 * track.getNumberOfClusters() - 5; + } else { + chi2 = static_cast(track.getTrackChi2()); + ndf = track.getNDF(); + } + result.push_back({.pt = pt, .eta = eta, .phi = phi, .chi2 = chi2, .ndf = ndf, .fitValid = std::isfinite(pt) && pt > 0.F && std::isfinite(eta) && std::isfinite(phi)}); + } + return result; +} + } // namespace TrackFile::TrackFile(const std::string& path, Chips chips) : mChips(chips), mFile(TFile::Open(path.c_str())) @@ -107,6 +138,16 @@ std::vector TrackFile::fakeFlags(long long entry) return result; } +std::vector TrackFile::kinematics(long long entry) +{ + if (entry >= mTree->GetEntries()) { + return {}; + } + return mChips == Chips::kITS + ? readKinematics(*mTree, entry, "ITSTrack") + : readKinematics(*mTree, entry, "MFTTrack"); +} + std::vector TrackFile::innerStates(long long entry) { if (mChips != Chips::kITS) { From a5fdc56b02cbac04d10bf81ba7cf4be4d295caad Mon Sep 17 00:00:00 2001 From: Maximiliano Puccio Date: Wed, 11 Nov 2026 12:27:34 +0100 Subject: [PATCH 10/11] GenericTracking: check material against aligned simulation geometry --- .../GenericTracking/study/tools/checkMatLUT.cxx | 13 ++++++++++--- 1 file changed, 10 insertions(+), 3 deletions(-) diff --git a/Detectors/GenericTracking/study/tools/checkMatLUT.cxx b/Detectors/GenericTracking/study/tools/checkMatLUT.cxx index ae63215338ee3..f1f78a5193216 100644 --- a/Detectors/GenericTracking/study/tools/checkMatLUT.cxx +++ b/Detectors/GenericTracking/study/tools/checkMatLUT.cxx @@ -534,7 +534,14 @@ void write(const std::string& path, const std::vector& steps) template int runFor(const Options& options) { - o2::base::GeometryManager::loadGeometry(options.geometry.empty() ? options.sim + "/sgn_geometry.root" : options.geometry); + // O2 exports the ideal geometry before applying alignment, then transports + // with the aligned geometry (O2MCApplicationBase::MisalignGeometry). + auto geometryPath = options.geometry; + if (geometryPath.empty()) { + const auto aligned = options.sim + "/sgn_geometry-aligned.root"; + geometryPath = std::filesystem::exists(aligned) ? aligned : options.sim + "/sgn_geometry.root"; + } + o2::base::GeometryManager::loadGeometry(geometryPath); const auto geometry = std::make_unique>(o2::gtrk::alice::placeChecked()); const std::unique_ptr lut{o2::base::MatLayerCylSet::loadFromFile(options.lut)}; if (!lut) { @@ -577,9 +584,9 @@ int run(int argc, char** argv) description.add_options() // ("help,h", "this message") // ("detector,d", bpo::value(&options.detector)->default_value(options.detector), "ITS, MFT or TRKFT3") // - ("sim", bpo::value(&options.sim), "the simulation's directory: sgn_Hits*.root, sgn_Kine.root, sgn_geometry.root") // + ("sim", bpo::value(&options.sim), "the simulation's directory: sgn_Hits*.root, sgn_Kine.root, sgn_geometry[-aligned].root") // ("lut", bpo::value(&options.lut), "the material LUT") // - ("geometry", bpo::value(&options.geometry), "the geometry, if not /sgn_geometry.root (the hits' own transport geometry)") // + ("geometry", bpo::value(&options.geometry), "override geometry; default /sgn_geometry-aligned.root if present, otherwise sgn_geometry.root") // ("output,o", bpo::value(&options.output), "write the steps to this ROOT file") // ("events,n", bpo::value(&options.events)->default_value(options.events), "the first MC events") // ("min-pt", bpo::value(&options.minPt)->default_value(options.minPt), "the particles' pT floor [GeV]") // From 7c06f0d778166a09ced87ffa0c6df3120d82f462 Mon Sep 17 00:00:00 2001 From: Maximiliano Puccio Date: Thu, 12 Nov 2026 09:54:45 +0100 Subject: [PATCH 11/11] GenericTracking: add independent references to track fitting Use the whole-chain Karimaki circle and longitudinal seed as an independent outward reference. Option 2 starts a new independent inward reference from the outward fitted parameters. Anchor each reference to cluster positions while retaining its transported direction and momentum within the pass. Use the reference in the initial road fit and reuse its outward fitted state for selected-track and beam-line inward fits. Follower refits use the same per-pass policy. No per-surface reference cache is needed. Keep ordinary fitting as the default and preserve the covariance reset, material model, mass and cuts. Expose reference configuration and fit diagnostics, and exercise the modes in the existing physical-truth and road regression tests. Validation: ten CPU regression modules pass; options 0 and 1 produce byte-identical first-timeframe outputs for ITS, MFT and TRK. CUDA fitter memcheck and initcheck are clean, including a partial final batch; full GPU tracking remains unvalidated. --- .../study/tools/CompleteAnatomy.h | 6 +- .../GenericTracking/study/tools/TrackBench.h | 59 +++ .../study/tools/trackBench.cxx | 1 + .../tracking/GPU/cuda/DeviceTracker.cuh | 2 +- .../Tracking/CovarianceTransport.h | 16 +- .../include/GenericTracking/Tracking/Follow.h | 10 +- .../GenericTracking/Tracking/Iteration.h | 4 + .../GenericTracking/Tracking/RoadFit.h | 2 +- .../GenericTracking/Tracking/TrackFitter.h | 194 ++++++-- .../GenericTracking/Tracking/TrackRefit.h | 14 +- .../test/device/testDeviceTracker.cxx | 74 ++- .../tracking/test/testRoadFit.cxx | 178 +++---- .../tracking/test/testTrackFitter.cxx | 437 +++++++++--------- .../GenericTracking/Workflow/GridConfig.h | 4 + .../GenericTracking/Workflow/RecoParam.h | 6 + 15 files changed, 637 insertions(+), 370 deletions(-) diff --git a/Detectors/GenericTracking/study/tools/CompleteAnatomy.h b/Detectors/GenericTracking/study/tools/CompleteAnatomy.h index 921306a73b6b2..13389c6b4a4ef 100644 --- a/Detectors/GenericTracking/study/tools/CompleteAnatomy.h +++ b/Detectors/GenericTracking/study/tools/CompleteAnatomy.h @@ -4342,7 +4342,11 @@ class CompleteAnatomy // The road fit, its last step nominal and allowed. const RowChain chain{.frame = &frame, .rows = road.data()}; const auto count = static_cast(road.size()); - const auto plain = fitOutward(propagator, geometry, material, chain, count, parameters.roadChi2MaxByNdf, nominalBudget); + auto plainParameters = parameters; + plainParameters.kinkScale2 = 1.F; + plainParameters.maxKinks = nominalBudget.max; + plainParameters.kinkChi2Min = nominalBudget.gate; + const auto plain = fitRoad(plainParameters, propagator, geometry, material, chain, count); const auto allowed = fitRoad(parameters, propagator, geometry, material, chain, count); if (plain.status != RoadFitStatus::kAccepted && allowed.status == RoadFitStatus::kAccepted) { need |= kLastStep; diff --git a/Detectors/GenericTracking/study/tools/TrackBench.h b/Detectors/GenericTracking/study/tools/TrackBench.h index 46bc3abc5bf73..999976c646eca 100644 --- a/Detectors/GenericTracking/study/tools/TrackBench.h +++ b/Detectors/GenericTracking/study/tools/TrackBench.h @@ -51,6 +51,7 @@ #include "DetectorsBase/MatLayerCylSet.h" #include +#include #include "FitInputs.h" #include @@ -139,6 +140,7 @@ struct Options { // NOLINT(clang-analyzer-optin.performance.Padding) std::string trace; std::string json; std::string tracksOut; + std::string fitResults; std::string clonesOut; ///< --clones-out std::string rowsOut; ///< --rows-out std::string fakesOut; ///< --fakes-out @@ -546,6 +548,57 @@ void writeTracks(const std::string& path, const typename TrackRefit::Tracks } } +/// Per-track fit output for comparisons, outside the measured tracking loop. +/// Particle IDs include source/event/track; the caller supplies the TF identity. +template +void writeFitResults(const std::string& path, const typename TrackRefit::Tracks& tracks, const TimeFrame& frame, const study::TrackingTruth* truth) +{ + std::ofstream out{path}; + if (!out) { + throw std::runtime_error("cannot write " + path); + } + out << "track,iteration,status,count,layerMask,chi2,ndf,kinks,good,particle,ptGen,etaGen,physicalPrimary"; + for (int i = 0; i < kNBoundParams; ++i) { + out << ",par" << i; + } + for (int i = 0; i < kBoundCovSize; ++i) { + out << ",cov" << i; + } + out << ",clusters\n" << std::setprecision(std::numeric_limits::max_digits10); + const auto view = tracks.view(); + const auto columns = frame.measurements(); + for (std::size_t track = 0; track < tracks.size(); ++track) { + const auto rows = std::span{view.clusters[track].index, view.count[track]}; + const auto good = truth ? study::goodParticle(truth->rows, rows) : study::kNoTrackParticle; + const auto particle = good != study::kNoTrackParticle ? good : truth ? study::majorityParticle(truth->rows, rows).particle : study::kNoTrackParticle; + const study::ParticleTruth* mc = nullptr; + std::uint64_t layerMask = 0U; + if (truth) { + const auto it = truth->particles.find(particle); + if (it != truth->particles.end()) { + mc = &it->second; + } + for (const auto row : rows) { + layerMask |= std::uint64_t{1} << truth->rows[row].layer; + } + } + out << track << ',' << static_cast(view.iteration[track]) << ',' << static_cast(view.status[track]) << ',' << rows.size() << ',' << layerMask + << ',' << view.chi2[track] << ',' << view.ndf[track] << ',' << static_cast(view.kinks[track]) << ',' << (good != study::kNoTrackParticle) << ',' << particle + << ',' << (mc ? mc->pt : 0.F) << ',' << (mc ? mc->eta : 0.F) << ',' << (mc && mc->origin == study::ParticleOrigin::kPhysicalPrimary); + for (const auto parameter : view.state[track].par) { + out << ',' << parameter; + } + for (const auto covariance : view.state[track].cov) { + out << ',' << covariance; + } + out << ','; + for (const auto row : rows) { + out << columns.clusterIndex[row] << ';'; + } + out << '\n'; + } +} + /// The checkpoint file of --save-state and --resume-state: a header, then /// each part raw -- the first iterations' settings, which --resume-state checks against its own, their summaries, the used mask and /// the tracks. The header's fingerprints hold what else the iterations read: @@ -1920,6 +1973,9 @@ int track(const Options& options) // NOLINT(readability-function-size): the benc if (!options.tracksOut.empty() && last) { writeTracks(options.tracksOut, deviceTracks, frame); } + if (!options.fitResults.empty() && last) { + writeFitResults(options.fitResults, deviceTracks, frame, truth ? &*truth : nullptr); + } std::cout << std::format("\n=== run {}{}: tracks, {}, {}: {} in {:.4f} s (plan {:.4f} s, upload {:.4f} s, iterations", repetition, repetition == 0U ? " (cold)" : "", label, options.backend, deviceTracks.size(), planSeconds + deviceSeconds, planSeconds, deviceRun.uploadSeconds); for (const double seconds : deviceRun.iterationSeconds) { std::cout << std::format(" {:.4f}", seconds); @@ -2042,6 +2098,9 @@ int track(const Options& options) // NOLINT(readability-function-size): the benc if (!options.tracksOut.empty() && last) { writeTracks(options.tracksOut, *result.tracks, frame); } + if (!options.fitResults.empty() && last) { + writeFitResults(options.fitResults, *result.tracks, frame, truth ? &*truth : nullptr); + } if (truth && last) { writeTruthStudies(options, *result.tracks, *truth, comparison, frame, tracking, *material, propagator, studies.followed, [&](std::size_t iteration) -> const IterationParameters& { return tracker.iteration(iteration); }, [&](std::uint64_t id, std::size_t iteration) { return funnel ? funnel->stages(id, iteration) : std::nullopt; }, [&](std::size_t iteration, std::span rows) { return funnel ? funnel->recordChi2(iteration, rows) : std::nullopt; }); } diff --git a/Detectors/GenericTracking/study/tools/trackBench.cxx b/Detectors/GenericTracking/study/tools/trackBench.cxx index 399304e0c0e79..97eb2e364239a 100644 --- a/Detectors/GenericTracking/study/tools/trackBench.cxx +++ b/Detectors/GenericTracking/study/tools/trackBench.cxx @@ -166,6 +166,7 @@ static int run(int argc, char** argv) ("fit-inputs", bpo::value(&options.fitInputs)->default_value(""), "write the last run's triplet and neighbour fit inputs here, for o2-gtrk-fit-bench") // ("trace", bpo::value(&options.trace)->default_value(""), "write the last run's zones and chunks here, as a Chrome trace") // ("json", bpo::value(&options.json)->default_value(""), "append one JSON line per run here") // + ("fit-results", bpo::value(&options.fitResults)->default_value(""), "write the last run's fit parameters, covariance, cluster identities and MC labels as CSV, outside timing") // ("tracks-out", bpo::value(&options.tracksOut)->default_value(""), "write the last run's tracks here, one line each: iteration, clusters, chi2") // ("follow-steps", bpo::value(&options.followSteps)->default_value(""), "write the follow's steps on correct roads against the truth here, replayed after each iteration") // diff --git a/Detectors/GenericTracking/tracking/GPU/cuda/DeviceTracker.cuh b/Detectors/GenericTracking/tracking/GPU/cuda/DeviceTracker.cuh index f816f226c8e81..fb3fed5735999 100644 --- a/Detectors/GenericTracking/tracking/GPU/cuda/DeviceTracker.cuh +++ b/Detectors/GenericTracking/tracking/GPU/cuda/DeviceTracker.cuh @@ -1114,7 +1114,7 @@ class DeviceTrackerImpl final : public DeviceTracker if (!tracks.allocate(mArena, trackCount) || !proposals.allocate(mArena, followCount) || !kinkSteps.allocate(mArena, kinkStepCount) || !followRows.allocate(mArena, followCount) || !followCounts.allocate(mArena, followCount) || !endStarts.allocate(mArena, hostParameters.followEnd ? trackCount : 0U)) { return StageStatus::kOutOfMemory; } - const TrackRefitElement refit{mFrame.data(), mGeometry.data(), mMaterial.data(), mPropagator, records.clusters, records.count, records.outwardState, records.chi2, std::as_const(selected.selected).data(), tracks.view(), 0U, static_cast(mIteration), kinkBudget(hostParameters)}; + const TrackRefitElement refit{mFrame.data(), mGeometry.data(), mMaterial.data(), mPropagator, records.clusters, records.count, records.outwardState, records.chi2, std::as_const(selected.selected).data(), tracks.view(), 0U, static_cast(mIteration), kinkBudget(hostParameters), hostParameters.fitReference}; algorithms::forEach(*mRuntime, mStream, trackCount, refit); // The follow, the host's sequence diff --git a/Detectors/GenericTracking/tracking/include/GenericTracking/Tracking/CovarianceTransport.h b/Detectors/GenericTracking/tracking/include/GenericTracking/Tracking/CovarianceTransport.h index 30a4ebcc86ee9..a1532ef255fe4 100644 --- a/Detectors/GenericTracking/tracking/include/GenericTracking/Tracking/CovarianceTransport.h +++ b/Detectors/GenericTracking/tracking/include/GenericTracking/Tracking/CovarianceTransport.h @@ -182,17 +182,11 @@ GPUhdi() bool transportFactors(const SolenoidJacobian& jac, const StepNoise& noi /// The reference, not the filtered state. A Kalman filter that linearises /// about its own estimate is an *extended* filter, and its linearisation point /// carries whatever error the filter currently has; linearising about a -/// reference trajectory - O2 does this with the `linRef` overload of -/// `PropagatorImpl::PropagateToXBxByBz` - keeps F and Q independent of the -/// measurements. For a two-pass fit that is a correctness property rather than -/// a refinement: the outward and inward passes then see *identical* F and Q -/// for the same physical step, which is what lets them be combined. -/// -/// Unlike O2 we do not also transport the state through the reference. O2 -/// writes `state = ref + F (state - ref)` because its exact propagation is a -/// stepped field integration and doing it once for the reference is a real -/// saving. Ours is a closed form, so linearising the state would trade -/// accuracy for a speed-up we do not need. +/// reference trajectory allows the fitter to evaluate transport independently +/// of its measurement updates. The caller chooses the linearisation state; +/// this routine returns the exact propagated reference and its Jacobian. +/// The reference filter uses them for the affine prediction +/// f(reference) + F (state - reference). /// /// The derivation, for a uniform Bz. With `b = Bz * B2C` and `kp = (q/p) b`, /// the turn angle is `psi = kp * s` - **independent of theta**, because the diff --git a/Detectors/GenericTracking/tracking/include/GenericTracking/Tracking/Follow.h b/Detectors/GenericTracking/tracking/include/GenericTracking/Tracking/Follow.h index 7e562f5b7ce62..2144df0a489bf 100644 --- a/Detectors/GenericTracking/tracking/include/GenericTracking/Tracking/Follow.h +++ b/Detectors/GenericTracking/tracking/include/GenericTracking/Tracking/Follow.h @@ -1429,7 +1429,7 @@ class FollowStage // Every lane holds the best so far, each passing row shared as // it comes. const auto one = follow::shared(group, candidate, from, branch.state); - if (++passing == 1U || follow::beats(one.row, one.chi2, {.row = child.row, .chi2 = child.chi2})) { + if (++passing == 1U || follow::beats(one.row, one.chi2, follow::Witness{.row = child.row, .chi2 = child.chi2})) { child = one; } continue; @@ -1513,7 +1513,7 @@ class FollowStage const bool confirmed = confirmation(branch, box, direction, taken, next, witness, WorkCounter::kFollowUpdates, group); const Scalar key = confirmed ? taken.chi2 + witness : taken.chi2; observe.lookaheadConfirmation(taken.row, taken.chi2, next, witness); - const bool better = best.child.row == ClusterList::kNoCluster || (confirmed != best.confirmed ? confirmed : follow::beats(taken.row, key, {.row = best.child.row, .chi2 = best.key})); + const bool better = best.child.row == ClusterList::kNoCluster || (confirmed != best.confirmed ? confirmed : follow::beats(taken.row, key, follow::Witness{.row = best.child.row, .chi2 = best.key})); if (better) { best.child = taken; best.confirmed = confirmed; @@ -1557,7 +1557,7 @@ class FollowStage continue; } const Scalar score = taken.chi2 + witness; - if (follow::beats(taken.row, score, {.row = child.row, .chi2 = bestScore})) { + if (follow::beats(taken.row, score, follow::Witness{.row = child.row, .chi2 = bestScore})) { child = taken; bestScore = score; onNext = next != ClusterList::kNoCluster; @@ -1909,7 +1909,7 @@ struct Refit { /// followed again with `followEarlier`. [[nodiscard]] GPUhdi() GPUAlwaysInline() bool accept(std::uint32_t track, const ClusterList& list, std::uint8_t n, std::uint8_t outward, std::uint8_t inward) const { - const auto fit = fitMeasurements(propagator, *geometry, *material, RowChain{.frame = frame, .rows = list.index}, static_cast(n), kinkBudget(*parameters)); + const auto fit = fitMeasurementsReference(propagator, *geometry, *material, RowChain{.frame = frame, .rows = list.index}, static_cast(n), parameters->fitReference, kinkBudget(*parameters)); if (!passesRoadChi2(*parameters, fit)) { return false; } @@ -2013,7 +2013,7 @@ struct TrimAmbiguous { rows.index[count - 1U] = ClusterList::kNoCluster; const auto shorter = static_cast(count - 1U); countWork(WorkCounter::kFollowTrimFits); - const auto fit = fitMeasurements(propagator, *geometry, *material, RowChain{.frame = frame, .rows = rows.index}, static_cast(shorter), kinkBudget(*parameters)); + const auto fit = fitMeasurementsReference(propagator, *geometry, *material, RowChain{.frame = frame, .rows = rows.index}, static_cast(shorter), parameters->fitReference, kinkBudget(*parameters)); if (!passesRoadChi2(*parameters, fit) || !keptAfterFollow(*parameters, shorter, fit.state)) { return; } diff --git a/Detectors/GenericTracking/tracking/include/GenericTracking/Tracking/Iteration.h b/Detectors/GenericTracking/tracking/include/GenericTracking/Tracking/Iteration.h index e2d8d711b4f5f..f4e4a546426cb 100644 --- a/Detectors/GenericTracking/tracking/include/GenericTracking/Tracking/Iteration.h +++ b/Detectors/GenericTracking/tracking/include/GenericTracking/Tracking/Iteration.h @@ -33,6 +33,7 @@ #include "GenericTracking/Base/LayerLinks.h" #include "GenericTracking/Base/LayerSpec.h" #include "GenericTracking/Tracking/Definitions.h" +#include "GenericTracking/Tracking/TrackFitter.h" #include "GenericTracking/Tracking/LinkHoles.h" #include "GenericTracking/Tracking/MathUtils.h" #include "GenericTracking/Tracking/Measurement.h" @@ -70,6 +71,7 @@ struct IterationSettings { int minClusters = 4; /// The kinks a fit pass may take; 0 for none. int maxKinks = 1; + FitReference fitReference = FitReference::kFilteredState; /// The triplet and neighbour cuts' kink allowance: the factor on theta0 at pT_min a triplet's middle /// point may have scattered by, beyond the nominal cuts, to be kept as a /// kinked triplet; 1 for none. @@ -294,6 +296,7 @@ struct IterationParameters { Scalar roadChi2MaxByNdf[(2 * maxChainLayers)-4]{}; /// The kinks a fit pass may take. std::uint8_t maxKinks = 0U; + FitReference fitReference = FitReference::kFilteredState; /// `IterationSettings::followEnd`, `innerTrim` /// and `crossingRule`. bool followEnd = false; @@ -568,6 +571,7 @@ template parameters.roadChi2MaxByNdf[ndf] = static_cast(chiSquareQuantile(ndf, settings.windows.nSigma)); } parameters.maxKinks = static_cast(std::max(settings.maxKinks, 0)); + parameters.fitReference = settings.fitReference; parameters.innerTrim = static_cast(std::clamp(settings.innerTrim, 0, static_cast(maxChainLayers))); parameters.crossingRule = settings.crossingRule; parameters.kinkChi2Min = static_cast(chiSquareQuantile(2, settings.windows.nSigma)); diff --git a/Detectors/GenericTracking/tracking/include/GenericTracking/Tracking/RoadFit.h b/Detectors/GenericTracking/tracking/include/GenericTracking/Tracking/RoadFit.h index fda09b03fb760..33504001cef66 100644 --- a/Detectors/GenericTracking/tracking/include/GenericTracking/Tracking/RoadFit.h +++ b/Detectors/GenericTracking/tracking/include/GenericTracking/Tracking/RoadFit.h @@ -133,7 +133,7 @@ template template [[nodiscard]] GPUdi() bool clearsBeamLine(const Propagator& propagator, const IterationParameters& parameters, const TimeFrame& frame, const TrackingGeometry& geometry, const Material& material, const std::uint32_t* rows, std::uint32_t count, const BoundTrackState& outwardState, Scalar& chi2) { - const auto inward = fitInward(propagator, geometry, material, RowChain{.frame = &frame, .rows = rows}, static_cast(count), outwardState, kinkBudget(parameters)); + const auto inward = fitInwardReference(propagator, geometry, material, RowChain{.frame = &frame, .rows = rows}, static_cast(count), outwardState, parameters.fitReference, kinkBudget(parameters)); if (inward.status != PropagationStatus::kOk) { return false; } diff --git a/Detectors/GenericTracking/tracking/include/GenericTracking/Tracking/TrackFitter.h b/Detectors/GenericTracking/tracking/include/GenericTracking/Tracking/TrackFitter.h index d9e5446e389d5..a1c6c4d6ceffb 100644 --- a/Detectors/GenericTracking/tracking/include/GenericTracking/Tracking/TrackFitter.h +++ b/Detectors/GenericTracking/tracking/include/GenericTracking/Tracking/TrackFitter.h @@ -10,12 +10,13 @@ // or submit itself to any jurisdiction. /// \file TrackFitter.h -/// \brief Independent outward and inward Kalman filters over a chain, each step linearised about its filtered state. +/// \brief Independent outward and inward Kalman filters over a chain. #ifndef GTRK_TRACKING_TRACKFITTER_H_ #define GTRK_TRACKING_TRACKFITTER_H_ #include "GenericTracking/Base/LayerSpec.h" +#include "GenericTracking/Base/LayerLinks.h" #include "GenericTracking/Tracking/CovarianceTransport.h" #include "GenericTracking/Tracking/Definitions.h" #include "GenericTracking/Tracking/Material.h" @@ -44,6 +45,18 @@ struct KinkBudget { Scalar gate = 0.F; ///< `chiSquareQuantile(2, nSigma)` }; +/// Linearisation of the road fits and final/follower refits. +/// The outward reference starts from the whole-chain circle seed. In +/// bidirectional mode the inward reference starts from the outward fitted +/// parameters. Each reference is transported independently and position-anchored +/// at every measurement, without Kalman updates to its direction or momentum. +/// Neither mode changes the reset covariance. +enum class FitReference : std::uint8_t { + kFilteredState = 0, + kCircleOutward = 1, + kIndependentBoth = 2 +}; + /// An iteration's kink budget, from its parameters' `maxKinks` and /// `kinkChi2Min`: the refit's and the extension's. template @@ -183,54 +196,84 @@ GPUdi() PropagationStatus transportStep(const Prop& prop, const StepMaterial& st return PropagationStatus::kOk; } -/// One step of a filter pass: \a state carried to \a measurement's plane -/// and updated with it, linearised about itself -- the filtered state the step -/// starts from, as ACTS's Kalman filter (a fixed whole-chain reference, never -/// corrected by the measurements, falls short of the outer planes of a low-momentum track near -/// its turning point). The covariance is \a factors, the pass's -/// (CovarianceFactors): `state.cov` is neither read nor written. \a Scaled, -/// the step's scattering theta0 is multiplied by \a noiseScale, both halves: -/// the kink allowance's last step (`allowanceStep`); the nominal step is its -/// own instantiation. \a Rule is the propagation's choice of crossing -/// (`CrossingRule`), and with it every pass below. Never inlined on the -/// device: one copy for each instantiation, not one in every fit of every -/// stage variant. -template -GPUdnii() PropagationStatus filterStep(const Prop& prop, const TrackingGeometry& geometry, const Material& material, const Measurement& measurement, PropagationDirection direction, BoundTrackState& state, CovarianceFactors& factors, Scalar& chi2, [[maybe_unused]] Scalar noiseScale = 1.F) +/// One filter step, with the pass covariance held in factors. Normally the +/// filtered state supplies the linearisation. With Reference, F and material +/// are evaluated on an independent trajectory and the mean is predicted as +/// f(reference) + F (state - reference). The propagated reference, including +/// mean energy loss, is optionally saved for the next step. After a successful +/// update its position is anchored to the measurement, as in native ITS; its +/// direction and momentum keep their propagated values. Fits accept an +/// intersection outside the sensor's active area. +/// Scaled multiplies both halves' scattering theta0 by noiseScale. +template +GPUdnii() PropagationStatus filterStep(const Prop& prop, const TrackingGeometry& geometry, const Material& material, const Measurement& measurement, PropagationDirection direction, BoundTrackState& state, CovarianceFactors& factors, Scalar& chi2, [[maybe_unused]] Scalar noiseScale = 1.F, [[maybe_unused]] const BoundTrackState* reference = nullptr, [[maybe_unused]] BoundTrackState* nextReference = nullptr) { - const auto& from = geometry.planeOf(state); + BoundTrackState prediction; + if constexpr (Reference) { + prediction = *reference; + } + auto& transportedState = Reference ? prediction : state; + const auto& from = geometry.planeOf(transportedState); const auto& to = geometry.planeOf(measurement); Scalar start[3]{}; - surfaceToGlobal(from, state.loc0(), state.loc1(), start); - const Scalar leaving = state.theta(); + surfaceToGlobal(from, transportedState.loc0(), transportedState.loc1(), start); + const Scalar leaving = transportedState.theta(); // An intersection just outside the sensor's active area is still the point // to update at: the cluster is there, so the track was. Bounds are for // searches, not fits (`jacobianAnalytic` takes kOutOfBounds). SolenoidJacobian jacobian; - const auto jacobianStatus = jacobianAnalytic(prop, state, from, to, direction, jacobian, &state); + const auto jacobianStatus = jacobianAnalytic(prop, transportedState, from, to, direction, jacobian, &transportedState); if (jacobianStatus != PropagationStatus::kOk) { return jacobianStatus; } Scalar end[3]{}; - surfaceToGlobal(to, state.loc0(), state.loc1(), end); + surfaceToGlobal(to, transportedState.loc0(), transportedState.loc1(), end); const auto step = stepMaterial(material, start, end); - const Scalar momentum = 1.F / o2::gpu::GPUCommonMath::Abs(state.qOverP()); + const Scalar momentum = 1.F / o2::gpu::GPUCommonMath::Abs(transportedState.qOverP()); const Scalar betaP = betaMomentum(momentum, prop.massGeV); // The nearer half before the transport (StepMaterial), at the direction on // the surface the step leaves; the farther half after it. const auto before = scatteringNoise(scaledTheta0(highlandTheta0(step.nearerX2X0, betaP), noiseScale), leaving); - const auto transported = transportStep(prop, step, betaP, before, direction, jacobian, state, factors, noiseScale); + const auto transported = transportStep(prop, step, betaP, before, direction, jacobian, transportedState, factors, noiseScale); if (transported != PropagationStatus::kOk) { return transported; } + if constexpr (Reference) { + if (nextReference != nullptr) { + *nextReference = prediction; + } + Scalar delta[kNBoundParams]; + for (int i = 0; i < kNBoundParams; ++i) { + delta[i] = state.par[i] - reference->par[i]; + } + o2::math_utils::detail::bringToPMPi(delta[kPhi]); + const Scalar* rows[3] = {jacobian.loc0, jacobian.loc1, jacobian.phi}; + for (int i = 0; i < 3; ++i) { + for (int j = 0; j < kNBoundParams; ++j) { + prediction.par[i] += rows[i][j] * delta[j]; + } + } + // These rows are the identity, as in the existing energy-loss transport. + prediction.par[kTheta] += delta[kTheta]; + prediction.par[kQOverP] += delta[kQOverP]; + state = prediction; + } const auto update = updateWithMeasurement(state, factors, measurement); if (!updateUsable(update, factors)) { return PropagationStatus::kNotFinite; } + if constexpr (Reference) { + if (nextReference != nullptr) { + // Re-anchor the reference to the observed position, without feeding the + // filtered direction, momentum or covariance into it (ITS ShiftRefToCluster). + nextReference->par[kLoc0] = measurement.loc0; + nextReference->par[kLoc1] = measurement.loc1; + } + } chi2 += update.chi2; return PropagationStatus::kOk; } @@ -263,16 +306,16 @@ GPUhdi() GPUAlwaysInline() void freeDirection(BoundTrackState& state, Covariance /// \a followed -- a measurement comes after this one in the pass -- and its /// increment is beyond the gate, taken again from the surface before with the /// direction freed, and counted in \a kinks. -template -GPUdi() PropagationStatus kinkStep(const Prop& prop, const TrackingGeometry& geometry, const Material& material, const Measurement& measurement, PropagationDirection direction, BoundTrackState& state, CovarianceFactors& factors, Scalar& chi2, const KinkBudget& budget, bool followed, int& kinks) +template +GPUdi() PropagationStatus kinkStep(const Prop& prop, const TrackingGeometry& geometry, const Material& material, const Measurement& measurement, PropagationDirection direction, BoundTrackState& state, CovarianceFactors& factors, Scalar& chi2, const KinkBudget& budget, bool followed, int& kinks, const BoundTrackState* reference = nullptr, BoundTrackState* nextReference = nullptr) { if (kinks >= budget.max || !followed) { - return filterStep(prop, geometry, material, measurement, direction, state, factors, chi2); + return filterStep(prop, geometry, material, measurement, direction, state, factors, chi2, 1.F, reference, nextReference); } const BoundTrackState keptState = state; const CovarianceFactors keptFactors = factors; const Scalar keptChi2 = chi2; - auto status = filterStep(prop, geometry, material, measurement, direction, state, factors, chi2); + auto status = filterStep(prop, geometry, material, measurement, direction, state, factors, chi2, 1.F, reference, nextReference); if (status != PropagationStatus::kOk || !(chi2 - keptChi2 > budget.gate)) { return status; } @@ -280,7 +323,7 @@ GPUdi() PropagationStatus kinkStep(const Prop& prop, const TrackingGeometry factors = keptFactors; chi2 = keptChi2; freeDirection(state, factors); - status = filterStep(prop, geometry, material, measurement, direction, state, factors, chi2); + status = filterStep(prop, geometry, material, measurement, direction, state, factors, chi2, 1.F, reference, nextReference); kinks += status == PropagationStatus::kOk ? 1 : 0; return status; } @@ -292,24 +335,24 @@ GPUdi() PropagationStatus kinkStep(const Prop& prop, const TrackingGeometry /// \a noiseScale, as the extension's kink step, and counted in \a kinks; /// \a refused if the increment is still beyond the gate. It frees no /// parameter. -template -GPUdi() PropagationStatus allowanceStep(const Propagator& prop, const TrackingGeometry& geometry, const Material& material, const Measurement& measurement, PropagationDirection direction, BoundTrackState& state, CovarianceFactors& factors, Scalar& chi2, const KinkBudget& budget, Scalar noiseScale, int& kinks, bool& refused) +template +GPUdi() PropagationStatus allowanceStep(const Propagator& prop, const TrackingGeometry& geometry, const Material& material, const Measurement& measurement, PropagationDirection direction, BoundTrackState& state, CovarianceFactors& factors, Scalar& chi2, const KinkBudget& budget, Scalar noiseScale, int& kinks, bool& refused, const BoundTrackState* reference = nullptr, BoundTrackState* nextReference = nullptr) { refused = false; if (kinks >= budget.max) { - return filterStep(prop, geometry, material, measurement, direction, state, factors, chi2); + return filterStep(prop, geometry, material, measurement, direction, state, factors, chi2, 1.F, reference, nextReference); } const BoundTrackState keptState = state; const CovarianceFactors keptFactors = factors; const Scalar keptChi2 = chi2; - auto status = filterStep(prop, geometry, material, measurement, direction, state, factors, chi2); + auto status = filterStep(prop, geometry, material, measurement, direction, state, factors, chi2, 1.F, reference, nextReference); if (status != PropagationStatus::kOk || !(chi2 - keptChi2 > budget.gate)) { return status; } state = keptState; factors = keptFactors; chi2 = keptChi2; - status = filterStep(prop, geometry, material, measurement, direction, state, factors, chi2, noiseScale); + status = filterStep(prop, geometry, material, measurement, direction, state, factors, chi2, noiseScale, reference, nextReference); if (status == PropagationStatus::kOk) { ++kinks; refused = !(chi2 - keptChi2 <= budget.gate); @@ -344,7 +387,7 @@ GPUhdi() GPUAlwaysInline() MeasurementUpdateResult startOnMeasurement(BoundTrack /// A filter pass from \a seed over the chain. Its covariance is carried as /// factors from the first update to the last (CovarianceFactors) and written /// to \a result's `cov` once, at the end. -template +template GPUdi() PropagationStatus runPass(const Prop& prop, const TrackingGeometry& geometry, const Material& material, const Chain& measurements, int n, PropagationDirection direction, BoundTrackState seed, BoundTrackState& result, Scalar& chi2, std::int16_t& failedMeasurement, const KinkBudget& passBudget, int& kinks, std::int16_t& kinkMeasurement, const Observer& observe) { const auto budget = kinkBudgetFor(passBudget, n); @@ -359,17 +402,27 @@ GPUdi() PropagationStatus runPass(const Prop& prop, const TrackingGeometry& return PropagationStatus::kNotFinite; } chi2 = firstUpdate.chi2; + BoundTrackState reference; + if constexpr (Reference) { + // Begin each pass from its own seed: circle outward, outward fit inward. + // startOnMeasurement has already anchored the seed's position. + reference = state; + } observe(direction, first, chi2); for (int step = 1; step < n; ++step) { const int index = first + step * increment; const Scalar before = chi2; const int kinksBefore = kinks; const auto& stepProp = stepPropagator(prop, geometry, measurements, index - increment, index); - const auto status = kinkStep(stepProp, geometry, material, measurements[index], direction, state, factors, chi2, budget, step < n - 1, kinks); + BoundTrackState nextReference; + const auto status = kinkStep(stepProp, geometry, material, measurements[index], direction, state, factors, chi2, budget, step < n - 1, kinks, &reference, &nextReference); if (status != PropagationStatus::kOk) { failedMeasurement = static_cast(index); return status; } + if constexpr (Reference) { + reference = nextReference; + } if (kinks != kinksBefore) { kinkMeasurement = static_cast(index); } @@ -382,7 +435,7 @@ GPUdi() PropagationStatus runPass(const Prop& prop, const TrackingGeometry& /// The inward pass of `fitMeasurements`: its seed, the outward pass's end /// state with a diffuse covariance, then the filter from there in. -template +template GPUdi() void inwardPass(const Prop& prop, const TrackingGeometry& geometry, const Material& material, const Chain& measurements, int n, Scalar qOverPSigma, const KinkBudget& kinkBudget, const Observer& observe, TrackFitResult& fit) { // From the outward pass's end state, its covariance made diffuse with the @@ -393,7 +446,7 @@ GPUdi() void inwardPass(const Prop& prop, const TrackingGeometry& geometry, setDiffuseCovariance(outerSeed, qOverPSigma); int kinks = 0; - const auto status = runPass(prop, geometry, material, measurements, n, PropagationDirection::kBackward, outerSeed, fit.state, fit.chi2, fit.failedMeasurement, kinkBudget, kinks, fit.kinkMeasurement, observe); + const auto status = runPass(prop, geometry, material, measurements, n, PropagationDirection::kBackward, outerSeed, fit.state, fit.chi2, fit.failedMeasurement, kinkBudget, kinks, fit.kinkMeasurement, observe); if (status != PropagationStatus::kOk) { fit.status = status; fit.failedPass = FitPass::kInward; @@ -441,6 +494,39 @@ GPUdi() TrackFitResult fitMeasurements(const Prop& prop, const TrackingGeometry< return fit; } +/// Refit using the whole-chain circle/longitudinal seed as an independent +/// linearisation trajectory outward. In bidirectional mode, an independent +/// reference starts inward from the outward fitted parameters. Both references +/// receive mean energy loss and position anchoring at each measurement; their +/// direction and momentum are not Kalman-updated within a pass. Outward-only +/// mode keeps the existing inward filter. No posterior covariance is recycled +/// as a prior. +/// Failed reference intersections are reported as failed fits, without fallback. +template +GPUdnii() TrackFitResult fitMeasurementsReference(const Prop& prop, const TrackingGeometry& geometry, const Material& material, const Chain& measurements, int n, FitReference referenceMode, const KinkBudget& kinkBudget = {}, const Observer& observe = {}) +{ + if (referenceMode == FitReference::kFilteredState) { + return fitMeasurements(prop, geometry, material, measurements, n, kinkBudget, observe); + } + TrackFitResult fit{.ndf = static_cast((2 * n) - kNBoundParams)}; + const auto seed = makeChainSeed(prop, geometry, measurements, n); + int outwardKinks = 0; + std::int16_t outwardKinkMeasurement = -1; + const auto status = detail::runPass(prop, geometry, material, measurements, n, PropagationDirection::kForward, seed.state, fit.outwardState, fit.outwardChi2, fit.failedMeasurement, kinkBudget, outwardKinks, outwardKinkMeasurement, observe); + fit.outwardKinks = static_cast(outwardKinks); + if (status != PropagationStatus::kOk) { + fit.status = status; + fit.failedPass = FitPass::kOutward; + return fit; + } + if (referenceMode == FitReference::kIndependentBoth) { + detail::inwardPass(prop, geometry, material, measurements, n, seed.qOverPSigma, kinkBudget, observe, fit); + } else { + detail::inwardPass(prop, geometry, material, measurements, n, seed.qOverPSigma, kinkBudget, observe, fit); + } + return fit; +} + /// `fitMeasurements` given its outward pass's end state, \a outwardState /// -- as the road fit keeps it for an accepted /// road (`RoadFitResult::state`): only the inward pass runs. The outward @@ -455,6 +541,22 @@ GPUdi() TrackFitResult fitInward(const Prop& prop, const TrackingGeometry& return fit; } +/// Finish a road fit without repeating its outward KF. Bidirectional mode +/// starts an independent inward reference from the outward fitted parameters; +/// outward-only and ordinary modes use the existing inward filter. +template +GPUdi() TrackFitResult fitInwardReference(const Prop& prop, const TrackingGeometry& geometry, const Material& material, const Chain& measurements, int n, const BoundTrackState& outwardState, FitReference mode, const KinkBudget& kinkBudget = {}, const Observer& observe = {}) +{ + if (mode != FitReference::kIndependentBoth) { + return fitInward(prop, geometry, material, measurements, n, outwardState, kinkBudget, observe); + } + TrackFitResult fit{.ndf = static_cast((2 * n) - kNBoundParams)}; + fit.outwardState = outwardState; + const auto seed = makeChainSeed(prop, geometry, measurements, n); + detail::inwardPass(prop, geometry, material, measurements, n, seed.qOverPSigma, kinkBudget, observe, fit); + return fit; +} + /// The road fit's candidate pass: the outward filter /// alone, as in `fitMeasurements`. struct RoadFitResult { @@ -480,7 +582,7 @@ struct RoadFitResult { /// beyond the gate, it frees no degree of freedom, and the result's kinks /// count it. A compile-time choice: without it the last step is a `kinkStep` /// like the others. -template +template GPUdi() RoadFitResult fitOutward(const Prop& prop, const TrackingGeometry& geometry, const Material& material, const Chain& measurements, int n, const Scalar* chi2MaxByNdf, const KinkBudget& kinkBudget = {}, [[maybe_unused]] Scalar lastStepScale = 1.F) { const auto budget = detail::kinkBudgetFor(kinkBudget, n); @@ -498,6 +600,10 @@ GPUdi() RoadFitResult fitOutward(const Prop& prop, const TrackingGeometry& return result; } result.chi2 = firstUpdate.chi2; + BoundTrackState reference; + if constexpr (Reference) { + reference = state; + } const int ndf = (2 * n) - kNBoundParams; int kinks = 0; // The allowance's last step, which frees no degree of freedom; counted by its @@ -515,11 +621,12 @@ GPUdi() RoadFitResult fitOutward(const Prop& prop, const TrackingGeometry& } countWork(WorkCounter::kRoadSteps); const auto& stepProp = stepPropagator(prop, geometry, measurements, index - 1, index); + BoundTrackState nextReference; if constexpr (LastStepAllowance) { if (index == n - 1) { int taken = kinks; bool refused = false; - const auto status = detail::allowanceStep(stepProp, geometry, material, measurements[index], PropagationDirection::kForward, state, factors, result.chi2, budget, lastStepScale, taken, refused); + const auto status = detail::allowanceStep(stepProp, geometry, material, measurements[index], PropagationDirection::kForward, state, factors, result.chi2, budget, lastStepScale, taken, refused, &reference, &nextReference); if (status != PropagationStatus::kOk) { countWork(WorkCounter::kRoadPropagation); result.status = RoadFitStatus::kPropagationFailed; @@ -531,15 +638,21 @@ GPUdi() RoadFitResult fitOutward(const Prop& prop, const TrackingGeometry& return result; } allowed = taken - kinks; + if constexpr (Reference) { + reference = nextReference; + } continue; } } - const auto status = detail::kinkStep(stepProp, geometry, material, measurements[index], PropagationDirection::kForward, state, factors, result.chi2, budget, index < n - 1, kinks); + const auto status = detail::kinkStep(stepProp, geometry, material, measurements[index], PropagationDirection::kForward, state, factors, result.chi2, budget, index < n - 1, kinks, &reference, &nextReference); if (status != PropagationStatus::kOk) { countWork(WorkCounter::kRoadPropagation); result.status = RoadFitStatus::kPropagationFailed; return result; } + if constexpr (Reference) { + reference = nextReference; + } } result.kinks = static_cast(kinks + allowed); covarianceFromFactors(factors, state.cov); @@ -559,6 +672,9 @@ GPUdi() RoadFitResult fitOutward(const Prop& prop, const TrackingGeometry& template [[nodiscard]] GPUdi() RoadFitResult fitRoad(const Parameters& parameters, const Prop& prop, const TrackingGeometry& geometry, const Material& material, const Chain& measurements, int n) { + if (parameters.fitReference != FitReference::kFilteredState) { + return fitOutward(prop, geometry, material, measurements, n, parameters.roadChi2MaxByNdf, kinkBudget(parameters), o2::gpu::CAMath::Sqrt(parameters.kinkScale2)); + } if constexpr (LastStepAllowance) { return fitOutward(prop, geometry, material, measurements, n, parameters.roadChi2MaxByNdf, kinkBudget(parameters), o2::gpu::CAMath::Sqrt(parameters.kinkScale2)); } else { diff --git a/Detectors/GenericTracking/tracking/include/GenericTracking/Tracking/TrackRefit.h b/Detectors/GenericTracking/tracking/include/GenericTracking/Tracking/TrackRefit.h index 32853986a3a14..98a7778d17b21 100644 --- a/Detectors/GenericTracking/tracking/include/GenericTracking/Tracking/TrackRefit.h +++ b/Detectors/GenericTracking/tracking/include/GenericTracking/Tracking/TrackRefit.h @@ -111,7 +111,8 @@ GPUhdi() GPUAlwaysInline() void writeFit(const Tracks& tracks, std::uint32_t tra /// One element per selected record: the inward pass on its rows from the road /// fit's outward state (`fitInward`), and the track's time from their boxes; -/// the pass may take the iteration's kinks. Its frame, geometry and material +/// the pass may take the iteration's kinks. Bidirectional reference mode starts +/// its inward reference from that outward state. Its frame, geometry and material /// are pointers, which may be a device's. template class TrackRefitElement @@ -120,7 +121,7 @@ class TrackRefitElement static constexpr std::size_t kMaxClusters = maxChainLayers; using View = typename BoundedColumns::template RecordT>::View; - TrackRefitElement(const TimeFrame* frame, const TrackingGeometry* geometry, const Material* material, const Propagator& propagator, const ClusterList* clusters, const std::uint8_t* count, const BoundTrackState* outward, const Scalar* outwardChi2, const std::uint32_t* selected, View tracks, std::uint32_t offset, std::uint8_t iteration, KinkBudget kinks) + TrackRefitElement(const TimeFrame* frame, const TrackingGeometry* geometry, const Material* material, const Propagator& propagator, const ClusterList* clusters, const std::uint8_t* count, const BoundTrackState* outward, const Scalar* outwardChi2, const std::uint32_t* selected, View tracks, std::uint32_t offset, std::uint8_t iteration, KinkBudget kinks, FitReference reference = FitReference::kFilteredState) : mFrame(frame), mGeometry(geometry), mMaterial(material), @@ -133,7 +134,8 @@ class TrackRefitElement mTracks(tracks), mOffset(offset), mIteration(iteration), - mKinks(kinks) + mKinks(kinks), + mReference(reference) { } @@ -144,7 +146,8 @@ class TrackRefitElement const auto& rows = mClusters[record]; const std::uint8_t count = mCount[record]; // The road fit ran the outward pass already and kept its end state. - const auto fit = fitInward(mPropagator, *mGeometry, *mMaterial, RowChain{.frame = mFrame, .rows = rows.index}, static_cast(count), mOutward[record], mKinks); + const auto chain = RowChain{.frame = mFrame, .rows = rows.index}; + const auto fit = fitInwardReference(mPropagator, *mGeometry, *mMaterial, chain, static_cast(count), mOutward[record], mReference, mKinks); writeFit(mTracks, track, fit, rows, count, trackTimeBox(mFrame->measurements(), rows, count)); mTracks.outwardChi2[track] = mOutwardChi2[record]; mTracks.iteration[track] = mIteration; @@ -166,6 +169,7 @@ class TrackRefitElement std::uint32_t mOffset = 0U; std::uint8_t mIteration = 0U; KinkBudget mKinks; + FitReference mReference; }; /// The timeframe's tracks: each iteration's appended after the ones before, @@ -200,7 +204,7 @@ class TrackRefit const auto offset = static_cast(mTracks.size()); mTracks.resize(offset + selected.size()); const auto records = selection.records(); - const TrackRefitElement element{&frame, &geometry, &material, propagator, records.clusters, records.count, records.outwardState, records.chi2, selected.data(), mTracks.view(), offset, iteration, kinkBudget(parameters)}; + const TrackRefitElement element{&frame, &geometry, &material, propagator, records.clusters, records.count, records.outwardState, records.chi2, selected.data(), mTracks.view(), offset, iteration, kinkBudget(parameters), parameters.fitReference}; forEach(static_cast(selected.size()), element); return StageStatus::kComplete; }); diff --git a/Detectors/GenericTracking/tracking/test/device/testDeviceTracker.cxx b/Detectors/GenericTracking/tracking/test/device/testDeviceTracker.cxx index e65c5913e89c2..a2944f3788c9f 100644 --- a/Detectors/GenericTracking/tracking/test/device/testDeviceTracker.cxx +++ b/Detectors/GenericTracking/tracking/test/device/testDeviceTracker.cxx @@ -172,6 +172,7 @@ struct Options { /// never both in one iteration. double kinkSigma = 1.; bool followEnd = false; + FitReference reference = FitReference::kFilteredState; }; std::string describe(const Options& options) { return std::string{options.source == VertexSource::kBeamLine ? "beam line" : "vertex list"} + ", width " + std::to_string(options.width) + ", shared " + std::to_string(options.sharedInnermost) + (options.extension ? ", follow" : ", no follow") + (options.kinkSigma > 1. ? ", kink allowance " + std::to_string(options.kinkSigma) + " sigma" : "") + (options.followEnd ? ", end again" : ""); } @@ -182,6 +183,7 @@ std::unique_ptr makeHost(const Options& options) static const UniformMaterial none{}; auto settings = iterations(options.source); for (auto& iteration : settings) { + iteration.fitReference = options.reference; iteration.sharedInnermost = options.sharedInnermost; iteration.follow = options.extension; } @@ -232,7 +234,44 @@ bool sameTracks(const View& host, std::size_t hostRows, const View& device, std: std::ranges::sort(result); return result; }; - return hostRows == deviceRows && rowsOf(host, hostRows) == rowsOf(device, deviceRows); + if (hostRows != deviceRows || rowsOf(host, hostRows) != rowsOf(device, deviceRows)) { + return false; + } + // Check states explicitly: a non-arithmetic memlayout column has no scalar + // maxAbsDifference in StageComparison, so matching identities is not a + // numerical check of its parameters and covariance. + const auto order = [](const View& view, std::size_t rows) { + std::vector, std::size_t>> sorted; + for (std::size_t i = 0; i < rows; ++i) { + sorted.emplace_back(std::vector{view.clusters[i].index, view.clusters[i].index + view.count[i]}, i); + } + std::ranges::sort(sorted); + return sorted; + }; + const auto lhs = order(host, hostRows); + const auto rhs = order(device, deviceRows); + const auto closeState = [](const BoundTrackState& a, const BoundTrackState& b) { + constexpr Scalar tolerance = 1.e-4F; + if (a.surface != b.surface || std::abs(a.surfacePhi - b.surfacePhi) > tolerance) { + return false; + } + for (int i = 0; i < kNBoundParams; ++i) { + Scalar delta = a.par[i] - b.par[i]; + if (i == kPhi) { + o2::math_utils::detail::bringToPMPi(delta); + } + if (!(std::abs(delta) <= tolerance)) { + return false; + } + } + return std::ranges::equal(a.cov, b.cov, [=](Scalar x, Scalar y) { return std::abs(x - y) <= tolerance; }); + }; + for (std::size_t i = 0; i < lhs.size(); ++i) { + if (!closeState(host.state[lhs[i].second], device.state[rhs[i].second]) || !closeState(host.outwardState[lhs[i].second], device.outwardState[rhs[i].second])) { + return false; + } + } + return true; #endif } @@ -397,22 +436,25 @@ BOOST_AUTO_TEST_CASE(EverySliceAndIterationMatchesTheHostTracker) auto* device = runtime(); const auto tracker = kernels().makeToyBarrel8Tracker(*device, block(kArenaBytes)); BOOST_TEST(tracker->arenaBytes() == kArenaBytes); - for (const auto& options : configurations()) { - BOOST_TEST_CONTEXT(describe(options)) - { - const auto totals = compare(*tracker, options); - BOOST_TEST(totals.triplets > 0U); - BOOST_TEST(totals.neighbours > 0U); - BOOST_TEST(totals.chained > 0U); - BOOST_TEST(totals.roads > 0U); - BOOST_TEST(totals.records > 0U); - BOOST_TEST(totals.gathered == totals.records); - BOOST_TEST(totals.tracks == kTracks); - BOOST_TEST(totals.maskedSlices > 0U); - if (options.kinkSigma > 1.) { - BOOST_TEST(totals.extended > 0U); + for (const auto reference : {FitReference::kFilteredState, FitReference::kCircleOutward, FitReference::kIndependentBoth}) { + for (auto options : configurations()) { + options.reference = reference; + BOOST_TEST_CONTEXT(describe(options) << ", reference " << static_cast(reference)) + { + const auto totals = compare(*tracker, options); + BOOST_TEST(totals.triplets > 0U); + BOOST_TEST(totals.neighbours > 0U); + BOOST_TEST(totals.chained > 0U); + BOOST_TEST(totals.roads > 0U); + BOOST_TEST(totals.records > 0U); + BOOST_TEST(totals.gathered == totals.records); + BOOST_TEST(totals.tracks == kTracks); + BOOST_TEST(totals.maskedSlices > 0U); + if (options.kinkSigma > 1.) { + BOOST_TEST(totals.extended > 0U); + } + BOOST_TEST_MESSAGE(device->backend() << " " << describe(options) << ": " << totals.slices << " slices, " << totals.tracklets << " tracklets, " << totals.triplets << " triplets, " << totals.neighbours << " neighbours, " << totals.roads << " roads, " << totals.records << " records, " << totals.rounds << " rounds, " << totals.tracks << " tracks, " << totals.extended << " extended"); } - BOOST_TEST_MESSAGE(device->backend() << " " << describe(options) << ": " << totals.slices << " slices, " << totals.tracklets << " tracklets, " << totals.triplets << " triplets, " << totals.neighbours << " neighbours, " << totals.roads << " roads, " << totals.records << " records, " << totals.rounds << " rounds, " << totals.tracks << " tracks, " << totals.extended << " extended"); } } } diff --git a/Detectors/GenericTracking/tracking/test/testRoadFit.cxx b/Detectors/GenericTracking/tracking/test/testRoadFit.cxx index 6631e5776d53b..afe14afb233cd 100644 --- a/Detectors/GenericTracking/tracking/test/testRoadFit.cxx +++ b/Detectors/GenericTracking/tracking/test/testRoadFit.cxx @@ -220,23 +220,28 @@ BOOST_AUTO_TEST_CASE(SmearedRoadsHaveTheirDegreesOfFreedom) hit.point.z = global[2]; } const Fixture fixture{hits, 100.}; - const auto cut = study::roadChi2Cut(parameters(), kLayers); - double sum = 0.; - std::size_t accepted = 0U; - for (std::size_t track = 0; track < kTracks; ++track) { - const auto rows = trackRows(fixture, track); - const RowChain chain{.frame = fixture.frame.get(), .rows = rows.data()}; - const auto fit = fitOutward(kPropagator, trackingGeometry(), kNoMaterial, chain, static_cast(kLayers), kNoCut.data()); - BOOST_REQUIRE((fit.status == RoadFitStatus::kAccepted)); - sum += static_cast(fit.chi2); - accepted += fit.chi2 <= cut ? 1U : 0U; + for (const auto mode : {FitReference::kFilteredState, FitReference::kCircleOutward, FitReference::kIndependentBoth}) { + const auto cut = study::roadChi2Cut(parameters(), kLayers); + auto iteration = parameters(); + iteration.fitReference = mode; + std::fill(std::begin(iteration.roadChi2MaxByNdf), std::end(iteration.roadChi2MaxByNdf), std::numeric_limits::infinity()); + double sum = 0.; + std::size_t accepted = 0U; + for (std::size_t track = 0; track < kTracks; ++track) { + const auto rows = trackRows(fixture, track); + const RowChain chain{.frame = fixture.frame.get(), .rows = rows.data()}; + const auto fit = fitRoad(iteration, kPropagator, trackingGeometry(), kNoMaterial, chain, static_cast(kLayers)); + BOOST_REQUIRE((fit.status == RoadFitStatus::kAccepted)); + sum += static_cast(fit.chi2); + accepted += fit.chi2 <= cut ? 1U : 0U; + } + const double mean = sum / static_cast(kTracks); + BOOST_TEST_MESSAGE("mean chi2 " << mean << " (11), within the cut " << accepted << "/" << kTracks); + // The mean's error over 400 fits is sqrt(2 * 11 / 400) = 0.23. + BOOST_TEST(mean > 10.F); + BOOST_TEST(mean < 12.F); + BOOST_TEST(accepted >= 0.99 * kTracks); } - const double mean = sum / static_cast(kTracks); - BOOST_TEST_MESSAGE("mean chi2 " << mean << " (11), within the cut " << accepted << "/" << kTracks); - // The mean's error over 400 fits is sqrt(2 * 11 / 400) = 0.23. - BOOST_TEST(mean > 10.F); - BOOST_TEST(mean < 12.F); - BOOST_TEST(accepted >= 0.99 * kTracks); } /// A helix whose direction turns by \a turn on the sensor of layer \a kink -- @@ -572,78 +577,81 @@ BOOST_AUTO_TEST_CASE(SelectedTracksAreRefittedWithBothPasses) // fit, the selection, the refit. constexpr std::size_t kTracks = 6U; const Fixture fixture{trackHits(kTracks)}; - const auto iteration = parameters(); - const RoadsRun roads{fixture}; - const auto& run = roads.result; - const auto& count = roads.count; - const SerialForEach serial{}; - const Stage stage{iteration, *fixture.frame, trackingGeometry(), kNoMaterial, kPropagator, *run.tracklets, *run.triplets, count.walk(*run.triplets, *run.neighbours)}; - RecordOutput records; - records.setMemoryPool(std::make_shared()); - records.prepare(count.roadCount()); - BOOST_REQUIRE(runSerial(records, stage) == StageStatus::kComplete); - GlobalSelection selection; - selection.setMemoryPool(std::make_shared()); - BOOST_REQUIRE(selection.append(records, SerialForEach{}) == StageStatus::kComplete); - std::vector mask(fixture.frame->size(), kFree); - BOOST_REQUIRE(selection.select(serial, [](std::span indices, const auto& less) { std::ranges::sort(indices, less); }, 0U, mask) == StageStatus::kComplete); - BOOST_REQUIRE(selection.selected().size() == kTracks); - - TrackRefit refit; - refit.setMemoryPool(std::make_shared()); - BOOST_REQUIRE(refit.run(serial, selection, iteration, *fixture.frame, trackingGeometry(), kNoMaterial, kPropagator, 0U) == StageStatus::kComplete); - BOOST_REQUIRE(refit.size() == kTracks); - const auto tracks = refit.tracks().view(); - for (std::uint32_t track = 0; track < kTracks; ++track) { - const auto& rows = tracks.clusters[track]; - const auto index = fixture.frame->measurement(rows.index[0]).clusterIndex / kLayers; - const auto h = helix(index); - BOOST_TEST_CONTEXT("track " << index) - { - // Both passes, as fitMeasurements gives them on the same rows. - BOOST_REQUIRE(tracks.status[track] == PropagationStatus::kOk); - BOOST_TEST(tracks.count[track] == kLayers); - std::array copied{}; - const RowChain chain{.frame = fixture.frame.get(), .rows = rows.index}; - for (std::size_t i = 0; i < kLayers; ++i) { - copied[i] = chain[static_cast(i)]; + for (const auto mode : {FitReference::kFilteredState, FitReference::kCircleOutward, FitReference::kIndependentBoth}) { + auto iteration = parameters(); + iteration.fitReference = mode; + const RoadsRun roads{fixture}; + const auto& run = roads.result; + const auto& count = roads.count; + const SerialForEach serial{}; + const Stage stage{iteration, *fixture.frame, trackingGeometry(), kNoMaterial, kPropagator, *run.tracklets, *run.triplets, count.walk(*run.triplets, *run.neighbours)}; + RecordOutput records; + records.setMemoryPool(std::make_shared()); + records.prepare(count.roadCount()); + BOOST_REQUIRE(runSerial(records, stage) == StageStatus::kComplete); + GlobalSelection selection; + selection.setMemoryPool(std::make_shared()); + BOOST_REQUIRE(selection.append(records, SerialForEach{}) == StageStatus::kComplete); + std::vector mask(fixture.frame->size(), kFree); + BOOST_REQUIRE(selection.select(serial, [](std::span indices, const auto& less) { std::ranges::sort(indices, less); }, 0U, mask) == StageStatus::kComplete); + BOOST_REQUIRE(selection.selected().size() == kTracks); + + TrackRefit refit; + refit.setMemoryPool(std::make_shared()); + BOOST_REQUIRE(refit.run(serial, selection, iteration, *fixture.frame, trackingGeometry(), kNoMaterial, kPropagator, 0U) == StageStatus::kComplete); + BOOST_REQUIRE(refit.size() == kTracks); + const auto tracks = refit.tracks().view(); + for (std::uint32_t track = 0; track < kTracks; ++track) { + const auto& rows = tracks.clusters[track]; + const auto index = fixture.frame->measurement(rows.index[0]).clusterIndex / kLayers; + const auto h = helix(index); + BOOST_TEST_CONTEXT("track " << index) + { + // Both passes, as fitMeasurements gives them on the same rows. + BOOST_REQUIRE(tracks.status[track] == PropagationStatus::kOk); + BOOST_TEST(tracks.count[track] == kLayers); + std::array copied{}; + const RowChain chain{.frame = fixture.frame.get(), .rows = rows.index}; + for (std::size_t i = 0; i < kLayers; ++i) { + copied[i] = chain[static_cast(i)]; + } + const auto full = fitMeasurementsReference(kPropagator, trackingGeometry(), kNoMaterial, copied.data(), static_cast(kLayers), mode); + BOOST_TEST(tracks.chi2[track] == full.chi2); + + BOOST_TEST(tracks.ndf[track] == (2 * static_cast(kLayers)) - 5); + // On the innermost surface, and the generating helix. + const auto& state = tracks.state[track]; + const auto first = fixture.frame->measurement(rows.index[0]); + BOOST_TEST(state.surface == first.surface); + BOOST_TEST(std::abs(state.loc0() - first.loc0) < 1.e-4F); + BOOST_TEST(std::abs(state.loc1() - first.loc1) < 1.e-4F); + BOOST_TEST(state.theta() == static_cast(std::atan2(1., h.slope)), boost::test_tools::tolerance(1.e-4F)); + const double momentum = h.ptGeV * std::hypot(1., h.slope); + BOOST_TEST(std::abs(state.qOverP()) == static_cast(1. / momentum), boost::test_tools::tolerance(1.e-3F)); } - const auto full = fitMeasurements(kPropagator, trackingGeometry(), kNoMaterial, copied.data(), static_cast(kLayers)); - BOOST_TEST(tracks.chi2[track] == full.chi2); - - BOOST_TEST(tracks.ndf[track] == (2 * static_cast(kLayers)) - 5); - // On the innermost surface, and the generating helix. - const auto& state = tracks.state[track]; - const auto first = fixture.frame->measurement(rows.index[0]); - BOOST_TEST(state.surface == first.surface); - BOOST_TEST(std::abs(state.loc0() - first.loc0) < 1.e-4F); - BOOST_TEST(std::abs(state.loc1() - first.loc1) < 1.e-4F); - BOOST_TEST(state.theta() == static_cast(std::atan2(1., h.slope)), boost::test_tools::tolerance(1.e-4F)); - const double momentum = h.ptGeV * std::hypot(1., h.slope); - BOOST_TEST(std::abs(state.qOverP()) == static_cast(1. / momentum), boost::test_tools::tolerance(1.e-3F)); } - } - for (const int threads : {1, 2, 4}) { - detail::TaskArena cpuArena{threads}; - CpuStage::template RecordT> cpu{cpuArena}; - cpu.setMemoryPool(std::make_shared()); - TrackRefit parallel; - parallel.setMemoryPool(std::make_shared()); - BOOST_REQUIRE(parallel.run([&](std::uint32_t elements, const auto& body) { cpu.forEachElement(elements, 1U, body); }, selection, iteration, *fixture.frame, trackingGeometry(), kNoMaterial, kPropagator, 0U) == StageStatus::kComplete); - // Bit for bit. - const auto other = parallel.tracks().view(); - const auto bits = [](Scalar value) { return std::bit_cast(value); }; - for (std::uint32_t track = 0; track < kTracks; ++track) { - const auto& a = other.state[track]; - const auto& b = tracks.state[track]; - BOOST_TEST(std::ranges::equal(a.par, b.par, {}, bits, bits)); - BOOST_TEST(std::ranges::equal(a.cov, b.cov, {}, bits, bits)); - BOOST_TEST(other.timePS[track] == tracks.timePS[track]); - BOOST_TEST(other.widthPS[track] == tracks.widthPS[track]); - BOOST_TEST(a.surface == b.surface); - BOOST_TEST(bits(other.chi2[track]) == bits(tracks.chi2[track])); - BOOST_TEST(std::ranges::equal(other.clusters[track].index, tracks.clusters[track].index)); + for (const int threads : {1, 2, 4}) { + detail::TaskArena cpuArena{threads}; + CpuStage::template RecordT> cpu{cpuArena}; + cpu.setMemoryPool(std::make_shared()); + TrackRefit parallel; + parallel.setMemoryPool(std::make_shared()); + BOOST_REQUIRE(parallel.run([&](std::uint32_t elements, const auto& body) { cpu.forEachElement(elements, 1U, body); }, selection, iteration, *fixture.frame, trackingGeometry(), kNoMaterial, kPropagator, 0U) == StageStatus::kComplete); + // Bit for bit. + const auto other = parallel.tracks().view(); + const auto bits = [](Scalar value) { return std::bit_cast(value); }; + for (std::uint32_t track = 0; track < kTracks; ++track) { + const auto& a = other.state[track]; + const auto& b = tracks.state[track]; + BOOST_TEST(std::ranges::equal(a.par, b.par, {}, bits, bits)); + BOOST_TEST(std::ranges::equal(a.cov, b.cov, {}, bits, bits)); + BOOST_TEST(other.timePS[track] == tracks.timePS[track]); + BOOST_TEST(other.widthPS[track] == tracks.widthPS[track]); + BOOST_TEST(a.surface == b.surface); + BOOST_TEST(bits(other.chi2[track]) == bits(tracks.chi2[track])); + BOOST_TEST(std::ranges::equal(other.clusters[track].index, tracks.clusters[track].index)); + } } } } diff --git a/Detectors/GenericTracking/tracking/test/testTrackFitter.cxx b/Detectors/GenericTracking/tracking/test/testTrackFitter.cxx index 1d49207d7c256..3569abfaf4c50 100644 --- a/Detectors/GenericTracking/tracking/test/testTrackFitter.cxx +++ b/Detectors/GenericTracking/tracking/test/testTrackFitter.cxx @@ -216,19 +216,24 @@ BOOST_AUTO_TEST_CASE(ThreePointSeedRemainsFiniteAsCurvatureDiverges) BOOST_AUTO_TEST_CASE(IndependentPassesFitTheSameNoiselessChain) { - const SyntheticChain chain; - const UniformMaterial material; - const auto fit = fitMeasurements(chain.propagator, chain.tracking, material, chain.measurements.data(), static_cast(chain.measurements.size())); - - BOOST_TEST(static_cast(fit.status) == static_cast(PropagationStatus::kOk)); - BOOST_TEST(fit.ndf == 1); - BOOST_TEST(fit.outwardState.surface == chain.truth[2].surface); - BOOST_TEST(fit.state.surface == chain.truth[0].surface); - BOOST_TEST(std::abs(fit.state.phi() - chain.truth[0].phi()) < 1e-5F); - BOOST_TEST(std::abs(fit.state.theta() - chain.truth[0].theta()) < 1e-5F); - BOOST_TEST(std::abs(fit.state.qOverP() - chain.truth[0].qOverP()) < 1e-3F); - BOOST_TEST(std::abs(fit.chi2) < 1e-4F); - BOOST_TEST(std::abs(fit.outwardChi2) < 1e-4F); + for (const auto referenceMode : {FitReference::kFilteredState, FitReference::kCircleOutward, FitReference::kIndependentBoth}) { + BOOST_TEST_CONTEXT("reference mode " << static_cast(referenceMode)) + { + const SyntheticChain chain; + const UniformMaterial material; + const auto fit = fitMeasurementsReference(chain.propagator, chain.tracking, material, chain.measurements.data(), static_cast(chain.measurements.size()), referenceMode); + + BOOST_TEST(static_cast(fit.status) == static_cast(PropagationStatus::kOk)); + BOOST_TEST(fit.ndf == 1); + BOOST_TEST(fit.outwardState.surface == chain.truth[2].surface); + BOOST_TEST(fit.state.surface == chain.truth[0].surface); + BOOST_TEST(std::abs(fit.state.phi() - chain.truth[0].phi()) < 1e-5F); + BOOST_TEST(std::abs(fit.state.theta() - chain.truth[0].theta()) < 1e-5F); + BOOST_TEST(std::abs(fit.state.qOverP() - chain.truth[0].qOverP()) < 1e-3F); + BOOST_TEST(std::abs(fit.chi2) < 1e-4F); + BOOST_TEST(std::abs(fit.outwardChi2) < 1e-4F); + } + } } /// A primary from a vertex in the diamond's tail crosses the far end of barrel @@ -239,105 +244,115 @@ BOOST_AUTO_TEST_CASE(IndependentPassesFitTheSameNoiselessChain) /// last. BOOST_AUTO_TEST_CASE(LayerZeroAtShallowIncidenceFitsBothWays) { - using Mixed = detectors::ToyBarrel4Disc8; - const auto detector = toyGeometry(); - const TrackingGeometry tracking{detector}; - const Propagator propagator; - const UniformMaterial material; - - SensorID start = firstSensorOfLayer(0); - for (SensorID sensor = start; sensor < firstSensorOfLayer(1); ++sensor) { - if (tracking.surface(sensor).origin[2] > tracking.surface(start).origin[2]) { - start = sensor; - } - } - const auto& first = tracking.surface(start); - std::array truth; - truth[0].par[kPhi] = std::atan2(first.origin[1], first.origin[0]); - truth[0].par[kTheta] = 0.048F; - truth[0].par[kQOverP] = 0.05F; - truth[0].surface = start; - const double cosIncidence = std::sin(0.048) * (std::cos(static_cast(truth[0].phi())) * static_cast(first.normal[0]) + std::sin(static_cast(truth[0].phi())) * static_cast(first.normal[1])); - BOOST_TEST_REQUIRE(std::abs(cosIncidence) < 0.05); - BOOST_TEST_REQUIRE(std::abs(cosIncidence) > static_cast(propagator.minAbsCosIncidence)); - - constexpr std::array kDiscs{8, 9, 10, 11}; - for (std::size_t i = 0; i < kDiscs.size(); ++i) { - const auto& from = tracking.surface(truth[i].surface); - for (SensorID sensor = firstSensorOfLayer(kDiscs[i]); sensor < firstSensorOfLayer(static_cast(kDiscs[i] + 1U)); ++sensor) { - BoundTrackState probe = truth[i]; - if (propagate(propagator, probe, from, tracking.surface(sensor), PropagationDirection::kForward).status == PropagationStatus::kOk) { - truth[i + 1U] = probe; - break; + for (const auto referenceMode : {FitReference::kFilteredState, FitReference::kCircleOutward, FitReference::kIndependentBoth}) { + BOOST_TEST_CONTEXT("reference mode " << static_cast(referenceMode)) + { + using Mixed = detectors::ToyBarrel4Disc8; + const auto detector = toyGeometry(); + const TrackingGeometry tracking{detector}; + const Propagator propagator; + const UniformMaterial material; + + SensorID start = firstSensorOfLayer(0); + for (SensorID sensor = start; sensor < firstSensorOfLayer(1); ++sensor) { + if (tracking.surface(sensor).origin[2] > tracking.surface(start).origin[2]) { + start = sensor; + } + } + const auto& first = tracking.surface(start); + std::array truth; + truth[0].par[kPhi] = std::atan2(first.origin[1], first.origin[0]); + truth[0].par[kTheta] = 0.048F; + truth[0].par[kQOverP] = 0.05F; + truth[0].surface = start; + const double cosIncidence = std::sin(0.048) * (std::cos(static_cast(truth[0].phi())) * static_cast(first.normal[0]) + std::sin(static_cast(truth[0].phi())) * static_cast(first.normal[1])); + BOOST_TEST_REQUIRE(std::abs(cosIncidence) < 0.05); + BOOST_TEST_REQUIRE(std::abs(cosIncidence) > static_cast(propagator.minAbsCosIncidence)); + + constexpr std::array kDiscs{8, 9, 10, 11}; + for (std::size_t i = 0; i < kDiscs.size(); ++i) { + const auto& from = tracking.surface(truth[i].surface); + for (SensorID sensor = firstSensorOfLayer(kDiscs[i]); sensor < firstSensorOfLayer(static_cast(kDiscs[i] + 1U)); ++sensor) { + BoundTrackState probe = truth[i]; + if (propagate(propagator, probe, from, tracking.surface(sensor), PropagationDirection::kForward).status == PropagationStatus::kOk) { + truth[i + 1U] = probe; + break; + } + } + BOOST_TEST_REQUIRE(truth[i + 1U].surface != truth[i].surface); + } + std::array measurements; + for (std::size_t i = 0; i < measurements.size(); ++i) { + measurements[i] = measurementFromState(truth[i], static_cast(i)); } - } - BOOST_TEST_REQUIRE(truth[i + 1U].surface != truth[i].surface); - } - std::array measurements; - for (std::size_t i = 0; i < measurements.size(); ++i) { - measurements[i] = measurementFromState(truth[i], static_cast(i)); - } - const auto fit = fitMeasurements(propagator, tracking, material, measurements.data(), static_cast(measurements.size())); - BOOST_TEST(static_cast(fit.status) == static_cast(PropagationStatus::kOk)); - BOOST_TEST(fit.state.surface == start); - BOOST_TEST(std::abs(fit.state.theta() - truth[0].theta()) < 1e-5F); + const auto fit = fitMeasurementsReference(propagator, tracking, material, measurements.data(), static_cast(measurements.size()), referenceMode); + BOOST_TEST(static_cast(fit.status) == static_cast(PropagationStatus::kOk)); + BOOST_TEST(fit.state.surface == start); + BOOST_TEST(std::abs(fit.state.theta() - truth[0].theta()) < 1e-5F); - const Propagator former{.minAbsCosIncidence = 0.05F}; - const auto failed = fitMeasurements(former, tracking, material, measurements.data(), static_cast(measurements.size())); - BOOST_TEST(static_cast(failed.status) == static_cast(PropagationStatus::kGrazing)); - BOOST_TEST(static_cast(failed.failedPass) == static_cast(FitPass::kInward)); - BOOST_TEST(failed.failedMeasurement == 0); + const Propagator former{.minAbsCosIncidence = 0.05F}; + const auto failed = fitMeasurementsReference(former, tracking, material, measurements.data(), static_cast(measurements.size()), referenceMode); + BOOST_TEST(static_cast(failed.status) == static_cast(PropagationStatus::kGrazing)); + BOOST_TEST(static_cast(failed.failedPass) == static_cast(FitPass::kInward)); + BOOST_TEST(failed.failedMeasurement == 0); + } + } } BOOST_AUTO_TEST_CASE(SmearedChainsGiveUnbiasedParametersAndUnitPulls) { - // A fit can be self-consistent and still wrong in - // exactly one way that matters here: if the inward pass inherited the - // outward pass's covariance it would consume every measurement twice, the - // estimate would be too confident by root two, and the pull widths would - // come out near 0.707 while every other test in this file still passed. - // That number is asserted against explicitly. - const SyntheticChain chain; - const UniformMaterial material; - - constexpr int kTrials = 3000; - const Scalar sigmaLoc = 1.5811388e-3F; // sqrt(2.5e-6) cm - - std::mt19937 rng(20260919); - std::normal_distribution gauss(0., 1.); - - std::array, kNBoundParams> pull; - - for (int trial = 0; trial < kTrials; ++trial) { - auto smeared = chain.measurements; - for (auto& m : smeared) { - m.loc0 += static_cast(gauss(rng) * static_cast(sigmaLoc)); - m.loc1 += static_cast(gauss(rng) * static_cast(sigmaLoc)); - } - - const auto fit = fitMeasurements(chain.propagator, chain.tracking, material, smeared.data(), static_cast(smeared.size())); - BOOST_REQUIRE(fit.status == PropagationStatus::kOk); - - collectPulls(trial, fit.state, chain.truth[0], pull); - } + for (const auto referenceMode : {FitReference::kFilteredState, FitReference::kCircleOutward, FitReference::kIndependentBoth}) { + BOOST_TEST_CONTEXT("reference mode " << static_cast(referenceMode)) + { + // A fit can be self-consistent and still wrong in + // exactly one way that matters here: if the inward pass inherited the + // outward pass's covariance it would consume every measurement twice, the + // estimate would be too confident by root two, and the pull widths would + // come out near 0.707 while every other test in this file still passed. + // That number is asserted against explicitly. + const SyntheticChain chain; + const UniformMaterial material; + + constexpr int kTrials = 3000; + const Scalar sigmaLoc = 1.5811388e-3F; // sqrt(2.5e-6) cm + + std::mt19937 rng(20260919); + std::normal_distribution gauss(0., 1.); + + std::array, kNBoundParams> pull; + + for (int trial = 0; trial < kTrials; ++trial) { + auto smeared = chain.measurements; + for (auto& m : smeared) { + m.loc0 += static_cast(gauss(rng) * static_cast(sigmaLoc)); + m.loc1 += static_cast(gauss(rng) * static_cast(sigmaLoc)); + } - const char* names[kNBoundParams] = {"loc0", "loc1", "phi", "theta", "qOverP"}; - for (int i = 0; i < kNBoundParams; ++i) { - const auto& p = pull[static_cast(i)]; - const auto n = static_cast(p.size()); + const auto fit = fitMeasurementsReference(chain.propagator, chain.tracking, material, smeared.data(), static_cast(smeared.size()), referenceMode); + BOOST_REQUIRE(fit.status == PropagationStatus::kOk); - const auto [meanPull, width] = moments(p); + collectPulls(trial, fit.state, chain.truth[0], pull); + } - BOOST_TEST_CONTEXT("parameter " << names[i] << " pull mean " << meanPull << " width " << width) - { - // The mean of kTrials unit-width pulls has standard error 1/sqrt(N), - // which is 0.018 here; four of those is a comfortable band that still - // catches any real bias. - BOOST_TEST(std::abs(meanPull) < 4. / std::sqrt(n)); - // The width itself is known to 1/sqrt(2N) = 1.3 %, so 10 % is generous - // -- and 0.707 is 29 % low, far outside it. - BOOST_TEST(std::abs(width - 1.) < 0.10); + const char* names[kNBoundParams] = {"loc0", "loc1", "phi", "theta", "qOverP"}; + for (int i = 0; i < kNBoundParams; ++i) { + const auto& p = pull[static_cast(i)]; + const auto n = static_cast(p.size()); + + const auto [meanPull, width] = moments(p); + + BOOST_TEST_CONTEXT("parameter " << names[i] << " pull mean " << meanPull << " width " << width) + { + // The mean of kTrials unit-width pulls has standard error 1/sqrt(N), + // which is 0.018 here; four of those is a comfortable band that still + // catches any real bias. + BOOST_TEST(std::abs(meanPull) < 4. / std::sqrt(n)); + // The width itself is known to 1/sqrt(2N) = 1.3 %, so 10 % is generous + // -- and 0.707 is 29 % low, far outside it. + BOOST_TEST(std::abs(width - 1.) < 0.10); + } + } } } } @@ -547,64 +562,69 @@ BOOST_AUTO_TEST_CASE(ChainSeedOfAStraightTrackHasNoCurvature) /// is. BOOST_AUTO_TEST_CASE(LowMomentumChainsFitWithUnitPulls) { - constexpr Scalar kQOverP = 2.5F; - constexpr Scalar kX2X0PerCm = 1.5e-3F; - constexpr std::size_t kLayers = 4; - BarrelChain chain(kQOverP); - BOOST_REQUIRE(chain.built); - const UniformMaterial material{.x2X0PerCm = kX2X0PerCm}; - - constexpr int kTrials = 1000; - const Scalar sigmaLoc = 1.5811388e-3F; // sqrt(2.5e-6) cm - - std::mt19937 rng(20260923); - std::normal_distribution gauss(0., 1.); - std::array, kNBoundParams> inwardPull; - std::array, kNBoundParams> outwardPull; - double inwardChi2PerDof = 0.; - double outwardChi2PerDof = 0.; - - for (int trial = 0; trial < kTrials; ++trial) { - // A kink can send the track through the 0.5 mm gap between two modules; - // draw again rather than fit a chain with a missing layer. - for (int attempt = 0; attempt < 10; ++attempt) { - chain.generate(kQOverP, kX2X0PerCm, &rng); - if (chain.built) { - break; + for (const auto referenceMode : {FitReference::kFilteredState, FitReference::kCircleOutward, FitReference::kIndependentBoth}) { + BOOST_TEST_CONTEXT("reference mode " << static_cast(referenceMode)) + { + constexpr Scalar kQOverP = 2.5F; + constexpr Scalar kX2X0PerCm = 1.5e-3F; + constexpr std::size_t kLayers = 4; + BarrelChain chain(kQOverP); + BOOST_REQUIRE(chain.built); + const UniformMaterial material{.x2X0PerCm = kX2X0PerCm}; + + constexpr int kTrials = 1000; + const Scalar sigmaLoc = 1.5811388e-3F; // sqrt(2.5e-6) cm + + std::mt19937 rng(20260923); + std::normal_distribution gauss(0., 1.); + std::array, kNBoundParams> inwardPull; + std::array, kNBoundParams> outwardPull; + double inwardChi2PerDof = 0.; + double outwardChi2PerDof = 0.; + + for (int trial = 0; trial < kTrials; ++trial) { + // A kink can send the track through the 0.5 mm gap between two modules; + // draw again rather than fit a chain with a missing layer. + for (int attempt = 0; attempt < 10; ++attempt) { + chain.generate(kQOverP, kX2X0PerCm, &rng); + if (chain.built) { + break; + } + } + BOOST_REQUIRE(chain.built); + auto smeared = chain.measurements; + for (auto& m : smeared) { + m.loc0 += static_cast(gauss(rng) * static_cast(sigmaLoc)); + m.loc1 += static_cast(gauss(rng) * static_cast(sigmaLoc)); + } + const auto fit = fitMeasurementsReference(chain.propagator, chain.tracking, material, smeared.data(), static_cast(smeared.size()), referenceMode); + BOOST_TEST_CONTEXT("trial " << trial << ", propagation status " << static_cast(fit.status)) { BOOST_REQUIRE(fit.status == PropagationStatus::kOk); } + collectPulls(trial, fit.state, chain.truth[0], inwardPull); + collectPulls(trial, fit.outwardState, chain.truth[kLayers - 1], outwardPull); + inwardChi2PerDof += static_cast(fit.chi2) / fit.ndf / kTrials; + outwardChi2PerDof += static_cast(fit.outwardChi2) / fit.ndf / kTrials; } - } - BOOST_REQUIRE(chain.built); - auto smeared = chain.measurements; - for (auto& m : smeared) { - m.loc0 += static_cast(gauss(rng) * static_cast(sigmaLoc)); - m.loc1 += static_cast(gauss(rng) * static_cast(sigmaLoc)); - } - const auto fit = fitMeasurements(chain.propagator, chain.tracking, material, smeared.data(), static_cast(smeared.size())); - BOOST_TEST_CONTEXT("trial " << trial << ", propagation status " << static_cast(fit.status)) { BOOST_REQUIRE(fit.status == PropagationStatus::kOk); } - collectPulls(trial, fit.state, chain.truth[0], inwardPull); - collectPulls(trial, fit.outwardState, chain.truth[kLayers - 1], outwardPull); - inwardChi2PerDof += static_cast(fit.chi2) / fit.ndf / kTrials; - outwardChi2PerDof += static_cast(fit.outwardChi2) / fit.ndf / kTrials; - } - const char* names[kNBoundParams] = {"loc0", "loc1", "phi", "theta", "qOverP"}; - for (const auto* pass : {&inwardPull, &outwardPull}) { - for (int i = 0; i < kNBoundParams; ++i) { - const auto& p = (*pass)[static_cast(i)]; - const auto n = static_cast(p.size()); - const auto [mean, width] = moments(p); - BOOST_TEST_CONTEXT((pass == &inwardPull ? "inward" : "outward") << " parameter " << names[i] << " pull mean " << mean << " width " << width) + const char* names[kNBoundParams] = {"loc0", "loc1", "phi", "theta", "qOverP"}; + for (const auto* pass : {&inwardPull, &outwardPull}) { + for (int i = 0; i < kNBoundParams; ++i) { + const auto& p = (*pass)[static_cast(i)]; + const auto n = static_cast(p.size()); + const auto [mean, width] = moments(p); + BOOST_TEST_CONTEXT((pass == &inwardPull ? "inward" : "outward") << " parameter " << names[i] << " pull mean " << mean << " width " << width) + { + BOOST_TEST(std::abs(mean) < 4. / std::sqrt(n)); + BOOST_TEST(std::abs(width - 1.) < 0.10); + } + } + } + BOOST_TEST_CONTEXT("inward chi2 per spatial dof " << inwardChi2PerDof << ", outward " << outwardChi2PerDof) { - BOOST_TEST(std::abs(mean) < 4. / std::sqrt(n)); - BOOST_TEST(std::abs(width - 1.) < 0.10); + BOOST_TEST(std::abs(inwardChi2PerDof - 1.) < 0.15); + BOOST_TEST(std::abs(outwardChi2PerDof - 1.) < 0.15); } } } - BOOST_TEST_CONTEXT("inward chi2 per spatial dof " << inwardChi2PerDof << ", outward " << outwardChi2PerDof) - { - BOOST_TEST(std::abs(inwardChi2PerDof - 1.) < 0.15); - BOOST_TEST(std::abs(outwardChi2PerDof - 1.) < 0.15); - } } /// The same chains losing energy: the chain seed is a circle through the @@ -617,61 +637,66 @@ BOOST_AUTO_TEST_CASE(LowMomentumChainsFitWithUnitPulls) /// the fit's model of where the loss is taken. BOOST_AUTO_TEST_CASE(LowMomentumChainsLosingEnergyFitWithUnitPulls) { - constexpr Scalar kQOverP = 2.5F; - constexpr Scalar kX2X0PerCm = 1.5e-3F; - constexpr Scalar kRho = 0.3F; // about 2.5 % of the momentum over the chain - constexpr std::size_t kLayers = 4; - BarrelChain chain(kQOverP); - BOOST_REQUIRE(chain.built); - const UniformMaterial material{.x2X0PerCm = kX2X0PerCm, .rho = kRho}; - - constexpr int kTrials = 1000; - const Scalar sigmaLoc = 1.5811388e-3F; - std::mt19937 rng(20261021); - std::normal_distribution gauss(0., 1.); - std::array, kNBoundParams> inwardPull; - double inwardChi2PerDof = 0.; - double outwardChi2PerDof = 0.; - double lostShare = 0.; - for (int trial = 0; trial < kTrials; ++trial) { - for (int attempt = 0; attempt < 10; ++attempt) { - chain.generate(kQOverP, kX2X0PerCm, &rng, kRho); - if (chain.built) { - break; - } - } - BOOST_REQUIRE(chain.built); - lostShare += (1. - static_cast(kQOverP / chain.truth[kLayers - 1].qOverP())) / kTrials; - auto smeared = chain.measurements; - for (auto& m : smeared) { - m.loc0 += static_cast(gauss(rng) * static_cast(sigmaLoc)); - m.loc1 += static_cast(gauss(rng) * static_cast(sigmaLoc)); - } - const auto fit = fitMeasurements(chain.propagator, chain.tracking, material, smeared.data(), static_cast(smeared.size())); - BOOST_TEST_CONTEXT("trial " << trial << ", propagation status " << static_cast(fit.status)) { BOOST_REQUIRE(fit.status == PropagationStatus::kOk); } - collectPulls(trial, fit.state, chain.truth[0], inwardPull); - inwardChi2PerDof += static_cast(fit.chi2) / fit.ndf / kTrials; - outwardChi2PerDof += static_cast(fit.outwardChi2) / fit.ndf / kTrials; - } - BOOST_TEST_MESSAGE("momentum lost over the chain " << lostShare << ", chi2 per dof inward " << inwardChi2PerDof << " outward " << outwardChi2PerDof); - BOOST_TEST(lostShare > 0.01); - for (int i = 0; i < kNBoundParams; ++i) { - const auto& p = inwardPull[static_cast(i)]; - const bool held = i != o2::gtrk::kQOverP; - const auto n = static_cast(p.size()); - const auto [mean, width] = moments(p); - BOOST_TEST_MESSAGE("inward parameter " << i << " pull mean " << mean << " width " << width); - BOOST_TEST_CONTEXT("inward parameter " << i << " pull mean " << mean << " width " << width) + for (const auto referenceMode : {FitReference::kFilteredState, FitReference::kCircleOutward, FitReference::kIndependentBoth}) { + BOOST_TEST_CONTEXT("reference mode " << static_cast(referenceMode)) { - BOOST_TEST((!held || std::abs(mean) < 4. / std::sqrt(n))); - BOOST_TEST((!held || std::abs(width - 1.) < 0.10)); + constexpr Scalar kQOverP = 2.5F; + constexpr Scalar kX2X0PerCm = 1.5e-3F; + constexpr Scalar kRho = 0.3F; // about 2.5 % of the momentum over the chain + constexpr std::size_t kLayers = 4; + BarrelChain chain(kQOverP); + BOOST_REQUIRE(chain.built); + const UniformMaterial material{.x2X0PerCm = kX2X0PerCm, .rho = kRho}; + + constexpr int kTrials = 1000; + const Scalar sigmaLoc = 1.5811388e-3F; + std::mt19937 rng(20261021); + std::normal_distribution gauss(0., 1.); + std::array, kNBoundParams> inwardPull; + double inwardChi2PerDof = 0.; + double outwardChi2PerDof = 0.; + double lostShare = 0.; + for (int trial = 0; trial < kTrials; ++trial) { + for (int attempt = 0; attempt < 10; ++attempt) { + chain.generate(kQOverP, kX2X0PerCm, &rng, kRho); + if (chain.built) { + break; + } + } + BOOST_REQUIRE(chain.built); + lostShare += (1. - static_cast(kQOverP / chain.truth[kLayers - 1].qOverP())) / kTrials; + auto smeared = chain.measurements; + for (auto& m : smeared) { + m.loc0 += static_cast(gauss(rng) * static_cast(sigmaLoc)); + m.loc1 += static_cast(gauss(rng) * static_cast(sigmaLoc)); + } + const auto fit = fitMeasurementsReference(chain.propagator, chain.tracking, material, smeared.data(), static_cast(smeared.size()), referenceMode); + BOOST_TEST_CONTEXT("trial " << trial << ", propagation status " << static_cast(fit.status)) { BOOST_REQUIRE(fit.status == PropagationStatus::kOk); } + collectPulls(trial, fit.state, chain.truth[0], inwardPull); + inwardChi2PerDof += static_cast(fit.chi2) / fit.ndf / kTrials; + outwardChi2PerDof += static_cast(fit.outwardChi2) / fit.ndf / kTrials; + } + BOOST_TEST_MESSAGE("momentum lost over the chain " << lostShare << ", chi2 per dof inward " << inwardChi2PerDof << " outward " << outwardChi2PerDof); + BOOST_TEST(lostShare > 0.01); + for (int i = 0; i < kNBoundParams; ++i) { + const auto& p = inwardPull[static_cast(i)]; + const bool held = i != o2::gtrk::kQOverP; + const auto n = static_cast(p.size()); + const auto [mean, width] = moments(p); + BOOST_TEST_MESSAGE("inward parameter " << i << " pull mean " << mean << " width " << width); + BOOST_TEST_CONTEXT("inward parameter " << i << " pull mean " << mean << " width " << width) + { + BOOST_TEST((!held || std::abs(mean) < 4. / std::sqrt(n))); + BOOST_TEST((!held || std::abs(width - 1.) < 0.10)); + } + } + BOOST_TEST_CONTEXT("inward chi2 per spatial dof " << inwardChi2PerDof << ", outward " << outwardChi2PerDof) + { + BOOST_TEST(std::abs(inwardChi2PerDof - 1.) < 0.15); + BOOST_TEST(std::abs(outwardChi2PerDof - 1.) < 0.15); + } } } - BOOST_TEST_CONTEXT("inward chi2 per spatial dof " << inwardChi2PerDof << ", outward " << outwardChi2PerDof) - { - BOOST_TEST(std::abs(inwardChi2PerDof - 1.) < 0.15); - BOOST_TEST(std::abs(outwardChi2PerDof - 1.) < 0.15); - } } BOOST_AUTO_TEST_CASE(PerigeePhiIsTheDirectionAtTheTransverseClosestApproach) diff --git a/Detectors/GenericTracking/workflow/include/GenericTracking/Workflow/GridConfig.h b/Detectors/GenericTracking/workflow/include/GenericTracking/Workflow/GridConfig.h index 1132b2e8cbe6d..91154eb704a8c 100644 --- a/Detectors/GenericTracking/workflow/include/GenericTracking/Workflow/GridConfig.h +++ b/Detectors/GenericTracking/workflow/include/GenericTracking/Workflow/GridConfig.h @@ -108,6 +108,9 @@ namespace o2::gtrk::workflow template inline void checkRecoParam(const RecoSettings& param) { + if (param.fitReference < 0 || param.fitReference > 2) { + throw std::invalid_argument("GTRKReco.fitReference must be 0 (filtered state), 1 (circle outward), or 2 (circle outward, outward-fit reference inward)"); + } const auto margins = param.frameMarginBC.size(); if (margins != 0U && margins != static_cast(Detector::nLayers)) { throw std::invalid_argument(std::format("{}.frameMarginBC has {} entries for {}'s {} layers", RecoParamOf::type::Instance().getName(), margins, Det::name, Detector::nLayers)); @@ -253,6 +256,7 @@ inline void checkRecoParam(const RecoSettings& param) .vertexSource = entry(param.iterationVertexSource, static_cast(defaults.vertexSource)) == 0 ? VertexSource::kBeamLine : VertexSource::kVertexList, .minClusters = entry(param.iterationMinClusters, defaults.minClusters), .maxKinks = entry(param.iterationMaxKinks, defaults.maxKinks), + .fitReference = static_cast(param.fitReference), .kinkSigma = kinkSigma, .graphKinkSigma = entryDouble(param.iterationGraphKinkSigma, defaults.graphKinkSigma), .skipLinks = flag(param.iterationSkipLinks, defaults.skipLinks), diff --git a/Detectors/GenericTracking/workflow/include/GenericTracking/Workflow/RecoParam.h b/Detectors/GenericTracking/workflow/include/GenericTracking/Workflow/RecoParam.h index 497b523d82e44..d4468644f8d33 100644 --- a/Detectors/GenericTracking/workflow/include/GenericTracking/Workflow/RecoParam.h +++ b/Detectors/GenericTracking/workflow/include/GenericTracking/Workflow/RecoParam.h @@ -64,6 +64,12 @@ struct RecoSettings { /// scattering, in the windows at pT_min and in the fit, and the fit's energy /// loss. The charged pion's by default. double massHypothesisGeV = o2::constants::physics::MassPionCharged; + /// Linearisation: 0 uses the evolving filtered state both ways (default). + /// 1 uses an independent circle reference outward and the filtered state inward. + /// 2 uses the circle reference outward and an independent reference seeded + /// from the outward fitted parameters inward. References are position-anchored + /// to clusters within each pass. Selected tracks reuse the road's outward fit. + int fitReference = 0; /// The uniform material per layer [X/X0] that sizes the windows and the /// grid, and the scattering when no material LUT is given; with a LUT the /// tracking takes its budgets.