ITS: debug dump

Author: f3sch

Detectors/ITSMFT/ITS/tracking/CMakeLists.txt

@@ -48,6 +48,7 @@ o2_add_library(ITSTrackingInterface
 o2_target_root_dictionary(ITStracking
                           HEADERS include/ITStracking/Tracklet.h
                                   include/ITStracking/Cluster.h
+                                  include/ITStracking/Seeding.h
                                   include/ITStracking/TrackingConfigParam.h
                           LINKDEF src/TrackingLinkDef.h)
 

Detectors/ITSMFT/ITS/tracking/include/ITStracking/Configuration.h

@@ -55,7 +55,7 @@ struct TrackingParameters {
   int MinTrackLength = 7;
   float NSigmaCut = 5;
   float NSigmaCutZ = 5;
-  float NSigmaCutPhi = 1;
+  float NSigmaCutPhi = 5;
   int NTimeSlices{1};
   /// Trackleting cuts
   float TrackletMinPt = 0.3f;

Detectors/ITSMFT/ITS/tracking/include/ITStracking/Definitions.h

@@ -30,7 +30,29 @@ using TimeEstBC = o2::dataformats::TimeStampWithError<uint32_t, uint16_t>;
 using Vertex = o2::dataformats::Vertex<TimeEstBC>;
 
 template <bool IsConst, typename T>
-using maybe_const = typename std::conditional<IsConst, const T, T>::type;
+using maybe_const = std::conditional_t<IsConst, const T, T>;
+
+// simple implemnetion of logging with exp. backoff
+struct LogLogThrottler {
+  uint32_t evCount{0};
+  uint32_t nextLog{1};
+  int32_t iteration{-1};
+  int32_t layer{-1};
+  bool needToLog(int32_t iter, int32_t lay)
+  {
+    if (iteration != iter || layer != lay) {
+      iteration = iter;
+      layer = lay;
+      evCount = 0;
+      nextLog = 1;
+    }
+    if (++evCount > nextLog) {
+      nextLog *= 2;
+      return true;
+    }
+    return false;
+  }
+};
 
 } // namespace o2::its
 

Detectors/ITSMFT/ITS/tracking/include/ITStracking/Seeding.h

@@ -26,6 +26,7 @@
 #include "ITStracking/Constants.h"
 #include "GPUCommonDef.h"
 #include "GPUCommonMath.h"
+#include "GPUCommonRtypes.h"
 #include "GPUCommonLogger.h"
 
 namespace o2::its
@@ -57,7 +58,7 @@ struct LinearizedTrack {
 
   float getChi2(const dataformats::VertexBase&) const;
 
-  int32_t cellIdx{-1}; /// < index of attached cell
+  int32_t cellIdx{-1}; ///< index of attached cell
   float x{0.f};        ///< reference X
   float y{0.f};        ///< Y at X
   float z{0.f};        ///< Z at X
@@ -68,7 +69,7 @@ struct LinearizedTrack {
   float tgL{0.f};      ///< tangent(lambda)
   float cosAlp{0.f};   ///< cos of alpha frame
   float sinAlp{0.f};   ///< sin of alpha frame
-  TimeEstBC time{};    ///< timestamp; exclusive range of BCs that are spanned [time, time+error)
+  TimeEstBC time;      ///< timestamp; exclusive range of BCs that are spanned [time, time+error)
   uint8_t status{0};   ///< status bits
 
   GPUhdi() void markDead() noexcept { status |= kDead; }
@@ -78,6 +79,8 @@ struct LinearizedTrack {
   {
     LOGP(info, "LT: cell={} x={} y={} z={} t={}+{} status={:b}", cellIdx, x, y, z, time.getTimeStamp(), time.getTimeStampError(), status);
   }
+
+  ClassDefNV(LinearizedTrack, 1);
 };
 
 struct VertexSeed : public Vertex {
@@ -116,6 +119,7 @@ constexpr int32_t NDim{2};
 struct DBSCANParams {
   std::array<float, NDim> eps; // Maximum distance per dimension (z,t)
   int32_t minPts;              // Minimum points to form a dense region
+  size_t avgNeigh{1024};       // average number of neighbors
 };
 
 // Clustering result
@@ -127,19 +131,8 @@ struct DBSCANResult {
   int32_t nNoise = 0;
 };
 
-// Flat neighbor list
-struct FlatNeighborList {
-  std::vector<size_t> indices;
-  std::vector<size_t> offsets;
-  [[nodiscard]] size_t getSize(size_t i) const
-  {
-    return offsets[i + 1] - offsets[i];
-  }
-  [[nodiscard]] std::span<const size_t> getNeighbors(size_t i) const
-  {
-    return {&indices[offsets[i]], getSize(i)};
-  }
-};
+// Neighbor list
+using NeighborList = std::vector<std::vector<size_t>>;
 
 // Point classification
 enum DBSCANLabel : int32_t {
@@ -149,6 +142,29 @@ enum DBSCANLabel : int32_t {
   DB_CORE = -(1 << 3),
 };
 
+// #define ITS_DB_MEASURE_TIMING
+#ifdef ITS_DB_MEASURE_TIMING
+class ScopedTimer
+{
+  std::string_view name;
+  std::chrono::high_resolution_clock::time_point start;
+
+ public:
+  explicit ScopedTimer(std::string_view name)
+    : name(name), start(std::chrono::high_resolution_clock::now()) {}
+
+  ~ScopedTimer()
+  {
+    auto end = std::chrono::high_resolution_clock::now();
+    double elapsed_ms = std::chrono::duration<double, std::milli>(end - start).count();
+    LOGP(info, "{} : {:.2f} ms", name, elapsed_ms);
+  }
+};
+#define SCOPED_TIMER(name) ScopedTimer _timer##__LINE__(name)
+#else
+#define SCOPED_TIMER(name) ((void)0)
+#endif
+
 class DBSCANDistance
 {
   using EPS = decltype(DBSCANParams::eps);
@@ -172,20 +188,21 @@ class DBSCANDistance
     const int end0 = l0.time.getTimeStamp() + l0.time.getTimeStampError();
     const int end1 = l1.time.getTimeStamp() + l1.time.getTimeStampError();
     const int diffT = o2::gpu::CAMath::Max(0, o2::gpu::CAMath::Max((int)l0.time.getTimeStamp(), (int)l1.time.getTimeStamp())) - o2::gpu::CAMath::Min(end0, end1);
-    return diffT <= mEps[1];
+    return diffT <= (int)mEps[1];
   }
 
-  // Batch compute
-  void computeNeighbors(const LinearizedTrack& query, const LinearizedTrack* tracks, const std::vector<size_t>& candidates, std::vector<size_t>& neighbors) const
+  float getDistance(const LinearizedTrack& l0, const LinearizedTrack& l1) const
   {
-    neighbors.clear();
-    neighbors.reserve(candidates.size());
-    for (auto idx : candidates) {
-      const auto& lt = tracks[idx];
-      if (areNeighbours(query, lt)) {
-        neighbors.push_back(idx);
-      }
+    if (l0.isDead() || l1.isDead()) {
+      return std::numeric_limits<float>::max();
+    }
+    const int end0 = l0.time.getTimeStamp() + l0.time.getTimeStampError();
+    const int end1 = l1.time.getTimeStamp() + l1.time.getTimeStampError();
+    const int diffT = o2::gpu::CAMath::Max(0, o2::gpu::CAMath::Max((int)l0.time.getTimeStamp(), (int)l1.time.getTimeStamp())) - o2::gpu::CAMath::Min(end0, end1);
+    if (diffT > (int)mEps[1]) {
+      return std::numeric_limits<float>::max();
     }
+    return o2::gpu::CAMath::Abs(l0.z - l1.z);
   }
 
  private:
@@ -199,11 +216,13 @@ class DBSCAN
   DBSCAN(const DBSCANParams& p);
 
   DBSCANResult cluster(const LinearizedTrack* tracks, size_t n);
+  const auto& getParams() { return mParams; }
+  const auto& getDistance() { return mDistance; }
 
  private:
-  void findNeighbors(const LinearizedTrack* tracks, size_t n, FlatNeighborList& neighbors);
-  void classify(size_t n, const FlatNeighborList& neighbors, std::vector<int32_t>& labels);
-  void expandCluster(size_t i, int32_t idx, const FlatNeighborList& neighbors, std::vector<int32_t>& labels) const;
+  void findNeighbors(const LinearizedTrack* tracks, size_t n, NeighborList& neighbors);
+  void classify(size_t n, const NeighborList& neighbors, std::vector<int32_t>& labels);
+  void expandCluster(size_t i, int32_t idx, const NeighborList& neighbors, std::vector<int32_t>& labels) const;
 
   DBSCANParams mParams;
   DBSCANDistance mDistance;
@@ -220,6 +239,7 @@ class Grid
   Grid(const LinearizedTrack* tracks, size_t n, const std::array<float, NDim>& cellSizes)
     : mTracks(tracks), mNPoints(n), mCellSizes(cellSizes)
   {
+    SCOPED_TIMER("Grid construction");
     computeBounds();
     computeGridDimensions();
     allocateCells();

Detectors/ITSMFT/ITS/tracking/include/ITStracking/TrackerTraits.h

@@ -16,7 +16,9 @@
 #ifndef TRACKINGITSU_INCLUDE_TRACKERTRAITS_H_
 #define TRACKINGITSU_INCLUDE_TRACKERTRAITS_H_
 
-#include <oneapi/tbb.h>
+#include <oneapi/tbb/task_arena.h>
+
+#include <type_traits>
 
 #include "DetectorsBase/Propagator.h"
 #include "ITStracking/Configuration.h"
@@ -27,8 +29,6 @@
 #include "ITStracking/Cell.h"
 #include "ITStracking/BoundedAllocator.h"
 
-// #define OPTIMISATION_OUTPUT
-
 namespace o2
 {
 namespace gpu
@@ -42,16 +42,20 @@ class TrackITSExt;
 template <int NLayers = 7>
 class TrackerTraits
 {
+  struct Neighbor {
+    int32_t cell{-1}, nextCell{-1}, level{-1};
+  };
+
  public:
   using IndexTableUtilsN = IndexTableUtils<NLayers>;
   using CellSeedN = CellSeed<NLayers>;
 
-  TrackerTraits() = default;
+  TrackerTraits();
   TrackerTraits(const TrackerTraits&) = delete;
   TrackerTraits(TrackerTraits&&) = delete;
   TrackerTraits& operator=(const TrackerTraits&) = delete;
   TrackerTraits& operator=(TrackerTraits&&) = delete;
-  virtual ~TrackerTraits() = default;
+  virtual ~TrackerTraits();
 
   virtual void adoptTimeFrame(TimeFrame<NLayers>* tf) { mTimeFrame = tf; }
   virtual void initialiseTimeFrame(const RecoIteration& reco) { mTimeFrame->initialise(reco); }
@@ -106,14 +110,19 @@ class TrackerTraits
   std::shared_ptr<BoundedMemoryResource> mMemoryPool;
   std::shared_ptr<tbb::task_arena> mTaskArena;
 
+  // Debug output
+  inline void debugComputeLayerTracklets(int iteration, int layer, const Cluster& currentCls, const Cluster& nextCls, const Vertex& pv, float sigmaZ, float sigmaPhi, bool accepted);
+  inline void debugComputeLayerCells(int iteration, int layer, int currentTrkl, int nextTrkl, bool accepted);
+  inline void debugFindCellsNeighbours(int iteration, int layer, int currentCell, int nextCell, bool accepted);
+
  protected:
   o2::gpu::GPUChainITS* mChain = nullptr;
   TimeFrame<NLayers>* mTimeFrame;
   std::vector<RecoIteration> mRecoParams;
 
   float mBz{-999.f};
   bool mIsZeroField{false};
-  dbscan::DBSCAN mDBScan;
+  dbscan::DBSCAN mDBScan{};
 };
 
 template <int NLayers>

Detectors/ITSMFT/ITS/tracking/include/ITStracking/TrackingConfigParam.h

@@ -66,7 +66,7 @@ struct TrackerParamConfig : public o2::conf::ConfigurableParamHelper<TrackerPara
   int seedingMaxFitIter{3};                // maximum iterations for fit
   int seedingMinTracksIter{50};            // minimum tracks needed for one iteration
   int seedingDBScanMinPt{100};             // DBSCAN: minimum number of cluster points
-  float seedingDBScanEpsZ{0.1f};           // DBSCAN: maximum epsilon z
+  float seedingDBScanEpsZ{0.2f};           // DBSCAN: maximum epsilon z
   float seedingDBScanEpsT{10.f};           // DBSCAN: maximum epsilon t (BC)
   float seedingVertexExtraErr2[6] = {0.f}; // impose additional errors on seeding vertices
   bool seedingUseMCTruth{false};           // skip seeding and impose MC event information

Detectors/ITSMFT/ITS/tracking/src/Seeding.cxx

@@ -150,11 +150,12 @@ DBSCANResult DBSCAN::cluster(const LinearizedTrack* tracks, size_t n)
   }
 
   // Step 1: Find neighbors for all points using grid
-  FlatNeighborList neighbors;
+  NeighborList neighbors;
   findNeighbors(tracks, n, neighbors);
   // Step 2: Classify points and form clusters
   classify(n, neighbors, result.labels);
   // Step 3: Count clusters and noise points
+  SCOPED_TIMER("finalizing");
   int32_t max_label = *std::ranges::max_element(result.labels);
   result.nClusters = max_label + 1;
   result.nNoise = static_cast<int32_t>(std::count(result.labels.begin(), result.labels.end(), DB_NOISE));
@@ -184,49 +185,39 @@ DBSCANResult DBSCAN::cluster(const LinearizedTrack* tracks, size_t n)
   return result;
 }
 
-void DBSCAN::findNeighbors(const LinearizedTrack* tracks, size_t n, FlatNeighborList& neighbors)
+void DBSCAN::findNeighbors(const LinearizedTrack* tracks, size_t n, NeighborList& neighbors)
 {
   Grid grid(tracks, n, mParams.eps);
-  std::vector<std::vector<size_t>> temp(n);
 
   // Parallel neighbor finding
+  SCOPED_TIMER("find neighbors");
+  neighbors.resize(n);
   tbb::parallel_for(
     tbb::blocked_range<size_t>(0, n), [&](const tbb::blocked_range<size_t>& range) {
-      std::vector<const GridCell*> neighbor_cells;
-      std::vector<size_t> candidates;
+      std::vector<const GridCell*> neighborCells;
 
       for (size_t i = range.begin(); i < range.end(); ++i) {
         const auto& query = tracks[i];
         if (query.isDead()) {
           continue;
         }
         auto coords = grid.getGridCoords(i);
-        grid.getNeighborCells(coords, neighbor_cells);
-        candidates.clear();
-        for (const GridCell* cell : neighbor_cells) {
+        grid.getNeighborCells(coords, neighborCells);
+        neighbors[i].reserve(mParams.avgNeigh);
+        for (const GridCell* cell : neighborCells) {
           for (auto idx : *cell) {
-            if (idx != i) {
-              candidates.push_back(idx);
+            if (idx != i && mDistance.areNeighbours(query, tracks[idx])) {
+              neighbors[i].push_back(idx);
             }
           }
         }
-        mDistance.computeNeighbors(query, tracks, candidates, temp[i]);
       }
     });
-  // flatten
-  neighbors.offsets.resize(n + 1);
-  neighbors.offsets[0] = 0;
-  for (size_t i = 0; i < n; ++i) {
-    neighbors.offsets[i + 1] = neighbors.offsets[i] + temp[i].size();
-  }
-  neighbors.indices.resize(neighbors.offsets[n]);
-  for (size_t i = 0; i < n; ++i) {
-    std::copy(temp[i].begin(), temp[i].end(), &neighbors.indices[neighbors.offsets[i]]);
-  }
 }
 
-void DBSCAN::classify(size_t n, const FlatNeighborList& neighbors, std::vector<int32_t>& labels)
+void DBSCAN::classify(size_t n, const NeighborList& neighbors, std::vector<int32_t>& labels)
 {
+  SCOPED_TIMER("classify");
   std::vector<bool> isCore(n, false);
   int32_t nextClsIdx{0};
 
@@ -235,7 +226,7 @@ void DBSCAN::classify(size_t n, const FlatNeighborList& neighbors, std::vector<i
     tbb::blocked_range<size_t>(0, n),
     [&](const tbb::blocked_range<size_t>& range) {
       for (size_t i = range.begin(); i < range.end(); ++i) {
-        if (static_cast<int32_t>(neighbors.getSize(i)) >= mParams.minPts) {
+        if (static_cast<int32_t>(neighbors[i].size()) >= mParams.minPts) {
           isCore[i] = true;
         }
       }
@@ -263,7 +254,7 @@ void DBSCAN::classify(size_t n, const FlatNeighborList& neighbors, std::vector<i
     });
 }
 
-void DBSCAN::expandCluster(size_t i, int32_t idx, const FlatNeighborList& neighbors, std::vector<int32_t>& labels) const
+void DBSCAN::expandCluster(size_t i, int32_t idx, const NeighborList& neighbors, std::vector<int32_t>& labels) const
 {
   std::queue<size_t> seeds;
   seeds.push(i);
@@ -274,8 +265,8 @@ void DBSCAN::expandCluster(size_t i, int32_t idx, const FlatNeighborList& neighb
     seeds.pop();
 
     // If current is a core point, add its neighbors to the cluster
-    if (static_cast<int32_t>(neighbors.getSize(current)) >= mParams.minPts) {
-      for (auto neighbor : neighbors.getNeighbors(current)) {
+    if (static_cast<int32_t>(neighbors[current].size()) >= mParams.minPts) {
+      for (const auto neighbor : neighbors[current]) {
         if (labels[neighbor] == DB_UNVISITED) {
           labels[neighbor] = idx;
           seeds.push(neighbor);