ITS: fix time assignments

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

Author: f3sch

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

@@ -23,64 +23,69 @@ namespace o2::its
 {
 // Time estimates are given in BC
 // error needs to cover maximum 1 orbit
-// this is an symmetric time error [t0-tE, t0+tE)
 using TimeStampType = uint32_t;
 using TimeStampErrorType = uint16_t;
+// this is an symmetric time error [t0-tE, t0+tE]
+using TimeStamp = o2::dataformats::TimeStampWithError<float, float>;
+// this is an asymmetric time interval [t0, t0+tE] used for internal calculations
 class TimeEstBC : public o2::dataformats::TimeStampWithError<TimeStampType, TimeStampErrorType>
 {
- public:
   using Base = o2::dataformats::TimeStampWithError<TimeStampType, TimeStampErrorType>;
+
+ public:
   GPUhdDefault() TimeEstBC() = default;
   GPUhdi() TimeEstBC(TimeStampType t, TimeStampErrorType e) : Base(t, e) {}
-
+  // convert to symmetric center+-half representation
+  GPUhdi() its::TimeStamp makeSymmetrical() const noexcept
+  {
+    const auto start = static_cast<float>(this->getTimeStamp());
+    const float half = (float)this->getTimeStampError() / 2.f;
+    return {start + half, half};
+  }
   // check if timestamps overlap within their interval
   GPUhdi() bool isCompatible(const TimeEstBC& o) const noexcept
   {
-    return !(upper() <= o.lower() || o.upper() <= lower());
+    return this->upper() > o.lower() && o.upper() > this->lower();
   }
-
   GPUhdi() TimeEstBC& operator+=(const TimeEstBC& o) noexcept
   {
     add(o);
     return *this;
   }
-
   GPUhdi() TimeEstBC operator+(const TimeEstBC& o) const noexcept
   {
     TimeEstBC res = *this;
     res += o;
     return res;
   }
-
+  // upper bound of interval t0+tE
   GPUhdi() TimeStampType upper() const noexcept
   {
     TimeStampType t = this->getTimeStamp();
     TimeStampType e = this->getTimeStampError();
     constexpr TimeStampType max = std::numeric_limits<TimeStampType>::max();
     return (t > (max - e)) ? max : t + e;
   }
-
+  // lower bound of interval t0
   GPUhdi() TimeStampType lower() const noexcept
   {
-    TimeStampType t = this->getTimeStamp();
-    TimeStampType e = this->getTimeStampError();
-    return (t > e) ? (t - e) : 0u;
+    return this->getTimeStamp();
   }
 
  private:
-  // add the other timestmap to this one
+  // intersect with the other timestamp
   // 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));
+    const TimeStampType lo = o2::gpu::CAMath::Max(this->lower(), o.lower());
+    const TimeStampType hi = o2::gpu::CAMath::Min(this->upper(), o.upper());
+    this->setTimeStamp(lo);
+    this->setTimeStampError(static_cast<TimeStampErrorType>(hi - lo));
   }
 
   ClassDefNV(TimeEstBC, 1);
 };
+
 } // namespace o2::its
 
 #endif

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

@@ -159,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
-  TimeEstBC mTime;                  ///< track time stamp with error in BC since start of TF, symmetrical
+  TimeStamp mTime;                  ///< track time stamp with error in BC since start of TF, symmetrical
 
   ClassDefNV(TrackITS, 7);
 };

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

@@ -22,6 +22,7 @@
 
 namespace o2::its
 {
+// NOTE: this uses the internal asymmetrical time reprenstation!
 using Vertex = o2::dataformats::Vertex<o2::its::TimeEstBC>;
 using VertexLabel = std::pair<o2::MCCompLabel, float>;
 } // namespace o2::its

Detectors/ITSMFT/ITS/tracking/GPU/cuda/TrackerTraitsGPU.cxx

@@ -399,7 +399,7 @@ void TrackerTraitsGPU<NLayers>::findRoads(const int iteration)
           ts += rofTS;
         }
       }
