Adding TPC scaler weights for 1D-distortion fluctuation correction

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

@@ -78,6 +78,7 @@ enum class CDBType {
                        ///
   CalTimeSeries,       ///< integrated DCAs for longer time interval
   CalScaler,           ///< Scaler from IDCs or combined estimator
+  CalScalerWeights,    ///< Weights for scalers
                        ///
   CorrMapParam,        ///< parameters for CorrectionMapsLoader configuration
                        ///
@@ -142,6 +143,7 @@ const std::unordered_map<CDBType, const std::string> CDBTypeMap{
   // time series
   {CDBType::CalTimeSeries, "TPC/Calib/TimeSeries"},
   {CDBType::CalScaler, "TPC/Calib/Scaler"},
+  {CDBType::CalScalerWeights, "TPC/Calib/ScalerWeights"},
   // correction maps loader params
   {CDBType::CorrMapParam, "TPC/Calib/CorrMapParam"},
   // distortion maps

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

@@ -27,6 +27,16 @@ 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
 */
 
+struct TPCScalerWeights {
+  float getWeight(float deltaTime) const;
+  bool isValid() const { return (mSamplingTimeMS > 0 && !mWeights.empty()); }
+  float mSamplingTimeMS = -1;  ///< sampling of the stored weights
+  float mFirstTimeStampMS = 0; ///< first timestamp
+  std::vector<float> mWeights; ///< stored weights
+
+  ClassDefNV(TPCScalerWeights, 1);
+};
+
 class TPCScaler
 {
  public:
@@ -62,6 +72,19 @@ class TPCScaler
   /// \param file output file
   void dumpToFile(const char* file, const char* name);
 
+  /// dump this object to a file for a given time range
+  /// \param file output file
+  /// \param startTime start time stamp in ms of the dumped object
+  /// \param endTime end time stamp in ms of the dumped object (if endTimeMS=-1 the end of the stored data is used)
+  void dumpToFile(const char* file, const char* name, double startTimeMS, double endTimeMS);
+
+  /// dump this object to several file each containing the scalers for n minutes
+  /// \param file output file
+  /// \param minutesPerObject minutes each output file is covering
+  /// \param marginMS additional margin to the beginning of each object (should be >= getIonDriftTimeMS())
+  /// \param marginCCDBMinutes margin for CCDB timestamp
+  void dumpToFileSlices(const char* file, const char* name, int minutesPerObject, float marginMS = 500, float marginCCDBMinutes = 1);
+
   /// load parameters from input file (which were written using the writeToFile method)
   /// \param inpf input file
   void loadFromFile(const char* inpf, const char* name);
@@ -87,6 +110,9 @@ class TPCScaler
   /// \return return first time in ms
   double getStartTimeStampMS() const { return mTimeStampMS; }
 
+  /// \return return end time in ms
+  double getEndTimeStampMS(o2::tpc::Side side) const { return mTimeStampMS + getDurationMS(side); }
+
   /// \return returns integration time for each scaler value
   float getIntegrationTimeMS() const { return mIntegrationTimeMS; }
 
@@ -97,16 +123,38 @@ class TPCScaler
   /// \param timestamp timestamp for which the last values are used to calculate the mean
   float getMeanScaler(double timestamp, o2::tpc::Side side) const;
 
+  /// \return returns duration in ms for which the scalers are defined
+  float getDurationMS(o2::tpc::Side side) const { return mIntegrationTimeMS * getNValues(side); }
+
+  /// clamp scaler values
+  /// \param minThreshold minimum accepted scaler value
+  /// \param maxThreshold maximum accepted scaler value
+  void clampScalers(float minThreshold, float maxThreshold, Side side);
+
+  /// setting the weights for the scalers
+  void setScalerWeights(const TPCScalerWeights& weights, Side side) { (side == o2::tpc::Side::A) ? mScalerWeightsA = weights : mScalerWeightsC = weights; }
+
+  /// \return returns stored weights for TPC scalers
+  const auto& getScalerWeights(o2::tpc::Side side) const { return (side == o2::tpc::Side::A) ? mScalerWeightsA : mScalerWeightsC; }
+
+  void useWeights(bool useWeights) { mUseWeights = useWeights; }
+
  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);
+  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
+  TPCScalerWeights mScalerWeightsA{}; ///< weights for the scalers for A-side
+  TPCScalerWeights mScalerWeightsC{}; ///< weights for the scalers for C-side
+  bool mUseWeights{false};            ///< use weights when calculating the mean scaler
+
+  // get index to data for given timestamp
+  int getDataIdx(double timestamp) const { return (timestamp - mTimeStampMS) / mIntegrationTimeMS + 0.5; }
+
+  ClassDefNV(TPCScaler, 2);
 };
 
 } // namespace o2::tpc

