TPC: Optimising scaling of distortions in MC

- Making simulation of scaled distortions in the digitizer fully
consistent with the scaled corrections in the reconstruction

- cleaning some old code

Author: matthias-kleiner

Detectors/TPC/simulation/src/Digitizer.cxx

@@ -257,7 +257,7 @@ void Digitizer::recalculateDistortions()
       mSpaceChargeDer->calcGlobalCorrWithGlobalDistIterative(side, nullptr, 0);
 
       LOGP(info, "Calculating scaled distortions with scaling factor {}", mLumiScaleFactor);
-      mSpaceCharge->calcGlobalDistWithGlobalCorrIterative(side, mSpaceChargeDer.get(), mLumiScaleFactor);
+      mSpaceCharge->calcGlobalDistWithGlobalCorrIterativeLinearCartesian(side, mSpaceChargeDer.get(), mLumiScaleFactor);
     }
     // set new lumi of avg map
     mSpaceCharge->setMeanLumi(CorrMapParam::Instance().lumiInst);

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

@@ -369,6 +369,20 @@ class SpaceCharge
   void calcGlobalCorrWithGlobalDistIterative(const Side side, const SpaceCharge<DataT>* scSCale = nullptr, float scale = 0, const int maxIter = 100, const DataT approachZ = 1, const DataT approachR = 1, const DataT approachPhi = 1, const DataT diffCorr = 50e-6);
   void calcGlobalCorrWithGlobalDistIterative(const SpaceCharge<DataT>* scSCale = nullptr, float scale = 0, const int maxIter = 100, const DataT approachZ = 1, const DataT approachR = 1, const DataT approachPhi = 1, const DataT diffCorr = 50e-6);
 
+  /// Calculate global distortions using the global corrections by scaling with second distortion object, which will be consistent with scaled corrections in cartesian coordinates (as done in the tracking)
+  /// \param side Side of the TPC
+  /// \param scSCale possible second sc object
+  /// \param scale scaling for second sc object
+  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);
+
+  /// 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
+  /// \param scale scaling for second sc object
+  void calcGlobalCorrWithGlobalDistIterativeLinearCartesian(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 calcGlobalCorrWithGlobalDistIterativeLinearCartesian(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);
+
   /// \return returns number of vertices in z direction
   unsigned short getNZVertices() const { return mParamGrid.NZVertices; }
 
@@ -1373,6 +1387,7 @@ class SpaceCharge
   void initRodAlignmentVoltages(const MisalignmentType misalignmentType, const FCType fcType, const int sector, const Side side, const float deltaPot);
 
   void calcGlobalDistCorrIterative(const DistCorrInterpolator<DataT>& globCorr, const int maxIter, const DataT approachZ, const DataT approachR, const DataT approachPhi, const DataT diffCorr, const SpaceCharge<DataT>* scSCale, float scale, const Type type);
+  void calcGlobalDistCorrIterativeLinearCartesian(const DistCorrInterpolator<DataT>& globCorr, const int maxIter, const DataT approachX, const DataT approachY, const DataT approachZ, const DataT diffCorr, const SpaceCharge<DataT>* scSCale, float scale, const Type type);
 
   ClassDefNV(SpaceCharge, 6);
 };

Detectors/TPC/spacecharge/src/SpaceCharge.cxx