-      track.getTimeStamp() = ts;
+      track.getTimeStamp() = ts.makeSymmetrical();
       track.setUserField(0);
       track.getParamOut().setUserField(0);
       mTimeFrameGPU->getTracks().emplace_back(track);

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

@@ -48,20 +48,12 @@ struct Line final {
 
   float originPoint[3] = {0, 0, 0};
   float cosinesDirector[3] = {0, 0, 0};
-  // float weightMatrix[6] = {1., 0., 0., 1., 0., 1.};
-  // weightMatrix is a symmetric matrix internally stored as
-  //    0 --> row = 0, col = 0
-  //    1 --> 0,1
-  //    2 --> 0,2
-  //    3 --> 1,1
-  //    4 --> 1,2
-  //    5 --> 2,2
   TimeEstBC mTime;
 
   ClassDefNV(Line, 1);
 };
 
-GPUhdi() Line::Line(const Tracklet& tracklet, const Cluster* innerClusters, const Cluster* outerClusters)
+GPUhdi() Line::Line(const Tracklet& tracklet, const Cluster* innerClusters, const Cluster* outerClusters) : mTime(tracklet.mTime)
 {
   originPoint[0] = innerClusters[tracklet.firstClusterIndex].xCoordinate;
   originPoint[1] = innerClusters[tracklet.firstClusterIndex].yCoordinate;
@@ -75,8 +67,6 @@ GPUhdi() Line::Line(const Tracklet& tracklet, const Cluster* innerClusters, cons
   cosinesDirector[0] *= inverseNorm;
   cosinesDirector[1] *= inverseNorm;
   cosinesDirector[2] *= inverseNorm;
-
-  mTime = tracklet.mTime;
 }
 
 // static functions:
@@ -140,9 +130,9 @@ GPUhdi() void Line::getDCAComponents(const Line& line, const float point[3], flo
   destArray[0] = line.originPoint[0] - point[0] + line.cosinesDirector[0] * cdelta;
   destArray[3] = line.originPoint[1] - point[1] + line.cosinesDirector[1] * cdelta;
   destArray[5] = line.originPoint[2] - point[2] + line.cosinesDirector[2] * cdelta;
-  destArray[1] = o2::gpu::CAMath::Sqrt(destArray[0] * destArray[0] + destArray[3] * destArray[3]);
-  destArray[2] = o2::gpu::CAMath::Sqrt(destArray[0] * destArray[0] + destArray[5] * destArray[5]);
-  destArray[4] = o2::gpu::CAMath::Sqrt(destArray[3] * destArray[3] + destArray[5] * destArray[5]);
+  destArray[1] = o2::gpu::CAMath::Hypot(destArray[0], destArray[3]);
+  destArray[2] = o2::gpu::CAMath::Hypot(destArray[0], destArray[5]);
+  destArray[4] = o2::gpu::CAMath::Hypot(destArray[3], destArray[5]);
 }
 
 inline bool Line::operator==(const Line& rhs) const
@@ -163,28 +153,27 @@ class ClusterLines final
 {
  public:
   ClusterLines() = default;
-  ClusterLines(const int firstLabel, const Line& firstLine, const int secondLabel, const Line& secondLine,
-               const bool weight = false);
-  ClusterLines(const Line& firstLine, const Line& secondLine);
-  void add(const int lineLabel, const Line& line, const bool weight = false);
+  ClusterLines(const int firstLabel, const Line& firstLine, const int secondLabel, const Line& secondLine);
+  void add(const int lineLabel, const Line& line);
   void computeClusterCentroid();
-  inline std::array<float, 3> getVertex() const { return mVertex; }
-  inline std::array<float, 6> getRMS2() const { return mRMS2; }
-  inline float getAvgDistance2() const { return mAvgDistance2; }
-  inline auto getSize() const noexcept { return mLabels.size(); }
+  auto const& getVertex() const { return mVertex; }
+  std::array<float, 6> getRMS2() const { return mRMS2; }
+  float getAvgDistance2() const { return mAvgDistance2; }
+  auto getSize() const noexcept { return mLabels.size(); }
   auto& getLabels() noexcept { return mLabels; }
   const auto& getTimeStamp() const noexcept { return mTime; }
   bool operator==(const ClusterLines& rhs) const noexcept;
+  float getR2() const noexcept { return (mVertex[0] * mVertex[0]) + (mVertex[1] * mVertex[1]); }
+  float getR() const noexcept { return std::sqrt(getR2()); }
 
  protected:
   std::array<float, 6> mAMatrix = {}; // AX=B
   std::array<float, 3> mBMatrix = {}; // AX=B
-  // std::array<float, 9> mWeightMatrix = {}; // weight matrix
-  std::array<float, 3> mVertex = {}; // cluster centroid position
-  std::array<float, 6> mRMS2 = {};   // symmetric matrix: diagonal is RMS2
-  float mAvgDistance2 = 0.f;         // substitute for chi2
-  TimeEstBC mTime;                   // time stamp
-  std::vector<int> mLabels;          // contributing labels
+  std::array<float, 3> mVertex = {};  // cluster centroid position
+  std::array<float, 6> mRMS2 = {};    // symmetric matrix: diagonal is RMS2
+  float mAvgDistance2 = 0.f;          // substitute for chi2
+  TimeEstBC mTime;                    // time stamp
+  std::vector<int> mLabels;           // contributing labels
 };
 
 } // namespace o2::its

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

