TPC SC: faster and more robust space-charge distortion calculation

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

@@ -140,6 +140,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();
 
@@ -192,7 +199,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, 4);
+  int mRecalcDistortionsMethod = 1;      ///< method of recalculateDistortions(): 0 = inversion of the scaled corrections, 1 = derivative corrections applied to the distortions
+  ClassDefNV(Digitizer, 5);
 };
 } // namespace tpc
 } // namespace o2

Detectors/TPC/simulation/src/Digitizer.cxx

@@ -351,6 +351,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) {
@@ -362,15 +370,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, 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, 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/CMakeLists.txt

@@ -17,7 +17,6 @@ o2_add_library(TPCSpaceCharge
                        src/DataContainer3D.cxx
                PUBLIC_LINK_LIBRARIES O2::TPCBaseRecSim
                                      O2::Field
-                                     Vc::Vc
                                      ROOT::Core
                                      O2::DataFormatsParameters)
 
@@ -35,7 +34,6 @@ o2_target_root_dictionary(TPCSpaceCharge
                                   include/TPCSpaceCharge/PoissonSolverHelpers.h
                                   include/TPCSpaceCharge/SpaceChargeHelpers.h
                                   include/TPCSpaceCharge/TriCubic.h
-                                  include/TPCSpaceCharge/Vector.h
                                   include/TPCSpaceCharge/Vector3D.h
                                   include/TPCSpaceCharge/SpaceChargeParameter.h
                           LINKDEF src/TPCSpacechargeLinkDef.h)
@@ -49,7 +47,7 @@ o2_add_test_root_macro(macro/createResidualDistortionObject.C
                       PUBLIC_LINK_LIBRARIES O2::TPCSpaceCharge
                                             O2::TPCSpaceChargeIO
                                             O2::CommonUtils
-                      LABELS tpc)
+                      LABELS tpc COMPILE_ONLY)
 
 
 install(FILES macro/createSCHistosFromHits.C
@@ -81,6 +79,14 @@ o2_add_test(TriCubic
             LABELS tpc
             CONFIGURATIONS RelWithDebInfo Release MinSizeRel)
 
+o2_add_test(SpaceCharge
+            COMPONENT_NAME spacecharge
+            PUBLIC_LINK_LIBRARIES O2::TPCSpaceCharge
+            SOURCES test/testO2TPCSpaceCharge.cxx
+            ENVIRONMENT O2_ROOT=${CMAKE_BINARY_DIR}/stage
+            LABELS tpc
+            CONFIGURATIONS RelWithDebInfo Release MinSizeRel)
+
 if (OpenMP_CXX_FOUND)
     target_compile_definitions(${targetName} PRIVATE WITH_OPENMP)
     target_link_libraries(${targetName} PRIVATE OpenMP::OpenMP_CXX)

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

@@ -20,6 +20,8 @@
 
 #include "Rtypes.h"
 #include <vector>
+#include <string_view>
+#include <utility>
 
 class TFile;
 class TTree;