Detectors/TPC/calibration/src/TPCCalibrationLinkDef.h

@@ -113,4 +113,5 @@
 
 #pragma link C++ class o2::tpc::CalculatedEdx + ;
 #pragma link C++ class o2::tpc::TPCScaler + ;
+#pragma link C++ struct o2::tpc::TPCScalerWeights + ;
 #endif

Detectors/TPC/calibration/src/TPCScaler.cxx

@@ -18,18 +18,81 @@
 #include <TFile.h>
 #include <TTree.h>
 #include "Framework/Logger.h"
+#include "CommonConstants/LHCConstants.h"
 
 using namespace o2::tpc;
 
 void TPCScaler::dumpToFile(const char* file, const char* name)
 {
   TFile out(file, "RECREATE");
-  TTree tree("TPCScaler", "TPCScaler");
+  TTree tree(name, name);
+  tree.SetAutoSave(0);
   tree.Branch("TPCScaler", this);
   tree.Fill();
   out.WriteObject(&tree, name);
 }
 
+void TPCScaler::dumpToFile(const char* file, const char* name, double startTimeMS, double endTimeMS)
+{
+  TPCScaler scaler;
+  scaler.mIonDriftTimeMS = mIonDriftTimeMS;
+  scaler.mRun = mRun;
+  scaler.mFirstTFOrbit = (startTimeMS <= 0) ? mFirstTFOrbit : (mFirstTFOrbit + (startTimeMS - mTimeStampMS) / (o2::constants::lhc::LHCOrbitMUS * 0.001));
+  scaler.mTimeStampMS = (startTimeMS <= 0) ? mTimeStampMS : startTimeMS;
+  scaler.mIntegrationTimeMS = mIntegrationTimeMS;
+
+  const int dataIdx = (startTimeMS <= 0) ? 0 : getDataIdx(startTimeMS);
+  LOGP(info, "dataIdx: {}   max: {}", dataIdx, getNValues(Side::A));
+  const int idxDataStart = (startTimeMS <= 0) ? 0 : dataIdx;
+  const int idxDataEndA = ((endTimeMS < 0) || (dataIdx > getNValues(Side::A))) ? getNValues(Side::A) : dataIdx;
+  const int idxDataEndC = ((endTimeMS < 0) || (dataIdx > getNValues(Side::C))) ? getNValues(Side::C) : dataIdx;
+  scaler.mScalerA = std::vector<float>(mScalerA.begin() + idxDataStart, mScalerA.begin() + idxDataEndA);
+  scaler.mScalerC = std::vector<float>(mScalerC.begin() + idxDataStart, mScalerC.begin() + idxDataEndC);
+  scaler.dumpToFile(file, name);
+}
+
+void TPCScaler::dumpToFileSlices(const char* file, const char* name, int minutesPerObject, float marginMS, float marginCCDBMinutes)
+{
+  if (getNValues(o2::tpc::Side::A) != getNValues(o2::tpc::Side::C)) {
+    LOGP(error, "Number of points stored for A-side and C-side is different");
+    return;
+  }
+  const long marginCCDBMS = marginCCDBMinutes * 60 * 1000;
+  const float msPerObjectTmp = minutesPerObject * 60 * 1000;
+  const int valuesPerSlice = static_cast<int>(msPerObjectTmp / mIntegrationTimeMS);
+  const int marginPerSlice = static_cast<int>(marginMS / mIntegrationTimeMS);
+  int nSlices = getNValues(Side::A) / valuesPerSlice; // number of output objects
+
+  LOGP(info, "Producing {} objects with a CCDB margin of {} ms and {} margin per slice", nSlices, marginCCDBMS, marginPerSlice);
+  if (nSlices == 0) {
+    nSlices = 1;
+  }
+
+  for (int i = 0; i < nSlices; ++i) {
+    const int idxDataStart = (i == 0) ? 0 : (i * valuesPerSlice - marginPerSlice);
+    int idxDataEnd = (i == nSlices - 1) ? getNValues(Side::A) : ((i + 1) * valuesPerSlice + marginPerSlice);
+    if (idxDataEnd > getNValues(Side::A)) {
+      idxDataEnd = getNValues(Side::A);
+    }
+
+    TPCScaler scaler;
+    scaler.mIonDriftTimeMS = mIonDriftTimeMS;
+    scaler.mRun = mRun;
+    scaler.mTimeStampMS = (i == 0) ? mTimeStampMS : (mTimeStampMS + idxDataStart * static_cast<double>(mIntegrationTimeMS));
+    scaler.mFirstTFOrbit = mFirstTFOrbit + (scaler.mTimeStampMS - mTimeStampMS) / (o2::constants::lhc::LHCOrbitMUS * 0.001);
+    scaler.mIntegrationTimeMS = mIntegrationTimeMS;
+    scaler.mScalerA = std::vector<float>(mScalerA.begin() + idxDataStart, mScalerA.begin() + idxDataEnd);
+    scaler.mScalerC = std::vector<float>(mScalerC.begin() + idxDataStart, mScalerC.begin() + idxDataEnd);
+
+    const long timePerSliceMS = valuesPerSlice * mIntegrationTimeMS;
+    const long tsCCDBStart = mTimeStampMS + i * timePerSliceMS;
+    const long tsCCDBStartMargin = (i == 0) ? (tsCCDBStart - marginCCDBMS) : tsCCDBStart;
+    const long tsCCDBEnd = (i == nSlices - 1) ? (getEndTimeStampMS(o2::tpc::Side::A) + marginCCDBMS) : (tsCCDBStart + timePerSliceMS);
+    const std::string fileOut = fmt::format("{}_{}_{}_{}.root", file, i, tsCCDBStartMargin, tsCCDBEnd - 1);
+    scaler.dumpToFile(fileOut.data(), name);
+  }
+}
+
 void TPCScaler::loadFromFile(const char* inpf, const char* name)
 {
   TFile out(inpf, "READ");
@@ -53,7 +116,7 @@ void TPCScaler::setFromTree(TTree& tpcScalerTree)
 float TPCScaler::getMeanScaler(double timestamp, o2::tpc::Side side) const
 {
   // index to data buffer
-  const int idxData = (timestamp - mTimeStampMS) / mIntegrationTimeMS + 0.5;
+  const int idxData = getDataIdx(timestamp);
   const int nVals = getNValuesIonDriftTime();
   const int nValues = getNValues(side);
   if ((nVals == 0) || (nVals > nValues)) {
@@ -67,8 +130,35 @@ float TPCScaler::getMeanScaler(double timestamp, o2::tpc::Side side) const
 
   // sump up values from last ion drift time
   float sum = 0;
+  float sumW = 0;
+  const bool useWeights = mUseWeights && getScalerWeights(side).isValid();
   for (int i = firstIdx; i < lastIdx; ++i) {
-    sum += getScalers(i, side);
+    float weight = 1;
+    if (useWeights) {
+      const double relTSMS = mTimeStampMS + i * mIntegrationTimeMS - timestamp;
+      weight = getScalerWeights(side).getWeight(relTSMS);
+    }
+    sum += getScalers(i, side) * weight;
+    sumW += weight;
   }
-  return (sum / nVals);
+  if (sumW != 0) {
+    return (sum / sumW);
+  }
+  return 0;
+}
+
+void TPCScaler::clampScalers(float minThreshold, float maxThreshold, Side side)
+{
+  auto& scaler = (side == o2::tpc::Side::A) ? mScalerA : mScalerC;
+  std::transform(std::begin(scaler), std::end(scaler), std::begin(scaler), [minThreshold, maxThreshold](auto val) { return std::clamp(val, minThreshold, maxThreshold); });
 }
+
+float TPCScalerWeights::getWeight(float deltaTime) const
+{
+  const int idx = (deltaTime - mFirstTimeStampMS) / mSamplingTimeMS + 0.5;
+  if ((idx < 0) || (idx > mWeights.size() - 1)) {
+    // LOGP(info, "Index out of range for deltaTime: {} mFirstTimeStampMS: {} mSamplingTimeMS: {}", deltaTime, mFirstTimeStampMS, mSamplingTimeMS);
+    return 1;
+  }
+  return mWeights[idx];
+}

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

@@ -19,7 +19,7 @@ namespace o2
 namespace tpc
 {
 
-o2::framework::DataProcessorSpec getTPCScalerSpec();
+o2::framework::DataProcessorSpec getTPCScalerSpec(bool enableWeights);
 
 } // end namespace tpc
 } // end namespace o2

Detectors/TPC/workflow/src/TPCScalerSpec.cxx

@@ -34,7 +34,7 @@ namespace tpc
 class TPCScalerSpec : public Task
 {
  public:
-  TPCScalerSpec(std::shared_ptr<o2::base::GRPGeomRequest> req) : mCCDBRequest(req){};
+  TPCScalerSpec(std::shared_ptr<o2::base::GRPGeomRequest> req, bool enableWeights) : mCCDBRequest(req), mEnableWeights(enableWeights){};
 
   void init(framework::InitContext& ic) final
   {
@@ -49,6 +49,12 @@ class TPCScalerSpec : public Task
       pc.inputs().get<TTree*>("tpcscaler");
     }
 
+    if (mEnableWeights) {
+      if (pc.inputs().isValid("tpcscalerw")) {
+        pc.inputs().get<TPCScalerWeights*>("tpcscalerw");
+      }
+    }
+
     if (pc.services().get<o2::framework::TimingInfo>().runNumber != mTPCScaler.getRun()) {
       LOGP(error, "Run number {} of processed data and run number {} of loaded TPC scaler doesnt match!", pc.services().get<o2::framework::TimingInfo>().runNumber, mTPCScaler.getRun());
     }
@@ -73,19 +79,37 @@ class TPCScalerSpec : public Task
         LOGP(info, "Setting ion drift time to: {}", mIonDriftTimeMS);
         mTPCScaler.setIonDriftTimeMS(mIonDriftTimeMS);
       }
+      if (mScalerWeights.isValid()) {
+        LOGP(info, "Setting TPC scaler weights");
+        mTPCScaler.setScalerWeights(mScalerWeights, Side::A);
+        mTPCScaler.setScalerWeights(mScalerWeights, Side::C);
+        mTPCScaler.useWeights(true);
+      }
+    }
+    if (matcher == ConcreteDataMatcher(o2::header::gDataOriginTPC, "TPCSCALERCCDBW", 0)) {
+      LOGP(info, "Updating TPC scaler weights");
+      mScalerWeights = *(TPCScalerWeights*)obj;
+      mTPCScaler.setScalerWeights(mScalerWeights, Side::A);
+      mTPCScaler.setScalerWeights(mScalerWeights, Side::C);
+      mTPCScaler.useWeights(true);
     }
   }
 
  private:
   std::shared_ptr<o2::base::GRPGeomRequest> mCCDBRequest; ///< info for CCDB request
