diff --git a/DataFormats/Reconstruction/include/ReconstructionDataFormats/Vertex.h b/DataFormats/Reconstruction/include/ReconstructionDataFormats/Vertex.h index bbc2c01359276..7f9332f9447d4 100644 --- a/DataFormats/Reconstruction/include/ReconstructionDataFormats/Vertex.h +++ b/DataFormats/Reconstruction/include/ReconstructionDataFormats/Vertex.h @@ -146,8 +146,8 @@ class Vertex : public VertexBase GPUd() void setChi2(float v) { mChi2 = v; } GPUd() float getChi2() const { return mChi2; } - GPUd() const Stamp& getTimeStamp() const { return mTimeStamp; } - GPUd() Stamp& getTimeStamp() { return mTimeStamp; } + GPUhd() const Stamp& getTimeStamp() const { return mTimeStamp; } + GPUhd() Stamp& getTimeStamp() { return mTimeStamp; } GPUd() void setTimeStamp(const Stamp& v) { mTimeStamp = v; } protected: diff --git a/DataFormats/common/include/CommonDataFormat/TimeStamp.h b/DataFormats/common/include/CommonDataFormat/TimeStamp.h index 90019de5122be..354c937c1a24b 100644 --- a/DataFormats/common/include/CommonDataFormat/TimeStamp.h +++ b/DataFormats/common/include/CommonDataFormat/TimeStamp.h @@ -28,7 +28,7 @@ class TimeStamp GPUhdDefault() TimeStamp() CON_DEFAULT; GPUhdDefault() ~TimeStamp() CON_DEFAULT; GPUdi() TimeStamp(T time) { mTimeStamp = time; } - GPUdi() T getTimeStamp() const { return mTimeStamp; } + GPUhdi() T getTimeStamp() const { return mTimeStamp; } GPUdi() void setTimeStamp(T t) { mTimeStamp = t; } GPUdi() bool operator==(const TimeStamp& t) const { return mTimeStamp == t.mTimeStamp; } diff --git a/Detectors/ITSMFT/ITS/tracking/GPU/ITStrackingGPU/TimeFrameGPU.h b/Detectors/ITSMFT/ITS/tracking/GPU/ITStrackingGPU/TimeFrameGPU.h index 2a2fadb21caeb..a67ef03e0ef22 100644 --- a/Detectors/ITSMFT/ITS/tracking/GPU/ITStrackingGPU/TimeFrameGPU.h +++ b/Detectors/ITSMFT/ITS/tracking/GPU/ITStrackingGPU/TimeFrameGPU.h @@ -210,7 +210,7 @@ class TimeFrameGPU : public TimeFrame /// interface int getNClustersInRofSpan(const int, const int, const int) const; IndexTableUtils* getDeviceIndexTableUtils() { return mIndexTableUtilsDevice; } - int* getDeviceROframesClusters(const int layer) { return mROframesClustersDevice[layer]; } + int* getDeviceROFramesClusters(const int layer) { return mROFramesClustersDevice[layer]; } std::vector>& getVerticesInChunks() { return mVerticesInChunks; } std::vector>& getNVerticesInChunks() { return mNVerticesInChunks; } std::vector& getTrackITSExt() { return mTrackITSExt; } @@ -218,7 +218,7 @@ class TimeFrameGPU : public TimeFrame int getNAllocatedROFs() const { return mNrof; } // Allocated means maximum nROF for each chunk while populated is the number of loaded ones. StaticTrackingParameters* getDeviceTrackingParameters() { return mTrackingParamsDevice; } Vertex* getDeviceVertices() { return mVerticesDevice; } - int* getDeviceROframesPV() { return mROframesPVDevice; } + int* getDeviceROFramesPV() { return mROFramesPVDevice; } unsigned char* getDeviceUsedClusters(const int); const o2::base::Propagator* getChainPropagator(); @@ -251,10 +251,10 @@ class TimeFrameGPU : public TimeFrame // Device pointers StaticTrackingParameters* mTrackingParamsDevice; IndexTableUtils* mIndexTableUtilsDevice; - std::array mROframesClustersDevice; + std::array mROFramesClustersDevice; std::array mUsedClustersDevice; Vertex* mVerticesDevice; - int* mROframesPVDevice; + int* mROFramesPVDevice; // Hybrid pref std::array mClustersDevice; @@ -314,7 +314,7 @@ size_t TimeFrameGPU::loadChunkData(const size_t chunk, const size_t off template inline int TimeFrameGPU::getNClustersInRofSpan(const int rofIdstart, const int rofSpanSize, const int layerId) const { - return static_cast(mROframesClusters[layerId][(rofIdstart + rofSpanSize) < mROframesClusters.size() ? rofIdstart + rofSpanSize : mROframesClusters.size() - 1] - mROframesClusters[layerId][rofIdstart]); + return static_cast(mROFramesClusters[layerId][(rofIdstart + rofSpanSize) < mROFramesClusters.size() ? rofIdstart + rofSpanSize : mROFramesClusters.size() - 1] - mROFramesClusters[layerId][rofIdstart]); } } // namespace gpu } // namespace its diff --git a/Detectors/ITSMFT/ITS/tracking/GPU/cuda/TimeFrameGPU.cu b/Detectors/ITSMFT/ITS/tracking/GPU/cuda/TimeFrameGPU.cu index e758bf4990c4f..1d0fca2585730 100644 --- a/Detectors/ITSMFT/ITS/tracking/GPU/cuda/TimeFrameGPU.cu +++ b/Detectors/ITSMFT/ITS/tracking/GPU/cuda/TimeFrameGPU.cu @@ -443,21 +443,21 @@ void TimeFrameGPU::initDevice(const int chunks, mMemChunks[iChunk].allocate(GpuTimeFrameChunk::computeRofPerChunk(mGpuParams, mAvailMemGB), mGpuStreams[iChunk]); } for (auto iLayer{0}; iLayer < nLayers; ++iLayer) { - checkGPUError(cudaMalloc(reinterpret_cast(&mROframesClustersDevice[iLayer]), mROframesClusters[iLayer].size() * sizeof(int))); + checkGPUError(cudaMalloc(reinterpret_cast(&mROFramesClustersDevice[iLayer]), mROFramesClusters[iLayer].size() * sizeof(int))); checkGPUError(cudaMalloc(reinterpret_cast(&(mUsedClustersDevice[iLayer])), sizeof(unsigned char) * mGpuParams.clustersPerROfCapacity * mNrof)); } checkGPUError(cudaMalloc(reinterpret_cast(&mVerticesDevice), sizeof(Vertex) * mGpuParams.maxVerticesCapacity)); - checkGPUError(cudaMalloc(reinterpret_cast(&mROframesPVDevice), sizeof(int) * (mNrof + 1))); + checkGPUError(cudaMalloc(reinterpret_cast(&mROFramesPVDevice), sizeof(int) * (mNrof + 1))); mFirstInit = false; } if (maxLayers < nLayers) { // Vertexer for (auto iLayer{0}; iLayer < nLayers; ++iLayer) { - checkGPUError(cudaMemcpy(mROframesClustersDevice[iLayer], mROframesClusters[iLayer].data(), mROframesClusters[iLayer].size() * sizeof(int), cudaMemcpyHostToDevice)); + checkGPUError(cudaMemcpy(mROFramesClustersDevice[iLayer], mROFramesClusters[iLayer].data(), mROFramesClusters[iLayer].size() * sizeof(int), cudaMemcpyHostToDevice)); } } else { // Tracker checkGPUError(cudaMemcpy(mVerticesDevice, mPrimaryVertices.data(), sizeof(Vertex) * mPrimaryVertices.size(), cudaMemcpyHostToDevice)); - checkGPUError(cudaMemcpy(mROframesPVDevice, mROframesPV.data(), sizeof(int) * mROframesPV.size(), cudaMemcpyHostToDevice)); + checkGPUError(cudaMemcpy(mROFramesPVDevice, mROFramesPV.data(), sizeof(int) * mROFramesPV.size(), cudaMemcpyHostToDevice)); if (!iteration) { for (auto iLayer{0}; iLayer < nLayers; ++iLayer) { checkGPUError(cudaMemset(mUsedClustersDevice[iLayer], 0, sizeof(unsigned char) * mGpuParams.clustersPerROfCapacity * mNrof)); diff --git a/Detectors/ITSMFT/ITS/tracking/GPU/cuda/VertexerTraitsGPU.cu b/Detectors/ITSMFT/ITS/tracking/GPU/cuda/VertexerTraitsGPU.cu index 9a1ed507ae5a4..f121609f9b019 100644 --- a/Detectors/ITSMFT/ITS/tracking/GPU/cuda/VertexerTraitsGPU.cu +++ b/Detectors/ITSMFT/ITS/tracking/GPU/cuda/VertexerTraitsGPU.cu @@ -627,8 +627,8 @@ void VertexerTraitsGPU::computeTracklets() gpu::trackleterKernelMultipleRof<<getStream(chunkId).get()>>>( mTimeFrameGPU->getChunk(chunkId).getDeviceClusters(0), // const Cluster* clustersNextLayer, // 0 2 mTimeFrameGPU->getChunk(chunkId).getDeviceClusters(1), // const Cluster* clustersCurrentLayer, // 1 1 - mTimeFrameGPU->getDeviceROframesClusters(0), // const int* sizeNextLClusters, - mTimeFrameGPU->getDeviceROframesClusters(1), // const int* sizeCurrentLClusters, + mTimeFrameGPU->getDeviceROFramesClusters(0), // const int* sizeNextLClusters, + mTimeFrameGPU->getDeviceROFramesClusters(1), // const int* sizeCurrentLClusters, mTimeFrameGPU->getChunk(chunkId).getDeviceIndexTables(0), // const int* nextIndexTables, mTimeFrameGPU->getChunk(chunkId).getDeviceTracklets(0), // Tracklet* Tracklets, mTimeFrameGPU->getChunk(chunkId).getDeviceNTrackletCluster(0), // int* foundTracklets, @@ -641,8 +641,8 @@ void VertexerTraitsGPU::computeTracklets() gpu::trackleterKernelMultipleRof<<getStream(chunkId).get()>>>( mTimeFrameGPU->getChunk(chunkId).getDeviceClusters(2), // const Cluster* clustersNextLayer, // 0 2 mTimeFrameGPU->getChunk(chunkId).getDeviceClusters(1), // const Cluster* clustersCurrentLayer, // 1 1 - mTimeFrameGPU->getDeviceROframesClusters(2), // const int* sizeNextLClusters, - mTimeFrameGPU->getDeviceROframesClusters(1), // const int* sizeCurrentLClusters, + mTimeFrameGPU->getDeviceROFramesClusters(2), // const int* sizeNextLClusters, + mTimeFrameGPU->getDeviceROFramesClusters(1), // const int* sizeCurrentLClusters, mTimeFrameGPU->getChunk(chunkId).getDeviceIndexTables(2), // const int* nextIndexTables, mTimeFrameGPU->getChunk(chunkId).getDeviceTracklets(1), // Tracklet* Tracklets, mTimeFrameGPU->getChunk(chunkId).getDeviceNTrackletCluster(1), // int* foundTracklets, @@ -655,8 +655,8 @@ void VertexerTraitsGPU::computeTracklets() gpu::trackletSelectionKernelMultipleRof<<getStream(chunkId).get()>>>( mTimeFrameGPU->getChunk(chunkId).getDeviceClusters(0), // const Cluster* clusters0, // Clusters on layer 0 mTimeFrameGPU->getChunk(chunkId).getDeviceClusters(1), // const Cluster* clusters1, // Clusters on layer 1 - mTimeFrameGPU->getDeviceROframesClusters(0), // const int* sizeClustersL0, // Number of clusters on layer 0 per ROF - mTimeFrameGPU->getDeviceROframesClusters(1), // const int* sizeClustersL1, // Number of clusters on layer 1 per ROF + mTimeFrameGPU->getDeviceROFramesClusters(0), // const int* sizeClustersL0, // Number of clusters on layer 0 per ROF + mTimeFrameGPU->getDeviceROFramesClusters(1), // const int* sizeClustersL1, // Number of clusters on layer 1 per ROF mTimeFrameGPU->getChunk(chunkId).getDeviceTracklets(0), // Tracklet* tracklets01, // Tracklets on layer 0-1 mTimeFrameGPU->getChunk(chunkId).getDeviceTracklets(1), // Tracklet* tracklets12, // Tracklets on layer 1-2 mTimeFrameGPU->getChunk(chunkId).getDeviceNTrackletCluster(0), // const int* nFoundTracklets01, // Number of tracklets found on layers 0-1 @@ -688,8 +688,8 @@ void VertexerTraitsGPU::computeTracklets() gpu::trackletSelectionKernelMultipleRof<<getStream(chunkId).get()>>>( mTimeFrameGPU->getChunk(chunkId).getDeviceClusters(0), // const Cluster* clusters0, // Clusters on layer 0 mTimeFrameGPU->getChunk(chunkId).getDeviceClusters(1), // const Cluster* clusters1, // Clusters on layer 1 - mTimeFrameGPU->getDeviceROframesClusters(0), // const int* sizeClustersL0, // Number of clusters on layer 0 per ROF - mTimeFrameGPU->getDeviceROframesClusters(1), // const int* sizeClustersL1, // Number of clusters on layer 1 per ROF + mTimeFrameGPU->getDeviceROFramesClusters(0), // const int* sizeClustersL0, // Number of clusters on layer 0 per ROF + mTimeFrameGPU->getDeviceROFramesClusters(1), // const int* sizeClustersL1, // Number of clusters on layer 1 per ROF mTimeFrameGPU->getChunk(chunkId).getDeviceTracklets(0), // Tracklet* tracklets01, // Tracklets on layer 0-1 mTimeFrameGPU->getChunk(chunkId).getDeviceTracklets(1), // Tracklet* tracklets12, // Tracklets on layer 1-2 mTimeFrameGPU->getChunk(chunkId).getDeviceNTrackletCluster(0), // const int* nFoundTracklets01, // Number of tracklets found on layers 0-1 @@ -723,8 +723,8 @@ void VertexerTraitsGPU::computeTracklets() std::vector usedLines; for (int rofId{0}; rofId < rofs; ++rofId) { auto rof = offset + rofId; - auto clustersL1offsetRof = mTimeFrameGPU->getROframeClusters(1)[rof] - mTimeFrameGPU->getROframeClusters(1)[offset]; // starting cluster offset for this ROF - auto nClustersL1Rof = mTimeFrameGPU->getROframeClusters(1)[rof + 1] - mTimeFrameGPU->getROframeClusters(1)[rof]; // number of clusters for this ROF + auto clustersL1offsetRof = mTimeFrameGPU->getROFrameClusters(1)[rof] - mTimeFrameGPU->getROFrameClusters(1)[offset]; // starting cluster offset for this ROF + auto nClustersL1Rof = mTimeFrameGPU->getROFrameClusters(1)[rof + 1] - mTimeFrameGPU->getROFrameClusters(1)[rof]; // number of clusters for this ROF auto linesOffsetRof = exclusiveFoundLinesHost[clustersL1offsetRof]; // starting line offset for this ROF auto nLinesRof = exclusiveFoundLinesHost[clustersL1offsetRof + nClustersL1Rof] - linesOffsetRof; gsl::span linesInRof(lines.data() + linesOffsetRof, static_cast::size_type>(nLinesRof)); @@ -756,7 +756,7 @@ void VertexerTraitsGPU::computeTracklets() int start{0}; for (int rofId{0}; rofId < mTimeFrameGPU->getNVerticesInChunks()[chunkId].size(); ++rofId) { gsl::span rofVerts{mTimeFrameGPU->getVerticesInChunks()[chunkId].data() + start, static_cast::size_type>(mTimeFrameGPU->getNVerticesInChunks()[chunkId][rofId])}; - mTimeFrameGPU->addPrimaryVertices(rofVerts); + mTimeFrameGPU->addPrimaryVertices(rofVerts, rofId); if (mTimeFrameGPU->hasMCinformation()) { mTimeFrameGPU->getVerticesLabels().emplace_back(); // TODO: add MC labels diff --git a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/ClusterLines.h b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/ClusterLines.h index 550b8405770dc..8a4a2193ca2ef 100644 --- a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/ClusterLines.h +++ b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/ClusterLines.h @@ -9,8 +9,8 @@ // granted to it by virtue of its status as an Intergovernmental Organization // or submit itself to any jurisdiction. -#ifndef O2_ITSMFT_TRACKING_LINE_H_ -#define O2_ITSMFT_TRACKING_LINE_H_ +#ifndef O2_ITS_CLUSTERLINES_H +#define O2_ITS_CLUSTERLINES_H #include #include @@ -19,11 +19,8 @@ #include "ITStracking/Tracklet.h" #include "GPUCommonMath.h" -namespace o2 +namespace o2::its { -namespace its -{ - struct Line final { GPUhd() Line(); GPUhd() Line(const Line&); @@ -31,7 +28,7 @@ struct Line final { GPUhd() Line(const float firstPoint[3], const float secondPoint[3]); GPUhd() Line(const Tracklet&, const Cluster*, const Cluster*); - inline static float getDistanceFromPoint(const Line& line, const std::array& point); + static float getDistanceFromPoint(const Line& line, const std::array& point); GPUhd() static float getDistanceFromPoint(const Line& line, const float point[3]); static std::array getDCAComponents(const Line& line, const std::array point); GPUhd() static void getDCAComponents(const Line& line, const float point[3], float destArray[6]); @@ -39,12 +36,14 @@ struct Line final { static bool areParallel(const Line&, const Line&, const float precision = 1e-14); GPUhd() unsigned char isEmpty() const { return (originPoint[0] == 0.f && originPoint[1] == 0.f && originPoint[2] == 0.f) && (cosinesDirector[0] == 0.f && cosinesDirector[1] == 0.f && cosinesDirector[2] == 0.f); } + GPUhdi() auto getDeltaROF() const { return rof[1] - rof[0]; } GPUhd() void print() const; bool operator==(const Line&) const; bool operator!=(const Line&) const; + short getMinROF() const { return rof[0] < rof[1] ? rof[0] : rof[1]; } - float originPoint[3], cosinesDirector[3]; // std::array originPoint, cosinesDirector; - float weightMatrix[6] = {1., 0., 0., 1., 0., 1.}; // std::array weightMatrix; + float originPoint[3], cosinesDirector[3]; + float weightMatrix[6] = {1., 0., 0., 1., 0., 1.}; // weightMatrix is a symmetric matrix internally stored as // 0 --> row = 0, col = 0 // 1 --> 0,1 @@ -52,11 +51,13 @@ struct Line final { // 3 --> 1,1 // 4 --> 1,2 // 5 --> 2,2 - // Debug quantities + short rof[2]; }; GPUhdi() Line::Line() : weightMatrix{1., 0., 0., 1., 0., 1.} { + rof[0] = 0; + rof[1] = 0; } GPUhdi() Line::Line(const Line& other) @@ -68,6 +69,9 @@ GPUhdi() Line::Line(const Line& other) for (int i{0}; i < 6; ++i) { weightMatrix[i] = other.weightMatrix[i]; } + for (int i{0}; i < 2; ++i) { + rof[i] = other.rof[i]; + } } GPUhdi() Line::Line(const float firstPoint[3], const float secondPoint[3]) @@ -83,6 +87,9 @@ GPUhdi() Line::Line(const float firstPoint[3], const float secondPoint[3]) for (int index{0}; index < 3; ++index) { cosinesDirector[index] *= inverseNorm; } + + rof[0] = 0; + rof[1] = 0; } GPUhdi() Line::Line(const Tracklet& tracklet, const Cluster* innerClusters, const Cluster* outerClusters) @@ -101,6 +108,9 @@ GPUhdi() Line::Line(const Tracklet& tracklet, const Cluster* innerClusters, cons for (int index{0}; index < 3; ++index) { cosinesDirector[index] *= inverseNorm; } + + rof[0] = tracklet.rof[0]; + rof[1] = tracklet.rof[1]; } // static functions: @@ -198,8 +208,8 @@ inline bool Line::operator!=(const Line& rhs) const GPUhdi() void Line::print() const { - printf("Line: originPoint = (%f, %f, %f), cosinesDirector = (%f, %f, %f)\n", - originPoint[0], originPoint[1], originPoint[2], cosinesDirector[0], cosinesDirector[1], cosinesDirector[2]); + printf("Line: originPoint = (%f, %f, %f), cosinesDirector = (%f, %f, %f), rofs = (%u, %u)\n", + originPoint[0], originPoint[1], originPoint[2], cosinesDirector[0], cosinesDirector[1], cosinesDirector[2], rof[0], rof[1]); } class ClusterLines final @@ -211,11 +221,13 @@ class ClusterLines final ClusterLines(const Line& firstLine, const Line& secondLine); void add(const int& lineLabel, const Line& line, const bool& weight = false); void computeClusterCentroid(); + void updateROFPoll(const Line&); inline std::vector& getLabels() { return mLabels; } inline int getSize() const { return mLabels.size(); } + inline short getROF() const { return mROF; } inline std::array getVertex() const { return mVertex; } inline std::array getRMS2() const { return mRMS2; } inline float getAvgDistance2() const { return mAvgDistance2; } @@ -230,8 +242,9 @@ class ClusterLines final std::array mVertex = {0.f}; // cluster centroid position std::array mRMS2 = {0.f}; // symmetric matrix: diagonal is RMS2 float mAvgDistance2 = 0.f; // substitute for chi2 + int mROFWeight = 0; // rof weight for voting + short mROF = -1; // rof }; -} // namespace its -} // namespace o2 -#endif /* O2_ITSMFT_TRACKING_LINE_H_ */ +} // namespace o2::its +#endif /* O2_ITS_CLUSTERLINES_H */ diff --git a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/Configuration.h b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/Configuration.h index a700dc1e806c0..8453dce2a11d5 100644 --- a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/Configuration.h +++ b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/Configuration.h @@ -111,7 +111,7 @@ struct VertexingParameters { std::vector LayerRadii = {2.33959f, 3.14076f, 3.91924f, 19.6213f, 24.5597f, 34.388f, 39.3329f}; int ZBins{1}; int PhiBins{128}; - + int deltaRof = 0; float zCut = 0.002f; float phiCut = 0.005f; float pairCut = 0.04f; diff --git a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/TimeFrame.h b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/TimeFrame.h index 3d99192d8f478..f3acb125e28ee 100644 --- a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/TimeFrame.h +++ b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/TimeFrame.h @@ -78,16 +78,16 @@ class TimeFrame friend class TimeFrameGPU; TimeFrame(int nLayers = 7); const Vertex& getPrimaryVertex(const int) const; - gsl::span getPrimaryVertices(int tf) const; + gsl::span getPrimaryVertices(int rofId) const; gsl::span getPrimaryVertices(int romin, int romax) const; - gsl::span> getPrimaryVerticesLabels(const int rof) const; - gsl::span> getPrimaryVerticesXAlpha(int tf) const; + gsl::span> getPrimaryVerticesLabels(const int rofId) const; + gsl::span> getPrimaryVerticesXAlpha(int rofId) const; void fillPrimaryVerticesXandAlpha(); - int getPrimaryVerticesNum(int rofID = -1) const; - void addPrimaryVertices(const std::vector& vertices); - void addPrimaryVertices(const gsl::span& vertices); - void addPrimaryVertices(const std::vector&); - void removePrimaryVerticesInROf(const int rofId); + int getPrimaryVerticesNum(int rofId = -1) const; + void addPrimaryVertices(const std::vector& vertices, const int rofId); + void addPrimaryVertices(const gsl::span& vertices, const int rofId); + // void addPrimaryVertices(const std::vector&); + void removePrimaryVerticesInROF(const int rofId); int loadROFrameData(const o2::itsmft::ROFRecord& rof, gsl::span clusters, const dataformats::MCTruthContainer* mcLabels = nullptr); @@ -100,7 +100,7 @@ class TimeFrame int getTotalClusters() const; bool empty() const; bool isGPU() const { return mIsGPU; } - int getSortedIndex(int rof, int layer, int i) const; + int getSortedIndex(int rofId, int layer, int i) const; int getSortedStartIndex(const int, const int) const; int getNrof() const; @@ -124,8 +124,10 @@ class TimeFrame gsl::span getClustersOnLayer(int rofId, int layerId) const; gsl::span getClustersPerROFrange(int rofMin, int range, int layerId) const; gsl::span getUnsortedClustersOnLayer(int rofId, int layerId) const; - gsl::span getROframesClustersPerROFrange(int rofMin, int range, int layerId) const; - gsl::span getROframeClusters(int layerId) const; + gsl::span getUsedClustersROF(int rofId, int layerId); + gsl::span getUsedClustersROF(int rofId, int layerId) const; + gsl::span getROFramesClustersPerROFrange(int rofMin, int range, int layerId) const; + gsl::span getROFrameClusters(int layerId) const; gsl::span getNClustersROFrange(int rofMin, int range, int layerId) const; gsl::span getIndexTablePerROFrange(int rofMin, int range, int layerId) const; gsl::span getIndexTable(int rofId, int layerId); @@ -145,12 +147,14 @@ class TimeFrame void initialise(const int iteration, const TrackingParameters& trkParam, const int maxLayers = 7); void resetRofPV() { - mPrimaryVertices.clear(); - mROframesPV.resize(1, 0); + // mPrimaryVertices.clear(); + deepVectorClear(mPrimaryVertices); + mROFramesPV.resize(1, 0); }; bool isClusterUsed(int layer, int clusterId) const; void markUsedCluster(int layer, int clusterId); + gsl::span getUsedClusters(const int layer); std::vector>& getTracklets(); std::vector>& getTrackletsLookupTable(); @@ -164,9 +168,9 @@ class TimeFrame std::vector>& getCellsNeighbours(); std::vector>& getCellsNeighboursLUT(); std::vector>& getRoads(); - std::vector& getTracks(int rof) { return mTracks[rof]; } - std::vector& getTracksLabel(const int rof) { return mTracksLabel[rof]; } - std::vector& getLinesLabel(const int rof) { return mLinesLabels[rof]; } + std::vector& getTracks(int rofId) { return mTracks[rofId]; } + std::vector& getTracksLabel(const int rofId) { return mTracksLabel[rofId]; } + std::vector& getLinesLabel(const int rofId) { return mLinesLabels[rofId]; } std::vector>& getVerticesLabels() { return mVerticesLabels; } int getNumberOfClusters() const; @@ -178,22 +182,29 @@ class TimeFrame bool checkMemory(unsigned long max) { return getArtefactsMemory() < max; } unsigned long getArtefactsMemory(); - int getROfCutClusterMult() const { return mCutClusterMult; }; - int getROfCutVertexMult() const { return mCutVertexMult; }; - int getROfCutAllMult() const { return mCutClusterMult + mCutVertexMult; } + int getROFCutClusterMult() const { return mCutClusterMult; }; + int getROFCutVertexMult() const { return mCutVertexMult; }; + int getROFCutAllMult() const { return mCutClusterMult + mCutVertexMult; } // Vertexer - void computeTrackletsScans(const int nThreads = 1); - int& getNTrackletsROf(int tf, int combId); - std::vector& getLines(int tf); - std::vector& getTrackletClusters(int tf); + void computeTrackletsPerROFScans(); + void computeTracletsPerClusterScans(); + int& getNTrackletsROF(int rofId, int combId); + std::vector& getLines(int rofId); + int getNLinesTotal() const + { + return std::accumulate(mLines.begin(), mLines.end(), 0, [](int sum, const auto& l) { return sum + l.size(); }); + } + std::vector& getTrackletClusters(int rofId); gsl::span getFoundTracklets(int rofId, int combId) const; gsl::span getFoundTracklets(int rofId, int combId); gsl::span getLabelsFoundTracklets(int rofId, int combId) const; gsl::span getNTrackletsCluster(int rofId, int combId); + gsl::span getExclusiveNTrackletsCluster(int rofId, int combId); uint32_t getTotalTrackletsTF(const int iLayer) { return mTotalTracklets[iLayer]; } int getTotalClustersPerROFrange(int rofMin, int range, int layerId) const; std::array& getBeamXY() { return mBeamPos; } + void insertLateVertex(const Vertex& vertex); // \Vertexer void initialiseRoadLabels(); @@ -219,7 +230,8 @@ class TimeFrame } } - virtual void setDevicePropagator(const o2::base::PropagatorImpl*){}; + virtual void setDevicePropagator(const o2::base::PropagatorImpl*) {}; + const o2::base::PropagatorImpl* getDevicePropagator() const { return mPropagatorDevice; } template @@ -249,15 +261,16 @@ class TimeFrame std::vector> mClusters; std::vector> mTrackingFrameInfo; std::vector> mClusterExternalIndices; - std::vector> mROframesClusters; + std::vector> mROFramesClusters; const dataformats::MCTruthContainer* mClusterLabels = nullptr; - std::array, 2> mNTrackletsPerCluster; // TODO: remove in favour of mNTrackletsPerROf + std::array, 2> mNTrackletsPerCluster; + std::array, 2> mNTrackletsPerClusterSum; std::vector> mNClustersPerROF; std::vector> mIndexTables; std::vector> mTrackletsLookupTable; std::vector> mUsedClusters; int mNrof = 0; - std::vector mROframesPV = {0}; + std::vector mROFramesPV = {0}; std::vector mPrimaryVertices; // State if memory will be externally managed. @@ -305,10 +318,10 @@ class TimeFrame int mCutVertexMult; // Vertexer - std::vector> mNTrackletsPerROf; + std::vector> mNTrackletsPerROF; std::vector> mLines; std::vector> mTrackletClusters; - std::vector> mTrackletsIndexROf; + std::vector> mTrackletsIndexROF; std::vector> mLinesLabels; std::vector> mVerticesLabels; std::array mTotalTracklets = {0, 0}; @@ -317,45 +330,45 @@ class TimeFrame inline const Vertex& TimeFrame::getPrimaryVertex(const int vertexIndex) const { return mPrimaryVertices[vertexIndex]; } -inline gsl::span TimeFrame::getPrimaryVertices(int rof) const +inline gsl::span TimeFrame::getPrimaryVertices(int rofId) const { - const int start = mROframesPV[rof]; - const int stop_idx = rof >= mNrof - 1 ? mNrof : rof + 1; - int delta = mMultiplicityCutMask[rof] ? mROframesPV[stop_idx] - start : 0; // return empty span if Rof is excluded + const int start = mROFramesPV[rofId]; + const int stop_idx = rofId >= mNrof - 1 ? mNrof : rofId + 1; + int delta = mMultiplicityCutMask[rofId] ? mROFramesPV[stop_idx] - start : 0; // return empty span if Rof is excluded return {&mPrimaryVertices[start], static_cast::size_type>(delta)}; } -inline gsl::span> TimeFrame::getPrimaryVerticesLabels(const int rof) const +inline gsl::span> TimeFrame::getPrimaryVerticesLabels(const int rofId) const { - const int start = mROframesPV[rof]; - const int stop_idx = rof >= mNrof - 1 ? mNrof : rof + 1; - int delta = mMultiplicityCutMask[rof] ? mROframesPV[stop_idx] - start : 0; // return empty span if Rof is excluded + const int start = mROFramesPV[rofId]; + const int stop_idx = rofId >= mNrof - 1 ? mNrof : rofId + 1; + int delta = mMultiplicityCutMask[rofId] ? mROFramesPV[stop_idx] - start : 0; // return empty span if Rof is excluded return {&mVerticesLabels[start], static_cast::size_type>(delta)}; } inline gsl::span TimeFrame::getPrimaryVertices(int romin, int romax) const { - return {&mPrimaryVertices[mROframesPV[romin]], static_cast::size_type>(mROframesPV[romax + 1] - mROframesPV[romin])}; + return {&mPrimaryVertices[mROFramesPV[romin]], static_cast::size_type>(mROFramesPV[romax + 1] - mROFramesPV[romin])}; } -inline gsl::span> TimeFrame::getPrimaryVerticesXAlpha(int rof) const +inline gsl::span> TimeFrame::getPrimaryVerticesXAlpha(int rofId) const { - const int start = mROframesPV[rof]; - const int stop_idx = rof >= mNrof - 1 ? mNrof : rof + 1; - int delta = mMultiplicityCutMask[rof] ? mROframesPV[stop_idx] - start : 0; // return empty span if Rof is excluded + const int start = mROFramesPV[rofId]; + const int stop_idx = rofId >= mNrof - 1 ? mNrof : rofId + 1; + int delta = mMultiplicityCutMask[rofId] ? mROFramesPV[stop_idx] - start : 0; // return empty span if Rof is excluded return {&(mPValphaX[start]), static_cast>::size_type>(delta)}; } -inline int TimeFrame::getPrimaryVerticesNum(int rofID) const +inline int TimeFrame::getPrimaryVerticesNum(int rofId) const { - return rofID < 0 ? mPrimaryVertices.size() : mROframesPV[rofID + 1] - mROframesPV[rofID]; + return rofId < 0 ? mPrimaryVertices.size() : mROFramesPV[rofId + 1] - mROFramesPV[rofId]; } inline bool TimeFrame::empty() const { return getTotalClusters() == 0; } -inline int TimeFrame::getSortedIndex(int rof, int layer, int index) const { return mROframesClusters[layer][rof] + index; } +inline int TimeFrame::getSortedIndex(int rofId, int layer, int index) const { return mROFramesClusters[layer][rofId] + index; } -inline int TimeFrame::getSortedStartIndex(const int rof, const int layer) const { return mROframesClusters[layer][rof]; } +inline int TimeFrame::getSortedStartIndex(const int rofId, const int layer) const { return mROFramesClusters[layer][rofId]; } inline int TimeFrame::getNrof() const { return mNrof; } @@ -370,9 +383,9 @@ inline float TimeFrame::getBeamX() const { return mBeamPos[0]; } inline float TimeFrame::getBeamY() const { return mBeamPos[1]; } -inline gsl::span TimeFrame::getROframeClusters(int layerId) const +inline gsl::span TimeFrame::getROFrameClusters(int layerId) const { - return {&mROframesClusters[layerId][0], static_cast::size_type>(mROframesClusters[layerId].size())}; + return {&mROFramesClusters[layerId][0], static_cast::size_type>(mROFramesClusters[layerId].size())}; } inline gsl::span TimeFrame::getClustersOnLayer(int rofId, int layerId) @@ -380,8 +393,8 @@ inline gsl::span TimeFrame::getClustersOnLayer(int rofId, int layerId) if (rofId < 0 || rofId >= mNrof) { return gsl::span(); } - int startIdx{mROframesClusters[layerId][rofId]}; - return {&mClusters[layerId][startIdx], static_cast::size_type>(mROframesClusters[layerId][rofId + 1] - startIdx)}; + int startIdx{mROFramesClusters[layerId][rofId]}; + return {&mClusters[layerId][startIdx], static_cast::size_type>(mROFramesClusters[layerId][rofId + 1] - startIdx)}; } inline gsl::span TimeFrame::getClustersOnLayer(int rofId, int layerId) const @@ -389,8 +402,26 @@ inline gsl::span TimeFrame::getClustersOnLayer(int rofId, int lay if (rofId < 0 || rofId >= mNrof) { return gsl::span(); } - int startIdx{mROframesClusters[layerId][rofId]}; - return {&mClusters[layerId][startIdx], static_cast::size_type>(mROframesClusters[layerId][rofId + 1] - startIdx)}; + int startIdx{mROFramesClusters[layerId][rofId]}; + return {&mClusters[layerId][startIdx], static_cast::size_type>(mROFramesClusters[layerId][rofId + 1] - startIdx)}; +} + +inline gsl::span TimeFrame::getUsedClustersROF(int rofId, int layerId) +{ + if (rofId < 0 || rofId >= mNrof) { + return gsl::span(); + } + int startIdx{mROFramesClusters[layerId][rofId]}; + return {&mUsedClusters[layerId][startIdx], static_cast::size_type>(mROFramesClusters[layerId][rofId + 1] - startIdx)}; +} + +inline gsl::span TimeFrame::getUsedClustersROF(int rofId, int layerId) const +{ + if (rofId < 0 || rofId >= mNrof) { + return gsl::span(); + } + int startIdx{mROFramesClusters[layerId][rofId]}; + return {&mUsedClusters[layerId][startIdx], static_cast::size_type>(mROFramesClusters[layerId][rofId + 1] - startIdx)}; } inline gsl::span TimeFrame::getClustersPerROFrange(int rofMin, int range, int layerId) const @@ -398,15 +429,15 @@ inline gsl::span TimeFrame::getClustersPerROFrange(int rofMin, in if (rofMin < 0 || rofMin >= mNrof) { return gsl::span(); } - int startIdx{mROframesClusters[layerId][rofMin]}; // First cluster of rofMin - int endIdx{mROframesClusters[layerId][std::min(rofMin + range, mNrof)]}; + int startIdx{mROFramesClusters[layerId][rofMin]}; // First cluster of rofMin + int endIdx{mROFramesClusters[layerId][std::min(rofMin + range, mNrof)]}; return {&mClusters[layerId][startIdx], static_cast::size_type>(endIdx - startIdx)}; } -inline gsl::span TimeFrame::getROframesClustersPerROFrange(int rofMin, int range, int layerId) const +inline gsl::span TimeFrame::getROFramesClustersPerROFrange(int rofMin, int range, int layerId) const { int chkdRange{std::min(range, mNrof - rofMin)}; - return {&mROframesClusters[layerId][rofMin], static_cast::size_type>(chkdRange)}; + return {&mROFramesClusters[layerId][rofMin], static_cast::size_type>(chkdRange)}; } inline gsl::span TimeFrame::getNClustersROFrange(int rofMin, int range, int layerId) const @@ -419,7 +450,7 @@ inline int TimeFrame::getTotalClustersPerROFrange(int rofMin, int range, int lay { int startIdx{rofMin}; // First cluster of rofMin int endIdx{std::min(rofMin + range, mNrof)}; - return mROframesClusters[layerId][endIdx] - mROframesClusters[layerId][startIdx]; + return mROFramesClusters[layerId][endIdx] - mROFramesClusters[layerId][startIdx]; } inline gsl::span TimeFrame::getIndexTablePerROFrange(int rofMin, int range, int layerId) const @@ -431,7 +462,7 @@ inline gsl::span TimeFrame::getIndexTablePerROFrange(int rofMin, int inline int TimeFrame::getClusterROF(int iLayer, int iCluster) { - return std::lower_bound(mROframesClusters[iLayer].begin(), mROframesClusters[iLayer].end(), iCluster + 1) - mROframesClusters[iLayer].begin() - 1; + return std::lower_bound(mROFramesClusters[iLayer].begin(), mROFramesClusters[iLayer].end(), iCluster + 1) - mROFramesClusters[iLayer].begin() - 1; } inline gsl::span TimeFrame::getUnsortedClustersOnLayer(int rofId, int layerId) const @@ -439,8 +470,8 @@ inline gsl::span TimeFrame::getUnsortedClustersOnLayer(int rofId, if (rofId < 0 || rofId >= mNrof) { return gsl::span(); } - int startIdx{mROframesClusters[layerId][rofId]}; - return {&mUnsortedClusters[layerId][startIdx], static_cast::size_type>(mROframesClusters[layerId][rofId + 1] - startIdx)}; + int startIdx{mROFramesClusters[layerId][rofId]}; + return {&mUnsortedClusters[layerId][startIdx], static_cast::size_type>(mROFramesClusters[layerId][rofId + 1] - startIdx)}; } inline const std::vector& TimeFrame::getTrackingFrameInfoOnLayer(int layerId) const @@ -482,14 +513,14 @@ inline gsl::span TimeFrame::getIndexTable(int rofId, int layer) static_cast::size_type>(mIndexTableUtils.getNphiBins() * mIndexTableUtils.getNzBins() + 1)}; } -inline std::vector& TimeFrame::getLines(int rof) +inline std::vector& TimeFrame::getLines(int rofId) { - return mLines[rof]; + return mLines[rofId]; } -inline std::vector& TimeFrame::getTrackletClusters(int rof) +inline std::vector& TimeFrame::getTrackletClusters(int rofId) { - return mTrackletClusters[rof]; + return mTrackletClusters[rofId]; } template @@ -519,6 +550,11 @@ inline bool TimeFrame::isClusterUsed(int layer, int clusterId) const return mUsedClusters[layer][clusterId]; } +inline gsl::span TimeFrame::getUsedClusters(const int layer) +{ + return {&mUsedClusters[layer][0], static_cast::size_type>(mUsedClusters[layer].size())}; +} + inline void TimeFrame::markUsedCluster(int layer, int clusterId) { mUsedClusters[layer][clusterId] = true; } inline std::vector>& TimeFrame::getTracklets() @@ -553,13 +589,23 @@ inline gsl::span TimeFrame::getNTrackletsCluster(int rofId, int combId) if (rofId < 0 || rofId >= mNrof) { return gsl::span(); } - auto startIdx{mROframesClusters[1][rofId]}; - return {&mNTrackletsPerCluster[combId][startIdx], static_cast::size_type>(mROframesClusters[1][rofId + 1] - startIdx)}; + auto startIdx{mROFramesClusters[1][rofId]}; + return {&mNTrackletsPerCluster[combId][startIdx], static_cast::size_type>(mROFramesClusters[1][rofId + 1] - startIdx)}; } -inline int& TimeFrame::getNTrackletsROf(int rof, int combId) +inline gsl::span TimeFrame::getExclusiveNTrackletsCluster(int rofId, int combId) { - return mNTrackletsPerROf[combId][rof]; + if (rofId < 0 || rofId >= mNrof) { + return gsl::span(); + } + auto clusStartIdx{mROFramesClusters[1][rofId]}; + + return {&mNTrackletsPerClusterSum[combId][clusStartIdx], static_cast::size_type>(mROFramesClusters[1][rofId + 1] - clusStartIdx)}; +} + +inline int& TimeFrame::getNTrackletsROF(int rofId, int combId) +{ + return mNTrackletsPerROF[combId][rofId]; } inline bool TimeFrame::isRoadFake(int i) const @@ -594,8 +640,8 @@ inline gsl::span TimeFrame::getFoundTracklets(int rofId, int combId) if (rofId < 0 || rofId >= mNrof) { return gsl::span(); } - auto startIdx{mNTrackletsPerROf[combId][rofId]}; - return {&mTracklets[combId][startIdx], static_cast::size_type>(mNTrackletsPerROf[combId][rofId + 1] - startIdx)}; + auto startIdx{mNTrackletsPerROF[combId][rofId]}; + return {&mTracklets[combId][startIdx], static_cast::size_type>(mNTrackletsPerROF[combId][rofId + 1] - startIdx)}; } inline gsl::span TimeFrame::getFoundTracklets(int rofId, int combId) const @@ -603,8 +649,8 @@ inline gsl::span TimeFrame::getFoundTracklets(int rofId, int com if (rofId < 0 || rofId >= mNrof) { return gsl::span(); } - auto startIdx{mNTrackletsPerROf[combId][rofId]}; - return {&mTracklets[combId][startIdx], static_cast::size_type>(mNTrackletsPerROf[combId][rofId + 1] - startIdx)}; + auto startIdx{mNTrackletsPerROF[combId][rofId]}; + return {&mTracklets[combId][startIdx], static_cast::size_type>(mNTrackletsPerROF[combId][rofId + 1] - startIdx)}; } inline gsl::span TimeFrame::getLabelsFoundTracklets(int rofId, int combId) const @@ -612,8 +658,8 @@ inline gsl::span TimeFrame::getLabelsFoundTracklets(int rofId if (rofId < 0 || rofId >= mNrof || !hasMCinformation()) { return gsl::span(); } - auto startIdx{mNTrackletsPerROf[combId][rofId]}; - return {&mTrackletLabels[combId][startIdx], static_cast::size_type>(mNTrackletsPerROf[combId][rofId + 1] - startIdx)}; + auto startIdx{mNTrackletsPerROF[combId][rofId]}; + return {&mTrackletLabels[combId][startIdx], static_cast::size_type>(mNTrackletsPerROF[combId][rofId + 1] - startIdx)}; } inline int TimeFrame::getNumberOfClusters() const @@ -670,6 +716,15 @@ inline size_t TimeFrame::getNumberOfUsedClusters() const return nClusters; } +inline void TimeFrame::insertLateVertex(const Vertex& vertex) +{ + int rofId = vertex.getTimeStamp().getTimeStamp(); + mPrimaryVertices.insert(mPrimaryVertices.begin() + mROFramesPV[rofId + 1], vertex); // inserted at the end of sequence of vertices belonging to rofId + for (int i = rofId + 1; i < mROFramesPV.size(); ++i) { + mROFramesPV[i] += 1; + } +} + } // namespace its } // namespace o2 diff --git a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/TrackingConfigParam.h b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/TrackingConfigParam.h index 36a5fd63b12d1..5ccb74babea6d 100644 --- a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/TrackingConfigParam.h +++ b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/TrackingConfigParam.h @@ -23,6 +23,8 @@ namespace its struct VertexerParamConfig : public o2::conf::ConfigurableParamHelper { bool allowSingleContribClusters = false; + // Number of ROFs to be considered for the vertexing + int deltaRof = 0; // geometrical cuts float zCut = 0.002f; diff --git a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/Tracklet.h b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/Tracklet.h index a7f7adbb8abad..37c02b90b7510 100644 --- a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/Tracklet.h +++ b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/Tracklet.h @@ -21,6 +21,10 @@ #include "GPUCommonMath.h" #include "GPUCommonDef.h" +#ifndef GPUCA_GPUCODE_DEVICE +#include +#endif + namespace o2 { namespace its @@ -36,11 +40,19 @@ struct Tracklet final { { return firstClusterIndex < 0 || secondClusterIndex < 0; } + GPUhdi() auto getDeltaRof() const { return rof[1] - rof[0]; } GPUhdi() void dump(); GPUhdi() void dump() const; GPUhdi() void dump(const int, const int); GPUhdi() void dump(const int, const int) const; GPUhdi() unsigned char operator<(const Tracklet&) const; +#ifndef GPUCA_GPUCODE_DEVICE + std::string asString() const + { + return "fClIdx: " + std::to_string(firstClusterIndex) + " sClIdx: " + std::to_string(secondClusterIndex) + + " rof1: " + std::to_string(rof[0]) + " rof2: " + std::to_string(rof[1]) + " delta: " + std::to_string(getDeltaRof()); + } +#endif int firstClusterIndex; int secondClusterIndex; @@ -104,26 +116,6 @@ GPUhdi() unsigned char Tracklet::operator<(const Tracklet& t) const return true; } -// GPUhdi() void Tracklet::dump() -// { -// printf("fClIdx: %d sClIdx: %d rof1: %hu rof2: %hu phi: %f tl: %f \n", firstClusterIndex, secondClusterIndex, rof[0], rof[1], phi, tanLambda); -// } - -// GPUhdi() void Tracklet::dump() const -// { -// printf("fClIdx: %d sClIdx: %d rof1: %hu rof2: %hu phi: %f tl: %f \n", firstClusterIndex, secondClusterIndex, rof[0], rof[1], phi, tanLambda); -// } - -// GPUhdi() void Tracklet::dump(const int offsetFirst, const int offsetSecond) -// { -// printf("fClIdx: %d sClIdx: %d rof1: %hu rof2: %hu phi: %f tl: %f \n", firstClusterIndex + offsetFirst, secondClusterIndex + offsetSecond, rof[0], rof[1], phi, tanLambda); -// } - -// GPUhdi() void Tracklet::dump(const int offsetFirst, const int offsetSecond) const -// { -// printf("fClIdx: %d sClIdx: %d rof1: %hu rof2: %hu phi: %f tl: %f \n", firstClusterIndex + offsetFirst, secondClusterIndex + offsetSecond, rof[0], rof[1], phi, tanLambda); -// } - GPUhdi() void Tracklet::dump(const int offsetFirst, const int offsetSecond) { printf("fClIdx: %d sClIdx: %d rof1: %hu rof2: %hu\n", firstClusterIndex + offsetFirst, secondClusterIndex + offsetSecond, rof[0], rof[1]); diff --git a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/Vertexer.h b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/Vertexer.h index eb46c60305bda..544eb872042f7 100644 --- a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/Vertexer.h +++ b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/Vertexer.h @@ -88,7 +88,10 @@ class Vertexer void dumpTraits(); template float evaluateTask(void (Vertexer::*)(T...), const char*, std::function logger, T&&... args); - void printEpilog(std::function logger, const float total); + void printEpilog(std::function logger, + bool isHybrid, + const unsigned int trackletN01, const unsigned int trackletN12, const unsigned selectedN, + const float initT, const float trackletT, const float selecT, const float vertexT); private: std::uint32_t mTimeFrameCounter = 0; diff --git a/Detectors/ITSMFT/ITS/tracking/src/ClusterLines.cxx b/Detectors/ITSMFT/ITS/tracking/src/ClusterLines.cxx index 8605bb0ecada2..32b78d4f2ef18 100644 --- a/Detectors/ITSMFT/ITS/tracking/src/ClusterLines.cxx +++ b/Detectors/ITSMFT/ITS/tracking/src/ClusterLines.cxx @@ -84,6 +84,9 @@ ClusterLines::ClusterLines(const int firstLabel, const Line& firstLine, const in const bool weight) { + updateROFPoll(firstLine); + updateROFPoll(secondLine); + mLabels.push_back(firstLabel); if (secondLabel > 0) { mLabels.push_back(secondLabel); // don't add info in case of beamline used @@ -188,7 +191,8 @@ ClusterLines::ClusterLines(const Line& firstLine, const Line& secondLine) std::array covarianceFirst{1., 1., 1.}; std::array covarianceSecond{1., 1., 1.}; - + updateROFPoll(firstLine); + updateROFPoll(secondLine); for (int i{0}; i < 6; ++i) { mWeightMatrix[i] = firstLine.weightMatrix[i] + secondLine.weightMatrix[i]; } @@ -274,6 +278,7 @@ ClusterLines::ClusterLines(const Line& firstLine, const Line& secondLine) void ClusterLines::add(const int& lineLabel, const Line& line, const bool& weight) { mLabels.push_back(lineLabel); + updateROFPoll(line); std::array covariance{1., 1., 1.}; for (int i{0}; i < 6; ++i) { @@ -361,5 +366,21 @@ bool ClusterLines::operator==(const ClusterLines& rhs) const } return retval && this->mAvgDistance2 == rhs.mAvgDistance2; } + +GPUhdi() void ClusterLines::updateROFPoll(const Line& line) +{ + // Boyer-Moore voting for rof label + if (mROFWeight == 0) { + mROF = line.getMinROF(); + mROFWeight = 1; + } else { + if (mROF == line.getMinROF()) { + mROFWeight++; + } else { + mROFWeight--; + } + } +} + } // namespace its } // namespace o2 diff --git a/Detectors/ITSMFT/ITS/tracking/src/TimeFrame.cxx b/Detectors/ITSMFT/ITS/tracking/src/TimeFrame.cxx index 6e7114d9ca54d..b75a475df3577 100644 --- a/Detectors/ITSMFT/ITS/tracking/src/TimeFrame.cxx +++ b/Detectors/ITSMFT/ITS/tracking/src/TimeFrame.cxx @@ -73,15 +73,19 @@ TimeFrame::TimeFrame(int nLayers) mTrackingFrameInfo.resize(nLayers); mClusterExternalIndices.resize(nLayers); mUsedClusters.resize(nLayers); - mROframesClusters.resize(nLayers, {0}); /// TBC: if resetting the timeframe is required, then this has to be done + mROFramesClusters.resize(nLayers, {0}); /// TBC: if resetting the timeframe is required, then this has to be done mNClustersPerROF.resize(nLayers); - mTrackletsIndexROf.resize(2, {0}); + mTrackletsIndexROF.resize(2, {0}); } -void TimeFrame::addPrimaryVertices(const std::vector& vertices) +void TimeFrame::addPrimaryVertices(const std::vector& vertices, const int rofId) { for (const auto& vertex : vertices) { - mPrimaryVertices.emplace_back(vertex); + if (vertex.getTimeStamp().getTimeStamp() != rofId) { + insertLateVertex(vertex); + } else { + mPrimaryVertices.emplace_back(vertex); + } if (!isBeamPositionOverridden) { const int w{vertex.getNContributors()}; mBeamPos[0] = (mBeamPos[0] * mBeamPosWeight + vertex.getX() * w) / (mBeamPosWeight + w); @@ -89,11 +93,12 @@ void TimeFrame::addPrimaryVertices(const std::vector& vertices) mBeamPosWeight += w; } } - mROframesPV.push_back(mPrimaryVertices.size()); + mROFramesPV.push_back(mPrimaryVertices.size()); } -void TimeFrame::addPrimaryVertices(const gsl::span& vertices) +void TimeFrame::addPrimaryVertices(const gsl::span& vertices, const int rofId) { + // auto rofId = mROFramesPV.size(); for (const auto& vertex : vertices) { mPrimaryVertices.emplace_back(vertex); if (!isBeamPositionOverridden) { @@ -103,18 +108,18 @@ void TimeFrame::addPrimaryVertices(const gsl::span& vertices) mBeamPosWeight += w; } } - mROframesPV.push_back(mPrimaryVertices.size()); + mROFramesPV.push_back(mPrimaryVertices.size()); } -void TimeFrame::addPrimaryVertices(const std::vector& lVertices) -{ - std::vector vertices; - for (auto& vertex : lVertices) { - vertices.emplace_back(o2::math_utils::Point3D(vertex.mX, vertex.mY, vertex.mZ), vertex.mRMS2, vertex.mContributors, vertex.mAvgDistance2); - vertices.back().setTimeStamp(vertex.mTimeStamp); - } - addPrimaryVertices(vertices); -} +// void TimeFrame::addPrimaryVertices(const std::vector& lVertices) +// { +// std::vector vertices; +// for (auto& vertex : lVertices) { +// vertices.emplace_back(o2::math_utils::Point3D(vertex.mX, vertex.mY, vertex.mZ), vertex.mRMS2, vertex.mContributors, vertex.mAvgDistance2); +// vertices.back().setTimeStamp(vertex.mTimeStamp); +// } +// addPrimaryVertices(vertices); +// } int TimeFrame::loadROFrameData(const o2::itsmft::ROFRecord& rof, gsl::span clusters, const dataformats::MCTruthContainer* mcLabels) @@ -141,10 +146,10 @@ int TimeFrame::loadROFrameData(const o2::itsmft::ROFRecord& rof, gsl::span rofs, deepVectorClear(mUnsortedClusters[iLayer]); deepVectorClear(mTrackingFrameInfo[iLayer]); deepVectorClear(mClusterExternalIndices[iLayer]); - mROframesClusters[iLayer].resize(1, 0); + mROFramesClusters[iLayer].resize(1, 0); if (iLayer < 2) { - deepVectorClear(mTrackletsIndexROf[iLayer]); + deepVectorClear(mTrackletsIndexROF[iLayer]); + deepVectorClear(mNTrackletsPerCluster[iLayer]); + deepVectorClear(mNTrackletsPerClusterSum[iLayer]); } } @@ -223,14 +230,15 @@ int TimeFrame::loadROFrameData(gsl::span rofs, addClusterExternalIndexToLayer(layer, clusterId); } for (unsigned int iL{0}; iL < mUnsortedClusters.size(); ++iL) { - // mNClustersPerROF[iL].push_back(mUnsortedClusters[iL].size() - mROframesClusters[iL].back()); - mROframesClusters[iL].push_back(mUnsortedClusters[iL].size()); + // mNClustersPerROF[iL].push_back(mUnsortedClusters[iL].size() - mROFramesClusters[iL].back()); + mROFramesClusters[iL].push_back(mUnsortedClusters[iL].size()); } mNrof++; } - for (auto& v : mNTrackletsPerCluster) { - v.resize(mUnsortedClusters[1].size()); + for (auto i = 0; i < mNTrackletsPerCluster.size(); ++i) { + mNTrackletsPerCluster[i].resize(mUnsortedClusters[1].size()); + mNTrackletsPerClusterSum[i].resize(mUnsortedClusters[1].size() + 1); // Exc sum "prepends" a 0 } if (mcLabels) { @@ -285,8 +293,8 @@ void TimeFrame::initialise(const int iteration, const TrackingParameters& trkPar mIndexTables.resize(mClusters.size(), std::vector(mNrof * (trkParam.ZBins * trkParam.PhiBins + 1), 0)); mLines.resize(mNrof); mTrackletClusters.resize(mNrof); - mNTrackletsPerROf.resize(2); - for (auto& v : mNTrackletsPerROf) { + mNTrackletsPerROF.resize(2); + for (auto& v : mNTrackletsPerROF) { v = std::vector(mNrof + 1, 0); } @@ -433,6 +441,15 @@ void TimeFrame::fillPrimaryVerticesXandAlpha() } } +void TimeFrame::computeTrackletsPerROFScans() +{ + for (ushort iLayer = 0; iLayer < 2; ++iLayer) { + mTotalTracklets[iLayer] = std::accumulate(mNTrackletsPerROF[iLayer].begin(), mNTrackletsPerROF[iLayer].end(), 0); + std::exclusive_scan(mNTrackletsPerROF[iLayer].begin(), mNTrackletsPerROF[iLayer].end(), mNTrackletsPerROF[iLayer].begin(), 0); + std::exclusive_scan(mNTrackletsPerCluster[iLayer].begin(), mNTrackletsPerCluster[iLayer].end(), mNTrackletsPerClusterSum[iLayer].begin(), 0); + } +} + void TimeFrame::checkTrackletLUTs() { for (uint32_t iLayer{0}; iLayer < getTracklets().size(); ++iLayer) { @@ -474,9 +491,9 @@ void TimeFrame::resizeVectors(int nLayers) mTrackingFrameInfo.resize(nLayers); mClusterExternalIndices.resize(nLayers); mUsedClusters.resize(nLayers); - mROframesClusters.resize(nLayers, {0}); + mROFramesClusters.resize(nLayers, {0}); mNClustersPerROF.resize(nLayers); - mTrackletsIndexROf.resize(2, {0}); + mTrackletsIndexROF.resize(2, {0}); } void TimeFrame::printTrackletLUTonLayer(int i) @@ -517,9 +534,9 @@ void TimeFrame::printCellLUTs() void TimeFrame::printVertices() { - std::cout << "Vertices in ROF (nROF = " << mNrof << ", lut size = " << mROframesPV.size() << ")" << std::endl; - for (unsigned int iR{0}; iR < mROframesPV.size(); ++iR) { - std::cout << mROframesPV[iR] << "\t"; + std::cout << "Vertices in ROF (nROF = " << mNrof << ", lut size = " << mROFramesPV.size() << ")" << std::endl; + for (unsigned int iR{0}; iR < mROFramesPV.size(); ++iR) { + std::cout << mROFramesPV[iR] << "\t"; } std::cout << "\n\n Vertices:" << std::endl; for (unsigned int iV{0}; iV < mPrimaryVertices.size(); ++iV) { @@ -531,9 +548,9 @@ void TimeFrame::printVertices() void TimeFrame::printROFoffsets() { std::cout << "--------" << std::endl; - for (unsigned int iLayer{0}; iLayer < mROframesClusters.size(); ++iLayer) { + for (unsigned int iLayer{0}; iLayer < mROFramesClusters.size(); ++iLayer) { std::cout << "Layer " << iLayer << std::endl; - for (auto value : mROframesClusters[iLayer]) { + for (auto value : mROFramesClusters[iLayer]) { std::cout << value << "\t"; } std::cout << std::endl; @@ -551,15 +568,5 @@ void TimeFrame::printNClsPerROF() std::cout << std::endl; } } - -void TimeFrame::computeTrackletsScans(const int nThreads) -{ - // #pragma omp parallel for num_threads(nThreads > 1 ? 2 : nThreads) - for (ushort iLayer = 0; iLayer < 2; ++iLayer) { - mTotalTracklets[iLayer] = std::accumulate(mNTrackletsPerROf[iLayer].begin(), mNTrackletsPerROf[iLayer].end(), 0); - std::exclusive_scan(mNTrackletsPerROf[iLayer].begin(), mNTrackletsPerROf[iLayer].end(), mNTrackletsPerROf[iLayer].begin(), 0); - } -} - } // namespace its } // namespace o2 diff --git a/Detectors/ITSMFT/ITS/tracking/src/Tracker.cxx b/Detectors/ITSMFT/ITS/tracking/src/Tracker.cxx index 42df15b24f052..c9d5dc9451f22 100644 --- a/Detectors/ITSMFT/ITS/tracking/src/Tracker.cxx +++ b/Detectors/ITSMFT/ITS/tracking/src/Tracker.cxx @@ -102,10 +102,10 @@ void Tracker::clustersToTracks(std::function logger, std::f } iVertex++; } while (iVertex < maxNvertices); - logger(fmt::format("- Tracklet finding: {} tracklets found in {:.2f} ms", nTracklets, timeTracklets)); - logger(fmt::format("- Cell finding: {} cells found in {:.2f} ms", nCells, timeCells)); - logger(fmt::format("- Neighbours finding: {} neighbours found in {:.2f} ms", nNeighbours, timeNeighbours)); - logger(fmt::format("- Track finding: {} tracks found in {:.2f} ms", nTracks + mTimeFrame->getNumberOfTracks(), timeRoads)); + logger(fmt::format(" - Tracklet finding: {} tracklets found in {:.2f} ms", nTracklets, timeTracklets)); + logger(fmt::format(" - Cell finding: {} cells found in {:.2f} ms", nCells, timeCells)); + logger(fmt::format(" - Neighbours finding: {} neighbours found in {:.2f} ms", nNeighbours, timeNeighbours)); + logger(fmt::format(" - Track finding: {} tracks found in {:.2f} ms", nTracks + mTimeFrame->getNumberOfTracks(), timeRoads)); total += timeTracklets + timeCells + timeNeighbours + timeRoads; total += evaluateTask(&Tracker::extendTracks, "Extending tracks", logger, iteration); } diff --git a/Detectors/ITSMFT/ITS/tracking/src/TrackingInterface.cxx b/Detectors/ITSMFT/ITS/tracking/src/TrackingInterface.cxx index 492976d164227..73056595c644c 100644 --- a/Detectors/ITSMFT/ITS/tracking/src/TrackingInterface.cxx +++ b/Detectors/ITSMFT/ITS/tracking/src/TrackingInterface.cxx @@ -208,7 +208,7 @@ void ITSTrackingInterface::run(framework::ProcessingContext& pc) for (auto& v : vtxVecLoc) { vertices.push_back(v); } - mTimeFrame->addPrimaryVertices(vtxVecLoc); + mTimeFrame->addPrimaryVertices(vtxVecLoc, iRof); } } LOG(info) << fmt::format(" - rejected {}/{} ROFs: random/mult.sel:{} (seed {}), vtx.sel:{}", cutRandomMult + cutVertexMult, rofspan.size(), cutRandomMult, multEst.lastRandomSeed, cutVertexMult); diff --git a/Detectors/ITSMFT/ITS/tracking/src/Vertexer.cxx b/Detectors/ITSMFT/ITS/tracking/src/Vertexer.cxx index acc59a66a7c2c..16777c66c0425 100644 --- a/Detectors/ITSMFT/ITS/tracking/src/Vertexer.cxx +++ b/Detectors/ITSMFT/ITS/tracking/src/Vertexer.cxx @@ -38,30 +38,39 @@ Vertexer::Vertexer(VertexerTraits* traits) float Vertexer::clustersToVertices(std::function logger) { - float total{0.f}; + float timeTracklet, timeSelection, timeVertexing, timeInit; + unsigned int nTracklets01, nTracklets12, nSelectedTracklets; TrackingParameters trkPars; trkPars.PhiBins = mTraits->getVertexingParameters().PhiBins; trkPars.ZBins = mTraits->getVertexingParameters().ZBins; - total += evaluateTask(&Vertexer::initialiseVertexer, "Vertexer initialisation", logger, trkPars); - total += evaluateTask(&Vertexer::findTracklets, "Vertexer tracklet finding", logger); - total += evaluateTask(&Vertexer::validateTracklets, "Vertexer adjacent tracklets validation", logger); - total += evaluateTask(&Vertexer::findVertices, "Vertexer vertex finding", logger); - printEpilog(logger, total); - return total; + timeInit += evaluateTask(&Vertexer::initialiseVertexer, "Vertexer initialisation", [](std::string) {}, trkPars); + timeTracklet = evaluateTask(&Vertexer::findTracklets, "Vertexer tracklet finding", [](std::string) {}); + nTracklets01 = mTimeFrame->getTotalTrackletsTF(0); + nTracklets12 = mTimeFrame->getTotalTrackletsTF(1); + timeSelection = evaluateTask(&Vertexer::validateTracklets, "Vertexer adjacent tracklets validation", [](std::string) {}); + nSelectedTracklets = mTimeFrame->getNLinesTotal(); + timeVertexing = evaluateTask(&Vertexer::findVertices, "Vertexer vertex finding", [](std::string) {}); + printEpilog(logger, false, nTracklets01, nTracklets12, nSelectedTracklets, timeInit, timeTracklet, timeSelection, timeVertexing); + return timeInit + timeTracklet + timeSelection + timeVertexing; } float Vertexer::clustersToVerticesHybrid(std::function logger) { - float total{0.f}; + float timeTracklet, timeSelection, timeVertexing, timeInit; + unsigned int nTracklets01, nTracklets12, nSelectedTracklets; TrackingParameters trkPars; trkPars.PhiBins = mTraits->getVertexingParameters().PhiBins; trkPars.ZBins = mTraits->getVertexingParameters().ZBins; - total += evaluateTask(&Vertexer::initialiseVertexerHybrid, "Hybrid Vertexer initialisation", logger, trkPars); - total += evaluateTask(&Vertexer::findTrackletsHybrid, "Hybrid Vertexer tracklet finding", logger); - total += evaluateTask(&Vertexer::validateTrackletsHybrid, "Hybrid Vertexer adjacent tracklets validation", logger); - total += evaluateTask(&Vertexer::findVerticesHybrid, "Hybrid Vertexer vertex finding", logger); - printEpilog(logger, total); - return total; + timeInit += evaluateTask( + &Vertexer::initialiseVertexerHybrid, "Hybrid Vertexer initialisation", [](std::string) {}, trkPars); + timeTracklet = evaluateTask(&Vertexer::findTrackletsHybrid, "Hybrid Vertexer tracklet finding", [](std::string) {}); + nTracklets01 = mTimeFrame->getTotalTrackletsTF(0); + nTracklets12 = mTimeFrame->getTotalTrackletsTF(1); + timeSelection = evaluateTask(&Vertexer::validateTrackletsHybrid, "Hybrid Vertexer adjacent tracklets validation", [](std::string) {}); + nSelectedTracklets = mTimeFrame->getNLinesTotal(); + timeVertexing = evaluateTask(&Vertexer::findVerticesHybrid, "Hybrid Vertexer vertex finding", [](std::string) {}); + printEpilog(logger, false, nTracklets01, nTracklets12, nSelectedTracklets, timeInit, timeTracklet, timeSelection, timeVertexing); + return timeInit + timeTracklet + timeSelection + timeVertexing; } void Vertexer::getGlobalConfiguration() @@ -71,6 +80,7 @@ void Vertexer::getGlobalConfiguration() auto& grc = o2::its::GpuRecoParamConfig::Instance(); VertexingParameters verPar; + verPar.deltaRof = vc.deltaRof; verPar.allowSingleContribClusters = vc.allowSingleContribClusters; verPar.zCut = vc.zCut; verPar.phiCut = vc.phiCut; @@ -102,8 +112,15 @@ void Vertexer::adoptTimeFrame(TimeFrame& tf) mTraits->adoptTimeFrame(&tf); } -void Vertexer::printEpilog(std::function logger, const float total) +void Vertexer::printEpilog(std::function logger, + bool isHybrid, + const unsigned int trackletN01, const unsigned int trackletN12, const unsigned selectedN, + const float initT, const float trackletT, const float selecT, const float vertexT) { + float total = initT + trackletT + selecT + vertexT; + logger(fmt::format(" - {}Vertexer: trackleting found {}|{} tracklets, completed in: {} ms", isHybrid ? "Hybrid" : "", trackletN01, trackletN12, trackletT)); + logger(fmt::format(" - {}Vertexer: selected {} tracklets, completed in: {} ms", isHybrid ? "Hybrid" : "", selectedN, selecT)); + logger(fmt::format(" - {}Vertexer: vertexing completed in: {} ms", isHybrid ? "Hybrid" : "", vertexT)); logger(fmt::format(" - Timeframe {} vertexing completed in: {} ms, using {} thread(s).", mTimeFrameCounter++, total, mTraits->getNThreads())); } diff --git a/Detectors/ITSMFT/ITS/tracking/src/VertexerTraits.cxx b/Detectors/ITSMFT/ITS/tracking/src/VertexerTraits.cxx index ceb5e4d0d44b7..1e8cea8f722b1 100644 --- a/Detectors/ITSMFT/ITS/tracking/src/VertexerTraits.cxx +++ b/Detectors/ITSMFT/ITS/tracking/src/VertexerTraits.cxx @@ -44,23 +44,24 @@ float smallestAngleDifference(float a, float b) return (diff < -constants::math::Pi) ? diff + constants::math::TwoPi : ((diff > constants::math::Pi) ? diff - constants::math::TwoPi : diff); } -template +template void trackleterKernelHost( - const gsl::span& clustersNextLayer, // 0 2 - const gsl::span& clustersCurrentLayer, // 1 1 + const gsl::span& clustersNextLayer, // 0 2 + const gsl::span& clustersCurrentLayer, // 1 1 + const gsl::span& usedClustersNextLayer, // 0 2 int* indexTableNext, const float phiCut, std::vector& tracklets, gsl::span foundTracklets, const IndexTableUtils& utils, - const int rof, - const int rofFoundTrackletsOffset, + const int pivotRof, + const int targetRof, + gsl::span rofFoundTrackletsOffsets, // we want to change those, to keep track of the offset in deltaRof>0 const int maxTrackletsPerCluster = static_cast(2e3)) { const int PhiBins{utils.getNphiBins()}; const int ZBins{utils.getNzBins()}; // loop on layer1 clusters - int cumulativeStoredTracklets{0}; for (int iCurrentLayerClusterIndex = 0; iCurrentLayerClusterIndex < clustersCurrentLayer.size(); ++iCurrentLayerClusterIndex) { int storedTracklets{0}; const Cluster& currentCluster{clustersCurrentLayer[iCurrentLayerClusterIndex]}; @@ -77,14 +78,18 @@ void trackleterKernelHost( const int maxRowClusterIndex{indexTableNext[firstBinIndex + ZBins]}; // loop on clusters next layer for (int iNextLayerClusterIndex{firstRowClusterIndex}; iNextLayerClusterIndex < maxRowClusterIndex && iNextLayerClusterIndex < static_cast(clustersNextLayer.size()); ++iNextLayerClusterIndex) { + if (usedClustersNextLayer[iNextLayerClusterIndex]) { + LOGP(warning, "skipping used cluster"); + continue; + } const Cluster& nextCluster{clustersNextLayer[iNextLayerClusterIndex]}; if (o2::gpu::GPUCommonMath::Abs(smallestAngleDifference(currentCluster.phi, nextCluster.phi)) < phiCut) { if (storedTracklets < maxTrackletsPerCluster) { - if constexpr (!DryRun) { + if constexpr (!EvalRun) { if constexpr (Mode == TrackletMode::Layer0Layer1) { - tracklets[rofFoundTrackletsOffset + cumulativeStoredTracklets + storedTracklets] = Tracklet{iNextLayerClusterIndex, iCurrentLayerClusterIndex, nextCluster, currentCluster, rof, rof}; + tracklets[rofFoundTrackletsOffsets[iCurrentLayerClusterIndex] + storedTracklets] = Tracklet{iNextLayerClusterIndex, iCurrentLayerClusterIndex, nextCluster, currentCluster, targetRof, pivotRof}; } else { - tracklets[rofFoundTrackletsOffset + cumulativeStoredTracklets + storedTracklets] = Tracklet{iCurrentLayerClusterIndex, iNextLayerClusterIndex, currentCluster, nextCluster, rof, rof}; + tracklets[rofFoundTrackletsOffsets[iCurrentLayerClusterIndex] + storedTracklets] = Tracklet{iCurrentLayerClusterIndex, iNextLayerClusterIndex, currentCluster, nextCluster, pivotRof, targetRof}; } } ++storedTracklets; @@ -93,10 +98,10 @@ void trackleterKernelHost( } } } - if constexpr (DryRun) { - foundTracklets[iCurrentLayerClusterIndex] = storedTracklets; + if constexpr (EvalRun) { + foundTracklets[iCurrentLayerClusterIndex] += storedTracklets; } else { - cumulativeStoredTracklets += storedTracklets; + rofFoundTrackletsOffsets[iCurrentLayerClusterIndex] += storedTracklets; } } } @@ -104,28 +109,39 @@ void trackleterKernelHost( void trackletSelectionKernelHost( const gsl::span clusters0, // 0 const gsl::span clusters1, // 1 + gsl::span usedClusters0, // Layer 0 + gsl::span usedClusters2, // Layer 2 const gsl::span& tracklets01, const gsl::span& tracklets12, + std::vector& usedTracklets, const gsl::span foundTracklets01, const gsl::span foundTracklets12, - std::vector& destTracklets, + std::vector& lines, const gsl::span& trackletLabels, std::vector& linesLabels, + const int pivotRofId, + const int targetRofId, const float tanLambdaCut = 0.025f, const float phiCut = 0.005f, const int maxTracklets = static_cast(1e2)) { int offset01{0}, offset12{0}; - std::vector usedTracklets(tracklets01.size(), false); for (unsigned int iCurrentLayerClusterIndex{0}; iCurrentLayerClusterIndex < clusters1.size(); ++iCurrentLayerClusterIndex) { int validTracklets{0}; for (int iTracklet12{offset12}; iTracklet12 < offset12 + foundTracklets12[iCurrentLayerClusterIndex]; ++iTracklet12) { for (int iTracklet01{offset01}; iTracklet01 < offset01 + foundTracklets01[iCurrentLayerClusterIndex]; ++iTracklet01) { - const float deltaTanLambda{o2::gpu::GPUCommonMath::Abs(tracklets01[iTracklet01].tanLambda - tracklets12[iTracklet12].tanLambda)}; - const float deltaPhi{o2::gpu::GPUCommonMath::Abs(smallestAngleDifference(tracklets01[iTracklet01].phi, tracklets12[iTracklet12].phi))}; + const auto& tracklet01{tracklets01[iTracklet01]}; + const auto& tracklet12{tracklets12[iTracklet12]}; + if (tracklet01.rof[0] != targetRofId || tracklet12.rof[1] != targetRofId) { + continue; + } + const float deltaTanLambda{o2::gpu::GPUCommonMath::Abs(tracklet01.tanLambda - tracklets12[iTracklet12].tanLambda)}; + const float deltaPhi{o2::gpu::GPUCommonMath::Abs(smallestAngleDifference(tracklet01.phi, tracklets12[iTracklet12].phi))}; if (!usedTracklets[iTracklet01] && deltaTanLambda < tanLambdaCut && deltaPhi < phiCut && validTracklets != maxTracklets) { + usedClusters0[tracklet01.firstClusterIndex] = true; + usedClusters2[tracklet12.secondClusterIndex] = true; usedTracklets[iTracklet01] = true; - destTracklets.emplace_back(tracklets01[iTracklet01], clusters0.data(), clusters1.data()); + lines.emplace_back(tracklet01, clusters0.data(), clusters1.data()); if (trackletLabels.size()) { linesLabels.emplace_back(trackletLabels[iTracklet01]); } @@ -173,63 +189,79 @@ void VertexerTraits::computeTracklets() #pragma omp parallel num_threads(mNThreads) { #pragma omp for schedule(dynamic) - for (int rofId = 0; rofId < mTimeFrame->getNrof(); ++rofId) { - trackleterKernelHost( - mTimeFrame->getClustersOnLayer(rofId, 0), - mTimeFrame->getClustersOnLayer(rofId, 1), - mTimeFrame->getIndexTable(rofId, 0).data(), - mVrtParams.phiCut, - mTimeFrame->getTracklets()[0], - mTimeFrame->getNTrackletsCluster(rofId, 0), - mIndexTableUtils, - rofId, - 0, - mVrtParams.maxTrackletsPerCluster); - trackleterKernelHost( - mTimeFrame->getClustersOnLayer(rofId, 2), - mTimeFrame->getClustersOnLayer(rofId, 1), - mTimeFrame->getIndexTable(rofId, 2).data(), - mVrtParams.phiCut, - mTimeFrame->getTracklets()[1], - mTimeFrame->getNTrackletsCluster(rofId, 1), - mIndexTableUtils, - rofId, - 0, - mVrtParams.maxTrackletsPerCluster); - mTimeFrame->getNTrackletsROf(rofId, 0) = std::accumulate(mTimeFrame->getNTrackletsCluster(rofId, 0).begin(), mTimeFrame->getNTrackletsCluster(rofId, 0).end(), 0); - mTimeFrame->getNTrackletsROf(rofId, 1) = std::accumulate(mTimeFrame->getNTrackletsCluster(rofId, 1).begin(), mTimeFrame->getNTrackletsCluster(rofId, 1).end(), 0); + for (int pivotRofId = 0; pivotRofId < mTimeFrame->getNrof(); ++pivotRofId) { // Pivot rofId: the rof for which the tracklets are computed + int startROF{std::max(0, pivotRofId - mVrtParams.deltaRof)}; + int endROF{std::min(mTimeFrame->getNrof(), pivotRofId + mVrtParams.deltaRof + 1)}; + for (auto targetRofId = startROF; targetRofId < endROF; ++targetRofId) { + trackleterKernelHost( + mTimeFrame->getClustersOnLayer(targetRofId, 0), // Clusters to be matched with the next layer in target rof + mTimeFrame->getClustersOnLayer(pivotRofId, 1), // Clusters to be matched with the current layer in pivot rof + mTimeFrame->getUsedClustersROF(targetRofId, 0), // Span of the used clusters in the target rof + mTimeFrame->getIndexTable(targetRofId, 0).data(), // Index table to access the data on the next layer in target rof + mVrtParams.phiCut, + mTimeFrame->getTracklets()[0], // Flat tracklet buffer + mTimeFrame->getNTrackletsCluster(pivotRofId, 0), // Span of the number of tracklets per each cluster in pivot rof + mIndexTableUtils, + pivotRofId, + targetRofId, + gsl::span(), // Offset in the tracklet buffer + mVrtParams.maxTrackletsPerCluster); + trackleterKernelHost( + mTimeFrame->getClustersOnLayer(targetRofId, 2), + mTimeFrame->getClustersOnLayer(pivotRofId, 1), + mTimeFrame->getUsedClustersROF(targetRofId, 2), + mTimeFrame->getIndexTable(targetRofId, 2).data(), + mVrtParams.phiCut, + mTimeFrame->getTracklets()[1], + mTimeFrame->getNTrackletsCluster(pivotRofId, 1), // Span of the number of tracklets per each cluster in pivot rof + mIndexTableUtils, + pivotRofId, + targetRofId, + gsl::span(), // Offset in the tracklet buffer + mVrtParams.maxTrackletsPerCluster); + } + mTimeFrame->getNTrackletsROF(pivotRofId, 0) = std::accumulate(mTimeFrame->getNTrackletsCluster(pivotRofId, 0).begin(), mTimeFrame->getNTrackletsCluster(pivotRofId, 0).end(), 0); + mTimeFrame->getNTrackletsROF(pivotRofId, 1) = std::accumulate(mTimeFrame->getNTrackletsCluster(pivotRofId, 1).begin(), mTimeFrame->getNTrackletsCluster(pivotRofId, 1).end(), 0); } #pragma omp single - mTimeFrame->computeTrackletsScans(mNThreads); + mTimeFrame->computeTrackletsPerROFScans(); #pragma omp single mTimeFrame->getTracklets()[0].resize(mTimeFrame->getTotalTrackletsTF(0)); #pragma omp single mTimeFrame->getTracklets()[1].resize(mTimeFrame->getTotalTrackletsTF(1)); #pragma omp for schedule(dynamic) - for (int rofId = 0; rofId < mTimeFrame->getNrof(); ++rofId) { - trackleterKernelHost( - mTimeFrame->getClustersOnLayer(rofId, 0), - mTimeFrame->getClustersOnLayer(rofId, 1), - mTimeFrame->getIndexTable(rofId, 0).data(), - mVrtParams.phiCut, - mTimeFrame->getTracklets()[0], - mTimeFrame->getNTrackletsCluster(rofId, 0), - mIndexTableUtils, - rofId, - mTimeFrame->getNTrackletsROf(rofId, 0), - mVrtParams.maxTrackletsPerCluster); - trackleterKernelHost( - mTimeFrame->getClustersOnLayer(rofId, 2), - mTimeFrame->getClustersOnLayer(rofId, 1), - mTimeFrame->getIndexTable(rofId, 2).data(), - mVrtParams.phiCut, - mTimeFrame->getTracklets()[1], - mTimeFrame->getNTrackletsCluster(rofId, 1), - mIndexTableUtils, - rofId, - mTimeFrame->getNTrackletsROf(rofId, 1), - mVrtParams.maxTrackletsPerCluster); + for (int pivotRofId = 0; pivotRofId < mTimeFrame->getNrof(); ++pivotRofId) { + int startROF{std::max(0, pivotRofId - mVrtParams.deltaRof)}; + int endROF{std::min(mTimeFrame->getNrof(), pivotRofId + mVrtParams.deltaRof + 1)}; + for (auto targetRofId = startROF; targetRofId < endROF; ++targetRofId) { + trackleterKernelHost( + mTimeFrame->getClustersOnLayer(targetRofId, 0), + mTimeFrame->getClustersOnLayer(pivotRofId, 1), + mTimeFrame->getUsedClustersROF(targetRofId, 0), + mTimeFrame->getIndexTable(targetRofId, 0).data(), + mVrtParams.phiCut, + mTimeFrame->getTracklets()[0], + mTimeFrame->getNTrackletsCluster(pivotRofId, 0), + mIndexTableUtils, + pivotRofId, + targetRofId, + mTimeFrame->getExclusiveNTrackletsCluster(pivotRofId, 0), + mVrtParams.maxTrackletsPerCluster); + trackleterKernelHost( + mTimeFrame->getClustersOnLayer(targetRofId, 2), + mTimeFrame->getClustersOnLayer(pivotRofId, 1), + mTimeFrame->getUsedClustersROF(targetRofId, 2), + mTimeFrame->getIndexTable(targetRofId, 2).data(), + mVrtParams.phiCut, + mTimeFrame->getTracklets()[1], + mTimeFrame->getNTrackletsCluster(pivotRofId, 1), + mIndexTableUtils, + pivotRofId, + targetRofId, + mTimeFrame->getExclusiveNTrackletsCluster(pivotRofId, 1), + mVrtParams.maxTrackletsPerCluster); + } } } @@ -238,7 +270,7 @@ void VertexerTraits::computeTracklets() for (auto& trk : mTimeFrame->getTracklets()[0]) { MCCompLabel label; int sortedId0{mTimeFrame->getSortedIndex(trk.rof[0], 0, trk.firstClusterIndex)}; - int sortedId1{mTimeFrame->getSortedIndex(trk.rof[0], 1, trk.secondClusterIndex)}; + int sortedId1{mTimeFrame->getSortedIndex(trk.rof[1], 1, trk.secondClusterIndex)}; for (auto& lab0 : mTimeFrame->getClusterLabels(0, mTimeFrame->getClusters()[0][sortedId0].clusterId)) { for (auto& lab1 : mTimeFrame->getClusterLabels(1, mTimeFrame->getClusters()[1][sortedId1].clusterId)) { if (lab0 == lab1 && lab0.isValid()) { @@ -282,8 +314,14 @@ void VertexerTraits::computeTracklets() std::ofstream out01("NTC01_cpu.txt"), out12("NTC12_cpu.txt"); for (int iRof{0}; iRof < mTimeFrame->getNrof(); ++iRof) { + out01 << "ROF: " << iRof << std::endl; + out12 << "ROF: " << iRof << std::endl; std::copy(mTimeFrame->getNTrackletsCluster(iRof, 0).begin(), mTimeFrame->getNTrackletsCluster(iRof, 0).end(), std::ostream_iterator(out01, "\t")); + out01 << std::endl; + std::copy(mTimeFrame->getExclusiveNTrackletsCluster(iRof, 0).begin(), mTimeFrame->getExclusiveNTrackletsCluster(iRof, 0).end(), std::ostream_iterator(out01, "\t")); std::copy(mTimeFrame->getNTrackletsCluster(iRof, 1).begin(), mTimeFrame->getNTrackletsCluster(iRof, 1).end(), std::ostream_iterator(out12, "\t")); + out12 << std::endl; + std::copy(mTimeFrame->getExclusiveNTrackletsCluster(iRof, 1).begin(), mTimeFrame->getExclusiveNTrackletsCluster(iRof, 1).end(), std::ostream_iterator(out12, "\t")); out01 << std::endl; out12 << std::endl; } @@ -295,20 +333,30 @@ void VertexerTraits::computeTracklets() void VertexerTraits::computeTrackletMatching() { #pragma omp parallel for num_threads(mNThreads) schedule(dynamic) - for (int rofId = 0; rofId < mTimeFrame->getNrof(); ++rofId) { - mTimeFrame->getLines(rofId).reserve(mTimeFrame->getNTrackletsCluster(rofId, 0).size()); - trackletSelectionKernelHost( - mTimeFrame->getClustersOnLayer(rofId, 0), - mTimeFrame->getClustersOnLayer(rofId, 1), - mTimeFrame->getFoundTracklets(rofId, 0), - mTimeFrame->getFoundTracklets(rofId, 1), - mTimeFrame->getNTrackletsCluster(rofId, 0), - mTimeFrame->getNTrackletsCluster(rofId, 1), - mTimeFrame->getLines(rofId), - mTimeFrame->getLabelsFoundTracklets(rofId, 0), - mTimeFrame->getLinesLabel(rofId), - mVrtParams.tanLambdaCut, - mVrtParams.phiCut); + for (int pivotRofId = 0; pivotRofId < mTimeFrame->getNrof(); ++pivotRofId) { + mTimeFrame->getLines(pivotRofId).reserve(mTimeFrame->getNTrackletsCluster(pivotRofId, 0).size()); + std::vector usedTracklets(mTimeFrame->getFoundTracklets(pivotRofId, 0).size(), false); + int startROF{std::max(0, pivotRofId - mVrtParams.deltaRof)}; + int endROF{std::min(mTimeFrame->getNrof(), pivotRofId + mVrtParams.deltaRof + 1)}; + for (auto targetRofId = startROF; targetRofId < endROF; ++targetRofId) { + trackletSelectionKernelHost( + mTimeFrame->getClustersOnLayer(targetRofId, 0), + mTimeFrame->getClustersOnLayer(pivotRofId, 1), + mTimeFrame->getUsedClustersROF(targetRofId, 0), + mTimeFrame->getUsedClustersROF(targetRofId, 2), + mTimeFrame->getFoundTracklets(pivotRofId, 0), + mTimeFrame->getFoundTracklets(pivotRofId, 1), + usedTracklets, + mTimeFrame->getNTrackletsCluster(pivotRofId, 0), + mTimeFrame->getNTrackletsCluster(pivotRofId, 1), + mTimeFrame->getLines(pivotRofId), + mTimeFrame->getLabelsFoundTracklets(pivotRofId, 0), + mTimeFrame->getLinesLabel(pivotRofId), + pivotRofId, + targetRofId, + mVrtParams.tanLambdaCut, + mVrtParams.phiCut); + } } #ifdef VTX_DEBUG @@ -408,7 +456,7 @@ void VertexerTraits::computeVertices() std::array vertex2{}; for (int iCluster2{iCluster1 + 1}; iCluster2 < noClustersVec[rofId]; ++iCluster2) { vertex2 = mTimeFrame->getTrackletClusters(rofId)[iCluster2].getVertex(); - if (std::abs(vertex1[2] - vertex2[2]) < mVrtParams.clusterCut) { + if (o2::gpu::GPUCommonMath::Abs(vertex1[2] - vertex2[2]) < mVrtParams.clusterCut) { float distance{(vertex1[0] - vertex2[0]) * (vertex1[0] - vertex2[0]) + (vertex1[1] - vertex2[1]) * (vertex1[1] - vertex2[1]) + (vertex1[2] - vertex2[2]) * (vertex1[2] - vertex2[2])}; @@ -448,7 +496,7 @@ void VertexerTraits::computeVertices() } } - if (beamDistance2 < nsigmaCut && std::abs(mTimeFrame->getTrackletClusters(rofId)[iCluster].getVertex()[2]) < mVrtParams.maxZPositionAllowed) { + if (beamDistance2 < nsigmaCut && o2::gpu::GPUCommonMath::Abs(mTimeFrame->getTrackletClusters(rofId)[iCluster].getVertex()[2]) < mVrtParams.maxZPositionAllowed) { atLeastOneFound = true; vertices.emplace_back(o2::math_utils::Point3D(mTimeFrame->getTrackletClusters(rofId)[iCluster].getVertex()[0], mTimeFrame->getTrackletClusters(rofId)[iCluster].getVertex()[1], @@ -457,8 +505,7 @@ void VertexerTraits::computeVertices() // off-diagonal: square mean of projections on planes. mTimeFrame->getTrackletClusters(rofId)[iCluster].getSize(), // Contributors mTimeFrame->getTrackletClusters(rofId)[iCluster].getAvgDistance2()); // In place of chi2 - - vertices.back().setTimeStamp(rofId); + vertices.back().setTimeStamp(mTimeFrame->getTrackletClusters(rofId)[iCluster].getROF()); if (mTimeFrame->hasMCinformation()) { mTimeFrame->getVerticesLabels().emplace_back(); for (auto& index : mTimeFrame->getTrackletClusters(rofId)[iCluster].getLabels()) { @@ -467,7 +514,7 @@ void VertexerTraits::computeVertices() } } } - mTimeFrame->addPrimaryVertices(vertices); + mTimeFrame->addPrimaryVertices(vertices, rofId); } #ifdef VTX_DEBUG TFile* dbg_file = TFile::Open("artefacts_tf.root", "update"); @@ -569,7 +616,7 @@ void VertexerTraits::computeVerticesInRof(int rofId, std::array vertex2{}; for (int iCluster2{iCluster1 + 1}; iCluster2 < nClusters; ++iCluster2) { vertex2 = clusterLines[iCluster2].getVertex(); - if (std::abs(vertex1[2] - vertex2[2]) < mVrtParams.clusterCut) { + if (o2::gpu::GPUCommonMath::Abs(vertex1[2] - vertex2[2]) < mVrtParams.clusterCut) { float distance{(vertex1[0] - vertex2[0]) * (vertex1[0] - vertex2[0]) + (vertex1[1] - vertex2[1]) * (vertex1[1] - vertex2[1]) + (vertex1[2] - vertex2[2]) * (vertex1[2] - vertex2[2])}; @@ -602,7 +649,7 @@ void VertexerTraits::computeVerticesInRof(int rofId, continue; } } - if (beamDistance2 < nsigmaCut && std::abs(clusterLines[iCluster].getVertex()[2]) < mVrtParams.maxZPositionAllowed) { + if (beamDistance2 < nsigmaCut && o2::gpu::GPUCommonMath::Abs(clusterLines[iCluster].getVertex()[2]) < mVrtParams.maxZPositionAllowed) { atLeastOneFound = true; ++foundVertices; vertices.emplace_back(o2::math_utils::Point3D(clusterLines[iCluster].getVertex()[0], @@ -612,7 +659,7 @@ void VertexerTraits::computeVerticesInRof(int rofId, // off-diagonal: square mean of projections on planes. clusterLines[iCluster].getSize(), // Contributors clusterLines[iCluster].getAvgDistance2()); // In place of chi2 - vertices.back().setTimeStamp(rofId); + vertices.back().setTimeStamp(clusterLines[iCluster].getROF()); if (labels) { for (auto& index : clusterLines[iCluster].getLabels()) { labels->push_back(tf->getLinesLabel(rofId)[index]); // then we can use nContributors from vertices to get the labels diff --git a/Detectors/Upgrades/ITS3/reconstruction/src/IOUtils.cxx b/Detectors/Upgrades/ITS3/reconstruction/src/IOUtils.cxx index fd8a7d333037e..f930ae00f22c7 100644 --- a/Detectors/Upgrades/ITS3/reconstruction/src/IOUtils.cxx +++ b/Detectors/Upgrades/ITS3/reconstruction/src/IOUtils.cxx @@ -91,7 +91,7 @@ int loadROFrameDataITS3(its::TimeFrame* tf, tf->addClusterExternalIndexToLayer(layer, clusterId); } for (unsigned int iL{0}; iL < tf->getUnsortedClusters().size(); ++iL) { - tf->mROframesClusters[iL].push_back(tf->getUnsortedClusters()[iL].size()); + tf->mROFramesClusters[iL].push_back(tf->getUnsortedClusters()[iL].size()); } tf->mNrof++; } diff --git a/Detectors/Upgrades/ITS3/workflow/src/TrackerSpec.cxx b/Detectors/Upgrades/ITS3/workflow/src/TrackerSpec.cxx index b3f5a26f9d8be..d94f985bc32ab 100644 --- a/Detectors/Upgrades/ITS3/workflow/src/TrackerSpec.cxx +++ b/Detectors/Upgrades/ITS3/workflow/src/TrackerSpec.cxx @@ -246,7 +246,7 @@ void TrackerDPL::run(ProcessingContext& pc) for (auto& v : vtxVecLoc) { vertices.push_back(v); } - timeFrame->addPrimaryVertices(vtxVecLoc); + timeFrame->addPrimaryVertices(vtxVecLoc, iRof); } } // LOG(info) << fmt::format(" - rejected {}/{} ROFs: random/mult.sel:{} (seed {}), vtx.sel:{}", cutRandomMult + cutVertexMult, rofspan.size(), cutRandomMult, multEst.lastRandomSeed, cutVertexMult);