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