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