From a8b31f71871a82f2441c7ffc5c91d0d37766acd4 Mon Sep 17 00:00:00 2001 From: Felix Schlepper Date: Wed, 22 Apr 2026 15:56:28 +0200 Subject: [PATCH 1/3] ITSMFT: ITS3: reuse clusterizer Signed-off-by: Felix Schlepper --- .../common/reconstruction/CMakeLists.txt | 2 - .../include/ITSMFTReconstruction/Clusterer.h | 576 ++++++++++++++++-- .../include/ITSMFTReconstruction/LookUp.h | 3 + .../common/reconstruction/src/Clusterer.cxx | 496 +-------------- .../src/ITSMFTReconstructionLinkDef.h | 1 - .../include/ITSMFTWorkflow/STFDecoderSpec.h | 2 +- .../ITS3/reconstruction/CMakeLists.txt | 3 +- .../include/ITS3Reconstruction/Clusterer.h | 313 +--------- .../include/ITS3Reconstruction/LookUp.h | 7 + .../ITS3/reconstruction/src/Clusterer.cxx | 463 +------------- .../src/ITS3ReconstructionLinkDef.h | 1 - .../include/ITS3Workflow/ClustererSpec.h | 1 - 12 files changed, 566 insertions(+), 1302 deletions(-) diff --git a/Detectors/ITSMFT/common/reconstruction/CMakeLists.txt b/Detectors/ITSMFT/common/reconstruction/CMakeLists.txt index fa6abd958cb89..43ceeaf343fd7 100644 --- a/Detectors/ITSMFT/common/reconstruction/CMakeLists.txt +++ b/Detectors/ITSMFT/common/reconstruction/CMakeLists.txt @@ -49,7 +49,6 @@ o2_target_root_dictionary( include/ITSMFTReconstruction/DigitPixelReader.h include/ITSMFTReconstruction/RawPixelReader.h include/ITSMFTReconstruction/PixelData.h - include/ITSMFTReconstruction/Clusterer.h include/ITSMFTReconstruction/ClustererParam.h include/ITSMFTReconstruction/BuildTopologyDictionary.h include/ITSMFTReconstruction/LookUp.h @@ -67,4 +66,3 @@ if (OpenMP_CXX_FOUND) target_compile_definitions(${targetName} PRIVATE WITH_OPENMP) target_link_libraries(${targetName} PRIVATE OpenMP::OpenMP_CXX) endif() - diff --git a/Detectors/ITSMFT/common/reconstruction/include/ITSMFTReconstruction/Clusterer.h b/Detectors/ITSMFT/common/reconstruction/include/ITSMFTReconstruction/Clusterer.h index 0bdbb701a9356..1b82832c73dae 100644 --- a/Detectors/ITSMFT/common/reconstruction/include/ITSMFTReconstruction/Clusterer.h +++ b/Detectors/ITSMFT/common/reconstruction/include/ITSMFTReconstruction/Clusterer.h @@ -20,20 +20,29 @@ // | *| #define _ALLOW_DIAGONAL_ALPIDE_CLUSTERS_ +#include +#include +#include #include #include #include #include #include +#include "CommonDataFormat/InteractionRecord.h" #include "ITSMFTBase/SegmentationAlpide.h" #include "DataFormatsITSMFT/CompCluster.h" #include "DataFormatsITSMFT/ROFRecord.h" +#include "Framework/Logger.h" #include "ITSMFTReconstruction/PixelReader.h" #include "ITSMFTReconstruction/PixelData.h" #include "ITSMFTReconstruction/LookUp.h" #include "SimulationDataFormat/MCCompLabel.h" +#include "SimulationDataFormat/MCTruthContainer.h" #include "CommonConstants/LHCConstants.h" -#include "Rtypes.h" + +#ifdef WITH_OPENMP +#include +#endif #ifdef _PERFORM_TIMING_ #include @@ -59,9 +68,8 @@ using CompClusCont = std::vector; using PatternCont = std::vector; using ROFRecCont = std::vector; -// template // container types (PMR or std::vectors) - -class Clusterer +template +class ClustererT { using PixelReader = o2::itsmft::PixelReader; using PixelData = o2::itsmft::PixelData; @@ -94,18 +102,10 @@ class Clusterer } void adjust(uint16_t row, uint16_t col) { - if (row < rowMin) { - rowMin = row; - } - if (row > rowMax) { - rowMax = row; - } - if (col < colMin) { - colMin = col; - } - if (col > colMax) { - colMax = col; - } + rowMin = std::min(row, rowMin); + rowMax = std::max(row, rowMax); + colMin = std::min(col, colMin); + colMax = std::max(col, colMax); } }; @@ -126,11 +126,11 @@ class Clusterer int index = 0; }; int id = -1; - Clusterer* parent = nullptr; // parent clusterer + ClustererT* parent = nullptr; // parent clusterer // buffers for entries in preClusterIndices in 2 columns, to avoid boundary checks, we reserve // extra elements in the beginning and the end - int column1[SegmentationAlpide::NRows + 2]; - int column2[SegmentationAlpide::NRows + 2]; + int column1[MaxRows + 2]{}; + int column2[MaxRows + 2]{}; int* curr = nullptr; // pointer on the 1st row of currently processed columnsX int* prev = nullptr; // pointer on the 1st row of previously processed columnsX // pixels[].first is the index of the next pixel of the same precluster in the pixels @@ -149,7 +149,7 @@ class Clusterer std::vector stats; // statistics for each thread results, used at merging /// ///< reset column buffer, for the performance reasons we use memset - void resetColumn(int* buff) { std::memset(buff, -1, sizeof(int) * SegmentationAlpide::NRows); } + void resetColumn(int* buff) { std::memset(buff, -1, sizeof(int) * MaxRows); } ///< swap current and previous column buffers void swapColumnBuffers() { std::swap(prev, curr); } @@ -183,7 +183,7 @@ class Clusterer void process(uint16_t chip, uint16_t nChips, CompClusCont* compClusPtr, PatternCont* patternsPtr, const ConstMCTruth* labelsDigPtr, MCTruth* labelsClPtr, const ROFRecord& rofPtr); - ClustererThread(Clusterer* par = nullptr, int _id = -1) : parent(par), id(_id), curr(column2 + 1), prev(column1 + 1) + ClustererThread(ClustererT* par = nullptr, int _id = -1) : parent(par), id(_id), curr(column2 + 1), prev(column1 + 1) { std::fill(std::begin(column1), std::end(column1), -1); std::fill(std::begin(column2), std::end(column2), -1); @@ -191,16 +191,18 @@ class Clusterer }; //========================================================= - Clusterer(); - ~Clusterer() = default; + ClustererT(); + ClustererT(ClustererT&&) = delete; + ClustererT& operator=(ClustererT&&) = delete; + ~ClustererT() = default; - Clusterer(const Clusterer&) = delete; - Clusterer& operator=(const Clusterer&) = delete; + ClustererT(const ClustererT&) = delete; + ClustererT& operator=(const ClustererT&) = delete; void process(int nThreads, PixelReader& r, CompClusCont* compClus, PatternCont* patterns, ROFRecCont* vecROFRec, MCTruth* labelsCl = nullptr); template - static void streamCluster(const std::vector& pixbuf, const std::array* lblBuff, const BBox& bbox, const LookUp& pattIdConverter, + static void streamCluster(const std::vector& pixbuf, const std::array* lblBuff, const BBox& bbox, const LookUpT& pattIdConverter, VCLUS* compClusPtr, VPAT* patternsPtr, MCTruth* labelsClusPtr, int nlab, bool isHuge = false); bool isContinuousReadOut() const { return mContinuousReadout; } @@ -235,7 +237,11 @@ class Clusterer ///< load the dictionary of cluster topologies void loadDictionary(const std::string& fileName) { mPattIdConverter.loadDictionary(fileName); } - void setDictionary(const TopologyDictionary* dict) { mPattIdConverter.setDictionary(dict); } + template + void setDictionary(const TD* dict) + { + mPattIdConverter.setDictionary(dict); + } TStopwatch& getTimer() { return mTimer; } // cannot be const TStopwatch& getTimerMerge() { return mTimerMerge; } // cannot be const @@ -248,13 +254,13 @@ class Clusterer bool mDropHugeClusters = false; ///< don't include clusters that would be split in more than one ///< mask continuously fired pixels in frames separated by less than this amount of BCs (fired from hit in prev. ROF) - int mMaxBCSeparationToMask = 6000. / o2::constants::lhc::LHCBunchSpacingNS + 10; + int mMaxBCSeparationToMask = (6000.f / o2::constants::lhc::LHCBunchSpacingNS) + 10; int mMaxRowColDiffToMask = 0; ///< provide their difference in col/row is <= than this int mNHugeClus = 0; ///< number of encountered huge clusters ///< Squashing options int mSquashingDepth = 0; ///< squashing is applied to next N rofs - int mMaxBCSeparationToSquash = 6000. / o2::constants::lhc::LHCBunchSpacingNS + 10; + int mMaxBCSeparationToSquash = (6000.f / o2::constants::lhc::LHCBunchSpacingNS) + 10; std::vector mSquashingLayerDepth; std::vector mMaxBCSeparationToSquashLayer; @@ -263,15 +269,16 @@ class Clusterer std::vector mChipsOld; // previously processed ROF's chips data (for masking) std::vector mFiredChipsPtr; // pointers on the fired chips data in the decoder cache - LookUp mPattIdConverter; //! Convert the cluster topology to the corresponding entry in the dictionary. + LookUpT mPattIdConverter; //! Convert the cluster topology to the corresponding entry in the dictionary. TStopwatch mTimer; TStopwatch mTimerMerge; }; +template template -void Clusterer::streamCluster(const std::vector& pixbuf, const std::array* lblBuff, const Clusterer::BBox& bbox, const LookUp& pattIdConverter, - VCLUS* compClusPtr, VPAT* patternsPtr, MCTruth* labelsClusPtr, int nlab, bool isHuge) +void ClustererT::streamCluster(const std::vector& pixbuf, const std::array* lblBuff, const ClustererT::BBox& bbox, const LookUpT& pattIdConverter, + VCLUS* compClusPtr, VPAT* patternsPtr, MCTruth* labelsClusPtr, int nlab, bool isHuge) { if (labelsClusPtr && lblBuff) { // MC labels were requested auto cnt = compClusPtr->size(); @@ -288,15 +295,15 @@ void Clusterer::streamCluster(const std::vector& pixbuf, const std::a int nbits = ir * colSpanW + ic; patt[nbits >> 3] |= (0x1 << (7 - (nbits % 8))); } - uint16_t pattID = (isHuge || pattIdConverter.size() == 0) ? CompCluster::InvalidPatternID : pattIdConverter.findGroupID(rowSpanW, colSpanW, patt.data()); + uint16_t pattID = (isHuge || pattIdConverter.size(bbox.chipID) == 0) ? CompCluster::InvalidPatternID : pattIdConverter.findGroupID(rowSpanW, colSpanW, bbox.chipID, patt.data()); uint16_t row = bbox.rowMin, col = bbox.colMin; - if (pattID == CompCluster::InvalidPatternID || pattIdConverter.isGroup(pattID)) { + if (pattID == CompCluster::InvalidPatternID || pattIdConverter.isGroup(pattID, bbox.chipID)) { if (pattID != CompCluster::InvalidPatternID) { // For grouped topologies, the reference pixel is the COG pixel float xCOG = 0., zCOG = 0.; ClusterPattern::getCOG(rowSpanW, colSpanW, patt.data(), xCOG, zCOG); - row += round(xCOG); - col += round(zCOG); + row += std::round(xCOG); + col += std::round(zCOG); } if (patternsPtr) { patternsPtr->emplace_back((unsigned char)rowSpanW); @@ -311,6 +318,503 @@ void Clusterer::streamCluster(const std::vector& pixbuf, const std::a compClusPtr->emplace_back(row, col, pattID, bbox.chipID); } +//__________________________________________________ +template +void ClustererT::process(int nThreads, PixelReader& reader, CompClusCont* compClus, + PatternCont* patterns, ROFRecCont* vecROFRec, MCTruth* labelsCl) +{ +#ifdef _PERFORM_TIMING_ + mTimer.Start(kFALSE); +#endif + nThreads = std::max(nThreads, 1); + auto autoDecode = reader.getDecodeNextAuto(); + o2::InteractionRecord lastIR{}; + do { + if (autoDecode) { + reader.setDecodeNextAuto(false); // internally do not autodecode + if (!reader.decodeNextTrigger()) { + break; // on the fly decoding was requested, but there were no data left + } + } + if (reader.getInteractionRecord().isDummy()) { + continue; // No IR info was found + } + if (!lastIR.isDummy() && lastIR >= reader.getInteractionRecord()) { + const int MaxErrLog = 2; + static int errLocCount = 0; + if (errLocCount++ < MaxErrLog) { + LOGP(warn, "Impossible ROF IR {}, does not exceed previous {}, discarding in clusterization", reader.getInteractionRecord().asString(), lastIR.asString()); + } + continue; + } + lastIR = reader.getInteractionRecord(); + // pre-fetch all non-empty chips of current ROF + ChipPixelData* curChipData = nullptr; + mFiredChipsPtr.clear(); + size_t nPix = 0; + while ((curChipData = reader.getNextChipData(mChips))) { + mFiredChipsPtr.push_back(curChipData); + nPix += curChipData->getData().size(); + } + + auto& rof = vecROFRec->emplace_back(reader.getInteractionRecord(), vecROFRec->size(), compClus->size(), 0); // create new ROF + + uint16_t nFired = mFiredChipsPtr.size(); + if (!nFired) { + if (autoDecode) { + continue; + } + break; // just 1 ROF was asked to be processed + } + nThreads = std::min(nFired, nThreads); +#ifndef WITH_OPENMP + nThreads = 1; +#endif + uint16_t chipStep = nThreads > 1 ? (nThreads == 2 ? 20 : 10) : nFired; + int dynGrp = std::min(4, std::max(1, nThreads / 2)); + if (nThreads > mThreads.size()) { + int oldSz = mThreads.size(); + mThreads.resize(nThreads); + for (int i = oldSz; i < nThreads; i++) { + mThreads[i] = std::make_unique(this, i); + } + } +#ifdef WITH_OPENMP +#pragma omp parallel for schedule(dynamic, dynGrp) num_threads(nThreads) + //>> start of MT region + for (uint16_t ic = 0; ic < nFired; ic += chipStep) { + auto ith = omp_get_thread_num(); + if (nThreads > 1) { + mThreads[ith]->process(ic, std::min(chipStep, uint16_t(nFired - ic)), + &mThreads[ith]->compClusters, + patterns ? &mThreads[ith]->patterns : nullptr, + labelsCl ? reader.getDigitsMCTruth() : nullptr, + labelsCl ? &mThreads[ith]->labels : nullptr, rof); + } else { // put directly to the destination + mThreads[0]->process(0, nFired, compClus, patterns, labelsCl ? reader.getDigitsMCTruth() : nullptr, labelsCl, rof); + } + } + //<< end of MT region +#else + mThreads[0]->process(0, nFired, compClus, patterns, labelsCl ? reader.getDigitsMCTruth() : nullptr, labelsCl, rof); +#endif + // copy data of all threads but the 1st one to final destination + if (nThreads > 1) { +#ifdef _PERFORM_TIMING_ + mTimerMerge.Start(false); +#endif + size_t nClTot = 0, nPattTot = 0; + int chid = 0; + std::vector thrStatIdx(nThreads); + for (int ith = 0; ith < nThreads; ith++) { + std::sort(mThreads[ith]->stats.begin(), mThreads[ith]->stats.end(), [](const ThreadStat& a, const ThreadStat& b) { return a.firstChip < b.firstChip; }); + thrStatIdx[ith] = 0; + nClTot += mThreads[ith]->compClusters.size(); + nPattTot += mThreads[ith]->patterns.size(); + } + compClus->reserve(nClTot); + if (patterns) { + patterns->reserve(nPattTot); + } + while (chid < nFired) { + for (int ith = 0; ith < nThreads; ith++) { + if (thrStatIdx[ith] >= mThreads[ith]->stats.size()) { + continue; + } + const auto& stat = mThreads[ith]->stats[thrStatIdx[ith]]; + if (stat.firstChip == chid) { + thrStatIdx[ith]++; + chid += stat.nChips; // next chip to look + if (stat.nClus > 0) { + const auto clbeg = mThreads[ith]->compClusters.begin() + stat.firstClus; + compClus->insert(compClus->end(), clbeg, clbeg + stat.nClus); + if (patterns) { + const auto ptbeg = mThreads[ith]->patterns.begin() + stat.firstPatt; + patterns->insert(patterns->end(), ptbeg, ptbeg + stat.nPatt); + } + if (labelsCl) { + labelsCl->mergeAtBack(mThreads[ith]->labels, stat.firstClus, stat.nClus); + } + } + } + } + } + for (int ith = 0; ith < nThreads; ith++) { + mThreads[ith]->patterns.clear(); + mThreads[ith]->compClusters.clear(); + mThreads[ith]->labels.clear(); + mThreads[ith]->stats.clear(); + } +#ifdef _PERFORM_TIMING_ + mTimerMerge.Stop(); +#endif + } else { + mThreads[0]->stats.clear(); + } + rof.setNEntries(compClus->size() - rof.getFirstEntry()); // update + } while (autoDecode); + reader.setDecodeNextAuto(autoDecode); // restore setting +#ifdef _PERFORM_TIMING_ + mTimer.Stop(); +#endif +} + +//__________________________________________________ +template +void ClustererT::ClustererThread::process(uint16_t chip, uint16_t nChips, CompClusCont* compClusPtr, PatternCont* patternsPtr, + const ConstMCTruth* labelsDigPtr, MCTruth* labelsClPtr, const ROFRecord& rofPtr) +{ + if (stats.empty() || stats.back().firstChip + stats.back().nChips != chip) { // there is a jump, register new block + stats.emplace_back(ThreadStat{.firstChip = chip, .nChips = 0, .firstClus = uint32_t(compClusPtr->size()), .firstPatt = patternsPtr ? uint32_t(patternsPtr->size()) : 0, .nClus = 0, .nPatt = 0}); + } + for (int ic = 0; ic < nChips; ic++) { + auto* curChipData = parent->mFiredChipsPtr[chip + ic]; + auto chipID = curChipData->getChipID(); + if (parent->mMaxBCSeparationToMask > 0) { // mask pixels fired from the previous ROF + const auto& chipInPrevROF = parent->mChipsOld[chipID]; + if (std::abs(rofPtr.getBCData().differenceInBC(chipInPrevROF.getInteractionRecord())) < parent->mMaxBCSeparationToMask) { + parent->mMaxRowColDiffToMask ? curChipData->maskFiredInSample(parent->mChipsOld[chipID], parent->mMaxRowColDiffToMask) : curChipData->maskFiredInSample(parent->mChipsOld[chipID]); + } + } + auto validPixID = curChipData->getFirstUnmasked(); + auto npix = curChipData->getData().size(); + if (validPixID < npix) { // chip data may have all of its pixels masked! + auto valp = validPixID++; + if (validPixID == npix) { // special case of a single pixel fired on the chip + finishChipSingleHitFast(valp, curChipData, compClusPtr, patternsPtr, labelsDigPtr, labelsClPtr); + } else { + initChip(curChipData, valp); + for (; validPixID < npix; validPixID++) { + if (!curChipData->getData()[validPixID].isMasked()) { + updateChip(curChipData, validPixID); + } + } + finishChip(curChipData, compClusPtr, patternsPtr, labelsDigPtr, labelsClPtr); + } + } + if (parent->mMaxBCSeparationToMask > 0) { // current chip data will be used in the next ROF to mask overflow pixels + parent->mChipsOld[chipID].swap(*curChipData); + } + } + auto& currStat = stats.back(); + currStat.nChips += nChips; + currStat.nClus = compClusPtr->size() - currStat.firstClus; + currStat.nPatt = patternsPtr ? (patternsPtr->size() - currStat.firstPatt) : 0; +} + +//__________________________________________________ +template +void ClustererT::ClustererThread::finishChip(ChipPixelData* curChipData, CompClusCont* compClusPtr, + PatternCont* patternsPtr, const ConstMCTruth* labelsDigPtr, MCTruth* labelsClusPtr) +{ + const auto& pixData = curChipData->getData(); + int nPreclusters = preClusters.size(); + // account for the eventual reindexing of preClusters: Id2 might have been reindexed to Id1, which later was reindexed to Id0 + for (int i = 1; i < nPreclusters; i++) { + if (preClusters[i].index != i) { // reindexing is always done towards smallest index + preClusters[i].index = preClusters[preClusters[i].index].index; + } + } + for (int i1 = 0; i1 < nPreclusters; ++i1) { + auto& preCluster = preClusters[i1]; + auto ci = preCluster.index; + if (ci < 0) { + continue; + } + BBox bbox(curChipData->getChipID()); + int nlab = 0; + int next = preCluster.head; + pixArrBuff.clear(); + while (next >= 0) { + const auto& pixEntry = pixels[next]; + const auto pix = pixData[pixEntry.second]; + pixArrBuff.push_back(pix); // needed for cluster topology + bbox.adjust(pix.getRowDirect(), pix.getCol()); + if (labelsClusPtr) { + if (parent->mSquashingDepth) { // the MCtruth for this pixel is stored in chip data: due to squashing we lose contiguity + fetchMCLabels(curChipData->getOrderedPixId(pixEntry.second), labelsDigPtr, nlab); + } else { // the MCtruth for this pixel is at curChipData->startID+pixEntry.second + fetchMCLabels(pixEntry.second + curChipData->getStartID(), labelsDigPtr, nlab); + } + } + next = pixEntry.first; + } + preCluster.index = -1; + for (int i2 = i1 + 1; i2 < nPreclusters; ++i2) { + auto& preCluster2 = preClusters[i2]; + if (preCluster2.index != ci) { + continue; + } + next = preCluster2.head; + while (next >= 0) { + const auto& pixEntry = pixels[next]; + const auto pix = pixData[pixEntry.second]; // PixelData + pixArrBuff.push_back(pix); // needed for cluster topology + bbox.adjust(pix.getRowDirect(), pix.getCol()); + if (labelsClusPtr) { + if (parent->mSquashingDepth) { // the MCtruth for this pixel is stored in chip data: due to squashing we lose contiguity + fetchMCLabels(curChipData->getOrderedPixId(pixEntry.second), labelsDigPtr, nlab); + } else { // the MCtruth for this pixel is at curChipData->startID+pixEntry.second + fetchMCLabels(pixEntry.second + curChipData->getStartID(), labelsDigPtr, nlab); + } + } + next = pixEntry.first; + } + preCluster2.index = -1; + } + if (bbox.isAcceptableSize()) { + parent->streamCluster(pixArrBuff, &labelsBuff, bbox, parent->mPattIdConverter, compClusPtr, patternsPtr, labelsClusPtr, nlab); + } else { + auto warnLeft = MaxHugeClusWarn - parent->mNHugeClus; + if (!parent->mDropHugeClusters) { + if (warnLeft > 0) { + LOGP(warn, "Splitting a huge cluster: chipID {}, rows {}:{} cols {}:{}{}", bbox.chipID, bbox.rowMin, bbox.rowMax, bbox.colMin, bbox.colMax, + warnLeft == 1 ? " (Further warnings will be muted)" : ""); +#ifdef WITH_OPENMP +#pragma omp critical +#endif + { + parent->mNHugeClus++; + } + } + BBox bboxT(bbox); // truncated box + std::vector pixbuf; + do { + bboxT.rowMin = bbox.rowMin; + bboxT.colMax = std::min(bbox.colMax, uint16_t(bboxT.colMin + o2::itsmft::ClusterPattern::MaxColSpan - 1)); + do { // Select a subset of pixels fitting the reduced bounding box + bboxT.rowMax = std::min(bbox.rowMax, uint16_t(bboxT.rowMin + o2::itsmft::ClusterPattern::MaxRowSpan - 1)); + for (const auto& pix : pixArrBuff) { + if (bboxT.isInside(pix.getRowDirect(), pix.getCol())) { + pixbuf.push_back(pix); + } + } + if (!pixbuf.empty()) { // Stream a piece of cluster only if the reduced bounding box is not empty + parent->streamCluster(pixbuf, &labelsBuff, bboxT, parent->mPattIdConverter, compClusPtr, patternsPtr, labelsClusPtr, nlab, true); + pixbuf.clear(); + } + bboxT.rowMin = bboxT.rowMax + 1; + } while (bboxT.rowMin < bbox.rowMax); + bboxT.colMin = bboxT.colMax + 1; + } while (bboxT.colMin < bbox.colMax); + } + } + } +} + +//__________________________________________________ +template +void ClustererT::ClustererThread::finishChipSingleHitFast(uint32_t hit, ChipPixelData* curChipData, CompClusCont* compClusPtr, + PatternCont* patternsPtr, const ConstMCTruth* labelsDigPtr, MCTruth* labelsClusPtr) +{ + auto pix = curChipData->getData()[hit]; + uint16_t row = pix.getRowDirect(), col = pix.getCol(); + + if (labelsClusPtr) { // MC labels were requested + int nlab = 0; + fetchMCLabels(curChipData->getStartID() + hit, labelsDigPtr, nlab); + auto cnt = compClusPtr->size(); + for (int i = nlab; i--;) { + labelsClusPtr->addElement(cnt, labelsBuff[i]); + } + } + + // add to compact clusters, which must be always filled + unsigned char patt[ClusterPattern::MaxPatternBytes]{0x1 << (7 - (0 % 8))}; // unrolled 1 hit version of full loop in finishChip + const auto chipID = curChipData->getChipID(); + uint16_t pattID = (parent->mPattIdConverter.size(chipID) == 0) ? CompCluster::InvalidPatternID : parent->mPattIdConverter.findGroupID(1, 1, chipID, patt); + if ((pattID == CompCluster::InvalidPatternID || parent->mPattIdConverter.isGroup(pattID, chipID)) && patternsPtr) { + patternsPtr->emplace_back(1); // rowspan + patternsPtr->emplace_back(1); // colspan + patternsPtr->insert(patternsPtr->end(), std::begin(patt), std::begin(patt) + 1); + } + compClusPtr->emplace_back(row, col, pattID, chipID); +} + +//__________________________________________________ +template +ClustererT::ClustererT() +{ +#ifdef _PERFORM_TIMING_ + mTimer.Stop(); + mTimer.Reset(); + mTimerMerge.Stop(); + mTimerMerge.Reset(); +#endif +} + +//__________________________________________________ +template +void ClustererT::ClustererThread::initChip(const ChipPixelData* curChipData, uint32_t first) +{ + // init chip with the 1st unmasked pixel (entry "from" in the mChipData) + prev = column1 + 1; + curr = column2 + 1; + resetColumn(curr); + pixels.clear(); + preClusters.clear(); + auto pix = curChipData->getData()[first]; + currCol = pix.getCol(); + curr[pix.getRowDirect()] = 0; // can use getRowDirect since the pixel is not masked + // start the first pre-cluster + preClusters.emplace_back(); + pixels.emplace_back(-1, first); // id of current pixel + noLeftCol = true; +} + +//__________________________________________________ +template +void ClustererT::ClustererThread::updateChip(const ChipPixelData* curChipData, uint32_t ip) +{ + const auto pix = curChipData->getData()[ip]; + uint16_t row = pix.getRowDirect(); // can use getRowDirect since the pixel is not masked + if (currCol != pix.getCol()) { // switch the buffers + swapColumnBuffers(); + resetColumn(curr); + noLeftCol = false; + if (pix.getCol() > currCol + 1) { + // no connection with previous column, this pixel cannot belong to any of the + // existing preclusters, create a new precluster and flag to check only the row above for next pixels of this column + currCol = pix.getCol(); + addNewPrecluster(ip, row); + noLeftCol = true; + return; + } + currCol = pix.getCol(); + } + + if (noLeftCol) { // check only the row above + if (curr[row - 1] >= 0) { + expandPreCluster(ip, row, curr[row - 1]); // attach to the precluster of the previous row + } else { + addNewPrecluster(ip, row); // start new precluster + } + } else { + // row above should be always checked + int nnb = 0, lowestIndex = curr[row - 1], lowestNb = 0, *nbrCol[4], nbrRow[4]; + if (lowestIndex >= 0) { + nbrCol[nnb] = curr; + nbrRow[nnb++] = row - 1; + } else { + lowestIndex = 0x7ffff; + lowestNb = -1; + } +#ifdef _ALLOW_DIAGONAL_ALPIDE_CLUSTERS_ + for (int i : {-1, 0, 1}) { + auto v = prev[row + i]; + if (v >= 0) { + nbrCol[nnb] = prev; + nbrRow[nnb] = row + i; + if (v < lowestIndex) { + lowestIndex = v; + lowestNb = nnb; + } + nnb++; + } + } +#else + if (prev[row] >= 0) { + nbrCol[nnb] = prev; + nbrRow[nnb] = row; + if (prev[row] < lowestIndex) { + lowestIndex = prev[row]; + lowestNb = nnb; + } + nnb++; + } +#endif + if (!nnb) { // no neighbours, create new precluster + addNewPrecluster(ip, row); // start new precluster + } else { + expandPreCluster(ip, row, lowestIndex); // attach to the adjascent precluster with smallest index + if (nnb > 1) { + for (int inb = 0; inb < nnb; inb++) { // reassign precluster index to smallest one, replicating updated values to columns caches + auto& prevIndex = (nbrCol[inb])[nbrRow[inb]]; + prevIndex = preClusters[prevIndex].index = lowestIndex; + } + } + } + } +} + +//__________________________________________________ +template +void ClustererT::ClustererThread::fetchMCLabels(int digID, const ConstMCTruth* labelsDig, int& nfilled) +{ + // transfer MC labels to cluster + if (nfilled >= MaxLabels) { + return; + } + const auto& lbls = labelsDig->getLabels(digID); + for (int i = lbls.size(); i--;) { + int ic = nfilled; + for (; ic--;) { // check if the label is already present + if (labelsBuff[ic] == lbls[i]) { + return; // label is found, do nothing + } + } + labelsBuff[nfilled++] = lbls[i]; + if (nfilled >= MaxLabels) { + break; + } + } + // +} + +//__________________________________________________ +template +void ClustererT::clear() +{ + // reset +#ifdef _PERFORM_TIMING_ + mTimer.Stop(); + mTimer.Reset(); + mTimerMerge.Stop(); + mTimerMerge.Reset(); +#endif +} + +//__________________________________________________ +template +void ClustererT::print(bool showsTiming) const +{ + // print settings + if (mSquashingLayerDepth.empty()) { + LOGP(info, "Clusterizer squashes overflow pixels separated by {} BC and <= {} in row/col seeking down to {} neighbour ROFs", mMaxBCSeparationToSquash, mMaxRowColDiffToMask, mSquashingDepth); + } else { + LOGP(info, "Clusterizer squashes overflow pixels <= {} in row/col", mMaxRowColDiffToMask); + for (size_t i{0}; i < mSquashingLayerDepth.size(); ++i) { + LOGP(info, "\tClusterizer on layer {} separated by {} BC seeking down to {} neighbour ROFs", i, mMaxBCSeparationToSquashLayer[i], mSquashingLayerDepth[i]); + } + } + LOGP(info, "Clusterizer masks overflow pixels separated by < {} BC and <= {} in row/col", mMaxBCSeparationToMask, mMaxRowColDiffToMask); + LOGP(info, "Clusterizer does {} drop huge clusters", mDropHugeClusters ? "" : "not"); + + if (showsTiming) { +#ifdef _PERFORM_TIMING_ + auto& tmr = const_cast(mTimer); // ugly but this is what root does internally + auto& tmrm = const_cast(mTimerMerge); + LOG(info) << "Inclusive clusterization timing (w/o disk IO): Cpu: " << tmr.CpuTime() + << " Real: " << tmr.RealTime() << " s in " << tmr.Counter() << " slots"; + LOG(info) << "Threads output merging timing : Cpu: " << tmrm.CpuTime() + << " Real: " << tmrm.RealTime() << " s in " << tmrm.Counter() << " slots"; + +#endif + } +} + +//__________________________________________________ +template +void ClustererT::reset() +{ + // reset for new run + clear(); + mNHugeClus = 0; +} + +using Clusterer = ClustererT; +extern template class ClustererT; + } // namespace itsmft } // namespace o2 #endif /* ALICEO2_ITS_CLUSTERER_H */ diff --git a/Detectors/ITSMFT/common/reconstruction/include/ITSMFTReconstruction/LookUp.h b/Detectors/ITSMFT/common/reconstruction/include/ITSMFTReconstruction/LookUp.h index 3537a1f408886..08499cc02d2bd 100644 --- a/Detectors/ITSMFT/common/reconstruction/include/ITSMFTReconstruction/LookUp.h +++ b/Detectors/ITSMFT/common/reconstruction/include/ITSMFTReconstruction/LookUp.h @@ -37,11 +37,14 @@ class LookUp LookUp(std::string fileName); static int groupFinder(int nRow, int nCol); int findGroupID(int nRow, int nCol, const unsigned char patt[ClusterPattern::MaxPatternBytes]) const; + int findGroupID(int nRow, int nCol, uint16_t, const unsigned char patt[ClusterPattern::MaxPatternBytes]) const { return findGroupID(nRow, nCol, patt); } int getTopologiesOverThreshold() const { return mTopologiesOverThreshold; } void loadDictionary(std::string fileName); void setDictionary(const TopologyDictionary* dict); bool isGroup(int id) const { return mDictionary.isGroup(id); } + bool isGroup(int id, uint16_t) const { return isGroup(id); } int size() const { return mDictionary.getSize(); } + int size(uint16_t) const { return size(); } auto getPattern(int id) const { return mDictionary.getPattern(id); } auto getDictionaty() const { return mDictionary; } diff --git a/Detectors/ITSMFT/common/reconstruction/src/Clusterer.cxx b/Detectors/ITSMFT/common/reconstruction/src/Clusterer.cxx index dcc268a4504a9..6093f07104e8c 100644 --- a/Detectors/ITSMFT/common/reconstruction/src/Clusterer.cxx +++ b/Detectors/ITSMFT/common/reconstruction/src/Clusterer.cxx @@ -10,498 +10,6 @@ // or submit itself to any jurisdiction. /// \file Clusterer.cxx -/// \brief Implementation of the ITS cluster finder -#include -#include -#include "Framework/Logger.h" +/// \brief Explicit instantiation of the ITS cluster finder template #include "ITSMFTReconstruction/Clusterer.h" -#include "SimulationDataFormat/MCTruthContainer.h" -#include "CommonDataFormat/InteractionRecord.h" - -#ifdef WITH_OPENMP -#include -#endif -using namespace o2::itsmft; - -//__________________________________________________ -void Clusterer::process(int nThreads, PixelReader& reader, CompClusCont* compClus, - PatternCont* patterns, ROFRecCont* vecROFRec, MCTruth* labelsCl) -{ -#ifdef _PERFORM_TIMING_ - mTimer.Start(kFALSE); -#endif - nThreads = std::max(nThreads, 1); - auto autoDecode = reader.getDecodeNextAuto(); - o2::InteractionRecord lastIR{}; - do { - if (autoDecode) { - reader.setDecodeNextAuto(false); // internally do not autodecode - if (!reader.decodeNextTrigger()) { - break; // on the fly decoding was requested, but there were no data left - } - } - if (reader.getInteractionRecord().isDummy()) { - continue; // No IR info was found - } - if (!lastIR.isDummy() && lastIR >= reader.getInteractionRecord()) { - const int MaxErrLog = 2; - static int errLocCount = 0; - if (errLocCount++ < MaxErrLog) { - LOGP(warn, "Impossible ROF IR {}, does not exceed previous {}, discarding in clusterization", reader.getInteractionRecord().asString(), lastIR.asString()); - } - continue; - } - lastIR = reader.getInteractionRecord(); - // pre-fetch all non-empty chips of current ROF - ChipPixelData* curChipData = nullptr; - mFiredChipsPtr.clear(); - size_t nPix = 0; - while ((curChipData = reader.getNextChipData(mChips))) { - mFiredChipsPtr.push_back(curChipData); - nPix += curChipData->getData().size(); - } - - auto& rof = vecROFRec->emplace_back(reader.getInteractionRecord(), vecROFRec->size(), compClus->size(), 0); // create new ROF - - uint16_t nFired = mFiredChipsPtr.size(); - if (!nFired) { - if (autoDecode) { - continue; - } - break; // just 1 ROF was asked to be processed - } - nThreads = std::min(nFired, nThreads); -#ifndef WITH_OPENMP - nThreads = 1; -#endif - uint16_t chipStep = nThreads > 1 ? (nThreads == 2 ? 20 : 10) : nFired; - int dynGrp = std::min(4, std::max(1, nThreads / 2)); - if (nThreads > mThreads.size()) { - int oldSz = mThreads.size(); - mThreads.resize(nThreads); - for (int i = oldSz; i < nThreads; i++) { - mThreads[i] = std::make_unique(this, i); - } - } -#ifdef WITH_OPENMP -#pragma omp parallel for schedule(dynamic, dynGrp) num_threads(nThreads) - //>> start of MT region - for (uint16_t ic = 0; ic < nFired; ic += chipStep) { - auto ith = omp_get_thread_num(); - if (nThreads > 1) { - mThreads[ith]->process(ic, std::min(chipStep, uint16_t(nFired - ic)), - &mThreads[ith]->compClusters, - patterns ? &mThreads[ith]->patterns : nullptr, - labelsCl ? reader.getDigitsMCTruth() : nullptr, - labelsCl ? &mThreads[ith]->labels : nullptr, rof); - } else { // put directly to the destination - mThreads[0]->process(0, nFired, compClus, patterns, labelsCl ? reader.getDigitsMCTruth() : nullptr, labelsCl, rof); - } - } - //<< end of MT region -#else - mThreads[0]->process(0, nFired, compClus, patterns, labelsCl ? reader.getDigitsMCTruth() : nullptr, labelsCl, rof); -#endif - // copy data of all threads but the 1st one to final destination - if (nThreads > 1) { -#ifdef _PERFORM_TIMING_ - mTimerMerge.Start(false); -#endif - size_t nClTot = 0, nPattTot = 0; - int chid = 0, thrStatIdx[nThreads]; - for (int ith = 0; ith < nThreads; ith++) { - std::sort(mThreads[ith]->stats.begin(), mThreads[ith]->stats.end(), [](const ThreadStat& a, const ThreadStat& b) { return a.firstChip < b.firstChip; }); - thrStatIdx[ith] = 0; - nClTot += mThreads[ith]->compClusters.size(); - nPattTot += mThreads[ith]->patterns.size(); - } - compClus->reserve(nClTot); - if (patterns) { - patterns->reserve(nPattTot); - } - while (chid < nFired) { - for (int ith = 0; ith < nThreads; ith++) { - if (thrStatIdx[ith] >= mThreads[ith]->stats.size()) { - continue; - } - const auto& stat = mThreads[ith]->stats[thrStatIdx[ith]]; - if (stat.firstChip == chid) { - thrStatIdx[ith]++; - chid += stat.nChips; // next chip to look - if (stat.nClus > 0) { - const auto clbeg = mThreads[ith]->compClusters.begin() + stat.firstClus; - auto szold = compClus->size(); - compClus->insert(compClus->end(), clbeg, clbeg + stat.nClus); - if (patterns) { - const auto ptbeg = mThreads[ith]->patterns.begin() + stat.firstPatt; - patterns->insert(patterns->end(), ptbeg, ptbeg + stat.nPatt); - } - if (labelsCl) { - labelsCl->mergeAtBack(mThreads[ith]->labels, stat.firstClus, stat.nClus); - } - } - } - } - } - for (int ith = 0; ith < nThreads; ith++) { - mThreads[ith]->patterns.clear(); - mThreads[ith]->compClusters.clear(); - mThreads[ith]->labels.clear(); - mThreads[ith]->stats.clear(); - } -#ifdef _PERFORM_TIMING_ - mTimerMerge.Stop(); -#endif - } else { - mThreads[0]->stats.clear(); - } - rof.setNEntries(compClus->size() - rof.getFirstEntry()); // update - } while (autoDecode); - reader.setDecodeNextAuto(autoDecode); // restore setting -#ifdef _PERFORM_TIMING_ - mTimer.Stop(); -#endif -} - -//__________________________________________________ -void Clusterer::ClustererThread::process(uint16_t chip, uint16_t nChips, CompClusCont* compClusPtr, PatternCont* patternsPtr, - const ConstMCTruth* labelsDigPtr, MCTruth* labelsClPtr, const ROFRecord& rofPtr) -{ - if (stats.empty() || stats.back().firstChip + stats.back().nChips != chip) { // there is a jump, register new block - stats.emplace_back(ThreadStat{.firstChip = chip, .nChips = 0, .firstClus = uint32_t(compClusPtr->size()), .firstPatt = patternsPtr ? uint32_t(patternsPtr->size()) : 0, .nClus = 0, .nPatt = 0}); - } - for (int ic = 0; ic < nChips; ic++) { - auto* curChipData = parent->mFiredChipsPtr[chip + ic]; - auto chipID = curChipData->getChipID(); - if (parent->mMaxBCSeparationToMask > 0) { // mask pixels fired from the previous ROF - const auto& chipInPrevROF = parent->mChipsOld[chipID]; - if (std::abs(rofPtr.getBCData().differenceInBC(chipInPrevROF.getInteractionRecord())) < parent->mMaxBCSeparationToMask) { - parent->mMaxRowColDiffToMask ? curChipData->maskFiredInSample(parent->mChipsOld[chipID], parent->mMaxRowColDiffToMask) : curChipData->maskFiredInSample(parent->mChipsOld[chipID]); - } - } - auto nclus0 = compClusPtr->size(); - auto validPixID = curChipData->getFirstUnmasked(); - auto npix = curChipData->getData().size(); - if (validPixID < npix) { // chip data may have all of its pixels masked! - auto valp = validPixID++; - if (validPixID == npix) { // special case of a single pixel fired on the chip - finishChipSingleHitFast(valp, curChipData, compClusPtr, patternsPtr, labelsDigPtr, labelsClPtr); - } else { - initChip(curChipData, valp); - for (; validPixID < npix; validPixID++) { - if (!curChipData->getData()[validPixID].isMasked()) { - updateChip(curChipData, validPixID); - } - } - finishChip(curChipData, compClusPtr, patternsPtr, labelsDigPtr, labelsClPtr); - } - } - if (parent->mMaxBCSeparationToMask > 0) { // current chip data will be used in the next ROF to mask overflow pixels - parent->mChipsOld[chipID].swap(*curChipData); - } - } - auto& currStat = stats.back(); - currStat.nChips += nChips; - currStat.nClus = compClusPtr->size() - currStat.firstClus; - currStat.nPatt = patternsPtr ? (patternsPtr->size() - currStat.firstPatt) : 0; -} - -//__________________________________________________ -void Clusterer::ClustererThread::finishChip(ChipPixelData* curChipData, CompClusCont* compClusPtr, - PatternCont* patternsPtr, const ConstMCTruth* labelsDigPtr, MCTruth* labelsClusPtr) -{ - const auto& pixData = curChipData->getData(); - int nPreclusters = preClusters.size(); - // account for the eventual reindexing of preClusters: Id2 might have been reindexed to Id1, which later was reindexed to Id0 - for (int i = 1; i < nPreclusters; i++) { - if (preClusters[i].index != i) { // reindexing is always done towards smallest index - preClusters[i].index = preClusters[preClusters[i].index].index; - } - } - for (int i1 = 0; i1 < nPreclusters; ++i1) { - auto& preCluster = preClusters[i1]; - auto ci = preCluster.index; - if (ci < 0) { - continue; - } - BBox bbox(curChipData->getChipID()); - int nlab = 0; - int next = preCluster.head; - pixArrBuff.clear(); - while (next >= 0) { - const auto& pixEntry = pixels[next]; - const auto pix = pixData[pixEntry.second]; - pixArrBuff.push_back(pix); // needed for cluster topology - bbox.adjust(pix.getRowDirect(), pix.getCol()); - if (labelsClusPtr) { - if (parent->mSquashingDepth) { // the MCtruth for this pixel is stored in chip data: due to squashing we lose contiguity - fetchMCLabels(curChipData->getOrderedPixId(pixEntry.second), labelsDigPtr, nlab); - } else { // the MCtruth for this pixel is at curChipData->startID+pixEntry.second - fetchMCLabels(pixEntry.second + curChipData->getStartID(), labelsDigPtr, nlab); - } - } - next = pixEntry.first; - } - preCluster.index = -1; - for (int i2 = i1 + 1; i2 < nPreclusters; ++i2) { - auto& preCluster2 = preClusters[i2]; - if (preCluster2.index != ci) { - continue; - } - next = preCluster2.head; - while (next >= 0) { - const auto& pixEntry = pixels[next]; - const auto pix = pixData[pixEntry.second]; // PixelData - pixArrBuff.push_back(pix); // needed for cluster topology - bbox.adjust(pix.getRowDirect(), pix.getCol()); - if (labelsClusPtr) { - if (parent->mSquashingDepth) { // the MCtruth for this pixel is stored in chip data: due to squashing we lose contiguity - fetchMCLabels(curChipData->getOrderedPixId(pixEntry.second), labelsDigPtr, nlab); - } else { // the MCtruth for this pixel is at curChipData->startID+pixEntry.second - fetchMCLabels(pixEntry.second + curChipData->getStartID(), labelsDigPtr, nlab); - } - } - next = pixEntry.first; - } - preCluster2.index = -1; - } - if (bbox.isAcceptableSize()) { - parent->streamCluster(pixArrBuff, &labelsBuff, bbox, parent->mPattIdConverter, compClusPtr, patternsPtr, labelsClusPtr, nlab); - } else { - auto warnLeft = MaxHugeClusWarn - parent->mNHugeClus; - if (!parent->mDropHugeClusters) { - if (warnLeft > 0) { - LOGP(warn, "Splitting a huge cluster: chipID {}, rows {}:{} cols {}:{}{}", bbox.chipID, bbox.rowMin, bbox.rowMax, bbox.colMin, bbox.colMax, - warnLeft == 1 ? " (Further warnings will be muted)" : ""); -#ifdef WITH_OPENMP -#pragma omp critical -#endif - { - parent->mNHugeClus++; - } - } - BBox bboxT(bbox); // truncated box - std::vector pixbuf; - do { - bboxT.rowMin = bbox.rowMin; - bboxT.colMax = std::min(bbox.colMax, uint16_t(bboxT.colMin + o2::itsmft::ClusterPattern::MaxColSpan - 1)); - do { // Select a subset of pixels fitting the reduced bounding box - bboxT.rowMax = std::min(bbox.rowMax, uint16_t(bboxT.rowMin + o2::itsmft::ClusterPattern::MaxRowSpan - 1)); - for (const auto& pix : pixArrBuff) { - if (bboxT.isInside(pix.getRowDirect(), pix.getCol())) { - pixbuf.push_back(pix); - } - } - if (!pixbuf.empty()) { // Stream a piece of cluster only if the reduced bounding box is not empty - parent->streamCluster(pixbuf, &labelsBuff, bboxT, parent->mPattIdConverter, compClusPtr, patternsPtr, labelsClusPtr, nlab, true); - pixbuf.clear(); - } - bboxT.rowMin = bboxT.rowMax + 1; - } while (bboxT.rowMin < bbox.rowMax); - bboxT.colMin = bboxT.colMax + 1; - } while (bboxT.colMin < bbox.colMax); - } - } - } -} - -//__________________________________________________ -void Clusterer::ClustererThread::finishChipSingleHitFast(uint32_t hit, ChipPixelData* curChipData, CompClusCont* compClusPtr, - PatternCont* patternsPtr, const ConstMCTruth* labelsDigPtr, MCTruth* labelsClusPtr) -{ - auto pix = curChipData->getData()[hit]; - uint16_t row = pix.getRowDirect(), col = pix.getCol(); - - if (labelsClusPtr) { // MC labels were requested - int nlab = 0; - fetchMCLabels(curChipData->getStartID() + hit, labelsDigPtr, nlab); - auto cnt = compClusPtr->size(); - for (int i = nlab; i--;) { - labelsClusPtr->addElement(cnt, labelsBuff[i]); - } - } - - // add to compact clusters, which must be always filled - unsigned char patt[ClusterPattern::MaxPatternBytes]{0x1 << (7 - (0 % 8))}; // unrolled 1 hit version of full loop in finishChip - uint16_t pattID = (parent->mPattIdConverter.size() == 0) ? CompCluster::InvalidPatternID : parent->mPattIdConverter.findGroupID(1, 1, patt); - if ((pattID == CompCluster::InvalidPatternID || parent->mPattIdConverter.isGroup(pattID)) && patternsPtr) { - patternsPtr->emplace_back(1); // rowspan - patternsPtr->emplace_back(1); // colspan - patternsPtr->insert(patternsPtr->end(), std::begin(patt), std::begin(patt) + 1); - } - compClusPtr->emplace_back(row, col, pattID, curChipData->getChipID()); -} - -//__________________________________________________ -Clusterer::Clusterer() : mPattIdConverter() -{ -#ifdef _PERFORM_TIMING_ - mTimer.Stop(); - mTimer.Reset(); - mTimerMerge.Stop(); - mTimerMerge.Reset(); -#endif -} - -//__________________________________________________ -void Clusterer::ClustererThread::initChip(const ChipPixelData* curChipData, uint32_t first) -{ - // init chip with the 1st unmasked pixel (entry "from" in the mChipData) - prev = column1 + 1; - curr = column2 + 1; - resetColumn(curr); - pixels.clear(); - preClusters.clear(); - auto pix = curChipData->getData()[first]; - currCol = pix.getCol(); - curr[pix.getRowDirect()] = 0; // can use getRowDirect since the pixel is not masked - // start the first pre-cluster - preClusters.emplace_back(); - pixels.emplace_back(-1, first); // id of current pixel - noLeftCol = true; -} - -//__________________________________________________ -void Clusterer::ClustererThread::updateChip(const ChipPixelData* curChipData, uint32_t ip) -{ - const auto pix = curChipData->getData()[ip]; - uint16_t row = pix.getRowDirect(); // can use getRowDirect since the pixel is not masked - if (currCol != pix.getCol()) { // switch the buffers - swapColumnBuffers(); - resetColumn(curr); - noLeftCol = false; - if (pix.getCol() > currCol + 1) { - // no connection with previous column, this pixel cannot belong to any of the - // existing preclusters, create a new precluster and flag to check only the row above for next pixels of this column - currCol = pix.getCol(); - addNewPrecluster(ip, row); - noLeftCol = true; - return; - } - currCol = pix.getCol(); - } - - if (noLeftCol) { // check only the row above - if (curr[row - 1] >= 0) { - expandPreCluster(ip, row, curr[row - 1]); // attach to the precluster of the previous row - } else { - addNewPrecluster(ip, row); // start new precluster - } - } else { - // row above should be always checked - int nnb = 0, lowestIndex = curr[row - 1], lowestNb = 0, *nbrCol[4], nbrRow[4]; - if (lowestIndex >= 0) { - nbrCol[nnb] = curr; - nbrRow[nnb++] = row - 1; - } else { - lowestIndex = 0x7ffff; - lowestNb = -1; - } -#ifdef _ALLOW_DIAGONAL_ALPIDE_CLUSTERS_ - for (int i : {-1, 0, 1}) { - auto v = prev[row + i]; - if (v >= 0) { - nbrCol[nnb] = prev; - nbrRow[nnb] = row + i; - if (v < lowestIndex) { - lowestIndex = v; - lowestNb = nnb; - } - nnb++; - } - } -#else - if (prev[row] >= 0) { - nbrCol[nnb] = prev; - nbrRow[nnb] = row; - if (prev[row] < lowestIndex) { - lowestIndex = v; - lowestNb = nnb; - } - nnb++; - } -#endif - if (!nnb) { // no neighbours, create new precluster - addNewPrecluster(ip, row); // start new precluster - } else { - expandPreCluster(ip, row, lowestIndex); // attach to the adjascent precluster with smallest index - if (nnb > 1) { - for (int inb = 0; inb < nnb; inb++) { // reassign precluster index to smallest one, replicating updated values to columns caches - auto& prevIndex = (nbrCol[inb])[nbrRow[inb]]; - prevIndex = preClusters[prevIndex].index = lowestIndex; - } - } - } - } -} - -//__________________________________________________ -void Clusterer::ClustererThread::fetchMCLabels(int digID, const ConstMCTruth* labelsDig, int& nfilled) -{ - // transfer MC labels to cluster - if (nfilled >= MaxLabels) { - return; - } - const auto& lbls = labelsDig->getLabels(digID); - for (int i = lbls.size(); i--;) { - int ic = nfilled; - for (; ic--;) { // check if the label is already present - if (labelsBuff[ic] == lbls[i]) { - return; // label is found, do nothing - } - } - labelsBuff[nfilled++] = lbls[i]; - if (nfilled >= MaxLabels) { - break; - } - } - // -} - -//__________________________________________________ -void Clusterer::clear() -{ - // reset -#ifdef _PERFORM_TIMING_ - mTimer.Stop(); - mTimer.Reset(); - mTimerMerge.Stop(); - mTimerMerge.Reset(); -#endif -} - -//__________________________________________________ -void Clusterer::print(bool showsTiming) const -{ - // print settings - if (mSquashingLayerDepth.empty()) { - LOGP(info, "Clusterizer squashes overflow pixels separated by {} BC and <= {} in row/col seeking down to {} neighbour ROFs", mMaxBCSeparationToSquash, mMaxRowColDiffToMask, mSquashingDepth); - } else { - LOGP(info, "Clusterizer squashes overflow pixels <= {} in row/col", mMaxRowColDiffToMask); - for (size_t i{0}; i < mSquashingLayerDepth.size(); ++i) { - LOGP(info, "\tClusterizer on layer {} separated by {} BC seeking down to {} neighbour ROFs", i, mMaxBCSeparationToSquashLayer[i], mSquashingLayerDepth[i]); - } - } - LOGP(info, "Clusterizer masks overflow pixels separated by < {} BC and <= {} in row/col", mMaxBCSeparationToMask, mMaxRowColDiffToMask); - LOGP(info, "Clusterizer does {} drop huge clusters", mDropHugeClusters ? "" : "not"); - - if (showsTiming) { -#ifdef _PERFORM_TIMING_ - auto& tmr = const_cast(mTimer); // ugly but this is what root does internally - auto& tmrm = const_cast(mTimerMerge); - LOG(info) << "Inclusive clusterization timing (w/o disk IO): Cpu: " << tmr.CpuTime() - << " Real: " << tmr.RealTime() << " s in " << tmr.Counter() << " slots"; - LOG(info) << "Threads output merging timing : Cpu: " << tmrm.CpuTime() - << " Real: " << tmrm.RealTime() << " s in " << tmrm.Counter() << " slots"; - -#endif - } -} - -//__________________________________________________ -void Clusterer::reset() -{ - // reset for new run - clear(); - mNHugeClus = 0; -} +template class o2::itsmft::ClustererT; diff --git a/Detectors/ITSMFT/common/reconstruction/src/ITSMFTReconstructionLinkDef.h b/Detectors/ITSMFT/common/reconstruction/src/ITSMFTReconstructionLinkDef.h index 19f4ca06d0220..87e540541388e 100644 --- a/Detectors/ITSMFT/common/reconstruction/src/ITSMFTReconstructionLinkDef.h +++ b/Detectors/ITSMFT/common/reconstruction/src/ITSMFTReconstructionLinkDef.h @@ -15,7 +15,6 @@ #pragma link off all classes; #pragma link off all functions; -#pragma link C++ class o2::itsmft::Clusterer + ; #pragma link C++ class o2::itsmft::PixelReader + ; #pragma link C++ class o2::itsmft::DigitPixelReader + ; #pragma link C++ class o2::itsmft::RawPixelReader < o2::itsmft::ChipMappingITS> + ; diff --git a/Detectors/ITSMFT/common/workflow/include/ITSMFTWorkflow/STFDecoderSpec.h b/Detectors/ITSMFT/common/workflow/include/ITSMFTWorkflow/STFDecoderSpec.h index 29b9f75bcbc4e..a9f6788b4b9fb 100644 --- a/Detectors/ITSMFT/common/workflow/include/ITSMFTWorkflow/STFDecoderSpec.h +++ b/Detectors/ITSMFT/common/workflow/include/ITSMFTWorkflow/STFDecoderSpec.h @@ -24,6 +24,7 @@ #include "Framework/Task.h" #include "DataFormatsITSMFT/DPLAlpideParam.h" #include "DataFormatsITSMFT/ROFRecord.h" +#include "ITSMFTReconstruction/Clusterer.h" #include "ITSMFTReconstruction/ChipMappingITS.h" #include "ITSMFTReconstruction/ChipMappingMFT.h" #include "ITSMFTReconstruction/RawPixelDecoder.h" @@ -38,7 +39,6 @@ class GRPGeomRequest; } namespace itsmft { -class Clusterer; struct STFDecoderInp { bool doClusters = true; diff --git a/Detectors/Upgrades/ITS3/reconstruction/CMakeLists.txt b/Detectors/Upgrades/ITS3/reconstruction/CMakeLists.txt index a666e15cc21e4..f3bd4be820510 100644 --- a/Detectors/Upgrades/ITS3/reconstruction/CMakeLists.txt +++ b/Detectors/Upgrades/ITS3/reconstruction/CMakeLists.txt @@ -33,8 +33,7 @@ o2_add_library(ITS3Reconstruction o2_target_root_dictionary( ITS3Reconstruction - HEADERS include/ITS3Reconstruction/Clusterer.h - include/ITS3Reconstruction/TopologyDictionary.h + HEADERS include/ITS3Reconstruction/TopologyDictionary.h include/ITS3Reconstruction/BuildTopologyDictionary.h include/ITS3Reconstruction/LookUp.h include/ITS3Reconstruction/IOUtils.h diff --git a/Detectors/Upgrades/ITS3/reconstruction/include/ITS3Reconstruction/Clusterer.h b/Detectors/Upgrades/ITS3/reconstruction/include/ITS3Reconstruction/Clusterer.h index a81db09217e9b..0e99c0f4596de 100644 --- a/Detectors/Upgrades/ITS3/reconstruction/include/ITS3Reconstruction/Clusterer.h +++ b/Detectors/Upgrades/ITS3/reconstruction/include/ITS3Reconstruction/Clusterer.h @@ -1,4 +1,4 @@ -// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// Copyright 2019-2026 CERN and copyright holders of ALICE O2. // See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. // All rights not expressly granted are reserved. // @@ -9,312 +9,19 @@ // granted to it by virtue of its status as an Intergovernmental Organization // or submit itself to any jurisdiction. -/// \file Clusterer.h -/// \brief Definition of the ITS cluster finder -#ifndef ALICEO2_ITS_CLUSTERER_H -#define ALICEO2_ITS_CLUSTERER_H - -#define _PERFORM_TIMING_ - -// uncomment this to not allow diagonal clusters, e.g. like |* | -// | *| -#define _ALLOW_DIAGONAL_ALPIDE_CLUSTERS_ - -#include -#include -#include -#include -#include +#ifndef ALICEO2_ITS3_CLUSTERER_H +#define ALICEO2_ITS3_CLUSTERER_H #include "ITSMFTBase/SegmentationAlpide.h" -#include "DataFormatsITSMFT/CompCluster.h" -#include "DataFormatsITSMFT/ROFRecord.h" -#include "ITSMFTReconstruction/PixelReader.h" -#include "ITSMFTReconstruction/PixelData.h" +#include "ITSMFTReconstruction/Clusterer.h" +#include "ITS3Base/SegmentationMosaix.h" #include "ITS3Reconstruction/LookUp.h" -#include "SimulationDataFormat/MCCompLabel.h" -#include "CommonConstants/LHCConstants.h" - -#ifdef _PERFORM_TIMING_ -#include -#endif - -class TTree; - -namespace o2 -{ -class MCCompLabel; -namespace dataformats -{ -template -class ConstMCTruthContainerView; -template -class MCTruthContainer; -} // namespace dataformats - -namespace its3 -{ -using CompClusCont = std::vector; -using PatternCont = std::vector; -using ROFRecCont = std::vector; - -// template // container types (PMR or std::vectors) - -class Clusterer +namespace o2::its3 { - using PixelReader = o2::itsmft::PixelReader; - using PixelData = o2::itsmft::PixelData; - using ChipPixelData = o2::itsmft::ChipPixelData; - using CompCluster = o2::itsmft::CompCluster; - using CompClusterExt = o2::itsmft::CompClusterExt; - using Label = o2::MCCompLabel; - using MCTruth = o2::dataformats::MCTruthContainer; - using ConstMCTruth = o2::dataformats::ConstMCTruthContainerView; - - public: - static constexpr int MaxLabels = 10; - static constexpr int MaxHugeClusWarn = 5; // max number of warnings for HugeCluster - - struct BBox { - uint16_t chipID = 0xffff; - uint16_t rowMin = 0xffff; - uint16_t colMin = 0xffff; - uint16_t rowMax = 0; - uint16_t colMax = 0; - BBox(uint16_t c) : chipID(c) {} - bool isInside(uint16_t row, uint16_t col) const { return row >= rowMin && row <= rowMax && col >= colMin && col <= colMax; } - [[nodiscard]] auto rowSpan() const { return rowMax - rowMin + 1; } - [[nodiscard]] auto colSpan() const { return colMax - colMin + 1; } - bool isAcceptableSize() const { return colMax - colMin < o2::itsmft::ClusterPattern::MaxColSpan && rowMax - rowMin < o2::itsmft::ClusterPattern::MaxRowSpan; } - void clear() - { - rowMin = colMin = 0xffff; - rowMax = colMax = 0; - } - void adjust(uint16_t row, uint16_t col) - { - if (row < rowMin) { - rowMin = row; - } - if (row > rowMax) { - rowMax = row; - } - if (col < colMin) { - colMin = col; - } - if (col > colMax) { - colMax = col; - } - } - }; - - //========================================================= - /// methods and transient data used within a thread - struct ThreadStat { - uint16_t firstChip = 0; - uint16_t nChips = 0; - uint32_t firstClus = 0; - uint32_t firstPatt = 0; - uint32_t nClus = 0; - uint32_t nPatt = 0; - }; - - struct ClustererThread { - Clusterer* parent = nullptr; // parent clusterer - int id = -1; - // buffers for entries in preClusterIndices in 2 columns, to avoid boundary checks, we reserve - // extra elements in the beginning and the end - int* column1 = nullptr; - int* column2 = nullptr; - int* curr = nullptr; // pointer on the 1st row of currently processed columnsX - int* prev = nullptr; // pointer on the 1st row of previously processed columnsX - int size = itsmft::SegmentationAlpide::NRows + 2; - // pixels[].first is the index of the next pixel of the same precluster in the pixels - // pixels[].second is the index of the referred pixel in the ChipPixelData (element of mChips) - std::vector> pixels; - std::vector preClusterHeads; // index of precluster head in the pixels - std::vector preClusterIndices; - uint16_t currCol = 0xffff; ///< Column being processed - bool noLeftCol = true; ///< flag that there is no column on the left to check - std::array labelsBuff; //! temporary buffer for building cluster labels - std::vector pixArrBuff; //! temporary buffer for pattern calc. - // - /// temporary storage for the thread output - CompClusCont compClusters; - PatternCont patterns; - MCTruth labels; - std::vector stats; // statistics for each thread results, used at merging - /// - ///< reset column buffer, for the performance reasons we use memset - void resetColumn(int* buff) const { std::memset(buff, -1, sizeof(int) * (size - 2)); } - - ///< swap current and previous column buffers - void swapColumnBuffers() { std::swap(prev, curr); } - - ///< add cluster at row (entry ip in the ChipPixeData) to the precluster with given index - void expandPreCluster(uint32_t ip, uint16_t row, int preClusIndex) - { - auto& firstIndex = preClusterHeads[preClusterIndices[preClusIndex]]; - pixels.emplace_back(firstIndex, ip); - firstIndex = pixels.size() - 1; - curr[row] = preClusIndex; - } - - ///< add new precluster at given row of current column for the fired pixel with index ip in the ChipPixelData - void addNewPrecluster(uint32_t ip, uint16_t row) - { - preClusterHeads.push_back(pixels.size()); - // new head does not point yet (-1) on other pixels, store just the entry of the pixel in the ChipPixelData - pixels.emplace_back(-1, ip); - int lastIndex = preClusterIndices.size(); - preClusterIndices.push_back(lastIndex); - curr[row] = lastIndex; // store index of the new precluster in the current column buffer - } - - void fetchMCLabels(int digID, const ConstMCTruth* labelsDig, int& nfilled); - void initChip(const ChipPixelData* curChipData, uint32_t first); - void updateChip(const ChipPixelData* curChipData, uint32_t ip); - void finishChip(ChipPixelData* curChipData, CompClusCont* compClus, PatternCont* patterns, - const ConstMCTruth* labelsDig, MCTruth* labelsClus); - void finishChipSingleHitFast(uint32_t hit, ChipPixelData* curChipData, CompClusCont* compClusPtr, - PatternCont* patternsPtr, const ConstMCTruth* labelsDigPtr, MCTruth* labelsClusPTr); - void process(uint16_t chip, uint16_t nChips, CompClusCont* compClusPtr, PatternCont* patternsPtr, - const ConstMCTruth* labelsDigPtr, MCTruth* labelsClPtr, const itsmft::ROFRecord& rofPtr); - - ~ClustererThread() - { - delete[] column1; - delete[] column2; - } +// Maximum number of rows is determined by Alpide +using Clusterer = o2::itsmft::ClustererT; +} // namespace o2::its3 +extern template class o2::itsmft::ClustererT; - ClustererThread(Clusterer* par = nullptr, int _id = -1) : parent(par), id(_id), curr(column2 + 1), prev(column1 + 1) {} - ClustererThread(const ClustererThread&) = delete; - ClustererThread(ClustererThread&&) = delete; - ClustererThread& operator=(const ClustererThread&) = delete; - ClustererThread& operator=(ClustererThread&&) = delete; - }; - //========================================================= - - Clusterer(); - Clusterer(Clusterer&&) = delete; - Clusterer& operator=(Clusterer&&) = delete; - ~Clusterer() = default; - Clusterer(const Clusterer&) = delete; - Clusterer& operator=(const Clusterer&) = delete; - - void process(int nThreads, PixelReader& r, CompClusCont* compClus, PatternCont* patterns, ROFRecCont* vecROFRec, MCTruth* labelsCl = nullptr); - - template - static void streamCluster(const std::vector& pixbuf, const std::array* lblBuff, const BBox& bbox, const its3::LookUp& pattIdConverter, - VCLUS* compClusPtr, VPAT* patternsPtr, MCTruth* labelsClusPtr, int nlab, bool isIB, bool isHuge = false); - - bool isContinuousReadOut() const { return mContinuousReadout; } - void setContinuousReadOut(bool v) { mContinuousReadout = v; } - - int getMaxBCSeparationToMask() const { return mMaxBCSeparationToMask; } - void setMaxBCSeparationToMask(int n) { mMaxBCSeparationToMask = n; } - - int getMaxRowColDiffToMask() const { return mMaxRowColDiffToMask; } - void setMaxRowColDiffToMask(int v) { mMaxRowColDiffToMask = v; } - - int getMaxROFDepthToSquash() const { return mSquashingDepth; } - void setMaxROFDepthToSquash(int v) { mSquashingDepth = v; } - - int getMaxBCSeparationToSquash() const { return mMaxBCSeparationToSquash; } - void setMaxBCSeparationToSquash(int n) { mMaxBCSeparationToSquash = n; } - - void print() const; - void clear(); - - ///< load the dictionary of cluster topologies - void setDictionary(const its3::TopologyDictionary* dict) - { - LOGP(info, "Setting TopologyDictionary with IB size={} & OB size={}", dict->getSize(true), dict->getSize(false)); - mPattIdConverter.setDictionary(dict); - // dict->print(); - } - - TStopwatch& getTimer() { return mTimer; } // cannot be const - TStopwatch& getTimerMerge() { return mTimerMerge; } // cannot be const - - void setNChips(int n) - { - mChips.resize(n); - mChipsOld.resize(n); - } - - private: - void flushClusters(CompClusCont* compClus, MCTruth* labels); - - // clusterization options - bool mContinuousReadout = true; ///< flag continuous readout - - ///< mask continuosly fired pixels in frames separated by less than this amount of BCs (fired from hit in prev. ROF) - int mMaxBCSeparationToMask = static_cast(6000. / o2::constants::lhc::LHCBunchSpacingNS + 10); - int mMaxRowColDiffToMask = 0; ///< provide their difference in col/row is <= than this - int mNHugeClus = 0; ///< number of encountered huge clusters - - ///< Squashing options - int mSquashingDepth = 0; ///< squashing is applied to next N rofs - int mMaxBCSeparationToSquash = 6000. / o2::constants::lhc::LHCBunchSpacingNS + 10; - - std::vector> mThreads; // buffers for threads - std::vector mChips; // currently processed ROF's chips data - std::vector mChipsOld; // previously processed ROF's chips data (for masking) - std::vector mFiredChipsPtr; // pointers on the fired chips data in the decoder cache - - its3::LookUp mPattIdConverter; //! Convert the cluster topology to the corresponding entry in the dictionary. - // -#ifdef _PERFORM_TIMING_ - TStopwatch mTimer; - TStopwatch mTimerMerge; #endif -}; - -template -void Clusterer::streamCluster(const std::vector& pixbuf, const std::array* lblBuff, const Clusterer::BBox& bbox, const its3::LookUp& pattIdConverter, - VCLUS* compClusPtr, VPAT* patternsPtr, MCTruth* labelsClusPtr, int nlab, bool isIB, bool isHuge) -{ - if (labelsClusPtr && lblBuff) { // MC labels were requested - auto cnt = compClusPtr->size(); - for (int i = nlab; i--;) { - labelsClusPtr->addElement(cnt, (*lblBuff)[i]); - } - } - auto colSpanW = bbox.colSpan(); - auto rowSpanW = bbox.rowSpan(); - // add to compact clusters, which must be always filled - std::array patt{}; - for (const auto& pix : pixbuf) { - uint32_t ir = pix.getRowDirect() - bbox.rowMin, ic = pix.getCol() - bbox.colMin; - int nbits = ir * colSpanW + ic; - patt[nbits >> 3] |= (0x1 << (7 - (nbits % 8))); - } - uint16_t pattID = (isHuge || pattIdConverter.size(isIB) == 0) ? CompCluster::InvalidPatternID : pattIdConverter.findGroupID(rowSpanW, colSpanW, isIB, patt.data()); - uint16_t row = bbox.rowMin, col = bbox.colMin; - LOGP(debug, "PattID: findGroupID({},{},{})={}", row, col, patt[0], pattID); - if (pattID == CompCluster::InvalidPatternID || pattIdConverter.isGroup(pattID, isIB)) { - if (pattID != CompCluster::InvalidPatternID) { - // For groupped topologies, the reference pixel is the COG pixel - float xCOG = 0., zCOG = 0.; - itsmft::ClusterPattern::getCOG(rowSpanW, colSpanW, patt.data(), xCOG, zCOG); - row += round(xCOG); - col += round(zCOG); - } - if (patternsPtr) { - patternsPtr->emplace_back((unsigned char)rowSpanW); - patternsPtr->emplace_back((unsigned char)colSpanW); - int nBytes = rowSpanW * colSpanW / 8; - if (((rowSpanW * colSpanW) % 8) != 0) { - nBytes++; - } - patternsPtr->insert(patternsPtr->end(), std::begin(patt), std::begin(patt) + nBytes); - } - } - compClusPtr->emplace_back(row, col, pattID, bbox.chipID); -} - -} // namespace its3 -} // namespace o2 -#endif /* ALICEO2_ITS_CLUSTERER_H */ diff --git a/Detectors/Upgrades/ITS3/reconstruction/include/ITS3Reconstruction/LookUp.h b/Detectors/Upgrades/ITS3/reconstruction/include/ITS3Reconstruction/LookUp.h index 809a129a0debf..37e1eadd8c31c 100644 --- a/Detectors/Upgrades/ITS3/reconstruction/include/ITS3Reconstruction/LookUp.h +++ b/Detectors/Upgrades/ITS3/reconstruction/include/ITS3Reconstruction/LookUp.h @@ -21,6 +21,7 @@ #ifndef ALICEO2_ITS3_LOOKUP_H #define ALICEO2_ITS3_LOOKUP_H +#include "ITS3Base/SpecsV2.h" #include "ITS3Reconstruction/TopologyDictionary.h" namespace o2::its3 @@ -32,12 +33,18 @@ class LookUp LookUp(std::string fileName); static int groupFinder(int nRow, int nCol); int findGroupID(int nRow, int nCol, bool IB, const unsigned char patt[itsmft::ClusterPattern::MaxPatternBytes]) const; + int findGroupID(int nRow, int nCol, uint16_t chipID, const unsigned char patt[itsmft::ClusterPattern::MaxPatternBytes]) const + { + return findGroupID(nRow, nCol, constants::detID::isDetITS3(chipID), patt); + } int getTopologiesOverThreshold(bool IB) const { return (IB) ? mTopologiesOverThresholdIB : mTopologiesOverThresholdOB; } void loadDictionary(std::string fileName); void setDictionary(const TopologyDictionary* dict); auto getDictionary() const { return mDictionary; } bool isGroup(int id, bool IB) const { return mDictionary.isGroup(id, IB); } + bool isGroup(int id, uint16_t chipID) const { return isGroup(id, constants::detID::isDetITS3(chipID)); } int size(bool IB) const { return mDictionary.getSize(IB); } + int size(uint16_t chipID) const { return size(constants::detID::isDetITS3(chipID)); } auto getPattern(int id, bool IB) const { return mDictionary.getPattern(id, IB); } private: diff --git a/Detectors/Upgrades/ITS3/reconstruction/src/Clusterer.cxx b/Detectors/Upgrades/ITS3/reconstruction/src/Clusterer.cxx index bce17b3759340..00c0006306fda 100644 --- a/Detectors/Upgrades/ITS3/reconstruction/src/Clusterer.cxx +++ b/Detectors/Upgrades/ITS3/reconstruction/src/Clusterer.cxx @@ -1,4 +1,4 @@ -// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// Copyright 2019-2026 CERN and copyright holders of ALICE O2. // See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. // All rights not expressly granted are reserved. // @@ -9,464 +9,5 @@ // granted to it by virtue of its status as an Intergovernmental Organization // or submit itself to any jurisdiction. -/// \file Clusterer.cxx -/// \brief Implementation of the ITS cluster finder - -#include - #include "ITS3Reconstruction/Clusterer.h" -#include "ITS3Base/SegmentationMosaix.h" -#include "SimulationDataFormat/MCTruthContainer.h" -#include "CommonDataFormat/InteractionRecord.h" - -#include "TTree.h" - -#ifdef WITH_OPENMP -#include -#endif -using namespace o2::itsmft; - -namespace o2::its3 -{ -void Clusterer::process(int nThreads, PixelReader& reader, CompClusCont* compClus, - PatternCont* patterns, ROFRecCont* vecROFRec, MCTruth* labelsCl) -{ -#ifdef _PERFORM_TIMING_ - mTimer.Start(false); -#endif - if (nThreads < 1) { - nThreads = 1; - } - auto autoDecode = reader.getDecodeNextAuto(); - do { - if (autoDecode) { - reader.setDecodeNextAuto(false); // internally do not autodecode - if (!reader.decodeNextTrigger()) { - break; // on the fly decoding was requested, but there were no data left - } - } - if (reader.getInteractionRecord().isDummy()) { - continue; // No IR info was found - } - // pre-fetch all non-empty chips of current ROF - ChipPixelData* curChipData = nullptr; - mFiredChipsPtr.clear(); - size_t nPix = 0; - while ((curChipData = reader.getNextChipData(mChips))) { - mFiredChipsPtr.push_back(curChipData); - nPix += curChipData->getData().size(); - } - - auto& rof = vecROFRec->emplace_back(reader.getInteractionRecord(), vecROFRec->size(), compClus->size(), 0); // create new ROF - - uint16_t nFired = mFiredChipsPtr.size(); - if (!nFired) { - if (autoDecode) { - continue; - } - break; // just 1 ROF was asked to be processed - } - if (nFired < nThreads) { - nThreads = nFired; - } -#ifndef WITH_OPENMP - nThreads = 1; -#endif - uint16_t chipStep = nThreads > 1 ? (nThreads == 2 ? 20 : 10) : nFired; - int dynGrp = std::min(4, std::max(1, nThreads / 2)); - if (nThreads > mThreads.size()) { - int oldSz = mThreads.size(); - mThreads.resize(nThreads); - for (size_t i = oldSz; i < nThreads; i++) { - mThreads[i] = std::make_unique(this, i); - } - } -#ifdef WITH_OPENMP -#pragma omp parallel for schedule(dynamic, dynGrp) num_threads(nThreads) - //>> start of MT region - for (uint16_t ic = 0; ic < nFired; ic += chipStep) { - auto ith = omp_get_thread_num(); - if (nThreads > 1) { - mThreads[ith]->process(ic, std::min(chipStep, uint16_t(nFired - ic)), - &mThreads[ith]->compClusters, - patterns ? &mThreads[ith]->patterns : nullptr, - labelsCl ? reader.getDigitsMCTruth() : nullptr, - labelsCl ? &mThreads[ith]->labels : nullptr, rof); - } else { // put directly to the destination - mThreads[0]->process(0, nFired, compClus, patterns, labelsCl ? reader.getDigitsMCTruth() : nullptr, labelsCl, rof); - } - } - //<< end of MT region -#else - mThreads[0]->process(0, nFired, compClus, patterns, labelsCl ? reader.getDigitsMCTruth() : nullptr, labelsCl, rof); -#endif - // copy data of all threads but the 1st one to final destination - if (nThreads > 1) { -#ifdef _PERFORM_TIMING_ - mTimerMerge.Start(false); -#endif - size_t nClTot = 0, nPattTot = 0; - int chid = 0; - std::vector thrStatIdx(nThreads, 0); - for (int ith = 0; ith < nThreads; ith++) { - std::sort(mThreads[ith]->stats.begin(), mThreads[ith]->stats.end(), [](const ThreadStat& a, const ThreadStat& b) { return a.firstChip < b.firstChip; }); - nClTot += mThreads[ith]->compClusters.size(); - nPattTot += mThreads[ith]->patterns.size(); - } - compClus->reserve(nClTot); - if (patterns) { - patterns->reserve(nPattTot); - } - while (chid < nFired) { - for (int ith = 0; ith < nThreads; ith++) { - if (thrStatIdx[ith] >= mThreads[ith]->stats.size()) { - continue; - } - const auto& stat = mThreads[ith]->stats[thrStatIdx[ith]]; - if (stat.firstChip == chid) { - thrStatIdx[ith]++; - chid += stat.nChips; // next chip to look - const auto clbeg = mThreads[ith]->compClusters.begin() + stat.firstClus; - compClus->insert(compClus->end(), clbeg, clbeg + stat.nClus); - if (patterns) { - const auto ptbeg = mThreads[ith]->patterns.begin() + stat.firstPatt; - patterns->insert(patterns->end(), ptbeg, ptbeg + stat.nPatt); - } - if (labelsCl) { - labelsCl->mergeAtBack(mThreads[ith]->labels, stat.firstClus, stat.nClus); - } - } - } - } - for (int ith = 0; ith < nThreads; ith++) { - mThreads[ith]->patterns.clear(); - mThreads[ith]->compClusters.clear(); - mThreads[ith]->labels.clear(); - mThreads[ith]->stats.clear(); - } -#ifdef _PERFORM_TIMING_ - mTimerMerge.Stop(); -#endif - } else { - mThreads[0]->stats.clear(); - } - rof.setNEntries(compClus->size() - rof.getFirstEntry()); // update - } while (autoDecode); - reader.setDecodeNextAuto(autoDecode); // restore setting -#ifdef _PERFORM_TIMING_ - mTimer.Stop(); -#endif -} - -//__________________________________________________ -void Clusterer::ClustererThread::process(uint16_t chip, uint16_t nChips, CompClusCont* compClusPtr, PatternCont* patternsPtr, - const ConstMCTruth* labelsDigPtr, MCTruth* labelsClPtr, const ROFRecord& rofPtr) -{ - if (stats.empty() || stats.back().firstChip + stats.back().nChips != chip) { // there is a jump, register new block - stats.emplace_back(ThreadStat{chip, 0, uint32_t(compClusPtr->size()), patternsPtr ? uint32_t(patternsPtr->size()) : 0, 0, 0}); - } - for (int ic = 0; ic < nChips; ic++) { - auto* curChipData = parent->mFiredChipsPtr[chip + ic]; - auto chipID = curChipData->getChipID(); - if (parent->mMaxBCSeparationToMask > 0) { // mask pixels fired from the previous ROF - const auto& chipInPrevROF = parent->mChipsOld[chipID]; - if (std::abs(rofPtr.getBCData().differenceInBC(chipInPrevROF.getInteractionRecord())) < parent->mMaxBCSeparationToMask) { - parent->mMaxRowColDiffToMask != 0 ? curChipData->maskFiredInSample(parent->mChipsOld[chipID], parent->mMaxRowColDiffToMask) : curChipData->maskFiredInSample(parent->mChipsOld[chipID]); - } - } - auto validPixID = curChipData->getFirstUnmasked(); - auto npix = curChipData->getData().size(); - if (validPixID < npix) { // chip data may have all of its pixels masked! - auto valp = validPixID++; - if (validPixID == npix) { // special case of a single pixel fired on the chip - finishChipSingleHitFast(valp, curChipData, compClusPtr, patternsPtr, labelsDigPtr, labelsClPtr); - } else { - initChip(curChipData, valp); - for (; validPixID < npix; validPixID++) { - if (!curChipData->getData()[validPixID].isMasked()) { - updateChip(curChipData, validPixID); - } - } - finishChip(curChipData, compClusPtr, patternsPtr, labelsDigPtr, labelsClPtr); - } - } - if (parent->mMaxBCSeparationToMask > 0) { // current chip data will be used in the next ROF to mask overflow pixels - parent->mChipsOld[chipID].swap(*curChipData); - } - } - auto& currStat = stats.back(); - currStat.nChips += nChips; - currStat.nClus = compClusPtr->size() - currStat.firstClus; - currStat.nPatt = patternsPtr ? (patternsPtr->size() - currStat.firstPatt) : 0; -} - -//__________________________________________________ -void Clusterer::ClustererThread::finishChip(ChipPixelData* curChipData, CompClusCont* compClusPtr, - PatternCont* patternsPtr, const ConstMCTruth* labelsDigPtr, MCTruth* labelsClusPtr) -{ - const auto& pixData = curChipData->getData(); - for (size_t i1 = 0; i1 < preClusterHeads.size(); ++i1) { - auto ci = preClusterIndices[i1]; - if (ci < 0) { - continue; - } - BBox bbox(curChipData->getChipID()); - int nlab = 0; - int next = preClusterHeads[i1]; - pixArrBuff.clear(); - while (next >= 0) { - const auto& pixEntry = pixels[next]; - const auto pix = pixData[pixEntry.second]; - pixArrBuff.push_back(pix); // needed for cluster topology - bbox.adjust(pix.getRowDirect(), pix.getCol()); - if (labelsClusPtr) { - if (parent->mSquashingDepth) { // the MCtruth for this pixel is stored in chip data: due to squashing we lose contiguity - fetchMCLabels(curChipData->getOrderedPixId(pixEntry.second), labelsDigPtr, nlab); - } else { // the MCtruth for this pixel is at curChipData->startID+pixEntry.second - fetchMCLabels(pixEntry.second + curChipData->getStartID(), labelsDigPtr, nlab); - } - } - next = pixEntry.first; - } - preClusterIndices[i1] = -1; - for (size_t i2 = i1 + 1; i2 < preClusterHeads.size(); ++i2) { - if (preClusterIndices[i2] != ci) { - continue; - } - next = preClusterHeads[i2]; - while (next >= 0) { - const auto& pixEntry = pixels[next]; - const auto pix = pixData[pixEntry.second]; // PixelData - pixArrBuff.push_back(pix); // needed for cluster topology - bbox.adjust(pix.getRowDirect(), pix.getCol()); - if (labelsClusPtr != nullptr) { - if (parent->mSquashingDepth) { // the MCtruth for this pixel is stored in chip data: due to squashing we lose contiguity - fetchMCLabels(curChipData->getOrderedPixId(pixEntry.second), labelsDigPtr, nlab); - } else { // the MCtruth for this pixel is at curChipData->startID+pixEntry.second - fetchMCLabels(pixEntry.second + curChipData->getStartID(), labelsDigPtr, nlab); - } - } - next = pixEntry.first; - } - preClusterIndices[i2] = -1; - } - if (bbox.isAcceptableSize()) { - parent->streamCluster(pixArrBuff, &labelsBuff, bbox, parent->mPattIdConverter, compClusPtr, patternsPtr, labelsClusPtr, nlab, constants::detID::isDetITS3(curChipData->getChipID())); - } else { - auto warnLeft = MaxHugeClusWarn - parent->mNHugeClus; - if (warnLeft > 0) { - LOGP(warn, "Splitting a huge cluster: chipID {}, rows {}:{} cols {}:{}{}", bbox.chipID, bbox.rowMin, bbox.rowMax, bbox.colMin, bbox.colMax, - warnLeft == 1 ? " (Further warnings will be muted)" : ""); -#ifdef WITH_OPENMP -#pragma omp critical -#endif - { - parent->mNHugeClus++; - } - } - BBox bboxT(bbox); // truncated box - std::vector pixbuf; - do { - bboxT.rowMin = bbox.rowMin; - bboxT.colMax = std::min(bbox.colMax, uint16_t(bboxT.colMin + o2::itsmft::ClusterPattern::MaxColSpan - 1)); - do { // Select a subset of pixels fitting the reduced bounding box - bboxT.rowMax = std::min(bbox.rowMax, uint16_t(bboxT.rowMin + o2::itsmft::ClusterPattern::MaxRowSpan - 1)); - for (const auto& pix : pixArrBuff) { - if (bboxT.isInside(pix.getRowDirect(), pix.getCol())) { - pixbuf.push_back(pix); - } - } - if (!pixbuf.empty()) { // Stream a piece of cluster only if the reduced bounding box is not empty - parent->streamCluster(pixbuf, &labelsBuff, bboxT, parent->mPattIdConverter, compClusPtr, patternsPtr, labelsClusPtr, nlab, constants::detID::isDetITS3(curChipData->getChipID()), true); - pixbuf.clear(); - } - bboxT.rowMin = bboxT.rowMax + 1; - } while (bboxT.rowMin < bbox.rowMax); - bboxT.colMin = bboxT.colMax + 1; - } while (bboxT.colMin < bbox.colMax); - } - } -} - -//__________________________________________________ -void Clusterer::ClustererThread::finishChipSingleHitFast(uint32_t hit, ChipPixelData* curChipData, CompClusCont* compClusPtr, - PatternCont* patternsPtr, const ConstMCTruth* labelsDigPtr, MCTruth* labelsClusPtr) -{ - auto pix = curChipData->getData()[hit]; - uint16_t row = pix.getRowDirect(), col = pix.getCol(); - - if (labelsClusPtr) { // MC labels were requested - int nlab = 0; - fetchMCLabels(curChipData->getStartID() + hit, labelsDigPtr, nlab); - auto cnt = compClusPtr->size(); - for (int i = nlab; i--;) { - labelsClusPtr->addElement(cnt, labelsBuff[i]); - } - } - - auto ib = constants::detID::isDetITS3(curChipData->getChipID()); - - // add to compact clusters, which must be always filled - unsigned char patt[ClusterPattern::MaxPatternBytes]{0x1 << (7 - (0 % 8))}; // unrolled 1 hit version of full loop in finishChip - uint16_t pattID = (parent->mPattIdConverter.size(ib) == 0) ? CompCluster::InvalidPatternID : parent->mPattIdConverter.findGroupID(1, 1, ib, patt); - if ((pattID == CompCluster::InvalidPatternID || parent->mPattIdConverter.isGroup(pattID, ib)) && patternsPtr) { - patternsPtr->emplace_back(1); // rowspan - patternsPtr->emplace_back(1); // colspan - patternsPtr->insert(patternsPtr->end(), std::begin(patt), std::begin(patt) + 1); - } - compClusPtr->emplace_back(row, col, pattID, curChipData->getChipID()); -} - -//__________________________________________________ -Clusterer::Clusterer() -{ -#ifdef _PERFORM_TIMING_ - mTimer.Stop(); - mTimer.Reset(); - mTimerMerge.Stop(); - mTimerMerge.Reset(); -#endif -} - -//__________________________________________________ -void Clusterer::ClustererThread::initChip(const ChipPixelData* curChipData, uint32_t first) -{ - // init chip with the 1st unmasked pixel (entry "from" in the mChipData) - size = itsmft::SegmentationAlpide::NRows + 2; - int chipId = curChipData->getChipID(); - if (its3::constants::detID::isDetITS3(chipId)) { - size = its3::SegmentationMosaix::NRows + 2; - } - - delete[] column1; - delete[] column2; - column1 = new int[size]; - column2 = new int[size]; - column1[0] = column1[size - 1] = -1; - column2[0] = column2[size - 1] = -1; - // init chip with the 1st unmasked pixel (entry "from" in the mChipData) - prev = column1 + 1; - curr = column2 + 1; - resetColumn(curr); - - pixels.clear(); - preClusterHeads.clear(); - preClusterIndices.clear(); - auto pix = curChipData->getData()[first]; - currCol = pix.getCol(); - curr[pix.getRowDirect()] = 0; // can use getRowDirect since the pixel is not masked - // start the first pre-cluster - preClusterHeads.push_back(0); - preClusterIndices.push_back(0); - pixels.emplace_back(-1, first); // id of current pixel - noLeftCol = true; // flag that there is no column on the left to check yet -} - -//__________________________________________________ -void Clusterer::ClustererThread::updateChip(const ChipPixelData* curChipData, uint32_t ip) -{ - const auto pix = curChipData->getData()[ip]; - uint16_t row = pix.getRowDirect(); // can use getRowDirect since the pixel is not masked - if (currCol != pix.getCol()) { // switch the buffers - swapColumnBuffers(); - resetColumn(curr); - noLeftCol = false; - if (pix.getCol() > currCol + 1) { - // no connection with previous column, this pixel cannot belong to any of the - // existing preclusters, create a new precluster and flag to check only the row above for next pixels of this column - currCol = pix.getCol(); - addNewPrecluster(ip, row); - noLeftCol = true; - return; - } - currCol = pix.getCol(); - } - - bool orphan = true; - - if (noLeftCol) { // check only the row above - if (curr[row - 1] >= 0) { - expandPreCluster(ip, row, curr[row - 1]); // attach to the precluster of the previous row - return; - } - } else { -#ifdef _ALLOW_DIAGONAL_ALPIDE_CLUSTERS_ - int neighbours[]{curr[row - 1], prev[row], prev[row + 1], prev[row - 1]}; -#else - int neighbours[]{curr[row - 1], prev[row]}; -#endif - for (auto pci : neighbours) { - if (pci < 0) { - continue; - } - if (orphan) { - expandPreCluster(ip, row, pci); // attach to the adjascent precluster - orphan = false; - continue; - } - // reassign precluster index to smallest one - if (preClusterIndices[pci] < preClusterIndices[curr[row]]) { - preClusterIndices[curr[row]] = preClusterIndices[pci]; - } else { - preClusterIndices[pci] = preClusterIndices[curr[row]]; - } - } - } - if (orphan) { - addNewPrecluster(ip, row); // start new precluster - } -} - -//__________________________________________________ -void Clusterer::ClustererThread::fetchMCLabels(int digID, const ConstMCTruth* labelsDig, int& nfilled) -{ - // transfer MC labels to cluster - if (nfilled >= MaxLabels) { - return; - } - const auto& lbls = labelsDig->getLabels(digID); - for (int i = lbls.size(); i--;) { - int ic = nfilled; - for (; ic--;) { // check if the label is already present - if (labelsBuff[ic] == lbls[i]) { - return; // label is found, do nothing - } - } - labelsBuff[nfilled++] = lbls[i]; - if (nfilled >= MaxLabels) { - break; - } - } - // -} - -//__________________________________________________ -void Clusterer::clear() -{ - // reset -#ifdef _PERFORM_TIMING_ - mTimer.Stop(); - mTimer.Reset(); - mTimerMerge.Stop(); - mTimerMerge.Reset(); -#endif -} - -//__________________________________________________ -void Clusterer::print() const -{ - LOGP(info, "Clusterizer has {} chips registered", mChips.size()); - LOGP(info, "Clusterizer squashes overflow pixels separated by {} BC and <= {} in row/col seeking down to {} neighbour ROFs", mMaxBCSeparationToSquash, mMaxRowColDiffToMask, mSquashingDepth); - LOGP(info, "Clusterizer masks overflow pixels separated by < {} BC and <= {} in row/col", mMaxBCSeparationToMask, mMaxRowColDiffToMask); - -#ifdef _PERFORM_TIMING_ - auto& tmr = const_cast(mTimer); // ugly but this is what root does internally - auto& tmrm = const_cast(mTimerMerge); - LOG(info) << "Inclusive clusterization timing (w/o disk IO): Cpu: " << tmr.CpuTime() - << " Real: " << tmr.RealTime() << " s in " << tmr.Counter() << " slots"; - LOG(info) << "Threads output merging timing : Cpu: " << tmrm.CpuTime() - << " Real: " << tmrm.RealTime() << " s in " << tmrm.Counter() << " slots"; - -#endif -} -} // namespace o2::its3 +template class o2::itsmft::ClustererT; diff --git a/Detectors/Upgrades/ITS3/reconstruction/src/ITS3ReconstructionLinkDef.h b/Detectors/Upgrades/ITS3/reconstruction/src/ITS3ReconstructionLinkDef.h index 2ebd89970d9a1..bfdd7c3815701 100644 --- a/Detectors/Upgrades/ITS3/reconstruction/src/ITS3ReconstructionLinkDef.h +++ b/Detectors/Upgrades/ITS3/reconstruction/src/ITS3ReconstructionLinkDef.h @@ -15,7 +15,6 @@ #pragma link off all classes; #pragma link off all functions; -#pragma link C++ class o2::its3::Clusterer + ; #pragma link C++ class o2::its3::TopologyDictionaryData + ; #pragma link C++ class o2::its3::TopologyDictionary + ; #pragma link C++ class o2::its3::BuildTopologyDictionary + ; diff --git a/Detectors/Upgrades/ITS3/workflow/include/ITS3Workflow/ClustererSpec.h b/Detectors/Upgrades/ITS3/workflow/include/ITS3Workflow/ClustererSpec.h index e6732551c60e2..b6ed982d79882 100644 --- a/Detectors/Upgrades/ITS3/workflow/include/ITS3Workflow/ClustererSpec.h +++ b/Detectors/Upgrades/ITS3/workflow/include/ITS3Workflow/ClustererSpec.h @@ -14,7 +14,6 @@ #ifndef O2_ITS3_CLUSTERERDPL #define O2_ITS3_CLUSTERERDPL -#include #include "DetectorsBase/GRPGeomHelper.h" #include "ITS3Reconstruction/Clusterer.h" #include "Framework/DataProcessorSpec.h" From be0442141326d8ac424db5a6e80ec12807de03a6 Mon Sep 17 00:00:00 2001 From: Felix Schlepper Date: Wed, 22 Apr 2026 16:42:57 +0200 Subject: [PATCH 2/3] ITS3: adapt staggering Signed-off-by: Felix Schlepper --- .../include/ITStracking/TrackingInterface.h | 5 +- .../ITS/tracking/src/TrackingInterface.cxx | 7 +- .../include/ITS3Reconstruction/IOUtils.h | 1 + .../ITS3Reconstruction/TrackingInterface.h | 2 +- .../ITS3/reconstruction/src/IOUtils.cxx | 32 +- .../reconstruction/src/TrackingInterface.cxx | 28 +- .../include/ITS3Simulation/Digitizer.h | 17 +- .../ITS3/simulation/src/Digitizer.cxx | 76 +++-- .../include/ITS3Workflow/ClustererSpec.h | 4 +- .../include/ITS3Workflow/DigitReaderSpec.h | 69 ++-- .../ITS3/workflow/src/ClustererSpec.cxx | 226 ++++++++++---- .../ITS3/workflow/src/DigitReaderSpec.cxx | 116 ++++--- .../ITS3/workflow/src/DigitWriterSpec.cxx | 72 +++-- .../ITS3/workflow/src/RecoWorkflow.cxx | 12 +- .../ITS3/workflow/src/TrackerSpec.cxx | 22 +- .../src/ITS3DigitizerSpec.cxx | 295 +++++++++++------- 16 files changed, 586 insertions(+), 398 deletions(-) diff --git a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/TrackingInterface.h b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/TrackingInterface.h index ac4b99a0a8cd8..16f75ba685797 100644 --- a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/TrackingInterface.h +++ b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/TrackingInterface.h @@ -34,6 +34,7 @@ namespace o2::its { class ITSTrackingInterface { + public: static constexpr int NLayers{7}; using VertexerN = Vertexer; using VertexerTraitsN = VertexerTraits; @@ -41,7 +42,6 @@ class ITSTrackingInterface using TrackerTraitsN = TrackerTraits; using TimeFrameN = TimeFrame; - public: ITSTrackingInterface(bool isMC, bool doStag, int trgType, @@ -80,6 +80,7 @@ class ITSTrackingInterface TimeFrameN* mTimeFrame = nullptr; protected: + virtual void requestTopologyDictionary(framework::ProcessingContext& pc); virtual void loadROF(gsl::span& trackROFspan, gsl::span clusters, gsl::span::iterator& pattIt, @@ -98,7 +99,7 @@ class ITSTrackingInterface const o2::itsmft::TopologyDictionary* mDict = nullptr; std::unique_ptr mTracker = nullptr; std::unique_ptr mVertexer = nullptr; - const o2::dataformats::MeanVertexObject* mMeanVertex; + const o2::dataformats::MeanVertexObject* mMeanVertex{}; std::shared_ptr mMemoryPool; std::shared_ptr mTaskArena; }; diff --git a/Detectors/ITSMFT/ITS/tracking/src/TrackingInterface.cxx b/Detectors/ITSMFT/ITS/tracking/src/TrackingInterface.cxx index fcd9024a74709..ce345177a1c94 100644 --- a/Detectors/ITSMFT/ITS/tracking/src/TrackingInterface.cxx +++ b/Detectors/ITSMFT/ITS/tracking/src/TrackingInterface.cxx @@ -362,7 +362,7 @@ void ITSTrackingInterface::updateTimeDependentParams(framework::ProcessingContex } if (!initOnceDone) { // this params need to be queried only once initOnceDone = true; - pc.inputs().get("itscldict"); // just to trigger the finaliseCCDB + requestTopologyDictionary(pc); pc.inputs().get*>("itsalppar"); if (pc.inputs().getPos("itsTGeo") >= 0) { pc.inputs().get("itsTGeo"); @@ -474,6 +474,11 @@ void ITSTrackingInterface::setTraitsFromProvider(VertexerTraitsN* vertexerTraits mVertexer->setMemoryPool(mMemoryPool); } +void ITSTrackingInterface::requestTopologyDictionary(framework::ProcessingContext& pc) +{ + pc.inputs().get("itscldict"); // just to trigger the finaliseCCDB +} + void ITSTrackingInterface::loadROF(gsl::span& trackROFspan, gsl::span clusters, gsl::span::iterator& pattIt, diff --git a/Detectors/Upgrades/ITS3/reconstruction/include/ITS3Reconstruction/IOUtils.h b/Detectors/Upgrades/ITS3/reconstruction/include/ITS3Reconstruction/IOUtils.h index d82cd26e63c56..05ec269a89219 100644 --- a/Detectors/Upgrades/ITS3/reconstruction/include/ITS3Reconstruction/IOUtils.h +++ b/Detectors/Upgrades/ITS3/reconstruction/include/ITS3Reconstruction/IOUtils.h @@ -83,6 +83,7 @@ int loadROFrameDataITS3(its::TimeFrame<7>* tf, gsl::span clusters, gsl::span::iterator& pattIt, const o2::its3::TopologyDictionary* dict, + int layer, const dataformats::MCTruthContainer* mcLabels = nullptr); } // namespace its3::ioutils } // namespace o2 diff --git a/Detectors/Upgrades/ITS3/reconstruction/include/ITS3Reconstruction/TrackingInterface.h b/Detectors/Upgrades/ITS3/reconstruction/include/ITS3Reconstruction/TrackingInterface.h index 3b743c59524d2..27cfa93bac70f 100644 --- a/Detectors/Upgrades/ITS3/reconstruction/include/ITS3Reconstruction/TrackingInterface.h +++ b/Detectors/Upgrades/ITS3/reconstruction/include/ITS3Reconstruction/TrackingInterface.h @@ -24,10 +24,10 @@ class ITS3TrackingInterface final : public its::ITSTrackingInterface using its::ITSTrackingInterface::ITSTrackingInterface; void setClusterDictionary(const o2::its3::TopologyDictionary* d) { mDict = d; } - void updateTimeDependentParams(framework::ProcessingContext& pc) final; void finaliseCCDB(framework::ConcreteDataMatcher& matcher, void* obj) final; protected: + void requestTopologyDictionary(framework::ProcessingContext& pc) final; void loadROF(gsl::span& trackROFspan, gsl::span clusters, gsl::span::iterator& pattIt, diff --git a/Detectors/Upgrades/ITS3/reconstruction/src/IOUtils.cxx b/Detectors/Upgrades/ITS3/reconstruction/src/IOUtils.cxx index 92e36cd2a4b84..146875804c625 100644 --- a/Detectors/Upgrades/ITS3/reconstruction/src/IOUtils.cxx +++ b/Detectors/Upgrades/ITS3/reconstruction/src/IOUtils.cxx @@ -61,13 +61,21 @@ int loadROFrameDataITS3(its::TimeFrame<7>* tf, gsl::span clusters, gsl::span::iterator& pattIt, const its3::TopologyDictionary* dict, + int layer, const dataformats::MCTruthContainer* mcLabels) { auto geom = its::GeometryTGeo::Instance(); geom->fillMatrixCache(o2::math_utils::bit2Mask(o2::math_utils::TransformType::T2L, o2::math_utils::TransformType::L2G)); - // tf->resetROFrameData(rofs.size()); // FIXME - // tf->prepareROFrameData(rofs, clusters); FIXME + tf->resetROFrameData(layer); + tf->prepareROFrameData(clusters, layer); + + // check for missing/empty/unset rofs + // the code requires consistent monotonically increasing input without gaps + const auto& timing = tf->getROFOverlapTableView().getLayer(layer); + if (timing.mNROFsTF != rofs.size()) { + LOGP(fatal, "Received inconsistent number of rofs on layer:{} expected:{} received:{}", layer, timing.mNROFsTF, rofs.size()); + } its::bounded_vector clusterSizeVec(clusters.size(), tf->getMemoryPool().get()); @@ -110,22 +118,22 @@ int loadROFrameDataITS3(its::TimeFrame<7>* tf, tf->addClusterToLayer(layer, gloXYZ.x(), gloXYZ.y(), gloXYZ.z(), tf->getUnsortedClusters()[layer].size()); tf->addClusterExternalIndexToLayer(layer, clusterId); } - for (unsigned int iL{0}; iL < tf->getUnsortedClusters().size(); ++iL) { - tf->mROFramesClusters[iL][iRof + 1] = (int)tf->getUnsortedClusters()[iL].size(); - } + tf->mROFramesClusters[layer][iRof + 1] = (int)tf->getUnsortedClusters()[layer].size(); } - // tf->setClusterSize(clusterSizeVec); FIXME + tf->setClusterSize(layer, clusterSizeVec); - for (auto& v : tf->mNTrackletsPerCluster) { - v.resize(tf->getUnsortedClusters()[1].size()); - } - for (auto& v : tf->mNTrackletsPerClusterSum) { - v.resize(tf->getUnsortedClusters()[1].size() + 1); + if (layer == 1) { + for (auto& v : tf->mNTrackletsPerCluster) { + v.resize(tf->getUnsortedClusters()[1].size()); + } + for (auto& v : tf->mNTrackletsPerClusterSum) { + v.resize(tf->getUnsortedClusters()[1].size() + 1); + } } if (mcLabels != nullptr) { - // tf->mClusterLabels = mcLabels; // FIXME + tf->mClusterLabels[layer] = mcLabels; } return 0; } diff --git a/Detectors/Upgrades/ITS3/reconstruction/src/TrackingInterface.cxx b/Detectors/Upgrades/ITS3/reconstruction/src/TrackingInterface.cxx index 9fe6f3735a845..90aee2a09d2d7 100644 --- a/Detectors/Upgrades/ITS3/reconstruction/src/TrackingInterface.cxx +++ b/Detectors/Upgrades/ITS3/reconstruction/src/TrackingInterface.cxx @@ -12,6 +12,7 @@ #include "ITS3Reconstruction/TrackingInterface.h" #include "ITS3Reconstruction/IOUtils.h" #include "ITSBase/GeometryTGeo.h" +#include "ITStracking/FastMultEstConfig.h" #include "ITStracking/TrackingConfigParam.h" #include "DataFormatsITSMFT/DPLAlpideParam.h" #include "DetectorsBase/GRPGeomHelper.h" @@ -20,30 +21,9 @@ namespace o2::its3 { -void ITS3TrackingInterface::updateTimeDependentParams(framework::ProcessingContext& pc) +void ITS3TrackingInterface::requestTopologyDictionary(framework::ProcessingContext& pc) { - o2::base::GRPGeomHelper::instance().checkUpdates(pc); - static bool initOnceDone = false; - if (!initOnceDone) { // this params need to be queried only once - initOnceDone = true; - pc.inputs().get("cldict"); // just to trigger the finaliseCCDB - pc.inputs().get*>("alppar"); - if (pc.inputs().getPos("itsTGeo") >= 0) { - pc.inputs().get("itsTGeo"); - } - auto geom = its::GeometryTGeo::Instance(); - geom->fillMatrixCache(o2::math_utils::bit2Mask(o2::math_utils::TransformType::T2L, o2::math_utils::TransformType::T2GRot, o2::math_utils::TransformType::T2G)); - initialise(); - if (pc.services().get().inputTimesliceId == 0) { // print settings only for the 1st pipeling - o2::its::VertexerParamConfig::Instance().printKeyValues(); - o2::its::TrackerParamConfig::Instance().printKeyValues(); - const auto& trParams = getTracker()->getParameters(); - for (size_t it = 0; it < trParams.size(); it++) { - const auto& par = trParams[it]; - LOGP(info, "recoIter#{} : {}", it, par.asString()); - } - } - } + pc.inputs().get("itscldict"); // just to trigger the finaliseCCDB } void ITS3TrackingInterface::finaliseCCDB(framework::ConcreteDataMatcher& matcher, void* obj) @@ -80,7 +60,7 @@ void ITS3TrackingInterface::loadROF(gsl::span& trackROF int layer, const dataformats::MCTruthContainer* mcLabels) { - // ioutils::loadROFrameDataITS3(mTimeFrame, trackROFspan, clusters, pattIt, mDict, mcLabels); + ioutils::loadROFrameDataITS3(mTimeFrame, trackROFspan, clusters, pattIt, mDict, layer, mcLabels); } } // namespace o2::its3 diff --git a/Detectors/Upgrades/ITS3/simulation/include/ITS3Simulation/Digitizer.h b/Detectors/Upgrades/ITS3/simulation/include/ITS3Simulation/Digitizer.h index 866973083983b..78bb9923dae97 100644 --- a/Detectors/Upgrades/ITS3/simulation/include/ITS3Simulation/Digitizer.h +++ b/Detectors/Upgrades/ITS3/simulation/include/ITS3Simulation/Digitizer.h @@ -54,8 +54,8 @@ class Digitizer : public TObject void init(); /// Steer conversion of hits to digits - void process(const std::vector* hits, int evID, int srcID); - void setEventTime(const o2::InteractionTimeRecord& irt); + void process(const std::vector* hits, int evID, int srcID, int layer); + void setEventTime(const o2::InteractionTimeRecord& irt, int layer); double getEndTimeOfROFMax() const { ///< return the time corresponding to end of the last reserved ROFrame : mROFrameMax @@ -64,7 +64,7 @@ class Digitizer : public TObject void setContinuous(bool v) { mParams.setContinuous(v); } bool isContinuous() const { return mParams.isContinuous(); } - void fillOutputContainer(uint32_t maxFrame = 0xffffffff); + void fillOutputContainer(uint32_t maxFrame = 0xffffffff, int layer = -1); // provide the common itsmft::GeometryTGeo to access matrices and segmentation void setGeometry(const o2::its::GeometryTGeo* gm) { mGeometry = gm; } @@ -76,13 +76,19 @@ class Digitizer : public TObject mEventROFrameMin = 0xffffffff; mEventROFrameMax = 0; } + void resetROFrameBounds() + { + mROFrameMin = 0; + mROFrameMax = 0; + mNewROFrame = 0; + } void setDeadChannelsMap(const o2::itsmft::NoiseMap* mp) { mDeadChanMap = mp; } private: - void processHit(const o2::itsmft::Hit& hit, uint32_t& maxFr, int evID, int srcID); + void processHit(const o2::itsmft::Hit& hit, uint32_t& maxFr, int evID, int srcID, int layer); void registerDigits(o2::its3::ChipDigitsContainer& chip, uint32_t roFrame, float tInROF, int nROF, - uint16_t row, uint16_t col, int nEle, o2::MCCompLabel& lbl); + uint16_t row, uint16_t col, int nEle, o2::MCCompLabel& lbl, int layer); ExtraDig* getExtraDigBuffer(uint32_t roFrame) { @@ -105,6 +111,7 @@ class Digitizer : public TObject uint32_t mROFrameMin = 0; ///< lowest RO frame of current digits uint32_t mROFrameMax = 0; ///< highest RO frame of current digits uint32_t mNewROFrame = 0; ///< ROFrame corresponding to provided time + bool mIsBeforeFirstRO = false; uint32_t mEventROFrameMin = 0xffffffff; ///< lowest RO frame for processed events (w/o automatic noise ROFs) uint32_t mEventROFrameMax = 0; ///< highest RO frame forfor processed events (w/o automatic noise ROFs) diff --git a/Detectors/Upgrades/ITS3/simulation/src/Digitizer.cxx b/Detectors/Upgrades/ITS3/simulation/src/Digitizer.cxx index 4b0374d925401..86282c2d8c525 100644 --- a/Detectors/Upgrades/ITS3/simulation/src/Digitizer.cxx +++ b/Detectors/Upgrades/ITS3/simulation/src/Digitizer.cxx @@ -25,6 +25,7 @@ #include #include #include +#include #include using o2::itsmft::Hit; @@ -72,18 +73,18 @@ void Digitizer::init() mIRFirstSampledTF = o2::raw::HBFUtils::Instance().getFirstSampledTFIR(); } -void Digitizer::process(const std::vector* hits, int evID, int srcID) +void Digitizer::process(const std::vector* hits, int evID, int srcID, int layer) { // digitize single event, the time must have been set beforehand - LOG(debug) << "Digitizing " << mGeometry->getName() << " hits of entry " << evID << " from source " + LOG(debug) << "Digitizing " << mGeometry->getName() << ":" << layer << " hits of entry " << evID << " from source " << srcID << " at time " << mEventTime << " ROFrame = " << mNewROFrame << ")" << " cont.mode: " << isContinuous() << " Min/Max ROFrames " << mROFrameMin << "/" << mROFrameMax; // is there something to flush ? if (mNewROFrame > mROFrameMin) { - fillOutputContainer(mNewROFrame - 1); // flush out all frame preceding the new one + fillOutputContainer(mNewROFrame - 1, layer); // flush out all frame preceding the new one } int nHits = hits->size(); @@ -94,18 +95,20 @@ void Digitizer::process(const std::vector* hits, int evID, int srcI [hits](auto lhs, auto rhs) { return (*hits)[lhs].GetDetectorID() < (*hits)[rhs].GetDetectorID(); }); - for (int i : hitIdx) { - processHit((*hits)[i], mROFrameMax, evID, srcID); + for (int i : hitIdx | std::views::filter([&](int idx) { + return mGeometry->getLayer((*hits)[idx].GetDetectorID()) == layer; + })) { + processHit((*hits)[i], mROFrameMax, evID, srcID, layer); } // in the triggered mode store digits after every MC event // TODO: in the real triggered mode this will not be needed, this is actually for the // single event processing only if (!mParams.isContinuous()) { - fillOutputContainer(mROFrameMax); + fillOutputContainer(mROFrameMax, layer); } } -void Digitizer::setEventTime(const o2::InteractionTimeRecord& irt) +void Digitizer::setEventTime(const o2::InteractionTimeRecord& irt, int layer) { // assign event time in ns mEventTime = irt; @@ -120,9 +123,20 @@ void Digitizer::setEventTime(const o2::InteractionTimeRecord& irt) if (mCollisionTimeWrtROF < 0 && nbc > 0) { nbc--; } - mNewROFrame = nbc / mParams.getROFrameLengthInBC(); + // we might get interactions to digitize from before + // the first sampled IR + if (nbc < 0) { + mNewROFrame = 0; + // this event is before the first RO + mIsBeforeFirstRO = true; + } else { + mNewROFrame = nbc / mParams.getROFrameLengthInBC(layer); + mIsBeforeFirstRO = false; + } + LOG(debug) << " NewROFrame " << mNewROFrame << " nbc " << nbc; + // in continuous mode depends on starts of periodic readout frame - mCollisionTimeWrtROF += (nbc % mParams.getROFrameLengthInBC()) * o2::constants::lhc::LHCBunchSpacingNS; + mCollisionTimeWrtROF += (nbc % mParams.getROFrameLengthInBC(layer)) * o2::constants::lhc::LHCBunchSpacingNS; } else { mNewROFrame = 0; } @@ -137,17 +151,14 @@ void Digitizer::setEventTime(const o2::InteractionTimeRecord& irt) } } -void Digitizer::fillOutputContainer(uint32_t frameLast) +void Digitizer::fillOutputContainer(uint32_t frameLast, int layer) { // fill output with digits from min.cached up to requested frame, generating the noise beforehand - if (frameLast > mROFrameMax) { - frameLast = mROFrameMax; - } + frameLast = std::min(frameLast, mROFrameMax); // make sure all buffers for extra digits are created up to the maxFrame getExtraDigBuffer(mROFrameMax); - LOG(info) << "Filling " << mGeometry->getName() << " digits output for RO frames " << mROFrameMin << ":" - << frameLast; + LOG(info) << "Filling IT3 digits output on layer " << layer << " for RO frames " << mROFrameMin << ":" << frameLast; o2::itsmft::ROFRecord rcROF; @@ -159,6 +170,9 @@ void Digitizer::fillOutputContainer(uint32_t frameLast) auto& extra = *(mExtraBuff.front().get()); for (size_t iChip{0}; iChip < mChips.size(); ++iChip) { auto& chip = mChips[iChip]; + if (chip.isDisabled() || mGeometry->getLayer(chip.getChipIndex()) != layer) { + continue; + } chip.addNoise(mROFrameMin, mROFrameMin, &mParams); auto& buffer = chip.getPreDigits(); if (buffer.empty()) { @@ -188,7 +202,7 @@ void Digitizer::fillOutputContainer(uint32_t frameLast) // finalize ROF record rcROF.setNEntries(mDigits->size() - rcROF.getFirstEntry()); // number of digits if (isContinuous()) { - rcROF.getBCData().setFromLong(mIRFirstSampledTF.toLong() + mROFrameMin * mParams.getROFrameLengthInBC()); + rcROF.getBCData().setFromLong(mIRFirstSampledTF.toLong() + mROFrameMin * mParams.getROFrameLengthInBC(layer)); } else { rcROF.getBCData() = mEventTime; // RS TODO do we need to add trigger delay? } @@ -202,7 +216,7 @@ void Digitizer::fillOutputContainer(uint32_t frameLast) } } -void Digitizer::processHit(const o2::itsmft::Hit& hit, uint32_t& maxFr, int evID, int srcID) +void Digitizer::processHit(const o2::itsmft::Hit& hit, uint32_t& maxFr, int evID, int srcID, int layer) { // convert single hit to digits auto chipID = hit.GetDetectorID(); @@ -223,15 +237,20 @@ void Digitizer::processHit(const o2::itsmft::Hit& hit, uint32_t& maxFr, int evID if (isContinuous()) { timeInROF += mCollisionTimeWrtROF; } + if (mIsBeforeFirstRO && timeInROF < 0) { + // disregard this hit because it comes from an event before readout starts and it does not effect this RO + return; + } + // calculate RO Frame for this hit if (timeInROF < 0) { timeInROF = 0.; } float tTot = mParams.getSignalShape().getMaxDuration(); // frame of the hit signal start wrt event ROFrame - int roFrameRel = int(timeInROF * mParams.getROFrameLengthInv()); + int roFrameRel = int(timeInROF * mParams.getROFrameLengthInv(layer)); // frame of the hit signal end wrt event ROFrame: in the triggered mode we read just 1 frame - uint32_t roFrameRelMax = mParams.isContinuous() ? (timeInROF + tTot) * mParams.getROFrameLengthInv() : roFrameRel; + uint32_t roFrameRelMax = mParams.isContinuous() ? (timeInROF + tTot) * mParams.getROFrameLengthInv(layer) : roFrameRel; int nFrames = roFrameRelMax + 1 - roFrameRel; uint32_t roFrameMax = mNewROFrame + roFrameRelMax; if (roFrameMax > maxFr) { @@ -240,7 +259,6 @@ void Digitizer::processHit(const o2::itsmft::Hit& hit, uint32_t& maxFr, int evID // here we start stepping in the depth of the sensor to generate charge diffision int detID{hit.GetDetectorID()}; - int layer = mGeometry->getLayer(detID); const auto& matrix = mGeometry->getMatrixL2G(detID); int nSteps = chip.isIB() ? mParams.getIBNSimSteps() : mParams.getNSimSteps(); float nStepsInv = chip.isIB() ? mParams.getIBNSimStepsInv() : mParams.getNSimStepsInv(); @@ -404,34 +422,30 @@ void Digitizer::processHit(const o2::itsmft::Hit& hit, uint32_t& maxFr, int evID continue; } uint16_t colIS = icol + colS; - registerDigits(chip, roFrameAbs, timeInROF, nFrames, rowIS, colIS, nEle, lbl); + registerDigits(chip, roFrameAbs, timeInROF, nFrames, rowIS, colIS, nEle, lbl, layer); } } } void Digitizer::registerDigits(o2::its3::ChipDigitsContainer& chip, uint32_t roFrame, float tInROF, int nROF, - uint16_t row, uint16_t col, int nEle, o2::MCCompLabel& lbl) + uint16_t row, uint16_t col, int nEle, o2::MCCompLabel& lbl, int layer) { // Register digits for given pixel, accounting for the possible signal contribution to // multiple ROFrame. The signal starts at time tInROF wrt the start of provided roFrame // In every ROFrame we check the collected signal during strobe - float tStrobe = mParams.getStrobeDelay() - tInROF; // strobe start wrt signal start + float tStrobe = mParams.getStrobeDelay(layer) - tInROF; // strobe start wrt signal start for (int i = 0; i < nROF; i++) { uint32_t roFr = roFrame + i; - int nEleROF = mParams.getSignalShape().getCollectedCharge(nEle, tStrobe, tStrobe + mParams.getStrobeLength()); - tStrobe += mParams.getROFrameLength(); // for the next ROF + int nEleROF = mParams.getSignalShape().getCollectedCharge(nEle, tStrobe, tStrobe + mParams.getStrobeLength(layer)); + tStrobe += mParams.getROFrameLength(layer); // for the next ROF // discard too small contributions, they have no chance to produce a digit if (nEleROF < (chip.isIB() ? mParams.getIBChargeThreshold() : mParams.getChargeThreshold())) { continue; } - if (roFr > mEventROFrameMax) { - mEventROFrameMax = roFr; - } - if (roFr < mEventROFrameMin) { - mEventROFrameMin = roFr; - } + mEventROFrameMax = std::max(roFr, mEventROFrameMax); + mEventROFrameMin = std::min(roFr, mEventROFrameMin); auto key = chip.getOrderingKey(roFr, row, col); PreDigit* pd = chip.findDigit(key); if (pd == nullptr) { diff --git a/Detectors/Upgrades/ITS3/workflow/include/ITS3Workflow/ClustererSpec.h b/Detectors/Upgrades/ITS3/workflow/include/ITS3Workflow/ClustererSpec.h index b6ed982d79882..a9213515e4ae9 100644 --- a/Detectors/Upgrades/ITS3/workflow/include/ITS3Workflow/ClustererSpec.h +++ b/Detectors/Upgrades/ITS3/workflow/include/ITS3Workflow/ClustererSpec.h @@ -27,6 +27,8 @@ namespace o2::its3 class ClustererDPL : public Task { public: + static constexpr int NLayers = 7; + ClustererDPL(const std::shared_ptr& gr, bool useMC) : mGGCCDBRequest(gr), mUseMC(useMC) {} void init(InitContext& ic) final; void run(ProcessingContext& pc) final; @@ -38,9 +40,9 @@ class ClustererDPL : public Task std::shared_ptr mGGCCDBRequest; bool mUseMC = true; - int mState = 0; int mNThreads = 1; bool mUseClusterDictionary = true; + std::vector mFilter; std::unique_ptr mClusterer = nullptr; }; diff --git a/Detectors/Upgrades/ITS3/workflow/include/ITS3Workflow/DigitReaderSpec.h b/Detectors/Upgrades/ITS3/workflow/include/ITS3Workflow/DigitReaderSpec.h index 9d0379aa77d78..3a167980d5048 100644 --- a/Detectors/Upgrades/ITS3/workflow/include/ITS3Workflow/DigitReaderSpec.h +++ b/Detectors/Upgrades/ITS3/workflow/include/ITS3Workflow/DigitReaderSpec.h @@ -14,42 +14,53 @@ #ifndef O2_ITS3_DIGITREADER #define O2_ITS3_DIGITREADER -#include "TFile.h" -#include "TTree.h" +#include + +#include +#include + #include "DataFormatsITSMFT/Digit.h" -#include "DataFormatsITSMFT/GBTCalibData.h" #include "DataFormatsITSMFT/ROFRecord.h" #include "Framework/DataProcessorSpec.h" #include "Framework/Task.h" -#include "Headers/DataHeader.h" #include "DataFormatsITSMFT/ROFRecord.h" #include "DetectorsCommonDataFormats/DetID.h" +#include "SimulationDataFormat/IOMCTruthContainerView.h" +#include "SimulationDataFormat/ConstMCTruthContainer.h" using namespace o2::framework; -namespace o2 -{ -namespace its3 +namespace o2::its3 { -class DigitReader : public Task +class ITS3DigitReader final : public Task { public: - DigitReader() = delete; - DigitReader(o2::detectors::DetID id, bool useMC, bool useCalib); - ~DigitReader() override = default; + static constexpr int NLayers = 7; + + ITS3DigitReader(bool useMC, bool useCalib); + ~ITS3DigitReader() override = default; void init(InitContext& ic) final; void run(ProcessingContext& pc) final; protected: void connectTree(const std::string& filename); - std::vector mDigits, *mDigitsPtr = &mDigits; - std::vector mCalib, *mCalibPtr = &mCalib; - std::vector mDigROFRec, *mDigROFRecPtr = &mDigROFRec; - std::vector mDigMC2ROFs, *mDigMC2ROFsPtr = &mDigMC2ROFs; + std::array*, NLayers> mDigits{}; + std::array*, NLayers> mDigROFRec{}; + std::array*, NLayers> mDigMC2ROFs{}; + std::vector> mConstLabels; + std::array mPLabels{}; + + const o2::header::DataOrigin mOrigin = o2::header::gDataOriginIT3; + + std::string getBranchName(const std::string& base, int index) + { + return base + "_" + std::to_string(index); + } - o2::header::DataOrigin mOrigin = o2::header::gDataOriginInvalid; + template + void setBranchAddress(const std::string& base, Ptr& addr, int layer); std::unique_ptr mFile; std::unique_ptr mTree; @@ -57,33 +68,19 @@ class DigitReader : public Task bool mUseMC = true; // use MC truth bool mUseCalib = true; // send calib data - std::string mDetName = ""; - std::string mDetNameLC = ""; - std::string mFileName = ""; + std::string mDetName; + std::string mDetNameLC; + std::string mFileName; std::string mDigTreeName = "o2sim"; - std::string mDigitBranchName = "Digit"; + std::string mDigBranchName = "Digit"; std::string mDigROFBranchName = "DigitROF"; - std::string mCalibBranchName = "Calib"; - - std::string mDigtMCTruthBranchName = "DigitMCTruth"; - std::string mDigtMC2ROFBranchName = "DigitMC2ROF"; -}; - -class ITS3DigitReader : public DigitReader -{ - public: - ITS3DigitReader(bool useMC = true, bool useCalib = false) - : DigitReader(o2::detectors::DetID::IT3, useMC, useCalib) - { - mOrigin = o2::header::gDataOriginIT3; - } + std::string mDigMCTruthBranchName = "DigitMCTruth"; }; /// create a processor spec /// read ITS/MFT Digit data from a root file framework::DataProcessorSpec getITS3DigitReaderSpec(bool useMC = true, bool useCalib = false, std::string defname = "it3digits.root"); -} // namespace its3 -} // namespace o2 +} // namespace o2::its3 #endif /* O2_ITS3_DigitREADER */ diff --git a/Detectors/Upgrades/ITS3/workflow/src/ClustererSpec.cxx b/Detectors/Upgrades/ITS3/workflow/src/ClustererSpec.cxx index 73b5f4650d02d..ae4071aded1d1 100644 --- a/Detectors/Upgrades/ITS3/workflow/src/ClustererSpec.cxx +++ b/Detectors/Upgrades/ITS3/workflow/src/ClustererSpec.cxx @@ -15,6 +15,7 @@ #include "Framework/ConfigParamRegistry.h" #include "Framework/CCDBParamSpec.h" +#include "Framework/InputRecordWalker.h" #include "ITS3Workflow/ClustererSpec.h" #include "DataFormatsITSMFT/Digit.h" #include "ITS3Base/SpecsV2.h" @@ -25,7 +26,6 @@ #include "SimulationDataFormat/MCCompLabel.h" #include "SimulationDataFormat/ConstMCTruthContainer.h" #include "DataFormatsITSMFT/ROFRecord.h" -#include "DataFormatsParameters/GRPObject.h" #include "ITSMFTReconstruction/DigitPixelReader.h" #include "DataFormatsITSMFT/DPLAlpideParam.h" #include "CommonConstants/LHCConstants.h" @@ -42,64 +42,145 @@ void ClustererDPL::init(InitContext& ic) mNThreads = std::max(1, ic.options().get("nthreads")); o2::base::GRPGeomHelper::instance().setRequest(mGGCCDBRequest); mClusterer->setNChips(its3::constants::detID::nChips + itsmft::ChipMappingITS::getNChips(itsmft::ChipMappingITS::MB) + itsmft::ChipMappingITS::getNChips(itsmft::ChipMappingITS::OB)); - mClusterer->print(); - mState = 1; + + // prepare data filter + for (int iLayer = 0; iLayer < NLayers; ++iLayer) { + mFilter.emplace_back("digits", "IT3", "DIGITS", iLayer, Lifetime::Timeframe); + mFilter.emplace_back("ROframe", "IT3", "DIGITSROF", iLayer, Lifetime::Timeframe); + if (mUseMC) { + mFilter.emplace_back("labels", "IT3", "DIGITSMCTR", iLayer, Lifetime::Timeframe); + } + } } void ClustererDPL::run(ProcessingContext& pc) { updateTimeDependentParams(pc); - auto digits = pc.inputs().get>("digits"); - auto rofs = pc.inputs().get>("ROframes"); - - gsl::span mc2rofs; - gsl::span labelbuffer; - if (mUseMC) { - labelbuffer = pc.inputs().get>("labels"); - mc2rofs = pc.inputs().get>("MC2ROframes"); + + // filter input and compose + std::array, NLayers> digits{}; + std::array, NLayers> rofs{}; + std::array, NLayers> labelsbuffer{}; + for (const DataRef& ref : InputRecordWalker{pc.inputs(), mFilter}) { + auto const* dh = DataRefUtils::getHeader(ref); + if (DataRefUtils::match(ref, {"digits", ConcreteDataTypeMatcher{"IT3", "DIGITS"}})) { + digits[dh->subSpecification] = pc.inputs().get>(ref); + } + if (DataRefUtils::match(ref, {"ROframe", ConcreteDataTypeMatcher{"IT3", "DIGITSROF"}})) { + rofs[dh->subSpecification] = pc.inputs().get>(ref); + } + if (DataRefUtils::match(ref, {"labels", ConcreteDataTypeMatcher{"IT3", "DIGITSMCTR"}})) { + labelsbuffer[dh->subSpecification] = pc.inputs().get>(ref); + } } - o2::dataformats::ConstMCTruthContainerView labels(labelbuffer); - LOGP(info, "ITS3Clusterer pulled {} digits in {} ROFs", digits.size(), rofs.size()); - LOGP(info, "ITS3Clusterer pulled {} labels", labels.getNElements()); + // query the first orbit in this TF + const auto firstTForbit = pc.services().get().firstTForbit; + const o2::InteractionRecord firstIR(0, firstTForbit); + const auto& par = o2::itsmft::DPLAlpideParam::Instance(); + uint64_t nClusters{0}; + TStopwatch sw; o2::itsmft::DigitPixelReader reader; - reader.setSquashingDepth(mClusterer->getMaxROFDepthToSquash()); - reader.setSquashingDist(mClusterer->getMaxRowColDiffToMask()); // Sharing same parameter/logic with masking - reader.setMaxBCSeparationToSquash(mClusterer->getMaxBCSeparationToSquash()); - reader.setDigits(digits); - reader.setROFRecords(rofs); - if (mUseMC) { - reader.setMC2ROFRecords(mc2rofs); - reader.setDigitsMCTruth(labels.getIndexedSize() > 0 ? &labels : nullptr); - } - reader.init(); - auto orig = o2::header::gDataOriginITS; - std::vector clusCompVec; - std::vector clusROFVec; - std::vector clusPattVec; - - std::unique_ptr> clusterLabels; - if (mUseMC) { - clusterLabels = std::make_unique>(); - } - mClusterer->process(mNThreads, reader, &clusCompVec, &clusPattVec, &clusROFVec, clusterLabels.get()); - pc.outputs().snapshot(Output{orig, "COMPCLUSTERS", 0}, clusCompVec); - pc.outputs().snapshot(Output{orig, "CLUSTERSROF", 0}, clusROFVec); - pc.outputs().snapshot(Output{orig, "PATTERNS", 0}, clusPattVec); - - if (mUseMC) { - pc.outputs().snapshot(Output{orig, "CLUSTERSMCTR", 0}, *clusterLabels); // at the moment requires snapshot - std::vector clusterMC2ROframes(mc2rofs.size()); - for (int i = mc2rofs.size(); i--;) { - clusterMC2ROframes[i] = mc2rofs[i]; // Simply, replicate it from digits ? + for (uint32_t iLayer{0}; iLayer < NLayers; ++iLayer) { + sw.Start(); + LOG(info) << "IT3Clusterer on layer " << iLayer << " pulled " << digits[iLayer].size() << " digits, in " << rofs[iLayer].size() << " RO frames"; + mClusterer->setMaxROFDepthToSquash(mClusterer->getMaxROFDepthToSquash(iLayer)); + o2::dataformats::ConstMCTruthContainerView labels(labelsbuffer[iLayer]); + reader.setSquashingDepth(mClusterer->getMaxROFDepthToSquash(iLayer)); + reader.setSquashingDist(mClusterer->getMaxRowColDiffToMask()); // Sharing same parameter/logic with masking + reader.setMaxBCSeparationToSquash(mClusterer->getMaxBCSeparationToSquash(iLayer)); + reader.setDigits(digits[iLayer]); + reader.setROFRecords(rofs[iLayer]); + if (mUseMC) { + LOG(info) << "IT3Clusterer on layer " << iLayer << " pulled " << labels.getNElements() << " labels "; + reader.setDigitsMCTruth(labels.getIndexedSize() > 0 ? &labels : nullptr); + } + reader.init(); + std::vector clusCompVec; + std::vector clusROFVec; + std::vector clusPattVec; + + std::unique_ptr> clusterLabels; + if (mUseMC) { + clusterLabels = std::make_unique>(); + } + mClusterer->process(mNThreads, reader, &clusCompVec, &clusPattVec, &clusROFVec, clusterLabels.get()); + + // ensure that the rof output is continuous + size_t nROFs = clusROFVec.size(); + const int nROFsPerOrbit = o2::constants::lhc::LHCMaxBunches / par.getROFLengthInBC(iLayer); + const int nROFsTF = nROFsPerOrbit * o2::base::GRPGeomHelper::getNHBFPerTF(); + // It can happen that in the digitization rofs without contributing hits are skipped or there are stray ROFs + // We will preserve the clusters as they are but the stray ROFs will be removed (leaving their clusters unaddressed). + std::vector expClusRofVec(nROFsTF); + for (int iROF{0}; iROF < nROFsTF; ++iROF) { + auto& rof = expClusRofVec[iROF]; + int orb = (iROF * par.getROFLengthInBC(iLayer) / o2::constants::lhc::LHCMaxBunches) + firstTForbit; + int bc = (iROF * par.getROFLengthInBC(iLayer) % o2::constants::lhc::LHCMaxBunches) + par.getROFDelayInBC(iLayer); + o2::InteractionRecord ir(bc, orb); + rof.setBCData(ir); + rof.setROFrame(iROF); + rof.setNEntries(0); + rof.setFirstEntry(-1); } - pc.outputs().snapshot(Output{orig, "CLUSTERSMC2ROF", 0}, clusterMC2ROframes); + uint32_t prevEntry{0}; + for (const auto& rof : clusROFVec) { + const auto& ir = rof.getBCData(); + if (ir < firstIR) { + LOGP(warn, "Discard ROF {} preceding TF 1st orbit {} on layer {}", ir.asString(), firstTForbit, iLayer); + continue; + } + auto irToFirst = ir - firstIR; + if (irToFirst.toLong() - par.getROFDelayInBC(iLayer) < 0) { + LOGP(warn, "Discard ROF {} preceding TF 1st orbit {} due to imposed ROF delay on layer {}", ir.asString(), firstTForbit, iLayer); + continue; + } + irToFirst -= par.getROFDelayInBC(iLayer); + const long irROF = irToFirst.toLong() / par.getROFLengthInBC(iLayer); + if (irROF >= nROFsTF) { + LOGP(warn, "Discard ROF {} exceeding TF orbit range on layer {}", ir.asString(), iLayer); + continue; + } + auto& expROF = expClusRofVec[irROF]; + if (expROF.getNEntries() == 0) { + expROF.setFirstEntry(rof.getFirstEntry()); + expROF.setNEntries(rof.getNEntries()); + } else { + if (expROF.getNEntries() < rof.getNEntries()) { + LOGP(warn, "Repeating {} with {} clusters, prefer to already processed instance with {} clusters on layer {}", rof.asString(), rof.getNEntries(), expROF.getNEntries(), iLayer); + expROF.setFirstEntry(rof.getFirstEntry()); + expROF.setNEntries(rof.getNEntries()); + } else { + LOGP(warn, "Repeating {} with {} clusters, discard preferring already processed instance with {} clusters on layer {}", rof.asString(), rof.getNEntries(), expROF.getNEntries(), iLayer); + } + } + } + int prevLast{0}; + for (auto& rof : expClusRofVec) { + if (rof.getFirstEntry() < 0) { + rof.setFirstEntry(prevLast); + } + prevLast = rof.getFirstEntry() + rof.getNEntries(); + } + nROFs = expClusRofVec.size(); + pc.outputs().snapshot(Output{"ITS", "CLUSTERSROF", iLayer}, expClusRofVec); + + pc.outputs().snapshot(Output{"ITS", "COMPCLUSTERS", iLayer}, clusCompVec); + pc.outputs().snapshot(Output{"ITS", "PATTERNS", iLayer}, clusPattVec); + + nClusters += clusCompVec.size(); + + if (mUseMC) { + pc.outputs().snapshot(Output{"ITS", "CLUSTERSMCTR", iLayer}, *clusterLabels); // at the moment requires snapshot + } + reader.reset(); + + sw.Stop(); + LOG(info) << "IT3Clusterer on layer " << iLayer << " pushed " << clusCompVec.size() << " clusters, in " << nROFs << " RO frames in " << sw.RealTime() << " s"; } - // TODO: in principle, after masking "overflow" pixels the MC2ROFRecord maxROF supposed to change, nominally to minROF - // -> consider recalculationg maxROF - LOG(info) << "ITS3Clusterer pushed " << clusCompVec.size() << " clusters, in " << clusROFVec.size() << " RO frames"; + LOG(info) << "IT3Clusterer produced " << nClusters << " clusters"; } ///_______________________________________ @@ -112,12 +193,14 @@ void ClustererDPL::updateTimeDependentParams(ProcessingContext& pc) pc.inputs().get("cldict"); // just to trigger the finaliseCCDB pc.inputs().get*>("alppar"); pc.inputs().get*>("cluspar"); + mClusterer->setContinuousReadOut(true); // settings for the fired pixel overflow masking const auto& alpParams = o2::itsmft::DPLAlpideParam::Instance(); const auto& clParams = o2::itsmft::ClustererParam::Instance(); if (clParams.maxBCDiffToMaskBias > 0 && clParams.maxBCDiffToSquashBias > 0) { - LOGP(fatal, "maxBCDiffToMaskBias = {} and maxBCDiffToMaskBias = {} cannot be set at the same time. Either set masking or squashing with a BCDiff > 0", clParams.maxBCDiffToMaskBias, clParams.maxBCDiffToSquashBias); + LOGP(fatal, "maxBCDiffToMaskBias = {} and maxBCDiffToSquashBias = {} cannot be set at the same time. Either set masking or squashing with a BCDiff > 0", clParams.maxBCDiffToMaskBias, clParams.maxBCDiffToSquashBias); } + mClusterer->setDropHugeClusters(clParams.dropHugeClusters); auto nbc = clParams.maxBCDiffToMaskBias; nbc += mClusterer->isContinuousReadOut() ? alpParams.roFrameLengthInBC : (alpParams.roFrameLengthTrig / o2::constants::lhc::LHCBunchSpacingNS); mClusterer->setMaxBCSeparationToMask(nbc); @@ -129,8 +212,14 @@ void ClustererDPL::updateTimeDependentParams(ProcessingContext& pc) if (clParams.maxSOTMUS > 0 && rofBC > 0) { nROFsToSquash = 2 + int(clParams.maxSOTMUS / (rofBC * o2::constants::lhc::LHCBunchSpacingMUS)); // use squashing } - mClusterer->setMaxROFDepthToSquash(clParams.maxBCDiffToSquashBias > 0 ? nROFsToSquash : 0); - mClusterer->print(); + mClusterer->setMaxROFDepthToSquash(nROFsToSquash); + if (mClusterer->isContinuousReadOut()) { // almost tautological + for (int iLayer{0}; iLayer < NLayers; ++iLayer) { + mClusterer->addMaxBCSeparationToSquash(alpParams.getROFLengthInBC(iLayer) + clParams.getMaxBCDiffToSquashBias(iLayer)); + mClusterer->addMaxROFDepthToSquash((clParams.getMaxBCDiffToSquashBias(iLayer) > 0) ? 2 + int(clParams.maxSOTMUS / (alpParams.getROFLengthInBC(iLayer) * o2::constants::lhc::LHCBunchSpacingMUS)) : 0); + } + } + mClusterer->print(false); } // we may have other params which need to be queried regularly } @@ -171,37 +260,36 @@ void ClustererDPL::endOfStream(o2::framework::EndOfStreamContext& /*ec*/) DataProcessorSpec getClustererSpec(bool useMC) { std::vector inputs; - inputs.emplace_back("digits", "IT3", "DIGITS", 0, Lifetime::Timeframe); - inputs.emplace_back("ROframes", "IT3", "DIGITSROF", 0, Lifetime::Timeframe); + std::vector outputs; + for (uint32_t iLayer = 0; iLayer < ClustererDPL::NLayers; ++iLayer) { + inputs.emplace_back("digits", "IT3", "DIGITS", iLayer, Lifetime::Timeframe); + inputs.emplace_back("ROframes", "IT3", "DIGITSROF", iLayer, Lifetime::Timeframe); + outputs.emplace_back("ITS", "COMPCLUSTERS", iLayer, Lifetime::Timeframe); + outputs.emplace_back("ITS", "PATTERNS", iLayer, Lifetime::Timeframe); + outputs.emplace_back("ITS", "CLUSTERSROF", iLayer, Lifetime::Timeframe); + if (useMC) { + inputs.emplace_back("labels", "IT3", "DIGITSMCTR", iLayer, Lifetime::Timeframe); + outputs.emplace_back("ITS", "CLUSTERSMCTR", iLayer, Lifetime::Timeframe); + } + } inputs.emplace_back("cldict", "IT3", "CLUSDICT", 0, Lifetime::Condition, ccdbParamSpec("IT3/Calib/ClusterDictionary")); inputs.emplace_back("cluspar", "ITS", "CLUSPARAM", 0, Lifetime::Condition, ccdbParamSpec("ITS/Config/ClustererParam")); inputs.emplace_back("alppar", "ITS", "ALPIDEPARAM", 0, Lifetime::Condition, ccdbParamSpec("ITS/Config/AlpideParam")); auto ggRequest = std::make_shared(false, // orbitResetTime - false, // GRPECS + true, // GRPECS false, // GRPLHCIF false, // GRPMagField false, // askMatLUT o2::base::GRPGeomRequest::None, // geometry inputs, true); - std::vector outputs; - outputs.emplace_back("ITS", "COMPCLUSTERS", 0, Lifetime::Timeframe); - outputs.emplace_back("ITS", "PATTERNS", 0, Lifetime::Timeframe); - outputs.emplace_back("ITS", "CLUSTERSROF", 0, Lifetime::Timeframe); - - if (useMC) { - inputs.emplace_back("labels", "IT3", "DIGITSMCTR", 0, Lifetime::Timeframe); - inputs.emplace_back("MC2ROframes", "IT3", "DIGITSMC2ROF", 0, Lifetime::Timeframe); - outputs.emplace_back("ITS", "CLUSTERSMCTR", 0, Lifetime::Timeframe); - outputs.emplace_back("ITS", "CLUSTERSMC2ROF", 0, Lifetime::Timeframe); - } return DataProcessorSpec{ - "its3-clusterer", - inputs, - outputs, - AlgorithmSpec{adaptFromTask(ggRequest, useMC)}, - Options{ + .name = "its3-clusterer", + .inputs = inputs, + .outputs = outputs, + .algorithm = AlgorithmSpec{adaptFromTask(ggRequest, useMC)}, + .options = Options{ {"ignore-cluster-dictionary", VariantType::Bool, false, {"do not use cluster dictionary, always store explicit patterns"}}, {"nthreads", VariantType::Int, 1, {"Number of clustering threads"}}}}; } diff --git a/Detectors/Upgrades/ITS3/workflow/src/DigitReaderSpec.cxx b/Detectors/Upgrades/ITS3/workflow/src/DigitReaderSpec.cxx index c2380be77f956..35068afadafad 100644 --- a/Detectors/Upgrades/ITS3/workflow/src/DigitReaderSpec.cxx +++ b/Detectors/Upgrades/ITS3/workflow/src/DigitReaderSpec.cxx @@ -13,8 +13,6 @@ #include -#include "TTree.h" - #include "Framework/ControlService.h" #include "Framework/ConfigParamRegistry.h" #include "Framework/Logger.h" @@ -23,65 +21,51 @@ #include "SimulationDataFormat/ConstMCTruthContainer.h" #include "SimulationDataFormat/IOMCTruthContainerView.h" #include +#include using namespace o2::framework; using namespace o2::itsmft; -namespace o2 -{ -namespace its3 +namespace o2::its3 { -DigitReader::DigitReader(o2::detectors::DetID id, bool useMC, bool useCalib) +ITS3DigitReader::ITS3DigitReader(bool useMC, bool useCalib) : mUseMC(useMC), mUseCalib(useCalib), mDetNameLC(mDetName = "IT3"), mDigTreeName("o2sim") { - assert(id == o2::detectors::DetID::IT3); - mDetNameLC = mDetName = id.getName(); - mDigTreeName = "o2sim"; - - mDigitBranchName = mDetName + mDigitBranchName; + mDigBranchName = mDetName + mDigBranchName; mDigROFBranchName = mDetName + mDigROFBranchName; - mCalibBranchName = mDetName + mCalibBranchName; + mDigMCTruthBranchName = mDetName + mDigMCTruthBranchName; - mDigtMCTruthBranchName = mDetName + mDigtMCTruthBranchName; - mDigtMC2ROFBranchName = mDetName + mDigtMC2ROFBranchName; - - mUseMC = useMC; - mUseCalib = useCalib; std::transform(mDetNameLC.begin(), mDetNameLC.end(), mDetNameLC.begin(), ::tolower); } -void DigitReader::init(InitContext& ic) +void ITS3DigitReader::init(InitContext& ic) { mFileName = ic.options().get((mDetNameLC + "-digit-infile").c_str()); connectTree(mFileName); } -void DigitReader::run(ProcessingContext& pc) +void ITS3DigitReader::run(ProcessingContext& pc) { auto ent = mTree->GetReadEntry() + 1; assert(ent < mTree->GetEntries()); // this should not happen - o2::dataformats::IOMCTruthContainerView* plabels = nullptr; - if (mUseMC) { - mTree->SetBranchAddress(mDigtMCTruthBranchName.c_str(), &plabels); - } mTree->GetEntry(ent); - LOG(info) << mDetName << "DigitReader pushes " << mDigROFRec.size() << " ROFRecords, " - << mDigits.size() << " digits at entry " << ent; - - // This is a very ugly way of providing DataDescription, which anyway does not need to contain detector name. - // To be fixed once the names-definition class is ready - pc.outputs().snapshot(Output{mOrigin, "DIGITSROF", 0}, mDigROFRec); - pc.outputs().snapshot(Output{mOrigin, "DIGITS", 0}, mDigits); - if (mUseCalib) { - pc.outputs().snapshot(Output{mOrigin, "GBTCALIB", 0}, mCalib); - } - - if (mUseMC) { - auto& sharedlabels = pc.outputs().make>(Output{mOrigin, "DIGITSMCTR", 0}); - plabels->copyandflatten(sharedlabels); - delete plabels; - pc.outputs().snapshot(Output{mOrigin, "DIGITSMC2ROF", 0}, mDigMC2ROFs); + for (uint32_t iLayer = 0; iLayer < NLayers; ++iLayer) { + if (!mDigROFRec[iLayer] || !mDigits[iLayer]) { + throw std::runtime_error("ITS3 digit reader requires all 7 layer branches to be present and populated in every entry"); + } + LOG(info) << mDetName << "DigitReader pushes " << mDigROFRec[iLayer]->size() << " ROFRecords, " << mDigits[iLayer]->size() << " digits at entry " << ent << " on layer " << iLayer; + pc.outputs().snapshot(Output{mOrigin, "DIGITSROF", iLayer}, *mDigROFRec[iLayer]); + pc.outputs().snapshot(Output{mOrigin, "DIGITS", iLayer}, *mDigits[iLayer]); + if (mUseMC) { + if (!mPLabels[iLayer]) { + throw std::runtime_error("ITS3 digit reader requires MC truth branches for all 7 layers to be present and populated in every entry"); + } + auto& sharedlabels = pc.outputs().make>(Output{mOrigin, "DIGITSMCTR", iLayer}); + mPLabels[iLayer]->copyandflatten(sharedlabels); + delete mPLabels[iLayer]; + mPLabels[iLayer] = nullptr; + } } if (mTree->GetReadEntry() + 1 >= mTree->GetEntries()) { @@ -90,27 +74,38 @@ void DigitReader::run(ProcessingContext& pc) } } -void DigitReader::connectTree(const std::string& filename) +template +void ITS3DigitReader::setBranchAddress(const std::string& base, Ptr& addr, int layer) +{ + const auto name = getBranchName(base, layer); + if (Int_t ret = mTree->SetBranchAddress(name.c_str(), &addr); ret != 0) { + LOGP(fatal, "failed to set branch address for {} ret={}", name, ret); + } +} + +void ITS3DigitReader::connectTree(const std::string& filename) { mTree.reset(nullptr); // in case it was already loaded mFile.reset(TFile::Open(filename.c_str())); assert(mFile && !mFile->IsZombie()); mTree.reset((TTree*)mFile->Get(mDigTreeName.c_str())); assert(mTree); - - mTree->SetBranchAddress(mDigROFBranchName.c_str(), &mDigROFRecPtr); - mTree->SetBranchAddress(mDigitBranchName.c_str(), &mDigitsPtr); - if (mUseCalib) { - if (!mTree->GetBranch(mCalibBranchName.c_str())) { - throw std::runtime_error("GBT calibration data requested but not found in the tree"); + for (int iLayer = 0; iLayer < NLayers; ++iLayer) { + const auto rofBranchName = getBranchName(mDigROFBranchName, iLayer); + const auto digBranchName = getBranchName(mDigBranchName, iLayer); + if (!mTree->GetBranch(rofBranchName.c_str()) || !mTree->GetBranch(digBranchName.c_str())) { + throw std::runtime_error("ITS3 digit reader requires all 7 layer branches in the input file"); } - mTree->SetBranchAddress(mCalibBranchName.c_str(), &mCalibPtr); - } - if (mUseMC) { - if (!mTree->GetBranch(mDigtMC2ROFBranchName.c_str()) || !mTree->GetBranch(mDigtMCTruthBranchName.c_str())) { - throw std::runtime_error("MC data requested but not found in the tree"); + setBranchAddress(mDigROFBranchName, mDigROFRec[iLayer], iLayer); + setBranchAddress(mDigBranchName, mDigits[iLayer], iLayer); + if (mUseMC) { + if (!mTree->GetBranch(getBranchName(mDigMCTruthBranchName, iLayer).c_str())) { + throw std::runtime_error("ITS3 digit reader requires MC truth branches for all 7 layers in the input file"); + } + if (!mPLabels[iLayer]) { + setBranchAddress(mDigMCTruthBranchName, mPLabels[iLayer], iLayer); + } } - mTree->SetBranchAddress(mDigtMC2ROFBranchName.c_str(), &mDigMC2ROFsPtr); } LOG(info) << "Loaded tree from " << filename << " with " << mTree->GetEntries() << " entries"; } @@ -118,14 +113,12 @@ void DigitReader::connectTree(const std::string& filename) DataProcessorSpec getITS3DigitReaderSpec(bool useMC, bool useCalib, std::string defname) { std::vector outputSpec; - outputSpec.emplace_back("IT3", "DIGITS", 0, Lifetime::Timeframe); - outputSpec.emplace_back("IT3", "DIGITSROF", 0, Lifetime::Timeframe); - if (useCalib) { - outputSpec.emplace_back("IT3", "GBTCALIB", 0, Lifetime::Timeframe); - } - if (useMC) { - outputSpec.emplace_back("IT3", "DIGITSMCTR", 0, Lifetime::Timeframe); - outputSpec.emplace_back("IT3", "DIGITSMC2ROF", 0, Lifetime::Timeframe); + for (uint32_t iLayer = 0; iLayer < ITS3DigitReader::NLayers; ++iLayer) { + outputSpec.emplace_back("IT3", "DIGITS", iLayer, Lifetime::Timeframe); + outputSpec.emplace_back("IT3", "DIGITSROF", iLayer, Lifetime::Timeframe); + if (useMC) { + outputSpec.emplace_back("IT3", "DIGITSMCTR", iLayer, Lifetime::Timeframe); + } } return DataProcessorSpec{ @@ -137,5 +130,4 @@ DataProcessorSpec getITS3DigitReaderSpec(bool useMC, bool useCalib, std::string {"it3-digit-infile", VariantType::String, defname, {"Name of the input digit file"}}}}; } -} // namespace its3 -} // namespace o2 +} // namespace o2::its3 diff --git a/Detectors/Upgrades/ITS3/workflow/src/DigitWriterSpec.cxx b/Detectors/Upgrades/ITS3/workflow/src/DigitWriterSpec.cxx index 8fa67b28e56c1..b9f97fa76964b 100644 --- a/Detectors/Upgrades/ITS3/workflow/src/DigitWriterSpec.cxx +++ b/Detectors/Upgrades/ITS3/workflow/src/DigitWriterSpec.cxx @@ -14,14 +14,13 @@ #include "ITS3Workflow/DigitWriterSpec.h" #include "DPLUtils/MakeRootTreeWriterSpec.h" #include "DataFormatsITSMFT/Digit.h" -#include "DataFormatsITSMFT/GBTCalibData.h" #include "Headers/DataHeader.h" #include "DetectorsCommonDataFormats/DetID.h" #include "DataFormatsITSMFT/ROFRecord.h" #include "SimulationDataFormat/ConstMCTruthContainer.h" #include "SimulationDataFormat/IOMCTruthContainerView.h" #include "SimulationDataFormat/MCCompLabel.h" -#include +#include #include #include @@ -45,10 +44,17 @@ DataProcessorSpec getDigitWriterSpec(bool mctruth, bool dec, bool calib, o2::hea std::string detStrL = dec ? "o2_" : ""; // for decoded digits prepend by o2 detStrL += detStr; std::transform(detStrL.begin(), detStrL.end(), detStrL.begin(), ::tolower); - auto logger = [](std::vector const& inDigits) { - LOG(info) << "RECEIVED DIGITS SIZE " << inDigits.size(); - }; + auto digitSizes = std::make_shared>(7, 0); + auto digitSizeGetter = [digitSizes](std::vector const& inDigits, DataRef const& ref) { + auto const* dh = DataRefUtils::getHeader(ref); + (*digitSizes)[dh->subSpecification] = inDigits.size(); + }; + auto rofSizes = std::make_shared>(7, 0); + auto rofSizeGetter = [rofSizes](std::vector const& inROFs, DataRef const& ref) { + auto const* dh = DataRefUtils::getHeader(ref); + (*rofSizes)[dh->subSpecification] = inROFs.size(); + }; // the callback to be set as hook for custom action when the writer is closed auto finishWriting = [](TFile* outputfile, TTree* outputtree) { const auto* brArr = outputtree->GetListOfBranches(); @@ -68,9 +74,11 @@ DataProcessorSpec getDigitWriterSpec(bool mctruth, bool dec, bool calib, o2::hea // handler for labels // This is necessary since we can't store the original label buffer in a ROOT entry -- as is -- if it exceeds a certain size. // We therefore convert it to a special split class. - auto fillLabels = [](TBranch& branch, std::vector const& labelbuffer, DataRef const& /*ref*/) { + auto fillLabels = [detStr, digitSizes, rofSizes](TBranch& branch, std::vector const& labelbuffer, DataRef const& ref) { o2::dataformats::ConstMCTruthContainerView labels(labelbuffer); - LOG(info) << "WRITING " << labels.getNElements() << " LABELS "; + auto const* dh = DataRefUtils::getHeader(ref); + auto layer = static_cast(dh->subSpecification); + LOG(info) << detStr << ": WRITING " << labels.getNElements() << " LABELS FOR LAYER " << layer << " WITH " << (*digitSizes)[layer] << " DIGITS IN " << (*rofSizes)[layer] << " ROFS"; o2::dataformats::IOMCTruthContainerView outputcontainer; auto ptr = &outputcontainer; @@ -80,25 +88,47 @@ DataProcessorSpec getDigitWriterSpec(bool mctruth, bool dec, bool calib, o2::hea br->ResetAddress(); }; + auto getIndex = [](DataRef const& ref) -> size_t { + auto const* dh = DataRefUtils::getHeader(ref); + return static_cast(dh->subSpecification); + }; + auto getName = [](std::string base, size_t index) -> std::string { + return base += "_" + std::to_string(index); + }; + + std::vector vecInpSpecDig, vecInpSpecROF, vecInpSpecLbl; + vecInpSpecDig.reserve(7); + vecInpSpecROF.reserve(7); + vecInpSpecLbl.reserve(7); + for (int iLayer = 0; iLayer < 7; iLayer++) { + vecInpSpecDig.emplace_back(getName("digits", iLayer), "IT3", "DIGITS", iLayer); + vecInpSpecROF.emplace_back(getName("digitsROF", iLayer), "IT3", "DIGITSROF", iLayer); + vecInpSpecLbl.emplace_back(getName("digitsMCTR", iLayer), "IT3", "DIGITSMCTR", iLayer); + } + return MakeRootTreeWriterSpec((detStr + "DigitWriter" + (dec ? "_dec" : "")).c_str(), (detStrL + "digits.root").c_str(), MakeRootTreeWriterSpec::TreeAttributes{"o2sim", "Digits tree"}, MakeRootTreeWriterSpec::CustomClose(finishWriting), // in case of labels we first read them as std::vector and process them correctly in the fillLabels hook - BranchDefinition>{InputSpec{"digitsMCTR", detOrig, "DIGITSMCTR", 0}, - (detStr + "DigitMCTruth").c_str(), - (mctruth ? 1 : 0), fillLabels}, - BranchDefinition>{InputSpec{"digitsMC2ROF", detOrig, "DIGITSMC2ROF", 0}, - (detStr + "DigitMC2ROF").c_str(), - (mctruth ? 1 : 0)}, - BranchDefinition>{InputSpec{"digits", detOrig, "DIGITS", 0}, - (detStr + "Digit").c_str(), - logger}, - BranchDefinition>{InputSpec{"calib", detOrig, "GBTCALIB", 0}, - (detStr + "Calib").c_str(), - (calib ? 1 : 0)}, - BranchDefinition>{InputSpec{"digitsROF", detOrig, "DIGITSROF", 0}, - (detStr + "DigitROF").c_str()})(); + BranchDefinition>{vecInpSpecDig, + detStr + "Digit", "digit-branch", + 7, + digitSizeGetter, + getIndex, + getName}, + BranchDefinition>{vecInpSpecROF, + detStr + "DigitROF", "digit-rof-branch", + 7, + rofSizeGetter, + getIndex, + getName}, + BranchDefinition>{vecInpSpecLbl, + "IT3DigitMCTruth", "digit-mctruth-branch", + (mctruth ? 7 : 0), + fillLabels, + getIndex, + getName})(); } DataProcessorSpec getITS3DigitWriterSpec(bool mctruth, bool dec, bool calib) diff --git a/Detectors/Upgrades/ITS3/workflow/src/RecoWorkflow.cxx b/Detectors/Upgrades/ITS3/workflow/src/RecoWorkflow.cxx index f27fda19fe00c..4f20c8ceea43d 100644 --- a/Detectors/Upgrades/ITS3/workflow/src/RecoWorkflow.cxx +++ b/Detectors/Upgrades/ITS3/workflow/src/RecoWorkflow.cxx @@ -40,7 +40,7 @@ framework::WorkflowSpec getWorkflow(bool useMC, its::TrackingMode::Type trmode, } if (!disableRootOutput) { - specs.emplace_back(o2::itsmft::getITSClusterWriterSpec(useMC, false)); + specs.emplace_back(o2::itsmft::getITSClusterWriterSpec(useMC, true)); } if (trmode != its::TrackingMode::Off) { @@ -66,11 +66,11 @@ framework::WorkflowSpec getWorkflow(bool useMC, its::TrackingMode::Type trmode, std::move(ggInputs.begin(), ggInputs.end(), std::back_inserter(taskInputs)); specs.emplace_back(DataProcessorSpec{ - "its3-gpu-tracker", - taskInputs, - task->outputs(), - AlgorithmSpec{adoptTask(task)}, - taskOptions}); + .name = "its3-gpu-tracker", + .inputs = taskInputs, + .outputs = task->outputs(), + .algorithm = AlgorithmSpec{adoptTask(task)}, + .options = taskOptions}); } else { specs.emplace_back(o2::its3::getTrackerSpec(useMC, useGeom, useTrig, trmode, overrideBeamPosition, dtype)); } diff --git a/Detectors/Upgrades/ITS3/workflow/src/TrackerSpec.cxx b/Detectors/Upgrades/ITS3/workflow/src/TrackerSpec.cxx index 8db02d7227e7f..c6a5405cadb94 100644 --- a/Detectors/Upgrades/ITS3/workflow/src/TrackerSpec.cxx +++ b/Detectors/Upgrades/ITS3/workflow/src/TrackerSpec.cxx @@ -46,7 +46,7 @@ TrackerDPL::TrackerDPL(std::shared_ptr gr, const bool overrBeamEst, o2::gpu::gpudatatypes::DeviceType dType) : mGGCCDBRequest(gr), mRecChain{o2::gpu::GPUReconstruction::CreateInstance(dType, true)}, - mITS3TrackingInterface{isMC, false, trgType, overrBeamEst} + mITS3TrackingInterface{isMC, true, trgType, overrBeamEst} { mITS3TrackingInterface.setTrackingMode(trMode); } @@ -91,16 +91,22 @@ void TrackerDPL::endOfStream(EndOfStreamContext& ec) DataProcessorSpec getTrackerSpec(bool useMC, bool useGeom, int trgType, its::TrackingMode::Type trMode, const bool overrBeamEst, o2::gpu::gpudatatypes::DeviceType dType) { std::vector inputs; - inputs.emplace_back("compClusters", "ITS", "COMPCLUSTERS", 0, Lifetime::Timeframe); - inputs.emplace_back("patterns", "ITS", "PATTERNS", 0, Lifetime::Timeframe); - inputs.emplace_back("ROframes", "ITS", "CLUSTERSROF", 0, Lifetime::Timeframe); + std::vector outputs; + for (uint32_t iLayer{0}; iLayer < 7; ++iLayer) { + inputs.emplace_back("compClusters", "ITS", "COMPCLUSTERS", iLayer, Lifetime::Timeframe); + inputs.emplace_back("patterns", "ITS", "PATTERNS", iLayer, Lifetime::Timeframe); + inputs.emplace_back("ROframes", "ITS", "CLUSTERSROF", iLayer, Lifetime::Timeframe); + if (useMC) { + inputs.emplace_back("itsmclabels", "ITS", "CLUSTERSMCTR", iLayer, Lifetime::Timeframe); + } + } if (trgType == 1) { inputs.emplace_back("phystrig", "ITS", "PHYSTRIG", 0, Lifetime::Timeframe); } else if (trgType == 2) { inputs.emplace_back("phystrig", "TRD", "TRKTRGRD", 0, Lifetime::Timeframe); } - inputs.emplace_back("cldict", "IT3", "CLUSDICT", 0, Lifetime::Condition, ccdbParamSpec("IT3/Calib/ClusterDictionary")); - inputs.emplace_back("alppar", "ITS", "ALPIDEPARAM", 0, Lifetime::Condition, ccdbParamSpec("ITS/Config/AlpideParam")); + inputs.emplace_back("itscldict", "IT3", "CLUSDICT", 0, Lifetime::Condition, ccdbParamSpec("IT3/Calib/ClusterDictionary")); + inputs.emplace_back("itsalppar", "ITS", "ALPIDEPARAM", 0, Lifetime::Condition, ccdbParamSpec("ITS/Config/AlpideParam")); auto ggRequest = std::make_shared(false, // orbitResetTime true, // GRPECS false, // GRPLHCIF @@ -117,7 +123,6 @@ DataProcessorSpec getTrackerSpec(bool useMC, bool useGeom, int trgType, its::Tra inputs.emplace_back("meanvtx", "GLO", "MEANVERTEX", 0, Lifetime::Condition, ccdbParamSpec("GLO/Calib/MeanVertex", {}, 1)); } - std::vector outputs; outputs.emplace_back("ITS", "TRACKS", 0, Lifetime::Timeframe); outputs.emplace_back("ITS", "TRACKCLSID", 0, Lifetime::Timeframe); outputs.emplace_back("ITS", "ITSTrackROF", 0, Lifetime::Timeframe); @@ -126,12 +131,9 @@ DataProcessorSpec getTrackerSpec(bool useMC, bool useGeom, int trgType, its::Tra outputs.emplace_back("ITS", "IRFRAMES", 0, Lifetime::Timeframe); if (useMC) { - inputs.emplace_back("itsmclabels", "ITS", "CLUSTERSMCTR", 0, Lifetime::Timeframe); - inputs.emplace_back("ITSMC2ROframes", "ITS", "CLUSTERSMC2ROF", 0, Lifetime::Timeframe); outputs.emplace_back("ITS", "VERTICESMCTR", 0, Lifetime::Timeframe); outputs.emplace_back("ITS", "VERTICESMCPUR", 0, Lifetime::Timeframe); outputs.emplace_back("ITS", "TRACKSMCTR", 0, Lifetime::Timeframe); - outputs.emplace_back("ITS", "ITSTrackMC2ROF", 0, Lifetime::Timeframe); } return DataProcessorSpec{ diff --git a/Steer/DigitizerWorkflow/src/ITS3DigitizerSpec.cxx b/Steer/DigitizerWorkflow/src/ITS3DigitizerSpec.cxx index 19552a407ec57..730e219da97a8 100644 --- a/Steer/DigitizerWorkflow/src/ITS3DigitizerSpec.cxx +++ b/Steer/DigitizerWorkflow/src/ITS3DigitizerSpec.cxx @@ -21,6 +21,7 @@ #include "DataFormatsITSMFT/Digit.h" #include "SimulationDataFormat/ConstMCTruthContainer.h" #include "DetectorsBase/BaseDPLDigitizer.h" +#include "DetectorsRaw/HBFUtils.h" #include "DetectorsCommonDataFormats/DetID.h" #include "DetectorsCommonDataFormats/SimTraits.h" #include "DataFormatsParameters/GRPObject.h" @@ -40,22 +41,6 @@ using namespace o2::framework; using SubSpecificationType = o2::framework::DataAllocator::SubSpecificationType; -namespace -{ -std::vector makeOutChannels(o2::header::DataOrigin detOrig, bool mctruth) -{ - std::vector outputs; - outputs.emplace_back(detOrig, "DIGITS", 0, Lifetime::Timeframe); - outputs.emplace_back(detOrig, "DIGITSROF", 0, Lifetime::Timeframe); - if (mctruth) { - outputs.emplace_back(detOrig, "DIGITSMC2ROF", 0, Lifetime::Timeframe); - outputs.emplace_back(detOrig, "DIGITSMCTR", 0, Lifetime::Timeframe); - } - outputs.emplace_back(detOrig, "ROMode", 0, Lifetime::Timeframe); - return outputs; -} -} // namespace - namespace o2::its3 { using namespace o2::base; @@ -76,13 +61,20 @@ class ITS3DPLDigitizerTask : BaseDPLDigitizer if (mFinished) { return; } + mFirstOrbitTF = pc.services().get().firstTForbit; + const o2::InteractionRecord firstIR(0, mFirstOrbitTF); updateTimeDependentParams(pc); + TStopwatch timer; + timer.Start(); + LOG(info) << " CALLING ITS3 DIGITIZATION "; + // read collision context from input auto context = pc.inputs().get("collisioncontext"); context->initSimChains(mID, mSimChains); const bool withQED = context->isQEDProvided() && !mDisableQED; auto& timesview = context->getEventRecords(withQED); + const auto& par = o2::itsmft::DPLAlpideParam::Instance(); LOG(info) << "GOT " << timesview.size() << " COLLISSION TIMES"; LOG(info) << "SIMCHAINS " << mSimChains.size(); @@ -90,107 +82,150 @@ class ITS3DPLDigitizerTask : BaseDPLDigitizer if (timesview.empty()) { return; } - TStopwatch timer; - timer.Start(); - LOG(info) << " CALLING ITS3 DIGITIZATION "; - mDigitizer.setDigits(&mDigits); - mDigitizer.setROFRecords(&mROFRecords); - mDigitizer.setMCLabels(&mLabels); + uint64_t nDigits{0}; + for (uint32_t iLayer = 0; iLayer < 7; ++iLayer) { + mDigitizer.setDigits(&mDigits[iLayer]); + mDigitizer.setROFRecords(&mROFRecords[iLayer]); + mDigitizer.setMCLabels(&mLabels[iLayer]); + mDigitizer.resetROFrameBounds(); + // digits are directly put into DPL owned resource + auto& digitsAccum = pc.outputs().make>(Output{mOrigin, "DIGITS", iLayer}); + + // rofs are accumulated first and the copied + const int nROFsPerOrbit = o2::constants::lhc::LHCMaxBunches / par.getROFLengthInBC(iLayer); + const int nROFsTF = nROFsPerOrbit * raw::HBFUtils::Instance().getNOrbitsPerTF(); + mROFRecordsAccum[iLayer].reserve(nROFsTF); + + auto accumulate = [this, &digitsAccum, &iLayer]() { + // accumulate result of single event processing on a specific layer, called after processing every event supplied + // AND after the final flushing via digitizer::fillOutputContainer + if (!mDigits[iLayer].size()) { + return; // no digits were flushed, nothing to accumulate + } + auto ndigAcc = digitsAccum.size(); + std::copy(mDigits[iLayer].begin(), mDigits[iLayer].end(), std::back_inserter(digitsAccum)); + + // fix ROFrecords references on ROF entries + auto nROFRecsOld = mROFRecordsAccum[iLayer].size(); + + for (int i = 0; i < mROFRecords[iLayer].size(); i++) { + auto& rof = mROFRecords[iLayer][i]; + rof.setFirstEntry(ndigAcc + rof.getFirstEntry()); + rof.print(); + } + + std::copy(mROFRecords[iLayer].begin(), mROFRecords[iLayer].end(), std::back_inserter(mROFRecordsAccum[iLayer])); + if (mWithMCTruth) { + mLabelsAccum[iLayer].mergeAtBack(mLabels[iLayer]); + } + LOG(info) << "Added " << mDigits[iLayer].size() << " digits on layer " << iLayer; + // clean containers from already accumulated stuff + mLabels[iLayer].clear(); + mDigits[iLayer].clear(); + mROFRecords[iLayer].clear(); + }; // and accumulate lambda + + const auto& eventParts = context->getEventParts(withQED); + const int64_t bcShift = mDigitizer.getParams().getROFrameBiasInBC(iLayer); // this accounts the misalignment and the opt. imposed rof delay + // loop over all composite collisions given from context (aka loop over all the interaction records) + for (int collID = 0; collID < timesview.size(); ++collID) { + auto irt = timesview[collID]; + if (irt.toLong() < bcShift) { // due to the ROF misalignment (+opt. delay) the collision would go to negative ROF ID, discard + continue; + } + irt -= bcShift; // account for the ROF start shift + + mDigitizer.setEventTime(irt, iLayer); + mDigitizer.resetEventROFrames(); // to estimate min/max ROF for this collID + // for each collision, loop over the constituents event and source IDs + // (background signal merging is basically taking place here) + for (const auto& part : eventParts[collID]) { - // digits are directly put into DPL owned resource - auto& digitsAccum = pc.outputs().make>(Output{mOrigin, "DIGITS", 0}); + // get the hits for this event and this source + mHits.clear(); + context->retrieveHits(mSimChains, o2::detectors::SimTraits::DETECTORBRANCHNAMES[o2::detectors::DetID::IT3][0].c_str(), part.sourceID, part.entryID, &mHits); - auto accumulate = [this, &digitsAccum]() { - // accumulate result of single event processing, called after processing every event supplied - // AND after the final flushing via digitizer::fillOutputContainer - if (mDigits.empty()) { - return; // no digits were flushed, nothing to accumulate + if (mHits.size() > 0) { + LOG(debug) << "For collision " << collID << " eventID " << part.entryID << " found " << mHits.size() << " hits "; + mDigitizer.process(&mHits, part.entryID, part.sourceID, iLayer); // call actual digitization procedure + } + } + accumulate(); } - auto ndigAcc = digitsAccum.size(); - std::copy(mDigits.begin(), mDigits.end(), std::back_inserter(digitsAccum)); - - // fix ROFrecords references on ROF entries - auto nROFRecsOld = mROFRecordsAccum.size(); - - for (int i = 0; i < mROFRecords.size(); i++) { - auto& rof = mROFRecords[i]; - rof.setFirstEntry(ndigAcc + rof.getFirstEntry()); - rof.print(); - - if (mFixMC2ROF < mMC2ROFRecordsAccum.size()) { // fix ROFRecord entry in MC2ROF records - for (int m2rid = mFixMC2ROF; m2rid < mMC2ROFRecordsAccum.size(); m2rid++) { - // need to register the ROFRecors entry for MC event starting from this entry - auto& mc2rof = mMC2ROFRecordsAccum[m2rid]; - if (rof.getROFrame() == mc2rof.minROF) { - mFixMC2ROF++; - mc2rof.rofRecordID = nROFRecsOld + i; - mc2rof.print(); - } + mDigitizer.fillOutputContainer(0xffffffff, iLayer); + accumulate(); + nDigits += digitsAccum.size(); + + // here we have all digits and labels and we can send them to consumer (aka snapshot it onto output) + // ensure that the rof output is continuous + if (nROFsTF != mROFRecordsAccum[iLayer].size()) { + // it can happen that in the digitization rofs without contributing hits are skipped + // however downstream consumers of the clusters cannot know apriori the time structure + // the cluster rofs do not account for the bias so it will start always at BC=0 + // also have to account for spillage into next TF + const size_t nROFsLayer = std::max((size_t)nROFsTF, mROFRecordsAccum[iLayer].size()); + std::vector expDigitRofVec(nROFsLayer); + for (int iROF{0}; iROF < nROFsLayer; ++iROF) { + auto& rof = expDigitRofVec[iROF]; + int orb = iROF * par.getROFLengthInBC(iLayer) / o2::constants::lhc::LHCMaxBunches + mFirstOrbitTF; + int bc = iROF * par.getROFLengthInBC(iLayer) % o2::constants::lhc::LHCMaxBunches + par.getROFDelayInBC(iLayer); + o2::InteractionRecord ir(bc, orb); + rof.setBCData(ir); + rof.setROFrame(iROF); + rof.setNEntries(0); + rof.setFirstEntry(-1); + } + uint32_t prevEntry{0}; + for (const auto& rof : mROFRecordsAccum[iLayer]) { + const auto& ir = rof.getBCData(); + if (ir < firstIR) { + LOGP(warn, "Discard ROF {} preceding TF 1st orbit {} on layer {}", ir.asString(), mFirstOrbitTF, iLayer); + continue; + } + auto irToFirst = ir - firstIR; + if (irToFirst.toLong() - par.getROFDelayInBC(iLayer) < 0) { + LOGP(warn, "Discard ROF {} preceding TF 1st orbit {} due to imposed ROF delay on layer {}", ir.asString(), mFirstOrbitTF, iLayer); + continue; + } + irToFirst -= par.getROFDelayInBC(iLayer); + const int irROF = irToFirst.toLong() / par.getROFLengthInBC(iLayer); + auto& expROF = expDigitRofVec[irROF]; + expROF.setFirstEntry(rof.getFirstEntry()); + expROF.setNEntries(rof.getNEntries()); + if (expROF.getBCData() != rof.getBCData()) { + LOGP(fatal, "detected mismatch between expected {} and received {}", expROF.asString(), rof.asString()); } } + int prevFirst{0}; + for (auto& rof : expDigitRofVec) { + if (rof.getFirstEntry() < 0) { + rof.setFirstEntry(prevFirst); + } + prevFirst = rof.getFirstEntry(); + } + // if more rofs where accumulated than ROFs possible in the TF, cut them away + // by construction expDigitRofVec is at least nROFsTF long + expDigitRofVec.resize(nROFsTF); + pc.outputs().snapshot(Output{mOrigin, "DIGITSROF", iLayer}, expDigitRofVec); + } else { + pc.outputs().snapshot(Output{mOrigin, "DIGITSROF", iLayer}, mROFRecordsAccum[iLayer]); } - - std::copy(mROFRecords.begin(), mROFRecords.end(), std::back_inserter(mROFRecordsAccum)); if (mWithMCTruth) { - mLabelsAccum.mergeAtBack(mLabels); + auto& sharedlabels = pc.outputs().make>(Output{mOrigin, "DIGITSMCTR", iLayer}); + mLabelsAccum[iLayer].flatten_to(sharedlabels); + // free space of existing label containers + mLabels[iLayer].clear_andfreememory(); + mLabelsAccum[iLayer].clear_andfreememory(); } - LOG(info) << "Added " << mDigits.size() << " digits "; - // clean containers from already accumulated stuff - mLabels.clear(); - mDigits.clear(); - mROFRecords.clear(); - }; // and accumulate lambda - - auto& eventParts = context->getEventParts(withQED); - // loop over all composite collisions given from context (aka loop over all the interaction records) - const int bcShift = mDigitizer.getParams().getROFrameBiasInBC(); - // loop over all composite collisions given from context (aka loop over all the interaction records) - for (size_t collID = 0; collID < timesview.size(); ++collID) { - auto irt = timesview[collID]; - if (irt.toLong() < bcShift) { // due to the ROF misalignment the collision would go to negative ROF ID, discard - continue; - } - irt -= bcShift; // account for the ROF start shift - - mDigitizer.setEventTime(irt); - mDigitizer.resetEventROFrames(); // to estimate min/max ROF for this collID - // for each collision, loop over the constituents event and source IDs - // (background signal merging is basically taking place here) - for (auto& part : eventParts[collID]) { - - // get the hits for this event and this source - mHits.clear(); - context->retrieveHits(mSimChains, o2::detectors::SimTraits::DETECTORBRANCHNAMES[mID][0].c_str(), part.sourceID, part.entryID, &mHits); - - if (!mHits.empty()) { - LOG(debug) << "For collision " << collID << " eventID " << part.entryID - << " found " << mHits.size() << " hits "; - mDigitizer.process(&mHits, part.entryID, part.sourceID); // call actual digitization procedure - } - } - mMC2ROFRecordsAccum.emplace_back(collID, -1, mDigitizer.getEventROFrameMin(), mDigitizer.getEventROFrameMax()); - accumulate(); } - mDigitizer.fillOutputContainer(); - accumulate(); - - // here we have all digits and labels and we can send them to consumer (aka snapshot it onto output) - - pc.outputs().snapshot(Output{mOrigin, "DIGITSROF", 0}, mROFRecordsAccum); - if (mWithMCTruth) { - pc.outputs().snapshot(Output{mOrigin, "DIGITSMC2ROF", 0}, mMC2ROFRecordsAccum); - auto& sharedlabels = pc.outputs().make>(Output{mOrigin, "DIGITSMCTR", 0}); - mLabelsAccum.flatten_to(sharedlabels); - // free space of existing label containers - mLabels.clear_andfreememory(); - mLabelsAccum.clear_andfreememory(); - } - LOG(info) << mID.getName() << ": Sending ROMode= " << mROMode << " to GRPUpdater"; + + LOG(info) << "IT3: Sending ROMode= " << mROMode << " to GRPUpdater"; pc.outputs().snapshot(Output{mOrigin, "ROMode", 0}, mROMode); timer.Stop(); LOG(info) << "Digitization took " << timer.CpuTime() << "s"; + LOG(info) << "Produced " << nDigits << " digits"; // we should be only called once; tell DPL that this process is ready to exit pc.services().get().readyToQuit(QuitRequest::Me); @@ -239,6 +274,17 @@ class ITS3DPLDigitizerTask : BaseDPLDigitizer digipar.setIBNSimSteps(doptIB.nIBSimSteps); digipar.setIBNoisePerPixel(doptIB.IBNoisePerPixel); + // staggering parameters + for (int iLayer{0}; iLayer < 7; ++iLayer) { + auto frameNS = aopt.getROFLengthInBC(iLayer) * o2::constants::lhc::LHCBunchSpacingNS; + digipar.addROFrameLayerLengthInBC(aopt.getROFLengthInBC(iLayer)); + // NOTE: the rof delay looks from the digitizer like an additional bias + digipar.addROFrameLayerBiasInBC(aopt.getROFBiasInBC(iLayer) + aopt.getROFDelayInBC(iLayer)); + digipar.addStrobeDelay(aopt.strobeDelay); + digipar.addStrobeLength(aopt.strobeLengthCont > 0 ? aopt.strobeLengthCont : frameNS - aopt.strobeDelay); + digipar.setROFrameLength(aopt.getROFLengthInBC(iLayer) * o2::constants::lhc::LHCBunchSpacingNS, iLayer); + } + mROMode = digipar.isContinuous() ? o2::parameters::GRPObject::CONTINUOUS : o2::parameters::GRPObject::PRESENT; LOG(info) << mID.getName() << " simulated in " << ((mROMode == o2::parameters::GRPObject::CONTINUOUS) ? "CONTINUOUS" : "TRIGGERED") @@ -286,26 +332,40 @@ class ITS3DPLDigitizerTask : BaseDPLDigitizer bool mWithMCTruth{true}; bool mFinished{false}; bool mDisableQED{false}; + unsigned long mFirstOrbitTF = 0x0; const o2::detectors::DetID mID{o2::detectors::DetID::IT3}; const o2::header::DataOrigin mOrigin{o2::header::gDataOriginIT3}; o2::its3::Digitizer mDigitizer{}; - std::vector mDigits{}; - std::vector mROFRecords{}; - std::vector mROFRecordsAccum{}; - std::vector mHits{}; + std::array, 7> mDigits{}; + std::array, 7> mROFRecords{}; + std::array, 7> mROFRecordsAccum{}; + std::vector mHits; std::vector* mHitsP{&mHits}; - o2::dataformats::MCTruthContainer mLabels{}; - o2::dataformats::MCTruthContainer mLabelsAccum{}; - std::vector mMC2ROFRecordsAccum{}; + std::array, 7> mLabels; + std::array, 7> mLabelsAccum; std::vector mSimChains{}; - - int mFixMC2ROF = 0; // 1st entry in mc2rofRecordsAccum to be fixed for ROFRecordID o2::parameters::GRPObject::ROMode mROMode = o2::parameters::GRPObject::PRESENT; // readout mode }; +namespace +{ +std::vector makeOutChannels(o2::header::DataOrigin detOrig, bool mctruth) +{ + std::vector outputs; + for (int iLayer{0}; iLayer < 7; ++iLayer) { + outputs.emplace_back(detOrig, "DIGITS", iLayer, Lifetime::Timeframe); + outputs.emplace_back(detOrig, "DIGITSROF", iLayer, Lifetime::Timeframe); + if (mctruth) { + outputs.emplace_back(detOrig, "DIGITSMCTR", iLayer, Lifetime::Timeframe); + } + } + outputs.emplace_back(detOrig, "ROMode", 0, Lifetime::Timeframe); + return outputs; +} +} // namespace + DataProcessorSpec getITS3DigitizerSpec(int channel, bool mctruth) { - std::string detStr = o2::detectors::DetID::getName(o2::detectors::DetID::IT3); auto detOrig = o2::header::gDataOriginIT3; std::vector inputs; inputs.emplace_back("collisioncontext", "SIM", "COLLISIONCONTEXT", static_cast(channel), Lifetime::Timeframe); @@ -316,10 +376,11 @@ DataProcessorSpec getITS3DigitizerSpec(int channel, bool mctruth) inputs.emplace_back("IT3_alpiderespvbb0", "IT3", "ALPIDERESPVbb0", 0, Lifetime::Condition, ccdbParamSpec("ITSMFT/Calib/ALPIDEResponseVbb0")); inputs.emplace_back("IT3_aptsresp", "IT3", "APTSRESP", 0, Lifetime::Condition, ccdbParamSpec("IT3/Calib/APTSResponse")); - return DataProcessorSpec{detStr + "Digitizer", - inputs, makeOutChannels(detOrig, mctruth), - AlgorithmSpec{adaptFromTask(mctruth)}, - Options{{"disable-qed", o2::framework::VariantType::Bool, false, {"disable QED handling"}}}}; + return DataProcessorSpec{.name = "IT3Digitizer", + .inputs = inputs, + .outputs = makeOutChannels(detOrig, mctruth), + .algorithm = AlgorithmSpec{adaptFromTask(mctruth)}, + .options = Options{{"disable-qed", o2::framework::VariantType::Bool, false, {"disable QED handling"}}}}; } } // namespace o2::its3 From 06000aab9ed91128baddbb4c99ac77d99a981d17 Mon Sep 17 00:00:00 2001 From: Felix Schlepper Date: Wed, 22 Apr 2026 21:02:20 +0200 Subject: [PATCH 3/3] ITS: fix truth seeding mode Signed-off-by: Felix Schlepper --- Detectors/ITSMFT/ITS/tracking/src/Vertexer.cxx | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/Detectors/ITSMFT/ITS/tracking/src/Vertexer.cxx b/Detectors/ITSMFT/ITS/tracking/src/Vertexer.cxx index cbff174634ec8..b827c134cf5b4 100644 --- a/Detectors/ITSMFT/ITS/tracking/src/Vertexer.cxx +++ b/Detectors/ITSMFT/ITS/tracking/src/Vertexer.cxx @@ -82,7 +82,6 @@ float Vertexer::clustersToVertices(LogFunc logger) // update LUT with all currently found vertices so in second iteration we can check vertPerROFThreshold sortVertices(); - mTimeFrame->updateROFVertexLookupTable(); } } catch (const BoundedMemoryResource::MemoryLimitExceeded& err) { handleException(err); @@ -126,6 +125,8 @@ void Vertexer::sortVertices() } mc.swap(sortedMC); } + // update LUT after sorting + mTimeFrame->updateROFVertexLookupTable(); } template