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/include/ITSReconstruction/FastMultEst.h b/Detectors/ITSMFT/ITS/reconstruction/include/ITSReconstruction/FastMultEst.h deleted file mode 100644 index 9e8299e89b404..0000000000000 --- a/Detectors/ITSMFT/ITS/reconstruction/include/ITSReconstruction/FastMultEst.h +++ /dev/null @@ -1,70 +0,0 @@ -// Copyright 2019-2020 CERN and copyright holders of ALICE O2. -// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. -// All rights not expressly granted are reserved. -// -// This software is distributed under the terms of the GNU General Public -// License v3 (GPL Version 3), copied verbatim in the file "COPYING". -// -// In applying this license CERN does not waive the privileges and immunities -// granted to it by virtue of its status as an Intergovernmental Organization -// or submit itself to any jurisdiction. - -/// \file FastMultEst.h -/// \brief Fast multiplicity estimator for ITS -/// \author ruben.shahoyan@cern.ch - -#ifndef ALICEO2_ITS_FASTMULTEST_ -#define ALICEO2_ITS_FASTMULTEST_ - -#include "ITSMFTReconstruction/ChipMappingITS.h" -#include "DataFormatsITSMFT/ROFRecord.h" -#include "DataFormatsITSMFT/CompCluster.h" -#include -#include "ITSReconstruction/FastMultEstConfig.h" -#include -#include - -namespace o2 -{ -namespace its -{ - -struct FastMultEst { - - static constexpr int NLayers = o2::itsmft::ChipMappingITS::NLayers; - - 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 - 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); - - void fillNClPerLayer(const gsl::span& clusters); - float process(const std::array ncl) - { - return FastMultEstConfig::Instance().imposeNoisePerChip > 0 ? processNoiseImposed(ncl) : processNoiseFree(ncl); - } - float processNoiseFree(const std::array ncl); - float processNoiseImposed(const std::array ncl); - float process(const gsl::span& clusters) - { - fillNClPerLayer(clusters); - return process(nClPerLayer); - } - static bool sSeedSet; - - ClassDefNV(FastMultEst, 1); -}; - -} // namespace its -} // namespace o2 - -#endif diff --git a/Detectors/ITSMFT/ITS/reconstruction/src/FastMultEst.cxx b/Detectors/ITSMFT/ITS/reconstruction/src/FastMultEst.cxx deleted file mode 100644 index b62cacc28d9a7..0000000000000 --- a/Detectors/ITSMFT/ITS/reconstruction/src/FastMultEst.cxx +++ /dev/null @@ -1,189 +0,0 @@ -// Copyright 2019-2020 CERN and copyright holders of ALICE O2. -// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. -// All rights not expressly granted are reserved. -// -// This software is distributed under the terms of the GNU General Public -// License v3 (GPL Version 3), copied verbatim in the file "COPYING". -// -// In applying this license CERN does not waive the privileges and immunities -// granted to it by virtue of its status as an Intergovernmental Organization -// or submit itself to any jurisdiction. - -/// \file FastMultEst.h -/// \brief Fast multiplicity estimator for ITS -/// \author ruben.shahoyan@cern.ch - -#include "ITSReconstruction/FastMultEst.h" -#include "DataFormatsITSMFT/DPLAlpideParam.h" -#include "Framework/Logger.h" -#include -#include -#include - -using namespace o2::its; - -bool FastMultEst::sSeedSet = false; - -///______________________________________________________ -FastMultEst::FastMultEst() -{ - if (!sSeedSet && FastMultEstConfig::Instance().cutRandomFraction > 0.f) { - sSeedSet = true; - if (FastMultEstConfig::Instance().randomSeed > 0) { - gRandom->SetSeed(FastMultEstConfig::Instance().randomSeed); - } else if (FastMultEstConfig::Instance().randomSeed < 0) { - gRandom->SetSeed(std::time(nullptr) % 0xffff); - } - } -} - -///______________________________________________________ -/// find multiplicity for given set of clusters -void FastMultEst::fillNClPerLayer(const gsl::span& clusters) -{ - int lr = FastMultEst::NLayers - 1, nchAcc = o2::itsmft::ChipMappingITS::getNChips() - o2::itsmft::ChipMappingITS::getNChipsPerLr(lr); - std::memset(&nClPerLayer[0], 0, sizeof(int) * FastMultEst::NLayers); - 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]++; - } -} - -///______________________________________________________ -/// find multiplicity for given number of clusters per layer -float FastMultEst::processNoiseFree(const std::array ncl) -{ - // 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 - 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; - } - 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.; - 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) -{ - // 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 - // - 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; - } - return mult; -} - -int FastMultEst::selectROFs(const gsl::span rofs, const gsl::span clus, - const gsl::span trig, std::vector& 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 - lastRandomSeed = gRandom->GetSeed(); - if (multEstConf.isMultCutRequested()) { - for (uint32_t irof = 0; irof < nrof; irof++) { - nsel += sel[irof] = multEstConf.isPassingMultCut(process(rofs[irof].getROFData(clus))); - } - } 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--; - } - 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++; - } - } - } - 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--; - } - } - } - } - } - LOGP(debug, "NSel = {} of {} rofs Seeds: before {} after {}", nsel, nrof, lastRandomSeed, gRandom->GetSeed()); - - return nsel; -} 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/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 new file mode 100644 index 0000000000000..46d6c0331f7d0 --- /dev/null +++ b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/FastMultEst.h @@ -0,0 +1,96 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +/// \file FastMultEst.h +/// \brief Fast multiplicity estimator for ITS +/// \author ruben.shahoyan@cern.ch + +#ifndef ALICEO2_ITS_FASTMULTEST_ +#define ALICEO2_ITS_FASTMULTEST_ + +#include "ITSMFTReconstruction/ChipMappingITS.h" +#include "DataFormatsITS/Vertex.h" +#include "DataFormatsITSMFT/ROFRecord.h" +#include "DataFormatsITSMFT/CompCluster.h" +#include "DataFormatsITSMFT/PhysTrigger.h" +#include "ITStracking/FastMultEstConfig.h" +#include "ITStracking/ROFLookupTables.h" +#include +#include + +namespace o2 +{ +namespace its +{ + +struct FastMultEst { + + static constexpr int NLayers = o2::itsmft::ChipMappingITS::NLayers; + 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 + FastMultEst(); + + static uint32_t getCurrentRandomSeed(); + 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); + } + } + } + } + } + + int countClustersOnLayer(const gsl::span& clusters) const; + float process(int nClusters) + { + return FastMultEstConfig::Instance().imposeNoisePerChip > 0 ? processNoiseImposed(nClusters) : processNoiseFree(nClusters); + } + float processNoiseFree(int nClusters); + float processNoiseImposed(int nClusters); + float process(const gsl::span& clusters) + { + return process(countClustersOnLayer(clusters)); + } + static bool sSeedSet; + + ClassDefNV(FastMultEst, 1); +}; + +} // namespace its +} // namespace o2 + +#endif diff --git a/Detectors/ITSMFT/ITS/reconstruction/include/ITSReconstruction/FastMultEstConfig.h b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/FastMultEstConfig.h similarity index 95% rename from Detectors/ITSMFT/ITS/reconstruction/include/ITSReconstruction/FastMultEstConfig.h rename to Detectors/ITSMFT/ITS/tracking/include/ITStracking/FastMultEstConfig.h index c6bce50995a4b..7d30bcba45844 100644 --- a/Detectors/ITSMFT/ITS/reconstruction/include/ITSReconstruction/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, 0); + } + + 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; + } + + 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 new file mode 100644 index 0000000000000..0547209129b70 --- /dev/null +++ b/Detectors/ITSMFT/ITS/tracking/src/FastMultEst.cxx @@ -0,0 +1,251 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +/// \file FastMultEst.h +/// \brief Fast multiplicity estimator for ITS +/// \author ruben.shahoyan@cern.ch + +#include "ITStracking/FastMultEst.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; + +///______________________________________________________ +FastMultEst::FastMultEst() +{ + if (!sSeedSet && FastMultEstConfig::Instance().cutRandomFraction > 0.f) { + sSeedSet = true; + if (FastMultEstConfig::Instance().randomSeed > 0) { + gRandom->SetSeed(FastMultEstConfig::Instance().randomSeed); + } else if (FastMultEstConfig::Instance().randomSeed < 0) { + gRandom->SetSeed(std::time(nullptr) % 0xffff); + } + } +} + +///______________________________________________________ +/// count clusters on the configured multiplicity layer +int FastMultEst::countClustersOnLayer(const gsl::span& clusters) const +{ + 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); + } + if (lr == targetLayer) { + ++count; + } + } + return count; +} + +///______________________________________________________ +/// find multiplicity for given number of clusters per layer +float FastMultEst::processNoiseFree(int nClusters) +{ + // Single-layer regime: estimate multiplicity from one configured layer only. + const auto& conf = FastMultEstConfig::Instance(); + 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; + } + mult = nClusters / acc; + return mult > 0 ? mult : 0; +} + +///______________________________________________________ +/// find multiplicity for given number of clusters per layer with mean noise imposed +float FastMultEst::processNoiseImposed(int nClusters) +{ + // Single-layer regime with imposed noise subtraction. + const auto& conf = FastMultEstConfig::Instance(); + 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 std::array, NLayers>& rofs, + const std::array, NLayers>& clus, + const gsl::span trig, + uint32_t firstTForbit, + bool doStaggering, + const ROFOverlapTableN::View& overlapView, + ROFMaskTableN& sel) +{ + const auto& multEstConf = FastMultEstConfig::Instance(); // parameters for mult estimation and cuts + 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(); + 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 { + 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])); + } + } 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 (!passes) { + continue; + } + } + if (multEstConf.cutRandomFraction > 0.f && gRandom->Rndm() < multEstConf.cutRandomFraction) { + continue; + } + enableCompatibleROFs(selectionLayer, selectionRof, overlapView, sel); + } + } + + 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/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/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 375a35e0d55be..75ab1e0990e2e 100644 --- a/Detectors/ITSMFT/ITS/tracking/src/TrackingInterface.cxx +++ b/Detectors/ITSMFT/ITS/tracking/src/TrackingInterface.cxx @@ -17,9 +17,10 @@ #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/ROFLookupTables.h" #include "ITStracking/TrackingConfigParam.h" #include "ITStracking/TrackingInterface.h" @@ -176,12 +177,12 @@ 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); float vertexerElapsedTime{0.f}; if (mRunVertexer) { // Run seeding vertexer @@ -192,63 +193,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 +255,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 +317,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/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 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();