TPC SCD: round instead of truncate when packing the unbinned residuals (#16046)

Co-authored-by: Claude Opus 5.5 <noreply@anthropic.com>

Author: matthias-kleiner

Detectors/TPC/calibration/SpacePoints/include/SpacePoints/TrackInterpolation.h

@@ -18,6 +18,8 @@
 #define ALICEO2_TPC_TRACKINTERPOLATION_H_
 
 #include <gsl/span>
+#include <algorithm>
+#include <cmath>
 #include "CommonDataFormat/EvIndex.h"
 #include "CommonDataFormat/RangeReference.h"
 #include "ReconstructionDataFormats/Track.h"
@@ -77,11 +79,11 @@ struct TPCClusterResiduals {
 struct UnbinnedResid {
   UnbinnedResid() = default;
   UnbinnedResid(float dyIn, float dzIn, float tgSlpIn, float yIn, float zIn,
-                unsigned char rowIn, unsigned char secIn, short chanIn = -1, bool rejFlag = false) : dy(static_cast<short>(dyIn * 0x7fff / param::MaxResid)),
-                                                                                                     dz(static_cast<short>(dzIn * 0x7fff / param::MaxResid)),
-                                                                                                     tgSlp(static_cast<short>(tgSlpIn * 0x7fff / param::MaxTgSlp)),
-                                                                                                     y(static_cast<short>(yIn * 0x7fff / param::MaxY)),
-                                                                                                     z(static_cast<short>(zIn * 0x7fff / param::MaxZ)),
+                unsigned char rowIn, unsigned char secIn, short chanIn = -1, bool rejFlag = false) : dy(packToShort(dyIn, param::MaxResid)),
+                                                                                                     dz(packToShort(dzIn, param::MaxResid)),
+                                                                                                     tgSlp(packToShort(tgSlpIn, param::MaxTgSlp, std::abs(tgSlpIn) < param::MaxTgSlp ? 0x7fff - 1 : 0x7fff)),
+                                                                                                     y(packToShort(yIn, param::MaxY)),
+                                                                                                     z(packToShort(zIn, param::MaxZ)),
                                                                                                      row(rowIn),
                                                                                                      sec(secIn),
                                                                                                      channel(chanIn),
@@ -96,6 +98,16 @@ struct UnbinnedResid {
   short channel{-1};    ///< extra channel info (ITS chip ID, TRD chamber, TOF main pad within the sector)
   bool rejected{false}; ///< residual is flagged as rejected in the validateTrack
 
+  /// value (|value| <= maximum) packed into a signed short, rounded to the nearest step. The scale is fixed at 0x7fff steps per maximum
+  /// (decoded everywhere as short * maximum / 0x7fff); limit only clips the result to +-limit (e.g. 0x7fff - 1 for an unclamped tgSlp, to
+  /// keep 0x7fff as the clamp marker), it does not change the scale. A plain cast would truncate towards zero, i.e. move every value
+  /// towards zero by up to one step: the decoded cluster position y + dy then depends on the sign of dy below the step size, which
+  /// biases residuals binned by the cluster position on sub-pad scales
+  static short packToShort(float value, float maximum, int limit = 0x7fff)
+  {
+    return static_cast<short>(std::clamp<long>(std::lround(value * 0x7fff / maximum), -limit, limit));
+  }
+
   /// true if tgSlp was saturated at +-param::MaxTgSlp (scdcalib.clampTgSlp): unclamped values have |tgSlp| <= 0x7fff - 1
   bool isTgSlpClamped() const { return tgSlp == 0x7fff || tgSlp == -0x7fff; }
   /// tgSlp marker of a position-only TPC cluster (scdcalib.keepClustersOnPropFail): no reference track at this cluster,

Detectors/TPC/calibration/SpacePoints/src/TrackInterpolation.cxx

@@ -913,7 +913,7 @@ void TrackInterpolation::interpolateTrack(int iSeed)
           auto dy = yv - trkWorkITS.getY();
           auto dz = zv - trkWorkITS.getZ();
           if ((std::abs(dy) < param::MaxResid) && (std::abs(dz) < param::MaxResid) && (std::abs(trkWorkITS.getY()) < param::MaxY) && (std::abs(trkWorkITS.getZ()) < param::MaxZ) && std::abs(xv) < param::MaxVtxX) {
-            short compXV = static_cast<short>(xv * 0x7fff / param::MaxVtxX);
+            short compXV = UnbinnedResid::packToShort(xv, param::MaxVtxX);
             mClRes.emplace_back(dy, dz, alpha / TMath::Pi(), trkWorkITS.getY(), trkWorkITS.getZ(), 190, -1, compXV);
             if (!gidTable[GTrackID::ITSTPC].isIndexSet()) {
               LOGP(fatal, "ITS-TPC seed index is not set for TOF track");
@@ -985,7 +985,7 @@ int TrackInterpolation::processTRDLayer(const o2::trd::TrackTRD& trkTRD, int iLa
   if (trkData) {
     auto slope = trdSP.getDy();
     if (std::abs(slope) < param::MaxTRDSlope) {
-      trkData->TRDTrkltSlope[iLayer] = slope * 0x7fff / param::MaxTRDSlope;
+      trkData->TRDTrkltSlope[iLayer] = UnbinnedResid::packToShort(slope, param::MaxTRDSlope);
     }
   }
   if (trk64) {
@@ -1308,7 +1308,7 @@ void TrackInterpolation::extrapolateTrack(int iSeed)
           auto dy = yv - trkWorkITS.getY();
           auto dz = zv - trkWorkITS.getZ();
           if ((std::abs(dy) < param::MaxResid) && (std::abs(dz) < param::MaxResid) && (std::abs(trkWorkITS.getY()) < param::MaxY) && (std::abs(trkWorkITS.getZ()) < param::MaxZ) && std::abs(xv) < param::MaxVtxX) {
-            short compXV = static_cast<short>(xv * 0x7fff / param::MaxVtxX);
+            short compXV = UnbinnedResid::packToShort(xv, param::MaxVtxX);
             mClRes.emplace_back(dy, dz, alpha / TMath::Pi(), trkWorkITS.getY(), trkWorkITS.getZ(), 190, -1, compXV);
             if (!gidTableFull[GTrackID::ITSTPC].isIndexSet()) {
               LOGP(fatal, "ITS-TPC seed index is not set for TOF track");