@@ -19,13 +19,13 @@
 #include <vector>
 #ifndef GPUCA_GPUCODE
 #include <format>
-#endif
 
+#include "Framework/Logger.h"
+#endif
 #include "CommonConstants/LHCConstants.h"
 #include "CommonDataFormat/RangeReference.h"
 #include "DataFormatsITS/TimeEstBC.h"
 #include "DataFormatsITS/Vertex.h"
-#include "GPUCommonLogger.h"
 #include "GPUCommonMath.h"
 #include "GPUCommonDef.h"
 
@@ -65,12 +65,9 @@ struct LayerTiming {
       int64_t end = getROFEndInBC(rofId);
       start = o2::gpu::CAMath::Max(start - mROFAddTimeErr, int64_t(0));
       end += mROFAddTimeErr;
-      const BCType half = (end - start + 1u) / 2u;
-      return {BCType(start) + half, static_cast<uint16_t>(half)};
+      return {static_cast<BCType>(start), static_cast<TimeStampErrorType>(end - start)};
     }
-    const BCType start = getROFStartInBC(rofId);
-    const BCType half = mROFLength / BCType(2);
-    return {start + half, static_cast<TimeStampErrorType>(half)};
+    return {getROFStartInBC(rofId), static_cast<TimeStampErrorType>(mROFLength)};
   }
 
   // return which ROF this BC belongs to
@@ -83,6 +80,17 @@ struct LayerTiming {
     return (bc - offset) / mROFLength;
   }
 
+  // return which ROF this timestamp belongs by its lower edge
+  GPUhi() BCType getROF(TimeStamp ts) const noexcept
+  {
+    const BCType offset = mROFDelay + mROFBias;
+    const BCType bc = (ts.getTimeStamp() < ts.getTimeStampError()) ? BCType(0) : static_cast<BCType>(o2::gpu::CAMath::Floor(ts.getTimeStamp() - ts.getTimeStampError()));
+    if (bc <= offset) {
+      return 0;
+    }
+    return (bc - offset) / mROFLength;
+  }
+
 #ifndef GPUCA_GPUCODE
   GPUh() std::string asString() const
   {
@@ -278,7 +286,7 @@ struct ROFOverlapTableView {
       mLayers[i].print();
     }
 
-    const uint32_t totalBytes = (flatTableSize * sizeof(TableEntry)) + (NLayers * NLayers * sizeof(TableIndex));
+    const uint32_t totalBytes = (flatTableSize * sizeof(TableEntry)) + (static_cast<unsigned long>(NLayers * NLayers) * sizeof(TableIndex));
     LOGF(info, "------------------------------------------------------------");
     LOGF(info, "Total overlap table size: %u entries", totalEntries);
     LOGF(info, "Flat table size: %zu entries", flatTableSize);
@@ -332,7 +340,7 @@ class ROFOverlapTable : public LayerTimingBase<NLayers>
   }
 
   GPUh() size_t getFlatTableSize() const noexcept { return mFlatTable.size(); }
