Skip to content
Closed
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
194 changes: 194 additions & 0 deletions Detectors/ITSMFT/ITS/tracking/GPU/ITStrackingGPU/ClusterLinesGPU.h
Original file line number Diff line number Diff line change
@@ -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<float>(-x0);
vertex[1] = static_cast<float>(-x1);
vertex[2] = static_cast<float>(-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<float>(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 */
134 changes: 133 additions & 1 deletion Detectors/ITSMFT/ITS/tracking/GPU/ITStrackingGPU/TimeFrameGPU.h
Original file line number Diff line number Diff line change
Expand Up @@ -21,6 +21,7 @@
#include "ITStracking/Configuration.h"
#include "ITStracking/TrackExtensionHypothesis.h"
#include "ITStrackingGPU/Utils.h"
#include "ITStrackingGPU/ClusterLinesGPU.h"

namespace o2::its::gpu
{
Expand Down Expand Up @@ -54,10 +55,14 @@ class TimeFrameGPU : public TimeFrame<NLayers>
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();
Expand Down Expand Up @@ -87,6 +92,35 @@ class TimeFrameGPU : public TimeFrame<NLayers>
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<NLayers>::prepareClusters(trkParam, maxLayers);
}
}
void allocateClusterSortStorage(const TrackingParameters& trkParam, const int maxLayers) override
{
if (maxLayers < NLayers) { // only if former seeding vertexer is run
TimeFrame<NLayers>::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; }
Expand All @@ -111,6 +145,19 @@ class TimeFrameGPU : public TimeFrame<NLayers>
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; }
Expand All @@ -119,6 +166,45 @@ class TimeFrameGPU : public TimeFrame<NLayers>
TrackExtensionHypothesis<NLayers>* 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<o2::MCCompLabel>& getLineLabelFlat() { return mLineLabelFlatHost; }
const std::vector<o2::MCCompLabel>& getLineLabelFlat() const { return mLineLabelFlatHost; }
Vertex* getDeviceDiamond() { return mDiamondDevice; }
std::array<CellNeighbour*, MaxCells>& getDeviceNeighboursAll() { return mNeighboursDevice; }
CellNeighbour* getDeviceNeighbours(const int layer) { return mNeighboursDevice[layer]; }
const TrackingFrameInfo** getDeviceArrayTrackingFrameInfo() const { return mTrackingFrameInfoDeviceArray; }
const Cluster** getDeviceArrayClusters() const { return mClustersDeviceArray; }
Expand Down Expand Up @@ -215,6 +301,11 @@ class TimeFrameGPU : public TimeFrame<NLayers>
const int** mClustersIndexTablesDeviceArray{nullptr};
uint8_t** mUsedClustersDeviceArray{nullptr};
const int** mROFramesClustersDeviceArray{nullptr};
int* mROFramesPVDevice;
std::array<int*, NLayers> mClusterSortKeysDevice{};
std::array<int*, NLayers> mClusterSortPermDevice{};
float* mClusterMinRDevice{nullptr};
float* mClusterMaxRDevice{nullptr};
std::array<Tracklet*, MaxLinks> mTrackletsDevice{};
std::array<int*, MaxLinks> mTrackletsLUTDevice{};
std::array<int*, MaxCells> mCellsLUTDevice{};
Expand All @@ -239,6 +330,47 @@ class TimeFrameGPU : public TimeFrame<NLayers>
CellNeighbour** mNeighboursDeviceArray{nullptr};
std::array<TrackingFrameInfo*, NLayers> mTrackingFrameInfoDevice{};
const TrackingFrameInfo** mTrackingFrameInfoDeviceArray{nullptr};
std::array<unsigned long long*, 3> 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<GPULine> mLinesHost;
std::vector<int> mLineRofHost;
std::vector<int> mLineClustersHost;
std::vector<VertexCand> mVertexCandsHost;
std::vector<int> mPeakOffsetsHost;
std::vector<int> mMemberOffsetsHost;
std::vector<int> mMemberLinesHost;
std::vector<o2::MCCompLabel> mLineLabelFlatHost;
int mNMembers{0};
Vertex* mDiamondDevice{nullptr};
bool mPersistentTablesLoaded{false};
std::bitset<NLayers> mUnsortedClustersUploaded{};
std::bitset<NLayers> mTrackingFrameInfoUploaded{};

// State
Streams mGpuStreams;
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -29,6 +29,9 @@ class TrackerTraitsGPU final : public TrackerTraits<NLayers>
void adoptTimeFrame(TimeFrame<NLayers>* 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;
Expand Down
Loading
Loading