18#include <oneapi/tbb/blocked_range.h>
19#include <oneapi/tbb/parallel_for.h>
20#include <oneapi/tbb/combinable.h>
35template <TrackletMode Mode,
bool EvalRun,
int nLayers>
36static void trackleterKernelHost(
37 const gsl::span<const Cluster>& clustersNextLayer,
38 const gsl::span<const Cluster>& clustersCurrentLayer,
39 const gsl::span<uint8_t>& usedClustersNextLayer,
43 gsl::span<int> foundTracklets,
44 const IndexTableUtils<nLayers>&
utils,
46 const short targetRof,
47 gsl::span<int> rofFoundTrackletsOffsets,
48 const int maxTrackletsPerCluster =
static_cast<int>(2e3))
50 const int PhiBins{
utils.getNphiBins()};
51 const int ZBins{
utils.getNzBins()};
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};
60 phiBinsNum += PhiBins;
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]};
68 for (
int iNextLayerClusterIndex{firstRowClusterIndex}; iNextLayerClusterIndex < maxRowClusterIndex && iNextLayerClusterIndex < static_cast<int>(clustersNextLayer.size()); ++iNextLayerClusterIndex) {
69 if (usedClustersNextLayer[iNextLayerClusterIndex]) {
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) {
77 tracklets[rofFoundTrackletsOffsets[iCurrentLayerClusterIndex] + storedTracklets] = Tracklet{iNextLayerClusterIndex, iCurrentLayerClusterIndex, nextCluster, currentCluster, targetRof, pivotRof};
79 tracklets[rofFoundTrackletsOffsets[iCurrentLayerClusterIndex] + storedTracklets] = Tracklet{iCurrentLayerClusterIndex, iNextLayerClusterIndex, currentCluster, nextCluster, pivotRof, targetRof};
88 if constexpr (EvalRun) {
89 foundTracklets[iCurrentLayerClusterIndex] += storedTracklets;
91 rofFoundTrackletsOffsets[iCurrentLayerClusterIndex] += storedTracklets;
96static void trackletSelectionKernelHost(
97 const gsl::span<const Cluster> clusters0,
98 const gsl::span<const Cluster> clusters1,
99 gsl::span<unsigned char> usedClusters0,
100 gsl::span<unsigned char> usedClusters2,
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))
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]) {
125 const auto& tracklet01{tracklets01[iTracklet01]};
126 const auto& tracklet12{tracklets12[iTracklet12]};
128 if (tracklet01.rof[0] != targetRofId0 || tracklet12.rof[1] != targetRofId2) {
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) {
136 __atomic_store_n(&usedClusters0[tracklet01.firstClusterIndex], 1, __ATOMIC_RELAXED);
137 __atomic_store_n(&usedClusters2[tracklet12.secondClusterIndex], 1, __ATOMIC_RELAXED);
139 usedClusters0[tracklet01.firstClusterIndex] = 1;
140 usedClusters2[tracklet12.secondClusterIndex] = 1;
142 usedTracklets[iTracklet01] =
true;
143 lines.emplace_back(tracklet01, clusters0.data(), clusters1.data());
144 if (!trackletLabels.empty()) {
145 linesLabels.emplace_back(trackletLabels[iTracklet01]);
151 offset01 += foundTracklets01[iCurrentLayerClusterIndex];
152 offset12 += foundTracklets12[iCurrentLayerClusterIndex];
156template <
int nLayers>
160 mIndexTableUtils.setTrackingParameters(vrtPar[0]);
161 for (
auto& par : mVrtParams) {
163 par.zSpan =
static_cast<int>(std::ceil(par.zCut * mIndexTableUtils.getInverseZCoordinate(0)));
168template <
int nLayers>
171 mTaskArena->execute([&] {
173 tbb::blocked_range<short>(0, (
short)mTimeFrame->getNrof()),
174 [&](
const tbb::blocked_range<short>& Rofs) {
175 for (short pivotRofId = Rofs.begin(); pivotRofId < Rofs.end(); ++pivotRofId) {
176 bool skipROF = iteration && (int)mTimeFrame->getPrimaryVertices(pivotRofId).size() > mVrtParams[iteration].vertPerRofThreshold;
177 short startROF{std::max((short)0, static_cast<short>(pivotRofId - mVrtParams[iteration].deltaRof))};
178 short endROF{std::min(static_cast<short>(mTimeFrame->getNrof()), static_cast<short>(pivotRofId + mVrtParams[iteration].deltaRof + 1))};
179 for (auto targetRofId = startROF; targetRofId < endROF; ++targetRofId) {
180 trackleterKernelHost<TrackletMode::Layer0Layer1, true>(
181 !skipROF ? mTimeFrame->getClustersOnLayer(targetRofId, 0) : gsl::span<Cluster>(),
182 !skipROF ? mTimeFrame->getClustersOnLayer(pivotRofId, 1) : gsl::span<Cluster>(),
183 mTimeFrame->getUsedClustersROF(targetRofId, 0),
184 mTimeFrame->getIndexTable(targetRofId, 0).data(),
185 mVrtParams[iteration].phiCut,
186 mTimeFrame->getTracklets()[0],
187 mTimeFrame->getNTrackletsCluster(pivotRofId, 0),
192 mVrtParams[iteration].maxTrackletsPerCluster);
193 trackleterKernelHost<TrackletMode::Layer1Layer2, true>(
194 !skipROF ? mTimeFrame->getClustersOnLayer(targetRofId, 2) : gsl::span<Cluster>(),
195 !skipROF ? mTimeFrame->getClustersOnLayer(pivotRofId, 1) : gsl::span<Cluster>(),
196 mTimeFrame->getUsedClustersROF(targetRofId, 2),
197 mTimeFrame->getIndexTable(targetRofId, 2).data(),
198 mVrtParams[iteration].phiCut,
199 mTimeFrame->getTracklets()[1],
200 mTimeFrame->getNTrackletsCluster(pivotRofId, 1),
205 mVrtParams[iteration].maxTrackletsPerCluster);
207 mTimeFrame->getNTrackletsROF(pivotRofId, 0) = std::accumulate(mTimeFrame->getNTrackletsCluster(pivotRofId, 0).begin(), mTimeFrame->getNTrackletsCluster(pivotRofId, 0).end(), 0);
208 mTimeFrame->getNTrackletsROF(pivotRofId, 1) = std::accumulate(mTimeFrame->getNTrackletsCluster(pivotRofId, 1).begin(), mTimeFrame->getNTrackletsCluster(pivotRofId, 1).end(), 0);
212 mTimeFrame->computeTrackletsPerROFScans();
213 if (
auto tot0 = mTimeFrame->getTotalTrackletsTF(0), tot1 = mTimeFrame->getTotalTrackletsTF(1);
214 tot0 == 0 || tot1 == 0) {
217 mTimeFrame->getTracklets()[0].resize(tot0);
218 mTimeFrame->getTracklets()[1].resize(tot1);
222 tbb::blocked_range<short>(0, (
short)mTimeFrame->getNrof()),
223 [&](
const tbb::blocked_range<short>& Rofs) {
224 for (short pivotRofId = Rofs.begin(); pivotRofId < Rofs.end(); ++pivotRofId) {
225 bool skipROF = iteration && (int)mTimeFrame->getPrimaryVertices(pivotRofId).size() > mVrtParams[iteration].vertPerRofThreshold;
226 short startROF{std::max((short)0, static_cast<short>(pivotRofId - mVrtParams[iteration].deltaRof))};
227 short endROF{std::min(static_cast<short>(mTimeFrame->getNrof()), static_cast<short>(pivotRofId + mVrtParams[iteration].deltaRof + 1))};
228 auto mobileOffset0 = mTimeFrame->getNTrackletsROF(pivotRofId, 0);
229 auto mobileOffset1 = mTimeFrame->getNTrackletsROF(pivotRofId, 1);
230 for (auto targetRofId = startROF; targetRofId < endROF; ++targetRofId) {
231 trackleterKernelHost<TrackletMode::Layer0Layer1, false>(
232 !skipROF ? mTimeFrame->getClustersOnLayer(targetRofId, 0) : gsl::span<Cluster>(),
233 !skipROF ? mTimeFrame->getClustersOnLayer(pivotRofId, 1) : gsl::span<Cluster>(),
234 mTimeFrame->getUsedClustersROF(targetRofId, 0),
235 mTimeFrame->getIndexTable(targetRofId, 0).data(),
236 mVrtParams[iteration].phiCut,
237 mTimeFrame->getTracklets()[0],
238 mTimeFrame->getNTrackletsCluster(pivotRofId, 0),
242 mTimeFrame->getExclusiveNTrackletsCluster(pivotRofId, 0),
243 mVrtParams[iteration].maxTrackletsPerCluster);
244 trackleterKernelHost<TrackletMode::Layer1Layer2, false>(
245 !skipROF ? mTimeFrame->getClustersOnLayer(targetRofId, 2) : gsl::span<Cluster>(),
246 !skipROF ? mTimeFrame->getClustersOnLayer(pivotRofId, 1) : gsl::span<Cluster>(),
247 mTimeFrame->getUsedClustersROF(targetRofId, 2),
248 mTimeFrame->getIndexTable(targetRofId, 2).data(),
249 mVrtParams[iteration].phiCut,
250 mTimeFrame->getTracklets()[1],
251 mTimeFrame->getNTrackletsCluster(pivotRofId, 1),
255 mTimeFrame->getExclusiveNTrackletsCluster(pivotRofId, 1),
256 mVrtParams[iteration].maxTrackletsPerCluster);
263 if (mTimeFrame->hasMCinformation()) {
264 for (
const auto& trk : mTimeFrame->getTracklets()[0]) {
266 if (!trk.isEmpty()) {
267 int sortedId0{mTimeFrame->getSortedIndex(trk.rof[0], 0, trk.firstClusterIndex)};
268 int sortedId1{mTimeFrame->getSortedIndex(trk.rof[1], 1, trk.secondClusterIndex)};
269 for (
const auto& lab0 : mTimeFrame->getClusterLabels(0, mTimeFrame->
getClusters()[0][sortedId0].clusterId)) {
270 for (
const auto& lab1 : mTimeFrame->getClusterLabels(1, mTimeFrame->
getClusters()[1][sortedId1].clusterId)) {
271 if (lab0 == lab1 && lab0.isValid()) {
276 if (
label.isValid()) {
281 mTimeFrame->getTrackletsLabel(0).emplace_back(
label);
286 debugComputeTracklets(iteration);
290template <
int nLayers>
293 mTaskArena->execute([&] {
294 tbb::combinable<int> totalLines{0};
296 tbb::blocked_range<short>(0, (
short)mTimeFrame->getNrof()),
297 [&](
const tbb::blocked_range<short>& Rofs) {
298 for (short pivotRofId = Rofs.begin(); pivotRofId < Rofs.end(); ++pivotRofId) {
299 if (iteration && (int)mTimeFrame->getPrimaryVertices(pivotRofId).size() > mVrtParams[iteration].vertPerRofThreshold) {
302 if (mTimeFrame->getFoundTracklets(pivotRofId, 0).empty()) {
305 mTimeFrame->getLines(pivotRofId).reserve(mTimeFrame->getNTrackletsCluster(pivotRofId, 0).size());
306 bounded_vector<bool> usedTracklets(mTimeFrame->getFoundTracklets(pivotRofId, 0).size(), false, mMemoryPool.get());
307 short startROF{std::max((short)0, static_cast<short>(pivotRofId - mVrtParams[iteration].deltaRof))};
308 short endROF{std::min(static_cast<short>(mTimeFrame->getNrof()), static_cast<short>(pivotRofId + mVrtParams[iteration].deltaRof + 1))};
311 bool safeWrite = mTaskArena->max_concurrency() > 1 && mVrtParams[iteration].deltaRof != 0 && ((Rofs.begin() - startROF < 0) || (endROF - Rofs.end() > 0));
313 for (short targetRofId0 = startROF; targetRofId0 < endROF; ++targetRofId0) {
314 for (short targetRofId2 = startROF; targetRofId2 < endROF; ++targetRofId2) {
315 if (std::abs(targetRofId0 - targetRofId2) > mVrtParams[iteration].deltaRof) {
318 trackletSelectionKernelHost(
319 mTimeFrame->getClustersOnLayer(targetRofId0, 0),
320 mTimeFrame->getClustersOnLayer(pivotRofId, 1),
321 mTimeFrame->getUsedClustersROF(targetRofId0, 0),
322 mTimeFrame->getUsedClustersROF(targetRofId2, 2),
323 mTimeFrame->getFoundTracklets(pivotRofId, 0),
324 mTimeFrame->getFoundTracklets(pivotRofId, 1),
326 mTimeFrame->getNTrackletsCluster(pivotRofId, 0),
327 mTimeFrame->getNTrackletsCluster(pivotRofId, 1),
328 mTimeFrame->getLines(pivotRofId),
329 mTimeFrame->getLabelsFoundTracklets(pivotRofId, 0),
330 mTimeFrame->getLinesLabel(pivotRofId),
334 mVrtParams[iteration].tanLambdaCut,
335 mVrtParams[iteration].phiCut);
338 totalLines.local() += mTimeFrame->getLines(pivotRofId).size();
341 mTimeFrame->setNLinesTotal(totalLines.combine(std::plus<int>()));
345 debugComputeTrackletMatching(iteration);
352template <
int nLayers>
355 auto nsigmaCut{std::min(mVrtParams[iteration].vertNsigmaCut * mVrtParams[iteration].vertNsigmaCut * (mVrtParams[iteration].vertRadiusSigma * mVrtParams[iteration].vertRadiusSigma + mVrtParams[iteration].trackletSigma * mVrtParams[iteration].trackletSigma), 1.98f)};
360 for (
int rofId{0}; rofId < mTimeFrame->getNrof(); ++rofId) {
361 if (iteration && (
int)mTimeFrame->getPrimaryVertices(rofId).size() > mVrtParams[iteration].vertPerRofThreshold) {
364 const int numTracklets{
static_cast<int>(mTimeFrame->getLines(rofId).size())};
367 for (
int line1{0}; line1 < numTracklets; ++line1) {
368 if (usedTracklets[line1]) {
371 for (
int line2{line1 + 1}; line2 < numTracklets; ++line2) {
372 if (usedTracklets[line2]) {
375 auto dca{Line::getDCA(mTimeFrame->getLines(rofId)[line1], mTimeFrame->getLines(rofId)[line2])};
376 if (dca < mVrtParams[iteration].pairCut) {
377 mTimeFrame->getTrackletClusters(rofId).emplace_back(line1, mTimeFrame->getLines(rofId)[line1], line2, mTimeFrame->getLines(rofId)[line2]);
378 std::array<float, 3> tmpVertex{mTimeFrame->getTrackletClusters(rofId).back().getVertex()};
379 if (tmpVertex[0] * tmpVertex[0] + tmpVertex[1] * tmpVertex[1] > 4.f) {
380 mTimeFrame->getTrackletClusters(rofId).pop_back();
383 usedTracklets[line1] =
true;
384 usedTracklets[line2] =
true;
385 for (
int tracklet3{0}; tracklet3 < numTracklets; ++tracklet3) {
386 if (usedTracklets[tracklet3]) {
389 if (Line::getDistanceFromPoint(mTimeFrame->getLines(rofId)[tracklet3], tmpVertex) < mVrtParams[iteration].pairCut) {
390 mTimeFrame->getTrackletClusters(rofId).back().add(tracklet3, mTimeFrame->getLines(rofId)[tracklet3]);
391 usedTracklets[tracklet3] =
true;
392 tmpVertex = mTimeFrame->getTrackletClusters(rofId).back().getVertex();
399 if (mVrtParams[iteration].allowSingleContribClusters) {
400 auto beamLine =
Line{{mTimeFrame->getBeamX(), mTimeFrame->getBeamY(), -50.f}, {mTimeFrame->getBeamX(), mTimeFrame->getBeamY(), 50.f}};
401 for (
size_t iLine{0}; iLine < numTracklets; ++iLine) {
402 if (!usedTracklets[iLine]) {
403 auto dca = Line::getDCA(mTimeFrame->getLines(rofId)[iLine], beamLine);
404 if (dca < mVrtParams[iteration].pairCut) {
405 mTimeFrame->getTrackletClusters(rofId).emplace_back(iLine, mTimeFrame->getLines(rofId)[iLine], -1, beamLine);
412 std::sort(mTimeFrame->getTrackletClusters(rofId).begin(), mTimeFrame->getTrackletClusters(rofId).end(),
413 [](ClusterLines& cluster1, ClusterLines& cluster2) { return cluster1.getSize() > cluster2.getSize(); });
414 noClustersVec[rofId] =
static_cast<int>(mTimeFrame->getTrackletClusters(rofId).size());
415 for (
int iCluster1{0}; iCluster1 < noClustersVec[rofId]; ++iCluster1) {
416 std::array<float, 3> vertex1{mTimeFrame->getTrackletClusters(rofId)[iCluster1].getVertex()};
417 std::array<float, 3> vertex2{};
418 for (
int iCluster2{iCluster1 + 1}; iCluster2 < noClustersVec[rofId]; ++iCluster2) {
419 vertex2 = mTimeFrame->getTrackletClusters(rofId)[iCluster2].getVertex();
420 if (o2::gpu::GPUCommonMath::Abs(vertex1[2] - vertex2[2]) < mVrtParams[iteration].clusterCut) {
421 float distance{(vertex1[0] - vertex2[0]) * (vertex1[0] - vertex2[0]) +
422 (vertex1[1] - vertex2[1]) * (vertex1[1] - vertex2[1]) +
423 (vertex1[2] - vertex2[2]) * (vertex1[2] - vertex2[2])};
424 if (
distance < mVrtParams[iteration].pairCut * mVrtParams[iteration].pairCut) {
425 for (
auto label : mTimeFrame->getTrackletClusters(rofId)[iCluster2].getLabels()) {
426 mTimeFrame->getTrackletClusters(rofId)[iCluster1].add(
label, mTimeFrame->getLines(rofId)[
label]);
427 vertex1 = mTimeFrame->getTrackletClusters(rofId)[iCluster1].getVertex();
429 mTimeFrame->getTrackletClusters(rofId).erase(mTimeFrame->getTrackletClusters(rofId).begin() + iCluster2);
431 --noClustersVec[rofId];
437 for (
int rofId{0}; rofId < mTimeFrame->getNrof(); ++rofId) {
438 std::sort(mTimeFrame->getTrackletClusters(rofId).begin(), mTimeFrame->getTrackletClusters(rofId).end(),
439 [](
const ClusterLines& cluster1,
const ClusterLines& cluster2) { return cluster1.getSize() > cluster2.getSize(); });
440 bool atLeastOneFound{
false};
441 for (
int iCluster{0}; iCluster < noClustersVec[rofId]; ++iCluster) {
442 bool lowMultCandidate{
false};
443 double beamDistance2{(mTimeFrame->getBeamX() - mTimeFrame->getTrackletClusters(rofId)[iCluster].getVertex()[0]) * (mTimeFrame->getBeamX() - mTimeFrame->getTrackletClusters(rofId)[iCluster].getVertex()[0]) +
444 (mTimeFrame->getBeamY() - mTimeFrame->getTrackletClusters(rofId)[iCluster].getVertex()[1]) * (mTimeFrame->getBeamY() - mTimeFrame->getTrackletClusters(rofId)[iCluster].getVertex()[1])};
445 if (atLeastOneFound && (lowMultCandidate = mTimeFrame->getTrackletClusters(rofId)[iCluster].getSize() < mVrtParams[iteration].clusterContributorsCut)) {
446 lowMultCandidate &= (beamDistance2 < mVrtParams[iteration].lowMultBeamDistCut * mVrtParams[iteration].lowMultBeamDistCut);
447 if (!lowMultCandidate) {
448 mTimeFrame->getTrackletClusters(rofId).erase(mTimeFrame->getTrackletClusters(rofId).begin() + iCluster);
449 noClustersVec[rofId]--;
454 if (beamDistance2 < nsigmaCut && o2::gpu::GPUCommonMath::Abs(mTimeFrame->getTrackletClusters(rofId)[iCluster].getVertex()[2]) < mVrtParams[iteration].maxZPositionAllowed) {
455 atLeastOneFound =
true;
457 mTimeFrame->getTrackletClusters(rofId)[iCluster].getVertex()[1],
458 mTimeFrame->getTrackletClusters(rofId)[iCluster].getVertex()[2]),
459 mTimeFrame->getTrackletClusters(rofId)[iCluster].getRMS2(),
461 mTimeFrame->getTrackletClusters(rofId)[iCluster].getSize(),
462 mTimeFrame->getTrackletClusters(rofId)[iCluster].getAvgDistance2());
465 vertex.setFlags(Vertex::UPCMode);
467 vertex.setTimeStamp(mTimeFrame->getTrackletClusters(rofId)[iCluster].getROF());
468 if (mTimeFrame->hasMCinformation()) {
470 for (
auto&
index : mTimeFrame->getTrackletClusters(rofId)[iCluster].getLabels()) {
471 labels.push_back(mTimeFrame->getLinesLabel(rofId)[
index]);
473 polls.push_back(computeMain(labels));
474 if (mVrtParams[iteration].outputContLabels) {
475 contLabels.insert(contLabels.end(), labels.begin(), labels.end());
481 mTimeFrame->addPrimaryVertices(vertices, iteration);
482 if (mTimeFrame->hasMCinformation()) {
483 mTimeFrame->addPrimaryVerticesLabels(polls);
484 if (mVrtParams[iteration].outputContLabels) {
485 mTimeFrame->addPrimaryVerticesContributorLabels(contLabels);
489 mTimeFrame->addPrimaryVerticesInROF(vertices, rofId, iteration);
490 if (mTimeFrame->hasMCinformation()) {
491 mTimeFrame->addPrimaryVerticesLabelsInROF(polls, rofId);
492 if (mVrtParams[iteration].outputContLabels) {
493 mTimeFrame->addPrimaryVerticesContributorLabelsInROF(contLabels, rofId);
497 if (vertices.empty() && !(iteration && (
int)mTimeFrame->getPrimaryVertices(rofId).size() > mVrtParams[iteration].vertPerRofThreshold)) {
498 mTimeFrame->getNoVertexROF()++;
505 debugComputeVertices(iteration);
509template <
int nLayers>
512 LOGP(info,
"Using truth seeds as vertices; will skip computations");
513 mTimeFrame->resetRofPV();
515 const auto irs = dc->getEventRecords();
524 std::map<int, VertInfo> vertices;
525 for (
int iSrc{0}; iSrc < mcReader.
getNSources(); ++iSrc) {
526 auto eveId2colId = dc->getCollisionIndicesForSource(iSrc);
527 for (
int iEve{0}; iEve < mcReader.
getNEvents(iSrc); ++iEve) {
528 const auto&
ir = irs[eveId2colId[iEve]];
531 int rofId = ((
ir - raw::HBFUtils::Instance().getFirstSampledTFIR()).toLong() - roFrameBiasInBC) / roFrameLengthInBC;
532 if (!vertices.contains(rofId)) {
540 vert.setTimeStamp(rofId);
541 vert.setNContributors(std::ranges::count_if(mcReader.
getTracks(iSrc, iEve), [](
const auto& trk) {
542 return trk.isPrimary() && trk.GetPt() > 0.2 && std::abs(trk.GetEta()) < 1.3;
544 vert.setXYZ((
float)eve.GetX(), (
float)eve.GetY(), (
float)eve.GetZ());
546 constexpr float cov = 50e-9;
547 vert.setCov(cov, cov, cov, cov, cov, cov);
548 vertices[rofId].vertices.push_back(vert);
549 vertices[rofId].srcs.push_back(iSrc);
550 vertices[rofId].events.push_back(iEve);
555 for (
int iROF{0}; iROF < mTimeFrame->getNrof(); ++iROF) {
558 if (vertices.contains(iROF)) {
559 const auto& vertInfo = vertices[iROF];
560 verts = vertInfo.vertices;
561 nVerts += verts.size();
562 for (
size_t i{0};
i < verts.size(); ++
i) {
564 polls.emplace_back(lbl, 1.f);
567 mTimeFrame->getNoVertexROF()++;
569 mTimeFrame->addPrimaryVertices(verts, 0);
570 mTimeFrame->addPrimaryVerticesLabels(polls);
572 LOGP(info,
"Found {}/{} ROFs with {} vertices -> <NV>={:.2f}", vertices.size(), mTimeFrame->getNrof(), nVerts, (
float)nVerts / (
float)vertices.size());
575template <
int nLayers>
578#if defined(VTX_DEBUG)
579 LOGP(info,
"Vertexer with debug output forcing single thread");
580 mTaskArena = std::make_shared<tbb::task_arena>(1);
582 if (arena ==
nullptr) {
583 mTaskArena = std::make_shared<tbb::task_arena>(std::abs(
n));
584 LOGP(info,
"Setting seeding vertexer with {} threads.",
n);
587 LOGP(info,
"Attaching vertexer to calling thread's arena");
592template <
int nLayers>
596 LOGP(info,
"writing debug output for computeTracklets");
597 for (
int rofId{0}; rofId < mTimeFrame->getNrof(); ++rofId) {
598 const auto& strk0 = mTimeFrame->getFoundTracklets(rofId, 0);
599 std::vector<Tracklet> trk0(strk0.begin(), strk0.end());
600 const auto& strk1 = mTimeFrame->getFoundTracklets(rofId, 1);
601 std::vector<Tracklet> trk1(strk1.begin(), strk1.end());
602 (*stream) <<
"tracklets"
603 <<
"Tracklets0=" << trk0
604 <<
"Tracklets1=" << trk1
605 <<
"iteration=" << iteration
612template <
int nLayers>
613void VertexerTraits<nLayers>::debugComputeTrackletMatching(
int iteration)
615 auto stream =
new utils::TreeStreamRedirector(
"artefacts_tf.root",
"update");
616 LOGP(info,
"writing debug output for computeTrackletMatching");
617 for (
int rofId{0}; rofId < mTimeFrame->getNrof(); ++rofId) {
619 <<
"Lines=" <<
toSTDVector(mTimeFrame->getLines(rofId))
620 <<
"NTrackletCluster01=" << mTimeFrame->getNTrackletsCluster(rofId, 0)
621 <<
"NTrackletCluster12=" << mTimeFrame->getNTrackletsCluster(rofId, 1)
622 <<
"iteration=" << iteration
626 if (mTimeFrame->hasMCinformation()) {
627 LOGP(info,
"\tdumping also MC information");
629 const auto irs = dc->getEventRecords();
634 std::map<int, int> eve2BcInROF, bcInRofNEve;
635 for (
int iSrc{0}; iSrc < mcReader.getNSources(); ++iSrc) {
636 auto eveId2colId = dc->getCollisionIndicesForSource(iSrc);
637 for (
int iEve{0}; iEve < mcReader.getNEvents(iSrc); ++iEve) {
638 const auto&
ir = irs[eveId2colId[iEve]];
640 const auto& eve = mcReader.getMCEventHeader(iSrc, iEve);
641 const int bcInROF = ((
ir - raw::HBFUtils::Instance().getFirstSampledTFIR()).toLong() - roFrameBiasInBC) % roFrameLengthInBC;
642 eve2BcInROF[iEve] = bcInROF;
643 ++bcInRofNEve[bcInROF];
648 std::unordered_map<int, int> bcROFNTracklets01, bcROFNTracklets12;
649 std::vector<std::vector<int>> tracklet01BC, tracklet12BC;
650 for (
int rofId{0}; rofId < mTimeFrame->getNrof(); ++rofId) {
652 const auto& tracklet01 = mTimeFrame->getFoundTracklets(rofId, 0);
653 const auto& lbls01 = mTimeFrame->getLabelsFoundTracklets(rofId, 0);
654 auto& trkls01 = tracklet01BC.emplace_back();
655 for (
int iTrklt{0}; iTrklt < (
int)tracklet01.size(); ++iTrklt) {
656 const auto& tracklet = tracklet01[iTrklt];
657 const auto& lbl = lbls01[iTrklt];
658 if (lbl.isCorrect()) {
659 ++bcROFNTracklets01[eve2BcInROF[lbl.getEventID()]];
660 trkls01.push_back(eve2BcInROF[lbl.getEventID()]);
662 trkls01.push_back(-1);
667 const auto& tracklet12 = mTimeFrame->getFoundTracklets(rofId, 1);
668 auto& trkls12 = tracklet12BC.emplace_back();
669 for (
int iTrklt{0}; iTrklt < (
int)tracklet12.size(); ++iTrklt) {
670 const auto& tracklet = tracklet12[iTrklt];
673 int sortedId1{mTimeFrame->getSortedIndex(tracklet.rof[0], 1, tracklet.firstClusterIndex)};
674 int sortedId2{mTimeFrame->getSortedIndex(tracklet.rof[1], 2, tracklet.secondClusterIndex)};
675 for (
const auto& lab1 : mTimeFrame->getClusterLabels(1, mTimeFrame->
getClusters()[1][sortedId1].clusterId)) {
676 for (
const auto& lab2 : mTimeFrame->getClusterLabels(2, mTimeFrame->
getClusters()[2][sortedId2].clusterId)) {
677 if (lab1 == lab2 && lab1.isValid()) {
682 if (
label.isValid()) {
687 if (
label.isCorrect()) {
688 ++bcROFNTracklets12[eve2BcInROF[
label.getEventID()]];
689 trkls12.push_back(eve2BcInROF[
label.getEventID()]);
691 trkls12.push_back(-1);
696 LOGP(info,
"\tdumping ntracklets/RofBC ({})", bcInRofNEve.size());
697 for (
const auto& [bcInRof, neve] : bcInRofNEve) {
698 (*stream) <<
"ntracklets"
699 <<
"bcInROF=" << bcInRof
700 <<
"ntrkl01=" << bcROFNTracklets01[bcInRof]
701 <<
"ntrkl12=" << bcROFNTracklets12[bcInRof]
703 <<
"iteration=" << iteration
707 std::unordered_map<int, int> bcROFNLines;
708 for (
int rofId{0}; rofId < mTimeFrame->getNrof(); ++rofId) {
709 const auto& lines = mTimeFrame->getLines(rofId);
710 const auto& lbls = mTimeFrame->getLinesLabel(rofId);
711 for (
int iLine{0}; iLine < (
int)lines.size(); ++iLine) {
712 const auto&
line = lines[iLine];
713 const auto& lbl = lbls[iLine];
714 if (lbl.isCorrect()) {
715 ++bcROFNLines[eve2BcInROF[lbl.getEventID()]];
720 LOGP(info,
"\tdumping nlines/RofBC");
721 for (
const auto& [bcInRof, neve] : bcInRofNEve) {
722 (*stream) <<
"nlines"
723 <<
"bcInROF=" << bcInRof
724 <<
"nline=" << bcROFNLines[bcInRof]
726 <<
"iteration=" << iteration
734template <
int nLayers>
735void VertexerTraits<nLayers>::debugComputeVertices(
int iteration)
737 auto stream =
new utils::TreeStreamRedirector(
"artefacts_tf.root",
"update");
738 LOGP(info,
"writing debug output for computeVertices");
739 for (
auto rofId{0}; rofId < mTimeFrame->getNrof(); ++rofId) {
740 (*stream) <<
"clusterlines"
741 <<
"clines_post=" <<
toSTDVector(mTimeFrame->getTrackletClusters(rofId))
742 <<
"iteration=" << iteration
746 if (mTimeFrame->hasMCinformation()) {
747 LOGP(info,
"\tdumping also MC information");
749 const auto irs = dc->getEventRecords();
754 std::map<int, int> eve2BcInROF, bcInRofNEve;
755 for (
int iSrc{0}; iSrc < mcReader.getNSources(); ++iSrc) {
756 auto eveId2colId = dc->getCollisionIndicesForSource(iSrc);
757 for (
int iEve{0}; iEve < mcReader.getNEvents(iSrc); ++iEve) {
758 const auto&
ir = irs[eveId2colId[iEve]];
760 const auto& eve = mcReader.getMCEventHeader(iSrc, iEve);
761 const int bcInROF = ((
ir - raw::HBFUtils::Instance().getFirstSampledTFIR()).toLong() - roFrameBiasInBC) % roFrameLengthInBC;
762 eve2BcInROF[iEve] = bcInROF;
763 ++bcInRofNEve[bcInROF];
768 std::unordered_map<int, int> bcROFNVtx;
769 std::unordered_map<int, float> bcROFNPur;
770 std::unordered_map<o2::MCCompLabel, size_t> uniqueVertices;
771 for (
int rofId{0}; rofId < mTimeFrame->getNrof(); ++rofId) {
772 const auto& pvs = mTimeFrame->getPrimaryVertices(rofId);
773 const auto& lblspv = mTimeFrame->getPrimaryVerticesMCRecInfo(rofId);
774 for (
int i{0};
i < (
int)pvs.size(); ++
i) {
775 const auto& pv = pvs[
i];
776 const auto& [lbl, pur] = lblspv[
i];
777 if (lbl.isCorrect()) {
778 ++uniqueVertices[lbl];
779 ++bcROFNVtx[eve2BcInROF[lbl.getEventID()]];
780 bcROFNPur[eve2BcInROF[lbl.getEventID()]] += pur;
785 std::unordered_map<int, int> bcROFNUVtx, bcROFNCVtx;
786 for (
const auto& [k, _] : eve2BcInROF) {
787 bcROFNUVtx[k] = bcROFNCVtx[k] = 0;
790 for (
const auto& [lbl,
c] : uniqueVertices) {
792 ++bcROFNUVtx[eve2BcInROF[lbl.getEventID()]];
794 ++bcROFNCVtx[eve2BcInROF[lbl.getEventID()]];
798 LOGP(info,
"\tdumping nvtx/RofBC");
799 for (
const auto& [bcInRof, neve] : bcInRofNEve) {
801 <<
"bcInROF=" << bcInRof
802 <<
"nvtx=" << bcROFNVtx[bcInRof]
803 <<
"nuvtx=" << bcROFNUVtx[bcInRof]
804 <<
"ncvtx=" << bcROFNCVtx[bcInRof]
805 <<
"npur=" << bcROFNPur[bcInRof]
807 <<
"iteration=" << iteration
812 std::unordered_map<o2::MCCompLabel, std::vector<Vertex>> cVtx;
813 for (
int rofId{0}; rofId < mTimeFrame->getNrof(); ++rofId) {
814 const auto& pvs = mTimeFrame->getPrimaryVertices(rofId);
815 const auto& lblspv = mTimeFrame->getPrimaryVerticesMCRecInfo(rofId);
816 for (
int i{0};
i < (
int)pvs.size(); ++
i) {
817 const auto& pv = pvs[
i];
818 const auto& [lbl, pur] = lblspv[
i];
819 if (lbl.isCorrect() && uniqueVertices.contains(lbl) && uniqueVertices[lbl] > 1) {
820 if (!cVtx.contains(lbl)) {
821 cVtx[lbl] = std::vector<Vertex>();
823 cVtx[lbl].push_back(pv);
828 for (
auto& [_, vertices] : cVtx) {
829 std::sort(vertices.begin(), vertices.end(), [](
const Vertex&
a,
const Vertex&
b) { return a.getNContributors() > b.getNContributors(); });
830 for (
int i{0};
i < (
int)vertices.size(); ++
i) {
831 const auto vtx = vertices[
i];
835 <<
"dx=" << vertices[0].getX() - vtx.getX()
836 <<
"dy=" << vertices[0].getY() - vtx.getY()
837 <<
"dz=" << vertices[0].getZ() - vtx.getZ()
838 <<
"drof=" << vertices[0].getTimeStamp().getTimeStamp() - vtx.getTimeStamp().getTimeStamp()
839 <<
"dnc=" << vertices[0].getNContributors() - vtx.getNContributors()
840 <<
"iteration=" << iteration
849template class VertexerTraits<7>;
Class to compute the primary vertex in ITS from tracklets.
static constexpr int maxTrackID()
static const DPLAlpideParam< N > & Instance()
HMPID cluster implementation.
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
size_t getNSources() const
Get number of sources.
std::vector< MCTrack > const & getTracks(int source, int event) const
variant returning all tracks for source and event at once
GLboolean GLboolean GLboolean b
GLsizei GLsizei GLfloat distance
GLuint GLsizei const GLchar * label
GLboolean GLboolean GLboolean GLboolean a
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.
o2::InteractionRecord ir(0, 0)
std::vector< Tracklet64 > tracklets