diff --git a/Detectors/ITSMFT/ITS/tracking/GPU/ITStrackingGPU/ClusterLinesGPU.h b/Detectors/ITSMFT/ITS/tracking/GPU/ITStrackingGPU/ClusterLinesGPU.h new file mode 100644 index 0000000000000..14914d3f7b600 --- /dev/null +++ b/Detectors/ITSMFT/ITS/tracking/GPU/ITStrackingGPU/ClusterLinesGPU.h @@ -0,0 +1,194 @@ +// 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 ClusterLinesGPU.h +/// \brief device-side line + N-line vertex fit for the GPU seeding vertexer. + +#ifndef O2_ITS_CLUSTERLINES_GPU_H +#define O2_ITS_CLUSTERLINES_GPU_H + +#include "DataFormatsITS/TimeEstBC.h" +#include "GPUCommonDef.h" +#include "GPUCommonMath.h" +#include "ITStracking/LineProjection.h" + +namespace o2::its::gpu +{ + +using LineTime = o2::its::LineTime; +using LineWindow = o2::its::LineWindow; + +struct LineProjSoA { + float* z{nullptr}; // projected z at the beamline; sort key and binary-search key, kept dense + LineTime* t{nullptr}; // time centre + half-width + int* idx{nullptr}; // sorted slot -> original line index + int* rof{nullptr}; // ROF of the line +}; + +struct VertexCand { + float x, y, z; + float rms2[6]; + float avgDist2; + int nGood; + float seed[3]; + o2::its::TimeEstBC time; + int size; + uint8_t ok; // 1 if the candidate passed the fit cuts + uint8_t keep; // 1 if it survived duplicate suppression (subset of ok) + uint8_t fine; +}; + +// Device-side line: origin point + unit direction, with a time stamp +struct GPULine { + GPUhdDefault() GPULine() = default; + + GPUhdi() GPULine(const float origin[3], const float direction[3], const o2::its::TimeEstBC& t) : mTime(t) + { + const float norm = o2::gpu::GPUCommonMath::Sqrt(direction[0] * direction[0] + + direction[1] * direction[1] + + direction[2] * direction[2]); + const float inv = norm > 0.f ? 1.f / norm : 0.f; + for (int i = 0; i < 3; ++i) { + originPoint[i] = origin[i]; + cosinesDirector[i] = direction[i] * inv; + } + } + + // Squared distance from a point to the (infinite) line: |delta - (delta.u) u|^2 + GPUhdi() static float getDistance2FromPoint(const GPULine& line, const float point[3]) + { + float delta[3]; + float proj = 0.f; + for (int i = 0; i < 3; ++i) { + delta[i] = point[i] - line.originPoint[i]; + proj += delta[i] * line.cosinesDirector[i]; + } + float d2 = 0.f; + for (int i = 0; i < 3; ++i) { + const float residual = delta[i] - proj * line.cosinesDirector[i]; + d2 += residual * residual; + } + return d2; + } + + GPUhdi() static void getDCAComponents(const GPULine& line, const float point[3], float out[6]) + { + float delta[3]; + float proj = 0.f; + for (int i = 0; i < 3; ++i) { + delta[i] = line.originPoint[i] - point[i]; + proj += delta[i] * line.cosinesDirector[i]; + } + float r[3]; + for (int i = 0; i < 3; ++i) { + r[i] = delta[i] - proj * line.cosinesDirector[i]; + } + out[0] = r[0]; // (0,0) XX + out[1] = o2::gpu::GPUCommonMath::Hypot(r[0], r[1]); // (0,1) XY + out[2] = r[1]; // (1,1) YY + out[3] = o2::gpu::GPUCommonMath::Hypot(r[0], r[2]); // (0,2) XZ + out[4] = o2::gpu::GPUCommonMath::Hypot(r[1], r[2]); // (1,2) YZ + out[5] = r[2]; // (2,2) ZZ + } + + float originPoint[3] = {0.f, 0.f, 0.f}; + float cosinesDirector[3] = {0.f, 0.f, 0.f}; + o2::its::TimeEstBC mTime; +}; + +class GPUClusterLinesFit +{ + public: + GPUhdDefault() GPUClusterLinesFit() = default; + + // Add one line's contribution: A_ij += (delta_ij*|d|^2 - d_i*d_j)/|d|^2, + // b_i += (d_i*(d.o) - |d|^2*o_i)/|d|^2. For a unit director |d|^2 == 1. + GPUhdi() void add(const GPULine& line) + { + const double d0 = line.cosinesDirector[0], d1 = line.cosinesDirector[1], d2 = line.cosinesDirector[2]; + const double o0 = line.originPoint[0], o1 = line.originPoint[1], o2 = line.originPoint[2]; + const double det = d0 * d0 + d1 * d1 + d2 * d2; // == 1 for a normalised director + if (det <= 0.) { + return; + } + if (mNContributors <= 0) { + mTime = line.mTime; + } else { + mTime += line.mTime; + } + mA[0] += (det - d0 * d0) / det; + mA[1] += (-d0 * d1) / det; + mA[2] += (-d0 * d2) / det; + mA[3] += (det - d1 * d1) / det; + mA[4] += (-d1 * d2) / det; + mA[5] += (det - d2 * d2) / det; + const double dDotO = d0 * o0 + d1 * o1 + d2 * o2; + mB[0] += (d0 * dDotO - det * o0) / det; + mB[1] += (d1 * dDotO - det * o1) / det; + mB[2] += (d2 * dDotO - det * o2) / det; + ++mNContributors; + } + + // Solve the symmetric system and write the vertex (= -A^-1 B) + GPUhdi() bool solve(float vertex[3]) const + { + const double a = mA[0], b = mA[1], c = mA[2], d = mA[3], e = mA[4], f = mA[5]; + const double c00 = d * f - e * e; + const double c01 = c * e - b * f; + const double c02 = b * e - c * d; + const double c11 = a * f - c * c; + const double c12 = b * c - a * e; + const double c22 = a * d - b * b; + const double det = a * c00 + b * c01 + c * c02; + if (o2::gpu::GPUCommonMath::Abs(det) < 1.e-12) { + return false; + } + const double invDet = 1. / det; + const double x0 = (c00 * mB[0] + c01 * mB[1] + c02 * mB[2]) * invDet; + const double x1 = (c01 * mB[0] + c11 * mB[1] + c12 * mB[2]) * invDet; + const double x2 = (c02 * mB[0] + c12 * mB[1] + c22 * mB[2]) * invDet; + vertex[0] = static_cast(-x0); + vertex[1] = static_cast(-x1); + vertex[2] = static_cast(-x2); + return true; + } + + GPUhdi() void addResidual(const GPULine& line, const float vertex[3]) + { + float dca[6]; + GPULine::getDCAComponents(line, vertex, dca); + const float d2 = GPULine::getDistance2FromPoint(line, vertex); + ++mResidualCount; + const float inv = 1.f / static_cast(mResidualCount); + for (int i = 0; i < 6; ++i) { + mRMS2[i] += (dca[i] - mRMS2[i]) * inv; + } + mAvgDistance2 += (d2 - mAvgDistance2) * inv; + } + + GPUhdi() int getNContributors() const { return mNContributors; } + GPUhdi() const float* getRMS2() const { return mRMS2; } // Packed symmetric covariance in {XX, XY, YY, XZ, YZ, ZZ} order + GPUhdi() float getAvgDistance2() const { return mAvgDistance2; } + GPUhdi() const o2::its::TimeEstBC& getTimeStamp() const { return mTime; } + + private: + double mA[6] = {0., 0., 0., 0., 0., 0.}; + double mB[3] = {0., 0., 0.}; + int mNContributors = 0; + float mRMS2[6] = {0.f, 0.f, 0.f, 0.f, 0.f, 0.f}; + float mAvgDistance2 = 0.f; + int mResidualCount = 0; + o2::its::TimeEstBC mTime; +}; + +} // namespace o2::its::gpu + +#endif /* O2_ITS_CLUSTERLINES_GPU_H */ diff --git a/Detectors/ITSMFT/ITS/tracking/GPU/ITStrackingGPU/TimeFrameGPU.h b/Detectors/ITSMFT/ITS/tracking/GPU/ITStrackingGPU/TimeFrameGPU.h index 0d70158b9bdb8..3f50f09186805 100644 --- a/Detectors/ITSMFT/ITS/tracking/GPU/ITStrackingGPU/TimeFrameGPU.h +++ b/Detectors/ITSMFT/ITS/tracking/GPU/ITStrackingGPU/TimeFrameGPU.h @@ -21,6 +21,7 @@ #include "ITStracking/Configuration.h" #include "ITStracking/TrackExtensionHypothesis.h" #include "ITStrackingGPU/Utils.h" +#include "ITStrackingGPU/ClusterLinesGPU.h" namespace o2::its::gpu { @@ -54,10 +55,14 @@ class TimeFrameGPU : public TimeFrame void createTrackingFrameInfoDeviceArray(const int = NLayers); void loadUnsortedClustersDevice(const int); void createUnsortedClustersDeviceArray(const int = NLayers); - void loadClustersDevice(const int); void createClustersDeviceArray(const int = NLayers); void loadClustersIndexTables(const int); void createClustersIndexTablesArray(const int = NLayers); + void createClustersDevice(const int); + void createClustersIndexTables(const int); + void createClusterRadiiDevice(); + void uploadClusterRadii(); + void sortClustersDevice(const int layer, const TrackingParameters& trkParam); void createUsedClustersDevice(const int); void createUsedClustersDeviceArray(const int = NLayers); void loadUsedClustersDevice(); @@ -87,6 +92,35 @@ class TimeFrameGPU : public TimeFrame void createTrackExtensionScratchDevice(const int nThreads, const int maxHypotheses); void downloadTrackITSExtDevice(); + // Seeding-vertexer + void createClusterOwnersDeviceArray(); + void createClusterOwnersDevice(); + void resetClusterOwnersDevice(); + void createClusterSortScratchDevice(const int layer); + + protected: + void prepareClusters(const TrackingParameters& trkParam, const int maxLayers) override + { + if (maxLayers < NLayers) { // only if former seeding vertexer is run + TimeFrame::prepareClusters(trkParam, maxLayers); + } + } + void allocateClusterSortStorage(const TrackingParameters& trkParam, const int maxLayers) override + { + if (maxLayers < NLayers) { // only if former seeding vertexer is run + TimeFrame::allocateClusterSortStorage(trkParam, maxLayers); + } + } + + public: + void createLinesDevice(const int nCells); + void createDiamondDevice(const Vertex& diamond); + unsigned int downloadLinesDevice(); + unsigned int getNLines(); + const auto& getHostLines() const { return mLinesHost; } + const auto& getHostLineRof() const { return mLineRofHost; } + const auto& getHostLineClusters() const { return mLineClustersHost; } + /// synchronization auto& getStream(const size_t stream) { return mGpuStreams[stream]; } auto& getStreams() { return mGpuStreams; } @@ -111,6 +145,19 @@ class TimeFrameGPU : public TimeFrame auto& getTrackITSExt() { return mTrackITSExt; } auto& getTrackIndices() { return mTrackIndices; } Vertex* getDeviceVertices() { return mPrimaryVerticesDevice; } + int* getDeviceROFramesClusters(const int layer) { return mROFramesClustersDevice[layer]; } + int* getDeviceClusterSortKeys(const int layer) { return mClusterSortKeysDevice[layer]; } + int* getDeviceClusterSortPerm(const int layer) { return mClusterSortPermDevice[layer]; } + Cluster* getDeviceUnsortedClusters(const int layer) { return mUnsortedClustersDevice[layer]; } + Cluster* getDeviceClusters(const int layer) { return mClustersDevice[layer]; } + int* getDeviceClustersIndexTable(const int layer) { return mClustersIndexTablesDevice[layer]; } + const float* getDeviceMinRs() const { return mClusterMinRDevice; } + const float* getDeviceMaxRs() const { return mClusterMaxRDevice; } + int* getDeviceROFramesPV() { return mROFramesPVDevice; } + unsigned char* getDeviceUsedClusters(const int); + const o2::base::Propagator* getChainPropagator(); + bool arePersistentTablesLoaded() { return mPersistentTablesLoaded; } + void setPersistentTablesLoaded(bool setValue) { mPersistentTablesLoaded = setValue; } // Hybrid TrackITSExt* getDeviceTrackITSExt() { return mTrackITSExtDevice; } @@ -119,6 +166,45 @@ class TimeFrameGPU : public TimeFrame TrackExtensionHypothesis* getDeviceNextTrackExtensionHypotheses() { return mNextTrackExtensionHypothesesDevice; } int* getDeviceNeighboursLUT(const int layer) { return mNeighboursLUTDevice[layer]; } CellNeighbour** getDeviceArrayNeighbours() { return mNeighboursDeviceArray; } + unsigned long long** getDeviceArrayClusterOwners() { return mClusterOwnersDeviceArray; } + GPULine* getDeviceLines() { return mLinesDevice; } + int* getDeviceLineSlots() { return mLineSlotsDevice; } + int* getDeviceLineRof() { return mLineRofDevice; } + int* getDeviceLineClusters() { return mLineClustersDevice; } + float* getDeviceLineChi2() { return mLineChi2Device; } + float* getDeviceLinePt() { return mLinePtDevice; } + float* getDeviceLineZs() { return mLineZsDevice; } + gpu::LineTime* getDeviceLineTimes() { return mLineTimesDevice; } + int* getDeviceLineSortedIdx() { return mLinesSortedIdx; } + LineProjSoA getLineProjSoA() { return {mLineZsDevice, mLineTimesDevice, mLinesSortedIdx, mLineRofDevice}; } + LineProjSoA getLineProjSortedSoA() { return {mLineZsSortedDevice, mLineTimesSortedDevice, mLinesSortedIdx, mLineRofSortedDevice}; } + int* getDeviceRofLineOffsets() { return mRofLineOffsetsDevice; } + int* getDeviceLineDensity() { return mLineDensityDevice; } + gpu::LineWindow* getDeviceLineWin() { return mLineWinDevice; } + uint8_t* getDeviceLineIsPeak() { return mLineIsPeakDevice; } + int* getDeviceLineDensityFine() { return mLineDensityFineDevice; } + gpu::LineWindow* getDeviceLineWinFine() { return mLineWinFineDevice; } + uint8_t* getDeviceLineIsPeakFine() { return mLineIsPeakFineDevice; } + int* getDevicePeakScan() { return mPeakScanDevice; } + int* getDevicePeakLineIdx() { return mPeakLineIdxDevice; } + int* getDevicePeakOffsets() { return mPeakOffsetsDevice; } + const int* getDeviceNPeaks() { return mPeakOffsetsDevice + this->getNrof(1); } + VertexCand* getDeviceVertexCands() { return mVertexCandsDevice; } + int* getDeviceMemberOffsets() { return mMemberOffsetsDevice; } + int* getDeviceMemberLines() { return mMemberLinesDevice; } + int downloadVertexCandsDevice(); + int getNMembers() const { return mNMembers; } + void downloadMemberOffsetsDevice(); // (MC only) + void createMemberLinesMCDevice(const int nMembers); // (MC only) + void downloadMemberLinesDevice(); // (MC only) + const auto& getHostVertexCands() const { return mVertexCandsHost; } + const auto& getHostPeakOffsets() const { return mPeakOffsetsHost; } + const auto& getHostMemberOffsets() const { return mMemberOffsetsHost; } + const auto& getHostMemberLines() const { return mMemberLinesHost; } + std::vector& getLineLabelFlat() { return mLineLabelFlatHost; } + const std::vector& getLineLabelFlat() const { return mLineLabelFlatHost; } + Vertex* getDeviceDiamond() { return mDiamondDevice; } + std::array& getDeviceNeighboursAll() { return mNeighboursDevice; } CellNeighbour* getDeviceNeighbours(const int layer) { return mNeighboursDevice[layer]; } const TrackingFrameInfo** getDeviceArrayTrackingFrameInfo() const { return mTrackingFrameInfoDeviceArray; } const Cluster** getDeviceArrayClusters() const { return mClustersDeviceArray; } @@ -215,6 +301,11 @@ class TimeFrameGPU : public TimeFrame const int** mClustersIndexTablesDeviceArray{nullptr}; uint8_t** mUsedClustersDeviceArray{nullptr}; const int** mROFramesClustersDeviceArray{nullptr}; + int* mROFramesPVDevice; + std::array mClusterSortKeysDevice{}; + std::array mClusterSortPermDevice{}; + float* mClusterMinRDevice{nullptr}; + float* mClusterMaxRDevice{nullptr}; std::array mTrackletsDevice{}; std::array mTrackletsLUTDevice{}; std::array mCellsLUTDevice{}; @@ -239,6 +330,47 @@ class TimeFrameGPU : public TimeFrame CellNeighbour** mNeighboursDeviceArray{nullptr}; std::array mTrackingFrameInfoDevice{}; const TrackingFrameInfo** mTrackingFrameInfoDeviceArray{nullptr}; + std::array mClusterOwnersDevice{}; + unsigned long long** mClusterOwnersDeviceArray{nullptr}; + int* mLineSlotsDevice{nullptr}; + GPULine* mLinesDevice{nullptr}; + int* mLineRofDevice{nullptr}; + int* mLineClustersDevice{nullptr}; + float* mLineChi2Device{nullptr}; + float* mLinePtDevice{nullptr}; + float* mLineZsDevice{nullptr}; + gpu::LineTime* mLineTimesDevice{nullptr}; + float* mLineZsSortedDevice{nullptr}; + gpu::LineTime* mLineTimesSortedDevice{nullptr}; + int* mLinesSortedIdx{nullptr}; + int* mLineRofSortedDevice{nullptr}; // per (sorted) line's ROF + int* mRofLineOffsetsDevice{nullptr}; // CSR offsets into the (rof,z)-sorted lines, size nRofs+1 + int* mLineDensityDevice{nullptr}; // per (sorted) line: count of time-compatible neighbours in its z-window + gpu::LineWindow* mLineWinDevice{nullptr}; // per (sorted) line: [lo,hi) bounds of its z-window (sorted coords) + uint8_t* mLineIsPeakDevice{nullptr}; // per (sorted) line: 1 if it is a local density peak (vertex candidate) + int* mLineDensityFineDevice{nullptr}; + gpu::LineWindow* mLineWinFineDevice{nullptr}; + uint8_t* mLineIsPeakFineDevice{nullptr}; + int* mPeakScanDevice{nullptr}; // per (sorted) line: number of peaks strictly before it + int* mPeakLineIdxDevice{nullptr}; // per peak slot: the sorted line index it came from + int* mPeakOffsetsDevice{nullptr}; // CSR offsets into the compacted peaks + VertexCand* mVertexCandsDevice{nullptr}; + int* mMemberOffsetsDevice{nullptr}; + int* mMemberLinesDevice{nullptr}; + int mNLinesCapacity{0}; // = nCells the line buffers were sized for + std::vector mLinesHost; + std::vector mLineRofHost; + std::vector mLineClustersHost; + std::vector mVertexCandsHost; + std::vector mPeakOffsetsHost; + std::vector mMemberOffsetsHost; + std::vector mMemberLinesHost; + std::vector mLineLabelFlatHost; + int mNMembers{0}; + Vertex* mDiamondDevice{nullptr}; + bool mPersistentTablesLoaded{false}; + std::bitset mUnsortedClustersUploaded{}; + std::bitset mTrackingFrameInfoUploaded{}; // State Streams mGpuStreams; diff --git a/Detectors/ITSMFT/ITS/tracking/GPU/ITStrackingGPU/TrackerTraitsGPU.h b/Detectors/ITSMFT/ITS/tracking/GPU/ITStrackingGPU/TrackerTraitsGPU.h index 0d84662666632..74ba41e27d49b 100644 --- a/Detectors/ITSMFT/ITS/tracking/GPU/ITStrackingGPU/TrackerTraitsGPU.h +++ b/Detectors/ITSMFT/ITS/tracking/GPU/ITStrackingGPU/TrackerTraitsGPU.h @@ -29,6 +29,9 @@ class TrackerTraitsGPU final : public TrackerTraits void adoptTimeFrame(TimeFrame* tf) final; void initialiseTimeFrame(const int iteration) final; + void computeVertexCandidates(const int iteration) final; + void computeVertices(const int iteration) final; + void computeLayerTracklets(const int iteration, int) final; void computeLayerCells(const int iteration) final; void findCellsNeighbours(const int iteration) final; diff --git a/Detectors/ITSMFT/ITS/tracking/GPU/ITStrackingGPU/TrackingKernels.h b/Detectors/ITSMFT/ITS/tracking/GPU/ITStrackingGPU/TrackingKernels.h index 94950c04877b9..135399afd8968 100644 --- a/Detectors/ITSMFT/ITS/tracking/GPU/ITStrackingGPU/TrackingKernels.h +++ b/Detectors/ITSMFT/ITS/tracking/GPU/ITStrackingGPU/TrackingKernels.h @@ -22,6 +22,7 @@ #include "ITStracking/TrackingTopology.h" #include "ITStracking/TrackExtensionHypothesis.h" #include "ITStrackingGPU/Utils.h" +#include "ITStrackingGPU/ClusterLinesGPU.h" #include "DetectorsBase/Propagator.h" namespace o2::its @@ -37,6 +38,10 @@ class IndexTableUtils; class Cluster; class TrackITSExt; class ExternalAllocator; +namespace gpu +{ +struct GPULine; +} template struct TrackingKernels { @@ -49,6 +54,7 @@ struct TrackingKernels { const typename ROFVertexLookupTable::View& vertexLUT, const int vertexId, const Vertex* vertices, + const bool vtxMode, const Cluster** clusters, const std::vector& nClusters, const int** ROFClusters, @@ -64,8 +70,8 @@ struct TrackingKernels { const typename TrackingTopology::View topology, bounded_vector& linkPhiCuts, const float resolutionPV, - std::array& minR, - std::array& maxR, + const float* minRs, + const float* maxRs, bounded_vector& resolutions, std::vector& radii, bounded_vector& linkMSAngles, @@ -86,6 +92,7 @@ struct TrackingKernels { const float bz, const float maxChi2ClusterAttachment, const float cellDeltaTanLambdaSigma, + const float cellDeltaPhiCut, const float nSigmaCut, const float* layerxX0, o2::its::ExternalAllocator* alloc, @@ -180,5 +187,135 @@ int finalizeCellNeighboursHandler(CellNeighbour* cellNeighbours, o2::its::ExternalAllocator* alloc, gpu::Stream& stream); +template +void sortClustersHandler(const Cluster* unsorted, + Cluster* sorted, + const int* clusterOffsets, + int* indexTable, + const IndexTableUtils* utils, + const typename ROFMaskTable::View& rofMask, + float beamX, float beamY, + int zBins, int phiBins, int nRofs, int nClustersLayer, int iLayer, + float* minRadiusLayer, float* maxRadiusLayer, + int* keys, + int* perm, + o2::its::ExternalAllocator* alloc, + gpu::Stream& stream); + +template +void registerClusterOwnershipHandler(const CellSeed* cellsLayersDevice, + const int nCells, + unsigned long long** clusterOwnersDeviceArray, + gpu::Stream& stream); + +template +void linearizeCellsToLinesHandler(const int nCells, + const CellSeed* cells, + const unsigned long long* const* clusterOwners, + const int* rofFramesClustersL1, + const int nRofsL1, + const int ownedClustersCut, + gpu::GPULine* lines, + int* lineRof, + int* lineClusters, + int* lineSlots, + const float beamX, + const float beamY, + const float maxZ, + const float minPt, + float* linesZs, + gpu::LineTime* lineTimes, + float* lineChi2, + float* linePt, + o2::its::ExternalAllocator* alloc, + gpu::Stream& stream); + +template +void sortLinesHandler(const int nLines, + const int nRofs, + const gpu::LineProjSoA soa, + const gpu::LineProjSoA sortedSoa, + const int* lineRof, + int* rofOffsets, + o2::its::ExternalAllocator* alloc, + gpu::Stream& stream); + +template +void scanDensityHandler(const int nLines, + const gpu::LineProjSoA sortedSoa, + const int* rofOffsets, + int* density, + gpu::LineWindow* win, + const float zWindow, + gpu::Stream& stream); + +template +void findPeaksHandler(const int nLines, + const int nRofs, + const gpu::LineProjSoA sortedSoa, + const int* rofOffsets, + const int* density, + const gpu::LineWindow* win, + uint8_t* isPeak, + const int* densityFine, + const gpu::LineWindow* winFine, + const int fineMinDensity, + uint8_t* isPeakFine, + int* peakScan, + int* peakLineIdx, + int* peakOffsets, + o2::its::ExternalAllocator* alloc, + gpu::Stream& stream); + +template +void fitPeaksHandler(const int* nPeaksDevice, + const int* peakLineIdx, + const gpu::LineWindow* win, + const gpu::LineProjSoA sortedSoa, + const gpu::GPULine* lines, + const float* lineChi2, + const float* linePt, + const float goodLineChi2Cut, + const float goodLinePtCut, + const float pairCut2, + const float nSigmaCut, + const int minContributors, + const float beamX, + const float beamY, + const uint8_t* isPeakFine, + const float fineMaxDrift, + gpu::VertexCand* cands, + gpu::Stream& stream); + +template +void dedupVertexCandidatesHandler(const int* nPeaksDevice, + const int* peakLineIdx, + const int* peakOffsets, + const gpu::LineProjSoA sortedSoa, + const float duplicateZCut, + const float duplicateZScale, + gpu::VertexCand* cands, + gpu::Stream& stream); + +// MC-only member collection (feeds the host majority-vote vertex labels). +template +void scanMemberOffsetsHandler(const gpu::VertexCand* cands, + int* memberOffsets, + const int nLines, + o2::its::ExternalAllocator* alloc, + gpu::Stream& stream); + +template +void collectLinesForMCHandler(const int* nPeaksDevice, + const int* peakLineIdx, + const gpu::LineWindow* win, + const gpu::LineProjSoA sortedSoa, + const gpu::GPULine* lines, + const float pairCut2, + const gpu::VertexCand* cands, + const int* memberOffsets, + int* memberLines, + gpu::Stream& stream); + } // namespace o2::its #endif // ITSTRACKINGGPU_TRACKINGKERNELS_H_ diff --git a/Detectors/ITSMFT/ITS/tracking/GPU/ITStrackingGPU/Utils.h b/Detectors/ITSMFT/ITS/tracking/GPU/ITStrackingGPU/Utils.h index e6909b28a687a..117a7d0d52932 100644 --- a/Detectors/ITSMFT/ITS/tracking/GPU/ITStrackingGPU/Utils.h +++ b/Detectors/ITSMFT/ITS/tracking/GPU/ITStrackingGPU/Utils.h @@ -38,10 +38,17 @@ #endif #ifdef ITS_GPU_LOG -#define GPULog(...) \ - do { \ - LOGP(info, __VA_ARGS__); \ - GPUChkErrS(cudaDeviceSynchronize()); \ +#if defined(__HIPCC__) +#define GPULogSync() GPUChkErrS(hipDeviceSynchronize()) +#elif defined(__CUDACC__) +#define GPULogSync() GPUChkErrS(cudaDeviceSynchronize()) +#else +#define GPULogSync() +#endif +#define GPULog(...) \ + do { \ + LOGP(info, __VA_ARGS__); \ + GPULogSync(); \ } while (0) #else #define GPULog(...) @@ -343,6 +350,36 @@ struct TypedAllocator { ExternalAllocator* mInternalAllocator; }; +// first i in [beg,end) with a[i] >= key +template +GPUdii() int deviceLowerBound(const T* a, int beg, int end, const T key) +{ + while (beg < end) { + const int mid = beg + (end - beg) / 2; + if (a[mid] < key) { + beg = mid + 1; + } else { + end = mid; + } + } + return beg; +} + +// first i in [beg,end) with a[i] > key +template +GPUdii() int deviceUpperBound(const T* a, int beg, int end, const T key) +{ + while (beg < end) { + const int mid = beg + (end - beg) / 2; + if (a[mid] <= key) { + beg = mid + 1; + } else { + end = mid; + } + } + return beg; +} + GPUdii() gpuSpan getClustersOnLayer(const int rof, const int totROFs, const int layer, diff --git a/Detectors/ITSMFT/ITS/tracking/GPU/cuda/TimeFrameGPU.cu b/Detectors/ITSMFT/ITS/tracking/GPU/cuda/TimeFrameGPU.cu index 593f5529cf8ea..ff5839f77b281 100644 --- a/Detectors/ITSMFT/ITS/tracking/GPU/cuda/TimeFrameGPU.cu +++ b/Detectors/ITSMFT/ITS/tracking/GPU/cuda/TimeFrameGPU.cu @@ -16,6 +16,7 @@ #include #include #include +#include #include #include @@ -23,6 +24,8 @@ #include "ITStracking/Constants.h" #include "ITStracking/BoundedAllocator.h" #include "ITStrackingGPU/Utils.h" +#include "ITStrackingGPU/ClusterLinesGPU.h" +#include "ITStrackingGPU/TrackingKernels.h" #include "GPUCommonDef.h" #include "GPUCommonMath.h" @@ -223,13 +226,6 @@ void TimeFrameGPU::createClustersDeviceArray(const int maxLayers) pinHostLayers(this->mClusters, mPinnedClusters, maxLayers); } -template -void TimeFrameGPU::loadClustersDevice(const int layer) -{ - GPUTimer timer(mGpuStreams[layer], "loading sorted clusters", layer); - uploadSlot(mClustersDevice, mClustersDeviceArray, layer, this->mClusters[layer], "sorted clusters"); -} - template void TimeFrameGPU::createClustersIndexTablesArray(const int maxLayers) { @@ -238,6 +234,22 @@ void TimeFrameGPU::createClustersIndexTablesArray(const int maxLayers) pinHostLayers(this->mIndexTables, mPinnedClustersIndexTables, maxLayers); } +template +void TimeFrameGPU::createClustersDevice(const int layer) +{ + GPUTimer timer(mGpuStreams[layer], "creating sorted clusters", layer); + createSlot(mClustersDevice, mClustersDeviceArray, layer, this->mUnsortedClusters[layer].size(), "sorted clusters"); +} + +template +void TimeFrameGPU::createClustersIndexTables(const int layer) +{ + GPUTimer timer(mGpuStreams[layer], "creating clusters index table", layer); + const int nBins = this->mIndexTableUtils.getNphiBins() * this->mIndexTableUtils.getNzBins(); + const size_t nEntries = static_cast(this->getNrof(layer)) * (nBins + 1); + createSlot(mClustersIndexTablesDevice, mClustersIndexTablesDeviceArray, layer, nEntries, "clusters index table entries"); +} + template void TimeFrameGPU::loadClustersIndexTables(const int layer) { @@ -245,6 +257,46 @@ void TimeFrameGPU::loadClustersIndexTables(const int layer) uploadSlot(mClustersIndexTablesDevice, mClustersIndexTablesDeviceArray, layer, this->mIndexTables[layer], "clusters indextable entries"); } +template +void TimeFrameGPU::createClusterRadiiDevice() +{ + GPUTimer timer("creating cluster radii bounds"); + mClusterMinRDevice = allocDevice(NLayers); + mClusterMaxRDevice = allocDevice(NLayers); +} + +template +void TimeFrameGPU::uploadClusterRadii() +{ + GPUTimer timer("loading cluster radii bounds"); + GPUChkErrS(cudaMemcpy(mClusterMinRDevice, this->getMinRs().data(), NLayers * sizeof(float), cudaMemcpyHostToDevice)); + GPUChkErrS(cudaMemcpy(mClusterMaxRDevice, this->getMaxRs().data(), NLayers * sizeof(float), cudaMemcpyHostToDevice)); +} + +template +void TimeFrameGPU::sortClustersDevice(const int layer, const TrackingParameters& trkParam) +{ + GPUTimer timer(mGpuStreams[layer], "sorting clusters on device", layer); + const int nClusters = static_cast(this->mUnsortedClusters[layer].size()); + const int nRofs = this->getNrof(layer); + const size_t nEntries = static_cast(nRofs) * (trkParam.ZBins * trkParam.PhiBins + 1); + createClustersDevice(layer); + createClustersIndexTables(layer); + if (!nClusters) { + GPUChkErrS(cudaMemsetAsync(mClustersIndexTablesDevice[layer], 0, nEntries * sizeof(int), mGpuStreams[layer].get())); + return; + } + createClusterSortScratchDevice(layer); + sortClustersHandler(mUnsortedClustersDevice[layer], mClustersDevice[layer], + mROFramesClustersDevice[layer], mClustersIndexTablesDevice[layer], + mIndexTableUtilsDevice, mDeviceROFMaskTableView, + this->mBeamPos[0], this->mBeamPos[1], + trkParam.ZBins, trkParam.PhiBins, nRofs, nClusters, layer, + mClusterMinRDevice, mClusterMaxRDevice, + mClusterSortKeysDevice[layer], mClusterSortPermDevice[layer], + this->getFrameworkAllocator(), mGpuStreams[layer]); +} + template void TimeFrameGPU::createUsedClustersDeviceArray(const int maxLayers) { @@ -375,6 +427,7 @@ template void TimeFrameGPU::uploadROFVertexLookupTable() { GPUTimer timer("updating device view of ROFVertexLookupTable"); + TimeFrame::updateROFVertexLookupTable(); const auto& hostTable = this->getROFVertexLookupTable(); const auto& hostView = this->getROFVertexLookupTableView(); using TableEntry = ROFVertexLookupTable::TableEntry; @@ -394,7 +447,7 @@ void TimeFrameGPU::createTrackletsLUTDevice(bool allocate, const int la { GPUTimer timer(mGpuStreams[layer], "creating tracklets LUTs", layer); const int fromLayer = this->mTrackingTopologyView.getLink(layer).fromLayer; - const size_t ncls = this->mClusters[fromLayer].size() + 1; + const size_t ncls = this->mUnsortedClusters[fromLayer].size() + 1; // host mClusters is empty: the sort runs on device if (allocate || mTrackletsLUTDevice[layer] == nullptr) { createSlot(mTrackletsLUTDevice, mTrackletsLUTDeviceArray, layer, ncls, "tracklets LUT"); } @@ -447,6 +500,169 @@ void TimeFrameGPU::createCellsBuffersArray() mNeighboursDeviceArray = allocSlotArray(MaxCells); } +template +void TimeFrameGPU::createClusterOwnersDeviceArray() +{ + GPUTimer timer("creating cluster owners array"); + mClusterOwnersDeviceArray = allocDevice(3); +} + +template +void TimeFrameGPU::createClusterOwnersDevice() +{ + const auto sizes = getClusterSizes(); + for (int l = 0; l < 3; ++l) { + mClusterOwnersDevice[l] = allocDeviceAsync(sizes[l], mGpuStreams[l], (o2::gpu::GPUMemoryResource::MEMORY_GPU | o2::gpu::GPUMemoryResource::MEMORY_STACK)); + GPUChkErrS(cudaMemcpyAsync(&mClusterOwnersDeviceArray[l], &mClusterOwnersDevice[l], sizeof(unsigned long long*), cudaMemcpyHostToDevice, mGpuStreams[l].get())); + } +} + +template +void TimeFrameGPU::resetClusterOwnersDevice() +{ + const auto sizes = getClusterSizes(); + for (int l = 0; l < 3; ++l) { + GPUChkErrS(cudaMemsetAsync(mClusterOwnersDevice[l], 0xFF, sizes[l] * sizeof(unsigned long long), mGpuStreams[l].get())); + } +} + +template +void TimeFrameGPU::createDiamondDevice(const Vertex& diamond) +{ + GPUTimer timer("Creating diamond device"); + mDiamondDevice = allocDevice(1, (o2::gpu::GPUMemoryResource::MEMORY_GPU | o2::gpu::GPUMemoryResource::MEMORY_STACK)); + GPUChkErrS(cudaMemcpyAsync(mDiamondDevice, &diamond, sizeof(Vertex), cudaMemcpyHostToDevice, mGpuStreams[0].get())); +} + +template +void TimeFrameGPU::createClusterSortScratchDevice(const int layer) +{ + GPUTimer timer("creating cluster-sort scratch device"); + const size_t nClusters = this->mUnsortedClusters[layer].size(); + mClusterSortKeysDevice[layer] = allocDeviceAsync(nClusters, mGpuStreams[layer]); + mClusterSortPermDevice[layer] = allocDeviceAsync(nClusters, mGpuStreams[layer]); +} + +template +void TimeFrameGPU::createLinesDevice(const int nCells) +{ + GPUTimer timer("creating lines device"); + constexpr auto kStack = (o2::gpu::GPUMemoryResource::MEMORY_GPU | o2::gpu::GPUMemoryResource::MEMORY_STACK); + mLinesDevice = allocDeviceAsync(nCells, mGpuStreams[0], kStack); + mLineRofDevice = allocDeviceAsync(nCells, mGpuStreams[0], kStack); + mLineClustersDevice = allocDeviceAsync(3 * nCells, mGpuStreams[0], kStack); + mLineChi2Device = allocDeviceAsync(nCells, mGpuStreams[0], kStack); + mLinePtDevice = allocDeviceAsync(nCells, mGpuStreams[0], kStack); + mLineSlotsDevice = allocDeviceAsync((nCells + 1), mGpuStreams[0], kStack); + mLineZsDevice = allocDeviceAsync(nCells, mGpuStreams[0], kStack); + mLineTimesDevice = allocDeviceAsync(nCells, mGpuStreams[0], kStack); + mLineZsSortedDevice = allocDeviceAsync(nCells, mGpuStreams[0], kStack); + mLineTimesSortedDevice = allocDeviceAsync(nCells, mGpuStreams[0], kStack); + mLinesSortedIdx = allocDeviceAsync(nCells, mGpuStreams[0], kStack); + mLineRofSortedDevice = allocDeviceAsync(nCells, mGpuStreams[0], kStack); + mRofLineOffsetsDevice = allocDeviceAsync((this->getNrof(1) + 1), mGpuStreams[0], kStack); + mLineDensityDevice = allocDeviceAsync(nCells, mGpuStreams[0], kStack); + mLineWinDevice = allocDeviceAsync(nCells, mGpuStreams[0], kStack); + mLineIsPeakDevice = allocDeviceAsync(nCells, mGpuStreams[0], kStack); + mLineDensityFineDevice = allocDeviceAsync(nCells, mGpuStreams[0], kStack); + mLineWinFineDevice = allocDeviceAsync(nCells, mGpuStreams[0], kStack); + mLineIsPeakFineDevice = allocDeviceAsync(nCells, mGpuStreams[0], kStack); + mPeakScanDevice = allocDeviceAsync((nCells + 1), mGpuStreams[0], kStack); + mPeakLineIdxDevice = allocDeviceAsync(nCells, mGpuStreams[0], kStack); + mPeakOffsetsDevice = allocDeviceAsync((this->getNrof(1) + 1), mGpuStreams[0], kStack); + mVertexCandsDevice = allocDeviceAsync(nCells, mGpuStreams[0], kStack); + mMemberOffsetsDevice = allocDeviceAsync((nCells + 1), mGpuStreams[0], kStack); + mNLinesCapacity = nCells; +} + +template +void TimeFrameGPU::createMemberLinesMCDevice(const int nMembers) +{ + constexpr auto kStack = (o2::gpu::GPUMemoryResource::MEMORY_GPU | o2::gpu::GPUMemoryResource::MEMORY_STACK); + // Guard against a 0-byte allocation when a pass produces no survivors. + mMemberLinesDevice = allocDeviceAsync(std::max(nMembers, 1), mGpuStreams[0], kStack); +} + +template +unsigned int TimeFrameGPU::getNLines() +{ + GPUTimer timer("getting number of lines"); + int nLinesSigned{0}; + GPUChkErrS(cudaMemcpyAsync(&nLinesSigned, mLineSlotsDevice + mNLinesCapacity, sizeof(int), cudaMemcpyDeviceToHost, mGpuStreams[0].get())); + GPUChkErrS(cudaStreamSynchronize(mGpuStreams[0].get())); // need the count before sizing the staging buffers + if (nLinesSigned < 0 || nLinesSigned > mNLinesCapacity) { + LOGP(fatal, "ITS GPU linearizer produced {} lines for a capacity of {}: bad compaction scan total.", nLinesSigned, mNLinesCapacity); + } + return static_cast(nLinesSigned); +} + +template +int TimeFrameGPU::downloadVertexCandsDevice() +{ + GPUTimer timer("downloading vertex candidates"); + mPeakOffsetsHost.resize(this->getNrof(1) + 1); + GPUChkErrS(cudaMemcpyAsync(mPeakOffsetsHost.data(), mPeakOffsetsDevice, mPeakOffsetsHost.size() * sizeof(int), cudaMemcpyDeviceToHost, mGpuStreams[0].get())); + GPUChkErrS(cudaStreamSynchronize(mGpuStreams[0].get())); // need nPeaks before sizing the staging buffers + const int nPeaks = mPeakOffsetsHost.empty() ? 0 : mPeakOffsetsHost.back(); + if (nPeaks < 0 || nPeaks > mNLinesCapacity) { + LOGP(fatal, "ITS GPU vertexer produced {} peaks for a capacity of {}: bad peak compaction total.", nPeaks, mNLinesCapacity); + } + + mVertexCandsHost.resize(nPeaks); + if (nPeaks) { + GPUChkErrS(cudaMemcpyAsync(mVertexCandsHost.data(), mVertexCandsDevice, nPeaks * sizeof(VertexCand), cudaMemcpyDeviceToHost, mGpuStreams[0].get())); + GPUChkErrS(cudaStreamSynchronize(mGpuStreams[0].get())); + } + return nPeaks; +} + +// MC only +template +void TimeFrameGPU::downloadMemberOffsetsDevice() +{ + const int nPeaks = mPeakOffsetsHost.empty() ? 0 : mPeakOffsetsHost.back(); + mMemberOffsetsHost.resize(nPeaks + 1); + if (nPeaks) { + GPUChkErrS(cudaMemcpyAsync(mMemberOffsetsHost.data(), mMemberOffsetsDevice, (nPeaks + 1) * sizeof(int), cudaMemcpyDeviceToHost, mGpuStreams[0].get())); + GPUChkErrS(cudaStreamSynchronize(mGpuStreams[0].get())); + } + mNMembers = (nPeaks && !mMemberOffsetsHost.empty()) ? mMemberOffsetsHost.back() : 0; +} + +template +void TimeFrameGPU::downloadMemberLinesDevice() +{ + mMemberLinesHost.resize(mNMembers); + if (mNMembers) { + GPUChkErrS(cudaMemcpyAsync(mMemberLinesHost.data(), mMemberLinesDevice, mNMembers * sizeof(int), cudaMemcpyDeviceToHost, mGpuStreams[0].get())); + GPUChkErrS(cudaStreamSynchronize(mGpuStreams[0].get())); + } +} + +template +unsigned int TimeFrameGPU::downloadLinesDevice() +{ + GPUTimer timer("downloading lines"); + int nLinesSigned{0}; + GPUChkErrS(cudaMemcpyAsync(&nLinesSigned, mLineSlotsDevice + mNLinesCapacity, sizeof(int), cudaMemcpyDeviceToHost, mGpuStreams[0].get())); + GPUChkErrS(cudaStreamSynchronize(mGpuStreams[0].get())); // need the count before sizing the staging buffers + if (nLinesSigned < 0 || nLinesSigned > mNLinesCapacity) { + LOGP(fatal, "ITS GPU linearizer produced {} lines for a capacity of {}: bad compaction scan total.", nLinesSigned, mNLinesCapacity); + } + const unsigned int nLines = static_cast(nLinesSigned); + GPULog("gpu-transfer: downloading {} lines, for {:.2f} MB.", nLines, nLines * (sizeof(GPULine) + sizeof(int)) / constants::MB); + mLinesHost.resize(nLines); + mLineRofHost.resize(nLines); + mLineClustersHost.resize(3 * nLines); + if (nLines) { + GPUChkErrS(cudaMemcpyAsync(mLinesHost.data(), mLinesDevice, nLines * sizeof(GPULine), cudaMemcpyDeviceToHost, mGpuStreams[0].get())); + GPUChkErrS(cudaMemcpyAsync(mLineRofHost.data(), mLineRofDevice, nLines * sizeof(int), cudaMemcpyDeviceToHost, mGpuStreams[0].get())); + GPUChkErrS(cudaMemcpyAsync(mLineClustersHost.data(), mLineClustersDevice, 3 * nLines * sizeof(int), cudaMemcpyDeviceToHost, mGpuStreams[0].get())); + GPUChkErrS(cudaStreamSynchronize(mGpuStreams[0].get())); + } + return nLines; +} + template void TimeFrameGPU::createCellsBuffers(const int layer, size_t capacity) { @@ -613,6 +829,9 @@ void TimeFrameGPU::wipe() { unregisterHostMemory(); o2::its::TimeFrame::wipe(); + mPersistentTablesLoaded = false; + mUnsortedClustersUploaded.reset(); + mTrackingFrameInfoUploaded.reset(); } template class TimeFrameGPU<7>; diff --git a/Detectors/ITSMFT/ITS/tracking/GPU/cuda/TrackerTraitsGPU.cxx b/Detectors/ITSMFT/ITS/tracking/GPU/cuda/TrackerTraitsGPU.cxx index b80fb6452f8a9..fa2441725dd7d 100644 --- a/Detectors/ITSMFT/ITS/tracking/GPU/cuda/TrackerTraitsGPU.cxx +++ b/Detectors/ITSMFT/ITS/tracking/GPU/cuda/TrackerTraitsGPU.cxx @@ -12,13 +12,47 @@ #include +#include +#include +#include +#include +#include +#include +#include + #include "ITStrackingGPU/TrackerTraitsGPU.h" #include "ITStrackingGPU/TrackingKernels.h" #include "ITStrackingGPU/LaunchGeometry.h" #include "ITStracking/Configuration.h" +#include "ITStracking/TrackingConfigParam.h" // VertexerParamConfig namespace o2::its { +namespace +{ +VertexLabel computeMainGPU(const std::vector& elements) +{ + auto composeVtxLabel = [](const o2::MCCompLabel& lbl) -> o2::MCCompLabel { + return {o2::MCCompLabel::maxTrackID(), lbl.getEventID(), lbl.getSourceID(), lbl.isFake()}; + }; + std::unordered_map frequency; + for (const auto& element : elements) { + ++frequency[composeVtxLabel(element)]; + } + o2::MCCompLabel elem{}; + size_t maxCount = 0; + for (const auto& [key, count] : frequency) { + if (count > maxCount) { + maxCount = count; + elem = key; + } + } + if (maxCount <= 1) { // need >50% + elem.setFakeFlag(); + } + return std::make_pair(elem, static_cast(maxCount) / static_cast(elements.size())); +} +} // namespace template void TrackerTraitsGPU::initialiseTimeFrame(const int iteration) @@ -32,8 +66,11 @@ void TrackerTraitsGPU::initialiseTimeFrame(const int iteration) if (this->mTrkParams[iteration].PassFlags[IterationStep::FirstPass]) { // on default stream mTimeFrameGPU->loadVertices(); - // TODO these tables can be put in persistent memory - mTimeFrameGPU->loadROFOverlapTable(); // this can be put in constant memory actually + if (!mTimeFrameGPU->arePersistentTablesLoaded()) { + mTimeFrameGPU->loadROFOverlapTable(); // this can be put in constant memory actually + mTimeFrameGPU->loadTrackingTopologies(); + mTimeFrameGPU->setPersistentTablesLoaded(true); + } mTimeFrameGPU->loadROFVertexLookupTable(); mTimeFrameGPU->loadTrackingTopologies(); // once the tables are in persistent memory just re-upload the vertex one @@ -47,10 +84,13 @@ void TrackerTraitsGPU::initialiseTimeFrame(const int iteration) mTimeFrameGPU->createTrackingFrameInfoDeviceArray(); mTimeFrameGPU->createROFrameClustersDeviceArray(); // device array + mTimeFrameGPU->createClusterRadiiDevice(); + mTimeFrameGPU->uploadClusterRadii(); mTimeFrameGPU->createTrackletsLUTDeviceArray(); mTimeFrameGPU->createTrackletsBuffersArray(); mTimeFrameGPU->createCellsBuffersArray(); mTimeFrameGPU->createCellsLUTDeviceArray(); + mTimeFrameGPU->createClusterOwnersDeviceArray(); } if (this->mTrkParams[iteration].PassFlags[IterationStep::FirstPass] || this->mTrkParams[iteration].PassFlags[IterationStep::UseUPCMask]) { mTimeFrameGPU->loadROFCutMask(iteration); @@ -73,9 +113,9 @@ void TrackerTraitsGPU::computeLayerTracklets(const int iteration, int i for (int iLayer{0}; iLayer < this->mTrkParams[iteration].NLayers; ++iLayer) { if (loadFirstPassData) { mTimeFrameGPU->createUsedClustersDevice(iLayer); - mTimeFrameGPU->loadClustersDevice(iLayer); - mTimeFrameGPU->loadClustersIndexTables(iLayer); + mTimeFrameGPU->loadUnsortedClustersDevice(iLayer); mTimeFrameGPU->loadROFrameClustersDevice(iLayer); + mTimeFrameGPU->sortClustersDevice(iLayer, this->mTrkParams[iteration]); } mTimeFrameGPU->recordEvent(iLayer); } @@ -89,10 +129,26 @@ void TrackerTraitsGPU::computeLayerTracklets(const int iteration, int i mTimeFrameGPU->pushMemoryStack(iteration); const auto nClusters = mTimeFrameGPU->getClusterSizes(); + const bool vtxMode = this->mTrkParams[iteration].PassFlags[IterationStep::SeedingVertexPass]; + const bool useDiamond = this->mTrkParams[iteration].UseDiamond; + if (useDiamond) { + const Vertex diamondVert(this->mTrkParams[iteration].Diamond, this->mTrkParams[iteration].DiamondCov, 1, 1.f); + mTimeFrameGPU->createDiamondDevice(diamondVert); + mTimeFrameGPU->recordEvent(0); + } + const Vertex* deviceVertices = useDiamond ? mTimeFrameGPU->getDeviceDiamond() : mTimeFrameGPU->getDeviceVertices(); + bounded_vector vtxPhiCuts(vtxMode ? hostTopology.nLinks : 0, + o2::its::VertexerParamConfig::Instance().phiCut, + this->getMemoryPool().get()); + auto& linkPhiCuts = vtxMode ? vtxPhiCuts : mTimeFrameGPU->getLinkPhiCuts(); + for (int linkId{0}; linkId < hostTopology.nLinks; ++linkId) { const auto link = hostTopology.getLink(linkId); mTimeFrameGPU->waitEvent(linkId, link.fromLayer); mTimeFrameGPU->waitEvent(linkId, link.toLayer); + if (useDiamond) { + mTimeFrameGPU->waitEvent(linkId, 0); // links not anchored on layer 0 must still wait for the diamond upload + } const auto key = CapacityEstimator::makeKey(SlabSite::Tracklets, iteration, iVertex + 1, linkId); const auto scale = static_cast(nClusters[link.fromLayer]); runOnSlab(mTimeFrameGPU->getCapacityEstimator(), key, scale, [&](const int capacity) { @@ -105,7 +161,8 @@ void TrackerTraitsGPU::computeLayerTracklets(const int iteration, int i mTimeFrameGPU->getDeviceROFOverlapTableView(), mTimeFrameGPU->getDeviceROFVertexLookupTableView(), iVertex, - mTimeFrameGPU->getDeviceVertices(), + deviceVertices, + vtxMode, mTimeFrameGPU->getDeviceArrayClusters(), nClusters, mTimeFrameGPU->getDeviceROFrameClusters(), @@ -119,10 +176,10 @@ void TrackerTraitsGPU::computeLayerTracklets(const int iteration, int i this->mTrkParams[iteration].PassFlags[IterationStep::SelectUPCVertices], this->mTrkParams[iteration].NSigmaCut, topology, - mTimeFrameGPU->getLinkPhiCuts(), + linkPhiCuts, this->mTrkParams[iteration].PVres, - mTimeFrameGPU->getMinRs(), - mTimeFrameGPU->getMaxRs(), + mTimeFrameGPU->getDeviceMinRs(), + mTimeFrameGPU->getDeviceMaxRs(), mTimeFrameGPU->getPositionResolutions(), this->mTrkParams[iteration].LayerRadii, mTimeFrameGPU->getLinkMSAngles(), @@ -140,7 +197,7 @@ void TrackerTraitsGPU::computeLayerCells(const int iteration) const auto hostTopology = mTimeFrameGPU->getTrackingTopologyView(); for (int iLayer{0}; iLayer < this->mTrkParams[iteration].NLayers; ++iLayer) { if (this->mTrkParams[iteration].PassFlags[IterationStep::FirstPass]) { - mTimeFrameGPU->loadUnsortedClustersDevice(iLayer); + mTimeFrameGPU->loadUnsortedClustersDevice(iLayer); // latched: a no-op if trackleting already did it mTimeFrameGPU->loadTrackingFrameInfoDevice(iLayer); } mTimeFrameGPU->recordEvent(iLayer); @@ -150,6 +207,15 @@ void TrackerTraitsGPU::computeLayerCells(const int iteration) const auto cellTopology = hostTopology.getCell(cellTopologyId); const auto first = hostTopology.getLink(cellTopology.firstLink); const auto second = hostTopology.getLink(cellTopology.secondLink); + const float cellDeltaPhiCut = this->mTrkParams[iteration].PassFlags[IterationStep::SeedingVertexPass] ? cellDeltaPhiBound(this->mBz, this->mTrkParams[iteration].CellDeltaPhiMinPt, + this->mTrkParams[iteration].LayerRadii[first.fromLayer], + this->mTrkParams[iteration].LayerRadii[first.toLayer], + this->mTrkParams[iteration].LayerRadii[second.toLayer], + mTimeFrameGPU->getLinkMSAngle(cellTopology.firstLink)) + : -1.f; + const float cellTanLNSigma = this->mTrkParams[iteration].CellDeltaTanLambdaNSigma > 0.f + ? this->mTrkParams[iteration].CellDeltaTanLambdaNSigma + : this->mTrkParams[iteration].NSigmaCut; const int currentLayerTrackletsNum{static_cast(mTimeFrameGPU->getNTracklets()[cellTopology.firstLink])}; if (!currentLayerTrackletsNum || !mTimeFrameGPU->getNTracklets()[cellTopology.secondLink]) { mTimeFrameGPU->getNCells()[cellTopologyId] = 0; @@ -180,7 +246,8 @@ void TrackerTraitsGPU::computeLayerCells(const int iteration) this->mBz, this->mTrkParams[iteration].MaxChi2ClusterAttachment, this->mTrkParams[iteration].CellDeltaTanLambdaSigma, - this->mTrkParams[iteration].NSigmaCut, + cellDeltaPhiCut, + cellTanLNSigma, mTimeFrameGPU->getDeviceLayerxX0(), mTimeFrameGPU->getFrameworkAllocator(), mTimeFrameGPU->getStreams()); @@ -191,6 +258,317 @@ void TrackerTraitsGPU::computeLayerCells(const int iteration) mTimeFrameGPU->syncStreams(false); } +template +void TrackerTraitsGPU::computeVertexCandidates(const int iteration) +{ + const int nCells = mTimeFrameGPU->getNCells()[0]; + if (!nCells) { + mTimeFrameGPU->setNLinesTotal(0); + return; + } + mTimeFrameGPU->createClusterOwnersDevice(); + mTimeFrameGPU->createLinesDevice(nCells); + mTimeFrameGPU->resetClusterOwnersDevice(); + mTimeFrameGPU->syncStreams(false); + + registerClusterOwnershipHandler(mTimeFrameGPU->getDeviceCells()[0], + nCells, + mTimeFrameGPU->getDeviceArrayClusterOwners(), + mTimeFrameGPU->getStream(0)); + + linearizeCellsToLinesHandler(nCells, + mTimeFrameGPU->getDeviceCells()[0], + mTimeFrameGPU->getDeviceArrayClusterOwners(), + mTimeFrameGPU->getDeviceROFramesClusters(1), + mTimeFrameGPU->getNrof(1), + this->mTrkParams[iteration].CellLineSharedClusterCut, + mTimeFrameGPU->getDeviceLines(), + mTimeFrameGPU->getDeviceLineRof(), + mTimeFrameGPU->getDeviceLineClusters(), + mTimeFrameGPU->getDeviceLineSlots(), + mTimeFrameGPU->getBeamX(), + mTimeFrameGPU->getBeamY(), + o2::its::VertexerParamConfig::Instance().maxZPositionAllowed, + o2::its::VertexerParamConfig::Instance().lineMinPt, + mTimeFrameGPU->getDeviceLineZs(), + mTimeFrameGPU->getDeviceLineTimes(), + mTimeFrameGPU->getDeviceLineChi2(), + mTimeFrameGPU->getDeviceLinePt(), + mTimeFrameGPU->getFrameworkAllocator(), + mTimeFrameGPU->getStream(0)); + const unsigned int nLines = mTimeFrameGPU->downloadLinesDevice(); + + sortLinesHandler(nLines, + mTimeFrameGPU->getNrof(1), + mTimeFrameGPU->getLineProjSoA(), + mTimeFrameGPU->getLineProjSortedSoA(), + mTimeFrameGPU->getDeviceLineRof(), + mTimeFrameGPU->getDeviceRofLineOffsets(), + mTimeFrameGPU->getFrameworkAllocator(), + mTimeFrameGPU->getStream(0)); + + const float zWindow = 0.5f * o2::its::VertexerParamConfig::Instance().clusterCut; + scanDensityHandler(static_cast(nLines), + mTimeFrameGPU->getLineProjSortedSoA(), + mTimeFrameGPU->getDeviceRofLineOffsets(), + mTimeFrameGPU->getDeviceLineDensity(), + mTimeFrameGPU->getDeviceLineWin(), + zWindow, + mTimeFrameGPU->getStream(0)); + + const float fineZWindow = o2::its::VertexerParamConfig::Instance().fineZWindow; + const bool doFine = fineZWindow > 0.f && fineZWindow < zWindow; + if (doFine) { + scanDensityHandler(static_cast(nLines), + mTimeFrameGPU->getLineProjSortedSoA(), + mTimeFrameGPU->getDeviceRofLineOffsets(), + mTimeFrameGPU->getDeviceLineDensityFine(), + mTimeFrameGPU->getDeviceLineWinFine(), + fineZWindow, + mTimeFrameGPU->getStream(0)); + } + + findPeaksHandler(nLines, + mTimeFrameGPU->getNrof(1), + mTimeFrameGPU->getLineProjSortedSoA(), + mTimeFrameGPU->getDeviceRofLineOffsets(), + mTimeFrameGPU->getDeviceLineDensity(), + mTimeFrameGPU->getDeviceLineWin(), + mTimeFrameGPU->getDeviceLineIsPeak(), + doFine ? mTimeFrameGPU->getDeviceLineDensityFine() : nullptr, + doFine ? mTimeFrameGPU->getDeviceLineWinFine() : nullptr, + o2::its::VertexerParamConfig::Instance().fineMinDensity, + doFine ? mTimeFrameGPU->getDeviceLineIsPeakFine() : nullptr, + mTimeFrameGPU->getDevicePeakScan(), + mTimeFrameGPU->getDevicePeakLineIdx(), + mTimeFrameGPU->getDevicePeakOffsets(), + mTimeFrameGPU->getFrameworkAllocator(), + mTimeFrameGPU->getStream(0)); + + const auto& vc = o2::its::VertexerParamConfig::Instance(); + const float goodLinePtCut = std::abs(mTimeFrameGPU->getBz()) > 0.01f ? vc.goodLinePtCut : -1.f; + fitPeaksHandler(mTimeFrameGPU->getDeviceNPeaks(), + mTimeFrameGPU->getDevicePeakLineIdx(), + mTimeFrameGPU->getDeviceLineWin(), + mTimeFrameGPU->getLineProjSortedSoA(), + mTimeFrameGPU->getDeviceLines(), + mTimeFrameGPU->getDeviceLineChi2(), + mTimeFrameGPU->getDeviceLinePt(), + vc.goodLineChi2Cut, + goodLinePtCut, + vc.pairCut * vc.pairCut, + vc.nSigmaCut, + vc.clusterContributorsCut, + mTimeFrameGPU->getBeamX(), + mTimeFrameGPU->getBeamY(), + doFine ? mTimeFrameGPU->getDeviceLineIsPeakFine() : nullptr, + vc.fineMaxDrift, + mTimeFrameGPU->getDeviceVertexCands(), + mTimeFrameGPU->getStream(0)); + + const float duplicateZCut = vc.duplicateZCut > 0.f + ? vc.duplicateZCut + : std::max(4.f * vc.pairCut, 0.5f * vc.clusterCut); + dedupVertexCandidatesHandler(mTimeFrameGPU->getDeviceNPeaks(), + mTimeFrameGPU->getDevicePeakLineIdx(), + mTimeFrameGPU->getDevicePeakOffsets(), + mTimeFrameGPU->getLineProjSortedSoA(), + duplicateZCut, + vc.duplicateZScale, + mTimeFrameGPU->getDeviceVertexCands(), + mTimeFrameGPU->getStream(0)); + + const bool withMC = mTimeFrameGPU->hasMCinformation() && this->mTrkParams[iteration].CreateArtefactLabels; + mTimeFrameGPU->downloadVertexCandsDevice(); // sets nPeaks + host candidate/peak-offset mirrors + if (withMC) { + scanMemberOffsetsHandler(mTimeFrameGPU->getDeviceVertexCands(), + mTimeFrameGPU->getDeviceMemberOffsets(), + static_cast(nLines), + mTimeFrameGPU->getFrameworkAllocator(), + mTimeFrameGPU->getStream(0)); + mTimeFrameGPU->downloadMemberOffsetsDevice(); // sets getNMembers() + mTimeFrameGPU->createMemberLinesMCDevice(mTimeFrameGPU->getNMembers()); + collectLinesForMCHandler(mTimeFrameGPU->getDeviceNPeaks(), + mTimeFrameGPU->getDevicePeakLineIdx(), + mTimeFrameGPU->getDeviceLineWin(), + mTimeFrameGPU->getLineProjSortedSoA(), + mTimeFrameGPU->getDeviceLines(), + vc.pairCut * vc.pairCut, + mTimeFrameGPU->getDeviceVertexCands(), + mTimeFrameGPU->getDeviceMemberOffsets(), + mTimeFrameGPU->getDeviceMemberLines(), + mTimeFrameGPU->getStream(0)); + mTimeFrameGPU->downloadMemberLinesDevice(); + } + const auto& lines = mTimeFrameGPU->getHostLines(); + const auto& lineRof = mTimeFrameGPU->getHostLineRof(); + const auto& lineClusters = mTimeFrameGPU->getHostLineClusters(); + const int nRofs = mTimeFrameGPU->getNrof(1); + if (withMC) { + mTimeFrameGPU->getLineLabelFlat().assign(nLines, o2::MCCompLabel()); // global-indexed, for the vertex-label vote + } + auto lineLabel = [&](const int* cl) -> o2::MCCompLabel { + const auto l0 = mTimeFrameGPU->getClusterLabels(0, cl[0]); + const auto l1 = mTimeFrameGPU->getClusterLabels(1, cl[1]); + const auto l2 = mTimeFrameGPU->getClusterLabels(2, cl[2]); + for (const auto& a : l0) { + if (!a.isValid()) { + continue; + } + bool in1{false}, in2{false}; + for (const auto& b : l1) { + if (b == a) { + in1 = true; + break; + } + } + for (const auto& c : l2) { + if (c == a) { + in2 = true; + break; + } + } + if (in1 && in2) { + return a; + } + } + return o2::MCCompLabel(); + }; + for (unsigned int i{0}; i < nLines; ++i) { + const int rof = lineRof[i]; + if (rof < 0 || rof >= nRofs) { + LOGP(fatal, "ITS GPU linearizer: line {} carries out-of-range ROF {} (nRofs={}).", i, rof, nRofs); + } + const auto& l = lines[i]; + mTimeFrameGPU->getLines(rof).emplace_back(std::array{l.originPoint[0], l.originPoint[1], l.originPoint[2]}, + std::array{l.cosinesDirector[0], l.cosinesDirector[1], l.cosinesDirector[2]}, + l.mTime); + if (withMC) { + const auto lbl = lineLabel(&lineClusters[3 * i]); + mTimeFrameGPU->getLinesLabel(rof).emplace_back(lbl); + mTimeFrameGPU->getLineLabelFlat()[i] = lbl; // same label, global-indexed for the vote + } + } + mTimeFrameGPU->setNLinesTotal(nLines); +} + +template +void TrackerTraitsGPU::computeVertices(const int iteration) +{ + const int nRofs = mTimeFrameGPU->getNrof(1); + const bool withMC = mTimeFrameGPU->hasMCinformation() && this->mTrkParams[iteration].CreateArtefactLabels; + const auto& vc = o2::its::VertexerParamConfig::Instance(); + const int suppressLowMultDebris = vc.suppressLowMultDebris; + const bool skipHighMultRofs = this->mTrkParams[iteration].PassFlags[IterationStep::SkipROFsAboveThreshold]; + + const auto& cands = mTimeFrameGPU->getHostVertexCands(); // per compacted peak slot [0, nPeaks) + const auto& peakOffsets = mTimeFrameGPU->getHostPeakOffsets(); // nRofs+1; per-ROF peak slices + const auto& memberOffsets = mTimeFrameGPU->getHostMemberOffsets(); // nPeaks+1; per-survivor member slices + const auto& memberLines = mTimeFrameGPU->getHostMemberLines(); // global line indices, grouped by survivor + const auto& lineLabelFlat = mTimeFrameGPU->getLineLabelFlat(); // nLines, global-indexed line labels + const int nLab = static_cast(lineLabelFlat.size()); + + std::vector> rofVertices(nRofs); + std::vector> rofLabels(nRofs); + const float goodSig = vc.goodContributorsSignificance; + + for (int rofId = 0; rofId < nRofs; ++rofId) { + if (skipHighMultRofs && + static_cast(mTimeFrameGPU->getROFVertexLookupTableView().getVertices(1, rofId).getEntries()) > vc.vertPerRofThreshold) { + continue; + } + // Survivors of this ROF, sorted by contributor count desc + std::vector accepted; + for (int p = peakOffsets[rofId]; p < peakOffsets[rofId + 1]; ++p) { + if (cands[p].keep) { + accepted.push_back(p); + } + } + std::sort(accepted.begin(), accepted.end(), [&](const int a, const int b) { return cands[a].size > cands[b].size; }); + + double rofLoad = 0.; + if (goodSig > 0.f) { // compute the number of contributors in this ROF + for (int p = peakOffsets[rofId]; p < peakOffsets[rofId + 1]; ++p) { + if (cands[p].ok && !cands[p].fine) { + rofLoad += cands[p].size; + } + } + } + const float sigThreshold = goodSig > 0.f ? goodSig * std::sqrt(static_cast(std::max(rofLoad, 1.))) : 0.f; + + for (const int p : accepted) { + const auto& c = cands[p]; + if (!rofVertices[rofId].empty()) { + if (goodSig > 0.f) { + if (c.nGood <= sigThreshold) { + continue; + } + } else if (c.size < suppressLowMultDebris) { + continue; + } + } + const float pos[3] = {c.x, c.y, c.z}; + Vertex vertex(pos, c.rms2, static_cast(c.size), c.avgDist2); + vertex.setTimeStamp(c.time); + rofVertices[rofId].push_back(vertex); + if (withMC) { + std::vector labels; + labels.reserve(memberOffsets[p + 1] - memberOffsets[p]); + for (int e = memberOffsets[p]; e < memberOffsets[p + 1]; ++e) { + const int gi = memberLines[e]; + if (gi < 0 || gi >= nLab) { + LOGP(error, "[seedDbg] GPU emit: member line {} out of range (nLines={}) rof={}", gi, nLab, rofId); + continue; + } + labels.push_back(lineLabelFlat[gi]); + } + rofLabels[rofId].push_back(computeMainGPU(labels)); + } + } + } + + for (int rofId = 0; rofId < nRofs; ++rofId) { + for (auto& vertex : rofVertices[rofId]) { + mTimeFrameGPU->addPrimaryVertex(vertex); + } + if (withMC) { + for (auto& label : rofLabels[rofId]) { + mTimeFrameGPU->addPrimaryVertexLabel(label); + } + } + } + + auto& pvs = mTimeFrameGPU->getPrimaryVertices(); + std::vector indices(pvs.size()); + std::iota(indices.begin(), indices.end(), 0); + std::sort(indices.begin(), indices.end(), [&pvs](const size_t i, const size_t j) { + const auto aLower = pvs[i].getTimeStamp().lower(); + const auto bLower = pvs[j].getTimeStamp().lower(); + if (aLower != bLower) { + return aLower < bLower; + } + return pvs[i].getNContributors() > pvs[j].getNContributors(); + }); + std::decay_t sortedVtx(pvs.get_allocator()); + sortedVtx.reserve(pvs.size()); + for (const size_t idx : indices) { + sortedVtx.push_back(pvs[idx]); + } + pvs.swap(sortedVtx); + if (withMC) { + auto& mc = mTimeFrameGPU->getPrimaryVerticesLabels(); + std::decay_t sortedMC(mc.get_allocator()); + sortedMC.reserve(mc.size()); + for (const size_t idx : indices) { + sortedMC.push_back(mc[idx]); + } + mc.swap(sortedMC); + } + mTimeFrameGPU->updateROFVertexLookupTable(); + + mTimeFrameGPU->popMemoryStack(iteration); // frees the whole seeding-pass stack frame +} + template void TrackerTraitsGPU::findCellsNeighbours(const int iteration) { diff --git a/Detectors/ITSMFT/ITS/tracking/GPU/cuda/TrackingKernels.cu b/Detectors/ITSMFT/ITS/tracking/GPU/cuda/TrackingKernels.cu index a179fd0962dd3..a7ff623d3fa55 100644 --- a/Detectors/ITSMFT/ITS/tracking/GPU/cuda/TrackingKernels.cu +++ b/Detectors/ITSMFT/ITS/tracking/GPU/cuda/TrackingKernels.cu @@ -14,6 +14,7 @@ #include #include #include +#include #include #include @@ -22,11 +23,19 @@ #include #include #include +#include #include #include #include #include #include +#include +#include +#include +#include +#include +#include +#include #include "DataFormatsITS/TrackITS.h" #include "ITStracking/Constants.h" @@ -43,6 +52,7 @@ #include "ITStrackingGPU/TrackingKernels.h" #include "ITStrackingGPU/Utils.h" #include "MathUtils/Utils.h" +#include "ITStrackingGPU/ClusterLinesGPU.h" #include "utils/strtag.h" // O2 track model @@ -272,6 +282,7 @@ GPUg() void __launch_bounds__(GPUThreads, MinBlocks.computeLayerCells) computeLa int* outputCounter, const int outputCapacity, const float cellDeltaTanLambdaSigma, + const float cellDeltaPhiCut, const float nSigmaCut) { const auto cellTopology = topology.getCell(cellTopologyId); @@ -291,6 +302,15 @@ GPUg() void __launch_bounds__(GPUThreads, MinBlocks.computeLayerCells) computeLa if (!currentTracklet.getTimeStamp().isCompatible(nextTracklet.getTimeStamp())) { continue; } + if (cellDeltaPhiCut > 0.f) { + float deltaPhi{o2::gpu::CAMath::Abs(currentTracklet.phi - nextTracklet.phi)}; + if (deltaPhi > o2::constants::math::PI) { + deltaPhi = o2::constants::math::TwoPI - deltaPhi; + } + if (deltaPhi > cellDeltaPhiCut) { + continue; + } + } const float deltaTanLambda{o2::gpu::CAMath::Abs(currentTracklet.tanLambda - nextTracklet.tanLambda)}; if (deltaTanLambda / cellDeltaTanLambdaSigma < nSigmaCut) { if constexpr (Emit) { @@ -397,6 +417,7 @@ GPUg() void __launch_bounds__(GPUThreads, MinBlocks.computeLayerTracklets) compu const typename ROFOverlapTable::View rofOverlaps, const typename ROFVertexLookupTable::View vertexLUT, const Vertex* vertices, + const bool vtxMode, const int vertexId, const Cluster** clusters, const int** ROFClusters, @@ -409,8 +430,8 @@ GPUg() void __launch_bounds__(GPUThreads, MinBlocks.computeLayerTracklets) compu const float NSigmaCut, const float phiCut, const float resolutionPV, - const float minR, - const float maxR, + const float* minRs, + const float* maxRs, const float positionResolution, const float meanDeltaR, const float MSAngle) @@ -418,6 +439,8 @@ GPUg() void __launch_bounds__(GPUThreads, MinBlocks.computeLayerTracklets) compu const auto link = topology.getLink(linkId); const int fromLayer = link.fromLayer; const int toLayer = link.toLayer; + const float minR = minRs[toLayer]; + const float maxR = maxRs[toLayer]; const int phiBins{utils->getNphiBins()}; const int zBins{utils->getNzBins()}; const int tableSize{phiBins * zBins + 1}; @@ -447,8 +470,14 @@ GPUg() void __launch_bounds__(GPUThreads, MinBlocks.computeLayerTracklets) compu continue; } - const auto& pvs = vertexLUT.getVertices(fromLayer, pivotROF); - auto primaryVertices = gpuSpan(&vertices[pvs.getFirstEntry()], pvs.getEntries()); + // The diamond is a single PV-independent vertex: the lookup table is not consulted at all, since during the seeding-vertex pass it holds no vertices yet + gpuSpan primaryVertices; + if (vtxMode) { + primaryVertices = gpuSpan(vertices, 1); + } else { + const auto& pvs = vertexLUT.getVertices(fromLayer, pivotROF); + primaryVertices = gpuSpan(&vertices[pvs.getFirstEntry()], pvs.getEntries()); + } if (primaryVertices.empty()) { continue; } @@ -472,7 +501,7 @@ GPUg() void __launch_bounds__(GPUThreads, MinBlocks.computeLayerTracklets) compu const float inverseR0{1.f / currentCluster.radius}; for (int iV{startVtx}; iV < endVtx; ++iV) { auto& primaryVertex{primaryVertices[iV]}; - if (!vertexLUT.isVertexCompatible(fromLayer, pivotROF, primaryVertex)) { + if (!vtxMode && !vertexLUT.isVertexCompatible(fromLayer, pivotROF, primaryVertex)) { continue; } if (primaryVertex.isFlagSet(Vertex::Flags::UPCMode) != selectUPCVertices) { @@ -504,7 +533,7 @@ GPUg() void __launch_bounds__(GPUThreads, MinBlocks.computeLayerTracklets) compu continue; } const auto ts = rofOverlaps.getTimeStamp(fromLayer, pivotROF, toLayer, targetROF); - if (!ts.isCompatible(primaryVertex.getTimeStamp())) { + if (!vtxMode && !ts.isCompatible(primaryVertex.getTimeStamp())) { continue; } for (int iPhiCount{0}; iPhiCount < phiBinsNum; iPhiCount++) { @@ -760,6 +789,446 @@ GPUg() void __launch_bounds__(GPUThreads, (std::is_same_v } } +GPUg() void vertexingRegisterCellClustersOwnership( + const CellSeed* cells, + const int nCells, + unsigned long long** clusterOwners) +{ + for (int k = blockIdx.x * blockDim.x + threadIdx.x; k < nCells; k += blockDim.x * gridDim.x) { + const CellSeed& cell = cells[k]; + if (o2::gpu::CAMath::Abs(cell.getQ2Pt()) < o2::constants::math::Almost0 || + o2::gpu::CAMath::Abs(cell.getSnp()) > o2::constants::math::Almost1) { + continue; + } + const float pt = cell.getPt(); + const float rank = pt > 1.e-6f ? 1.f / pt : 1.e9f; + const unsigned long long key = (static_cast(__float_as_uint(rank)) << 32) | static_cast(k); + o2::gpu::GPUCommonMath::AtomicMin(&clusterOwners[0][cell.getFirstClusterIndex()], key); + o2::gpu::GPUCommonMath::AtomicMin(&clusterOwners[1][cell.getSecondClusterIndex()], key); + o2::gpu::GPUCommonMath::AtomicMin(&clusterOwners[2][cell.getThirdClusterIndex()], key); + } +} + +GPUdi() int clusterROF(const int* rofArr, const int nRofs, const int clusterIdx) +{ + const int key = clusterIdx + 1; + int lo = 0, hi = nRofs + 1; + while (lo < hi) { + const int mid = (lo + hi) >> 1; + if (rofArr[mid] < key) { + lo = mid + 1; + } else { + hi = mid; + } + } + return lo - 1; +} + +template +GPUg() void dedupCellsKernel( + const int nCells, + const CellSeed* cells, + const unsigned long long* const* clusterOwners, + const int ownedClustersCut, + const float beamX, + const float beamY, + const float maxZ, + const float minPt, + int* cellAccepted) +{ + for (int i = blockIdx.x * blockDim.x + threadIdx.x; i < nCells; i += blockDim.x * gridDim.x) { + const CellSeed& cell = cells[i]; + std::array origin, direction; + if (!cell.getPxPyPzGlo(direction)) { + cellAccepted[i] = 0; + continue; + } + const bool owned0 = static_cast(clusterOwners[0][cell.getFirstClusterIndex()]) == static_cast(i); + const bool owned1 = static_cast(clusterOwners[1][cell.getSecondClusterIndex()]) == static_cast(i); + const bool owned2 = static_cast(clusterOwners[2][cell.getThirdClusterIndex()]) == static_cast(i); + const bool keepCell = (static_cast(owned0) + static_cast(owned1) + static_cast(owned2)) >= 3 - ownedClustersCut; + cell.getXYZGlo(origin); + const float dx = origin[0] - beamX; + const float dy = origin[1] - beamY; + const float den = direction[0] * direction[0] + direction[1] * direction[1]; + const bool projOk = den >= constants::Tolerance && o2::gpu::CAMath::Abs(origin[2] - (dx * direction[0] + dy * direction[1]) / den * direction[2]) < maxZ; + const bool ptOk = minPt <= 0.f || cell.getPt() >= minPt; + cellAccepted[i] = keepCell && projOk && ptOk ? 1 : 0; + } +} + +template +GPUg() void linearizeCellsKernel( + const int nCells, + const CellSeed* cells, + const int* rofFramesClustersL1, // layer-1 ROF boundaries, size nRofsL1 + 1 + const int nRofsL1, + const int* lineSlots, // exclusive-scanned accept flags, size nCells + 1 + GPULine* lines, + int* lineRof, + const float beamX, + const float beamY, + float* lineZs, + gpu::LineTime* lineTimes, + int* lineClusters, // 3 per line (L0,L1,L2 cluster ids), for the host-side MC label derivation + float* lineChi2, + float* linePt) +{ + for (int i = blockIdx.x * blockDim.x + threadIdx.x; i < nCells; i += blockDim.x * gridDim.x) { + const int slot = lineSlots[i]; + if (slot == lineSlots[i + 1]) { + continue; + } + const CellSeed& cell = cells[i]; + std::array origin, direction; + cell.getXYZGlo(origin); + cell.getPxPyPzGlo(direction); + lines[slot] = GPULine{origin.data(), direction.data(), cell.getTimeStamp()}; + lineRof[slot] = clusterROF(rofFramesClustersL1, nRofsL1, cell.getSecondClusterIndex()); + const float dx = origin[0] - beamX; + const float dy = origin[1] - beamY; + const float den = direction[0] * direction[0] + direction[1] * direction[1]; + const float s0 = -(dx * direction[0] + dy * direction[1]) / den; + lineZs[slot] = origin[2] + s0 * direction[2]; + const auto sym = cell.getTimeStamp().makeSymmetrical(); + lineTimes[slot] = LineTime{sym.getTimeStamp(), sym.getTimeStampError()}; + lineClusters[3 * slot + 0] = cell.getFirstClusterIndex(); + lineClusters[3 * slot + 1] = cell.getSecondClusterIndex(); + lineClusters[3 * slot + 2] = cell.getThirdClusterIndex(); + lineChi2[slot] = cell.getChi2(); + linePt[slot] = cell.getPt(); + } +} + +template +GPUg() void gatherSortedLinesKernel(const int nLines, LineProjSoA lineProj, LineProjSoA lineProjSorted) +{ + const int* sortedIdx = lineProj.idx; + for (int i = blockIdx.x * blockDim.x + threadIdx.x; i < nLines; i += blockDim.x * gridDim.x) { + lineProjSorted.z[i] = lineProj.z[sortedIdx[i]]; + lineProjSorted.t[i] = lineProj.t[sortedIdx[i]]; + lineProjSorted.rof[i] = lineProj.rof[sortedIdx[i]]; + } +} + +template +GPUg() void scanDensityKernel(int* zDensity, LineWindow* win, const int nLines, const int* offsets, const LineProjSoA lineProjSorted, const float zWindow) +{ + const float* z = lineProjSorted.z; + const LineTime* t = lineProjSorted.t; + const int* rof = lineProjSorted.rof; + for (int iLine = blockIdx.x * blockDim.x + threadIdx.x; iLine < nLines; iLine += blockDim.x * gridDim.x) { + const int rofId = rof[iLine]; + const int rofOffset = offsets[rofId]; + const int nextRofOffset = offsets[rofId + 1]; + const float zk = z[iLine]; + const int lo = deviceLowerBound(z, rofOffset, nextRofOffset, zk - zWindow); // first line with z >= zk - zWindow + const int hi = deviceUpperBound(z, rofOffset, nextRofOffset, zk + zWindow); // first line with z > zk + zWindow + win[iLine] = LineWindow{lo, hi}; + const LineTime ti = t[iLine]; + int count = 0; + for (int j = lo; j < hi; ++j) { + const LineTime tj = t[j]; + if (o2::gpu::GPUCommonMath::Abs(ti.tc - tj.tc) <= (ti.th + tj.th)) { // count only if time compatible (includes self) + ++count; + } + } + zDensity[iLine] = count; + } +} + +template +GPUg() void fitPeaksKernel(const int* nPeaksDevice, + const int* peakLineIdx, + const gpu::LineWindow* win, + const LineProjSoA lineProjSorted, + const GPULine* lines, + const float* lineChi2, // global-indexed, same indexing as lines[] + const float* linePt, // idem + const float goodLineChi2Cut, // a contributor counts towards nGood only if its own + const float goodLinePtCut, // cell passes both; <= 0 disables that half + const float pairCut2, + const float nSigmaCut, + const int minContributors, + const float beamX, + const float beamY, + const uint8_t* isZPeakFine, // null when the fine pass is off + const float fineMaxDrift, // <= 0 disables; see VertexerParamConfig::fineMaxDrift + VertexCand* cands) +{ + const int nPeaks = *nPeaksDevice; + const LineTime* t = lineProjSorted.t; + const int* idx = lineProjSorted.idx; + for (int p = blockIdx.x * blockDim.x + threadIdx.x; p < nPeaks; p += blockDim.x * gridDim.x) { + cands[p].ok = 0; + cands[p].nGood = 0; + const int k = peakLineIdx[p]; + cands[p].fine = isZPeakFine != nullptr ? isZPeakFine[k] : 0; + const LineTime tk = t[k]; + const LineWindow wk = win[k]; + + GPUClusterLinesFit seed; + int nMembers = 0; + for (int j = wk.lo; j < wk.hi; ++j) { + const LineTime tj = t[j]; + if (o2::gpu::GPUCommonMath::Abs(tk.tc - tj.tc) <= tk.th + tj.th) { + seed.add(lines[idx[j]]); + ++nMembers; + } + } + float seedVertex[3]; + if (nMembers < 2 || !seed.solve(seedVertex)) { + continue; + } + + GPUClusterLinesFit fit; + int nKept = 0; + int nGood = 0; + for (int j = wk.lo; j < wk.hi; ++j) { + const LineTime tj = t[j]; + if (o2::gpu::GPUCommonMath::Abs(tk.tc - tj.tc) <= tk.th + tj.th) { + const GPULine& line = lines[idx[j]]; + if (GPULine::getDistance2FromPoint(line, seedVertex) < pairCut2) { + fit.add(line); + const float c = lineChi2[idx[j]]; // the kept set is exactly what collectLinesForMCKernel re-walks + const float pt = linePt[idx[j]]; + const bool okChi2 = (goodLineChi2Cut <= 0.f || c <= goodLineChi2Cut); + const bool okPt = (goodLinePtCut <= 0.f || pt >= goodLinePtCut); + nGood += okChi2 && okPt; + ++nKept; + } + } + } + float vertex[3]; + if (nKept < 2 || !fit.solve(vertex)) { + continue; + } + cands[p].seed[0] = seedVertex[0]; + cands[p].seed[1] = seedVertex[1]; + cands[p].seed[2] = seedVertex[2]; + const float bd2 = (beamX - vertex[0]) * (beamX - vertex[0]) + (beamY - vertex[1]) * (beamY - vertex[1]); + if (nKept < minContributors || !(bd2 < nSigmaCut)) { + continue; + } + if (fineMaxDrift > 0.f && cands[p].fine && + o2::gpu::GPUCommonMath::Abs(vertex[2] - seedVertex[2]) > fineMaxDrift) { + continue; + } + + for (int j = wk.lo; j < wk.hi; ++j) { + const LineTime tj = t[j]; + if (o2::gpu::GPUCommonMath::Abs(tk.tc - tj.tc) <= tk.th + tj.th) { + const GPULine& line = lines[idx[j]]; + if (GPULine::getDistance2FromPoint(line, seedVertex) < pairCut2) { + fit.addResidual(line, vertex); + } + } + } + + cands[p].x = vertex[0]; + cands[p].y = vertex[1]; + cands[p].z = vertex[2]; + for (int i = 0; i < 6; ++i) { + cands[p].rms2[i] = fit.getRMS2()[i]; + } + cands[p].avgDist2 = fit.getAvgDistance2(); + cands[p].nGood = nGood; + cands[p].time = fit.getTimeStamp(); + cands[p].size = nKept; + cands[p].ok = 1; + } +} + +// Strict local maximum of the density over a line's own z-window, ties broken by smaller z +GPUdi() bool isDensityPeak(const int* density, const float* z, const LineWindow w, const int iLine) +{ + const int di = density[iLine]; + const float zi = z[iLine]; + for (int j = w.lo; j < w.hi; ++j) { + const int dj = density[j]; + if (dj > di || (dj == di && z[j] < zi)) { + return false; + } + } + return true; +} + +template +GPUg() void findPeaksKernel(const int* zDensity, const LineWindow* win, const int nLines, const LineProjSoA lineProjSorted, uint8_t* isZPeak, + const int* zDensityFine, const LineWindow* winFine, + const int fineMinDensity, uint8_t* isZPeakFine) +{ + const float* z = lineProjSorted.z; + for (int iLine = blockIdx.x * blockDim.x + threadIdx.x; iLine < nLines; iLine += blockDim.x * gridDim.x) { + uint8_t peak = zDensity[iLine] >= 2 && isDensityPeak(zDensity, z, win[iLine], iLine); + + // fine pass: if the coarse pass did not find a peak, check if the fine density is above threshold and is a peak + uint8_t fine = 0; + if (!peak && zDensityFine != nullptr) { + fine = zDensityFine[iLine] >= fineMinDensity && isDensityPeak(zDensityFine, z, winFine[iLine], iLine); + peak = fine; + } + isZPeak[iLine] = peak; + if (isZPeakFine != nullptr) { + isZPeakFine[iLine] = fine; + } + } +} + +template +GPUg() void dedupVertexCandidatesKernel(const int* nPeaksDevice, + const int* peakLineIdx, + const int* peakOffsets, + const LineProjSoA lineProjSorted, + const float duplicateZCut, + const float duplicateZScale, + VertexCand* cands) +{ + const int nPeaks = *nPeaksDevice; + for (int p = blockIdx.x * blockDim.x + threadIdx.x; p < nPeaks; p += blockDim.x * gridDim.x) { + cands[p].keep = 0; // every visited slot must be written: this array is never memset + if (!cands[p].ok) { + continue; + } + const int r = lineProjSorted.rof[peakLineIdx[p]]; + const float zp = cands[p].z; + const int sp = cands[p].size; + float radius = duplicateZCut; + if (duplicateZScale > 0.f && sp > 0) { + radius = duplicateZScale / o2::gpu::GPUCommonMath::Sqrt((float)sp); + } + const auto tp = cands[p].time; + uint8_t survive = 1; + for (int q = peakOffsets[r]; q < peakOffsets[r + 1] && survive; ++q) { + if (q == p || !cands[q].ok) { + continue; + } + if (!tp.isCompatible(cands[q].time)) { + continue; + } + if (o2::gpu::GPUCommonMath::Abs(zp - cands[q].z) >= radius) { + continue; + } + const int sq = cands[q].size; + if (sq > sp || (sq == sp && q < p)) { + survive = 0; + } + } + cands[p].keep = survive; + } +} + +template +GPUg() void emitKeysForClusterSortingKernel(const Cluster* unsorted, + const int* clusterOffsets, // this layer, size nRofs+1 + const IndexTableUtils* utils, + const typename ROFMaskTable::View rofMask, + float beamX, float beamY, + int zBins, int phiBins, int nRofs, int iLayer, + float* minRadiusLayer, float* maxRadiusLayer, + int* keys) +{ + const int numBins = zBins * phiBins; + for (int iROF = blockIdx.x; iROF < nRofs; iROF += gridDim.x) { + const bool enabled = rofMask.isROFEnabled(iLayer, iROF); + const int start = clusterOffsets[iROF]; + const int n = clusterOffsets[iROF + 1] - start; + for (int i = threadIdx.x; i < n; i += blockDim.x) { + const Cluster& c = unsorted[start + i]; + const float x = c.xCoordinate - beamX, y = c.yCoordinate - beamY; + const float phi = math_utils::computePhi(x, y); + int zBin = utils->getZBinIndex(iLayer, c.zCoordinate); + zBin = o2::gpu::GPUCommonMath::Max(0, o2::gpu::GPUCommonMath::Min(zBin, zBins - 1)); // TODO: count bogus (clamped) if hasBogusClusters() is ever needed + const int bin = utils->getBinIndex(zBin, utils->getPhiBinIndex(phi)); + if (enabled) { + const float r = math_utils::hypot(x, y); + o2::gpu::GPUCommonMath::AtomicMin(&minRadiusLayer[iLayer], r); + o2::gpu::GPUCommonMath::AtomicMax(&maxRadiusLayer[iLayer], r); + } + keys[start + i] = iROF * numBins + bin; + } + } +} + +template +GPUg() void gatherSortedClustersKernel(const Cluster* unsorted, + Cluster* sorted, + const int* perm, + const IndexTableUtils* utils, + float beamX, float beamY, + int zBins, int nClustersLayer, int iLayer) +{ + for (int j = blockIdx.x * blockDim.x + threadIdx.x; j < nClustersLayer; j += blockDim.x * gridDim.x) { + Cluster c = unsorted[perm[j]]; + const float x = c.xCoordinate - beamX, y = c.yCoordinate - beamY; + const float phi = math_utils::computePhi(x, y); + int zBin = utils->getZBinIndex(iLayer, c.zCoordinate); + zBin = o2::gpu::GPUCommonMath::Max(0, o2::gpu::GPUCommonMath::Min(zBin, zBins - 1)); + c.phi = phi; + c.radius = math_utils::hypot(x, y); + c.indexTableBinIndex = utils->getBinIndex(zBin, utils->getPhiBinIndex(phi)); + sorted[j] = c; + } +} + +template +GPUg() void buildClusterIndexTableKernel(const int* sortedKeys, + const int* clusterOffsets, // this layer, size nRofs+1 + int* indexTable, // output, size nRofs*(numBins+1) + int numBins, int nRofs) +{ + const int stride = numBins + 1; + for (int iROF = blockIdx.x; iROF < nRofs; iROF += gridDim.x) { + const int rofStart = clusterOffsets[iROF]; + const int rofEnd = clusterOffsets[iROF + 1]; + int* base = indexTable + iROF * stride; + for (int b = threadIdx.x; b <= numBins; b += blockDim.x) { + const int keyB = iROF * numBins + b; + base[b] = deviceLowerBound(sortedKeys, rofStart, rofEnd, keyB) - rofStart; // ROF-local + } + } +} + +struct MemberCount { + const VertexCand* cands; + GPUhdi() int operator()(const int p) const { return cands[p].keep ? cands[p].size : 0; } +}; + +template +GPUg() void collectLinesForMCKernel(const int* nPeaksDevice, + const int* peakLineIdx, + const gpu::LineWindow* win, + const LineProjSoA lineProjSorted, + const GPULine* lines, + const float pairCut2, + const VertexCand* cands, + const int* memberOffsets, + int* memberLines) +{ + const int nPeaks = *nPeaksDevice; + const LineTime* t = lineProjSorted.t; + const int* idx = lineProjSorted.idx; + for (int p = blockIdx.x * blockDim.x + threadIdx.x; p < nPeaks; p += blockDim.x * gridDim.x) { + if (!cands[p].keep) { + continue; + } + const int k = peakLineIdx[p]; + const LineTime tk = t[k]; + const LineWindow wk = win[k]; + const float* seedVertex = cands[p].seed; + int localIdx = 0; + for (int j = wk.lo; j < wk.hi; ++j) { + const LineTime tj = t[j]; + if (o2::gpu::GPUCommonMath::Abs(tk.tc - tj.tc) <= tk.th + tj.th) { + const GPULine& line = lines[idx[j]]; + if (GPULine::getDistance2FromPoint(line, seedVertex) < pairCut2) { + memberLines[memberOffsets[p] + localIdx++] = idx[j]; + } + } + } + } +} + } // namespace gpu template @@ -772,6 +1241,7 @@ int TrackingKernels::computeTrackletsInROFsHandler(const IndexTableUtil const typename ROFVertexLookupTable::View& vertexLUT, const int vertexId, const Vertex* vertices, + const bool vtxMode, const Cluster** clusters, const std::vector& nClusters, const int** ROFClusters, @@ -787,8 +1257,8 @@ int TrackingKernels::computeTrackletsInROFsHandler(const IndexTableUtil const typename TrackingTopology::View topology, bounded_vector& linkPhiCuts, const float resolutionPV, - std::array& minRs, - std::array& maxRs, + const float* minRs, + const float* maxRs, bounded_vector& resolutions, std::vector& radii, bounded_vector& linkMSAngles, @@ -806,6 +1276,7 @@ int TrackingKernels::computeTrackletsInROFsHandler(const IndexTableUtil rofOverlaps, vertexLUT, vertices, + vtxMode, vertexId, clusters, ROFClusters, @@ -818,8 +1289,8 @@ int TrackingKernels::computeTrackletsInROFsHandler(const IndexTableUtil NSigmaCut, linkPhiCuts[linkId], resolutionPV, - minRs[toLayer], - maxRs[toLayer], + minRs, + maxRs, resolutions[fromLayer], radii[toLayer] - radii[fromLayer], linkMSAngles[linkId]); @@ -872,6 +1343,7 @@ int TrackingKernels::computeCellsHandler( const float bz, const float maxChi2ClusterAttachment, const float cellDeltaTanLambdaSigma, + const float cellDeltaPhiCut, const float nSigmaCut, const float* layerxX0, o2::its::ExternalAllocator* alloc, @@ -890,7 +1362,7 @@ int TrackingKernels::computeCellsHandler( gpu::computeLayerCellCandidatesKernel<<>>( tracklets, trackletsLUT, nTracklets, cellTopologyId, topology, sortedClusters, unsortedClusters, tfInfo, layerxX0, bz, - nullptr, nullptr, outputCounter, 0, cellDeltaTanLambdaSigma, nSigmaCut); + nullptr, nullptr, outputCounter, 0, cellDeltaTanLambdaSigma, cellDeltaPhiCut, nSigmaCut); int nCandidates = 0; GPUChkErrS(cudaMemcpyAsync(&nCandidates, outputCounter, sizeof(int), cudaMemcpyDeviceToHost, stream.get())); stream.sync(); @@ -910,7 +1382,7 @@ int TrackingKernels::computeCellsHandler( tracklets, trackletsLUT, nTracklets, cellTopologyId, topology, sortedClusters, unsortedClusters, tfInfo, layerxX0, bz, thrust::raw_pointer_cast(candidates), thrust::raw_pointer_cast(candidateKeys), - outputCounter, nCandidates, cellDeltaTanLambdaSigma, nSigmaCut); + outputCounter, nCandidates, cellDeltaTanLambdaSigma, cellDeltaPhiCut, nSigmaCut); // order the candidates by momentum before fitting them, so that the ELoss iteration count inside is uniform { @@ -1027,6 +1499,253 @@ void TrackingKernels::computeCellNeighboursHandler(CellSeed** cellsLaye alloc->popTagOffStack(CandidateTag); } +template +void sortClustersHandler(const Cluster* unsorted, // this layer (resident unsorted) + Cluster* sorted, // this layer (output) + const int* clusterOffsets, // this layer ROF boundaries, size nRofs+1 + int* indexTable, // this layer (output), size nRofs*(zBins*phiBins+1) + const IndexTableUtils* utils, + const typename ROFMaskTable::View& rofMask, + float beamX, float beamY, + int zBins, int phiBins, int nRofs, int nClustersLayer, int iLayer, + float* minRadiusLayer, float* maxRadiusLayer, // per-layer device arrays + int* keys, // scratch, size nClustersLayer + int* perm, // scratch, size nClustersLayer + o2::its::ExternalAllocator* alloc, + gpu::Stream& stream) +{ + if (nClustersLayer == 0) { + return; + } + const int numBins = zBins * phiBins; + auto policy = THRUST_NAMESPACE::par_nosync(gpu::TypedAllocator(alloc)).on(stream.get()); + thrust::fill_n(policy, minRadiusLayer + iLayer, 1, std::numeric_limits::max()); + thrust::fill_n(policy, maxRadiusLayer + iLayer, 1, std::numeric_limits::min()); + + gpu::emitKeysForClusterSortingKernel<<>>( + unsorted, clusterOffsets, utils, rofMask, beamX, beamY, zBins, phiBins, nRofs, iLayer, + minRadiusLayer, maxRadiusLayer, keys); + + thrust::sequence(policy, perm, perm + nClustersLayer); + thrust::stable_sort_by_key(policy, keys, keys + nClustersLayer, perm); + + gpu::gatherSortedClustersKernel<<>>( + unsorted, sorted, perm, utils, beamX, beamY, zBins, nClustersLayer, iLayer); + + gpu::buildClusterIndexTableKernel<<>>( + keys, clusterOffsets, indexTable, numBins, nRofs); +} + +template +void registerClusterOwnershipHandler(const CellSeed* cells, + const int nCells, + unsigned long long** clusterOwnersDeviceArray, + gpu::Stream& stream) +{ + + gpu::vertexingRegisterCellClustersOwnership<<>>( + cells, + nCells, + clusterOwnersDeviceArray); +} + +template +void linearizeCellsToLinesHandler(const int nCells, + const CellSeed* cells, + const unsigned long long* const* clusterOwners, + const int* rofFramesClustersL1, + const int nRofsL1, + const int ownedClustersCut, + gpu::GPULine* lines, + int* lineRof, + int* lineClusters, + int* lineSlots, // nCells + 1 scratch: accept flags, scanned in place into slots + const float beamX, + const float beamY, + const float maxZ, + const float minPt, + float* linesZs, + gpu::LineTime* lineTimes, + float* lineChi2, + float* linePt, + o2::its::ExternalAllocator* alloc, + gpu::Stream& stream) +{ + gpu::dedupCellsKernel<<>>( + nCells, + cells, + clusterOwners, + ownedClustersCut, + beamX, + beamY, + maxZ, + minPt, + lineSlots); + auto nosync_policy = THRUST_NAMESPACE::par_nosync(gpu::TypedAllocator(alloc)).on(stream.get()); + thrust::exclusive_scan(nosync_policy, lineSlots, lineSlots + nCells + 1, lineSlots); + gpu::linearizeCellsKernel<<>>( + nCells, + cells, + rofFramesClustersL1, + nRofsL1, + lineSlots, + lines, + lineRof, + beamX, + beamY, + linesZs, + lineTimes, + lineClusters, + lineChi2, + linePt); +} + +// Orders lines by (ROF, z): primary key the ROF, secondary the projected z within the ROF. +struct RofZLess { + const int* rof; + const float* z; + GPUhdi() bool operator()(const int a, const int b) const + { + return rof[a] != rof[b] ? rof[a] < rof[b] : z[a] < z[b]; + } +}; + +template +void sortLinesHandler(const int nLines, + const int nRofs, + const gpu::LineProjSoA soa, + const gpu::LineProjSoA sortedSoa, + const int* lineRof, + int* rofOffsets, + o2::its::ExternalAllocator* alloc, + gpu::Stream& stream) +{ + if (nLines < 2) { + return; + } + auto policy = THRUST_NAMESPACE::par_nosync(gpu::TypedAllocator(alloc)).on(stream.get()); + thrust::sequence(policy, soa.idx, soa.idx + nLines); + thrust::sort(policy, soa.idx, soa.idx + nLines, RofZLess{lineRof, soa.z}); + gpu::gatherSortedLinesKernel<<>>(nLines, soa, sortedSoa); + auto rofSorted = thrust::make_permutation_iterator(lineRof, soa.idx); + thrust::lower_bound(policy, rofSorted, rofSorted + nLines, + thrust::make_counting_iterator(0), thrust::make_counting_iterator(nRofs + 1), + rofOffsets); +} + +template +void scanDensityHandler(const int nLines, + const gpu::LineProjSoA sortedSoa, + const int* rofOffsets, + int* density, + gpu::LineWindow* win, + const float zWindow, + gpu::Stream& stream) +{ + if (nLines < 2) { + return; + } + gpu::scanDensityKernel<<>>(density, win, nLines, rofOffsets, sortedSoa, zWindow); +} + +template +void findPeaksHandler(const int nLines, + const int nRofs, + const gpu::LineProjSoA sortedSoa, + const int* rofOffsets, + const int* density, + const gpu::LineWindow* win, + uint8_t* isPeak, + const int* densityFine, + const gpu::LineWindow* winFine, + const int fineMinDensity, + uint8_t* isPeakFine, + int* peakScan, + int* peakLineIdx, + int* peakOffsets, + o2::its::ExternalAllocator* alloc, + gpu::Stream& stream) +{ + if (nLines < 2) { + return; + } + gpu::findPeaksKernel<<>>(density, win, nLines, sortedSoa, isPeak, + densityFine, winFine, fineMinDensity, isPeakFine); + auto nosync_policy = THRUST_NAMESPACE::par_nosync(gpu::TypedAllocator(alloc)).on(stream.get()); + thrust::exclusive_scan(nosync_policy, isPeak, isPeak + nLines + 1, peakScan, 0, thrust::plus()); + thrust::scatter_if(nosync_policy, thrust::make_counting_iterator(0), thrust::make_counting_iterator(nLines), + peakScan, isPeak, peakLineIdx); + thrust::gather(nosync_policy, rofOffsets, rofOffsets + nRofs + 1, peakScan, peakOffsets); +} + +template +void fitPeaksHandler(const int* nPeaksDevice, + const int* peakLineIdx, + const gpu::LineWindow* win, + const gpu::LineProjSoA sortedSoa, + const gpu::GPULine* lines, + const float* lineChi2, + const float* linePt, + const float goodLineChi2Cut, + const float goodLinePtCut, + const float pairCut2, + const float nSigmaCut, + const int minContributors, + const float beamX, + const float beamY, + const uint8_t* isPeakFine, + const float fineMaxDrift, + gpu::VertexCand* cands, + gpu::Stream& stream) +{ + gpu::fitPeaksKernel<<>>( + nPeaksDevice, peakLineIdx, win, sortedSoa, lines, lineChi2, linePt, goodLineChi2Cut, goodLinePtCut, pairCut2, nSigmaCut, minContributors, beamX, beamY, isPeakFine, fineMaxDrift, cands); +} + +template +void dedupVertexCandidatesHandler(const int* nPeaksDevice, + const int* peakLineIdx, + const int* peakOffsets, + const gpu::LineProjSoA sortedSoa, + const float duplicateZCut, + const float duplicateZScale, + gpu::VertexCand* cands, + gpu::Stream& stream) +{ + gpu::dedupVertexCandidatesKernel<<>>( + nPeaksDevice, peakLineIdx, peakOffsets, sortedSoa, duplicateZCut, duplicateZScale, cands); +} + +// MC only +template +void scanMemberOffsetsHandler(const gpu::VertexCand* cands, + int* memberOffsets, + const int nLines, + o2::its::ExternalAllocator* alloc, + gpu::Stream& stream) +{ + auto nosync_policy = THRUST_NAMESPACE::par_nosync(gpu::TypedAllocator(alloc)).on(stream.get()); + auto memberCount = thrust::make_transform_iterator(thrust::make_counting_iterator(0), gpu::MemberCount{cands}); + thrust::exclusive_scan(nosync_policy, memberCount, memberCount + nLines + 1, memberOffsets, 0, thrust::plus()); +} + +// MC only +template +void collectLinesForMCHandler(const int* nPeaksDevice, + const int* peakLineIdx, + const gpu::LineWindow* win, + const gpu::LineProjSoA sortedSoa, + const gpu::GPULine* lines, + const float pairCut2, + const gpu::VertexCand* cands, + const int* memberOffsets, + int* memberLines, + gpu::Stream& stream) +{ + gpu::collectLinesForMCKernel<<>>( + nPeaksDevice, peakLineIdx, win, sortedSoa, lines, pairCut2, cands, memberOffsets, memberLines); +} + int finalizeCellNeighboursHandler(CellNeighbour* cellNeighbours, int* neighboursLUT, const int nTargetCells, @@ -1379,4 +2098,342 @@ template struct TrackingKernels<11>; template struct TrackingKernels<13>; #endif +template void registerClusterOwnershipHandler<7>(const CellSeed* cells, + const int nCells, + unsigned long long** clusterOwnersDeviceArray, + gpu::Stream& stream); + +template void sortClustersHandler<7>(const Cluster* unsorted, Cluster* sorted, const int* clusterOffsets, + int* indexTable, const IndexTableUtils<7>* utils, + const typename ROFMaskTable<7>::View& rofMask, float beamX, float beamY, + int zBins, int phiBins, int nRofs, int nClustersLayer, int iLayer, + float* minRadiusLayer, float* maxRadiusLayer, int* keys, int* perm, + o2::its::ExternalAllocator* alloc, gpu::Stream& stream); + +template void linearizeCellsToLinesHandler<7>(const int nCells, + const CellSeed* cells, + const unsigned long long* const* clusterOwners, + const int* rofFramesClustersL1, + const int nRofsL1, + const int ownedClustersCut, + gpu::GPULine* lines, + int* lineRof, + int* lineClusters, + int* lineSlots, + const float beamX, + const float beamY, + const float maxZ, + const float minPt, + float* linesZs, + gpu::LineTime* lineTimes, + float* lineChi2, + float* linePt, + o2::its::ExternalAllocator* alloc, + gpu::Stream& stream); + +template void sortLinesHandler<7>(const int nLines, + const int nRofs, + const gpu::LineProjSoA soa, + const gpu::LineProjSoA sortedSoa, + const int* lineRof, + int* rofOffsets, + o2::its::ExternalAllocator* alloc, + gpu::Stream& stream); + +template void scanDensityHandler<7>(const int nLines, + const gpu::LineProjSoA sortedSoa, + const int* rofOffsets, + int* density, + gpu::LineWindow* win, + const float zWindow, + gpu::Stream& stream); + +template void findPeaksHandler<7>(const int nLines, + const int nRofs, + const gpu::LineProjSoA sortedSoa, + const int* rofOffsets, + const int* density, + const gpu::LineWindow* win, + uint8_t* isPeak, + const int* densityFine, + const gpu::LineWindow* winFine, + const int fineMinDensity, + uint8_t* isPeakFine, + int* peakScan, + int* peakLineIdx, + int* peakOffsets, + o2::its::ExternalAllocator* alloc, + gpu::Stream& stream); + +template void fitPeaksHandler<7>(const int* nPeaksDevice, + const int* peakLineIdx, + const gpu::LineWindow* win, + const gpu::LineProjSoA sortedSoa, + const gpu::GPULine* lines, + const float* lineChi2, + const float* linePt, + const float goodLineChi2Cut, + const float goodLinePtCut, + const float pairCut2, + const float nSigmaCut, + const int minContributors, + const float beamX, + const float beamY, + const uint8_t* isPeakFine, + const float fineMaxDrift, + gpu::VertexCand* cands, + gpu::Stream& stream); + +template void dedupVertexCandidatesHandler<7>(const int* nPeaksDevice, + const int* peakLineIdx, + const int* peakOffsets, + const gpu::LineProjSoA sortedSoa, + const float duplicateZCut, + const float duplicateZScale, + gpu::VertexCand* cands, + gpu::Stream& stream); + +template void scanMemberOffsetsHandler<7>(const gpu::VertexCand* cands, + int* memberOffsets, + const int nLines, + o2::its::ExternalAllocator* alloc, + gpu::Stream& stream); + +template void collectLinesForMCHandler<7>(const int* nPeaksDevice, + const int* peakLineIdx, + const gpu::LineWindow* win, + const gpu::LineProjSoA sortedSoa, + const gpu::GPULine* lines, + const float pairCut2, + const gpu::VertexCand* cands, + const int* memberOffsets, + int* memberLines, + gpu::Stream& stream); + +#ifdef ENABLE_UPGRADES +template void registerClusterOwnershipHandler<11>(const CellSeed* cells, + const int nCells, + unsigned long long** clusterOwnersDeviceArray, + gpu::Stream& stream); + +template void sortClustersHandler<11>(const Cluster* unsorted, Cluster* sorted, const int* clusterOffsets, + int* indexTable, const IndexTableUtils<11>* utils, + const typename ROFMaskTable<11>::View& rofMask, float beamX, float beamY, + int zBins, int phiBins, int nRofs, int nClustersLayer, int iLayer, + float* minRadiusLayer, float* maxRadiusLayer, int* keys, int* perm, + o2::its::ExternalAllocator* alloc, gpu::Stream& stream); + +template void linearizeCellsToLinesHandler<11>(const int nCells, + const CellSeed* cells, + const unsigned long long* const* clusterOwners, + const int* rofFramesClustersL1, + const int nRofsL1, + const int ownedClustersCut, + gpu::GPULine* lines, + int* lineRof, + int* lineClusters, + int* lineSlots, + const float beamX, + const float beamY, + const float maxZ, + const float minPt, + float* linesZs, + gpu::LineTime* lineTimes, + float* lineChi2, + float* linePt, + o2::its::ExternalAllocator* alloc, + gpu::Stream& stream); + +template void sortLinesHandler<11>(const int nLines, + const int nRofs, + const gpu::LineProjSoA soa, + const gpu::LineProjSoA sortedSoa, + const int* lineRof, + int* rofOffsets, + o2::its::ExternalAllocator* alloc, + gpu::Stream& stream); + +template void scanDensityHandler<11>(const int nLines, + const gpu::LineProjSoA sortedSoa, + const int* rofOffsets, + int* density, + gpu::LineWindow* win, + const float zWindow, + gpu::Stream& stream); + +template void findPeaksHandler<11>(const int nLines, + const int nRofs, + const gpu::LineProjSoA sortedSoa, + const int* rofOffsets, + const int* density, + const gpu::LineWindow* win, + uint8_t* isPeak, + const int* densityFine, + const gpu::LineWindow* winFine, + const int fineMinDensity, + uint8_t* isPeakFine, + int* peakScan, + int* peakLineIdx, + int* peakOffsets, + o2::its::ExternalAllocator* alloc, + gpu::Stream& stream); + +template void fitPeaksHandler<11>(const int* nPeaksDevice, + const int* peakLineIdx, + const gpu::LineWindow* win, + const gpu::LineProjSoA sortedSoa, + const gpu::GPULine* lines, + const float* lineChi2, + const float* linePt, + const float goodLineChi2Cut, + const float goodLinePtCut, + const float pairCut2, + const float nSigmaCut, + const int minContributors, + const float beamX, + const float beamY, + const uint8_t* isPeakFine, + const float fineMaxDrift, + gpu::VertexCand* cands, + gpu::Stream& stream); + +template void dedupVertexCandidatesHandler<11>(const int* nPeaksDevice, + const int* peakLineIdx, + const int* peakOffsets, + const gpu::LineProjSoA sortedSoa, + const float duplicateZCut, + const float duplicateZScale, + gpu::VertexCand* cands, + gpu::Stream& stream); + +template void scanMemberOffsetsHandler<11>(const gpu::VertexCand* cands, + int* memberOffsets, + const int nLines, + o2::its::ExternalAllocator* alloc, + gpu::Stream& stream); + +template void collectLinesForMCHandler<11>(const int* nPeaksDevice, + const int* peakLineIdx, + const gpu::LineWindow* win, + const gpu::LineProjSoA sortedSoa, + const gpu::GPULine* lines, + const float pairCut2, + const gpu::VertexCand* cands, + const int* memberOffsets, + int* memberLines, + gpu::Stream& stream); + +template void registerClusterOwnershipHandler<13>(const CellSeed* cells, + const int nCells, + unsigned long long** clusterOwnersDeviceArray, + gpu::Stream& stream); + +template void sortClustersHandler<13>(const Cluster* unsorted, Cluster* sorted, const int* clusterOffsets, + int* indexTable, const IndexTableUtils<13>* utils, + const typename ROFMaskTable<13>::View& rofMask, float beamX, float beamY, + int zBins, int phiBins, int nRofs, int nClustersLayer, int iLayer, + float* minRadiusLayer, float* maxRadiusLayer, int* keys, int* perm, + o2::its::ExternalAllocator* alloc, gpu::Stream& stream); + +template void linearizeCellsToLinesHandler<13>(const int nCells, + const CellSeed* cells, + const unsigned long long* const* clusterOwners, + const int* rofFramesClustersL1, + const int nRofsL1, + const int ownedClustersCut, + gpu::GPULine* lines, + int* lineRof, + int* lineClusters, + int* lineSlots, + const float beamX, + const float beamY, + const float maxZ, + const float minPt, + float* linesZs, + gpu::LineTime* lineTimes, + float* lineChi2, + float* linePt, + o2::its::ExternalAllocator* alloc, + gpu::Stream& stream); + +template void sortLinesHandler<13>(const int nLines, + const int nRofs, + const gpu::LineProjSoA soa, + const gpu::LineProjSoA sortedSoa, + const int* lineRof, + int* rofOffsets, + o2::its::ExternalAllocator* alloc, + gpu::Stream& stream); + +template void scanDensityHandler<13>(const int nLines, + const gpu::LineProjSoA sortedSoa, + const int* rofOffsets, + int* density, + gpu::LineWindow* win, + const float zWindow, + gpu::Stream& stream); + +template void findPeaksHandler<13>(const int nLines, + const int nRofs, + const gpu::LineProjSoA sortedSoa, + const int* rofOffsets, + const int* density, + const gpu::LineWindow* win, + uint8_t* isPeak, + const int* densityFine, + const gpu::LineWindow* winFine, + const int fineMinDensity, + uint8_t* isPeakFine, + int* peakScan, + int* peakLineIdx, + int* peakOffsets, + o2::its::ExternalAllocator* alloc, + gpu::Stream& stream); + +template void fitPeaksHandler<13>(const int* nPeaksDevice, + const int* peakLineIdx, + const gpu::LineWindow* win, + const gpu::LineProjSoA sortedSoa, + const gpu::GPULine* lines, + const float* lineChi2, + const float* linePt, + const float goodLineChi2Cut, + const float goodLinePtCut, + const float pairCut2, + const float nSigmaCut, + const int minContributors, + const float beamX, + const float beamY, + const uint8_t* isPeakFine, + const float fineMaxDrift, + gpu::VertexCand* cands, + gpu::Stream& stream); + +template void dedupVertexCandidatesHandler<13>(const int* nPeaksDevice, + const int* peakLineIdx, + const int* peakOffsets, + const gpu::LineProjSoA sortedSoa, + const float duplicateZCut, + const float duplicateZScale, + gpu::VertexCand* cands, + gpu::Stream& stream); + +template void scanMemberOffsetsHandler<13>(const gpu::VertexCand* cands, + int* memberOffsets, + const int nLines, + o2::its::ExternalAllocator* alloc, + gpu::Stream& stream); + +template void collectLinesForMCHandler<13>(const int* nPeaksDevice, + const int* peakLineIdx, + const gpu::LineWindow* win, + const gpu::LineProjSoA sortedSoa, + const gpu::GPULine* lines, + const float pairCut2, + const gpu::VertexCand* cands, + const int* memberOffsets, + int* memberLines, + gpu::Stream& stream); +#endif + } // namespace o2::its diff --git a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/ClusterLines.h b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/ClusterLines.h index bcb8a98a62cab..2e0960964b527 100644 --- a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/ClusterLines.h +++ b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/ClusterLines.h @@ -32,6 +32,7 @@ struct Line final { Line() = default; Line(const Tracklet&, const Cluster*, const Cluster*); + Line(const std::array& origin, const std::array& direction, const TimeEstBC& time); bool operator==(const Line&) const = default; static float getDistance2FromPoint(const Line& line, const std::array& point); diff --git a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/Configuration.h b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/Configuration.h index 8e1c68fc31c6c..ecfe64794fbae 100644 --- a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/Configuration.h +++ b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/Configuration.h @@ -16,6 +16,8 @@ #ifndef TRACKINGITSU_INCLUDE_CONFIGURATION_H_ #define TRACKINGITSU_INCLUDE_CONFIGURATION_H_ +#include +#include #include #ifndef GPUCA_GPUCODE_DEVICE #include @@ -25,6 +27,7 @@ #include "CommonUtils/EnumFlags.h" #include "DetectorsBase/Propagator.h" +#include "CommonConstants/MathConstants.h" #include "ITStracking/Constants.h" #include "ITStracking/LayerMask.h" @@ -42,6 +45,7 @@ enum class IterationStep : uint16_t { MarkVerticesAsUPC, TrackFollowerTop, TrackFollowerBot, + SeedingVertexPass, // this iteration runs the seeding-vertex step instead of track finding }; using IterationSteps = o2::utils::EnumFlags; @@ -111,6 +115,8 @@ struct TrackingParameters { float TrackletMinPt = 0.3f; /// Cell finding cuts float CellDeltaTanLambdaSigma = 0.007f; + float CellDeltaTanLambdaNSigma = -1.f; + float CellDeltaPhiMinPt = -1.f; /// Fitter parameters o2::base::PropagatorImpl::MatCorrType CorrType = o2::base::PropagatorImpl::MatCorrType::USEMatCorrNONE; float MaxChi2ClusterAttachment = 60.f; @@ -138,8 +144,24 @@ struct TrackingParameters { float SharedClusterMaxDeltaEta = 0.03f; // For tracks sharing clusters, maximum allowed delta eta at the cluster position bool SharedClusterOppositeSign = false; // For tracks sharing clusters, require opposite sign of the tracklets int SharedMaxClusters = 0; // Maximal allowed shared clusters (excluding first cluster) + + int CellLineSharedClusterCut = 2; // Seeding-vertex pass: max clusters a cell->Line may share with already-accepted Lines before it is dropped }; +inline float cellDeltaPhiBound(const float bz, const float ptMin, + const float rIn, const float rMid, const float rOut, + const float msAngle) +{ + if (ptMin <= 0.f) { + return -1.f; + } + const float oneOverR = 0.001f * 0.3f * std::abs(bz) / ptMin; // 1 / curvature radius [1/cm] + const float sA = std::min(0.25f * (rIn + rMid) * oneOverR, 1.f - 1.e-6f); + const float sB = std::min(0.25f * (rMid + rOut) * oneOverR, 1.f - 1.e-6f); + const float bound = 2.f * (std::asin(sB) - std::asin(sA)) + 2.f * msAngle; + return std::min(bound, static_cast(o2::constants::math::PI)); +} + struct VertexingParameters { std::string asString() const; @@ -171,7 +193,8 @@ struct VertexingParameters { int zSpan = -1; bool SaveTimeBenchmarks = false; - bool useTruthSeeding = false; // overwrite found vertices with MC events + bool useTruthSeeding = false; // overwrite found vertices with MC events + bool useParallelSeeding = false; // use the GPU-oriented parallel seeding int nThreads = 1; bool PrintMemory = false; // print allocator usage in epilog report diff --git a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/LineProjection.h b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/LineProjection.h new file mode 100644 index 0000000000000..951a6b413eb57 --- /dev/null +++ b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/LineProjection.h @@ -0,0 +1,36 @@ +// Copyright 2019-2026 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 LineProjection.h +/// \brief small types shared by the host and device seeding vertexers, describing lines +/// projected onto the beam line: their time interval and their z-window bounds. + +#ifndef O2_ITS_TRACKING_LINE_PROJECTION_H_ +#define O2_ITS_TRACKING_LINE_PROJECTION_H_ + +namespace o2::its +{ + +// Symmetrised time interval of a line: centre +/- half-width. +struct LineTime { + float tc{0.f}; // time centre + float th{0.f}; // time half-width +}; + +// Half-open [lo, hi) range of sorted-line slots falling inside one z-window. +struct LineWindow { + int lo{0}; + int hi{0}; +}; + +} // namespace o2::its + +#endif diff --git a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/TimeFrame.h b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/TimeFrame.h index 11246fa0ee3b0..89a2a0b5a5c33 100644 --- a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/TimeFrame.h +++ b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/TimeFrame.h @@ -239,6 +239,12 @@ struct TimeFrame { void computeTracletsPerClusterScans(); int& getNTrackletsROF(int rofId, int combId) { return mNTrackletsPerROF[combId][rofId]; } auto& getLines(int rofId) { return mLines[rofId]; } + struct LineQuality { + float chi2{-1.f}; + float pt{-1.f}; + }; + auto& getLinesQuality(int rofId) { return mLinesQuality[rofId]; } + const auto& getLinesQuality(int rofId) const { return mLinesQuality[rofId]; } int getNLinesTotal() const noexcept { return mTotalLines; } void setNLinesTotal(uint32_t a) noexcept { mTotalLines = a; } auto& getTrackletClusters(int rofId) { return mTrackletClusters[rofId]; } @@ -306,12 +312,14 @@ struct TimeFrame { virtual const char* getName() const noexcept { return "CPU"; } protected: - void prepareClusters(const TrackingParameters& trkParam, const int maxLayers = NLayers); + virtual void prepareClusters(const TrackingParameters& trkParam, const int maxLayers = NLayers); + virtual void allocateClusterSortStorage(const TrackingParameters& trkParam, const int maxLayers); float mBz = 5.; unsigned int mNTotalLowPtVertices = 0; int mBeamPosWeight = 0; std::array mBeamPos = {0.f, 0.f}; bool isBeamPositionOverridden = false; + bool mSystErrorsApplied = false; std::array mMinR; std::array mMaxR; bounded_vector mLinkPhiCuts; @@ -332,6 +340,7 @@ struct TimeFrame { bounded_vector mPrimaryVerticesLabels; std::vector> mNTrackletsPerROF; std::vector> mLines; + std::vector> mLinesQuality; // lockstep with mLines, see getLinesQuality() std::vector> mTrackletClusters; std::array, 2> mTrackletsIndexROF; std::vector> mLinesLabels; diff --git a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/Tracker.h b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/Tracker.h index c60f280de5307..b8fe6d0149f39 100644 --- a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/Tracker.h +++ b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/Tracker.h @@ -59,6 +59,10 @@ class Tracker const LogFunc& = [](const std::string& s) { std::cout << s << '\n'; }, const LogFunc& = [](const std::string& s) { std::cerr << s << '\n'; }); + // Self-contained seeding-vertex phase (own TimeFrame init + own timer), run before clustersToTracks. + float clustersToVertices( + const LogFunc& = [](const std::string& s) { std::cout << s << '\n'; }); + void setParameters(const std::vector& p) { mTrkParams = p; } void setMemoryPool(std::shared_ptr pool) { mMemoryPool = pool; } std::vector& getParameters() { return mTrkParams; } @@ -72,6 +76,9 @@ class Tracker void initialiseTimeFrame(int iteration) { mTraits->initialiseTimeFrame(iteration); } void computeTracklets(int iteration, int iVertex) { mTraits->computeLayerTracklets(iteration, iVertex); } void computeCells(int iteration) { mTraits->computeLayerCells(iteration); } + void computeVertexCandidates(int iteration) { mTraits->computeVertexCandidates(iteration); } + void computeBeamFromVertices(int iteration) { mTraits->computeBeamFromVertices(iteration); } + void computeVertices(int iteration) { mTraits->computeVertices(iteration); } void findCellsNeighbours(int iteration) { mTraits->findCellsNeighbours(iteration); } void findRoads(int iteration) { mTraits->findRoads(iteration); } @@ -100,10 +107,13 @@ class Tracker Neighbouring, Roading, Extending, + CellLinearising, + BeamPositioning, + SeedingVertices, NSteps, }; Steps mCurStep{TFInit}; - static constexpr std::array StateNames{"TimeFrame initialisation", "Tracklet finding", "Cell finding", "Neighbour finding", "Road finding", "Track extending"}; + static constexpr std::array StateNames{"TimeFrame initialisation", "Tracklet finding", "Cell finding", "Neighbour finding", "Road finding", "Track extending", "Cell linearisation", "Beam positioning", "Seeding vertices"}; 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 1e3b9e7d7fbf7..cea5e71d6fa4f 100644 --- a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/TrackerTraits.h +++ b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/TrackerTraits.h @@ -71,6 +71,10 @@ class TrackerTraits virtual void findCellsNeighbours(const int iteration); virtual void findRoads(const int iteration); + virtual void computeVertexCandidates(const int iteration); + virtual void computeBeamFromVertices(const int iteration); + virtual void computeVertices(const int iteration); + template void processNeighbours(int iteration, int defaultCellTopologyId, int iLevel, uint64_t capacityKey, const bounded_vector& currentSeeds, bounded_vector& updatedSeeds); @@ -99,6 +103,8 @@ class TrackerTraits virtual int getTFNumberOfCells() const { return mTimeFrame->getNumberOfCells(); } private: + bool skipROF(int iteration, int rof) const; // seeding-vertex pass: skip ROFs above the per-ROF vertex threshold + std::shared_ptr mMemoryPool; protected: diff --git a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/TrackingConfigParam.h b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/TrackingConfigParam.h index ee67b7234450d..3aa5e822e9e0d 100644 --- a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/TrackingConfigParam.h +++ b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/TrackingConfigParam.h @@ -44,6 +44,14 @@ struct VertexerParamConfig : public o2::conf::ConfigurableParamHelpercell->line->parallel seeding) to the tracker passes + float diamondTrackletingNSigmaCut = 6.6136f; // z-window n-sigma for the diamond trackleting (gates the cell tanLambda too unless the next one is set) + float diamondCellTanLambdaNSigma = -1.f; + float diamondTrackletingPVres = 0.9676f; // effective PV resolution [cm] feeding the trackleting z-sigma + float diamondTrackletingCellDeltaTanLambdaSigma = 0.007f; // cell tanLambda sigma for the diamond celling + float diamondTrackletingCellDeltaPhiMinPt = -1.f; + int cellLineSharedClusterCut = 1; // Max clusters a cell->Line may share with already-accepted Lines before being dropped. + + bool saveTimeBenchmarks = false; // dump metrics on file + bool overrideBeamEstimation = false; // use beam position from meanVertex CCDB object + int trackingMode = -1; // -1: unset, 0=sync, 1=async, 2=cosmics used by gpuwf only + bool doUPCIteration = false; // Perform an additional iteration for UPC events on tagged vertices. You want to combine this config with VertexerParamConfig.nIterations=2 + int nIterations = constants::MaxIter; // overwrite the number of iterations + int reseedIfShorter = 6; // for the final refit reseed the track with circle if they are shorter than this value + bool shiftRefToCluster{true}; // TrackFit: after update shift the linearization reference to cluster + bool repeatRefitOut{false}; // repeat outward refit using inward refit as a seed + bool createArtefactLabels{false}; // create on-the-fly labels for the artefacts bool trackFollowerTop[constants::MaxIter] = {}; bool trackFollowerBot[constants::MaxIter] = {}; float trackFollowerNSigmaCutZ = 1.f; diff --git a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/TrackingInterface.h b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/TrackingInterface.h index 14c5d6a62e0ad..5a3381769e3b3 100644 --- a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/TrackingInterface.h +++ b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/TrackingInterface.h @@ -29,6 +29,8 @@ #include "GPUChainITS.h" #include +#include +#include namespace o2::its { @@ -103,6 +105,7 @@ class ITSTrackingInterface const o2::dataformats::MeanVertexObject* mMeanVertex{}; std::shared_ptr mMemoryPool; std::shared_ptr mTaskArena; + // MC performance / vertex-dump state }; } // namespace o2::its diff --git a/Detectors/ITSMFT/ITS/tracking/src/ClusterLines.cxx b/Detectors/ITSMFT/ITS/tracking/src/ClusterLines.cxx index 3e3e1b8b46338..6af528a10fbac 100644 --- a/Detectors/ITSMFT/ITS/tracking/src/ClusterLines.cxx +++ b/Detectors/ITSMFT/ITS/tracking/src/ClusterLines.cxx @@ -28,6 +28,13 @@ Line::Line(const Tracklet& tracklet, const Cluster* innerClusters, const Cluster cosinesDirector /= std::sqrt(ROOT::Math::Dot(cosinesDirector, cosinesDirector)); } +Line::Line(const std::array& origin, const std::array& direction, const TimeEstBC& time) : mTime(time) +{ + originPoint = SVector3f(origin[0], origin[1], origin[2]); + cosinesDirector = SVector3f(direction[0], direction[1], direction[2]); + cosinesDirector /= std::sqrt(ROOT::Math::Dot(cosinesDirector, cosinesDirector)); +} + float Line::getDistance2FromPoint(const Line& line, const std::array& point) { const SVector3f p(point.data(), 3); diff --git a/Detectors/ITSMFT/ITS/tracking/src/Configuration.cxx b/Detectors/ITSMFT/ITS/tracking/src/Configuration.cxx index eb4d90ee9d15f..fe93e9a5c9db1 100644 --- a/Detectors/ITSMFT/ITS/tracking/src/Configuration.cxx +++ b/Detectors/ITSMFT/ITS/tracking/src/Configuration.cxx @@ -315,6 +315,7 @@ std::vector TrackingMode::getVertexingParameters(TrackingMo p.ZBins = vc.ZBins; p.PhiBins = vc.PhiBins; p.useTruthSeeding = vc.useTruthSeeding; + p.useParallelSeeding = vc.useParallelSeeding; p.maxTrackletsPerCluster = vc.maxTrackletsPerCluster; p.zCut = vc.zCut; p.phiCut = vc.phiCut; diff --git a/Detectors/ITSMFT/ITS/tracking/src/TimeFrame.cxx b/Detectors/ITSMFT/ITS/tracking/src/TimeFrame.cxx index 08c0164288388..d92a02dab9f20 100644 --- a/Detectors/ITSMFT/ITS/tracking/src/TimeFrame.cxx +++ b/Detectors/ITSMFT/ITS/tracking/src/TimeFrame.cxx @@ -187,6 +187,17 @@ void TimeFrame::prepareROFrameData(gsl::span +void TimeFrame::allocateClusterSortStorage(const TrackingParameters& trkParam, const int maxLayers) +{ + for (unsigned int iLayer{0}; iLayer < std::min((int)mClusters.size(), maxLayers); ++iLayer) { + clearResizeBoundedVector(mClusters[iLayer], mUnsortedClusters[iLayer].size(), getMaybeFrameworkHostResource(maxLayers != NLayers)); + } + for (int iLayer{0}; iLayer < NLayers; ++iLayer) { + clearResizeBoundedVector(mIndexTables[iLayer], getNrof(iLayer) * ((trkParam.ZBins * trkParam.PhiBins) + 1), getMaybeFrameworkHostResource()); + } +} + template void TimeFrame::prepareClusters(const TrackingParameters& trkParam, const int maxLayers) { @@ -301,6 +312,7 @@ void TimeFrame::initialise(const TrackingParameters& trkParam, const in deepVectorClear(mTracks); deepVectorClear(mTracksLabel); deepVectorClear(mLines); + deepVectorClear(mLinesQuality); deepVectorClear(mLinesLabels); if (trkParam.PassFlags[IterationStep::ResetVertices]) { deepVectorClear(mPrimaryVertices); @@ -308,28 +320,27 @@ void TimeFrame::initialise(const TrackingParameters& trkParam, const in } clearResizeBoundedVector(mLinesLabels, getNrof(1), mMemoryPool.get()); mIndexTableUtils.setTrackingParameters(trkParam); - clearResizeBoundedVector(mPositionResolution, trkParam.NLayers, mMemoryPool.get()); - clearResizeBoundedVector(mBogusClusters, trkParam.NLayers, mMemoryPool.get()); + clearResizeBoundedVector(mPositionResolution, NLayers, mMemoryPool.get()); + clearResizeBoundedVector(mBogusClusters, NLayers, mMemoryPool.get()); deepVectorClear(mTrackletClusters); for (unsigned int iLayer{0}; iLayer < std::min((int)mClusters.size(), maxLayers); ++iLayer) { - clearResizeBoundedVector(mClusters[iLayer], mUnsortedClusters[iLayer].size(), getMaybeFrameworkHostResource(maxLayers != NLayers)); clearResizeBoundedVector(mUsedClusters[iLayer], mUnsortedClusters[iLayer].size(), getMaybeFrameworkHostResource(maxLayers != NLayers)); mPositionResolution[iLayer] = o2::gpu::CAMath::Sqrt((0.5f * (trkParam.SystErrorZ2[iLayer] + trkParam.SystErrorY2[iLayer])) + (trkParam.LayerResolution[iLayer] * trkParam.LayerResolution[iLayer])); } clearResizeBoundedVector(mLines, getNrof(1), mMemoryPool.get()); + clearResizeBoundedVector(mLinesQuality, getNrof(1), mMemoryPool.get()); clearResizeBoundedVector(mTrackletClusters, getNrof(1), mMemoryPool.get()); - - for (int iLayer{0}; iLayer < NLayers; ++iLayer) { - clearResizeBoundedVector(mIndexTables[iLayer], getNrof(iLayer) * ((trkParam.ZBins * trkParam.PhiBins) + 1), getMaybeFrameworkHostResource()); - } - for (int iLayer{0}; iLayer < trkParam.NLayers; ++iLayer) { - if (trkParam.SystErrorY2[iLayer] > 0.f || trkParam.SystErrorZ2[iLayer] > 0.f) { - for (auto& tfInfo : mTrackingFrameInfo[iLayer]) { - /// Account for alignment systematics in the cluster covariance matrix - tfInfo.covarianceTrackingFrame[0] += trkParam.SystErrorY2[iLayer]; - tfInfo.covarianceTrackingFrame[2] += trkParam.SystErrorZ2[iLayer]; + allocateClusterSortStorage(trkParam, maxLayers); + if (!mSystErrorsApplied) { + for (int iLayer{0}; iLayer < trkParam.NLayers; ++iLayer) { + if (trkParam.SystErrorY2[iLayer] > 0.f || trkParam.SystErrorZ2[iLayer] > 0.f) { + for (auto& tfInfo : mTrackingFrameInfo[iLayer]) { + tfInfo.covarianceTrackingFrame[0] += trkParam.SystErrorY2[iLayer]; + tfInfo.covarianceTrackingFrame[2] += trkParam.SystErrorZ2[iLayer]; + } } } + mSystErrorsApplied = true; } mMinR.fill(std::numeric_limits::max()); @@ -505,6 +516,7 @@ template void TimeFrame::wipe() { resetTrackExtensionCounters(); + mSystErrorsApplied = false; deepVectorClear(mTracks); deepVectorClear(mTracklets); deepVectorClear(mCells); @@ -526,6 +538,7 @@ void TimeFrame::wipe() deepVectorClear(mTrackletsIndexROF); deepVectorClear(mTrackletClusters); deepVectorClear(mLines); + deepVectorClear(mLinesQuality); // if we use the external host allocator then the assumption is that we // don't clear the memory ourself if (!hasFrameworkAllocator()) { diff --git a/Detectors/ITSMFT/ITS/tracking/src/Tracker.cxx b/Detectors/ITSMFT/ITS/tracking/src/Tracker.cxx index d8ff8442f908f..a2145bb240dc9 100644 --- a/Detectors/ITSMFT/ITS/tracking/src/Tracker.cxx +++ b/Detectors/ITSMFT/ITS/tracking/src/Tracker.cxx @@ -35,6 +35,48 @@ Tracker::Tracker(TrackerTraits* traits) : mTraits(traits) { } +template +float Tracker::clustersToVertices(const LogFunc& logger) +{ + mTraits->updateTrackingParameters(mTrkParams); + int it{-1}; + for (int i = 0; i < (int)mTrkParams.size(); ++i) { + if (mTrkParams[i].PassFlags[IterationStep::SeedingVertexPass]) { + it = i; + break; + } + } + if (it < 0) { + return 0.f; + } + mMemoryPool->setMaxMemory(mTrkParams[it].MaxMemory); + logger(std::format("==== ITS {} Seeding-vertex pass ====", mTraits->getName())); + float total{0.f}; + constexpr int kMaxBootstrapPasses = 5; + constexpr float kBeamConvergence2 = 5.e-3f * 5.e-3f; // (50 um)^2 + for (int pass = 0; pass < kMaxBootstrapPasses; ++pass) { + const float prevBeamX = mTimeFrame->getBeamX(); + const float prevBeamY = mTimeFrame->getBeamY(); + total += evaluateTask(&Tracker::initialiseTimeFrame, StateNames[mCurStep = TFInit], it, logger, it); + total += evaluateTask(&Tracker::computeTracklets, StateNames[mCurStep = Trackleting], it, logger, it, -1); + const int nTracklets = mTraits->getTFNumberOfTracklets(); + total += evaluateTask(&Tracker::computeCells, StateNames[mCurStep = Celling], it, logger, it); + total += evaluateTask(&Tracker::computeVertexCandidates, StateNames[mCurStep = CellLinearising], it, logger, it); + logger(std::format(" - Seeding pass {}: {} tracklets, {} cells, {} lines", pass, + nTracklets, mTraits->getTFNumberOfCells(), + mTimeFrame->getNLinesTotal())); + total += evaluateTask(&Tracker::computeVertices, StateNames[mCurStep = SeedingVertices], it, logger, it); + total += evaluateTask(&Tracker::computeBeamFromVertices, StateNames[mCurStep = BeamPositioning], it, logger, it); + const float dx = mTimeFrame->getBeamX() - prevBeamX; + const float dy = mTimeFrame->getBeamY() - prevBeamY; + if (dx * dx + dy * dy < kBeamConvergence2) { + logger(std::format(" - Beam bootstrap converged after pass {} (beam shift < 50 um)", pass)); + break; + } + } + return total; +} + template float Tracker::clustersToTracks(const LogFunc& logger, const LogFunc& error) { @@ -44,7 +86,11 @@ float Tracker::clustersToTracks(const LogFunc& logger, const LogFunc& e mTraits->updateTrackingParameters(mTrkParams); int maxNvertices{-1}; - if (mTrkParams[0].PerPrimaryVertexProcessing) { + int firstTrackingIteration{0}; + while (firstTrackingIteration < (int)mTrkParams.size() && mTrkParams[firstTrackingIteration].PassFlags[IterationStep::SeedingVertexPass]) { + ++firstTrackingIteration; + } + if (firstTrackingIteration < (int)mTrkParams.size() && mTrkParams[firstTrackingIteration].PerPrimaryVertexProcessing) { maxNvertices = mTimeFrame->getROFVertexLookupTableView().getMaxVerticesPerROF(); } @@ -80,6 +126,9 @@ float Tracker::clustersToTracks(const LogFunc& logger, const LogFunc& e if (mTrkParams[iteration].PassFlags[IterationStep::UseUPCMask]) { mTimeFrame->useUPCMask(); } + if (mTrkParams[iteration].PassFlags[IterationStep::SeedingVertexPass]) { + continue; // the seeding-vertex pass runs as its own phase (clustersToVertices), not as a tracking iteration + } 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())}; diff --git a/Detectors/ITSMFT/ITS/tracking/src/TrackerTraits.cxx b/Detectors/ITSMFT/ITS/tracking/src/TrackerTraits.cxx index 79511e6a9add5..3f3ba6dd21303 100644 --- a/Detectors/ITSMFT/ITS/tracking/src/TrackerTraits.cxx +++ b/Detectors/ITSMFT/ITS/tracking/src/TrackerTraits.cxx @@ -15,7 +15,10 @@ #include #include +#include +#include #include +#include #include #include #include @@ -35,12 +38,16 @@ #include "ITStracking/Constants.h" #include "ITStracking/IndexTableUtils.h" #include "ITStracking/LayerMask.h" +#include "ITStracking/LineProjection.h" #include "ITStracking/ROFLookupTables.h" #include "ITStracking/SlabBumpAllocator.h" #include "ITStracking/TrackerTraits.h" #include "ITStracking/TrackFollower.h" #include "ITStracking/TrackHelpers.h" #include "ITStracking/Tracklet.h" +#include "ITStracking/TrackingConfigParam.h" + +#include namespace o2::its { @@ -57,6 +64,8 @@ template void TrackerTraits::computeLayerTracklets(const int iteration, int iVertex) { const auto topology = mTimeFrame->getTrackingTopologyView(); + const bool vtxMode = mTrkParams[iteration].PassFlags[IterationStep::SeedingVertexPass]; + const Vertex diamondVert(mTrkParams[iteration].Diamond, mTrkParams[iteration].DiamondCov, 1, 1.f); gsl::span diamondSpan(&diamondVert, 1); @@ -94,7 +103,7 @@ void TrackerTraits::computeLayerTracklets(const int iteration, int iVer } const float meanDeltaR = mTrkParams[iteration].LayerRadii[link.toLayer] - mTrkParams[iteration].LayerRadii[link.fromLayer]; - const float phiCut = mTimeFrame->getLinkPhiCut(linkId); + const float phiCut = vtxMode ? o2::its::VertexerParamConfig::Instance().phiCut : mTimeFrame->getLinkPhiCut(linkId); const float msAngle = mTimeFrame->getLinkMSAngle(linkId); for (int iCluster = 0; iCluster < int(layer0.size()); ++iCluster) { @@ -107,7 +116,7 @@ void TrackerTraits::computeLayerTracklets(const int iteration, int iVer for (int iV = startVtx; iV < endVtx; ++iV) { const auto& pv = primaryVertices[iV]; - if (!mTimeFrame->getROFVertexLookupTableView().isVertexCompatible(link.fromLayer, pivotROF, pv)) { + if (!vtxMode && !mTimeFrame->getROFVertexLookupTableView().isVertexCompatible(link.fromLayer, pivotROF, pv)) { continue; } if (pv.isFlagSet(Vertex::Flags::UPCMode) != mTrkParams[iteration].PassFlags[IterationStep::SelectUPCVertices]) { @@ -139,7 +148,7 @@ void TrackerTraits::computeLayerTracklets(const int iteration, int iVer continue; } const auto ts = mTimeFrame->getROFOverlapTableView().getTimeStamp(link.fromLayer, pivotROF, link.toLayer, targetROF); - if (!ts.isCompatible(pv.getTimeStamp())) { + if (!vtxMode && !ts.isCompatible(pv.getTimeStamp())) { continue; } const auto& targetIndexTable = mTimeFrame->getIndexTable(targetROF, link.toLayer); @@ -275,6 +284,561 @@ void TrackerTraits::computeLayerTracklets(const int iteration, int iVer }); } +template +void TrackerTraits::computeVertexCandidates(const int iteration) +{ + const auto& cells = mTimeFrame->getCells()[0]; + const bool withMC = mTimeFrame->hasMCinformation() && mTrkParams[iteration].CreateArtefactLabels; + const auto& cellLabels = mTimeFrame->getCellsLabel(0); + const auto nClusterLayer0 = mTimeFrame->getUnsortedClusters()[0].size(); + const auto nClusterLayer1 = mTimeFrame->getUnsortedClusters()[1].size(); + const auto nClusterLayer2 = mTimeFrame->getUnsortedClusters()[2].size(); + const int nCells = static_cast(cells.size()); + const int keepThreshold = 3 - mTrkParams[iteration].CellLineSharedClusterCut; + const auto& vtxConf = o2::its::VertexerParamConfig::Instance(); + const float maxZ = vtxConf.maxZPositionAllowed; + const float lineMinPt = vtxConf.lineMinPt; + const float beamX = mTimeFrame->getBeamX(); + const float beamY = mTimeFrame->getBeamY(); + auto makeKey = [](float attribute, int cellIdx) -> size_t { + const uint32_t attributeInt = std::bit_cast(attribute); + return (static_cast(attributeInt) << 32) | static_cast(cellIdx); + }; + + // Per-cell precompute (geometry, line, pivot ROF, label), indexed directly by cell id k. + bounded_vector kCl0(nCells, 0, mMemoryPool.get()); + bounded_vector kCl1(nCells, 0, mMemoryPool.get()); + bounded_vector kCl2(nCells, 0, mMemoryPool.get()); + bounded_vector kPivotRof(nCells, 0, mMemoryPool.get()); + bounded_vector kGeomOk(nCells, 0, mMemoryPool.get()); + bounded_vector kProjOk(nCells, 0, mMemoryPool.get()); + bounded_vector kLine(nCells, mMemoryPool.get()); + bounded_vector kLabel(withMC ? nCells : 0, o2::MCCompLabel{}, mMemoryPool.get()); + mTaskArena->execute([&] { + tbb::parallel_for(0, nCells, [&](const int k) { + const auto& cell = cells[k]; + const int c1 = cell.getSecondClusterIndex(); + kCl0[k] = cell.getFirstClusterIndex(); + kCl1[k] = c1; + kCl2[k] = cell.getThirdClusterIndex(); + std::array origin, direction; + cell.getXYZGlo(origin); + if (!cell.getPxPyPzGlo(direction)) { + return; + } + kGeomOk[k] = 1; + kLine[k] = Line(origin, direction, cell.getTimeStamp()); + kPivotRof[k] = mTimeFrame->getClusterROF(1, c1); + const float dx = origin[0] - beamX; + const float dy = origin[1] - beamY; + const float den = direction[0] * direction[0] + direction[1] * direction[1]; + kProjOk[k] = den >= constants::Tolerance && + o2::gpu::CAMath::Abs(origin[2] - (dx * direction[0] + dy * direction[1]) / den * direction[2]) < maxZ; + if (withMC) { + kLabel[k] = cellLabels[k]; + } + }); + }); + + bounded_vector clusterOwnershipLayer0(nClusterLayer0, std::numeric_limits::max(), mMemoryPool.get()); + bounded_vector clusterOwnershipLayer1(nClusterLayer1, std::numeric_limits::max(), mMemoryPool.get()); + bounded_vector clusterOwnershipLayer2(nClusterLayer2, std::numeric_limits::max(), mMemoryPool.get()); + mTaskArena->execute([&] { + tbb::parallel_for(0, nCells, [&](const int k) { + if (!kGeomOk[k]) { + return; + } + const float pt = cells[k].getPt(); + const size_t key = makeKey(pt > 1.e-6f ? 1.f / pt : 1.e9f, k); + o2::gpu::GPUCommonMath::AtomicMin(&clusterOwnershipLayer0[kCl0[k]], key); + o2::gpu::GPUCommonMath::AtomicMin(&clusterOwnershipLayer1[kCl1[k]], key); + o2::gpu::GPUCommonMath::AtomicMin(&clusterOwnershipLayer2[kCl2[k]], key); + }); + }); + + bounded_vector cellAccepted(nCells, 0, mMemoryPool.get()); + mTaskArena->execute([&] { + tbb::parallel_for(0, nCells, [&](const int k) { + if (!kGeomOk[k]) { + return; + } + int nOwned = 0; + nOwned += static_cast(clusterOwnershipLayer0[kCl0[k]]) == static_cast(k); + nOwned += static_cast(clusterOwnershipLayer1[kCl1[k]]) == static_cast(k); + nOwned += static_cast(clusterOwnershipLayer2[kCl2[k]]) == static_cast(k); + const bool ptOk = lineMinPt <= 0.f || cells[k].getPt() >= lineMinPt; + cellAccepted[k] = (nOwned >= keepThreshold) && kProjOk[k] && ptOk ? 1 : 0; + }); + }); + + int total{0}; + for (int k = 0; k < nCells; ++k) { + if (!cellAccepted[k]) { + continue; + } + const int pivotRof = kPivotRof[k]; + mTimeFrame->getLines(pivotRof).emplace_back(kLine[k]); + mTimeFrame->getLinesQuality(pivotRof).emplace_back( + typename TimeFrame::LineQuality{cells[k].getChi2(), cells[k].getPt()}); + if (withMC) { + mTimeFrame->getLinesLabel(pivotRof).emplace_back(kLabel[k]); + } + ++total; + } + mTimeFrame->setNLinesTotal(total); +} + +namespace +{ +// Majority MC label of a set of line labels (ported from VertexerTraits): we only care about the +// source&event of the contributing tracks, not their trackId; flag as fake if no label has a strict +// majority (>50%). +VertexLabel computeMain(const bounded_vector& elements) +{ + auto composeVtxLabel = [](const o2::MCCompLabel& lbl) -> o2::MCCompLabel { + return {o2::MCCompLabel::maxTrackID(), lbl.getEventID(), lbl.getSourceID(), lbl.isFake()}; + }; + std::unordered_map frequency; + for (const auto& element : elements) { + ++frequency[composeVtxLabel(element)]; + } + o2::MCCompLabel elem{}; + size_t maxCount = 0; + for (const auto& [key, count] : frequency) { + if (count > maxCount) { + maxCount = count; + elem = key; + } + } + if (maxCount <= 1) { // need >50% + elem.setFakeFlag(); + } + return std::make_pair(elem, static_cast(maxCount) / static_cast(elements.size())); +} + +inline bool timeCompatible(const LineTime& a, const LineTime& b) +{ + return o2::gpu::GPUCommonMath::Abs(a.tc - b.tc) <= (a.th + b.th); +} + +void computeDensities(const bounded_vector& Z, + const bounded_vector& T, + const float zWindow, + const bool incremental, + bounded_vector& win, + bounded_vector& density) +{ + const int m = static_cast(Z.size()); + int lo{0}, hi{0}; + for (int k = 0; k < m; ++k) { + const float zk = Z[k]; + if (incremental) { + while (lo < m && Z[lo] < zk - zWindow) { + ++lo; + } + while (hi < m && Z[hi] <= zk + zWindow) { + ++hi; + } + } else { + lo = static_cast(std::lower_bound(Z.begin(), Z.end(), zk - zWindow) - Z.begin()); + hi = static_cast(std::upper_bound(Z.begin(), Z.end(), zk + zWindow) - Z.begin()); + } + win[k] = LineWindow{lo, hi}; + int count = 0; + for (int j = lo; j < hi; ++j) { + count += timeCompatible(T[k], T[j]); + } + density[k] = count; + } +} + +void markLeftmostMaxima(const bounded_vector& density, + const bounded_vector& Z, + const bounded_vector& win, + const bool incremental, + std::pmr::memory_resource* pool, + bounded_vector& isMax) +{ + const int m = static_cast(density.size()); + if (incremental) { + bounded_vector maxDeque(m, pool); + int front{0}, back{0}, next{0}; + for (int k = 0; k < m; ++k) { + while (next < win[k].hi) { + while (back > front && density[maxDeque[back - 1]] < density[next]) { + --back; + } + maxDeque[back++] = next++; + } + while (front < back && maxDeque[front] < win[k].lo) { + ++front; + } + isMax[k] = maxDeque[front] == k; // the window always contains k itself, so the deque is non-empty + } + return; + } + for (int k = 0; k < m; ++k) { + uint8_t best = 1; + for (int j = win[k].lo; j < win[k].hi; ++j) { + if (j == k) { + continue; + } + if (density[j] > density[k] || (density[j] == density[k] && Z[j] < Z[k])) { + best = 0; + break; + } + } + isMax[k] = best; + } +} +} // namespace + +template +bool TrackerTraits::skipROF(int iteration, int rof) const +{ + return mTrkParams[iteration].PassFlags[IterationStep::SkipROFsAboveThreshold] && + (int)mTimeFrame->getROFVertexLookupTableView().getVertices(1, rof).getEntries() > + o2::its::VertexerParamConfig::Instance().vertPerRofThreshold; +} + +template +void TrackerTraits::computeBeamFromVertices(const int) +{ + const auto& vertices = mTimeFrame->getPrimaryVertices(); + double sx{0.}, sy{0.}, sw{0.}; + for (const auto& v : vertices) { + const double w = static_cast(v.getNContributors()); + sx += w * v.getX(); + sy += w * v.getY(); + sw += w; + } + if (sw <= 0.) { + return; + } + const float bx = static_cast(sx / sw); + const float by = static_cast(sy / sw); + mTimeFrame->resetBeamXY(bx, by, static_cast(sw)); + LOGP(info, "[beamFit] beam from {} vertices = ({:.4f}, {:.4f}) cm (sum nContrib={:.0f})", vertices.size(), bx, by, sw); +} + +template +void TrackerTraits::computeVertices(const int iteration) +{ + const int nRofs = mTimeFrame->getNrof(1); + const bool withMC = mTimeFrame->hasMCinformation() && mTrkParams[iteration].CreateArtefactLabels; + std::vector> rofVertices(nRofs); + std::vector> rofLabels(nRofs); + + const float beamX = mTimeFrame->getBeamX(); + const float beamY = mTimeFrame->getBeamY(); + const auto& vc = o2::its::VertexerParamConfig::Instance(); + const float maxZ = vc.maxZPositionAllowed; + const float clusterCut = vc.clusterCut; + const float pairCut2 = vc.pairCut * vc.pairCut; + const float zWindow = 0.5f * clusterCut; + const float nSigmaCut = vc.nSigmaCut; + const float duplicateZCut = vc.duplicateZCut > 0.f + ? vc.duplicateZCut + : std::max(4.f * vc.pairCut, 0.5f * clusterCut); + const int minContributors = vc.clusterContributorsCut; + const int suppressLowMultDebris = vc.suppressLowMultDebris; + const bool incremental = vc.incrementalSeeding; + const float fineZWindow = vc.fineZWindow; + const bool doFine = fineZWindow > 0.f && fineZWindow < zWindow; + const int fineMinDensity = vc.fineMinDensity; + const float fineMaxDrift = vc.fineMaxDrift; + const float goodSig = vc.goodContributorsSignificance; + const float duplicateZScale = vc.duplicateZScale; + const float goodLineChi2Cut = vc.goodLineChi2Cut; + const float goodLinePtCut = std::abs(mBz) > 0.01f ? vc.goodLinePtCut : -1.f; + + const auto processROF = [&](const int rofId) { + if (skipROF(iteration, rofId)) { + return; + } + auto& lines = mTimeFrame->getLines(rofId); + const int n = static_cast(lines.size()); + if (n < 2) { + return; + } + const auto lineSpan = gsl::span(lines.data(), lines.size()); + + // project every line onto the beam line + bounded_vector zc(mMemoryPool.get()); + bounded_vector tc(mMemoryPool.get()); + bounded_vector lic(mMemoryPool.get()); + zc.reserve(n); + tc.reserve(n); + lic.reserve(n); + for (int i = 0; i < n; ++i) { + const auto& line = lines[i]; + const float dx = line.originPoint(0) - beamX; + const float dy = line.originPoint(1) - beamY; + const float ux = line.cosinesDirector(0); + const float uy = line.cosinesDirector(1); + const float uz = line.cosinesDirector(2); + const float den = ux * ux + uy * uy; + if (den < constants::Tolerance) { + continue; + } + const float s0 = -(dx * ux + dy * uy) / den; + const float z = line.originPoint(2) + s0 * uz; + if (!(o2::gpu::CAMath::Abs(z) < maxZ)) { + continue; + } + const auto sym = line.mTime.makeSymmetrical(); + zc.push_back(z); + tc.push_back(LineTime{sym.getTimeStamp(), sym.getTimeStampError()}); + lic.push_back(i); + } + const int m = static_cast(zc.size()); + if (m < 2) { + return; + } + + // sort by z + bounded_vector order(m, mMemoryPool.get()); + std::iota(order.begin(), order.end(), 0); + std::sort(order.begin(), order.end(), [&](const int a, const int b) { return zc[a] < zc[b]; }); + bounded_vector Z(m, mMemoryPool.get()); + bounded_vector T(m, mMemoryPool.get()); + bounded_vector LI(m, mMemoryPool.get()); + for (int i = 0; i < m; ++i) { + Z[i] = zc[order[i]]; + T[i] = tc[order[i]]; + LI[i] = lic[order[i]]; + } + + // density of time-compatible neighbours in a z window + bounded_vector density(m, mMemoryPool.get()); + bounded_vector win(m, mMemoryPool.get()); + computeDensities(Z, T, zWindow, incremental, win, density); + bounded_vector densityFine(doFine ? m : 0, mMemoryPool.get()); + bounded_vector winFine(doFine ? m : 0, mMemoryPool.get()); + if (doFine) { + computeDensities(Z, T, fineZWindow, incremental, winFine, densityFine); + } + + // peaks: leftmost density maxima + bounded_vector isMax(m, mMemoryPool.get()); + markLeftmostMaxima(density, Z, win, incremental, mMemoryPool.get(), isMax); + bounded_vector isMaxFine(doFine ? m : 0, mMemoryPool.get()); + if (doFine) { + markLeftmostMaxima(densityFine, Z, winFine, incremental, mMemoryPool.get(), isMaxFine); + } + bounded_vector peaks(mMemoryPool.get()); + bounded_vector isPeakFine(m, 0, mMemoryPool.get()); + for (int k = 0; k < m; ++k) { + bool peak = density[k] >= 2 && isMax[k]; + if (!peak && doFine && densityFine[k] >= fineMinDensity && isMaxFine[k]) { + peak = true; + isPeakFine[k] = 1; + } + if (peak) { + peaks.push_back(k); + } + } + const int np = static_cast(peaks.size()); + if (np == 0) { + return; + } + + // one candidate per peak: seed fit over the window, prune outliers, refit, apply the cuts + bounded_vector cand(np, mMemoryPool.get()); + bounded_vector ok(np, mMemoryPool.get()); + bounded_vector nGoodCand(np, 0, mMemoryPool.get()); + bounded_vector fineCand(np, 0, mMemoryPool.get()); + const auto& linesQuality = mTimeFrame->getLinesQuality(rofId); + const int nQual = static_cast(linesQuality.size()); + for (int p = 0; p < np; p++) { + const int k = peaks[p]; + fineCand[p] = isPeakFine[k]; + bounded_vector members(mMemoryPool.get()); + members.reserve(density[k]); + for (int j = win[k].lo; j < win[k].hi; ++j) { + if (timeCompatible(T[k], T[j])) { + members.push_back(LI[j]); + } + } + if (members.size() < 2) { + ok[p] = 0; + continue; + } + ClusterLines seed{gsl::span{members.data(), members.size()}, lineSpan}; + if (!seed.isValid()) { + ok[p] = 0; + continue; + } + bounded_vector kept(mMemoryPool.get()); + kept.reserve(members.size()); + for (const int idx : members) { + if (Line::getDistance2FromPoint(lines[idx], seed.getVertex()) < pairCut2) { + kept.push_back(idx); + } + } + if (kept.size() < 2) { + ok[p] = 0; + continue; + } + int nGood = 0; + if (nQual == 0) { + nGood = static_cast(kept.size()); + } else { + for (const int idx : kept) { + const float chi2 = idx < nQual ? linesQuality[idx].chi2 : -1.f; + const float pt = idx < nQual ? linesQuality[idx].pt : -1.f; + const bool okChi2 = goodLineChi2Cut <= 0.f || chi2 <= goodLineChi2Cut; + const bool okPt = goodLinePtCut <= 0.f || pt >= goodLinePtCut; + nGood += okChi2 && okPt; + } + } + nGoodCand[p] = nGood; + ClusterLines fit{gsl::span{kept.data(), kept.size()}, lineSpan}; + const float bd2 = (beamX - fit.getVertex()[0]) * (beamX - fit.getVertex()[0]) + + (beamY - fit.getVertex()[1]) * (beamY - fit.getVertex()[1]); + if (!fit.isValid() || static_cast(fit.getSize()) < minContributors || !(bd2 < nSigmaCut)) { + ok[p] = 0; + continue; + } + if (fineMaxDrift > 0.f && fineCand[p] && std::abs(fit.getVertex()[2] - seed.getVertex()[2]) > fineMaxDrift) { + ok[p] = 0; + continue; + } + cand[p] = std::move(fit); + ok[p] = 1; + } + + // duplicate suppression + bounded_vector keep(np, mMemoryPool.get()); + for (int p = 0; p < np; ++p) { + if (!ok[p]) { + keep[p] = 0; + continue; + } + uint8_t survive = 1; + const float zp = cand[p].getVertex()[2]; + const auto sp = cand[p].getSize(); + const float radius = (duplicateZScale > 0.f && sp > 0) + ? duplicateZScale / std::sqrt(static_cast(sp)) + : duplicateZCut; + for (int q = 0; q < np && survive; ++q) { + if (q == p || !ok[q]) { + continue; + } + if (!cand[p].getTimeStamp().isCompatible(cand[q].getTimeStamp())) { + continue; + } + if (o2::gpu::GPUCommonMath::Abs(zp - cand[q].getVertex()[2]) < radius) { + const auto sq = cand[q].getSize(); + if (sq > sp || (sq == sp && q < p)) { + survive = 0; + } + } + } + keep[p] = survive; + } + + // emit + bounded_vector accepted(mMemoryPool.get()); + for (int p = 0; p < np; ++p) { + if (keep[p]) { + accepted.push_back(p); + } + } + std::sort(accepted.begin(), accepted.end(), + [&](int a, int b) { return cand[a].getSize() > cand[b].getSize(); }); + double rofLoad = 0.; + if (goodSig > 0.f) { + for (int p = 0; p < np; ++p) { + if (ok[p] && !fineCand[p]) { + rofLoad += cand[p].getSize(); + } + } + } + const float sigThreshold = goodSig > 0.f ? goodSig * std::sqrt(static_cast(std::max(rofLoad, 1.))) : 0.f; + for (const int p : accepted) { + if (!rofVertices[rofId].empty()) { + if (goodSig > 0.f) { + if (nGoodCand[p] <= sigThreshold) { + continue; + } + } else if (static_cast(cand[p].getSize()) < suppressLowMultDebris) { + continue; + } + } + Vertex vertex{cand[p].getVertex().data(), + cand[p].getRMS2(), + (ushort)cand[p].getSize(), + cand[p].getAvgDistance2()}; + vertex.setTimeStamp(cand[p].getTimeStamp()); + rofVertices[rofId].push_back(vertex); + if (withMC) { + auto& lineLabels = mTimeFrame->getLinesLabel(rofId); + const int nLab = static_cast(lineLabels.size()); + bounded_vector labels(mMemoryPool.get()); + for (const auto idx : cand[p].getLabels()) { + if (idx < 0 || idx >= nLab) { + LOGP(error, "[seedDbg] OOB lineLabel idx={} nLab={} rof={}", idx, nLab, rofId); + continue; + } + labels.push_back(lineLabels[idx]); + } + rofLabels[rofId].push_back(computeMain(labels)); + } + } + }; + + if (mTaskArena->max_concurrency() <= 1) { + for (int rofId{0}; rofId < nRofs; ++rofId) { + processROF(rofId); + } + } else { + mTaskArena->execute([&] { + tbb::parallel_for(0, nRofs, [&](const int rofId) { + processROF(rofId); + }); + }); + } + for (int rofId{0}; rofId < nRofs; ++rofId) { + for (auto& vertex : rofVertices[rofId]) { + mTimeFrame->addPrimaryVertex(vertex); + } + if (withMC) { + for (auto& label : rofLabels[rofId]) { + mTimeFrame->addPrimaryVertexLabel(label); + } + } + } + + auto& pvs = mTimeFrame->getPrimaryVertices(); + bounded_vector indices(pvs.size(), mMemoryPool.get()); + std::iota(indices.begin(), indices.end(), 0); + std::sort(indices.begin(), indices.end(), [&pvs](const size_t i, const size_t j) { + const auto aLower = pvs[i].getTimeStamp().lower(); + const auto bLower = pvs[j].getTimeStamp().lower(); + if (aLower != bLower) { + return aLower < bLower; + } + return pvs[i].getNContributors() > pvs[j].getNContributors(); + }); + bounded_vector sortedVtx(pvs.get_allocator()); + sortedVtx.reserve(pvs.size()); + for (const size_t idx : indices) { + sortedVtx.push_back(pvs[idx]); + } + pvs.swap(sortedVtx); + if (withMC) { + auto& mc = mTimeFrame->getPrimaryVerticesLabels(); + bounded_vector sortedMC(mc.get_allocator()); + sortedMC.reserve(mc.size()); + for (const size_t idx : indices) { + sortedMC.push_back(mc[idx]); + } + mc.swap(sortedMC); + } + mTimeFrame->updateROFVertexLookupTable(); +} + template void TrackerTraits::computeLayerCells(const int iteration) { @@ -302,6 +866,14 @@ void TrackerTraits::computeLayerCells(const int iteration) const auto& cellTopology = topology.getCell(cellTopologyId); const auto& firstLink = topology.getLink(cellTopology.firstLink); const auto& secondLink = topology.getLink(cellTopology.secondLink); + const float cellDeltaPhiCut = + mTrkParams[iteration].PassFlags[IterationStep::SeedingVertexPass] + ? cellDeltaPhiBound(mBz, mTrkParams[iteration].CellDeltaPhiMinPt, + mTrkParams[iteration].LayerRadii[firstLink.fromLayer], + mTrkParams[iteration].LayerRadii[firstLink.toLayer], + mTrkParams[iteration].LayerRadii[secondLink.toLayer], + mTimeFrame->getLinkMSAngle(cellTopology.firstLink)) + : -1.f; const Tracklet& currentTracklet{mTimeFrame->getTracklets()[cellTopology.firstLink][iTracklet]}; const int nextLayerClusterIndex{currentTracklet.secondClusterIndex}; const int nextLayerFirstTrackletIndex{mTimeFrame->getTrackletsLookupTable()[cellTopology.secondLink][nextLayerClusterIndex]}; @@ -314,9 +886,21 @@ void TrackerTraits::computeLayerCells(const int iteration) if (!currentTracklet.getTimeStamp().isCompatible(nextTracklet.getTimeStamp())) { continue; } + if (cellDeltaPhiCut > 0.f) { + float deltaPhi{std::abs(currentTracklet.phi - nextTracklet.phi)}; + if (deltaPhi > o2::constants::math::PI) { // the two directions straddle the +-pi wrap + deltaPhi = o2::constants::math::TwoPI - deltaPhi; + } + if (deltaPhi > cellDeltaPhiCut) { + continue; + } + } const float deltaTanLambdaSigma = std::abs(currentTracklet.tanLambda - nextTracklet.tanLambda) / mTrkParams[iteration].CellDeltaTanLambdaSigma; - if (deltaTanLambdaSigma < mTrkParams[iteration].NSigmaCut) { + const float cellTanLNSigma = mTrkParams[iteration].CellDeltaTanLambdaNSigma > 0.f + ? mTrkParams[iteration].CellDeltaTanLambdaNSigma + : mTrkParams[iteration].NSigmaCut; + if (deltaTanLambdaSigma < cellTanLNSigma) { /// Track seed preparation. Clusters are numbered progressively from the innermost going outward. const int clusId[3]{ diff --git a/Detectors/ITSMFT/ITS/tracking/src/TrackingInterface.cxx b/Detectors/ITSMFT/ITS/tracking/src/TrackingInterface.cxx index 83a1086ec5263..331bc45242860 100644 --- a/Detectors/ITSMFT/ITS/tracking/src/TrackingInterface.cxx +++ b/Detectors/ITSMFT/ITS/tracking/src/TrackingInterface.cxx @@ -11,8 +11,14 @@ #include #include +#include #include #include +#include +#include +#include +#include +#include #include @@ -35,6 +41,7 @@ #include "Framework/InputRecordWalker.h" #include "Framework/DataRefUtils.h" #include "Framework/DeviceSpec.h" +#include "SimulationDataFormat/MCTrack.h" using namespace o2::framework; using namespace o2::its; @@ -51,6 +58,33 @@ void ITSTrackingInterface::initialise() auto trackParams = TrackingMode::getTrackingParameters(mMode); auto vertParams = TrackingMode::getVertexingParameters(mMode); overrideParameters(trackParams, vertParams); + if (trackConf.seedingVertexIteration && !trackParams.empty()) { + + TrackingParameters seedingPass = trackParams.front(); + seedingPass.PassFlags = IterationSteps{IterationStep::FirstPass, IterationStep::RebuildClusterLUT, + IterationStep::ResetVertices, IterationStep::SeedingVertexPass}; + seedingPass.PerPrimaryVertexProcessing = false; + seedingPass.UseDiamond = true; // PV-independent + seedingPass.NLayers = 3; // 3-layer {0,1,2} vertexing topology + seedingPass.MinTrackLength = 3; // only trackleting+celling on the 3 inner layers + seedingPass.CreateArtefactLabels = mIsMC; + seedingPass.ZBins = vertConf.ZBins; + seedingPass.PhiBins = vertConf.PhiBins; + seedingPass.Diamond[0] = trackConf.diamondPos[0]; + seedingPass.Diamond[1] = trackConf.diamondPos[1]; + seedingPass.Diamond[2] = trackConf.diamondPos[2]; + seedingPass.NSigmaCut = trackConf.diamondTrackletingNSigmaCut; + seedingPass.PVres = trackConf.diamondTrackletingPVres; + seedingPass.CellDeltaTanLambdaSigma = trackConf.diamondTrackletingCellDeltaTanLambdaSigma; + seedingPass.CellDeltaTanLambdaNSigma = trackConf.diamondCellTanLambdaNSigma; + seedingPass.CellDeltaPhiMinPt = trackConf.diamondTrackletingCellDeltaPhiMinPt; + seedingPass.CellLineSharedClusterCut = trackConf.cellLineSharedClusterCut; + trackParams.push_back(seedingPass); + LOGP(info, "Appended a seeding-vertex pass slot (stub sub-steps; mVertexer still active): TrackletMinPt={:.4f} NSigmaCut={:.4f} PVres={:.4f} CellDeltaTanLambdaSigma={:.6f} CellDeltaTanLambdaNSigma={:.4f} CellDeltaPhiMinPt={:.4f} CellLineSharedClusterCut={} ZBins={} PhiBins={}", + seedingPass.TrackletMinPt, seedingPass.NSigmaCut, seedingPass.PVres, + seedingPass.CellDeltaTanLambdaSigma, seedingPass.CellDeltaTanLambdaNSigma, seedingPass.CellDeltaPhiMinPt, seedingPass.CellLineSharedClusterCut, + seedingPass.ZBins, seedingPass.PhiBins); + } LOGP(info, "Initializing tracker in {} phase reconstruction with {} passes for tracking and {}/{} for vertexing", TrackingMode::toString(mMode), trackParams.size(), o2::its::VertexerParamConfig::Instance().nIterations, vertParams.size()); mTracker->setParameters(trackParams); mVertexer->setParameters(vertParams); @@ -221,8 +255,12 @@ void ITSTrackingInterface::run(framework::ProcessingContext& pc) float vertexerElapsedTime{0.f}, trackerElapsedTime{0.f}; if (mRunVertexer) { - // Run seeding vertexer - vertexerElapsedTime = mVertexer->clustersToVertices(logger); + // Run seeding vertexer. With seedingVertexIteration the tracker-owned seeding phase replaces the + // standalone vertexer (Option A): same call site, so the consumer block below is unchanged. Until + // the seeding sub-steps are implemented it produces no vertices, so keep the flag off in production. + vertexerElapsedTime = o2::its::TrackerParamConfig::Instance().seedingVertexIteration + ? mTracker->clustersToVertices(logger) + : mVertexer->clustersToVertices(logger); const auto& vtx = mTimeFrame->getPrimaryVertices(); vertices.insert(vertices.begin(), vtx.begin(), vtx.end()); if (mIsMC) {