ITS: debug dump

Author: f3sch

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/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
@@ -46,12 +46,12 @@ class TrackerTraits
   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 +106,17 @@ 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);
+
  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/src/Configuration.cxx

@@ -280,6 +280,7 @@ std::vector<RecoIteration> TrackingMode::getRecoIterations(TrackingMode::Type mo
     recoIterations[0].params.MaxChi2NDF = 40.;
   }
 
+  recoIterations.resize(1);
   return recoIterations;
 }
 

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

@@ -15,25 +15,19 @@
 
 #include <oneapi/tbb/parallel_for.h>
 #include <oneapi/tbb/combinable.h>
+#include <oneapi/tbb/blocked_range.h>
+#include <oneapi/tbb/parallel_sort.h>
+
 #include <algorithm>
-#include <iostream>
 #include <iterator>
 #include <limits>
 #include <ranges>
-#include <type_traits>
-
-#ifdef OPTIMISATION_OUTPUT
-#include <format>
-#include <fstream>
-#endif
-
-#include <oneapi/tbb/blocked_range.h>
-#include <oneapi/tbb/parallel_sort.h>
 
 #include "DetectorsRaw/HBFUtils.h"
 #include "Steer/MCKinematicsReader.h"
 #include "SimulationDataFormat/O2DatabasePDG.h"
 #include "CommonConstants/MathConstants.h"
+#include "CommonUtils/TreeStreamRedirector.h"
 #include "DetectorsBase/Propagator.h"
 #include "ITSMFTBase/DPLAlpideParam.h"
 #include "GPUCommonMath.h"
@@ -46,12 +40,20 @@
 #include "ITStracking/IndexTableUtils.h"
 #include "ITStracking/Tracklet.h"
 #include "ReconstructionDataFormats/Track.h"
-#include "CommonUtils/TreeStreamRedirector.h"
 
