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