TPC: Adding scaler weights for 1D-distortion fluctuation correction

Author: matthias-kleiner

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,17 @@ 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 getDurationMS() const { return mSamplingTimeMS * mWeights.size(); }
+  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:
@@ -87,6 +98,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 +111,36 @@ 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); }
+
+  /// setting the weights for the scalers
+  void setScalerWeights(const TPCScalerWeights& weights) { mScalerWeights = weights; }
+
+  /// \return returns stored weights for TPC scalers
+  const auto& getScalerWeights() const { return mScalerWeights; }
+
+  /// enable usage of weights
+  void useWeights(bool useWeights) { mUseWeights = useWeights; }
+
+  /// return if weights are used
+  bool weightsUsed() const { return mUseWeights; }
+
  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 mScalerWeights{}; ///< weights for the scalers for A-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

@@ -24,7 +24,8 @@ 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);
@@ -53,7 +54,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 +68,44 @@ 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().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().getWeight(relTSMS);
+    }
+    sum += getScalers(i, side) * weight;
+    sumW += weight;
+  }
+  if (sumW != 0) {
+    return (sum / sumW);
+  }
+  return 0;
+}
+
+float TPCScalerWeights::getWeight(float deltaTime) const
+{
+  const float idxF = (deltaTime - mFirstTimeStampMS) / mSamplingTimeMS;
+  const int idx = idxF;
+  if ((idx < 0) || (idx > mWeights.size() - 1)) {
+    LOGP(error, "Index out of range for deltaTime: {} mFirstTimeStampMS: {} mSamplingTimeMS: {}", deltaTime, mFirstTimeStampMS, mSamplingTimeMS);
+    // set weight 1 to in case it is out of range. This can only happen if the TPC scaler is not valid for given time
+    return 1;
+  }
+
+  if ((idxF == idx) || (idx == mWeights.size() - 1)) {
+    // no interpolation required
+    return mWeights[idx];
+  } else {
+    // interpolate scaler
+    const float y0 = mWeights[idx];
+    const float y1 = mWeights[idx + 1];
+    const float x0 = idx;
+    const float x1 = idx + 1;
+    const float x = idxF;
+    const float y = ((y0 * (x1 - x)) + y1 * (x - x0)) / (x1 - x0);
+    return y;
   }
-  return (sum / nVals);
 }

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,35 @@ 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);
+        mTPCScaler.useWeights(true);
+      }
+    }
+    if (matcher == ConcreteDataMatcher(o2::header::gDataOriginTPC, "TPCSCALERWCCDB", 0)) {
+      LOGP(info, "Updating TPC scaler weights");
+      mScalerWeights = *(TPCScalerWeights*)obj;
+      mTPCScaler.setScalerWeights(mScalerWeights);
+      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 +124,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;
 }