@@ -757,15 +757,6 @@ void SpaceCharge<DataT>::calcGlobalDistCorrIterative(const DistCorrInterpolator<
       const DataT radius = getRVertex(iR, side);
       for (unsigned int iZ = 1; iZ < mParamGrid.NZVertices; ++iZ) {
         const DataT z = getZVertex(iZ, side);
-
-        unsigned int nearestiZ = iZ;
-        unsigned int nearestiR = iR;
-        unsigned int nearestiPhi = iPhi;
-
-        DataT nearestZ = getZVertex(nearestiZ, side);
-        DataT nearestR = getRVertex(nearestiR, side);
-        DataT nearestPhi = getPhiVertex(nearestiPhi, side);
-
         //
         //==========================================================================================
         //==== start algorithm: use tricubic upsampling to numerically approach the query point ====
@@ -774,9 +765,9 @@ void SpaceCharge<DataT>::calcGlobalDistCorrIterative(const DistCorrInterpolator<
         // 1. calculate difference from nearest point to query point with stepwidth factor x
         // and approach the new point
         //
-        DataT stepR = (radius - nearestR) * approachR;
-        DataT stepZ = (z - nearestZ) * approachZ;
-        DataT stepPhi = (phi - nearestPhi) * approachPhi;
+        DataT stepR = 0;
+        DataT stepZ = 0;
+        DataT stepPhi = 0;
 
         // needed to check for convergence
         DataT lastCorrdR = std::numeric_limits<DataT>::max();
@@ -790,9 +781,9 @@ void SpaceCharge<DataT>::calcGlobalDistCorrIterative(const DistCorrInterpolator<
 
         for (int iter = 0; iter < maxIter; ++iter) {
           // 2. get new point coordinates
-          const DataT rCurrPos = getRVertex(nearestiR, side) + stepR;
-          const DataT zCurrPos = getZVertex(nearestiZ, side) + stepZ;
-          const DataT phiCurrPos = getPhiVertex(nearestiPhi, side) + stepPhi;
+          const DataT rCurrPos = radius + stepR;
+          const DataT zCurrPos = z + stepZ;
+          const DataT phiCurrPos = phi + stepPhi;
 
           // abort calculation of drift path if electron reached inner/outer field cage or central electrode
           if (rCurrPos <= getRMinSim(side) || rCurrPos >= getRMaxSim(side) || getSide(zCurrPos) != side) {
@@ -867,6 +858,177 @@ void SpaceCharge<DataT>::calcGlobalDistCorrIterative(const DistCorrInterpolator<
   }
 }
 
+template <typename DataT>
+void SpaceCharge<DataT>::calcGlobalDistWithGlobalCorrIterativeLinearCartesian(const Side side, const SpaceCharge<DataT>* scSCale, float scale, const int maxIter, const DataT approachX, const DataT approachY, const DataT approachZ, const DataT diffCorr)
+{
+  calcGlobalDistCorrIterativeLinearCartesian(getGlobalCorrInterpolator(side), maxIter, approachX, approachY, approachZ, diffCorr, scSCale, scale, Type::Distortions);
+}
+
+template <typename DataT>
+void SpaceCharge<DataT>::calcGlobalDistWithGlobalCorrIterativeLinearCartesian(const SpaceCharge<DataT>* scSCale, float scale, const int maxIter, const DataT approachX, const DataT approachY, const DataT approachZ, const DataT diffCorr)
+{
+#pragma omp parallel for num_threads(sNThreads)
+  for (int iside = 0; iside < FNSIDES; ++iside) {
+    const o2::tpc::Side side = (iside == 0) ? Side::A : Side::C;
+    calcGlobalDistWithGlobalCorrIterativeLinearCartesian(side, scSCale, scale, maxIter, approachX, approachY, approachZ, diffCorr);
+  }
+}
+
+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)
+{
+  calcGlobalDistCorrIterativeLinearCartesian(getGlobalDistInterpolator(side), maxIter, approachX, approachY, approachZ, diffCorr, scSCale, scale, Type::Corrections);
+}
+
+template <typename DataT>
+void SpaceCharge<DataT>::calcGlobalCorrWithGlobalDistIterativeLinearCartesian(const SpaceCharge<DataT>* scSCale, float scale, const int maxIter, const DataT approachX, const DataT approachY, const DataT approachZ, const DataT diffCorr)
+{
+#pragma omp parallel for num_threads(sNThreads)
+  for (int iside = 0; iside < FNSIDES; ++iside) {
+    const o2::tpc::Side side = (iside == 0) ? Side::A : Side::C;
+    calcGlobalCorrWithGlobalDistIterativeLinearCartesian(side, scSCale, scale, maxIter, approachX, approachY, approachZ, diffCorr);
+  }
+}
+
+template <typename DataT>
+void SpaceCharge<DataT>::calcGlobalDistCorrIterativeLinearCartesian(const DistCorrInterpolator<DataT>& globCorr, const int maxIter, const DataT approachX, const DataT approachY, const DataT approachZ, const DataT diffCorr, const SpaceCharge<DataT>* scSCale, float scale, const Type type)
+{
+  const Side side = globCorr.getSide();
+  if (type == Type::Distortions) {
+    initContainer(mGlobalDistdR[side], true);
+    initContainer(mGlobalDistdZ[side], true);
+    initContainer(mGlobalDistdRPhi[side], true);
+  } else {
+    initContainer(mGlobalCorrdR[side], true);
+    initContainer(mGlobalCorrdZ[side], true);
+    initContainer(mGlobalCorrdRPhi[side], true);
+  }
+
+  const auto& scSCaleInterpolator = (type == Type::Distortions) ? scSCale->mInterpolatorGlobalCorr[side] : scSCale->mInterpolatorGlobalDist[side];
+
+#pragma omp parallel for num_threads(sNThreads)
+  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 x = getXFromPolar(radius, phi);
+      const DataT y = getYFromPolar(radius, phi);
+
+      for (unsigned int iZ = 1; iZ < mParamGrid.NZVertices; ++iZ) {
+        const DataT z = getZVertex(iZ, side);
+
+        DataT stepX = 0;
+        DataT stepY = 0;
+        DataT stepZ = 0;
+
+        // needed to check for convergence
+        DataT lastCorrX = std::numeric_limits<DataT>::max();
+        DataT lastCorrY = std::numeric_limits<DataT>::max();
+        DataT lastCorrZ = std::numeric_limits<DataT>::max();
+        DataT lastX = std::numeric_limits<DataT>::max();
+        DataT lastY = std::numeric_limits<DataT>::max();
+        DataT lastZ = std::numeric_limits<DataT>::max();
+
+        for (int iter = 0; iter < maxIter; ++iter) {
+          const DataT xCurrPos = x + stepX;
+          const DataT yCurrPos = y + stepY;
+          const DataT zCurrPos = z + stepZ;
+
+          // abort calculation of drift path if electron reached inner/outer field cage or central electrode
+          const DataT rCurrPos = getRadiusFromCartesian(xCurrPos, yCurrPos);
+          if (rCurrPos <= getRMinSim(side) || rCurrPos >= getRMaxSim(side) || getSide(zCurrPos) != side) {
+            break;
+          }
+
+          // interpolate global correction at new point and calculate position of global correction
+          DataT corrX = 0;
+          DataT corrY = 0;
+          DataT corrZ = 0;
+          (type == Type::Distortions) ? getCorrections(xCurrPos, yCurrPos, zCurrPos, side, corrX, corrY, corrZ) : getDistortions(xCurrPos, yCurrPos, zCurrPos, side, corrX, corrY, corrZ);
+
+          if (scSCale && scale != 0) {
+            DataT corrXScale = 0;
+            DataT corrYScale = 0;
+            DataT corrZScale = 0;
+            (type == Type::Distortions) ? scSCale->getCorrections(xCurrPos, yCurrPos, zCurrPos, side, corrXScale, corrYScale, corrZScale) : getDistortions(xCurrPos, yCurrPos, zCurrPos, side, corrX, corrY, corrZ);
+            corrX += scale * corrXScale;
+            corrY += scale * corrYScale;
+            corrZ += scale * corrZScale;
+          }
+
+          // check for convergence
+          const DataT diffCorrX = std::abs(corrX - lastCorrX);
+          const DataT diffCorrY = std::abs(corrY - lastCorrY);
+          const DataT diffCorrZ = std::abs(corrZ - lastCorrZ);
+
+          lastCorrX = corrX;
+          lastCorrY = corrY;
+          lastCorrZ = corrZ;
+          lastX = xCurrPos;
+          lastY = yCurrPos;
+          lastZ = zCurrPos;
+
+          // stop algorithm if converged
+          if ((diffCorrX < diffCorr) && (diffCorrY < diffCorr) && (diffCorrZ < diffCorr)) {
+            break;
+          }
+
+          const DataT xNewPos = xCurrPos + corrX;
+          const DataT yNewPos = yCurrPos + corrY;
+          const DataT zNewPos = zCurrPos + corrZ;
+
+          // approach desired coordinate
+          stepX += (x - xNewPos) * approachX;
+          stepY += (y - yNewPos) * approachY;
+          stepZ += (z - zNewPos) * approachZ;
+        }
+
+        const DataT xNew = lastX + lastCorrX;
+        const DataT yNew = lastY + lastCorrY;
+        const DataT radiusNew = getRadiusFromCartesian(xNew, yNew);
+        const DataT corrdR = -(radiusNew - getRadiusFromCartesian(lastX, lastY));
+
+        float phiNew = getPhiFromCartesian(xNew, yNew);
+        o2::math_utils::bringTo02PiGen(phiNew);
+
+        float phiLast = getPhiFromCartesian(lastX, lastY);
+        o2::math_utils::bringTo02PiGen(phiLast);
+
+        DataT deltaPhi = (phiNew - phiLast);
+        // handle edge cases
+        if (deltaPhi > PI) {
+          deltaPhi -= 2 * PI;
+        } else if (deltaPhi < -PI) {
+          deltaPhi += 2 * PI;
+        }
+        const DataT corrdRPhi = -deltaPhi * radiusNew;
+
+        // set global distortions if algorithm converged or iterations exceed max numbers of iterations
+        if (type == Type::Distortions) {
+          mGlobalDistdR[side](iZ, iR, iPhi) = corrdR;
+          mGlobalDistdRPhi[side](iZ, iR, iPhi) = corrdRPhi;
+          mGlobalDistdZ[side](iZ, iR, iPhi) = -lastCorrZ;
+        } else {
+          mGlobalCorrdR[side](iZ, iR, iPhi) = corrdR;
+          mGlobalCorrdRPhi[side](iZ, iR, iPhi) = corrdRPhi;
+          mGlobalCorrdZ[side](iZ, iR, iPhi) = -lastCorrZ;
+        }
+      }
+    }
+    for (unsigned int iR = 0; iR < mParamGrid.NRVertices; ++iR) {
+      if (type == Type::Distortions) {
+        mGlobalDistdR[side](0, iR, iPhi) = 3 * (mGlobalDistdR[side](1, iR, iPhi) - mGlobalDistdR[side](2, iR, iPhi)) + mGlobalDistdR[side](3, iR, iPhi);
+        mGlobalDistdRPhi[side](0, iR, iPhi) = 3 * (mGlobalDistdRPhi[side](1, iR, iPhi) - mGlobalDistdRPhi[side](2, iR, iPhi)) + mGlobalDistdRPhi[side](3, iR, iPhi);
+        mGlobalDistdZ[side](0, iR, iPhi) = 3 * (mGlobalDistdZ[side](1, iR, iPhi) - mGlobalDistdZ[side](2, iR, iPhi)) + mGlobalDistdZ[side](3, iR, iPhi);
+      } else {
+        mGlobalCorrdR[side](0, iR, iPhi) = 3 * (mGlobalCorrdR[side](1, iR, iPhi) - mGlobalCorrdR[side](2, iR, iPhi)) + mGlobalCorrdR[side](3, iR, iPhi);
+        mGlobalCorrdRPhi[side](0, iR, iPhi) = 3 * (mGlobalCorrdRPhi[side](1, iR, iPhi) - mGlobalCorrdRPhi[side](2, iR, iPhi)) + mGlobalCorrdRPhi[side](3, iR, iPhi);
+        mGlobalCorrdZ[side](0, iR, iPhi) = 3 * (mGlobalCorrdZ[side](1, iR, iPhi) - mGlobalCorrdZ[side](2, iR, iPhi)) + mGlobalCorrdZ[side](3, iR, iPhi);
+      }
+    }
+  }
+}
+
 template <typename DataT>
 NumericalFields<DataT> SpaceCharge<DataT>::getElectricFieldsInterpolator(const Side side) const
 {