TPC: Add sector edge fluctuation correction infrastructure

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,14 +48,19 @@ 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);
   static void addInput(std::vector<o2::framework::InputSpec>& inputs, o2::framework::InputSpec&& isp);
 
-  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
+  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/include/TPCCalibration/TPCFastSpaceChargeCorrectionHelper.h

@@ -41,6 +41,9 @@ using namespace o2::gpu;
 
 class TPCFastSpaceChargeCorrectionHelper
 {
+ public:
+  using SectorScales = std::array<double, TPCFastTransformGeo::getNumberOfSectors()>;
+
  public:
   /// _____________  Constructors / destructors __________________________
 
@@ -115,15 +118,32 @@ class TPCFastSpaceChargeCorrectionHelper
   /// initialise inverse transformation from linear combination of several input corrections
   void initInverse(std::vector<o2::gpu::TPCFastSpaceChargeCorrection*>& corrections, const std::vector<float>& scaling, bool prn);
 
-  /// merge several corrections
+  /// weighted add of several corrections
   /// \param mainCorrection main correction
   /// \param scale scaling factor for the main correction
   /// \param additionalCorrections vector of pairs of additional corrections and their scaling factors
-  /// \param prn printout flag
   /// \return main correction merged with additional corrections
+  void addCorrections(
+    o2::gpu::TPCFastSpaceChargeCorrection& mainCorrection, double scale,
+    const std::vector<std::pair<const o2::gpu::TPCFastSpaceChargeCorrection*, double>>& additionalCorrections);
+
+  /// weighted add of several corrections with sector-dependent scaling factors
+  /// \param mainCorrection main correction
+  /// \param scale scaling factor for the main correction
+  /// \param additionalCorrections vector of pairs of additional corrections and their sector-dependent scaling factors
+  /// \return main correction merged with additional corrections
+  void addCorrections(
+    o2::gpu::TPCFastSpaceChargeCorrection& mainCorrection, SectorScales scale,
+    const std::vector<std::pair<const o2::gpu::TPCFastSpaceChargeCorrection*, SectorScales>>& additionalCorrections);
+
+  /// merge of two corrections sector-wise
+  /// \param destinationCorrection main correction to which the source correction will be added
+  /// \param sourceCorrection correction to be added to the main correction
+  /// \param sectors vector of sector indices for which the correction will be added
+  /// \return main correction merged with the source correction
   void mergeCorrections(
-    o2::gpu::TPCFastSpaceChargeCorrection& mainCorrection, float scale,
-    const std::vector<std::pair<const o2::gpu::TPCFastSpaceChargeCorrection*, float>>& additionalCorrections, bool prn);
+    o2::gpu::TPCFastSpaceChargeCorrection& destinationCorrection, const o2::gpu::TPCFastSpaceChargeCorrection& sourceCorrection,
+    const std::vector<int>& sectors);
 
   /// how far the voxel mean is allowed to be outside of the voxel (1.1 means 10%)
   void setVoxelMeanValidityRange(double range)

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;
 }