TPC: Digitizer recalculates the distortions with the derivative corrections without inverting the average map (new default)

Digitizer::recalculateDistortions() (distortion type 2) inverted the average
distortions to corrections, added the scaled corrections of the derivative
map and inverted the sum back to distortions. The first inversion is the
slowest and the least precise step.

SpaceCharge::scaleGlobalDistWithDerivativeCorr() applies the scaled
derivative corrections directly to the average distortions. With
T(x) = x + D(x) the new distortion at the vertex p is D_new(p) = T(x) - p,
where x is the fixed point of x = p - scale * C_der(T(x)) (at most 20
iterations, stop at a change of the start point < 2 um). The distorted point
y = T(x) then fulfils y + C(y) + scale * C_der(y) = p for the inverse C of
D, i.e. the distortions are consistent with the scaled corrections at the
distorted point as they are applied in the reconstruction (C_avg + s * C_der
at the measured point). Only the derivative distortions are inverted.
The r-phi distortion is calculated from the cross and dot product of the
vertex and the distorted point: the difference of two angles wrapped with
the single precision 2 pi gave an offset of 0.4 um for phi > pi.
The function returns false if the global distortions of the object or the
global corrections of the derivative map are not set (distortions
unchanged); the Digitizer then stops with a fatal instead of continuing with
the unscaled distortions.

- Digitizer::setRecalcDistortionsMethod() and option
  --recalculate-distortions-method of the TPC digitizer: 0 = inversion of
  the scaled corrections (previous method), 1 = new method (default). Other
  values are a fatal error. The distortions of a digitisation with lumi
  scaling (distortion type 2) change.
- ClassDefNV of the Digitizer 3 -> 4 (new member)

Validation with run 559917 (distortion maps and correction splines of the
same MC map production): closure with the correction splines of the
reconstruction, |c + D(c) - q| for the measured point q and the corrected
point c = q + C_avg(q) + s * C_der(q), 99% quantile in local X / Y / Z
(5 < |z| < 245 cm, pad rows >= 10):

  scale s   new method           previous method      no scaling
  -0.32     0.62/0.94/2.00 mm    1.02/1.39/2.73 mm    0.58/0.96/1.98 mm
  -2        0.83/1.01/2.17 mm    1.07/1.39/2.74 mm

The new method closes as well as the unscaled average map with its own
correction splines. The previous method has larger tails (99.9% in r-phi at
s = -2: 4.4 mm vs 2.0 mm) and 10x more points with a closure > 1 cm at
s = -2 (2396 vs 239). The remaining large residuals of the new method are at
the innermost pad rows (IFC).

Time of the recalculation (CCDB maps of run 559917, 65x257x720, both
sides): 85.9 s -> 27.1 s with 1 thread, 22.9 s -> 7.3 s with 4 threads
(default of n-threads-distortions).

Tests: SpaceCharge_scaleWithDerivative_test (analytical fields, scale -3, -1
and 2: closure with the scaled corrections, agreement with the inversion of
the scaled corrections, return value without distortions or derivative
corrections), SpaceCharge_scaleWithDerivativeFloat_test (single precision:
r-phi distortion of radial distortions and radial derivative corrections for
phi <= pi and phi > pi).

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

Author: matthias-kleiner

Detectors/TPC/simulation/include/TPCSimulation/Digitizer.h

@@ -138,6 +138,13 @@ class Digitizer
   void setMeanLumiDistortions(float meanLumi);
   void setMeanLumiDistortionsDerivative(float meanLumi);
 
+  /// method of recalculateDistortions()
+  /// \param method 0: invert the average and the derivative distortions to corrections and invert the scaled corrections C_avg + s * C_der to distortions,
+  ///               1: invert only the derivative distortions and apply the scaled derivative corrections directly to the average distortions (SpaceCharge::scaleGlobalDistWithDerivativeCorr).
+  ///               Other values are a fatal error
+  void setRecalcDistortionsMethod(int method);
+  int getRecalcDistortionsMethod() const { return mRecalcDistortionsMethod; }
+
   /// in case of scaled distortions, the distortions can be recalculated to ensure consistent distortions and corrections
   void recalculateDistortions();
 
