From 4c28bdc483b467a153e821d94bf3a0217b107928 Mon Sep 17 00:00:00 2001 From: Gabriele Cimador Date: Fri, 9 Oct 2026 09:42:13 +0200 Subject: [PATCH] ITS GPU seeding vertexer: add clamps to adaptive low mult threshold, and propate cells towards beamline --- .../GPU/ITStrackingGPU/TrackingKernels.h | 2 + .../tracking/GPU/cuda/TrackerTraitsGPU.cxx | 9 ++- .../ITS/tracking/GPU/cuda/TrackingKernels.cu | 16 +++-- .../include/ITStracking/ClusterLines.h | 20 ++++++ .../include/ITStracking/Configuration.h | 2 + .../include/ITStracking/VertexUtils.h | 64 +++++++++++++++++++ .../ITSMFT/ITS/tracking/src/Configuration.cxx | 2 + .../ITSMFT/ITS/tracking/src/TrackerTraits.cxx | 11 ++-- .../include/ITSMFTTracking/Constants.h | 6 +- .../ITSMFTTracking/ITSTrackingConfigParam.h | 12 ++-- 10 files changed, 127 insertions(+), 17 deletions(-) diff --git a/Detectors/ITSMFT/ITS/tracking/GPU/ITStrackingGPU/TrackingKernels.h b/Detectors/ITSMFT/ITS/tracking/GPU/ITStrackingGPU/TrackingKernels.h index 45617bc00b502..5c25bdc833890 100644 --- a/Detectors/ITSMFT/ITS/tracking/GPU/ITStrackingGPU/TrackingKernels.h +++ b/Detectors/ITSMFT/ITS/tracking/GPU/ITStrackingGPU/TrackingKernels.h @@ -211,6 +211,8 @@ struct TrackingKernels { int* lineSlots, const float beamX, const float beamY, + const float bz, + const float curvatureScale, const float maxZ, const float minPt, float* linesZs, diff --git a/Detectors/ITSMFT/ITS/tracking/GPU/cuda/TrackerTraitsGPU.cxx b/Detectors/ITSMFT/ITS/tracking/GPU/cuda/TrackerTraitsGPU.cxx index cda6894b0356a..c28ba6630808f 100644 --- a/Detectors/ITSMFT/ITS/tracking/GPU/cuda/TrackerTraitsGPU.cxx +++ b/Detectors/ITSMFT/ITS/tracking/GPU/cuda/TrackerTraitsGPU.cxx @@ -266,6 +266,8 @@ void TrackerTraitsGPU::computeVertexCandidates(const int iteration) mTimeFrameGPU->getDeviceLineSlots(), mTimeFrameGPU->getBeamX(), mTimeFrameGPU->getBeamY(), + this->getBz(), + this->mTrkParams[iteration].VtxLineCurvatureScale, this->mTrkParams[iteration].VtxMaxZPositionAllowed, this->mTrkParams[iteration].VtxLineMinPt, mTimeFrameGPU->getDeviceLineZs(), @@ -451,13 +453,13 @@ void TrackerTraitsGPU::computeVertices(const int iteration) } } } - const float sigThreshold = goodSig > 0.f ? goodSig * std::sqrt(static_cast(std::max(rofLoad, 1.))) : 0.f; + const float debrisThreshold = goodSig > 0.f ? getDebrisThreshold(goodSig, rofLoad, suppressLowMultDebris) : 0.f; for (const int p : accepted) { const auto& c = cands[p]; if (!rofVertices[rofId].empty()) { if (goodSig > 0.f) { - if (c.nGood <= sigThreshold) { + if (c.nGood < debrisThreshold) { continue; } } else if (c.size < suppressLowMultDebris) { @@ -503,6 +505,9 @@ void TrackerTraitsGPU::computeVertices(const int iteration) } } } + if (!this->mTrkParams[iteration].PassFlags[IterationStep::MarkVerticesAsUPC]) { // UPC ROFs are near-empty by construction + pruneOverpopulatedRofs(rofVertices, rofLabels, this->mTrkParams[iteration].VtxOverpopulatedRofNSigma, suppressLowMultDebris, constants::VtxOverpopulatedRofTrimFraction); + } for (int rofId = 0; rofId < nRofs; ++rofId) { for (auto& vertex : rofVertices[rofId]) { diff --git a/Detectors/ITSMFT/ITS/tracking/GPU/cuda/TrackingKernels.cu b/Detectors/ITSMFT/ITS/tracking/GPU/cuda/TrackingKernels.cu index 83607580c7370..43ee8436892d4 100644 --- a/Detectors/ITSMFT/ITS/tracking/GPU/cuda/TrackingKernels.cu +++ b/Detectors/ITSMFT/ITS/tracking/GPU/cuda/TrackingKernels.cu @@ -874,6 +874,8 @@ GPUg() void dedupCellsKernel( const int ownedClustersCut, const float beamX, const float beamY, + const float bz, + const float curvatureScale, const float maxZ, const float minPt, int* cellAccepted) @@ -881,7 +883,7 @@ GPUg() void dedupCellsKernel( for (int i = blockIdx.x * blockDim.x + threadIdx.x; i < nCells; i += blockDim.x * gridDim.x) { const CellSeed& cell = cells[i]; std::array origin, direction; - if (!cell.getPxPyPzGlo(direction)) { + if (!getCellLineAtBeam(cell, beamX, beamY, bz, curvatureScale, origin, direction)) { cellAccepted[i] = 0; continue; } @@ -889,7 +891,6 @@ GPUg() void dedupCellsKernel( const bool owned1 = static_cast(clusterOwners[1][cell.getSecondClusterIndex()]) == static_cast(i); const bool owned2 = static_cast(clusterOwners[2][cell.getThirdClusterIndex()]) == static_cast(i); const bool keepCell = (static_cast(owned0) + static_cast(owned1) + static_cast(owned2)) >= 3 - ownedClustersCut; - cell.getXYZGlo(origin); const float dx = origin[0] - beamX; const float dy = origin[1] - beamY; const float den = direction[0] * direction[0] + direction[1] * direction[1]; @@ -910,6 +911,8 @@ GPUg() void linearizeCellsKernel( int* lineRof, const float beamX, const float beamY, + const float bz, + const float curvatureScale, float* lineZs, o2::its::TimeEstBC* lineTimes, int* lineClusters, // 3 per line (L0,L1,L2 cluster ids), for the host-side MC label derivation @@ -923,8 +926,7 @@ GPUg() void linearizeCellsKernel( } const CellSeed& cell = cells[i]; std::array origin, direction; - cell.getXYZGlo(origin); - cell.getPxPyPzGlo(direction); + getCellLineAtBeam(cell, beamX, beamY, bz, curvatureScale, origin, direction); // accepted by dedupCellsKernel: succeeds lines[slot] = o2::its::Line{origin.data(), direction.data(), cell.getTimeStamp()}; lineRof[slot] = deviceUpperBound(rofFramesClustersL1, 0, nRofsL1 + 1, cell.getSecondClusterIndex()) - 1; float zAtBeam; @@ -1577,6 +1579,8 @@ void TrackingKernels::linearizeCellsToLinesHandler(const int nCells, int* lineSlots, // nCells + 1 scratch: accept flags, scanned in place into slots const float beamX, const float beamY, + const float bz, + const float curvatureScale, const float maxZ, const float minPt, float* linesZs, @@ -1593,6 +1597,8 @@ void TrackingKernels::linearizeCellsToLinesHandler(const int nCells, ownedClustersCut, beamX, beamY, + bz, + curvatureScale, maxZ, minPt, lineSlots); @@ -1608,6 +1614,8 @@ void TrackingKernels::linearizeCellsToLinesHandler(const int nCells, lineRof, beamX, beamY, + bz, + curvatureScale, linesZs, lineTimes, lineClusters, diff --git a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/ClusterLines.h b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/ClusterLines.h index d63c42b54736c..f1255a4ee9d1f 100644 --- a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/ClusterLines.h +++ b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/ClusterLines.h @@ -19,6 +19,7 @@ #include #endif #include "ITStracking/Cluster.h" +#include "ITStracking/Cell.h" #include "ITSMFTTracking/Constants.h" #include "ITStracking/Tracklet.h" #include "GPUCommonDef.h" @@ -121,6 +122,25 @@ struct Line final { TimeEstBC mTime; }; +GPUdi() bool getCellLineAtBeam(const CellSeed& cell, const float beamX, const float beamY, const float bz, const float curvatureScale, std::array& origin, std::array& direction) +{ + cell.getXYZGlo(origin); + if (!cell.getPxPyPzGlo(direction)) { + return false; + } + o2::track::TrackParametrization par{cell}; + par.setQ2Pt(par.getQ2Pt() * curvatureScale); + std::array dcaOrigin, dcaDirection; + if (par.propagateParamToDCA({beamX, beamY, 0.f}, bz)) { + par.getXYZGlo(dcaOrigin); + if (par.getPxPyPzGlo(dcaDirection)) { + origin = dcaOrigin; + direction = dcaDirection; + } + } + return true; +} + /// Least-squares vertex fit over a set of lines (the normal equations AX = -B). class ClusterLines final { diff --git a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/Configuration.h b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/Configuration.h index a253be93e54c3..61a681bc30cfe 100644 --- a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/Configuration.h +++ b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/Configuration.h @@ -151,6 +151,7 @@ struct TrackingParameters { int VertPerRofThreshold = 0; // max vertices in a ROF for the UPC pass to still run on it float VtxPhiCut = -1.f; float VtxLineMinPt = -1.f; + float VtxLineCurvatureScale = 0.5f; float VtxMaxZPositionAllowed = -1.f; float VtxClusterCut = -1.f; float VtxPairCut = -1.f; @@ -165,6 +166,7 @@ struct TrackingParameters { float VtxGoodContributorsSignificance = -1.f; int VtxClusterContributorsCut = -1; int VtxSuppressLowMultDebris = -1; + float VtxOverpopulatedRofNSigma = -1.f; }; struct VertexingParameters { diff --git a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/VertexUtils.h b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/VertexUtils.h index eb0b843b667cc..b8b5174c60455 100644 --- a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/VertexUtils.h +++ b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/VertexUtils.h @@ -19,9 +19,14 @@ #include "SimulationDataFormat/MCCompLabel.h" #include "ITStracking/Configuration.h" +#include "Framework/Logger.h" + +#include +#include #include #include #include +#include namespace o2::its { @@ -60,6 +65,65 @@ inline Vertex makeDiamondVertex(const TrackingParameters& trkParam) return diamond; } +/// Good lines a further vertex needs to not count as debris in a ROF that already has one: goodSig * sqrt(ROF load), clamped to +/// [constants::VtxMinGoodThreshold, suppressLowMultDebris] unless the debris cut is off (UPC pass). +inline float getDebrisThreshold(const float goodSig, const double rofLoad, const int suppressLowMultDebris) +{ + const float threshold = goodSig * std::sqrt(static_cast(std::max(rofLoad, 1.))); + if (suppressLowMultDebris < constants::VtxMinGoodThreshold) { + return threshold; + } + return std::clamp(threshold, constants::VtxMinGoodThreshold, static_cast(suppressLowMultDebris)); +} + +/// Caps ROFs whose vertex count is an outlier of the TF +template +int pruneOverpopulatedRofs(std::vector& rofVertices, std::vector& rofLabels, const float nSigma, const int minContributors, const float trimFraction) +{ + const int nRofs = static_cast(rofVertices.size()); + if (nSigma <= 0.f || nRofs == 0) { + return 0; + } + std::vector counts(nRofs); + for (int r = 0; r < nRofs; ++r) { + counts[r] = static_cast(rofVertices[r].size()); + } + std::vector sorted(counts); + std::sort(sorted.begin(), sorted.end()); + const int nUsed = std::max(1, nRofs - std::max(1, static_cast(trimFraction * nRofs))); + double sum = 0.; + for (int i = 0; i < nUsed; ++i) { + sum += sorted[i]; + } + const double mean = sum / nUsed; + const double threshold = mean + nSigma * std::sqrt(mean + 1.); + int removed = 0; + for (int r = 0; r < nRofs; ++r) { + if (counts[r] <= threshold) { + continue; + } + auto& vtx = rofVertices[r]; + const bool withLabels = static_cast(rofLabels.size()) == nRofs && rofLabels[r].size() == vtx.size(); + size_t out = 1; // the largest vertex always stays + for (size_t i = 1; i < vtx.size(); ++i) { + if (vtx[i].getNContributors() >= minContributors) { + vtx[out] = vtx[i]; + if (withLabels) { + rofLabels[r][out] = rofLabels[r][i]; + } + ++out; + } + } + LOGP(info, "Seeding vertexer: overpopulated ROF {} pruned {} -> {} vertices (threshold {:.1f}, mean {:.2f} per ROF)", r, vtx.size(), out, threshold, mean); + removed += static_cast(vtx.size() - out); + vtx.erase(vtx.begin() + out, vtx.end()); + if (withLabels) { + rofLabels[r].erase(rofLabels[r].begin() + out, rofLabels[r].end()); + } + } + return removed; +} + } // namespace o2::its #endif /* O2_ITS_TRACKING_VERTEXUTILS_H_ */ diff --git a/Detectors/ITSMFT/ITS/tracking/src/Configuration.cxx b/Detectors/ITSMFT/ITS/tracking/src/Configuration.cxx index 149c19ba83343..302969b4c31f2 100644 --- a/Detectors/ITSMFT/ITS/tracking/src/Configuration.cxx +++ b/Detectors/ITSMFT/ITS/tracking/src/Configuration.cxx @@ -308,6 +308,7 @@ std::vector TrackingMode::getTrackingParameters(TrackingMode seedingPass.VertPerRofThreshold = vc.vertPerRofThreshold; seedingPass.VtxPhiCut = vc.phiCut; seedingPass.VtxLineMinPt = vc.lineMinPt; + seedingPass.VtxLineCurvatureScale = vc.lineCurvatureScale; seedingPass.VtxMaxZPositionAllowed = vc.maxZPositionAllowed; seedingPass.VtxClusterCut = vc.clusterCut; seedingPass.VtxPairCut = vc.pairCut; @@ -322,6 +323,7 @@ std::vector TrackingMode::getTrackingParameters(TrackingMode seedingPass.VtxGoodContributorsSignificance = vc.goodContributorsSignificance; seedingPass.VtxClusterContributorsCut = vc.clusterContributorsCut; seedingPass.VtxSuppressLowMultDebris = vc.suppressLowMultDebris; + seedingPass.VtxOverpopulatedRofNSigma = vc.overpopulatedRofNSigma; std::vector seedingPasses; seedingPasses.push_back(seedingPass); diff --git a/Detectors/ITSMFT/ITS/tracking/src/TrackerTraits.cxx b/Detectors/ITSMFT/ITS/tracking/src/TrackerTraits.cxx index ea459e0bd311a..8365fab5f74c6 100644 --- a/Detectors/ITSMFT/ITS/tracking/src/TrackerTraits.cxx +++ b/Detectors/ITSMFT/ITS/tracking/src/TrackerTraits.cxx @@ -305,6 +305,7 @@ void TrackerTraits::computeVertexCandidates(const int iteration) const float lineMinPt = mTrkParams[iteration].VtxLineMinPt; const float beamX = mTimeFrame->getBeamX(); const float beamY = mTimeFrame->getBeamY(); + const float lineCurvatureScale = mTrkParams[iteration].VtxLineCurvatureScale; auto makeKey = [](float attribute, int cellIdx) -> size_t { const uint32_t attributeInt = std::bit_cast(attribute); return (static_cast(attributeInt) << 32) | static_cast(cellIdx); @@ -327,8 +328,7 @@ void TrackerTraits::computeVertexCandidates(const int iteration) kCl1[k] = c1; kCl2[k] = cell.getThirdClusterIndex(); std::array origin, direction; - cell.getXYZGlo(origin); - if (!cell.getPxPyPzGlo(direction)) { + if (!getCellLineAtBeam(cell, beamX, beamY, getBz(), lineCurvatureScale, origin, direction)) { return; } kGeomOk[k] = 1; @@ -711,11 +711,11 @@ void TrackerTraits::computeVertices(const int iteration) } } } - const float sigThreshold = goodSig > 0.f ? goodSig * std::sqrt(static_cast(std::max(rofLoad, 1.))) : 0.f; + const float debrisThreshold = goodSig > 0.f ? getDebrisThreshold(goodSig, rofLoad, suppressLowMultDebris) : 0.f; for (const int p : accepted) { if (!rofVertices[rofId].empty()) { if (goodSig > 0.f) { - if (nGoodCand[p] <= sigThreshold) { + if (nGoodCand[p] < debrisThreshold) { continue; } } else if (static_cast(cand[p].getSize()) < suppressLowMultDebris) { @@ -758,6 +758,9 @@ void TrackerTraits::computeVertices(const int iteration) }); }); } + if (!tp.PassFlags[IterationStep::MarkVerticesAsUPC]) { // UPC ROFs are near-empty by construction + pruneOverpopulatedRofs(rofVertices, rofLabels, tp.VtxOverpopulatedRofNSigma, suppressLowMultDebris, constants::VtxOverpopulatedRofTrimFraction); + } for (int rofId{0}; rofId < nRofs; ++rofId) { for (auto& vertex : rofVertices[rofId]) { mTimeFrame->addPrimaryVertex(vertex); diff --git a/Detectors/ITSMFT/common/tracking/include/ITSMFTTracking/Constants.h b/Detectors/ITSMFT/common/tracking/include/ITSMFTTracking/Constants.h index 26f7e187c13eb..e3ff1d8cf87b2 100644 --- a/Detectors/ITSMFT/common/tracking/include/ITSMFTTracking/Constants.h +++ b/Detectors/ITSMFT/common/tracking/include/ITSMFTTracking/Constants.h @@ -41,8 +41,10 @@ constexpr int MinNumberOfConcurrentSeeds = (1 << 8); // minimum chunk size for constexpr int MaxNumberOfConcurrentSeeds = (1 << 12); // maximum chunk size for a worker for the final track fit/extraploation step constexpr float MaxTrackSeedQ2Pt = 1.e3f; // maximum q/pt for track seeds -constexpr int MaxBootstrapPasses = 5; // beam bootstrap: cap on the re-trackleting passes -constexpr float BeamConvergence2 = 5.e-3f * 5.e-3f; // beam bootstrap: stop below a (50 um)^2 beam shift +constexpr int MaxBootstrapPasses = 5; // beam bootstrap: cap on the re-trackleting passes +constexpr float BeamConvergence2 = 5.e-3f * 5.e-3f; // beam bootstrap: stop below a (50 um)^2 beam shift +constexpr float VtxMinGoodThreshold = 2.f; // seeding emit: floor of the k*sqrt(ROF load) debris threshold (bounds included) +constexpr float VtxOverpopulatedRofTrimFraction = 0.02f; // overpopulated-ROF pruning: busiest fraction of ROFs left out of the per-ROF mean namespace helpers { diff --git a/Detectors/ITSMFT/common/tracking/include/ITSMFTTracking/ITSTrackingConfigParam.h b/Detectors/ITSMFT/common/tracking/include/ITSMFTTracking/ITSTrackingConfigParam.h index 589fd7632a4c3..ae27d432039d7 100644 --- a/Detectors/ITSMFT/common/tracking/include/ITSMFTTracking/ITSTrackingConfigParam.h +++ b/Detectors/ITSMFT/common/tracking/include/ITSMFTTracking/ITSTrackingConfigParam.h @@ -42,15 +42,17 @@ struct VertexerParamConfig : public o2::conf::ConfigurableParamHelper= suppressLowMultDebris contributors (not in the UPC pass); <=0 disables float duplicateZScale = 0.7f; // per-candidate dedup radius scale/sqrt(size); <=0 uses duplicateZCut int seedMemberRadiusTime = 0; int seedMemberRadiusZ = 2;