TPC: adding option to use custom scalers (IDCs) for SC correction map scaling (#12205)

* TPC: adding option to use custom scalers (IDCs) for SC correction map scaling

Adding a generic method to provide scaling of the space-charge
correction maps based on the IDCs (and cluster currents).
In the `o2-tpc-scaler-workflow` for each TF the last n values (dependent
on the specified ion drift time and integration length of the scalers)
are averaged and provided as scaling value.

- A- and C-side values are averaged for now!

Example how to run with scalersfrom local file:
```
  scaler="file:///data/mkleiner/run3/LHC23zo/539399_lowIR/=TPC/Calib/Scaler"

  o2-ctf-reader-workflow --copy-cmd no-copy --onlyDet TPC --ctf-input ${INPUT} --severity warning --condition-remap "${scaler}" |
  o2-tpc-scaler-workflow |
  o2-tpc-reco-workflow --lumi-type 2 --input-type compressed-clusters-ctf --output-type "tracks,disable-writer" --disable-mc -b
```

* Fixing setting of lumi from TPC scalers

* Adding check for run number, adding option to set ion drift time

- Adding mInstLumiFactor as optional scaling factor

* Parse lumi-type from TPC_CORR_SCALING to ASK_CTP_LUMI_GPU

* Avoid copy when loading object from file

Author: matthias-kleiner

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

@@ -88,6 +88,7 @@ enum class CDBType {
   CalCorrDerivMap,     ///< Cluster correction map (derivative map)
                        ///
   CalTimeSeries,       ///< integrated DCAs for longer time interval
+  CalScaler,           ///< Scaler from IDCs or combined estimator
 };
 
 /// Upload intervention type
@@ -149,6 +150,7 @@ const std::unordered_map<CDBType, const std::string> CDBTypeMap{
   {CDBType::CalCorrDerivMap, "TPC/Calib/CorrectionMapDerivativeV2"},
   // time series
   {CDBType::CalTimeSeries, "TPC/Calib/TimeSeries"},
+  {CDBType::CalScaler, "TPC/Calib/Scaler"},
 };
 
 /// Poor enum reflection ...

Detectors/TPC/calibration/CMakeLists.txt

@@ -52,6 +52,7 @@ o2_add_library(TPCCalibration
                        src/SACCCDBHelper.cxx
                        src/TPCFastSpaceChargeCorrectionHelper.cxx
                        src/CalculatedEdx.cxx
+                       src/TPCScaler.cxx
                PUBLIC_LINK_LIBRARIES O2::DataFormatsTPC O2::TPCBase
                                      O2::TPCReconstruction ROOT::Minuit
                                      Microsoft.GSL::GSL
@@ -103,7 +104,8 @@ o2_target_root_dictionary(TPCCalibration
                                   include/TPCCalibration/VDriftHelper.h
                                   include/TPCCalibration/SACCCDBHelper.h
                                   include/TPCCalibration/TPCFastSpaceChargeCorrectionHelper.h
-                                  include/TPCCalibration/CalculatedEdx.h)
+                                  include/TPCCalibration/CalculatedEdx.h
+                                  include/TPCCalibration/TPCScaler.h)
 
 o2_add_test_root_macro(macro/comparePedestalsAndNoise.C
                        PUBLIC_LINK_LIBRARIES O2::TPCBase

Detectors/TPC/calibration/include/TPCCalibration/CorrectionMapsLoader.h

@@ -50,7 +50,7 @@ class CorrectionMapsLoader : public o2::gpu::CorrectionMapsHelper
   void init(o2::framework::InitContext& ic);
   void copySettings(const CorrectionMapsLoader& src);
 
-  static void requestCCDBInputs(std::vector<o2::framework::InputSpec>& inputs, std::vector<o2::framework::ConfigParamSpec>& options, bool requestCTPLumi = false, int lumiScaleMode = 0);
+  static void requestCCDBInputs(std::vector<o2::framework::InputSpec>& inputs, std::vector<o2::framework::ConfigParamSpec>& options, int lumiScaleType = 0, int lumiScaleMode = 0);
   static void addOptions(std::vector<o2::framework::ConfigParamSpec>& options);
 
  protected:

Detectors/TPC/calibration/include/TPCCalibration/TPCScaler.h

@@ -0,0 +1,113 @@
+// Copyright 2019-2020 CERN and copyright holders of ALICE O2.
+// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
+// All rights not expressly granted are reserved.
+//
+// This software is distributed under the terms of the GNU General Public
+// License v3 (GPL Version 3), copied verbatim in the file "COPYING".
+//
+// In applying this license CERN does not waive the privileges and immunities
+// granted to it by virtue of its status as an Intergovernmental Organization
+// or submit itself to any jurisdiction.
+
+/// \file TPCScaler.h
+/// \author Matthias Kleiner <mkleiner@ikf.uni-frankfurt.de>
+
+#ifndef ALICEO2_TPC_TPCSCALER
+#define ALICEO2_TPC_TPCSCALER
+
+#include "DataFormatsTPC/Defs.h"
+#include <vector>
+
+class TTree;
+
+namespace o2::tpc
+{
+
+/*
+Class for storing the scalers which are used to calculate an estimate for the mean space-charge density for the last ion drift time
+*/
+
+class TPCScaler
+{
+ public:
+  /// default constructor
+  TPCScaler() = default;
+
+  /// default move assignment
+  TPCScaler& operator=(TPCScaler&& other) = default;
+
+  /// \return returns number of stored TPC scaler values
+  int getNValues(o2::tpc::Side side) const { return (side == o2::tpc::Side::A ? mScalerA.size() : mScalerC.size()); }
+
+  /// set the parameters for the coefficients of the polynomial
+  /// \param params parameter for the coefficients
+  void setScaler(const std::vector<float>& values, const o2::tpc::Side side) { (side == o2::tpc::Side::A ? (mScalerA = values) : (mScalerC = values)); }
+
+  /// \return returns ion drift time in ms
+  void setIonDriftTimeMS(float ionDriftTimeMS) { mIonDriftTimeMS = ionDriftTimeMS; }
+
+  /// \return returns run number for which this object is valid
+  void setRun(int run) { mRun = run; }
+
+  /// \param firstTFOrbit first TF orbit of first data
+  void setFirstTFOrbit(unsigned int firstTFOrbit) { mFirstTFOrbit = firstTFOrbit; }
+
+  /// \param timeStampMS first time in ms
+  void setStartTimeStampMS(double timeStampMS) { mTimeStampMS = timeStampMS; }
+
+  /// \param integrationTimeMS integration time for each scaler value
+  void setIntegrationTimeMS(float integrationTimeMS) { mIntegrationTimeMS = integrationTimeMS; }
+
+  /// dump this object to a file
+  /// \param file output file
+  void dumpToFile(const char* file, const char* name);
+
+  /// load parameters from input file (which were written using the writeToFile method)
+  /// \param inpf input file
+  void loadFromFile(const char* inpf, const char* name);
+
+  /// set this object from input tree
+  void setFromTree(TTree& tpcScalerTree);
+
+  /// \return returns stored scalers for given side and data index
+  float getScalers(unsigned int idx, o2::tpc::Side side) const { return (side == o2::tpc::Side::A) ? mScalerA[idx] : mScalerC[idx]; }
+
+  /// \return returns stored scalers for given side
+  const auto& getScalers(o2::tpc::Side side) const { return (side == o2::tpc::Side::A) ? mScalerA : mScalerC; }
+
+  /// \return returns ion drift time in ms
+  float getIonDriftTimeMS() const { return mIonDriftTimeMS; }
+
+  /// \return returns run number for which this object is valid
+  int getRun() const { return mRun; }
+
+  /// \return returns first TF orbit of first data
+  unsigned int getFirstTFOrbit() const { return mFirstTFOrbit; }
+
+  /// \return return first time in ms
+  double getStartTimeStampMS() const { return mTimeStampMS; }
+
+  /// \return returns integration time for each scaler value
+  float getIntegrationTimeMS() const { return mIntegrationTimeMS; }
+
+  /// \return returns numbers of scalers for one ion drift time
+  int getNValuesIonDriftTime() const { return mIonDriftTimeMS / mIntegrationTimeMS + 0.5; }
+
+  /// \return returns mean scaler value for last ion drift time
+  /// \param timestamp timestamp for which the last values are used to calculate the mean
+  float getMeanScaler(double timestamp, o2::tpc::Side side) const;
+
+ private:
+  float mIonDriftTimeMS{200};    ///< ion drift time in ms
+  int mRun{};                    ///< run for which this object is valid
+  unsigned int mFirstTFOrbit{};  ///< first TF orbit of the stored scalers
+  double mTimeStampMS{};         ///< time stamp of the first stored values
+  float mIntegrationTimeMS{1.};  ///< integration time for each stored value in ms
+  std::vector<float> mScalerA{}; ///< TPC scaler for A-side
+  std::vector<float> mScalerC{}; ///< TPC scaler for C-side
+
+  ClassDefNV(TPCScaler, 1);
+};
+
+} // namespace o2::tpc
+#endif

Detectors/TPC/calibration/src/CorrectionMapsLoader.cxx

@@ -46,7 +46,7 @@ void CorrectionMapsLoader::extractCCDBInputs(ProcessingContext& pc)
   int dumRep = 0;
   o2::ctp::LumiInfo lumiObj;
   static o2::ctp::LumiInfo lumiPrev;
-  if (getUseCTPLumi() && mInstLumiOverride <= 0.) {
+  if (getLumiScaleType() == 1 && mInstLumiOverride <= 0.) {
     if (pc.inputs().get<gsl::span<char>>("CTPLumi").size() == sizeof(o2::ctp::LumiInfo)) {
       lumiPrev = lumiObj = pc.inputs().get<o2::ctp::LumiInfo>("CTPLumi");
     } else {
@@ -56,13 +56,16 @@ void CorrectionMapsLoader::extractCCDBInputs(ProcessingContext& pc)
       lumiObj = lumiPrev;
     }
     setInstLumi(mInstLumiFactor * (mCTPLumiSource == 0 ? lumiObj.getLumi() : lumiObj.getLumiAlt()));
+  } else if (getLumiScaleType() == 2 && mInstLumiOverride <= 0.) {
+    float tpcScaler = pc.inputs().get<float>("tpcscaler");
+    setInstLumi(mInstLumiFactor * tpcScaler);
   }
 }
 
 //________________________________________________________
-void CorrectionMapsLoader::requestCCDBInputs(std::vector<InputSpec>& inputs, std::vector<o2::framework::ConfigParamSpec>& options, bool requestCTPLumi, int lumiScaleMode)
+void CorrectionMapsLoader::requestCCDBInputs(std::vector<InputSpec>& inputs, std::vector<o2::framework::ConfigParamSpec>& options, int lumiScaleType, int lumiScaleMode)
 {
-  addInput(inputs, {"tpcCorrMap", "TPC", "CorrMap", 0, Lifetime::Condition, ccdbParamSpec(CDBTypeMap.at(CDBType::CalCorrMap), {}, 1)});          // time-dependent
+  addInput(inputs, {"tpcCorrMap", "TPC", "CorrMap", 0, Lifetime::Condition, ccdbParamSpec(CDBTypeMap.at(CDBType::CalCorrMap), {}, 1)}); // time-dependent
   if (lumiScaleMode == 0) {
     addInput(inputs, {"tpcCorrMapRef", "TPC", "CorrMapRef", 0, Lifetime::Condition, ccdbParamSpec(CDBTypeMap.at(CDBType::CalCorrMapRef), {}, 0)}); // load once
   } else if (lumiScaleMode == 1) {
@@ -71,8 +74,10 @@ void CorrectionMapsLoader::requestCCDBInputs(std::vector<InputSpec>& inputs, std
     LOG(fatal) << "Correction mode unknown! Choose either 0 (default) or 1 (derivative map) for flag corrmap-lumi-mode.";
   }
 
-  if (requestCTPLumi) {
+  if (lumiScaleType == 1) {
     addInput(inputs, {"CTPLumi", "CTP", "LUMI", 0, Lifetime::Timeframe});
+  } else if (lumiScaleType == 2) {
+    addInput(inputs, {"tpcscaler", o2::header::gDataOriginTPC, "TPCSCALER", 0, Lifetime::Timeframe});
   }
   addOptions(options);
 }
@@ -131,8 +136,10 @@ void CorrectionMapsLoader::init(o2::framework::InitContext& ic)
   const auto& inputRouts = ic.services().get<const o2::framework::DeviceSpec>().inputs;
   for (const auto& route : inputRouts) {
     if (route.matcher == InputSpec{"CTPLumi", "CTP", "LUMI", 0, Lifetime::Timeframe}) {
-      setUseCTPLumi(true);
+      setLumiScaleType(1);
       break;
+    } else if (route.matcher == InputSpec{"tpcscaler", o2::header::gDataOriginTPC, "TPCSCALER", 0, Lifetime::Timeframe}) {
+      setLumiScaleType(2);
     }
   }
   mMeanLumiOverride = ic.options().get<float>("corrmap-lumi-mean");
@@ -146,16 +153,22 @@ void CorrectionMapsLoader::init(o2::framework::InitContext& ic)
   if (mInstLumiOverride != 0.) {
     setInstLumi(mInstLumiOverride);
   }
-  LOGP(info, "CTP Lumi request for TPC corr.map scaling={}, override values: lumiMean={} lumiInst={} lumiScaleMode={}, LumiInst scale={}, CTP Lumi source={}",
-       getUseCTPLumi() ? "ON" : "OFF", mMeanLumiOverride, mInstLumiOverride, mLumiScaleMode, mInstLumiFactor, mCTPLumiSource);
+  const std::array<std::string, 3> lumiS{"OFF", "CTP", "TPC scaler"};
+  int scaleType = getLumiScaleType();
+  if (scaleType >= lumiS.size()) {
+    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);
 }
 
 //________________________________________________________
 void CorrectionMapsLoader::copySettings(const CorrectionMapsLoader& src)
 {
   setInstLumi(src.getInstLumi(), false);
   setMeanLumi(src.getMeanLumi(), false);
-  setUseCTPLumi(src.getUseCTPLumi());
+  setLumiScaleType(src.getLumiScaleType());
   setMeanLumiOverride(src.getMeanLumiOverride());
   setInstLumiOverride(src.getInstLumiOverride());
   setLumiScaleMode(src.getLumiScaleMode());

Detectors/TPC/calibration/src/TPCCalibrationLinkDef.h

@@ -110,4 +110,5 @@
 #pragma link C++ class o2::tpc::TPCFastSpaceChargeCorrectionHelper + ;
 
 #pragma link C++ class o2::tpc::CalculatedEdx + ;
+#pragma link C++ class o2::tpc::TPCScaler + ;
 #endif

Detectors/TPC/calibration/src/TPCScaler.cxx

@@ -0,0 +1,74 @@
+// Copyright 2019-2020 CERN and copyright holders of ALICE O2.
+// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
+// All rights not expressly granted are reserved.
+//
+// This software is distributed under the terms of the GNU General Public
+// License v3 (GPL Version 3), copied verbatim in the file "COPYING".
+//
+// In applying this license CERN does not waive the privileges and immunities
+// granted to it by virtue of its status as an Intergovernmental Organization
+// or submit itself to any jurisdiction.
+
+/// \file  TPCScaler.cxx
+/// \brief Definition of TPCScaler class
+///
+/// \author Matthias Kleiner <mkleiner@ikf.uni-frankfurt.de>
+
+#include "TPCCalibration/TPCScaler.h"
+#include <TFile.h>
+#include <TTree.h>
+#include "Framework/Logger.h"
+
+using namespace o2::tpc;
+
+void TPCScaler::dumpToFile(const char* file, const char* name)
+{
+  TFile out(file, "RECREATE");
+  TTree tree("TPCScaler", "TPCScaler");
+  tree.Branch("TPCScaler", this);
+  tree.Fill();
+  out.WriteObject(&tree, name);
+}
+
+void TPCScaler::loadFromFile(const char* inpf, const char* name)
+{
+  TFile out(inpf, "READ");
+  TTree* tree = (TTree*)out.Get(name);
+  setFromTree(*tree);
+}
+
+void TPCScaler::setFromTree(TTree& tpcScalerTree)
+{
+  TPCScaler* scalerTmp = this;
+  tpcScalerTree.SetBranchAddress("TPCScaler", &scalerTmp);
+  const int entries = tpcScalerTree.GetEntries();
+  if (entries > 0) {
+    tpcScalerTree.GetEntry(0);
+  } else {
+    LOGP(error, "TPCScaler not found in input file");
+  }
+  tpcScalerTree.SetBranchAddress("TPCScaler", nullptr);
+}
+
+float TPCScaler::getMeanScaler(double timestamp, o2::tpc::Side side) const
+{
+  // index to data buffer
+  const int idxData = (timestamp - mTimeStampMS) / mIntegrationTimeMS + 0.5;
+  const int nVals = getNValuesIonDriftTime();
+  const int nValues = getNValues(side);
+  if ((nVals == 0) || (nVals > nValues)) {
+    return -1;
+    LOGP(error, "Empty data provided {}", nValues);
+  }
+
+  // clamp indices to min and max
+  const int lastIdx = std::clamp(idxData, nVals, nValues);
+  const int firstIdx = (lastIdx == nValues) ? (nValues - nVals) : std::clamp(idxData - nVals, 0, nValues);
+
+  // sump up values from last ion drift time
+  float sum = 0;
+  for (int i = firstIdx; i < lastIdx; ++i) {
+    sum += getScalers(i, side);
+  }
+  return (sum / nVals);
+}

Detectors/TPC/workflow/CMakeLists.txt

@@ -45,6 +45,7 @@ o2_add_library(TPCWorkflow
                        src/TPCTimeSeriesWriterSpec.cxx
                        src/TPCTimeSeriesReaderSpec.cxx
                        src/TPCMergeTimeSeriesSpec.cxx
+                       src/TPCScalerSpec.cxx
                TARGETVARNAME targetName
                PUBLIC_LINK_LIBRARIES O2::Framework O2::DataFormatsTPC
                                      O2::DPLUtils O2::TPCReconstruction
@@ -246,6 +247,11 @@ o2_add_executable(merge-time-series-workflow
                   COMPONENT_NAME tpc
                   PUBLIC_LINK_LIBRARIES O2::TPCWorkflow)
 
+o2_add_executable(scaler-workflow
+                  COMPONENT_NAME tpc
+                  SOURCES src/tpc-scaler.cxx
+                  PUBLIC_LINK_LIBRARIES O2::TPCWorkflow)
+
 o2_add_test(workflow
             COMPONENT_NAME tpc
             LABELS tpc workflow

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

@@ -84,7 +84,7 @@ framework::WorkflowSpec getWorkflow(CompletionPolicyData* policyData,
                                     bool askDISTSTF = true,
                                     bool selIR = false,
                                     bool filteredInp = false,
-                                    bool requireCTPLumi = false,
+                                    int lumiScaleType = 0,
                                     int deadMapSources = -1);
 
 void cleanupCallback();

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

@@ -240,7 +240,7 @@ class TPCCalibPadGainTracksDevice : public o2::framework::Task
   }
 };
 
-DataProcessorSpec getTPCCalibPadGainTracksSpec(const uint32_t publishAfterTFs, const bool debug, const bool useLastExtractedMapAsReference, const std::string polynomialsFile, bool disablePolynomialsCCDB, bool requestCTPLumi)
+DataProcessorSpec getTPCCalibPadGainTracksSpec(const uint32_t publishAfterTFs, const bool debug, const bool useLastExtractedMapAsReference, const std::string polynomialsFile, bool disablePolynomialsCCDB, int lumiType)
 {
   std::vector<InputSpec> inputs;
   inputs.emplace_back("trackTPC", gDataOriginTPC, "TRACKS", 0, Lifetime::Timeframe);
@@ -282,7 +282,7 @@ DataProcessorSpec getTPCCalibPadGainTracksSpec(const uint32_t publishAfterTFs, c
     {"useEveryNthTF", VariantType::Int, 10, {"Using only a fraction of the data: 1: Use every TF, 10: Use only every tenth TF."}},
     {"maxTracksPerTF", VariantType::Int, 10000, {"Maximum number of processed tracks per TF (-1 for processing all tracks)"}},
   };
-  o2::tpc::CorrectionMapsLoader::requestCCDBInputs(inputs, opts, requestCTPLumi);
+  o2::tpc::CorrectionMapsLoader::requestCCDBInputs(inputs, opts, lumiType);
 
   auto ccdbRequest = std::make_shared<o2::base::GRPGeomRequest>(false,                          // orbitResetTime
                                                                 false,                          // GRPECS=true