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};
52template <int nLayers>
53void TrackerTraits<nLayers>::computeLayerTracklets(const int iteration, int iROFslice, int iVertex)
54{
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 int minROF = o2::gpu::CAMath::Max(startROF, pivotROF - mTrkParams[iteration].DeltaROF);
77 int maxROF = o2::gpu::CAMath::Min(endROF - 1, pivotROF + mTrkParams[iteration].DeltaROF);
78 gsl::span<const Vertex> primaryVertices = mTrkParams[iteration].UseDiamond ? diamondSpan : mTimeFrame->getPrimaryVertices(minROF, maxROF);
79 if (primaryVertices.empty()) {
80 return 0;
81 }
82 const int startVtx = iVertex >= 0 ? iVertex : 0;
83 const int endVtx = iVertex >= 0 ? o2::gpu::CAMath::Min(iVertex + 1, int(primaryVertices.size())) : int(primaryVertices.size());
84 if (endVtx <= startVtx) {
85 return 0;
86 }
87
88 int localCount = 0;
89 auto& tracklets = mTimeFrame->getTracklets()[iLayer];
90 for (int targetROF0{minROF}; targetROF0 <= maxROF; ++targetROF0) {
91 if (!mTimeFrame->mMultiplicityCutMask[targetROF0]) {
92 continue;
93 }
94 auto layer0 = mTimeFrame->getClustersOnLayer(targetROF0, iLayer);
95 if (layer0.empty()) {
96 continue;
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(targetROF0, 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 targetROF1{minROF}; targetROF1 <= maxROF; ++targetROF1) {
132 if (!mTimeFrame->mMultiplicityCutMask[targetROF1] || std::abs(targetROF0 - targetROF1) > mTrkParams[iteration].DeltaROF) {
133 continue;
134 }
135 auto layer1 = mTimeFrame->getClustersOnLayer(targetROF1, iLayer + 1);
136 if (layer1.empty()) {
137 continue;
138 }
139 for (int iPhi = 0; iPhi < phiBinsNum; ++iPhi) {
140 int iPhiBin = (bins.y + iPhi) % mTrkParams[iteration].PhiBins;
141 int firstBinIdx = mTimeFrame->mIndexTableUtils.getBinIndex(bins.x, iPhiBin);
142 int maxBinIdx = firstBinIdx + (bins.z - bins.x) + 1;
143 int firstRow = mTimeFrame->getIndexTable(targetROF1, iLayer + 1)[firstBinIdx];
144 int lastRow = mTimeFrame->getIndexTable(targetROF1, 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) ||
176 o2::gpu::GPUCommonMath::Abs(deltaPhi - o2::constants::math::TwoPI) < mTimeFrame->getPhiCut(iLayer))) {
177 float phi = o2::gpu::GPUCommonMath::ATan2(currentCluster.yCoordinate - nextCluster.yCoordinate, currentCluster.xCoordinate - nextCluster.xCoordinate);
178 float tanL = (currentCluster.zCoordinate - nextCluster.zCoordinate) / (currentCluster.radius - nextCluster.radius);
179 if constexpr (decltype(Tag)::value == PassMode::OnePass::value) {
180 tracklets.emplace_back(currentSortedIndex, mTimeFrame->getSortedIndex(targetROF1, iLayer + 1, iNext), tanL, phi, targetROF0, targetROF1);
181 } else if constexpr (decltype(Tag)::value == PassMode::TwoPassCount::value) {
182 ++localCount;
183 } else if constexpr (decltype(Tag)::value == PassMode::TwoPassInsert::value) {
184 const int idx = base + offset++;
185 tracklets[idx] = Tracklet(currentSortedIndex, mTimeFrame->getSortedIndex(targetROF1, iLayer + 1, iNext), tanL, phi, targetROF0, targetROF1);
186 }
187 }
188 }
189 }
190 }
191 }
192 }
193 }
194 return localCount;
195 };
196
197 int dummy{0};
198 if (mTaskArena->max_concurrency() <= 1) {
199 for (int pivotROF{startROF}; pivotROF < endROF; ++pivotROF) {
200 for (int iLayer{0}; iLayer < mTrkParams[iteration].TrackletsPerRoad(); ++iLayer) {
201 forTracklets(PassMode::OnePass{}, iLayer, pivotROF, 0, dummy);
202 }
203 }
204 } else {
205 bounded_vector<bounded_vector<int>> perROFCount(mTrkParams[iteration].TrackletsPerRoad(), bounded_vector<int>(endROF - startROF + 1, 0, mMemoryPool.get()), mMemoryPool.get());
206 tbb::parallel_for(
207 tbb::blocked_range2d<int, int>(0, mTrkParams[iteration].TrackletsPerRoad(), 1,
208 startROF, endROF, 1),
209 [&](auto const& Range) {
210 for (int iLayer{Range.rows().begin()}; iLayer < Range.rows().end(); ++iLayer) {
211 for (int pivotROF = Range.cols().begin(); pivotROF < Range.cols().end(); ++pivotROF) {
212 perROFCount[iLayer][pivotROF - startROF] = forTracklets(PassMode::TwoPassCount{}, iLayer, pivotROF, 0, dummy);
213 }
214 }
215 });
216
217 tbb::parallel_for(
218 tbb::blocked_range<int>(0, mTrkParams[iteration].TrackletsPerRoad()),
219 [&](auto const& Layers) {
220 for (int iLayer{Layers.begin()}; iLayer < Layers.end(); ++iLayer) {
221 std::exclusive_scan(perROFCount[iLayer].begin(), perROFCount[iLayer].end(), perROFCount[iLayer].begin(), 0);
222 mTimeFrame->getTracklets()[iLayer].resize(perROFCount[iLayer].back());
223 }
224 });
225
226 tbb::parallel_for(
227 tbb::blocked_range2d<int, int>(0, mTrkParams[iteration].TrackletsPerRoad(), 1,
228 startROF, endROF, 1),
229 [&](auto const& Range) {
230 for (int iLayer{Range.rows().begin()}; iLayer < Range.rows().end(); ++iLayer) {
231 if (perROFCount[iLayer].back() == 0) {
232 continue;
233 }
234 for (int pivotROF = Range.cols().begin(); pivotROF < Range.cols().end(); ++pivotROF) {
235 int baseIdx = perROFCount[iLayer][pivotROF - startROF];
236 if (baseIdx == perROFCount[iLayer][pivotROF - startROF + 1]) {
237 continue;
238 }
239 int localIdx = 0;
240 forTracklets(PassMode::TwoPassInsert{}, iLayer, pivotROF, baseIdx, localIdx);
241 }
242 }
243 });
244 }
245
246 tbb::parallel_for(
247 tbb::blocked_range<int>(0, mTrkParams[iteration].TrackletsPerRoad()),
248 [&](const tbb::blocked_range<int>& Layers) {
249 for (int iLayer = Layers.begin(); iLayer < Layers.end(); ++iLayer) {
251 auto& trkl{mTimeFrame->getTracklets()[iLayer]};
252 tbb::parallel_sort(trkl.begin(), trkl.end(), [](const Tracklet& a, const Tracklet& b) -> bool {
253 return a.firstClusterIndex < b.firstClusterIndex || (a.firstClusterIndex == b.firstClusterIndex && a.secondClusterIndex < b.secondClusterIndex);
254 });
256 trkl.erase(std::unique(trkl.begin(), trkl.end(), [](const Tracklet& a, const Tracklet& b) -> bool {
257 return a.firstClusterIndex == b.firstClusterIndex && a.secondClusterIndex == b.secondClusterIndex;
258 }),
259 trkl.end());
260 trkl.shrink_to_fit();
261 if (iLayer > 0) {
262 auto& lut{mTimeFrame->getTrackletsLookupTable()[iLayer - 1]};
263 if (!trkl.empty()) {
264 for (const auto& tkl : trkl) {
265 lut[tkl.firstClusterIndex + 1]++;
266 }
267 std::inclusive_scan(lut.begin(), lut.end(), lut.begin());
268 }
269 }
270 }
271 });
272
274 if (mTimeFrame->hasMCinformation()) {
275 tbb::parallel_for(
276 tbb::blocked_range<int>(0, mTrkParams[iteration].TrackletsPerRoad()),
277 [&](const tbb::blocked_range<int>& Layers) {
278 for (int iLayer = Layers.begin(); iLayer < Layers.end(); ++iLayer) {
279 for (auto& trk : mTimeFrame->getTracklets()[iLayer]) {
280 MCCompLabel label;
281 int currentId{mTimeFrame->getClusters()[iLayer][trk.firstClusterIndex].clusterId};
282 int nextId{mTimeFrame->getClusters()[iLayer + 1][trk.secondClusterIndex].clusterId};
283 for (const auto& lab1 : mTimeFrame->getClusterLabels(iLayer, currentId)) {
284 for (const auto& lab2 : mTimeFrame->getClusterLabels(iLayer + 1, nextId)) {
285 if (lab1 == lab2 && lab1.isValid()) {
286 label = lab1;
287 break;
288 }
289 }
290 if (label.isValid()) {
291 break;
292 }
293 }
294 mTimeFrame->getTrackletsLabel(iLayer).emplace_back(label);
295 }
296 }
297 });
298 }
299 });
300}
301
302template <int nLayers>
304{
305#ifdef OPTIMISATION_OUTPUT
306 static int iter{0};
307 std::ofstream off(std::format("cells{}.txt", iter++));
308#endif
309
310 for (int iLayer = 0; iLayer < mTrkParams[iteration].CellsPerRoad(); ++iLayer) {
311 deepVectorClear(mTimeFrame->getCells()[iLayer]);
312 if (iLayer > 0) {
313 deepVectorClear(mTimeFrame->getCellsLookupTable()[iLayer - 1]);
314 }
315 if (mTimeFrame->hasMCinformation()) {
316 deepVectorClear(mTimeFrame->getCellsLabel(iLayer));
317 }
318 }
319
320 mTaskArena->execute([&] {
321 auto forTrackletCells = [&](auto Tag, int iLayer, bounded_vector<CellSeed>& layerCells, int iTracklet, int offset = 0) -> int {
322 const Tracklet& currentTracklet{mTimeFrame->getTracklets()[iLayer][iTracklet]};
323 const int nextLayerClusterIndex{currentTracklet.secondClusterIndex};
324 const int nextLayerFirstTrackletIndex{mTimeFrame->getTrackletsLookupTable()[iLayer][nextLayerClusterIndex]};
325 const int nextLayerLastTrackletIndex{mTimeFrame->getTrackletsLookupTable()[iLayer][nextLayerClusterIndex + 1]};
326 int foundCells{0};
327 for (int iNextTracklet{nextLayerFirstTrackletIndex}; iNextTracklet < nextLayerLastTrackletIndex; ++iNextTracklet) {
328 const Tracklet& nextTracklet{mTimeFrame->getTracklets()[iLayer + 1][iNextTracklet]};
329 const auto& nextLbl = mTimeFrame->getTrackletsLabel(iLayer + 1)[iNextTracklet];
330 bool print = false;
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++] = CellSeed(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()) {
462 for (int iLayer{0}; iLayer < mTrkParams[iteration].CellsPerRoad(); ++iLayer) {
463 for (const auto& cell : mTimeFrame->getCells()[iLayer]) {
464 MCCompLabel currentLab{mTimeFrame->getTrackletsLabel(iLayer)[cell.getFirstTrackletIndex()]};
465 MCCompLabel nextLab{mTimeFrame->getTrackletsLabel(iLayer + 1)[cell.getSecondTrackletIndex()]};
466 mTimeFrame->getCellsLabel(iLayer).emplace_back(currentLab == nextLab ? currentLab : MCCompLabel());
467 }
468 }
469 }
470}
471
472template <int nLayers>
474{
475#ifdef OPTIMISATION_OUTPUT
476 std::ofstream off(std::format("cellneighs{}.txt", iteration));
477#endif
478
479 struct Neighbor {
480 int cell{-1}, nextCell{-1}, level{-1};
481 };
482
483 mTaskArena->execute([&] {
484 for (int iLayer{0}; iLayer < mTrkParams[iteration].CellsPerRoad() - 1; ++iLayer) {
485 deepVectorClear(mTimeFrame->getCellsNeighbours()[iLayer]);
486 deepVectorClear(mTimeFrame->getCellsNeighboursLUT()[iLayer]);
487 if (mTimeFrame->getCells()[iLayer + 1].empty() ||
488 mTimeFrame->getCellsLookupTable()[iLayer].empty()) {
489 continue;
490 }
491
492 int nCells{static_cast<int>(mTimeFrame->getCells()[iLayer].size())};
493 bounded_vector<Neighbor> cellsNeighbours(mMemoryPool.get());
494
495 auto forCellNeighbour = [&](auto Tag, int iCell, int offset = 0) -> int {
496 const auto& currentCellSeed{mTimeFrame->getCells()[iLayer][iCell]};
497 const int nextLayerTrackletIndex{currentCellSeed.getSecondTrackletIndex()};
498 const int nextLayerFirstCellIndex{mTimeFrame->getCellsLookupTable()[iLayer][nextLayerTrackletIndex]};
499 const int nextLayerLastCellIndex{mTimeFrame->getCellsLookupTable()[iLayer][nextLayerTrackletIndex + 1]};
500 int foundNextCells{0};
501 for (int iNextCell{nextLayerFirstCellIndex}; iNextCell < nextLayerLastCellIndex; ++iNextCell) {
502 auto nextCellSeed{mTimeFrame->getCells()[iLayer + 1][iNextCell]};
503 if (nextCellSeed.getFirstTrackletIndex() != nextLayerTrackletIndex) {
504 break;
505 }
506
507 if (mTrkParams[iteration].DeltaROF) { // TODO this has to be improved for the staggering
508 const auto& trkl00 = mTimeFrame->getTracklets()[iLayer][currentCellSeed.getFirstTrackletIndex()];
509 const auto& trkl01 = mTimeFrame->getTracklets()[iLayer + 1][currentCellSeed.getSecondTrackletIndex()];
510 const auto& trkl10 = mTimeFrame->getTracklets()[iLayer + 1][nextCellSeed.getFirstTrackletIndex()];
511 const auto& trkl11 = mTimeFrame->getTracklets()[iLayer + 2][nextCellSeed.getSecondTrackletIndex()];
512 if ((std::max({trkl00.getMaxRof(), trkl01.getMaxRof(), trkl10.getMaxRof(), trkl11.getMaxRof()}) - std::min({trkl00.getMinRof(), trkl01.getMinRof(), trkl10.getMinRof(), trkl10.getMinRof()})) > mTrkParams[0].DeltaROF) {
513 continue;
514 }
515 }
516
517 if (!nextCellSeed.rotate(currentCellSeed.getAlpha()) ||
518 !nextCellSeed.propagateTo(currentCellSeed.getX(), getBz())) {
519 continue;
520 }
521 float chi2 = currentCellSeed.getPredictedChi2(nextCellSeed);
522
523#ifdef OPTIMISATION_OUTPUT
524 bool good{mTimeFrame->getCellsLabel(iLayer)[iCell] == mTimeFrame->getCellsLabel(iLayer + 1)[iNextCell]};
525 off << std::format("{}\t{:d}\t{}", iLayer, good, chi2) << std::endl;
526#endif
527
528 if (chi2 > mTrkParams[0].MaxChi2ClusterAttachment) {
529 continue;
530 }
531
532 if constexpr (decltype(Tag)::value == PassMode::OnePass::value) {
533 cellsNeighbours.emplace_back(iCell, iNextCell, currentCellSeed.getLevel() + 1);
534 } else if constexpr (decltype(Tag)::value == PassMode::TwoPassCount::value) {
535 ++foundNextCells;
536 } else if constexpr (decltype(Tag)::value == PassMode::TwoPassInsert::value) {
537 cellsNeighbours[offset++] = {iCell, iNextCell, currentCellSeed.getLevel() + 1};
538 } else {
539 static_assert(false, "Unknown mode!");
540 }
541 }
542 return foundNextCells;
543 };
544
545 if (mTaskArena->max_concurrency() <= 1) {
546 for (int iCell{0}; iCell < nCells; ++iCell) {
547 forCellNeighbour(PassMode::OnePass{}, iCell);
548 }
549 } else {
550 bounded_vector<int> perCellCount(nCells + 1, 0, mMemoryPool.get());
551 tbb::parallel_for(
552 tbb::blocked_range<int>(0, nCells),
553 [&](const tbb::blocked_range<int>& Cells) {
554 for (int iCell = Cells.begin(); iCell < Cells.end(); ++iCell) {
555 perCellCount[iCell] = forCellNeighbour(PassMode::TwoPassCount{}, iCell);
556 }
557 });
558
559 std::exclusive_scan(perCellCount.begin(), perCellCount.end(), perCellCount.begin(), 0);
560 int totalCellNeighbours = perCellCount.back();
561 if (totalCellNeighbours == 0) {
562 deepVectorClear(mTimeFrame->getCellsNeighbours()[iLayer]);
563 continue;
564 }
565 cellsNeighbours.resize(totalCellNeighbours);
566
567 tbb::parallel_for(
568 tbb::blocked_range<int>(0, nCells),
569 [&](const tbb::blocked_range<int>& Cells) {
570 for (int iCell = Cells.begin(); iCell < Cells.end(); ++iCell) {
571 int offset = perCellCount[iCell];
572 if (offset == perCellCount[iCell + 1]) {
573 continue;
574 }
575 forCellNeighbour(PassMode::TwoPassInsert{}, iCell, offset);
576 }
577 });
578 }
579
580 if (cellsNeighbours.empty()) {
581 continue;
582 }
583
584 tbb::parallel_sort(cellsNeighbours.begin(), cellsNeighbours.end(), [](const auto& a, const auto& b) {
585 return a.nextCell < b.nextCell;
586 });
587
588 auto& cellsNeighbourLUT = mTimeFrame->getCellsNeighboursLUT()[iLayer];
589 cellsNeighbourLUT.assign(mTimeFrame->getCells()[iLayer + 1].size(), 0);
590 for (const auto& neigh : cellsNeighbours) {
591 ++cellsNeighbourLUT[neigh.nextCell];
592 }
593 std::inclusive_scan(cellsNeighbourLUT.begin(), cellsNeighbourLUT.end(), cellsNeighbourLUT.begin());
594
595 mTimeFrame->getCellsNeighbours()[iLayer].reserve(cellsNeighbours.size());
596 std::ranges::transform(cellsNeighbours, std::back_inserter(mTimeFrame->getCellsNeighbours()[iLayer]), [](const auto& neigh) { return neigh.cell; });
597
598 auto it = cellsNeighbours.begin();
599 int current = it->nextCell;
600 int maxLvl = it->level;
601 ++it;
602 for (; it != cellsNeighbours.end(); ++it) {
603 if (it->nextCell == current) {
604 maxLvl = std::max(maxLvl, it->level);
605 } else {
606 mTimeFrame->getCells()[iLayer + 1][current].setLevel(maxLvl);
607 current = it->nextCell;
608 maxLvl = it->level;
609 }
610 }
611 mTimeFrame->getCells()[iLayer + 1][current].setLevel(maxLvl);
612 }
613 });
614}
615
616template <int nLayers>
617void TrackerTraits<nLayers>::processNeighbours(int iLayer, int iLevel, const bounded_vector<CellSeed>& currentCellSeed, const bounded_vector<int>& currentCellId, bounded_vector<CellSeed>& 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 CellSeed& 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 CellSeed& 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 if (mTrkParams[0].DeltaROF) { // TODO this has to be improved for the staggering
661 const auto& trklNeigh = mTimeFrame->getTracklets()[iLayer - 1][neighbourCell.getFirstTrackletIndex()];
662 short minRof{std::numeric_limits<short>::max()}, maxRof{std::numeric_limits<short>::min()};
663 for (int iLayer{0}; iLayer < mTrkParams[0].NLayers; ++iLayer) {
664 if (const auto clsId = currentCell.getCluster(iLayer); clsId != constants::UnusedIndex) {
665 const short clsROF = mTimeFrame->getClusterROF(iLayer, clsId);
666 minRof = std::min(minRof, clsROF);
667 maxRof = std::max(maxRof, clsROF);
668 }
669 }
670 if ((std::max(trklNeigh.getMaxRof(), maxRof) - std::min(trklNeigh.getMinRof(), minRof)) > mTrkParams[0].DeltaROF) {
671 continue;
672 }
673 }
674
676 CellSeed seed{currentCell};
677 const auto& trHit = mTimeFrame->getTrackingFrameInfoOnLayer(iLayer - 1)[neighbourCell.getFirstClusterIndex()];
678
679 if (!seed.rotate(trHit.alphaTrackingFrame)) {
680 CA_DEBUGGER(failed[1]++);
681 continue;
682 }
683
684 if (!propagator->propagateToX(seed, trHit.xTrackingFrame, getBz(), o2::base::PropagatorImpl<float>::MAX_SIN_PHI, o2::base::PropagatorImpl<float>::MAX_STEP, mTrkParams[0].CorrType)) {
685 CA_DEBUGGER(failed[2]++);
686 continue;
687 }
688
689 if (mTrkParams[0].CorrType == o2::base::PropagatorF::MatCorrType::USEMatCorrNONE) {
690 if (!seed.correctForMaterial(mTrkParams[0].LayerxX0[iLayer - 1], mTrkParams[0].LayerxX0[iLayer - 1] * constants::Radl * constants::Rho, true)) {
691 continue;
692 }
693 }
694
695 auto predChi2{seed.getPredictedChi2Quiet(trHit.positionTrackingFrame, trHit.covarianceTrackingFrame)};
696 if ((predChi2 > mTrkParams[0].MaxChi2ClusterAttachment) || predChi2 < 0.f) {
697 CA_DEBUGGER(failed[3]++);
698 continue;
699 }
700 seed.setChi2(seed.getChi2() + predChi2);
701 if (!seed.o2::track::TrackParCov::update(trHit.positionTrackingFrame, trHit.covarianceTrackingFrame)) {
702 CA_DEBUGGER(failed[4]++);
703 continue;
704 }
705
706 if constexpr (decltype(Tag)::value != PassMode::TwoPassCount::value) {
707 seed.getClusters()[iLayer - 1] = neighbourCell.getFirstClusterIndex();
708 seed.setLevel(neighbourCell.getLevel());
709 seed.setFirstTrackletIndex(neighbourCell.getFirstTrackletIndex());
710 seed.setSecondTrackletIndex(neighbourCell.getSecondTrackletIndex());
711 }
712
713 if constexpr (decltype(Tag)::value == PassMode::OnePass::value) {
714 updatedCellSeeds.push_back(seed);
715 updatedCellsIds.push_back(neighbourCellId);
716 } else if constexpr (decltype(Tag)::value == PassMode::TwoPassCount::value) {
717 ++foundSeeds;
718 } else if constexpr (decltype(Tag)::value == PassMode::TwoPassInsert::value) {
719 updatedCellSeeds[offset] = seed;
720 updatedCellsIds[offset++] = neighbourCellId;
721 } else {
722 static_assert(false, "Unknown mode!");
723 }
724 }
725 return foundSeeds;
726 };
727
728 const int nCells = static_cast<int>(currentCellSeed.size());
729 if (mTaskArena->max_concurrency() <= 1) {
730 for (int iCell{0}; iCell < nCells; ++iCell) {
731 forCellNeighbours(PassMode::OnePass{}, iCell);
732 }
733 } else {
734 bounded_vector<int> perCellCount(nCells + 1, 0, mMemoryPool.get());
735 tbb::parallel_for(
736 tbb::blocked_range<int>(0, nCells),
737 [&](const tbb::blocked_range<int>& Cells) {
738 for (int iCell = Cells.begin(); iCell < Cells.end(); ++iCell) {
739 perCellCount[iCell] = forCellNeighbours(PassMode::TwoPassCount{}, iCell);
740 }
741 });
742
743 std::exclusive_scan(perCellCount.begin(), perCellCount.end(), perCellCount.begin(), 0);
744 auto totalNeighbours{perCellCount.back()};
745 if (totalNeighbours == 0) {
746 return;
747 }
748 updatedCellSeeds.resize(totalNeighbours);
749 updatedCellsIds.resize(totalNeighbours);
750
751 tbb::parallel_for(
752 tbb::blocked_range<int>(0, nCells),
753 [&](const tbb::blocked_range<int>& Cells) {
754 for (int iCell = Cells.begin(); iCell < Cells.end(); ++iCell) {
755 int offset = perCellCount[iCell];
756 if (offset == perCellCount[iCell + 1]) {
757 continue;
758 }
759 forCellNeighbours(PassMode::TwoPassInsert{}, iCell, offset);
760 }
761 });
762 }
763 });
764
765#ifdef CA_DEBUG
766 std::cout << "\t\t- Found " << updatedCellSeeds.size() << " cell seeds out of " << attempts << " attempts" << std::endl;
767 std::cout << "\t\t\t> " << failed[0] << " failed because of level" << std::endl;
768 std::cout << "\t\t\t> " << failed[1] << " failed because of rotation" << std::endl;
769 std::cout << "\t\t\t> " << failed[2] << " failed because of propagation" << std::endl;
770 std::cout << "\t\t\t> " << failed[3] << " failed because of chi2 cut" << std::endl;
771 std::cout << "\t\t\t> " << failed[4] << " failed because of update" << std::endl;
772 std::cout << "\t\t\t> " << failedByMismatch << " failed because of mismatch" << std::endl;
773#endif
774}
775
776template <int nLayers>
777void TrackerTraits<nLayers>::findRoads(const int iteration)
778{
779 CA_DEBUGGER(std::cout << "Finding roads, iteration " << iteration << std::endl);
780
781 for (int startLevel{mTrkParams[iteration].CellsPerRoad()}; startLevel >= mTrkParams[iteration].CellMinimumLevel(); --startLevel) {
782 CA_DEBUGGER(std::cout << "\t > Processing level " << startLevel << std::endl);
783 auto seedFilter = [&](const CellSeed& seed) {
784 return seed.getQ2Pt() <= 1.e3 && seed.getChi2() <= mTrkParams[0].MaxChi2NDF * ((startLevel + 2) * 2 - 5);
785 };
786 bounded_vector<CellSeed> trackSeeds(mMemoryPool.get());
787 for (int startLayer{mTrkParams[iteration].CellsPerRoad() - 1}; startLayer >= startLevel - 1; --startLayer) {
788 if ((mTrkParams[iteration].StartLayerMask & (1 << (startLayer + 2))) == 0) {
789 continue;
790 }
791 CA_DEBUGGER(std::cout << "\t\t > Starting processing layer " << startLayer << std::endl);
792 bounded_vector<int> lastCellId(mMemoryPool.get()), updatedCellId(mMemoryPool.get());
793 bounded_vector<CellSeed> lastCellSeed(mMemoryPool.get()), updatedCellSeed(mMemoryPool.get());
794
795 processNeighbours(startLayer, startLevel, mTimeFrame->getCells()[startLayer], lastCellId, updatedCellSeed, updatedCellId);
796
797 int level = startLevel;
798 for (int iLayer{startLayer - 1}; iLayer > 0 && level > 2; --iLayer) {
799 lastCellSeed.swap(updatedCellSeed);
800 lastCellId.swap(updatedCellId);
801 deepVectorClear(updatedCellSeed);
802 deepVectorClear(updatedCellId);
803 processNeighbours(iLayer, --level, lastCellSeed, lastCellId, updatedCellSeed, updatedCellId);
804 }
805 deepVectorClear(lastCellId);
806 deepVectorClear(lastCellSeed);
807
808 if (!updatedCellSeed.empty()) {
809 trackSeeds.reserve(trackSeeds.size() + std::count_if(updatedCellSeed.begin(), updatedCellSeed.end(), seedFilter));
810 std::copy_if(updatedCellSeed.begin(), updatedCellSeed.end(), std::back_inserter(trackSeeds), seedFilter);
811 }
812 }
813
814 if (trackSeeds.empty()) {
815 continue;
816 }
817
818 bounded_vector<TrackITSExt> tracks(mMemoryPool.get());
819 mTaskArena->execute([&] {
820 auto forSeed = [&](auto Tag, int iSeed, int offset = 0) {
821 const CellSeed& seed{trackSeeds[iSeed]};
822 TrackITSExt temporaryTrack{seed};
823 temporaryTrack.resetCovariance();
824 temporaryTrack.setChi2(0);
825 for (int iL{0}; iL < 7; ++iL) {
826 temporaryTrack.setExternalClusterIndex(iL, seed.getCluster(iL), seed.getCluster(iL) != constants::UnusedIndex);
827 }
828
829 bool fitSuccess = fitTrack(temporaryTrack, 0, mTrkParams[0].NLayers, 1, mTrkParams[0].MaxChi2ClusterAttachment, mTrkParams[0].MaxChi2NDF);
830 if (!fitSuccess) {
831 return 0;
832 }
833
834 temporaryTrack.getParamOut() = temporaryTrack.getParamIn();
835 temporaryTrack.resetCovariance();
836 temporaryTrack.setChi2(0);
837 fitSuccess = fitTrack(temporaryTrack, mTrkParams[0].NLayers - 1, -1, -1, mTrkParams[0].MaxChi2ClusterAttachment, mTrkParams[0].MaxChi2NDF, 50.f);
838 if (!fitSuccess || temporaryTrack.getPt() < mTrkParams[iteration].MinPt[mTrkParams[iteration].NLayers - temporaryTrack.getNClusters()]) {
839 return 0;
840 }
841
842 if constexpr (decltype(Tag)::value == PassMode::OnePass::value) {
843 tracks.push_back(temporaryTrack);
844 } else if constexpr (decltype(Tag)::value == PassMode::TwoPassCount::value) {
845 // nothing to do
846 } else if constexpr (decltype(Tag)::value == PassMode::TwoPassInsert::value) {
847 tracks[offset] = temporaryTrack;
848 } else {
849 static_assert(false, "Unknown mode!");
850 }
851 return 1;
852 };
853
854 const int nSeeds = static_cast<int>(trackSeeds.size());
855 if (mTaskArena->max_concurrency() <= 1) {
856 for (int iSeed{0}; iSeed < nSeeds; ++iSeed) {
857 forSeed(PassMode::OnePass{}, iSeed);
858 }
859 } else {
860 bounded_vector<int> perSeedCount(nSeeds + 1, 0, mMemoryPool.get());
861 tbb::parallel_for(
862 tbb::blocked_range<int>(0, nSeeds),
863 [&](const tbb::blocked_range<int>& Seeds) {
864 for (int iSeed = Seeds.begin(); iSeed < Seeds.end(); ++iSeed) {
865 perSeedCount[iSeed] = forSeed(PassMode::TwoPassCount{}, iSeed);
866 }
867 });
868
869 std::exclusive_scan(perSeedCount.begin(), perSeedCount.end(), perSeedCount.begin(), 0);
870 auto totalTracks{perSeedCount.back()};
871 if (totalTracks == 0) {
872 return;
873 }
874 tracks.resize(totalTracks);
875
876 tbb::parallel_for(
877 tbb::blocked_range<int>(0, nSeeds),
878 [&](const tbb::blocked_range<int>& Seeds) {
879 for (int iSeed = Seeds.begin(); iSeed < Seeds.end(); ++iSeed) {
880 if (perSeedCount[iSeed] == perSeedCount[iSeed + 1]) {
881 continue;
882 }
883 forSeed(PassMode::TwoPassInsert{}, iSeed, perSeedCount[iSeed]);
884 }
885 });
886 }
887
888 deepVectorClear(trackSeeds);
889 tbb::parallel_sort(tracks.begin(), tracks.end(), [](const auto& a, const auto& b) {
890 return a.getChi2() < b.getChi2();
891 });
892 });
893
894 for (auto& track : tracks) {
895 int nShared = 0;
896 bool isFirstShared{false};
897 for (int iLayer{0}; iLayer < mTrkParams[0].NLayers; ++iLayer) {
898 if (track.getClusterIndex(iLayer) == constants::UnusedIndex) {
899 continue;
900 }
901 nShared += int(mTimeFrame->isClusterUsed(iLayer, track.getClusterIndex(iLayer)));
902 isFirstShared |= !iLayer && mTimeFrame->isClusterUsed(iLayer, track.getClusterIndex(iLayer));
903 }
904
905 if (nShared > mTrkParams[0].ClusterSharing) {
906 continue;
907 }
908
909 std::array<int, 3> rofs{INT_MAX, INT_MAX, INT_MAX};
910 for (int iLayer{0}; iLayer < mTrkParams[0].NLayers; ++iLayer) {
911 if (track.getClusterIndex(iLayer) == constants::UnusedIndex) {
912 continue;
913 }
914 mTimeFrame->markUsedCluster(iLayer, track.getClusterIndex(iLayer));
915 int currentROF = mTimeFrame->getClusterROF(iLayer, track.getClusterIndex(iLayer));
916 for (int iR{0}; iR < 3; ++iR) {
917 if (rofs[iR] == INT_MAX) {
918 rofs[iR] = currentROF;
919 }
920 if (rofs[iR] == currentROF) {
921 break;
922 }
923 }
924 }
925 if (rofs[2] != INT_MAX) {
926 continue;
927 }
928 track.setUserField(0);
929 track.getParamOut().setUserField(0);
930 if (rofs[1] != INT_MAX) {
931 track.setNextROFbit();
932 }
933 mTimeFrame->getTracks(o2::gpu::CAMath::Min(rofs[0], rofs[1])).emplace_back(track);
934 }
935 }
936}
937
938template <int nLayers>
940{
941 for (int rof{0}; rof < mTimeFrame->getNrof(); ++rof) {
942 for (auto& track : mTimeFrame->getTracks(rof)) {
943 auto backup{track};
944 bool success{false};
945 // the order here biases towards top extension, tracks should probably be fitted separately in the directions and then compared.
946 if ((mTrkParams[iteration].UseTrackFollowerMix || mTrkParams[iteration].UseTrackFollowerTop) && track.getLastClusterLayer() != mTrkParams[iteration].NLayers - 1) {
947 success = success || trackFollowing(&track, rof, true, iteration);
948 }
949 if ((mTrkParams[iteration].UseTrackFollowerMix || (mTrkParams[iteration].UseTrackFollowerBot && !success)) && track.getFirstClusterLayer() != 0) {
950 success = success || trackFollowing(&track, rof, false, iteration);
951 }
952 if (success) {
954 track.resetCovariance();
955 track.setChi2(0);
956 bool fitSuccess = fitTrack(track, 0, mTrkParams[iteration].NLayers, 1, mTrkParams[iteration].MaxChi2ClusterAttachment, mTrkParams[0].MaxChi2NDF);
957 if (!fitSuccess) {
958 track = backup;
959 continue;
960 }
961 track.getParamOut() = track;
962 track.resetCovariance();
963 track.setChi2(0);
964 fitSuccess = fitTrack(track, mTrkParams[iteration].NLayers - 1, -1, -1, mTrkParams[iteration].MaxChi2ClusterAttachment, mTrkParams[0].MaxChi2NDF, 50.);
965 if (!fitSuccess) {
966 track = backup;
967 continue;
968 }
969 mTimeFrame->mNExtendedTracks++;
970 mTimeFrame->mNExtendedUsedClusters += track.getNClusters() - backup.getNClusters();
971 auto pattern = track.getPattern();
972 auto diff = (pattern & ~backup.getPattern()) & 0xff;
973 pattern |= (diff << 24);
974 track.setPattern(pattern);
976 for (int iLayer{0}; iLayer < mTrkParams[iteration].NLayers; ++iLayer) {
977 if (track.getClusterIndex(iLayer) == constants::UnusedIndex) {
978 continue;
979 }
980 mTimeFrame->markUsedCluster(iLayer, track.getClusterIndex(iLayer));
981 }
982 }
983 }
984 }
985}
986
987template <int nLayers>
989{
990 const auto propagator = o2::base::Propagator::Instance();
991 mTimeFrame->fillPrimaryVerticesXandAlpha();
992
993 for (auto& cell : mTimeFrame->getCells()[0]) {
994 auto& cluster3_glo = mTimeFrame->getClusters()[2][cell.getThirdClusterIndex()];
995 auto& cluster2_glo = mTimeFrame->getClusters()[1][cell.getSecondClusterIndex()];
996 auto& cluster1_glo = mTimeFrame->getClusters()[0][cell.getFirstClusterIndex()];
997 if (mTimeFrame->isClusterUsed(2, cluster1_glo.clusterId) ||
998 mTimeFrame->isClusterUsed(1, cluster2_glo.clusterId) ||
999 mTimeFrame->isClusterUsed(0, cluster3_glo.clusterId)) {
1000 continue;
1001 }
1002
1003 std::array<int, 3> rofs{
1004 mTimeFrame->getClusterROF(2, cluster3_glo.clusterId),
1005 mTimeFrame->getClusterROF(1, cluster2_glo.clusterId),
1006 mTimeFrame->getClusterROF(0, cluster1_glo.clusterId)};
1007 if (rofs[0] != rofs[1] && rofs[1] != rofs[2] && rofs[0] != rofs[2]) {
1008 continue;
1009 }
1010
1011 int rof{rofs[0]};
1012 if (rofs[1] == rofs[2]) {
1013 rof = rofs[2];
1014 }
1015
1016 auto pvs{mTimeFrame->getPrimaryVertices(rof)};
1017 auto pvsXAlpha{mTimeFrame->getPrimaryVerticesXAlpha(rof)};
1018
1019 const auto& cluster3_tf = mTimeFrame->getTrackingFrameInfoOnLayer(2)[cluster3_glo.clusterId];
1020 TrackITSExt temporaryTrack{buildTrackSeed(cluster1_glo, cluster2_glo, cluster3_tf)};
1021 temporaryTrack.setExternalClusterIndex(0, cluster1_glo.clusterId, true);
1022 temporaryTrack.setExternalClusterIndex(1, cluster2_glo.clusterId, true);
1023 temporaryTrack.setExternalClusterIndex(2, cluster3_glo.clusterId, true);
1024
1026 bool fitSuccess = fitTrack(temporaryTrack, 1, -1, -1);
1027 if (!fitSuccess) {
1028 continue;
1029 }
1030 fitSuccess = false;
1031
1032 TrackITSExt bestTrack{temporaryTrack}, backup{temporaryTrack};
1033 float bestChi2{std::numeric_limits<float>::max()};
1034 for (int iV{0}; iV < (int)pvs.size(); ++iV) {
1035 temporaryTrack = backup;
1036 if (!temporaryTrack.rotate(pvsXAlpha[iV][1])) {
1037 continue;
1038 }
1039 if (!propagator->propagateTo(temporaryTrack, pvsXAlpha[iV][0], true)) {
1040 continue;
1041 }
1042
1043 float pvRes{mTrkParams[0].PVres / o2::gpu::CAMath::Sqrt(float(pvs[iV].getNContributors()))};
1044 const float posVtx[2]{0.f, pvs[iV].getZ()};
1045 const float covVtx[3]{pvRes, 0.f, pvRes};
1046 float chi2 = temporaryTrack.getPredictedChi2Quiet(posVtx, covVtx);
1047 if (chi2 < bestChi2) {
1048 if (!temporaryTrack.track::TrackParCov::update(posVtx, covVtx)) {
1049 continue;
1050 }
1051 bestTrack = temporaryTrack;
1052 bestChi2 = chi2;
1053 }
1054 }
1055
1056 bestTrack.resetCovariance();
1057 bestTrack.setChi2(0.f);
1058 fitSuccess = fitTrack(bestTrack, 0, mTrkParams[0].NLayers, 1, mTrkParams[0].MaxChi2ClusterAttachment, mTrkParams[0].MaxChi2NDF);
1059 if (!fitSuccess) {
1060 continue;
1061 }
1062 bestTrack.getParamOut() = bestTrack;
1063 bestTrack.resetCovariance();
1064 bestTrack.setChi2(0.f);
1065 fitSuccess = fitTrack(bestTrack, mTrkParams[0].NLayers - 1, -1, -1, mTrkParams[0].MaxChi2ClusterAttachment, mTrkParams[0].MaxChi2NDF, 50.);
1066 if (!fitSuccess) {
1067 continue;
1068 }
1069 mTimeFrame->markUsedCluster(0, bestTrack.getClusterIndex(0));
1070 mTimeFrame->markUsedCluster(1, bestTrack.getClusterIndex(1));
1071 mTimeFrame->markUsedCluster(2, bestTrack.getClusterIndex(2));
1072 mTimeFrame->getTracks(rof).emplace_back(bestTrack);
1073 }
1074}
1075
1076template <int nLayers>
1077bool TrackerTraits<nLayers>::fitTrack(TrackITSExt& track, int start, int end, int step, float chi2clcut, float chi2ndfcut, float maxQoverPt, int nCl)
1078{
1079 auto propInstance = o2::base::Propagator::Instance();
1080
1081 for (int iLayer{start}; iLayer != end; iLayer += step) {
1082 if (track.getClusterIndex(iLayer) == constants::UnusedIndex) {
1083 continue;
1084 }
1085 const TrackingFrameInfo& trackingHit = mTimeFrame->getTrackingFrameInfoOnLayer(iLayer)[track.getClusterIndex(iLayer)];
1086
1087 if (!track.rotate(trackingHit.alphaTrackingFrame)) {
1088 return false;
1089 }
1090
1091 if (!propInstance->propagateToX(track, trackingHit.xTrackingFrame, getBz(), o2::base::PropagatorImpl<float>::MAX_SIN_PHI, o2::base::PropagatorImpl<float>::MAX_STEP, mTrkParams[0].CorrType)) {
1092 return false;
1093 }
1094
1095 if (mTrkParams[0].CorrType == o2::base::PropagatorF::MatCorrType::USEMatCorrNONE) {
1096 if (!track.correctForMaterial(mTrkParams[0].LayerxX0[iLayer], mTrkParams[0].LayerxX0[iLayer] * constants::Radl * constants::Rho, true)) {
1097 continue;
1098 }
1099 }
1100
1101 auto predChi2{track.getPredictedChi2Quiet(trackingHit.positionTrackingFrame, trackingHit.covarianceTrackingFrame)};
1102 if ((nCl >= 3 && predChi2 > chi2clcut) || predChi2 < 0.f) {
1103 return false;
1104 }
1105 track.setChi2(track.getChi2() + predChi2);
1106 if (!track.o2::track::TrackParCov::update(trackingHit.positionTrackingFrame, trackingHit.covarianceTrackingFrame)) {
1107 return false;
1108 }
1109 nCl++;
1110 }
1111 return std::abs(track.getQ2Pt()) < maxQoverPt && track.getChi2() < chi2ndfcut * (nCl * 2 - 5);
1112}
1113
1114template <int nLayers>
1115bool TrackerTraits<nLayers>::trackFollowing(TrackITSExt* track, int rof, bool outward, const int iteration)
1116{
1117 auto propInstance = o2::base::Propagator::Instance();
1118 const int step = -1 + outward * 2;
1119 const int end = outward ? mTrkParams[iteration].NLayers - 1 : 0;
1120 bounded_vector<TrackITSExt> hypotheses(1, *track, mMemoryPool.get()); // possibly avoid reallocation
1121 for (size_t iHypo{0}; iHypo < hypotheses.size(); ++iHypo) {
1122 auto hypo{hypotheses[iHypo]};
1123 int iLayer = static_cast<int>(outward ? hypo.getLastClusterLayer() : hypo.getFirstClusterLayer());
1124 // per layer we add new hypotheses
1125 while (iLayer != end) {
1126 iLayer += step; // step through all layers until we reach the end, this allows for skipping on empty layers
1127 const float r = mTrkParams[iteration].LayerRadii[iLayer];
1128 // get an estimate of the trackinf-frame x for the next step
1129 float x{-999};
1130 if (!hypo.getXatLabR(r, x, mTimeFrame->getBz(), o2::track::DirAuto) || x <= 0.f) {
1131 continue;
1132 }
1133 // estimate hypo's trk parameters at that x
1134 auto& hypoParam{outward ? hypo.getParamOut() : hypo.getParamIn()};
1135 if (!propInstance->propagateToX(hypoParam, x, mTimeFrame->getBz(), PropagatorF::MAX_SIN_PHI,
1136 PropagatorF::MAX_STEP, mTrkParams[iteration].CorrType)) {
1137 continue;
1138 }
1139
1140 if (mTrkParams[iteration].CorrType == PropagatorF::MatCorrType::USEMatCorrNONE) { // account for material affects if propagator does not
1141 if (!hypoParam.correctForMaterial(mTrkParams[iteration].LayerxX0[iLayer], mTrkParams[iteration].LayerxX0[iLayer] * constants::Radl * constants::Rho, true)) {
1142 continue;
1143 }
1144 }
1145
1146 // calculate the search window on this layer
1147 const float phi{hypoParam.getPhi()};
1148 const float ePhi{o2::gpu::CAMath::Sqrt(hypoParam.getSigmaSnp2() / hypoParam.getCsp2())};
1149 const float z{hypoParam.getZ()};
1150 const float eZ{o2::gpu::CAMath::Sqrt(hypoParam.getSigmaZ2())};
1151 const int4 selectedBinsRect{getBinsRect(iLayer, phi, mTrkParams[iteration].TrackFollowerNSigmaCutPhi * ePhi, z, mTrkParams[iteration].TrackFollowerNSigmaCutZ * eZ)};
1152 if (selectedBinsRect.x == 0 && selectedBinsRect.y == 0 && selectedBinsRect.z == 0 && selectedBinsRect.w == 0) {
1153 continue;
1154 }
1155
1156 int phiBinsNum{selectedBinsRect.w - selectedBinsRect.y + 1};
1157
1158 if (phiBinsNum < 0) {
1159 phiBinsNum += mTrkParams[iteration].PhiBins;
1160 }
1161
1162 gsl::span<const Cluster> layer1 = mTimeFrame->getClustersOnLayer(rof, iLayer);
1163 if (layer1.empty()) {
1164 continue;
1165 }
1166
1167 // check all clusters in search windows for possible new hypotheses
1168 for (int iPhiCount = 0; iPhiCount < phiBinsNum; iPhiCount++) {
1169 int iPhiBin = (selectedBinsRect.y + iPhiCount) % mTrkParams[iteration].PhiBins;
1170 const int firstBinIndex{mTimeFrame->mIndexTableUtils.getBinIndex(selectedBinsRect.x, iPhiBin)};
1171 const int maxBinIndex{firstBinIndex + selectedBinsRect.z - selectedBinsRect.x + 1};
1172 const int firstRowClusterIndex = mTimeFrame->getIndexTable(rof, iLayer)[firstBinIndex];
1173 const int maxRowClusterIndex = mTimeFrame->getIndexTable(rof, iLayer)[maxBinIndex];
1174
1175 for (int iNextCluster{firstRowClusterIndex}; iNextCluster < maxRowClusterIndex; ++iNextCluster) {
1176 if (iNextCluster >= (int)layer1.size()) {
1177 break;
1178 }
1179 const Cluster& nextCluster{layer1[iNextCluster]};
1180
1181 if (mTimeFrame->isClusterUsed(iLayer, nextCluster.clusterId)) {
1182 continue;
1183 }
1184
1185 const TrackingFrameInfo& trackingHit = mTimeFrame->getTrackingFrameInfoOnLayer(iLayer)[nextCluster.clusterId];
1186
1187 auto tbupdated{hypo};
1188 auto& tbuParams = outward ? tbupdated.getParamOut() : tbupdated.getParamIn();
1189 if (!tbuParams.rotate(trackingHit.alphaTrackingFrame)) {
1190 continue;
1191 }
1192
1193 if (!propInstance->propagateToX(tbuParams, trackingHit.xTrackingFrame, mTimeFrame->getBz(),
1194 PropagatorF::MAX_SIN_PHI, PropagatorF::MAX_STEP, PropagatorF::MatCorrType::USEMatCorrNONE)) {
1195 continue;
1196 }
1197
1198 auto predChi2{tbuParams.getPredictedChi2Quiet(trackingHit.positionTrackingFrame, trackingHit.covarianceTrackingFrame)};
1199 if (predChi2 >= track->getChi2() * mTrkParams[iteration].NSigmaCut) {
1200 continue;
1201 }
1202
1203 if (!tbuParams.o2::track::TrackParCov::update(trackingHit.positionTrackingFrame, trackingHit.covarianceTrackingFrame)) {
1204 continue;
1205 }
1206 tbupdated.setChi2(tbupdated.getChi2() + predChi2);
1207 tbupdated.setExternalClusterIndex(iLayer, nextCluster.clusterId, true);
1208 hypotheses.emplace_back(tbupdated);
1209 }
1210 }
1211 }
1212 }
1213
1214 TrackITSExt* bestHypo{track};
1215 bool swapped{false};
1216 for (auto& hypo : hypotheses) {
1217 if (hypo.isBetter(*bestHypo, track->getChi2() * mTrkParams[iteration].NSigmaCut)) {
1218 bestHypo = &hypo;
1219 swapped = true;
1220 }
1221 }
1222 *track = *bestHypo;
1223 return swapped;
1224}
1225
1228template <int nLayers>
1230{
1231 float ca{-999.f}, sa{-999.f};
1232 o2::gpu::CAMath::SinCos(tf3.alphaTrackingFrame, sa, ca);
1233 const float x1 = cluster1.xCoordinate * ca + cluster1.yCoordinate * sa;
1234 const float y1 = -cluster1.xCoordinate * sa + cluster1.yCoordinate * ca;
1235 const float z1 = cluster1.zCoordinate;
1236 const float x2 = cluster2.xCoordinate * ca + cluster2.yCoordinate * sa;
1237 const float y2 = -cluster2.xCoordinate * sa + cluster2.yCoordinate * ca;
1238 const float z2 = cluster2.zCoordinate;
1239 const float x3 = tf3.xTrackingFrame;
1240 const float y3 = tf3.positionTrackingFrame[0];
1241 const float z3 = tf3.positionTrackingFrame[1];
1242 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};
1243 if (mIsZeroField) {
1244 tgp = o2::gpu::CAMath::ATan2(y3 - y1, x3 - x1);
1245 snp = tgp / o2::gpu::CAMath::Sqrt(1.f + tgp * tgp);
1246 } else {
1247 crv = math_utils::computeCurvature(x3, y3, x2, y2, x1, y1);
1248 snp = crv * (x3 - math_utils::computeCurvatureCentreX(x3, y3, x2, y2, x1, y1));
1249 q2pt = crv / (mBz * o2::constants::math::B2C);
1250 q2pt2 = crv * crv;
1251 }
1252 tgl12 = math_utils::computeTanDipAngle(x1, y1, x2, y2, z1, z2);
1253 tgl23 = math_utils::computeTanDipAngle(x2, y2, x3, y3, z2, z3);
1254 sg2q2pt = track::kC1Pt2max * (q2pt2 > 0.0005f ? (q2pt2 < 1.f ? q2pt2 : 1.f) : 0.0005f);
1255 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}};
1256}
1257
1258template <int nLayers>
1260{
1261 mBz = bz;
1262 mIsZeroField = std::abs(mBz) < 0.01;
1263 mTimeFrame->setBz(bz);
1264}
1265
1266template <int nLayers>
1268{
1269 return o2::base::Propagator::Instance()->getMatLUT() && (mTrkParams[0].CorrType == o2::base::PropagatorImpl<float>::MatCorrType::USEMatCorrLUT);
1270}
1271
1272template <int nLayers>
1273void TrackerTraits<nLayers>::setNThreads(int n, std::shared_ptr<tbb::task_arena>& arena)
1274{
1275#if defined(OPTIMISATION_OUTPUT) || defined(CA_DEBUG)
1276 mTaskArena = std::make_shared<tbb::task_arena>(1);
1277#else
1278 if (arena == nullptr) {
1279 mTaskArena = std::make_shared<tbb::task_arena>(std::abs(n));
1280 LOGP(info, "Setting tracker with {} threads.", n);
1281 } else {
1282 mTaskArena = arena;
1283 LOGP(info, "Attaching tracker to calling thread's arena");
1284 }
1285#endif
1286}
1287
1288template class TrackerTraits<7>;
1289
1290} // namespace o2::its
Base track model for the Barrel, params only, w/o covariance.
#define CA_DEBUGGER(x)
Definition Definitions.h:25
const auto bins
Definition PID.cxx:49
void print() const
#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 findRoads(const int iteration)
void setNThreads(int n, std::shared_ptr< tbb::task_arena > &arena)
virtual void processNeighbours(int iLayer, int iLevel, const bounded_vector< CellSeed > &currentCellSeed, const bounded_vector< int > &currentCellId, bounded_vector< CellSeed > &updatedCellSeed, bounded_vector< int > &updatedCellId)
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
constexpr int UnusedIndex
Definition Constants.h:30
constexpr float Radl
Definition Constants.h:32
constexpr float Tolerance
Definition Constants.h:28
constexpr float Rho
Definition Constants.h:33
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:55
float zCoordinate
Definition Cluster.h:56
float xCoordinate
Definition Cluster.h:54
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:82
std::array< float, 3 > covarianceTrackingFrame
Definition Cluster.h:83
std::array< uint16_t, 5 > pattern
std::vector< Tracklet64 > tracklets