TPC: Fix usage of derivative map for scaling of space-charge corrections

- add option to override lumi of reference map (reference map can be the
derivative map)

- changing definition of scaling when using the derivative map: The
derivative map is defined as "map_high_IR - map_low_IR".
The stored lumi of the derivative map should correspond to
"<Lumi_high_IR> - <Lumi_low_IR>"

- Changing condition when scaling is used: for the derivative map the
scaling is negative, when the inst_lumi < mean_lumi

Detectors/TPC/calibration/src/CorrectionMapsLoader.cxx

@@ -87,7 +87,7 @@ void CorrectionMapsLoader::addOptions(std::vector<ConfigParamSpec>& options)
 {
   addOption(options, ConfigParamSpec{"corrmap-lumi-mean", VariantType::Float, 0.f, {"override TPC corr.map mean lumi (if > 0), disable corrections if < 0"}});
   addOption(options, ConfigParamSpec{"corrmap-lumi-inst", VariantType::Float, 0.f, {"override instantaneous CTP lumi (if > 0) for TPC corr.map scaling, disable corrections if < 0"}});
-  addOption(options, ConfigParamSpec{"corrmap-lumi-mode", VariantType::Int, 0, {"scaling mode: (default) 0 = static + scale * full; 1 = full + scale * derivative"}});
+  addOption(options, ConfigParamSpec{"corrmap-lumi-ref", VariantType::Float, 0.f, {"override TPC corr.mapRef mean lumi (if > 0)"}});
   addOption(options, ConfigParamSpec{"ctp-lumi-factor", VariantType::Float, 1.0f, {"scaling to apply to instantaneous lumi from CTP (but not corrmap-lumi-inst)"}});
   addOption(options, ConfigParamSpec{"ctp-lumi-source", VariantType::Int, 0, {"CTP lumi source: 0 = LumiInfo.getLumi(), 1 = LumiInfo.getLumiAlt()"}});
 }
@@ -124,6 +124,10 @@ bool CorrectionMapsLoader::accountCCDBInputs(const ConcreteDataMatcher& matcher,
   if (matcher == ConcreteDataMatcher("TPC", "CorrMapRef", 0)) {
     setCorrMapRef((o2::gpu::TPCFastTransform*)obj);
     mCorrMapRef->rectifyAfterReadingFromFile();
+    if (getMeanLumiRefOverride() == 0 && mCorrMapRef->getLumi() > 0.) {
+      setMeanLumiRef(mCorrMapRef->getLumi());
+    }
+    LOGP(debug, "MeanLumiRefOverride={} MeanLumiMap={} -> meanLumi = {}", getMeanLumiRefOverride(), mCorrMapRef->getLumi(), getMeanLumiRef());
     setUpdatedMapRef();
     return true;
   }
@@ -143,13 +147,16 @@ void CorrectionMapsLoader::init(o2::framework::InitContext& ic)
     }
   }
   mMeanLumiOverride = ic.options().get<float>("corrmap-lumi-mean");
+  mMeanLumiRefOverride = ic.options().get<float>("corrmap-lumi-ref");
   mInstLumiOverride = ic.options().get<float>("corrmap-lumi-inst");
-  mLumiScaleMode = ic.options().get<int>("corrmap-lumi-mode");
   mInstLumiFactor = ic.options().get<float>("ctp-lumi-factor");
   mCTPLumiSource = ic.options().get<int>("ctp-lumi-source");
   if (mMeanLumiOverride != 0.) {
     setMeanLumi(mMeanLumiOverride);
   }
+  if (mMeanLumiRefOverride != 0.) {
+    setMeanLumiRef(mMeanLumiRefOverride);
+  }
   if (mInstLumiOverride != 0.) {
     setInstLumi(mInstLumiOverride);
   }
@@ -159,21 +166,24 @@ void CorrectionMapsLoader::init(o2::framework::InitContext& ic)
     LOGP(fatal, "Wrong lumi-scale-type provided!");
   }
 
