ITS: tracking same as dev

Signed-off-by: Felix Schlepper <felix.schlepper@cern.ch>

Author: f3sch

DataFormats/Detectors/ITSMFT/ITS/include/DataFormatsITS/TimeEstBC.h

@@ -52,18 +52,6 @@ class TimeEstBC : public o2::dataformats::TimeStampWithError<TimeStampType, Time
     return res;
   }
 
- private:
-  // add the other timestmap to this one
-  // this assumes already that both overlap
-  GPUhdi() void add(const TimeEstBC& o) noexcept
-  {
-    const TimeStampType lo = o2::gpu::CAMath::Max(lower(), o.lower());
-    const TimeStampType hi = o2::gpu::CAMath::Min(upper(), o.upper());
-    const TimeStampType half = (hi - lo) / 2u;
-    this->setTimeStamp(lo + half);
-    this->setTimeStampError(static_cast<TimeStampErrorType>(half));
-  }
-
   GPUhdi() TimeStampType upper() const noexcept
   {
     TimeStampType t = this->getTimeStamp();
@@ -79,8 +67,20 @@ class TimeEstBC : public o2::dataformats::TimeStampWithError<TimeStampType, Time
     return (t > e) ? (t - e) : 0u;
   }
 
+ private:
+  // add the other timestmap to this one
+  // this assumes already that both overlap
+  GPUhdi() void add(const TimeEstBC& o) noexcept
+  {
+    const TimeStampType lo = o2::gpu::CAMath::Max(lower(), o.lower());
+    const TimeStampType hi = o2::gpu::CAMath::Min(upper(), o.upper());
+    const TimeStampType half = (hi - lo) / 2u;
+    this->setTimeStamp(lo + half);
+    this->setTimeStampError(static_cast<TimeStampErrorType>(half));
+  }
+
   ClassDefNV(TimeEstBC, 1);
 };
 } // namespace o2::its
 
-#endif
+#endif

DataFormats/Detectors/ITSMFT/ITS/include/DataFormatsITS/TrackITS.h

@@ -1,4 +1,4 @@
-// Copyright 2019-2020 CERN and copyright holders of ALICE O2.
+// Copyright 2019-2026 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.
 //
@@ -21,7 +21,7 @@
 #include "GPUCommonDef.h"
 #include "ReconstructionDataFormats/Track.h"
 #include "CommonDataFormat/RangeReference.h"
-#include "CommonDataFormat/TimeStamp.h"
+#include "DataFormatsITS/TimeEstBC.h"
 
 namespace o2
 {
@@ -41,7 +41,6 @@ class TrackITS : public o2::track::TrackParCov
 
   using Cluster = o2::itsmft::Cluster;
   using ClusRefs = o2::dataformats::RangeRefComp<4>;
-  using Timestamp = o2::dataformats::TimeStampWithError<uint32_t, uint16_t>;
 
  public:
   using o2::track::TrackParCov::TrackParCov; // inherit base constructors
@@ -160,7 +159,7 @@ class TrackITS : public o2::track::TrackParCov
   float mChi2 = 0.;                 ///< Chi2 for this track
   uint32_t mPattern = 0;            ///< layers pattern
   uint32_t mClusterSizes = 0u;      ///< 4bit packed cluster sizes
-  Timestamp mTime;                  ///< track time stamp with error in BC since start of TF, symmetrical
+  TimeEstBC mTime;                  ///< track time stamp with error in BC since start of TF, symmetrical
 
   ClassDefNV(TrackITS, 7);
 };

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

@@ -17,7 +17,6 @@
 #include <limits>
 #include <format>
 #include <string>
-#include <numeric>
 #include <vector>
 
 #include "CommonConstants/LHCConstants.h"
@@ -33,7 +32,7 @@ namespace o2::its
 
 // Layer timing definition
 struct LayerTiming {
-  using BCType = uint32_t;
+  using BCType = TimeStampType;
   BCType mNROFsTF{0};       // number of ROFs per timeframe
   BCType mROFLength{0};     // ROF length in BC
   BCType mROFDelay{0};      // delay of ROFs wrt start of first orbit in TF in BC
@@ -69,7 +68,14 @@ struct LayerTiming {
     }
     const BCType start = getROFStartInBC(rofId);
     const BCType half = mROFLength / BCType(2);
-    return {start + half, static_cast<uint16_t>(half)};
+    return {start + half, static_cast<TimeStampErrorType>(half)};
+  }
+
+  // return which ROF this BC belongs to
+  GPUhi() BCType getROF(BCType bc) const noexcept
+  {
+    BCType rof = (bc - mROFDelay - mROFBias) / mROFLength;
+    return rof >= 0 ? rof : 0;
   }
 
   GPUh() std::string asString() const
