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);
91 float meanDeltaR{mTrkParams[iteration].LayerRadii[iLayer + 1] - mTrkParams[iteration].LayerRadii[iLayer]};
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)};
98 if (mTimeFrame->isClusterUsed(iLayer, currentCluster.clusterId)) {
101 const float inverseR0{1.f / currentCluster.
radius};
103 for (
int iV{startVtx}; iV < endVtx; ++iV) {
104 const auto& primaryVertex{primaryVertices[iV]};
105 if (primaryVertex.isFlagSet(2) && iteration != 3) {
108 const float resolution = o2::gpu::CAMath::Sqrt(math_utils::Sq(mTrkParams[iteration].PVres) / primaryVertex.getNContributors() + math_utils::Sq(mTimeFrame->getPositionResolution(iLayer)));
110 const float tanLambda{(currentCluster.zCoordinate - primaryVertex.getZ()) * inverseR0};
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};
115 const float sqInverseDeltaZ0{1.f / (math_utils::Sq(currentCluster.zCoordinate - primaryVertex.getZ()) + 2.e-8f)};
116 const float sigmaZ{o2::gpu::CAMath::Sqrt(math_utils::Sq(resolution) * math_utils::Sq(tanLambda) * ((math_utils::Sq(inverseR0) + sqInverseDeltaZ0) * math_utils::Sq(meanDeltaR) + 1.f) + math_utils::Sq(meanDeltaR * mTimeFrame->getMSangle(iLayer)))};
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) {
123 int phiBinsNum{selectedBinsRect.
w - selectedBinsRect.y + 1};
125 if (phiBinsNum < 0) {
126 phiBinsNum += mTrkParams[iteration].PhiBins;
129 for (
int rof1{minRof}; rof1 <= maxRof; ++rof1) {
130 auto layer1 = mTimeFrame->getClustersOnLayer(rof1, iLayer + 1);
131 if (layer1.empty()) {
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;
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()) {
155 const Cluster& nextCluster{layer1[iNextCluster]};
156 if (mTimeFrame->isClusterUsed(iLayer + 1, nextCluster.clusterId)) {
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)};
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()) {
175 if (
label.isValid()) {
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;
182 if (deltaZ / sigmaZ < mTrkParams[iteration].NSigmaCut &&
183 (deltaPhi < mTimeFrame->getPhiCut(iLayer) ||
186 mTimeFrame->getTrackletsLookupTable()[iLayer - 1][currentSortedIndex]++;
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);
204 tbb::blocked_range<int>(0, mTrkParams[iteration].TrackletsPerRoad()),
205 [&](
const tbb::blocked_range<int>& Layers) {
206 for (
int iLayer = Layers.begin(); iLayer < Layers.end(); ++iLayer) {
208 auto& trkl{mTimeFrame->getTracklets()[iLayer]};
209 tbb::parallel_sort(trkl.begin(), trkl.end(), [](
const Tracklet&
a,
const Tracklet&
b) ->
bool {
210 return a.firstClusterIndex < b.firstClusterIndex || (a.firstClusterIndex == b.firstClusterIndex && a.secondClusterIndex < b.secondClusterIndex);
213 trkl.erase(std::unique(trkl.begin(), trkl.end(), [](
const Tracklet&
a,
const Tracklet&
b) ->
bool {
214 return a.firstClusterIndex == b.firstClusterIndex && a.secondClusterIndex == b.secondClusterIndex;
217 trkl.shrink_to_fit();
219 auto& lut{mTimeFrame->getTrackletsLookupTable()[iLayer - 1]};
220 std::fill(lut.begin(), lut.end(), 0);
224 for (
const auto& tkl : trkl) {
225 lut[tkl.firstClusterIndex]++;
227 std::exclusive_scan(lut.begin(), lut.end(), lut.begin(), 0);
228 lut.push_back(trkl.size());
235 if (mTimeFrame->hasMCinformation()) {
236 for (
int iLayer{0}; iLayer < mTrkParams[iteration].TrackletsPerRoad(); ++iLayer) {
237 for (
auto& trk : mTimeFrame->getTracklets()[iLayer]) {
239 int currentId{mTimeFrame->getClusters()[iLayer][trk.firstClusterIndex].clusterId};
240 int nextId{mTimeFrame->getClusters()[iLayer + 1][trk.secondClusterIndex].clusterId};
241 for (
auto& lab1 : mTimeFrame->getClusterLabels(iLayer, currentId)) {
242 for (
auto& lab2 : mTimeFrame->getClusterLabels(iLayer + 1, nextId)) {
243 if (lab1 == lab2 && lab1.isValid()) {
248 if (
label.isValid()) {
252 mTimeFrame->getTrackletsLabel(iLayer).emplace_back(
label);
258template <
int nLayers>
261#ifdef OPTIMISATION_OUTPUT
263 std::ofstream off(std::format(
"cells{}.txt", iter++));
266 for (
int iLayer = 0; iLayer < mTrkParams[iteration].CellsPerRoad(); ++iLayer) {
271 if (mTimeFrame->hasMCinformation()) {
276 mTaskArena->execute([&] {
278 const Tracklet& currentTracklet{mTimeFrame->getTracklets()[iLayer][iTracklet]};
280 const int nextLayerFirstTrackletIndex{
281 mTimeFrame->getTrackletsLookupTable()[iLayer][nextLayerClusterIndex]};
282 const int nextLayerLastTrackletIndex{
283 mTimeFrame->getTrackletsLookupTable()[iLayer][nextLayerClusterIndex + 1]};
286 for (
int iNextTracklet{nextLayerFirstTrackletIndex}; iNextTracklet < nextLayerLastTrackletIndex; ++iNextTracklet) {
287 if (mTimeFrame->getTracklets()[iLayer + 1][iNextTracklet].firstClusterIndex != nextLayerClusterIndex) {
290 const Tracklet& nextTracklet{mTimeFrame->getTracklets()[iLayer + 1][iNextTracklet]};
291 const float deltaTanLambda{std::abs(currentTracklet.tanLambda - nextTracklet.tanLambda)};
293#ifdef OPTIMISATION_OUTPUT
294 float resolution{o2::gpu::CAMath::Sqrt(0.5f * (mTrkParams[iteration].SystErrorZ2[iLayer] + mTrkParams[iteration].SystErrorZ2[iLayer + 1] + mTrkParams[iteration].SystErrorZ2[iLayer + 2] + mTrkParams[iteration].SystErrorY2[iLayer] + mTrkParams[iteration].SystErrorY2[iLayer + 1] + mTrkParams[iteration].SystErrorY2[iLayer + 2])) / mTrkParams[iteration].LayerResolution[iLayer]};
295 resolution = resolution > 1.e-12 ? resolution : 1.f;
296 bool good{mTimeFrame->getTrackletsLabel(iLayer)[iTracklet] == mTimeFrame->getTrackletsLabel(iLayer + 1)[iNextTracklet]};
297 float signedDelta{currentTracklet.
tanLambda - nextTracklet.tanLambda};
298 off << std::format(
"{}\t{:d}\t{}\t{}\t{}\t{}", iLayer, good, signedDelta, signedDelta / (mTrkParams[iteration].CellDeltaTanLambdaSigma), tanLambda, resolution) << std::endl;
301 if (deltaTanLambda / mTrkParams[iteration].CellDeltaTanLambdaSigma < mTrkParams[iteration].NSigmaCut) {
305 mTimeFrame->getClusters()[iLayer][currentTracklet.firstClusterIndex].clusterId,
306 mTimeFrame->getClusters()[iLayer + 1][nextTracklet.firstClusterIndex].clusterId,
307 mTimeFrame->getClusters()[iLayer + 2][nextTracklet.secondClusterIndex].clusterId};
308 const auto& cluster1_glo = mTimeFrame->getUnsortedClusters()[iLayer][clusId[0]];
309 const auto& cluster2_glo = mTimeFrame->getUnsortedClusters()[iLayer + 1][clusId[1]];
310 const auto& cluster3_tf = mTimeFrame->getTrackingFrameInfoOnLayer(iLayer + 2)[clusId[2]];
311 auto track{buildTrackSeed(cluster1_glo, cluster2_glo, cluster3_tf)};
315 for (
int iC{2}; iC--;) {
316 const TrackingFrameInfo& trackingHit = mTimeFrame->getTrackingFrameInfoOnLayer(iLayer + iC)[clusId[iC]];
326 if (!track.correctForMaterial(mTrkParams[0].LayerxX0[iLayer + iC], mTrkParams[0].LayerxX0[iLayer] *
constants::Radl *
constants::Rho,
true)) {
331 if (!iC && predChi2 > mTrkParams[iteration].MaxChi2ClusterAttachment) {
343 if constexpr (
decltype(Tag)
::value == PassMode::OnePass::value) {
344 layerCells.emplace_back(iLayer, clusId[0], clusId[1], clusId[2], iTracklet, iNextTracklet, track, chi2);
346 }
else if constexpr (
decltype(Tag)
::value == PassMode::TwoPassCount::value) {
348 }
else if constexpr (
decltype(Tag)
::value == PassMode::TwoPassInsert::value) {
349 layerCells[
offset++] =
CellSeed(iLayer, clusId[0], clusId[1], clusId[2], iTracklet, iNextTracklet, track, chi2);
351 static_assert(
false,
"Unknown mode!");
360 tbb::blocked_range<int>(0, mTrkParams[iteration].CellsPerRoad()),
361 [&](
const tbb::blocked_range<int>& Layers) {
362 for (
int iLayer = Layers.begin(); iLayer < Layers.end(); ++iLayer) {
363 if (mTimeFrame->getTracklets()[iLayer + 1].empty() ||
364 mTimeFrame->getTracklets()[iLayer].empty()) {
368 auto& layerCells = mTimeFrame->getCells()[iLayer];
369 const int currentLayerTrackletsNum{
static_cast<int>(mTimeFrame->getTracklets()[iLayer].size())};
371 if (mTaskArena->max_concurrency() <= 1) {
372 for (
int iTracklet{0}; iTracklet < currentLayerTrackletsNum; ++iTracklet) {
373 perTrackletCount[iTracklet] = forTrackletCells(
PassMode::OnePass{}, iLayer, layerCells, iTracklet);
375 std::exclusive_scan(perTrackletCount.begin(), perTrackletCount.end(), perTrackletCount.begin(), 0);
378 tbb::blocked_range<int>(0, currentLayerTrackletsNum),
379 [&](
const tbb::blocked_range<int>& Tracklets) {
380 for (
int iTracklet = Tracklets.begin(); iTracklet < Tracklets.end(); ++iTracklet) {
385 std::exclusive_scan(perTrackletCount.begin(), perTrackletCount.end(), perTrackletCount.begin(), 0);
386 auto totalCells{perTrackletCount.back()};
387 if (totalCells == 0) {
390 layerCells.resize(totalCells);
393 tbb::blocked_range<int>(0, currentLayerTrackletsNum),
394 [&](
const tbb::blocked_range<int>& Tracklets) {
395 for (
int iTracklet = Tracklets.begin(); iTracklet < Tracklets.end(); ++iTracklet) {
396 int offset = perTrackletCount[iTracklet];
397 if (
offset == perTrackletCount[iTracklet + 1]) {
406 auto& lut = mTimeFrame->getCellsLookupTable()[iLayer - 1];
407 lut.resize(currentLayerTrackletsNum + 1);
408 std::copy_n(perTrackletCount.begin(), currentLayerTrackletsNum + 1, lut.begin());
415 if (mTimeFrame->hasMCinformation()) {
416 for (
int iLayer{0}; iLayer < mTrkParams[iteration].CellsPerRoad(); ++iLayer) {
417 for (
auto& cell : mTimeFrame->getCells()[iLayer]) {
418 MCCompLabel currentLab{mTimeFrame->getTrackletsLabel(iLayer)[cell.getFirstTrackletIndex()]};
419 MCCompLabel nextLab{mTimeFrame->getTrackletsLabel(iLayer + 1)[cell.getSecondTrackletIndex()]};
420 mTimeFrame->getCellsLabel(iLayer).emplace_back(currentLab == nextLab ? currentLab :
MCCompLabel());
425 if constexpr (debugLevel) {
426 for (
int iLayer{0}; iLayer < mTrkParams[iteration].CellsPerRoad(); ++iLayer) {
427 std::cout <<
"Cells on layer " << iLayer <<
" " << mTimeFrame->getCells()[iLayer].size() << std::endl;
432template <
int nLayers>
435#ifdef OPTIMISATION_OUTPUT
436 std::ofstream off(std::format(
"cellneighs{}.txt", iteration));
440 int cell{-1}, nextCell{-1},
level{-1};
443 mTaskArena->execute([&] {
444 for (
int iLayer{0}; iLayer < mTrkParams[iteration].CellsPerRoad() - 1; ++iLayer) {
447 if (mTimeFrame->getCells()[iLayer + 1].empty() ||
448 mTimeFrame->getCellsLookupTable()[iLayer].empty()) {
452 int nCells{
static_cast<int>(mTimeFrame->getCells()[iLayer].size())};
455 auto forCellNeighbour = [&](
auto Tag,
int iCell,
int offset = 0) ->
int {
456 const auto& currentCellSeed{mTimeFrame->getCells()[iLayer][iCell]};
457 const int nextLayerTrackletIndex{currentCellSeed.getSecondTrackletIndex()};
458 const int nextLayerFirstCellIndex{mTimeFrame->getCellsLookupTable()[iLayer][nextLayerTrackletIndex]};
459 const int nextLayerLastCellIndex{mTimeFrame->getCellsLookupTable()[iLayer][nextLayerTrackletIndex + 1]};
460 int foundNextCells{0};
461 for (
int iNextCell{nextLayerFirstCellIndex}; iNextCell < nextLayerLastCellIndex; ++iNextCell) {
462 auto nextCellSeed{mTimeFrame->getCells()[iLayer + 1][iNextCell]};
463 if (nextCellSeed.getFirstTrackletIndex() != nextLayerTrackletIndex) {
467 if (!nextCellSeed.rotate(currentCellSeed.getAlpha()) ||
468 !nextCellSeed.propagateTo(currentCellSeed.getX(), getBz())) {
471 float chi2 = currentCellSeed.getPredictedChi2(nextCellSeed);
473#ifdef OPTIMISATION_OUTPUT
474 bool good{mTimeFrame->getCellsLabel(iLayer)[iCell] == mTimeFrame->getCellsLabel(iLayer + 1)[iNextCell]};
475 off << std::format(
"{}\t{:d}\t{}", iLayer, good, chi2) << std::endl;
478 if (chi2 > mTrkParams[0].MaxChi2ClusterAttachment) {
482 if constexpr (
decltype(Tag)
::value == PassMode::OnePass::value) {
483 cellsNeighbours.emplace_back(iCell, iNextCell, currentCellSeed.getLevel() + 1);
484 }
else if constexpr (
decltype(Tag)
::value == PassMode::TwoPassCount::value) {
486 }
else if constexpr (
decltype(Tag)
::value == PassMode::TwoPassInsert::value) {
487 cellsNeighbours[
offset++] = {iCell, iNextCell, currentCellSeed.getLevel() + 1};
489 static_assert(
false,
"Unknown mode!");
492 return foundNextCells;
495 if (mTaskArena->max_concurrency() <= 1) {
496 for (
int iCell{0}; iCell < nCells; ++iCell) {
502 tbb::blocked_range<int>(0, nCells),
503 [&](
const tbb::blocked_range<int>& Cells) {
504 for (
int iCell = Cells.begin(); iCell < Cells.end(); ++iCell) {
509 std::exclusive_scan(perCellCount.begin(), perCellCount.end(), perCellCount.begin(), 0);
510 int totalCellNeighbours = perCellCount.back();
511 if (totalCellNeighbours == 0) {
515 cellsNeighbours.resize(totalCellNeighbours);
518 tbb::blocked_range<int>(0, nCells),
519 [&](
const tbb::blocked_range<int>& Cells) {
520 for (
int iCell = Cells.begin(); iCell < Cells.end(); ++iCell) {
521 int offset = perCellCount[iCell];
522 if (
offset == perCellCount[iCell + 1]) {
530 if (cellsNeighbours.empty()) {
534 tbb::parallel_sort(cellsNeighbours.begin(), cellsNeighbours.end(), [](
const auto&
a,
const auto&
b) {
535 return a.nextCell < b.nextCell;
538 auto& cellsNeighbourLUT = mTimeFrame->getCellsNeighboursLUT()[iLayer];
539 cellsNeighbourLUT.assign(mTimeFrame->getCells()[iLayer + 1].size(), 0);
540 for (
const auto& neigh : cellsNeighbours) {
541 ++cellsNeighbourLUT[neigh.nextCell];
543 std::inclusive_scan(cellsNeighbourLUT.begin(), cellsNeighbourLUT.end(), cellsNeighbourLUT.begin());
545 mTimeFrame->getCellsNeighbours()[iLayer].reserve(cellsNeighbours.size());
546 std::ranges::transform(cellsNeighbours, std::back_inserter(mTimeFrame->getCellsNeighbours()[iLayer]), [](
const auto& neigh) { return neigh.cell; });
548 auto it = cellsNeighbours.begin();
549 int current = it->nextCell;
550 int maxLvl = it->level;
552 for (; it != cellsNeighbours.end(); ++it) {
553 if (it->nextCell == current) {
554 maxLvl = std::max(maxLvl, it->level);
556 mTimeFrame->getCells()[iLayer + 1][current].setLevel(maxLvl);
557 current = it->nextCell;
561 mTimeFrame->getCells()[iLayer + 1][current].setLevel(maxLvl);
566template <
int nLayers>
569 CA_DEBUGGER(std::cout <<
"Processing neighbours layer " << iLayer <<
" level " << iLevel <<
", size of the cell seeds: " << currentCellSeed.size() << std::endl);
573 int failed[5]{0, 0, 0, 0, 0}, attempts{0}, failedByMismatch{0};
576 mTaskArena->execute([&] {
577 auto forCellNeighbours = [&](
auto Tag,
int iCell,
int offset = 0) ->
int {
578 const CellSeed& currentCell{currentCellSeed[iCell]};
580 if constexpr (
decltype(Tag)
::value != PassMode::TwoPassInsert::value) {
581 if (currentCell.getLevel() != iLevel) {
584 if (currentCellId.empty() && (mTimeFrame->isClusterUsed(iLayer, currentCell.getFirstClusterIndex()) ||
585 mTimeFrame->isClusterUsed(iLayer + 1, currentCell.getSecondClusterIndex()) ||
586 mTimeFrame->isClusterUsed(iLayer + 2, currentCell.getThirdClusterIndex()))) {
591 const int cellId = currentCellId.empty() ? iCell : currentCellId[iCell];
592 const int startNeighbourId{cellId ? mTimeFrame->getCellsNeighboursLUT()[iLayer - 1][cellId - 1] : 0};
593 const int endNeighbourId{mTimeFrame->getCellsNeighboursLUT()[iLayer - 1][cellId]};
595 for (
int iNeighbourCell{startNeighbourId}; iNeighbourCell < endNeighbourId; ++iNeighbourCell) {
597 const int neighbourCellId = mTimeFrame->getCellsNeighbours()[iLayer - 1][iNeighbourCell];
598 const CellSeed& neighbourCell = mTimeFrame->getCells()[iLayer - 1][neighbourCellId];
599 if (neighbourCell.getSecondTrackletIndex() != currentCell.getFirstTrackletIndex()) {
603 if (mTimeFrame->isClusterUsed(iLayer - 1, neighbourCell.getFirstClusterIndex())) {
606 if (currentCell.getLevel() - 1 != neighbourCell.getLevel()) {
612 const auto& trHit = mTimeFrame->getTrackingFrameInfoOnLayer(iLayer - 1)[neighbourCell.getFirstClusterIndex()];
614 if (!seed.rotate(trHit.alphaTrackingFrame)) {
624 if (mCorrType == o2::base::PropagatorF::MatCorrType::USEMatCorrNONE) {
625 if (!seed.correctForMaterial(mTrkParams[0].LayerxX0[iLayer - 1], mTrkParams[0].LayerxX0[iLayer - 1] *
constants::Radl *
constants::Rho,
true)) {
630 auto predChi2{seed.getPredictedChi2Quiet(trHit.positionTrackingFrame, trHit.covarianceTrackingFrame)};
631 if ((predChi2 > mTrkParams[0].MaxChi2ClusterAttachment) || predChi2 < 0.f) {
635 seed.setChi2(seed.getChi2() + predChi2);
636 if (!seed.o2::track::TrackParCov::update(trHit.positionTrackingFrame, trHit.covarianceTrackingFrame)) {
641 if constexpr (
decltype(Tag)
::value != PassMode::TwoPassCount::value) {
642 seed.getClusters()[iLayer - 1] = neighbourCell.getFirstClusterIndex();
643 seed.setLevel(neighbourCell.getLevel());
644 seed.setFirstTrackletIndex(neighbourCell.getFirstTrackletIndex());
645 seed.setSecondTrackletIndex(neighbourCell.getSecondTrackletIndex());
648 if constexpr (
decltype(Tag)
::value == PassMode::OnePass::value) {
649 updatedCellSeeds.push_back(seed);
650 updatedCellsIds.push_back(neighbourCellId);
651 }
else if constexpr (
decltype(Tag)
::value == PassMode::TwoPassCount::value) {
653 }
else if constexpr (
decltype(Tag)
::value == PassMode::TwoPassInsert::value) {
654 updatedCellSeeds[
offset] = seed;
655 updatedCellsIds[
offset++] = neighbourCellId;
657 static_assert(
false,
"Unknown mode!");
663 const int nCells =
static_cast<int>(currentCellSeed.size());
664 if (mTaskArena->max_concurrency() <= 1) {
665 for (
int iCell{0}; iCell < nCells; ++iCell) {
671 tbb::blocked_range<int>(0, nCells),
672 [&](
const tbb::blocked_range<int>& Cells) {
673 for (
int iCell = Cells.begin(); iCell < Cells.end(); ++iCell) {
678 std::exclusive_scan(perCellCount.begin(), perCellCount.end(), perCellCount.begin(), 0);
679 auto totalNeighbours{perCellCount.back()};
680 if (totalNeighbours == 0) {
683 updatedCellSeeds.resize(totalNeighbours);
684 updatedCellsIds.resize(totalNeighbours);
687 tbb::blocked_range<int>(0, nCells),
688 [&](
const tbb::blocked_range<int>& Cells) {
689 for (
int iCell = Cells.begin(); iCell < Cells.end(); ++iCell) {
690 int offset = perCellCount[iCell];
691 if (
offset == perCellCount[iCell + 1]) {
701 std::cout <<
"\t\t- Found " << updatedCellSeeds.size() <<
" cell seeds out of " << attempts <<
" attempts" << std::endl;
702 std::cout <<
"\t\t\t> " <<
failed[0] <<
" failed because of level" << std::endl;
703 std::cout <<
"\t\t\t> " <<
failed[1] <<
" failed because of rotation" << std::endl;
704 std::cout <<
"\t\t\t> " <<
failed[2] <<
" failed because of propagation" << std::endl;
705 std::cout <<
"\t\t\t> " <<
failed[3] <<
" failed because of chi2 cut" << std::endl;
706 std::cout <<
"\t\t\t> " <<
failed[4] <<
" failed because of update" << std::endl;
707 std::cout <<
"\t\t\t> " << failedByMismatch <<
" failed because of mismatch" << std::endl;
711template <
int nLayers>
714 CA_DEBUGGER(std::cout <<
"Finding roads, iteration " << iteration << std::endl);
716 for (
int startLevel{mTrkParams[iteration].CellsPerRoad()}; startLevel >= mTrkParams[iteration].CellMinimumLevel(); --startLevel) {
717 CA_DEBUGGER(std::cout <<
"\t > Processing level " << startLevel << std::endl);
718 auto seedFilter = [&](
const CellSeed& seed) {
719 return seed.getQ2Pt() <= 1.e3 && seed.getChi2() <= mTrkParams[0].MaxChi2NDF * ((startLevel + 2) * 2 - 5);
722 for (
int startLayer{mTrkParams[iteration].CellsPerRoad() - 1}; startLayer >= startLevel - 1; --startLayer) {
723 if ((mTrkParams[iteration].StartLayerMask & (1 << (startLayer + 2))) == 0) {
726 CA_DEBUGGER(std::cout <<
"\t\t > Starting processing layer " << startLayer << std::endl);
730 processNeighbours(startLayer, startLevel, mTimeFrame->getCells()[startLayer], lastCellId, updatedCellSeed, updatedCellId);
732 int level = startLevel;
733 for (
int iLayer{startLayer - 1}; iLayer > 0 &&
level > 2; --iLayer) {
734 lastCellSeed.swap(updatedCellSeed);
735 lastCellId.swap(updatedCellId);
738 processNeighbours(iLayer, --
level, lastCellSeed, lastCellId, updatedCellSeed, updatedCellId);
743 if (!updatedCellSeed.empty()) {
744 trackSeeds.reserve(trackSeeds.size() + std::count_if(updatedCellSeed.begin(), updatedCellSeed.end(), seedFilter));
745 std::copy_if(updatedCellSeed.begin(), updatedCellSeed.end(), std::back_inserter(trackSeeds), seedFilter);
749 if (trackSeeds.empty()) {
754 mTaskArena->execute([&] {
755 auto forSeed = [&](
auto Tag,
int iSeed,
int offset = 0) {
756 const CellSeed& seed{trackSeeds[iSeed]};
758 temporaryTrack.resetCovariance();
759 temporaryTrack.setChi2(0);
760 for (
int iL{0}; iL < 7; ++iL) {
761 temporaryTrack.setExternalClusterIndex(iL, seed.getCluster(iL), seed.getCluster(iL) !=
constants::UnusedIndex);
764 bool fitSuccess = fitTrack(temporaryTrack, 0, mTrkParams[0].NLayers, 1, mTrkParams[0].MaxChi2ClusterAttachment, mTrkParams[0].MaxChi2NDF);
769 temporaryTrack.getParamOut() = temporaryTrack.getParamIn();
770 temporaryTrack.resetCovariance();
771 temporaryTrack.setChi2(0);
772 fitSuccess = fitTrack(temporaryTrack, mTrkParams[0].NLayers - 1, -1, -1, mTrkParams[0].MaxChi2ClusterAttachment, mTrkParams[0].MaxChi2NDF, 50.f);
773 if (!fitSuccess || temporaryTrack.getPt() < mTrkParams[iteration].MinPt[mTrkParams[iteration].NLayers - temporaryTrack.getNClusters()]) {
777 if constexpr (
decltype(Tag)
::value == PassMode::OnePass::value) {
778 tracks.push_back(temporaryTrack);
779 }
else if constexpr (
decltype(Tag)
::value == PassMode::TwoPassCount::value) {
781 }
else if constexpr (
decltype(Tag)
::value == PassMode::TwoPassInsert::value) {
782 tracks[
offset] = temporaryTrack;
784 static_assert(
false,
"Unknown mode!");
789 const int nSeeds =
static_cast<int>(trackSeeds.size());
790 if (mTaskArena->max_concurrency() <= 1) {
791 for (
int iSeed{0}; iSeed < nSeeds; ++iSeed) {
797 tbb::blocked_range<int>(0, nSeeds),
798 [&](
const tbb::blocked_range<int>& Seeds) {
799 for (
int iSeed = Seeds.begin(); iSeed < Seeds.end(); ++iSeed) {
804 std::exclusive_scan(perSeedCount.begin(), perSeedCount.end(), perSeedCount.begin(), 0);
805 auto totalTracks{perSeedCount.back()};
806 if (totalTracks == 0) {
809 tracks.resize(totalTracks);
812 tbb::blocked_range<int>(0, nSeeds),
813 [&](
const tbb::blocked_range<int>& Seeds) {
814 for (
int iSeed = Seeds.begin(); iSeed < Seeds.end(); ++iSeed) {
815 if (perSeedCount[iSeed] == perSeedCount[iSeed + 1]) {
824 tbb::parallel_sort(tracks.begin(), tracks.end(), [](
const auto&
a,
const auto&
b) {
825 return a.getChi2() < b.getChi2();
829 for (
auto& track : tracks) {
831 bool isFirstShared{
false};
832 for (
int iLayer{0}; iLayer < mTrkParams[0].NLayers; ++iLayer) {
836 nShared +=
int(mTimeFrame->isClusterUsed(iLayer, track.getClusterIndex(iLayer)));
837 isFirstShared |= !iLayer && mTimeFrame->isClusterUsed(iLayer, track.getClusterIndex(iLayer));
840 if (nShared > mTrkParams[0].ClusterSharing) {
844 std::array<int, 3> rofs{INT_MAX, INT_MAX, INT_MAX};
845 for (
int iLayer{0}; iLayer < mTrkParams[0].NLayers; ++iLayer) {
849 mTimeFrame->markUsedCluster(iLayer, track.getClusterIndex(iLayer));
850 int currentROF = mTimeFrame->getClusterROF(iLayer, track.getClusterIndex(iLayer));
851 for (
int iR{0}; iR < 3; ++iR) {
852 if (rofs[iR] == INT_MAX) {
853 rofs[iR] = currentROF;
855 if (rofs[iR] == currentROF) {
860 if (rofs[2] != INT_MAX) {
863 track.setUserField(0);
864 track.getParamOut().setUserField(0);
865 if (rofs[1] != INT_MAX) {
866 track.setNextROFbit();
868 mTimeFrame->getTracks(o2::gpu::CAMath::Min(rofs[0], rofs[1])).emplace_back(track);
873template <
int nLayers>
876 for (
int rof{0}; rof < mTimeFrame->getNrof(); ++rof) {
877 for (
auto& track : mTimeFrame->getTracks(rof)) {
881 if ((mTrkParams[iteration].UseTrackFollowerMix || mTrkParams[iteration].UseTrackFollowerTop) && track.getLastClusterLayer() != mTrkParams[iteration].NLayers - 1) {
882 success = success || trackFollowing(&track, rof,
true, iteration);
884 if ((mTrkParams[iteration].UseTrackFollowerMix || (mTrkParams[iteration].UseTrackFollowerBot && !success)) && track.getFirstClusterLayer() != 0) {
885 success = success || trackFollowing(&track, rof,
false, iteration);
889 track.resetCovariance();
891 bool fitSuccess = fitTrack(track, 0, mTrkParams[iteration].NLayers, 1, mTrkParams[iteration].MaxChi2ClusterAttachment, mTrkParams[0].MaxChi2NDF);
896 track.getParamOut() = track;
897 track.resetCovariance();
899 fitSuccess = fitTrack(track, mTrkParams[iteration].NLayers - 1, -1, -1, mTrkParams[iteration].MaxChi2ClusterAttachment, mTrkParams[0].MaxChi2NDF, 50.);
904 mTimeFrame->mNExtendedTracks++;
905 mTimeFrame->mNExtendedUsedClusters += track.getNClusters() - backup.getNClusters();
906 auto pattern = track.getPattern();
907 auto diff = (
pattern & ~backup.getPattern()) & 0xff;
911 for (
int iLayer{0}; iLayer < mTrkParams[iteration].NLayers; ++iLayer) {
915 mTimeFrame->markUsedCluster(iLayer, track.getClusterIndex(iLayer));
922template <
int nLayers>
926 mTimeFrame->fillPrimaryVerticesXandAlpha();
928 for (
auto& cell : mTimeFrame->getCells()[0]) {
929 auto& cluster3_glo = mTimeFrame->getClusters()[2][cell.getThirdClusterIndex()];
930 auto& cluster2_glo = mTimeFrame->getClusters()[1][cell.getSecondClusterIndex()];
931 auto& cluster1_glo = mTimeFrame->getClusters()[0][cell.getFirstClusterIndex()];
932 if (mTimeFrame->isClusterUsed(2, cluster1_glo.clusterId) ||
933 mTimeFrame->isClusterUsed(1, cluster2_glo.clusterId) ||
934 mTimeFrame->isClusterUsed(0, cluster3_glo.clusterId)) {
938 std::array<int, 3> rofs{
939 mTimeFrame->getClusterROF(2, cluster3_glo.clusterId),
940 mTimeFrame->getClusterROF(1, cluster2_glo.clusterId),
941 mTimeFrame->getClusterROF(0, cluster1_glo.clusterId)};
942 if (rofs[0] != rofs[1] && rofs[1] != rofs[2] && rofs[0] != rofs[2]) {
947 if (rofs[1] == rofs[2]) {
951 auto pvs{mTimeFrame->getPrimaryVertices(rof)};
952 auto pvsXAlpha{mTimeFrame->getPrimaryVerticesXAlpha(rof)};
954 const auto& cluster3_tf = mTimeFrame->getTrackingFrameInfoOnLayer(2)[cluster3_glo.clusterId];
955 TrackITSExt temporaryTrack{buildTrackSeed(cluster1_glo, cluster2_glo, cluster3_tf)};
956 temporaryTrack.setExternalClusterIndex(0, cluster1_glo.clusterId,
true);
957 temporaryTrack.setExternalClusterIndex(1, cluster2_glo.clusterId,
true);
958 temporaryTrack.setExternalClusterIndex(2, cluster3_glo.clusterId,
true);
961 bool fitSuccess = fitTrack(temporaryTrack, 1, -1, -1);
967 TrackITSExt bestTrack{temporaryTrack}, backup{temporaryTrack};
968 float bestChi2{std::numeric_limits<float>::max()};
969 for (
int iV{0}; iV < (
int)pvs.size(); ++iV) {
970 temporaryTrack = backup;
971 if (!temporaryTrack.rotate(pvsXAlpha[iV][1])) {
974 if (!propagator->propagateTo(temporaryTrack, pvsXAlpha[iV][0],
true)) {
978 float pvRes{mTrkParams[0].PVres / o2::gpu::CAMath::Sqrt(
float(pvs[iV].getNContributors()))};
979 const float posVtx[2]{0.f, pvs[iV].getZ()};
980 const float covVtx[3]{pvRes, 0.f, pvRes};
981 float chi2 = temporaryTrack.getPredictedChi2Quiet(posVtx, covVtx);
982 if (chi2 < bestChi2) {
983 if (!temporaryTrack.track::TrackParCov::update(posVtx, covVtx)) {
986 bestTrack = temporaryTrack;
991 bestTrack.resetCovariance();
992 bestTrack.setChi2(0.f);
993 fitSuccess = fitTrack(bestTrack, 0, mTrkParams[0].NLayers, 1, mTrkParams[0].MaxChi2ClusterAttachment, mTrkParams[0].MaxChi2NDF);
997 bestTrack.getParamOut() = bestTrack;
998 bestTrack.resetCovariance();
999 bestTrack.setChi2(0.f);
1000 fitSuccess = fitTrack(bestTrack, mTrkParams[0].NLayers - 1, -1, -1, mTrkParams[0].MaxChi2ClusterAttachment, mTrkParams[0].MaxChi2NDF, 50.);
1004 mTimeFrame->markUsedCluster(0, bestTrack.getClusterIndex(0));
1005 mTimeFrame->markUsedCluster(1, bestTrack.getClusterIndex(1));
1006 mTimeFrame->markUsedCluster(2, bestTrack.getClusterIndex(2));
1007 mTimeFrame->getTracks(rof).emplace_back(bestTrack);
1011template <
int nLayers>
1016 for (
int iLayer{
start}; iLayer !=
end; iLayer += step) {
1020 const TrackingFrameInfo& trackingHit = mTimeFrame->getTrackingFrameInfoOnLayer(iLayer)[track.getClusterIndex(iLayer)];
1022 if (!track.rotate(trackingHit.alphaTrackingFrame)) {
1030 if (mCorrType == o2::base::PropagatorF::MatCorrType::USEMatCorrNONE) {
1036 auto predChi2{track.getPredictedChi2Quiet(trackingHit.positionTrackingFrame, trackingHit.covarianceTrackingFrame)};
1037 if ((nCl >= 3 && predChi2 > chi2clcut) || predChi2 < 0.f) {
1040 track.setChi2(track.getChi2() + predChi2);
1041 if (!track.o2::track::TrackParCov::update(trackingHit.positionTrackingFrame, trackingHit.covarianceTrackingFrame)) {
1046 return std::abs(track.getQ2Pt()) < maxQoverPt && track.getChi2() < chi2ndfcut * (nCl * 2 - 5);
1049template <
int nLayers>
1053 const int step = -1 + outward * 2;
1054 const int end = outward ? mTrkParams[iteration].NLayers - 1 : 0;
1056 for (
size_t iHypo{0}; iHypo < hypotheses.size(); ++iHypo) {
1057 auto hypo{hypotheses[iHypo]};
1058 int iLayer =
static_cast<int>(outward ? hypo.getLastClusterLayer() : hypo.getFirstClusterLayer());
1060 while (iLayer !=
end) {
1062 const float r = mTrkParams[iteration].LayerRadii[iLayer];
1069 auto& hypoParam{outward ? hypo.getParamOut() : hypo.getParamIn()};
1070 if (!propInstance->propagateToX(hypoParam,
x, mTimeFrame->getBz(), PropagatorF::MAX_SIN_PHI,
1071 PropagatorF::MAX_STEP, mTrkParams[iteration].CorrType)) {
1075 if (mTrkParams[iteration].CorrType == PropagatorF::MatCorrType::USEMatCorrNONE) {
1076 if (!hypoParam.correctForMaterial(mTrkParams[iteration].LayerxX0[iLayer], mTrkParams[iteration].LayerxX0[iLayer] *
constants::Radl *
constants::Rho,
true)) {
1082 const float phi{hypoParam.getPhi()};
1083 const float ePhi{o2::gpu::CAMath::Sqrt(hypoParam.getSigmaSnp2() / hypoParam.getCsp2())};
1084 const float z{hypoParam.getZ()};
1085 const float eZ{o2::gpu::CAMath::Sqrt(hypoParam.getSigmaZ2())};
1086 const int4 selectedBinsRect{getBinsRect(iLayer, phi, mTrkParams[iteration].TrackFollowerNSigmaCutPhi * ePhi,
z, mTrkParams[iteration].TrackFollowerNSigmaCutZ * eZ)};
1087 if (selectedBinsRect.x == 0 && selectedBinsRect.y == 0 && selectedBinsRect.z == 0 && selectedBinsRect.w == 0) {
1091 int phiBinsNum{selectedBinsRect.
w - selectedBinsRect.y + 1};
1093 if (phiBinsNum < 0) {
1094 phiBinsNum += mTrkParams[iteration].PhiBins;
1097 gsl::span<const Cluster> layer1 = mTimeFrame->getClustersOnLayer(rof, iLayer);
1098 if (layer1.empty()) {
1103 for (
int iPhiCount = 0; iPhiCount < phiBinsNum; iPhiCount++) {
1104 int iPhiBin = (selectedBinsRect.y + iPhiCount) % mTrkParams[iteration].PhiBins;
1105 const int firstBinIndex{mTimeFrame->mIndexTableUtils.getBinIndex(selectedBinsRect.x, iPhiBin)};
1106 const int maxBinIndex{firstBinIndex + selectedBinsRect.z - selectedBinsRect.x + 1};
1107 const int firstRowClusterIndex = mTimeFrame->getIndexTable(rof, iLayer)[firstBinIndex];
1108 const int maxRowClusterIndex = mTimeFrame->getIndexTable(rof, iLayer)[maxBinIndex];
1110 for (
int iNextCluster{firstRowClusterIndex}; iNextCluster < maxRowClusterIndex; ++iNextCluster) {
1111 if (iNextCluster >= (
int)layer1.size()) {
1114 const Cluster& nextCluster{layer1[iNextCluster]};
1116 if (mTimeFrame->isClusterUsed(iLayer, nextCluster.clusterId)) {
1120 const TrackingFrameInfo& trackingHit = mTimeFrame->getTrackingFrameInfoOnLayer(iLayer)[nextCluster.clusterId];
1122 auto tbupdated{hypo};
1123 auto& tbuParams = outward ? tbupdated.getParamOut() : tbupdated.getParamIn();
1128 if (!propInstance->propagateToX(tbuParams, trackingHit.
xTrackingFrame, mTimeFrame->getBz(),
1129 PropagatorF::MAX_SIN_PHI, PropagatorF::MAX_STEP, PropagatorF::MatCorrType::USEMatCorrNONE)) {
1134 if (predChi2 >= track->getChi2() * mTrkParams[iteration].NSigmaCut) {
1141 tbupdated.setChi2(tbupdated.getChi2() + predChi2);
1142 tbupdated.setExternalClusterIndex(iLayer, nextCluster.clusterId,
true);
1143 hypotheses.emplace_back(tbupdated);
1150 bool swapped{
false};
1151 for (
auto& hypo : hypotheses) {
1152 if (hypo.isBetter(*bestHypo, track->getChi2() * mTrkParams[iteration].NSigmaCut)) {
1163template <
int nLayers>
1166 float ca{-999.f}, sa{-999.f};
1177 float tgp{1.f}, crv{1.f}, snp{-999.f}, tgl12{-999.f}, tgl23{-999.f}, q2pt{1.f /
track::kMostProbablePt}, q2pt2{1.f}, sg2q2pt{-999.f};
1179 tgp = o2::gpu::CAMath::ATan2(y3 -
y1, x3 -
x1);
1180 snp = tgp / o2::gpu::CAMath::Sqrt(1.f + tgp * tgp);
1182 crv = math_utils::computeCurvature(x3, y3, x2, y2,
x1,
y1);
1183 snp = crv * (
x3 - math_utils::computeCurvatureCentreX(x3, y3, x2, y2,
x1,
y1));
1187 tgl12 = math_utils::computeTanDipAngle(
x1,
y1, x2, y2, z1, z2);
1188 tgl23 = math_utils::computeTanDipAngle(x2, y2, x3, y3, z2, z3);
1189 sg2q2pt =
track::kC1Pt2max * (q2pt2 > 0.0005f ? (q2pt2 < 1.f ? q2pt2 : 1.f) : 0.0005f);
1190 return {tf3.
xTrackingFrame, tf3.
alphaTrackingFrame, {
y3, z3, snp, 0.5f * (tgl12 + tgl23), q2pt}, {tf3.
covarianceTrackingFrame[0], tf3.
covarianceTrackingFrame[1], tf3.
covarianceTrackingFrame[2], 0.f, 0.f,
track::kCSnp2max, 0.f, 0.f, 0.f,
track::kCTgl2max, 0.f, 0.f, 0.f, 0.f, sg2q2pt}};
1193template <
int nLayers>
1197 mIsZeroField = std::abs(mBz) < 0.01;
1198 mTimeFrame->setBz(bz);
1201template <
int nLayers>
1207template <
int nLayers>
1210#if defined(OPTIMISATION_OUTPUT) || defined(CA_DEBUG)
1211 mTaskArena = std::make_shared<tbb::task_arena>(1);
1213 if (arena ==
nullptr) {
1214 mTaskArena = std::make_shared<tbb::task_arena>(std::abs(
n));
1215 LOGP(info,
"Setting tracker with {} threads.",
n);
1218 LOGP(info,
"Attaching tracker to calling thread's arena");