TPC: Add sector edge fluctuation correction methods

Add option to correct for sector edge fluctuations.
The time intervals and sectors with optional scaling factors are parsed from CSV table.
Two new CCDB types are added:
1. CalSecEdgeCorrection: the correction map with per sector delta fluctuations
2. CalSecEdgeInfo: the time intervals with optional per sector scaling factors for which the CalSecEdgeCorrection will be applied

The applied scaling factors are added to the TPC timeseries for visualisation.

Author: matthias-kleiner

Detectors/Calibration/include/DetectorsCalibration/IntegratedClusterCalibrator.h

@@ -331,17 +331,18 @@ struct TimeSeriesdEdx {
 };
 
 struct TimeSeriesITSTPC {
-  float mVDrift = 0;                          ///< drift velocity in cm/us
-  float mPressure = 0;                        ///< pressure
-  float mTemperature = 0;                     ///< temperature
-  TimeSeries mTSTPC;                          ///< TPC standalone DCAs
-  TimeSeries mTSITSTPC;                       ///< ITS-TPC standalone DCAs
-  ITSTPC_Matching mITSTPCAll;                 ///< ITS-TPC matching efficiency for ITS standalone + afterburner
-  ITSTPC_Matching mITSTPCStandalone;          ///< ITS-TPC matching efficiency for ITS standalone
-  ITSTPC_Matching mITSTPCAfterburner;         ///< ITS-TPC matchin efficiency  fir ITS afterburner
-  TimeSeriesdEdx mdEdxQTot;                   ///< time series for dE/dx qTot monitoring
-  TimeSeriesdEdx mdEdxQMax;                   ///< time series for dE/dx qMax monitoring
-  std::vector<unsigned int> mOccupancyMapTPC; ///< cluster occupancy map
+  float mVDrift = 0;                                    ///< drift velocity in cm/us
+  float mPressure = 0;                                  ///< pressure
+  float mTemperature = 0;                               ///< temperature
+  TimeSeries mTSTPC;                                    ///< TPC standalone DCAs
+  TimeSeries mTSITSTPC;                                 ///< ITS-TPC standalone DCAs
+  ITSTPC_Matching mITSTPCAll;                           ///< ITS-TPC matching efficiency for ITS standalone + afterburner
+  ITSTPC_Matching mITSTPCStandalone;                    ///< ITS-TPC matching efficiency for ITS standalone
+  ITSTPC_Matching mITSTPCAfterburner;                   ///< ITS-TPC matchin efficiency  fir ITS afterburner
+  TimeSeriesdEdx mdEdxQTot;                             ///< time series for dE/dx qTot monitoring
+  TimeSeriesdEdx mdEdxQMax;                             ///< time series for dE/dx qMax monitoring
+  std::vector<unsigned int> mOccupancyMapTPC;           ///< cluster occupancy map
+  std::vector<std::pair<int, float>> mSecEdgeFlucCorr;  ///< applied sector edge fluctuation correction
 
   std::vector<float> nPrimVertices;                  ///< number of primary vertices
   std::vector<float> nPrimVertices_ITS;              ///< number of primary vertices selected with ITS cut 0.2<nContributorsITS/nContributors<0.8
@@ -502,7 +503,7 @@ struct TimeSeriesITSTPC {
     nVertexContributors_Quantiles.resize(nTotalQ);
   }
 
-  ClassDefNV(TimeSeriesITSTPC, 6);
+  ClassDefNV(TimeSeriesITSTPC, 7);
 };
 
 } // end namespace tpc

Detectors/TPC/baserecsim/include/TPCBaseRecSim/CDBTypes.h

