TPC Splines: fix the inverse correction

Author: sgorbuno

Detectors/TPC/calibration/src/TPCFastSpaceChargeCorrectionHelper.cxx

@@ -146,7 +146,6 @@ void TPCFastSpaceChargeCorrectionHelper::fillSpaceChargeCorrectionFromMap(TPCFas
         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);
@@ -160,7 +159,6 @@ void TPCFastSpaceChargeCorrectionHelper::fillSpaceChargeCorrectionFromMap(TPCFas
             pointCorr[3 * i + 1] = du;
             pointCorr[3 * i + 2] = dv;
             info.updateMaxValues(20. * dx, 20. * du, 20. * dv);
-            info.updateMaxValuesInv(-20. * dx, -20. * du, -20. * dv);
           }
           helper.approximateDataPoints(spline, splineParameters, 0., spline.getGridX1().getUmax(), 0., spline.getGridX2().getUmax(), &pointSU[0],
                                        &pointSV[0], &pointCorr[0], nDataPoints);
@@ -908,46 +906,60 @@ void TPCFastSpaceChargeCorrectionHelper::initInverse(std::vector<o2::gpu::TPCFas
 
   for (int slice = 0; slice < mGeo.getNumberOfSlices(); slice++) {
     // LOG(info) << "inverse transform for slice " << slice ;
-    double vLength = (slice < mGeo.getNumberOfSlicesA()) ? mGeo.getTPCzLengthA() : mGeo.getTPCzLengthC();
 
     auto myThread = [&](int iThread) {
       Spline2DHelper<float> helper;
       std::vector<float> splineParameters;
-      ChebyshevFit1D chebFitterX, chebFitterU, chebFitterV;
 
       for (int row = iThread; row < mGeo.getNumberOfRows(); row += mNthreads) {
         TPCFastSpaceChargeCorrection::SplineType spline = correction.getSpline(slice, row);
         helper.setSpline(spline, 10, 10);
-        std::vector<double> dataPointCU, dataPointCV, dataPointF;
-
-        float u0, u1, v0, v1;
-        mGeo.convScaledUVtoUV(slice, row, 0., 0., u0, v0);
-        mGeo.convScaledUVtoUV(slice, row, 1., 1., u1, v1);
 
         double x = mGeo.getRowInfo(row).x;
-        int nPointsU = (spline.getGridX1().getNumberOfKnots() - 1) * 10;
-        int nPointsV = (spline.getGridX2().getNumberOfKnots() - 1) * 10;
-
-        double stepU = (u1 - u0) / (nPointsU - 1);
-        double stepV = (v1 - v0) / (nPointsV - 1);
+        auto& sliceRowInfo = correction.getSliceRowInfo(slice, row);
 
-        if (prn) {
-          LOG(info) << "u0 " << u0 << " u1 " << u1 << " v0 " << v0 << " v1 " << v1;
+        std::vector<double> gridU;
+        {
+          const auto& grid = spline.getGridX1();
+          for (int i = 0; i < grid.getNumberOfKnots(); i++) {
+            if (i == grid.getNumberOfKnots() - 1) {
+              gridU.push_back(grid.getKnot(i).u);
+              break;
+            }
+            for (double s = 1.; s > 0.; s -= 0.1) {
+              gridU.push_back(s * grid.getKnot(i).u + (1. - s) * grid.getKnot(i + 1).u);
+            }
+          }
+        }
+        std::vector<double> gridV;
+        {
+          const auto& grid = spline.getGridX2();
+          for (int i = 0; i < grid.getNumberOfKnots(); i++) {
+            if (i == grid.getNumberOfKnots() - 1) {
+              gridV.push_back(grid.getKnot(i).u);
+              break;
+            }
+            for (double s = 1.; s > 0.; s -= 0.1) {
+              gridV.push_back(s * grid.getKnot(i).u + (1. - s) * grid.getKnot(i + 1).u);
+            }
+          }
         }
-        TPCFastSpaceChargeCorrection::RowActiveArea& area = correction.getSliceRowInfo(slice, row).activeArea;
+
+        std::vector<double> dataPointCU, dataPointCV, dataPointF;
+        dataPointCU.reserve(gridU.size() * gridV.size());
+        dataPointCV.reserve(gridU.size() * gridV.size());
+        dataPointF.reserve(gridU.size() * gridV.size());
+
+        TPCFastSpaceChargeCorrection::RowActiveArea& area = sliceRowInfo.activeArea;
         area.cuMin = 1.e10;
         area.cuMax = -1.e10;
+        double cvMin = 1.e10;
 
-        /*
-        v1 = area.vMax;
-        stepV = (v1 - v0) / (nPointsU - 1);
-        if (stepV < 1.f) {
-          stepV = 1.f;
-        }
-        */
+        for (int iu = 0; iu < gridU.size(); iu++) {
+          for (int iv = 0; iv < gridV.size(); iv++) {
+            float u, v;
+            correction.convGridToUV(slice, row, gridU[iu], gridV[iv], u, v);
 
-        for (double u = u0; u < u1 + stepU; u += stepU) {
-          for (double v = v0; v < v1 + stepV; v += stepV) {
             float dx, du, dv;
             correction.getCorrection(slice, row, u, v, dx, du, dv);
             dx *= scaling[0];
@@ -976,39 +988,41 @@ void TPCFastSpaceChargeCorrectionHelper::initInverse(std::vector<o2::gpu::TPCFas
             dataPointF.push_back(dx);
             dataPointF.push_back(du);
             dataPointF.push_back(dv);
-
-            if (prn) {
-              LOG(info) << "measurement cu " << cu << " cv " << cv << " dx " << dx << " du " << du << " dv " << dv;
-            }
-          } // v
-        }   // u
+          }
+        }
 
         if (area.cuMax - area.cuMin < 0.2) {
           area.cuMax = .1;
           area.cuMin = -.1;
         }
-        if (area.cvMax < 0.1) {
+        if (area.cvMax - cvMin < 0.2) {
           area.cvMax = .1;
+          cvMin = -.1;
         }
+
         if (prn) {
           LOG(info) << "slice " << slice << " row " << row << " max drift L = " << correction.getMaxDriftLength(slice, row)
                     << " active area: cuMin " << area.cuMin << " cuMax " << area.cuMax << " vMax " << area.vMax << " cvMax " << area.cvMax;
         }
 
-        TPCFastSpaceChargeCorrection::SliceRowInfo& info = correction.getSliceRowInfo(slice, row);
-        info.gridCorrU0 = area.cuMin;
-        info.scaleCorrUtoGrid = spline.getGridX1().getUmax() / (area.cuMax - area.cuMin);
-        info.scaleCorrVtoGrid = spline.getGridX2().getUmax() / area.cvMax;
+        // define the grid for the inverse correction
 
-        info.gridCorrU0 = u0;
-        info.gridCorrV0 = info.gridV0;
-        info.scaleCorrUtoGrid = spline.getGridX1().getUmax() / (u1 - info.gridCorrU0);
-        info.scaleCorrVtoGrid = spline.getGridX2().getUmax() / (v1 - info.gridCorrV0);
+        sliceRowInfo.gridCorrU0 = area.cuMin;
+        sliceRowInfo.gridCorrV0 = cvMin;
+        sliceRowInfo.scaleCorrUtoGrid = spline.getGridX1().getUmax() / (area.cuMax - area.cuMin);
+        sliceRowInfo.scaleCorrVtoGrid = spline.getGridX2().getUmax() / area.cvMax;
+
+        /*
+        sliceRowInfo.gridCorrU0 = sliceRowInfo.gridU0;
+        sliceRowInfo.gridCorrV0 = sliceRowInfo.gridV0;
+        sliceRowInfo.scaleCorrUtoGrid = sliceRowInfo.scaleUtoGrid;
+        sliceRowInfo.scaleCorrVtoGrid = sliceRowInfo.scaleVtoGrid;
+        */
 
         int nDataPoints = dataPointCU.size();
         for (int i = 0; i < nDataPoints; i++) {
-          dataPointCU[i] = (dataPointCU[i] - info.gridCorrU0) * info.scaleCorrUtoGrid;
-          dataPointCV[i] = (dataPointCV[i] - info.gridCorrV0) * info.scaleCorrVtoGrid;
+          dataPointCU[i] = (dataPointCU[i] - sliceRowInfo.gridCorrU0) * sliceRowInfo.scaleCorrUtoGrid;
+          dataPointCV[i] = (dataPointCV[i] - sliceRowInfo.gridCorrV0) * sliceRowInfo.scaleCorrVtoGrid;
         }
 
         splineParameters.resize(spline.getNumberOfParameters());

GPU/TPCFastTransformation/TPCFastSpaceChargeCorrection.cxx

@@ -19,6 +19,7 @@
 
 #if !defined(GPUCA_GPUCODE)
 #include <iostream>
+#include <string>
 #include <cmath>
 #include "Spline2DHelper.h"
 #endif
@@ -514,25 +515,54 @@ double TPCFastSpaceChargeCorrection::testInverse(bool prn)
   tpcR2max = tpcR2max / cos(2 * M_PI / mGeo.getNumberOfSlicesA() / 2) + 1.;
   tpcR2max = tpcR2max * tpcR2max;
 
-  double maxDtpc[3] = {0, 0, 0};
-  double maxD = 0;
+  struct MaxValue {
+    double V{0.};
+    int Roc{-1};
+    int Row{-1};
+
+    void update(double v, int roc, int row)
+    {
+      if (fabs(v) > fabs(V)) {
+        V = v;
+        Roc = roc;
+        Row = row;
+      }
+    }
+    void update(const MaxValue& other)
+    {
+      update(other.V, other.Roc, other.Row);
+    }
+
+    std::string toString()
+    {
+      std::stringstream ss;
+      ss << V << "(" << Roc << "," << Row << ")";
+      return ss.str();
+    }
+  };
+
+  MaxValue maxDtpc[3];
+  MaxValue maxD;
 
   for (int32_t slice = 0; slice < mGeo.getNumberOfSlices(); slice++) {
     if (prn) {
       LOG(info) << "check inverse transform for slice " << slice;
     }
-    double vLength = (slice < mGeo.getNumberOfSlicesA()) ? mGeo.getTPCzLengthA() : mGeo.getTPCzLengthC();
-    double maxDslice[3] = {0, 0, 0};
+    double vLength = mGeo.getTPCzLength(slice);
+    MaxValue maxDslice[3];
     for (int32_t row = 0; row < mGeo.getNumberOfRows(); row++) {
       float u0, u1, v0, v1;
       mGeo.convScaledUVtoUV(slice, row, 0., 0., u0, v0);
       mGeo.convScaledUVtoUV(slice, row, 1., 1., u1, v1);
       double x = mGeo.getRowInfo(row).x;
       double stepU = (u1 - u0) / 100.;
       double stepV = (v1 - v0) / 100.;
-      double maxDrow[3] = {0, 0, 0};
+      MaxValue maxDrow[3];
       for (double u = u0; u < u1; u += stepU) {
         for (double v = v0; v < v1; v += stepV) {
+          if (v < getSliceRowInfo(slice, row).gridV0) {
+            continue;
+          }
           float dx, du, dv;
           getCorrection(slice, row, u, v, dx, du, dv);
           double cx = x + dx;
@@ -545,11 +575,9 @@ double TPCFastSpaceChargeCorrection::testInverse(bool prn)
           float nx, nu, nv;
           getCorrectionInvCorrectedX(slice, row, cu, cv, nx);
           getCorrectionInvUV(slice, row, cu, cv, nu, nv);
-          double d[3] = {nx - cx, nu - u, nv - v};
+          double d[3] = {(cx - nx) - dx, (cu - nu) - du, (cv - nv) - dv};
           for (int32_t i = 0; i < 3; i++) {
-            if (fabs(d[i]) > fabs(maxDrow[i])) {
-              maxDrow[i] = d[i];
-            }
+            maxDrow[i].update(d[i], slice, row);
           }
 
           if (0 && prn && fabs(d[0]) + fabs(d[1]) + fabs(d[2]) > 0.1) {
@@ -560,32 +588,26 @@ double TPCFastSpaceChargeCorrection::testInverse(bool prn)
           }
         }
       }
-      if (0 && prn) {
+      if (1 && prn) {
         LOG(info) << "slice " << slice << " row " << row
-                  << " dx " << maxDrow[0] << " du " << maxDrow[1] << " dv " << maxDrow[2];
+                  << " dx " << maxDrow[0].V << " du " << maxDrow[1].V << " dv " << maxDrow[2].V;
       }
       for (int32_t i = 0; i < 3; i++) {
-        if (fabs(maxDslice[i]) < fabs(maxDrow[i])) {
-          maxDslice[i] = maxDrow[i];
-        }
-        if (fabs(maxDtpc[i]) < fabs(maxDrow[i])) {
-          maxDtpc[i] = maxDrow[i];
-        }
-        if (fabs(maxD) < fabs(maxDrow[i])) {
-          maxD = maxDrow[i];
-        }
+        maxDslice[i].update(maxDrow[i]);
+        maxDtpc[i].update(maxDrow[i]);
+        maxD.update(maxDrow[i]);
       }
     }
     if (prn) {
-      LOG(info) << "inverse correction: slice " << slice
-                << " dx " << maxDslice[0] << " du " << maxDslice[1] << " dv " << maxDslice[2];
+      LOG(info) << "inverse correction: slice " << slice << ". Max deviations: "
+                << " dx " << maxDslice[0].toString() << " du " << maxDslice[1].toString() << " dv " << maxDslice[2].toString();
     }
   } // slice
 
   LOG(info) << "Test inverse TPC correction. max deviations: "
-            << " dx " << maxDtpc[0] << " du " << maxDtpc[1] << " dv " << maxDtpc[2] << " cm";
+            << " dx " << maxDtpc[0].toString() << " du " << maxDtpc[1].toString() << " dv " << maxDtpc[2].toString() << " cm";
 
-  return maxD;
+  return maxD.V;
 }
 
 #endif // GPUCA_GPUCODE

GPU/TPCFastTransformation/TPCFastSpaceChargeCorrection.h

@@ -68,8 +68,6 @@ class TPCFastSpaceChargeCorrection : public FlatObject
     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()
@@ -94,28 +92,6 @@ class TPCFastSpaceChargeCorrection : public FlatObject
       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);
   };
 
@@ -494,7 +470,15 @@ GPUdi() void TPCFastSpaceChargeCorrection::getCorrectionInvCorrectedX(
   float dx = 0;
   spline.interpolateU(splineData, gridU, gridV, &dx);
   const auto& info = getSliceRowInfo(slice, row);
-  dx = GPUCommonMath::Max(info.minInvCorr[0], GPUCommonMath::Min(info.maxInvCorr[0], dx));
+
+  float s = corrV / info.gridCorrV0;
+  if (s < 0.) {
+    s = 0.;
+  }
+  if (s > 1.) {
+    s = 1.;
+  }
+  dx = GPUCommonMath::Clamp(s * dx, info.minCorr[0], info.maxCorr[0]);
   x = mGeo.getRowInfo(row).x + dx;
 }
 
@@ -510,8 +494,15 @@ GPUdi() void TPCFastSpaceChargeCorrection::getCorrectionInvUV(
   float duv[2];
   spline.interpolateU(splineData, gridU, gridV, duv);
   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]));
+  float s = corrV / info.gridCorrV0;
+  if (s < 0.) {
+    s = 0.;
+  }
+  if (s > 1.) {
+    s = 1.;
+  }
+  duv[0] = GPUCommonMath::Clamp(s * duv[0], info.minCorr[1], info.maxCorr[1]);
+  duv[1] = GPUCommonMath::Clamp(s * duv[1], info.minCorr[2], info.maxCorr[2]);
   nomU = corrU - duv[0];
   nomV = corrV - duv[1];
 }