TPC: Adding loading of distortions maps from CCDB for MC (#12300)

* TPC: Adding loading of distortions maps from CCDB for MC

Distortions from pre-calculated maps can be loaded from CCDB with the option
--tpc-simulate-distortions in the o2-sim-digitizer-workflow:

o2-sim-digitizer-workflow $ARGS_ALL -b --tpc-lanes 1 --onlyDet TPC --tpc-simulate-distortions

* Implement scaling of distortions in MC

- Scaling of the distortions can be enabled with `o2-sim-digitizer-workflow --tpc-distortion-type 2`, where the lumi is taken from `TPCCorrMap.lumiInst`.
- space charge: adding meta data for scaling
- In reconstruction the --corrmap-lumi-mode 2 should be used for MC and
luminosity set via TPCCorrMap.lumiInst

- Fixing inverse transform for scalMode 1 and 2!

- other small changes

* renaming CCDB path for MC correction map to be more consistent

* Adding lumi source as meta data

Author: matthias-kleiner

Detectors/TPC/base/include/TPCBase/CDBTypes.h

@@ -71,13 +71,18 @@ enum class CDBType {
                        ///
   CalCorrMap,          ///< Cluster correction map (high IR rate distortions)
   CalCorrMapRef,       ///< Cluster correction reference map (static distortions)
+  CalCorrMapMC,        ///< Cluster correction map (high IR rate distortions) for MC
+  CalCorrDerivMapMC,   ///< Cluster correction reference map (static distortions) for MC
                        ///
   CalCorrDerivMap,     ///< Cluster correction map (derivative map)
                        ///
   CalTimeSeries,       ///< integrated DCAs for longer time interval
   CalScaler,           ///< Scaler from IDCs or combined estimator
                        ///
   CorrMapParam,        ///< parameters for CorrectionMapsLoader configuration
+                       ///
+  DistortionMapMC,     ///< full distortions (static + IR dependant) for MC used in the digitizer
+  DistortionMapDerivMC ///< derivative distortions for MC used in the digitizer for scaling
 };
 
 /// Storage name in CCDB for each calibration and parameter type
@@ -129,13 +134,19 @@ const std::unordered_map<CDBType, const std::string> CDBTypeMap{
   // correction maps
   {CDBType::CalCorrMap, "TPC/Calib/CorrectionMapV2"},
   {CDBType::CalCorrMapRef, "TPC/Calib/CorrectionMapRefV2"},
+  // correction maps for MC
+  {CDBType::CalCorrMapMC, "TPC/Calib/CorrectionMapMCV2"},
+  {CDBType::CalCorrDerivMapMC, "TPC/Calib/CorrectionMapDerivativeMCV2"},
   // derivative map correction
   {CDBType::CalCorrDerivMap, "TPC/Calib/CorrectionMapDerivativeV2"},
   // time series
   {CDBType::CalTimeSeries, "TPC/Calib/TimeSeries"},
   {CDBType::CalScaler, "TPC/Calib/Scaler"},
   // correction maps loader params
   {CDBType::CorrMapParam, "TPC/Calib/CorrMapParam"},
+  // distortion maps
+  {CDBType::DistortionMapMC, "TPC/Calib/DistortionMapMC"},
+  {CDBType::DistortionMapDerivMC, "TPC/Calib/DistortionMapDerivativeMC"},
 };
 
 } // namespace o2::tpc

Detectors/TPC/calibration/src/CorrectionMapsLoader.cxx