@@ -155,7 +162,8 @@ class Digitizer
   int mDistortionScaleType = 0;        ///< type=0: no scaling of distortions, type=1 distortions without any scaling, type=2 distortions scaling with lumi
   float mLumiScaleFactor = 0;          ///< value used to scale the derivative map
   bool mUseScaledDistortions = false;  ///< whether the distortions are already scaled
-  ClassDefNV(Digitizer, 3);
+  int mRecalcDistortionsMethod = 1;    ///< method of recalculateDistortions(): 0 = inversion of the scaled corrections, 1 = derivative corrections applied to the distortions
+  ClassDefNV(Digitizer, 4);
 };
 } // namespace tpc
 } // namespace o2

Detectors/TPC/simulation/src/Digitizer.cxx

@@ -249,6 +249,14 @@ void Digitizer::setMeanLumiDistortionsDerivative(float meanLumi)
   mSpaceChargeDer->setMeanLumi(meanLumi);
 }
 
+void Digitizer::setRecalcDistortionsMethod(int method)
+{
+  if ((method != 0) && (method != 1)) {
+    LOGP(fatal, "Unknown method {} for the recalculation of the distortions (0: inversion of the scaled corrections, 1: scaled derivative corrections applied to the distortions)", method);
+  }
+  mRecalcDistortionsMethod = method;
+}
+
 void Digitizer::recalculateDistortions()
 {
   if (!mSpaceCharge || !mSpaceChargeDer) {
@@ -260,15 +268,24 @@ void Digitizer::recalculateDistortions()
   if (mSpaceCharge->getMeanLumi() != CorrMapParam::Instance().lumiInst) {
     for (int iside = 0; iside < 2; ++iside) {
       const o2::tpc::Side side = (iside == 0) ? Side::A : Side::C;
-      // this needs to be done only once
-      LOGP(info, "Calculating corrections for average distortions");
-      mSpaceCharge->calcGlobalCorrWithGlobalDistIterative(side, nullptr, 0);
-
       LOGP(info, "Calculating corrections for derivative distortions");
       mSpaceChargeDer->calcGlobalCorrWithGlobalDistIterative(side, nullptr, 0);
 
-      LOGP(info, "Calculating scaled distortions with scaling factor {}", mLumiScaleFactor);
-      mSpaceCharge->calcGlobalDistWithGlobalCorrIterativeLinearCartesian(side, mSpaceChargeDer.get(), mLumiScaleFactor);
+      if (mRecalcDistortionsMethod == 0) {
+        // invert the average distortions, then invert the scaled corrections C_avg + s * C_der
+        LOGP(info, "Calculating corrections for average distortions");
+        mSpaceCharge->calcGlobalCorrWithGlobalDistIterative(side, nullptr, 0);
+
+        LOGP(info, "Calculating scaled distortions with scaling factor {}", mLumiScaleFactor);
+        mSpaceCharge->calcGlobalDistWithGlobalCorrIterativeLinearCartesian(side, mSpaceChargeDer.get(), mLumiScaleFactor);
+      } else {
+        // apply the scaled derivative corrections directly to the average distortions (no inversion of the average distortions)
+        LOGP(info, "Scaling distortions with the derivative corrections and scaling factor {}", mLumiScaleFactor);
+        if (!mSpaceCharge->scaleGlobalDistWithDerivativeCorr(side, *mSpaceChargeDer, mLumiScaleFactor)) {
+          // continuing would use the unscaled average distortions for the rest of the job
+          LOGP(fatal, "Distortions of side {} could not be scaled with the derivative map", static_cast<int>(side));
+        }
+      }
     }
     // set new lumi of avg map
     mSpaceCharge->setMeanLumi(CorrMapParam::Instance().lumiInst);

Detectors/TPC/spacecharge/include/TPCSpaceCharge/SpaceCharge.h

@@ -409,6 +409,24 @@ class SpaceCharge
   void calcGlobalDistWithGlobalCorrIterativeLinearCartesian(const Side side, const SpaceCharge<DataT>* scSCale = nullptr, float scale = 0, const int maxIter = 100, const DataT approachX = 1, const DataT approachY = 1, const DataT approachZ = 1, const DataT diffCorr = 50e-6);
   void calcGlobalDistWithGlobalCorrIterativeLinearCartesian(const SpaceCharge<DataT>* scSCale = nullptr, float scale = 0, const int maxIter = 100, const DataT approachX = 1, const DataT approachY = 1, const DataT approachZ = 1, const DataT diffCorr = 50e-6);
 
+  /// Scale the global distortions with the global corrections of a derivative map, consistent with the scaled corrections C(y) + scale * C_der(y) at the
+  /// distorted point y (as done in the tracking), without inverting the global distortions of this object.
+  /// With T(x) = x + D(x) (the global distortions of this object), the new distortion at the vertex p is D_new(p) = T(x) - p, where x is the fixed point of
+  /// x = p - scale * C_der(T(x)) (iterated starting at x = p). Then y = T(x) fulfils y + C(y) = x for the inverse C of D, and y + C(y) + scale * C_der(y) = p.
+  /// The iteration converges fast since C_der is small and smooth (each step reduces the error by about |scale| * |dC_der/dy|).
+  /// \param side side of the TPC
+  /// \param scDer derivative map: its global corrections must be set (e.g. with calcGlobalCorrWithGlobalDistIterative())
+  /// \param scale scaling factor of the derivative map
+  /// \param maxIter maximum number of fixed-point iterations per vertex
+  /// \param diffCorr the iteration stops when the start point x changes by less than this value (cm) in x, y and z. The default (2 um) is above the
+  ///                 rounding noise of the cartesian corrections in single precision (O(0.5 um) at the outer radius) times the scaling factor
+  /// \return true if the distortions were scaled, false if the global distortions of this object or the global corrections of the derivative map are not set (the distortions are not changed)
+  bool scaleGlobalDistWithDerivativeCorr(const Side side, const SpaceCharge<DataT>& scDer, const float scale, const int maxIter = 20, const DataT diffCorr = 2e-4);
+
+  /// Scale the global distortions of both sides with the global corrections of a derivative map (see the function for one side)
+  /// \return true if the distortions of both sides were scaled, false if at least one side was not scaled
+  bool scaleGlobalDistWithDerivativeCorr(const SpaceCharge<DataT>& scDer, const float scale, const int maxIter = 20, const DataT diffCorr = 2e-4);
+
   /// Calculate global corrections using the global distortions by scaling with second distortion object, which will be consistent with scaled corrections in cartesian coordinates
   /// \param side Side of the TPC
   /// \param scSCale possible second sc object

Detectors/TPC/spacecharge/src/SpaceCharge.cxx

@@ -943,6 +943,105 @@ void SpaceCharge<DataT>::calcGlobalDistWithGlobalCorrIterativeLinearCartesian(co
   }
 }
 
+template <typename DataT>
+bool SpaceCharge<DataT>::scaleGlobalDistWithDerivativeCorr(const Side side, const SpaceCharge<DataT>& scDer, const float scale, const int maxIter, const DataT diffCorr)
+{
+  if (!mGlobalDistdR[side].getNDataPoints() || !scDer.mGlobalCorrdR[side].getNDataPoints()) {
+    LOGP(error, "Global distortions or global corrections of the derivative map not set for side {}: distortions are not scaled", static_cast<int>(side));
+    return false;
+  }
+
+  // the new distortions are written to separate containers, since the current distortions are interpolated during the calculation
+  DataContainer distdR(mParamGrid.NZVertices, mParamGrid.NRVertices, mParamGrid.NPhiVertices);
+  DataContainer distdZ(mParamGrid.NZVertices, mParamGrid.NRVertices, mParamGrid.NPhiVertices);
+  DataContainer distdRPhi(mParamGrid.NZVertices, mParamGrid.NRVertices, mParamGrid.NPhiVertices);
+  unsigned long nNotConverged = 0;
+
+#pragma omp parallel for num_threads(sNThreads) reduction(+ : nNotConverged)
+  for (unsigned int iPhi = 0; iPhi < mParamGrid.NPhiVertices; ++iPhi) {
+    const DataT phi = getPhiVertex(iPhi, side);
+    for (unsigned int iR = 0; iR < mParamGrid.NRVertices; ++iR) {
+      const DataT radius = getRVertex(iR, side);
+      const DataT px = getXFromPolar(radius, phi);
+      const DataT py = getYFromPolar(radius, phi);
+      for (unsigned int iZ = 0; iZ < mParamGrid.NZVertices; ++iZ) {
+        const DataT pz = getZVertex(iZ, side);
+
+        // distorted point y = T(p) of the vertex p, from the distortions at the vertex
+        const DataT radiusDist = radius + mGlobalDistdR[side](iZ, iR, iPhi);
+        const DataT phiDist = phi + mGlobalDistdRPhi[side](iZ, iR, iPhi) / radius;
+        DataT yX = getXFromPolar(radiusDist, phiDist);
+        DataT yY = getYFromPolar(radiusDist, phiDist);
+        DataT yZ = pz + mGlobalDistdZ[side](iZ, iR, iPhi);
+
+        // fixed-point iteration for the start point x = p - scale * C_der(T(x)), starting at x = p
+        DataT xX = px;
+        DataT xY = py;
+        DataT xZ = pz;
+        bool converged = false;
+        for (int iter = 0; iter < maxIter; ++iter) {
+          DataT corrX = 0;
+          DataT corrY = 0;
+          DataT corrZ = 0;
+          scDer.getCorrections(yX, yY, yZ, side, corrX, corrY, corrZ);
+          const DataT xXNew = px - scale * corrX;
+          const DataT xYNew = py - scale * corrY;
+          const DataT xZNew = pz - scale * corrZ;
+          const DataT diff = std::max(std::abs(xXNew - xX), std::max(std::abs(xYNew - xY), std::abs(xZNew - xZ)));
+          xX = xXNew;
+          xY = xYNew;
+          xZ = xZNew;
+
+          // distorted point of the new start point
+          DataT distX = 0;
+          DataT distY = 0;
+          DataT distZ = 0;
+          getDistortions(xX, xY, xZ, side, distX, distY, distZ);
+          yX = xX + distX;
+          yY = xY + distY;
+          yZ = xZ + distZ;
+          if (diff < diffCorr) {
+            converged = true;
+            break;
+          }
+        }
+        nNotConverged += !converged;
+
+        // new distortion of the vertex: y - p in r, r phi (at the radius of the vertex) and z. The angle between p and y is calculated from the cross
+        // and the dot product: the difference of the two angles is close to -2 pi for phi > pi and has to be wrapped, which in single precision gives
+        // an offset (2 pi in single precision is larger than 2 pi by 1.7e-7) and a larger rounding error
+        const DataT deltaPhi = std::atan2(px * yY - py * yX, px * yX + py * yY);
+        distdR(iZ, iR, iPhi) = getRadiusFromCartesian(yX, yY) - radius;
+        distdRPhi(iZ, iR, iPhi) = deltaPhi * radius;
+        distdZ(iZ, iR, iPhi) = yZ - pz;
+      }
+    }
+  }
+
+  // swap the data: the interpolators of the global distortions keep the address of the containers
+  mGlobalDistdR[side].getData().swap(distdR.getData());
+  mGlobalDistdZ[side].getData().swap(distdZ.getData());
+  mGlobalDistdRPhi[side].getData().swap(distdRPhi.getData());
+
+  if (nNotConverged) {
+    LOGP(warning, "Scaling of the global distortions with the derivative map (scale {}): {} of {} vertices on side {} not converged after {} iterations", scale, nNotConverged, mGlobalDistdR[side].getNDataPoints(), static_cast<int>(side), maxIter);
+  }
+  warnNonFinite(mGlobalDistdR[side], mGlobalDistdRPhi[side], mGlobalDistdZ[side], "global distortions scaled with the derivative map", side);
+  return true;
+}
+
+template <typename DataT>
+bool SpaceCharge<DataT>::scaleGlobalDistWithDerivativeCorr(const SpaceCharge<DataT>& scDer, const float scale, const int maxIter, const DataT diffCorr)
+{
+  // the sides are processed one after the other: the calculation for each side is parallelised
+  bool scaled = true;
+  for (int iside = 0; iside < FNSIDES; ++iside) {
+    const o2::tpc::Side side = (iside == 0) ? Side::A : Side::C;
+    scaled &= scaleGlobalDistWithDerivativeCorr(side, scDer, scale, maxIter, diffCorr);
+  }
+  return scaled;
+}
+
 template <typename DataT>
 void SpaceCharge<DataT>::calcGlobalCorrWithGlobalDistIterativeLinearCartesian(const Side side, const SpaceCharge<DataT>* scSCale, float scale, const int maxIter, const DataT approachX, const DataT approachY, const DataT approachZ, const DataT diffCorr)
 {