@@ -136,7 +142,7 @@ struct ROFOverlapTableView {
     return mLayers[layer];
   }
 
-  GPUh() const LayerTiming& getClockLayer() const noexcept
+  GPUh() int getClock() const noexcept
   {
     // we take the fastest layer as clock
     int fastest = 0;
@@ -145,10 +151,14 @@ struct ROFOverlapTableView {
       const auto& layer = getLayer(iL);
       if (layer.mROFLength < shortestROF) {
         fastest = iL;
-        shortestROF = layer.mROFLength;
       }
     }
-    return mLayers[fastest];
+    return fastest;
+  }
+
+  GPUh() const LayerTiming& getClockLayer() const noexcept
+  {
+    return mLayers[getClock()];
   }
 
   GPUhdi() const TableEntry& getOverlap(int32_t from, int32_t to, size_t rofIdx) const noexcept
@@ -246,7 +256,7 @@ struct ROFOverlapTableView {
     for (int32_t i = 0; i < NLayers; ++i) {
       for (int32_t j = 0; j < NLayers; ++j) {
         if (i != j) {
-          const size_t linearIdx = i * NLayers + j;
+          const size_t linearIdx = (i * NLayers) + j;
           const auto& idx = mIndices[linearIdx];
           totalEntries += idx.getEntries();
           flatTableSize += idx.getEntries();
@@ -580,8 +590,9 @@ class ROFVertexLookupTable : public LayerTimingBase<NLayers>
       while (firstVertex < lastVertex) {
         int64_t vUpper = (int64_t)vertices[firstVertex].getTimeStamp().getTimeStamp() +
                          (int64_t)vertices[firstVertex].getTimeStamp().getTimeStampError();
-        if (vUpper > rofLower)
+        if (vUpper > rofLower) {
           break;
+        }
         ++firstVertex;
       }
       size_t count = (lastVertex > firstVertex) ? (lastVertex - firstVertex) : 0;
@@ -601,8 +612,9 @@ class ROFVertexLookupTable : public LayerTimingBase<NLayers>
     while (firstVertex < lastVertex) {
       int64_t vUpper = (int64_t)vertices[firstVertex].getTimeStamp().getTimeStamp() +
                        (int64_t)vertices[firstVertex].getTimeStamp().getTimeStampError();
-      if (vUpper > rofLower)
+      if (vUpper > rofLower) {
         break;
+      }
       ++firstVertex;
     }
     size_t count = (lastVertex > firstVertex) ? (lastVertex - firstVertex) : 0;

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

@@ -73,7 +73,7 @@ void TimeFrame<NLayers>::loadROFrameData(gsl::span<const o2::itsmft::ROFRecord>
 
   // check for missing/empty/unset rofs
   // the code requires consistent monotonically increasing input without gaps
-  const auto& timing = mROFOverlapTableView.getLayer(layer);
+  const auto& timing = mROFOverlapTableView.getLayer(layer >= 0 ? layer : 0);
   if (timing.mNROFsTF != rofs.size()) {
     LOGP(fatal, "Received inconsistent number of rofs on layer:{} expected:{} received:{}", layer, timing.mNROFsTF, rofs.size());
   }
@@ -103,52 +103,82 @@ void TimeFrame<NLayers>::loadROFrameData(gsl::span<const o2::itsmft::ROFRecord>
         locXYZ = dict->getClusterCoordinates(c, patt, false);
         clusterSize = patt.getNPixels();
       }
-      mClusterSize[layer][clusterId] = std::clamp(clusterSize, 0u, 255u);
+      mClusterSize[layer >= 0 ? layer : 0][clusterId] = std::clamp(clusterSize, 0u, 255u);
       auto sensorID = c.getSensorID();
       // Inverse transformation to the local --> tracking
       auto trkXYZ = geom->getMatrixT2L(sensorID) ^ locXYZ;
       // Transformation to the local --> global
       auto gloXYZ = geom->getMatrixL2G(sensorID) * locXYZ;
-      addTrackingFrameInfoToLayer(layer, gloXYZ.x(), gloXYZ.y(), gloXYZ.z(), trkXYZ.x(), geom->getSensorRefAlpha(sensorID),
+      addTrackingFrameInfoToLayer(lay, gloXYZ.x(), gloXYZ.y(), gloXYZ.z(), trkXYZ.x(), geom->getSensorRefAlpha(sensorID),
                                   std::array<float, 2>{trkXYZ.y(), trkXYZ.z()},
                                   std::array<float, 3>{sigmaY2, sigmaYZ, sigmaZ2});
       /// Rotate to the global frame
-      addClusterToLayer(layer, gloXYZ.x(), gloXYZ.y(), gloXYZ.z(), mUnsortedClusters[layer].size());
-      addClusterExternalIndexToLayer(layer, clusterId);
+      addClusterToLayer(lay, gloXYZ.x(), gloXYZ.y(), gloXYZ.z(), mUnsortedClusters[lay].size());
+      addClusterExternalIndexToLayer(lay, clusterId);
+    }
+    // effectively calculating an exclusive sum
+    if (layer >= 0) {
+      mROFramesClusters[layer][iRof + 1] = mUnsortedClusters[layer].size();
+    } else {
+      for (unsigned int iL{0}; iL < mUnsortedClusters.size(); ++iL) {
+        mROFramesClusters[iL][iRof + 1] = mUnsortedClusters[iL].size();
+      }
     }
-    mROFramesClusters[layer][iRof + 1] = mUnsortedClusters[layer].size(); // effectively calculating an exclusive sum
   }
 
-  if (layer == 1) {
+  if (layer == 1 || layer == -1) {
     for (auto i = 0; i < mNTrackletsPerCluster.size(); ++i) {
       mNTrackletsPerCluster[i].resize(mUnsortedClusters[1].size());
       mNTrackletsPerClusterSum[i].resize(mUnsortedClusters[1].size() + 1);
     }
   }
 
   if (mcLabels != nullptr) {
-    mClusterLabels[layer] = mcLabels;
+    mClusterLabels[layer >= 0 ? layer : 0] = mcLabels;
   } else {
-    mClusterLabels[layer] = nullptr;
+    mClusterLabels[layer >= 0 ? layer : 0] = nullptr;
   }
 }
 
 template <int NLayers>
 void TimeFrame<NLayers>::resetROFrameData(int layer)
 {
-  deepVectorClear(mUnsortedClusters[layer], getMaybeFrameworkHostResource());
-  deepVectorClear(mTrackingFrameInfo[layer], getMaybeFrameworkHostResource());
-  deepVectorClear(mClusterExternalIndices[layer], mMemoryPool.get());
-  clearResizeBoundedVector(mROFramesClusters[layer], mROFOverlapTableView.getLayer(layer).mNROFsTF + 1, getMaybeFrameworkHostResource());
+  if (layer >= 0) {
+    deepVectorClear(mUnsortedClusters[layer], getMaybeFrameworkHostResource());
+    deepVectorClear(mTrackingFrameInfo[layer], getMaybeFrameworkHostResource());
+    deepVectorClear(mClusterExternalIndices[layer], mMemoryPool.get());
+    clearResizeBoundedVector(mROFramesClusters[layer], mROFOverlapTableView.getLayer(layer).mNROFsTF + 1, getMaybeFrameworkHostResource());
+  } else {
+    for (int iLayer{0}; iLayer < NLayers; ++iLayer) {
+      deepVectorClear(mUnsortedClusters[iLayer], getMaybeFrameworkHostResource());
+      deepVectorClear(mTrackingFrameInfo[iLayer], getMaybeFrameworkHostResource());
+      deepVectorClear(mClusterExternalIndices[iLayer], mMemoryPool.get());
+      clearResizeBoundedVector(mROFramesClusters[iLayer], mROFOverlapTableView.getLayer(iLayer).mNROFsTF + 1, getMaybeFrameworkHostResource());
+    }
+  }
 }
 
 template <int NLayers>
 void TimeFrame<NLayers>::prepareROFrameData(gsl::span<const itsmft::CompClusterExt> clusters, int layer)
 {
-  mUnsortedClusters[layer].reserve(clusters.size());
-  mTrackingFrameInfo[layer].reserve(clusters.size());
-  mClusterExternalIndices[layer].reserve(clusters.size());
-  clearResizeBoundedVector(mClusterSize[layer], clusters.size(), mMemoryPool.get());
+  if (layer >= 0) {
+    mUnsortedClusters[layer].reserve(clusters.size());
+    mTrackingFrameInfo[layer].reserve(clusters.size());
+    mClusterExternalIndices[layer].reserve(clusters.size());
+    clearResizeBoundedVector(mClusterSize[layer], clusters.size(), mMemoryPool.get());
+  } else {
+    auto* geom = GeometryTGeo::Instance();
+    clearResizeBoundedVector(mClusterSize[0], clusters.size(), mMemoryPool.get());
+    std::array<size_t, NLayers> clusterCountPerLayer{0};
+    for (const auto& cls : clusters) {
+      ++clusterCountPerLayer[geom->getLayer(cls.getChipID())];
+    }
+    for (int iLayer{0}; iLayer < NLayers; ++iLayer) {
+      mUnsortedClusters[iLayer].reserve(clusterCountPerLayer[iLayer]);
+      mTrackingFrameInfo[iLayer].reserve(clusterCountPerLayer[iLayer]);
+      mClusterExternalIndices[iLayer].reserve(clusterCountPerLayer[iLayer]);
+    }
+  }
 }
 
 template <int NLayers>

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

@@ -278,16 +278,15 @@ void TrackerTraits<NLayers>::computeLayerCells(const int iteration)
       int foundCells{0};
       for (int iNextTracklet{nextLayerFirstTrackletIndex}; iNextTracklet < nextLayerLastTrackletIndex; ++iNextTracklet) {
         const Tracklet& nextTracklet{mTimeFrame->getTracklets()[iLayer + 1][iNextTracklet]};
-        const auto& nextLbl = mTimeFrame->getTrackletsLabel(iLayer + 1)[iNextTracklet];
         if (mTimeFrame->getTracklets()[iLayer + 1][iNextTracklet].firstClusterIndex != nextLayerClusterIndex) {
           break;
         }
         if (!currentTracklet.getTimeStamp().isCompatible(nextTracklet.getTimeStamp())) {
           continue;
         }
-        const float deltaTanLambda{std::abs(currentTracklet.tanLambda - nextTracklet.tanLambda)};
 
-        if (deltaTanLambda / mTrkParams[iteration].CellDeltaTanLambdaSigma < mTrkParams[iteration].NSigmaCut) {
+        const float deltaTanLambdaSigma = std::abs(currentTracklet.tanLambda - nextTracklet.tanLambda) / mTrkParams[iteration].CellDeltaTanLambdaSigma;
+        if (deltaTanLambdaSigma < mTrkParams[iteration].NSigmaCut) {
 
           /// Track seed preparation. Clusters are numbered progressively from the innermost going outward.
           const int clusId[3]{

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

@@ -9,6 +9,7 @@
 // granted to it by virtue of its status as an Intergovernmental Organization
 // or submit itself to any jurisdiction.
 
+#include <algorithm>
 #include <memory>
 
 #include <oneapi/tbb/task_arena.h>
@@ -129,9 +130,9 @@ void ITSTrackingInterface::run(framework::ProcessingContext& pc)
   }
   const auto& alpParams = o2::itsmft::DPLAlpideParam<o2::detectors::DetID::ITS>::Instance();
   for (int iLayer = 0; iLayer < ((mDoStaggering) ? NLayers : 1); ++iLayer) {
-    LOGP(info, "ITSTracker:{} pulled {} clusters, {} RO frames", iLayer, compClusters[iLayer].size(), rofsinput[iLayer].size());
+    LOGP(info, "ITSTracker{} pulled {} clusters, {} RO frames", ((mDoStaggering) ? std::format(":{}", iLayer) : ""), compClusters[iLayer].size(), rofsinput[iLayer].size());
     if (compClusters[iLayer].empty()) {
-      LOGP(warn, " -> received no processable data on layer {}", iLayer);
+      LOGP(warn, " -> received no processable data{}", (mDoStaggering) ? std::format(" on layer {}", iLayer) : "");
     }
     if (mIsMC) {
       LOG(info) << " -> " << labels[iLayer]->getIndexedSize() << " MC label objects";
@@ -186,9 +187,9 @@ void ITSTrackingInterface::run(framework::ProcessingContext& pc)
 
   mTracker->setBz(o2::base::Propagator::Instance()->getNominalBz());
 
-  for (int iLayer = 0; iLayer < NLayers; ++iLayer) {
+  for (int iLayer = 0; iLayer < ((mDoStaggering) ? NLayers : 1); ++iLayer) {
     gsl::span<const unsigned char>::iterator pattIt = patterns[iLayer].begin();
-    loadROF(rofsinput[iLayer], compClusters[iLayer], pattIt, iLayer, labels[iLayer]);
+    loadROF(rofsinput[iLayer], compClusters[iLayer], pattIt, ((mDoStaggering) ? iLayer : -1), labels[iLayer]);
   }
 
   auto logger = [&](const std::string& s) { LOG(info) << s; };
@@ -299,30 +300,46 @@ void ITSTrackingInterface::run(framework::ProcessingContext& pc)
     auto& tracks = mTimeFrame->getTracks();
     allTrackLabels.reserve(mTimeFrame->getTracksLabel().size()); // should be 0 if not MC
     std::copy(mTimeFrame->getTracksLabel().begin(), mTimeFrame->getTracksLabel().end(), std::back_inserter(allTrackLabels));
-    // Some conversions that needs to be moved in the tracker internals
-    // also we create the track to clock ROF association here
-    // the clock ROF is just the fastest ROF (the number of ROFs does not necessarily reflect the actual ROFs due to
-    // possible delay of other layers)
+
+    // create the track to clock ROF association here
+    // the clock ROF is just the fastest ROF
+    // the number of ROFs does not necessarily reflect the actual ROFs
+    // due to possible delay of other layers, however it is guaranteed to be >=0
     // tracks are guaranteed to be sorted here by their lower edge
-    const auto& clockROF = mTimeFrame->getROFOverlapTableView().getClockLayer();
-    // TODO:
+    // NOTE: we are not setting the BCData of these ROFs (should we?)
+    const auto& clockLayer = mTimeFrame->getROFOverlapTableView().getClockLayer();
+    int highestROF{0};
+    for (const auto& trc : tracks) {
+      highestROF = std::max(highestROF, (int)clockLayer.getROF(trc.getTimeStamp().lower()));
+    }
+    allTrackROFs.resize(highestROF);
 
+    // Some conversions that needs to be moved in the tracker internals
+    std::vector<int> rofEntries(highestROF + 1, 0);
     for (unsigned int iTrk{0}; iTrk < tracks.size(); ++iTrk) {
       auto& trc{tracks[iTrk]};
-      trc.setFirstClusterEntry(allClusIdx.size()); // before adding tracks, create final cluster indices
+      trc.setFirstClusterEntry((int)allClusIdx.size()); // before adding tracks, create final cluster indices
       int ncl = trc.getNumberOfClusters(), nclf = 0;
       for (int ic = TrackITSExt::MaxClusters; ic--;) { // track internally keeps in->out cluster indices, but we want to store the references as out->in!!!
         auto clid = trc.getClusterIndex(ic);
         if (clid >= 0) {
-          trc.setClusterSize(ic, mTimeFrame->getClusterSize(ic, clid));
+          trc.setClusterSize(ic, mTimeFrame->getClusterSize((mDoStaggering) ? ic : 0, clid));
           allClusIdx.push_back(clid);
           nclf++;
         }
       }
       assert(ncl == nclf);
       allTracks.emplace_back(trc);
+      auto rof = clockLayer.getROF(trc.getTimeStamp().lower());
+      ++rofEntries[rof];
+    }
+    std::exclusive_scan(rofEntries.begin(), rofEntries.end(), rofEntries.begin(), 0);
+    for (size_t iROF{0}; iROF < allTrackROFs.size(); ++iROF) {
+      allTrackROFs[iROF].setFirstEntry(rofEntries[iROF]);
+      allTrackROFs[iROF].setNEntries(rofEntries[iROF + 1] - rofEntries[iROF]);
     }
   }
+
   LOGP(info, "ITSTracker pushed {} tracks in {} rofs and {} vertices", allTracks.size(), allTrackROFs.size(), vertices.size());
   if (mIsMC) {
     LOGP(info, "ITSTracker pushed {} track labels", allTrackLabels.size());