diff --git a/DataFormats/Detectors/ITSMFT/ITS/include/DataFormatsITS/TrackITS.h b/DataFormats/Detectors/ITSMFT/ITS/include/DataFormatsITS/TrackITS.h index 20fb7c63ebacd..9b63509cc9424 100644 --- a/DataFormats/Detectors/ITSMFT/ITS/include/DataFormatsITS/TrackITS.h +++ b/DataFormats/Detectors/ITSMFT/ITS/include/DataFormatsITS/TrackITS.h @@ -35,6 +35,11 @@ namespace its class TrackITS : public o2::track::TrackParCov { + public: + static constexpr unsigned int ExtendedPatternShift = 24; + static constexpr int MaxLayersInTrackPattern = 8; + + private: enum UserBits { kSharedClusters = 1 << 28 }; @@ -106,8 +111,39 @@ class TrackITS : public o2::track::TrackParCov GPUhdi() uint32_t getPattern() const { return mPattern; } bool hasHitOnLayer(uint32_t i) const { return mPattern & (0x1 << i); } bool isFakeOnLayer(uint32_t i) const { return !(mPattern & (0x1 << (16 + i))); } - bool isExtendedOnLayer(uint32_t i) const { return (mPattern & (0x1 << (24 + i))); } // only correct if getNClusters <= 8 on layers <= 8 - uint32_t getLastClusterLayer() const + bool isExtendedOnLayer(uint32_t i) const { return (mPattern & (0x1 << (ExtendedPatternShift + i))); } // only correct if getNClusters <= 8 on layers <= 8 + template + GPUhdi() static constexpr uint32_t getLayerPatternMask() + { + return (NLayers >= 32) ? 0xffffffffu : ((1u << NLayers) - 1u); + } + template + GPUhdi() void setExtendedLayerPattern(uint32_t pattern) + { + pattern &= getLayerPatternMask(); + setUserField(static_cast(pattern)); + if constexpr (NLayers <= MaxLayersInTrackPattern) { + setPattern(getPattern() | (pattern << ExtendedPatternShift)); + } + } + template + GPUhdi() uint32_t getExtendedLayerPattern() const + { + const auto mask = getLayerPatternMask(); + if constexpr (NLayers <= MaxLayersInTrackPattern) { + const auto pattern = (getPattern() >> ExtendedPatternShift) & mask; + if (pattern) { + return pattern; + } + } + return getUserField() & mask; + } + GPUhdi() void clearExtendedLayerPattern() + { + setUserField(0); + getParamOut().setUserField(0); + } + GPUhdi() uint32_t getLastClusterLayer() const { uint32_t r{0}, v{mPattern & ((1 << 16) - 1)}; while (v >>= 1) { @@ -115,7 +151,7 @@ class TrackITS : public o2::track::TrackParCov } return r; } - uint32_t getFirstClusterLayer() const + GPUhdi() uint32_t getFirstClusterLayer() const { int s{0}; while (!(mPattern & (1 << s))) { diff --git a/Detectors/ITSMFT/ITS/tracking/GPU/ITStrackingGPU/TimeFrameGPU.h b/Detectors/ITSMFT/ITS/tracking/GPU/ITStrackingGPU/TimeFrameGPU.h index 5f56e3f272473..e2542d841c8bf 100644 --- a/Detectors/ITSMFT/ITS/tracking/GPU/ITStrackingGPU/TimeFrameGPU.h +++ b/Detectors/ITSMFT/ITS/tracking/GPU/ITStrackingGPU/TimeFrameGPU.h @@ -19,6 +19,7 @@ #include "ITStracking/BoundedAllocator.h" #include "ITStracking/TimeFrame.h" #include "ITStracking/Configuration.h" +#include "ITStracking/TrackExtensionHypothesis.h" #include "ITStrackingGPU/Utils.h" namespace o2::its::gpu @@ -90,6 +91,7 @@ class TimeFrameGPU : public TimeFrame void createNeighboursDevice(const unsigned int layer); void createNeighboursLUTDevice(const int, const unsigned int); void createTrackITSExtDevice(const size_t); + void createTrackExtensionScratchDevice(const int nThreads, const int maxHypotheses); void downloadTrackITSExtDevice(); void downloadCellsNeighboursDevice(std::vector>&, const int); void downloadNeighboursLUTDevice(bounded_vector&, const int); @@ -125,6 +127,8 @@ class TimeFrameGPU : public TimeFrame // Hybrid TrackITSExt* getDeviceTrackITSExt() { return mTrackITSExtDevice; } + TrackExtensionHypothesis* getDeviceActiveTrackExtensionHypotheses() { return mActiveTrackExtensionHypothesesDevice; } + TrackExtensionHypothesis* getDeviceNextTrackExtensionHypotheses() { return mNextTrackExtensionHypothesesDevice; } int* getDeviceNeighboursLUT(const int layer) { return mNeighboursLUTDevice[layer]; } gsl::span getDeviceNeighboursLUTs() { return mNeighboursLUTDevice; } CellNeighbour** getDeviceArrayNeighbours() { return mNeighboursDeviceArray; } @@ -222,6 +226,8 @@ class TimeFrameGPU : public TimeFrame float** mCellSeedsChi2DeviceArray; TrackITSExt* mTrackITSExtDevice; + TrackExtensionHypothesis* mActiveTrackExtensionHypothesesDevice{nullptr}; + TrackExtensionHypothesis* mNextTrackExtensionHypothesesDevice{nullptr}; std::array mNeighboursDevice{}; CellNeighbour** mNeighboursDeviceArray{nullptr}; std::array mTrackingFrameInfoDevice; diff --git a/Detectors/ITSMFT/ITS/tracking/GPU/ITStrackingGPU/TrackingKernels.h b/Detectors/ITSMFT/ITS/tracking/GPU/ITStrackingGPU/TrackingKernels.h index 161283db2a2bc..634ac5217a089 100644 --- a/Detectors/ITSMFT/ITS/tracking/GPU/ITStrackingGPU/TrackingKernels.h +++ b/Detectors/ITSMFT/ITS/tracking/GPU/ITStrackingGPU/TrackingKernels.h @@ -13,11 +13,13 @@ #ifndef ITSTRACKINGGPU_TRACKINGKERNELS_H_ #define ITSTRACKINGGPU_TRACKINGKERNELS_H_ +#include #include #include "ITStracking/BoundedAllocator.h" #include "ITStracking/ROFLookupTables.h" #include "ITStracking/TrackingTopology.h" +#include "ITStracking/TrackExtensionHypothesis.h" #include "ITStrackingGPU/Utils.h" #include "DetectorsBase/Propagator.h" @@ -208,7 +210,6 @@ void countTrackSeedHandler(TrackSeed* trackSeeds, const std::vector& layerxX0Host, const unsigned int nSeeds, const float Bz, - const int startLevel, const float maxChi2ClusterAttachment, const float maxChi2NDF, const int reseedIfShorter, @@ -222,20 +223,35 @@ template void computeTrackSeedHandler(TrackSeed* trackSeeds, const TrackingFrameInfo** foundTrackingFrameInfo, const Cluster** unsortedClusters, + const IndexTableUtils* utils, + const typename ROFMaskTable::View& rofMask, + const typename ROFOverlapTable::View& rofOverlaps, + const Cluster** clusters, + const unsigned char** usedClusters, + const int** clustersIndexTables, + const int** ROFClusters, o2::its::TrackITSExt* tracks, const int* seedLUT, + TrackExtensionHypothesis* activeHypotheses, + TrackExtensionHypothesis* nextHypotheses, const std::vector& layerRadiiHost, const std::vector& minPtsHost, const std::vector& layerxX0Host, const unsigned int nSeeds, const unsigned int nTracks, const float Bz, - const int startLevel, const float maxChi2ClusterAttachment, const float maxChi2NDF, const int reseedIfShorter, const bool repeatRefitOut, const bool shiftRefToCluster, + const int nLayers, + const int phiBins, + const int maxHypotheses, + const bool extendTop, + const bool extendBot, + const float nSigmaCutPhi, + const float nSigmaCutZ, const o2::base::Propagator* propagator, const o2::base::PropagatorF::MatCorrType matCorrType, o2::its::ExternalAllocator* alloc); diff --git a/Detectors/ITSMFT/ITS/tracking/GPU/cuda/TimeFrameGPU.cu b/Detectors/ITSMFT/ITS/tracking/GPU/cuda/TimeFrameGPU.cu index 5fff30f5162b1..bc2d90165ee96 100644 --- a/Detectors/ITSMFT/ITS/tracking/GPU/cuda/TimeFrameGPU.cu +++ b/Detectors/ITSMFT/ITS/tracking/GPU/cuda/TimeFrameGPU.cu @@ -12,6 +12,7 @@ #include +#include #include #include @@ -581,6 +582,21 @@ void TimeFrameGPU::createTrackITSExtDevice(const size_t nSeeds) GPUChkErrS(cudaMemset(mTrackITSExtDevice, 0, mNTracks * sizeof(o2::its::TrackITSExt))); } +template +void TimeFrameGPU::createTrackExtensionScratchDevice(const int nThreads, const int maxHypotheses) +{ + GPUTimer timer("reserving track extension scratch"); + const size_t nHypotheses = static_cast(std::max(1, nThreads)) * std::max(1, maxHypotheses); + GPULog("gpu-allocation: reserving {} track extension hypotheses per scratch buffer, for {:.2f} MB each.", nHypotheses, nHypotheses * sizeof(o2::its::TrackExtensionHypothesis) / constants::MB); + mActiveTrackExtensionHypothesesDevice = nullptr; + mNextTrackExtensionHypothesesDevice = nullptr; + if (nHypotheses == 0) { + return; + } + allocMem(reinterpret_cast(&mActiveTrackExtensionHypothesesDevice), nHypotheses * sizeof(o2::its::TrackExtensionHypothesis), this->hasFrameworkAllocator(), (o2::gpu::GPUMemoryResource::MEMORY_GPU | o2::gpu::GPUMemoryResource::MEMORY_STACK)); + allocMem(reinterpret_cast(&mNextTrackExtensionHypothesesDevice), nHypotheses * sizeof(o2::its::TrackExtensionHypothesis), this->hasFrameworkAllocator(), (o2::gpu::GPUMemoryResource::MEMORY_GPU | o2::gpu::GPUMemoryResource::MEMORY_STACK)); +} + template void TimeFrameGPU::downloadCellsDevice() { diff --git a/Detectors/ITSMFT/ITS/tracking/GPU/cuda/TrackerTraitsGPU.cxx b/Detectors/ITSMFT/ITS/tracking/GPU/cuda/TrackerTraitsGPU.cxx index 141d558712e6d..7ee18d69ceee0 100644 --- a/Detectors/ITSMFT/ITS/tracking/GPU/cuda/TrackerTraitsGPU.cxx +++ b/Detectors/ITSMFT/ITS/tracking/GPU/cuda/TrackerTraitsGPU.cxx @@ -12,12 +12,19 @@ #include +#include +#include + #include "ITStrackingGPU/TrackerTraitsGPU.h" #include "ITStrackingGPU/TrackingKernels.h" #include "ITStracking/Configuration.h" namespace o2::its { +namespace +{ +constexpr int trackExtensionLaunchThreads = 60 * 256; +} template void TrackerTraitsGPU::initialiseTimeFrame(const int iteration) @@ -301,10 +308,11 @@ template void TrackerTraitsGPU::findRoads(const int iteration) { bounded_vector> firstClusters(this->mTrkParams[iteration].NLayers, bounded_vector(this->getMemoryPool().get()), this->getMemoryPool().get()); - bounded_vector> sharedFirstClusters(this->mTrkParams[iteration].NLayers, bounded_vector(this->getMemoryPool().get()), this->getMemoryPool().get()); firstClusters.resize(this->mTrkParams[iteration].NLayers); - sharedFirstClusters.resize(this->mTrkParams[iteration].NLayers); const auto hostTopology = mTimeFrameGPU->getTrackingTopologyView(); + const bool extendTop = this->mTrkParams[iteration].PassFlags[IterationStep::TrackFollowerTop]; + const bool extendBot = this->mTrkParams[iteration].PassFlags[IterationStep::TrackFollowerBot]; + const bool extendTracks = extendTop || extendBot; for (int startLevel{this->mTrkParams[iteration].CellsPerRoad()}; startLevel >= this->mTrkParams[iteration].CellMinimumLevel(); --startLevel) { bounded_vector> trackSeeds(this->getMemoryPool().get()); for (int startCellTopologyId{0}; startCellTopologyId < hostTopology.nCells; ++startCellTopologyId) { @@ -353,7 +361,6 @@ void TrackerTraitsGPU::findRoads(const int iteration) this->mTrkParams[iteration].LayerxX0, trackSeeds.size(), this->mBz, - startLevel, this->mTrkParams[iteration].MaxChi2ClusterAttachment, this->mTrkParams[iteration].MaxChi2NDF, this->mTrkParams[iteration].ReseedIfShorter, @@ -363,23 +370,41 @@ void TrackerTraitsGPU::findRoads(const int iteration) this->mTrkParams[iteration].CorrType, mTimeFrameGPU->getFrameworkAllocator()); mTimeFrameGPU->createTrackITSExtDevice(trackSeeds.size()); + if (extendTracks) { + mTimeFrameGPU->createTrackExtensionScratchDevice(trackExtensionLaunchThreads, this->mTrkParams[iteration].TrackFollowerMaxHypotheses); + } computeTrackSeedHandler(mTimeFrameGPU->getDeviceTrackSeeds(), mTimeFrameGPU->getDeviceArrayTrackingFrameInfo(), mTimeFrameGPU->getDeviceArrayUnsortedClusters(), + mTimeFrameGPU->getDeviceIndexTableUtils(), + mTimeFrameGPU->getDeviceROFMaskTableView(), + mTimeFrameGPU->getDeviceROFOverlapTableView(), + mTimeFrameGPU->getDeviceArrayClusters(), + (const unsigned char**)mTimeFrameGPU->getDeviceArrayUsedClusters(), + mTimeFrameGPU->getDeviceArrayClustersIndexTables(), + mTimeFrameGPU->getDeviceROFrameClusters(), mTimeFrameGPU->getDeviceTrackITSExt(), mTimeFrameGPU->getDeviceTrackSeedsLUT(), + extendTracks ? mTimeFrameGPU->getDeviceActiveTrackExtensionHypotheses() : nullptr, + extendTracks ? mTimeFrameGPU->getDeviceNextTrackExtensionHypotheses() : nullptr, this->mTrkParams[iteration].LayerRadii, this->mTrkParams[iteration].MinPt, this->mTrkParams[iteration].LayerxX0, trackSeeds.size(), mTimeFrameGPU->getNTrackSeeds(), this->mBz, - startLevel, this->mTrkParams[iteration].MaxChi2ClusterAttachment, this->mTrkParams[iteration].MaxChi2NDF, this->mTrkParams[iteration].ReseedIfShorter, this->mTrkParams[iteration].RepeatRefitOut, this->mTrkParams[iteration].ShiftRefToCluster, + this->mTrkParams[iteration].NLayers, + this->mTrkParams[iteration].PhiBins, + this->mTrkParams[iteration].TrackFollowerMaxHypotheses, + extendTop, + extendBot, + this->mTrkParams[iteration].TrackFollowerNSigmaCutPhi, + this->mTrkParams[iteration].TrackFollowerNSigmaCutZ, mTimeFrameGPU->getDevicePropagator(), this->mTrkParams[iteration].CorrType, mTimeFrameGPU->getFrameworkAllocator()); diff --git a/Detectors/ITSMFT/ITS/tracking/GPU/cuda/TrackingKernels.cu b/Detectors/ITSMFT/ITS/tracking/GPU/cuda/TrackingKernels.cu index 571afe08fc209..6d171da82f274 100644 --- a/Detectors/ITSMFT/ITS/tracking/GPU/cuda/TrackingKernels.cu +++ b/Detectors/ITSMFT/ITS/tracking/GPU/cuda/TrackingKernels.cu @@ -20,6 +20,7 @@ #include #include #include +#include #include #include @@ -31,6 +32,7 @@ #include "ITStracking/Tracklet.h" #include "ITStracking/Cluster.h" #include "ITStracking/Cell.h" +#include "ITStracking/TrackFollower.h" #include "ITStracking/TrackHelpers.h" #include "DataFormatsITS/TrackITS.h" #include "ITStrackingGPU/TrackingKernels.h" @@ -108,35 +110,115 @@ struct compare_track_chi2 { } }; -template +template +GPUdi() void finaliseTrackExtensionTrial(const uint32_t backupPattern, + TrackITSInternal& trial, + const TrackingFrameInfo* const* trackingFrameInfo, + const float* layerxX0, + const int nLayers, + const float bz, + const float maxChi2ClusterAttachment, + const float maxChi2NDF, + const o2::base::Propagator* propagator, + const o2::base::PropagatorF::MatCorrType matCorrType, + const bool shiftRefToCluster, + const bool repeatRefitOut, + TrackITSInternal& best, + uint32_t& bestDiff) +{ + const auto diff = (trial.getPattern() & ~backupPattern) & TrackITS::getLayerPatternMask(); + if (!diff || + !o2::its::track::refitTrack(trial, trackingFrameInfo, layerxX0, nLayers, bz, maxChi2ClusterAttachment, maxChi2NDF, propagator, matCorrType, shiftRefToCluster, repeatRefitOut)) { + return; + } + if (o2::its::track::isBetter(trial, best)) { + best = trial; + bestDiff = diff; + } +} + +template +GPUg() void __launch_bounds__(256, 1) countTrackSeedsKernel( + TrackSeed* trackSeeds, + const TrackingFrameInfo** foundTrackingFrameInfo, + const Cluster** unsortedClusters, + int* seedLUT, + const float* layerRadii, + const float* minPts, + const float* layerxX0, + const unsigned int nSeeds, + const float bz, + const float maxChi2ClusterAttachment, + const float maxChi2NDF, + const int reseedIfShorter, + const bool repeatRefitOut, + const bool shiftRefToCluster, + const o2::base::Propagator* propagator, + const o2::base::PropagatorF::MatCorrType matCorrType) +{ + for (int iCurrentTrackSeedIndex = blockIdx.x * blockDim.x + threadIdx.x; iCurrentTrackSeedIndex < nSeeds; iCurrentTrackSeedIndex += blockDim.x * gridDim.x) { + TrackITSInternal temporaryTrack; + if (o2::its::track::refitTrack(trackSeeds[iCurrentTrackSeedIndex], + temporaryTrack, + maxChi2ClusterAttachment, + maxChi2NDF, + bz, + foundTrackingFrameInfo, + unsortedClusters, + layerxX0, + layerRadii, + minPts, + propagator, + matCorrType, + reseedIfShorter, + shiftRefToCluster, + repeatRefitOut)) { + seedLUT[iCurrentTrackSeedIndex] = 1; + } + } +} + +template GPUg() void __launch_bounds__(256, 1) fitTrackSeedsKernel( TrackSeed* trackSeeds, const TrackingFrameInfo** foundTrackingFrameInfo, const Cluster** unsortedClusters, + const IndexTableUtils* utils, + const typename ROFMaskTable::View rofMask, + const typename ROFOverlapTable::View rofOverlaps, + const Cluster** clusters, + const unsigned char** usedClusters, + const int** clustersIndexTables, + const int** ROFClusters, o2::its::TrackITSExt* tracks, - maybe_const* seedLUT, + const int* seedLUT, + TrackExtensionHypothesis* activeHypothesesScratch, + TrackExtensionHypothesis* nextHypothesesScratch, const float* layerRadii, const float* minPts, const float* layerxX0, const unsigned int nSeeds, const float bz, - const int startLevel, const float maxChi2ClusterAttachment, const float maxChi2NDF, const int reseedIfShorter, const bool repeatRefitOut, const bool shiftRefToCluster, + const int nLayers, + const int phiBins, + const int maxHypothesesConfig, + const bool extendTop, + const bool extendBot, + const float nSigmaCutPhi, + const float nSigmaCutZ, const o2::base::Propagator* propagator, const o2::base::PropagatorF::MatCorrType matCorrType) { for (int iCurrentTrackSeedIndex = blockIdx.x * blockDim.x + threadIdx.x; iCurrentTrackSeedIndex < nSeeds; iCurrentTrackSeedIndex += blockDim.x * gridDim.x) { - - if constexpr (!initRun) { - if (seedLUT[iCurrentTrackSeedIndex] == seedLUT[iCurrentTrackSeedIndex + 1]) { - continue; - } + if (seedLUT[iCurrentTrackSeedIndex] == seedLUT[iCurrentTrackSeedIndex + 1]) { + continue; } - TrackITSExt temporaryTrack; + TrackITSInternal temporaryTrack; bool refitSuccess = o2::its::track::refitTrack(trackSeeds[iCurrentTrackSeedIndex], temporaryTrack, maxChi2ClusterAttachment, @@ -153,11 +235,166 @@ GPUg() void __launch_bounds__(256, 1) fitTrackSeedsKernel( shiftRefToCluster, repeatRefitOut); if (refitSuccess) { - if constexpr (initRun) { - seedLUT[iCurrentTrackSeedIndex] = 1; - } else { - tracks[seedLUT[iCurrentTrackSeedIndex]] = temporaryTrack; + if ((extendTop || extendBot) && activeHypothesesScratch && nextHypothesesScratch) { + const int maxHypotheses = o2::gpu::CAMath::Max(maxHypothesesConfig, 1); + const int threadIndex = blockIdx.x * blockDim.x + threadIdx.x; + auto* activeHypotheses = activeHypothesesScratch + threadIndex * maxHypotheses; + auto* nextHypotheses = nextHypothesesScratch + threadIndex * maxHypotheses; + const auto backupPattern = temporaryTrack.getPattern(); + auto best = temporaryTrack; + uint32_t bestDiff{0}; + TrackITSInternal topResult; + TrackITSInternal botResult; + bool hasTopResult{false}; + bool hasBotResult{false}; + const uint32_t lastLayer = static_cast(nLayers - 1); + + if (extendTop && temporaryTrack.getLastClusterLayer() != lastLayer) { + auto candidate = temporaryTrack; + const auto startHypothesis = TrackExtensionHypothesis{temporaryTrack, true}; + TrackExtensionHypothesis bestHypothesis; + if (followTrackExtensionDirection(startHypothesis, + *utils, + rofMask, + rofOverlaps, + clusters, + usedClusters, + clustersIndexTables, + ROFClusters, + foundTrackingFrameInfo, + layerRadii, + layerxX0, + nLayers, + phiBins, + maxHypotheses, + bz, + maxChi2ClusterAttachment, + maxChi2NDF, + nSigmaCutPhi, + nSigmaCutZ, + true, + propagator, + matCorrType, + activeHypotheses, + nextHypotheses, + bestHypothesis)) { + updateTrackFromExtensionHypothesis(bestHypothesis, true, nLayers, candidate); + topResult = candidate; + hasTopResult = true; + finaliseTrackExtensionTrial(backupPattern, candidate, foundTrackingFrameInfo, layerxX0, nLayers, bz, maxChi2ClusterAttachment, maxChi2NDF, propagator, matCorrType, shiftRefToCluster, repeatRefitOut, best, bestDiff); + } + } + if (extendBot && temporaryTrack.getFirstClusterLayer() != 0) { + auto candidate = temporaryTrack; + const auto startHypothesis = TrackExtensionHypothesis{temporaryTrack, false}; + TrackExtensionHypothesis bestHypothesis; + if (followTrackExtensionDirection(startHypothesis, + *utils, + rofMask, + rofOverlaps, + clusters, + usedClusters, + clustersIndexTables, + ROFClusters, + foundTrackingFrameInfo, + layerRadii, + layerxX0, + nLayers, + phiBins, + maxHypotheses, + bz, + maxChi2ClusterAttachment, + maxChi2NDF, + nSigmaCutPhi, + nSigmaCutZ, + false, + propagator, + matCorrType, + activeHypotheses, + nextHypotheses, + bestHypothesis)) { + updateTrackFromExtensionHypothesis(bestHypothesis, false, nLayers, candidate); + botResult = candidate; + hasBotResult = true; + finaliseTrackExtensionTrial(backupPattern, candidate, foundTrackingFrameInfo, layerxX0, nLayers, bz, maxChi2ClusterAttachment, maxChi2NDF, propagator, matCorrType, shiftRefToCluster, repeatRefitOut, best, bestDiff); + } + } + if (extendTop && extendBot) { + if (hasTopResult && topResult.getFirstClusterLayer() != 0) { + auto candidate = topResult; + const auto startHypothesis = TrackExtensionHypothesis{topResult, false}; + TrackExtensionHypothesis bestHypothesis; + if (followTrackExtensionDirection(startHypothesis, + *utils, + rofMask, + rofOverlaps, + clusters, + usedClusters, + clustersIndexTables, + ROFClusters, + foundTrackingFrameInfo, + layerRadii, + layerxX0, + nLayers, + phiBins, + maxHypotheses, + bz, + maxChi2ClusterAttachment, + maxChi2NDF, + nSigmaCutPhi, + nSigmaCutZ, + false, + propagator, + matCorrType, + activeHypotheses, + nextHypotheses, + bestHypothesis)) { + updateTrackFromExtensionHypothesis(bestHypothesis, false, nLayers, candidate); + finaliseTrackExtensionTrial(backupPattern, candidate, foundTrackingFrameInfo, layerxX0, nLayers, bz, maxChi2ClusterAttachment, maxChi2NDF, propagator, matCorrType, shiftRefToCluster, repeatRefitOut, best, bestDiff); + } + } + if (hasBotResult && botResult.getLastClusterLayer() != lastLayer) { + auto candidate = botResult; + const auto startHypothesis = TrackExtensionHypothesis{botResult, true}; + TrackExtensionHypothesis bestHypothesis; + if (followTrackExtensionDirection(startHypothesis, + *utils, + rofMask, + rofOverlaps, + clusters, + usedClusters, + clustersIndexTables, + ROFClusters, + foundTrackingFrameInfo, + layerRadii, + layerxX0, + nLayers, + phiBins, + maxHypotheses, + bz, + maxChi2ClusterAttachment, + maxChi2NDF, + nSigmaCutPhi, + nSigmaCutZ, + true, + propagator, + matCorrType, + activeHypotheses, + nextHypotheses, + bestHypothesis)) { + updateTrackFromExtensionHypothesis(bestHypothesis, true, nLayers, candidate); + finaliseTrackExtensionTrial(backupPattern, candidate, foundTrackingFrameInfo, layerxX0, nLayers, bz, maxChi2ClusterAttachment, maxChi2NDF, propagator, matCorrType, shiftRefToCluster, repeatRefitOut, best, bestDiff); + } + } + } + temporaryTrack = best; + tracks[seedLUT[iCurrentTrackSeedIndex]] = makeTrackITSExt(temporaryTrack); + if (bestDiff) { + tracks[seedLUT[iCurrentTrackSeedIndex]].setExtendedLayerPattern(bestDiff); + } + continue; } + tracks[seedLUT[iCurrentTrackSeedIndex]] = makeTrackITSExt(temporaryTrack); } } } @@ -1042,7 +1279,6 @@ void countTrackSeedHandler(TrackSeed* trackSeeds, const std::vector& layerxX0Host, const unsigned int nSeeds, const float bz, - const int startLevel, const float maxChi2ClusterAttachment, const float maxChi2NDF, const int reseedIfShorter, @@ -1058,18 +1294,16 @@ void countTrackSeedHandler(TrackSeed* trackSeeds, thrust::device_vector minPts(minPtsHost); thrust::device_vector layerRadii(layerRadiiHost); thrust::device_vector layerxX0(layerxX0Host); - gpu::fitTrackSeedsKernel<<<60, 256>>>( + gpu::countTrackSeedsKernel<<<60, 256>>>( trackSeeds, // CellSeed* foundTrackingFrameInfo, // TrackingFrameInfo** unsortedClusters, // Cluster** - nullptr, // TrackITSExt* seedLUT, // int* thrust::raw_pointer_cast(&layerRadii[0]), // const float* thrust::raw_pointer_cast(&minPts[0]), // const float* thrust::raw_pointer_cast(&layerxX0[0]), // const float* nSeeds, // const unsigned int bz, // const float - startLevel, // const int maxChi2ClusterAttachment, // float maxChi2NDF, // float reseedIfShorter, // int @@ -1085,20 +1319,35 @@ template void computeTrackSeedHandler(TrackSeed* trackSeeds, const TrackingFrameInfo** foundTrackingFrameInfo, const Cluster** unsortedClusters, + const IndexTableUtils* utils, + const typename ROFMaskTable::View& rofMask, + const typename ROFOverlapTable::View& rofOverlaps, + const Cluster** clusters, + const unsigned char** usedClusters, + const int** clustersIndexTables, + const int** ROFClusters, o2::its::TrackITSExt* tracks, const int* seedLUT, + TrackExtensionHypothesis* activeHypotheses, + TrackExtensionHypothesis* nextHypotheses, const std::vector& layerRadiiHost, const std::vector& minPtsHost, const std::vector& layerxX0Host, const unsigned int nSeeds, const unsigned int nTracks, const float bz, - const int startLevel, const float maxChi2ClusterAttachment, const float maxChi2NDF, const int reseedIfShorter, const bool repeatRefitOut, const bool shiftRefToCluster, + const int nLayers, + const int phiBins, + const int maxHypotheses, + const bool extendTop, + const bool extendBot, + const float nSigmaCutPhi, + const float nSigmaCutZ, const o2::base::Propagator* propagator, const o2::base::PropagatorF::MatCorrType matCorrType, o2::its::ExternalAllocator* alloc) @@ -1106,23 +1355,38 @@ void computeTrackSeedHandler(TrackSeed* trackSeeds, thrust::device_vector minPts(minPtsHost); thrust::device_vector layerRadii(layerRadiiHost); thrust::device_vector layerxX0(layerxX0Host); - gpu::fitTrackSeedsKernel<<<60, 256>>>( + gpu::fitTrackSeedsKernel<<<60, 256>>>( trackSeeds, // CellSeed* foundTrackingFrameInfo, // TrackingFrameInfo** unsortedClusters, // Cluster** + utils, // IndexTableUtils* + rofMask, // ROFMaskTable::View + rofOverlaps, // ROFOverlapTable::View + clusters, // Cluster** + usedClusters, // unsigned char** + clustersIndexTables, // int** + ROFClusters, // int** tracks, // TrackITSExt* seedLUT, // const int* + activeHypotheses, // TrackExtensionHypothesis* + nextHypotheses, // TrackExtensionHypothesis* thrust::raw_pointer_cast(&layerRadii[0]), // const float* thrust::raw_pointer_cast(&minPts[0]), // const float* thrust::raw_pointer_cast(&layerxX0[0]), // const float* nSeeds, // const unsigned int bz, // const float - startLevel, // const int maxChi2ClusterAttachment, // float maxChi2NDF, // float reseedIfShorter, // int repeatRefitOut, // bool shiftRefToCluster, // bool + nLayers, // int + phiBins, // int + maxHypotheses, // int + extendTop, // bool + extendBot, // bool + nSigmaCutPhi, // float + nSigmaCutZ, // float propagator, // const o2::base::Propagator* matCorrType); // o2::base::PropagatorF::MatCorrType auto sync_policy = THRUST_NAMESPACE::par(gpu::TypedAllocator(alloc)); @@ -1284,7 +1548,6 @@ template void countTrackSeedHandler(TrackSeed<7>* trackSeeds, const std::vector& layerxX0Host, const unsigned int nSeeds, const float bz, - const int startLevel, const float maxChi2ClusterAttachment, const float maxChi2NDF, const int reseedIfShorter, @@ -1297,20 +1560,35 @@ template void countTrackSeedHandler(TrackSeed<7>* trackSeeds, template void computeTrackSeedHandler(TrackSeed<7>* trackSeeds, const TrackingFrameInfo** foundTrackingFrameInfo, const Cluster** unsortedClusters, + const IndexTableUtils<7>* utils, + const ROFMaskTable<7>::View& rofMask, + const ROFOverlapTable<7>::View& rofOverlaps, + const Cluster** clusters, + const unsigned char** usedClusters, + const int** clustersIndexTables, + const int** ROFClusters, o2::its::TrackITSExt* tracks, const int* seedLUT, + TrackExtensionHypothesis<7>* activeHypotheses, + TrackExtensionHypothesis<7>* nextHypotheses, const std::vector& layerRadiiHost, const std::vector& minPtsHost, const std::vector& layerxX0Host, const unsigned int nSeeds, const unsigned int nTracks, const float bz, - const int startLevel, const float maxChi2ClusterAttachment, const float maxChi2NDF, const int reseedIfShorter, const bool repeatRefitOut, const bool shiftRefToCluster, + const int nLayers, + const int phiBins, + const int maxHypotheses, + const bool extendTop, + const bool extendBot, + const float nSigmaCutPhi, + const float nSigmaCutZ, const o2::base::Propagator* propagator, const o2::base::PropagatorF::MatCorrType matCorrType, o2::its::ExternalAllocator* alloc); @@ -1470,7 +1748,6 @@ template void countTrackSeedHandler(TrackSeed<11>* trackSeeds, const std::vector& layerxX0Host, const unsigned int nSeeds, const float bz, - const int startLevel, const float maxChi2ClusterAttachment, const float maxChi2NDF, const int reseedIfShorter, @@ -1483,20 +1760,35 @@ template void countTrackSeedHandler(TrackSeed<11>* trackSeeds, template void computeTrackSeedHandler(TrackSeed<11>* trackSeeds, const TrackingFrameInfo** foundTrackingFrameInfo, const Cluster** unsortedClusters, + const IndexTableUtils<11>* utils, + const ROFMaskTable<11>::View& rofMask, + const ROFOverlapTable<11>::View& rofOverlaps, + const Cluster** clusters, + const unsigned char** usedClusters, + const int** clustersIndexTables, + const int** ROFClusters, o2::its::TrackITSExt* tracks, const int* seedLUT, + TrackExtensionHypothesis<11>* activeHypotheses, + TrackExtensionHypothesis<11>* nextHypotheses, const std::vector& layerRadiiHost, const std::vector& minPtsHost, const std::vector& layerxX0Host, const unsigned int nSeeds, const unsigned int nTracks, const float bz, - const int startLevel, const float maxChi2ClusterAttachment, const float maxChi2NDF, const int reseedIfShorter, const bool repeatRefitOut, const bool shiftRefToCluster, + const int nLayers, + const int phiBins, + const int maxHypotheses, + const bool extendTop, + const bool extendBot, + const float nSigmaCutPhi, + const float nSigmaCutZ, const o2::base::Propagator* propagator, const o2::base::PropagatorF::MatCorrType matCorrType, o2::its::ExternalAllocator* alloc); diff --git a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/Configuration.h b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/Configuration.h index 275752854665b..ee8baf20e66c9 100644 --- a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/Configuration.h +++ b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/Configuration.h @@ -32,7 +32,7 @@ namespace o2::its { // Steering of dedicated steps in an iteration -enum class IterationStep : uint8_t { +enum class IterationStep : uint16_t { FirstPass = 0, RebuildClusterLUT, UseUPCMask, @@ -40,6 +40,8 @@ enum class IterationStep : uint8_t { ResetVertices, SkipROFsAboveThreshold, MarkVerticesAsUPC, + TrackFollowerTop, + TrackFollowerBot, }; using IterationSteps = o2::utils::EnumFlags; @@ -94,6 +96,9 @@ struct TrackingParameters { bool DoUPCIteration = false; bool FataliseUponFailure = true; bool CreateArtefactLabels{false}; + float TrackFollowerNSigmaCutZ = 1.f; + float TrackFollowerNSigmaCutPhi = 1.f; + int TrackFollowerMaxHypotheses = 1; bool PrintMemory = false; // print allocator usage in epilog report size_t MaxMemory = std::numeric_limits::max(); bool DropTFUponFailure = false; diff --git a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/IndexTableUtils.h b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/IndexTableUtils.h index 4e8d5bcfea42a..a7b44e91d9093 100644 --- a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/IndexTableUtils.h +++ b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/IndexTableUtils.h @@ -113,14 +113,17 @@ GPUhdi() void IndexTableUtils::print() const } template -GPUhdi() int4 getBinsRect(const Cluster& currentCluster, const int layerIndex, - const float z1, const float z2, const float maxdeltaz, const float maxdeltaphi, +GPUhdi() int4 getBinsRect(const int layerIndex, + const float phi, + const float z, + const float maxdeltaz, + const float maxdeltaphi, const IndexTableUtils& utils) { - const float zRangeMin = o2::gpu::GPUCommonMath::Min(z1, z2) - maxdeltaz; - const float phiRangeMin = (maxdeltaphi > o2::constants::math::PI) ? 0.f : currentCluster.phi - maxdeltaphi; - const float zRangeMax = o2::gpu::GPUCommonMath::Max(z1, z2) + maxdeltaz; - const float phiRangeMax = (maxdeltaphi > o2::constants::math::PI) ? o2::constants::math::TwoPI : currentCluster.phi + maxdeltaphi; + const float zRangeMin = z - maxdeltaz; + const float phiRangeMin = (maxdeltaphi > o2::constants::math::PI) ? 0.f : phi - maxdeltaphi; + const float zRangeMax = z + maxdeltaz; + const float phiRangeMax = (maxdeltaphi > o2::constants::math::PI) ? o2::constants::math::TwoPI : phi + maxdeltaphi; if (zRangeMax < -utils.getLayerZ(layerIndex) || zRangeMin > utils.getLayerZ(layerIndex) || zRangeMin > zRangeMax) { @@ -133,5 +136,15 @@ GPUhdi() int4 getBinsRect(const Cluster& currentCluster, const int layerIndex, utils.getPhiBinIndex(math_utils::getNormalizedPhi(phiRangeMax))}; } +template +GPUhdi() int4 getBinsRect(const Cluster& currentCluster, const int layerIndex, + const float z1, const float z2, const float maxdeltaz, const float maxdeltaphi, + const IndexTableUtils& utils) +{ + const float zMean = 0.5f * (z1 + z2); + const float zDelta = 0.5f * o2::gpu::GPUCommonMath::Abs(z1 - z2) + maxdeltaz; + return getBinsRect(layerIndex, currentCluster.phi, zMean, zDelta, maxdeltaphi, utils); +} + } // namespace o2::its #endif /* TRACKINGITSU_INCLUDE_INDEXTABLEUTILS_H_ */ diff --git a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/ROFLookupTables.h b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/ROFLookupTables.h index ce20169e36c64..1584b5dcd79aa 100644 --- a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/ROFLookupTables.h +++ b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/ROFLookupTables.h @@ -37,6 +37,7 @@ namespace o2::its // Layer timing definition struct LayerTiming { using BCType = TimeStampType; + using BCRange = dataformats::RangeReference; BCType mNROFsTF{0}; // number of ROFs per timeframe BCType mROFLength{0}; // ROF length in BC BCType mROFDelay{0}; // delay of ROFs wrt start of first orbit in TF in BC @@ -73,7 +74,7 @@ struct LayerTiming { } // return which ROF this BC belongs to - GPUhi() BCType getROF(BCType bc) const noexcept + GPUhdi() BCType getROF(BCType bc) const noexcept { const BCType offset = mROFDelay + mROFBias; if (bc <= offset) { @@ -83,7 +84,7 @@ struct LayerTiming { } // return which ROF this timestamp belongs by its lower edge - GPUhi() BCType getROF(TimeStamp ts) const noexcept + GPUhdi() BCType getROF(TimeStamp ts) const noexcept { const BCType offset = mROFDelay + mROFBias; const BCType bc = (ts.getTimeStamp() < ts.getTimeStampError()) ? BCType(0) : static_cast(o2::gpu::CAMath::Floor(ts.getTimeStamp() - ts.getTimeStampError())); @@ -93,6 +94,50 @@ struct LayerTiming { return (bc - offset) / mROFLength; } + // return which ROF this floating point (number of BCs) time belongs + GPUhdi() BCType getROF(float time) const noexcept + { + const float offset = static_cast(mROFDelay + mROFBias); + if (time <= offset) { + return 0; + } + return static_cast((time - offset) / mROFLength); + } + + GPUhdi() bool intersectROF(BCType rof, float lower, float upper) const noexcept + { + const auto rofTS = getROFTimeBounds(rof, true); + return static_cast(rofTS.upper()) > lower && upper > static_cast(rofTS.lower()); + } + + // return clamped ROF range with strictly positive overlap with timestamp interval + GPUhdi() BCRange getROFRange(TimeStamp ts) const noexcept + { + const float lower = ts.getTimeStamp() - ts.getTimeStampError(); + const float upper = ts.getTimeStamp() + ts.getTimeStampError(); + return getROFRange(lower, upper); + } + + GPUhdi() BCRange getROFRange(TimeEstBC ts) const noexcept + { + return getROFRange(static_cast(ts.lower()), static_cast(ts.upper())); + } + + GPUhdi() BCRange getROFRange(float lower, float upper) const noexcept + { + const BCType maxROF = mNROFsTF - 1; + BCType first = o2::gpu::CAMath::Clamp(getROF(lower - mROFAddTimeErr), BCType{0}, maxROF); + BCType last = o2::gpu::CAMath::Clamp(getROF(upper + mROFAddTimeErr), BCType{0}, maxROF); + + if (first <= last && !intersectROF(first, lower, upper)) { + ++first; + } + if (last >= first && !intersectROF(last, lower, upper)) { + --last; + } + return {first, first <= last ? static_cast(last - first + 1) : BCType{0}}; + } + #ifndef GPUCA_GPUCODE GPUh() std::string asString() const { diff --git a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/TimeFrame.h b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/TimeFrame.h index 3fef2dc640cbc..9a31e4014bff5 100644 --- a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/TimeFrame.h +++ b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/TimeFrame.h @@ -212,6 +212,10 @@ struct TimeFrame { virtual size_t getNumberOfNeighbours() const; size_t getNumberOfTracks() const; size_t getNumberOfUsedClusters() const; + void resetTrackExtensionCounters(); + void addTrackExtensionCounters(size_t nTracks, size_t nClusters); + size_t getNExtendedTracks() const { return mNExtendedTracks; } + size_t getNExtendedClusters() const { return mNExtendedClusters; } /// memory management void setMemoryPool(std::shared_ptr pool); @@ -280,6 +284,8 @@ struct TimeFrame { std::vector> mCells; bounded_vector mTracks; bounded_vector mTracksLabel; + size_t mNExtendedTracks = 0; + size_t mNExtendedClusters = 0; std::vector> mCellsNeighbours; std::vector> mCellsNeighboursTopology; std::vector> mCellsLookupTable; @@ -604,6 +610,20 @@ inline size_t TimeFrame::getNumberOfUsedClusters() const return nClusters; } +template +inline void TimeFrame::resetTrackExtensionCounters() +{ + mNExtendedTracks = 0; + mNExtendedClusters = 0; +} + +template +inline void TimeFrame::addTrackExtensionCounters(size_t nTracks, size_t nClusters) +{ + mNExtendedTracks += nTracks; + mNExtendedClusters += nClusters; +} + } // namespace its } // namespace o2 diff --git a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/TrackExtensionHypothesis.h b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/TrackExtensionHypothesis.h new file mode 100644 index 0000000000000..ddfc88c239f36 --- /dev/null +++ b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/TrackExtensionHypothesis.h @@ -0,0 +1,56 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +#ifndef TRACKINGITSU_INCLUDE_TRACKEXTENSIONHYPOTHESIS_H_ +#define TRACKINGITSU_INCLUDE_TRACKEXTENSIONHYPOTHESIS_H_ + +#include + +#include "GPUCommonDef.h" +#include "DataFormatsITS/TimeEstBC.h" +#include "ITStracking/Constants.h" +#include "ITStracking/TrackITSInternal.h" +#include "ReconstructionDataFormats/Track.h" + +namespace o2::its +{ + +template +struct TrackExtensionHypothesis { + TrackExtensionHypothesis() = default; + GPUhdi() TrackExtensionHypothesis(const TrackITSInternal& track, bool outward) + { + initialiseFromTrack(track, outward); + } + + GPUhdi() void initialiseFromTrack(const TrackITSInternal& track, bool outward) + { + param = outward ? track.paramOut : track.paramIn; + time = track.time; + chi2 = track.getChi2(); + nClusters = track.getNClusters(); + edgeLayer = outward ? track.getLastClusterLayer() : track.getFirstClusterLayer(); + for (int iLayer{0}; iLayer < NLayers; ++iLayer) { + clusters[iLayer] = track.getClusterIndex(iLayer); + } + } + + o2::track::TrackParCov param; + std::array clusters{}; + TimeEstBC time; + float chi2{0.f}; + int nClusters{0}; + int edgeLayer{constants::UnusedIndex}; +}; + +} // namespace o2::its + +#endif /* TRACKINGITSU_INCLUDE_TRACKEXTENSIONHYPOTHESIS_H_ */ diff --git a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/TrackFollower.h b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/TrackFollower.h new file mode 100644 index 0000000000000..4d3b197074909 --- /dev/null +++ b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/TrackFollower.h @@ -0,0 +1,233 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +/// \file TrackFollower.h +/// \brief Hypothesis search used by CPU and GPU track extension. + +#ifndef TRACKINGITSU_INCLUDE_TRACKFOLLOWER_H_ +#define TRACKINGITSU_INCLUDE_TRACKFOLLOWER_H_ + +#include "GPUCommonDef.h" +#include "GPUCommonMath.h" +#include "DetectorsBase/Propagator.h" + +#include "ITStracking/Cluster.h" +#include "ITStracking/Constants.h" +#include "ITStracking/IndexTableUtils.h" +#include "ITStracking/MathUtils.h" +#include "ITStracking/ROFLookupTables.h" +#include "ITStracking/TrackExtensionHypothesis.h" +#include "ITStracking/TrackHelpers.h" + +namespace o2::its +{ + +template +GPUhdi() void keepTrackExtensionHypothesis(const TrackExtensionHypothesis& hypo, + TrackExtensionHypothesis* keptHypotheses, + int& nKeptHypotheses, + const int maxHypotheses) +{ + if (nKeptHypotheses < maxHypotheses) { + keptHypotheses[nKeptHypotheses++] = hypo; + return; + } + + int worst{0}; + for (int i{1}; i < nKeptHypotheses; ++i) { + if (track::isBetter(keptHypotheses[worst].nClusters, keptHypotheses[worst].chi2, keptHypotheses[i].nClusters, keptHypotheses[i].chi2)) { + worst = i; + } + } + if (track::isBetter(hypo.nClusters, hypo.chi2, keptHypotheses[worst].nClusters, keptHypotheses[worst].chi2)) { + keptHypotheses[worst] = hypo; + } +} + +template +GPUhdi() void updateTrackFromExtensionHypothesis(const TrackExtensionHypothesis& hypo, + const bool outward, + const int nLayers, + TrackITSInternal& track) +{ + if (outward) { + track.paramOut = hypo.param; + } else { + track.paramIn = hypo.param; + } + track.time = hypo.time; + track.setChi2(hypo.chi2); + for (int iLayer{0}; iLayer < nLayers; ++iLayer) { + if (track.getClusterIndex(iLayer) == constants::UnusedIndex && hypo.clusters[iLayer] != constants::UnusedIndex) { + track.setClusterIndex(iLayer, hypo.clusters[iLayer]); + } + } +} + +template +GPUhdi() bool followTrackExtensionDirection(const TrackExtensionHypothesis& startHypothesis, + const IndexTableUtils& utils, + const typename ROFMaskTable::View& rofMask, + const typename ROFOverlapTable::View& rofOverlaps, + const Cluster* const* clusters, + const unsigned char* const* usedClusters, + const int* const* clustersIndexTables, + const int* const* ROFClusters, + const TrackingFrameInfo* const* trackingFrameInfo, + const float* layerRadii, + const float* layerxX0, + const int nLayers, + const int phiBins, + const int maxHypothesesConfig, + const float bz, + const float maxChi2ClusterAttachment, + const float maxChi2NDF, + const float nSigmaCutPhi, + const float nSigmaCutZ, + const bool outward, + const o2::base::Propagator* propagator, + const o2::base::PropagatorF::MatCorrType matCorrType, + TrackExtensionHypothesis* activeHypotheses, + TrackExtensionHypothesis* nextHypotheses, + TrackExtensionHypothesis& bestHypothesis) +{ + const int step = outward ? 1 : -1; + const int end = outward ? nLayers - 1 : 0; + const int maxHypotheses = o2::gpu::CAMath::Max(maxHypothesesConfig, 1); + int nActive{1}; + int nNext{0}; + activeHypotheses[0] = startHypothesis; + + const int tableSize = utils.getNphiBins() * utils.getNzBins() + 1; + for (int iLayer = activeHypotheses[0].edgeLayer + step; nActive > 0; iLayer += step) { + if ((step > 0 && iLayer > end) || (step < 0 && iLayer < end)) { + break; + } + nNext = 0; + for (int iHypo{0}; iHypo < nActive; ++iHypo) { + auto hypo = activeHypotheses[iHypo]; + const float r = layerRadii[iLayer]; + float x{-999.f}; + if (!hypo.param.getXatLabR(r, x, bz, o2::track::DirAuto) || x <= 0.f) { + continue; + } + + if (!propagator->propagateToX(hypo.param, x, bz, o2::base::PropagatorF::MAX_SIN_PHI, + o2::base::PropagatorF::MAX_STEP, matCorrType)) { + continue; + } + if (matCorrType == o2::base::PropagatorF::MatCorrType::USEMatCorrNONE && + !hypo.param.correctForMaterial(layerxX0[iLayer], layerxX0[iLayer] * constants::Radl * constants::Rho, true)) { + continue; + } + + const float ePhi{o2::gpu::CAMath::Sqrt(hypo.param.getSigmaSnp2() / hypo.param.getCsp2())}; + const float eZ{o2::gpu::CAMath::Sqrt(hypo.param.getSigmaZ2())}; + const int4 selectedBins = getBinsRect(iLayer, hypo.param.getPhi(), hypo.param.getZ(), nSigmaCutZ * eZ, nSigmaCutPhi * ePhi, utils); + if (selectedBins.x < 0) { + continue; + } + + int phiBinsNum = selectedBins.w - selectedBins.y + 1; + if (phiBinsNum < 0) { + phiBinsNum += phiBins; + } + + const auto rofRange = rofOverlaps.getLayer(iLayer).getROFRange(hypo.time); + for (int rof = rofRange.getFirstEntry(); rof < rofRange.getEntriesBound(); ++rof) { + if (!rofMask.isROFEnabled(iLayer, rof)) { + continue; + } + const int rofStart = ROFClusters[iLayer][rof]; + const int nLayerClusters = ROFClusters[iLayer][rof + 1] - rofStart; + if (nLayerClusters <= 0) { + continue; + } + const Cluster* layerClusters = clusters[iLayer] + rofStart; + const int* indexTable = clustersIndexTables[iLayer] + rof * tableSize; + const int zBinRange = selectedBins.z - selectedBins.x + 1; + for (int iPhiCount = 0; iPhiCount < phiBinsNum; ++iPhiCount) { + const int iPhiBin = (selectedBins.y + iPhiCount) % phiBins; + const int firstBinIndex = utils.getBinIndex(selectedBins.x, iPhiBin); + const int maxBinIndex = firstBinIndex + zBinRange; + const int firstRowClusterIndex = indexTable[firstBinIndex]; + const int maxRowClusterIndex = indexTable[maxBinIndex]; + for (int iNextCluster{firstRowClusterIndex}; iNextCluster < maxRowClusterIndex; ++iNextCluster) { + if (iNextCluster >= nLayerClusters) { + break; + } + const Cluster& nextCluster = layerClusters[iNextCluster]; + if (usedClusters[iLayer][nextCluster.clusterId]) { + continue; + } + + const TrackingFrameInfo& trackingHit = trackingFrameInfo[iLayer][nextCluster.clusterId]; + auto updated = hypo; + if (!updated.param.rotate(trackingHit.alphaTrackingFrame) || + !propagator->propagateToX(updated.param, trackingHit.xTrackingFrame, bz, + o2::base::PropagatorF::MAX_SIN_PHI, + o2::base::PropagatorF::MAX_STEP, + matCorrType)) { + continue; + } + + const auto predChi2 = updated.param.getPredictedChi2Quiet(trackingHit.positionTrackingFrame, trackingHit.covarianceTrackingFrame); + if (predChi2 < 0.f || predChi2 > maxChi2ClusterAttachment) { + continue; + } + if (!updated.param.o2::track::TrackParCov::update(trackingHit.positionTrackingFrame, trackingHit.covarianceTrackingFrame)) { + continue; + } + updated.chi2 += predChi2; + updated.clusters[iLayer] = nextCluster.clusterId; + ++updated.nClusters; + updated.edgeLayer = iLayer; + updated.time += rofOverlaps.getLayer(iLayer).getROFTimeBounds(rof, true); + keepTrackExtensionHypothesis(updated, nextHypotheses, nNext, maxHypotheses); + } + } + } + keepTrackExtensionHypothesis(hypo, nextHypotheses, nNext, maxHypotheses); + } + if (nNext == 0) { + break; + } + for (int iHypo{0}; iHypo < nNext; ++iHypo) { + activeHypotheses[iHypo] = nextHypotheses[iHypo]; + } + nActive = nNext; + } + + const TrackExtensionHypothesis* bestHypo{nullptr}; + for (int iHypo{0}; iHypo < nActive; ++iHypo) { + const auto& hypo = activeHypotheses[iHypo]; + if (hypo.nClusters == startHypothesis.nClusters) { + continue; + } + const float maxChi2 = maxChi2NDF * static_cast(hypo.nClusters * 2 - 5); + if (hypo.chi2 >= maxChi2) { + continue; + } + if (!bestHypo || track::isBetter(hypo.nClusters, hypo.chi2, bestHypo->nClusters, bestHypo->chi2)) { + bestHypo = &hypo; + } + } + if (!bestHypo) { + return false; + } + + bestHypothesis = *bestHypo; + return true; +} + +} // namespace o2::its + +#endif // TRACKINGITSU_INCLUDE_TRACKFOLLOWER_H_ diff --git a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/TrackHelpers.h b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/TrackHelpers.h index d244b39ff9d11..4077e57f7d571 100644 --- a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/TrackHelpers.h +++ b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/TrackHelpers.h @@ -16,25 +16,33 @@ #ifndef O2_ITS_TRACKING_TRACKHELPERS_H_ #define O2_ITS_TRACKING_TRACKHELPERS_H_ +#include "CommonConstants/MathConstants.h" #include "DataFormatsITS/TrackITS.h" #include "ITStracking/Cell.h" #include "ITStracking/Cluster.h" #include "ITStracking/Constants.h" #include "ITStracking/MathUtils.h" +#include "ITStracking/TrackITSInternal.h" #include "DetectorsBase/Propagator.h" #include "ReconstructionDataFormats/Track.h" namespace o2::its::track { -// Prefer 1) longer track 2) sorted in chi2 +GPUhdi() bool isBetter(const int nClustersA, const float chi2A, const int nClustersB, const float chi2B) +{ + return (nClustersA > nClustersB) || (nClustersA == nClustersB && chi2A < chi2B); +} + GPUhdi() bool isBetter(const o2::its::TrackITS& a, const o2::its::TrackITS& b) { - const auto ncla = a.getNumberOfClusters(); - const auto nclb = b.getNumberOfClusters(); - // is a as long as b ? then decide on chi2 - // otherwise prefer longer - return (ncla == nclb) ? (a.getChi2() < b.getChi2()) : ncla > nclb; + return isBetter(a.getNumberOfClusters(), a.getChi2(), b.getNumberOfClusters(), b.getChi2()); +} + +template +GPUhdi() bool isBetter(const o2::its::TrackITSInternal& a, const o2::its::TrackITSInternal& b) +{ + return isBetter(a.getNumberOfClusters(), a.getChi2(), b.getNumberOfClusters(), b.getChi2()); } // Find the populated interior layer closest to the radial midpoint. @@ -58,7 +66,7 @@ GPUdi() int selectReseedMidLayer(int minLayer, int maxLayer, const float* layerR return midLayer; } -GPUdi() void resetTrackCovariance(TrackITSExt& track) +GPUdi() void resetTrackCovariance(o2::track::TrackParCov& track) { track.resetCovariance(); track.setCov(track.getQ2Pt() * track.getQ2Pt() * track.getCov()[o2::track::CovLabels::kSigQ2Pt2], o2::track::CovLabels::kSigQ2Pt2); @@ -97,19 +105,20 @@ GPUdi() o2::track::TrackParCov buildTrackSeed(const Cluster& cluster1, } template -GPUdi() TrackITSExt seedTrackForRefit(const TrackSeed& seed, - const TrackingFrameInfo* const* foundTrackingFrameInfo, - const Cluster* const* unsortedClusters, - const float* layerRadii, - const float bz, - const int reseedIfShorter) +GPUdi() TrackITSInternal seedTrackForRefit(const TrackSeed& seed, + const TrackingFrameInfo* const* foundTrackingFrameInfo, + const Cluster* const* unsortedClusters, + const float* layerRadii, + const float bz, + const int reseedIfShorter) { - TrackITSExt temporaryTrack(seed); + TrackITSInternal temporaryTrack; + temporaryTrack.paramIn = static_cast(seed); int lrMin = NLayers; int lrMax = 0; for (int iL{0}; iL < NLayers; ++iL) { const int idx = seed.getCluster(iL); - temporaryTrack.setExternalClusterIndex(iL, idx, idx != constants::UnusedIndex); + temporaryTrack.setClusterIndex(iL, idx); if (idx != constants::UnusedIndex) { lrMin = o2::gpu::CAMath::Min(lrMin, iL); lrMax = o2::gpu::CAMath::Max(lrMax, iL); @@ -123,15 +132,17 @@ GPUdi() TrackITSExt seedTrackForRefit(const TrackSeed& seed, const auto& cluster0TF = foundTrackingFrameInfo[lrMin][seed.getCluster(lrMin)]; const auto& cluster1GL = unsortedClusters[lrMid][seed.getCluster(lrMid)]; const auto& cluster2GL = unsortedClusters[lrMax][seed.getCluster(lrMax)]; - temporaryTrack.getParamIn() = buildTrackSeed(cluster2GL, cluster1GL, cluster0TF, bz, true); + temporaryTrack.paramIn = buildTrackSeed(cluster2GL, cluster1GL, cluster0TF, bz, true); } } - resetTrackCovariance(temporaryTrack); + resetTrackCovariance(temporaryTrack.paramIn); return temporaryTrack; } -GPUdi() bool fitTrack(TrackITSExt& trk, +template +GPUdi() bool fitTrack(TrackITSInternal& trk, + o2::track::TrackParCov& param, int start, int end, int step, @@ -154,43 +165,43 @@ GPUdi() bool fitTrack(TrackITSExt& trk, const TrackingFrameInfo& trackingHit = tfInfos[iLayer][trk.getClusterIndex(iLayer)]; if (linRef) { - if (!trk.o2::track::TrackParCovF::rotate(trackingHit.alphaTrackingFrame, *linRef, bz)) { + if (!param.o2::track::TrackParCovF::rotate(trackingHit.alphaTrackingFrame, *linRef, bz)) { return false; } - if (!propagator->propagateToX(trk, *linRef, trackingHit.xTrackingFrame, bz, + if (!propagator->propagateToX(param, *linRef, trackingHit.xTrackingFrame, bz, o2::base::PropagatorImpl::MAX_SIN_PHI, o2::base::PropagatorImpl::MAX_STEP, matCorrType)) { return false; } if (matCorrType == o2::base::PropagatorF::MatCorrType::USEMatCorrNONE) { - if (!trk.correctForMaterial(*linRef, layerxX0[iLayer], layerxX0[iLayer] * constants::Radl * constants::Rho, true)) { + if (!param.correctForMaterial(*linRef, layerxX0[iLayer], layerxX0[iLayer] * constants::Radl * constants::Rho, true)) { continue; } } } else { - if (!trk.o2::track::TrackParCovF::rotate(trackingHit.alphaTrackingFrame)) { + if (!param.o2::track::TrackParCovF::rotate(trackingHit.alphaTrackingFrame)) { return false; } - if (!propagator->propagateToX(trk, trackingHit.xTrackingFrame, bz, + if (!propagator->propagateToX(param, trackingHit.xTrackingFrame, bz, o2::base::PropagatorImpl::MAX_SIN_PHI, o2::base::PropagatorImpl::MAX_STEP, matCorrType)) { return false; } if (matCorrType == o2::base::PropagatorF::MatCorrType::USEMatCorrNONE) { - if (!trk.correctForMaterial(layerxX0[iLayer], layerxX0[iLayer] * constants::Radl * constants::Rho, true)) { + if (!param.correctForMaterial(layerxX0[iLayer], layerxX0[iLayer] * constants::Radl * constants::Rho, true)) { continue; } } } - const auto predChi2{trk.getPredictedChi2Quiet(trackingHit.positionTrackingFrame, trackingHit.covarianceTrackingFrame)}; + const auto predChi2{param.getPredictedChi2Quiet(trackingHit.positionTrackingFrame, trackingHit.covarianceTrackingFrame)}; if ((nCl >= 3 && predChi2 > chi2clcut) || predChi2 < 0.f) { return false; } trk.setChi2(trk.getChi2() + predChi2); - if (!trk.o2::track::TrackParCov::update(trackingHit.positionTrackingFrame, trackingHit.covarianceTrackingFrame)) { + if (!param.o2::track::TrackParCov::update(trackingHit.positionTrackingFrame, trackingHit.covarianceTrackingFrame)) { return false; } if (linRef && shiftRefToCluster) { @@ -200,36 +211,30 @@ GPUdi() bool fitTrack(TrackITSExt& trk, nCl++; } - return o2::gpu::CAMath::Abs(trk.getQ2Pt()) < maxQoverPt && trk.getChi2() < chi2ndfcut * (float)((nCl * 2) - 5); + return o2::gpu::CAMath::Abs(param.getQ2Pt()) < maxQoverPt && trk.getChi2() < chi2ndfcut * (float)((nCl * 2) - 5); } template -GPUdi() bool refitTrack(const TrackSeed& trackSeed, - TrackITSExt& temporaryTrack, - float chi2clcut, - float chi2ndfcut, - const float bz, +GPUdi() bool refitTrack(TrackITSInternal& track, const TrackingFrameInfo* const* tfInfos, - const Cluster* const* clusters, const float* layerxX0, - const float* layerRadii, - const float* minPt, + const int nLayers, + const float bz, + const float chi2clcut, + const float chi2ndfcut, const o2::base::Propagator* propagator, const o2::base::PropagatorF::MatCorrType matCorrType, - const int reseedIfShorter, const bool shiftRefToCluster, - const bool repeatRefitOut) + const bool repeatRefitOut, + const float minPt = -1.f) { - temporaryTrack = seedTrackForRefit(trackSeed, - tfInfos, - clusters, - layerRadii, - bz, - reseedIfShorter); - o2::track::TrackPar linRef{temporaryTrack}; - bool fitSuccess = fitTrack(temporaryTrack, + o2::track::TrackPar linRef{track.paramIn}; + resetTrackCovariance(track.paramIn); + track.setChi2(0); + bool fitSuccess = fitTrack(track, + track.paramIn, 0, - NLayers, + nLayers, 1, chi2clcut, chi2ndfcut, @@ -245,12 +250,14 @@ GPUdi() bool refitTrack(const TrackSeed& trackSeed, if (!fitSuccess) { return false; } - temporaryTrack.getParamOut() = temporaryTrack.getParamIn(); - linRef = temporaryTrack.getParamOut(); // use refitted track as lin.reference - resetTrackCovariance(temporaryTrack); - temporaryTrack.setChi2(0); - fitSuccess = fitTrack(temporaryTrack, - NLayers - 1, + + track.paramOut = track.paramIn; + linRef = track.paramOut; + resetTrackCovariance(track.paramIn); + track.setChi2(0); + fitSuccess = fitTrack(track, + track.paramIn, + nLayers - 1, -1, -1, chi2clcut, @@ -264,36 +271,81 @@ GPUdi() bool refitTrack(const TrackSeed& trackSeed, matCorrType, &linRef, shiftRefToCluster); - if (!fitSuccess || temporaryTrack.getPt() < minPt[NLayers - temporaryTrack.getNClusters()]) { + if (!fitSuccess) { + return false; + } + if (minPt > 0.f && track.getPt() < minPt) { return false; } if (repeatRefitOut) { // repeat outward refit seeding and linearizing with the stable inward fit result - o2::track::TrackParCov saveInw{temporaryTrack}; + o2::track::TrackParCov saveInw{track.paramIn}; linRef = saveInw; // use refitted track as lin.reference - float saveChi2 = temporaryTrack.getChi2(); - track::resetTrackCovariance(temporaryTrack); - temporaryTrack.setChi2(0); - fitSuccess = o2::its::track::fitTrack(temporaryTrack, - 0, - NLayers, - 1, - chi2clcut, - chi2ndfcut, - o2::constants::math::VeryBig, - 0, - bz, - tfInfos, - layerxX0, - propagator, - matCorrType, - &linRef, - shiftRefToCluster); + float saveChi2 = track.getChi2(); + track.paramOut = saveInw; + track::resetTrackCovariance(track.paramOut); + track.setChi2(0); + fitSuccess = fitTrack(track, + track.paramOut, + 0, + nLayers, + 1, + chi2clcut, + chi2ndfcut, + o2::constants::math::VeryBig, + 0, + bz, + tfInfos, + layerxX0, + propagator, + matCorrType, + &linRef, + shiftRefToCluster); if (!fitSuccess) { return false; } - temporaryTrack.getParamOut() = temporaryTrack.getParamIn(); - temporaryTrack.getParamIn() = saveInw; - temporaryTrack.setChi2(saveChi2); + track.paramIn = saveInw; + track.setChi2(saveChi2); + } + return true; +} + +template +GPUdi() bool refitTrack(const TrackSeed& trackSeed, + TrackITSInternal& temporaryTrack, + float chi2clcut, + float chi2ndfcut, + const float bz, + const TrackingFrameInfo* const* tfInfos, + const Cluster* const* clusters, + const float* layerxX0, + const float* layerRadii, + const float* minPt, + const o2::base::Propagator* propagator, + const o2::base::PropagatorF::MatCorrType matCorrType, + const int reseedIfShorter, + const bool shiftRefToCluster, + const bool repeatRefitOut) +{ + temporaryTrack = seedTrackForRefit(trackSeed, + tfInfos, + clusters, + layerRadii, + bz, + reseedIfShorter); + bool fitSuccess = refitTrack(temporaryTrack, + tfInfos, + layerxX0, + NLayers, + bz, + chi2clcut, + chi2ndfcut, + propagator, + matCorrType, + shiftRefToCluster, + repeatRefitOut, + minPt[NLayers - temporaryTrack.getNClusters()]); + if (!fitSuccess) { + return false; } return true; } diff --git a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/TrackITSInternal.h b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/TrackITSInternal.h new file mode 100644 index 0000000000000..1f56afb1468c8 --- /dev/null +++ b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/TrackITSInternal.h @@ -0,0 +1,113 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +#ifndef TRACKINGITSU_INCLUDE_TRACKITSINTERNAL_H_ +#define TRACKINGITSU_INCLUDE_TRACKITSINTERNAL_H_ + +#include + +#include "GPUCommonDef.h" +#include "DataFormatsITS/TrackITS.h" +#include "DataFormatsITS/TimeEstBC.h" +#include "ITStracking/Constants.h" +#include "ReconstructionDataFormats/Track.h" + +namespace o2::its +{ + +template +struct TrackITSInternal { + GPUhdi() TrackITSInternal() { resetClusters(); } + + GPUhdi() void resetClusters() + { + for (int iLayer{0}; iLayer < NLayers; ++iLayer) { + clusters[iLayer] = constants::UnusedIndex; + } + nClusters = 0; + } + + GPUhdi() int getClusterIndex(int layer) const { return clusters[layer]; } + + GPUhdi() void setClusterIndex(int layer, int cluster) + { + if (clusters[layer] == constants::UnusedIndex && cluster != constants::UnusedIndex) { + ++nClusters; + } else if (clusters[layer] != constants::UnusedIndex && cluster == constants::UnusedIndex) { + --nClusters; + } + clusters[layer] = cluster; + } + + GPUhdi() int getNClusters() const { return nClusters; } + GPUhdi() int getNumberOfClusters() const { return nClusters; } + GPUhdi() float getChi2() const { return chi2; } + GPUhdi() void setChi2(float value) { chi2 = value; } + GPUdi() float getPt() const { return paramIn.getPt(); } + + GPUhdi() uint32_t getPattern() const + { + uint32_t pattern{0}; + for (int iLayer{0}; iLayer < NLayers; ++iLayer) { + if (clusters[iLayer] != constants::UnusedIndex) { + pattern |= (0x1u << iLayer); + } + } + return pattern; + } + + GPUhdi() int getFirstClusterLayer() const + { + for (int iLayer{0}; iLayer < NLayers; ++iLayer) { + if (clusters[iLayer] != constants::UnusedIndex) { + return iLayer; + } + } + return constants::UnusedIndex; + } + + GPUhdi() int getLastClusterLayer() const + { + for (int iLayer{NLayers - 1}; iLayer >= 0; --iLayer) { + if (clusters[iLayer] != constants::UnusedIndex) { + return iLayer; + } + } + return constants::UnusedIndex; + } + + o2::track::TrackParCov paramIn; + o2::track::TrackParCov paramOut; + std::array clusters{}; + TimeEstBC time; + float chi2{0.f}; + int nClusters{0}; +}; + +template +GPUhdi() TrackITSExt makeTrackITSExt(const TrackITSInternal& track) +{ + TrackITSExt out; + out.getParamIn() = track.paramIn; + out.getParamOut() = track.paramOut; + out.setChi2(track.chi2); + out.getTimeStamp() = track.time.makeSymmetrical(); + for (int iLayer{0}; iLayer < NLayers; ++iLayer) { + if (track.clusters[iLayer] != constants::UnusedIndex) { + out.setExternalClusterIndex(iLayer, track.clusters[iLayer], true); + } + } + return out; +} + +} // namespace o2::its + +#endif /* TRACKINGITSU_INCLUDE_TRACKITSINTERNAL_H_ */ diff --git a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/Tracker.h b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/Tracker.h index 240b0eb1e2f63..2362b6b2d9816 100644 --- a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/Tracker.h +++ b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/Tracker.h @@ -99,10 +99,11 @@ class Tracker Celling, Neighbouring, Roading, + Extending, NSteps, }; Steps mCurStep{TFInit}; - static constexpr std::array StateNames{"TimeFrame initialisation", "Tracklet finding", "Cell finding", "Neighbour finding", "Road finding"}; + static constexpr std::array StateNames{"TimeFrame initialisation", "Tracklet finding", "Cell finding", "Neighbour finding", "Road finding", "Track extending"}; std::vector> mTimingStats; void addTimingStatCurStep(int iteration, double timeMs); }; diff --git a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/TrackerTraits.h b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/TrackerTraits.h index f536e86fe95d5..aa28355c429a6 100644 --- a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/TrackerTraits.h +++ b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/TrackerTraits.h @@ -17,12 +17,16 @@ #define TRACKINGITSU_INCLUDE_TRACKERTRAITS_H_ #include +#include +#include "DetectorsBase/Propagator.h" #include "ITStracking/Configuration.h" #include "ITStracking/IndexTableUtils.h" #include "ITStracking/TimeFrame.h" #include "ITStracking/Cell.h" #include "ITStracking/BoundedAllocator.h" +#include "ITStracking/TrackExtensionHypothesis.h" +#include "ITStracking/TrackITSInternal.h" // #define OPTIMISATION_OUTPUT @@ -85,6 +89,24 @@ class TrackerTraits std::shared_ptr mTaskArena; protected: + struct TrackFollowerScratch { + explicit TrackFollowerScratch(std::pmr::memory_resource* memoryResource) + : activeHypotheses(memoryResource), nextHypotheses(memoryResource) + { + } + + bounded_vector> activeHypotheses; + bounded_vector> nextHypotheses; + }; + + bool finaliseTrackSeed(const TrackSeedN& seed, + TrackITSExt& track, + const int iteration, + const TrackingFrameInfo* const* tfInfos, + const Cluster* const* unsortedClusters, + const o2::base::Propagator* propagator); + bool trackFollowing(TrackITSInternal* track, bool outward, const int iteration, TrackFollowerScratch& scratch); + o2::gpu::GPUChainITS* mChain = nullptr; TimeFrame* mTimeFrame; std::vector mTrkParams; diff --git a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/TrackingConfigParam.h b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/TrackingConfigParam.h index 69aa3c5fdaf06..d967736168c0b 100644 --- a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/TrackingConfigParam.h +++ b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/TrackingConfigParam.h @@ -96,15 +96,20 @@ struct TrackerParamConfig : public o2::conf::ConfigurableParamHelper::max(); bool dropTFUponFailure = false; - bool fataliseUponFailure = true; // granular management of the fatalisation in async mode + bool fataliseUponFailure = true; // granular management of the fatalisation in async mode // Selections on tracks sharing clusters - bool allowSharingFirstCluster = false; // allow first cluster sharing among tracks + bool allowSharingFirstCluster = false; // allow first cluster sharing among tracks float sharedClusterMaxDeltaPhi = 0.05f; // Maximum allowed delta phi at the cluster position float sharedClusterMaxDeltaEta = 0.03f; // Maximum allowed delta eta at the cluster position bool sharedClusterOppositeSign = false; // Require opposite sign of the tracklets diff --git a/Detectors/ITSMFT/ITS/tracking/src/Configuration.cxx b/Detectors/ITSMFT/ITS/tracking/src/Configuration.cxx index f07a8f3394c05..eb4fc399adc07 100644 --- a/Detectors/ITSMFT/ITS/tracking/src/Configuration.cxx +++ b/Detectors/ITSMFT/ITS/tracking/src/Configuration.cxx @@ -59,6 +59,16 @@ std::string TrackingParameters::asString() const if (MaxHoles) { str += std::format(" MaxHoles:{} HoleMask:{}", MaxHoles, HoleLayerMask.asString()); } + if (PassFlags[IterationStep::TrackFollowerTop] || PassFlags[IterationStep::TrackFollowerBot]) { + const bool top = PassFlags[IterationStep::TrackFollowerTop], bot = PassFlags[IterationStep::TrackFollowerBot]; + str += std::format(" TrackFollower:{} NSigmaZ/Phi:{:.2f}/{:.2f}", + top && bot ? "mix" : (top ? "top" : "bot"), + TrackFollowerNSigmaCutZ, + TrackFollowerNSigmaCutPhi); + if (TrackFollowerMaxHypotheses > 1) { + str += std::format(" MaxHypotheses:{}", TrackFollowerMaxHypotheses); + } + } if (std::numeric_limits::max() != MaxMemory) { str += std::format(" MemLimit {:.2f} GB", double(MaxMemory) / constants::GB); } @@ -191,7 +201,6 @@ std::vector TrackingMode::getTrackingParameters(TrackingMode if (trackParams.size() > 3 && tc.doUPCIteration) { trackParams[3].PassFlags.set(IterationStep::UseUPCMask, IterationStep::RebuildClusterLUT, IterationStep::SelectUPCVertices); } - float bFactor = std::abs(o2::base::Propagator::Instance()->getNominalBz()) / 5.0066791f; float bFactorTracklets = bFactor < 0.01f ? 1.f : bFactor; // for tracklets only @@ -207,6 +216,9 @@ std::vector TrackingMode::getTrackingParameters(TrackingMode p.RepeatRefitOut = tc.repeatRefitOut; p.ShiftRefToCluster = tc.shiftRefToCluster; p.CreateArtefactLabels = tc.createArtefactLabels; + p.TrackFollowerNSigmaCutZ = tc.trackFollowerNSigmaCutZ; + p.TrackFollowerNSigmaCutPhi = tc.trackFollowerNSigmaCutPhi; + p.TrackFollowerMaxHypotheses = std::max(1, tc.trackFollowerMaxHypotheses); p.PrintMemory = tc.printMemory; p.MaxMemory = tc.maxMemory; @@ -221,6 +233,12 @@ std::vector TrackingMode::getTrackingParameters(TrackingMode if (iter < constants::MaxIter) { p.MaxHoles = tc.maxHolesIter[iter]; p.HoleLayerMask = tc.holeLayerMaskIter[iter]; + if (tc.trackFollowerTop[iter]) { + p.PassFlags.set(IterationStep::TrackFollowerTop); + } + if (tc.trackFollowerBot[iter]) { + p.PassFlags.set(IterationStep::TrackFollowerBot); + } } if (tc.useMatCorrTGeo) { diff --git a/Detectors/ITSMFT/ITS/tracking/src/TimeFrame.cxx b/Detectors/ITSMFT/ITS/tracking/src/TimeFrame.cxx index 8375004cbfbad..232fc0178a30d 100644 --- a/Detectors/ITSMFT/ITS/tracking/src/TimeFrame.cxx +++ b/Detectors/ITSMFT/ITS/tracking/src/TimeFrame.cxx @@ -474,6 +474,7 @@ void TimeFrame::setFrameworkAllocator(ExternalAllocator* ext) template void TimeFrame::wipe() { + resetTrackExtensionCounters(); deepVectorClear(mTracks); deepVectorClear(mTracklets); deepVectorClear(mCells); diff --git a/Detectors/ITSMFT/ITS/tracking/src/Tracker.cxx b/Detectors/ITSMFT/ITS/tracking/src/Tracker.cxx index f17d961fc7bb7..95c2de41bfa97 100644 --- a/Detectors/ITSMFT/ITS/tracking/src/Tracker.cxx +++ b/Detectors/ITSMFT/ITS/tracking/src/Tracker.cxx @@ -20,6 +20,7 @@ #include "ITStracking/TrackingConfigParam.h" #include +#include #include #include #include @@ -74,6 +75,7 @@ float Tracker::clustersToTracks(const LogFunc& logger, const LogFunc& e float timeFrame{0.}, timeTracklets{0.}, timeCells{0.}, timeNeighbours{0.}, timeRoads{0.}; size_t nTracklets{0}, nCells{0}, nNeighbours{0}; int nTracks{-static_cast(mTimeFrame->getNumberOfTracks())}; + mTimeFrame->resetTrackExtensionCounters(); iVertex = std::min(maxNvertices, 0); logger(std::format("==== ITS {} Tracking iteration {} summary ====", mTraits->getName(), iteration)); total += timeFrame = evaluateTask(&Tracker::initialiseTimeFrame, StateNames[mCurStep = TFInit], iteration, evalLog, iteration); @@ -91,6 +93,9 @@ float Tracker::clustersToTracks(const LogFunc& logger, const LogFunc& e logger(std::format(" - Cell finding: {} cells found in {:.2f} ms", nCells, timeCells)); logger(std::format(" - Neighbours finding: {} neighbours found in {:.2f} ms", nNeighbours, timeNeighbours)); logger(std::format(" - Track finding: {} tracks found in {:.2f} ms", nTracks + mTimeFrame->getNumberOfTracks(), timeRoads)); + if (mTrkParams[iteration].PassFlags[IterationStep::TrackFollowerTop] || mTrkParams[iteration].PassFlags[IterationStep::TrackFollowerBot]) { + logger(std::format(" - Integrated track extension: {} tracks accepted using {} clusters", mTimeFrame->getNExtendedTracks(), mTimeFrame->getNExtendedClusters())); + } total += timeTracklets + timeCells + timeNeighbours + timeRoads; } } catch (const BoundedMemoryResource::MemoryLimitExceeded& err) { diff --git a/Detectors/ITSMFT/ITS/tracking/src/TrackerTraits.cxx b/Detectors/ITSMFT/ITS/tracking/src/TrackerTraits.cxx index c4439dc74d29e..93fcb78bb7379 100644 --- a/Detectors/ITSMFT/ITS/tracking/src/TrackerTraits.cxx +++ b/Detectors/ITSMFT/ITS/tracking/src/TrackerTraits.cxx @@ -14,9 +14,15 @@ /// #include +#include +#include #include +#include +#include +#include #include #include +#include #include #include @@ -30,6 +36,7 @@ #include "ITStracking/LayerMask.h" #include "ITStracking/ROFLookupTables.h" #include "ITStracking/TrackerTraits.h" +#include "ITStracking/TrackFollower.h" #include "ITStracking/TrackHelpers.h" #include "ITStracking/Tracklet.h" @@ -657,6 +664,105 @@ void TrackerTraits::processNeighbours(int iteration, int defaultCellTop }); } +template +bool TrackerTraits::finaliseTrackSeed(const TrackSeedN& seed, + TrackITSExt& track, + const int iteration, + const TrackingFrameInfo* const* tfInfos, + const Cluster* const* unsortedClusters, + const o2::base::Propagator* propagator) +{ + TrackITSInternal internalTrack; + if (!track::refitTrack(seed, + internalTrack, + mTrkParams[iteration].MaxChi2ClusterAttachment, + mTrkParams[iteration].MaxChi2NDF, + mBz, + tfInfos, + unsortedClusters, + mTrkParams[iteration].LayerxX0.data(), + mTrkParams[iteration].LayerRadii.data(), + mTrkParams[iteration].MinPt.data(), + propagator, + mTrkParams[iteration].CorrType, + mTrkParams[iteration].ReseedIfShorter, + mTrkParams[iteration].ShiftRefToCluster, + mTrkParams[iteration].RepeatRefitOut)) { + return false; + } + + const bool extendTop = mTrkParams[iteration].PassFlags[IterationStep::TrackFollowerTop]; + const bool extendBot = mTrkParams[iteration].PassFlags[IterationStep::TrackFollowerBot]; + if (!extendTop && !extendBot) { + track = makeTrackITSExt(internalTrack); + return true; + } + + const auto backup = internalTrack; + auto best = internalTrack; + uint32_t bestDiff{0}; + TrackFollowerScratch scratch{mMemoryPool.get()}; + const uint32_t lastLayer = static_cast(mTrkParams[iteration].NLayers - 1); + + auto finaliseExtensionCandidate = [&](TrackITSInternal& candidate) { + const auto diff = (candidate.getPattern() & ~backup.getPattern()) & TrackITS::getLayerPatternMask(); + if (!diff || + !track::refitTrack(candidate, + tfInfos, + mTrkParams[iteration].LayerxX0.data(), + mTrkParams[iteration].NLayers, + mBz, + mTrkParams[iteration].MaxChi2ClusterAttachment, + mTrkParams[iteration].MaxChi2NDF, + propagator, + mTrkParams[iteration].CorrType, + mTrkParams[iteration].ShiftRefToCluster, + mTrkParams[iteration].RepeatRefitOut)) { + return; + } + if (track::isBetter(candidate, best)) { + best = candidate; + bestDiff = diff; + } + }; + + std::optional> topResult, botResult; + if (extendTop && backup.getLastClusterLayer() != lastLayer) { + auto candidate = backup; + if (trackFollowing(&candidate, true, iteration, scratch)) { + topResult = candidate; + finaliseExtensionCandidate(candidate); + } + } + if (extendBot && backup.getFirstClusterLayer() != 0) { + auto candidate = backup; + if (trackFollowing(&candidate, false, iteration, scratch)) { + botResult = candidate; + finaliseExtensionCandidate(candidate); + } + } + if (extendTop && extendBot) { + if (topResult && topResult->getFirstClusterLayer() != 0) { + auto candidate = *topResult; + if (trackFollowing(&candidate, false, iteration, scratch)) { + finaliseExtensionCandidate(candidate); + } + } + if (botResult && botResult->getLastClusterLayer() != lastLayer) { + auto candidate = *botResult; + if (trackFollowing(&candidate, true, iteration, scratch)) { + finaliseExtensionCandidate(candidate); + } + } + } + + track = makeTrackITSExt(best); + if (bestDiff) { + track.setExtendedLayerPattern(bestDiff); + } + return true; +} + template void TrackerTraits::findRoads(const int iteration) { @@ -717,65 +823,34 @@ void TrackerTraits::findRoads(const int iteration) bounded_vector tracks(mMemoryPool.get()); mTaskArena->execute([&] { - auto forSeed = [&](auto Tag, int iSeed, int offset = 0) { - TrackITSExt temporaryTrack; - bool refitSuccess = track::refitTrack(trackSeeds[iSeed], - temporaryTrack, - mTrkParams[iteration].MaxChi2ClusterAttachment, - mTrkParams[iteration].MaxChi2NDF, - mBz, - tfInfos, - unsortedClusters, - mTrkParams[iteration].LayerxX0.data(), - mTrkParams[iteration].LayerRadii.data(), - mTrkParams[iteration].MinPt.data(), - propagator, - mTrkParams[iteration].CorrType, - mTrkParams[iteration].ReseedIfShorter, - mTrkParams[iteration].ShiftRefToCluster, - mTrkParams[iteration].RepeatRefitOut); - - if (refitSuccess) { - if constexpr (decltype(Tag)::value == PassMode::OnePass::value) { - tracks.push_back(temporaryTrack); - } else if constexpr (decltype(Tag)::value == PassMode::TwoPassCount::value) { - // nothing to do - } else if constexpr (decltype(Tag)::value == PassMode::TwoPassInsert::value) { - tracks[offset] = temporaryTrack; - } else { - static_assert(false, "Unknown mode!"); - } - return 1; - } - return 0; - }; - const int nSeeds = static_cast(trackSeeds.size()); - if (mTaskArena->max_concurrency() <= 1) { - for (int iSeed{0}; iSeed < nSeeds; ++iSeed) { - forSeed(PassMode::OnePass{}, iSeed); - } - } else { - // The double-pass allows us to avoid sizeable memory spikes - bounded_vector perSeedCount(nSeeds + 1, 0, mMemoryPool.get()); - tbb::parallel_for(0, nSeeds, [&](const int iSeed) { - perSeedCount[iSeed] = forSeed(PassMode::TwoPassCount{}, iSeed); - }); - - std::exclusive_scan(perSeedCount.begin(), perSeedCount.end(), perSeedCount.begin(), 0); - auto totalTracks{perSeedCount.back()}; - if (totalTracks == 0) { - return; - } - tracks.resize(totalTracks); - - tbb::parallel_for(0, nSeeds, [&](const int iSeed) { - if (perSeedCount[iSeed] == perSeedCount[iSeed + 1]) { - return; + const int nWorkers = std::min(static_cast(mTaskArena->max_concurrency()), nSeeds); + const int chunkSize = std::min(nSeeds, std::clamp(nSeeds / (16 * nWorkers), 256, 4096)); + std::atomic nextSeed{0}; + std::mutex tracksMutex; + tbb::parallel_for(0, nWorkers, [&](const int) { + bounded_vector localTracks(mMemoryPool.get()); + localTracks.reserve(chunkSize); + while (true) { + const int firstSeed = nextSeed.fetch_add(chunkSize, std::memory_order_relaxed); + if (firstSeed >= nSeeds) { + break; } - forSeed(PassMode::TwoPassInsert{}, iSeed, perSeedCount[iSeed]); - }); - } + const int lastSeed = std::min(firstSeed + chunkSize, nSeeds); + for (int iSeed{firstSeed}; iSeed < lastSeed; ++iSeed) { + TrackITSExt temporaryTrack; + if (finaliseTrackSeed(trackSeeds[iSeed], temporaryTrack, iteration, tfInfos, unsortedClusters, propagator)) { + localTracks.push_back(temporaryTrack); + } + } + if (!localTracks.empty()) { + std::lock_guard lock{tracksMutex}; + tracks.insert(tracks.end(), std::make_move_iterator(localTracks.begin()), std::make_move_iterator(localTracks.end())); + localTracks.clear(); + } + } + deepVectorClear(localTracks); + }); deepVectorClear(trackSeeds); }); @@ -790,7 +865,9 @@ void TrackerTraits::findRoads(const int iteration) } template -void TrackerTraits::acceptTracks(int iteration, bounded_vector& tracks, bounded_vector>& firstClusters) +void TrackerTraits::acceptTracks(int iteration, + bounded_vector& tracks, + bounded_vector>& firstClusters) { auto& trks = mTimeFrame->getTracks(); trks.reserve(trks.size() + tracks.size()); @@ -851,8 +928,15 @@ void TrackerTraits::acceptTracks(int iteration, bounded_vector smallestROFHalf) { track.getTimeStamp().setTimeStampError(smallestROFHalf); } - track.setUserField(0); - track.getParamOut().setUserField(0); + const auto diff = track.getExtendedLayerPattern(); + if (diff) { + size_t nExtendedClusters = 0; + for (int iLayer{0}; iLayer < mTrkParams[iteration].NLayers; ++iLayer) { + nExtendedClusters += static_cast(diff & (0x1u << iLayer)); + } + mTimeFrame->addTrackExtensionCounters(1, nExtendedClusters); + } + track.clearExtendedLayerPattern(); trks.emplace_back(track); if (mTrkParams[iteration].AllowSharingFirstCluster) { @@ -907,6 +991,66 @@ void TrackerTraits::markTracks(int iteration) } } +template +bool TrackerTraits::trackFollowing(TrackITSInternal* track, bool outward, const int iteration, TrackFollowerScratch& scratch) +{ + const int maxHypotheses = std::max(1, mTrkParams[iteration].TrackFollowerMaxHypotheses); + if (static_cast(scratch.activeHypotheses.size()) < maxHypotheses) { + scratch.activeHypotheses.resize(maxHypotheses); + } + if (static_cast(scratch.nextHypotheses.size()) < maxHypotheses) { + scratch.nextHypotheses.resize(maxHypotheses); + } + + const Cluster* clustersPtrs[NLayers]{}; + const unsigned char* usedClustersPtrs[NLayers]{}; + const int* clustersIndexTablesPtrs[NLayers]{}; + const int* rofClustersPtrs[NLayers]{}; + const TrackingFrameInfo* tfInfoPtrs[NLayers]{}; + for (int iLayer{0}; iLayer < NLayers; ++iLayer) { + clustersPtrs[iLayer] = mTimeFrame->getClusters()[iLayer].data(); + usedClustersPtrs[iLayer] = mTimeFrame->getUsedClusters(iLayer).data(); + clustersIndexTablesPtrs[iLayer] = mTimeFrame->getIndexTable(0, iLayer).data(); + rofClustersPtrs[iLayer] = mTimeFrame->getROFrameClusters(iLayer).data(); + tfInfoPtrs[iLayer] = mTimeFrame->getTrackingFrameInfoOnLayer(iLayer).data(); + } + + auto startHypothesis = TrackExtensionHypothesis{*track, outward}; + TrackExtensionHypothesis bestHypothesis; + const bool ok = followTrackExtensionDirection( + startHypothesis, + mTimeFrame->getIndexTableUtils(), + mTimeFrame->getROFMaskView(), + mTimeFrame->getROFOverlapTableView(), + clustersPtrs, + usedClustersPtrs, + clustersIndexTablesPtrs, + rofClustersPtrs, + tfInfoPtrs, + mTrkParams[iteration].LayerRadii.data(), + mTrkParams[iteration].LayerxX0.data(), + mTrkParams[iteration].NLayers, + mTrkParams[iteration].PhiBins, + maxHypotheses, + mBz, + mTrkParams[iteration].MaxChi2ClusterAttachment, + mTrkParams[iteration].MaxChi2NDF, + mTrkParams[iteration].TrackFollowerNSigmaCutPhi, + mTrkParams[iteration].TrackFollowerNSigmaCutZ, + outward, + o2::base::Propagator::Instance(), + mTrkParams[iteration].CorrType, + scratch.activeHypotheses.data(), + scratch.nextHypotheses.data(), + bestHypothesis); + if (!ok) { + return false; + } + + updateTrackFromExtensionHypothesis(bestHypothesis, outward, mTrkParams[iteration].NLayers, *track); + return true; +} + template void TrackerTraits::setBz(float bz) {