From 36a7d639c72b719c5a195889c79deeee632f48e3 Mon Sep 17 00:00:00 2001 From: Felix Schlepper Date: Fri, 15 May 2026 11:10:02 +0200 Subject: [PATCH] ITS3: opt. staggering and fixes to alignment Signed-off-by: Felix Schlepper --- .../ITS3/alignment/src/alignment-workflow.cxx | 2 + .../include/ITS3Reconstruction/Clusterer.h | 47 +-- .../include/ITS3Reconstruction/IOUtils.h | 1 + .../ITS3Reconstruction/TrackingInterface.h | 3 +- .../ITS3/reconstruction/src/Clusterer.cxx | 231 +++++++------ .../ITS3/reconstruction/src/IOUtils.cxx | 49 ++- .../reconstruction/src/TrackingInterface.cxx | 43 ++- .../include/ITS3Simulation/Digitizer.h | 17 +- .../ITS3/simulation/src/Digitizer.cxx | 87 ++--- .../ITS3TrackingStudyParam.h | 2 +- .../ITS3/study/macros/PlotMisalignment.C | 44 ++- .../Upgrades/ITS3/study/macros/PlotPulls.C | 28 +- .../Upgrades/ITS3/study/src/TrackingStudy.cxx | 144 ++++++--- .../src/its3-tracking-study-workflow.cxx | 2 + .../include/ITS3Workflow/ClustererSpec.h | 9 +- .../include/ITS3Workflow/DigitReaderSpec.h | 75 ++--- .../include/ITS3Workflow/DigitWriterSpec.h | 10 +- .../include/ITS3Workflow/RecoWorkflow.h | 1 + .../include/ITS3Workflow/TrackerSpec.h | 3 +- .../ITS3/workflow/src/ClustererSpec.cxx | 235 ++++++++++---- .../ITS3/workflow/src/DigitReaderSpec.cxx | 128 ++++---- .../ITS3/workflow/src/DigitWriterSpec.cxx | 97 +++--- .../ITS3/workflow/src/RecoWorkflow.cxx | 21 +- .../ITS3/workflow/src/TrackerSpec.cxx | 35 +- .../workflow/src/digit-writer-workflow.cxx | 7 +- .../ITS3/workflow/src/its3-reco-workflow.cxx | 5 +- .../src/ITS3DigitizerSpec.cxx | 303 +++++++++++------- .../DigitizerWorkflow/src/ITS3DigitizerSpec.h | 2 +- .../src/SimpleDigitizerWorkflow.cxx | 7 +- 29 files changed, 1004 insertions(+), 634 deletions(-) diff --git a/Detectors/Upgrades/ITS3/alignment/src/alignment-workflow.cxx b/Detectors/Upgrades/ITS3/alignment/src/alignment-workflow.cxx index 6ec7615885556..fd9495fb3f206 100644 --- a/Detectors/Upgrades/ITS3/alignment/src/alignment-workflow.cxx +++ b/Detectors/Upgrades/ITS3/alignment/src/alignment-workflow.cxx @@ -15,6 +15,7 @@ #include "GlobalTrackingWorkflowHelpers/InputHelper.h" #include "GlobalTrackingWorkflowHelpers/NoInpDummyOutSpec.h" #include "DetectorsRaw/HBFUtilsInitializer.h" +#include "DataFormatsITSMFT/DPLAlpideParamInitializer.h" #include "ITS3Align/AlignmentSpec.h" using namespace o2::framework; @@ -38,6 +39,7 @@ void customize(std::vector& workflowOptions) {"disable-root-input", VariantType::Bool, false, {"disable root-files input reader"}}, {"configKeyValues", VariantType::String, "", {"Semicolon separated key=value strings ..."}}}; o2::raw::HBFUtilsInitializer::addConfigOption(options); + o2::itsmft::DPLAlpideParamInitializer::addITSConfigOption(options); std::swap(workflowOptions, options); } #include "Framework/runDataProcessing.h" diff --git a/Detectors/Upgrades/ITS3/reconstruction/include/ITS3Reconstruction/Clusterer.h b/Detectors/Upgrades/ITS3/reconstruction/include/ITS3Reconstruction/Clusterer.h index a81db09217e9b..ab75c3fb1047b 100644 --- a/Detectors/Upgrades/ITS3/reconstruction/include/ITS3Reconstruction/Clusterer.h +++ b/Detectors/Upgrades/ITS3/reconstruction/include/ITS3Reconstruction/Clusterer.h @@ -20,6 +20,7 @@ // | *| #define _ALLOW_DIAGONAL_ALPIDE_CLUSTERS_ +#include #include #include #include @@ -94,18 +95,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); } }; @@ -121,6 +114,10 @@ class Clusterer }; struct ClustererThread { + struct PreCluster { + int head = 0; // index of precluster head in the pixels + int index = 0; + }; Clusterer* parent = nullptr; // parent clusterer int id = -1; // buffers for entries in preClusterIndices in 2 columns, to avoid boundary checks, we reserve @@ -133,12 +130,11 @@ class Clusterer // 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. + std::vector preClusters; //! preclusters info // /// temporary storage for the thread output CompClusCont compClusters; @@ -152,10 +148,11 @@ class Clusterer ///< 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 ///< 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]]; + auto& firstIndex = preClusters[preClusters[preClusIndex].index].head; pixels.emplace_back(firstIndex, ip); firstIndex = pixels.size() - 1; curr[row] = preClusIndex; @@ -164,11 +161,10 @@ class Clusterer ///< 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()); + int lastIndex = preClusters.size(); + preClusters.emplace_back(pixels.size(), lastIndex); // 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 } @@ -212,19 +208,25 @@ class Clusterer bool isContinuousReadOut() const { return mContinuousReadout; } void setContinuousReadOut(bool v) { mContinuousReadout = v; } + bool isDropHugeClusters() const { return mDropHugeClusters; } + void setDropHugeClusters(bool v) { mDropHugeClusters = 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; } + int getMaxROFDepthToSquash(int layer = -1) const { return (layer < 0) ? mSquashingDepth : mSquashingLayerDepth[layer]; } void setMaxROFDepthToSquash(int v) { mSquashingDepth = v; } + void addMaxROFDepthToSquash(int v) { mSquashingLayerDepth.push_back(v); } - int getMaxBCSeparationToSquash() const { return mMaxBCSeparationToSquash; } + int getMaxBCSeparationToSquash(int layer = -1) const { return (layer < 0) ? mMaxBCSeparationToSquash : mMaxBCSeparationToSquashLayer[layer]; } void setMaxBCSeparationToSquash(int n) { mMaxBCSeparationToSquash = n; } + void addMaxBCSeparationToSquash(int n) { mMaxBCSeparationToSquashLayer.push_back(n); } - void print() const; + void print(bool showsTiming) const; + void reset(); void clear(); ///< load the dictionary of cluster topologies @@ -249,6 +251,7 @@ class Clusterer // clusterization options bool mContinuousReadout = true; ///< flag continuous readout + bool mDropHugeClusters = false; ///< don't include clusters that would be split in more than one ///< 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); @@ -258,6 +261,8 @@ class Clusterer ///< Squashing options int mSquashingDepth = 0; ///< squashing is applied to next N rofs int mMaxBCSeparationToSquash = 6000. / o2::constants::lhc::LHCBunchSpacingNS + 10; + std::vector mSquashingLayerDepth; + std::vector mMaxBCSeparationToSquashLayer; std::vector> mThreads; // buffers for threads std::vector mChips; // currently processed ROF's chips data 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..23cbb0aad487c 100644 --- a/Detectors/Upgrades/ITS3/reconstruction/include/ITS3Reconstruction/TrackingInterface.h +++ b/Detectors/Upgrades/ITS3/reconstruction/include/ITS3Reconstruction/TrackingInterface.h @@ -24,10 +24,11 @@ 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 overrideParameters(std::vector& t, std::vector& v) final; + 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/Clusterer.cxx b/Detectors/Upgrades/ITS3/reconstruction/src/Clusterer.cxx index bce17b3759340..0b43a7cfea693 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. // @@ -10,7 +10,7 @@ // or submit itself to any jurisdiction. /// \file Clusterer.cxx -/// \brief Implementation of the ITS cluster finder +/// \brief Implementation of the ITS3 cluster finder #include @@ -19,8 +19,6 @@ #include "SimulationDataFormat/MCTruthContainer.h" #include "CommonDataFormat/InteractionRecord.h" -#include "TTree.h" - #ifdef WITH_OPENMP #include #endif @@ -32,12 +30,11 @@ void Clusterer::process(int nThreads, PixelReader& reader, CompClusCont* compClu PatternCont* patterns, ROFRecCont* vecROFRec, MCTruth* labelsCl) { #ifdef _PERFORM_TIMING_ - mTimer.Start(false); + mTimer.Start(kFALSE); #endif - if (nThreads < 1) { - nThreads = 1; - } + nThreads = std::max(nThreads, 1); auto autoDecode = reader.getDecodeNextAuto(); + o2::InteractionRecord lastIR{}; do { if (autoDecode) { reader.setDecodeNextAuto(false); // internally do not autodecode @@ -48,6 +45,15 @@ void Clusterer::process(int nThreads, PixelReader& reader, CompClusCont* compClu 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(); @@ -66,9 +72,7 @@ void Clusterer::process(int nThreads, PixelReader& reader, CompClusCont* compClu } break; // just 1 ROF was asked to be processed } - if (nFired < nThreads) { - nThreads = nFired; - } + nThreads = std::min(nFired, nThreads); #ifndef WITH_OPENMP nThreads = 1; #endif @@ -77,7 +81,7 @@ void Clusterer::process(int nThreads, PixelReader& reader, CompClusCont* compClu if (nThreads > mThreads.size()) { int oldSz = mThreads.size(); mThreads.resize(nThreads); - for (size_t i = oldSz; i < nThreads; i++) { + for (int i = oldSz; i < nThreads; i++) { mThreads[i] = std::make_unique(this, i); } } @@ -106,10 +110,10 @@ void Clusterer::process(int nThreads, PixelReader& reader, CompClusCont* compClu mTimerMerge.Start(false); #endif size_t nClTot = 0, nPattTot = 0; - int chid = 0; - std::vector thrStatIdx(nThreads, 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(); } @@ -126,14 +130,16 @@ void Clusterer::process(int nThreads, PixelReader& reader, CompClusCont* compClu 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); + 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); + } } } } @@ -163,7 +169,7 @@ void Clusterer::ClustererThread::process(uint16_t chip, uint16_t nChips, CompClu 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}); + 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]; @@ -171,7 +177,7 @@ void Clusterer::ClustererThread::process(uint16_t chip, uint16_t nChips, CompClu 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]); + parent->mMaxRowColDiffToMask ? curChipData->maskFiredInSample(parent->mChipsOld[chipID], parent->mMaxRowColDiffToMask) : curChipData->maskFiredInSample(parent->mChipsOld[chipID]); } } auto validPixID = curChipData->getFirstUnmasked(); @@ -205,14 +211,22 @@ void Clusterer::ClustererThread::finishChip(ChipPixelData* curChipData, CompClus 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]; + 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 = preClusterHeads[i1]; + int next = preCluster.head; pixArrBuff.clear(); while (next >= 0) { const auto& pixEntry = pixels[next]; @@ -228,18 +242,19 @@ void Clusterer::ClustererThread::finishChip(ChipPixelData* curChipData, CompClus } next = pixEntry.first; } - preClusterIndices[i1] = -1; - for (size_t i2 = i1 + 1; i2 < preClusterHeads.size(); ++i2) { - if (preClusterIndices[i2] != ci) { + preCluster.index = -1; + for (int i2 = i1 + 1; i2 < nPreclusters; ++i2) { + auto& preCluster2 = preClusters[i2]; + if (preCluster2.index != ci) { continue; } - next = preClusterHeads[i2]; + 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 != nullptr) { + 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 @@ -248,42 +263,46 @@ void Clusterer::ClustererThread::finishChip(ChipPixelData* curChipData, CompClus } next = pixEntry.first; } - preClusterIndices[i2] = -1; + preCluster2.index = -1; } if (bbox.isAcceptableSize()) { - parent->streamCluster(pixArrBuff, &labelsBuff, bbox, parent->mPattIdConverter, compClusPtr, patternsPtr, labelsClusPtr, nlab, constants::detID::isDetITS3(curChipData->getChipID())); + const bool isIB = constants::detID::isDetITS3(curChipData->getChipID()); + parent->streamCluster(pixArrBuff, &labelsBuff, bbox, parent->mPattIdConverter, compClusPtr, patternsPtr, labelsClusPtr, nlab, isIB); } 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)" : ""); + 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++; + { + 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); + 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); + if (!pixbuf.empty()) { // Stream a piece of cluster only if the reduced bounding box is not empty + const bool isIB = constants::detID::isDetITS3(curChipData->getChipID()); + parent->streamCluster(pixbuf, &labelsBuff, bboxT, parent->mPattIdConverter, compClusPtr, patternsPtr, labelsClusPtr, nlab, isIB, true); + pixbuf.clear(); + } + bboxT.rowMin = bboxT.rowMax + 1; + } while (bboxT.rowMin < bbox.rowMax); + bboxT.colMin = bboxT.colMax + 1; + } while (bboxT.colMin < bbox.colMax); + } } } } @@ -350,14 +369,12 @@ void Clusterer::ClustererThread::initChip(const ChipPixelData* curChipData, uint resetColumn(curr); pixels.clear(); - preClusterHeads.clear(); - preClusterIndices.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 - preClusterHeads.push_back(0); - preClusterIndices.push_back(0); + preClusters.emplace_back(); pixels.emplace_back(-1, first); // id of current pixel noLeftCol = true; // flag that there is no column on the left to check yet } @@ -382,39 +399,55 @@ void Clusterer::ClustererThread::updateChip(const ChipPixelData* curChipData, ui 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 { + addNewPrecluster(ip, row); // start new precluster } } else { + // row above should be always checked + int nnb = 0, lowestIndex = curr[row - 1], *nbrCol[4], nbrRow[4]; + if (lowestIndex >= 0) { + nbrCol[nnb] = curr; + nbrRow[nnb++] = row - 1; + } else { + lowestIndex = 0x7ffff; + } #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; + 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; + } + nnb++; } - if (orphan) { - expandPreCluster(ip, row, pci); // attach to the adjascent precluster - orphan = false; - continue; + } +#else + if (prev[row] >= 0) { + nbrCol[nnb] = prev; + nbrRow[nnb] = row; + if (prev[row] < lowestIndex) { + lowestIndex = prev[row]; } - // reassign precluster index to smallest one - if (preClusterIndices[pci] < preClusterIndices[curr[row]]) { - preClusterIndices[curr[row]] = preClusterIndices[pci]; - } else { - preClusterIndices[pci] = preClusterIndices[curr[row]]; + 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; + } } } } - if (orphan) { - addNewPrecluster(ip, row); // start new precluster - } } //__________________________________________________ @@ -453,20 +486,30 @@ void Clusterer::clear() } //__________________________________________________ -void Clusterer::print() const +void Clusterer::print(bool showsTiming) 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); + // 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"; + 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 diff --git a/Detectors/Upgrades/ITS3/reconstruction/src/IOUtils.cxx b/Detectors/Upgrades/ITS3/reconstruction/src/IOUtils.cxx index 92e36cd2a4b84..1a7bb8adb1d58 100644 --- a/Detectors/Upgrades/ITS3/reconstruction/src/IOUtils.cxx +++ b/Detectors/Upgrades/ITS3/reconstruction/src/IOUtils.cxx @@ -61,27 +61,35 @@ 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); - its::bounded_vector clusterSizeVec(clusters.size(), tf->getMemoryPool().get()); + // check for missing/empty/unset rofs + // the code requires consistent monotonically increasing input without gaps + const auto& timing = tf->getROFOverlapTableView().getLayer(layer >= 0 ? layer : 0); + 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(), 0, tf->getMemoryPool().get()); for (size_t iRof{0}; iRof < rofs.size(); ++iRof) { const auto& rof = rofs[iRof]; for (int clusterId{rof.getFirstEntry()}; clusterId < rof.getFirstEntry() + rof.getNEntries(); ++clusterId) { const auto& c = clusters[clusterId]; const auto sensorID = c.getSensorID(); - const auto layer = geom->getLayer(sensorID); + const auto lay = geom->getLayer(sensorID); float sigmaY2{0}, sigmaZ2{0}, sigmaYZ{0}; uint8_t clusterSize{0}; const auto locXYZ = extractClusterData(c, pattIt, dict, sigmaY2, sigmaZ2, clusterSize); - clusterSizeVec.push_back(clusterSize); + clusterSizeVec[clusterId] = clusterSize; // Transformation to the local --> global const auto gloXYZ = geom->getMatrixL2G(sensorID) * locXYZ; @@ -102,30 +110,37 @@ int loadROFrameDataITS3(its::TimeFrame<7>* tf, } math_utils::detail::bringToPMPi(alpha); // alpha is defined on -Pi,Pi - tf->addTrackingFrameInfoToLayer(layer, gloXYZ.x(), gloXYZ.y(), gloXYZ.z(), x, alpha, + tf->addTrackingFrameInfoToLayer(lay, gloXYZ.x(), gloXYZ.y(), gloXYZ.z(), x, alpha, std::array{y, trkXYZ.z()}, std::array{sigmaY2, sigmaYZ, sigmaZ2}); /// Rotate to the global frame - tf->addClusterToLayer(layer, gloXYZ.x(), gloXYZ.y(), gloXYZ.z(), tf->getUnsortedClusters()[layer].size()); - tf->addClusterExternalIndexToLayer(layer, clusterId); + tf->addClusterToLayer(lay, gloXYZ.x(), gloXYZ.y(), gloXYZ.z(), tf->getUnsortedClusters()[lay].size()); + tf->addClusterExternalIndexToLayer(lay, clusterId); } - for (unsigned int iL{0}; iL < tf->getUnsortedClusters().size(); ++iL) { - tf->mROFramesClusters[iL][iRof + 1] = (int)tf->getUnsortedClusters()[iL].size(); + // effectively calculating an exclusive sum + if (layer >= 0) { + tf->mROFramesClusters[layer][iRof + 1] = tf->mUnsortedClusters[layer].size(); + } else { + for (unsigned int iL{0}; iL < tf->mUnsortedClusters.size(); ++iL) { + tf->mROFramesClusters[iL][iRof + 1] = tf->mUnsortedClusters[iL].size(); + } } } - // tf->setClusterSize(clusterSizeVec); FIXME + tf->setClusterSize(layer >= 0 ? layer : 0, 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 || layer == -1) { + for (auto i = 0; i < tf->mNTrackletsPerCluster.size(); ++i) { + tf->mNTrackletsPerCluster[i].resize(tf->mUnsortedClusters[1].size()); + tf->mNTrackletsPerClusterSum[i].resize(tf->mUnsortedClusters[1].size() + 1); + } } if (mcLabels != nullptr) { - // tf->mClusterLabels = mcLabels; // FIXME + tf->mClusterLabels[layer >= 0 ? layer : 0] = mcLabels; + } else { + tf->mClusterLabels[layer >= 0 ? layer : 0] = nullptr; } return 0; } diff --git a/Detectors/Upgrades/ITS3/reconstruction/src/TrackingInterface.cxx b/Detectors/Upgrades/ITS3/reconstruction/src/TrackingInterface.cxx index 9fe6f3735a845..1ebaf3971c534 100644 --- a/Detectors/Upgrades/ITS3/reconstruction/src/TrackingInterface.cxx +++ b/Detectors/Upgrades/ITS3/reconstruction/src/TrackingInterface.cxx @@ -9,41 +9,34 @@ // granted to it by virtue of its status as an Intergovernmental Organization // or submit itself to any jurisdiction. +#include "ITS3Base/SpecsV2.h" #include "ITS3Reconstruction/TrackingInterface.h" #include "ITS3Reconstruction/IOUtils.h" #include "ITSBase/GeometryTGeo.h" -#include "ITStracking/TrackingConfigParam.h" #include "DataFormatsITSMFT/DPLAlpideParam.h" #include "DetectorsBase/GRPGeomHelper.h" -#include "Framework/DeviceSpec.h" namespace o2::its3 { -void ITS3TrackingInterface::updateTimeDependentParams(framework::ProcessingContext& pc) +void ITS3TrackingInterface::overrideParameters(std::vector& t, std::vector& v) { - 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()); - } - } + // only override IB radii + for (auto& tt : t) { + tt.LayerRadii[0] = constants::radii[0]; + tt.LayerRadii[1] = constants::radii[1]; + tt.LayerRadii[2] = constants::radii[2]; } + for (auto& vv : v) { + vv.LayerRadii[0] = constants::radii[0]; + vv.LayerRadii[1] = constants::radii[1]; + vv.LayerRadii[2] = constants::radii[2]; + } +} + +void ITS3TrackingInterface::requestTopologyDictionary(framework::ProcessingContext& pc) +{ + pc.inputs().get("itscldict"); // just to trigger the finaliseCCDB } void ITS3TrackingInterface::finaliseCCDB(framework::ConcreteDataMatcher& matcher, void* obj) @@ -80,7 +73,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..6fbed92bb1400 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,23 @@ 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) { + if (layer < 0) { + return true; + } + 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 +126,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 +154,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 +173,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() || (layer >= 0 && mGeometry->getLayer(chip.getChipIndex()) != layer)) { + continue; + } chip.addNoise(mROFrameMin, mROFrameMin, &mParams); auto& buffer = chip.getPreDigits(); if (buffer.empty()) { @@ -188,7 +205,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 +219,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 lay) { // convert single hit to digits auto chipID = hit.GetDetectorID(); @@ -223,25 +240,27 @@ 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(lay)); // 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(lay) : roFrameRel; int nFrames = roFrameRelMax + 1 - roFrameRel; uint32_t roFrameMax = mNewROFrame + roFrameRelMax; - if (roFrameMax > maxFr) { - maxFr = roFrameMax; // if signal extends beyond current maxFrame, increase the latter - } + maxFr = std::max(roFrameMax, maxFr); // if signal extends beyond current maxFrame, increase the latter // 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); + const int layer = mGeometry->getLayer(chipID); + const auto& matrix = mGeometry->getMatrixL2G(chipID); int nSteps = chip.isIB() ? mParams.getIBNSimSteps() : mParams.getNSimSteps(); float nStepsInv = chip.isIB() ? mParams.getIBNSimStepsInv() : mParams.getNSimStepsInv(); math_utils::Vector3D xyzLocS, xyzLocE; @@ -404,34 +423,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, lay); } } } 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/study/include/ITS3TrackingStudy/ITS3TrackingStudyParam.h b/Detectors/Upgrades/ITS3/study/include/ITS3TrackingStudy/ITS3TrackingStudyParam.h index 1150dc19a1162..db68ff7ec1b98 100644 --- a/Detectors/Upgrades/ITS3/study/include/ITS3TrackingStudy/ITS3TrackingStudyParam.h +++ b/Detectors/Upgrades/ITS3/study/include/ITS3TrackingStudy/ITS3TrackingStudyParam.h @@ -36,7 +36,7 @@ struct ITS3TrackingStudyParam : o2::conf::ConfigurableParamHelper::MatCorrType CorrType = o2::base::PropagatorImpl::MatCorrType::USEMatCorrLUT; + o2::base::PropagatorImpl::MatCorrType CorrType = o2::base::PropagatorImpl::MatCorrType::USEMatCorrTGeo; /// studies bool doDCA = false; diff --git a/Detectors/Upgrades/ITS3/study/macros/PlotMisalignment.C b/Detectors/Upgrades/ITS3/study/macros/PlotMisalignment.C index 5b9372bc3402e..3f2170776e1fd 100644 --- a/Detectors/Upgrades/ITS3/study/macros/PlotMisalignment.C +++ b/Detectors/Upgrades/ITS3/study/macros/PlotMisalignment.C @@ -57,7 +57,7 @@ void processTree(TFile* f, const char* treeName) } // branch variables - float dY, dZ, phi, eta; + float dY, dZ, phi, eta, dcaXY, dcaZ; int lay; auto* trk = new o2::track::TrackParCov; auto* mcTrk = new o2::track::TrackPar; @@ -67,8 +67,15 @@ void processTree(TFile* f, const char* treeName) tree->SetBranchAddress("eta", &eta); tree->SetBranchAddress("lay", &lay); tree->SetBranchAddress("trk", &trk); + tree->SetBranchAddress("dcaXY", &dcaXY); + tree->SetBranchAddress("dcaZ", &dcaZ); tree->SetBranchAddress("mcTrk", &mcTrk); + const int nPtBins = 35; + const double ptLimits[nPtBins] = {0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5, 0.6, 0.7, 0.8, 0.9, 1., 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2., 2.2, 2.5, 3., 4., 5., 6., 8., 10., 15., 20.}; + const int yDCABins{1000}; + const float yDCARange{500}; + // --- book histograms --- // [ptBin][lay] for each plot type // dY/dZ vs phi @@ -79,6 +86,12 @@ void processTree(TFile* f, const char* treeName) TProfile2D* hProf[kNVar][kNPtBins][kNLay]; // pulls TH1F* hPull[kNPar][kNPtBins][kNLay]; + // DCAxy + TH2F* hDCAxyVsPt = new TH2F(Form("%s_hDCAxyVsPt", treeName), ";#it{p}_{T,MC} (GeV/#it{c});DCA_{#it{xy}} (#mum);entries", nPtBins - 1, ptLimits, yDCABins, -yDCARange, yDCARange); + TH2F* hDCAxyVsPhi = new TH2F(Form("%s_hDCAxyVsPhi", treeName), ";#phi (rad);DCA_{#it{xy}} (#mum);entries", 100, 0, 2 * TMath::Pi(), yDCABins, -yDCARange, yDCARange); + // DCAz + TH2F* hDCAzVsPt = new TH2F(Form("%s_hDCAzVsPt", treeName), ";#it{p}_{T,MC} (GeV/#it{c});DCA_{#it{z}} (#mum);entries", nPtBins - 1, ptLimits, yDCABins, -yDCARange, yDCARange); + TH2F* hDCAzVsPhi = new TH2F(Form("%s_hDCAzVsPhi", treeName), ";#phi (rad);DCA_{#it{z}} (#mum);entries", 100, 0, 2 * TMath::Pi(), yDCABins, -yDCARange, yDCARange); for (int ipt = 0; ipt < kNPtBins; ipt++) { for (int ilay = 0; ilay < kNLay; ilay++) { @@ -118,6 +131,15 @@ void processTree(TFile* f, const char* treeName) float dYum = dY * 10000.f; float dZum = dZ * 10000.f; + if (lay == -1) { + hDCAxyVsPt->Fill(pt, dcaXY * 10000.); + hDCAzVsPt->Fill(pt, dcaZ * 10000.); + if (pt >= 1.0 && pt <= 2.0) { + hDCAxyVsPhi->Fill(phi, dcaXY * 10000.); + hDCAzVsPhi->Fill(phi, dcaZ * 10000.); + } + } + // integrated (ipt=0) + differential int iptDiff = getPtBin(pt); for (int ipt : {0, iptDiff}) { @@ -216,6 +238,26 @@ void processTree(TFile* f, const char* treeName) c->SaveAs(Form("%s.png", c->GetName())); } } + + // write file out + auto oFile = TFile::Open(Form("plotMisalignment_%s.root", treeName), "RECREATE"); + hDCAxyVsPt->Write(); + hDCAzVsPt->Write(); + hDCAxyVsPhi->Write(); + hDCAzVsPhi->Write(); + for (int ipt = 0; ipt < kNPtBins; ipt++) { + for (int ilay = 0; ilay < kNLay; ilay++) { + for (int iv = 0; iv < kNVar; iv++) { + hVsPhi[iv][ipt][ilay]->Write(); + hVsEta[iv][ipt][ilay]->Write(); + hProf[iv][ipt][ilay]->Write(); + } + for (auto& ip : hPull) { + ip[ipt][ilay]->Write(); + } + } + } + oFile->Close(); } void PlotMisalignment(const char* fname = "its3TrackStudy.root") diff --git a/Detectors/Upgrades/ITS3/study/macros/PlotPulls.C b/Detectors/Upgrades/ITS3/study/macros/PlotPulls.C index 371a94cda0e70..bc59261199065 100644 --- a/Detectors/Upgrades/ITS3/study/macros/PlotPulls.C +++ b/Detectors/Upgrades/ITS3/study/macros/PlotPulls.C @@ -9,11 +9,14 @@ // granted to it by virtue of its status as an Intergovernmental Organization // or submit itself to any jurisdiction. -/// \file PlotDCA.C -/// \brief Simple macro to plot ITS3 impact parameter resolution +/// \file PlotPulls.C +/// \brief Simple macro to plot ITS3 pulls #if !defined(__CLING__) || defined(__ROOTCLING__) +#include +#include #include +#include #include #include @@ -27,7 +30,6 @@ #include "ReconstructionDataFormats/GlobalTrackID.h" #include "ReconstructionDataFormats/Track.h" -#include "ReconstructionDataFormats/DCA.h" #include "SimulationDataFormat/MCTrack.h" #endif @@ -49,7 +51,13 @@ void PlotPulls(const char* fName = "its3TrackStudy.root") { TH1::SetDefaultSumw2(); std::unique_ptr inFile(TFile::Open(fName)); + if (!inFile || inFile->IsZombie()) { + return; + } auto tree = inFile->Get("pull"); + if (!tree) { + return; + } uint8_t src; // track type tree->SetBranchAddress("src", &src); @@ -133,13 +141,13 @@ void PlotPulls(const char* fName = "its3TrackStudy.root") std::vector projs; const char* fitOpt{"QWMERSB"}; for (int i{0}; i < o2::track::kNParams; ++i) { - for (auto iPt{0}; iPt < nPtBins; ++iPt) { + for (auto iPt{0}; iPt < nPtBins - 1; ++iPt) { auto hProj = pulls[i]->ProjectionY(Form("%s_%d", pulls[i]->GetName(), iPt), iPt + 1, iPt + 1); hProj->SetName(Form("p%s_pt%d", pNames[i], iPt)); hProj->SetTitle(Form("Pull %s #it{p}_{T}#in[%.2f, %.2f)", pNames[i], ptLimits[iPt], ptLimits[iPt + 1])); projs.push_back(hProj); if (hProj->GetEntries() < 100) { - return; + continue; } fGaus->SetParameter(1, 0); fGaus->SetParameter(2, 1); @@ -153,14 +161,18 @@ void PlotPulls(const char* fName = "its3TrackStudy.root") } } - hMahDist2->Scale(1. / hMahDist2->Integral("width")); TF1* fchi2Fit = new TF1("fchi2_fit", chi2_pdf, 0.1, 6, 3); fchi2Fit->SetParNames("A", "k", "s"); fchi2Fit->SetParameter(0, 1); fchi2Fit->SetParameter(1, 5); fchi2Fit->SetParameter(2, 1); - auto fitres = hMahDist2->Fit(fchi2Fit, "RMQS"); - fitres->Print(); + if (hMahDist2->Integral("width") > 0.) { + hMahDist2->Scale(1. / hMahDist2->Integral("width")); + auto fitres = hMahDist2->Fit(fchi2Fit, "RMQS"); + if (fitres.Get()) { + fitres->Print(); + } + } TFile outFile("plotPulls.root", "RECREATE"); for (int i{0}; i < o2::track::kNParams; ++i) { diff --git a/Detectors/Upgrades/ITS3/study/src/TrackingStudy.cxx b/Detectors/Upgrades/ITS3/study/src/TrackingStudy.cxx index 63a08d1c9c6ab..9bdcda8b77a4c 100644 --- a/Detectors/Upgrades/ITS3/study/src/TrackingStudy.cxx +++ b/Detectors/Upgrades/ITS3/study/src/TrackingStudy.cxx @@ -66,8 +66,8 @@ class TrackingStudySpec : public Task TrackingStudySpec(TrackingStudySpec&&) = delete; TrackingStudySpec& operator=(const TrackingStudySpec&) = delete; TrackingStudySpec& operator=(TrackingStudySpec&&) = delete; - TrackingStudySpec(std::shared_ptr dr, std::shared_ptr gr, GTrackID::mask_t src, bool useMC) - : mDataRequest(dr), mGGCCDBRequest(gr), mTracksSrc(src), mUseMC(useMC) {} + TrackingStudySpec(std::shared_ptr dr, std::shared_ptr gr, GTrackID::mask_t src, bool useMC, bool withPV) + : mDataRequest(dr), mGGCCDBRequest(gr), mTracksSrc(src), mUseMC(useMC), mWithPV(withPV) {} ~TrackingStudySpec() final = default; void init(InitContext& ic) final; void run(ProcessingContext& pc) final; @@ -80,7 +80,8 @@ class TrackingStudySpec : public Task void prepareITSClusters(); bool selectTrack(GTrackID trkID, bool checkMCTruth = true) const; T2VMap buildT2V(bool includeCont = false, bool requireMCMatch = true) const; - bool refitITSPVTrack(o2::track::TrackParCov& trFit, GTrackID gidx); + bool refitITSPVTrack(o2::track::TrackParCov& trFit, GTrackID gidx, const o2::dataformats::VertexBase& pv); + void getImpactParams(const o2::track::TrackParCov& trk, const o2::dataformats::VertexBase& pv, float ip[2], float bz); void doDCAStudy(); void doDCARefitStudy(); @@ -146,6 +147,7 @@ class TrackingStudySpec : public Task std::shared_ptr mDataRequest; std::shared_ptr mGGCCDBRequest; bool mUseMC{false}; + bool mWithPV{false}; GTrackID::mask_t mTracksSrc; o2::vertexing::PVertexer mVertexer; o2::steer::MCKinematicsReader mMCReader; // reader of MC information @@ -183,6 +185,7 @@ void TrackingStudySpec::updateTimeDependentParams(ProcessingContext& pc) mVertexer.init(); o2::its::GeometryTGeo::Instance()->fillMatrixCache(o2::math_utils::bit2Mask(o2::math_utils::TransformType::T2L, o2::math_utils::TransformType::L2G, o2::math_utils::TransformType::T2G)); mParams = &ITS3TrackingStudyParam::Instance(); + mParams->printKeyValues(true, true); if (mParams->doMisalignment) { mMisalignment = {}; if (!mParams->misAlgJson.empty()) { @@ -381,36 +384,19 @@ T2VMap TrackingStudySpec::buildT2V(bool includeCont, bool requireMCMatch) const return std::move(t2v); } -bool TrackingStudySpec::refitITSPVTrack(o2::track::TrackParCov& trFit, GTrackID gidx) +bool TrackingStudySpec::refitITSPVTrack(o2::track::TrackParCov& trFit, GTrackID gidx, const o2::dataformats::VertexBase& pv) { if (gidx.getSource() != GTrackID::ITS) { return false; } - static auto pvvec = mRecoData.getPrimaryVertices(); - static auto t2v = buildT2V(true, true); - static std::vector itsTracksROF; - if (static bool done{false}; !done) { - done = true; - const auto& itsTracksROFRec = mRecoData.getITSTracksROFRecords(); - itsTracksROF.resize(mRecoData.getITSTracks().size()); - for (unsigned irf = 0, cnt = 0; irf < itsTracksROFRec.size(); irf++) { - int ntr = itsTracksROFRec[irf].getNEntries(); - for (int itr = 0; itr < ntr; itr++) { - itsTracksROF[cnt++] = irf; - } - } - } - auto prop = o2::base::Propagator::Instance(); + + const auto geom = o2::its::GeometryTGeo::Instance(); + const auto prop = o2::base::Propagator::Instance(); std::array clArr{nullptr}; const auto trkIn = mRecoData.getTrackParam(gidx); - const auto trkOut = mRecoData.getTrackParamOut(gidx); const auto& itsTrOrig = mRecoData.getITSTrack(gidx); - int ncl = itsTrOrig.getNumberOfClusters(), rof = itsTracksROF[gidx.getIndex()]; - const auto& itsTrackClusRefs = mRecoData.getITSTracksClusterRefs(); - int clEntry = itsTrOrig.getFirstClusterEntry(); - const auto propagator = o2::base::Propagator::Instance(); + // convert PV to a fake cluster in the track DCA frame - const auto& pv = pvvec[t2v[gidx]]; auto trkPV = trkIn; if (!prop->propagateToDCA(pv, trkPV, prop->getNominalBz(), 2.0, mParams->CorrType)) { mTrackCounter -= gidx.getSource(); @@ -421,18 +407,29 @@ bool TrackingStudySpec::refitITSPVTrack(o2::track::TrackParCov& trFit, GTrackID o2::math_utils::sincos(trkPV.getAlpha(), sinAlp, cosAlp); // vertex position rotated to track frame TrackingCluster pvCls; + pvCls.alpha = trkPV.getAlpha(); pvCls.setXYZ((pv.getX() * cosAlp) + (pv.getY() * sinAlp), (-pv.getX() * sinAlp) + (pv.getY() * cosAlp), pv.getZ()); - pvCls.setSigmaY2(0.5f * (pv.getSigmaX2() + pv.getSigmaY2())); - pvCls.setSigmaZ2(pv.getSigmaZ2()); + pvCls.setErrors(0.5f * (pv.getSigmaX2() + pv.getSigmaY2()), pv.getSigmaZ2(), 0.f); + pvCls.setSensorID(-1); clArr[0] = &pvCls; - for (int icl = 0; icl < ncl; ++icl) { // ITS clusters are referred in layer decreasing order - clArr[ncl - icl] = &mITScl[itsTrackClusRefs[clEntry + icl]]; + + const int ncl = itsTrOrig.getNumberOfClusters(); + for (int icl = 0; icl < ncl; ++icl) { + const auto& curClu = mITScl[mITSclRef[itsTrOrig.getClusterEntry(icl)]]; + const int llr = geom->getLayer(curClu.getSensorID()); + if (clArr[1 + llr]) { + LOGP(fatal, "Cluster at lr {} was already assigned, old sens {}, new sens {}", llr, clArr[1 + llr]->getSensorID(), curClu.getSensorID()); + } + clArr[1 + llr] = &curClu; } - // start refit - trFit = trkOut; + + trFit = mRecoData.getTrackParamOut(gidx); trFit.resetCovariance(1'000); float chi2{0}; - for (int icl = ncl; icl >= 0; --icl) { // go backwards + for (int icl = clArr.size() - 1; icl >= 0; --icl) { // go backwards + if (!clArr[icl]) { + continue; + } if (!trFit.rotate(clArr[icl]->alpha) || !prop->propagateToX(trFit, clArr[icl]->getX(), prop->getNominalBz(), 0.85, 2.0, mParams->CorrType)) { mTrackCounter -= gidx.getSource(); return false; @@ -443,7 +440,6 @@ bool TrackingStudySpec::refitITSPVTrack(o2::track::TrackParCov& trFit, GTrackID return false; } } - // chi2 < conf.maxChi2; should I cut here? return true; }; @@ -496,10 +492,10 @@ void TrackingStudySpec::doDCAStudy() const auto& trk = mRecoData.getTrackParam(contributorsGID[cis]); auto trkRefit = trk; // for ITS standalone tracks instead of having the trk at the pv we refit with the pv - if (mParams->refitITS && cis == GTrackID::ITS && !refitITSPVTrack(trkRefit, contributorsGID[cis])) { + if (mWithPV && mParams->refitITS && cis == GTrackID::ITS && !refitITSPVTrack(trkRefit, contributorsGID[cis], pv)) { mTrackCounter -= cis; continue; - } else { + } else if (!(mWithPV && mParams->refitITS && cis == GTrackID::ITS)) { trkRefit.invalidate(); }; @@ -684,17 +680,30 @@ void TrackingStudySpec::doPullStudy() mTrackCounter &= trkID.getSource(); return; } - const auto mcTrk = mMCReader.getTrack(mRecoData.getTrackMCLabel(trkID)); + const auto& lbl = mRecoData.getTrackMCLabel(trkID); + const auto mcTrk = mMCReader.getTrack(lbl); if (!mcTrk) { return; } + if (!mcTrk->isPrimary()) { + mTrackCounter &= trkID.getSource(); + return; + } auto trk = mRecoData.getTrackParam(trkID); - // for ITS standalone tracks we add the PV as an additional measurement point - if (mParams->refitITS && trkID.getSource() == GTrackID::ITS && !refitITSPVTrack(trk, trkID)) { - mTrackCounter -= trkID.getSource(); - ++nPullsFail; - return; + // for ITS standalone tracks we add the MC event vertex as an additional measurement point + if (mParams->refitITS && trkID.getSource() == GTrackID::ITS) { + const auto& eve = mMCReader.getMCEventHeader(lbl.getSourceID(), lbl.getEventID()); + o2::dataformats::VertexBase mcEve; + mcEve.setXYZ((float)eve.GetX(), (float)eve.GetY(), (float)eve.GetZ()); + mcEve.setSigmaX(20e-4f); + mcEve.setSigmaY(20e-4f); + mcEve.setSigmaZ(20e-4f); + if (!refitITSPVTrack(trk, trkID, mcEve)) { + mTrackCounter -= trkID.getSource(); + ++nPullsFail; + return; + } } std::array xyz{(float)mcTrk->GetStartVertexCoordinatesX(), (float)mcTrk->GetStartVertexCoordinatesY(), (float)mcTrk->GetStartVertexCoordinatesZ()}, @@ -870,9 +879,10 @@ void TrackingStudySpec::doResidStudy() auto doRefits = [&](const o2::its::TrackITS& iTrack, const o2::MCCompLabel& lbl) { std::array cl; std::array clArr{nullptr}; + dataformats::VertexBase pv; + float ip[2]; if (mParams->addPVAsCluster) { const auto& eve = mMCReader.getMCEventHeader(lbl.getSourceID(), lbl.getEventID()); - dataformats::VertexBase pv; auto trFitOut = iTrack.getParamIn(); pv.setXYZ(eve.GetX(), eve.GetY(), eve.GetZ()); if (!prop->propagateToDCA(pv, trFitOut, bz, base::Propagator::MAX_STEP, mParams->CorrType)) { @@ -885,8 +895,7 @@ void TrackingStudySpec::doResidStudy() o2::math_utils::sincos(trFitOut.getAlpha(), sinAlp, cosAlp); cl[0].alpha = trFitOut.getAlpha(); cl[0].setXYZ((pv.getX() * cosAlp) + (pv.getY() * sinAlp), (-pv.getX() * sinAlp) + (pv.getY() * cosAlp), pv.getZ()); - cl[0].setSigmaY2(0.5f * (pv.getSigmaX2() + pv.getSigmaY2())); - cl[0].setSigmaZ2(pv.getSigmaZ2()); + cl[0].setErrors(0.5f * (pv.getSigmaX2() + pv.getSigmaY2()), pv.getSigmaZ2(), 0.f); cl[0].setSensorID(-1); clArr[0] = &cl[0]; } @@ -918,6 +927,9 @@ void TrackingStudySpec::doResidStudy() } auto phi = i == 0 ? tInt.getPhi() : tInt.getPhiPos(); o2::math_utils::bringTo02Pi(phi); + if (clArr[0]) { + getImpactParams(tInt, pv, ip, bz); + } (*mDBGOut) << "res" << "dYInt=" << clArr[i]->getY() - tInt.getY() << "dZInt=" << clArr[i]->getZ() - tInt.getZ() @@ -932,6 +944,8 @@ void TrackingStudySpec::doResidStudy() << "alpha=" << clArr[i]->alpha << "sens=" << clArr[i]->getSensorID() << "phi=" << phi + << "dcaXY=" << ip[0] + << "dcaZ=" << ip[1] << "pt=" << tInt.getPt() << "chip=" << constants::detID::getSensorID(clArr[i]->getSensorID()) << "lay=" << i - 1 @@ -964,7 +978,8 @@ void TrackingStudySpec::doMisalignmentStudy() const auto geom = o2::its::GeometryTGeo::Instance(); int goodRefit{0}, notPassedSel{0}, fitFail{0}, fitFailMis{0}; - + o2::dataformats::VertexBase pv; + float ip[2]; float chi2{0}; auto writeTree = [&](const char* treeName, const std::array& clArr, @@ -994,12 +1009,17 @@ void TrackingStudySpec::doMisalignmentStudy() if (mcTrkAtX.rotate(tInt.getAlpha()) && prop->PropagateToXBxByBz(mcTrkAtX, tInt.getX())) { auto phi = i == 0 ? tInt.getPhi() : tInt.getPhiPos(); o2::math_utils::bringTo02Pi(phi); + if (clArr[0]) { + getImpactParams(tInt, pv, ip, prop->getNominalBz()); + } (*mDBGOut) << treeName << "trk=" << tInt << "mcTrk=" << mcTrkAtX << "chi2=" << chi2 << "dY=" << dY << "dZ=" << dZ + << "dcaXY=" << ip[0] + << "dcaZ=" << ip[1] << "phi=" << phi << "eta=" << tInt.getEta() << "lay=" << i - 1 @@ -1034,7 +1054,6 @@ void TrackingStudySpec::doMisalignmentStudy() std::array clArr{nullptr}; if (mParams->addPVAsCluster) { const auto& eve = mMCReader.getMCEventHeader(lbl.getSourceID(), lbl.getEventID()); - dataformats::VertexBase pv; auto trFitOut = iTrack.getParamIn(); pv.setXYZ(eve.GetX(), eve.GetY(), eve.GetZ()); if (!prop->propagateToDCA(pv, trFitOut, prop->getNominalBz(), base::Propagator::MAX_STEP, mParams->CorrType)) { @@ -1047,8 +1066,7 @@ void TrackingStudySpec::doMisalignmentStudy() o2::math_utils::sincos(trFitOut.getAlpha(), sinAlp, cosAlp); cl[0].alpha = trFitOut.getAlpha(); cl[0].setXYZ((pv.getX() * cosAlp) + (pv.getY() * sinAlp), (-pv.getX() * sinAlp) + (pv.getY() * cosAlp), pv.getZ()); - cl[0].setSigmaY2(0.5f * (pv.getSigmaX2() + pv.getSigmaY2())); - cl[0].setSigmaZ2(pv.getSigmaZ2()); + cl[0].setErrors(0.5f * (pv.getSigmaX2() + pv.getSigmaY2()), pv.getSigmaZ2(), 0.f); cl[0].setSensorID(-1); clArr[0] = &cl[0]; } @@ -1155,6 +1173,34 @@ void TrackingStudySpec::doMisalignmentStudy() LOGP(info, "\tdoMisalignmentStudy: refitted {} out of {} tracks ({} !sel, {} !fit, {} !fitMis)", goodRefit, itsTracks.size(), notPassedSel, fitFail, fitFailMis); } +// get IP for current PV +// copied from TrackITS::getImpactParams +void TrackingStudySpec::getImpactParams(const o2::track::TrackParCov& trk, const o2::dataformats::VertexBase& pv, float ip[2], float bz) +{ + float x = pv.getX(), y = pv.getY(), z = pv.getZ(); + float f1 = trk.getSnp(), r1 = std::sqrt((1.f - f1) * (1.f + f1)); + float xt = trk.getX(), yt = trk.getY(); + float sn = std::sin(trk.getAlpha()), cs = std::cos(trk.getAlpha()); + float a = (x * cs) + (y * sn); + y = (-x * sn) + (y * cs); + x = a; + xt -= x; + yt -= y; + float rp4 = trk.getCurvature(bz); + if ((std::abs(bz) < o2::constants::math::Almost0) || (std::abs(rp4) < o2::constants::math::Almost0)) { + ip[0] = -((xt * f1) - (yt * r1)); + ip[1] = trk.getZ() + ((ip[0] * f1 - xt) / r1 * trk.getTgl()) - z; + return; + } + sn = (rp4 * xt) - f1; + cs = (rp4 * yt) + r1; + a = (2 * (xt * f1 - yt * r1)) - (rp4 * (xt * xt + yt * yt)); + float rr = std::sqrt((sn * sn) + (cs * cs)); + ip[0] = -a / (1 + rr); + float f2 = -sn / rr, r2 = std::sqrt((1.f - f2) * (1.f + f2)); + ip[1] = trk.getZ() + (trk.getTgl() / rp4 * std::asin((f2 * r1) - (f1 * r2))) - z; +} + DataProcessorSpec getTrackingStudySpec(GTrackID::mask_t srcTracks, GTrackID::mask_t srcClusters, bool useMC, bool withPV) { std::vector outputs; @@ -1179,7 +1225,7 @@ DataProcessorSpec getTrackingStudySpec(GTrackID::mask_t srcTracks, GTrackID::mas .name = "its3-track-study", .inputs = dataRequest->inputs, .outputs = outputs, - .algorithm = AlgorithmSpec{adaptFromTask(dataRequest, ggRequest, srcTracks, useMC)}, + .algorithm = AlgorithmSpec{adaptFromTask(dataRequest, ggRequest, srcTracks, useMC, withPV)}, .options = {}}; } diff --git a/Detectors/Upgrades/ITS3/study/src/its3-tracking-study-workflow.cxx b/Detectors/Upgrades/ITS3/study/src/its3-tracking-study-workflow.cxx index 482ef2bb71e1d..c26ce34dea403 100644 --- a/Detectors/Upgrades/ITS3/study/src/its3-tracking-study-workflow.cxx +++ b/Detectors/Upgrades/ITS3/study/src/its3-tracking-study-workflow.cxx @@ -19,6 +19,7 @@ #include "Framework/CallbacksPolicy.h" #include "DetectorsBase/DPLWorkflowUtils.h" #include "GlobalTrackingWorkflowHelpers/InputHelper.h" +#include "DataFormatsITSMFT/DPLAlpideParamInitializer.h" #include "DetectorsRaw/HBFUtilsInitializer.h" using namespace o2::framework; @@ -43,6 +44,7 @@ void customize(std::vector& workflowOptions) {"disable-root-input", VariantType::Bool, false, {"disable root-files input reader"}}, {"configKeyValues", VariantType::String, "", {"Semicolon separated key=value strings ..."}}}; o2::raw::HBFUtilsInitializer::addConfigOption(options); + o2::itsmft::DPLAlpideParamInitializer::addITSConfigOption(options); std::swap(workflowOptions, options); } diff --git a/Detectors/Upgrades/ITS3/workflow/include/ITS3Workflow/ClustererSpec.h b/Detectors/Upgrades/ITS3/workflow/include/ITS3Workflow/ClustererSpec.h index e6732551c60e2..71efa54db94f0 100644 --- a/Detectors/Upgrades/ITS3/workflow/include/ITS3Workflow/ClustererSpec.h +++ b/Detectors/Upgrades/ITS3/workflow/include/ITS3Workflow/ClustererSpec.h @@ -28,7 +28,9 @@ namespace o2::its3 class ClustererDPL : public Task { public: - ClustererDPL(const std::shared_ptr& gr, bool useMC) : mGGCCDBRequest(gr), mUseMC(useMC) {} + static constexpr int NLayers = 7; + + ClustererDPL(const std::shared_ptr& gr, bool useMC, bool doStag) : mGGCCDBRequest(gr), mUseMC(useMC), mDoStaggering(doStag) {} void init(InitContext& ic) final; void run(ProcessingContext& pc) final; void finaliseCCDB(ConcreteDataMatcher& matcher, void* obj) final; @@ -39,15 +41,16 @@ class ClustererDPL : public Task std::shared_ptr mGGCCDBRequest; bool mUseMC = true; - int mState = 0; int mNThreads = 1; bool mUseClusterDictionary = true; + bool mDoStaggering = false; + std::vector mFilter; std::unique_ptr mClusterer = nullptr; }; /// create a processor spec /// run ITS cluster finder -framework::DataProcessorSpec getClustererSpec(bool useMC); +framework::DataProcessorSpec getClustererSpec(bool useMC, bool doStag); } // namespace o2::its3 diff --git a/Detectors/Upgrades/ITS3/workflow/include/ITS3Workflow/DigitReaderSpec.h b/Detectors/Upgrades/ITS3/workflow/include/ITS3Workflow/DigitReaderSpec.h index 9d0379aa77d78..e22f634c976cf 100644 --- a/Detectors/Upgrades/ITS3/workflow/include/ITS3Workflow/DigitReaderSpec.h +++ b/Detectors/Upgrades/ITS3/workflow/include/ITS3Workflow/DigitReaderSpec.h @@ -14,76 +14,77 @@ #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 doStag, 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) const + { + if (mDoStaggering) { + return base + "_" + std::to_string(index); + } + return base; + } - 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; bool mUseMC = true; // use MC truth bool mUseCalib = true; // send calib data + bool mDoStaggering = false; - 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"); +framework::DataProcessorSpec getITS3DigitReaderSpec(bool useMC = true, bool doStag = false, 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/include/ITS3Workflow/DigitWriterSpec.h b/Detectors/Upgrades/ITS3/workflow/include/ITS3Workflow/DigitWriterSpec.h index 309c515d92b93..97fb9481a236c 100644 --- a/Detectors/Upgrades/ITS3/workflow/include/ITS3Workflow/DigitWriterSpec.h +++ b/Detectors/Upgrades/ITS3/workflow/include/ITS3Workflow/DigitWriterSpec.h @@ -14,13 +14,9 @@ #include "Framework/DataProcessorSpec.h" -namespace o2 +namespace o2::its3 { -namespace its3 -{ - -o2::framework::DataProcessorSpec getITS3DigitWriterSpec(bool mctruth = true, bool dec = false, bool calib = false); -} // namespace its3 -} // end namespace o2 +o2::framework::DataProcessorSpec getITS3DigitWriterSpec(bool mctruth = true, bool doStag = false, bool dec = false, bool calib = false); +} // namespace o2::its3 #endif /* STEER_ITSMFTDIGITWRITER_H_ */ diff --git a/Detectors/Upgrades/ITS3/workflow/include/ITS3Workflow/RecoWorkflow.h b/Detectors/Upgrades/ITS3/workflow/include/ITS3Workflow/RecoWorkflow.h index 010e1cd0a8127..164eb778b5093 100644 --- a/Detectors/Upgrades/ITS3/workflow/include/ITS3Workflow/RecoWorkflow.h +++ b/Detectors/Upgrades/ITS3/workflow/include/ITS3Workflow/RecoWorkflow.h @@ -22,6 +22,7 @@ namespace o2::its3::reco_workflow { framework::WorkflowSpec getWorkflow(bool useMC, + bool doStag, its::TrackingMode::Type trmode, o2::gpu::gpudatatypes::DeviceType dtype, bool useGPUWorkflow, diff --git a/Detectors/Upgrades/ITS3/workflow/include/ITS3Workflow/TrackerSpec.h b/Detectors/Upgrades/ITS3/workflow/include/ITS3Workflow/TrackerSpec.h index 66d58a5f5a92c..4b8c286f6c9fb 100644 --- a/Detectors/Upgrades/ITS3/workflow/include/ITS3Workflow/TrackerSpec.h +++ b/Detectors/Upgrades/ITS3/workflow/include/ITS3Workflow/TrackerSpec.h @@ -42,6 +42,7 @@ class TrackerDPL : public framework::Task public: TrackerDPL(std::shared_ptr gr, bool isMC, + bool doStag, int trgType, its::TrackingMode::Type trmode = its::TrackingMode::Unset, const bool overrBeamEst = false, @@ -69,7 +70,7 @@ class TrackerDPL : public framework::Task /// create a processor spec /// run ITS CA tracker -framework::DataProcessorSpec getTrackerSpec(bool useMC, bool useGeom, int useTrig, its::TrackingMode::Type trMode, const bool overrBeamEst = false, gpu::gpudatatypes::DeviceType dType = gpu::gpudatatypes::DeviceType::CPU); +framework::DataProcessorSpec getTrackerSpec(bool useMC, bool doStag, bool useGeom, int useTrig, its::TrackingMode::Type trMode, const bool overrBeamEst = false, gpu::gpudatatypes::DeviceType dType = gpu::gpudatatypes::DeviceType::CPU); } // namespace o2::its3 diff --git a/Detectors/Upgrades/ITS3/workflow/src/ClustererSpec.cxx b/Detectors/Upgrades/ITS3/workflow/src/ClustererSpec.cxx index 73b5f4650d02d..94803466d70b4 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,149 @@ 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 < (mDoStaggering ? NLayers : 1); ++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 < (mDoStaggering ? NLayers : 1); ++iLayer) { + int layer = (mDoStaggering) ? iLayer : -1; + sw.Start(); + LOG(info) << "Clusterer" << ((mDoStaggering) ? std::format(" on layer {}", layer) : "") << " pulled " << digits[iLayer].size() << " digits, in " << rofs[iLayer].size() << " RO frames"; + mClusterer->setMaxROFDepthToSquash(mClusterer->getMaxROFDepthToSquash(layer)); + o2::dataformats::ConstMCTruthContainerView labels(labelsbuffer[iLayer]); + reader.setSquashingDepth(mClusterer->getMaxROFDepthToSquash(layer)); + reader.setSquashingDist(mClusterer->getMaxRowColDiffToMask()); // Sharing same parameter/logic with masking + reader.setMaxBCSeparationToSquash(mClusterer->getMaxBCSeparationToSquash(layer)); + reader.setDigits(digits[iLayer]); + reader.setROFRecords(rofs[iLayer]); + if (mUseMC) { + LOG(info) << "Clusterer" << ((mDoStaggering) ? std::format(" on layer {}", layer) : "") << " 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 {}{}", ir.asString(), firstTForbit, ((mDoStaggering) ? std::format(" on layer {}", layer) : "")); + 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{}", ir.asString(), firstTForbit, ((mDoStaggering) ? std::format(" on layer {}", 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{}", ir.asString(), ((mDoStaggering) ? std::format(" on layer {}", layer) : "")); + 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{}", rof.asString(), rof.getNEntries(), expROF.getNEntries(), ((mDoStaggering) ? std::format(" on layer {}", layer) : "")); + expROF.setFirstEntry(rof.getFirstEntry()); + expROF.setNEntries(rof.getNEntries()); + } else { + LOGP(warn, "Repeating {} with {} clusters, discard preferring already processed instance with {} clusters{}", rof.asString(), rof.getNEntries(), expROF.getNEntries(), ((mDoStaggering) ? std::format(" on layer {}", layer) : "")); + } + } + } + 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 + // write dummy MC2ROF vector to keep writer/readers backward compatible + static std::vector dummyMC2ROF; + pc.outputs().snapshot(Output{"ITS", "CLUSTERSMC2ROF", iLayer}, dummyMC2ROF); + } + 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 +197,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 +216,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 } @@ -165,43 +258,43 @@ void ClustererDPL::finaliseCCDB(ConcreteDataMatcher& matcher, void* obj) void ClustererDPL::endOfStream(o2::framework::EndOfStreamContext& /*ec*/) { - mClusterer->print(); + mClusterer->print(true); } -DataProcessorSpec getClustererSpec(bool useMC) +DataProcessorSpec getClustererSpec(bool useMC, bool doStag) { 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 < (doStag ? ClustererDPL::NLayers : 1); ++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); + outputs.emplace_back("ITS", "CLUSTERSMC2ROF", 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, doStag)}, + .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..141457c319b9b 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 doStag, bool useCalib) : mUseMC(useMC), mDoStaggering(doStag), 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 < (mDoStaggering ? NLayers : 1); ++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,52 +74,60 @@ 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 < (mDoStaggering ? NLayers : 1); ++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 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"; } -DataProcessorSpec getITS3DigitReaderSpec(bool useMC, bool useCalib, std::string defname) +DataProcessorSpec getITS3DigitReaderSpec(bool useMC, bool doStag, 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 < (doStag ? ITS3DigitReader::NLayers : 1); ++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{ - "its3-digit-reader", - Inputs{}, - outputSpec, - AlgorithmSpec{adaptFromTask(useMC, useCalib)}, - Options{ + .name = "its3-digit-reader", + .inputs = Inputs{}, + .outputs = outputSpec, + .algorithm = AlgorithmSpec{adaptFromTask(useMC, doStag, useCalib)}, + .options = Options{ {"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..dcc28b5745f5c 100644 --- a/Detectors/Upgrades/ITS3/workflow/src/DigitWriterSpec.cxx +++ b/Detectors/Upgrades/ITS3/workflow/src/DigitWriterSpec.cxx @@ -14,41 +14,44 @@ #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 using namespace o2::framework; using SubSpecificationType = o2::framework::DataAllocator::SubSpecificationType; -namespace o2 -{ -namespace its3 +namespace o2::its3 { template using BranchDefinition = MakeRootTreeWriterSpec::BranchDefinition; using MCCont = o2::dataformats::ConstMCTruthContainer; -/// create the processor spec -/// describing a processor receiving digits for ITS/MFT and writing them to file -DataProcessorSpec getDigitWriterSpec(bool mctruth, bool dec, bool calib, o2::header::DataOrigin detOrig, o2::detectors::DetID detId) +DataProcessorSpec getITS3DigitWriterSpec(bool mctruth, bool doStag, bool dec, bool calib) { - std::string detStr = o2::detectors::DetID::getName(detId); + std::string detStr = o2::detectors::DetID::getName(o2::detectors::DetID::IT3); std::string detStrL = dec ? "o2_" : ""; // for decoded digits prepend by o2 + const int mLayers = doStag ? 7 : 1; 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>(mLayers, 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>(mLayers, 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 +71,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,31 +85,51 @@ 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 = [&doStag](std::string base, size_t index) -> std::string { + if (doStag) { + return base += "_" + std::to_string(index); + } + return base; + }; + + std::vector vecInpSpecDig, vecInpSpecROF, vecInpSpecLbl; + vecInpSpecDig.reserve(mLayers); + vecInpSpecROF.reserve(mLayers); + vecInpSpecLbl.reserve(mLayers); + for (int iLayer = 0; iLayer < mLayers; 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::TreeAttributes{.name = "o2sim", .title = "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()})(); -} - -DataProcessorSpec getITS3DigitWriterSpec(bool mctruth, bool dec, bool calib) -{ - return getDigitWriterSpec(mctruth, dec, calib, o2::header::gDataOriginIT3, o2::detectors::DetID::IT3); + BranchDefinition>{vecInpSpecDig, + detStr + "Digit", "digit-branch", + mLayers, + digitSizeGetter, + getIndex, + getName}, + BranchDefinition>{vecInpSpecROF, + detStr + "DigitROF", "digit-rof-branch", + mLayers, + rofSizeGetter, + getIndex, + getName}, + BranchDefinition>{vecInpSpecLbl, + "IT3DigitMCTruth", "digit-mctruth-branch", + (mctruth ? mLayers : 0), + fillLabels, + getIndex, + getName})(); } -} // end namespace its3 -} // end namespace o2 +} // namespace o2::its3 + // end namespace o2 diff --git a/Detectors/Upgrades/ITS3/workflow/src/RecoWorkflow.cxx b/Detectors/Upgrades/ITS3/workflow/src/RecoWorkflow.cxx index f27fda19fe00c..87643772a9b2f 100644 --- a/Detectors/Upgrades/ITS3/workflow/src/RecoWorkflow.cxx +++ b/Detectors/Upgrades/ITS3/workflow/src/RecoWorkflow.cxx @@ -26,27 +26,28 @@ static std::shared_ptr gTask; namespace o2::its3::reco_workflow { -framework::WorkflowSpec getWorkflow(bool useMC, its::TrackingMode::Type trmode, o2::gpu::gpudatatypes::DeviceType dtype, bool useGPUWorkflow, +framework::WorkflowSpec getWorkflow(bool useMC, bool doStag, its::TrackingMode::Type trmode, o2::gpu::gpudatatypes::DeviceType dtype, bool useGPUWorkflow, bool upstreamDigits, bool upstreamClusters, bool disableRootOutput, bool useGeom, int useTrig, bool overrideBeamPosition) { framework::WorkflowSpec specs; if (!(upstreamDigits || upstreamClusters)) { - specs.emplace_back(o2::its3::getITS3DigitReaderSpec(useMC, false, "it3digits.root")); + specs.emplace_back(o2::its3::getITS3DigitReaderSpec(useMC, doStag, false, "it3digits.root")); } if (!upstreamClusters) { - specs.emplace_back(o2::its3::getClustererSpec(useMC)); + specs.emplace_back(o2::its3::getClustererSpec(useMC, doStag)); } if (!disableRootOutput) { - specs.emplace_back(o2::itsmft::getITSClusterWriterSpec(useMC, false)); + specs.emplace_back(o2::itsmft::getITSClusterWriterSpec(useMC, doStag, false)); } if (trmode != its::TrackingMode::Off) { if (useGPUWorkflow) { o2::gpu::GPURecoWorkflowSpec::Config cfg; cfg.runITSTracking = true; + cfg.itsStaggered = doStag, cfg.isITS3 = true; cfg.itsTriggerType = useTrig; cfg.itsOverrBeamEst = overrideBeamPosition; @@ -66,13 +67,13 @@ 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)); + specs.emplace_back(o2::its3::getTrackerSpec(useMC, doStag, useGeom, useTrig, trmode, overrideBeamPosition, dtype)); } if (!disableRootOutput) { specs.emplace_back(o2::its::getTrackWriterSpec(useMC)); diff --git a/Detectors/Upgrades/ITS3/workflow/src/TrackerSpec.cxx b/Detectors/Upgrades/ITS3/workflow/src/TrackerSpec.cxx index 8db02d7227e7f..94e711a05a2d6 100644 --- a/Detectors/Upgrades/ITS3/workflow/src/TrackerSpec.cxx +++ b/Detectors/Upgrades/ITS3/workflow/src/TrackerSpec.cxx @@ -41,12 +41,13 @@ using Vertex = o2::dataformats::Vertex>; TrackerDPL::TrackerDPL(std::shared_ptr gr, bool isMC, + bool doStag, int trgType, its::TrackingMode::Type trMode, const bool overrBeamEst, o2::gpu::gpudatatypes::DeviceType dType) : mGGCCDBRequest(gr), mRecChain{o2::gpu::GPUReconstruction::CreateInstance(dType, true)}, - mITS3TrackingInterface{isMC, false, trgType, overrBeamEst} + mITS3TrackingInterface{isMC, doStag, trgType, overrBeamEst} { mITS3TrackingInterface.setTrackingMode(trMode); } @@ -88,19 +89,25 @@ void TrackerDPL::endOfStream(EndOfStreamContext& ec) LOGF(info, "ITS3 CA-Tracker total timing: Cpu: %.3e Real: %.3e s in %d slots", mTimer.CpuTime(), mTimer.RealTime(), mTimer.Counter() - 1); } -DataProcessorSpec getTrackerSpec(bool useMC, bool useGeom, int trgType, its::TrackingMode::Type trMode, const bool overrBeamEst, o2::gpu::gpudatatypes::DeviceType dType) +DataProcessorSpec getTrackerSpec(bool useMC, bool doStag, 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 < (doStag ? 7 : 1); ++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 +124,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,20 +132,17 @@ 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{ - "its3-tracker", - inputs, - outputs, - AlgorithmSpec{adaptFromTask(ggRequest, useMC, trgType, trMode, overrBeamEst, dType)}, - Options{}}; + .name = "its3-tracker", + .inputs = inputs, + .outputs = outputs, + .algorithm = AlgorithmSpec{adaptFromTask(ggRequest, useMC, doStag, trgType, trMode, overrBeamEst, dType)}, + .options = Options{}}; } } // namespace its3 diff --git a/Detectors/Upgrades/ITS3/workflow/src/digit-writer-workflow.cxx b/Detectors/Upgrades/ITS3/workflow/src/digit-writer-workflow.cxx index f37fce71dc2d2..0f36c2c826f9a 100644 --- a/Detectors/Upgrades/ITS3/workflow/src/digit-writer-workflow.cxx +++ b/Detectors/Upgrades/ITS3/workflow/src/digit-writer-workflow.cxx @@ -13,6 +13,7 @@ #include "CommonUtils/ConfigurableParam.h" #include "Framework/ConfigParamSpec.h" #include "Framework/CompletionPolicyHelpers.h" +#include "DataFormatsITSMFT/DPLAlpideParamInitializer.h" using namespace o2::framework; @@ -31,7 +32,7 @@ void customize(std::vector& workflowOptions) ConfigParamSpec{"disable-mc", VariantType::Bool, false, {"disable mc truth"}}, ConfigParamSpec{"enable-calib-data", VariantType::Bool, false, {"enable writing GBT calibration data"}}, ConfigParamSpec{"configKeyValues", VariantType::String, "", {"semicolon separated key=value strings"}}}; - + o2::itsmft::DPLAlpideParamInitializer::addITSConfigOption(options); std::swap(workflowOptions, options); } @@ -46,7 +47,7 @@ WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) bool calib = cfgc.options().get("enable-calib-data"); // Update the (declared) parameters if changed from the command line o2::conf::ConfigurableParam::updateFromString(cfgc.options().get("configKeyValues")); - - wf.emplace_back(o2::its3::getITS3DigitWriterSpec(useMC, true, calib)); + bool doStag = o2::itsmft::DPLAlpideParamInitializer::isITSStaggeringEnabled(cfgc); + wf.emplace_back(o2::its3::getITS3DigitWriterSpec(useMC, doStag, true, calib)); return wf; } diff --git a/Detectors/Upgrades/ITS3/workflow/src/its3-reco-workflow.cxx b/Detectors/Upgrades/ITS3/workflow/src/its3-reco-workflow.cxx index 018ad807a974a..beb6815bc9bcd 100644 --- a/Detectors/Upgrades/ITS3/workflow/src/its3-reco-workflow.cxx +++ b/Detectors/Upgrades/ITS3/workflow/src/its3-reco-workflow.cxx @@ -16,6 +16,7 @@ #include "DetectorsRaw/HBFUtilsInitializer.h" #include "Framework/CallbacksPolicy.h" #include "Framework/CompletionPolicyHelpers.h" +#include "DataFormatsITSMFT/DPLAlpideParamInitializer.h" #include "GPUO2Interface.h" #include "GPUReconstruction.h" @@ -51,6 +52,7 @@ void customize(std::vector& workflowOptions) {"use-gpu-workflow", o2::framework::VariantType::Bool, false, {"use GPU workflow (default: false)"}}, {"gpu-device", o2::framework::VariantType::Int, 1, {"use gpu device: CPU=1,CUDA=2,HIP=3 (default: CPU)"}}}; o2::raw::HBFUtilsInitializer::addConfigOption(options); + o2::itsmft::DPLAlpideParamInitializer::addITSConfigOption(options); std::swap(workflowOptions, options); } @@ -69,6 +71,7 @@ WorkflowSpec defineDataProcessing(ConfigContext const& configcontext) auto disableRootOutput = configcontext.options().get("disable-root-output"); auto useGeom = configcontext.options().get("use-full-geometry"); auto useGPUWfx = configcontext.options().get("use-gpu-workflow"); + bool doStag = o2::itsmft::DPLAlpideParamInitializer::isITSStaggeringEnabled(configcontext); o2::conf::ConfigurableParam::updateFromString(configcontext.options().get("configKeyValues")); int trType = 0; @@ -81,7 +84,7 @@ WorkflowSpec defineDataProcessing(ConfigContext const& configcontext) LOG(fatal) << "Unknown trigger type requested for events prescaling: " << selTrig; } } - auto wf = o2::its3::reco_workflow::getWorkflow(useMC, o2::its::TrackingMode::fromString(trmode), gpuDevice, useGPUWfx, extDigits, extClusters, disableRootOutput, useGeom, trType, beamPosOVerride); + auto wf = o2::its3::reco_workflow::getWorkflow(useMC, doStag, o2::its::TrackingMode::fromString(trmode), gpuDevice, useGPUWfx, extDigits, extClusters, disableRootOutput, useGeom, trType, beamPosOVerride); // configure dpl timer to inject correct firstTForbit: start from the 1st orbit of TF containing 1st sampled orbit o2::raw::HBFUtilsInitializer hbfIni(configcontext, wf); diff --git a/Steer/DigitizerWorkflow/src/ITS3DigitizerSpec.cxx b/Steer/DigitizerWorkflow/src/ITS3DigitizerSpec.cxx index 19552a407ec57..760c1e14ef400 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; @@ -64,7 +49,7 @@ class ITS3DPLDigitizerTask : BaseDPLDigitizer public: using BaseDPLDigitizer::init; - ITS3DPLDigitizerTask(bool mctruth = true) : BaseDPLDigitizer(InitServices::FIELD | InitServices::GEOM), mWithMCTruth(mctruth) {} + ITS3DPLDigitizerTask(bool mctruth = true, bool doStag = false) : BaseDPLDigitizer(InitServices::FIELD | InitServices::GEOM), mWithMCTruth(mctruth), mDoStaggering(doStag) {} void initDigitizerTask(framework::InitContext& ic) override { @@ -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,151 @@ 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 < (mDoStaggering ? 7 : 1); ++iLayer) { + const int layer = (mDoStaggering) ? iLayer : -1; + 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" << ((mDoStaggering) ? std::format(" 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(layer); // 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, layer); + 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, layer); // 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, layer); + 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 {}{}", ir.asString(), mFirstOrbitTF, ((mDoStaggering) ? std::format(" on layer {}", 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{}", ir.asString(), mFirstOrbitTF, ((mDoStaggering) ? std::format(" on layer {}", 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 +275,19 @@ class ITS3DPLDigitizerTask : BaseDPLDigitizer digipar.setIBNSimSteps(doptIB.nIBSimSteps); digipar.setIBNoisePerPixel(doptIB.IBNoisePerPixel); + // staggering parameters + if (mDoStaggering) { + 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 +335,41 @@ 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 + bool mDoStaggering = false; }; -DataProcessorSpec getITS3DigitizerSpec(int channel, bool mctruth) +namespace +{ +std::vector makeOutChannels(o2::header::DataOrigin detOrig, bool mctruth, bool doStag) +{ + std::vector outputs; + for (int iLayer{0}; iLayer < (doStag ? 7 : 1); ++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, bool doStag) { - 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 +380,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, doStag), + .algorithm = AlgorithmSpec{adaptFromTask(mctruth)}, + .options = Options{{"disable-qed", o2::framework::VariantType::Bool, false, {"disable QED handling"}}}}; } } // namespace o2::its3 diff --git a/Steer/DigitizerWorkflow/src/ITS3DigitizerSpec.h b/Steer/DigitizerWorkflow/src/ITS3DigitizerSpec.h index 33fe2383ec21f..3dd566e576673 100644 --- a/Steer/DigitizerWorkflow/src/ITS3DigitizerSpec.h +++ b/Steer/DigitizerWorkflow/src/ITS3DigitizerSpec.h @@ -17,7 +17,7 @@ namespace o2::its3 { -o2::framework::DataProcessorSpec getITS3DigitizerSpec(int channel, bool mctruth = true); +o2::framework::DataProcessorSpec getITS3DigitizerSpec(int channel, bool mctruth = true, bool doStag = false); } // namespace o2::its3 // end namespace o2 diff --git a/Steer/DigitizerWorkflow/src/SimpleDigitizerWorkflow.cxx b/Steer/DigitizerWorkflow/src/SimpleDigitizerWorkflow.cxx index b4f9c1643d150..3b7bf3088a8f9 100644 --- a/Steer/DigitizerWorkflow/src/SimpleDigitizerWorkflow.cxx +++ b/Steer/DigitizerWorkflow/src/SimpleDigitizerWorkflow.cxx @@ -639,7 +639,7 @@ WorkflowSpec defineDataProcessing(ConfigContext const& configcontext) // the ITS part if (isEnabled(o2::detectors::DetID::ITS)) { detList.emplace_back(o2::detectors::DetID::ITS); - bool doStag = o2::itsmft::DPLAlpideParamInitializer::isMFTStaggeringEnabled(configcontext); + bool doStag = o2::itsmft::DPLAlpideParamInitializer::isITSStaggeringEnabled(configcontext); // connect the ITS digitization digitizerSpecs.emplace_back(o2::itsmft::getITSDigitizerSpec(fanoutsize++, mctruth, doStag)); // connect ITS digit writer @@ -650,10 +650,11 @@ WorkflowSpec defineDataProcessing(ConfigContext const& configcontext) // the ITS3 part if (isEnabled(o2::detectors::DetID::IT3)) { detList.emplace_back(o2::detectors::DetID::IT3); + bool doStag = o2::itsmft::DPLAlpideParamInitializer::isITSStaggeringEnabled(configcontext); // connect the ITS digitization - specs.emplace_back(o2::its3::getITS3DigitizerSpec(fanoutsize++, mctruth)); + specs.emplace_back(o2::its3::getITS3DigitizerSpec(fanoutsize++, mctruth, doStag)); // // connect ITS digit writer - specs.emplace_back(o2::its3::getITS3DigitWriterSpec(mctruth)); + specs.emplace_back(o2::its3::getITS3DigitWriterSpec(mctruth, doStag)); } // the ALICE 3 TRK part