TPC Splines: add limits for SP correction values per TPC row

sgorbuno · cbmsw · commit e32971881fc2 · 2025-10-27T14:21:46.000Z
diff --git a/Detectors/TPC/calibration/src/TPCFastSpaceChargeCorrectionHelper.cxx b/Detectors/TPC/calibration/src/TPCFastSpaceChargeCorrectionHelper.cxx
@@ -144,6 +144,9 @@ void TPCFastSpaceChargeCorrectionHelper::fillSpaceChargeCorrectionFromMap(TPCFas
         float* splineParameters = correction.getSplineData(slice, row);
         const std::vector<o2::gpu::TPCFastSpaceChargeCorrectionMap::CorrectionPoint>& data = mCorrectionMap.getPoints(slice, row);
         int nDataPoints = data.size();
+        auto& info = correction.getSliceRowInfo(slice, row);
+        info.resetMaxValues();
+        info.resetMaxValuesInv();
         if (nDataPoints >= 4) {
           std::vector<double> pointSU(nDataPoints);
           std::vector<double> pointSV(nDataPoints);
@@ -156,6 +159,8 @@ void TPCFastSpaceChargeCorrectionHelper::fillSpaceChargeCorrectionFromMap(TPCFas
             pointCorr[3 * i + 0] = dx;
             pointCorr[3 * i + 1] = du;
             pointCorr[3 * i + 2] = dv;
+            info.updateMaxValues(2. * dx, 2. * du, 2. * dv);
+            info.updateMaxValuesInv(-2. * dx, -2. * du, -2. * dv);
           }
           helper.approximateDataPoints(spline, splineParameters, 0., spline.getGridX1().getUmax(), 0., spline.getGridX2().getUmax(), &pointSU[0],
                                        &pointSV[0], &pointCorr[0], nDataPoints);
@@ -767,9 +772,20 @@ std::unique_ptr<o2::gpu::TPCFastSpaceChargeCorrection> TPCFastSpaceChargeCorrect
 
             double yStep = (yLast - yFirst) / 2;
 
+            double zFirst = z - dz / 2.;
+            double zLast = z + dz / 2.;
+            double zStep = (zLast - zFirst) / 2.;
+
+            if (0) { // no smoothing
+              yFirst = y;
+              yLast = y;
+              zFirst = z;
+              zLast = z;
+            }
+
             for (double py = yFirst; py <= yLast + yStep / 2.; py += yStep) {
 
-              for (double pz = z - dz / 2.; pz <= z + dz / 2. + 1.e-4; pz += dz / 2.) {
+              for (double pz = zFirst; pz <= zLast + zStep / 2.; pz += zStep) {
                 map.addCorrectionPoint(iRoc, iRow, py, pz, correctionX, correctionY,
                                        correctionZ);
               }
diff --git a/GPU/TPCFastTransformation/TPCFastSpaceChargeCorrection.h b/GPU/TPCFastTransformation/TPCFastSpaceChargeCorrection.h
@@ -58,15 +58,64 @@ class TPCFastSpaceChargeCorrection : public FlatObject
   };
 
   struct SliceRowInfo {
-    float gridU0{0.f};           //< U coordinate of the U-grid start
-    float scaleUtoGrid{0.f};     //< scale U to U-grid coordinate
-    float gridV0{0.f};           ///< V coordinate of the V-grid start
-    float scaleVtoGrid{0.f};     //< scale V to V-grid coordinate
-    float gridCorrU0{0.f};       ///< U coordinate of the U-grid start for corrected U
-    float scaleCorrUtoGrid{0.f}; ///< scale corrected U to U-grid coordinate
-    float gridCorrV0{0.f};       ///< V coordinate of the V-grid start for corrected V
-    float scaleCorrVtoGrid{0.f}; ///< scale corrected V to V-grid coordinate
+    float gridU0{0.f};                        //< U coordinate of the U-grid start
+    float scaleUtoGrid{0.f};                  //< scale U to U-grid coordinate
+    float gridV0{0.f};                        ///< V coordinate of the V-grid start
+    float scaleVtoGrid{0.f};                  //< scale V to V-grid coordinate
+    float gridCorrU0{0.f};                    ///< U coordinate of the U-grid start for corrected U
+    float scaleCorrUtoGrid{0.f};              ///< scale corrected U to U-grid coordinate
+    float gridCorrV0{0.f};                    ///< V coordinate of the V-grid start for corrected V
+    float scaleCorrVtoGrid{0.f};              ///< scale corrected V to V-grid coordinate
+    float maxCorr[3]{10.f, 10.f, 10.f};       ///< max correction for dX, dU, dV
+    float minCorr[3]{-10.f, -10.f, -10.f};    ///< min correction for dX, dU, dV
+    float maxInvCorr[3]{10.f, 10.f, 10.f};    ///< max inverse correction for dX, dU, dV
+    float minInvCorr[3]{-10.f, -10.f, -10.f}; ///< min inverse correction for dX, dU, dV
     RowActiveArea activeArea;
+
+    void resetMaxValues()
+    {
+      maxCorr[0] = 1.f;
+      minCorr[0] = -1.f;
+      maxCorr[1] = 1.f;
+      minCorr[1] = -1.f;
+      maxCorr[2] = 1.f;
+      minCorr[2] = -1.f;
+    }
+
+    void updateMaxValues(float dx, float du, float dv)
+    {
+      maxCorr[0] = GPUCommonMath::Max(maxCorr[0], dx);
+      minCorr[0] = GPUCommonMath::Min(minCorr[0], dx);
+
+      maxCorr[1] = GPUCommonMath::Max(maxCorr[1], du);
+      minCorr[1] = GPUCommonMath::Min(minCorr[1], du);
+
+      maxCorr[2] = GPUCommonMath::Max(maxCorr[2], dv);
+      minCorr[2] = GPUCommonMath::Min(minCorr[2], dv);
+    }
+
+    void resetMaxValuesInv()
+    {
+      maxInvCorr[0] = 1.f;
+      minInvCorr[0] = -1.f;
+      maxInvCorr[1] = 1.f;
+      minInvCorr[1] = -1.f;
+      maxInvCorr[2] = 1.f;
+      minInvCorr[2] = -1.f;
+    }
+
+    void updateMaxValuesInv(float dx, float du, float dv)
+    {
+      maxInvCorr[0] = GPUCommonMath::Max(maxInvCorr[0], dx);
+      minInvCorr[0] = GPUCommonMath::Min(minInvCorr[0], dx);
+
+      maxInvCorr[1] = GPUCommonMath::Max(maxInvCorr[1], du);
+      minInvCorr[1] = GPUCommonMath::Min(minInvCorr[1], du);
+
+      maxInvCorr[2] = GPUCommonMath::Max(maxInvCorr[2], dv);
+      minInvCorr[2] = GPUCommonMath::Min(minInvCorr[2], dv);
+    }
+
     ClassDefNV(SliceRowInfo, 2);
   };
 
@@ -397,12 +446,10 @@ GPUdi() int32_t TPCFastSpaceChargeCorrection::getCorrection(int32_t slice, int32
   convUVtoGrid(slice, row, u, v, gridU, gridV);
   float dxuv[3];
   spline.interpolateU(splineData, gridU, gridV, dxuv);
-  if (CAMath::Abs(dxuv[0]) > 100 || CAMath::Abs(dxuv[1]) > 100 || CAMath::Abs(dxuv[2]) > 100) {
-    dxuv[0] = dxuv[1] = dxuv[2] = 0;
-  }
-  dx = dxuv[0];
-  du = dxuv[1];
-  dv = dxuv[2];
+  const auto& info = getSliceRowInfo(slice, row);
+  dx = GPUCommonMath::Max(info.minCorr[0], GPUCommonMath::Min(info.maxCorr[0], dxuv[0]));
+  du = GPUCommonMath::Max(info.minCorr[1], GPUCommonMath::Min(info.maxCorr[1], dxuv[1]));
+  dv = GPUCommonMath::Max(info.minCorr[2], GPUCommonMath::Min(info.maxCorr[2], dxuv[2]));
   return 0;
 }
 
@@ -414,12 +461,10 @@ GPUdi() int32_t TPCFastSpaceChargeCorrection::getCorrectionOld(int32_t slice, in
   convUVtoGrid(slice, row, u, v, gridU, gridV);
   float dxuv[3];
   spline.interpolateUold(splineData, gridU, gridV, dxuv);
-  if (CAMath::Abs(dxuv[0]) > 100 || CAMath::Abs(dxuv[1]) > 100 || CAMath::Abs(dxuv[2]) > 100) {
-    dxuv[0] = dxuv[1] = dxuv[2] = 0;
-  }
-  dx = dxuv[0];
-  du = dxuv[1];
-  dv = dxuv[2];
+  const auto& info = getSliceRowInfo(slice, row);
+  dx = GPUCommonMath::Max(info.minCorr[0], GPUCommonMath::Min(info.maxCorr[0], dxuv[0]));
+  du = GPUCommonMath::Max(info.minCorr[1], GPUCommonMath::Min(info.maxCorr[1], dxuv[1]));
+  dv = GPUCommonMath::Max(info.minCorr[2], GPUCommonMath::Min(info.maxCorr[2], dxuv[2]));
   return 0;
 }
 
@@ -433,9 +478,8 @@ GPUdi() void TPCFastSpaceChargeCorrection::getCorrectionInvCorrectedX(
   const float* splineData = getSplineData(slice, row, 1);
   float dx = 0;
   spline.interpolateU(splineData, gridU, gridV, &dx);
-  if (CAMath::Abs(dx) > 100) {
-    dx = 0;
-  }
+  const auto& info = getSliceRowInfo(slice, row);
+  dx = GPUCommonMath::Max(info.minInvCorr[0], GPUCommonMath::Min(info.maxInvCorr[0], dx));
   x = mGeo.getRowInfo(row).x + dx;
 }
 
@@ -450,9 +494,9 @@ GPUdi() void TPCFastSpaceChargeCorrection::getCorrectionInvUV(
 
   float duv[2];
   spline.interpolateU(splineData, gridU, gridV, duv);
-  if (CAMath::Abs(duv[0]) > 100 || CAMath::Abs(duv[1]) > 100) {
-    duv[0] = duv[1] = 0;
-  }
+  const auto& info = getSliceRowInfo(slice, row);
+  duv[0] = GPUCommonMath::Max(info.minInvCorr[1], GPUCommonMath::Min(info.maxInvCorr[1], duv[0]));
+  duv[1] = GPUCommonMath::Max(info.minInvCorr[2], GPUCommonMath::Min(info.maxInvCorr[2], duv[1]));
   nomU = corrU - duv[0];
   nomV = corrV - duv[1];
 }
diff --git a/GPU/TPCFastTransformation/macro/TPCFastTransformInit.C b/GPU/TPCFastTransformation/macro/TPCFastTransformInit.C
@@ -21,8 +21,6 @@
 /// root -l TPCFastTransformInit.C'("debugVoxRes.root")'
 ///
 
-#include "Algorithm/RangeTokenizer.h"
-
 #if !defined(__CLING__) || defined(__ROOTCLING__)
 
 #include <filesystem>
@@ -41,6 +39,8 @@
 #include "TPCCalibration/TPCFastSpaceChargeCorrectionHelper.h"
 #endif
 
+#include "Algorithm/RangeTokenizer.h"
+
 using namespace o2::tpc;
 using namespace o2::gpu;
 
@@ -99,8 +99,9 @@ void TPCFastTransformInit(const char* fileName = "debugVoxRes.root",
   trackResiduals.setZ2XBinning(z2xBins);
   trackResiduals.init();
 
-  {
-    std::cout << "input track residuals: " << std::endl;
+  { // debug output
+
+    std::cout << " ===== input track residuals ==== " << std::endl;
     std::cout << "voxel tree y2xBins: " << y2xBins.size() << std::endl;
 
     for (auto y2x : y2xBins) {
@@ -127,6 +128,7 @@ void TPCFastTransformInit(const char* fileName = "debugVoxRes.root",
     for (int i = 0; i < nZ2Xbins; i++) {
       std::cout << "getZ2X(bin) : " << trackResiduals.getZ2X(i) << std::endl;
     }
+    std::cout << " ==================================== " << std::endl;
   }
 
   std::cout << "create fast transformation ... " << std::endl;
@@ -310,15 +312,16 @@ void TPCFastTransformInit(const char* fileName = "debugVoxRes.root",
     geo.convUVtoLocal(iRoc, u + cu, v + cv, cy, cz);
     cy -= y;
     cz -= z;
+
     double d[3] = {cx - correctionX, cy - correctionY, cz - correctionZ};
     if (voxEntries >= 1.) {
       for (int i = 0; i < 3; i++) {
         if (fabs(maxDiff[i]) < fabs(d[i])) {
           maxDiff[i] = d[i];
           maxDiffRoc[i] = iRoc;
           maxDiffRow[i] = iRow;
-          std::cout << " roc " << iRoc << " row " << iRow << " xyz " << i
-                    << " diff " << d[i] << " entries " << voxEntries << " y " << y2xBin << " z " << z2xBin << std::endl;
+          // std::cout << " roc " << iRoc << " row " << iRow << " xyz " << i
+          //  << " diff " << d[i] << " entries " << voxEntries << " y " << y2xBin << " z " << z2xBin << std::endl;
         }
         sumDiff[i] += d[i] * d[i];
       }