-using o2::base::PropagatorF;
+// optimization output
+#define OPTIMISATION_TRACKLETS (1 << 0)
+#define OPTIMISATION (OPTIMISATION_TRACKLETS)
+// #define OPTIMISATION 0
+#define OPTIMISATION_ANY (OPTIMISATION != 0)
+#define OPTIMISATION_SET(flag) ((OPTIMISATION & (flag)) != 0)
 
 namespace o2::its
 {
+namespace
+{
+utils::TreeStreamRedirector* sDBGOut{nullptr};
+}
 
 struct PassMode {
   using OnePass = std::integral_constant<int, 0>;
@@ -60,13 +62,26 @@ struct PassMode {
 };
 
 template <int NLayers>
-void TrackerTraits<NLayers>::computeLayerTracklets(const int iteration, const int iSlice, const int iVertex)
+TrackerTraits<NLayers>::TrackerTraits()
 {
-#ifdef OPTIMISATION_OUTPUT
-  static int iter{0};
-  std::ofstream off(std::format("tracklets{}.txt", iter++));
+#if OPTIMISATION_ANY
+  sDBGOut = new utils::TreeStreamRedirector("its_debug.root");
 #endif
+}
 
+template <int NLayers>
+TrackerTraits<NLayers>::~TrackerTraits()
+{
+#if OPTIMISATION_ANY
+  sDBGOut->Close();
+  delete sDBGOut;
+  sDBGOut = nullptr;
+#endif
+}
+
+template <int NLayers>
+void TrackerTraits<NLayers>::computeLayerTracklets(const int iteration, const int iSlice, const int iVertex)
+{
   for (int iLayer = 0; iLayer < mRecoParams[iteration].params.TrackletsPerRoad(); ++iLayer) {
     mTimeFrame->getTracklets()[iLayer].clear();
     mTimeFrame->getTrackletsLabel(iLayer).clear();
@@ -133,8 +148,7 @@ void TrackerTraits<NLayers>::computeLayerTracklets(const int iteration, const in
           const float tanLambda = deltaZ * inverseR0;
           const float zAtRmin = tanLambda * (mTimeFrame->getMinR(iLayer + 1) - currentCluster.radius) + currentCluster.zCoordinate;
           const float zAtRmax = tanLambda * (mTimeFrame->getMaxR(iLayer + 1) - currentCluster.radius) + currentCluster.zCoordinate;
-          const float sqInvDeltaZ0 = (deltaZ2 < res2) ? res2 : 1.f / (deltaZ2 + constants::Tolerance);
-          const float sigmaZ = o2::gpu::CAMath::Sqrt((res2 + pvZ2) * math_utils::Sq(tanLambda) * ((math_utils::Sq(inverseR0) + sqInvDeltaZ0) * math_utils::Sq(meanDeltaR) + 1.f) + math_utils::Sq(meanDeltaR * mTimeFrame->getMSangle(iLayer)));
+          const float sigmaZ = o2::gpu::CAMath::Sqrt((res2 + pvZ2) * math_utils::Sq(tanLambda) * math_utils::Sq((inverseR0 * math_utils::Sq(meanDeltaR)) + 1.f) + math_utils::Sq(meanDeltaR * mTimeFrame->getMSangle(iLayer)));
           const float nSigmaZ = sigmaZ * mRecoParams[iteration].params.NSigmaCutZ;
           // phi-angle resolution
           const float sigmaPhi = o2::gpu::CAMath::Sqrt(pvR2 + math_utils::Sq(mTimeFrame->getPhiCut(iLayer)));
@@ -174,28 +188,14 @@ void TrackerTraits<NLayers>::computeLayerTracklets(const int iteration, const in
                 if (iteration && mTimeFrame->isClusterUsed(iLayer + 1, nextCluster.clusterId)) {
                   continue;
                 }
-                float deltaPhi = o2::gpu::GPUCommonMath::Abs(currentCluster.phi - nextCluster.phi);
-                float deltaZ = o2::gpu::GPUCommonMath::Abs((tanLambda * (nextCluster.radius - currentCluster.radius)) + currentCluster.zCoordinate - nextCluster.zCoordinate);
 
-#ifdef OPTIMISATION_OUTPUT
-                MCCompLabel label;
-                int currentId{currentCluster.clusterId};
-                int nextId{nextCluster.clusterId};
-                for (auto& lab1 : mTimeFrame->getClusterLabels(iLayer, currentId)) {
-                  for (auto& lab2 : mTimeFrame->getClusterLabels(iLayer + 1, nextId)) {
-                    if (lab1 == lab2 && lab1.isValid()) {
-                      label = lab1;
-                      break;
-                    }
-                  }
-                  if (label.isValid()) {
-                    break;
-                  }
-                }
-                off << std::format("{}\t{:d}\t{}\t{}\t{}\t{}", iLayer, label.isValid(), (tanLambda * (nextCluster.radius - currentCluster.radius) + currentCluster.zCoordinate - nextCluster.zCoordinate) / sigmaZ, tanLambda, resolution, sigmaZ) << std::endl;
-#endif
+                const float deltaPhi = o2::gpu::CAMath::Abs(o2::math_utils::toPMPi(currentCluster.phi - nextCluster.phi));
+                const float deltaZ = o2::gpu::CAMath::Abs((tanLambda * (nextCluster.radius - currentCluster.radius)) + currentCluster.zCoordinate - nextCluster.zCoordinate);
 
-                if (deltaZ < nSigmaZ && ((deltaPhi < nSigmaPhi || o2::gpu::GPUCommonMath::Abs(deltaPhi - o2::constants::math::TwoPI) < nSigmaPhi))) {
+                const bool accept = deltaZ < nSigmaZ && o2::gpu::GPUCommonMath::Abs(deltaPhi) < nSigmaPhi;
+                debugComputeLayerTracklets(iteration, iLayer, currentCluster, nextCluster, pv, sigmaZ, sigmaPhi, accept);
+
+                if (accept) {
                   const float phi{o2::gpu::CAMath::ATan2(currentCluster.yCoordinate - nextCluster.yCoordinate, currentCluster.xCoordinate - nextCluster.xCoordinate)};
                   const float tanL = (currentCluster.zCoordinate - nextCluster.zCoordinate) / (currentCluster.radius - nextCluster.radius);
 
@@ -346,7 +346,7 @@ void TrackerTraits<NLayers>::computeLayerCells(const int iteration)
 
         const float deltaTanLambda{std::abs(currentTracklet.tanLambda - nextTracklet.tanLambda)};
 
-#ifdef OPTIMISATION_OUTPUT
+#ifdef OPTIMISATION_OUTPUT_
         float resolution{o2::gpu::CAMath::Sqrt(0.5f * (mRecoParams[iteration].params.SystErrorZ2[iLayer] + mRecoParams[iteration].params.SystErrorZ2[iLayer + 1] + mRecoParams[iteration].params.SystErrorZ2[iLayer + 2] + mRecoParams[iteration].params.SystErrorY2[iLayer] + mRecoParams[iteration].params.SystErrorY2[iLayer + 1] + mRecoParams[iteration].params.SystErrorY2[iLayer + 2])) / mRecoParams[iteration].params.LayerResolution[iLayer]};
         resolution = resolution > 1.e-12 ? resolution : 1.f;
         bool good{mTimeFrame->getTrackletsLabel(iLayer)[iTracklet] == mTimeFrame->getTrackletsLabel(iLayer + 1)[iNextTracklet]};
@@ -786,7 +786,7 @@ void TrackerTraits<NLayers>::findCellSeeds(const int iteration)
 template <int NLayers>
 void TrackerTraits<NLayers>::findCellsNeighbours(const int iteration)
 {
-#ifdef OPTIMISATION_OUTPUT
+#ifdef OPTIMISATION_OUTPUT_
   std::ofstream off(std::format("cellneighs{}.txt", iteration));
 #endif
 
@@ -833,7 +833,7 @@ void TrackerTraits<NLayers>::findCellsNeighbours(const int iteration)
           }
           float chi2 = currentCellSeed.getPredictedChi2(nextCellSeed); /// TODO: switch to the chi2 wrt cluster to avoid correlation
 
-#ifdef OPTIMISATION_OUTPUT
+#ifdef OPTIMISATION_OUTPUT_
           bool good{mTimeFrame->getCellsLabel(iLayer)[iCell] == mTimeFrame->getCellsLabel(iLayer + 1)[iNextCell]};
           off << std::format("{}\t{:d}\t{}", iLayer, good, chi2) << std::endl;
 #endif
@@ -1437,12 +1437,12 @@ bool TrackerTraits<NLayers>::trackFollowing(TrackITSExt* track, int rof, bool ou
       }
       // estimate hypo's trk parameters at that x
       auto& hypoParam{outward ? hypo.getParamOut() : hypo.getParamIn()};
-      if (!propInstance->propagateToX(hypoParam, x, mTimeFrame->getBz(), PropagatorF::MAX_SIN_PHI,
-                                      PropagatorF::MAX_STEP, mRecoParams[iteration].params.CorrType)) {
+      if (!propInstance->propagateToX(hypoParam, x, mTimeFrame->getBz(), base::PropagatorF::MAX_SIN_PHI,
+                                      base::PropagatorF::MAX_STEP, mRecoParams[iteration].params.CorrType)) {
         continue;
       }
 
-      if (mRecoParams[iteration].params.CorrType == PropagatorF::MatCorrType::USEMatCorrNONE) { // account for material affects if propagator does not
+      if (mRecoParams[iteration].params.CorrType == base::PropagatorF::MatCorrType::USEMatCorrNONE) { // account for material affects if propagator does not
         if (!hypoParam.correctForMaterial(mRecoParams[iteration].params.LayerxX0[iLayer], mRecoParams[iteration].params.LayerxX0[iLayer] * constants::Radl * constants::Rho, true)) {
           continue;
         }
@@ -1496,7 +1496,7 @@ bool TrackerTraits<NLayers>::trackFollowing(TrackITSExt* track, int rof, bool ou
           }
 
           if (!propInstance->propagateToX(tbuParams, trackingHit.xTrackingFrame, mTimeFrame->getBz(),
-                                          PropagatorF::MAX_SIN_PHI, PropagatorF::MAX_STEP, PropagatorF::MatCorrType::USEMatCorrNONE)) {
+                                          base::PropagatorF::MAX_SIN_PHI, base::PropagatorF::MAX_STEP, base::PropagatorF::MatCorrType::USEMatCorrNONE)) {
             continue;
           }
 
@@ -1641,7 +1641,8 @@ void TrackerTraits<NLayers>::updateTrackingParameters(const std::vector<RecoIter
 template <int NLayers>
 void TrackerTraits<NLayers>::setNThreads(int n, std::shared_ptr<tbb::task_arena>& arena)
 {
-#if defined(OPTIMISATION_OUTPUT) || defined(CA_DEBUG)
+#if OPTIMISATION_ANY
+  LOGP(info, "Debug output enabled; enforcing single-threaded");
   mTaskArena = std::make_shared<tbb::task_arena>(1);
 #else
   if (arena == nullptr) {
@@ -1654,6 +1655,62 @@ void TrackerTraits<NLayers>::setNThreads(int n, std::shared_ptr<tbb::task_arena>
 #endif
 }
 
+/// Debug dumps
+template <int NLayers>
+inline void TrackerTraits<NLayers>::debugComputeLayerTracklets(int iteration, int layer, const Cluster& currentCls, const Cluster& nextCls, const Vertex& pv, float sigmaZ, float sigmaPhi, bool accepted)
+{
+#if !(OPTIMISATION_SET(OPTIMISATION_TRACKLETS))
+  return; // no-op
+#endif
+
+  static LogLogThrottler logger;
+  LOG_IF(info, logger.needToLog(iteration, layer)) << "debugTree: LayerTracklets:" << iteration << ":" << layer << " dumped entries " << logger.evCount;
+
+  static uint32_t cnt{0};
+  if ((++cnt % 50) != 0) {
+    return;
+  }
+
+  MCCompLabel label;
+  if (mTimeFrame->hasMCinformation()) {
+    int currentId{currentCls.clusterId};
+    int nextId{nextCls.clusterId};
+    for (auto& lab1 : mTimeFrame->getClusterLabels(layer, currentId)) {
+      for (auto& lab2 : mTimeFrame->getClusterLabels(layer + 1, nextId)) {
+        if (lab1 == lab2 && lab1.isValid()) {
+          label = lab1;
+          break;
+        }
+      }
+      if (label.isValid()) {
+        break;
+      }
+    }
+  }
+
+  const float deltaPvZ = currentCls.zCoordinate - pv.getZ();
+  const float tanLambda = deltaPvZ / currentCls.radius;
+  const float deltaZ = o2::gpu::GPUCommonMath::Abs((tanLambda * (nextCls.radius - currentCls.radius)) + currentCls.zCoordinate - nextCls.zCoordinate);
+  const float tglDiff = ((nextCls.zCoordinate - currentCls.zCoordinate) / (nextCls.radius - currentCls.radius)) - tanLambda;
+  const float deltaPhi = o2::gpu::CAMath::Abs(o2::math_utils::toPMPi(currentCls.phi - nextCls.phi));
+  (*sDBGOut) << "tracklets"
+             << "iter=" << iteration
+             << "lay=" << layer
+             << "lbl=" << label
+             << "pv=" << pv
+             << "curCls=" << currentCls
+             << "nextCls=" << nextCls
+             << "sigZ=" << sigmaZ
+             << "delZ=" << deltaZ
+             << "sigPhi=" << sigmaPhi
+             << "delPhi=" << deltaPhi
+             << "tgl=" << tanLambda
+             << "tglDiff=" << tglDiff
+             << "acc=" << accepted
+             << "\n";
+}
+
+// explicitly instaniate the ITS2/ITS3 tracker functions
 template class TrackerTraits<7>;
 
 } // namespace o2::its