ITS: fix rare valid lookup at the edge of acceptance

Detectors/ITSMFT/ITS/macros/test/CheckStaggering.C

@@ -80,9 +80,9 @@ void CheckStaggering(int runNumber, int max = -1, const std::string& dir = "")
   auto& ccdbmgr = o2::ccdb::BasicCCDBManager::instance();
   ccdbmgr.setURL("https://alice-ccdb.cern.ch");
   auto runDuration = ccdbmgr.getRunDuration(runNumber);
-  auto tRun = runDuration.first + (runDuration.second - runDuration.first) / 2; // time stamp for the middle of the run duration
+  auto tRun = runDuration.first + ((runDuration.second - runDuration.first) / 2); // time stamp for the middle of the run duration
   ccdbmgr.setTimestamp(tRun);
-  printf("Run %d has TS %lld", runNumber, tRun);
+  printf("Run %d has TS %ld", runNumber, tRun);
   auto geoAligned = ccdbmgr.get<TGeoManager>("GLO/Config/GeometryAligned");
   auto magField = ccdbmgr.get<o2::parameters::GRPMagField>("GLO/Config/GRPMagField");
   auto grpLHC = ccdbmgr.get<o2::parameters::GRPLHCIFData>("GLO/Config/GRPLHCIF");
@@ -134,7 +134,7 @@ void CheckStaggering(int runNumber, int max = -1, const std::string& dir = "")
   auto hVtxZ = new TH1F("hVtxZ", "seeding vertices Z", 200, -16, 16);
   auto hVtxNCont = new TH1F("hVtxNCont", "seeding vertices contributors", 100, 0, 100);
   auto hVtxZNCont = new TProfile("hVtxZNCont", "seeding vertices z-contributors", 200, -16, 16);
-  auto hVtxCls = new TProfile("hVtxCls", ";Cls/TF;Cls/Vtx", 400, 20000, 60000);
+  auto hVtxCls = new TProfile("hVtxCls", ";Cls/TF;Cls/Vtx", 2000, 600000, 900000);
   auto hVtxTS = new TH1D("hVtxTS", "vtx time t0;t0 (BC)", o2::constants::lhc::LHCMaxBunches, 0, o2::constants::lhc::LHCMaxBunches);
 
   const float minVtxWeight{5};
@@ -183,8 +183,12 @@ void CheckStaggering(int runNumber, int max = -1, const std::string& dir = "")
 
     tTrks->SetBranchAddress("ITSTrack", &trkArrPtr);
     tTrks->SetBranchAddress("Vertices", &vtxArrPtr);
