Project
Loading...
Searching...
No Matches
VertexerTraits.cxx
Go to the documentation of this file.
1// Copyright 2019-2020 CERN and copyright holders of ALICE O2.
2// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
3// All rights not expressly granted are reserved.
4//
5// This software is distributed under the terms of the GNU General Public
6// License v3 (GPL Version 3), copied verbatim in the file "COPYING".
7//
8// In applying this license CERN does not waive the privileges and immunities
9// granted to it by virtue of its status as an Intergovernmental Organization
10// or submit itself to any jurisdiction.
12
13#include <memory>
14#include <ranges>
15#include <map>
16#include <algorithm>
17
18#include <oneapi/tbb/blocked_range.h>
19#include <oneapi/tbb/parallel_for.h>
20#include <oneapi/tbb/combinable.h>
21
31
32namespace o2::its
33{
34
35template <TrackletMode Mode, bool EvalRun, int nLayers>
36static void trackleterKernelHost(
37 const gsl::span<const Cluster>& clustersNextLayer, // 0 2
38 const gsl::span<const Cluster>& clustersCurrentLayer, // 1 1
39 const gsl::span<uint8_t>& usedClustersNextLayer, // 0 2
40 int* indexTableNext,
41 const float phiCut,
42 bounded_vector<Tracklet>& tracklets,
43 gsl::span<int> foundTracklets,
44 const IndexTableUtils<nLayers>& utils,
45 const short pivotRof,
46 const short targetRof,
47 gsl::span<int> rofFoundTrackletsOffsets, // we want to change those, to keep track of the offset in deltaRof>0
48 const int maxTrackletsPerCluster = static_cast<int>(2e3))
49{
50 const int PhiBins{utils.getNphiBins()};
51 const int ZBins{utils.getNzBins()};
52 // loop on layer1 clusters
53 for (int iCurrentLayerClusterIndex = 0; iCurrentLayerClusterIndex < clustersCurrentLayer.size(); ++iCurrentLayerClusterIndex) {
54 int storedTracklets{0};
55 const Cluster& currentCluster{clustersCurrentLayer[iCurrentLayerClusterIndex]};
56 const int4 selectedBinsRect{VertexerTraits<nLayers>::getBinsRect(currentCluster, (int)Mode, 0.f, 50.f, phiCut / 2, utils)};
57 if (selectedBinsRect.x != 0 || selectedBinsRect.y != 0 || selectedBinsRect.z != 0 || selectedBinsRect.w != 0) {
58 int phiBinsNum{selectedBinsRect.w - selectedBinsRect.y + 1};
59 if (phiBinsNum < 0) {
60 phiBinsNum += PhiBins;
61 }
62 // loop on phi bins next layer
63 for (int iPhiBin{selectedBinsRect.y}, iPhiCount{0}; iPhiCount < phiBinsNum; iPhiBin = ++iPhiBin == PhiBins ? 0 : iPhiBin, iPhiCount++) {
64 const int firstBinIndex{utils.getBinIndex(selectedBinsRect.x, iPhiBin)};
65 const int firstRowClusterIndex{indexTableNext[firstBinIndex]};
66 const int maxRowClusterIndex{indexTableNext[firstBinIndex + ZBins]};
67 // loop on clusters next layer
68 for (int iNextLayerClusterIndex{firstRowClusterIndex}; iNextLayerClusterIndex < maxRowClusterIndex && iNextLayerClusterIndex < static_cast<int>(clustersNextLayer.size()); ++iNextLayerClusterIndex) {
69 if (usedClustersNextLayer[iNextLayerClusterIndex]) {
70 continue;
71 }
72 const Cluster& nextCluster{clustersNextLayer[iNextLayerClusterIndex]};
73 if (o2::gpu::GPUCommonMath::Abs(math_utils::smallestAngleDifference(currentCluster.phi, nextCluster.phi)) < phiCut) {
74 if (storedTracklets < maxTrackletsPerCluster) {
75 if constexpr (!EvalRun) {
76 if constexpr (Mode == TrackletMode::Layer0Layer1) {
77 tracklets[rofFoundTrackletsOffsets[iCurrentLayerClusterIndex] + storedTracklets] = Tracklet{iNextLayerClusterIndex, iCurrentLayerClusterIndex, nextCluster, currentCluster, targetRof, pivotRof};
78 } else {
79 tracklets[rofFoundTrackletsOffsets[iCurrentLayerClusterIndex] + storedTracklets] = Tracklet{iCurrentLayerClusterIndex, iNextLayerClusterIndex, currentCluster, nextCluster, pivotRof, targetRof};
80 }
81 }
82 ++storedTracklets;
83 }
84 }
85 }
86 }
87 }
88 if constexpr (EvalRun) {
89 foundTracklets[iCurrentLayerClusterIndex] += storedTracklets;
90 } else {
91 rofFoundTrackletsOffsets[iCurrentLayerClusterIndex] += storedTracklets;
92 }
93 }
94}
95
96static void trackletSelectionKernelHost(
97 const gsl::span<const Cluster> clusters0, // 0
98 const gsl::span<const Cluster> clusters1, // 1
99 gsl::span<unsigned char> usedClusters0, // Layer 0
100 gsl::span<unsigned char> usedClusters2, // Layer 2
101 const gsl::span<const Tracklet>& tracklets01,
102 const gsl::span<const Tracklet>& tracklets12,
103 bounded_vector<bool>& usedTracklets,
104 const gsl::span<int> foundTracklets01,
105 const gsl::span<int> foundTracklets12,
106 bounded_vector<Line>& lines,
107 const gsl::span<const o2::MCCompLabel>& trackletLabels,
108 bounded_vector<o2::MCCompLabel>& linesLabels,
109 const short targetRofId0,
110 const short targetRofId2,
111 bool safeWrites = false,
112 const float tanLambdaCut = 0.025f,
113 const float phiCut = 0.005f,
114 const int maxTracklets = static_cast<int>(1e2))
115{
116 int offset01{0}, offset12{0};
117 for (unsigned int iCurrentLayerClusterIndex{0}; iCurrentLayerClusterIndex < clusters1.size(); ++iCurrentLayerClusterIndex) {
118 int validTracklets{0};
119 for (int iTracklet12{offset12}; iTracklet12 < offset12 + foundTracklets12[iCurrentLayerClusterIndex]; ++iTracklet12) {
120 for (int iTracklet01{offset01}; iTracklet01 < offset01 + foundTracklets01[iCurrentLayerClusterIndex]; ++iTracklet01) {
121 if (usedTracklets[iTracklet01]) {
122 continue;
123 }
124
125 const auto& tracklet01{tracklets01[iTracklet01]};
126 const auto& tracklet12{tracklets12[iTracklet12]};
127
128 if (tracklet01.rof[0] != targetRofId0 || tracklet12.rof[1] != targetRofId2) {
129 continue;
130 }
131
132 const float deltaTanLambda{o2::gpu::GPUCommonMath::Abs(tracklet01.tanLambda - tracklet12.tanLambda)};
133 const float deltaPhi{o2::gpu::GPUCommonMath::Abs(math_utils::smallestAngleDifference(tracklet01.phi, tracklet12.phi))};
134 if (deltaTanLambda < tanLambdaCut && deltaPhi < phiCut && validTracklets != maxTracklets) {
135 if (safeWrites) {
136 __atomic_store_n(&usedClusters0[tracklet01.firstClusterIndex], 1, __ATOMIC_RELAXED);
137 __atomic_store_n(&usedClusters2[tracklet12.secondClusterIndex], 1, __ATOMIC_RELAXED);
138 } else {
139 usedClusters0[tracklet01.firstClusterIndex] = 1;
140 usedClusters2[tracklet12.secondClusterIndex] = 1;
141 }
142 usedTracklets[iTracklet01] = true;
143 lines.emplace_back(tracklet01, clusters0.data(), clusters1.data());
144 if (!trackletLabels.empty()) {
145 linesLabels.emplace_back(trackletLabels[iTracklet01]);
146 }
147 ++validTracklets;
148 }
149 }
150 }
151 offset01 += foundTracklets01[iCurrentLayerClusterIndex];
152 offset12 += foundTracklets12[iCurrentLayerClusterIndex];
153 }
154}
155
156template <int nLayers>
157void VertexerTraits<nLayers>::updateVertexingParameters(const std::vector<VertexingParameters>& vrtPar, const TimeFrameGPUParameters& tfPar)
158{
159 mVrtParams = vrtPar;
160 mIndexTableUtils.setTrackingParameters(vrtPar[0]);
161 for (auto& par : mVrtParams) {
162 par.phiSpan = static_cast<int>(std::ceil(mIndexTableUtils.getNphiBins() * par.phiCut / o2::constants::math::TwoPI));
163 par.zSpan = static_cast<int>(std::ceil(par.zCut * mIndexTableUtils.getInverseZCoordinate(0)));
164 }
165}
166
167// Main functions
168template <int nLayers>
170{
171 mTaskArena->execute([&] {
172 tbb::parallel_for(0, mTimeFrame->getNrof(), [&](const short pivotRofId) {
173 bool skipROF = iteration && (int)mTimeFrame->getPrimaryVertices(pivotRofId).size() > mVrtParams[iteration].vertPerRofThreshold;
174 short startROF{std::max((short)0, static_cast<short>(pivotRofId - mVrtParams[iteration].deltaRof))};
175 short endROF{std::min(static_cast<short>(mTimeFrame->getNrof()), static_cast<short>(pivotRofId + mVrtParams[iteration].deltaRof + 1))};
176 for (auto targetRofId = startROF; targetRofId < endROF; ++targetRofId) {
177 trackleterKernelHost<TrackletMode::Layer0Layer1, true>(
178 !skipROF ? mTimeFrame->getClustersOnLayer(targetRofId, 0) : gsl::span<Cluster>(), // Clusters to be matched with the next layer in target rof
179 !skipROF ? mTimeFrame->getClustersOnLayer(pivotRofId, 1) : gsl::span<Cluster>(), // Clusters to be matched with the current layer in pivot rof
180 mTimeFrame->getUsedClustersROF(targetRofId, 0), // Span of the used clusters in the target rof
181 mTimeFrame->getIndexTable(targetRofId, 0).data(), // Index table to access the data on the next layer in target rof
182 mVrtParams[iteration].phiCut,
183 mTimeFrame->getTracklets()[0], // Flat tracklet buffer
184 mTimeFrame->getNTrackletsCluster(pivotRofId, 0), // Span of the number of tracklets per each cluster in pivot rof
185 mIndexTableUtils,
186 pivotRofId,
187 targetRofId,
188 gsl::span<int>(), // Offset in the tracklet buffer
189 mVrtParams[iteration].maxTrackletsPerCluster);
190 trackleterKernelHost<TrackletMode::Layer1Layer2, true>(
191 !skipROF ? mTimeFrame->getClustersOnLayer(targetRofId, 2) : gsl::span<Cluster>(),
192 !skipROF ? mTimeFrame->getClustersOnLayer(pivotRofId, 1) : gsl::span<Cluster>(),
193 mTimeFrame->getUsedClustersROF(targetRofId, 2),
194 mTimeFrame->getIndexTable(targetRofId, 2).data(),
195 mVrtParams[iteration].phiCut,
196 mTimeFrame->getTracklets()[1],
197 mTimeFrame->getNTrackletsCluster(pivotRofId, 1), // Span of the number of tracklets per each cluster in pivot rof
198 mIndexTableUtils,
199 pivotRofId,
200 targetRofId,
201 gsl::span<int>(), // Offset in the tracklet buffer
202 mVrtParams[iteration].maxTrackletsPerCluster);
203 }
204 mTimeFrame->getNTrackletsROF(pivotRofId, 0) = std::accumulate(mTimeFrame->getNTrackletsCluster(pivotRofId, 0).begin(), mTimeFrame->getNTrackletsCluster(pivotRofId, 0).end(), 0);
205 mTimeFrame->getNTrackletsROF(pivotRofId, 1) = std::accumulate(mTimeFrame->getNTrackletsCluster(pivotRofId, 1).begin(), mTimeFrame->getNTrackletsCluster(pivotRofId, 1).end(), 0);
206 });
207
208 mTimeFrame->computeTrackletsPerROFScans();
209 if (auto tot0 = mTimeFrame->getTotalTrackletsTF(0), tot1 = mTimeFrame->getTotalTrackletsTF(1);
210 tot0 == 0 || tot1 == 0) {
211 return;
212 } else {
213 mTimeFrame->getTracklets()[0].resize(tot0);
214 mTimeFrame->getTracklets()[1].resize(tot1);
215 }
216
217 tbb::parallel_for(0, mTimeFrame->getNrof(), [&](const short pivotRofId) {
218 bool skipROF = iteration && (int)mTimeFrame->getPrimaryVertices(pivotRofId).size() > mVrtParams[iteration].vertPerRofThreshold;
219 short startROF{std::max((short)0, static_cast<short>(pivotRofId - mVrtParams[iteration].deltaRof))};
220 short endROF{std::min(static_cast<short>(mTimeFrame->getNrof()), static_cast<short>(pivotRofId + mVrtParams[iteration].deltaRof + 1))};
221 auto mobileOffset0 = mTimeFrame->getNTrackletsROF(pivotRofId, 0);
222 auto mobileOffset1 = mTimeFrame->getNTrackletsROF(pivotRofId, 1);
223 for (auto targetRofId = startROF; targetRofId < endROF; ++targetRofId) {
224 trackleterKernelHost<TrackletMode::Layer0Layer1, false>(
225 !skipROF ? mTimeFrame->getClustersOnLayer(targetRofId, 0) : gsl::span<Cluster>(),
226 !skipROF ? mTimeFrame->getClustersOnLayer(pivotRofId, 1) : gsl::span<Cluster>(),
227 mTimeFrame->getUsedClustersROF(targetRofId, 0),
228 mTimeFrame->getIndexTable(targetRofId, 0).data(),
229 mVrtParams[iteration].phiCut,
230 mTimeFrame->getTracklets()[0],
231 mTimeFrame->getNTrackletsCluster(pivotRofId, 0),
232 mIndexTableUtils,
233 pivotRofId,
234 targetRofId,
235 mTimeFrame->getExclusiveNTrackletsCluster(pivotRofId, 0),
236 mVrtParams[iteration].maxTrackletsPerCluster);
237 trackleterKernelHost<TrackletMode::Layer1Layer2, false>(
238 !skipROF ? mTimeFrame->getClustersOnLayer(targetRofId, 2) : gsl::span<Cluster>(),
239 !skipROF ? mTimeFrame->getClustersOnLayer(pivotRofId, 1) : gsl::span<Cluster>(),
240 mTimeFrame->getUsedClustersROF(targetRofId, 2),
241 mTimeFrame->getIndexTable(targetRofId, 2).data(),
242 mVrtParams[iteration].phiCut,
243 mTimeFrame->getTracklets()[1],
244 mTimeFrame->getNTrackletsCluster(pivotRofId, 1),
245 mIndexTableUtils,
246 pivotRofId,
247 targetRofId,
248 mTimeFrame->getExclusiveNTrackletsCluster(pivotRofId, 1),
249 mVrtParams[iteration].maxTrackletsPerCluster);
250 }
251 });
252 });
253
255 if (mTimeFrame->hasMCinformation()) {
256 for (const auto& trk : mTimeFrame->getTracklets()[0]) {
258 if (!trk.isEmpty()) {
259 int sortedId0{mTimeFrame->getSortedIndex(trk.rof[0], 0, trk.firstClusterIndex)};
260 int sortedId1{mTimeFrame->getSortedIndex(trk.rof[1], 1, trk.secondClusterIndex)};
261 for (const auto& lab0 : mTimeFrame->getClusterLabels(0, mTimeFrame->getClusters()[0][sortedId0].clusterId)) {
262 for (const auto& lab1 : mTimeFrame->getClusterLabels(1, mTimeFrame->getClusters()[1][sortedId1].clusterId)) {
263 if (lab0 == lab1 && lab0.isValid()) {
264 label = lab0;
265 break;
266 }
267 }
268 if (label.isValid()) {
269 break;
270 }
271 }
272 }
273 mTimeFrame->getTrackletsLabel(0).emplace_back(label);
274 }
275 }
276
277#ifdef VTX_DEBUG
278 debugComputeTracklets(iteration);
279#endif
280}
281
282template <int nLayers>
284{
285 mTaskArena->execute([&] {
286 tbb::combinable<int> totalLines{0};
287 tbb::parallel_for(
288 tbb::blocked_range<short>(0, (short)mTimeFrame->getNrof()),
289 [&](const tbb::blocked_range<short>& Rofs) {
290 for (short pivotRofId = Rofs.begin(); pivotRofId < Rofs.end(); ++pivotRofId) {
291 if (iteration && (int)mTimeFrame->getPrimaryVertices(pivotRofId).size() > mVrtParams[iteration].vertPerRofThreshold) {
292 continue;
293 }
294 if (mTimeFrame->getFoundTracklets(pivotRofId, 0).empty()) {
295 continue;
296 }
297 mTimeFrame->getLines(pivotRofId).reserve(mTimeFrame->getNTrackletsCluster(pivotRofId, 0).size());
298 bounded_vector<bool> usedTracklets(mTimeFrame->getFoundTracklets(pivotRofId, 0).size(), false, mMemoryPool.get());
299 short startROF{std::max((short)0, static_cast<short>(pivotRofId - mVrtParams[iteration].deltaRof))};
300 short endROF{std::min(static_cast<short>(mTimeFrame->getNrof()), static_cast<short>(pivotRofId + mVrtParams[iteration].deltaRof + 1))};
301
302 // needed only if multi-threaded using deltaRof and only at the overlap edges of the ranges
303 bool safeWrite = mTaskArena->max_concurrency() > 1 && mVrtParams[iteration].deltaRof != 0 && ((Rofs.begin() - startROF < 0) || (endROF - Rofs.end() > 0));
304
305 for (short targetRofId0 = startROF; targetRofId0 < endROF; ++targetRofId0) {
306 for (short targetRofId2 = startROF; targetRofId2 < endROF; ++targetRofId2) {
307 if (std::abs(targetRofId0 - targetRofId2) > mVrtParams[iteration].deltaRof) { // do not allow over 3 ROFs
308 continue;
309 }
310 trackletSelectionKernelHost(
311 mTimeFrame->getClustersOnLayer(targetRofId0, 0),
312 mTimeFrame->getClustersOnLayer(pivotRofId, 1),
313 mTimeFrame->getUsedClustersROF(targetRofId0, 0),
314 mTimeFrame->getUsedClustersROF(targetRofId2, 2),
315 mTimeFrame->getFoundTracklets(pivotRofId, 0),
316 mTimeFrame->getFoundTracklets(pivotRofId, 1),
317 usedTracklets,
318 mTimeFrame->getNTrackletsCluster(pivotRofId, 0),
319 mTimeFrame->getNTrackletsCluster(pivotRofId, 1),
320 mTimeFrame->getLines(pivotRofId),
321 mTimeFrame->getLabelsFoundTracklets(pivotRofId, 0),
322 mTimeFrame->getLinesLabel(pivotRofId),
323 targetRofId0,
324 targetRofId2,
325 safeWrite,
326 mVrtParams[iteration].tanLambdaCut,
327 mVrtParams[iteration].phiCut);
328 }
329 }
330 totalLines.local() += mTimeFrame->getLines(pivotRofId).size();
331 }
332 });
333 mTimeFrame->setNLinesTotal(totalLines.combine(std::plus<int>()));
334 });
335
336#ifdef VTX_DEBUG
337 debugComputeTrackletMatching(iteration);
338#endif
339
340 // from here on we do not use tracklets from L1-2 anymore, so let's free them
341 deepVectorClear(mTimeFrame->getTracklets()[1]);
342}
343
344template <int nLayers>
346{
347 auto nsigmaCut{std::min(mVrtParams[iteration].vertNsigmaCut * mVrtParams[iteration].vertNsigmaCut * (mVrtParams[iteration].vertRadiusSigma * mVrtParams[iteration].vertRadiusSigma + mVrtParams[iteration].trackletSigma * mVrtParams[iteration].trackletSigma), 1.98f)};
348 bounded_vector<Vertex> vertices(mMemoryPool.get());
349 bounded_vector<std::pair<o2::MCCompLabel, float>> polls(mMemoryPool.get());
350 bounded_vector<o2::MCCompLabel> contLabels(mMemoryPool.get());
351 bounded_vector<int> noClustersVec(mTimeFrame->getNrof(), 0, mMemoryPool.get());
352 for (int rofId{0}; rofId < mTimeFrame->getNrof(); ++rofId) {
353 if (iteration && (int)mTimeFrame->getPrimaryVertices(rofId).size() > mVrtParams[iteration].vertPerRofThreshold) {
354 continue;
355 }
356 const int numTracklets{static_cast<int>(mTimeFrame->getLines(rofId).size())};
357
358 bounded_vector<bool> usedTracklets(numTracklets, false, mMemoryPool.get());
359 for (int line1{0}; line1 < numTracklets; ++line1) {
360 if (usedTracklets[line1]) {
361 continue;
362 }
363 for (int line2{line1 + 1}; line2 < numTracklets; ++line2) {
364 if (usedTracklets[line2]) {
365 continue;
366 }
367 auto dca{Line::getDCA(mTimeFrame->getLines(rofId)[line1], mTimeFrame->getLines(rofId)[line2])};
368 if (dca < mVrtParams[iteration].pairCut) {
369 mTimeFrame->getTrackletClusters(rofId).emplace_back(line1, mTimeFrame->getLines(rofId)[line1], line2, mTimeFrame->getLines(rofId)[line2]);
370 std::array<float, 3> tmpVertex{mTimeFrame->getTrackletClusters(rofId).back().getVertex()};
371 if (tmpVertex[0] * tmpVertex[0] + tmpVertex[1] * tmpVertex[1] > 4.f) {
372 mTimeFrame->getTrackletClusters(rofId).pop_back();
373 break;
374 }
375 usedTracklets[line1] = true;
376 usedTracklets[line2] = true;
377 for (int tracklet3{0}; tracklet3 < numTracklets; ++tracklet3) {
378 if (usedTracklets[tracklet3]) {
379 continue;
380 }
381 if (Line::getDistanceFromPoint(mTimeFrame->getLines(rofId)[tracklet3], tmpVertex) < mVrtParams[iteration].pairCut) {
382 mTimeFrame->getTrackletClusters(rofId).back().add(tracklet3, mTimeFrame->getLines(rofId)[tracklet3]);
383 usedTracklets[tracklet3] = true;
384 tmpVertex = mTimeFrame->getTrackletClusters(rofId).back().getVertex();
385 }
386 }
387 break;
388 }
389 }
390 }
391 if (mVrtParams[iteration].allowSingleContribClusters) {
392 auto beamLine = Line{{mTimeFrame->getBeamX(), mTimeFrame->getBeamY(), -50.f}, {mTimeFrame->getBeamX(), mTimeFrame->getBeamY(), 50.f}}; // use beam position as contributor
393 for (size_t iLine{0}; iLine < numTracklets; ++iLine) {
394 if (!usedTracklets[iLine]) {
395 auto dca = Line::getDCA(mTimeFrame->getLines(rofId)[iLine], beamLine);
396 if (dca < mVrtParams[iteration].pairCut) {
397 mTimeFrame->getTrackletClusters(rofId).emplace_back(iLine, mTimeFrame->getLines(rofId)[iLine], -1, beamLine); // beamline must be passed as second line argument
398 }
399 }
400 }
401 }
402
403 // Cluster merging
404 std::sort(mTimeFrame->getTrackletClusters(rofId).begin(), mTimeFrame->getTrackletClusters(rofId).end(),
405 [](ClusterLines& cluster1, ClusterLines& cluster2) { return cluster1.getSize() > cluster2.getSize(); });
406 noClustersVec[rofId] = static_cast<int>(mTimeFrame->getTrackletClusters(rofId).size());
407 for (int iCluster1{0}; iCluster1 < noClustersVec[rofId]; ++iCluster1) {
408 std::array<float, 3> vertex1{mTimeFrame->getTrackletClusters(rofId)[iCluster1].getVertex()};
409 std::array<float, 3> vertex2{};
410 for (int iCluster2{iCluster1 + 1}; iCluster2 < noClustersVec[rofId]; ++iCluster2) {
411 vertex2 = mTimeFrame->getTrackletClusters(rofId)[iCluster2].getVertex();
412 if (o2::gpu::GPUCommonMath::Abs(vertex1[2] - vertex2[2]) < mVrtParams[iteration].clusterCut) {
413 float distance{(vertex1[0] - vertex2[0]) * (vertex1[0] - vertex2[0]) +
414 (vertex1[1] - vertex2[1]) * (vertex1[1] - vertex2[1]) +
415 (vertex1[2] - vertex2[2]) * (vertex1[2] - vertex2[2])};
416 if (distance < mVrtParams[iteration].pairCut * mVrtParams[iteration].pairCut) {
417 for (auto label : mTimeFrame->getTrackletClusters(rofId)[iCluster2].getLabels()) {
418 mTimeFrame->getTrackletClusters(rofId)[iCluster1].add(label, mTimeFrame->getLines(rofId)[label]);
419 vertex1 = mTimeFrame->getTrackletClusters(rofId)[iCluster1].getVertex();
420 }
421 mTimeFrame->getTrackletClusters(rofId).erase(mTimeFrame->getTrackletClusters(rofId).begin() + iCluster2);
422 --iCluster2;
423 --noClustersVec[rofId];
424 }
425 }
426 }
427 }
428 }
429 for (int rofId{0}; rofId < mTimeFrame->getNrof(); ++rofId) {
430 std::sort(mTimeFrame->getTrackletClusters(rofId).begin(), mTimeFrame->getTrackletClusters(rofId).end(),
431 [](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.
432 bool atLeastOneFound{false};
433 for (int iCluster{0}; iCluster < noClustersVec[rofId]; ++iCluster) {
434 bool lowMultCandidate{false};
435 double beamDistance2{(mTimeFrame->getBeamX() - mTimeFrame->getTrackletClusters(rofId)[iCluster].getVertex()[0]) * (mTimeFrame->getBeamX() - mTimeFrame->getTrackletClusters(rofId)[iCluster].getVertex()[0]) +
436 (mTimeFrame->getBeamY() - mTimeFrame->getTrackletClusters(rofId)[iCluster].getVertex()[1]) * (mTimeFrame->getBeamY() - mTimeFrame->getTrackletClusters(rofId)[iCluster].getVertex()[1])};
437 if (atLeastOneFound && (lowMultCandidate = mTimeFrame->getTrackletClusters(rofId)[iCluster].getSize() < mVrtParams[iteration].clusterContributorsCut)) { // We might have pile up with nContr > cut.
438 lowMultCandidate &= (beamDistance2 < mVrtParams[iteration].lowMultBeamDistCut * mVrtParams[iteration].lowMultBeamDistCut);
439 if (!lowMultCandidate) { // Not the first cluster and not a low multiplicity candidate, we can remove it
440 mTimeFrame->getTrackletClusters(rofId).erase(mTimeFrame->getTrackletClusters(rofId).begin() + iCluster);
441 noClustersVec[rofId]--;
442 continue;
443 }
444 }
445
446 if (beamDistance2 < nsigmaCut && o2::gpu::GPUCommonMath::Abs(mTimeFrame->getTrackletClusters(rofId)[iCluster].getVertex()[2]) < mVrtParams[iteration].maxZPositionAllowed) {
447 atLeastOneFound = true;
448 auto& vertex = vertices.emplace_back(o2::math_utils::Point3D<float>(mTimeFrame->getTrackletClusters(rofId)[iCluster].getVertex()[0],
449 mTimeFrame->getTrackletClusters(rofId)[iCluster].getVertex()[1],
450 mTimeFrame->getTrackletClusters(rofId)[iCluster].getVertex()[2]),
451 mTimeFrame->getTrackletClusters(rofId)[iCluster].getRMS2(), // Symm matrix. Diagonal: RMS2 components,
452 // off-diagonal: square mean of projections on planes.
453 mTimeFrame->getTrackletClusters(rofId)[iCluster].getSize(), // Contributors
454 mTimeFrame->getTrackletClusters(rofId)[iCluster].getAvgDistance2()); // In place of chi2
455
456 if (iteration) {
457 vertex.setFlags(Vertex::UPCMode);
458 }
459 vertex.setTimeStamp(mTimeFrame->getTrackletClusters(rofId)[iCluster].getROF());
460 if (mTimeFrame->hasMCinformation()) {
461 bounded_vector<o2::MCCompLabel> labels(mMemoryPool.get());
462 for (auto& index : mTimeFrame->getTrackletClusters(rofId)[iCluster].getLabels()) {
463 labels.push_back(mTimeFrame->getLinesLabel(rofId)[index]); // then we can use nContributors from vertices to get the labels
464 }
465 polls.push_back(computeMain(labels));
466 if (mVrtParams[iteration].outputContLabels) {
467 contLabels.insert(contLabels.end(), labels.begin(), labels.end());
468 }
469 }
470 }
471 }
472 if (!iteration) {
473 mTimeFrame->addPrimaryVertices(vertices, iteration);
474 if (mTimeFrame->hasMCinformation()) {
475 mTimeFrame->addPrimaryVerticesLabels(polls);
476 if (mVrtParams[iteration].outputContLabels) {
477 mTimeFrame->addPrimaryVerticesContributorLabels(contLabels);
478 }
479 }
480 } else {
481 mTimeFrame->addPrimaryVerticesInROF(vertices, rofId, iteration);
482 if (mTimeFrame->hasMCinformation()) {
483 mTimeFrame->addPrimaryVerticesLabelsInROF(polls, rofId);
484 if (mVrtParams[iteration].outputContLabels) {
485 mTimeFrame->addPrimaryVerticesContributorLabelsInROF(contLabels, rofId);
486 }
487 }
488 }
489 if (vertices.empty() && !(iteration && (int)mTimeFrame->getPrimaryVertices(rofId).size() > mVrtParams[iteration].vertPerRofThreshold)) {
490 mTimeFrame->getNoVertexROF()++;
491 }
492 vertices.clear();
493 polls.clear();
494 }
495
496#ifdef VTX_DEBUG
497 debugComputeVertices(iteration);
498#endif
499}
500
501template <int nLayers>
503{
504 LOGP(info, "Using truth seeds as vertices; will skip computations");
505 mTimeFrame->resetRofPV();
506 const auto dc = o2::steer::DigitizationContext::loadFromFile("collisioncontext.root");
507 const auto irs = dc->getEventRecords();
508 int64_t roFrameBiasInBC = o2::itsmft::DPLAlpideParam<o2::detectors::DetID::ITS>::Instance().roFrameBiasInBC;
509 int64_t roFrameLengthInBC = o2::itsmft::DPLAlpideParam<o2::detectors::DetID::ITS>::Instance().roFrameLengthInBC;
511 struct VertInfo {
512 bounded_vector<Vertex> vertices;
514 bounded_vector<int> events;
515 };
516 std::map<int, VertInfo> vertices;
517 const int iSrc = 0; // take only events from collision generator
518 auto eveId2colId = dc->getCollisionIndicesForSource(iSrc);
519 for (int iEve{0}; iEve < mcReader.getNEvents(iSrc); ++iEve) {
520 const auto& ir = irs[eveId2colId[iEve]];
521 if (!ir.isDummy()) { // do we need this, is this for diffractive events?
522 const auto& eve = mcReader.getMCEventHeader(iSrc, iEve);
523 int rofId = ((ir - raw::HBFUtils::Instance().getFirstSampledTFIR()).toLong() - roFrameBiasInBC) / roFrameLengthInBC;
524 if (!vertices.contains(rofId)) {
525 vertices[rofId] = {
526 .vertices = bounded_vector<Vertex>(mMemoryPool.get()),
527 .srcs = bounded_vector<int>(mMemoryPool.get()),
528 .events = bounded_vector<int>(mMemoryPool.get()),
529 };
530 }
531 Vertex vert;
532 vert.setTimeStamp(rofId);
533 // set minimum to 1 sometimes for diffractive events there is nothing acceptance
534 vert.setNContributors(std::max(1L, std::ranges::count_if(mcReader.getTracks(iSrc, iEve), [](const auto& trk) {
535 return trk.isPrimary() && trk.GetPt() > 0.05 && std::abs(trk.GetEta()) < 1.1;
536 })));
537 vert.setXYZ((float)eve.GetX(), (float)eve.GetY(), (float)eve.GetZ());
538 vert.setChi2(1); // not used as constraint
539 constexpr float cov = 50e-9;
540 vert.setCov(cov, cov, cov, cov, cov, cov);
541 vertices[rofId].vertices.push_back(vert);
542 vertices[rofId].srcs.push_back(iSrc);
543 vertices[rofId].events.push_back(iEve);
544 }
545 mcReader.releaseTracksForSourceAndEvent(iSrc, iEve);
546 }
547 size_t nVerts{0};
548 for (int iROF{0}; iROF < mTimeFrame->getNrof(); ++iROF) {
549 bounded_vector<Vertex> verts(mMemoryPool.get());
550 bounded_vector<std::pair<o2::MCCompLabel, float>> polls(mMemoryPool.get());
551 if (vertices.contains(iROF)) {
552 const auto& vertInfo = vertices[iROF];
553 verts = vertInfo.vertices;
554 nVerts += verts.size();
555 for (size_t i{0}; i < verts.size(); ++i) {
556 o2::MCCompLabel lbl(o2::MCCompLabel::maxTrackID(), vertInfo.events[i], vertInfo.srcs[i], false);
557 polls.emplace_back(lbl, 1.f);
558 }
559 } else {
560 mTimeFrame->getNoVertexROF()++;
561 }
562 mTimeFrame->addPrimaryVertices(verts, 0);
563 mTimeFrame->addPrimaryVerticesLabels(polls);
564 }
565 LOGP(info, "Found {}/{} ROFs with {} vertices -> <NV>={:.2f}", vertices.size(), mTimeFrame->getNrof(), nVerts, (float)nVerts / (float)vertices.size());
566}
567
568template <int nLayers>
569void VertexerTraits<nLayers>::setNThreads(int n, std::shared_ptr<tbb::task_arena>& arena)
570{
571#if defined(VTX_DEBUG)
572 LOGP(info, "Vertexer with debug output forcing single thread");
573 mTaskArena = std::make_shared<tbb::task_arena>(1);
574#else
575 if (arena == nullptr) {
576 mTaskArena = std::make_shared<tbb::task_arena>(std::abs(n));
577 LOGP(info, "Setting seeding vertexer with {} threads.", n);
578 } else {
579 mTaskArena = arena;
580 LOGP(info, "Attaching vertexer to calling thread's arena");
581 }
582#endif
583}
584
585template <int nLayers>
587{
588 auto stream = new utils::TreeStreamRedirector("artefacts_tf.root", "recreate");
589 LOGP(info, "writing debug output for computeTracklets");
590 for (int rofId{0}; rofId < mTimeFrame->getNrof(); ++rofId) {
591 const auto& strk0 = mTimeFrame->getFoundTracklets(rofId, 0);
592 std::vector<Tracklet> trk0(strk0.begin(), strk0.end());
593 const auto& strk1 = mTimeFrame->getFoundTracklets(rofId, 1);
594 std::vector<Tracklet> trk1(strk1.begin(), strk1.end());
595 (*stream) << "tracklets"
596 << "Tracklets0=" << trk0
597 << "Tracklets1=" << trk1
598 << "iteration=" << iteration
599 << "\n";
600 }
601 stream->Close();
602 delete stream;
603}
604
605template <int nLayers>
606void VertexerTraits<nLayers>::debugComputeTrackletMatching(int iteration)
607{
608 auto stream = new utils::TreeStreamRedirector("artefacts_tf.root", "update");
609 LOGP(info, "writing debug output for computeTrackletMatching");
610 for (int rofId{0}; rofId < mTimeFrame->getNrof(); ++rofId) {
611 (*stream) << "lines"
612 << "Lines=" << toSTDVector(mTimeFrame->getLines(rofId))
613 << "NTrackletCluster01=" << mTimeFrame->getNTrackletsCluster(rofId, 0)
614 << "NTrackletCluster12=" << mTimeFrame->getNTrackletsCluster(rofId, 1)
615 << "iteration=" << iteration
616 << "\n";
617 }
618
619 if (mTimeFrame->hasMCinformation()) {
620 LOGP(info, "\tdumping also MC information");
621 const auto dc = o2::steer::DigitizationContext::loadFromFile("collisioncontext.root");
622 const auto irs = dc->getEventRecords();
623 int64_t roFrameBiasInBC = o2::itsmft::DPLAlpideParam<o2::detectors::DetID::ITS>::Instance().roFrameBiasInBC;
624 int64_t roFrameLengthInBC = o2::itsmft::DPLAlpideParam<o2::detectors::DetID::ITS>::Instance().roFrameLengthInBC;
626
627 std::map<int, int> eve2BcInROF, bcInRofNEve;
628 for (int iSrc{0}; iSrc < mcReader.getNSources(); ++iSrc) {
629 auto eveId2colId = dc->getCollisionIndicesForSource(iSrc);
630 for (int iEve{0}; iEve < mcReader.getNEvents(iSrc); ++iEve) {
631 const auto& ir = irs[eveId2colId[iEve]];
632 if (!ir.isDummy()) { // do we need this, is this for diffractive events?
633 const auto& eve = mcReader.getMCEventHeader(iSrc, iEve);
634 const int bcInROF = ((ir - raw::HBFUtils::Instance().getFirstSampledTFIR()).toLong() - roFrameBiasInBC) % roFrameLengthInBC;
635 eve2BcInROF[iEve] = bcInROF;
636 ++bcInRofNEve[bcInROF];
637 }
638 }
639 }
640
641 std::unordered_map<int, int> bcROFNTracklets01, bcROFNTracklets12;
642 std::vector<std::vector<int>> tracklet01BC, tracklet12BC;
643 for (int rofId{0}; rofId < mTimeFrame->getNrof(); ++rofId) {
644 { // 0-1
645 const auto& tracklet01 = mTimeFrame->getFoundTracklets(rofId, 0);
646 const auto& lbls01 = mTimeFrame->getLabelsFoundTracklets(rofId, 0);
647 auto& trkls01 = tracklet01BC.emplace_back();
648 for (int iTrklt{0}; iTrklt < (int)tracklet01.size(); ++iTrklt) {
649 const auto& tracklet = tracklet01[iTrklt];
650 const auto& lbl = lbls01[iTrklt];
651 if (lbl.isCorrect()) {
652 ++bcROFNTracklets01[eve2BcInROF[lbl.getEventID()]];
653 trkls01.push_back(eve2BcInROF[lbl.getEventID()]);
654 } else {
655 trkls01.push_back(-1);
656 }
657 }
658 }
659 { // 1-2 computed on the fly!
660 const auto& tracklet12 = mTimeFrame->getFoundTracklets(rofId, 1);
661 auto& trkls12 = tracklet12BC.emplace_back();
662 for (int iTrklt{0}; iTrklt < (int)tracklet12.size(); ++iTrklt) {
663 const auto& tracklet = tracklet12[iTrklt];
665
666 int sortedId1{mTimeFrame->getSortedIndex(tracklet.rof[0], 1, tracklet.firstClusterIndex)};
667 int sortedId2{mTimeFrame->getSortedIndex(tracklet.rof[1], 2, tracklet.secondClusterIndex)};
668 for (const auto& lab1 : mTimeFrame->getClusterLabels(1, mTimeFrame->getClusters()[1][sortedId1].clusterId)) {
669 for (const auto& lab2 : mTimeFrame->getClusterLabels(2, mTimeFrame->getClusters()[2][sortedId2].clusterId)) {
670 if (lab1 == lab2 && lab1.isValid()) {
671 label = lab1;
672 break;
673 }
674 }
675 if (label.isValid()) {
676 break;
677 }
678 }
679
680 if (label.isCorrect()) {
681 ++bcROFNTracklets12[eve2BcInROF[label.getEventID()]];
682 trkls12.push_back(eve2BcInROF[label.getEventID()]);
683 } else {
684 trkls12.push_back(-1);
685 }
686 }
687 }
688 }
689 LOGP(info, "\tdumping ntracklets/RofBC ({})", bcInRofNEve.size());
690 for (const auto& [bcInRof, neve] : bcInRofNEve) {
691 (*stream) << "ntracklets"
692 << "bcInROF=" << bcInRof
693 << "ntrkl01=" << bcROFNTracklets01[bcInRof]
694 << "ntrkl12=" << bcROFNTracklets12[bcInRof]
695 << "neve=" << neve
696 << "iteration=" << iteration
697 << "\n";
698 }
699
700 std::unordered_map<int, int> bcROFNLines;
701 for (int rofId{0}; rofId < mTimeFrame->getNrof(); ++rofId) {
702 const auto& lines = mTimeFrame->getLines(rofId);
703 const auto& lbls = mTimeFrame->getLinesLabel(rofId);
704 for (int iLine{0}; iLine < (int)lines.size(); ++iLine) {
705 const auto& line = lines[iLine];
706 const auto& lbl = lbls[iLine];
707 if (lbl.isCorrect()) {
708 ++bcROFNLines[eve2BcInROF[lbl.getEventID()]];
709 }
710 }
711 }
712
713 LOGP(info, "\tdumping nlines/RofBC");
714 for (const auto& [bcInRof, neve] : bcInRofNEve) {
715 (*stream) << "nlines"
716 << "bcInROF=" << bcInRof
717 << "nline=" << bcROFNLines[bcInRof]
718 << "neve=" << neve
719 << "iteration=" << iteration
720 << "\n";
721 }
722 }
723 stream->Close();
724 delete stream;
725}
726
727template <int nLayers>
728void VertexerTraits<nLayers>::debugComputeVertices(int iteration)
729{
730 auto stream = new utils::TreeStreamRedirector("artefacts_tf.root", "update");
731 LOGP(info, "writing debug output for computeVertices");
732 for (auto rofId{0}; rofId < mTimeFrame->getNrof(); ++rofId) {
733 (*stream) << "clusterlines"
734 << "clines_post=" << toSTDVector(mTimeFrame->getTrackletClusters(rofId))
735 << "iteration=" << iteration
736 << "\n";
737 }
738
739 if (mTimeFrame->hasMCinformation()) {
740 LOGP(info, "\tdumping also MC information");
741 const auto dc = o2::steer::DigitizationContext::loadFromFile("collisioncontext.root");
742 const auto irs = dc->getEventRecords();
743 int64_t roFrameBiasInBC = o2::itsmft::DPLAlpideParam<o2::detectors::DetID::ITS>::Instance().roFrameBiasInBC;
744 int64_t roFrameLengthInBC = o2::itsmft::DPLAlpideParam<o2::detectors::DetID::ITS>::Instance().roFrameLengthInBC;
746
747 std::map<int, int> eve2BcInROF, bcInRofNEve;
748 for (int iSrc{0}; iSrc < mcReader.getNSources(); ++iSrc) {
749 auto eveId2colId = dc->getCollisionIndicesForSource(iSrc);
750 for (int iEve{0}; iEve < mcReader.getNEvents(iSrc); ++iEve) {
751 const auto& ir = irs[eveId2colId[iEve]];
752 if (!ir.isDummy()) { // do we need this, is this for diffractive events?
753 const auto& eve = mcReader.getMCEventHeader(iSrc, iEve);
754 const int bcInROF = ((ir - raw::HBFUtils::Instance().getFirstSampledTFIR()).toLong() - roFrameBiasInBC) % roFrameLengthInBC;
755 eve2BcInROF[iEve] = bcInROF;
756 ++bcInRofNEve[bcInROF];
757 }
758 }
759 }
760
761 std::unordered_map<int, int> bcROFNVtx;
762 std::unordered_map<int, float> bcROFNPur;
763 std::unordered_map<o2::MCCompLabel, size_t> uniqueVertices;
764 for (int rofId{0}; rofId < mTimeFrame->getNrof(); ++rofId) {
765 const auto& pvs = mTimeFrame->getPrimaryVertices(rofId);
766 const auto& lblspv = mTimeFrame->getPrimaryVerticesMCRecInfo(rofId);
767 for (int i{0}; i < (int)pvs.size(); ++i) {
768 const auto& pv = pvs[i];
769 const auto& [lbl, pur] = lblspv[i];
770 if (lbl.isCorrect()) {
771 ++uniqueVertices[lbl];
772 ++bcROFNVtx[eve2BcInROF[lbl.getEventID()]];
773 bcROFNPur[eve2BcInROF[lbl.getEventID()]] += pur;
774 }
775 }
776 }
777
778 std::unordered_map<int, int> bcROFNUVtx, bcROFNCVtx;
779 for (const auto& [k, _] : eve2BcInROF) {
780 bcROFNUVtx[k] = bcROFNCVtx[k] = 0;
781 }
782
783 for (const auto& [lbl, c] : uniqueVertices) {
784 if (c <= 1) {
785 ++bcROFNUVtx[eve2BcInROF[lbl.getEventID()]];
786 } else {
787 ++bcROFNCVtx[eve2BcInROF[lbl.getEventID()]];
788 }
789 }
790
791 LOGP(info, "\tdumping nvtx/RofBC");
792 for (const auto& [bcInRof, neve] : bcInRofNEve) {
793 (*stream) << "nvtx"
794 << "bcInROF=" << bcInRof
795 << "nvtx=" << bcROFNVtx[bcInRof] // all vertices
796 << "nuvtx=" << bcROFNUVtx[bcInRof] // unique vertices
797 << "ncvtx=" << bcROFNCVtx[bcInRof] // cloned vertices
798 << "npur=" << bcROFNPur[bcInRof]
799 << "neve=" << neve
800 << "iteration=" << iteration
801 << "\n";
802 }
803
804 // check dist of clones
805 std::unordered_map<o2::MCCompLabel, std::vector<Vertex>> cVtx;
806 for (int rofId{0}; rofId < mTimeFrame->getNrof(); ++rofId) {
807 const auto& pvs = mTimeFrame->getPrimaryVertices(rofId);
808 const auto& lblspv = mTimeFrame->getPrimaryVerticesMCRecInfo(rofId);
809 for (int i{0}; i < (int)pvs.size(); ++i) {
810 const auto& pv = pvs[i];
811 const auto& [lbl, pur] = lblspv[i];
812 if (lbl.isCorrect() && uniqueVertices.contains(lbl) && uniqueVertices[lbl] > 1) {
813 if (!cVtx.contains(lbl)) {
814 cVtx[lbl] = std::vector<Vertex>();
815 }
816 cVtx[lbl].push_back(pv);
817 }
818 }
819 }
820
821 for (auto& [_, vertices] : cVtx) {
822 std::sort(vertices.begin(), vertices.end(), [](const Vertex& a, const Vertex& b) { return a.getNContributors() > b.getNContributors(); });
823 for (int i{0}; i < (int)vertices.size(); ++i) {
824 const auto vtx = vertices[i];
825 (*stream) << "cvtx"
826 << "vertex=" << vtx
827 << "i=" << i
828 << "dx=" << vertices[0].getX() - vtx.getX()
829 << "dy=" << vertices[0].getY() - vtx.getY()
830 << "dz=" << vertices[0].getZ() - vtx.getZ()
831 << "drof=" << vertices[0].getTimeStamp().getTimeStamp() - vtx.getTimeStamp().getTimeStamp()
832 << "dnc=" << vertices[0].getNContributors() - vtx.getNContributors()
833 << "iteration=" << iteration
834 << "\n";
835 }
836 }
837 }
838 stream->Close();
839 delete stream;
840}
841
842template class VertexerTraits<7>;
843} // namespace o2::its
uint64_t vertex
Definition RawEventData.h:9
int32_t i
Mode
Definition Utils.h:89
uint32_t c
Definition RawData.h:2
Class to compute the primary vertex in ITS from tracklets.
static constexpr int maxTrackID()
HMPID cluster implementation.
Definition Cluster.h:27
virtual void computeTracklets(const int iteration=0)
virtual void updateVertexingParameters(const std::vector< VertexingParameters > &vrtPar, const TimeFrameGPUParameters &gpuTfPar)
static DigitizationContext * loadFromFile(std::string_view filename="")
size_t getNEvents(int source) const
Get number of events.
o2::dataformats::MCEventHeader const & getMCEventHeader(int source, int event) const
retrieves the MCEventHeader for a given eventID and sourceID
void releaseTracksForSourceAndEvent(int source, int event)
API to ask releasing tracks (freeing memory) for source + event.
std::vector< MCTrack > const & getTracks(int source, int event) const
variant returning all tracks for source and event at once
GLdouble n
Definition glcorearb.h:1982
GLuint index
Definition glcorearb.h:781
GLboolean GLboolean GLboolean b
Definition glcorearb.h:1233
GLsizei GLsizei GLfloat distance
Definition glcorearb.h:5506
GLuint GLsizei const GLchar * label
Definition glcorearb.h:2519
GLboolean GLboolean GLboolean GLboolean a
Definition glcorearb.h:1233
GLuint GLuint stream
Definition glcorearb.h:1806
constexpr float TwoPI
int32_t const char int32_t line
void deepVectorClear(std::vector< T > &vec)
std::pmr::vector< T > bounded_vector
std::vector< T > toSTDVector(const bounded_vector< T > &b)
std::vector< Cluster > getClusters(int event)
Common utility functions.
int32_t w
o2::InteractionRecord ir(0, 0)
std::vector< Tracklet64 > tracklets