From e48c4bd87d933b160ff1a942d2cd0af8ab18b5e3 Mon Sep 17 00:00:00 2001 From: Maximiliano Puccio Date: Wed, 25 Mar 2026 15:20:17 +0100 Subject: [PATCH 1/3] Move FastMultEstimation to ITS tracking library --- Detectors/ITSMFT/ITS/reconstruction/CMakeLists.txt | 6 +----- .../ITS/reconstruction/src/ITSReconstructionLinkDef.h | 3 --- Detectors/ITSMFT/ITS/tracking/CMakeLists.txt | 5 +++++ .../include/ITStracking}/FastMultEst.h | 2 +- .../include/ITStracking}/FastMultEstConfig.h | 0 .../ITS/{reconstruction => tracking}/src/FastMultEst.cxx | 2 +- .../{reconstruction => tracking}/src/FastMultEstConfig.cxx | 2 +- Detectors/ITSMFT/ITS/tracking/src/TrackingInterface.cxx | 4 ++-- Detectors/ITSMFT/ITS/tracking/src/TrackingLinkDef.h | 4 ++++ 9 files changed, 15 insertions(+), 13 deletions(-) rename Detectors/ITSMFT/ITS/{reconstruction/include/ITSReconstruction => tracking/include/ITStracking}/FastMultEst.h (98%) rename Detectors/ITSMFT/ITS/{reconstruction/include/ITSReconstruction => tracking/include/ITStracking}/FastMultEstConfig.h (100%) rename Detectors/ITSMFT/ITS/{reconstruction => tracking}/src/FastMultEst.cxx (99%) rename Detectors/ITSMFT/ITS/{reconstruction => tracking}/src/FastMultEstConfig.cxx (94%) diff --git a/Detectors/ITSMFT/ITS/reconstruction/CMakeLists.txt b/Detectors/ITSMFT/ITS/reconstruction/CMakeLists.txt index d2126be1da2c6..be42015b95795 100644 --- a/Detectors/ITSMFT/ITS/reconstruction/CMakeLists.txt +++ b/Detectors/ITSMFT/ITS/reconstruction/CMakeLists.txt @@ -11,8 +11,6 @@ o2_add_library(ITSReconstruction SOURCES src/RecoGeomHelper.cxx - src/FastMultEstConfig.cxx - src/FastMultEst.cxx PUBLIC_LINK_LIBRARIES O2::ITSBase O2::ITSMFTReconstruction O2::DataFormatsITS @@ -20,6 +18,4 @@ o2_add_library(ITSReconstruction o2_target_root_dictionary( ITSReconstruction - HEADERS include/ITSReconstruction/RecoGeomHelper.h - include/ITSReconstruction/FastMultEst.h - include/ITSReconstruction/FastMultEstConfig.h) + HEADERS include/ITSReconstruction/RecoGeomHelper.h) diff --git a/Detectors/ITSMFT/ITS/reconstruction/src/ITSReconstructionLinkDef.h b/Detectors/ITSMFT/ITS/reconstruction/src/ITSReconstructionLinkDef.h index 67622303fc840..3bc8ee0f5403b 100644 --- a/Detectors/ITSMFT/ITS/reconstruction/src/ITSReconstructionLinkDef.h +++ b/Detectors/ITSMFT/ITS/reconstruction/src/ITSReconstructionLinkDef.h @@ -16,8 +16,5 @@ #pragma link off all functions; #pragma link C++ class o2::its::RecoGeomHelper + ; -#pragma link C++ class o2::its::FastMultEst + ; -#pragma link C++ class o2::its::FastMultEstConfig + ; -#pragma link C++ class o2::conf::ConfigurableParamHelper < o2::its::FastMultEstConfig> + ; #endif diff --git a/Detectors/ITSMFT/ITS/tracking/CMakeLists.txt b/Detectors/ITSMFT/ITS/tracking/CMakeLists.txt index d341043bf7e37..8d8304d16764f 100644 --- a/Detectors/ITSMFT/ITS/tracking/CMakeLists.txt +++ b/Detectors/ITSMFT/ITS/tracking/CMakeLists.txt @@ -14,6 +14,8 @@ o2_add_library(ITStracking SOURCES src/ClusterLines.cxx src/Cluster.cxx src/Configuration.cxx + src/FastMultEstConfig.cxx + src/FastMultEst.cxx src/TimeFrame.cxx src/IOUtils.cxx src/Tracker.cxx @@ -28,6 +30,7 @@ o2_add_library(ITStracking O2::DataFormatsITSMFT O2::SimulationDataFormat O2::ITSBase + O2::CommonUtils O2::ITSReconstruction O2::ITSMFTReconstruction O2::DataFormatsITS @@ -49,6 +52,8 @@ o2_target_root_dictionary(ITStracking include/ITStracking/Tracklet.h include/ITStracking/Cluster.h include/ITStracking/Definitions.h + include/ITStracking/FastMultEst.h + include/ITStracking/FastMultEstConfig.h include/ITStracking/TrackingConfigParam.h LINKDEF src/TrackingLinkDef.h) diff --git a/Detectors/ITSMFT/ITS/reconstruction/include/ITSReconstruction/FastMultEst.h b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/FastMultEst.h similarity index 98% rename from Detectors/ITSMFT/ITS/reconstruction/include/ITSReconstruction/FastMultEst.h rename to Detectors/ITSMFT/ITS/tracking/include/ITStracking/FastMultEst.h index 9e8299e89b404..1ae8713e75a05 100644 --- a/Detectors/ITSMFT/ITS/reconstruction/include/ITSReconstruction/FastMultEst.h +++ b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/FastMultEst.h @@ -20,7 +20,7 @@ #include "DataFormatsITSMFT/ROFRecord.h" #include "DataFormatsITSMFT/CompCluster.h" #include -#include "ITSReconstruction/FastMultEstConfig.h" +#include "ITStracking/FastMultEstConfig.h" #include #include diff --git a/Detectors/ITSMFT/ITS/reconstruction/include/ITSReconstruction/FastMultEstConfig.h b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/FastMultEstConfig.h similarity index 100% rename from Detectors/ITSMFT/ITS/reconstruction/include/ITSReconstruction/FastMultEstConfig.h rename to Detectors/ITSMFT/ITS/tracking/include/ITStracking/FastMultEstConfig.h diff --git a/Detectors/ITSMFT/ITS/reconstruction/src/FastMultEst.cxx b/Detectors/ITSMFT/ITS/tracking/src/FastMultEst.cxx similarity index 99% rename from Detectors/ITSMFT/ITS/reconstruction/src/FastMultEst.cxx rename to Detectors/ITSMFT/ITS/tracking/src/FastMultEst.cxx index b62cacc28d9a7..4b854be552f80 100644 --- a/Detectors/ITSMFT/ITS/reconstruction/src/FastMultEst.cxx +++ b/Detectors/ITSMFT/ITS/tracking/src/FastMultEst.cxx @@ -13,7 +13,7 @@ /// \brief Fast multiplicity estimator for ITS /// \author ruben.shahoyan@cern.ch -#include "ITSReconstruction/FastMultEst.h" +#include "ITStracking/FastMultEst.h" #include "DataFormatsITSMFT/DPLAlpideParam.h" #include "Framework/Logger.h" #include diff --git a/Detectors/ITSMFT/ITS/reconstruction/src/FastMultEstConfig.cxx b/Detectors/ITSMFT/ITS/tracking/src/FastMultEstConfig.cxx similarity index 94% rename from Detectors/ITSMFT/ITS/reconstruction/src/FastMultEstConfig.cxx rename to Detectors/ITSMFT/ITS/tracking/src/FastMultEstConfig.cxx index 63c43cf26ba15..1568d8ed9f9fb 100644 --- a/Detectors/ITSMFT/ITS/reconstruction/src/FastMultEstConfig.cxx +++ b/Detectors/ITSMFT/ITS/tracking/src/FastMultEstConfig.cxx @@ -9,7 +9,7 @@ // granted to it by virtue of its status as an Intergovernmental Organization // or submit itself to any jurisdiction. -#include "ITSReconstruction/FastMultEstConfig.h" +#include "ITStracking/FastMultEstConfig.h" #include "TRandom.h" O2ParamImpl(o2::its::FastMultEstConfig); diff --git a/Detectors/ITSMFT/ITS/tracking/src/TrackingInterface.cxx b/Detectors/ITSMFT/ITS/tracking/src/TrackingInterface.cxx index 375a35e0d55be..2d46de523e493 100644 --- a/Detectors/ITSMFT/ITS/tracking/src/TrackingInterface.cxx +++ b/Detectors/ITSMFT/ITS/tracking/src/TrackingInterface.cxx @@ -17,8 +17,8 @@ #include "DataFormatsITSMFT/DPLAlpideParam.h" #include "ITSBase/GeometryTGeo.h" -#include "ITSReconstruction/FastMultEstConfig.h" -#include "ITSReconstruction/FastMultEst.h" +#include "ITStracking/FastMultEstConfig.h" +#include "ITStracking/FastMultEst.h" #include "ITStracking/TrackingConfigParam.h" #include "ITStracking/TrackingInterface.h" diff --git a/Detectors/ITSMFT/ITS/tracking/src/TrackingLinkDef.h b/Detectors/ITSMFT/ITS/tracking/src/TrackingLinkDef.h index 8861f1a255cd8..9efd6dde0176d 100644 --- a/Detectors/ITSMFT/ITS/tracking/src/TrackingLinkDef.h +++ b/Detectors/ITSMFT/ITS/tracking/src/TrackingLinkDef.h @@ -42,4 +42,8 @@ #pragma link C++ class o2::its::ITSGpuTrackingParamConfig + ; #pragma link C++ class o2::conf::ConfigurableParamHelper < o2::its::ITSGpuTrackingParamConfig> + ; +#pragma link C++ class o2::its::FastMultEst + ; +#pragma link C++ class o2::its::FastMultEstConfig + ; +#pragma link C++ class o2::conf::ConfigurableParamHelper < o2::its::FastMultEstConfig> + ; + #endif From dd608d651c9ab2592bfcf1a6d1683b303632fd50 Mon Sep 17 00:00:00 2001 From: Maximiliano Puccio Date: Fri, 27 Mar 2026 09:47:16 +0100 Subject: [PATCH 2/3] Implement new kind of multiplicity mask --- .../GPU/ITStrackingGPU/TimeFrameGPU.h | 4 + .../ITS/tracking/GPU/cuda/TimeFrameGPU.cu | 15 +- .../include/ITStracking/FastMultEst.h | 60 ++-- .../include/ITStracking/FastMultEstConfig.h | 1 + .../include/ITStracking/ROFLookupTables.h | 161 ++++++++++ .../tracking/include/ITStracking/TimeFrame.h | 9 +- .../ITSMFT/ITS/tracking/src/FastMultEst.cxx | 282 +++++++++++------- .../ITSMFT/ITS/tracking/src/TimeFrame.cxx | 14 +- .../ITSMFT/ITS/tracking/src/TrackerTraits.cxx | 19 +- .../ITS/tracking/src/TrackingInterface.cxx | 104 +++---- .../ITS/tracking/test/testROFLookupTables.cxx | 34 +++ .../TRK/reconstruction/src/TimeFrame.cxx | 2 +- 12 files changed, 490 insertions(+), 215 deletions(-) diff --git a/Detectors/ITSMFT/ITS/tracking/GPU/ITStrackingGPU/TimeFrameGPU.h b/Detectors/ITSMFT/ITS/tracking/GPU/ITStrackingGPU/TimeFrameGPU.h index 00bb25e67e354..12be726d11435 100644 --- a/Detectors/ITSMFT/ITS/tracking/GPU/ITStrackingGPU/TimeFrameGPU.h +++ b/Detectors/ITSMFT/ITS/tracking/GPU/ITStrackingGPU/TimeFrameGPU.h @@ -31,6 +31,7 @@ class TimeFrameGPU final : public TimeFrame using typename TimeFrame::IndexTableUtilsN; using typename TimeFrame::ROFOverlapTableN; using typename TimeFrame::ROFVertexLookupTableN; + using typename TimeFrame::ROFMaskTableN; public: TimeFrameGPU() = default; @@ -107,6 +108,7 @@ class TimeFrameGPU final : public TimeFrame IndexTableUtilsN* getDeviceIndexTableUtils() { return mIndexTableUtilsDevice; } const auto getDeviceROFOverlapTableView() { return mDeviceROFOverlapTableView; } const auto getDeviceROFVertexLookupTableView() { return mDeviceROFVertexLookupTableView; } + const auto getDeviceROFMaskTableView() { return mDeviceROFMaskTableView; } int* getDeviceROFramesClusters(const int layer) { return mROFramesClustersDevice[layer]; } auto& getTrackITSExt() { return mTrackITSExt; } Vertex* getDeviceVertices() { return mPrimaryVerticesDevice; } @@ -177,9 +179,11 @@ class TimeFrameGPU final : public TimeFrame // device navigation views ROFOverlapTableN::View mDeviceROFOverlapTableView; ROFVertexLookupTableN::View mDeviceROFVertexLookupTableView; + ROFMaskTableN::View mDeviceROFMaskTableView; // Hybrid pref uint8_t* mMultMaskDevice; + int32_t* mMultMaskOffsetsDevice; Vertex* mPrimaryVerticesDevice; int* mROFramesPVDevice; std::array mClustersDevice; diff --git a/Detectors/ITSMFT/ITS/tracking/GPU/cuda/TimeFrameGPU.cu b/Detectors/ITSMFT/ITS/tracking/GPU/cuda/TimeFrameGPU.cu index 4ff4636963da5..d2e4efc88a61c 100644 --- a/Detectors/ITSMFT/ITS/tracking/GPU/cuda/TimeFrameGPU.cu +++ b/Detectors/ITSMFT/ITS/tracking/GPU/cuda/TimeFrameGPU.cu @@ -252,11 +252,18 @@ void TimeFrameGPU::loadMultiplicityCutMask(const int iteration) { if (!iteration || iteration == 3) { // we need to re-load the swapped mult-mask in upc iteration GPUTimer timer("loading multiplicity cut mask"); - GPULog("gpu-transfer: iteration {} loading multiplicity cut mask with {} elements, for {:.2f} MB.", iteration, this->mMultiplicityCutMask.size(), this->mMultiplicityCutMask.size() * sizeof(uint8_t) / constants::MB); - if (!iteration) { // only allocate on first call - allocMem(reinterpret_cast(&mMultMaskDevice), this->mMultiplicityCutMask.size() * sizeof(uint8_t), this->hasFrameworkAllocator()); + const auto& hostTable = this->mMultiplicityCutMask; + const auto hostView = hostTable.getView(); + GPULog("gpu-transfer: iteration {} loading multiplicity cut mask with {} elements, for {:.2f} MB.", + iteration, hostTable.getFlatMaskSize(), hostTable.getFlatMaskSize() * sizeof(uint8_t) / constants::MB); + if (!iteration) { // only allocate on first call; offsets are stable across iterations + allocMem(reinterpret_cast(&mMultMaskDevice), hostTable.getFlatMaskSize() * sizeof(uint8_t), this->hasFrameworkAllocator()); + allocMem(reinterpret_cast(&mMultMaskOffsetsDevice), NLayers * sizeof(int32_t), this->hasFrameworkAllocator()); + GPUChkErrS(cudaMemcpy(mMultMaskOffsetsDevice, hostView.mLayerROFOffsets, NLayers * sizeof(int32_t), cudaMemcpyHostToDevice)); } - GPUChkErrS(cudaMemcpy(mMultMaskDevice, this->mMultiplicityCutMask.data(), this->mMultiplicityCutMask.size() * sizeof(uint8_t), cudaMemcpyHostToDevice)); + // Re-copy the flat mask on every qualifying iteration (e.g. after swapMasks() for UPC) + GPUChkErrS(cudaMemcpy(mMultMaskDevice, hostView.mFlatMask, hostTable.getFlatMaskSize() * sizeof(uint8_t), cudaMemcpyHostToDevice)); + mDeviceROFMaskTableView = hostTable.getDeviceView(mMultMaskDevice, mMultMaskOffsetsDevice); } } diff --git a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/FastMultEst.h b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/FastMultEst.h index 1ae8713e75a05..46d6c0331f7d0 100644 --- a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/FastMultEst.h +++ b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/FastMultEst.h @@ -17,10 +17,12 @@ #define ALICEO2_ITS_FASTMULTEST_ #include "ITSMFTReconstruction/ChipMappingITS.h" +#include "DataFormatsITS/Vertex.h" #include "DataFormatsITSMFT/ROFRecord.h" #include "DataFormatsITSMFT/CompCluster.h" -#include +#include "DataFormatsITSMFT/PhysTrigger.h" #include "ITStracking/FastMultEstConfig.h" +#include "ITStracking/ROFLookupTables.h" #include #include @@ -32,32 +34,56 @@ namespace its struct FastMultEst { static constexpr int NLayers = o2::itsmft::ChipMappingITS::NLayers; + using ROFOverlapTableN = ROFOverlapTable; + using ROFMaskTableN = ROFMaskTable; - float mult = 0.; /// estimated signal clusters multipliciy at reference (1st?) layer - float noisePerChip = 0.; /// estimated or imposed noise per chip - float cov[3] = {0.}; /// covariance matrix of estimation - float chi2 = 0.; /// chi2 - int nLayersUsed = 0; /// number of layers actually used + float mult = 0.; /// estimated signal clusters multiplicity on the selected multiplicity layer + float noisePerChip = 0.; /// imposed noise per chip (when enabled by configuration) + float cov[3] = {0.}; /// retained for compatibility; set to zero in single-layer mode + float chi2 = 0.; /// retained for compatibility; set to zero in single-layer mode + int nLayersUsed = 0; /// number of layers used by estimator (0/1 in single-layer mode) uint32_t lastRandomSeed = 0; /// state of the gRandom before - - std::array nClPerLayer{0}; // measured N Cl per layer selectROFs FastMultEst(); static uint32_t getCurrentRandomSeed(); - int selectROFs(const gsl::span rofs, const gsl::span clus, - const gsl::span trig, std::vector& sel); + int selectROFs(const std::array, NLayers>& rofs, + const std::array, NLayers>& clus, + const gsl::span trig, + uint32_t firstTForbit, + bool doStaggering, + const ROFOverlapTableN::View& overlapView, + ROFMaskTableN& sel); + void selectROFsWithVertices(const auto& vertices, const ROFOverlapTableN::View& overlapView, ROFMaskTableN& sel) + { + const auto& multEstConf = FastMultEstConfig::Instance(); + if (!multEstConf.isVtxMultCutRequested()) { + return; + } + + for (const auto& vertex : vertices) { + 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); + } + } + } + } + } - void fillNClPerLayer(const gsl::span& clusters); - float process(const std::array ncl) + int countClustersOnLayer(const gsl::span& clusters) const; + float process(int nClusters) { - return FastMultEstConfig::Instance().imposeNoisePerChip > 0 ? processNoiseImposed(ncl) : processNoiseFree(ncl); + return FastMultEstConfig::Instance().imposeNoisePerChip > 0 ? processNoiseImposed(nClusters) : processNoiseFree(nClusters); } - float processNoiseFree(const std::array ncl); - float processNoiseImposed(const std::array ncl); + float processNoiseFree(int nClusters); + float processNoiseImposed(int nClusters); float process(const gsl::span& clusters) { - fillNClPerLayer(clusters); - return process(nClPerLayer); + return process(countClustersOnLayer(clusters)); } static bool sSeedSet; diff --git a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/FastMultEstConfig.h b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/FastMultEstConfig.h index c6bce50995a4b..7d30bcba45844 100644 --- a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/FastMultEstConfig.h +++ b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/FastMultEstConfig.h @@ -34,6 +34,7 @@ struct FastMultEstConfig : public o2::conf::ConfigurableParamHelper std::vector mFlatTable; }; +// GPU-friendly view of the ROF mask table +template +struct ROFMaskTableView { + const uint8_t* mFlatMask{nullptr}; + const int32_t* mLayerROFOffsets{nullptr}; // size NLayers+1 + + GPUhdi() bool isROFEnabled(int32_t layer, int32_t rofId) const noexcept + { + assert(layer >= 0 && layer < NLayers); + return mFlatMask[mLayerROFOffsets[layer] + rofId] != 0; + } + +#ifndef GPUCA_GPUCODE + GPUh() void printAll() const + { + for (int32_t i = 0; i < NLayers; ++i) { + printLayer(i); + } + } + + GPUh() void printLayer(int32_t layer) const + { + constexpr int w_rof = 10; + constexpr int w_active = 10; + int32_t nROFs = mLayerROFOffsets[layer + 1] - mLayerROFOffsets[layer]; + LOGF(info, "Mask table: Layer %d", layer); + LOGF(info, "%*s | %*s", w_rof, "ROF", w_active, "Enabled"); + LOGF(info, "%.*s-+-%.*s", w_rof, "----------", w_active, "----------"); + for (int32_t i = 0; i < nROFs; ++i) { + LOGF(info, "%*d | %*d", w_rof, i, w_active, (int)isROFEnabled(layer, i)); + } + } +#endif +}; + +// Per-ROF per-layer boolean mask (uint8_t for GPU compatibility). +template +class ROFMaskTable : public LayerTimingBase +{ + public: + using BCRange = dataformats::RangeReference; + using View = ROFMaskTableView; + + GPUdDefault() ROFMaskTable() = default; + GPUh() explicit ROFMaskTable(const LayerTimingBase& timingBase) : LayerTimingBase(timingBase) { init(); } + + GPUh() void init() + { + int32_t totalROFs = 0; + for (int32_t layer{0}; layer < NLayers; ++layer) { + mLayerROFOffsets[layer] = totalROFs; + totalROFs += this->getLayer(layer).mNROFsTF; + } + mLayerROFOffsets[NLayers] = totalROFs; // sentinel + mFlatMask.resize(totalROFs, 1); + } + + GPUh() size_t getFlatMaskSize() const noexcept { return mFlatMask.size(); } + + GPUh() bool isROFEnabled(int32_t layer, int32_t rofId) const noexcept { return mFlatMask[mLayerROFOffsets[layer] + rofId] != 0; } + + GPUh() void setROFEnabled(int32_t layer, int32_t rofId, uint8_t state = 1) noexcept + { + assert(layer >= 0 && layer < NLayers); + assert(rofId >= 0 && rofId < mLayerROFOffsets[layer + 1] - mLayerROFOffsets[layer]); + mFlatMask[mLayerROFOffsets[layer] + rofId] = state; + } + + GPUh() void setROFsEnabled(int32_t layer, int32_t firstRof, int32_t nRofs, uint8_t state = 1) noexcept + { + assert(layer >= 0 && layer < NLayers); + assert(firstRof >= 0); + assert(firstRof + nRofs <= mLayerROFOffsets[layer + 1] - mLayerROFOffsets[layer]); + std::memset(mFlatMask.data() + mLayerROFOffsets[layer] + firstRof, state, nRofs); + } + + // Enable all ROFs in all layers that are time-compatible with the given BC range + GPUh() void selectROF(const BCRange& t) + { + const int32_t bcStart = t.getFirstEntry(); + const int32_t bcEnd = t.getEntriesBound(); + for (int32_t layer{0}; layer < NLayers; ++layer) { + const auto& lay = this->getLayer(layer); + const int32_t offset = mLayerROFOffsets[layer]; + for (int32_t rofId{0}; rofId < lay.mNROFsTF; ++rofId) { + if (static_cast(lay.getROFStartInBC(rofId)) < bcEnd && + static_cast(lay.getROFEndInBC(rofId)) > bcStart) { + mFlatMask[offset + rofId] = 1; + } + } + } + } + + // Reset mask to 0, then enable all ROFs compatible with any of the given BC ranges + GPUh() void selectROFs(const std::vector& ts) + { + resetMask(); + for (const auto& t : ts) { + selectROF(t); + } + } + + GPUh() void resetMask(uint8_t s = 0) + { + std::memset(mFlatMask.data(), s, mFlatMask.size()); + } + + GPUh() void invertMask() + { + std::ranges::transform(mFlatMask, mFlatMask.begin(), [](uint8_t x) { return 1 - x; }); + } + + GPUh() void swap(ROFMaskTable& other) noexcept + { + std::swap(mFlatMask, other.mFlatMask); + std::swap(mLayerROFOffsets, other.mLayerROFOffsets); + } + + GPUh() View getView() const + { + View view; + view.mFlatMask = mFlatMask.data(); + view.mLayerROFOffsets = mLayerROFOffsets; + return view; + } + + GPUh() View getDeviceView(const uint8_t* deviceFlatMaskPtr, const int32_t* deviceOffsetPtr) const + { + View view; + view.mFlatMask = deviceFlatMaskPtr; + view.mLayerROFOffsets = deviceOffsetPtr; + return view; + } +#ifndef GPUCA_GPUCODE + GPUh() std::string asString() const + { + std::string mask_str; + for (int32_t i = 0; i < NLayers; ++i) { + int32_t nROFs = mLayerROFOffsets[i + 1] - mLayerROFOffsets[i]; + int32_t enabledROFs = 0; + for (int32_t j = 0; j < nROFs; ++j) { + if (isROFEnabled(i, j)) { + ++enabledROFs; + } + } + mask_str += std::format("Layer {} ROFs enabled: {}/{} | ", i, enabledROFs, nROFs); + } + return mask_str; + } + + GPUh() void print() const + { + LOG(info) << asString(); + } +#endif + + private: + int32_t mLayerROFOffsets[NLayers + 1]{}; // NLayers entries + 1 sentinel + std::vector mFlatMask; +}; + } // namespace o2::its #endif diff --git a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/TimeFrame.h b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/TimeFrame.h index cc0b73f1eac76..c8058319638d9 100644 --- a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/TimeFrame.h +++ b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/TimeFrame.h @@ -67,6 +67,7 @@ struct TimeFrame { using IndexTableUtilsN = IndexTableUtils; using ROFOverlapTableN = ROFOverlapTable; using ROFVertexLookupTableN = ROFVertexLookupTable; + using ROFMaskTableN = ROFMaskTable; using CellSeedN = CellSeed; friend class gpu::TimeFrameGPU; @@ -220,8 +221,8 @@ struct TimeFrame { std::array& getBeamXY() { return mBeamPos; } // \Vertexer - void setMultiplicityCutMask(const std::vector& cutMask) { mMultiplicityCutMask = cutMask; } - void setROFMask(const std::vector& rofMask) { mROFMask = rofMask; } + void setMultiplicityCutMask(const ROFMaskTableN& cutMask) { mMultiplicityCutMask = cutMask; } + void setROFMask(const ROFMaskTableN& rofMask) { mROFMask = rofMask; } void swapMasks() { mMultiplicityCutMask.swap(mROFMask); } int hasBogusClusters() const { return std::accumulate(mBogusClusters.begin(), mBogusClusters.end(), 0); } @@ -266,7 +267,7 @@ struct TimeFrame { bounded_vector mTracksLabel; std::vector> mCellsNeighbours; std::vector> mCellsLookupTable; - std::vector mMultiplicityCutMask; + ROFMaskTableN mMultiplicityCutMask; const o2::base::PropagatorImpl* mPropagatorDevice = nullptr; // Needed only for GPU @@ -290,7 +291,7 @@ struct TimeFrame { bounded_vector mPositionResolution; std::array, NLayers> mClusterSize; - std::vector mROFMask; + ROFMaskTableN mROFMask; bounded_vector> mPValphaX; /// PV x and alpha for track propagation std::vector> mTrackletLabels; std::vector> mCellLabels; diff --git a/Detectors/ITSMFT/ITS/tracking/src/FastMultEst.cxx b/Detectors/ITSMFT/ITS/tracking/src/FastMultEst.cxx index 4b854be552f80..0547209129b70 100644 --- a/Detectors/ITSMFT/ITS/tracking/src/FastMultEst.cxx +++ b/Detectors/ITSMFT/ITS/tracking/src/FastMultEst.cxx @@ -14,14 +14,76 @@ /// \author ruben.shahoyan@cern.ch #include "ITStracking/FastMultEst.h" -#include "DataFormatsITSMFT/DPLAlpideParam.h" #include "Framework/Logger.h" #include #include +#include #include using namespace o2::its; +namespace +{ + +// Convert trigger IR to ROF index on a given layer using LayerTiming +int findROFForIR(const o2::InteractionRecord& ir, + const o2::InteractionRecord& tfStartIR, + const LayerTiming& layerTiming) +{ + // Convert IR to BC-from-TF-start, which is the time base expected by LayerTiming. + const int64_t bcFromTFStart = ir.differenceInBC(tfStartIR); + if (bcFromTFStart < 0) { + return -1; + } + return layerTiming.getROF(static_cast(bcFromTFStart)); +} + +template +void enableCompatibleROFs(int baseLayer, + int baseRof, + const typename o2::its::ROFOverlapTable::View& overlapView, + o2::its::ROFMaskTable& sel) +{ + sel.setROFEnabled(baseLayer, baseRof); + for (int layer = 0; layer < NLayers; ++layer) { + if (layer == baseLayer) { + continue; + } + const auto& overlap = overlapView.getOverlap(baseLayer, layer, baseRof); + if (overlap.getEntries() > 0) { + sel.setROFsEnabled(layer, overlap.getFirstEntry(), overlap.getEntries()); + } + } +} + +template +std::vector buildMultiplicityCounts(const std::array, NLayers>& rofs, + const std::array, NLayers>& clus, + bool doStaggering, + int multLayer) +{ + std::vector multCounts; + if (doStaggering) { + multCounts.resize(rofs[multLayer].size()); + for (size_t iRof = 0; iRof < rofs[multLayer].size(); ++iRof) { + multCounts[iRof] = rofs[multLayer][iRof].getNEntries(); + } + return multCounts; + } + + static const o2::itsmft::ChipMappingITS chipMapping; + multCounts.resize(rofs[0].size(), 0); + for (size_t iRof = 0; iRof < rofs[0].size(); ++iRof) { + for (const auto& cluster : rofs[0][iRof].getROFData(clus[0])) { + if (chipMapping.getLayer(cluster.getSensorID()) == multLayer) { + ++multCounts[iRof]; + } + } + } + return multCounts; +} +} // namespace + bool FastMultEst::sSeedSet = false; ///______________________________________________________ @@ -38,152 +100,152 @@ FastMultEst::FastMultEst() } ///______________________________________________________ -/// find multiplicity for given set of clusters -void FastMultEst::fillNClPerLayer(const gsl::span& clusters) +/// count clusters on the configured multiplicity layer +int FastMultEst::countClustersOnLayer(const gsl::span& clusters) const { - int lr = FastMultEst::NLayers - 1, nchAcc = o2::itsmft::ChipMappingITS::getNChips() - o2::itsmft::ChipMappingITS::getNChipsPerLr(lr); - std::memset(&nClPerLayer[0], 0, sizeof(int) * FastMultEst::NLayers); + const int targetLayer = std::clamp(FastMultEstConfig::Instance().cutMultClusLayer, 0, NLayers - 1); + int count = 0; + int lr = FastMultEst::NLayers - 1; + int nchAcc = o2::itsmft::ChipMappingITS::getNChips() - o2::itsmft::ChipMappingITS::getNChipsPerLr(lr); for (int i = clusters.size(); i--;) { // profit from clusters being ordered in chip increasing order while (clusters[i].getSensorID() < nchAcc) { assert(lr >= 0); nchAcc -= o2::itsmft::ChipMappingITS::getNChipsPerLr(--lr); } - nClPerLayer[lr]++; + if (lr == targetLayer) { + ++count; + } } + return count; } ///______________________________________________________ /// find multiplicity for given number of clusters per layer -float FastMultEst::processNoiseFree(const std::array ncl) +float FastMultEst::processNoiseFree(int nClusters) { - // we assume that on the used layers the observed number of clusters is defined by the - // the noise ~ nu * Nchips and contribution from the signal tracks Ntr*mAccCorr + // Single-layer regime: estimate multiplicity from one configured layer only. const auto& conf = FastMultEstConfig::Instance(); - - float mat[3] = {0}, b[2] = {0}; - nLayersUsed = 0; - for (int il = conf.firstLayer; il <= conf.lastLayer; il++) { - if (ncl[il] > 0) { - int nch = o2::itsmft::ChipMappingITS::getNChipsPerLr(il); - float err2i = 1. / ncl[il]; - float m2n = nch * err2i; - mat[0] += err2i * conf.accCorr[il] * conf.accCorr[il]; - mat[2] += nch * m2n; - mat[1] += conf.accCorr[il] * m2n; // non-diagonal element - b[0] += conf.accCorr[il]; - b[1] += nch; - nLayersUsed++; - } - } - mult = noisePerChip = chi2 = -1; - float det = mat[0] * mat[2] - mat[1] * mat[1]; - if (nLayersUsed < 2 || std::abs(det) < 1e-15) { - return -1; + const int layer = std::clamp(conf.cutMultClusLayer, 0, NLayers - 1); + const float acc = conf.accCorr[layer]; + nLayersUsed = nClusters > 0 ? 1 : 0; + noisePerChip = 0.f; + chi2 = 0.f; + cov[0] = cov[1] = cov[2] = 0.f; + if (nLayersUsed == 0 || acc <= 0.f) { + mult = -1.f; + return -1.f; } - float detI = 1. / det; - mult = detI * (b[0] * mat[2] - b[1] * mat[1]); - noisePerChip = detI * (b[1] * mat[0] - b[0] * mat[1]); - cov[0] = mat[2] * detI; - cov[2] = mat[0] * detI; - cov[1] = -mat[1] * detI; - chi2 = 0.; - for (int il = conf.firstLayer; il <= conf.lastLayer; il++) { - if (ncl[il] > 0) { - int nch = o2::itsmft::ChipMappingITS::getNChipsPerLr(il); - float diff = mult * conf.accCorr[il] + nch * noisePerChip - ncl[il]; - chi2 += diff * diff / ncl[il]; - } - } - chi2 = nLayersUsed > 2 ? chi2 / (nLayersUsed - 2) : 0.; + mult = nClusters / acc; return mult > 0 ? mult : 0; } ///______________________________________________________ /// find multiplicity for given number of clusters per layer with mean noise imposed -float FastMultEst::processNoiseImposed(const std::array ncl) +float FastMultEst::processNoiseImposed(int nClusters) { - // we assume that on the used layers the observed number of clusters is defined by the - // the noise ~ nu * Nchips and contribution from the signal tracks Ntr*conf.accCorr - // + // Single-layer regime with imposed noise subtraction. const auto& conf = FastMultEstConfig::Instance(); - float accSum = 0., accWSum = 0., noiseSum = 0.; - nLayersUsed = 0; - for (int il = conf.firstLayer; il <= conf.lastLayer; il++) { - if (ncl[il] > 0) { - float err = 1. / ncl[il]; - accSum += conf.accCorr[il]; - accWSum += conf.accCorr[il] * conf.accCorr[il] * err; - noiseSum += o2::itsmft::ChipMappingITS::getNChipsPerLr(il) * conf.accCorr[il] * err; - nLayersUsed++; - } - } - mult = 0; - if (nLayersUsed) { - mult = (accSum - noisePerChip * noiseSum) / accWSum; + const int layer = std::clamp(conf.cutMultClusLayer, 0, NLayers - 1); + const float acc = conf.accCorr[layer]; + const float nch = static_cast(o2::itsmft::ChipMappingITS::getNChipsPerLr(layer)); + nLayersUsed = nClusters > 0 ? 1 : 0; + chi2 = 0.f; + cov[0] = cov[1] = cov[2] = 0.f; + if (nLayersUsed == 0 || acc <= 0.f) { + mult = -1.f; + return -1.f; } + mult = (nClusters - noisePerChip * nch) / acc; return mult; } -int FastMultEst::selectROFs(const gsl::span rofs, const gsl::span clus, - const gsl::span trig, std::vector& sel) +int FastMultEst::selectROFs(const std::array, NLayers>& rofs, + const std::array, NLayers>& clus, + const gsl::span trig, + uint32_t firstTForbit, + bool doStaggering, + const ROFOverlapTableN::View& overlapView, + ROFMaskTableN& sel) { - int nrof = rofs.size(), nsel = 0; const auto& multEstConf = FastMultEstConfig::Instance(); // parameters for mult estimation and cuts - sel.clear(); - sel.resize(nrof, true); // by default select all + const int selectionLayer = overlapView.getClock(); + int multLayer = std::clamp(multEstConf.cutMultClusLayer, 0, NLayers - 1); + if (doStaggering && rofs[multLayer].empty()) { + LOGP(warning, "FastMultEst multiplicity layer {} has no ROFs, falling back to selection layer {}", multLayer, selectionLayer); + multLayer = selectionLayer; + } + + const auto multCounts = buildMultiplicityCounts(rofs, clus, doStaggering, multLayer); + const int selectionRofCount = doStaggering ? static_cast(rofs[selectionLayer].size()) : static_cast(rofs[0].size()); + + sel.resetMask(); lastRandomSeed = gRandom->GetSeed(); - if (multEstConf.isMultCutRequested()) { - for (uint32_t irof = 0; irof < nrof; irof++) { - nsel += sel[irof] = multEstConf.isPassingMultCut(process(rofs[irof].getROFData(clus))); + const o2::InteractionRecord tfStartIR{0, firstTForbit}; + + if (!trig.empty()) { + const auto& selectionLayerTiming = overlapView.getLayer(selectionLayer); + const auto& multLayerTiming = overlapView.getLayer(multLayer); + + for (const auto& trigger : trig) { + const int selectionRof = findROFForIR(trigger.ir, tfStartIR, selectionLayerTiming); + if (selectionRof < 0) { + continue; + } + if (multEstConf.cutRandomFraction > 0.f && gRandom->Rndm() < multEstConf.cutRandomFraction) { + continue; + } + if (multEstConf.isMultCutRequested()) { + const int triggerMultRof = doStaggering ? findROFForIR(trigger.ir, tfStartIR, multLayerTiming) : selectionRof; + if (triggerMultRof < 0 || triggerMultRof >= static_cast(multCounts.size())) { + continue; + } + if (!multEstConf.isPassingMultCut(process(multCounts[triggerMultRof]))) { + continue; + } + } + enableCompatibleROFs(selectionLayer, selectionRof, overlapView, sel); } } else { - nsel = nrof; - } - using IdNT = std::pair; - if (multEstConf.cutRandomFraction > 0.) { - int ntrig = trig.size(), currTrig = 0; - if (multEstConf.preferTriggered) { - const auto& alpParams = o2::itsmft::DPLAlpideParam::Instance(); - std::vector nTrigROF; - nTrigROF.reserve(nrof); - for (uint32_t irof = 0; irof < nrof; irof++) { - if (sel[irof]) { - if (nsel && gRandom->Rndm() < multEstConf.cutRandomFraction) { - nsel--; + LOGP(warning, "FastMultEst received no physics/TRD triggers, falling back to ROF-driven filtering on layer {}", selectionLayer); + for (int selectionRof = 0; selectionRof < selectionRofCount; ++selectionRof) { + if (multEstConf.isMultCutRequested()) { + bool passes = false; + if (!doStaggering || selectionLayer == multLayer) { + if (selectionRof < static_cast(multCounts.size())) { + passes = multEstConf.isPassingMultCut(process(multCounts[selectionRof])); } - auto irROF = rofs[irof].getBCData(); - while (currTrig < ntrig && trig[currTrig].ir < irROF) { // triggers are sorted, jump to 1st one not less than current ROF - currTrig++; - } - auto& trof = nTrigROF.emplace_back(irof, 0); - irROF += alpParams.roFrameLengthInBC; - while (currTrig < ntrig && trig[currTrig].ir < irROF) { - trof.second++; - currTrig++; + } else { + const auto& overlap = overlapView.getOverlap(selectionLayer, multLayer, selectionRof); + for (int rof = overlap.getFirstEntry(); rof < overlap.getEntriesBound(); ++rof) { + if (rof < static_cast(multCounts.size())) { + if (multEstConf.isPassingMultCut(process(multCounts[rof]))) { + passes = true; + break; + } + } } } - } - if (nsel > 0) { - sort(nTrigROF.begin(), nTrigROF.end(), [](const IdNT& a, const IdNT& b) { return a.second > b.second; }); // order in number of triggers - auto last = nTrigROF.begin() + nsel; - sort(nTrigROF.begin(), last, [](const IdNT& a, const IdNT& b) { return a.first < b.first; }); // order in ROF ID first nsel ROFs - } - for (int i = nsel; i < int(nTrigROF.size()); i++) { // reject ROFs in the tail - sel[nTrigROF[i].first] = false; - } - } else { // dummy random rejection - for (int irof = 0; irof < nrof; irof++) { - if (sel[irof]) { - float sr = gRandom->Rndm(); - if (gRandom->Rndm() < multEstConf.cutRandomFraction) { - sel[irof] = false; - nsel--; - } + if (!passes) { + continue; } } + if (multEstConf.cutRandomFraction > 0.f && gRandom->Rndm() < multEstConf.cutRandomFraction) { + continue; + } + enableCompatibleROFs(selectionLayer, selectionRof, overlapView, sel); } } - LOGP(debug, "NSel = {} of {} rofs Seeds: before {} after {}", nsel, nrof, lastRandomSeed, gRandom->GetSeed()); + + int nsel = 0; + for (int irof = 0; irof < selectionRofCount; ++irof) { + nsel += sel.isROFEnabled(selectionLayer, irof); + } + + if (!trig.empty() && multEstConf.preferTriggered) { + LOGP(debug, "FastMultEst preferTriggered is ignored in trigger-driven mask mode"); + } + + LOGP(debug, "NSel = {} of {} rofs on layer {} Seeds: before {} after {}", nsel, selectionRofCount, selectionLayer, lastRandomSeed, gRandom->GetSeed()); return nsel; } diff --git a/Detectors/ITSMFT/ITS/tracking/src/TimeFrame.cxx b/Detectors/ITSMFT/ITS/tracking/src/TimeFrame.cxx index e7e3290d2d7b4..8f5ae0cee25e5 100644 --- a/Detectors/ITSMFT/ITS/tracking/src/TimeFrame.cxx +++ b/Detectors/ITSMFT/ITS/tracking/src/TimeFrame.cxx @@ -191,10 +191,9 @@ void TimeFrame::prepareClusters(const TrackingParameters& trkParam, con bounded_vector lutPerBin(numBins, 0, mMemoryPool.get()); for (int iLayer{0}, stopLayer = std::min(trkParam.NLayers, maxLayers); iLayer < stopLayer; ++iLayer) { for (int rof{0}; rof < getNrof(iLayer); ++rof) { - // FIXME - // if (!iLayer && !mMultiplicityCutMask[rof]) { - // continue; - // } + if (!mMultiplicityCutMask.isROFEnabled(iLayer, rof)) { + continue; + } const auto& unsortedClusters{getUnsortedClustersOnLayer(rof, iLayer)}; const int clustersNum{static_cast(unsortedClusters.size())}; auto* tableBase = mIndexTables[iLayer].data() + rof * stride; @@ -374,10 +373,9 @@ void TimeFrame::computeTrackletsPerROFScans() { for (ushort iLayer = 0; iLayer < 2; ++iLayer) { for (unsigned int iRof{0}; iRof < getNrof(1); ++iRof) { - // FIXME? - // if (mMultiplicityCutMask[iRof]) { - mTotalTracklets[iLayer] += mNTrackletsPerROF[iLayer][iRof]; - // } + if (mMultiplicityCutMask.isROFEnabled(iLayer, iRof)) { + mTotalTracklets[iLayer] += mNTrackletsPerROF[iLayer][iRof]; + } } std::exclusive_scan(mNTrackletsPerROF[iLayer].begin(), mNTrackletsPerROF[iLayer].end(), mNTrackletsPerROF[iLayer].begin(), 0); std::exclusive_scan(mNTrackletsPerCluster[iLayer].begin(), mNTrackletsPerCluster[iLayer].end(), mNTrackletsPerClusterSum[iLayer].begin(), 0); diff --git a/Detectors/ITSMFT/ITS/tracking/src/TrackerTraits.cxx b/Detectors/ITSMFT/ITS/tracking/src/TrackerTraits.cxx index a1db58027ca00..19671282b976d 100644 --- a/Detectors/ITSMFT/ITS/tracking/src/TrackerTraits.cxx +++ b/Detectors/ITSMFT/ITS/tracking/src/TrackerTraits.cxx @@ -24,11 +24,12 @@ #include "CommonConstants/MathConstants.h" #include "DetectorsBase/Propagator.h" #include "GPUCommonMath.h" +#include "ITStracking/BoundedAllocator.h" #include "ITStracking/Cell.h" #include "ITStracking/Constants.h" -#include "ITStracking/TrackerTraits.h" -#include "ITStracking/BoundedAllocator.h" #include "ITStracking/IndexTableUtils.h" +#include "ITStracking/ROFLookupTables.h" +#include "ITStracking/TrackerTraits.h" #include "ITStracking/Tracklet.h" #include "ReconstructionDataFormats/Track.h" @@ -59,10 +60,9 @@ void TrackerTraits::computeLayerTracklets(const int iteration, int iVer mTaskArena->execute([&] { auto forTracklets = [&](auto Tag, int iLayer, int pivotROF, int base, int& offset) -> int { - // FIXME - // if (!mTimeFrame->mMultiplicityCutMask[pivotROF]) { - // return 0; - // } + if (!mTimeFrame->mMultiplicityCutMask.isROFEnabled(iLayer, pivotROF)) { + return 0; + } gsl::span primaryVertices = mTrkParams[iteration].UseDiamond ? diamondSpan : mTimeFrame->getPrimaryVertices(iLayer, pivotROF); if (primaryVertices.empty()) { return 0; @@ -121,10 +121,9 @@ void TrackerTraits::computeLayerTracklets(const int iteration, int iVer } for (int targetROF = rofOverlap.getFirstEntry(); targetROF < rofOverlap.getEntriesBound(); ++targetROF) { - // FIXME - // if (!mTimeFrame->mMultiplicityCutMask[targetROF]) { - // continue; - // } + if (!mTimeFrame->mMultiplicityCutMask.isROFEnabled(iLayer + 1, targetROF)) { + continue; + } auto layer1 = mTimeFrame->getClustersOnLayer(targetROF, iLayer + 1); if (layer1.empty()) { continue; diff --git a/Detectors/ITSMFT/ITS/tracking/src/TrackingInterface.cxx b/Detectors/ITSMFT/ITS/tracking/src/TrackingInterface.cxx index 2d46de523e493..3c153d8ff6895 100644 --- a/Detectors/ITSMFT/ITS/tracking/src/TrackingInterface.cxx +++ b/Detectors/ITSMFT/ITS/tracking/src/TrackingInterface.cxx @@ -20,6 +20,7 @@ #include "ITStracking/FastMultEstConfig.h" #include "ITStracking/FastMultEst.h" +#include "ITStracking/ROFLookupTables.h" #include "ITStracking/TrackingConfigParam.h" #include "ITStracking/TrackingInterface.h" @@ -176,12 +177,13 @@ void ITSTrackingInterface::run(framework::ProcessingContext& pc) auto logger = [&](const std::string& s) { LOG(info) << s; }; auto fatalLogger = [&](const std::string& s) { LOG(fatal) << s; }; auto errorLogger = [&](const std::string& s) { LOG(error) << s; }; + const auto firstTForbit = pc.services().get().firstTForbit; FastMultEst multEst; // mult estimator - std::vector processingMask, processUPCMask; - // int cutVertexMult{0}, cutUPCVertex{0}, cutRandomMult = int(trackROFvec.size()) - multEst.selectROFs(trackROFvec, compClusters, physTriggers, processingMask); - // processUPCMask.resize(processingMask.size(), false); - // mTimeFrame->setMultiplicityCutMask(processingMask); + o2::its::ROFMaskTable processingMask{mTimeFrame->getROFOverlapTable()}, processUPCMask{mTimeFrame->getROFOverlapTable()}; + multEst.selectROFs(rofsinput, compClusters, physTriggers, firstTForbit, mDoStaggering, mTimeFrame->getROFOverlapTableView(), processingMask); + mTimeFrame->setMultiplicityCutMask(processingMask); + logger(processingMask.asString()); float vertexerElapsedTime{0.f}; if (mRunVertexer) { // Run seeding vertexer @@ -192,63 +194,43 @@ void ITSTrackingInterface::run(framework::ProcessingContext& pc) vertices.insert(vertices.begin(), vtx.begin(), vtx.end()); } } - // const auto& multEstConf = FastMultEstConfig::Instance(); // parameters for mult estimation and cuts - // gsl::span vMCRecInfo; - // gsl::span vMCContLabels; - // for (auto iRof{0}; iRof < trackROFspan.size(); ++iRof) { - // bounded_vector vtxVecLoc; - // auto& vtxROF = vertROFvec.emplace_back(trackROFspan[iRof]); - // vtxROF.setFirstEntry(vertices.size()); - // if (mRunVertexer) { - // auto vtxSpan = mTimeFrame->getPrimaryVertices(iRof); - // if (mIsMC) { - // vMCRecInfo = mTimeFrame->getPrimaryVerticesMCRecInfo(iRof); - // } - // if (o2::its::TrackerParamConfig::Instance().doUPCIteration) { - // if (!vtxSpan.empty()) { - // if (vtxSpan[0].isFlagSet(Vertex::UPCMode) == 1) { // at least one vertex in this ROF and it is from second vertex iteration - // LOGP(debug, "ROF {} rejected as vertices are from the UPC iteration", iRof); - // processUPCMask[iRof] = true; - // cutUPCVertex++; - // vtxROF.setFlag(o2::itsmft::ROFRecord::VtxUPCMode); - // } else { // in all cases except if as standard mode vertex was found, the ROF was processed with UPC settings - // vtxROF.setFlag(o2::itsmft::ROFRecord::VtxStdMode); - // } - // } else { - // vtxROF.setFlag(o2::itsmft::ROFRecord::VtxUPCMode); - // } - // } else { - // vtxROF.setFlag(o2::itsmft::ROFRecord::VtxStdMode); - // } - // vtxROF.setNEntries(vtxSpan.size()); - // bool selROF = vtxSpan.empty(); - // for (int iV{0}, iVC{0}; iV < vtxSpan.size(); ++iV) { - // const auto& v = vtxSpan[iV]; - // if (multEstConf.isVtxMultCutRequested() && !multEstConf.isPassingVtxMultCut(v.getNContributors())) { - // iVC += v.getNContributors(); - // continue; // skip vertex of unwanted multiplicity - // } - // selROF = true; - // vertices.push_back(v); - // if (mIsMC && !VertexerParamConfig::Instance().useTruthSeeding) { - // allVerticesLabels.push_back(vMCRecInfo[iV].first); - // allVerticesPurities.push_back(vMCRecInfo[iV].second); - // } - // iVC += v.getNContributors(); - // } - // if (processingMask[iRof] && !selROF) { // passed selection in clusters and not in vertex multiplicity - // LOGP(info, "ROF {} rejected by the vertex multiplicity selection [{},{}]", iRof, multEstConf.cutMultVtxLow, multEstConf.cutMultVtxHigh); - // processingMask[iRof] = selROF; - // cutVertexMult++; - // } - // } - // } + multEst.selectROFsWithVertices(vertices, mTimeFrame->getROFOverlapTableView(), processingMask); + + const auto& multEstConf = FastMultEstConfig::Instance(); // parameters for mult estimation and cuts + gsl::span vMCRecInfo; + gsl::span vMCContLabels; + const int clockLayerId{mDoStaggering ? mTimeFrame->getROFOverlapTableView().getClock() : 0}; + auto clockROFspan = rofsinput[clockLayerId]; + auto clockTiming = mTimeFrame->getROFOverlapTableView().getClockLayer(); + for (auto iRof{0}; iRof < clockROFspan.size(); ++iRof) { + bounded_vector vtxVecLoc; + auto& vtxROF = vertROFvec.emplace_back(clockROFspan[iRof]); + vtxROF.setFirstEntry(vertices.size()); + + if (mRunVertexer) { + auto vtxSpan = mTimeFrame->getPrimaryVertices(clockLayerId, iRof); + if (o2::its::TrackerParamConfig::Instance().doUPCIteration) { + if (!vtxSpan.empty()) { + if (vtxSpan[0].isFlagSet(Vertex::UPCMode) == 1) { // at least one vertex in this ROF and it is from second vertex iteration + LOGP(debug, "ROF {} rejected as vertices are from the UPC iteration", iRof); + processUPCMask.selectROF({clockTiming.getROFStartInBC(iRof), clockTiming.getROFEndInBC(iRof)}); + vtxROF.setFlag(o2::itsmft::ROFRecord::VtxUPCMode); + } else { // in all cases except if as standard mode vertex was found, the ROF was processed with UPC settings + vtxROF.setFlag(o2::itsmft::ROFRecord::VtxStdMode); + } + } else { + vtxROF.setFlag(o2::itsmft::ROFRecord::VtxUPCMode); + } + } else { + vtxROF.setFlag(o2::itsmft::ROFRecord::VtxStdMode); + } + vtxROF.setNEntries(vtxSpan.size()); + } + } if (mRunVertexer && !compClusters.empty()) { LOG(info) << fmt::format(" - Vertex seeding total elapsed time: {} ms for {} vertices found", vertexerElapsedTime, mTimeFrame->getPrimaryVerticesNum()); - // FIXME - // LOG(info) << fmt::format(" - FastMultEst: rejected {}/{} ROFs: random/mult.sel:{} (seed {}), vtx.sel:{}", cutRandomMult + cutVertexMult, trackROFspan.size(), cutRandomMult, multEst.lastRandomSeed, cutVertexMult); } if (mOverrideBeamEstimation) { LOG(info) << fmt::format(" - Beam position set to: {}, {} from meanvertex object", mTimeFrame->getBeamX(), mTimeFrame->getBeamY()); @@ -274,10 +256,6 @@ void ITSTrackingInterface::run(framework::ProcessingContext& pc) LOG(warning) << fmt::format(" - The processed timeframe had {} clusters with wild z coordinates, check the dictionaries", mTimeFrame->hasBogusClusters()); } - // FIXME - // if (processingMask[iROF]) { - // irFrames.emplace_back(tracksROF.getBCData(), tracksROF.getBCData() + nBCPerTF - 1).info = tracks.size(); - // } auto& tracks = mTimeFrame->getTracks(); allTrackLabels.reserve(mTimeFrame->getTracksLabel().size()); // should be 0 if not MC std::copy(mTimeFrame->getTracksLabel().begin(), mTimeFrame->getTracksLabel().end(), std::back_inserter(allTrackLabels)); @@ -340,6 +318,10 @@ void ITSTrackingInterface::run(framework::ProcessingContext& pc) for (size_t iROF{0}; iROF < allTrackROFs.size(); ++iROF) { allTrackROFs[iROF].setFirstEntry(rofEntries[iROF]); allTrackROFs[iROF].setNEntries(rofEntries[iROF + 1] - rofEntries[iROF]); + if (processingMask.isROFEnabled(clockLayerId, iROF)) { + auto& irFrame = irFrames.emplace_back(allTrackROFs[iROF].getBCData(), allTrackROFs[iROF].getBCData() + clockLayer.mROFLength - 1); + irFrame.info = allTrackROFs[iROF].getNEntries(); + } } // same thing for vertices rofs std::fill(rofEntries.begin(), rofEntries.end(), 0); diff --git a/Detectors/ITSMFT/ITS/tracking/test/testROFLookupTables.cxx b/Detectors/ITSMFT/ITS/tracking/test/testROFLookupTables.cxx index de02c77ff9acc..5e22005f2ce00 100644 --- a/Detectors/ITSMFT/ITS/tracking/test/testROFLookupTables.cxx +++ b/Detectors/ITSMFT/ITS/tracking/test/testROFLookupTables.cxx @@ -48,6 +48,40 @@ BOOST_AUTO_TEST_CASE(layertiming_base) BOOST_CHECK_EQUAL(layer1.mROFLength, 600); } +BOOST_AUTO_TEST_CASE(rofmask_construct_from_timing) +{ + o2::its::ROFOverlapTable<2> timing; + timing.defineLayer(0, 3, 100, 0, 0, 0); + timing.defineLayer(1, 4, 50, 25, 0, 0); + + o2::its::ROFMaskTable<2> mask{timing}; + const auto view = mask.getView(); + + BOOST_REQUIRE(view.mFlatMask != nullptr); + BOOST_REQUIRE(view.mLayerROFOffsets != nullptr); + BOOST_CHECK_EQUAL(view.mLayerROFOffsets[0], 0); + BOOST_CHECK_EQUAL(view.mLayerROFOffsets[1], 3); + BOOST_CHECK_EQUAL(view.mLayerROFOffsets[2], 7); + + for (int rof{0}; rof < 3; ++rof) { + BOOST_CHECK(!view.isROFEnabled(0, rof)); + } + for (int rof{0}; rof < 4; ++rof) { + BOOST_CHECK(!view.isROFEnabled(1, rof)); + } + + mask.selectROF({110, 20}); + + BOOST_CHECK(!view.isROFEnabled(0, 0)); + BOOST_CHECK(view.isROFEnabled(0, 1)); + BOOST_CHECK(!view.isROFEnabled(0, 2)); + + BOOST_CHECK(!view.isROFEnabled(1, 0)); + BOOST_CHECK(view.isROFEnabled(1, 1)); + BOOST_CHECK(view.isROFEnabled(1, 2)); + BOOST_CHECK(!view.isROFEnabled(1, 3)); +} + // ROFOverlapTable BOOST_AUTO_TEST_CASE(rofoverlap_basic) { diff --git a/Detectors/Upgrades/ALICE3/TRK/reconstruction/src/TimeFrame.cxx b/Detectors/Upgrades/ALICE3/TRK/reconstruction/src/TimeFrame.cxx index 6e8876f609b39..957560aea8cae 100644 --- a/Detectors/Upgrades/ALICE3/TRK/reconstruction/src/TimeFrame.cxx +++ b/Detectors/Upgrades/ALICE3/TRK/reconstruction/src/TimeFrame.cxx @@ -213,7 +213,7 @@ void TimeFrame::getPrimaryVerticesFromMC(TTree* mcHeaderTree, int nRofs iRof++; } } - this->mMultiplicityCutMask.resize(nRofs, true); /// all ROFs are valid with MC primary vertices. + this->mMultiplicityCutMask.resetMask(1u); /// all ROFs are valid with MC primary vertices. // Update the vertex lookup table with the newly added vertices this->updateROFVertexLookupTable(); From 32114cb099e617ae1f72c115d200e8ea92ad6597 Mon Sep 17 00:00:00 2001 From: Maximiliano Puccio Date: Thu, 2 Apr 2026 14:42:41 +0200 Subject: [PATCH 3/3] Adapt GPU code to the new mult mask --- .../tracking/GPU/ITStrackingGPU/TimeFrameGPU.h | 1 - .../GPU/ITStrackingGPU/TrackingKernels.h | 4 ++-- .../ITS/tracking/GPU/cuda/TrackerTraitsGPU.cxx | 4 ++-- .../ITS/tracking/GPU/cuda/TrackingKernels.cu | 14 +++++++++----- .../tracking/include/ITStracking/FastMultEst.h | 12 ++++++------ .../include/ITStracking/FastMultEstConfig.h | 16 ++++++++-------- .../ITS/tracking/src/TrackingInterface.cxx | 2 +- 7 files changed, 28 insertions(+), 25 deletions(-) diff --git a/Detectors/ITSMFT/ITS/tracking/GPU/ITStrackingGPU/TimeFrameGPU.h b/Detectors/ITSMFT/ITS/tracking/GPU/ITStrackingGPU/TimeFrameGPU.h index 12be726d11435..9bbb190dd3aca 100644 --- a/Detectors/ITSMFT/ITS/tracking/GPU/ITStrackingGPU/TimeFrameGPU.h +++ b/Detectors/ITSMFT/ITS/tracking/GPU/ITStrackingGPU/TimeFrameGPU.h @@ -143,7 +143,6 @@ class TimeFrameGPU final : public TimeFrame o2::track::TrackParCovF** getDeviceArrayTrackSeeds() { return mCellSeedsDeviceArray; } float** getDeviceArrayTrackSeedsChi2() { return mCellSeedsChi2DeviceArray; } int* getDeviceNeighboursIndexTables(const int layer) { return mNeighboursIndexTablesDevice[layer]; } - uint8_t* getDeviceMultCutMask() { return mMultMaskDevice; } void setDevicePropagator(const o2::base::PropagatorImpl* p) final { this->mPropagatorDevice = p; } diff --git a/Detectors/ITSMFT/ITS/tracking/GPU/ITStrackingGPU/TrackingKernels.h b/Detectors/ITSMFT/ITS/tracking/GPU/ITStrackingGPU/TrackingKernels.h index 043357743372c..8663f22f02d98 100644 --- a/Detectors/ITSMFT/ITS/tracking/GPU/ITStrackingGPU/TrackingKernels.h +++ b/Detectors/ITSMFT/ITS/tracking/GPU/ITStrackingGPU/TrackingKernels.h @@ -36,7 +36,7 @@ class ExternalAllocator; template void countTrackletsInROFsHandler(const IndexTableUtils* utils, - const uint8_t* multMask, + const typename ROFMaskTable::View& multMask, const int layer, const typename ROFOverlapTable::View& rofOverlaps, const typename ROFVertexLookupTable::View& vertexLUT, @@ -66,7 +66,7 @@ void countTrackletsInROFsHandler(const IndexTableUtils* utils, template void computeTrackletsInROFsHandler(const IndexTableUtils* utils, - const uint8_t* multMask, + const typename ROFMaskTable::View& multMask, const int layer, const typename ROFOverlapTable::View& rofOverlaps, const typename ROFVertexLookupTable::View& vertexLUT, diff --git a/Detectors/ITSMFT/ITS/tracking/GPU/cuda/TrackerTraitsGPU.cxx b/Detectors/ITSMFT/ITS/tracking/GPU/cuda/TrackerTraitsGPU.cxx index c5ef06e864de0..8f41fb2e619a4 100644 --- a/Detectors/ITSMFT/ITS/tracking/GPU/cuda/TrackerTraitsGPU.cxx +++ b/Detectors/ITSMFT/ITS/tracking/GPU/cuda/TrackerTraitsGPU.cxx @@ -85,7 +85,7 @@ void TrackerTraitsGPU::computeLayerTracklets(const int iteration, int i mTimeFrameGPU->createTrackletsLUTDevice(iteration, iLayer); mTimeFrameGPU->waitEvent(iLayer, iLayer + 1); // wait stream until all data is available countTrackletsInROFsHandler(mTimeFrameGPU->getDeviceIndexTableUtils(), - mTimeFrameGPU->getDeviceMultCutMask(), + mTimeFrameGPU->getDeviceROFMaskTableView(), iLayer, mTimeFrameGPU->getDeviceROFOverlapTableView(), mTimeFrameGPU->getDeviceROFVertexLookupTableView(), @@ -117,7 +117,7 @@ void TrackerTraitsGPU::computeLayerTracklets(const int iteration, int i continue; } computeTrackletsInROFsHandler(mTimeFrameGPU->getDeviceIndexTableUtils(), - mTimeFrameGPU->getDeviceMultCutMask(), + mTimeFrameGPU->getDeviceROFMaskTableView(), iLayer, mTimeFrameGPU->getDeviceROFOverlapTableView(), mTimeFrameGPU->getDeviceROFVertexLookupTableView(), diff --git a/Detectors/ITSMFT/ITS/tracking/GPU/cuda/TrackingKernels.cu b/Detectors/ITSMFT/ITS/tracking/GPU/cuda/TrackingKernels.cu index 0d2d4f116561a..69f76e53b8479 100644 --- a/Detectors/ITSMFT/ITS/tracking/GPU/cuda/TrackingKernels.cu +++ b/Detectors/ITSMFT/ITS/tracking/GPU/cuda/TrackingKernels.cu @@ -536,7 +536,7 @@ GPUg() void __launch_bounds__(256, 1) computeLayerCellsKernel( template GPUg() void __launch_bounds__(256, 1) computeLayerTrackletsMultiROFKernel( const IndexTableUtils* utils, - const uint8_t* multMask, + const typename ROFMaskTable::View multMask, const int layerIndex, const typename ROFOverlapTable::View rofOverlaps, const typename ROFVertexLookupTable::View vertexLUT, @@ -565,6 +565,10 @@ GPUg() void __launch_bounds__(256, 1) computeLayerTrackletsMultiROFKernel( const int totalROFs0 = rofOverlaps.getLayer(layerIndex).mNROFsTF; const int totalROFs1 = rofOverlaps.getLayer(layerIndex + 1).mNROFsTF; for (unsigned int pivotROF{blockIdx.x}; pivotROF < totalROFs0; pivotROF += gridDim.x) { + if (!multMask.isROFEnabled(layerIndex, pivotROF)) { + continue; + } + const auto& pvs = vertexLUT.getVertices(layerIndex, pivotROF); auto primaryVertices = gpuSpan(&vertices[pvs.getFirstEntry()], pvs.getEntries()); if (primaryVertices.empty()) { @@ -782,7 +786,7 @@ GPUg() void __launch_bounds__(256, 1) processNeighboursKernel( template void countTrackletsInROFsHandler(const IndexTableUtils* utils, - const uint8_t* multMask, + const typename ROFMaskTable::View& multMask, const int layer, const typename ROFOverlapTable::View& rofOverlaps, const typename ROFVertexLookupTable::View& vertexLUT, @@ -840,7 +844,7 @@ void countTrackletsInROFsHandler(const IndexTableUtils* utils, template void computeTrackletsInROFsHandler(const IndexTableUtils* utils, - const uint8_t* multMask, + const typename ROFMaskTable::View& multMask, const int layer, const typename ROFOverlapTable::View& rofOverlaps, const typename ROFVertexLookupTable::View& vertexLUT, @@ -1300,7 +1304,7 @@ void computeTrackSeedHandler(CellSeed* trackSeeds, /// Explicit instantiation of ITS2 handlers template void countTrackletsInROFsHandler<7>(const IndexTableUtils<7>* utils, - const uint8_t* multMask, + const ROFMaskTable<7>::View& multMask, const int layer, const ROFOverlapTable<7>::View& rofOverlaps, const ROFVertexLookupTable<7>::View& vertexLUT, @@ -1329,7 +1333,7 @@ template void countTrackletsInROFsHandler<7>(const IndexTableUtils<7>* utils, gpu::Streams& streams); template void computeTrackletsInROFsHandler<7>(const IndexTableUtils<7>* utils, - const uint8_t* multMask, + const ROFMaskTable<7>::View& multMask, const int layer, const ROFOverlapTable<7>::View& rofOverlaps, const ROFVertexLookupTable<7>::View& vertexLUT, diff --git a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/FastMultEst.h b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/FastMultEst.h index 46d6c0331f7d0..f71dda666efc7 100644 --- a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/FastMultEst.h +++ b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/FastMultEst.h @@ -37,12 +37,12 @@ struct FastMultEst { using ROFOverlapTableN = ROFOverlapTable; using ROFMaskTableN = ROFMaskTable; - float mult = 0.; /// estimated signal clusters multiplicity on the selected multiplicity layer - float noisePerChip = 0.; /// imposed noise per chip (when enabled by configuration) - float cov[3] = {0.}; /// retained for compatibility; set to zero in single-layer mode - float chi2 = 0.; /// retained for compatibility; set to zero in single-layer mode - int nLayersUsed = 0; /// number of layers used by estimator (0/1 in single-layer mode) - uint32_t lastRandomSeed = 0; /// state of the gRandom before + float mult = 0.; /// estimated signal clusters multiplicity on the selected multiplicity layer + float noisePerChip = 0.; /// imposed noise per chip (when enabled by configuration) + float cov[3] = {0.}; /// retained for compatibility; set to zero in single-layer mode + float chi2 = 0.; /// retained for compatibility; set to zero in single-layer mode + int nLayersUsed = 0; /// number of layers used by estimator (0/1 in single-layer mode) + uint32_t lastRandomSeed = 0; /// state of the gRandom before FastMultEst(); static uint32_t getCurrentRandomSeed(); diff --git a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/FastMultEstConfig.h b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/FastMultEstConfig.h index 7d30bcba45844..a76cd66d4df53 100644 --- a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/FastMultEstConfig.h +++ b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/FastMultEstConfig.h @@ -34,14 +34,14 @@ struct FastMultEstConfig : public o2::conf::ConfigurableParamHelper0 : set as is, <0 : use current time - bool preferTriggered = true; /// prefer ROFs with highest number of physics triggers + int cutMultClusLayer = NLayers - 1; /// layer used for cluster multiplicity selection (by default the outermost one) + float cutMultClusLow = 0; /// reject ROF with estimated cluster mult. below this value (no cut if <0) + float cutMultClusHigh = -1; /// reject ROF with estimated cluster mult. above this value (no cut if <0) + float cutMultVtxLow = -1; /// reject seed vertex if its multiplicity below this value (no cut if <0) + float cutMultVtxHigh = -1; /// reject seed vertex if its multiplicity above this value (no cut if <0) + float cutRandomFraction = -1.; /// apply random cut rejecting requested fraction + int randomSeed = 0; /// 0 - do not seet seed, >0 : set as is, <0 : use current time + bool preferTriggered = true; /// prefer ROFs with highest number of physics triggers bool isMultCutRequested() const { return cutMultClusLow >= 0.f && cutMultClusHigh > 0.f; }; bool isVtxMultCutRequested() const { return cutMultVtxLow >= 0.f && cutMultVtxHigh > 0.f; }; diff --git a/Detectors/ITSMFT/ITS/tracking/src/TrackingInterface.cxx b/Detectors/ITSMFT/ITS/tracking/src/TrackingInterface.cxx index 3c153d8ff6895..0bc6809f15028 100644 --- a/Detectors/ITSMFT/ITS/tracking/src/TrackingInterface.cxx +++ b/Detectors/ITSMFT/ITS/tracking/src/TrackingInterface.cxx @@ -214,7 +214,7 @@ void ITSTrackingInterface::run(framework::ProcessingContext& pc) if (vtxSpan[0].isFlagSet(Vertex::UPCMode) == 1) { // at least one vertex in this ROF and it is from second vertex iteration LOGP(debug, "ROF {} rejected as vertices are from the UPC iteration", iRof); processUPCMask.selectROF({clockTiming.getROFStartInBC(iRof), clockTiming.getROFEndInBC(iRof)}); - vtxROF.setFlag(o2::itsmft::ROFRecord::VtxUPCMode); + vtxROF.setFlag(o2::itsmft::ROFRecord::VtxUPCMode); } else { // in all cases except if as standard mode vertex was found, the ROF was processed with UPC settings vtxROF.setFlag(o2::itsmft::ROFRecord::VtxStdMode); }