-  LOGP(info, "Scaling for TPC corr.map scaling={}, override values: lumiMean={} lumiInst={} lumiScaleMode={}, LumiInst scale={}, CTP Lumi source={}",
-       lumiS[scaleType], mMeanLumiOverride, mInstLumiOverride, mLumiScaleMode, mInstLumiFactor, mCTPLumiSource);
+  LOGP(info, "Scaling for TPC corr.map scaling={}, override values: lumiMean={} lumiRefMean={} lumiInst={} lumiScaleMode={}, LumiInst scale={}, CTP Lumi source={}",
+       lumiS[scaleType], mMeanLumiOverride, mMeanLumiRefOverride, mInstLumiOverride, mLumiScaleMode, mInstLumiFactor, mCTPLumiSource);
 }
 
 //________________________________________________________
 void CorrectionMapsLoader::copySettings(const CorrectionMapsLoader& src)
 {
   setInstLumi(src.getInstLumi(), false);
   setMeanLumi(src.getMeanLumi(), false);
+  setMeanLumiRef(src.getMeanLumiRef());
   setLumiScaleType(src.getLumiScaleType());
   setMeanLumiOverride(src.getMeanLumiOverride());
+  setMeanLumiRefOverride(src.getMeanLumiRefOverride());
   setInstLumiOverride(src.getInstLumiOverride());
   setLumiScaleMode(src.getLumiScaleMode());
   mInstLumiFactor = src.mInstLumiFactor;
   mCTPLumiSource = src.mCTPLumiSource;
+  mLumiScaleMode = src.mLumiScaleMode;
 }
 
 #endif // #ifndef GPUCA_GPUCODE_DEVICE

Detectors/TPC/workflow/include/TPCWorkflow/RecoWorkflow.h

@@ -85,7 +85,8 @@ framework::WorkflowSpec getWorkflow(CompletionPolicyData* policyData,
                                     bool selIR = false,
                                     bool filteredInp = false,
                                     int lumiScaleType = 0,
-                                    int deadMapSources = -1);
+                                    int deadMapSources = -1,
+                                    int lumiScaleMode = 0);
 
 void cleanupCallback();
 

Detectors/TPC/workflow/src/RecoWorkflow.cxx

@@ -98,7 +98,7 @@ const std::unordered_map<std::string, OutputType> OutputMap{
 
 framework::WorkflowSpec getWorkflow(CompletionPolicyData* policyData, std::vector<int> const& tpcSectors, unsigned long tpcSectorMask, std::vector<int> const& laneConfiguration,
                                     bool propagateMC, unsigned nLanes, std::string const& cfgInput, std::string const& cfgOutput, bool disableRootInput,
-                                    int caClusterer, int zsOnTheFly, bool askDISTSTF, bool selIR, bool filteredInp, int lumiScaleType, int deadMapSources)
+                                    int caClusterer, int zsOnTheFly, bool askDISTSTF, bool selIR, bool filteredInp, int lumiScaleType, int deadMapSources, int lumiScaleMode)
 {
   InputType inputType;
   try {
@@ -441,6 +441,7 @@ framework::WorkflowSpec getWorkflow(CompletionPolicyData* policyData, std::vecto
     o2::gpu::GPURecoWorkflowSpec::Config cfg;
     cfg.runTPCTracking = true;
     cfg.lumiScaleType = lumiScaleType;
+    cfg.lumiScaleMode = lumiScaleMode;
     cfg.decompressTPC = decompressTPC;
     cfg.decompressTPCFromROOT = decompressTPC && inputType == InputType::CompClusters;
     cfg.caClusterer = caClusterer;

Detectors/TPC/workflow/src/tpc-reco-workflow.cxx

@@ -71,6 +71,7 @@ void customize(std::vector<o2::framework::ConfigParamSpec>& workflowOptions)
     {"lumi-type", o2::framework::VariantType::Int, 0, {"1 = require CTP lumi for TPC correction scaling, 2 = require TPC scalers for TPC correction scaling"}},
     {"select-ir-frames", VariantType::Bool, false, {"Subscribe and filter according to external IR Frames"}},
     {"tpc-deadMap-sources", VariantType::Int, -1, {"Sources to consider for TPC dead channel map creation; -1=all, 0=deactivated"}},
+    {"corrmap-lumi-mode", o2::framework::VariantType::Int, 0, {"scaling mode: (default) 0 = static + scale * full; 1 = full + scale * derivative"}},
   };
   o2::raw::HBFUtilsInitializer::addConfigOption(options);
   std::swap(workflowOptions, options);
@@ -182,7 +183,8 @@ WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)
                                                 cfgc.options().get<bool>("select-ir-frames"),
                                                 cfgc.options().get<bool>("filtered-input"),
                                                 lumiType,
-                                                cfgc.options().get<int>("tpc-deadMap-sources"));
+                                                cfgc.options().get<int>("tpc-deadMap-sources"),
+                                                cfgc.options().get<int>("corrmap-lumi-mode"));
 
   // configure dpl timer to inject correct firstTForbit: start from the 1st orbit of TF containing 1st sampled orbit
   o2::raw::HBFUtilsInitializer hbfIni(cfgc, wf);