-  static GPUh() constexpr size_t getIndicesSize() { return NLayers * NLayers; }
+  static GPUh() constexpr size_t getIndicesSize() { return static_cast<size_t>(NLayers * NLayers); }
 
  private:
   GPUh() void buildMapping(int32_t from, int32_t to, std::vector<TableEntry>& table)
@@ -346,7 +354,7 @@ class ROFOverlapTable : public LayerTimingBase<NLayers>
       int64_t fromEnd = (int64_t)layerFrom.getROFEndInBC(iROF) + layerFrom.mROFAddTimeErr;
 
       int32_t firstROFTo = o2::gpu::CAMath::Max(0, (int32_t)((fromStart - (int64_t)layerTo.mROFAddTimeErr - (int64_t)layerTo.mROFDelay - (int64_t)layerTo.mROFBias) / (int64_t)layerTo.mROFLength));
-      int32_t lastROFTo = (int32_t)((fromEnd + (int64_t)layerTo.mROFAddTimeErr - (int64_t)layerTo.mROFDelay - (int64_t)layerTo.mROFBias - 1) / (int64_t)layerTo.mROFLength);
+      auto lastROFTo = (int32_t)((fromEnd + (int64_t)layerTo.mROFAddTimeErr - (int64_t)layerTo.mROFDelay - (int64_t)layerTo.mROFBias - 1) / (int64_t)layerTo.mROFLength);
       firstROFTo = o2::gpu::CAMath::Max(0, firstROFTo);
       lastROFTo = o2::gpu::CAMath::Min((int32_t)layerTo.mNROFsTF - 1, lastROFTo);
 
@@ -443,8 +451,8 @@ struct ROFVertexLookupTableView {
     const auto& layerDef = mLayers[layer];
     int64_t rofLower = o2::gpu::CAMath::Max((int64_t)layerDef.getROFStartInBC(rofIdx) - (int64_t)layerDef.mROFAddTimeErr, int64_t(0));
     int64_t rofUpper = (int64_t)layerDef.getROFEndInBC(rofIdx) + layerDef.mROFAddTimeErr;
-    int64_t vLower = (int64_t)vertex.getTimeStamp().getTimeStamp() - (int64_t)vertex.getTimeStamp().getTimeStampError();
-    int64_t vUpper = (int64_t)vertex.getTimeStamp().getTimeStamp() + (int64_t)vertex.getTimeStamp().getTimeStampError();
+    auto vLower = (int64_t)vertex.getTimeStamp().lower();
+    auto vUpper = (int64_t)vertex.getTimeStamp().upper();
     return vUpper >= rofLower && vLower < rofUpper;
   }
 
@@ -601,8 +609,7 @@ class ROFVertexLookupTable : public LayerTimingBase<NLayers>
       size_t lastVertex = binarySearchFirst(vertices, nVertices, vertexSearchStart, rofUpper);
       size_t firstVertex = vertexSearchStart;
       while (firstVertex < lastVertex) {
-        int64_t vUpper = (int64_t)vertices[firstVertex].getTimeStamp().getTimeStamp() +
-                         (int64_t)vertices[firstVertex].getTimeStamp().getTimeStampError();
+        auto vUpper = (int64_t)vertices[firstVertex].getTimeStamp().upper();
         if (vUpper > rofLower) {
           break;
         }