From ad4352bf2ee982f83a2128f67adbc950ab05ffe5 Mon Sep 17 00:00:00 2001 From: Felix Schlepper Date: Thu, 19 Jun 2025 11:38:46 +0200 Subject: [PATCH 1/8] ITS: for deltaROF consider neighbouring vertices Signed-off-by: Felix Schlepper --- .../ITS/tracking/include/ITStracking/TimeFrame.h | 3 ++- Detectors/ITSMFT/ITS/tracking/src/TrackerTraits.cxx | 12 +++++++++--- 2 files changed, 11 insertions(+), 4 deletions(-) diff --git a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/TimeFrame.h b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/TimeFrame.h index 9434fc2292750..22f67c6783bea 100644 --- a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/TimeFrame.h +++ b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/TimeFrame.h @@ -397,7 +397,8 @@ inline gsl::span TimeFrame::getPrimaryVertices(int romin, if (mPrimaryVertices.empty()) { return {}; } - return {&mPrimaryVertices[mROFramesPV[romin]], static_cast::size_type>(mROFramesPV[romax + 1] - mROFramesPV[romin])}; + const int stop_idx = romax >= mNrof - 1 ? mNrof : romax + 1; + return {&mPrimaryVertices[mROFramesPV[romin]], static_cast::size_type>(mROFramesPV[stop_idx] - mROFramesPV[romin])}; } template diff --git a/Detectors/ITSMFT/ITS/tracking/src/TrackerTraits.cxx b/Detectors/ITSMFT/ITS/tracking/src/TrackerTraits.cxx index da7f31bd678b5..014031aaedccb 100644 --- a/Detectors/ITSMFT/ITS/tracking/src/TrackerTraits.cxx +++ b/Detectors/ITSMFT/ITS/tracking/src/TrackerTraits.cxx @@ -74,11 +74,14 @@ void TrackerTraits::computeLayerTracklets(const int iteration, int iROF mTaskArena->execute([&] { for (int rof0{startROF}; rof0 < endROF; ++rof0) { - gsl::span primaryVertices = mTrkParams[iteration].UseDiamond ? diamondSpan : mTimeFrame->getPrimaryVertices(rof0); - const int startVtx{iVertex >= 0 ? iVertex : 0}; - const int endVtx{iVertex >= 0 ? o2::gpu::CAMath::Min(iVertex + 1, static_cast(primaryVertices.size())) : static_cast(primaryVertices.size())}; int minRof = o2::gpu::CAMath::Max(startROF, rof0 - mTrkParams[iteration].DeltaROF); int maxRof = o2::gpu::CAMath::Min(endROF - 1, rof0 + mTrkParams[iteration].DeltaROF); + gsl::span primaryVertices = mTrkParams[iteration].UseDiamond ? diamondSpan : mTimeFrame->getPrimaryVertices(minRof, maxRof); + if (primaryVertices.empty()) { + continue; + } + const int startVtx{iVertex >= 0 ? iVertex : 0}; + const int endVtx{iVertex >= 0 ? o2::gpu::CAMath::Min(iVertex + 1, static_cast(primaryVertices.size())) : static_cast(primaryVertices.size())}; tbb::parallel_for( tbb::blocked_range(0, mTrkParams[iteration].TrackletsPerRoad()), @@ -127,6 +130,9 @@ void TrackerTraits::computeLayerTracklets(const int iteration, int iROF } for (int rof1{minRof}; rof1 <= maxRof; ++rof1) { + if (!mTimeFrame->mMultiplicityCutMask[rof1]) { + continue; + } auto layer1 = mTimeFrame->getClustersOnLayer(rof1, iLayer + 1); if (layer1.empty()) { continue; From ccd8f0c27b74923774483e203b460bc5da731982 Mon Sep 17 00:00:00 2001 From: Felix Schlepper Date: Thu, 26 Jun 2025 17:42:17 +0200 Subject: [PATCH 2/8] ITS: remove copying vertices addPrimaryVertices Signed-off-by: Felix Schlepper --- .../tracking/include/ITStracking/TimeFrame.h | 4 +- .../ITSMFT/ITS/tracking/src/TimeFrame.cxx | 43 +++---------------- .../ITS/tracking/src/TrackingInterface.cxx | 4 +- .../ITS/tracking/src/VertexerTraits.cxx | 4 +- 4 files changed, 10 insertions(+), 45 deletions(-) diff --git a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/TimeFrame.h b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/TimeFrame.h index 22f67c6783bea..bc129020e0710 100644 --- a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/TimeFrame.h +++ b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/TimeFrame.h @@ -77,11 +77,9 @@ struct TimeFrame { gsl::span> getPrimaryVerticesXAlpha(int rofId) const; void fillPrimaryVerticesXandAlpha(); int getPrimaryVerticesNum(int rofId = -1) const; - void addPrimaryVertices(const bounded_vector& vertices); void addPrimaryVerticesLabels(bounded_vector>& labels); void addPrimaryVerticesContributorLabels(bounded_vector& labels); - void addPrimaryVertices(const bounded_vector& vertices, const int rofId, const int iteration); - void addPrimaryVertices(const gsl::span& vertices, const int rofId, const int iteration); + void addPrimaryVertices(const bounded_vector& vertices, const int iteration); void addPrimaryVerticesInROF(const bounded_vector& vertices, const int rofId, const int iteration); void addPrimaryVerticesLabelsInROF(const bounded_vector>& labels, const int rofId); void addPrimaryVerticesContributorLabelsInROF(const bounded_vector& labels, const int rofId); diff --git a/Detectors/ITSMFT/ITS/tracking/src/TimeFrame.cxx b/Detectors/ITSMFT/ITS/tracking/src/TimeFrame.cxx index 4115726756e73..694bfb244e9aa 100644 --- a/Detectors/ITSMFT/ITS/tracking/src/TimeFrame.cxx +++ b/Detectors/ITSMFT/ITS/tracking/src/TimeFrame.cxx @@ -59,24 +59,19 @@ TimeFrame::~TimeFrame() } template -void TimeFrame::addPrimaryVertices(const bounded_vector& vertices) +void TimeFrame::addPrimaryVertices(const bounded_vector& vertices, const int iteration) { for (const auto& vertex : vertices) { - mPrimaryVertices.emplace_back(vertex); - if (!isBeamPositionOverridden) { + mPrimaryVertices.emplace_back(vertex); // put a copy in the present + mTotVertPerIteration[iteration]++; + if (!isBeamPositionOverridden) { // beam position is updated only at first occurrence of the vertex. A bit sketchy if we have past/future vertices, it should not impact too much. const float w = vertex.getNContributors(); mBeamPos[0] = (mBeamPos[0] * mBeamPosWeight + vertex.getX() * w) / (mBeamPosWeight + w); mBeamPos[1] = (mBeamPos[1] * mBeamPosWeight + vertex.getY() * w) / (mBeamPosWeight + w); mBeamPosWeight += w; } } - mROFramesPV.push_back(mPrimaryVertices.size()); -} - -template -void TimeFrame::addPrimaryVertices(const bounded_vector& vertices, const int rofId, const int iteration) -{ - addPrimaryVertices(gsl::span(vertices), rofId, iteration); + mROFramesPV.push_back(mPrimaryVertices.size()); // current rof must have number of vertices up to present } template @@ -119,34 +114,6 @@ void TimeFrame::addPrimaryVerticesContributorLabelsInROF(const bounded_ mVerticesContributorLabels.insert(mVerticesContributorLabels.begin() + n, labels.begin(), labels.end()); } -template -void TimeFrame::addPrimaryVertices(const gsl::span& vertices, const int rofId, const int iteration) -{ - bounded_vector futureVertices(mMemoryPool.get()); - for (const auto& vertex : vertices) { - if (vertex.getTimeStamp().getTimeStamp() < rofId) { // put a copy in the past - insertPastVertex(vertex, iteration); - } else { - if (vertex.getTimeStamp().getTimeStamp() > rofId) { // or put a copy in the future - futureVertices.emplace_back(vertex); - } - } - mPrimaryVertices.emplace_back(vertex); // put a copy in the present - mTotVertPerIteration[iteration]++; - if (!isBeamPositionOverridden) { // beam position is updated only at first occurrence of the vertex. A bit sketchy if we have past/future vertices, it should not impact too much. - const float w = vertex.getNContributors(); - mBeamPos[0] = (mBeamPos[0] * mBeamPosWeight + vertex.getX() * w) / (mBeamPosWeight + w); - mBeamPos[1] = (mBeamPos[1] * mBeamPosWeight + vertex.getY() * w) / (mBeamPosWeight + w); - mBeamPosWeight += w; - } - } - mROFramesPV.push_back(mPrimaryVertices.size()); // current rof must have number of vertices up to present - for (auto& vertex : futureVertices) { - mPrimaryVertices.emplace_back(vertex); - mTotVertPerIteration[iteration]++; - } -} - template int TimeFrame::loadROFrameData(gsl::span rofs, gsl::span clusters, diff --git a/Detectors/ITSMFT/ITS/tracking/src/TrackingInterface.cxx b/Detectors/ITSMFT/ITS/tracking/src/TrackingInterface.cxx index a908f8b2a1f1e..00a69a37cb51a 100644 --- a/Detectors/ITSMFT/ITS/tracking/src/TrackingInterface.cxx +++ b/Detectors/ITSMFT/ITS/tracking/src/TrackingInterface.cxx @@ -165,7 +165,7 @@ void ITSTrackingInterface::run(framework::ProcessingContext& pc) gsl::span> vMCRecInfo; gsl::span vMCContLabels; for (auto iRof{0}; iRof < trackROFspan.size(); ++iRof) { - std::vector vtxVecLoc; + bounded_vector vtxVecLoc; auto& vtxROF = vertROFvec.emplace_back(trackROFspan[iRof]); vtxROF.setFirstEntry(vertices.size()); if (mRunVertexer) { @@ -223,7 +223,7 @@ void ITSTrackingInterface::run(framework::ProcessingContext& pc) for (auto& v : vtxVecLoc) { vertices.push_back(v); } - mTimeFrame->addPrimaryVertices(vtxVecLoc, iRof, 0); + mTimeFrame->addPrimaryVertices(vtxVecLoc, 0); } } if (mRunVertexer) { diff --git a/Detectors/ITSMFT/ITS/tracking/src/VertexerTraits.cxx b/Detectors/ITSMFT/ITS/tracking/src/VertexerTraits.cxx index 1069f1808fb2a..fdfc951b5cece 100644 --- a/Detectors/ITSMFT/ITS/tracking/src/VertexerTraits.cxx +++ b/Detectors/ITSMFT/ITS/tracking/src/VertexerTraits.cxx @@ -520,7 +520,7 @@ void VertexerTraits::computeVertices(const int iteration) } } if (!iteration) { - mTimeFrame->addPrimaryVertices(vertices, rofId, iteration); + mTimeFrame->addPrimaryVertices(vertices, iteration); if (mTimeFrame->hasMCinformation()) { mTimeFrame->addPrimaryVerticesLabels(polls); if (mVrtParams[iteration].outputContLabels) { @@ -625,7 +625,7 @@ void VertexerTraits::addTruthSeedingVertices() } else { mTimeFrame->getNoVertexROF()++; } - mTimeFrame->addPrimaryVertices(verts, iROF, 0); + mTimeFrame->addPrimaryVertices(verts, 0); mTimeFrame->addPrimaryVerticesLabels(polls); } LOGP(info, "Found {}/{} ROFs with {} vertices -> ={:.2f}", vertices.size(), mTimeFrame->getNrof(), nVerts, (float)nVerts / (float)vertices.size()); From ed109e2dbb0d4fabba149c1a0281f11b65dd66a9 Mon Sep 17 00:00:00 2001 From: Felix Schlepper Date: Tue, 1 Jul 2025 17:13:22 +0200 Subject: [PATCH 3/8] ITS: add extra protection if no tracklets were produced in ROF Signed-off-by: Felix Schlepper --- Detectors/ITSMFT/ITS/tracking/include/ITStracking/TimeFrame.h | 2 +- Detectors/ITSMFT/ITS/tracking/src/VertexerTraits.cxx | 3 +++ 2 files changed, 4 insertions(+), 1 deletion(-) diff --git a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/TimeFrame.h b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/TimeFrame.h index bc129020e0710..3f0d291d5e51d 100644 --- a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/TimeFrame.h +++ b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/TimeFrame.h @@ -593,7 +593,7 @@ inline gsl::span TimeFrame::getExclusiveNTrackletsCluster(int rofI template inline gsl::span TimeFrame::getFoundTracklets(int rofId, int combId) { - if (rofId < 0 || rofId >= mNrof) { + if (rofId < 0 || rofId >= mNrof || mTracklets[combId].empty()) { return {}; } auto startIdx{mNTrackletsPerROF[combId][rofId]}; diff --git a/Detectors/ITSMFT/ITS/tracking/src/VertexerTraits.cxx b/Detectors/ITSMFT/ITS/tracking/src/VertexerTraits.cxx index fdfc951b5cece..8ff4dffef5cdc 100644 --- a/Detectors/ITSMFT/ITS/tracking/src/VertexerTraits.cxx +++ b/Detectors/ITSMFT/ITS/tracking/src/VertexerTraits.cxx @@ -326,6 +326,9 @@ void VertexerTraits::computeTrackletMatching(const int iteration) if (iteration && (int)mTimeFrame->getPrimaryVertices(pivotRofId).size() > mVrtParams[iteration].vertPerRofThreshold) { continue; } + if (mTimeFrame->getFoundTracklets(pivotRofId, 0).empty()) { + continue; + } mTimeFrame->getLines(pivotRofId).reserve(mTimeFrame->getNTrackletsCluster(pivotRofId, 0).size()); bounded_vector usedTracklets(mTimeFrame->getFoundTracklets(pivotRofId, 0).size(), false, mMemoryPool.get()); short startROF{std::max((short)0, static_cast(pivotRofId - mVrtParams[iteration].deltaRof))}; From a67d748f65302d210401629de83031c3b245ecab Mon Sep 17 00:00:00 2001 From: Felix Schlepper Date: Mon, 7 Jul 2025 18:27:33 +0200 Subject: [PATCH 4/8] ITS: deltaROF fix line validation Signed-off-by: Felix Schlepper --- .../ITS/tracking/src/VertexerTraits.cxx | 76 ++++++++++++------- 1 file changed, 49 insertions(+), 27 deletions(-) diff --git a/Detectors/ITSMFT/ITS/tracking/src/VertexerTraits.cxx b/Detectors/ITSMFT/ITS/tracking/src/VertexerTraits.cxx index 8ff4dffef5cdc..5af7dc2c4c44b 100644 --- a/Detectors/ITSMFT/ITS/tracking/src/VertexerTraits.cxx +++ b/Detectors/ITSMFT/ITS/tracking/src/VertexerTraits.cxx @@ -97,21 +97,22 @@ void trackleterKernelHost( } } -void trackletSelectionKernelHost( +static void trackletSelectionKernelHost( const gsl::span clusters0, // 0 const gsl::span clusters1, // 1 gsl::span usedClusters0, // Layer 0 gsl::span usedClusters2, // Layer 2 const gsl::span& tracklets01, const gsl::span& tracklets12, - bounded_vector& usedTracklets, + bounded_vector& usedTracklets, const gsl::span foundTracklets01, const gsl::span foundTracklets12, bounded_vector& lines, const gsl::span& trackletLabels, bounded_vector& linesLabels, - const short pivotRofId, - const short targetRofId, + const short targetRofId0, + const short targetRofId2, + bool safeWrites = false, const float tanLambdaCut = 0.025f, const float phiCut = 0.005f, const int maxTracklets = static_cast(1e2)) @@ -121,16 +122,27 @@ void trackletSelectionKernelHost( int validTracklets{0}; for (int iTracklet12{offset12}; iTracklet12 < offset12 + foundTracklets12[iCurrentLayerClusterIndex]; ++iTracklet12) { for (int iTracklet01{offset01}; iTracklet01 < offset01 + foundTracklets01[iCurrentLayerClusterIndex]; ++iTracklet01) { + if (usedTracklets[iTracklet01]) { + continue; + } + const auto& tracklet01{tracklets01[iTracklet01]}; const auto& tracklet12{tracklets12[iTracklet12]}; - if (tracklet01.rof[0] != targetRofId || tracklet12.rof[1] != targetRofId) { + + if (tracklet01.rof[0] != targetRofId0 || tracklet12.rof[1] != targetRofId2) { continue; } + const float deltaTanLambda{o2::gpu::GPUCommonMath::Abs(tracklet01.tanLambda - tracklet12.tanLambda)}; const float deltaPhi{o2::gpu::GPUCommonMath::Abs(math_utils::smallestAngleDifference(tracklet01.phi, tracklet12.phi))}; - if (!usedTracklets[iTracklet01] && deltaTanLambda < tanLambdaCut && deltaPhi < phiCut && validTracklets != maxTracklets) { - usedClusters0[tracklet01.firstClusterIndex] = true; - usedClusters2[tracklet12.secondClusterIndex] = true; + if (deltaTanLambda < tanLambdaCut && deltaPhi < phiCut && validTracklets != maxTracklets) { + if (safeWrites) { + __atomic_store_n(&usedClusters0[tracklet01.firstClusterIndex], 1, __ATOMIC_RELAXED); + __atomic_store_n(&usedClusters2[tracklet12.secondClusterIndex], 1, __ATOMIC_RELAXED); + } else { + usedClusters0[tracklet01.firstClusterIndex] = 1; + usedClusters2[tracklet12.secondClusterIndex] = 1; + } usedTracklets[iTracklet01] = true; lines.emplace_back(tracklet01, clusters0.data(), clusters1.data()); if (!trackletLabels.empty()) { @@ -330,27 +342,37 @@ void VertexerTraits::computeTrackletMatching(const int iteration) continue; } mTimeFrame->getLines(pivotRofId).reserve(mTimeFrame->getNTrackletsCluster(pivotRofId, 0).size()); - bounded_vector usedTracklets(mTimeFrame->getFoundTracklets(pivotRofId, 0).size(), false, mMemoryPool.get()); + bounded_vector usedTracklets(mTimeFrame->getFoundTracklets(pivotRofId, 0).size(), false, mMemoryPool.get()); short startROF{std::max((short)0, static_cast(pivotRofId - mVrtParams[iteration].deltaRof))}; short endROF{std::min(static_cast(mTimeFrame->getNrof()), static_cast(pivotRofId + mVrtParams[iteration].deltaRof + 1))}; - for (short targetRofId = startROF; targetRofId < endROF; ++targetRofId) { - trackletSelectionKernelHost( - mTimeFrame->getClustersOnLayer(targetRofId, 0), - mTimeFrame->getClustersOnLayer(pivotRofId, 1), - mTimeFrame->getUsedClustersROF(targetRofId, 0), - mTimeFrame->getUsedClustersROF(targetRofId, 2), - mTimeFrame->getFoundTracklets(pivotRofId, 0), - mTimeFrame->getFoundTracklets(pivotRofId, 1), - usedTracklets, - mTimeFrame->getNTrackletsCluster(pivotRofId, 0), - mTimeFrame->getNTrackletsCluster(pivotRofId, 1), - mTimeFrame->getLines(pivotRofId), - mTimeFrame->getLabelsFoundTracklets(pivotRofId, 0), - mTimeFrame->getLinesLabel(pivotRofId), - pivotRofId, - targetRofId, - mVrtParams[iteration].tanLambdaCut, - mVrtParams[iteration].phiCut); + + // needed only if multi-threaded using deltaRof and only at the overlap edges of the ranges + bool safeWrite = mTaskArena->max_concurrency() > 1 && mVrtParams[iteration].deltaRof != 0 && ((Rofs.begin() - startROF < 0) || (endROF - Rofs.end() > 0)); + + for (short targetRofId0 = startROF; targetRofId0 < endROF; ++targetRofId0) { + for (short targetRofId2 = startROF; targetRofId2 < endROF; ++targetRofId2) { + if (std::abs(targetRofId0 - targetRofId2) > mVrtParams[iteration].deltaRof) { // do not allow over 3 ROFs + continue; + } + trackletSelectionKernelHost( + mTimeFrame->getClustersOnLayer(targetRofId0, 0), + mTimeFrame->getClustersOnLayer(pivotRofId, 1), + mTimeFrame->getUsedClustersROF(targetRofId0, 0), + mTimeFrame->getUsedClustersROF(targetRofId2, 2), + mTimeFrame->getFoundTracklets(pivotRofId, 0), + mTimeFrame->getFoundTracklets(pivotRofId, 1), + usedTracklets, + mTimeFrame->getNTrackletsCluster(pivotRofId, 0), + mTimeFrame->getNTrackletsCluster(pivotRofId, 1), + mTimeFrame->getLines(pivotRofId), + mTimeFrame->getLabelsFoundTracklets(pivotRofId, 0), + mTimeFrame->getLinesLabel(pivotRofId), + targetRofId0, + targetRofId2, + safeWrite, + mVrtParams[iteration].tanLambdaCut, + mVrtParams[iteration].phiCut); + } } } }); From f1a2d9eba10fa1997960276e0cb8accba0ceef85 Mon Sep 17 00:00:00 2001 From: Felix Schlepper Date: Sun, 13 Jul 2025 17:01:09 +0200 Subject: [PATCH 5/8] ITS: Vertexer free L12 trklts to save memory Signed-off-by: Felix Schlepper --- Detectors/ITSMFT/ITS/tracking/src/VertexerTraits.cxx | 3 +++ 1 file changed, 3 insertions(+) diff --git a/Detectors/ITSMFT/ITS/tracking/src/VertexerTraits.cxx b/Detectors/ITSMFT/ITS/tracking/src/VertexerTraits.cxx index 5af7dc2c4c44b..414e715d7eacf 100644 --- a/Detectors/ITSMFT/ITS/tracking/src/VertexerTraits.cxx +++ b/Detectors/ITSMFT/ITS/tracking/src/VertexerTraits.cxx @@ -407,6 +407,9 @@ void VertexerTraits::computeTrackletMatching(const int iteration) ln_tre->Write(); trackletFile->Close(); #endif + + // from here on we do not use tracklets from L1-2 anymore, so let's free them + deepVectorClear(mTimeFrame->getTracklets()[1]); } void VertexerTraits::computeVertices(const int iteration) From 1529e7128b02ea50a9c6ad9122f9c3e843622a5f Mon Sep 17 00:00:00 2001 From: Felix Schlepper Date: Tue, 15 Jul 2025 10:55:10 +0200 Subject: [PATCH 6/8] ITS: Tracker early exit if no vertex + disallow trklt > dROF Signed-off-by: Felix Schlepper --- .../ITSMFT/ITS/tracking/src/TimeFrame.cxx | 4 +- Detectors/ITSMFT/ITS/tracking/src/Tracker.cxx | 4 +- .../ITSMFT/ITS/tracking/src/TrackerTraits.cxx | 360 ++++++++++-------- 3 files changed, 214 insertions(+), 154 deletions(-) diff --git a/Detectors/ITSMFT/ITS/tracking/src/TimeFrame.cxx b/Detectors/ITSMFT/ITS/tracking/src/TimeFrame.cxx index 694bfb244e9aa..2e9ce23719f90 100644 --- a/Detectors/ITSMFT/ITS/tracking/src/TimeFrame.cxx +++ b/Detectors/ITSMFT/ITS/tracking/src/TimeFrame.cxx @@ -362,7 +362,7 @@ void TimeFrame::initialise(const int iteration, const TrackingParameter if (iLayer < (int)mCells.size()) { deepVectorClear(mCells[iLayer]); deepVectorClear(mTrackletsLookupTable[iLayer]); - mTrackletsLookupTable[iLayer].resize(mClusters[iLayer + 1].size(), 0); + mTrackletsLookupTable[iLayer].resize(mClusters[iLayer + 1].size() + 1, 0); deepVectorClear(mCellLabels[iLayer]); } @@ -653,8 +653,8 @@ void TimeFrame::wipe() deepVectorClear(mCellsLookupTable); deepVectorClear(mTotVertPerIteration); deepVectorClear(mPrimaryVertices); - deepVectorClear(mROFramesPV); deepVectorClear(mClusters); + deepVectorClear(mTrackletsLookupTable); deepVectorClear(mTrackingFrameInfo); deepVectorClear(mClusterExternalIndices); deepVectorClear(mROFramesClusters); diff --git a/Detectors/ITSMFT/ITS/tracking/src/Tracker.cxx b/Detectors/ITSMFT/ITS/tracking/src/Tracker.cxx index aa93e32e1db9c..ba722c410f95c 100644 --- a/Detectors/ITSMFT/ITS/tracking/src/Tracker.cxx +++ b/Detectors/ITSMFT/ITS/tracking/src/Tracker.cxx @@ -52,7 +52,9 @@ void Tracker::clustersToTracks(const LogFunc& logger, const LogFunc& error) int maxNvertices{-1}; if (mTrkParams[0].PerPrimaryVertexProcessing) { for (int iROF{0}; iROF < mTimeFrame->getNrof(); ++iROF) { - maxNvertices = std::max(maxNvertices, (int)mTimeFrame->getPrimaryVertices(iROF).size()); + int minRof = o2::gpu::CAMath::Max(0, iROF - mTrkParams[0].DeltaROF); + int maxRof = o2::gpu::CAMath::Min(mTimeFrame->getNrof(), iROF + mTrkParams[0].DeltaROF); + maxNvertices = std::max(maxNvertices, (int)mTimeFrame->getPrimaryVertices(minRof, maxRof).size()); } } diff --git a/Detectors/ITSMFT/ITS/tracking/src/TrackerTraits.cxx b/Detectors/ITSMFT/ITS/tracking/src/TrackerTraits.cxx index 014031aaedccb..136ebc647cc38 100644 --- a/Detectors/ITSMFT/ITS/tracking/src/TrackerTraits.cxx +++ b/Detectors/ITSMFT/ITS/tracking/src/TrackerTraits.cxx @@ -21,6 +21,7 @@ #ifdef OPTIMISATION_OUTPUT #include +#include #endif #include @@ -42,8 +43,6 @@ using o2::base::PropagatorF; namespace o2::its { -static constexpr int debugLevel{0}; - struct PassMode { using OnePass = std::integral_constant; using TwoPassCount = std::integral_constant; @@ -73,136 +72,173 @@ void TrackerTraits::computeLayerTracklets(const int iteration, int iROF int endROF{o2::gpu::GPUCommonMath::Min(mTrkParams[iteration].nROFsPerIterations > 0 ? (iROFslice + 1) * mTrkParams[iteration].nROFsPerIterations + mTrkParams[iteration].DeltaROF : mTimeFrame->getNrof(), mTimeFrame->getNrof())}; mTaskArena->execute([&] { - for (int rof0{startROF}; rof0 < endROF; ++rof0) { - int minRof = o2::gpu::CAMath::Max(startROF, rof0 - mTrkParams[iteration].DeltaROF); - int maxRof = o2::gpu::CAMath::Min(endROF - 1, rof0 + mTrkParams[iteration].DeltaROF); - gsl::span primaryVertices = mTrkParams[iteration].UseDiamond ? diamondSpan : mTimeFrame->getPrimaryVertices(minRof, maxRof); + auto forTracklets = [&](auto Tag, int iLayer, int pivotROF, int base, int& offset) -> int { + int minROF = o2::gpu::CAMath::Max(startROF, pivotROF - mTrkParams[iteration].DeltaROF); + int maxROF = o2::gpu::CAMath::Min(endROF - 1, pivotROF + mTrkParams[iteration].DeltaROF); + gsl::span primaryVertices = mTrkParams[iteration].UseDiamond ? diamondSpan : mTimeFrame->getPrimaryVertices(minROF, maxROF); if (primaryVertices.empty()) { - continue; + return 0; + } + const int startVtx = iVertex >= 0 ? iVertex : 0; + const int endVtx = iVertex >= 0 ? o2::gpu::CAMath::Min(iVertex + 1, int(primaryVertices.size())) : int(primaryVertices.size()); + if (endVtx <= startVtx) { + return 0; } - const int startVtx{iVertex >= 0 ? iVertex : 0}; - const int endVtx{iVertex >= 0 ? o2::gpu::CAMath::Min(iVertex + 1, static_cast(primaryVertices.size())) : static_cast(primaryVertices.size())}; - tbb::parallel_for( - tbb::blocked_range(0, mTrkParams[iteration].TrackletsPerRoad()), - [&](const tbb::blocked_range& Layers) { - for (int iLayer = Layers.begin(); iLayer < Layers.end(); ++iLayer) { - gsl::span layer0 = mTimeFrame->getClustersOnLayer(rof0, iLayer); - if (layer0.empty()) { + int localCount = 0; + auto& tracklets = mTimeFrame->getTracklets()[iLayer]; + for (int targetROF0{minROF}; targetROF0 <= maxROF; ++targetROF0) { + if (!mTimeFrame->mMultiplicityCutMask[targetROF0]) { + continue; + } + auto layer0 = mTimeFrame->getClustersOnLayer(targetROF0, iLayer); + if (layer0.empty()) { + continue; + } + const float meanDeltaR = mTrkParams[iteration].LayerRadii[iLayer + 1] - mTrkParams[iteration].LayerRadii[iLayer]; + + for (int iCluster = 0; iCluster < int(layer0.size()); ++iCluster) { + const Cluster& currentCluster = layer0[iCluster]; + const int currentSortedIndex = mTimeFrame->getSortedIndex(targetROF0, iLayer, iCluster); + if (mTimeFrame->isClusterUsed(iLayer, currentCluster.clusterId)) { + continue; + } + const float inverseR0 = 1.f / currentCluster.radius; + + for (int iV = startVtx; iV < endVtx; ++iV) { + const auto& pv = primaryVertices[iV]; + if (pv.isFlagSet(Vertex::Flags::UPCMode) && iteration != 3) { continue; } - float meanDeltaR{mTrkParams[iteration].LayerRadii[iLayer + 1] - mTrkParams[iteration].LayerRadii[iLayer]}; - - const int currentLayerClustersNum{static_cast(layer0.size())}; - for (int iCluster{0}; iCluster < currentLayerClustersNum; ++iCluster) { - const Cluster& currentCluster{layer0[iCluster]}; - const int currentSortedIndex{mTimeFrame->getSortedIndex(rof0, iLayer, iCluster)}; + const float resolution = o2::gpu::CAMath::Sqrt(math_utils::Sq(mTimeFrame->getPositionResolution(iLayer)) + math_utils::Sq(mTrkParams[iteration].PVres) / float(pv.getNContributors())); + const float tanLambda = (currentCluster.zCoordinate - pv.getZ()) * inverseR0; + const float zAtRmin = tanLambda * (mTimeFrame->getMinR(iLayer + 1) - currentCluster.radius) + currentCluster.zCoordinate; + const float zAtRmax = tanLambda * (mTimeFrame->getMaxR(iLayer + 1) - currentCluster.radius) + currentCluster.zCoordinate; + const float sqInvDeltaZ0 = 1.f / (math_utils::Sq(currentCluster.zCoordinate - pv.getZ()) + constants::Tolerance); + const float sigmaZ = o2::gpu::CAMath::Sqrt( + math_utils::Sq(resolution) * math_utils::Sq(tanLambda) * ((math_utils::Sq(inverseR0) + sqInvDeltaZ0) * math_utils::Sq(meanDeltaR) + 1.f) + math_utils::Sq(meanDeltaR * mTimeFrame->getMSangle(iLayer))); + + auto bins = getBinsRect(currentCluster, iLayer + 1, zAtRmin, zAtRmax, sigmaZ * mTrkParams[iteration].NSigmaCut, mTimeFrame->getPhiCut(iLayer)); + if (bins.x == 0 && bins.y == 0 && bins.z == 0 && bins.w == 0) { + continue; + } + int phiBinsNum = bins.w - bins.y + 1; + if (phiBinsNum < 0) { + phiBinsNum += mTrkParams[iteration].PhiBins; + } - if (mTimeFrame->isClusterUsed(iLayer, currentCluster.clusterId)) { + for (int targetROF1{minROF}; targetROF1 <= maxROF; ++targetROF1) { + if (!mTimeFrame->mMultiplicityCutMask[targetROF1] || std::abs(targetROF0 - targetROF1) > mTrkParams[iteration].DeltaROF) { continue; } - const float inverseR0{1.f / currentCluster.radius}; - - for (int iV{startVtx}; iV < endVtx; ++iV) { - const auto& primaryVertex{primaryVertices[iV]}; - if (primaryVertex.isFlagSet(2) && iteration != 3) { - continue; - } - const float resolution = o2::gpu::CAMath::Sqrt(math_utils::Sq(mTrkParams[iteration].PVres) / primaryVertex.getNContributors() + math_utils::Sq(mTimeFrame->getPositionResolution(iLayer))); - - const float tanLambda{(currentCluster.zCoordinate - primaryVertex.getZ()) * inverseR0}; - - const float zAtRmin{tanLambda * (mTimeFrame->getMinR(iLayer + 1) - currentCluster.radius) + currentCluster.zCoordinate}; - const float zAtRmax{tanLambda * (mTimeFrame->getMaxR(iLayer + 1) - currentCluster.radius) + currentCluster.zCoordinate}; - - const float sqInverseDeltaZ0{1.f / (math_utils::Sq(currentCluster.zCoordinate - primaryVertex.getZ()) + 2.e-8f)}; /// protecting from overflows adding the detector resolution - const float sigmaZ{o2::gpu::CAMath::Sqrt(math_utils::Sq(resolution) * math_utils::Sq(tanLambda) * ((math_utils::Sq(inverseR0) + sqInverseDeltaZ0) * math_utils::Sq(meanDeltaR) + 1.f) + math_utils::Sq(meanDeltaR * mTimeFrame->getMSangle(iLayer)))}; - - const int4 selectedBinsRect{getBinsRect(currentCluster, iLayer + 1, zAtRmin, zAtRmax, sigmaZ * mTrkParams[iteration].NSigmaCut, mTimeFrame->getPhiCut(iLayer))}; - if (selectedBinsRect.x == 0 && selectedBinsRect.y == 0 && selectedBinsRect.z == 0 && selectedBinsRect.w == 0) { - continue; - } - - int phiBinsNum{selectedBinsRect.w - selectedBinsRect.y + 1}; - - if (phiBinsNum < 0) { - phiBinsNum += mTrkParams[iteration].PhiBins; - } - - for (int rof1{minRof}; rof1 <= maxRof; ++rof1) { - if (!mTimeFrame->mMultiplicityCutMask[rof1]) { - continue; + auto layer1 = mTimeFrame->getClustersOnLayer(targetROF1, iLayer + 1); + if (layer1.empty()) { + continue; + } + for (int iPhi = 0; iPhi < phiBinsNum; ++iPhi) { + int iPhiBin = (bins.y + iPhi) % mTrkParams[iteration].PhiBins; + int firstBinIdx = mTimeFrame->mIndexTableUtils.getBinIndex(bins.x, iPhiBin); + int maxBinIdx = firstBinIdx + (bins.z - bins.x) + 1; + int firstRow = mTimeFrame->getIndexTable(targetROF1, iLayer + 1)[firstBinIdx]; + int lastRow = mTimeFrame->getIndexTable(targetROF1, iLayer + 1)[maxBinIdx]; + for (int iNext = firstRow; iNext < lastRow; ++iNext) { + if (iNext >= int(layer1.size())) { + break; } - auto layer1 = mTimeFrame->getClustersOnLayer(rof1, iLayer + 1); - if (layer1.empty()) { + const Cluster& nextCluster = layer1[iNext]; + if (mTimeFrame->isClusterUsed(iLayer + 1, nextCluster.clusterId)) { continue; } - for (int iPhiCount{0}; iPhiCount < phiBinsNum; iPhiCount++) { - int iPhiBin = (selectedBinsRect.y + iPhiCount) % mTrkParams[iteration].PhiBins; - const int firstBinIndex{mTimeFrame->mIndexTableUtils.getBinIndex(selectedBinsRect.x, iPhiBin)}; - const int maxBinIndex{firstBinIndex + selectedBinsRect.z - selectedBinsRect.x + 1}; - if constexpr (debugLevel) { - if (firstBinIndex < 0 || firstBinIndex > mTimeFrame->getIndexTable(rof1, iLayer + 1).size() || - maxBinIndex < 0 || maxBinIndex > mTimeFrame->getIndexTable(rof1, iLayer + 1).size()) { - std::cout << iLayer << "\t" << iCluster << "\t" << zAtRmin << "\t" << zAtRmax << "\t" << sigmaZ * mTrkParams[iteration].NSigmaCut << "\t" << mTimeFrame->getPhiCut(iLayer) << std::endl; - std::cout << currentCluster.zCoordinate << "\t" << primaryVertex.getZ() << "\t" << currentCluster.radius << std::endl; - std::cout << mTimeFrame->getMinR(iLayer + 1) << "\t" << currentCluster.radius << "\t" << currentCluster.zCoordinate << std::endl; - std::cout << "Illegal access to IndexTable " << firstBinIndex << "\t" << maxBinIndex << "\t" << selectedBinsRect.z << "\t" << selectedBinsRect.x << std::endl; - exit(1); - } - } - const int firstRowClusterIndex = mTimeFrame->getIndexTable(rof1, iLayer + 1)[firstBinIndex]; - const int maxRowClusterIndex = mTimeFrame->getIndexTable(rof1, iLayer + 1)[maxBinIndex]; - for (int iNextCluster{firstRowClusterIndex}; iNextCluster < maxRowClusterIndex; ++iNextCluster) { - if (iNextCluster >= (int)layer1.size()) { - break; - } - - const Cluster& nextCluster{layer1[iNextCluster]}; - if (mTimeFrame->isClusterUsed(iLayer + 1, nextCluster.clusterId)) { - continue; - } - - const float deltaPhi{o2::gpu::GPUCommonMath::Abs(currentCluster.phi - nextCluster.phi)}; - const float deltaZ{o2::gpu::GPUCommonMath::Abs(tanLambda * (nextCluster.radius - currentCluster.radius) + - currentCluster.zCoordinate - nextCluster.zCoordinate)}; + float deltaPhi = o2::gpu::GPUCommonMath::Abs(currentCluster.phi - nextCluster.phi); + float deltaZ = o2::gpu::GPUCommonMath::Abs((tanLambda * (nextCluster.radius - currentCluster.radius)) + currentCluster.zCoordinate - nextCluster.zCoordinate); #ifdef OPTIMISATION_OUTPUT - MCCompLabel label; - int currentId{currentCluster.clusterId}; - int nextId{nextCluster.clusterId}; - for (auto& lab1 : mTimeFrame->getClusterLabels(iLayer, currentId)) { - for (auto& lab2 : mTimeFrame->getClusterLabels(iLayer + 1, nextId)) { - if (lab1 == lab2 && lab1.isValid()) { - label = lab1; - break; - } - } - if (label.isValid()) { - break; - } + MCCompLabel label; + int currentId{currentCluster.clusterId}; + int nextId{nextCluster.clusterId}; + for (auto& lab1 : mTimeFrame->getClusterLabels(iLayer, currentId)) { + for (auto& lab2 : mTimeFrame->getClusterLabels(iLayer + 1, nextId)) { + if (lab1 == lab2 && lab1.isValid()) { + label = lab1; + break; } - off << std::format("{}\t{:d}\t{}\t{}\t{}\t{}", iLayer, label.isValid(), (tanLambda * (nextCluster.radius - currentCluster.radius) + currentCluster.zCoordinate - nextCluster.zCoordinate) / sigmaZ, tanLambda, resolution, sigmaZ) << std::endl; + } + if (label.isValid()) { + break; + } + } + off << std::format("{}\t{:d}\t{}\t{}\t{}\t{}", iLayer, label.isValid(), (tanLambda * (nextCluster.radius - currentCluster.radius) + currentCluster.zCoordinate - nextCluster.zCoordinate) / sigmaZ, tanLambda, resolution, sigmaZ) << std::endl; #endif - if (deltaZ / sigmaZ < mTrkParams[iteration].NSigmaCut && - (deltaPhi < mTimeFrame->getPhiCut(iLayer) || - o2::gpu::GPUCommonMath::Abs(deltaPhi - o2::constants::math::TwoPI) < mTimeFrame->getPhiCut(iLayer))) { - if (iLayer > 0) { - mTimeFrame->getTrackletsLookupTable()[iLayer - 1][currentSortedIndex]++; - } - const float phi{o2::gpu::GPUCommonMath::ATan2(currentCluster.yCoordinate - nextCluster.yCoordinate, - currentCluster.xCoordinate - nextCluster.xCoordinate)}; - const float tanL{(currentCluster.zCoordinate - nextCluster.zCoordinate) / - (currentCluster.radius - nextCluster.radius)}; - mTimeFrame->getTracklets()[iLayer].emplace_back(currentSortedIndex, mTimeFrame->getSortedIndex(rof1, iLayer + 1, iNextCluster), tanL, phi, rof0, rof1); - } + if (deltaZ / sigmaZ < mTrkParams[iteration].NSigmaCut && + (deltaPhi < mTimeFrame->getPhiCut(iLayer) || + o2::gpu::GPUCommonMath::Abs(deltaPhi - o2::constants::math::TwoPI) < mTimeFrame->getPhiCut(iLayer))) { + float phi = o2::gpu::GPUCommonMath::ATan2(currentCluster.yCoordinate - nextCluster.yCoordinate, currentCluster.xCoordinate - nextCluster.xCoordinate); + float tanL = (currentCluster.zCoordinate - nextCluster.zCoordinate) / (currentCluster.radius - nextCluster.radius); + if constexpr (decltype(Tag)::value == PassMode::OnePass::value) { + tracklets.emplace_back(currentSortedIndex, mTimeFrame->getSortedIndex(targetROF1, iLayer + 1, iNext), tanL, phi, targetROF0, targetROF1); + } else if constexpr (decltype(Tag)::value == PassMode::TwoPassCount::value) { + ++localCount; + } else if constexpr (decltype(Tag)::value == PassMode::TwoPassInsert::value) { + const int idx = base + offset++; + tracklets[idx] = Tracklet(currentSortedIndex, mTimeFrame->getSortedIndex(targetROF1, iLayer + 1, iNext), tanL, phi, targetROF0, targetROF1); } } } } } } + } + } + return localCount; + }; + + int dummy{0}; + if (mTaskArena->max_concurrency() <= 1) { + for (int pivotROF{startROF}; pivotROF < endROF; ++pivotROF) { + for (int iLayer{0}; iLayer < mTrkParams[iteration].TrackletsPerRoad(); ++iLayer) { + forTracklets(PassMode::OnePass{}, iLayer, pivotROF, 0, dummy); + } + } + } else { + bounded_vector> perROFCount(mTrkParams[iteration].TrackletsPerRoad(), bounded_vector(endROF - startROF + 1, 0, mMemoryPool.get()), mMemoryPool.get()); + tbb::parallel_for( + tbb::blocked_range2d(0, mTrkParams[iteration].TrackletsPerRoad(), 1, + startROF, endROF, 1), + [&](auto const& Range) { + for (int iLayer{Range.rows().begin()}; iLayer < Range.rows().end(); ++iLayer) { + for (int pivotROF = Range.cols().begin(); pivotROF < Range.cols().end(); ++pivotROF) { + perROFCount[iLayer][pivotROF - startROF] = forTracklets(PassMode::TwoPassCount{}, iLayer, pivotROF, 0, dummy); + } + } + }); + + tbb::parallel_for( + tbb::blocked_range(0, mTrkParams[iteration].TrackletsPerRoad()), + [&](auto const& Layers) { + for (int iLayer{Layers.begin()}; iLayer < Layers.end(); ++iLayer) { + std::exclusive_scan(perROFCount[iLayer].begin(), perROFCount[iLayer].end(), perROFCount[iLayer].begin(), 0); + mTimeFrame->getTracklets()[iLayer].resize(perROFCount[iLayer].back()); + } + }); + + tbb::parallel_for( + tbb::blocked_range2d(0, mTrkParams[iteration].TrackletsPerRoad(), 1, + startROF, endROF, 1), + [&](auto const& Range) { + for (int iLayer{Range.rows().begin()}; iLayer < Range.rows().end(); ++iLayer) { + if (perROFCount[iLayer].back() == 0) { + continue; + } + for (int pivotROF = Range.cols().begin(); pivotROF < Range.cols().end(); ++pivotROF) { + int baseIdx = perROFCount[iLayer][pivotROF - startROF]; + if (baseIdx == perROFCount[iLayer][pivotROF - startROF + 1]) { + continue; + } + int localIdx = 0; + forTracklets(PassMode::TwoPassInsert{}, iLayer, pivotROF, baseIdx, localIdx); + } + } }); } @@ -223,42 +259,43 @@ void TrackerTraits::computeLayerTracklets(const int iteration, int iROF trkl.shrink_to_fit(); if (iLayer > 0) { /// recalculate lut auto& lut{mTimeFrame->getTrackletsLookupTable()[iLayer - 1]}; - std::fill(lut.begin(), lut.end(), 0); - if (trkl.empty()) { - return; - } - for (const auto& tkl : trkl) { - lut[tkl.firstClusterIndex]++; + if (!trkl.empty()) { + for (const auto& tkl : trkl) { + lut[tkl.firstClusterIndex + 1]++; + } + std::inclusive_scan(lut.begin(), lut.end(), lut.begin()); } - std::exclusive_scan(lut.begin(), lut.end(), lut.begin(), 0); - lut.push_back(trkl.size()); } } }); - }); - /// Create tracklets labels - if (mTimeFrame->hasMCinformation()) { - for (int iLayer{0}; iLayer < mTrkParams[iteration].TrackletsPerRoad(); ++iLayer) { - for (auto& trk : mTimeFrame->getTracklets()[iLayer]) { - MCCompLabel label; - int currentId{mTimeFrame->getClusters()[iLayer][trk.firstClusterIndex].clusterId}; - int nextId{mTimeFrame->getClusters()[iLayer + 1][trk.secondClusterIndex].clusterId}; - for (auto& lab1 : mTimeFrame->getClusterLabels(iLayer, currentId)) { - for (auto& lab2 : mTimeFrame->getClusterLabels(iLayer + 1, nextId)) { - if (lab1 == lab2 && lab1.isValid()) { - label = lab1; - break; + /// Create tracklets labels + if (mTimeFrame->hasMCinformation()) { + tbb::parallel_for( + tbb::blocked_range(0, mTrkParams[iteration].TrackletsPerRoad()), + [&](const tbb::blocked_range& Layers) { + for (int iLayer = Layers.begin(); iLayer < Layers.end(); ++iLayer) { + for (auto& trk : mTimeFrame->getTracklets()[iLayer]) { + MCCompLabel label; + int currentId{mTimeFrame->getClusters()[iLayer][trk.firstClusterIndex].clusterId}; + int nextId{mTimeFrame->getClusters()[iLayer + 1][trk.secondClusterIndex].clusterId}; + for (const auto& lab1 : mTimeFrame->getClusterLabels(iLayer, currentId)) { + for (const auto& lab2 : mTimeFrame->getClusterLabels(iLayer + 1, nextId)) { + if (lab1 == lab2 && lab1.isValid()) { + label = lab1; + break; + } + } + if (label.isValid()) { + break; + } + } + mTimeFrame->getTrackletsLabel(iLayer).emplace_back(label); } } - if (label.isValid()) { - break; - } - } - mTimeFrame->getTrackletsLabel(iLayer).emplace_back(label); - } + }); } - } + }); } template @@ -283,17 +320,19 @@ void TrackerTraits::computeLayerCells(const int iteration) auto forTrackletCells = [&](auto Tag, int iLayer, bounded_vector& layerCells, int iTracklet, int offset = 0) -> int { const Tracklet& currentTracklet{mTimeFrame->getTracklets()[iLayer][iTracklet]}; const int nextLayerClusterIndex{currentTracklet.secondClusterIndex}; - const int nextLayerFirstTrackletIndex{ - mTimeFrame->getTrackletsLookupTable()[iLayer][nextLayerClusterIndex]}; - const int nextLayerLastTrackletIndex{ - mTimeFrame->getTrackletsLookupTable()[iLayer][nextLayerClusterIndex + 1]}; - + const int nextLayerFirstTrackletIndex{mTimeFrame->getTrackletsLookupTable()[iLayer][nextLayerClusterIndex]}; + const int nextLayerLastTrackletIndex{mTimeFrame->getTrackletsLookupTable()[iLayer][nextLayerClusterIndex + 1]}; int foundCells{0}; for (int iNextTracklet{nextLayerFirstTrackletIndex}; iNextTracklet < nextLayerLastTrackletIndex; ++iNextTracklet) { + const Tracklet& nextTracklet{mTimeFrame->getTracklets()[iLayer + 1][iNextTracklet]}; + const auto& nextLbl = mTimeFrame->getTrackletsLabel(iLayer + 1)[iNextTracklet]; + bool print = false; if (mTimeFrame->getTracklets()[iLayer + 1][iNextTracklet].firstClusterIndex != nextLayerClusterIndex) { break; } - const Tracklet& nextTracklet{mTimeFrame->getTracklets()[iLayer + 1][iNextTracklet]}; + if (mTrkParams[iteration].DeltaROF && currentTracklet.getSpanRof(nextTracklet) > mTrkParams[iteration].DeltaROF) { // TODO this has to be improved for the staggering + continue; + } const float deltaTanLambda{std::abs(currentTracklet.tanLambda - nextTracklet.tanLambda)}; #ifdef OPTIMISATION_OUTPUT @@ -420,19 +459,13 @@ void TrackerTraits::computeLayerCells(const int iteration) /// Create cells labels if (mTimeFrame->hasMCinformation()) { for (int iLayer{0}; iLayer < mTrkParams[iteration].CellsPerRoad(); ++iLayer) { - for (auto& cell : mTimeFrame->getCells()[iLayer]) { + for (const auto& cell : mTimeFrame->getCells()[iLayer]) { MCCompLabel currentLab{mTimeFrame->getTrackletsLabel(iLayer)[cell.getFirstTrackletIndex()]}; MCCompLabel nextLab{mTimeFrame->getTrackletsLabel(iLayer + 1)[cell.getSecondTrackletIndex()]}; mTimeFrame->getCellsLabel(iLayer).emplace_back(currentLab == nextLab ? currentLab : MCCompLabel()); } } } - - if constexpr (debugLevel) { - for (int iLayer{0}; iLayer < mTrkParams[iteration].CellsPerRoad(); ++iLayer) { - std::cout << "Cells on layer " << iLayer << " " << mTimeFrame->getCells()[iLayer].size() << std::endl; - } - } } template @@ -470,6 +503,16 @@ void TrackerTraits::findCellsNeighbours(const int iteration) break; } + if (mTrkParams[iteration].DeltaROF) { // TODO this has to be improved for the staggering + const auto& trkl00 = mTimeFrame->getTracklets()[iLayer][currentCellSeed.getFirstTrackletIndex()]; + const auto& trkl01 = mTimeFrame->getTracklets()[iLayer + 1][currentCellSeed.getSecondTrackletIndex()]; + const auto& trkl10 = mTimeFrame->getTracklets()[iLayer + 1][nextCellSeed.getFirstTrackletIndex()]; + const auto& trkl11 = mTimeFrame->getTracklets()[iLayer + 2][nextCellSeed.getSecondTrackletIndex()]; + if ((std::max({trkl00.getMaxRof(), trkl01.getMaxRof(), trkl10.getMaxRof(), trkl11.getMaxRof()}) - std::min({trkl00.getMinRof(), trkl01.getMinRof(), trkl10.getMinRof(), trkl10.getMinRof()})) > mTrkParams[0].DeltaROF) { + continue; + } + } + if (!nextCellSeed.rotate(currentCellSeed.getAlpha()) || !nextCellSeed.propagateTo(currentCellSeed.getX(), getBz())) { continue; @@ -613,6 +656,21 @@ void TrackerTraits::processNeighbours(int iLayer, int iLevel, const bou CA_DEBUGGER(failed[0]++); continue; } + if (mTrkParams[0].DeltaROF) { // TODO this has to be improved for the staggering + const auto& trklNeigh = mTimeFrame->getTracklets()[iLayer - 1][neighbourCell.getFirstTrackletIndex()]; + short minRof{std::numeric_limits::max()}, maxRof{std::numeric_limits::min()}; + for (int iLayer{0}; iLayer < mTrkParams[0].NLayers; ++iLayer) { + if (const auto clsId = currentCell.getCluster(iLayer); clsId != constants::UnusedIndex) { + const short clsROF = mTimeFrame->getClusterROF(iLayer, clsId); + minRof = std::min(minRof, clsROF); + maxRof = std::max(maxRof, clsROF); + } + } + if ((std::max(trklNeigh.getMaxRof(), maxRof) - std::min(trklNeigh.getMinRof(), minRof)) > mTrkParams[0].DeltaROF) { + continue; + } + } + /// Let's start the fitting procedure CellSeed seed{currentCell}; const auto& trHit = mTimeFrame->getTrackingFrameInfoOnLayer(iLayer - 1)[neighbourCell.getFirstClusterIndex()]; From e76eb64c24641e5b1e76fb77c1cc8af8908ba0ca Mon Sep 17 00:00:00 2001 From: Felix Schlepper Date: Sat, 19 Jul 2025 14:57:00 +0200 Subject: [PATCH 7/8] ITS: drof class mods Signed-off-by: Felix Schlepper --- .../ITS/tracking/include/ITStracking/Cell.h | 28 ++++++++----------- .../include/ITStracking/ClusterLines.h | 8 +++--- .../tracking/include/ITStracking/Tracklet.h | 22 +++++---------- 3 files changed, 23 insertions(+), 35 deletions(-) diff --git a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/Cell.h b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/Cell.h index d81ba4426ca55..fcea96abbfa82 100644 --- a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/Cell.h +++ b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/Cell.h @@ -15,12 +15,8 @@ #ifndef TRACKINGITSU_INCLUDE_CACELL_H_ #define TRACKINGITSU_INCLUDE_CACELL_H_ -#ifndef GPUCA_GPUCODE_DEVICE -#include -#include -#include -#endif +#include "ITStracking/Constants.h" #include "GPUCommonDef.h" namespace o2::its @@ -39,12 +35,12 @@ class Cell final GPUhd() int* getLevelPtr() { return &mLevel; } private: - int mFirstClusterIndex{0}; - int mSecondClusterIndex{0}; - int mThirdClusterIndex{0}; - int mFirstTrackletIndex{0}; - int mSecondTrackletIndex{0}; - int mLevel{0}; + int mFirstClusterIndex{constants::UnusedIndex}; + int mSecondClusterIndex{constants::UnusedIndex}; + int mThirdClusterIndex{constants::UnusedIndex}; + int mFirstTrackletIndex{constants::UnusedIndex}; + int mSecondTrackletIndex{constants::UnusedIndex}; + int mLevel{constants::UnusedIndex}; }; class CellSeed final : public o2::track::TrackParCovF @@ -82,14 +78,14 @@ class CellSeed final : public o2::track::TrackParCovF GPUhd() int getCluster(int i) const { return mClusters[i]; } GPUhd() void printCell() const { - printf("trkl: %d, %d\t lvl: %d\t chi2: %f\n", mTracklets[0], mTracklets[1], mLevel, mChi2); + printf("trkl: %d, %d\t lvl: %d\t chi2: %f\tcls: [%d | %d | %d | %d | %d | %d | %d]\n", mTracklets[0], mTracklets[1], mLevel, mChi2, mClusters[0], mClusters[1], mClusters[2], mClusters[3], mClusters[4], mClusters[5], mClusters[6]); } private: - float mChi2 = 0.f; - int mLevel = 0; - int mTracklets[2] = {-1, -1}; - int mClusters[7] = {-1, -1, -1, -1, -1, -1, -1}; + float mChi2 = -999.f; + int mLevel = constants::UnusedIndex; + int mTracklets[2] = {constants::UnusedIndex, constants::UnusedIndex}; + int mClusters[7] = {constants::UnusedIndex, constants::UnusedIndex, constants::UnusedIndex, constants::UnusedIndex, constants::UnusedIndex, constants::UnusedIndex, constants::UnusedIndex}; }; } // namespace o2::its diff --git a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/ClusterLines.h b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/ClusterLines.h index 3ffeda9adcfd5..0e7ad474ae455 100644 --- a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/ClusterLines.h +++ b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/ClusterLines.h @@ -42,8 +42,8 @@ struct Line final { bool operator!=(const Line&) const; short getMinROF() const { return rof[0] < rof[1] ? rof[0] : rof[1]; } - float originPoint[3] = {0}; - float cosinesDirector[3] = {0}; + float originPoint[3] = {0, 0, 0}; + float cosinesDirector[3] = {0, 0, 0}; // float weightMatrix[6] = {1., 0., 0., 1., 0., 1.}; // weightMatrix is a symmetric matrix internally stored as // 0 --> row = 0, col = 0 @@ -52,7 +52,7 @@ struct Line final { // 3 --> 1,1 // 4 --> 1,2 // 5 --> 2,2 - short rof[2] = {-1, -1}; + short rof[2] = {constants::UnusedIndex, constants::UnusedIndex}; ClassDefNV(Line, 1); }; @@ -207,7 +207,7 @@ class ClusterLines final std::array mRMS2 = {0.f}; // symmetric matrix: diagonal is RMS2 float mAvgDistance2 = 0.f; // substitute for chi2 int mROFWeight = 0; // rof weight for voting - short mROF = -1; // rof + short mROF = constants::UnusedIndex; // rof }; } // namespace o2::its diff --git a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/Tracklet.h b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/Tracklet.h index ba3474e6e86c6..5741a9fc65947 100644 --- a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/Tracklet.h +++ b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/Tracklet.h @@ -16,6 +16,7 @@ #ifndef TRACKINGITS_INCLUDE_TRACKLET_H_ #define TRACKINGITS_INCLUDE_TRACKLET_H_ +#include "ITStracking/Constants.h" #include "ITStracking/Cluster.h" #include "GPUCommonRtypes.h" #include "GPUCommonMath.h" @@ -41,9 +42,10 @@ struct Tracklet final { { return firstClusterIndex < 0 || secondClusterIndex < 0; } + GPUhdi() auto getMinRof() const noexcept { return o2::gpu::CAMath::Min(rof[0], rof[1]); } + GPUhdi() auto getMaxRof() const noexcept { return o2::gpu::CAMath::Max(rof[0], rof[1]); } GPUhdi() auto getDeltaRof() const { return rof[1] - rof[0]; } - GPUhdi() void dump() const; - GPUhdi() void dump(const int, const int) const; + GPUhdi() auto getSpanRof(const Tracklet& o) const noexcept { return o2::gpu::CAMath::Max(getMaxRof(), o.getMaxRof()) - o2::gpu::CAMath::Min(getMinRof(), o.getMinRof()); } GPUhdi() unsigned char operator<(const Tracklet&) const; #if !defined(GPUCA_NO_FMT) && !defined(GPUCA_GPUCODE_DEVICE) std::string asString() const @@ -53,11 +55,11 @@ struct Tracklet final { void print() const { LOG(info) << asString(); } #endif - int firstClusterIndex{-1}; - int secondClusterIndex{-1}; + int firstClusterIndex{constants::UnusedIndex}; + int secondClusterIndex{constants::UnusedIndex}; float tanLambda{-999}; float phi{-999}; - short rof[2] = {-1, -1}; + short rof[2] = {constants::UnusedIndex, constants::UnusedIndex}; ClassDefNV(Tracklet, 1); }; @@ -93,16 +95,6 @@ GPUhdi() unsigned char Tracklet::operator<(const Tracklet& t) const return true; } -GPUhdi() void Tracklet::dump(const int offsetFirst, const int offsetSecond) const -{ - printf("fClIdx: %d sClIdx: %d rof1: %hu rof2: %hu\n", firstClusterIndex + offsetFirst, secondClusterIndex + offsetSecond, rof[0], rof[1]); -} - -GPUhdi() void Tracklet::dump() const -{ - printf("fClIdx: %d sClIdx: %d rof1: %hu rof2: %hu\n", firstClusterIndex, secondClusterIndex, rof[0], rof[1]); -} - } // namespace o2::its #endif /* TRACKINGITS_INCLUDE_TRACKLET_H_ */ From ddf1fd0aef447df01d53b2aa0af2ec6ac79f03de Mon Sep 17 00:00:00 2001 From: Felix Schlepper Date: Wed, 23 Jul 2025 14:05:09 +0200 Subject: [PATCH 8/8] ITS: Vertexer add debug dump Signed-off-by: Felix Schlepper --- .../include/ITStracking/BoundedAllocator.h | 8 + .../include/ITStracking/VertexerTraits.h | 5 + .../ITS/tracking/src/VertexerTraits.cxx | 366 +++++++++++++----- 3 files changed, 275 insertions(+), 104 deletions(-) diff --git a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/BoundedAllocator.h b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/BoundedAllocator.h index eced0c64c73a5..ac9f72089602d 100644 --- a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/BoundedAllocator.h +++ b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/BoundedAllocator.h @@ -174,6 +174,14 @@ void clearResizeBoundedArray(std::array, S>& arr, size_t size, } } +template +std::vector toSTDVector(const bounded_vector& b) +{ + std::vector t(b.size()); + std::copy(b.cbegin(), b.cend(), t.begin()); + return t; +} + } // namespace o2::its #endif diff --git a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/VertexerTraits.h b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/VertexerTraits.h index 54424136fcfe1..1213ad0a423b8 100644 --- a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/VertexerTraits.h +++ b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/VertexerTraits.h @@ -119,6 +119,11 @@ class VertexerTraits private: std::shared_ptr mMemoryPool; std::shared_ptr mTaskArena; + + // debug output + void debugComputeTracklets(int iteration); + void debugComputeTrackletMatching(int iteration); + void debugComputeVertices(int iteration); }; inline void VertexerTraits::initialise(const TrackingParameters& trackingParams, const int iteration) diff --git a/Detectors/ITSMFT/ITS/tracking/src/VertexerTraits.cxx b/Detectors/ITSMFT/ITS/tracking/src/VertexerTraits.cxx index 414e715d7eacf..bcafa98972d78 100644 --- a/Detectors/ITSMFT/ITS/tracking/src/VertexerTraits.cxx +++ b/Detectors/ITSMFT/ITS/tracking/src/VertexerTraits.cxx @@ -26,18 +26,12 @@ #include "Steer/MCKinematicsReader.h" #include "ITSMFTBase/DPLAlpideParam.h" #include "DetectorsRaw/HBFUtils.h" - -#ifdef VTX_DEBUG -#include "TTree.h" -#include "TFile.h" -#include -#include -#endif +#include "CommonUtils/TreeStreamRedirector.h" using namespace o2::its; template -void trackleterKernelHost( +static void trackleterKernelHost( const gsl::span& clustersNextLayer, // 0 2 const gsl::span& clustersCurrentLayer, // 1 1 const gsl::span& usedClustersNextLayer, // 0 2 @@ -285,46 +279,7 @@ void VertexerTraits::computeTracklets(const int iteration) } #ifdef VTX_DEBUG - // Dump on file - TFile* trackletFile = TFile::Open("artefacts_tf.root", "recreate"); - TTree* tr_tre = new TTree("tracklets", "tf"); - std::vector trkl_vec_0(0); - std::vector trkl_vec_1(0); - std::vector clus0(0); - std::vector clus1(0); - std::vector clus2(0); - tr_tre->Branch("Tracklets0", &trkl_vec_0); - tr_tre->Branch("Tracklets1", &trkl_vec_1); - for (int rofId{0}; rofId < mTimeFrame->getNrof(); ++rofId) { - trkl_vec_0.clear(); - trkl_vec_1.clear(); - for (auto& tr : mTimeFrame->getFoundTracklets(rofId, 0)) { - trkl_vec_0.push_back(tr); - } - for (auto& tr : mTimeFrame->getFoundTracklets(rofId, 1)) { - trkl_vec_1.push_back(tr); - } - tr_tre->Fill(); - } - trackletFile->cd(); - tr_tre->Write(); - trackletFile->Close(); - - std::ofstream out01("NTC01_cpu.txt"), out12("NTC12_cpu.txt"); - for (int iRof{0}; iRof < mTimeFrame->getNrof(); ++iRof) { - out01 << "ROF: " << iRof << std::endl; - out12 << "ROF: " << iRof << std::endl; - std::copy(mTimeFrame->getNTrackletsCluster(iRof, 0).begin(), mTimeFrame->getNTrackletsCluster(iRof, 0).end(), std::ostream_iterator(out01, "\t")); - out01 << std::endl; - std::copy(mTimeFrame->getExclusiveNTrackletsCluster(iRof, 0).begin(), mTimeFrame->getExclusiveNTrackletsCluster(iRof, 0).end(), std::ostream_iterator(out01, "\t")); - std::copy(mTimeFrame->getNTrackletsCluster(iRof, 1).begin(), mTimeFrame->getNTrackletsCluster(iRof, 1).end(), std::ostream_iterator(out12, "\t")); - out12 << std::endl; - std::copy(mTimeFrame->getExclusiveNTrackletsCluster(iRof, 1).begin(), mTimeFrame->getExclusiveNTrackletsCluster(iRof, 1).end(), std::ostream_iterator(out12, "\t")); - out01 << std::endl; - out12 << std::endl; - } - out01.close(); - out12.close(); + debugComputeTracklets(iteration); #endif } @@ -379,33 +334,7 @@ void VertexerTraits::computeTrackletMatching(const int iteration) }); #ifdef VTX_DEBUG - TFile* trackletFile = TFile::Open("artefacts_tf.root", "update"); - TTree* ln_tre = new TTree("lines", "tf"); - std::vector lines_vec(0); - std::vector nTrackl01(0); - std::vector nTrackl12(0); - ln_tre->Branch("Lines", &lines_vec); - ln_tre->Branch("NTrackletCluster01", &nTrackl01); - ln_tre->Branch("NTrackletCluster12", &nTrackl12); - for (int rofId{0}; rofId < mTimeFrame->getNrof(); ++rofId) { - lines_vec.clear(); - nTrackl01.clear(); - nTrackl12.clear(); - for (auto& ln : mTimeFrame->getLines(rofId)) { - lines_vec.push_back(ln); - } - for (auto& n : mTimeFrame->getNTrackletsCluster(rofId, 0)) { - nTrackl01.push_back(n); - } - for (auto& n : mTimeFrame->getNTrackletsCluster(rofId, 1)) { - nTrackl12.push_back(n); - } - - ln_tre->Fill(); - } - trackletFile->cd(); - ln_tre->Write(); - trackletFile->Close(); + debugComputeTrackletMatching(iteration); #endif // from here on we do not use tracklets from L1-2 anymore, so let's free them @@ -418,9 +347,6 @@ void VertexerTraits::computeVertices(const int iteration) bounded_vector vertices(mMemoryPool.get()); bounded_vector> polls(mMemoryPool.get()); bounded_vector contLabels(mMemoryPool.get()); -#ifdef VTX_DEBUG - std::vector> dbg_clusLines(mTimeFrame->getNrof()); -#endif bounded_vector noClustersVec(mTimeFrame->getNrof(), 0, mMemoryPool.get()); for (int rofId{0}; rofId < mTimeFrame->getNrof(); ++rofId) { if (iteration && (int)mTimeFrame->getPrimaryVertices(rofId).size() > mVrtParams[iteration].vertPerRofThreshold) { @@ -501,12 +427,7 @@ void VertexerTraits::computeVertices(const int iteration) } for (int rofId{0}; rofId < mTimeFrame->getNrof(); ++rofId) { std::sort(mTimeFrame->getTrackletClusters(rofId).begin(), mTimeFrame->getTrackletClusters(rofId).end(), - [](ClusterLines& cluster1, ClusterLines& cluster2) { return cluster1.getSize() > cluster2.getSize(); }); // ensure clusters are ordered by contributors, so that we can cat after the first. -#ifdef VTX_DEBUG - for (auto& cl : mTimeFrame->getTrackletClusters(rofId)) { - dbg_clusLines[rofId].push_back(cl); - } -#endif + [](const ClusterLines& cluster1, const ClusterLines& cluster2) { return cluster1.getSize() > cluster2.getSize(); }); // ensure clusters are ordered by contributors, so that we can cat after the first. bool atLeastOneFound{false}; for (int iCluster{0}; iCluster < noClustersVec[rofId]; ++iCluster) { bool lowMultCandidate{false}; @@ -570,27 +491,9 @@ void VertexerTraits::computeVertices(const int iteration) vertices.clear(); polls.clear(); } + #ifdef VTX_DEBUG - TFile* dbg_file = TFile::Open("artefacts_tf.root", "update"); - TTree* ln_clus_lines_tree = new TTree("clusterlines", "tf"); - std::vector cl_lines_vec_pre(0); - std::vector cl_lines_vec_post(0); - ln_clus_lines_tree->Branch("cllines_pre", &cl_lines_vec_pre); - ln_clus_lines_tree->Branch("cllines_post", &cl_lines_vec_post); - for (auto rofId{0}; rofId < mTimeFrame->getNrof(); ++rofId) { - cl_lines_vec_pre.clear(); - cl_lines_vec_post.clear(); - for (auto& clln : mTimeFrame->getTrackletClusters(rofId)) { - cl_lines_vec_post.push_back(clln); - } - for (auto& cl : dbg_clusLines[rofId]) { - cl_lines_vec_pre.push_back(cl); - } - ln_clus_lines_tree->Fill(); - } - dbg_file->cd(); - ln_clus_lines_tree->Write(); - dbg_file->Close(); + debugComputeVertices(iteration); #endif } @@ -662,6 +565,7 @@ void VertexerTraits::addTruthSeedingVertices() void VertexerTraits::setNThreads(int n, std::shared_ptr& arena) { #if defined(VTX_DEBUG) + LOGP(info, "Vertexer with debug output forcing single thread"); mTaskArena = std::make_shared(1); #else if (arena == nullptr) { @@ -673,3 +577,257 @@ void VertexerTraits::setNThreads(int n, std::shared_ptr& arena) } #endif } + +void VertexerTraits::debugComputeTracklets(int iteration) +{ + auto stream = new utils::TreeStreamRedirector("artefacts_tf.root", "recreate"); + LOGP(info, "writing debug output for computeTracklets"); + for (int rofId{0}; rofId < mTimeFrame->getNrof(); ++rofId) { + const auto& strk0 = mTimeFrame->getFoundTracklets(rofId, 0); + std::vector trk0(strk0.begin(), strk0.end()); + const auto& strk1 = mTimeFrame->getFoundTracklets(rofId, 1); + std::vector trk1(strk1.begin(), strk1.end()); + (*stream) << "tracklets" + << "Tracklets0=" << trk0 + << "Tracklets1=" << trk1 + << "iteration=" << iteration + << "\n"; + } + stream->Close(); + delete stream; +} + +void VertexerTraits::debugComputeTrackletMatching(int iteration) +{ + auto stream = new utils::TreeStreamRedirector("artefacts_tf.root", "update"); + LOGP(info, "writing debug output for computeTrackletMatching"); + for (int rofId{0}; rofId < mTimeFrame->getNrof(); ++rofId) { + (*stream) << "lines" + << "Lines=" << toSTDVector(mTimeFrame->getLines(rofId)) + << "NTrackletCluster01=" << mTimeFrame->getNTrackletsCluster(rofId, 0) + << "NTrackletCluster12=" << mTimeFrame->getNTrackletsCluster(rofId, 1) + << "iteration=" << iteration + << "\n"; + } + + if (mTimeFrame->hasMCinformation()) { + LOGP(info, "\tdumping also MC information"); + const auto dc = o2::steer::DigitizationContext::loadFromFile("collisioncontext.root"); + const auto irs = dc->getEventRecords(); + int64_t roFrameBiasInBC = o2::itsmft::DPLAlpideParam::Instance().roFrameBiasInBC; + int64_t roFrameLengthInBC = o2::itsmft::DPLAlpideParam::Instance().roFrameLengthInBC; + o2::steer::MCKinematicsReader mcReader(dc); + + std::map eve2BcInROF, bcInRofNEve; + for (int iSrc{0}; iSrc < mcReader.getNSources(); ++iSrc) { + auto eveId2colId = dc->getCollisionIndicesForSource(iSrc); + for (int iEve{0}; iEve < mcReader.getNEvents(iSrc); ++iEve) { + const auto& ir = irs[eveId2colId[iEve]]; + if (!ir.isDummy()) { // do we need this, is this for diffractive events? + const auto& eve = mcReader.getMCEventHeader(iSrc, iEve); + const int bcInROF = ((ir - raw::HBFUtils::Instance().getFirstSampledTFIR()).toLong() - roFrameBiasInBC) % roFrameLengthInBC; + eve2BcInROF[iEve] = bcInROF; + ++bcInRofNEve[bcInROF]; + } + } + } + + std::unordered_map bcROFNTracklets01, bcROFNTracklets12; + std::vector> tracklet01BC, tracklet12BC; + for (int rofId{0}; rofId < mTimeFrame->getNrof(); ++rofId) { + { // 0-1 + const auto& tracklet01 = mTimeFrame->getFoundTracklets(rofId, 0); + const auto& lbls01 = mTimeFrame->getLabelsFoundTracklets(rofId, 0); + auto& trkls01 = tracklet01BC.emplace_back(); + for (int iTrklt{0}; iTrklt < (int)tracklet01.size(); ++iTrklt) { + const auto& tracklet = tracklet01[iTrklt]; + const auto& lbl = lbls01[iTrklt]; + if (lbl.isCorrect()) { + ++bcROFNTracklets01[eve2BcInROF[lbl.getEventID()]]; + trkls01.push_back(eve2BcInROF[lbl.getEventID()]); + } else { + trkls01.push_back(-1); + } + } + } + { // 1-2 computed on the fly! + const auto& tracklet12 = mTimeFrame->getFoundTracklets(rofId, 1); + auto& trkls12 = tracklet12BC.emplace_back(); + for (int iTrklt{0}; iTrklt < (int)tracklet12.size(); ++iTrklt) { + const auto& tracklet = tracklet12[iTrklt]; + o2::MCCompLabel label; + + int sortedId1{mTimeFrame->getSortedIndex(tracklet.rof[0], 1, tracklet.firstClusterIndex)}; + int sortedId2{mTimeFrame->getSortedIndex(tracklet.rof[1], 2, tracklet.secondClusterIndex)}; + for (const auto& lab1 : mTimeFrame->getClusterLabels(1, mTimeFrame->getClusters()[1][sortedId1].clusterId)) { + for (const auto& lab2 : mTimeFrame->getClusterLabels(2, mTimeFrame->getClusters()[2][sortedId2].clusterId)) { + if (lab1 == lab2 && lab1.isValid()) { + label = lab1; + break; + } + } + if (label.isValid()) { + break; + } + } + + if (label.isCorrect()) { + ++bcROFNTracklets12[eve2BcInROF[label.getEventID()]]; + trkls12.push_back(eve2BcInROF[label.getEventID()]); + } else { + trkls12.push_back(-1); + } + } + } + } + LOGP(info, "\tdumping ntracklets/RofBC ({})", bcInRofNEve.size()); + for (const auto& [bcInRof, neve] : bcInRofNEve) { + (*stream) << "ntracklets" + << "bcInROF=" << bcInRof + << "ntrkl01=" << bcROFNTracklets01[bcInRof] + << "ntrkl12=" << bcROFNTracklets12[bcInRof] + << "neve=" << neve + << "iteration=" << iteration + << "\n"; + } + + std::unordered_map bcROFNLines; + for (int rofId{0}; rofId < mTimeFrame->getNrof(); ++rofId) { + const auto& lines = mTimeFrame->getLines(rofId); + const auto& lbls = mTimeFrame->getLinesLabel(rofId); + for (int iLine{0}; iLine < (int)lines.size(); ++iLine) { + const auto& line = lines[iLine]; + const auto& lbl = lbls[iLine]; + if (lbl.isCorrect()) { + ++bcROFNLines[eve2BcInROF[lbl.getEventID()]]; + } + } + } + + LOGP(info, "\tdumping nlines/RofBC"); + for (const auto& [bcInRof, neve] : bcInRofNEve) { + (*stream) << "nlines" + << "bcInROF=" << bcInRof + << "nline=" << bcROFNLines[bcInRof] + << "neve=" << neve + << "iteration=" << iteration + << "\n"; + } + } + stream->Close(); + delete stream; +} + +void VertexerTraits::debugComputeVertices(int iteration) +{ + auto stream = new utils::TreeStreamRedirector("artefacts_tf.root", "update"); + LOGP(info, "writing debug output for computeVertices"); + for (auto rofId{0}; rofId < mTimeFrame->getNrof(); ++rofId) { + (*stream) << "clusterlines" + << "clines_post=" << toSTDVector(mTimeFrame->getTrackletClusters(rofId)) + << "iteration=" << iteration + << "\n"; + } + + if (mTimeFrame->hasMCinformation()) { + LOGP(info, "\tdumping also MC information"); + const auto dc = o2::steer::DigitizationContext::loadFromFile("collisioncontext.root"); + const auto irs = dc->getEventRecords(); + int64_t roFrameBiasInBC = o2::itsmft::DPLAlpideParam::Instance().roFrameBiasInBC; + int64_t roFrameLengthInBC = o2::itsmft::DPLAlpideParam::Instance().roFrameLengthInBC; + o2::steer::MCKinematicsReader mcReader(dc); + + std::map eve2BcInROF, bcInRofNEve; + for (int iSrc{0}; iSrc < mcReader.getNSources(); ++iSrc) { + auto eveId2colId = dc->getCollisionIndicesForSource(iSrc); + for (int iEve{0}; iEve < mcReader.getNEvents(iSrc); ++iEve) { + const auto& ir = irs[eveId2colId[iEve]]; + if (!ir.isDummy()) { // do we need this, is this for diffractive events? + const auto& eve = mcReader.getMCEventHeader(iSrc, iEve); + const int bcInROF = ((ir - raw::HBFUtils::Instance().getFirstSampledTFIR()).toLong() - roFrameBiasInBC) % roFrameLengthInBC; + eve2BcInROF[iEve] = bcInROF; + ++bcInRofNEve[bcInROF]; + } + } + } + + std::unordered_map bcROFNVtx; + std::unordered_map bcROFNPur; + std::unordered_map uniqueVertices; + for (int rofId{0}; rofId < mTimeFrame->getNrof(); ++rofId) { + const auto& pvs = mTimeFrame->getPrimaryVertices(rofId); + const auto& lblspv = mTimeFrame->getPrimaryVerticesMCRecInfo(rofId); + for (int i{0}; i < (int)pvs.size(); ++i) { + const auto& pv = pvs[i]; + const auto& [lbl, pur] = lblspv[i]; + if (lbl.isCorrect()) { + ++uniqueVertices[lbl]; + ++bcROFNVtx[eve2BcInROF[lbl.getEventID()]]; + bcROFNPur[eve2BcInROF[lbl.getEventID()]] += pur; + } + } + } + + std::unordered_map bcROFNUVtx, bcROFNCVtx; + for (const auto& [k, _] : eve2BcInROF) { + bcROFNUVtx[k] = bcROFNCVtx[k] = 0; + } + + for (const auto& [lbl, c] : uniqueVertices) { + if (c <= 1) { + ++bcROFNUVtx[eve2BcInROF[lbl.getEventID()]]; + } else { + ++bcROFNCVtx[eve2BcInROF[lbl.getEventID()]]; + } + } + + LOGP(info, "\tdumping nvtx/RofBC"); + for (const auto& [bcInRof, neve] : bcInRofNEve) { + (*stream) << "nvtx" + << "bcInROF=" << bcInRof + << "nvtx=" << bcROFNVtx[bcInRof] // all vertices + << "nuvtx=" << bcROFNUVtx[bcInRof] // unique vertices + << "ncvtx=" << bcROFNCVtx[bcInRof] // cloned vertices + << "npur=" << bcROFNPur[bcInRof] + << "neve=" << neve + << "iteration=" << iteration + << "\n"; + } + + // check dist of clones + std::unordered_map> cVtx; + for (int rofId{0}; rofId < mTimeFrame->getNrof(); ++rofId) { + const auto& pvs = mTimeFrame->getPrimaryVertices(rofId); + const auto& lblspv = mTimeFrame->getPrimaryVerticesMCRecInfo(rofId); + for (int i{0}; i < (int)pvs.size(); ++i) { + const auto& pv = pvs[i]; + const auto& [lbl, pur] = lblspv[i]; + if (lbl.isCorrect() && uniqueVertices.contains(lbl) && uniqueVertices[lbl] > 1) { + if (!cVtx.contains(lbl)) { + cVtx[lbl] = std::vector(); + } + cVtx[lbl].push_back(pv); + } + } + } + + for (auto& [_, vertices] : cVtx) { + std::sort(vertices.begin(), vertices.end(), [](const Vertex& a, const Vertex& b) { return a.getNContributors() > b.getNContributors(); }); + for (int i{0}; i < (int)vertices.size(); ++i) { + const auto vtx = vertices[i]; + (*stream) << "cvtx" + << "vertex=" << vtx + << "i=" << i + << "dx=" << vertices[0].getX() - vtx.getX() + << "dy=" << vertices[0].getY() - vtx.getY() + << "dz=" << vertices[0].getZ() - vtx.getZ() + << "drof=" << vertices[0].getTimeStamp().getTimeStamp() - vtx.getTimeStamp().getTimeStamp() + << "dnc=" << vertices[0].getNContributors() - vtx.getNContributors() + << "iteration=" << iteration + << "\n"; + } + } + } + stream->Close(); + delete stream; +}