@@ -83,6 +83,8 @@ enum class CDBType {
   CalScaler,            ///< Scaler from IDCs or combined estimator
   CalScalerWeights,     ///< Weights for scalers
   CalMShape,            ///< calibration object for M-shape distortions
+  CalSecEdgeCorrection, ///< calibration object for sector edge distortions
+  CalSecEdgeInfo,       ///< time slots and scaling for sector edge distortions
                         ///
   CorrMapParam,         ///< parameters for CorrectionMapsLoader configuration
                         ///
@@ -154,6 +156,8 @@ const std::unordered_map<CDBType, const std::string> CDBTypeMap{
   {CDBType::CalScaler, "TPC/Calib/Scaler"},
   {CDBType::CalScalerWeights, "TPC/Calib/ScalerWeights"},
   {CDBType::CalMShape, "TPC/Calib/MShapePotential"},
+  {CDBType::CalSecEdgeCorrection, "TPC/Calib/CorrectionMapSecEdgeFluc"},
+  {CDBType::CalSecEdgeInfo, "TPC/Calib/SecEdgeFlucInfo"},
   // correction maps loader params
   {CDBType::CorrMapParam, "TPC/Calib/CorrMapParam"},
   // distortion maps

Detectors/TPC/calibration/CMakeLists.txt

@@ -61,6 +61,7 @@ o2_add_library(TPCCalibration
                        src/CMVContainer.cxx
                        src/CorrectionMapsLoader.cxx
                        src/CMVHelper.cxx
+                       src/SectorEdgeFluctuations.cxx
                PUBLIC_LINK_LIBRARIES O2::DataFormatsTPC O2::TPCBaseRecSim
                                      O2::TPCReconstruction ROOT::Minuit
                                      Microsoft.GSL::GSL
@@ -121,7 +122,8 @@ o2_target_root_dictionary(TPCCalibration
                                   include/TPCCalibration/PressureTemperatureHelper.h
                                   include/TPCCalibration/CMVContainer.h
                                   include/TPCCalibration/CorrectionMapsLoader.h
-                                  include/TPCCalibration/CMVHelper.h)
+                                  include/TPCCalibration/CMVHelper.h
+                                  include/TPCCalibration/SectorEdgeFluctuations.h)
 
 o2_add_test_root_macro(macro/comparePedestalsAndNoise.C
                        PUBLIC_LINK_LIBRARIES O2::TPCBaseRecSim

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

@@ -19,6 +19,7 @@
 #include <vector>
 #include "CorrectionMapsHelper.h"
 #include "CorrectionMapsTypes.h"
+#include "TPCCalibration/SectorEdgeFluctuations.h"
 
 namespace o2
 {
@@ -47,6 +48,9 @@ class CorrectionMapsLoader : public o2::gpu::CorrectionMapsHelper
   void checkMeanScaleConsistency(float meanLumi, float threshold) const;
 
   static void requestCCDBInputs(std::vector<o2::framework::InputSpec>& inputs, const o2::tpc::CorrectionMapsGloOpts& gloOpts);
+  void enableSecEdgeFlucCorrection(const bool enable = true) { mApplySecEdgeFlucCorr = enable; }
+  bool applySecEdgeFlucCorrection() const { return mApplySecEdgeFlucCorr; }
+  const auto& getSectorEdgeFlucInfo() const { return mSecEdgeFlucInfo; }
 
  protected:
   static void addOption(std::vector<o2::framework::ConfigParamSpec>& options, o2::framework::ConfigParamSpec&& osp);
@@ -55,6 +59,8 @@ class CorrectionMapsLoader : public o2::gpu::CorrectionMapsHelper
   float mInstLumiCTPFactor = 1.0;      // multiplicative factor for inst. lumi
   int mLumiCTPSource = 0;              // 0: main, 1: alternative CTP lumi source
   bool mIDC2CTPFallbackActive = false; // flag indicating that fallback from IDC to CTP scaling is active
+  o2::tpc::SectorEdgeFluctuations mSecEdgeFlucInfo; // definition of sector edge fluctuation distortion map scaling
+  bool mApplySecEdgeFlucCorr = true;   // flag indicating if sector edge fluctuation correction is enabled
 };
 
 } // namespace tpc

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

@@ -0,0 +1,108 @@
+// 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 SectorEdgeFluctuations.h
+/// \brief Class to parse and query time-dependent TPC sector edge fluctuation intervals
+/// \author Matthias Kleiner <matthias.kleiner@cern.ch>
+
+#ifndef ALICEO2_TPC_SECTOREDGEFLUCTUATIONS_H
+#define ALICEO2_TPC_SECTOREDGEFLUCTUATIONS_H
+
+#include <map>
+#include <vector>
+#include <string>
+#include <utility>
+#include "Rtypes.h"
+
+class TTree;
+
+namespace o2::tpc
+{
+
+/// One time interval during which a set of TPC sectors has edge fluctuations.
+/// Each sector carries an optional scaling factor (default 1.0).
+struct SectorEdgeInterval {
+  Long64_t startTimeMS{0};                    ///< interval start, Unix time in ms
+  Long64_t endTimeMS{0};                      ///< interval end,   Unix time in ms
+  std::vector<std::pair<int, float>> sectors; ///< {o2SectorId, scalingFactor}
+  ClassDefNV(SectorEdgeInterval, 1);
+};
+
+/// Parses and queries TPC sector edge fluctuation intervals from a CSV text file.
+///
+/// Expected line format (comma-separated, whitespace around tokens is ignored):
+///   runNumber, startMS, endMS, durationMS, label[, SectorID[=scale], ...]
+///
+/// The sector list is optional. If omitted (or all tokens fail to parse), the
+/// interval is applied to all 36 sectors (A0-A17 and C0-C17) with scale 1.0.
+///
+/// Examples:
+///   560352,1732244770094,1732244771094,1000,edge distortions,A3
+///   560352,1732244771344,1732244776344,5000,edge distortions,A3=1.2,C0,C1=0.6,C2,C3
+///
+/// If the same sector appears in multiple overlapping intervals at a queried
+/// timestamp, the scale from the interval with the latest end-time is used.
+class SectorEdgeFluctuations
+{
+ public:
+  /// Load intervals from file. Clears any previously loaded data.
+  /// Throws std::runtime_error if the file cannot be opened.
+  bool loadFromCSVFile(const std::string& filename);
+
+  /// dump this object to a file
+  /// \param file output file
+  /// \param name name of the output object
+  void dumpToFile(const char* file, const char* name = "ccdb_object", const char* brName = "SectorEdgeFluctuation");
+
+  /// load from input file (which were written using the dumpToFile method)
+  /// \param inpf input file
+  /// \param name name of the object in the file
+  void loadFromFile(const char* inpf, const char* name = "ccdb_object", const int iEntry = 0, const char* brName = "SectorEdgeFluctuation");
+
+  /// set this object from input tree
+  void setFromTree(TTree& tree, const int iEntry = 0, const char* brName = "SectorEdgeFluctuation");
+
+  /// Returns all {o2SectorId, scalingFactor} pairs active for the given run
+  /// at the given Unix timestamp (milliseconds). Returns empty if none are active
+  /// or if the run is not known.
+  std::vector<std::pair<int, float>> getSectorsAtTime(int run, Long64_t timestampMS) const;
+
+  /// Convert a sector string such as "A3" or "C14" to the O2 integer sector
+  /// index (0-35). Returns -1 on parse error.
+  static int parseSectorId(const std::string& sectorStr);
+
+  /// Total number of intervals across all runs.
+  size_t size() const
+  {
+    size_t n = 0;
+    for (const auto& [run, v] : mIntervals) {
+      n += v.size();
+    }
+    return n;
+  }
+
+  /// number of total runs stored
+  size_t getNRuns() const { return mIntervals.size(); }
+  bool empty() const { return mIntervals.empty(); }
+
+  /// get stored data
+  const auto& getIntervals() const { return mIntervals; }
+
+ private:
+  /// Per-run intervals, each sorted by startTimeMS.
+  std::map<int, std::vector<SectorEdgeInterval>> mIntervals;
+
+  ClassDefNV(SectorEdgeFluctuations, 1);
+};
+
+} // namespace o2::tpc
+
+#endif // ALICEO2_TPC_SECTOREDGEFLUCTUATIONS_H

Detectors/TPC/calibration/src/CorrectionMapsLoader.cxx

@@ -35,6 +35,12 @@ void CorrectionMapsLoader::extractCCDBInputs(ProcessingContext& pc, float tpcSca
   if (lumiMode != LumiScaleMode::NoCorrection) {
     pc.inputs().get<o2::gpu::TPCFastTransform*>("tpcCorrMapRef");
   }
+
+  if(mApplySecEdgeFlucCorr) {
+    pc.inputs().get<o2::gpu::TPCFastTransform*>("tpcCorrMapSecFluc");
+    pc.inputs().get<TTree*>("tpSecFlucInfo");
+  }
+
   const int maxDumRep = 5;
   int dumRep = 0;
   o2::ctp::LumiInfo lumiObj;
@@ -93,11 +99,11 @@ void CorrectionMapsLoader::requestCCDBInputs(std::vector<InputSpec>& inputs, con
 {
   LOGP(info, "Requesting CCDB inputs for TPC correction maps with lumiType={} and lumiMode={}", static_cast<int>(gloOpts.lumiType), static_cast<int>(gloOpts.lumiMode));
   if (gloOpts.lumiMode == LumiScaleMode::Linear) {
-    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
+    addInput(inputs, {"tpcCorrMap", o2::header::gDataOriginTPC, "CorrMap", 0, Lifetime::Condition, ccdbParamSpec(CDBTypeMap.at(CDBType::CalCorrMap), {}, 1)});          // time-dependent
+    addInput(inputs, {"tpcCorrMapRef", o2::header::gDataOriginTPC, "CorrMapRef", 0, Lifetime::Condition, ccdbParamSpec(CDBTypeMap.at(CDBType::CalCorrMapRef), {}, 0)}); // load once
   } else if (gloOpts.lumiMode == LumiScaleMode::DerivativeMap) {
-    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
+    addInput(inputs, {"tpcCorrMap", o2::header::gDataOriginTPC, "CorrMap", 0, Lifetime::Condition, ccdbParamSpec(CDBTypeMap.at(CDBType::CalCorrMap), {}, 1)});            // time-dependent
+    addInput(inputs, {"tpcCorrMapRef", o2::header::gDataOriginTPC, "CorrMapRef", 0, Lifetime::Condition, ccdbParamSpec(CDBTypeMap.at(CDBType::CalCorrDerivMap), {}, 1)}); // time-dependent
   } else if (gloOpts.lumiMode == LumiScaleMode::DerivativeMapMC) {
     // for MC corrections
     addInput(inputs, {"tpcCorrMap", "TPC", "CorrMap", 0, Lifetime::Condition, ccdbParamSpec(CDBTypeMap.at(CDBType::CalCorrMapMC), {}, 1)});            // time-dependent
@@ -110,11 +116,17 @@ void CorrectionMapsLoader::requestCCDBInputs(std::vector<InputSpec>& inputs, con
     LOG(fatal) << "Correction mode unknown! Choose either 0 (default) or 1 (derivative map) for flag corrmap-lumi-mode.";
   }
 
+  // load sector edge fluctuation correction only for data
+  if (gloOpts.enableSecEdgeFlucCorrection) {
+    addInput(inputs, {"tpcCorrMapSecFluc", o2::header::gDataOriginTPC, "CorrMapSecFluc", 0, Lifetime::Condition, ccdbParamSpec(CDBTypeMap.at(CDBType::CalSecEdgeCorrection), {}, 1)}); // time-dependent
+    addInput(inputs, {"tpSecFlucInfo", o2::header::gDataOriginTPC, "InfoMapSecFluc", 0, Lifetime::Condition, ccdbParamSpec(CDBTypeMap.at(CDBType::CalSecEdgeInfo), {}, 1)}); // time-dependent
+  }
+
   if (gloOpts.requestCTPLumi) {
     addInput(inputs, {"CTPLumi", "CTP", "LUMI", 0, Lifetime::Timeframe});
   }
 
-  addInput(inputs, {"tpcCorrPar", "TPC", "CorrMapParam", 0, Lifetime::Condition, ccdbParamSpec(CDBTypeMap.at(CDBType::CorrMapParam), {}, 0)}); // load once
+  addInput(inputs, {"tpcCorrPar", o2::header::gDataOriginTPC, "CorrMapParam", 0, Lifetime::Condition, ccdbParamSpec(CDBTypeMap.at(CDBType::CorrMapParam), {}, 0)}); // load once
 }
 
 //________________________________________________________
@@ -136,7 +148,7 @@ void CorrectionMapsLoader::addOption(std::vector<ConfigParamSpec>& options, Conf
 //________________________________________________________
 bool CorrectionMapsLoader::accountCCDBInputs(const ConcreteDataMatcher& matcher, void* obj)
 {
-  if (matcher == ConcreteDataMatcher("TPC", "CorrMap", 0)) {
+  if (matcher == ConcreteDataMatcher(o2::header::gDataOriginTPC, "CorrMap", 0)) {
     setCorrMap((o2::gpu::TPCFastTransform*)obj);
     mCorrMap->rectifyAfterReadingFromFile();
     mCorrMap->setCTP2IDCFallBackThreshold(o2::tpc::CorrMapParam::Instance().CTP2IDCFallBackThreshold);
@@ -165,7 +177,7 @@ bool CorrectionMapsLoader::accountCCDBInputs(const ConcreteDataMatcher& matcher,
     setUpdatedMap();
     return true;
   }
-  if (matcher == ConcreteDataMatcher("TPC", "CorrMapRef", 0)) {
+  if (matcher == ConcreteDataMatcher(o2::header::gDataOriginTPC, "CorrMapRef", 0)) {
     setCorrMapRef((o2::gpu::TPCFastTransform*)obj);
     mCorrMapRef->rectifyAfterReadingFromFile();
     mCorrMapRef->setCTP2IDCFallBackThreshold(o2::tpc::CorrMapParam::Instance().CTP2IDCFallBackThreshold);
@@ -194,7 +206,7 @@ bool CorrectionMapsLoader::accountCCDBInputs(const ConcreteDataMatcher& matcher,
     setUpdatedMapRef();
     return true;
   }
-  if (matcher == ConcreteDataMatcher("TPC", "CorrMapParam", 0)) {
+  if (matcher == ConcreteDataMatcher(o2::header::gDataOriginTPC, "CorrMapParam", 0)) {
     const auto& par = o2::tpc::CorrMapParam::Instance();
     mMeanLumiOverride = par.lumiMean; // negative value switches off corrections !!!
     mMeanLumiRefOverride = par.lumiMeanRef;
@@ -225,6 +237,18 @@ bool CorrectionMapsLoader::accountCCDBInputs(const ConcreteDataMatcher& matcher,
          canUseCorrections() ? "ON" : "OFF",
          lumiS[scaleType], mMeanLumiOverride, mMeanLumiRefOverride, static_cast<int>(getLumiScaleMode()), mLumiCTPSource, mInstCTPLumiOverride, mInstLumiCTPFactor);
   }
+  if (matcher == ConcreteDataMatcher(o2::header::gDataOriginTPC, "CorrMapSecFluc", 0)) {
+    setCorrMapSecEdgeFluc((o2::gpu::TPCFastTransform*)obj);
+    mCorrMapSecEdgeFluc->rectifyAfterReadingFromFile();
+    setUpdatedMapSecEdgeFluc();
+    return true;
+  }
+  if (matcher == ConcreteDataMatcher(o2::header::gDataOriginTPC, "InfoMapSecFluc", 0)) {
+    LOGP(info, "Updating TPC sector edge fluctuation info");
+    mSecEdgeFlucInfo.setFromTree(*((TTree*)obj));
+    LOGP(info, "Loaded sector edge fluctuation information with {} intervals for {} runs", mSecEdgeFlucInfo.size(), mSecEdgeFlucInfo.getNRuns());
+    return true;
+  }
   return false;
 }
 

Detectors/TPC/calibration/src/CorrectionMapsOptions.cxx

@@ -33,6 +33,7 @@ CorrectionMapsGloOpts CorrectionMapsOptions::parseGlobalOptions(const o2::framew
   tpcopt.lumiMode = static_cast<LumiScaleMode>(lumiModeVal);
 
   tpcopt.enableMShapeCorrection = opts.get<bool>("enable-M-shape-correction");
+  tpcopt.enableSecEdgeFlucCorrection = opts.get<bool>("enable-sec-edge-fluc-correction");
   tpcopt.requestCTPLumi = !opts.get<bool>("disable-ctp-lumi-request");
   tpcopt.checkCTPIDCconsistency = !opts.get<bool>("disable-lumi-type-consistency-check");
   if (!tpcopt.requestCTPLumi && tpcopt.lumiType == LumiScaleType::CTPLumi) {
@@ -49,6 +50,7 @@ void CorrectionMapsOptions::addGlobalOptions(std::vector<ConfigParamSpec>& optio
   addOption(options, ConfigParamSpec{"enable-M-shape-correction", o2::framework::VariantType::Bool, false, {"Enable M-shape distortion correction"}});
   addOption(options, ConfigParamSpec{"disable-ctp-lumi-request", o2::framework::VariantType::Bool, false, {"do not request CTP lumi (regardless what is used for corrections)"}});
   addOption(options, ConfigParamSpec{"disable-lumi-type-consistency-check", o2::framework::VariantType::Bool, false, {"disable check of selected CTP or IDC scaling source being consistent with the map"}});
+  addOption(options, ConfigParamSpec{"enable-sec-edge-fluc-correction", o2::framework::VariantType::Bool, false, {"Enable sector edge fluctuation correction"}});
 }
 
 void CorrectionMapsOptions::addOption(std::vector<ConfigParamSpec>& options, ConfigParamSpec&& osp)