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 <iostream>
14#include <memory>
15#include <string>
16#include <chrono>
17
18#include <oneapi/tbb/blocked_range.h>
19#include <oneapi/tbb/parallel_for.h>
20
21#include "ITStracking/VertexerTraits.h"
22#include "ITStracking/BoundedAllocator.h"
23#include "ITStracking/ClusterLines.h"
24#include "ITStracking/Tracklet.h"
25
26#ifdef VTX_DEBUG
27#include "TTree.h"
28#include "TFile.h"
29#include <fstream>
30#include <ostream>
31#endif
32
33using namespace o2::its;
34
35template <TrackletMode Mode, bool EvalRun>
36void 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& 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::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
96void 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<uint8_t>& 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 pivotRofId,
110 const short targetRofId,
111 const float tanLambdaCut = 0.025f,
112 const float phiCut = 0.005f,
113 const int maxTracklets = static_cast<int>(1e2))
114{
115 int offset01{0}, offset12{0};
116 for (unsigned int iCurrentLayerClusterIndex{0}; iCurrentLayerClusterIndex < clusters1.size(); ++iCurrentLayerClusterIndex) {
117 int validTracklets{0};
118 for (int iTracklet12{offset12}; iTracklet12 < offset12 + foundTracklets12[iCurrentLayerClusterIndex]; ++iTracklet12) {
119 for (int iTracklet01{offset01}; iTracklet01 < offset01 + foundTracklets01[iCurrentLayerClusterIndex]; ++iTracklet01) {
120 const auto& tracklet01{tracklets01[iTracklet01]};
121 const auto& tracklet12{tracklets12[iTracklet12]};
122 if (tracklet01.rof[0] != targetRofId || tracklet12.rof[1] != targetRofId) {
123 continue;
124 }
125 const float deltaTanLambda{o2::gpu::GPUCommonMath::Abs(tracklet01.tanLambda - tracklet12.tanLambda)};
126 const float deltaPhi{o2::gpu::GPUCommonMath::Abs(math_utils::smallestAngleDifference(tracklet01.phi, tracklet12.phi))};
127 if (!usedTracklets[iTracklet01] && deltaTanLambda < tanLambdaCut && deltaPhi < phiCut && validTracklets != maxTracklets) {
128 usedClusters0[tracklet01.firstClusterIndex] = true;
129 usedClusters2[tracklet12.secondClusterIndex] = true;
130 usedTracklets[iTracklet01] = true;
131 lines.emplace_back(tracklet01, clusters0.data(), clusters1.data());
132 if (!trackletLabels.empty()) {
133 linesLabels.emplace_back(trackletLabels[iTracklet01]);
134 }
135 ++validTracklets;
136 }
137 }
138 }
139 offset01 += foundTracklets01[iCurrentLayerClusterIndex];
140 offset12 += foundTracklets12[iCurrentLayerClusterIndex];
141 }
142}
143
144bounded_vector<std::pair<int, int>> VertexerTraits::selectClusters(const int* indexTable,
145 const std::array<int, 4>& selectedBinsRect,
146 const IndexTableUtils& utils)
147{
148 bounded_vector<std::pair<int, int>> filteredBins{mMemoryPool.get()};
149 int phiBinsNum{selectedBinsRect[3] - selectedBinsRect[1] + 1};
150 if (phiBinsNum < 0) {
151 phiBinsNum += utils.getNphiBins();
152 }
153 filteredBins.reserve(phiBinsNum);
154 for (int iPhiBin{selectedBinsRect[1]}, iPhiCount{0}; iPhiCount < phiBinsNum;
155 iPhiBin = ++iPhiBin == utils.getNphiBins() ? 0 : iPhiBin, iPhiCount++) {
156 const int firstBinIndex{utils.getBinIndex(selectedBinsRect[0], iPhiBin)};
157 filteredBins.emplace_back(
158 indexTable[firstBinIndex],
159 utils.countRowSelectedBins(indexTable, iPhiBin, selectedBinsRect[0], selectedBinsRect[2]));
160 }
161 return filteredBins;
162}
163
164void VertexerTraits::updateVertexingParameters(const std::vector<VertexingParameters>& vrtPar, const TimeFrameGPUParameters& tfPar)
165{
166 mVrtParams = vrtPar;
167 mIndexTableUtils.setTrackingParameters(vrtPar[0]);
168 for (auto& par : mVrtParams) {
169 par.phiSpan = static_cast<int>(std::ceil(mIndexTableUtils.getNphiBins() * par.phiCut / o2::constants::math::TwoPI));
170 par.zSpan = static_cast<int>(std::ceil(par.zCut * mIndexTableUtils.getInverseZCoordinate(0)));
171 }
172}
173
174// Main functions
175void VertexerTraits::computeTracklets(const int iteration)
176{
177 mTaskArena->execute([&] {
178 tbb::parallel_for(
179 tbb::blocked_range<short>(0, (short)mTimeFrame->getNrof()),
180 [&](const tbb::blocked_range<short>& Rofs) {
181 for (short pivotRofId = Rofs.begin(); pivotRofId < Rofs.end(); ++pivotRofId) {
182 bool skipROF = iteration && (int)mTimeFrame->getPrimaryVertices(pivotRofId).size() > mVrtParams[iteration].vertPerRofThreshold;
183 short startROF{std::max((short)0, static_cast<short>(pivotRofId - mVrtParams[iteration].deltaRof))};
184 short endROF{std::min(static_cast<short>(mTimeFrame->getNrof()), static_cast<short>(pivotRofId + mVrtParams[iteration].deltaRof + 1))};
185 for (auto targetRofId = startROF; targetRofId < endROF; ++targetRofId) {
186 trackleterKernelHost<TrackletMode::Layer0Layer1, true>(
187 !skipROF ? mTimeFrame->getClustersOnLayer(targetRofId, 0) : gsl::span<Cluster>(), // Clusters to be matched with the next layer in target rof
188 !skipROF ? mTimeFrame->getClustersOnLayer(pivotRofId, 1) : gsl::span<Cluster>(), // Clusters to be matched with the current layer in pivot rof
189 mTimeFrame->getUsedClustersROF(targetRofId, 0), // Span of the used clusters in the target rof
190 mTimeFrame->getIndexTable(targetRofId, 0).data(), // Index table to access the data on the next layer in target rof
191 mVrtParams[iteration].phiCut,
192 mTimeFrame->getTracklets()[0], // Flat tracklet buffer
193 mTimeFrame->getNTrackletsCluster(pivotRofId, 0), // Span of the number of tracklets per each cluster in pivot rof
194 mIndexTableUtils,
195 pivotRofId,
196 targetRofId,
197 gsl::span<int>(), // Offset in the tracklet buffer
198 mVrtParams[iteration].maxTrackletsPerCluster);
199 trackleterKernelHost<TrackletMode::Layer1Layer2, true>(
200 !skipROF ? mTimeFrame->getClustersOnLayer(targetRofId, 2) : gsl::span<Cluster>(),
201 !skipROF ? mTimeFrame->getClustersOnLayer(pivotRofId, 1) : gsl::span<Cluster>(),
202 mTimeFrame->getUsedClustersROF(targetRofId, 2),
203 mTimeFrame->getIndexTable(targetRofId, 2).data(),
204 mVrtParams[iteration].phiCut,
205 mTimeFrame->getTracklets()[1],
206 mTimeFrame->getNTrackletsCluster(pivotRofId, 1), // Span of the number of tracklets per each cluster in pivot rof
207 mIndexTableUtils,
208 pivotRofId,
209 targetRofId,
210 gsl::span<int>(), // Offset in the tracklet buffer
211 mVrtParams[iteration].maxTrackletsPerCluster);
212 }
213 mTimeFrame->getNTrackletsROF(pivotRofId, 0) = std::accumulate(mTimeFrame->getNTrackletsCluster(pivotRofId, 0).begin(), mTimeFrame->getNTrackletsCluster(pivotRofId, 0).end(), 0);
214 mTimeFrame->getNTrackletsROF(pivotRofId, 1) = std::accumulate(mTimeFrame->getNTrackletsCluster(pivotRofId, 1).begin(), mTimeFrame->getNTrackletsCluster(pivotRofId, 1).end(), 0);
215 }
216 });
217 });
218
219 mTimeFrame->computeTrackletsPerROFScans();
220 mTimeFrame->getTracklets()[0].resize(mTimeFrame->getTotalTrackletsTF(0));
221 mTimeFrame->getTracklets()[1].resize(mTimeFrame->getTotalTrackletsTF(1));
222
223 mTaskArena->execute([&] {
224 tbb::parallel_for(
225 tbb::blocked_range<short>(0, (short)mTimeFrame->getNrof()),
226 [&](const tbb::blocked_range<short>& Rofs) {
227 for (short pivotRofId = Rofs.begin(); pivotRofId < Rofs.end(); ++pivotRofId) {
228 bool skipROF = iteration && (int)mTimeFrame->getPrimaryVertices(pivotRofId).size() > mVrtParams[iteration].vertPerRofThreshold;
229 short startROF{std::max((short)0, static_cast<short>(pivotRofId - mVrtParams[iteration].deltaRof))};
230 short endROF{std::min(static_cast<short>(mTimeFrame->getNrof()), static_cast<short>(pivotRofId + mVrtParams[iteration].deltaRof + 1))};
231 auto mobileOffset0 = mTimeFrame->getNTrackletsROF(pivotRofId, 0);
232 auto mobileOffset1 = mTimeFrame->getNTrackletsROF(pivotRofId, 1);
233 for (auto targetRofId = startROF; targetRofId < endROF; ++targetRofId) {
234 trackleterKernelHost<TrackletMode::Layer0Layer1, false>(
235 !skipROF ? mTimeFrame->getClustersOnLayer(targetRofId, 0) : gsl::span<Cluster>(),
236 !skipROF ? mTimeFrame->getClustersOnLayer(pivotRofId, 1) : gsl::span<Cluster>(),
237 mTimeFrame->getUsedClustersROF(targetRofId, 0),
238 mTimeFrame->getIndexTable(targetRofId, 0).data(),
239 mVrtParams[iteration].phiCut,
240 mTimeFrame->getTracklets()[0],
241 mTimeFrame->getNTrackletsCluster(pivotRofId, 0),
242 mIndexTableUtils,
243 pivotRofId,
244 targetRofId,
245 mTimeFrame->getExclusiveNTrackletsCluster(pivotRofId, 0),
246 mVrtParams[iteration].maxTrackletsPerCluster);
247 trackleterKernelHost<TrackletMode::Layer1Layer2, false>(
248 !skipROF ? mTimeFrame->getClustersOnLayer(targetRofId, 2) : gsl::span<Cluster>(),
249 !skipROF ? mTimeFrame->getClustersOnLayer(pivotRofId, 1) : gsl::span<Cluster>(),
250 mTimeFrame->getUsedClustersROF(targetRofId, 2),
251 mTimeFrame->getIndexTable(targetRofId, 2).data(),
252 mVrtParams[iteration].phiCut,
253 mTimeFrame->getTracklets()[1],
254 mTimeFrame->getNTrackletsCluster(pivotRofId, 1),
255 mIndexTableUtils,
256 pivotRofId,
257 targetRofId,
258 mTimeFrame->getExclusiveNTrackletsCluster(pivotRofId, 1),
259 mVrtParams[iteration].maxTrackletsPerCluster);
260 }
261 }
262 });
263 });
264
266 if (mTimeFrame->hasMCinformation()) {
267 for (const auto& trk : mTimeFrame->getTracklets()[0]) {
268 o2::MCCompLabel label;
269 int sortedId0{mTimeFrame->getSortedIndex(trk.rof[0], 0, trk.firstClusterIndex)};
270 int sortedId1{mTimeFrame->getSortedIndex(trk.rof[1], 1, trk.secondClusterIndex)};
271 for (const auto& lab0 : mTimeFrame->getClusterLabels(0, mTimeFrame->getClusters()[0][sortedId0].clusterId)) {
272 for (const auto& lab1 : mTimeFrame->getClusterLabels(1, mTimeFrame->getClusters()[1][sortedId1].clusterId)) {
273 if (lab0 == lab1 && lab0.isValid()) {
274 label = lab0;
275 break;
276 }
277 }
278 if (label.isValid()) {
279 break;
280 }
281 }
282 mTimeFrame->getTrackletsLabel(0).emplace_back(label);
283 }
284 }
285
286#ifdef VTX_DEBUG
287 // Dump on file
288 TFile* trackletFile = TFile::Open("artefacts_tf.root", "recreate");
289 TTree* tr_tre = new TTree("tracklets", "tf");
290 std::vector<o2::its::Tracklet> trkl_vec_0(0);
291 std::vector<o2::its::Tracklet> trkl_vec_1(0);
292 std::vector<o2::its::Cluster> clus0(0);
293 std::vector<o2::its::Cluster> clus1(0);
294 std::vector<o2::its::Cluster> clus2(0);
295 tr_tre->Branch("Tracklets0", &trkl_vec_0);
296 tr_tre->Branch("Tracklets1", &trkl_vec_1);
297 for (int rofId{0}; rofId < mTimeFrame->getNrof(); ++rofId) {
298 trkl_vec_0.clear();
299 trkl_vec_1.clear();
300 for (auto& tr : mTimeFrame->getFoundTracklets(rofId, 0)) {
301 trkl_vec_0.push_back(tr);
302 }
303 for (auto& tr : mTimeFrame->getFoundTracklets(rofId, 1)) {
304 trkl_vec_1.push_back(tr);
305 }
306 tr_tre->Fill();
307 }
308 trackletFile->cd();
309 tr_tre->Write();
310 trackletFile->Close();
311
312 std::ofstream out01("NTC01_cpu.txt"), out12("NTC12_cpu.txt");
313 for (int iRof{0}; iRof < mTimeFrame->getNrof(); ++iRof) {
314 out01 << "ROF: " << iRof << std::endl;
315 out12 << "ROF: " << iRof << std::endl;
316 std::copy(mTimeFrame->getNTrackletsCluster(iRof, 0).begin(), mTimeFrame->getNTrackletsCluster(iRof, 0).end(), std::ostream_iterator<double>(out01, "\t"));
317 out01 << std::endl;
318 std::copy(mTimeFrame->getExclusiveNTrackletsCluster(iRof, 0).begin(), mTimeFrame->getExclusiveNTrackletsCluster(iRof, 0).end(), std::ostream_iterator<double>(out01, "\t"));
319 std::copy(mTimeFrame->getNTrackletsCluster(iRof, 1).begin(), mTimeFrame->getNTrackletsCluster(iRof, 1).end(), std::ostream_iterator<double>(out12, "\t"));
320 out12 << std::endl;
321 std::copy(mTimeFrame->getExclusiveNTrackletsCluster(iRof, 1).begin(), mTimeFrame->getExclusiveNTrackletsCluster(iRof, 1).end(), std::ostream_iterator<double>(out12, "\t"));
322 out01 << std::endl;
323 out12 << std::endl;
324 }
325 out01.close();
326 out12.close();
327#endif
328}
329
330void VertexerTraits::computeTrackletMatching(const int iteration)
331{
332 mTaskArena->execute([&] {
333 tbb::parallel_for(
334 tbb::blocked_range<short>(0, (short)mTimeFrame->getNrof()),
335 [&](const tbb::blocked_range<short>& Rofs) {
336 for (short pivotRofId = Rofs.begin(); pivotRofId < Rofs.end(); ++pivotRofId) {
337 if (iteration && (int)mTimeFrame->getPrimaryVertices(pivotRofId).size() > mVrtParams[iteration].vertPerRofThreshold) {
338 continue;
339 }
340 mTimeFrame->getLines(pivotRofId).reserve(mTimeFrame->getNTrackletsCluster(pivotRofId, 0).size());
341 bounded_vector<uint8_t> usedTracklets(mTimeFrame->getFoundTracklets(pivotRofId, 0).size(), false, mMemoryPool.get());
342 short startROF{std::max((short)0, static_cast<short>(pivotRofId - mVrtParams[iteration].deltaRof))};
343 short endROF{std::min(static_cast<short>(mTimeFrame->getNrof()), static_cast<short>(pivotRofId + mVrtParams[iteration].deltaRof + 1))};
344 for (short targetRofId = startROF; targetRofId < endROF; ++targetRofId) {
345 trackletSelectionKernelHost(
346 mTimeFrame->getClustersOnLayer(targetRofId, 0),
347 mTimeFrame->getClustersOnLayer(pivotRofId, 1),
348 mTimeFrame->getUsedClustersROF(targetRofId, 0),
349 mTimeFrame->getUsedClustersROF(targetRofId, 2),
350 mTimeFrame->getFoundTracklets(pivotRofId, 0),
351 mTimeFrame->getFoundTracklets(pivotRofId, 1),
352 usedTracklets,
353 mTimeFrame->getNTrackletsCluster(pivotRofId, 0),
354 mTimeFrame->getNTrackletsCluster(pivotRofId, 1),
355 mTimeFrame->getLines(pivotRofId),
356 mTimeFrame->getLabelsFoundTracklets(pivotRofId, 0),
357 mTimeFrame->getLinesLabel(pivotRofId),
358 pivotRofId,
359 targetRofId,
360 mVrtParams[iteration].tanLambdaCut,
361 mVrtParams[iteration].phiCut);
362 }
363 }
364 });
365 });
366
367#ifdef VTX_DEBUG
368 TFile* trackletFile = TFile::Open("artefacts_tf.root", "update");
369 TTree* ln_tre = new TTree("lines", "tf");
370 std::vector<o2::its::Line> lines_vec(0);
371 std::vector<int> nTrackl01(0);
372 std::vector<int> nTrackl12(0);
373 ln_tre->Branch("Lines", &lines_vec);
374 ln_tre->Branch("NTrackletCluster01", &nTrackl01);
375 ln_tre->Branch("NTrackletCluster12", &nTrackl12);
376 for (int rofId{0}; rofId < mTimeFrame->getNrof(); ++rofId) {
377 lines_vec.clear();
378 nTrackl01.clear();
379 nTrackl12.clear();
380 for (auto& ln : mTimeFrame->getLines(rofId)) {
381 lines_vec.push_back(ln);
382 }
383 for (auto& n : mTimeFrame->getNTrackletsCluster(rofId, 0)) {
384 nTrackl01.push_back(n);
385 }
386 for (auto& n : mTimeFrame->getNTrackletsCluster(rofId, 1)) {
387 nTrackl12.push_back(n);
388 }
389
390 ln_tre->Fill();
391 }
392 trackletFile->cd();
393 ln_tre->Write();
394 trackletFile->Close();
395#endif
396}
397
398void VertexerTraits::computeVertices(const int iteration)
399{
400 auto nsigmaCut{std::min(mVrtParams[iteration].vertNsigmaCut * mVrtParams[iteration].vertNsigmaCut * (mVrtParams[iteration].vertRadiusSigma * mVrtParams[iteration].vertRadiusSigma + mVrtParams[iteration].trackletSigma * mVrtParams[iteration].trackletSigma), 1.98f)};
401 bounded_vector<Vertex> vertices(mMemoryPool.get());
402 bounded_vector<std::pair<o2::MCCompLabel, float>> polls(mMemoryPool.get());
403#ifdef VTX_DEBUG
404 std::vector<std::vector<ClusterLines>> dbg_clusLines(mTimeFrame->getNrof());
405#endif
406 bounded_vector<int> noClustersVec(mTimeFrame->getNrof(), 0, mMemoryPool.get());
407 for (int rofId{0}; rofId < mTimeFrame->getNrof(); ++rofId) {
408 if (iteration && (int)mTimeFrame->getPrimaryVertices(rofId).size() > mVrtParams[iteration].vertPerRofThreshold) {
409 continue;
410 }
411 const int numTracklets{static_cast<int>(mTimeFrame->getLines(rofId).size())};
412
413 bounded_vector<bool> usedTracklets(numTracklets, false, mMemoryPool.get());
414 for (int line1{0}; line1 < numTracklets; ++line1) {
415 if (usedTracklets[line1]) {
416 continue;
417 }
418 for (int line2{line1 + 1}; line2 < numTracklets; ++line2) {
419 if (usedTracklets[line2]) {
420 continue;
421 }
422 auto dca{Line::getDCA(mTimeFrame->getLines(rofId)[line1], mTimeFrame->getLines(rofId)[line2])};
423 if (dca < mVrtParams[iteration].pairCut) {
424 mTimeFrame->getTrackletClusters(rofId).emplace_back(line1, mTimeFrame->getLines(rofId)[line1], line2, mTimeFrame->getLines(rofId)[line2]);
425 std::array<float, 3> tmpVertex{mTimeFrame->getTrackletClusters(rofId).back().getVertex()};
426 if (tmpVertex[0] * tmpVertex[0] + tmpVertex[1] * tmpVertex[1] > 4.f) {
427 mTimeFrame->getTrackletClusters(rofId).pop_back();
428 break;
429 }
430 usedTracklets[line1] = true;
431 usedTracklets[line2] = true;
432 for (int tracklet3{0}; tracklet3 < numTracklets; ++tracklet3) {
433 if (usedTracklets[tracklet3]) {
434 continue;
435 }
436 if (Line::getDistanceFromPoint(mTimeFrame->getLines(rofId)[tracklet3], tmpVertex) < mVrtParams[iteration].pairCut) {
437 mTimeFrame->getTrackletClusters(rofId).back().add(tracklet3, mTimeFrame->getLines(rofId)[tracklet3]);
438 usedTracklets[tracklet3] = true;
439 tmpVertex = mTimeFrame->getTrackletClusters(rofId).back().getVertex();
440 }
441 }
442 break;
443 }
444 }
445 }
446 if (mVrtParams[iteration].allowSingleContribClusters) {
447 auto beamLine = Line{{mTimeFrame->getBeamX(), mTimeFrame->getBeamY(), -50.f}, {mTimeFrame->getBeamX(), mTimeFrame->getBeamY(), 50.f}}; // use beam position as contributor
448 for (size_t iLine{0}; iLine < numTracklets; ++iLine) {
449 if (!usedTracklets[iLine]) {
450 auto dca = Line::getDCA(mTimeFrame->getLines(rofId)[iLine], beamLine);
451 if (dca < mVrtParams[iteration].pairCut) {
452 mTimeFrame->getTrackletClusters(rofId).emplace_back(iLine, mTimeFrame->getLines(rofId)[iLine], -1, beamLine); // beamline must be passed as second line argument
453 }
454 }
455 }
456 }
457
458 // Cluster merging
459 std::sort(mTimeFrame->getTrackletClusters(rofId).begin(), mTimeFrame->getTrackletClusters(rofId).end(),
460 [](ClusterLines& cluster1, ClusterLines& cluster2) { return cluster1.getSize() > cluster2.getSize(); });
461 noClustersVec[rofId] = static_cast<int>(mTimeFrame->getTrackletClusters(rofId).size());
462 for (int iCluster1{0}; iCluster1 < noClustersVec[rofId]; ++iCluster1) {
463 std::array<float, 3> vertex1{mTimeFrame->getTrackletClusters(rofId)[iCluster1].getVertex()};
464 std::array<float, 3> vertex2{};
465 for (int iCluster2{iCluster1 + 1}; iCluster2 < noClustersVec[rofId]; ++iCluster2) {
466 vertex2 = mTimeFrame->getTrackletClusters(rofId)[iCluster2].getVertex();
467 if (o2::gpu::GPUCommonMath::Abs(vertex1[2] - vertex2[2]) < mVrtParams[iteration].clusterCut) {
468 float distance{(vertex1[0] - vertex2[0]) * (vertex1[0] - vertex2[0]) +
469 (vertex1[1] - vertex2[1]) * (vertex1[1] - vertex2[1]) +
470 (vertex1[2] - vertex2[2]) * (vertex1[2] - vertex2[2])};
471 if (distance < mVrtParams[iteration].pairCut * mVrtParams[iteration].pairCut) {
472 for (auto label : mTimeFrame->getTrackletClusters(rofId)[iCluster2].getLabels()) {
473 mTimeFrame->getTrackletClusters(rofId)[iCluster1].add(label, mTimeFrame->getLines(rofId)[label]);
474 vertex1 = mTimeFrame->getTrackletClusters(rofId)[iCluster1].getVertex();
475 }
476 mTimeFrame->getTrackletClusters(rofId).erase(mTimeFrame->getTrackletClusters(rofId).begin() + iCluster2);
477 --iCluster2;
478 --noClustersVec[rofId];
479 }
480 }
481 }
482 }
483 }
484 for (int rofId{0}; rofId < mTimeFrame->getNrof(); ++rofId) {
485 vertices.clear();
486 std::sort(mTimeFrame->getTrackletClusters(rofId).begin(), mTimeFrame->getTrackletClusters(rofId).end(),
487 [](ClusterLines& cluster1, ClusterLines& cluster2) { return cluster1.getSize() > cluster2.getSize(); }); // ensure clusters are ordered by contributors, so that we can cat after the first.
488#ifdef VTX_DEBUG
489 for (auto& cl : mTimeFrame->getTrackletClusters(rofId)) {
490 dbg_clusLines[rofId].push_back(cl);
491 }
492#endif
493 bool atLeastOneFound{false};
494 for (int iCluster{0}; iCluster < noClustersVec[rofId]; ++iCluster) {
495 bool lowMultCandidate{false};
496 double beamDistance2{(mTimeFrame->getBeamX() - mTimeFrame->getTrackletClusters(rofId)[iCluster].getVertex()[0]) * (mTimeFrame->getBeamX() - mTimeFrame->getTrackletClusters(rofId)[iCluster].getVertex()[0]) +
497 (mTimeFrame->getBeamY() - mTimeFrame->getTrackletClusters(rofId)[iCluster].getVertex()[1]) * (mTimeFrame->getBeamY() - mTimeFrame->getTrackletClusters(rofId)[iCluster].getVertex()[1])};
498 if (atLeastOneFound && (lowMultCandidate = mTimeFrame->getTrackletClusters(rofId)[iCluster].getSize() < mVrtParams[iteration].clusterContributorsCut)) { // We might have pile up with nContr > cut.
499 lowMultCandidate &= (beamDistance2 < mVrtParams[iteration].lowMultBeamDistCut * mVrtParams[iteration].lowMultBeamDistCut);
500 if (!lowMultCandidate) { // Not the first cluster and not a low multiplicity candidate, we can remove it
501 mTimeFrame->getTrackletClusters(rofId).erase(mTimeFrame->getTrackletClusters(rofId).begin() + iCluster);
502 noClustersVec[rofId]--;
503 continue;
504 }
505 }
506
507 if (beamDistance2 < nsigmaCut && o2::gpu::GPUCommonMath::Abs(mTimeFrame->getTrackletClusters(rofId)[iCluster].getVertex()[2]) < mVrtParams[iteration].maxZPositionAllowed) {
508 atLeastOneFound = true;
509 auto& vertex = vertices.emplace_back(o2::math_utils::Point3D<float>(mTimeFrame->getTrackletClusters(rofId)[iCluster].getVertex()[0],
510 mTimeFrame->getTrackletClusters(rofId)[iCluster].getVertex()[1],
511 mTimeFrame->getTrackletClusters(rofId)[iCluster].getVertex()[2]),
512 mTimeFrame->getTrackletClusters(rofId)[iCluster].getRMS2(), // Symm matrix. Diagonal: RMS2 components,
513 // off-diagonal: square mean of projections on planes.
514 mTimeFrame->getTrackletClusters(rofId)[iCluster].getSize(), // Contributors
515 mTimeFrame->getTrackletClusters(rofId)[iCluster].getAvgDistance2()); // In place of chi2
516
517 if (iteration) {
518 vertex.setFlags(Vertex::UPCMode);
519 }
520 vertex.setTimeStamp(mTimeFrame->getTrackletClusters(rofId)[iCluster].getROF());
521 if (mTimeFrame->hasMCinformation()) {
522 bounded_vector<o2::MCCompLabel> labels(mMemoryPool.get());
523 for (auto& index : mTimeFrame->getTrackletClusters(rofId)[iCluster].getLabels()) {
524 labels.push_back(mTimeFrame->getLinesLabel(rofId)[index]); // then we can use nContributors from vertices to get the labels
525 }
526 polls.push_back(computeMain(labels));
527 }
528 }
529 }
530 if (!iteration) {
531 mTimeFrame->addPrimaryVertices(vertices, rofId, iteration);
532 if (mTimeFrame->hasMCinformation()) {
533 mTimeFrame->addPrimaryVerticesLabels(polls);
534 }
535 } else {
536 mTimeFrame->addPrimaryVerticesInROF(vertices, rofId, iteration);
537 if (mTimeFrame->hasMCinformation()) {
538 mTimeFrame->addPrimaryVerticesLabelsInROF(polls, rofId);
539 }
540 }
541 if (vertices.empty() && !(iteration && (int)mTimeFrame->getPrimaryVertices(rofId).size() > mVrtParams[iteration].vertPerRofThreshold)) {
542 mTimeFrame->getNoVertexROF()++;
543 }
544 }
545#ifdef VTX_DEBUG
546 TFile* dbg_file = TFile::Open("artefacts_tf.root", "update");
547 TTree* ln_clus_lines_tree = new TTree("clusterlines", "tf");
548 std::vector<o2::its::ClusterLines> cl_lines_vec_pre(0);
549 std::vector<o2::its::ClusterLines> cl_lines_vec_post(0);
550 ln_clus_lines_tree->Branch("cllines_pre", &cl_lines_vec_pre);
551 ln_clus_lines_tree->Branch("cllines_post", &cl_lines_vec_post);
552 for (auto rofId{0}; rofId < mTimeFrame->getNrof(); ++rofId) {
553 cl_lines_vec_pre.clear();
554 cl_lines_vec_post.clear();
555 for (auto& clln : mTimeFrame->getTrackletClusters(rofId)) {
556 cl_lines_vec_post.push_back(clln);
557 }
558 for (auto& cl : dbg_clusLines[rofId]) {
559 cl_lines_vec_pre.push_back(cl);
560 }
561 ln_clus_lines_tree->Fill();
562 }
563 dbg_file->cd();
564 ln_clus_lines_tree->Write();
565 dbg_file->Close();
566#endif
567}
568
569void VertexerTraits::computeVerticesInRof(int rofId,
570 gsl::span<const o2::its::Line>& lines,
571 bounded_vector<bool>& usedLines,
572 bounded_vector<o2::its::ClusterLines>& clusterLines,
573 std::array<float, 2>& beamPosXY,
574 bounded_vector<Vertex>& vertices,
575 bounded_vector<int>& verticesInRof,
576 TimeFrame7* tf,
577 bounded_vector<o2::MCCompLabel>* labels,
578 const int iteration)
579{
580 int foundVertices{0};
581 auto nsigmaCut{std::min(mVrtParams[iteration].vertNsigmaCut * mVrtParams[iteration].vertNsigmaCut * (mVrtParams[iteration].vertRadiusSigma * mVrtParams[iteration].vertRadiusSigma + mVrtParams[iteration].trackletSigma * mVrtParams[iteration].trackletSigma), 1.98f)};
582 const int numTracklets{static_cast<int>(lines.size())};
583 for (int line1{0}; line1 < numTracklets; ++line1) {
584 if (usedLines[line1]) {
585 continue;
586 }
587 for (int line2{line1 + 1}; line2 < numTracklets; ++line2) {
588 if (usedLines[line2]) {
589 continue;
590 }
591 auto dca{Line::getDCA(lines[line1], lines[line2])};
592 if (dca < mVrtParams[iteration].pairCut) {
593 clusterLines.emplace_back(line1, lines[line1], line2, lines[line2]);
594 std::array<float, 3> tmpVertex{clusterLines.back().getVertex()};
595 if (tmpVertex[0] * tmpVertex[0] + tmpVertex[1] * tmpVertex[1] > 4.f) {
596 clusterLines.pop_back();
597 break;
598 }
599 usedLines[line1] = true;
600 usedLines[line2] = true;
601 for (int tracklet3{0}; tracklet3 < numTracklets; ++tracklet3) {
602 if (usedLines[tracklet3]) {
603 continue;
604 }
605 if (Line::getDistanceFromPoint(lines[tracklet3], tmpVertex) < mVrtParams[iteration].pairCut) {
606 clusterLines.back().add(tracklet3, lines[tracklet3]);
607 usedLines[tracklet3] = true;
608 tmpVertex = clusterLines.back().getVertex();
609 }
610 }
611 break;
612 }
613 }
614 }
615
616 if (mVrtParams[iteration].allowSingleContribClusters) {
617 auto beamLine = Line{{tf->getBeamX(), tf->getBeamY(), -50.f}, {tf->getBeamX(), tf->getBeamY(), 50.f}}; // use beam position as contributor
618 for (size_t iLine{0}; iLine < numTracklets; ++iLine) {
619 if (!usedLines[iLine]) {
620 auto dca = Line::getDCA(lines[iLine], beamLine);
621 if (dca < mVrtParams[iteration].pairCut) {
622 clusterLines.emplace_back(iLine, lines[iLine], -1, beamLine); // beamline must be passed as second line argument
623 }
624 }
625 }
626 }
627
628 // Cluster merging
629 std::sort(clusterLines.begin(), clusterLines.end(), [](ClusterLines& cluster1, ClusterLines& cluster2) { return cluster1.getSize() > cluster2.getSize(); });
630 size_t nClusters{clusterLines.size()};
631 for (int iCluster1{0}; iCluster1 < nClusters; ++iCluster1) {
632 std::array<float, 3> vertex1{clusterLines[iCluster1].getVertex()};
633 std::array<float, 3> vertex2{};
634 for (int iCluster2{iCluster1 + 1}; iCluster2 < nClusters; ++iCluster2) {
635 vertex2 = clusterLines[iCluster2].getVertex();
636 if (o2::gpu::GPUCommonMath::Abs(vertex1[2] - vertex2[2]) < mVrtParams[iteration].clusterCut) {
637 float distance{(vertex1[0] - vertex2[0]) * (vertex1[0] - vertex2[0]) +
638 (vertex1[1] - vertex2[1]) * (vertex1[1] - vertex2[1]) +
639 (vertex1[2] - vertex2[2]) * (vertex1[2] - vertex2[2])};
640 if (distance < mVrtParams[iteration].pairCut * mVrtParams[iteration].pairCut) {
641 for (auto label : clusterLines[iCluster2].getLabels()) {
642 clusterLines[iCluster1].add(label, lines[label]);
643 vertex1 = clusterLines[iCluster1].getVertex();
644 }
645 clusterLines.erase(clusterLines.begin() + iCluster2);
646 --iCluster2;
647 --nClusters;
648 }
649 }
650 }
651 }
652
653 std::sort(clusterLines.begin(), clusterLines.end(),
654 [](ClusterLines& cluster1, ClusterLines& cluster2) { return cluster1.getSize() > cluster2.getSize(); }); // ensure clusters are ordered by contributors, so that we can cut after the first.
655 bool atLeastOneFound{false};
656 for (int iCluster{0}; iCluster < nClusters; ++iCluster) {
657 bool lowMultCandidate{false};
658 double beamDistance2{(tf->getBeamX() - clusterLines[iCluster].getVertex()[0]) * (tf->getBeamX() - clusterLines[iCluster].getVertex()[0]) +
659 (tf->getBeamY() - clusterLines[iCluster].getVertex()[1]) * (tf->getBeamY() - clusterLines[iCluster].getVertex()[1])};
660
661 if (atLeastOneFound && (lowMultCandidate = clusterLines[iCluster].getSize() < mVrtParams[iteration].clusterContributorsCut)) { // We might have pile up with nContr > cut.
662 lowMultCandidate &= (beamDistance2 < mVrtParams[iteration].lowMultBeamDistCut * mVrtParams[iteration].lowMultBeamDistCut);
663 if (!lowMultCandidate) { // Not the first cluster and not a low multiplicity candidate, we can remove it
664 clusterLines.erase(clusterLines.begin() + iCluster);
665 nClusters--;
666 continue;
667 }
668 }
669 if (beamDistance2 < nsigmaCut && o2::gpu::GPUCommonMath::Abs(clusterLines[iCluster].getVertex()[2]) < mVrtParams[iteration].maxZPositionAllowed) {
670 atLeastOneFound = true;
671 ++foundVertices;
672 vertices.emplace_back(o2::math_utils::Point3D<float>(clusterLines[iCluster].getVertex()[0],
673 clusterLines[iCluster].getVertex()[1],
674 clusterLines[iCluster].getVertex()[2]),
675 clusterLines[iCluster].getRMS2(), // Symm matrix. Diagonal: RMS2 components,
676 // off-diagonal: square mean of projections on planes.
677 clusterLines[iCluster].getSize(), // Contributors
678 clusterLines[iCluster].getAvgDistance2()); // In place of chi2
679 vertices.back().setTimeStamp(clusterLines[iCluster].getROF());
680 if (labels) {
681 for (auto& index : clusterLines[iCluster].getLabels()) {
682 labels->push_back(tf->getLinesLabel(rofId)[index]); // then we can use nContributors from vertices to get the labels
683 }
684 }
685 }
686 }
687 verticesInRof.push_back(foundVertices);
688}
689
690void VertexerTraits::setNThreads(int n, std::shared_ptr<tbb::task_arena>& arena)
691{
692#if defined(VTX_DEBUG)
693 mTaskArena = std::make_shared<tbb::task_arena>(1);
694#else
695 if (arena == nullptr) {
696 mTaskArena = std::make_shared<tbb::task_arena>(std::abs(n));
697 LOGP(info, "Setting seeding vertexer with {} threads.", n);
698 } else {
699 mTaskArena = arena;
700 LOGP(info, "Attaching vertexer to calling thread's arena");
701 }
702#endif
703}
vertex
uint64_t vertex
Definition RawEventData.h:9
Mode
Mode
Definition Utils.h:89
trackletSelectionKernelHost
void trackletSelectionKernelHost(const gsl::span< const Cluster > clusters0, const gsl::span< const Cluster > clusters1, gsl::span< unsigned char > usedClusters0, gsl::span< unsigned char > usedClusters2, const gsl::span< const Tracklet > &tracklets01, const gsl::span< const Tracklet > &tracklets12, bounded_vector< uint8_t > &usedTracklets, const gsl::span< int > foundTracklets01, const gsl::span< int > foundTracklets12, bounded_vector< Line > &lines, const gsl::span< const o2::MCCompLabel > &trackletLabels, bounded_vector< o2::MCCompLabel > &linesLabels, const short pivotRofId, const short targetRofId, const float tanLambdaCut=0.025f, const float phiCut=0.005f, const int maxTracklets=static_cast< int >(1e2))
Definition VertexerTraits.cxx:96
trackleterKernelHost
void trackleterKernelHost(const gsl::span< const Cluster > &clustersNextLayer, const gsl::span< const Cluster > &clustersCurrentLayer, const gsl::span< uint8_t > &usedClustersNextLayer, int *indexTableNext, const float phiCut, bounded_vector< Tracklet > &tracklets, gsl::span< int > foundTracklets, const IndexTableUtils &utils, const short pivotRof, const short targetRof, gsl::span< int > rofFoundTrackletsOffsets, const int maxTrackletsPerCluster=static_cast< int >(2e3))
Definition VertexerTraits.cxx:36
VertexerTraits.h
Class to compute the primary vertex in ITS from tracklets.
nClusters
int nClusters
Definition bench_Clusterizer.cxx:120
o2::its::IndexTableUtils::setTrackingParameters
void setTrackingParameters(const T &params)
Definition IndexTableUtils.h:56
o2::its::IndexTableUtils::getInverseZCoordinate
float getInverseZCoordinate(const int layerIndex) const
Definition IndexTableUtils.h:69
o2::its::VertexerTraits::selectClusters
bounded_vector< std::pair< int, int > > selectClusters(const int *indexTable, const std::array< int, 4 > &selectedBinsRect, const IndexTableUtils &utils)
Definition VertexerTraits.cxx:144
o2::its::VertexerTraits::computeTrackletMatching
virtual void computeTrackletMatching(const int iteration=0)
Definition VertexerTraits.cxx:330
o2::its::VertexerTraits::computeTracklets
virtual void computeTracklets(const int iteration=0)
Definition VertexerTraits.cxx:175
o2::its::VertexerTraits::setNThreads
void setNThreads(int n, std::shared_ptr< tbb::task_arena > &arena)
Definition VertexerTraits.cxx:690
o2::its::VertexerTraits::computeVertices
virtual void computeVertices(const int iteration=0)
Definition VertexerTraits.cxx:398
o2::its::VertexerTraits::mVrtParams
std::vector< VertexingParameters > mVrtParams
Definition VertexerTraits.h:126
o2::its::VertexerTraits::mIndexTableUtils
IndexTableUtils mIndexTableUtils
Definition VertexerTraits.h:127
o2::its::VertexerTraits::computeVerticesInRof
void computeVerticesInRof(int, gsl::span< const o2::its::Line > &, bounded_vector< bool > &, bounded_vector< o2::its::ClusterLines > &, std::array< float, 2 > &, bounded_vector< Vertex > &, bounded_vector< int > &, TimeFrame7 *, bounded_vector< o2::MCCompLabel > *, const int iteration=0)
Definition VertexerTraits.cxx:569
o2::its::VertexerTraits::updateVertexingParameters
virtual void updateVertexingParameters(const std::vector< VertexingParameters > &vrtPar, const TimeFrameGPUParameters &gpuTfPar)
Definition VertexerTraits.cxx:164
o2::its::VertexerTraits::computeMain
static std::pair< o2::MCCompLabel, float > computeMain(const bounded_vector< o2::MCCompLabel > &elements)
Definition VertexerTraits.h:104
n
GLdouble n
Definition glcorearb.h:1982
index
GLuint index
Definition glcorearb.h:781
distance
GLsizei GLsizei GLfloat distance
Definition glcorearb.h:5506
label
GLuint GLsizei const GLchar * label
Definition glcorearb.h:2519
o2::constants::math::TwoPI
constexpr float TwoPI
Definition MathConstants.h:31
o2::its::bounded_vector
std::pmr::vector< T > bounded_vector
Definition BoundedAllocator.h:116
o2::mid::getClusters
std::vector< Cluster > getClusters(int event)
Definition testClusterizer.cxx:91
utils
Common utility functions.
tf
std::unique_ptr< GPUReconstructionTimeframe > tf
Definition standalone.cxx:84
int4::w
int32_t w
Definition GPUCommonDefAPI.h:70
o2::its::Line::getDistanceFromPoint
const Cluster const Cluster *static float getDistanceFromPoint(const Line &line, const std::array< float, 3 > &point)
Definition ClusterLines.h:117
o2::its::TimeFrame::getNrof
int getNrof() const
Definition TimeFrame.h:101
o2::its::TimeFrame::getBeamX
float getBeamX() const
Definition TimeFrame.h:110
o2::its::TimeFrame::hasMCinformation
bool hasMCinformation() const
Definition TimeFrame.h:147
o2::its::TimeFrame::getNoVertexROF
unsigned int & getNoVertexROF()
Definition TimeFrame.h:216
o2::its::TimeFrame::addPrimaryVerticesInROF
void addPrimaryVerticesInROF(const bounded_vector< Vertex > &vertices, const int rofId, const int iteration)
Definition TimeFrame.cxx:89
o2::its::TimeFrame::getTrackletClusters
auto & getTrackletClusters(int rofId)
Definition TimeFrame.h:207
o2::its::TimeFrame::addPrimaryVerticesLabelsInROF
void addPrimaryVerticesLabelsInROF(const bounded_vector< std::pair< MCCompLabel, float > > &labels, const int rofId)
Definition TimeFrame.cxx:99
o2::its::TimeFrame::addPrimaryVerticesLabels
void addPrimaryVerticesLabels(bounded_vector< std::pair< MCCompLabel, float > > &labels)
Definition TimeFrame.cxx:83
o2::its::TimeFrame::getPrimaryVertices
gsl::span< const Vertex > getPrimaryVertices(int rofId) const
Definition TimeFrame.h:354
o2::its::TimeFrame::getLinesLabel
auto & getLinesLabel(const int rofId)
Definition TimeFrame.h:174
o2::its::TimeFrame::getBeamY
float getBeamY() const
Definition TimeFrame.h:111
o2::its::TimeFrame::getLines
auto & getLines(int rofId)
Definition TimeFrame.h:202
o2::its::TimeFrame::addPrimaryVertices
void addPrimaryVertices(const bounded_vector< Vertex > &vertices)
Definition TimeFrame.cxx:62
tracklets
std::vector< Tracklet64 > tracklets
Definition test_ctf_io_trd.cxx:40