Project
Loading...
Searching...
No Matches
TrackerTraits.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.
15
16#include <algorithm>
17#include <iostream>
18#include <iterator>
19#include <ranges>
20#include <type_traits>
21
22#ifdef OPTIMISATION_OUTPUT
23#include <format>
24#include <fstream>
25#endif
26
27#include <oneapi/tbb/blocked_range.h>
28#include <oneapi/tbb/parallel_sort.h>
29
32#include "GPUCommonMath.h"
33#include "ITStracking/Cell.h"
40
42
43namespace o2::its
44{
45
46struct PassMode {
47 using OnePass = std::integral_constant<int, 0>;
48 using TwoPassCount = std::integral_constant<int, 1>;
49 using TwoPassInsert = std::integral_constant<int, 2>;
50};
51
52template <int nLayers>
53void TrackerTraits<nLayers>::computeLayerTracklets(const int iteration, int iROFslice, int iVertex)
55#ifdef OPTIMISATION_OUTPUT
56 static int iter{0};
57 std::ofstream off(std::format("tracklets{}.txt", iter++));
58#endif
59
60 for (int iLayer = 0; iLayer < mTrkParams[iteration].TrackletsPerRoad(); ++iLayer) {
61 mTimeFrame->getTracklets()[iLayer].clear();
62 mTimeFrame->getTrackletsLabel(iLayer).clear();
63 if (iLayer > 0) {
64 std::fill(mTimeFrame->getTrackletsLookupTable()[iLayer - 1].begin(),
65 mTimeFrame->getTrackletsLookupTable()[iLayer - 1].end(), 0);
66 }
67 }
69 const Vertex diamondVert({mTrkParams[iteration].Diamond[0], mTrkParams[iteration].Diamond[1], mTrkParams[iteration].Diamond[2]}, {25.e-6f, 0.f, 0.f, 25.e-6f, 0.f, 36.f}, 1, 1.f);
70 gsl::span<const Vertex> diamondSpan(&diamondVert, 1);
71 int startROF{mTrkParams[iteration].nROFsPerIterations > 0 ? iROFslice * mTrkParams[iteration].nROFsPerIterations : 0};
72 int endROF{o2::gpu::GPUCommonMath::Min(mTrkParams[iteration].nROFsPerIterations > 0 ? (iROFslice + 1) * mTrkParams[iteration].nROFsPerIterations + mTrkParams[iteration].DeltaROF : mTimeFrame->getNrof(), mTimeFrame->getNrof())};
73
74 mTaskArena->execute([&] {
75 auto forTracklets = [&](auto Tag, int iLayer, int pivotROF, int base, int& offset) -> int {
76 if (!mTimeFrame->mMultiplicityCutMask[pivotROF]) {
77 return 0;
78 }
79 int minROF = o2::gpu::CAMath::Max(startROF, pivotROF - mTrkParams[iteration].DeltaROF);
80 int maxROF = o2::gpu::CAMath::Min(endROF - 1, pivotROF + mTrkParams[iteration].DeltaROF);
81 gsl::span<const Vertex> primaryVertices = mTrkParams[iteration].UseDiamond ? diamondSpan : mTimeFrame->getPrimaryVertices(minROF, maxROF);
82 if (primaryVertices.empty()) {
83 return 0;
84 }
85 const int startVtx = iVertex >= 0 ? iVertex : 0;
86 const int endVtx = iVertex >= 0 ? o2::gpu::CAMath::Min(iVertex + 1, int(primaryVertices.size())) : int(primaryVertices.size());
87 if (endVtx <= startVtx) {
88 return 0;
89 }
90
91 int localCount = 0;
92 auto& tracklets = mTimeFrame->getTracklets()[iLayer];
93 auto layer0 = mTimeFrame->getClustersOnLayer(pivotROF, iLayer);
94 if (layer0.empty()) {
95 return 0;
96 }
97
98 const float meanDeltaR = mTrkParams[iteration].LayerRadii[iLayer + 1] - mTrkParams[iteration].LayerRadii[iLayer];
99
100 for (int iCluster = 0; iCluster < int(layer0.size()); ++iCluster) {
101 const Cluster& currentCluster = layer0[iCluster];
102 const int currentSortedIndex = mTimeFrame->getSortedIndex(pivotROF, iLayer, iCluster);
103 if (mTimeFrame->isClusterUsed(iLayer, currentCluster.clusterId)) {
104 continue;
105 }
106 const float inverseR0 = 1.f / currentCluster.radius;
107
108 for (int iV = startVtx; iV < endVtx; ++iV) {
109 const auto& pv = primaryVertices[iV];
110 if ((pv.isFlagSet(Vertex::Flags::UPCMode) && iteration != 3) || (iteration == 3 && !pv.isFlagSet(Vertex::Flags::UPCMode))) {
111 continue;
112 }
113
114 const float resolution = o2::gpu::CAMath::Sqrt(math_utils::Sq(mTimeFrame->getPositionResolution(iLayer)) + math_utils::Sq(mTrkParams[iteration].PVres) / float(pv.getNContributors()));
115 const float tanLambda = (currentCluster.zCoordinate - pv.getZ()) * inverseR0;
116 const float zAtRmin = tanLambda * (mTimeFrame->getMinR(iLayer + 1) - currentCluster.radius) + currentCluster.zCoordinate;
117 const float zAtRmax = tanLambda * (mTimeFrame->getMaxR(iLayer + 1) - currentCluster.radius) + currentCluster.zCoordinate;
118 const float sqInvDeltaZ0 = 1.f / (math_utils::Sq(currentCluster.zCoordinate - pv.getZ()) + constants::Tolerance);
119 const float sigmaZ = o2::gpu::CAMath::Sqrt(
120 math_utils::Sq(resolution) * math_utils::Sq(tanLambda) * ((math_utils::Sq(inverseR0) + sqInvDeltaZ0) * math_utils::Sq(meanDeltaR) + 1.f) + math_utils::Sq(meanDeltaR * mTimeFrame->getMSangle(iLayer)));
121
122 auto bins = getBinsRect(currentCluster, iLayer + 1, zAtRmin, zAtRmax, sigmaZ * mTrkParams[iteration].NSigmaCut, mTimeFrame->getPhiCut(iLayer));
123 if (bins.x == 0 && bins.y == 0 && bins.z == 0 && bins.w == 0) {
124 continue;
125 }
126 int phiBinsNum = bins.w - bins.y + 1;
127 if (phiBinsNum < 0) {
128 phiBinsNum += mTrkParams[iteration].PhiBins;
129 }
130
131 for (int targetROF{minROF}; targetROF <= maxROF; ++targetROF) {
132 if (!mTimeFrame->mMultiplicityCutMask[targetROF]) {
133 continue;
134 }
135 auto layer1 = mTimeFrame->getClustersOnLayer(targetROF, iLayer + 1);
136 if (layer1.empty()) {
137 continue;
138 }
139 for (int iPhi = 0; iPhi < phiBinsNum; ++iPhi) {
140 const int iPhiBin = (bins.y + iPhi) % mTrkParams[iteration].PhiBins;
141 const int firstBinIdx = mTimeFrame->mIndexTableUtils.getBinIndex(bins.x, iPhiBin);
142 const int maxBinIdx = firstBinIdx + (bins.z - bins.x) + 1;
143 const int firstRow = mTimeFrame->getIndexTable(targetROF, iLayer + 1)[firstBinIdx];
144 const int lastRow = mTimeFrame->getIndexTable(targetROF, iLayer + 1)[maxBinIdx];
145 for (int iNext = firstRow; iNext < lastRow; ++iNext) {
146 if (iNext >= int(layer1.size())) {
147 break;
148 }
149 const Cluster& nextCluster = layer1[iNext];
150 if (mTimeFrame->isClusterUsed(iLayer + 1, nextCluster.clusterId)) {
151 continue;
152 }
153 float deltaPhi = o2::gpu::GPUCommonMath::Abs(currentCluster.phi - nextCluster.phi);
154 float deltaZ = o2::gpu::GPUCommonMath::Abs((tanLambda * (nextCluster.radius - currentCluster.radius)) + currentCluster.zCoordinate - nextCluster.zCoordinate);
155
156#ifdef OPTIMISATION_OUTPUT
157 MCCompLabel label;
158 int currentId{currentCluster.clusterId};
159 int nextId{nextCluster.clusterId};
160 for (auto& lab1 : mTimeFrame->getClusterLabels(iLayer, currentId)) {
161 for (auto& lab2 : mTimeFrame->getClusterLabels(iLayer + 1, nextId)) {
162 if (lab1 == lab2 && lab1.isValid()) {
163 label = lab1;
164 break;
165 }
166 }
167 if (label.isValid()) {
168 break;
169 }
170 }
171 off << std::format("{}\t{:d}\t{}\t{}\t{}\t{}", iLayer, label.isValid(), (tanLambda * (nextCluster.radius - currentCluster.radius) + currentCluster.zCoordinate - nextCluster.zCoordinate) / sigmaZ, tanLambda, resolution, sigmaZ) << std::endl;
172#endif
173
174 if (deltaZ / sigmaZ < mTrkParams[iteration].NSigmaCut &&
175 ((deltaPhi < mTimeFrame->getPhiCut(iLayer) || o2::gpu::GPUCommonMath::Abs(deltaPhi - o2::constants::math::TwoPI) < mTimeFrame->getPhiCut(iLayer)))) {
176 const float phi{o2::gpu::CAMath::ATan2(currentCluster.yCoordinate - nextCluster.yCoordinate, currentCluster.xCoordinate - nextCluster.xCoordinate)};
177 const float tanL = (currentCluster.zCoordinate - nextCluster.zCoordinate) / (currentCluster.radius - nextCluster.radius);
178 if constexpr (decltype(Tag)::value == PassMode::OnePass::value) {
179 tracklets.emplace_back(currentSortedIndex, mTimeFrame->getSortedIndex(targetROF, iLayer + 1, iNext), tanL, phi, pivotROF, targetROF);
180 } else if constexpr (decltype(Tag)::value == PassMode::TwoPassCount::value) {
181 ++localCount;
182 } else if constexpr (decltype(Tag)::value == PassMode::TwoPassInsert::value) {
183 const int idx = base + offset++;
184 tracklets[idx] = Tracklet(currentSortedIndex, mTimeFrame->getSortedIndex(targetROF, iLayer + 1, iNext), tanL, phi, pivotROF, targetROF);
185 }
186 }
187 }
188 }
189 }
190 }
191 }
192 return localCount;
193 };
194
195 int dummy{0};
196 if (mTaskArena->max_concurrency() <= 1) {
197 for (int pivotROF{startROF}; pivotROF < endROF; ++pivotROF) {
198 for (int iLayer{0}; iLayer < mTrkParams[iteration].TrackletsPerRoad(); ++iLayer) {
199 forTracklets(PassMode::OnePass{}, iLayer, pivotROF, 0, dummy);
200 }
201 }
202 } else {
203 bounded_vector<bounded_vector<int>> perROFCount(mTrkParams[iteration].TrackletsPerRoad(), bounded_vector<int>(endROF - startROF + 1, 0, mMemoryPool.get()), mMemoryPool.get());
204 tbb::parallel_for(
205 tbb::blocked_range2d<int, int>(0, mTrkParams[iteration].TrackletsPerRoad(), 1,
206 startROF, endROF, 1),
207 [&](auto const& Range) {
208 for (int iLayer{Range.rows().begin()}; iLayer < Range.rows().end(); ++iLayer) {
209 for (int pivotROF = Range.cols().begin(); pivotROF < Range.cols().end(); ++pivotROF) {
210 perROFCount[iLayer][pivotROF - startROF] = forTracklets(PassMode::TwoPassCount{}, iLayer, pivotROF, 0, dummy);
211 }
212 }
213 });
214
215 tbb::parallel_for(
216 tbb::blocked_range<int>(0, mTrkParams[iteration].TrackletsPerRoad()),
217 [&](auto const& Layers) {
218 for (int iLayer{Layers.begin()}; iLayer < Layers.end(); ++iLayer) {
219 std::exclusive_scan(perROFCount[iLayer].begin(), perROFCount[iLayer].end(), perROFCount[iLayer].begin(), 0);
220 mTimeFrame->getTracklets()[iLayer].resize(perROFCount[iLayer].back());
221 }
222 });
223
224 tbb::parallel_for(
225 tbb::blocked_range2d<int, int>(0, mTrkParams[iteration].TrackletsPerRoad(), 1,
226 startROF, endROF, 1),
227 [&](auto const& Range) {
228 for (int iLayer{Range.rows().begin()}; iLayer < Range.rows().end(); ++iLayer) {
229 if (perROFCount[iLayer].back() == 0) {
230 continue;
231 }
232 for (int pivotROF = Range.cols().begin(); pivotROF < Range.cols().end(); ++pivotROF) {
233 int baseIdx = perROFCount[iLayer][pivotROF - startROF];
234 if (baseIdx == perROFCount[iLayer][pivotROF - startROF + 1]) {
235 continue;
236 }
237 int localIdx = 0;
238 forTracklets(PassMode::TwoPassInsert{}, iLayer, pivotROF, baseIdx, localIdx);
239 }
240 }
241 });
242 }
243
244 tbb::parallel_for(
245 tbb::blocked_range<int>(0, mTrkParams[iteration].TrackletsPerRoad()),
246 [&](const tbb::blocked_range<int>& Layers) {
247 for (int iLayer = Layers.begin(); iLayer < Layers.end(); ++iLayer) {
249 auto& trkl{mTimeFrame->getTracklets()[iLayer]};
250 tbb::parallel_sort(trkl.begin(), trkl.end(), [](const Tracklet& a, const Tracklet& b) -> bool {
251 if (a.firstClusterIndex != b.firstClusterIndex) {
252 return a.firstClusterIndex < b.firstClusterIndex;
253 }
254 return a.secondClusterIndex < b.secondClusterIndex;
255 });
257 trkl.erase(std::unique(trkl.begin(), trkl.end(), [](const Tracklet& a, const Tracklet& b) -> bool {
258 return a.firstClusterIndex == b.firstClusterIndex && a.secondClusterIndex == b.secondClusterIndex;
259 }),
260 trkl.end());
261 trkl.shrink_to_fit();
262 if (iLayer > 0) {
263 auto& lut{mTimeFrame->getTrackletsLookupTable()[iLayer - 1]};
264 if (!trkl.empty()) {
265 for (const auto& tkl : trkl) {
266 lut[tkl.firstClusterIndex + 1]++;
267 }
268 std::inclusive_scan(lut.begin(), lut.end(), lut.begin());
269 }
270 }
271 }
272 });
273
275 if (mTimeFrame->hasMCinformation() && mTrkParams[iteration].createArtefactLabels) {
276 tbb::parallel_for(
277 tbb::blocked_range<int>(0, mTrkParams[iteration].TrackletsPerRoad()),
278 [&](const tbb::blocked_range<int>& Layers) {
279 for (int iLayer = Layers.begin(); iLayer < Layers.end(); ++iLayer) {
280 for (auto& trk : mTimeFrame->getTracklets()[iLayer]) {
281 MCCompLabel label;
282 int currentId{mTimeFrame->getClusters()[iLayer][trk.firstClusterIndex].clusterId};
283 int nextId{mTimeFrame->getClusters()[iLayer + 1][trk.secondClusterIndex].clusterId};
284 for (const auto& lab1 : mTimeFrame->getClusterLabels(iLayer, currentId)) {
285 for (const auto& lab2 : mTimeFrame->getClusterLabels(iLayer + 1, nextId)) {
286 if (lab1 == lab2 && lab1.isValid()) {
287 label = lab1;
288 break;
289 }
290 }
291 if (label.isValid()) {
292 break;
293 }
294 }
295 mTimeFrame->getTrackletsLabel(iLayer).emplace_back(label);
296 }
297 }
298 });
299 }
300 });
301} // namespace o2::its
302
303template <int nLayers>
305{
306#ifdef OPTIMISATION_OUTPUT
307 static int iter{0};
308 std::ofstream off(std::format("cells{}.txt", iter++));
309#endif
310
311 for (int iLayer = 0; iLayer < mTrkParams[iteration].CellsPerRoad(); ++iLayer) {
312 deepVectorClear(mTimeFrame->getCells()[iLayer]);
313 if (iLayer > 0) {
314 deepVectorClear(mTimeFrame->getCellsLookupTable()[iLayer - 1]);
315 }
316 if (mTimeFrame->hasMCinformation() && mTrkParams[iteration].createArtefactLabels) {
317 deepVectorClear(mTimeFrame->getCellsLabel(iLayer));
318 }
319 }
320
321 mTaskArena->execute([&] {
322 auto forTrackletCells = [&](auto Tag, int iLayer, bounded_vector<CellSeedN>& layerCells, int iTracklet, int offset = 0) -> int {
323 const Tracklet& currentTracklet{mTimeFrame->getTracklets()[iLayer][iTracklet]};
324 const int nextLayerClusterIndex{currentTracklet.secondClusterIndex};
325 const int nextLayerFirstTrackletIndex{mTimeFrame->getTrackletsLookupTable()[iLayer][nextLayerClusterIndex]};
326 const int nextLayerLastTrackletIndex{mTimeFrame->getTrackletsLookupTable()[iLayer][nextLayerClusterIndex + 1]};
327 int foundCells{0};
328 for (int iNextTracklet{nextLayerFirstTrackletIndex}; iNextTracklet < nextLayerLastTrackletIndex; ++iNextTracklet) {
329 const Tracklet& nextTracklet{mTimeFrame->getTracklets()[iLayer + 1][iNextTracklet]};
330 const auto& nextLbl = mTimeFrame->getTrackletsLabel(iLayer + 1)[iNextTracklet];
331 if (mTimeFrame->getTracklets()[iLayer + 1][iNextTracklet].firstClusterIndex != nextLayerClusterIndex) {
332 break;
333 }
334 if (mTrkParams[iteration].DeltaROF && currentTracklet.getSpanRof(nextTracklet) > mTrkParams[iteration].DeltaROF) { // TODO this has to be improved for the staggering
335 continue;
336 }
337 const float deltaTanLambda{std::abs(currentTracklet.tanLambda - nextTracklet.tanLambda)};
338
339#ifdef OPTIMISATION_OUTPUT
340 float resolution{o2::gpu::CAMath::Sqrt(0.5f * (mTrkParams[iteration].SystErrorZ2[iLayer] + mTrkParams[iteration].SystErrorZ2[iLayer + 1] + mTrkParams[iteration].SystErrorZ2[iLayer + 2] + mTrkParams[iteration].SystErrorY2[iLayer] + mTrkParams[iteration].SystErrorY2[iLayer + 1] + mTrkParams[iteration].SystErrorY2[iLayer + 2])) / mTrkParams[iteration].LayerResolution[iLayer]};
341 resolution = resolution > 1.e-12 ? resolution : 1.f;
342 bool good{mTimeFrame->getTrackletsLabel(iLayer)[iTracklet] == mTimeFrame->getTrackletsLabel(iLayer + 1)[iNextTracklet]};
343 float signedDelta{currentTracklet.tanLambda - nextTracklet.tanLambda};
344 off << std::format("{}\t{:d}\t{}\t{}\t{}\t{}", iLayer, good, signedDelta, signedDelta / (mTrkParams[iteration].CellDeltaTanLambdaSigma), tanLambda, resolution) << std::endl;
345#endif
346
347 if (deltaTanLambda / mTrkParams[iteration].CellDeltaTanLambdaSigma < mTrkParams[iteration].NSigmaCut) {
348
350 const int clusId[3]{
351 mTimeFrame->getClusters()[iLayer][currentTracklet.firstClusterIndex].clusterId,
352 mTimeFrame->getClusters()[iLayer + 1][nextTracklet.firstClusterIndex].clusterId,
353 mTimeFrame->getClusters()[iLayer + 2][nextTracklet.secondClusterIndex].clusterId};
354 const auto& cluster1_glo = mTimeFrame->getUnsortedClusters()[iLayer][clusId[0]];
355 const auto& cluster2_glo = mTimeFrame->getUnsortedClusters()[iLayer + 1][clusId[1]];
356 const auto& cluster3_tf = mTimeFrame->getTrackingFrameInfoOnLayer(iLayer + 2)[clusId[2]];
357 auto track{buildTrackSeed(cluster1_glo, cluster2_glo, cluster3_tf)};
358
359 float chi2{0.f};
360 bool good{false};
361 for (int iC{2}; iC--;) {
362 const TrackingFrameInfo& trackingHit = mTimeFrame->getTrackingFrameInfoOnLayer(iLayer + iC)[clusId[iC]];
363
364 if (!track.rotate(trackingHit.alphaTrackingFrame)) {
365 break;
366 }
367
368 if (!track.propagateTo(trackingHit.xTrackingFrame, getBz())) {
369 break;
370 }
371
372 if (!track.correctForMaterial(mTrkParams[0].LayerxX0[iLayer + iC], mTrkParams[0].LayerxX0[iLayer] * constants::Radl * constants::Rho, true)) {
373 break;
374 }
375
376 const auto predChi2{track.getPredictedChi2Quiet(trackingHit.positionTrackingFrame, trackingHit.covarianceTrackingFrame)};
377 if (!iC && predChi2 > mTrkParams[iteration].MaxChi2ClusterAttachment) {
378 break;
379 }
380
381 if (!track.o2::track::TrackParCov::update(trackingHit.positionTrackingFrame, trackingHit.covarianceTrackingFrame)) {
382 break;
383 }
384
385 good = !iC;
386 chi2 += predChi2;
387 }
388 if (good) {
389 if constexpr (decltype(Tag)::value == PassMode::OnePass::value) {
390 layerCells.emplace_back(iLayer, clusId[0], clusId[1], clusId[2], iTracklet, iNextTracklet, track, chi2);
391 ++foundCells;
392 } else if constexpr (decltype(Tag)::value == PassMode::TwoPassCount::value) {
393 ++foundCells;
394 } else if constexpr (decltype(Tag)::value == PassMode::TwoPassInsert::value) {
395 layerCells[offset++] = CellSeedN(iLayer, clusId[0], clusId[1], clusId[2], iTracklet, iNextTracklet, track, chi2);
396 } else {
397 static_assert(false, "Unknown mode!");
398 }
399 }
400 }
401 }
402 return foundCells;
403 };
404
405 tbb::parallel_for(
406 tbb::blocked_range<int>(0, mTrkParams[iteration].CellsPerRoad()),
407 [&](const tbb::blocked_range<int>& Layers) {
408 for (int iLayer = Layers.begin(); iLayer < Layers.end(); ++iLayer) {
409 if (mTimeFrame->getTracklets()[iLayer + 1].empty() ||
410 mTimeFrame->getTracklets()[iLayer].empty()) {
411 continue;
412 }
413
414 auto& layerCells = mTimeFrame->getCells()[iLayer];
415 const int currentLayerTrackletsNum{static_cast<int>(mTimeFrame->getTracklets()[iLayer].size())};
416 bounded_vector<int> perTrackletCount(currentLayerTrackletsNum + 1, 0, mMemoryPool.get());
417 if (mTaskArena->max_concurrency() <= 1) {
418 for (int iTracklet{0}; iTracklet < currentLayerTrackletsNum; ++iTracklet) {
419 perTrackletCount[iTracklet] = forTrackletCells(PassMode::OnePass{}, iLayer, layerCells, iTracklet);
420 }
421 std::exclusive_scan(perTrackletCount.begin(), perTrackletCount.end(), perTrackletCount.begin(), 0);
422 } else {
423 tbb::parallel_for(
424 tbb::blocked_range<int>(0, currentLayerTrackletsNum),
425 [&](const tbb::blocked_range<int>& Tracklets) {
426 for (int iTracklet = Tracklets.begin(); iTracklet < Tracklets.end(); ++iTracklet) {
427 perTrackletCount[iTracklet] = forTrackletCells(PassMode::TwoPassCount{}, iLayer, layerCells, iTracklet);
428 }
429 });
430
431 std::exclusive_scan(perTrackletCount.begin(), perTrackletCount.end(), perTrackletCount.begin(), 0);
432 auto totalCells{perTrackletCount.back()};
433 if (totalCells == 0) {
434 continue;
435 }
436 layerCells.resize(totalCells);
437
438 tbb::parallel_for(
439 tbb::blocked_range<int>(0, currentLayerTrackletsNum),
440 [&](const tbb::blocked_range<int>& Tracklets) {
441 for (int iTracklet = Tracklets.begin(); iTracklet < Tracklets.end(); ++iTracklet) {
442 int offset = perTrackletCount[iTracklet];
443 if (offset == perTrackletCount[iTracklet + 1]) {
444 continue;
445 }
446 forTrackletCells(PassMode::TwoPassInsert{}, iLayer, layerCells, iTracklet, offset);
447 }
448 });
449 }
450
451 if (iLayer > 0) {
452 auto& lut = mTimeFrame->getCellsLookupTable()[iLayer - 1];
453 lut.resize(currentLayerTrackletsNum + 1);
454 std::copy_n(perTrackletCount.begin(), currentLayerTrackletsNum + 1, lut.begin());
455 }
456 }
457 });
458 });
459
461 if (mTimeFrame->hasMCinformation() && mTrkParams[iteration].createArtefactLabels) {
462 tbb::parallel_for(
463 tbb::blocked_range<int>(0, mTrkParams[iteration].CellsPerRoad()),
464 [&](const tbb::blocked_range<int>& Layers) {
465 for (int iLayer = Layers.begin(); iLayer < Layers.end(); ++iLayer) {
466 mTimeFrame->getCellsLabel(iLayer).reserve(mTimeFrame->getCells()[iLayer].size());
467 for (const auto& cell : mTimeFrame->getCells()[iLayer]) {
468 MCCompLabel currentLab{mTimeFrame->getTrackletsLabel(iLayer)[cell.getFirstTrackletIndex()]};
469 MCCompLabel nextLab{mTimeFrame->getTrackletsLabel(iLayer + 1)[cell.getSecondTrackletIndex()]};
470 mTimeFrame->getCellsLabel(iLayer).emplace_back(currentLab == nextLab ? currentLab : MCCompLabel());
471 }
472 }
473 });
474 }
475}
476
477template <int nLayers>
479{
480#ifdef OPTIMISATION_OUTPUT
481 std::ofstream off(std::format("cellneighs{}.txt", iteration));
482#endif
483
484 struct Neighbor {
485 int cell{-1}, nextCell{-1}, level{-1};
486 };
487
488 mTaskArena->execute([&] {
489 for (int iLayer{0}; iLayer < mTrkParams[iteration].NeighboursPerRoad(); ++iLayer) {
490 deepVectorClear(mTimeFrame->getCellsNeighbours()[iLayer]);
491 deepVectorClear(mTimeFrame->getCellsNeighboursLUT()[iLayer]);
492 if (mTimeFrame->getCells()[iLayer + 1].empty() ||
493 mTimeFrame->getCellsLookupTable()[iLayer].empty()) {
494 continue;
495 }
496
497 int nCells{static_cast<int>(mTimeFrame->getCells()[iLayer].size())};
498 bounded_vector<Neighbor> cellsNeighbours(mMemoryPool.get());
499
500 auto forCellNeighbour = [&](auto Tag, int iCell, int offset = 0) -> int {
501 const auto& currentCellSeed{mTimeFrame->getCells()[iLayer][iCell]};
502 const int nextLayerTrackletIndex{currentCellSeed.getSecondTrackletIndex()};
503 const int nextLayerFirstCellIndex{mTimeFrame->getCellsLookupTable()[iLayer][nextLayerTrackletIndex]};
504 const int nextLayerLastCellIndex{mTimeFrame->getCellsLookupTable()[iLayer][nextLayerTrackletIndex + 1]};
505 int foundNextCells{0};
506 for (int iNextCell{nextLayerFirstCellIndex}; iNextCell < nextLayerLastCellIndex; ++iNextCell) {
507 auto nextCellSeed{mTimeFrame->getCells()[iLayer + 1][iNextCell]};
508 if (nextCellSeed.getFirstTrackletIndex() != nextLayerTrackletIndex) {
509 break;
510 }
511
512 if (mTrkParams[iteration].DeltaROF) { // TODO this has to be improved for the staggering
513 const auto& trkl00 = mTimeFrame->getTracklets()[iLayer][currentCellSeed.getFirstTrackletIndex()];
514 const auto& trkl01 = mTimeFrame->getTracklets()[iLayer + 1][currentCellSeed.getSecondTrackletIndex()];
515 const auto& trkl10 = mTimeFrame->getTracklets()[iLayer + 1][nextCellSeed.getFirstTrackletIndex()];
516 const auto& trkl11 = mTimeFrame->getTracklets()[iLayer + 2][nextCellSeed.getSecondTrackletIndex()];
517 if ((std::max({trkl00.getMaxRof(), trkl01.getMaxRof(), trkl10.getMaxRof(), trkl11.getMaxRof()}) -
518 std::min({trkl00.getMinRof(), trkl01.getMinRof(), trkl10.getMinRof(), trkl11.getMinRof()})) > mTrkParams[0].DeltaROF) {
519 continue;
520 }
521 }
522
523 if (!nextCellSeed.rotate(currentCellSeed.getAlpha()) ||
524 !nextCellSeed.propagateTo(currentCellSeed.getX(), getBz())) {
525 continue;
526 }
527 float chi2 = currentCellSeed.getPredictedChi2(nextCellSeed);
528
529#ifdef OPTIMISATION_OUTPUT
530 bool good{mTimeFrame->getCellsLabel(iLayer)[iCell] == mTimeFrame->getCellsLabel(iLayer + 1)[iNextCell]};
531 off << std::format("{}\t{:d}\t{}", iLayer, good, chi2) << std::endl;
532#endif
533
534 if (chi2 > mTrkParams[0].MaxChi2ClusterAttachment) {
535 continue;
536 }
537
538 if constexpr (decltype(Tag)::value == PassMode::OnePass::value) {
539 cellsNeighbours.emplace_back(iCell, iNextCell, currentCellSeed.getLevel() + 1);
540 } else if constexpr (decltype(Tag)::value == PassMode::TwoPassCount::value) {
541 ++foundNextCells;
542 } else if constexpr (decltype(Tag)::value == PassMode::TwoPassInsert::value) {
543 cellsNeighbours[offset++] = {iCell, iNextCell, currentCellSeed.getLevel() + 1};
544 } else {
545 static_assert(false, "Unknown mode!");
546 }
547 }
548 return foundNextCells;
549 };
550
551 if (mTaskArena->max_concurrency() <= 1) {
552 for (int iCell{0}; iCell < nCells; ++iCell) {
553 forCellNeighbour(PassMode::OnePass{}, iCell);
554 }
555 } else {
556 bounded_vector<int> perCellCount(nCells + 1, 0, mMemoryPool.get());
557 tbb::parallel_for(
558 tbb::blocked_range<int>(0, nCells),
559 [&](const tbb::blocked_range<int>& Cells) {
560 for (int iCell = Cells.begin(); iCell < Cells.end(); ++iCell) {
561 perCellCount[iCell] = forCellNeighbour(PassMode::TwoPassCount{}, iCell);
562 }
563 });
564
565 std::exclusive_scan(perCellCount.begin(), perCellCount.end(), perCellCount.begin(), 0);
566 int totalCellNeighbours = perCellCount.back();
567 if (totalCellNeighbours == 0) {
568 deepVectorClear(mTimeFrame->getCellsNeighbours()[iLayer]);
569 continue;
570 }
571 cellsNeighbours.resize(totalCellNeighbours);
572
573 tbb::parallel_for(
574 tbb::blocked_range<int>(0, nCells),
575 [&](const tbb::blocked_range<int>& Cells) {
576 for (int iCell = Cells.begin(); iCell < Cells.end(); ++iCell) {
577 int offset = perCellCount[iCell];
578 if (offset == perCellCount[iCell + 1]) {
579 continue;
580 }
581 forCellNeighbour(PassMode::TwoPassInsert{}, iCell, offset);
582 }
583 });
584 }
585
586 if (cellsNeighbours.empty()) {
587 continue;
588 }
589
590 tbb::parallel_sort(cellsNeighbours.begin(), cellsNeighbours.end(), [](const auto& a, const auto& b) {
591 return a.nextCell < b.nextCell;
592 });
593
594 auto& cellsNeighbourLUT = mTimeFrame->getCellsNeighboursLUT()[iLayer];
595 cellsNeighbourLUT.assign(mTimeFrame->getCells()[iLayer + 1].size(), 0);
596 for (const auto& neigh : cellsNeighbours) {
597 ++cellsNeighbourLUT[neigh.nextCell];
598 }
599 std::inclusive_scan(cellsNeighbourLUT.begin(), cellsNeighbourLUT.end(), cellsNeighbourLUT.begin());
600
601 mTimeFrame->getCellsNeighbours()[iLayer].reserve(cellsNeighbours.size());
602 std::ranges::transform(cellsNeighbours, std::back_inserter(mTimeFrame->getCellsNeighbours()[iLayer]), [](const auto& neigh) { return neigh.cell; });
603
604 for (auto it = cellsNeighbours.begin(); it != cellsNeighbours.end();) {
605 int cellIdx = it->nextCell;
606 int maxLvl = it->level;
607 while (++it != cellsNeighbours.end() && it->nextCell == cellIdx) {
608 maxLvl = std::max(maxLvl, it->level);
609 }
610 o2::gpu::CAMath::AtomicMax(mTimeFrame->getCells()[iLayer + 1][cellIdx].getLevelPtr(), maxLvl);
611 }
612 }
613 });
614}
615
616template <int nLayers>
617void TrackerTraits<nLayers>::processNeighbours(int iLayer, int iLevel, const bounded_vector<CellSeedN>& currentCellSeed, const bounded_vector<int>& currentCellId, bounded_vector<CellSeedN>& updatedCellSeeds, bounded_vector<int>& updatedCellsIds)
618{
619 CA_DEBUGGER(std::cout << "Processing neighbours layer " << iLayer << " level " << iLevel << ", size of the cell seeds: " << currentCellSeed.size() << std::endl);
620 auto propagator = o2::base::Propagator::Instance();
621
622#ifdef CA_DEBUG
623 int failed[5]{0, 0, 0, 0, 0}, attempts{0}, failedByMismatch{0};
624#endif
625
626 mTaskArena->execute([&] {
627 auto forCellNeighbours = [&](auto Tag, int iCell, int offset = 0) -> int {
628 const auto& currentCell{currentCellSeed[iCell]};
629
630 if constexpr (decltype(Tag)::value != PassMode::TwoPassInsert::value) {
631 if (currentCell.getLevel() != iLevel) {
632 return 0;
633 }
634 if (currentCellId.empty() && (mTimeFrame->isClusterUsed(iLayer, currentCell.getFirstClusterIndex()) ||
635 mTimeFrame->isClusterUsed(iLayer + 1, currentCell.getSecondClusterIndex()) ||
636 mTimeFrame->isClusterUsed(iLayer + 2, currentCell.getThirdClusterIndex()))) {
637 return 0;
638 }
639 }
640
641 const int cellId = currentCellId.empty() ? iCell : currentCellId[iCell];
642 const int startNeighbourId{cellId ? mTimeFrame->getCellsNeighboursLUT()[iLayer - 1][cellId - 1] : 0};
643 const int endNeighbourId{mTimeFrame->getCellsNeighboursLUT()[iLayer - 1][cellId]};
644 int foundSeeds{0};
645 for (int iNeighbourCell{startNeighbourId}; iNeighbourCell < endNeighbourId; ++iNeighbourCell) {
646 CA_DEBUGGER(attempts++);
647 const int neighbourCellId = mTimeFrame->getCellsNeighbours()[iLayer - 1][iNeighbourCell];
648 const auto& neighbourCell = mTimeFrame->getCells()[iLayer - 1][neighbourCellId];
649 if (neighbourCell.getSecondTrackletIndex() != currentCell.getFirstTrackletIndex()) {
650 CA_DEBUGGER(failedByMismatch++);
651 continue;
652 }
653 if (mTimeFrame->isClusterUsed(iLayer - 1, neighbourCell.getFirstClusterIndex())) {
654 continue;
655 }
656 if (currentCell.getLevel() - 1 != neighbourCell.getLevel()) {
657 CA_DEBUGGER(failed[0]++);
658 continue;
659 }
660
662 CellSeedN seed{currentCell};
663 const auto& trHit = mTimeFrame->getTrackingFrameInfoOnLayer(iLayer - 1)[neighbourCell.getFirstClusterIndex()];
664
665 if (!seed.rotate(trHit.alphaTrackingFrame)) {
666 CA_DEBUGGER(failed[1]++);
667 continue;
668 }
669
670 if (!propagator->propagateToX(seed, trHit.xTrackingFrame, getBz(), o2::base::PropagatorImpl<float>::MAX_SIN_PHI, o2::base::PropagatorImpl<float>::MAX_STEP, mTrkParams[0].CorrType)) {
671 CA_DEBUGGER(failed[2]++);
672 continue;
673 }
674
675 if (mTrkParams[0].CorrType == o2::base::PropagatorF::MatCorrType::USEMatCorrNONE) {
676 if (!seed.correctForMaterial(mTrkParams[0].LayerxX0[iLayer - 1], mTrkParams[0].LayerxX0[iLayer - 1] * constants::Radl * constants::Rho, true)) {
677 continue;
678 }
679 }
680
681 auto predChi2{seed.getPredictedChi2Quiet(trHit.positionTrackingFrame, trHit.covarianceTrackingFrame)};
682 if ((predChi2 > mTrkParams[0].MaxChi2ClusterAttachment) || predChi2 < 0.f) {
683 CA_DEBUGGER(failed[3]++);
684 continue;
685 }
686 seed.setChi2(seed.getChi2() + predChi2);
687 if (!seed.o2::track::TrackParCov::update(trHit.positionTrackingFrame, trHit.covarianceTrackingFrame)) {
688 CA_DEBUGGER(failed[4]++);
689 continue;
690 }
691
692 if constexpr (decltype(Tag)::value != PassMode::TwoPassCount::value) {
693 seed.getClusters()[iLayer - 1] = neighbourCell.getFirstClusterIndex();
694 seed.setLevel(neighbourCell.getLevel());
695 seed.setFirstTrackletIndex(neighbourCell.getFirstTrackletIndex());
696 seed.setSecondTrackletIndex(neighbourCell.getSecondTrackletIndex());
697 }
698
699 if constexpr (decltype(Tag)::value == PassMode::OnePass::value) {
700 updatedCellSeeds.push_back(seed);
701 updatedCellsIds.push_back(neighbourCellId);
702 } else if constexpr (decltype(Tag)::value == PassMode::TwoPassCount::value) {
703 ++foundSeeds;
704 } else if constexpr (decltype(Tag)::value == PassMode::TwoPassInsert::value) {
705 updatedCellSeeds[offset] = seed;
706 updatedCellsIds[offset++] = neighbourCellId;
707 } else {
708 static_assert(false, "Unknown mode!");
709 }
710 }
711 return foundSeeds;
712 };
713
714 const int nCells = static_cast<int>(currentCellSeed.size());
715 if (mTaskArena->max_concurrency() <= 1) {
716 for (int iCell{0}; iCell < nCells; ++iCell) {
717 forCellNeighbours(PassMode::OnePass{}, iCell);
718 }
719 } else {
720 bounded_vector<int> perCellCount(nCells + 1, 0, mMemoryPool.get());
721 tbb::parallel_for(
722 tbb::blocked_range<int>(0, nCells),
723 [&](const tbb::blocked_range<int>& Cells) {
724 for (int iCell = Cells.begin(); iCell < Cells.end(); ++iCell) {
725 perCellCount[iCell] = forCellNeighbours(PassMode::TwoPassCount{}, iCell);
726 }
727 });
728
729 std::exclusive_scan(perCellCount.begin(), perCellCount.end(), perCellCount.begin(), 0);
730 auto totalNeighbours{perCellCount.back()};
731 if (totalNeighbours == 0) {
732 return;
733 }
734 updatedCellSeeds.resize(totalNeighbours);
735 updatedCellsIds.resize(totalNeighbours);
736
737 tbb::parallel_for(
738 tbb::blocked_range<int>(0, nCells),
739 [&](const tbb::blocked_range<int>& Cells) {
740 for (int iCell = Cells.begin(); iCell < Cells.end(); ++iCell) {
741 int offset = perCellCount[iCell];
742 if (offset == perCellCount[iCell + 1]) {
743 continue;
744 }
745 forCellNeighbours(PassMode::TwoPassInsert{}, iCell, offset);
746 }
747 });
748 }
749 });
750
751#ifdef CA_DEBUG
752 std::cout << "\t\t- Found " << updatedCellSeeds.size() << " cell seeds out of " << attempts << " attempts" << std::endl;
753 std::cout << "\t\t\t> " << failed[0] << " failed because of level" << std::endl;
754 std::cout << "\t\t\t> " << failed[1] << " failed because of rotation" << std::endl;
755 std::cout << "\t\t\t> " << failed[2] << " failed because of propagation" << std::endl;
756 std::cout << "\t\t\t> " << failed[3] << " failed because of chi2 cut" << std::endl;
757 std::cout << "\t\t\t> " << failed[4] << " failed because of update" << std::endl;
758 std::cout << "\t\t\t> " << failedByMismatch << " failed because of mismatch" << std::endl;
759#endif
760}
761
762template <int nLayers>
763void TrackerTraits<nLayers>::findRoads(const int iteration)
764{
765 CA_DEBUGGER(std::cout << "Finding roads, iteration " << iteration << std::endl);
766
767 for (int startLevel{mTrkParams[iteration].CellsPerRoad()}; startLevel >= mTrkParams[iteration].CellMinimumLevel(); --startLevel) {
768 CA_DEBUGGER(std::cout << "\t > Processing level " << startLevel << std::endl);
769 auto seedFilter = [&](const auto& seed) {
770 return seed.getQ2Pt() <= 1.e3 && seed.getChi2() <= mTrkParams[0].MaxChi2NDF * ((startLevel + 2) * 2 - 5);
771 };
772 bounded_vector<CellSeedN> trackSeeds(mMemoryPool.get());
773 for (int startLayer{mTrkParams[iteration].CellsPerRoad() - 1}; startLayer >= startLevel - 1; --startLayer) {
774 if ((mTrkParams[iteration].StartLayerMask & (1 << (startLayer + 2))) == 0) {
775 continue;
776 }
777 CA_DEBUGGER(std::cout << "\t\t > Starting processing layer " << startLayer << std::endl);
778 bounded_vector<int> lastCellId(mMemoryPool.get()), updatedCellId(mMemoryPool.get());
779 bounded_vector<CellSeedN> lastCellSeed(mMemoryPool.get()), updatedCellSeed(mMemoryPool.get());
780
781 processNeighbours(startLayer, startLevel, mTimeFrame->getCells()[startLayer], lastCellId, updatedCellSeed, updatedCellId);
782
783 int level = startLevel;
784 for (int iLayer{startLayer - 1}; iLayer > 0 && level > 2; --iLayer) {
785 lastCellSeed.swap(updatedCellSeed);
786 lastCellId.swap(updatedCellId);
787 deepVectorClear(updatedCellSeed);
788 deepVectorClear(updatedCellId);
789 processNeighbours(iLayer, --level, lastCellSeed, lastCellId, updatedCellSeed, updatedCellId);
790 }
791 deepVectorClear(lastCellId);
792 deepVectorClear(lastCellSeed);
793
794 if (!updatedCellSeed.empty()) {
795 trackSeeds.reserve(trackSeeds.size() + std::count_if(updatedCellSeed.begin(), updatedCellSeed.end(), seedFilter));
796 std::copy_if(updatedCellSeed.begin(), updatedCellSeed.end(), std::back_inserter(trackSeeds), seedFilter);
797 }
798 }
799
800 if (trackSeeds.empty()) {
801 continue;
802 }
803
804 bounded_vector<TrackITSExt> tracks(mMemoryPool.get());
805 mTaskArena->execute([&] {
806 auto forSeed = [&](auto Tag, int iSeed, int offset = 0) {
807 const auto& seed{trackSeeds[iSeed]};
808 TrackITSExt temporaryTrack{seed};
809 temporaryTrack.resetCovariance();
810 temporaryTrack.setChi2(0);
811 for (int iL{0}; iL < nLayers; ++iL) {
812 temporaryTrack.setExternalClusterIndex(iL, seed.getCluster(iL), seed.getCluster(iL) != constants::UnusedIndex);
813 }
814
815 bool fitSuccess = fitTrack(temporaryTrack, 0, mTrkParams[0].NLayers, 1, mTrkParams[0].MaxChi2ClusterAttachment, mTrkParams[0].MaxChi2NDF);
816 if (!fitSuccess) {
817 return 0;
818 }
819
820 temporaryTrack.getParamOut() = temporaryTrack.getParamIn();
821 temporaryTrack.resetCovariance();
822 temporaryTrack.setChi2(0);
823 fitSuccess = fitTrack(temporaryTrack, mTrkParams[0].NLayers - 1, -1, -1, mTrkParams[0].MaxChi2ClusterAttachment, mTrkParams[0].MaxChi2NDF, 50.f);
824 if (!fitSuccess || temporaryTrack.getPt() < mTrkParams[iteration].MinPt[mTrkParams[iteration].NLayers - temporaryTrack.getNClusters()]) {
825 return 0;
826 }
827
828 if constexpr (decltype(Tag)::value == PassMode::OnePass::value) {
829 tracks.push_back(temporaryTrack);
830 } else if constexpr (decltype(Tag)::value == PassMode::TwoPassCount::value) {
831 // nothing to do
832 } else if constexpr (decltype(Tag)::value == PassMode::TwoPassInsert::value) {
833 tracks[offset] = temporaryTrack;
834 } else {
835 static_assert(false, "Unknown mode!");
836 }
837 return 1;
838 };
839
840 const int nSeeds = static_cast<int>(trackSeeds.size());
841 if (mTaskArena->max_concurrency() <= 1) {
842 for (int iSeed{0}; iSeed < nSeeds; ++iSeed) {
843 forSeed(PassMode::OnePass{}, iSeed);
844 }
845 } else {
846 bounded_vector<int> perSeedCount(nSeeds + 1, 0, mMemoryPool.get());
847 tbb::parallel_for(
848 tbb::blocked_range<int>(0, nSeeds),
849 [&](const tbb::blocked_range<int>& Seeds) {
850 for (int iSeed = Seeds.begin(); iSeed < Seeds.end(); ++iSeed) {
851 perSeedCount[iSeed] = forSeed(PassMode::TwoPassCount{}, iSeed);
852 }
853 });
854
855 std::exclusive_scan(perSeedCount.begin(), perSeedCount.end(), perSeedCount.begin(), 0);
856 auto totalTracks{perSeedCount.back()};
857 if (totalTracks == 0) {
858 return;
859 }
860 tracks.resize(totalTracks);
861
862 tbb::parallel_for(
863 tbb::blocked_range<int>(0, nSeeds),
864 [&](const tbb::blocked_range<int>& Seeds) {
865 for (int iSeed = Seeds.begin(); iSeed < Seeds.end(); ++iSeed) {
866 if (perSeedCount[iSeed] == perSeedCount[iSeed + 1]) {
867 continue;
868 }
869 forSeed(PassMode::TwoPassInsert{}, iSeed, perSeedCount[iSeed]);
870 }
871 });
872 }
873
874 deepVectorClear(trackSeeds);
875 tbb::parallel_sort(tracks.begin(), tracks.end(), [](const auto& a, const auto& b) {
876 return a.getChi2() < b.getChi2();
877 });
878 });
879
880 for (auto& track : tracks) {
881 int nShared = 0;
882 bool isFirstShared{false};
883 for (int iLayer{0}; iLayer < mTrkParams[0].NLayers; ++iLayer) {
884 if (track.getClusterIndex(iLayer) == constants::UnusedIndex) {
885 continue;
886 }
887 nShared += int(mTimeFrame->isClusterUsed(iLayer, track.getClusterIndex(iLayer)));
888 isFirstShared |= !iLayer && mTimeFrame->isClusterUsed(iLayer, track.getClusterIndex(iLayer));
889 }
890
891 if (nShared > mTrkParams[0].ClusterSharing) {
892 continue;
893 }
894
895 std::array<int, 3> rofs{INT_MAX, INT_MAX, INT_MAX};
896 for (int iLayer{0}; iLayer < mTrkParams[0].NLayers; ++iLayer) {
897 if (track.getClusterIndex(iLayer) == constants::UnusedIndex) {
898 continue;
899 }
900 mTimeFrame->markUsedCluster(iLayer, track.getClusterIndex(iLayer));
901 int currentROF = mTimeFrame->getClusterROF(iLayer, track.getClusterIndex(iLayer));
902 for (int iR{0}; iR < 3; ++iR) {
903 if (rofs[iR] == INT_MAX) {
904 rofs[iR] = currentROF;
905 }
906 if (rofs[iR] == currentROF) {
907 break;
908 }
909 }
910 }
911 if (rofs[2] != INT_MAX) {
912 continue;
913 }
914 track.setUserField(0);
915 track.getParamOut().setUserField(0);
916 if (rofs[1] != INT_MAX) {
917 track.setNextROFbit();
918 }
919 mTimeFrame->getTracks(o2::gpu::CAMath::Min(rofs[0], rofs[1])).emplace_back(track);
920 }
921 }
922}
923
924template <int nLayers>
926{
927 for (int rof{0}; rof < mTimeFrame->getNrof(); ++rof) {
928 for (auto& track : mTimeFrame->getTracks(rof)) {
929 auto backup{track};
930 bool success{false};
931 // the order here biases towards top extension, tracks should probably be fitted separately in the directions and then compared.
932 if ((mTrkParams[iteration].UseTrackFollowerMix || mTrkParams[iteration].UseTrackFollowerTop) && track.getLastClusterLayer() != mTrkParams[iteration].NLayers - 1) {
933 success = success || trackFollowing(&track, rof, true, iteration);
934 }
935 if ((mTrkParams[iteration].UseTrackFollowerMix || (mTrkParams[iteration].UseTrackFollowerBot && !success)) && track.getFirstClusterLayer() != 0) {
936 success = success || trackFollowing(&track, rof, false, iteration);
937 }
938 if (success) {
940 track.resetCovariance();
941 track.setChi2(0);
942 bool fitSuccess = fitTrack(track, 0, mTrkParams[iteration].NLayers, 1, mTrkParams[iteration].MaxChi2ClusterAttachment, mTrkParams[0].MaxChi2NDF);
943 if (!fitSuccess) {
944 track = backup;
945 continue;
946 }
947 track.getParamOut() = track;
948 track.resetCovariance();
949 track.setChi2(0);
950 fitSuccess = fitTrack(track, mTrkParams[iteration].NLayers - 1, -1, -1, mTrkParams[iteration].MaxChi2ClusterAttachment, mTrkParams[0].MaxChi2NDF, 50.);
951 if (!fitSuccess) {
952 track = backup;
953 continue;
954 }
955 mTimeFrame->mNExtendedTracks++;
956 mTimeFrame->mNExtendedUsedClusters += track.getNClusters() - backup.getNClusters();
957 auto pattern = track.getPattern();
958 auto diff = (pattern & ~backup.getPattern()) & 0xff;
959 pattern |= (diff << 24);
960 track.setPattern(pattern);
962 for (int iLayer{0}; iLayer < mTrkParams[iteration].NLayers; ++iLayer) {
963 if (track.getClusterIndex(iLayer) == constants::UnusedIndex) {
964 continue;
965 }
966 mTimeFrame->markUsedCluster(iLayer, track.getClusterIndex(iLayer));
967 }
968 }
969 }
970 }
971}
972
973template <int nLayers>
975{
976 const auto propagator = o2::base::Propagator::Instance();
977 mTimeFrame->fillPrimaryVerticesXandAlpha();
978
979 for (auto& cell : mTimeFrame->getCells()[0]) {
980 auto& cluster3_glo = mTimeFrame->getClusters()[2][cell.getThirdClusterIndex()];
981 auto& cluster2_glo = mTimeFrame->getClusters()[1][cell.getSecondClusterIndex()];
982 auto& cluster1_glo = mTimeFrame->getClusters()[0][cell.getFirstClusterIndex()];
983 if (mTimeFrame->isClusterUsed(2, cluster1_glo.clusterId) ||
984 mTimeFrame->isClusterUsed(1, cluster2_glo.clusterId) ||
985 mTimeFrame->isClusterUsed(0, cluster3_glo.clusterId)) {
986 continue;
987 }
988
989 std::array<int, 3> rofs{
990 mTimeFrame->getClusterROF(2, cluster3_glo.clusterId),
991 mTimeFrame->getClusterROF(1, cluster2_glo.clusterId),
992 mTimeFrame->getClusterROF(0, cluster1_glo.clusterId)};
993 if (rofs[0] != rofs[1] && rofs[1] != rofs[2] && rofs[0] != rofs[2]) {
994 continue;
995 }
996
997 int rof{rofs[0]};
998 if (rofs[1] == rofs[2]) {
999 rof = rofs[2];
1000 }
1001
1002 auto pvs{mTimeFrame->getPrimaryVertices(rof)};
1003 auto pvsXAlpha{mTimeFrame->getPrimaryVerticesXAlpha(rof)};
1004
1005 const auto& cluster3_tf = mTimeFrame->getTrackingFrameInfoOnLayer(2)[cluster3_glo.clusterId];
1006 TrackITSExt temporaryTrack{buildTrackSeed(cluster1_glo, cluster2_glo, cluster3_tf)};
1007 temporaryTrack.setExternalClusterIndex(0, cluster1_glo.clusterId, true);
1008 temporaryTrack.setExternalClusterIndex(1, cluster2_glo.clusterId, true);
1009 temporaryTrack.setExternalClusterIndex(2, cluster3_glo.clusterId, true);
1010
1012 bool fitSuccess = fitTrack(temporaryTrack, 1, -1, -1);
1013 if (!fitSuccess) {
1014 continue;
1015 }
1016 fitSuccess = false;
1017
1018 TrackITSExt bestTrack{temporaryTrack}, backup{temporaryTrack};
1019 float bestChi2{std::numeric_limits<float>::max()};
1020 for (int iV{0}; iV < (int)pvs.size(); ++iV) {
1021 temporaryTrack = backup;
1022 if (!temporaryTrack.rotate(pvsXAlpha[iV][1])) {
1023 continue;
1024 }
1025 if (!propagator->propagateTo(temporaryTrack, pvsXAlpha[iV][0], true)) {
1026 continue;
1027 }
1028
1029 float pvRes{mTrkParams[0].PVres / o2::gpu::CAMath::Sqrt(float(pvs[iV].getNContributors()))};
1030 const float posVtx[2]{0.f, pvs[iV].getZ()};
1031 const float covVtx[3]{pvRes, 0.f, pvRes};
1032 float chi2 = temporaryTrack.getPredictedChi2Quiet(posVtx, covVtx);
1033 if (chi2 < bestChi2) {
1034 if (!temporaryTrack.track::TrackParCov::update(posVtx, covVtx)) {
1035 continue;
1036 }
1037 bestTrack = temporaryTrack;
1038 bestChi2 = chi2;
1039 }
1040 }
1041
1042 bestTrack.resetCovariance();
1043 bestTrack.setChi2(0.f);
1044 fitSuccess = fitTrack(bestTrack, 0, mTrkParams[0].NLayers, 1, mTrkParams[0].MaxChi2ClusterAttachment, mTrkParams[0].MaxChi2NDF);
1045 if (!fitSuccess) {
1046 continue;
1047 }
1048 bestTrack.getParamOut() = bestTrack;
1049 bestTrack.resetCovariance();
1050 bestTrack.setChi2(0.f);
1051 fitSuccess = fitTrack(bestTrack, mTrkParams[0].NLayers - 1, -1, -1, mTrkParams[0].MaxChi2ClusterAttachment, mTrkParams[0].MaxChi2NDF, 50.);
1052 if (!fitSuccess) {
1053 continue;
1054 }
1055 mTimeFrame->markUsedCluster(0, bestTrack.getClusterIndex(0));
1056 mTimeFrame->markUsedCluster(1, bestTrack.getClusterIndex(1));
1057 mTimeFrame->markUsedCluster(2, bestTrack.getClusterIndex(2));
1058 mTimeFrame->getTracks(rof).emplace_back(bestTrack);
1059 }
1060}
1061
1062template <int nLayers>
1063bool TrackerTraits<nLayers>::fitTrack(TrackITSExt& track, int start, int end, int step, float chi2clcut, float chi2ndfcut, float maxQoverPt, int nCl)
1064{
1065 auto propInstance = o2::base::Propagator::Instance();
1066
1067 for (int iLayer{start}; iLayer != end; iLayer += step) {
1068 if (track.getClusterIndex(iLayer) == constants::UnusedIndex) {
1069 continue;
1070 }
1071 const TrackingFrameInfo& trackingHit = mTimeFrame->getTrackingFrameInfoOnLayer(iLayer)[track.getClusterIndex(iLayer)];
1072
1073 if (!track.rotate(trackingHit.alphaTrackingFrame)) {
1074 return false;
1075 }
1076
1077 if (!propInstance->propagateToX(track, trackingHit.xTrackingFrame, getBz(), o2::base::PropagatorImpl<float>::MAX_SIN_PHI, o2::base::PropagatorImpl<float>::MAX_STEP, mTrkParams[0].CorrType)) {
1078 return false;
1079 }
1080
1081 if (mTrkParams[0].CorrType == o2::base::PropagatorF::MatCorrType::USEMatCorrNONE) {
1082 if (!track.correctForMaterial(mTrkParams[0].LayerxX0[iLayer], mTrkParams[0].LayerxX0[iLayer] * constants::Radl * constants::Rho, true)) {
1083 continue;
1084 }
1085 }
1086
1087 auto predChi2{track.getPredictedChi2Quiet(trackingHit.positionTrackingFrame, trackingHit.covarianceTrackingFrame)};
1088 if ((nCl >= 3 && predChi2 > chi2clcut) || predChi2 < 0.f) {
1089 return false;
1090 }
1091 track.setChi2(track.getChi2() + predChi2);
1092 if (!track.o2::track::TrackParCov::update(trackingHit.positionTrackingFrame, trackingHit.covarianceTrackingFrame)) {
1093 return false;
1094 }
1095 nCl++;
1096 }
1097 return std::abs(track.getQ2Pt()) < maxQoverPt && track.getChi2() < chi2ndfcut * (nCl * 2 - 5);
1098}
1099
1100template <int nLayers>
1101bool TrackerTraits<nLayers>::trackFollowing(TrackITSExt* track, int rof, bool outward, const int iteration)
1102{
1103 auto propInstance = o2::base::Propagator::Instance();
1104 const int step = -1 + outward * 2;
1105 const int end = outward ? mTrkParams[iteration].NLayers - 1 : 0;
1106 bounded_vector<TrackITSExt> hypotheses(1, *track, mMemoryPool.get()); // possibly avoid reallocation
1107 for (size_t iHypo{0}; iHypo < hypotheses.size(); ++iHypo) {
1108 auto hypo{hypotheses[iHypo]};
1109 int iLayer = static_cast<int>(outward ? hypo.getLastClusterLayer() : hypo.getFirstClusterLayer());
1110 // per layer we add new hypotheses
1111 while (iLayer != end) {
1112 iLayer += step; // step through all layers until we reach the end, this allows for skipping on empty layers
1113 const float r = mTrkParams[iteration].LayerRadii[iLayer];
1114 // get an estimate of the trackinf-frame x for the next step
1115 float x{-999};
1116 if (!hypo.getXatLabR(r, x, mTimeFrame->getBz(), o2::track::DirAuto) || x <= 0.f) {
1117 continue;
1118 }
1119 // estimate hypo's trk parameters at that x
1120 auto& hypoParam{outward ? hypo.getParamOut() : hypo.getParamIn()};
1121 if (!propInstance->propagateToX(hypoParam, x, mTimeFrame->getBz(), PropagatorF::MAX_SIN_PHI,
1122 PropagatorF::MAX_STEP, mTrkParams[iteration].CorrType)) {
1123 continue;
1124 }
1125
1126 if (mTrkParams[iteration].CorrType == PropagatorF::MatCorrType::USEMatCorrNONE) { // account for material affects if propagator does not
1127 if (!hypoParam.correctForMaterial(mTrkParams[iteration].LayerxX0[iLayer], mTrkParams[iteration].LayerxX0[iLayer] * constants::Radl * constants::Rho, true)) {
1128 continue;
1129 }
1130 }
1131
1132 // calculate the search window on this layer
1133 const float phi{hypoParam.getPhi()};
1134 const float ePhi{o2::gpu::CAMath::Sqrt(hypoParam.getSigmaSnp2() / hypoParam.getCsp2())};
1135 const float z{hypoParam.getZ()};
1136 const float eZ{o2::gpu::CAMath::Sqrt(hypoParam.getSigmaZ2())};
1137 const int4 selectedBinsRect{getBinsRect(iLayer, phi, mTrkParams[iteration].TrackFollowerNSigmaCutPhi * ePhi, z, mTrkParams[iteration].TrackFollowerNSigmaCutZ * eZ)};
1138 if (selectedBinsRect.x == 0 && selectedBinsRect.y == 0 && selectedBinsRect.z == 0 && selectedBinsRect.w == 0) {
1139 continue;
1140 }
1141
1142 int phiBinsNum{selectedBinsRect.w - selectedBinsRect.y + 1};
1143
1144 if (phiBinsNum < 0) {
1145 phiBinsNum += mTrkParams[iteration].PhiBins;
1146 }
1147
1148 gsl::span<const Cluster> layer1 = mTimeFrame->getClustersOnLayer(rof, iLayer);
1149 if (layer1.empty()) {
1150 continue;
1151 }
1152
1153 // check all clusters in search windows for possible new hypotheses
1154 for (int iPhiCount = 0; iPhiCount < phiBinsNum; iPhiCount++) {
1155 int iPhiBin = (selectedBinsRect.y + iPhiCount) % mTrkParams[iteration].PhiBins;
1156 const int firstBinIndex{mTimeFrame->mIndexTableUtils.getBinIndex(selectedBinsRect.x, iPhiBin)};
1157 const int maxBinIndex{firstBinIndex + selectedBinsRect.z - selectedBinsRect.x + 1};
1158 const int firstRowClusterIndex = mTimeFrame->getIndexTable(rof, iLayer)[firstBinIndex];
1159 const int maxRowClusterIndex = mTimeFrame->getIndexTable(rof, iLayer)[maxBinIndex];
1160
1161 for (int iNextCluster{firstRowClusterIndex}; iNextCluster < maxRowClusterIndex; ++iNextCluster) {
1162 if (iNextCluster >= (int)layer1.size()) {
1163 break;
1164 }
1165 const Cluster& nextCluster{layer1[iNextCluster]};
1166
1167 if (mTimeFrame->isClusterUsed(iLayer, nextCluster.clusterId)) {
1168 continue;
1169 }
1170
1171 const TrackingFrameInfo& trackingHit = mTimeFrame->getTrackingFrameInfoOnLayer(iLayer)[nextCluster.clusterId];
1172
1173 auto tbupdated{hypo};
1174 auto& tbuParams = outward ? tbupdated.getParamOut() : tbupdated.getParamIn();
1175 if (!tbuParams.rotate(trackingHit.alphaTrackingFrame)) {
1176 continue;
1177 }
1178
1179 if (!propInstance->propagateToX(tbuParams, trackingHit.xTrackingFrame, mTimeFrame->getBz(),
1180 PropagatorF::MAX_SIN_PHI, PropagatorF::MAX_STEP, PropagatorF::MatCorrType::USEMatCorrNONE)) {
1181 continue;
1182 }
1183
1184 auto predChi2{tbuParams.getPredictedChi2Quiet(trackingHit.positionTrackingFrame, trackingHit.covarianceTrackingFrame)};
1185 if (predChi2 >= track->getChi2() * mTrkParams[iteration].NSigmaCut) {
1186 continue;
1187 }
1188
1189 if (!tbuParams.o2::track::TrackParCov::update(trackingHit.positionTrackingFrame, trackingHit.covarianceTrackingFrame)) {
1190 continue;
1191 }
1192 tbupdated.setChi2(tbupdated.getChi2() + predChi2);
1193 tbupdated.setExternalClusterIndex(iLayer, nextCluster.clusterId, true);
1194 hypotheses.emplace_back(tbupdated);
1195 }
1196 }
1197 }
1198 }
1199
1200 TrackITSExt* bestHypo{track};
1201 bool swapped{false};
1202 for (auto& hypo : hypotheses) {
1203 if (hypo.isBetter(*bestHypo, track->getChi2() * mTrkParams[iteration].NSigmaCut)) {
1204 bestHypo = &hypo;
1205 swapped = true;
1206 }
1207 }
1208 *track = *bestHypo;
1209 return swapped;
1210}
1211
1214template <int nLayers>
1216{
1217 float ca{-999.f}, sa{-999.f};
1218 o2::gpu::CAMath::SinCos(tf3.alphaTrackingFrame, sa, ca);
1219 const float x1 = cluster1.xCoordinate * ca + cluster1.yCoordinate * sa;
1220 const float y1 = -cluster1.xCoordinate * sa + cluster1.yCoordinate * ca;
1221 const float z1 = cluster1.zCoordinate;
1222 const float x2 = cluster2.xCoordinate * ca + cluster2.yCoordinate * sa;
1223 const float y2 = -cluster2.xCoordinate * sa + cluster2.yCoordinate * ca;
1224 const float z2 = cluster2.zCoordinate;
1225 const float x3 = tf3.xTrackingFrame;
1226 const float y3 = tf3.positionTrackingFrame[0];
1227 const float z3 = tf3.positionTrackingFrame[1];
1228 float tgp{1.f}, crv{1.f}, snp{-999.f}, tgl12{-999.f}, tgl23{-999.f}, q2pt{1.f / track::kMostProbablePt}, q2pt2{1.f}, sg2q2pt{-999.f};
1229 if (mIsZeroField) {
1230 tgp = o2::gpu::CAMath::ATan2(y3 - y1, x3 - x1);
1231 snp = tgp / o2::gpu::CAMath::Sqrt(1.f + tgp * tgp);
1232 } else {
1233 crv = math_utils::computeCurvature(x3, y3, x2, y2, x1, y1);
1234 snp = crv * (x3 - math_utils::computeCurvatureCentreX(x3, y3, x2, y2, x1, y1));
1235 q2pt = crv / (mBz * o2::constants::math::B2C);
1236 q2pt2 = crv * crv;
1237 }
1238 tgl12 = math_utils::computeTanDipAngle(x1, y1, x2, y2, z1, z2);
1239 tgl23 = math_utils::computeTanDipAngle(x2, y2, x3, y3, z2, z3);
1240 sg2q2pt = track::kC1Pt2max * (q2pt2 > 0.0005f ? (q2pt2 < 1.f ? q2pt2 : 1.f) : 0.0005f);
1241 return {tf3.xTrackingFrame, tf3.alphaTrackingFrame, {y3, z3, snp, 0.5f * (tgl12 + tgl23), q2pt}, {tf3.covarianceTrackingFrame[0], tf3.covarianceTrackingFrame[1], tf3.covarianceTrackingFrame[2], 0.f, 0.f, track::kCSnp2max, 0.f, 0.f, 0.f, track::kCTgl2max, 0.f, 0.f, 0.f, 0.f, sg2q2pt}};
1242}
1243
1244template <int nLayers>
1246{
1247 mBz = bz;
1248 mIsZeroField = std::abs(mBz) < 0.01;
1249 mTimeFrame->setBz(bz);
1250}
1251
1252template <int nLayers>
1254{
1255 return o2::base::Propagator::Instance()->getMatLUT() && (mTrkParams[0].CorrType == o2::base::PropagatorImpl<float>::MatCorrType::USEMatCorrLUT);
1256}
1257
1258template <int nLayers>
1259void TrackerTraits<nLayers>::setNThreads(int n, std::shared_ptr<tbb::task_arena>& arena)
1260{
1261#if defined(OPTIMISATION_OUTPUT) || defined(CA_DEBUG)
1262 mTaskArena = std::make_shared<tbb::task_arena>(1);
1263#else
1264 if (arena == nullptr) {
1265 mTaskArena = std::make_shared<tbb::task_arena>(std::abs(n));
1266 LOGP(info, "Setting tracker with {} threads.", n);
1267 } else {
1268 mTaskArena = arena;
1269 LOGP(info, "Attaching tracker to calling thread's arena");
1270 }
1271#endif
1272}
1273
1274template class TrackerTraits<7>;
1275
1276} // namespace o2::its
Base track model for the Barrel, params only, w/o covariance.
#define CA_DEBUGGER(x)
Definition Definitions.h:23
const auto bins
Definition PID.cxx:49
#define failed(...)
Definition Utils.h:42
useful math constants
GPUd() value_type estimateLTFast(o2 static GPUd() float estimateLTIncrement(const o2 PropagatorImpl * Instance(bool uninitialized=false)
Definition Propagator.h:143
HMPID cluster implementation.
Definition Cluster.h:27
virtual void processNeighbours(int iLayer, int iLevel, const bounded_vector< CellSeedN > &currentCellSeed, const bounded_vector< int > &currentCellId, bounded_vector< CellSeedN > &updatedCellSeed, bounded_vector< int > &updatedCellId)
virtual void findRoads(const int iteration)
void setNThreads(int n, std::shared_ptr< tbb::task_arena > &arena)
virtual bool trackFollowing(TrackITSExt *track, int rof, bool outward, const int iteration)
virtual void computeLayerTracklets(const int iteration, int iROFslice, int iVertex)
virtual void computeLayerCells(const int iteration)
virtual void findShortPrimaries()
virtual void setBz(float bz)
virtual void findCellsNeighbours(const int iteration)
virtual void extendTracks(const int iteration)
GLdouble n
Definition glcorearb.h:1982
GLint GLenum GLint x
Definition glcorearb.h:403
GLuint GLfloat GLfloat GLfloat GLfloat y1
Definition glcorearb.h:5034
GLsizeiptr size
Definition glcorearb.h:659
GLuint GLuint end
Definition glcorearb.h:469
GLuint GLfloat GLfloat GLfloat x1
Definition glcorearb.h:5034
GLboolean GLboolean GLboolean b
Definition glcorearb.h:1233
GLsizei const GLfloat * value
Definition glcorearb.h:819
GLintptr offset
Definition glcorearb.h:660
GLuint GLsizei const GLchar * label
Definition glcorearb.h:2519
GLboolean r
Definition glcorearb.h:1233
GLuint start
Definition glcorearb.h:469
GLboolean GLboolean GLboolean GLboolean a
Definition glcorearb.h:1233
GLdouble GLdouble GLdouble z
Definition glcorearb.h:843
constexpr float B2C
constexpr float TwoPI
float float float float float y3
Definition MathUtils.h:45
float float float float x3
Definition MathUtils.h:44
void deepVectorClear(std::vector< T > &vec)
std::pmr::vector< T > bounded_vector
constexpr Int_t TrackletsPerRoad
Definition Constants.h:39
Enum< T >::Iterator begin(Enum< T >)
Definition Defs.h:173
constexpr float kCTgl2max
TrackParCovF TrackParCov
Definition Track.h:33
constexpr float kCSnp2max
constexpr float kMostProbablePt
constexpr float kC1Pt2max
int32_t w
float yCoordinate
Definition Cluster.h:46
float zCoordinate
Definition Cluster.h:47
float xCoordinate
Definition Cluster.h:45
std::integral_constant< int, 0 > OnePass
std::integral_constant< int, 1 > TwoPassCount
std::integral_constant< int, 2 > TwoPassInsert
std::array< float, 2 > positionTrackingFrame
Definition Cluster.h:73
std::array< float, 3 > covarianceTrackingFrame
Definition Cluster.h:74
std::array< uint16_t, 5 > pattern
std::vector< Tracklet64 > tracklets