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