-    for (int i{0}; i < 7; ++i) {
-      tCls->SetBranchAddress(Form("ITSClusterComp_%d", i), &clsArr[i]);
+    if (tCls->GetBranchStatus("ITSClusterComp")) {
+      tCls->SetBranchAddress("ITSClusterComp", &clsArr[0]);
+    } else {
+      for (int i{0}; i < 7; ++i) {
+        tCls->SetBranchAddress(Form("ITSClusterComp_%d", i), &clsArr[i]);
+      }
     }
 
     for (int iTF{0}; tTrks->LoadTree(iTF) >= 0; ++iTF) {
@@ -193,7 +197,9 @@ void CheckStaggering(int runNumber, int max = -1, const std::string& dir = "")
 
       size_t ncls = 0;
       for (int i{0}; i < 7; ++i) {
-        ncls += clsArr[i]->size();
+        if (clsArr[i]) {
+          ncls += clsArr[i]->size();
+        }
       }
 
       // for each TF built pool of positive and negaitve tracks

Detectors/ITSMFT/ITS/tracking/GPU/cuda/TrackingKernels.cu

@@ -47,20 +47,6 @@ namespace o2::its
 namespace gpu
 {
 
-struct sort_tracklets {
-  GPUhd() bool operator()(const Tracklet& a, const Tracklet& b)
-  {
-    if (a.firstClusterIndex != b.firstClusterIndex) {
-      return a.firstClusterIndex < b.firstClusterIndex;
-    }
-    return a.secondClusterIndex < b.secondClusterIndex;
-  }
-};
-
-struct equal_tracklets {
-  GPUhd() bool operator()(const Tracklet& a, const Tracklet& b) { return a.firstClusterIndex == b.firstClusterIndex && a.secondClusterIndex == b.secondClusterIndex; }
-};
-
 template <typename T1, typename T2>
 struct sort_by_second {
   GPUhd() bool operator()(const gpuPair<T1, T2>& a, const gpuPair<T1, T2>& b) const { return a.second < b.second; }
@@ -408,7 +394,7 @@ GPUg() void __launch_bounds__(256, 1) computeLayerTrackletsMultiROFKernel(
         const float sqInverseDeltaZ0{1.f / (math_utils::Sq(currentCluster.zCoordinate - primaryVertex.getZ()) + constants::Tolerance)}; /// protecting from overflows adding the detector resolution
         const float sigmaZ{o2::gpu::CAMath::Sqrt(math_utils::Sq(resolution) * math_utils::Sq(tanLambda) * ((math_utils::Sq(inverseR0) + sqInverseDeltaZ0) * math_utils::Sq(meanDeltaR) + 1.f) + math_utils::Sq(meanDeltaR * MSAngle))};
         const int4 selectedBinsRect{o2::its::getBinsRect(currentCluster, layerIndex + 1, zAtRmin, zAtRmax, sigmaZ * NSigmaCut, phiCut, *utils)};
-        if (selectedBinsRect.x == 0 && selectedBinsRect.y == 0 && selectedBinsRect.z == 0 && selectedBinsRect.w == 0) {
+        if (selectedBinsRect.x < 0) {
           continue;
         }
         int phiBinsNum{selectedBinsRect.w - selectedBinsRect.y + 1};
@@ -687,8 +673,8 @@ void computeTrackletsInROFsHandler(const IndexTableUtils<NLayers>* utils,
     mulScatAng[layer]);
   thrust::device_ptr<Tracklet> tracklets_ptr(spanTracklets[layer]);
   auto nosync_policy = THRUST_NAMESPACE::par_nosync(gpu::TypedAllocator<char>(alloc)).on(streams[layer].get());
-  thrust::sort(nosync_policy, tracklets_ptr, tracklets_ptr + nTracklets[layer], gpu::sort_tracklets());
-  auto unique_end = thrust::unique(nosync_policy, tracklets_ptr, tracklets_ptr + nTracklets[layer], gpu::equal_tracklets());
+  thrust::sort(nosync_policy, tracklets_ptr, tracklets_ptr + nTracklets[layer]);
+  auto unique_end = thrust::unique(nosync_policy, tracklets_ptr, tracklets_ptr + nTracklets[layer]);
   nTracklets[layer] = unique_end - tracklets_ptr;
   if (layer) {
     GPUChkErrS(cudaMemsetAsync(trackletsLUTsHost[layer], 0, (nClusters[layer] + 1) * sizeof(int), streams[layer].get()));

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

@@ -56,18 +56,17 @@ class SeedBase : public o2::track::TrackParCovF
   GPUhd() const auto& clustersRaw() const { return mClusters; }
 
  private:
-  float mChi2 = -999.f;
-  int mLevel = constants::UnusedIndex;
+  float mChi2{constants::UnsetValue};
+  int mLevel{constants::UnusedIndex};
   std::array<int, 2> mTracklets = constants::helpers::initArray<int, 2, constants::UnusedIndex>();
   std::array<int, NClusters> mClusters = constants::helpers::initArray<int, NClusters, constants::UnusedIndex>();
   TimeEstBC mTime;
 };
 
 /// CellSeed: connections of three clusters
-class CellSeed final : public SeedBase<3>
+class CellSeed final : public SeedBase<constants::ClustersPerCell>
 {
-  static constexpr int NStoredClusters = 3;
-  using Base = SeedBase<NStoredClusters>;
+  using Base = SeedBase<constants::ClustersPerCell>;
 
  public:
   GPUhdDefault() CellSeed() = default;
@@ -98,7 +97,7 @@ class CellSeed final : public SeedBase<3>
   GPUhd() int getCluster(int layer) const
   {
     const int rel = layer - getInnerLayer();
-    return (rel >= 0 && rel < NStoredClusters) ? this->clustersRaw()[rel] : constants::UnusedIndex;
+    return (rel >= 0 && rel < constants::ClustersPerCell) ? this->clustersRaw()[rel] : constants::UnusedIndex;
   }
 };
 

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

@@ -19,9 +19,6 @@
 #include <array>
 #include <utility>
 
-#include "GPUCommonDef.h"
-#include "GPUCommonDefAPI.h"
-
 namespace o2::its::constants
 {
 
@@ -30,12 +27,14 @@ constexpr float MB = KB * KB;
 constexpr float GB = MB * KB;
 constexpr bool DoTimeBenchmarks = true;
 constexpr bool SaveTimeBenchmarks = false;
-constexpr float Tolerance = 1e-12; // numerical tolerance
-constexpr int ClustersPerCell = 3;
-constexpr int UnusedIndex = -1;
-constexpr float Radl = 9.36f; // Radiation length of Si [cm]
-constexpr float Rho = 2.33f;  // Density of Si [g/cm^3]
-constexpr int MaxIter = 4;    // Max. supported iterations
+constexpr float Tolerance = 1e-12;                  // numerical tolerance
+constexpr int ClustersPerCell = 3;                  // number of clusters for a cell
+constexpr int UnusedIndex = -1;                     // global unused flag
+constexpr float UnsetValue = -999.f;                // global unset value
+constexpr float Radl = 9.36f;                       // Radiation length of Si [cm]
+constexpr float Rho = 2.33f;                        // Density of Si [g/cm^3]
+constexpr int MaxIter = 4;                          // Max. supported iterations
+constexpr int MaxSelectedTrackletsPerCluster = 100; // vertexer: max lines per cluster
 
 namespace helpers
 {

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

@@ -50,6 +50,19 @@ struct LogLogThrottler {
     return false;
   }
 };
+
+struct TimingStats {
+  std::uint64_t calls = 0;
+  double totalTimeMs = 0.;
+
+  void add(double timeMs)
+  {
+    ++calls;
+    totalTimeMs += timeMs;
+  }
+  double averageTimeMs() const { return calls ? totalTimeMs / static_cast<double>(calls) : 0.; }
+};
+
 } // namespace o2::its
 
-#endif
+#endif

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

@@ -124,7 +124,7 @@ GPUhdi() int4 getBinsRect(const Cluster& currentCluster, const int layerIndex,
 
   if (zRangeMax < -utils.getLayerZ(layerIndex) ||
       zRangeMin > utils.getLayerZ(layerIndex) || zRangeMin > zRangeMax) {
-    return int4{0, 0, 0, 0};
+    return int4{-1, -1, -1, -1};
   }
 
   return int4{o2::gpu::GPUCommonMath::Max(0, utils.getZBinIndex(layerIndex, zRangeMin)),

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

@@ -89,6 +89,12 @@ GPUhdi() float smallestAngleDifference(float a, float b)
   return o2::gpu::CAMath::Remainderf(b - a, o2::constants::math::TwoPI);
 }
 
+GPUhdi() bool isPhiDifferenceBelow(const float phiA, const float phiB, const float phiCut)
+{
+  const float deltaPhi = o2::gpu::CAMath::Abs(phiA - phiB);
+  return deltaPhi < phiCut || deltaPhi > o2::constants::math::TwoPI - phiCut;
+}
+
 GPUhdi() constexpr float Sq(float v)
 {
   return v * v;

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

@@ -26,9 +26,7 @@
 #include "ITStracking/Cell.h"
 #include "ITStracking/Cluster.h"
 #include "ITStracking/Configuration.h"
-#include "ITStracking/Constants.h"
 #include "ITStracking/ClusterLines.h"
-#include "ITStracking/Definitions.h"
 #include "ITStracking/Tracklet.h"
 #include "ITStracking/IndexTableUtils.h"
 #include "ITStracking/ExternalAllocator.h"
@@ -103,7 +101,7 @@ struct TimeFrame {
   void setBeamPosition(const float x, const float y, const float s2, const float base = 50.f, const float systematic = 0.f)
   {
     isBeamPositionOverridden = true;
-    resetBeamXY(x, y, s2 / o2::gpu::CAMath::Sqrt(base * base + systematic));
+    resetBeamXY(x, y, s2 / o2::gpu::CAMath::Sqrt((base * base) + systematic));
   }
 
   float getBeamX() const { return mBeamPos[0]; }
@@ -249,7 +247,7 @@ struct TimeFrame {
 
   // Propagator
   const o2::base::PropagatorImpl<float>* getDevicePropagator() const { return mPropagatorDevice; }
-  virtual void setDevicePropagator(const o2::base::PropagatorImpl<float>*) {};
+  virtual void setDevicePropagator(const o2::base::PropagatorImpl<float>* /*unused*/) {};
 
   template <typename... T>
   void addClusterToLayer(int layer, T&&... args);

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

@@ -25,10 +25,12 @@
 #include <string_view>
 #include <utility>
 #include <sstream>
+#include <vector>
 
 #include <oneapi/tbb/task_arena.h>
 
 #include "ITStracking/Configuration.h"
+#include "ITStracking/Definitions.h"
 #include "ITStracking/TimeFrame.h"
 #include "ITStracking/TrackerTraits.h"
 #include "ITStracking/BoundedAllocator.h"
@@ -91,16 +93,18 @@ class Tracker
   double mTotalTime{0};
   std::shared_ptr<BoundedMemoryResource> mMemoryPool;
 
-  enum State {
+  enum Steps {
     TFInit = 0,
     Trackleting,
     Celling,
     Neighbouring,
     Roading,
-    NStates,
+    NSteps,
   };
-  State mCurState{TFInit};
-  static constexpr std::array<const char*, NStates> StateNames{"TimeFrame initialisation", "Tracklet finding", "Cell finding", "Neighbour finding", "Road finding"};
+  Steps mCurStep{TFInit};
+  static constexpr std::array<const char*, NSteps> StateNames{"TimeFrame initialisation", "Tracklet finding", "Cell finding", "Neighbour finding", "Road finding"};
+  std::vector<std::array<TimingStats, NSteps>> mTimingStats;
+  void addTimingStatCurStep(int iteration, double timeMs);
 };
 
 template <int NLayers>
@@ -125,22 +129,23 @@ float Tracker<NLayers>::evaluateTask(void (Tracker<NLayers>::*task)(T...), std::
     }
     logger(sstream.str());
 
-    if (mTrkParams[0].SaveTimeBenchmarks) {
+    if (mTrkParams[iteration].SaveTimeBenchmarks) {
       std::string taskNameStr(taskName);
       std::transform(taskNameStr.begin(), taskNameStr.end(), taskNameStr.begin(),
                      [](unsigned char c) { return std::tolower(c); });
       std::replace(taskNameStr.begin(), taskNameStr.end(), ' ', '_');
       if (std::ofstream file{"its_time_benchmarks.txt", std::ios::app}) {
         file << "trk:" << iteration << '\t' << taskNameStr << '\t' << diff << '\n';
       }
+      addTimingStatCurStep(iteration, diff);
     }
 
   } else {
     (this->*task)(std::forward<F>(args)...);
   }
 
   if (mTrkParams[iteration].PrintMemory) {
-    LOGP(info, "iter:{}:{}: {}", iteration, StateNames[mCurState], mMemoryPool->asString());
+    LOGP(info, "iter:{}:{}: {}", iteration, StateNames[mCurStep], mMemoryPool->asString());
   }
 
   return diff;

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

@@ -24,27 +24,35 @@
 #include "GPUCommonDef.h"
 #include "GPUCommonLogger.h"
 
-#ifndef GPUCA_GPUCODE_DEVICE
-#ifndef GPU_NO_FMT
-#include <string>
-#include <fmt/format.h>
-#endif
-#endif
-
 namespace o2::its
 {
 
+// tracklets are entirely determined by their two cluster idx
 struct Tracklet final {
   GPUhdDefault() Tracklet() = default;
-  GPUhdi() Tracklet(const int, const int, const Cluster&, const Cluster&, const TimeEstBC& t);
-  GPUhdi() Tracklet(const int, const int, float tanL, float phi, const TimeEstBC& t);
-  GPUhdDefault() bool operator==(const Tracklet&) const = default;
-  GPUhdi() unsigned char isEmpty() const
+  GPUhdi() Tracklet(const int firstClusterOrderingIndex, const int secondClusterOrderingIndex,
+                    const Cluster& firstCluster, const Cluster& secondCluster, const TimeEstBC& t)
+    : firstClusterIndex(firstClusterOrderingIndex),
+      secondClusterIndex(secondClusterOrderingIndex),
+      tanLambda((firstCluster.zCoordinate - secondCluster.zCoordinate) / (firstCluster.radius - secondCluster.radius)),
+      phi(o2::gpu::GPUCommonMath::ATan2(firstCluster.yCoordinate - secondCluster.yCoordinate, firstCluster.xCoordinate - secondCluster.xCoordinate)),
+      mTime(t) {}
+
+  GPUhdi() Tracklet(const int idx0, const int idx1, float tanL, float phi, const TimeEstBC& t)
+    : firstClusterIndex(idx0),
+      secondClusterIndex(idx1),
+      tanLambda(tanL),
+      phi(phi),
+      mTime(t) {}
+  GPUhdi() bool operator<(const Tracklet& o) const noexcept
+  {
+    return (firstClusterIndex != o.firstClusterIndex) ? firstClusterIndex < o.firstClusterIndex : secondClusterIndex < o.secondClusterIndex;
+  }
+  GPUhdi() bool operator==(const Tracklet& o) const noexcept
   {
-    return firstClusterIndex < 0 || secondClusterIndex < 0;
+    return firstClusterIndex == o.firstClusterIndex && secondClusterIndex == o.secondClusterIndex;
   }
   GPUhdi() bool isCompatible(const Tracklet& o) const { return mTime.isCompatible(o.mTime); }
-  GPUhdi() unsigned char operator<(const Tracklet&) const;
   GPUhd() void print() const
   {
     LOGP(info, "TRKLT: fClIdx:{} sClIdx:{} ts:{}+/-{} TgL={} Phi={}", firstClusterIndex, secondClusterIndex, mTime.getTimeStamp(), mTime.getTimeStampError(), tanLambda, phi);
@@ -54,44 +62,13 @@ struct Tracklet final {
 
   int firstClusterIndex{constants::UnusedIndex};
   int secondClusterIndex{constants::UnusedIndex};
-  float tanLambda{-999};
-  float phi{-999};
+  float tanLambda{constants::UnsetValue};
+  float phi{constants::UnsetValue};
   TimeEstBC mTime;
 
   ClassDefNV(Tracklet, 1);
 };
 
-GPUhdi() Tracklet::Tracklet(const int firstClusterOrderingIndex, const int secondClusterOrderingIndex,
-                            const Cluster& firstCluster, const Cluster& secondCluster, const TimeEstBC& t)
-  : firstClusterIndex(firstClusterOrderingIndex),
-    secondClusterIndex(secondClusterOrderingIndex),
-    tanLambda((firstCluster.zCoordinate - secondCluster.zCoordinate) /
-              (firstCluster.radius - secondCluster.radius)),
-    phi(o2::gpu::GPUCommonMath::ATan2(firstCluster.yCoordinate - secondCluster.yCoordinate,
-                                      firstCluster.xCoordinate - secondCluster.xCoordinate)),
-    mTime(t)
-{
-  // Nothing to do
-}
-
-GPUhdi() Tracklet::Tracklet(const int idx0, const int idx1, float tanL, float phi, const TimeEstBC& t)
-  : firstClusterIndex(idx0),
-    secondClusterIndex(idx1),
-    tanLambda(tanL),
-    phi(phi),
-    mTime(t)
-{
-  // Nothing to do
-}
-
-GPUhdi() unsigned char Tracklet::operator<(const Tracklet& t) const
-{
-  if (isEmpty()) {
-    return false;
-  }
-  return true;
-}
-
 } // namespace o2::its
 
 #endif /* TRACKINGITS_INCLUDE_TRACKLET_H_ */