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
17
18#include <algorithm>
19#include <cassert>
20#include <iostream>
21
22#ifdef OPTIMISATION_OUTPUT
23#include <format>
24#endif
25
28#include "GPUCommonMath.h"
29#include "ITStracking/Cell.h"
34
35#ifdef WITH_OPENMP
36#include <omp.h>
37#endif
38
40
41namespace
42{
43inline float Sq(float q)
44{
45 return q * q;
46}
47} // namespace
48
49namespace o2
50{
51namespace its
52{
53
54constexpr int debugLevel{0};
55
56void TrackerTraits::computeLayerTracklets(const int iteration, int iROFslice, int iVertex)
57{
59
60#ifdef OPTIMISATION_OUTPUT
61 static int iter{0};
62 std::ofstream off(std::format("tracklets{}.txt", iter++));
63#endif
64
65 for (int iLayer = 0; iLayer < mTrkParams[iteration].TrackletsPerRoad(); ++iLayer) {
66 tf->getTracklets()[iLayer].clear();
67 tf->getTrackletsLabel(iLayer).clear();
68 if (iLayer > 0) {
69 std::fill(tf->getTrackletsLookupTable()[iLayer - 1].begin(), tf->getTrackletsLookupTable()[iLayer - 1].end(), 0);
70 }
71 }
72
73 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);
74 gsl::span<const Vertex> diamondSpan(&diamondVert, 1);
75 int startROF{mTrkParams[iteration].nROFsPerIterations > 0 ? iROFslice * mTrkParams[iteration].nROFsPerIterations : 0};
76 int endROF{gpu::GPUCommonMath::Min(mTrkParams[iteration].nROFsPerIterations > 0 ? (iROFslice + 1) * mTrkParams[iteration].nROFsPerIterations + mTrkParams[iteration].DeltaROF : tf->getNrof(), tf->getNrof())};
77 for (int rof0{startROF}; rof0 < endROF; ++rof0) {
78 gsl::span<const Vertex> primaryVertices = mTrkParams[iteration].UseDiamond ? diamondSpan : tf->getPrimaryVertices(rof0);
79 const int startVtx{iVertex >= 0 ? iVertex : 0};
80 const int endVtx{iVertex >= 0 ? o2::gpu::CAMath::Min(iVertex + 1, static_cast<int>(primaryVertices.size())) : static_cast<int>(primaryVertices.size())};
81 int minRof = o2::gpu::CAMath::Max(startROF, rof0 - mTrkParams[iteration].DeltaROF);
82 int maxRof = o2::gpu::CAMath::Min(endROF - 1, rof0 + mTrkParams[iteration].DeltaROF);
83#pragma omp parallel for num_threads(mNThreads)
84 for (int iLayer = 0; iLayer < mTrkParams[iteration].TrackletsPerRoad(); ++iLayer) {
85 gsl::span<const Cluster> layer0 = tf->getClustersOnLayer(rof0, iLayer);
86 if (layer0.empty()) {
87 continue;
88 }
89 float meanDeltaR{mTrkParams[iteration].LayerRadii[iLayer + 1] - mTrkParams[iteration].LayerRadii[iLayer]};
90
91 const int currentLayerClustersNum{static_cast<int>(layer0.size())};
92 for (int iCluster{0}; iCluster < currentLayerClustersNum; ++iCluster) {
93 const Cluster& currentCluster{layer0[iCluster]};
94 const int currentSortedIndex{tf->getSortedIndex(rof0, iLayer, iCluster)};
95
96 if (tf->isClusterUsed(iLayer, currentCluster.clusterId)) {
97 continue;
98 }
99 const float inverseR0{1.f / currentCluster.radius};
100
101 for (int iV{startVtx}; iV < endVtx; ++iV) {
102 auto& primaryVertex{primaryVertices[iV]};
103 if (primaryVertex.isFlagSet(2) && iteration != 3) {
104 continue;
105 }
106 const float resolution = o2::gpu::CAMath::Sqrt(Sq(mTrkParams[iteration].PVres) / primaryVertex.getNContributors() + Sq(tf->getPositionResolution(iLayer)));
107
108 const float tanLambda{(currentCluster.zCoordinate - primaryVertex.getZ()) * inverseR0};
109
110 const float zAtRmin{tanLambda * (tf->getMinR(iLayer + 1) - currentCluster.radius) + currentCluster.zCoordinate};
111 const float zAtRmax{tanLambda * (tf->getMaxR(iLayer + 1) - currentCluster.radius) + currentCluster.zCoordinate};
112
113 const float sqInverseDeltaZ0{1.f / (Sq(currentCluster.zCoordinate - primaryVertex.getZ()) + 2.e-8f)};
114 const float sigmaZ{o2::gpu::CAMath::Sqrt(Sq(resolution) * Sq(tanLambda) * ((Sq(inverseR0) + sqInverseDeltaZ0) * Sq(meanDeltaR) + 1.f) + Sq(meanDeltaR * tf->getMSangle(iLayer)))};
115
116 const int4 selectedBinsRect{getBinsRect(currentCluster, iLayer + 1, zAtRmin, zAtRmax,
117 sigmaZ * mTrkParams[iteration].NSigmaCut, tf->getPhiCut(iLayer))};
118 if (selectedBinsRect.x == 0 && selectedBinsRect.y == 0 && selectedBinsRect.z == 0 && selectedBinsRect.w == 0) {
119 continue;
120 }
121
122 int phiBinsNum{selectedBinsRect.w - selectedBinsRect.y + 1};
123
124 if (phiBinsNum < 0) {
125 phiBinsNum += mTrkParams[iteration].PhiBins;
126 }
127
128 for (int rof1{minRof}; rof1 <= maxRof; ++rof1) {
129 gsl::span<const Cluster> layer1 = tf->getClustersOnLayer(rof1, iLayer + 1);
130 if (layer1.empty()) {
131 continue;
132 }
133 for (int iPhiCount{0}; iPhiCount < phiBinsNum; iPhiCount++) {
134 int iPhiBin = (selectedBinsRect.y + iPhiCount) % mTrkParams[iteration].PhiBins;
135 const int firstBinIndex{tf->mIndexTableUtils.getBinIndex(selectedBinsRect.x, iPhiBin)};
136 const int maxBinIndex{firstBinIndex + selectedBinsRect.z - selectedBinsRect.x + 1};
137 if constexpr (debugLevel) {
138 if (firstBinIndex < 0 || firstBinIndex > tf->getIndexTable(rof1, iLayer + 1).size() ||
139 maxBinIndex < 0 || maxBinIndex > tf->getIndexTable(rof1, iLayer + 1).size()) {
140 std::cout << iLayer << "\t" << iCluster << "\t" << zAtRmin << "\t" << zAtRmax << "\t" << sigmaZ * mTrkParams[iteration].NSigmaCut << "\t" << tf->getPhiCut(iLayer) << std::endl;
141 std::cout << currentCluster.zCoordinate << "\t" << primaryVertex.getZ() << "\t" << currentCluster.radius << std::endl;
142 std::cout << tf->getMinR(iLayer + 1) << "\t" << currentCluster.radius << "\t" << currentCluster.zCoordinate << std::endl;
143 std::cout << "Illegal access to IndexTable " << firstBinIndex << "\t" << maxBinIndex << "\t" << selectedBinsRect.z << "\t" << selectedBinsRect.x << std::endl;
144 exit(1);
145 }
146 }
147 const int firstRowClusterIndex = tf->getIndexTable(rof1, iLayer + 1)[firstBinIndex];
148 const int maxRowClusterIndex = tf->getIndexTable(rof1, iLayer + 1)[maxBinIndex];
149 for (int iNextCluster{firstRowClusterIndex}; iNextCluster < maxRowClusterIndex; ++iNextCluster) {
150 if (iNextCluster >= (int)layer1.size()) {
151 break;
152 }
153
154 const Cluster& nextCluster{layer1[iNextCluster]};
155 if (tf->isClusterUsed(iLayer + 1, nextCluster.clusterId)) {
156 continue;
157 }
158
159 const float deltaPhi{gpu::GPUCommonMath::Abs(currentCluster.phi - nextCluster.phi)};
160 const float deltaZ{gpu::GPUCommonMath::Abs(tanLambda * (nextCluster.radius - currentCluster.radius) +
161 currentCluster.zCoordinate - nextCluster.zCoordinate)};
162
163#ifdef OPTIMISATION_OUTPUT
165 int currentId{currentCluster.clusterId};
166 int nextId{nextCluster.clusterId};
167 for (auto& lab1 : tf->getClusterLabels(iLayer, currentId)) {
168 for (auto& lab2 : tf->getClusterLabels(iLayer + 1, nextId)) {
169 if (lab1 == lab2 && lab1.isValid()) {
170 label = lab1;
171 break;
172 }
173 }
174 if (label.isValid()) {
175 break;
176 }
177 }
178 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;
179#endif
180
181 if (deltaZ / sigmaZ < mTrkParams[iteration].NSigmaCut &&
182 (deltaPhi < tf->getPhiCut(iLayer) ||
183 gpu::GPUCommonMath::Abs(deltaPhi - constants::math::TwoPi) < tf->getPhiCut(iLayer))) {
184 if (iLayer > 0) {
185 tf->getTrackletsLookupTable()[iLayer - 1][currentSortedIndex]++;
186 }
187 const float phi{o2::gpu::GPUCommonMath::ATan2(currentCluster.yCoordinate - nextCluster.yCoordinate,
188 currentCluster.xCoordinate - nextCluster.xCoordinate)};
189 const float tanL{(currentCluster.zCoordinate - nextCluster.zCoordinate) /
190 (currentCluster.radius - nextCluster.radius)};
191 tf->getTracklets()[iLayer].emplace_back(currentSortedIndex, tf->getSortedIndex(rof1, iLayer + 1, iNextCluster), tanL, phi, rof0, rof1);
192 }
193 }
194 }
195 }
196 }
197 }
198 }
199 }
200 if (!tf->checkMemory(mTrkParams[iteration].MaxMemory)) {
201 return;
202 }
203
204#pragma omp parallel for num_threads(mNThreads)
205 for (int iLayer = 0; iLayer < mTrkParams[iteration].CellsPerRoad(); ++iLayer) {
207 auto& trkl{tf->getTracklets()[iLayer + 1]};
208 std::sort(trkl.begin(), trkl.end(), [](const Tracklet& a, const Tracklet& b) {
209 return a.firstClusterIndex < b.firstClusterIndex || (a.firstClusterIndex == b.firstClusterIndex && a.secondClusterIndex < b.secondClusterIndex);
210 });
212 auto& lut{tf->getTrackletsLookupTable()[iLayer]};
213 int id0{-1}, id1{-1};
214 std::vector<Tracklet> newTrk;
215 newTrk.reserve(trkl.size());
216 for (auto& trk : trkl) {
217 if (trk.firstClusterIndex == id0 && trk.secondClusterIndex == id1) {
218 lut[id0]--;
219 } else {
220 id0 = trk.firstClusterIndex;
221 id1 = trk.secondClusterIndex;
222 newTrk.push_back(trk);
223 }
224 }
225 trkl.swap(newTrk);
226
228 std::exclusive_scan(lut.begin(), lut.end(), lut.begin(), 0);
229 lut.push_back(trkl.size());
230 }
232 std::sort(tf->getTracklets()[0].begin(), tf->getTracklets()[0].end(), [](const Tracklet& a, const Tracklet& b) {
233 return a.firstClusterIndex < b.firstClusterIndex || (a.firstClusterIndex == b.firstClusterIndex && a.secondClusterIndex < b.secondClusterIndex);
234 });
235 int id0{-1}, id1{-1};
236 std::vector<Tracklet> newTrk;
237 newTrk.reserve(tf->getTracklets()[0].size());
238 for (auto& trk : tf->getTracklets()[0]) {
239 if (trk.firstClusterIndex != id0 || trk.secondClusterIndex != id1) {
240 id0 = trk.firstClusterIndex;
241 id1 = trk.secondClusterIndex;
242 newTrk.push_back(trk);
243 }
244 }
245 tf->getTracklets()[0].swap(newTrk);
246
248 if (tf->hasMCinformation()) {
249 for (int iLayer{0}; iLayer < mTrkParams[iteration].TrackletsPerRoad(); ++iLayer) {
250 for (auto& trk : tf->getTracklets()[iLayer]) {
252 int currentId{tf->getClusters()[iLayer][trk.firstClusterIndex].clusterId};
253 int nextId{tf->getClusters()[iLayer + 1][trk.secondClusterIndex].clusterId};
254 for (auto& lab1 : tf->getClusterLabels(iLayer, currentId)) {
255 for (auto& lab2 : tf->getClusterLabels(iLayer + 1, nextId)) {
256 if (lab1 == lab2 && lab1.isValid()) {
257 label = lab1;
258 break;
259 }
260 }
261 if (label.isValid()) {
262 break;
263 }
264 }
265 tf->getTrackletsLabel(iLayer).emplace_back(label);
266 }
267 }
268 }
269}
270
271void TrackerTraits::computeLayerCells(const int iteration)
272{
273#ifdef OPTIMISATION_OUTPUT
274 static int iter{0};
275 std::ofstream off(std::format("cells{}.txt", iter++));
276#endif
277
278 for (int iLayer = 0; iLayer < mTrkParams[iteration].CellsPerRoad(); ++iLayer) {
279 mTimeFrame->getCells()[iLayer].clear();
280 mTimeFrame->getCellsLabel(iLayer).clear();
281 if (iLayer > 0) {
282 mTimeFrame->getCellsLookupTable()[iLayer - 1].clear();
283 }
284 }
285
287#pragma omp parallel for num_threads(mNThreads)
288 for (int iLayer = 0; iLayer < mTrkParams[iteration].CellsPerRoad(); ++iLayer) {
289
290 if (tf->getTracklets()[iLayer + 1].empty() ||
291 tf->getTracklets()[iLayer].empty()) {
292 continue;
293 }
294
295#ifdef OPTIMISATION_OUTPUT
296 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]};
297 resolution = resolution > 1.e-12 ? resolution : 1.f;
298#endif
299 const int currentLayerTrackletsNum{static_cast<int>(tf->getTracklets()[iLayer].size())};
300 for (int iTracklet{0}; iTracklet < currentLayerTrackletsNum; ++iTracklet) {
301
302 const Tracklet& currentTracklet{tf->getTracklets()[iLayer][iTracklet]};
303 const int nextLayerClusterIndex{currentTracklet.secondClusterIndex};
304 const int nextLayerFirstTrackletIndex{
305 tf->getTrackletsLookupTable()[iLayer][nextLayerClusterIndex]};
306 const int nextLayerLastTrackletIndex{
307 tf->getTrackletsLookupTable()[iLayer][nextLayerClusterIndex + 1]};
308
309 if (nextLayerFirstTrackletIndex == nextLayerLastTrackletIndex) {
310 continue;
311 }
312
313 for (int iNextTracklet{nextLayerFirstTrackletIndex}; iNextTracklet < nextLayerLastTrackletIndex; ++iNextTracklet) {
314 if (tf->getTracklets()[iLayer + 1][iNextTracklet].firstClusterIndex != nextLayerClusterIndex) {
315 break;
316 }
317 const Tracklet& nextTracklet{tf->getTracklets()[iLayer + 1][iNextTracklet]};
318 const float deltaTanLambda{std::abs(currentTracklet.tanLambda - nextTracklet.tanLambda)};
319
320#ifdef OPTIMISATION_OUTPUT
321 bool good{tf->getTrackletsLabel(iLayer)[iTracklet] == tf->getTrackletsLabel(iLayer + 1)[iNextTracklet]};
322 float signedDelta{currentTracklet.tanLambda - nextTracklet.tanLambda};
323 off << std::format("{}\t{:d}\t{}\t{}\t{}\t{}", iLayer, good, signedDelta, signedDelta / (mTrkParams[iteration].CellDeltaTanLambdaSigma), tanLambda, resolution) << std::endl;
324#endif
325
326 if (deltaTanLambda / mTrkParams[iteration].CellDeltaTanLambdaSigma < mTrkParams[iteration].NSigmaCut) {
327
329 const int clusId[3]{
330 mTimeFrame->getClusters()[iLayer][currentTracklet.firstClusterIndex].clusterId,
331 mTimeFrame->getClusters()[iLayer + 1][nextTracklet.firstClusterIndex].clusterId,
332 mTimeFrame->getClusters()[iLayer + 2][nextTracklet.secondClusterIndex].clusterId};
333 const auto& cluster1_glo = mTimeFrame->getUnsortedClusters()[iLayer].at(clusId[0]);
334 const auto& cluster2_glo = mTimeFrame->getUnsortedClusters()[iLayer + 1].at(clusId[1]);
335 const auto& cluster3_tf = mTimeFrame->getTrackingFrameInfoOnLayer(iLayer + 2).at(clusId[2]);
336 auto track{buildTrackSeed(cluster1_glo, cluster2_glo, cluster3_tf)};
337
338 float chi2{0.f};
339 bool good{false};
340 for (int iC{2}; iC--;) {
341 const TrackingFrameInfo& trackingHit = mTimeFrame->getTrackingFrameInfoOnLayer(iLayer + iC).at(clusId[iC]);
342
343 if (!track.rotate(trackingHit.alphaTrackingFrame)) {
344 break;
345 }
346
347 if (!track.propagateTo(trackingHit.xTrackingFrame, getBz())) {
348 break;
349 }
350
351 constexpr float radl = 9.36f; // Radiation length of Si [cm]
352 constexpr float rho = 2.33f; // Density of Si [g/cm^3]
353 if (!track.correctForMaterial(mTrkParams[0].LayerxX0[iLayer + iC], mTrkParams[0].LayerxX0[iLayer] * radl * rho, true)) {
354 break;
355 }
356
357 auto predChi2{track.getPredictedChi2Quiet(trackingHit.positionTrackingFrame, trackingHit.covarianceTrackingFrame)};
358 if (!track.o2::track::TrackParCov::update(trackingHit.positionTrackingFrame, trackingHit.covarianceTrackingFrame)) {
359 break;
360 }
361 if (!iC && predChi2 > mTrkParams[iteration].MaxChi2ClusterAttachment) {
362 break;
363 }
364 good = !iC;
365 chi2 += predChi2;
366 }
367 if (!good) {
368 continue;
369 }
370 if (iLayer > 0 && (int)tf->getCellsLookupTable()[iLayer - 1].size() <= iTracklet) {
371 tf->getCellsLookupTable()[iLayer - 1].resize(iTracklet + 1, tf->getCells()[iLayer].size());
372 }
373 tf->getCells()[iLayer].emplace_back(iLayer, clusId[0], clusId[1], clusId[2],
374 iTracklet, iNextTracklet, track, chi2);
375 }
376 }
377 }
378 if (iLayer > 0) {
379 tf->getCellsLookupTable()[iLayer - 1].resize(currentLayerTrackletsNum + 1, tf->getCells()[iLayer].size());
380 }
381 }
382 if (!tf->checkMemory(mTrkParams[iteration].MaxMemory)) {
383 return;
384 }
385
387 if (tf->hasMCinformation()) {
388 for (int iLayer{0}; iLayer < mTrkParams[iteration].CellsPerRoad(); ++iLayer) {
389 for (auto& cell : tf->getCells()[iLayer]) {
390 MCCompLabel currentLab{tf->getTrackletsLabel(iLayer)[cell.getFirstTrackletIndex()]};
391 MCCompLabel nextLab{tf->getTrackletsLabel(iLayer + 1)[cell.getSecondTrackletIndex()]};
392 tf->getCellsLabel(iLayer).emplace_back(currentLab == nextLab ? currentLab : MCCompLabel());
393 }
394 }
395 }
396
397 if constexpr (debugLevel) {
398 for (int iLayer{0}; iLayer < mTrkParams[iteration].CellsPerRoad(); ++iLayer) {
399 std::cout << "Cells on layer " << iLayer << " " << tf->getCells()[iLayer].size() << std::endl;
400 }
401 }
402}
403
404void TrackerTraits::findCellsNeighbours(const int iteration)
405{
406#ifdef OPTIMISATION_OUTPUT
407 std::ofstream off(std::format("cellneighs{}.txt", iteration));
408#endif
409 for (int iLayer{0}; iLayer < mTrkParams[iteration].CellsPerRoad() - 1; ++iLayer) {
410 const int nextLayerCellsNum{static_cast<int>(mTimeFrame->getCells()[iLayer + 1].size())};
411 mTimeFrame->getCellsNeighboursLUT()[iLayer].clear();
412 mTimeFrame->getCellsNeighboursLUT()[iLayer].resize(nextLayerCellsNum, 0);
413 if (mTimeFrame->getCells()[iLayer + 1].empty() ||
414 mTimeFrame->getCellsLookupTable()[iLayer].empty()) {
415 mTimeFrame->getCellsNeighbours()[iLayer].clear();
416 continue;
417 }
418
419 int layerCellsNum{static_cast<int>(mTimeFrame->getCells()[iLayer].size())};
420 std::vector<std::pair<int, int>> cellsNeighbours;
421 cellsNeighbours.reserve(nextLayerCellsNum);
422
423 for (int iCell{0}; iCell < layerCellsNum; ++iCell) {
424
425 const auto& currentCellSeed{mTimeFrame->getCells()[iLayer][iCell]};
426 const int nextLayerTrackletIndex{currentCellSeed.getSecondTrackletIndex()};
427 const int nextLayerFirstCellIndex{mTimeFrame->getCellsLookupTable()[iLayer][nextLayerTrackletIndex]};
428 const int nextLayerLastCellIndex{mTimeFrame->getCellsLookupTable()[iLayer][nextLayerTrackletIndex + 1]};
429 for (int iNextCell{nextLayerFirstCellIndex}; iNextCell < nextLayerLastCellIndex; ++iNextCell) {
430
431 auto nextCellSeed{mTimeFrame->getCells()[iLayer + 1][iNextCell]};
432 if (nextCellSeed.getFirstTrackletIndex() != nextLayerTrackletIndex) {
433 break;
434 }
435
436 if (!nextCellSeed.rotate(currentCellSeed.getAlpha()) ||
437 !nextCellSeed.propagateTo(currentCellSeed.getX(), getBz())) {
438 continue;
439 }
440 float chi2 = currentCellSeed.getPredictedChi2(nextCellSeed);
441
442#ifdef OPTIMISATION_OUTPUT
443 bool good{mTimeFrame->getCellsLabel(iLayer)[iCell] == mTimeFrame->getCellsLabel(iLayer + 1)[iNextCell]};
444 off << std::format("{}\t{:d}\t{}", iLayer, good, chi2) << std::endl;
445#endif
446
447 if (chi2 > mTrkParams[0].MaxChi2ClusterAttachment) {
448 continue;
449 }
450
451 mTimeFrame->getCellsNeighboursLUT()[iLayer][iNextCell]++;
452 cellsNeighbours.push_back(std::make_pair(iCell, iNextCell));
453 const int currentCellLevel{currentCellSeed.getLevel()};
454
455 if (currentCellLevel >= nextCellSeed.getLevel()) {
456 mTimeFrame->getCells()[iLayer + 1][iNextCell].setLevel(currentCellLevel + 1);
457 }
458 }
459 }
460 std::sort(cellsNeighbours.begin(), cellsNeighbours.end(), [](const std::pair<int, int>& a, const std::pair<int, int>& b) {
461 return a.second < b.second;
462 });
463 mTimeFrame->getCellsNeighbours()[iLayer].clear();
464 mTimeFrame->getCellsNeighbours()[iLayer].reserve(cellsNeighbours.size());
465 for (auto& cellNeighboursIndex : cellsNeighbours) {
466 mTimeFrame->getCellsNeighbours()[iLayer].push_back(cellNeighboursIndex.first);
467 }
468 std::inclusive_scan(mTimeFrame->getCellsNeighboursLUT()[iLayer].begin(), mTimeFrame->getCellsNeighboursLUT()[iLayer].end(), mTimeFrame->getCellsNeighboursLUT()[iLayer].begin());
469 }
470}
471
472void TrackerTraits::processNeighbours(int iLayer, int iLevel, const std::vector<CellSeed>& currentCellSeed, const std::vector<int>& currentCellId, std::vector<CellSeed>& updatedCellSeeds, std::vector<int>& updatedCellsIds)
473{
474 bool print = iLayer == 3 && iLevel == 2;
475 if (iLevel < 2 || iLayer < 1) {
476 std::cout << "Error: layer " << iLayer << " or level " << iLevel << " cannot be processed by processNeighbours" << std::endl;
477 exit(1);
478 }
479 CA_DEBUGGER(std::cout << "Processing neighbours layer " << iLayer << " level " << iLevel << ", size of the cell seeds: " << currentCellSeed.size() << std::endl);
480 updatedCellSeeds.reserve(mTimeFrame->getCellsNeighboursLUT()[iLayer - 1].size());
481 updatedCellsIds.reserve(updatedCellSeeds.size());
482 auto propagator = o2::base::Propagator::Instance();
483#ifdef CA_DEBUG
484 int failed[5]{0, 0, 0, 0, 0}, attempts{0}, failedByMismatch{0};
485#endif
486
487#pragma omp parallel for num_threads(mNThreads)
488 for (unsigned int iCell = 0; iCell < currentCellSeed.size(); ++iCell) {
489 const CellSeed& currentCell{currentCellSeed[iCell]};
490 if (currentCell.getLevel() != iLevel) {
491 continue;
492 }
493 if (currentCellId.empty() && (mTimeFrame->isClusterUsed(iLayer, currentCell.getFirstClusterIndex()) ||
494 mTimeFrame->isClusterUsed(iLayer + 1, currentCell.getSecondClusterIndex()) ||
495 mTimeFrame->isClusterUsed(iLayer + 2, currentCell.getThirdClusterIndex()))) {
496 continue;
497 }
498 const int cellId = currentCellId.empty() ? iCell : currentCellId[iCell];
499 const int startNeighbourId{cellId ? mTimeFrame->getCellsNeighboursLUT()[iLayer - 1][cellId - 1] : 0};
500 const int endNeighbourId{mTimeFrame->getCellsNeighboursLUT()[iLayer - 1][cellId]};
501
502 for (int iNeighbourCell{startNeighbourId}; iNeighbourCell < endNeighbourId; ++iNeighbourCell) {
503 CA_DEBUGGER(attempts++);
504 const int neighbourCellId = mTimeFrame->getCellsNeighbours()[iLayer - 1][iNeighbourCell];
505 const CellSeed& neighbourCell = mTimeFrame->getCells()[iLayer - 1][neighbourCellId];
506 if (neighbourCell.getSecondTrackletIndex() != currentCell.getFirstTrackletIndex()) {
507 CA_DEBUGGER(failedByMismatch++);
508 continue;
509 }
510 if (mTimeFrame->isClusterUsed(iLayer - 1, neighbourCell.getFirstClusterIndex())) {
511 continue;
512 }
513 if (currentCell.getLevel() - 1 != neighbourCell.getLevel()) {
514 CA_DEBUGGER(failed[0]++);
515 continue;
516 }
518 CellSeed seed{currentCell};
519 auto& trHit = mTimeFrame->getTrackingFrameInfoOnLayer(iLayer - 1).at(neighbourCell.getFirstClusterIndex());
520
521 if (!seed.rotate(trHit.alphaTrackingFrame)) {
522 CA_DEBUGGER(failed[1]++);
523 continue;
524 }
525
526 if (!propagator->propagateToX(seed, trHit.xTrackingFrame, getBz(), o2::base::PropagatorImpl<float>::MAX_SIN_PHI, o2::base::PropagatorImpl<float>::MAX_STEP, mCorrType)) {
527 CA_DEBUGGER(failed[2]++);
528 continue;
529 }
530
531 if (mCorrType == o2::base::PropagatorF::MatCorrType::USEMatCorrNONE) {
532 float radl = 9.36f; // Radiation length of Si [cm]
533 float rho = 2.33f; // Density of Si [g/cm^3]
534 if (!seed.correctForMaterial(mTrkParams[0].LayerxX0[iLayer - 1], mTrkParams[0].LayerxX0[iLayer - 1] * radl * rho, true)) {
535 continue;
536 }
537 }
538
539 auto predChi2{seed.getPredictedChi2Quiet(trHit.positionTrackingFrame, trHit.covarianceTrackingFrame)};
540 if ((predChi2 > mTrkParams[0].MaxChi2ClusterAttachment) || predChi2 < 0.f) {
541 CA_DEBUGGER(failed[3]++);
542 continue;
543 }
544 seed.setChi2(seed.getChi2() + predChi2);
545 if (!seed.o2::track::TrackParCov::update(trHit.positionTrackingFrame, trHit.covarianceTrackingFrame)) {
546 CA_DEBUGGER(failed[4]++);
547 continue;
548 }
549 seed.getClusters()[iLayer - 1] = neighbourCell.getFirstClusterIndex();
550 seed.setLevel(neighbourCell.getLevel());
551 seed.setFirstTrackletIndex(neighbourCell.getFirstTrackletIndex());
552 seed.setSecondTrackletIndex(neighbourCell.getSecondTrackletIndex());
553#pragma omp critical
554 {
555 updatedCellsIds.push_back(neighbourCellId);
556 updatedCellSeeds.push_back(seed);
557 }
558 }
559 }
560#ifdef CA_DEBUG
561 std::cout << "\t\t- Found " << updatedCellSeeds.size() << " cell seeds out of " << attempts << " attempts" << std::endl;
562 std::cout << "\t\t\t> " << failed[0] << " failed because of level" << std::endl;
563 std::cout << "\t\t\t> " << failed[1] << " failed because of rotation" << std::endl;
564 std::cout << "\t\t\t> " << failed[2] << " failed because of propagation" << std::endl;
565 std::cout << "\t\t\t> " << failed[3] << " failed because of chi2 cut" << std::endl;
566 std::cout << "\t\t\t> " << failed[4] << " failed because of update" << std::endl;
567 std::cout << "\t\t\t> " << failedByMismatch << " failed because of mismatch" << std::endl;
568#endif
569}
570
571void TrackerTraits::findRoads(const int iteration)
572{
573 CA_DEBUGGER(std::cout << "Finding roads, iteration " << iteration << std::endl);
574 for (int startLevel{mTrkParams[iteration].CellsPerRoad()}; startLevel >= mTrkParams[iteration].CellMinimumLevel(); --startLevel) {
575 CA_DEBUGGER(std::cout << "\t > Processing level " << startLevel << std::endl);
576 const int minimumLayer{startLevel - 1};
577 std::vector<CellSeed> trackSeeds;
578 for (int startLayer{mTrkParams[iteration].CellsPerRoad() - 1}; startLayer >= minimumLayer; --startLayer) {
579 if ((mTrkParams[iteration].StartLayerMask & (1 << (startLayer + 2))) == 0) {
580 continue;
581 }
582 CA_DEBUGGER(std::cout << "\t\t > Starting processing layer " << startLayer << std::endl);
583 std::vector<int> lastCellId, updatedCellId;
584 std::vector<CellSeed> lastCellSeed, updatedCellSeed;
585
586 processNeighbours(startLayer, startLevel, mTimeFrame->getCells()[startLayer], lastCellId, updatedCellSeed, updatedCellId);
587
588 int level = startLevel;
589 for (int iLayer{startLayer - 1}; iLayer > 0 && level > 2; --iLayer) {
590 lastCellSeed.swap(updatedCellSeed);
591 lastCellId.swap(updatedCellId);
592 std::vector<CellSeed>().swap(updatedCellSeed);
593 updatedCellId.clear();
594 processNeighbours(iLayer, --level, lastCellSeed, lastCellId, updatedCellSeed, updatedCellId);
595 }
596 for (auto& seed : updatedCellSeed) {
597 if (seed.getQ2Pt() > 1.e3 || seed.getChi2() > mTrkParams[0].MaxChi2NDF * ((startLevel + 2) * 2 - 5)) {
598 continue;
599 }
600 trackSeeds.push_back(seed);
601 }
602 }
603
604 std::vector<TrackITSExt> tracks(trackSeeds.size());
605 std::atomic<size_t> trackIndex{0};
606#pragma omp parallel for num_threads(mNThreads)
607 for (size_t seedId = 0; seedId < trackSeeds.size(); ++seedId) {
608 const CellSeed& seed{trackSeeds[seedId]};
609 TrackITSExt temporaryTrack{seed};
610 temporaryTrack.resetCovariance();
611 temporaryTrack.setChi2(0);
612 for (int iL{0}; iL < 7; ++iL) {
613 temporaryTrack.setExternalClusterIndex(iL, seed.getCluster(iL), seed.getCluster(iL) != constants::its::UnusedIndex);
614 }
615
616 bool fitSuccess = fitTrack(temporaryTrack, 0, mTrkParams[0].NLayers, 1, mTrkParams[0].MaxChi2ClusterAttachment, mTrkParams[0].MaxChi2NDF);
617 if (!fitSuccess) {
618 continue;
619 }
620 temporaryTrack.getParamOut() = temporaryTrack.getParamIn();
621 temporaryTrack.resetCovariance();
622 temporaryTrack.setChi2(0);
623 fitSuccess = fitTrack(temporaryTrack, mTrkParams[0].NLayers - 1, -1, -1, mTrkParams[0].MaxChi2ClusterAttachment, mTrkParams[0].MaxChi2NDF, 50.f);
624 if (!fitSuccess || temporaryTrack.getPt() < mTrkParams[iteration].MinPt[mTrkParams[iteration].NLayers - temporaryTrack.getNClusters()]) {
625 continue;
626 }
627 tracks[trackIndex++] = temporaryTrack;
628 }
629
630 tracks.resize(trackIndex);
631 std::sort(tracks.begin(), tracks.end(), [](const TrackITSExt& a, const TrackITSExt& b) {
632 return a.getChi2() < b.getChi2();
633 });
634
635 for (auto& track : tracks) {
636 int nShared = 0;
637 bool isFirstShared{false};
638 for (int iLayer{0}; iLayer < mTrkParams[0].NLayers; ++iLayer) {
639 if (track.getClusterIndex(iLayer) == constants::its::UnusedIndex) {
640 continue;
641 }
642 nShared += int(mTimeFrame->isClusterUsed(iLayer, track.getClusterIndex(iLayer)));
643 isFirstShared |= !iLayer && mTimeFrame->isClusterUsed(iLayer, track.getClusterIndex(iLayer));
644 }
645
646 if (nShared > mTrkParams[0].ClusterSharing) {
647 continue;
648 }
649
650 std::array<int, 3> rofs{INT_MAX, INT_MAX, INT_MAX};
651 for (int iLayer{0}; iLayer < mTrkParams[0].NLayers; ++iLayer) {
652 if (track.getClusterIndex(iLayer) == constants::its::UnusedIndex) {
653 continue;
654 }
655 mTimeFrame->markUsedCluster(iLayer, track.getClusterIndex(iLayer));
656 int currentROF = mTimeFrame->getClusterROF(iLayer, track.getClusterIndex(iLayer));
657 for (int iR{0}; iR < 3; ++iR) {
658 if (rofs[iR] == INT_MAX) {
659 rofs[iR] = currentROF;
660 }
661 if (rofs[iR] == currentROF) {
662 break;
663 }
664 }
665 }
666 if (rofs[2] != INT_MAX) {
667 continue;
668 }
669 track.setUserField(0);
670 track.getParamOut().setUserField(0);
671 if (rofs[1] != INT_MAX) {
672 track.setNextROFbit();
673 }
674 mTimeFrame->getTracks(o2::gpu::CAMath::Min(rofs[0], rofs[1])).emplace_back(track);
675 }
676 }
677}
678
679void TrackerTraits::extendTracks(const int iteration)
680{
681 for (int rof{0}; rof < mTimeFrame->getNrof(); ++rof) {
682 for (auto& track : mTimeFrame->getTracks(rof)) {
683 auto backup{track};
684 bool success{false};
685 // the order here biases towards top extension, tracks should probably be fitted separately in the directions and then compared.
686 if ((mTrkParams[iteration].UseTrackFollowerMix || mTrkParams[iteration].UseTrackFollowerTop) && track.getLastClusterLayer() != mTrkParams[iteration].NLayers - 1) {
687 success = success || trackFollowing(&track, rof, true, iteration);
688 }
689 if ((mTrkParams[iteration].UseTrackFollowerMix || (mTrkParams[iteration].UseTrackFollowerBot && !success)) && track.getFirstClusterLayer() != 0) {
690 success = success || trackFollowing(&track, rof, false, iteration);
691 }
692 if (success) {
694 track.resetCovariance();
695 track.setChi2(0);
696 bool fitSuccess = fitTrack(track, 0, mTrkParams[iteration].NLayers, 1, mTrkParams[iteration].MaxChi2ClusterAttachment, mTrkParams[0].MaxChi2NDF);
697 if (!fitSuccess) {
698 track = backup;
699 continue;
700 }
701 track.getParamOut() = track;
702 track.resetCovariance();
703 track.setChi2(0);
704 fitSuccess = fitTrack(track, mTrkParams[iteration].NLayers - 1, -1, -1, mTrkParams[iteration].MaxChi2ClusterAttachment, mTrkParams[0].MaxChi2NDF, 50.);
705 if (!fitSuccess) {
706 track = backup;
707 continue;
708 }
710 mTimeFrame->mNExtendedUsedClusters += track.getNClusters() - backup.getNClusters();
711 auto pattern = track.getPattern();
712 auto diff = (pattern & ~backup.getPattern()) & 0xff;
713 pattern |= (diff << 24);
714 track.setPattern(pattern);
716 for (int iLayer{0}; iLayer < mTrkParams[iteration].NLayers; ++iLayer) {
717 if (track.getClusterIndex(iLayer) == constants::its::UnusedIndex) {
718 continue;
719 }
720 mTimeFrame->markUsedCluster(iLayer, track.getClusterIndex(iLayer));
721 }
722 }
723 }
724 }
725}
726
728{
729 const auto propagator = o2::base::Propagator::Instance();
731
732 for (auto& cell : mTimeFrame->getCells()[0]) {
733 auto& cluster3_glo = mTimeFrame->getClusters()[2][cell.getThirdClusterIndex()];
734 auto& cluster2_glo = mTimeFrame->getClusters()[1][cell.getSecondClusterIndex()];
735 auto& cluster1_glo = mTimeFrame->getClusters()[0][cell.getFirstClusterIndex()];
736 if (mTimeFrame->isClusterUsed(2, cluster1_glo.clusterId) ||
737 mTimeFrame->isClusterUsed(1, cluster2_glo.clusterId) ||
738 mTimeFrame->isClusterUsed(0, cluster3_glo.clusterId)) {
739 continue;
740 }
741
742 std::array<int, 3> rofs{
743 mTimeFrame->getClusterROF(2, cluster3_glo.clusterId),
744 mTimeFrame->getClusterROF(1, cluster2_glo.clusterId),
745 mTimeFrame->getClusterROF(0, cluster1_glo.clusterId)};
746 if (rofs[0] != rofs[1] && rofs[1] != rofs[2] && rofs[0] != rofs[2]) {
747 continue;
748 }
749
750 int rof{rofs[0]};
751 if (rofs[1] == rofs[2]) {
752 rof = rofs[2];
753 }
754
755 auto pvs{mTimeFrame->getPrimaryVertices(rof)};
756 auto pvsXAlpha{mTimeFrame->getPrimaryVerticesXAlpha(rof)};
757
758 const auto& cluster3_tf = mTimeFrame->getTrackingFrameInfoOnLayer(2).at(cluster3_glo.clusterId);
759 TrackITSExt temporaryTrack{buildTrackSeed(cluster1_glo, cluster2_glo, cluster3_tf)};
760 temporaryTrack.setExternalClusterIndex(0, cluster1_glo.clusterId, true);
761 temporaryTrack.setExternalClusterIndex(1, cluster2_glo.clusterId, true);
762 temporaryTrack.setExternalClusterIndex(2, cluster3_glo.clusterId, true);
763
765 bool fitSuccess = fitTrack(temporaryTrack, 1, -1, -1);
766 if (!fitSuccess) {
767 continue;
768 }
769 fitSuccess = false;
770
771 TrackITSExt bestTrack{temporaryTrack}, backup{temporaryTrack};
772 float bestChi2{std::numeric_limits<float>::max()};
773 for (int iV{0}; iV < (int)pvs.size(); ++iV) {
774 temporaryTrack = backup;
775 if (!temporaryTrack.rotate(pvsXAlpha[iV][1])) {
776 continue;
777 }
778 if (!propagator->propagateTo(temporaryTrack, pvsXAlpha[iV][0], true)) {
779 continue;
780 }
781
782 float pvRes{mTrkParams[0].PVres / o2::gpu::CAMath::Sqrt(float(pvs[iV].getNContributors()))};
783 const float posVtx[2]{0.f, pvs[iV].getZ()};
784 const float covVtx[3]{pvRes, 0.f, pvRes};
785 float chi2 = temporaryTrack.getPredictedChi2Quiet(posVtx, covVtx);
786 if (chi2 < bestChi2) {
787 if (!temporaryTrack.track::TrackParCov::update(posVtx, covVtx)) {
788 continue;
789 }
790 bestTrack = temporaryTrack;
791 bestChi2 = chi2;
792 }
793 }
794
795 bestTrack.resetCovariance();
796 bestTrack.setChi2(0.f);
797 fitSuccess = fitTrack(bestTrack, 0, mTrkParams[0].NLayers, 1, mTrkParams[0].MaxChi2ClusterAttachment, mTrkParams[0].MaxChi2NDF);
798 if (!fitSuccess) {
799 continue;
800 }
801 bestTrack.getParamOut() = bestTrack;
802 bestTrack.resetCovariance();
803 bestTrack.setChi2(0.f);
804 fitSuccess = fitTrack(bestTrack, mTrkParams[0].NLayers - 1, -1, -1, mTrkParams[0].MaxChi2ClusterAttachment, mTrkParams[0].MaxChi2NDF, 50.);
805 if (!fitSuccess) {
806 continue;
807 }
808 mTimeFrame->markUsedCluster(0, bestTrack.getClusterIndex(0));
809 mTimeFrame->markUsedCluster(1, bestTrack.getClusterIndex(1));
810 mTimeFrame->markUsedCluster(2, bestTrack.getClusterIndex(2));
811 mTimeFrame->getTracks(rof).emplace_back(bestTrack);
812 }
813}
814
815bool TrackerTraits::fitTrack(TrackITSExt& track, int start, int end, int step, float chi2clcut, float chi2ndfcut, float maxQoverPt, int nCl)
816{
817 auto propInstance = o2::base::Propagator::Instance();
818
819 for (int iLayer{start}; iLayer != end; iLayer += step) {
820 if (track.getClusterIndex(iLayer) == constants::its::UnusedIndex) {
821 continue;
822 }
823 const TrackingFrameInfo& trackingHit = mTimeFrame->getTrackingFrameInfoOnLayer(iLayer).at(track.getClusterIndex(iLayer));
824
825 if (!track.rotate(trackingHit.alphaTrackingFrame)) {
826 return false;
827 }
828
829 if (!propInstance->propagateToX(track, trackingHit.xTrackingFrame, getBz(), o2::base::PropagatorImpl<float>::MAX_SIN_PHI, o2::base::PropagatorImpl<float>::MAX_STEP, mCorrType)) {
830 return false;
831 }
832
833 if (mCorrType == o2::base::PropagatorF::MatCorrType::USEMatCorrNONE) {
834 constexpr float radl = 9.36f; // Radiation length of Si [cm]
835 constexpr float rho = 2.33f; // Density of Si [g/cm^3]
836 if (!track.correctForMaterial(mTrkParams[0].LayerxX0[iLayer], mTrkParams[0].LayerxX0[iLayer] * radl * rho, true)) {
837 continue;
838 }
839 }
840
841 auto predChi2{track.getPredictedChi2Quiet(trackingHit.positionTrackingFrame, trackingHit.covarianceTrackingFrame)};
842 if ((nCl >= 3 && predChi2 > chi2clcut) || predChi2 < 0.f) {
843 return false;
844 }
845 track.setChi2(track.getChi2() + predChi2);
846 if (!track.o2::track::TrackParCov::update(trackingHit.positionTrackingFrame, trackingHit.covarianceTrackingFrame)) {
847 return false;
848 }
849 nCl++;
850 }
851 return std::abs(track.getQ2Pt()) < maxQoverPt && track.getChi2() < chi2ndfcut * (nCl * 2 - 5);
852}
853
854bool TrackerTraits::trackFollowing(TrackITSExt* track, int rof, bool outward, const int iteration)
855{
856 auto propInstance = o2::base::Propagator::Instance();
857 const int step = -1 + outward * 2;
858 const int end = outward ? mTrkParams[iteration].NLayers - 1 : 0;
859 std::vector<TrackITSExt> hypotheses(1, *track); // possibly avoid reallocation
860 for (size_t iHypo{0}; iHypo < hypotheses.size(); ++iHypo) {
861 auto hypo{hypotheses[iHypo]};
862 int iLayer = static_cast<int>(outward ? hypo.getLastClusterLayer() : hypo.getFirstClusterLayer());
863 // per layer we add new hypotheses
864 while (iLayer != end) {
865 iLayer += step; // step through all layers until we reach the end, this allows for skipping on empty layers
866 const float r = mTrkParams[iteration].LayerRadii[iLayer];
867 // get an estimate of the trackinf-frame x for the next step
868 float x{-999};
869 if (!hypo.getXatLabR(r, x, mTimeFrame->getBz(), o2::track::DirAuto) || x <= 0.f) {
870 continue;
871 }
872 // estimate hypo's trk parameters at that x
873 auto& hypoParam{outward ? hypo.getParamOut() : hypo.getParamIn()};
874 if (!propInstance->propagateToX(hypoParam, x, mTimeFrame->getBz(), PropagatorF::MAX_SIN_PHI,
875 PropagatorF::MAX_STEP, mTrkParams[iteration].CorrType)) {
876 continue;
877 }
878
879 if (mTrkParams[iteration].CorrType == PropagatorF::MatCorrType::USEMatCorrNONE) { // account for material affects if propagator does not
880 constexpr float radl = 9.36f; // Radiation length of Si [cm]
881 constexpr float rho = 2.33f; // Density of Si [g/cm^3]
882 if (!hypoParam.correctForMaterial(mTrkParams[iteration].LayerxX0[iLayer], mTrkParams[iteration].LayerxX0[iLayer] * radl * rho, true)) {
883 continue;
884 }
885 }
886
887 // calculate the search window on this layer
888 const float phi{hypoParam.getPhi()};
889 const float ePhi{o2::gpu::CAMath::Sqrt(hypoParam.getSigmaSnp2() / hypoParam.getCsp2())};
890 const float z{hypoParam.getZ()};
891 const float eZ{o2::gpu::CAMath::Sqrt(hypoParam.getSigmaZ2())};
892 const int4 selectedBinsRect{getBinsRect(iLayer, phi, mTrkParams[iteration].TrackFollowerNSigmaCutPhi * ePhi, z, mTrkParams[iteration].TrackFollowerNSigmaCutZ * eZ)};
893 if (selectedBinsRect.x == 0 && selectedBinsRect.y == 0 && selectedBinsRect.z == 0 && selectedBinsRect.w == 0) {
894 continue;
895 }
896
897 int phiBinsNum{selectedBinsRect.w - selectedBinsRect.y + 1};
898
899 if (phiBinsNum < 0) {
900 phiBinsNum += mTrkParams[iteration].PhiBins;
901 }
902
903 gsl::span<const Cluster> layer1 = mTimeFrame->getClustersOnLayer(rof, iLayer);
904 if (layer1.empty()) {
905 continue;
906 }
907
908 // check all clusters in search windows for possible new hypotheses
909 for (int iPhiCount = 0; iPhiCount < phiBinsNum; iPhiCount++) {
910 int iPhiBin = (selectedBinsRect.y + iPhiCount) % mTrkParams[iteration].PhiBins;
911 const int firstBinIndex{mTimeFrame->mIndexTableUtils.getBinIndex(selectedBinsRect.x, iPhiBin)};
912 const int maxBinIndex{firstBinIndex + selectedBinsRect.z - selectedBinsRect.x + 1};
913 const int firstRowClusterIndex = mTimeFrame->getIndexTable(rof, iLayer)[firstBinIndex];
914 const int maxRowClusterIndex = mTimeFrame->getIndexTable(rof, iLayer)[maxBinIndex];
915
916 for (int iNextCluster{firstRowClusterIndex}; iNextCluster < maxRowClusterIndex; ++iNextCluster) {
917 if (iNextCluster >= (int)layer1.size()) {
918 break;
919 }
920 const Cluster& nextCluster{layer1[iNextCluster]};
921
922 if (mTimeFrame->isClusterUsed(iLayer, nextCluster.clusterId)) {
923 continue;
924 }
925
926 const TrackingFrameInfo& trackingHit = mTimeFrame->getTrackingFrameInfoOnLayer(iLayer).at(nextCluster.clusterId);
927
928 auto tbupdated{hypo};
929 auto& tbuParams = outward ? tbupdated.getParamOut() : tbupdated.getParamIn();
930 if (!tbuParams.rotate(trackingHit.alphaTrackingFrame)) {
931 continue;
932 }
933
934 if (!propInstance->propagateToX(tbuParams, trackingHit.xTrackingFrame, mTimeFrame->getBz(),
935 PropagatorF::MAX_SIN_PHI, PropagatorF::MAX_STEP, PropagatorF::MatCorrType::USEMatCorrNONE)) {
936 continue;
937 }
938
939 auto predChi2{tbuParams.getPredictedChi2Quiet(trackingHit.positionTrackingFrame, trackingHit.covarianceTrackingFrame)};
940 if (predChi2 >= track->getChi2() * mTrkParams[iteration].NSigmaCut) {
941 continue;
942 }
943
944 if (!tbuParams.o2::track::TrackParCov::update(trackingHit.positionTrackingFrame, trackingHit.covarianceTrackingFrame)) {
945 continue;
946 }
947 tbupdated.setChi2(tbupdated.getChi2() + predChi2);
948 tbupdated.setExternalClusterIndex(iLayer, nextCluster.clusterId, true);
949 hypotheses.emplace_back(tbupdated);
950 }
951 }
952 }
953 }
954
955 TrackITSExt* bestHypo{track};
956 bool swapped{false};
957 for (auto& hypo : hypotheses) {
958 if (hypo.isBetter(*bestHypo, track->getChi2() * mTrkParams[iteration].NSigmaCut)) {
959 bestHypo = &hypo;
960 swapped = true;
961 }
962 }
963 *track = *bestHypo;
964 return swapped;
965}
966
969track::TrackParCov TrackerTraits::buildTrackSeed(const Cluster& cluster1, const Cluster& cluster2, const TrackingFrameInfo& tf3)
970{
971 const float ca = o2::gpu::CAMath::Cos(tf3.alphaTrackingFrame), sa = o2::gpu::CAMath::Sin(tf3.alphaTrackingFrame);
972 const float x1 = cluster1.xCoordinate * ca + cluster1.yCoordinate * sa;
973 const float y1 = -cluster1.xCoordinate * sa + cluster1.yCoordinate * ca;
974 const float z1 = cluster1.zCoordinate;
975 const float x2 = cluster2.xCoordinate * ca + cluster2.yCoordinate * sa;
976 const float y2 = -cluster2.xCoordinate * sa + cluster2.yCoordinate * ca;
977 const float z2 = cluster2.zCoordinate;
978 const float x3 = tf3.xTrackingFrame;
979 const float y3 = tf3.positionTrackingFrame[0];
980 const float z3 = tf3.positionTrackingFrame[1];
981
982 const bool zeroField{std::abs(getBz()) < o2::constants::math::Almost0};
983 const float tgp = zeroField ? o2::gpu::CAMath::ATan2(y3 - y1, x3 - x1) : 1.f;
984 const float crv = zeroField ? 1.f : math_utils::computeCurvature(x3, y3, x2, y2, x1, y1);
985 const float snp = zeroField ? tgp / o2::gpu::CAMath::Sqrt(1.f + tgp * tgp) : crv * (x3 - math_utils::computeCurvatureCentreX(x3, y3, x2, y2, x1, y1));
986 const float tgl12 = math_utils::computeTanDipAngle(x1, y1, x2, y2, z1, z2);
987 const float tgl23 = math_utils::computeTanDipAngle(x2, y2, x3, y3, z2, z3);
988 const float q2pt = zeroField ? 1.f / o2::track::kMostProbablePt : crv / (getBz() * o2::constants::math::B2C);
989 const float q2pt2 = crv * crv;
990 const float sg2q2pt = track::kC1Pt2max * (q2pt2 > 0.0005 ? (q2pt2 < 1 ? q2pt2 : 1) : 0.0005);
991 return track::TrackParCov(tf3.xTrackingFrame, tf3.alphaTrackingFrame,
992 {y3, z3, snp, 0.5f * (tgl12 + tgl23), q2pt},
993 {tf3.covarianceTrackingFrame[0],
994 tf3.covarianceTrackingFrame[1], tf3.covarianceTrackingFrame[2],
995 0.f, 0.f, track::kCSnp2max,
996 0.f, 0.f, 0.f, track::kCTgl2max,
997 0.f, 0.f, 0.f, 0.f, sg2q2pt});
998}
999
1001{
1002 mBz = bz;
1003 mTimeFrame->setBz(bz);
1004}
1005
1007
1009{
1010#ifdef WITH_OPENMP
1011 mNThreads = n > 0 ? n : 1;
1012#else
1013 mNThreads = 1;
1014#endif
1015}
1016
1021
1026
1028{
1029 return mTimeFrame->getNumberOfCells();
1030}
1031
1036
1037} // namespace its
1038} // namespace o2
Base track model for the Barrel, params only, w/o covariance.
#define CA_DEBUGGER(x)
Definition Definitions.h:25
void print() const
#define failed(...)
Definition Utils.h:42
useful math constants
GPUd() value_type estimateLTFast(o2 static GPUd() float estimateLTIncrement(const o2 PropagatorImpl * Instance(bool uninitialized=false)
Definition Propagator.h:143
gsl::span< const Vertex > getPrimaryVertices(int rofId) const
Definition TimeFrame.h:355
void fillPrimaryVerticesXandAlpha()
bool isClusterUsed(int layer, int clusterId) const
Definition TimeFrame.h:570
IndexTableUtils mIndexTableUtils
Definition TimeFrame.h:264
std::vector< std::vector< int > > & getCellsLookupTable()
Definition TimeFrame.h:650
void markUsedCluster(int layer, int clusterId)
Definition TimeFrame.h:580
std::vector< std::vector< CellSeed > > & getCells()
Definition TimeFrame.h:648
std::vector< std::vector< Cluster > > & getClusters()
Definition TimeFrame.h:638
gsl::span< int > getIndexTable(int rofId, int layerId)
Definition TimeFrame.h:529
float getBz() const
Definition TimeFrame.h:223
std::vector< std::vector< Cluster > > & getUnsortedClusters()
Definition TimeFrame.h:643
int getNumberOfTracklets() const
Definition TimeFrame.h:705
int getNumberOfClusters() const
Definition TimeFrame.h:687
int getClusterROF(int iLayer, int iCluster)
Definition TimeFrame.h:485
std::vector< MCCompLabel > & getCellsLabel(int layer)
Definition TimeFrame.h:142
std::vector< TrackITSExt > & getTracks(int rofId)
Definition TimeFrame.h:169
std::vector< std::vector< int > > & getCellsNeighboursLUT()
Definition TimeFrame.h:656
gsl::span< Cluster > getClustersOnLayer(int rofId, int layerId)
Definition TimeFrame.h:413
int getNrof() const
Definition TimeFrame.h:395
std::vector< std::vector< int > > & getCellsNeighbours()
Definition TimeFrame.h:655
const std::vector< TrackingFrameInfo > & getTrackingFrameInfoOnLayer(int layerId) const
Definition TimeFrame.h:499
gsl::span< const std::array< float, 2 > > getPrimaryVerticesXAlpha(int rofId) const
Definition TimeFrame.h:376
int getNumberOfCells() const
Definition TimeFrame.h:696
void setBz(float bz)
Definition TimeFrame.h:222
virtual void findRoads(const int iteration)
virtual void processNeighbours(int iLayer, int iLevel, const std::vector< CellSeed > &currentCellSeed, const std::vector< int > &currentCellId, std::vector< CellSeed > &updatedCellSeed, std::vector< int > &updatedCellId)
virtual void extendTracks(const int iteration)
virtual int getTFNumberOfCells() const
virtual void adoptTimeFrame(TimeFrame *tf)
virtual void computeLayerTracklets(const int iteration, int iROFslice, int iVertex)
virtual void setBz(float bz)
virtual int getTFNumberOfTracklets() const
virtual void findCellsNeighbours(const int iteration)
virtual int getTFNumberOfClusters() const
virtual void findShortPrimaries()
const int4 getBinsRect(const Cluster &, int layer, float z1, float z2, float maxdeltaz, float maxdeltaphi) const noexcept
std::vector< TrackingParameters > mTrkParams
o2::base::PropagatorImpl< float >::MatCorrType mCorrType
virtual bool trackFollowing(TrackITSExt *track, int rof, bool outward, const int iteration)
virtual void computeLayerCells(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
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 Almost0
constexpr float B2C
constexpr int UnusedIndex
Definition Constants.h:54
constexpr float TwoPi
Definition Constants.h:46
constexpr int debugLevel
TrackParCovF TrackParCov
Definition Track.h:33
constexpr float kMostProbablePt
value_T rho
Definition TrackUtils.h:36
constexpr float kC1Pt2max
a couple of static helper functions to create timestamp values for CCDB queries or override obsolete ...
std::unique_ptr< GPUReconstructionTimeframe > tf
int32_t w
float yCoordinate
Definition Cluster.h:45
float zCoordinate
Definition Cluster.h:46
float xCoordinate
Definition Cluster.h:44
std::array< uint16_t, 5 > pattern