@@ -67,11 +67,16 @@ void CorrectionMapsLoader::extractCCDBInputs(ProcessingContext& pc)
 //________________________________________________________
 void CorrectionMapsLoader::requestCCDBInputs(std::vector<InputSpec>& inputs, std::vector<o2::framework::ConfigParamSpec>& options, const CorrectionMapsLoaderGloOpts& gloOpts)
 {
-  addInput(inputs, {"tpcCorrMap", "TPC", "CorrMap", 0, Lifetime::Condition, ccdbParamSpec(CDBTypeMap.at(CDBType::CalCorrMap), {}, 1)}); // time-dependent
   if (gloOpts.lumiMode == 0) {
+    addInput(inputs, {"tpcCorrMap", "TPC", "CorrMap", 0, Lifetime::Condition, ccdbParamSpec(CDBTypeMap.at(CDBType::CalCorrMap), {}, 1)});          // time-dependent
     addInput(inputs, {"tpcCorrMapRef", "TPC", "CorrMapRef", 0, Lifetime::Condition, ccdbParamSpec(CDBTypeMap.at(CDBType::CalCorrMapRef), {}, 0)}); // load once
   } else if (gloOpts.lumiMode == 1) {
+    addInput(inputs, {"tpcCorrMap", "TPC", "CorrMap", 0, Lifetime::Condition, ccdbParamSpec(CDBTypeMap.at(CDBType::CalCorrMap), {}, 1)});            // time-dependent
     addInput(inputs, {"tpcCorrMapRef", "TPC", "CorrMapRef", 0, Lifetime::Condition, ccdbParamSpec(CDBTypeMap.at(CDBType::CalCorrDerivMap), {}, 1)}); // time-dependent
+  } else if (gloOpts.lumiMode == 2) {
+    // for MC corrections
+    addInput(inputs, {"tpcCorrMap", "TPC", "CorrMap", 0, Lifetime::Condition, ccdbParamSpec(CDBTypeMap.at(CDBType::CalCorrMapMC), {}, 1)});            // time-dependent
+    addInput(inputs, {"tpcCorrMapRef", "TPC", "CorrMapRef", 0, Lifetime::Condition, ccdbParamSpec(CDBTypeMap.at(CDBType::CalCorrDerivMapMC), {}, 1)}); // time-dependent
   } else {
     LOG(fatal) << "Correction mode unknown! Choose either 0 (default) or 1 (derivative map) for flag corrmap-lumi-mode.";
   }
@@ -99,7 +104,7 @@ void CorrectionMapsLoader::addGlobalOptions(std::vector<ConfigParamSpec>& option
 {
   // these are options which should be added at the workflow level, since they modify the inputs of the devices
   addOption(options, ConfigParamSpec{"lumi-type", o2::framework::VariantType::Int, 0, {"1 = require CTP lumi for TPC correction scaling, 2 = require TPC scalers for TPC correction scaling"}});
-  addOption(options, ConfigParamSpec{"corrmap-lumi-mode", o2::framework::VariantType::Int, 0, {"scaling mode: (default) 0 = static + scale * full; 1 = full + scale * derivative"}});
+  addOption(options, ConfigParamSpec{"corrmap-lumi-mode", o2::framework::VariantType::Int, 0, {"scaling mode: (default) 0 = static + scale * full; 1 = full + scale * derivative; 2 = full + scale * derivative (for MC)"}});
 }
 
 //________________________________________________________

Detectors/TPC/simulation/CMakeLists.txt

@@ -22,7 +22,7 @@ o2_add_library(TPCSimulation
                        src/SAMPAProcessing.cxx
                        src/IDCSim.cxx
                PUBLIC_LINK_LIBRARIES O2::DetectorsBase O2::SimulationDataFormat
-                                     O2::TPCBase O2::TPCSpaceCharge
+                                     O2::TPCBase O2::TPCSpaceCharge O2::TPCCalibration
                                      ROOT::Physics)
 
 o2_target_root_dictionary(TPCSimulation

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

@@ -53,7 +53,7 @@ enum class SCDistortionType : int;
 class Digitizer
 {
  public:
-  using SC = SpaceCharge<double>;
+  using SC = SpaceCharge<float>;
 
   /// Default constructor
   Digitizer();
@@ -121,6 +121,9 @@ class Digitizer
   /// \param spaceCharge unique pointer to spaceCharge object
   void setUseSCDistortions(SC* spaceCharge);
 
+  /// \param spaceCharge unique pointer to spaceCharge object
+  void setSCDistortionsDerivative(SC* spaceCharge);
+
   /// Enable the use of space-charge distortions by providing global distortions and global corrections stored in a ROOT file
   /// The storage of the values should be done by the methods provided in the SpaceCharge class
   /// \param file containing distortions
@@ -129,17 +132,26 @@ class Digitizer
   void setVDrift(float v) { mVDrift = v; }
   void setTDriftOffset(float t) { mTDriftOffset = t; }
 
+  void setDistortionScaleType(int distortionScaleType) { mDistortionScaleType = distortionScaleType; }
+  int getDistortionScaleType() const { return mDistortionScaleType; }
+  void setLumiScaleFactor();
+  void setMeanLumiDistortions(float meanLumi);
+  void setMeanLumiDistortionsDerivative(float meanLumi);
+
  private:
-  DigitContainer mDigitContainer;    ///< Container for the Digits
-  std::unique_ptr<SC> mSpaceCharge;  ///< Handler of space-charge distortions
-  Sector mSector = -1;               ///< ID of the currently processed sector
-  double mEventTime = 0.f;           ///< Time of the currently processed event
-  double mOutputDigitTimeOffset = 0; ///< Time of the first IR sampled in the digitizer
-  float mVDrift = 0;                 ///< VDrift for current timestamp
-  float mTDriftOffset = 0;           ///< drift time additive offset in \mus
-  bool mIsContinuous;                ///< Switch for continuous readout
-  bool mUseSCDistortions = false;    ///< Flag to switch on the use of space-charge distortions
-  ClassDefNV(Digitizer, 1);
+  DigitContainer mDigitContainer;      ///< Container for the Digits
+  std::unique_ptr<SC> mSpaceCharge;    ///< Handler of full distortions (static + IR dependant)
+  std::unique_ptr<SC> mSpaceChargeDer; ///< Handler of reference static distortions
+  Sector mSector = -1;                 ///< ID of the currently processed sector
+  double mEventTime = 0.f;             ///< Time of the currently processed event
+  double mOutputDigitTimeOffset = 0;   ///< Time of the first IR sampled in the digitizer
+  float mVDrift = 0;                   ///< VDrift for current timestamp
+  float mTDriftOffset = 0;             ///< drift time additive offset in \mus
+  bool mIsContinuous;                  ///< Switch for continuous readout
+  bool mUseSCDistortions = false;      ///< Flag to switch on the use of space-charge distortions
+  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
+  ClassDefNV(Digitizer, 2);
 };
 } // namespace tpc
 } // namespace o2

Detectors/TPC/simulation/src/Digitizer.cxx

@@ -27,6 +27,7 @@
 #include "TPCBase/CDBInterface.h"
 #include "TPCSpaceCharge/SpaceCharge.h"
 #include "TPCBase/Mapper.h"
+#include "TPCCalibration/CorrMapParam.h"
 
 #include <fairlogger/Logger.h>
 
@@ -40,10 +41,6 @@ Digitizer::Digitizer() = default;
 
 void Digitizer::init()
 {
-  // Calculate distortion lookup tables if initial space-charge density is provided
-  if (mUseSCDistortions) {
-    mSpaceCharge->init();
-  }
   auto& gemAmplification = GEMAmplification::instance();
   gemAmplification.updateParameters();
   auto& electronTransport = ElectronTransport::instance();
@@ -83,8 +80,10 @@ void Digitizer::process(const std::vector<o2::tpc::HitGroup>& hits,
       GlobalPosition3D posEle(eh.GetX(), eh.GetY(), eh.GetZ());
 
       // Distort the electron position in case space-charge distortions are used
-      if (mUseSCDistortions) {
+      if (mDistortionScaleType == 1) {
         mSpaceCharge->distortElectron(posEle);
+      } else if (mDistortionScaleType == 2) {
+        mSpaceCharge->distortElectron(posEle, mSpaceChargeDer.get(), mLumiScaleFactor);
       }
 
       /// Remove electrons that end up more than three sigma of the hit's average diffusion away from the current sector
@@ -190,6 +189,15 @@ void Digitizer::setUseSCDistortions(SC* spaceCharge)
 {
   mUseSCDistortions = true;
   mSpaceCharge.reset(spaceCharge);
+  mSpaceCharge->initAfterReadingFromFile();
+  mSpaceCharge->printMetaData();
+}
+
+void Digitizer::setSCDistortionsDerivative(SC* spaceCharge)
+{
+  mSpaceChargeDer.reset(spaceCharge);
+  mSpaceChargeDer->initAfterReadingFromFile();
+  mSpaceChargeDer->printMetaData();
 }
 
 void Digitizer::setUseSCDistortions(std::string_view finp)
@@ -213,3 +221,19 @@ void Digitizer::setStartTime(double time)
   sampaProcessing.updateParameters(mVDrift);
   mDigitContainer.setStartTime(sampaProcessing.getTimeBinFromTime(time - mOutputDigitTimeOffset));
 }
+
+void Digitizer::setLumiScaleFactor()
+{
+  mLumiScaleFactor = (CorrMapParam::Instance().lumiInst - mSpaceCharge->getMeanLumi()) / mSpaceChargeDer->getMeanLumi();
+  LOGP(info, "Setting Lumi scale factor: lumiInst: {}  lumi mean: {} lumi mean derivative: {} lumi scale factor: {}", CorrMapParam::Instance().lumiInst, mSpaceCharge->getMeanLumi(), mSpaceChargeDer->getMeanLumi(), mLumiScaleFactor);
+}
+
+void Digitizer::setMeanLumiDistortions(float meanLumi)
+{
+  mSpaceCharge->setMeanLumi(meanLumi);
+}
+
+void Digitizer::setMeanLumiDistortionsDerivative(float meanLumi)
+{
+  mSpaceChargeDer->setMeanLumi(meanLumi);
+}

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

@@ -572,7 +572,9 @@ class SpaceCharge
 
   /// Distort electron position using distortion lookup tables
   /// \param point 3D coordinates of the electron
-  void distortElectron(GlobalPosition3D& point) const;
+  /// \param scSCale other sc object which is used for scaling of the distortions
+  /// \param scale scaling value
+  void distortElectron(GlobalPosition3D& point, const SpaceCharge<DataT>* scSCale = nullptr, float scale = 0) const;
 
   /// set the distortions directly from a look up table
   /// \param distdZ distortions in z direction
@@ -1170,12 +1172,24 @@ class SpaceCharge
   /// substract global corrections from other sc object (global corrections -= other.global corrections)
   /// can be used to calculate the derivative: (this - other)/normalization
   /// for normalization see scaleCorrections()
-  void substractGlobalCorrections(const SpaceCharge<DataT>& otherSC, const Side side);
+  void subtractGlobalCorrections(const SpaceCharge<DataT>& otherSC, const Side side);
+
+  /// substract global distortions from other sc object (global distortions -= other.global distortions)
+  /// can be used to calculate the derivative: (this - other)/normalization
+  void subtractGlobalDistortions(const SpaceCharge<DataT>& otherSC, const Side side);
 
   /// scale corrections by factor
   /// \param scaleFac global corrections are multiplied by this factor
   void scaleCorrections(const float scaleFac, const Side side);
 
+  /// setting meta data for this object
+  void setMetaData(const SCMetaData& meta) { mMeta = meta; }
+  const auto& getMetaData() const { return mMeta; }
+  void printMetaData() const { mMeta.print(); }
+  float getMeanLumi() const { return mMeta.meanLumi; }
+  void setMeanLumi(float lumi) { mMeta.meanLumi = lumi; }
+  void initAfterReadingFromFile();
+
  private:
   ParamSpaceCharge mParamGrid{};                                                                          ///< parameters of the grid on which the calculations are performed
   inline static int sNThreads{getOMPMaxThreads()};                                                        ///<! number of threads which are used during the calculations
@@ -1229,6 +1243,7 @@ class SpaceCharge
   AnalyticalDistCorr<DataT> mAnaDistCorr;                                                                                                                                                                                                                                                     ///< analytical distortions and corrections
   bool mUseAnaDistCorr{false};                                                                                                                                                                                                                                                                ///< flag if analytical distortions will be used in the distortElectron() and getCorrections() function
   BField mBField{};                                                                                                                                                                                                                                                                           ///<! B-Field                                                                                                                                                                                                                                                                                                                                                                                                                                                                                     ///<! B field
+  SCMetaData mMeta{};                                                                                                                                                                                                                                                                         ///< meta data
 
   /// check if the addition of two values are close to zero.
   /// This avoids errors during the integration of the electric fields when the sum of the nominal electric with the electric field from the space charge is close to 0 (usually this is not the case!).
@@ -1330,7 +1345,7 @@ class SpaceCharge
   /// set potentialsdue to ROD misalignment
   void initRodAlignmentVoltages(const MisalignmentType misalignmentType, const FCType fcType, const int sector, const Side side, const float deltaPot);
 
-  ClassDefNV(SpaceCharge, 5);
+  ClassDefNV(SpaceCharge, 6);
 };
 
 } // namespace tpc

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

@@ -23,12 +23,41 @@
 #include "TPCSpaceCharge/TriCubic.h"
 #include "DataFormatsTPC/Defs.h"
 #include "TFormula.h"
+#include "Framework/Logger.h"
 
 namespace o2
 {
 namespace tpc
 {
 
+struct SCMetaData {
+
+  enum LumiType {
+    CTP = 0,
+    IDC = 1
+  };
+
+  void print() const
+  {
+    const std::array<std::string, 2> collisionTypes{"PP", "Pb-Pb"};
+    const std::array<std::string, 2> slumiSource{"CTP", "IDC"};
+    if (collisionType < collisionTypes.size()) {
+      LOGP(info, "meanLumi: {}, source of lumi: {}, IR: {}kHz, run: {}, collisionType {}", meanLumi, slumiSource[lumiSource], ir, run, collisionTypes[collisionType]);
+    } else {
+      LOGP(info, "Specified collision type {} not allowed", collisionType);
+    }
+  }
+
+  float meanLumi = 0;    ///< mean lumi the sc object corresponds to
+  float ir = 0;          ///< IR
+  int run = 0;           ///< run number this object anchored to to
+  int collisionType = 0; ///< 0=PP, 1-Pb-Pb
+  LumiType lumiSource{}; ///< source of luminosity
+
+ private:
+  ClassDefNV(SCMetaData, 1);
+};
+
 ///
 /// this class contains an analytical description of the space charge, potential and the electric fields.
 /// The analytical functions can be used to test the poisson solver and the caluclation of distortions/corrections.