+  const bool mEnableWeights{false};                       ///< use weights for TPC scalers
+  TPCScalerWeights mScalerWeights{};                      ///< scaler weights
   float mIonDriftTimeMS{-1};                              ///< ion drift time
   TPCScaler mTPCScaler;                                   ///< tpc scaler
 };
 
-o2::framework::DataProcessorSpec getTPCScalerSpec()
+o2::framework::DataProcessorSpec getTPCScalerSpec(bool enableWeights)
 {
   std::vector<InputSpec> inputs;
   inputs.emplace_back("tpcscaler", o2::header::gDataOriginTPC, "TPCSCALERCCDB", 0, Lifetime::Condition, ccdbParamSpec(o2::tpc::CDBTypeMap.at(o2::tpc::CDBType::CalScaler), {}, 1)); // time-dependent
+  if (enableWeights) {
+    inputs.emplace_back("tpcscalerw", o2::header::gDataOriginTPC, "TPCSCALERWCCDB", 0, Lifetime::Condition, ccdbParamSpec(o2::tpc::CDBTypeMap.at(o2::tpc::CDBType::CalScalerWeights), {}, 0)); // non time-dependent
+  }
 
   auto ccdbRequest = std::make_shared<o2::base::GRPGeomRequest>(true,                           // orbitResetTime
                                                                 false,                          // GRPECS=true for nHBF per TF
@@ -102,7 +126,7 @@ o2::framework::DataProcessorSpec getTPCScalerSpec()
     "tpc-scaler",
     inputs,
     outputs,
-    AlgorithmSpec{adaptFromTask<TPCScalerSpec>(ccdbRequest)},
+    AlgorithmSpec{adaptFromTask<TPCScalerSpec>(ccdbRequest, enableWeights)},
     Options{
       {"ion-drift-time", VariantType::Float, -1.f, {"Overwrite ion drift time if a value >0 is provided"}}}};
 }

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

@@ -21,7 +21,9 @@ using namespace o2::framework;
 void customize(std::vector<o2::framework::ConfigParamSpec>& workflowOptions)
 {
   // option allowing to set parameters
-  std::vector<ConfigParamSpec> options{ConfigParamSpec{"configKeyValues", VariantType::String, "", {"Semicolon separated key=value strings"}}};
+  std::vector<ConfigParamSpec> options{
+    ConfigParamSpec{"configKeyValues", VariantType::String, "", {"Semicolon separated key=value strings"}},
+    {"enableWeights", VariantType::Bool, false, {"Enable weights for TPC scalers"}}};
   std::swap(workflowOptions, options);
 }
 
@@ -31,6 +33,7 @@ WorkflowSpec defineDataProcessing(ConfigContext const& config)
 {
   WorkflowSpec workflow;
   o2::conf::ConfigurableParam::updateFromString(config.options().get<std::string>("configKeyValues"));
-  workflow.emplace_back(o2::tpc::getTPCScalerSpec());
+  const bool enableWeights = config.options().get<bool>("enableWeights");
+  workflow.emplace_back(o2::tpc::getTPCScalerSpec(enableWeights));
   return workflow;
 }