GPU/TPCFastTransformation/CorrectionMapsHelper.cxx

@@ -26,6 +26,7 @@ void CorrectionMapsHelper::clear()
   mUpdatedFlags = 0;
   mInstLumi = 0.f;
   mMeanLumi = 0.f;
+  mMeanLumiRef = 0.f;
 }
 
 void CorrectionMapsHelper::setOwner(bool v)

GPU/TPCFastTransformation/CorrectionMapsHelper.h

@@ -80,6 +80,13 @@ class CorrectionMapsHelper
     }
   }
 
+  void setMeanLumiRef(float v)
+  {
+    if (v != mMeanLumi) {
+      mMeanLumiRef = v;
+    }
+  }
+
   void setLumiScaleMode(int v)
   {
     if (v != mLumiScaleMode) {
@@ -93,7 +100,8 @@ class CorrectionMapsHelper
     if (mMeanLumi < 0.f || mInstLumi < 0.f) {
       mLumiScale = -1.f;
     } else if (mLumiScaleMode == 1) {
-      mLumiScale = mMeanLumi ? mInstLumi / mMeanLumi - 1. : 0.f;
+      mLumiScale = mMeanLumiRef ? (mInstLumi - mMeanLumi) / mMeanLumiRef : 0.f;
+      LOGP(info, "mInstLumi: {}  mMeanLumi: {} mMeanLumiRef: {}", mInstLumi, mMeanLumi, mMeanLumiRef);
     } else {
       mLumiScale = mMeanLumi ? mInstLumi / mMeanLumi : 0.f;
     }
@@ -105,6 +113,8 @@ class CorrectionMapsHelper
 
   GPUd() float getInstLumi() const { return mInstLumi; }
   GPUd() float getMeanLumi() const { return mMeanLumi; }
+  GPUd() float getMeanLumiRef() const { return mMeanLumiRef; }
+
   GPUd() float getLumiScale() const { return mLumiScale; }
   GPUd() int getLumiScaleMode() const { return mLumiScaleMode; }
 
@@ -127,7 +137,9 @@ class CorrectionMapsHelper
   int getLumiScaleType() const { return mLumiScaleType; }
 
   void setMeanLumiOverride(float f) { mMeanLumiOverride = f; }
+  void setMeanLumiRefOverride(float f) { mMeanLumiRefOverride = f; }
   float getMeanLumiOverride() const { return mMeanLumiOverride; }
+  float getMeanLumiRefOverride() const { return mMeanLumiRefOverride; }
 
   void setInstLumiOverride(float f) { mInstLumiOverride = f; }
   float getInstLumiOverride() const { return mInstLumiOverride; }
@@ -143,14 +155,16 @@ class CorrectionMapsHelper
   int mUpdatedFlags = 0;
   float mInstLumi = 0.;                                         // instanteneous luminosity (a.u)
   float mMeanLumi = 0.;                                         // mean luminosity of the map (a.u)
+  float mMeanLumiRef = 0.;                                      // mean luminosity of the ref map (a.u)
   float mLumiScale = 0.;                                        // precalculated mInstLumi/mMeanLumi
   int mLumiScaleMode = 0;                                       // scaling-mode of the correciton maps
   float mMeanLumiOverride = -1.f;                               // optional value to override mean lumi
+  float mMeanLumiRefOverride = -1.f;                            // optional value to override ref mean lumi
   float mInstLumiOverride = -1.f;                               // optional value to override inst lumi
   GPUCA_NAMESPACE::gpu::TPCFastTransform* mCorrMap{nullptr};    // current transform
   GPUCA_NAMESPACE::gpu::TPCFastTransform* mCorrMapRef{nullptr}; // reference transform
 #ifndef GPUCA_ALIROOT_LIB
-  ClassDefNV(CorrectionMapsHelper, 3);
+  ClassDefNV(CorrectionMapsHelper, 4);
 #endif
 };
 

GPU/TPCFastTransformation/TPCFastTransform.h

@@ -442,7 +442,7 @@ GPUdi() void TPCFastTransform::TransformInternal(int slice, int row, float& u, f
 {
   if (mApplyCorrection) {
     float dx = 0.f, du = 0.f, dv = 0.f;
-    if (scale >= 0.f) {
+    if ((scale >= 0.f) || (scaleMode == 1)) {
 #ifndef GPUCA_GPUCODE
       if (mCorrectionSlow) {
         float ly, lz;
@@ -463,14 +463,14 @@ GPUdi() void TPCFastTransform::TransformInternal(int slice, int row, float& u, f
 #endif // GPUCA_GPUCODE
       {
         mCorrection.getCorrection(slice, row, u, v, dx, du, dv);
-        if (ref && scale > 0.f) { // scaling was requested
-          if (scaleMode == 0) {
+        if (ref) {
+          if ((scale > 0.f) && (scaleMode == 0)) { // scaling was requested
             float dxRef, duRef, dvRef;
             ref->mCorrection.getCorrection(slice, row, u, v, dxRef, duRef, dvRef);
             dx = (dx - dxRef) * scale + dxRef;
             du = (du - duRef) * scale + duRef;
             dv = (dv - dvRef) * scale + dvRef;
-          } else if (scaleMode == 1) {
+          } else if ((scale != 0.f) && (scaleMode == 1)) {
             float dxRef, duRef, dvRef;
             ref->mCorrection.getCorrection(slice, row, u, v, dxRef, duRef, dvRef);
             dx = dxRef * scale + dx;
@@ -499,8 +499,21 @@ GPUdi() void TPCFastTransform::TransformInternal(int slice, int row, float& u, f
       float YZtoNominalY;
       float YZtoNominalZ;
       InverseTransformYZtoNominalYZ(slice, row, ly, lz, YZtoNominalY, YZtoNominalZ);
+
+      float dxRef, duRef, dvRef;
+      ref->mCorrection.getCorrection(slice, row, u, v, dxRef, duRef, dvRef);
+
+      float dxOrig, duOrig, dvOrig;
+      mCorrection.getCorrection(slice, row, u, v, dxOrig, duOrig, dvOrig);
+
       o2::utils::DebugStreamer::instance()->getStreamer("debug_fasttransform", "UPDATE") << o2::utils::DebugStreamer::instance()->getUniqueTreeName("tree_Transform").data()
                                                                                          // corrections in x, u, v
+                                                                                         << "dxOrig=" << dxOrig
+                                                                                         << "duOrig=" << duOrig
+                                                                                         << "dvOrig=" << dvOrig
+                                                                                         << "dxRef=" << dxRef
+                                                                                         << "duRef=" << duRef
+                                                                                         << "dvRef=" << dvRef
                                                                                          << "dx=" << dx
                                                                                          << "du=" << du
                                                                                          << "dv=" << dv
@@ -525,6 +538,7 @@ GPUdi() void TPCFastTransform::TransformInternal(int slice, int row, float& u, f
                                                                                          << "invYZtoX=" << invYZtoX
                                                                                          << "YZtoNominalY=" << YZtoNominalY
                                                                                          << "YZtoNominalZ=" << YZtoNominalZ
+                                                                                         << "scaleMode=" << scaleMode
                                                                                          << "\n";
     })
 

GPU/Workflow/src/gpu-reco-workflow.cxx

@@ -53,6 +53,7 @@ void customize(std::vector<ConfigParamSpec>& workflowOptions)
     {"input-type", VariantType::String, "digits", {"digitizer, digits, zsraw, zsonthefly, clustersnative, compressed-clusters-root, compressed-clusters-ctf, trd-tracklets"}},
     {"output-type", VariantType::String, "tracks", {"clustersnative, tracks, compressed-clusters-ctf, qa, no-shared-cluster-map, send-clusters-per-sector, trd-tracks, error-qa, tpc-triggers"}},
     {"lumi-type", VariantType::Int, 0, {"1 = require CTP lumi for TPC correction scaling, 2 = require TPC scalers for TPC correction scaling"}},
+    {"corrmap-lumi-mode", VariantType::Int, 0, {"scaling mode: (default) 0 = static + scale * full; 1 = full + scale * derivative"}},
     {"disable-root-input", VariantType::Bool, true, {"disable root-files input reader"}},
     {"disable-mc", VariantType::Bool, false, {"disable sending of MC information"}},
     {"ignore-dist-stf", VariantType::Bool, false, {"do not subscribe to FLP/DISTSUBTIMEFRAME/0 message (no lost TF recovery)"}},
@@ -137,6 +138,7 @@ WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)
 
   auto inputType = cfgc.options().get<std::string>("input-type");
   auto lumiType = cfgc.options().get<int>("lumi-type");
+  auto lumiMode = cfgc.options().get<int>("corrmap-lumi-mode");
   bool doMC = !cfgc.options().get<bool>("disable-mc");
 
   o2::conf::ConfigurableParam::updateFromFile(cfgc.options().get<std::string>("configFile"));
@@ -158,6 +160,7 @@ WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)
   GPURecoWorkflowSpec::Config cfg;
   cfg.runTPCTracking = true;
   cfg.lumiScaleType = lumiType;
+  cfg.lumiScaleMode = lumiMode;
   cfg.decompressTPC = isEnabled(inputTypes, ioType::CompClustCTF);
   cfg.decompressTPCFromROOT = isEnabled(inputTypes, ioType::CompClustROOT);
   cfg.zsDecoder = isEnabled(inputTypes, ioType::ZSRaw);

prodtests/full-system-test/dpl-workflow.sh

@@ -273,8 +273,10 @@ ASK_CTP_LUMI_GPU=
 local restOpt=
 while [[ $# -gt 0 ]]; do
   case "$1" in
-    --lumi-type=*) ASK_CTP_LUMI_GPU=" --lumi-type ${1#*=}"; [[ ${1#*=} == "2" ]] && NEED_TPC_SCALERS_WF=1; shift 1;;
-    --lumi-type) ASK_CTP_LUMI_GPU=" --lumi-type ${2}"; [[ ${2} == "2" ]] && NEED_TPC_SCALERS_WF=1; shift 2;;
+    --lumi-type=*) ASK_CTP_LUMI_GPU="${ASK_CTP_LUMI_GPU} --lumi-type ${1#*=}"; shift 1;;
+    --lumi-type) ASK_CTP_LUMI_GPU="${ASK_CTP_LUMI_GPU} --lumi-type ${2}"; shift 2;;
+    --corrmap-lumi-mode=*) ASK_CTP_LUMI_GPU="${ASK_CTP_LUMI_GPU} --corrmap-lumi-mode ${1#*=}"; shift 1;;
+    --corrmap-lumi-mode) ASK_CTP_LUMI_GPU="${ASK_CTP_LUMI_GPU} --corrmap-lumi-mode ${2}"; shift 2;;
     *) restOpt+=" $1"; shift 1;;
   esac
 done