Merge a8b53cc into sapling-pr-archive-ktf

Author: ktf

DataFormats/Detectors/TPC/include/DataFormatsTPC/ClusterNative.h

@@ -15,6 +15,7 @@
 #ifndef ALICEO2_DATAFORMATSTPC_CLUSTERNATIVE_H
 #define ALICEO2_DATAFORMATSTPC_CLUSTERNATIVE_H
 #ifndef GPUCA_GPUCODE_DEVICE
+#include <climits>
 #include <cstdint>
 #include <cstddef> // for size_t
 #include <utility>
@@ -62,6 +63,9 @@ struct ClusterNative {
   static constexpr int scalePadPacked = 64;       //< ~60 is needed for 0.1mm precision, but power of two avoids rounding
   static constexpr int scaleSigmaTimePacked = 32; // 1/32nd of pad/timebin precision for cluster size
   static constexpr int scaleSigmaPadPacked = 32;
+  static constexpr int scaleSaturatedQtot = 8;
+  static constexpr int maxRegularQtot = 25 * 1024;
+  static constexpr int maxSaturatedQtot = (USHRT_MAX - maxRegularQtot) * scaleSaturatedQtot;
 
   uint32_t timeFlagsPacked; //< Contains the time in the lower 24 bits in a packed format, contains the flags in the
                             // upper 8 bits
@@ -83,7 +87,14 @@ struct ClusterNative {
   }
 
   GPUd() uint16_t getQmax() const { return qMax; }
-  GPUd() uint16_t getQtot() const { return qTot; }
+  GPUd() uint16_t getQtot() const
+  {
+    if (isSaturated()) [[unlikely]] {
+      auto sQtot = getSaturatedQtot();
+      return sQtot < USHRT_MAX ? sQtot : USHRT_MAX;
+    }
+    return qTot;
+  }
   GPUd() uint8_t getFlags() const { return timeFlagsPacked >> 24; }
   GPUd() uint32_t getTimePacked() const { return timeFlagsPacked & 0xFFFFFF; }
   GPUd() void setTimePackedFlags(uint32_t timePacked, uint8_t flags)
@@ -119,7 +130,13 @@ struct ClusterNative {
   /// Y = (12.4 - 0.5 * (66 - 1)) * 4.16mm = -83.616mm
   GPUd() float getPad() const { return unpackPad(padPacked); }
   GPUd() void setPad(float pad) { padPacked = packPad(pad); }
-  GPUd() float getSigmaTime() const { return float(sigmaTimePacked) * (1.f / scaleSigmaTimePacked); }
+  GPUd() float getSigmaTime() const
+  {
+    if (isSaturated()) [[unlikely]] {
+      return 0;
+    }
+    return float(sigmaTimePacked) * (1.f / scaleSigmaTimePacked);
+  }
   GPUd() void setSigmaTime(float sigmaTime)
   {
     uint32_t tmp = sigmaTime * scaleSigmaTimePacked + 0.5;
@@ -138,6 +155,31 @@ struct ClusterNative {
     sigmaPadPacked = tmp;
   }
 
+  GPUd() bool isSaturated() const { return qTot > maxRegularQtot; }
+
+  GPUd() void setSaturatedQtot(uint32_t qtot)
+  {
+    this->qTot = USHRT_MAX;
+    if (qtot < maxSaturatedQtot) {
+      this->qTot = ((qtot + scaleSaturatedQtot / 2) / scaleSaturatedQtot) + maxRegularQtot;
+    }
+  }
+
+  GPUd() uint32_t getSaturatedQtot() const
+  {
+    return uint32_t(qTot - maxRegularQtot) * scaleSaturatedQtot;
+  }
+
+  GPUd() void setSaturatedTailLength(uint32_t tail)
+  {
+    sigmaTimePacked = encodeTailLength(tail);
+  }
+
+  GPUd() uint32_t getSaturatedTailLength() const
+  {
+    return decodeTailLength(sigmaTimePacked);
+  }
+
   GPUd() bool operator<(const ClusterNative& rhs) const
   {
     if (this->getTimePacked() != rhs.getTimePacked()) {
@@ -167,6 +209,93 @@ struct ClusterNative {
            this->qTot == rhs.qTot &&
            this->getFlags() == rhs.getFlags();
   }
+
+ private:
+  static constexpr GPUd() uint32_t decodeTailLength(uint8_t code)
+  {
+    // Quantize tail length into 8bits.
+    // Max expected length is 1500 tbs.
+    // But allow outliers up to 8000 tbs.
+    //
+    // Full code layout is:
+    //
+    // | Code range | Decoded values |  Step | Codes |
+    // | ---------: | -------------: | ----: | ----: |
+    // |    `0..63` |        `0..63` |   `1` |  `64` |
+    // |   `64..95` |      `64..126` |   `2` |  `32` |
+    // |  `96..127` |     `128..252` |   `4` |  `32` |
+    // | `128..159` |     `256..504` |   `8` |  `32` |
+    // | `160..223` |    `512..1520` |  `16` |  `64` |
+    // | `224..239` |   `1552..2032` |  `32` |  `16` |
+    // | `240..255` |   `2048..8048` | `400` |  `16` |
+    //
+
+    if (code < 64) {
+      return code;
+    }
+
+    if (code < 160) {
+      uint32_t q = (uint32_t)code - 64u;
+      uint32_t exponent = (q >> 5) + 1u; // 1, 2, 3
+      uint32_t mantissa = q & 31u;       // 0..31
+
+      return (32u + mantissa) << exponent;
+    }
+
+    if (code < 224) {
+      return 512u + 16u * ((uint32_t)code - 160u);
+    }
+
+    if (code < 240) {
+      return 1552u + 32u * ((uint32_t)code - 224u);
+    }
+
+    return 2048u + 400u * ((uint32_t)code - 240u);
+  }
+
+  static constexpr GPUd() uint8_t encodeTailLength(uint32_t value)
+  {
+    // Saturate above representable range.
+    if (value >= decodeTailLength(255)) [[unlikely]] {
+      return 255;
+    }
+
+    // Binary search for the first code whose decoded value >= value.
+    uint8_t lo = 0;
+    uint8_t hi = 255;
+
+    while (lo < hi) {
+      uint8_t mid = lo + ((hi - lo) >> 1);
+      uint32_t decoded = decodeTailLength(mid);
+
+      if (decoded < value) {
+        lo = mid + 1;
+      } else {
+        hi = mid;
+      }
+    }
+
+    // lo is now the first code with decoded >= value.
+    if (lo == 0) [[unlikely]] {
+      return 0;
+    }
+
+    uint8_t above_code = lo;
+    uint8_t below_code = lo - 1;
+
+    uint32_t above_value = decodeTailLength(above_code);
+    uint32_t below_value = decodeTailLength(below_code);
+
+    uint32_t above_error = above_value - value;
+    uint32_t below_error = value - below_value;
+
+    // Tie-break downward.
+    if (below_error <= above_error) {
+      return below_code;
+    } else {
+      return above_code;
+    }
+  }
 };
 
 // This is an index struct to access TPC clusters inside sectors and rows. It shall not own the data, but just point to

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

@@ -18,6 +18,7 @@
 
 #include <array>
 #include "ITStracking/Constants.h"
+#include "ITStracking/MathUtils.h"
 #include "GPUCommonRtypes.h"
 #include "GPUCommonDef.h"
 
@@ -29,26 +30,21 @@ class IndexTableUtils;
 
 struct Cluster final {
   GPUhdDefault() Cluster() = default;
-  GPUhd() Cluster(const float x, const float y, const float z, const int idx);
-  template <int nLayers>
-  GPUhd() Cluster(const int, const IndexTableUtils<nLayers>& utils, const Cluster&);
-  template <int nLayers>
-  GPUhd() Cluster(const int, const float3&, const IndexTableUtils<nLayers>& utils, const Cluster&);
-  GPUhdDefault() Cluster(const Cluster&) = default;
-  GPUhdDefault() Cluster(Cluster&&) noexcept = default;
-  GPUhdDefault() ~Cluster() = default;
-
-  GPUhdDefault() Cluster& operator=(const Cluster&) = default;
-  GPUhdDefault() Cluster& operator=(Cluster&&) noexcept = default;
-  GPUhdDefault() bool operator==(const Cluster&) const = default;
-
+  GPUhd() Cluster(const float x, const float y, const float z, const int idx)
+    : xCoordinate(x),
+      yCoordinate(y),
+      zCoordinate(z),
+      phi(math_utils::computeNormalizedPhi(x, y)),
+      radius(math_utils::hypot(x, y)),
+      clusterId(idx),
+      indexTableBinIndex(0) {}
   GPUhd() void print() const;
 
-  float xCoordinate{-999.f};
-  float yCoordinate{-999.f};
-  float zCoordinate{-999.f};
-  float phi{-999.f};
-  float radius{-999.f};
+  float xCoordinate{constants::UnsetValue};
+  float yCoordinate{constants::UnsetValue};
+  float zCoordinate{constants::UnsetValue};
+  float phi{constants::UnsetValue};
+  float radius{constants::UnsetValue};
   int clusterId{constants::UnusedIndex};
   int indexTableBinIndex{constants::UnusedIndex};
 
@@ -57,23 +53,18 @@ struct Cluster final {
 
 struct TrackingFrameInfo final {
   GPUhdDefault() TrackingFrameInfo() = default;
-  GPUhd() TrackingFrameInfo(float x, float y, float z, float xTF, float alpha, std::array<float, 2>&& posTF, std::array<float, 3>&& covTF);
-  GPUhdDefault() TrackingFrameInfo(const TrackingFrameInfo&) = default;
-  GPUhdDefault() TrackingFrameInfo(TrackingFrameInfo&&) noexcept = default;
-  GPUhdDefault() ~TrackingFrameInfo() = default;
-
-  GPUhdDefault() TrackingFrameInfo& operator=(const TrackingFrameInfo&) = default;
-  GPUhdDefault() TrackingFrameInfo& operator=(TrackingFrameInfo&&) = default;
+  GPUhd() TrackingFrameInfo(float x, float y, float z, float xTF, float alpha, const std::array<float, 2>& posTF, const std::array<float, 3>& covTF)
+    : xCoordinate(x), yCoordinate(y), zCoordinate(z), xTrackingFrame(xTF), alphaTrackingFrame(alpha), positionTrackingFrame(posTF), covarianceTrackingFrame(covTF) {}
 
   GPUhd() void print() const;
 
-  float xCoordinate{-999.f};
-  float yCoordinate{-999.f};
-  float zCoordinate{-999.f};
-  float xTrackingFrame{-999.f};
-  float alphaTrackingFrame{-999.f};
-  std::array<float, 2> positionTrackingFrame = {-999.f, -999.f};
-  std::array<float, 3> covarianceTrackingFrame = {-999.f, -999.f, -999.f};
+  float xCoordinate{constants::UnsetValue};
+  float yCoordinate{constants::UnsetValue};
+  float zCoordinate{constants::UnsetValue};
+  float xTrackingFrame{constants::UnsetValue};
+  float alphaTrackingFrame{constants::UnsetValue};
+  std::array<float, 2> positionTrackingFrame{constants::UnsetValue, constants::UnsetValue};
+  std::array<float, 3> covarianceTrackingFrame{constants::UnsetValue, constants::UnsetValue, constants::UnsetValue};
 
   ClassDefNV(TrackingFrameInfo, 1);
 };

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

@@ -41,6 +41,11 @@ GPUhdi() constexpr float getNormalizedPhi(float phi)
   return phi;
 }
 
+GPUhdi() constexpr float computeNormalizedPhi(float x, float y)
+{
+  return getNormalizedPhi(computePhi(x, y));
+}
+
 GPUhdi() float computeCurvature(float x1, float y1, float x2, float y2, float x3, float y3)
 {
   // in case the triangle is degenerate we return infinite curvature.

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

@@ -647,15 +647,14 @@ class ROFVertexLookupTable : public LayerTimingBase<NLayers>
     mFlatTable[flatTableIdx].setEntries(static_cast<T>(count));
   }
 
-  // Binary search for first vertex where maxBC >= targetBC
+  // Binary search for first vertex where lowerBC >= targetBC
   GPUh() size_t binarySearchFirst(const Vertex* vertices, size_t nVertices, size_t searchStart, BCType targetBC) const
   {
     size_t left = searchStart;
     size_t right = nVertices;
     while (left < right) {
       size_t mid = left + ((right - left) / 2);
-      int64_t lower = (int64_t)vertices[mid].getTimeStamp().getTimeStamp() -
-                      (int64_t)vertices[mid].getTimeStamp().getTimeStampError();
+      int64_t lower = (int64_t)vertices[mid].getTimeStamp().lower();
       if (lower < targetBC) {
         left = mid + 1;
       } else {

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

@@ -8,77 +8,16 @@
 // In applying this license CERN does not waive the privileges and immunities
 // granted to it by virtue of its status as an Intergovernmental Organization
 // or submit itself to any jurisdiction.
-///
-/// \file Cluster.cxx
-/// \brief
-///
-#include "GPUCommonMath.h"
-#include "GPUCommonArray.h"
 
 #include "ITStracking/Cluster.h"
-#include "ITStracking/Definitions.h"
-#include "ITStracking/MathUtils.h"
-#include "ITStracking/IndexTableUtils.h"
 
 using namespace o2::its;
 
-using math_utils::computePhi;
-using math_utils::getNormalizedPhi;
-
-Cluster::Cluster(const float x, const float y, const float z, const int index)
-  : xCoordinate{x},
-    yCoordinate{y},
-    zCoordinate{z},
-    phi{getNormalizedPhi(computePhi(x, y))},
-    radius{o2::gpu::GPUCommonMath::Hypot(x, y)},
-    clusterId{index},
-    indexTableBinIndex{0}
-{
-  // Nothing to do
-}
-
-template <int nLayers>
-Cluster::Cluster(const int layerIndex, const IndexTableUtils<nLayers>& utils, const Cluster& other)
-  : xCoordinate{other.xCoordinate},
-    yCoordinate{other.yCoordinate},
-    zCoordinate{other.zCoordinate},
-    phi{getNormalizedPhi(computePhi(other.xCoordinate, other.yCoordinate))},
-    radius{o2::gpu::GPUCommonMath::Hypot(other.xCoordinate, other.yCoordinate)},
-    clusterId{other.clusterId},
-    indexTableBinIndex{utils.getBinIndex(utils.getZBinIndex(layerIndex, zCoordinate),
-                                         utils.getPhiBinIndex(phi))}
-//, montecarloId{ other.montecarloId }
-{
-  // Nothing to do
-}
-
-template <int nLayers>
-Cluster::Cluster(const int layerIndex, const float3& primaryVertex, const IndexTableUtils<nLayers>& utils, const Cluster& other)
-  : xCoordinate{other.xCoordinate},
-    yCoordinate{other.yCoordinate},
-    zCoordinate{other.zCoordinate},
-    phi{getNormalizedPhi(
-      computePhi(xCoordinate - primaryVertex.x, yCoordinate - primaryVertex.y))},
-    radius{o2::gpu::GPUCommonMath::Hypot(xCoordinate - primaryVertex.x, yCoordinate - primaryVertex.y)},
-    clusterId{other.clusterId},
-    indexTableBinIndex{utils.getBinIndex(utils.getZBinIndex(layerIndex, zCoordinate),
-                                         utils.getPhiBinIndex(phi))}
-{
-  // Nothing to do
-}
-
 GPUhd() void Cluster::print() const
 {
   printf("Cluster: %f %f %f %f %f %d %d\n", xCoordinate, yCoordinate, zCoordinate, phi, radius, clusterId, indexTableBinIndex);
 }
 
-TrackingFrameInfo::TrackingFrameInfo(float x, float y, float z, float xTF, float alpha, std::array<float, 2>&& posTF,
-                                     std::array<float, 3>&& covTF)
-  : xCoordinate{x}, yCoordinate{y}, zCoordinate{z}, xTrackingFrame{xTF}, alphaTrackingFrame{alpha}, positionTrackingFrame{posTF}, covarianceTrackingFrame{covTF}
-{
-  // Nothing to do
-}
-
 GPUhd() void TrackingFrameInfo::print() const
 {
   printf("x: %f y: %f z: %f xTF: %f alphaTF: %f posTF: %f %f covTF: %f %f %f\n",

Detectors/ITSMFT/common/workflow/src/EntropyDecoderSpec.cxx

@@ -71,26 +71,28 @@ void EntropyDecoderSpec<N>::run(ProcessingContext& pc)
     auto buff = pc.inputs().get<gsl::span<o2::ctf::BufferType>>(getBinding(nm + "CTF", iLayer));
     // since the buff is const, we cannot use EncodedBlocks::relocate directly, instead we wrap its data to another flat object
     // const auto ctfImage = o2::itsmft::CTF::getImage(buff.data());
-    const auto& ctf = o2::itsmft::CTF::getImage(buff.data());
-    if (ctf.getHeader().maxStreams != nLayers) {
-      LOGP(fatal, "Number of streams {} in the CTF header is not equal to NLayers {} from AlpideParam in {}staggered mode",
-           ctf.getHeader().maxStreams, nLayers, mDoStaggering ? "" : "non-");
+    if (buff.size()) {
+      const auto& ctf = o2::itsmft::CTF::getImage(buff.data());
+      if (ctf.getHeader().maxStreams != nLayers) {
+        LOGP(fatal, "Number of streams {} in the CTF header is not equal to NLayers {} from AlpideParam in {}staggered mode",
+             ctf.getHeader().maxStreams, nLayers, mDoStaggering ? "" : "non-");
+      }
     }
     // this produces weird memory problems in unrelated devices, to be understood
     // auto& trigs = pc.outputs().make<std::vector<o2::itsmft::PhysTrigger>>(OutputRef{"phystrig"}); // dummy output
     auto& rofs = pc.outputs().make<std::vector<o2::itsmft::ROFRecord>>(OutputRef{nm + "ROframes", iLayer});
     if (mGetDigits) {
       auto& digits = pc.outputs().make<std::vector<o2::itsmft::Digit>>(OutputRef{nm + "Digits", iLayer});
       if (buff.size()) {
-        iosize += mCTFCoder.decode(ctf, rofs, digits, mNoiseMap, mPattIdConverter);
+        iosize += mCTFCoder.decode(o2::itsmft::CTF::getImage(buff.data()), rofs, digits, mNoiseMap, mPattIdConverter);
       }
       ndigcl += digits.size();
       nrofs += rofs.size();
     } else {
       auto& compcl = pc.outputs().make<std::vector<o2::itsmft::CompClusterExt>>(OutputRef{nm + "compClusters", iLayer});
       auto& patterns = pc.outputs().make<std::vector<unsigned char>>(OutputRef{nm + "patterns", iLayer});
       if (buff.size()) {
-        iosize += mCTFCoder.decode(ctf, rofs, compcl, patterns, mNoiseMap, mPattIdConverter);
+        iosize += mCTFCoder.decode(o2::itsmft::CTF::getImage(buff.data()), rofs, compcl, patterns, mNoiseMap, mPattIdConverter);
       }
       ndigcl += compcl.size();
     }