TPC: Adding v-shape distortions

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
+  CalVShape,           ///< Scaling value for shape distortions on A-side
                        ///
   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::CalVShape, "TPC/Calib/VShape"},
   // correction maps loader params
   {CDBType::CorrMapParam, "TPC/Calib/CorrMapParam"},
   // distortion maps

Detectors/TPC/calibration/CMakeLists.txt

@@ -54,6 +54,7 @@ o2_add_library(TPCCalibration
                        src/CalculatedEdx.cxx
                        src/TPCScaler.cxx
                        src/CorrMapParam.cxx
+                       src/TPCVShapeScaler.cxx
                PUBLIC_LINK_LIBRARIES O2::DataFormatsTPC O2::TPCBase
                                      O2::TPCReconstruction ROOT::Minuit
                                      Microsoft.GSL::GSL
@@ -107,7 +108,8 @@ o2_target_root_dictionary(TPCCalibration
                                   include/TPCCalibration/TPCFastSpaceChargeCorrectionHelper.h
                                   include/TPCCalibration/CalculatedEdx.h
                                   include/TPCCalibration/TPCScaler.h
-                                  include/TPCCalibration/CorrMapParam.h)
+                                  include/TPCCalibration/CorrMapParam.h
+                                  include/TPCCalibration/TPCVShapeScaler.h)
 
 o2_add_test_root_macro(macro/comparePedestalsAndNoise.C
                        PUBLIC_LINK_LIBRARIES O2::TPCBase

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

@@ -40,6 +40,7 @@ namespace tpc
 struct CorrectionMapsLoaderGloOpts {
   int lumiType = 0;
   int lumiMode = 0;
+  bool enableVShapeCorrection = false;
 };
 
 class CorrectionMapsLoader : public o2::gpu::CorrectionMapsHelper

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

@@ -0,0 +1,90 @@
+// 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 TPCVShapeScaler.h
+/// \author Matthias Kleiner <mkleiner@ikf.uni-frankfurt.de>
+
+#ifndef ALICEO2_TPC_TPCVShapeScaler
+#define ALICEO2_TPC_TPCVShapeScaler
+
+#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 TPCVShapeScaler
+{
+ public:
+  /// default constructor
+  TPCVShapeScaler() = default;
+
+  /// default move assignment
+  TPCVShapeScaler& operator=(TPCVShapeScaler&& other) = default;
+
+  /// \return returns number of stored TPC scaler values
+  int getNValues() const { return mScalerA.size(); }
+
+  /// setting the scalers <time, value>
+  void setScaler(const std::vector<std::pair<double, float>>& values) { mScalerA = values; }
+
+  /// \return returns run number for which this object is valid
+  void setRun(int run) { mRun = run; }
+
+  /// 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);
+
+  /// set sampling time of the stored values
+  float setSamplingTimeMS(float t) { return mSamplingTimeMS = t; }
+
+  /// \return returns stored scalers for given side and data index
+  std::pair<double, float> getScalers(unsigned int idx) const { return mScalerA[idx]; }
+
+  /// \return returns stored scalers for given side
+  const auto& getScalers() const { return mScalerA; }
+
+  /// \return returns run number for which this object is valid
+  int getRun() const { return mRun; }
+
+  /// \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 getScaler(const double timestamp) const;
+
+  /// \return returns sampling time of the stored values
+  float getSamplingTimeMS() const { return mSamplingTimeMS; }
+
+ private:
+  int mRun{};                                       ///< run for which this object is valid
+  std::vector<std::pair<double, float>> mScalerA{}; ///< TPC scaler for A-side
+  float mSamplingTimeMS = 1;                        ///< sampling time of the V-shape values
+
+  /// if distance to neighbouring data is larger than sampling time return 0 for the scaling
+  bool checkDeltaTime(double deltaTime) const { return deltaTime < 1.5 * mSamplingTimeMS; }
+
+  ClassDefNV(TPCVShapeScaler, 1);
+};
+
+} // namespace o2::tpc
+#endif

Detectors/TPC/calibration/src/CorrectionMapsLoader.cxx

@@ -105,6 +105,7 @@ void CorrectionMapsLoader::addGlobalOptions(std::vector<ConfigParamSpec>& option
   // these are options which should be added at the workflow level, since they modify the inputs of the devices
   addOption(options, ConfigParamSpec{"lumi-type", o2::framework::VariantType::Int, 0, {"1 = require CTP lumi for TPC correction scaling, 2 = require TPC scalers for TPC correction scaling"}});
   addOption(options, ConfigParamSpec{"corrmap-lumi-mode", o2::framework::VariantType::Int, 0, {"scaling mode: (default) 0 = static + scale * full; 1 = full + scale * derivative; 2 = full + scale * derivative (for MC)"}});
+  addOption(options, ConfigParamSpec{"enable-V-shape-correction", o2::framework::VariantType::Bool, false, {"Enable V-shape distortion correction"}});
 }
 
 //________________________________________________________
@@ -113,6 +114,7 @@ CorrectionMapsLoaderGloOpts CorrectionMapsLoader::parseGlobalOptions(const o2::f
   CorrectionMapsLoaderGloOpts tpcopt;
   tpcopt.lumiType = opts.get<int>("lumi-type");
   tpcopt.lumiMode = opts.get<int>("corrmap-lumi-mode");
+  tpcopt.enableVShapeCorrection = opts.get<bool>("enable-V-shape-correction");
   return tpcopt;
 }
 

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++ class o2::tpc::TPCVShapeScaler + ;
 #endif

Detectors/TPC/calibration/src/TPCVShapeScaler.cxx

@@ -0,0 +1,84 @@
+// 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  TPCVShapeScaler.cxx
+/// \brief Definition of TPCVShapeScaler class
+///
+/// \author Matthias Kleiner <mkleiner@ikf.uni-frankfurt.de>
+
+#include "TPCCalibration/TPCVShapeScaler.h"
+#include <TFile.h>
+#include <TTree.h>
+#include "Framework/Logger.h"
+#include "CommonConstants/LHCConstants.h"
+
+using namespace o2::tpc;
+
+void TPCVShapeScaler::dumpToFile(const char* file, const char* name)
+{
+  TFile out(file, "RECREATE");
+  TTree tree(name, name);
+  tree.SetAutoSave(0);
+  tree.Branch("TPCVShapeScaler", this);
+  tree.Fill();
+  out.WriteObject(&tree, name);
+}
+
+void TPCVShapeScaler::loadFromFile(const char* inpf, const char* name)
+{
+  TFile out(inpf, "READ");
+  TTree* tree = (TTree*)out.Get(name);
+  setFromTree(*tree);
+}
+
+void TPCVShapeScaler::setFromTree(TTree& TPCVShapeScalerTree)
+{
+  TPCVShapeScaler* scalerTmp = this;
+  TPCVShapeScalerTree.SetBranchAddress("TPCVShapeScaler", &scalerTmp);
+  const int entries = TPCVShapeScalerTree.GetEntries();
+  if (entries > 0) {
+    TPCVShapeScalerTree.GetEntry(0);
+  } else {
+    LOGP(error, "TPCVShapeScaler not found in input file");
+  }
+  TPCVShapeScalerTree.SetBranchAddress("TPCVShapeScaler", nullptr);
+}
+
+float TPCVShapeScaler::getScaler(const double timestamp) const
+{
+  auto idx = std::distance(mScalerA.begin(), std::lower_bound(mScalerA.begin(), mScalerA.end(), std::pair<double, float>(timestamp, std::numeric_limits<float>::min())));
+
+  if (idx == mScalerA.size()) {
+    // check end range
+    const double deltaTime = std::abs(mScalerA.back().first - timestamp);
+    if (checkDeltaTime(deltaTime)) {
+      return mScalerA.back().second;
+    }
+  } else if (idx == 0) {
+    const double deltaTime = std::abs(mScalerA.front().first - timestamp);
+    if (checkDeltaTime(deltaTime)) {
+      return mScalerA.front().second;
+    }
+  } else {
+    // check if upper or lower value is closer
+    const double deltaTimeLow = std::abs(mScalerA[idx - 1].first - timestamp);
+    const double deltaTimeUp = std::abs(mScalerA[idx].first - timestamp);
+
+    // return closes value
+    if ((deltaTimeLow < deltaTimeUp) && checkDeltaTime(deltaTimeLow)) {
+      return mScalerA[idx - 1].second;
+    } else if (checkDeltaTime(deltaTimeUp)) {
+      return mScalerA[idx - 1].second;
+    }
+  }
+
+  return 0;
+}

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

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

Detectors/TPC/workflow/src/TPCScalerSpec.cxx

@@ -22,6 +22,7 @@
 #include "TPCBase/CDBInterface.h"
 #include "DetectorsBase/GRPGeomHelper.h"
 #include "TPCCalibration/TPCScaler.h"
+#include "TPCCalibration/TPCVShapeScaler.h"
 #include "TTree.h"
 
 using namespace o2::framework;
@@ -34,7 +35,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 enableVShape) : mCCDBRequest(req), mEnableVShape(enableVShape){};
 
   void init(framework::InitContext& ic) final
   {
@@ -49,6 +50,12 @@ class TPCScalerSpec : public Task
       pc.inputs().get<TTree*>("tpcscaler");
     }
 
+    if (mEnableVShape) {
+      if (pc.inputs().isValid("vshape")) {
+        pc.inputs().get<TTree*>("vshape");
+      }
+    }
+
     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());
     }
@@ -61,6 +68,11 @@ class TPCScalerSpec : public Task
     float meanScaler = (scalerA + scalerC) / 2;
     LOGP(info, "Publishing TPC scaler: {}", meanScaler);
     pc.outputs().snapshot(Output{header::gDataOriginTPC, "TPCSCALER"}, meanScaler);
+
+    if (mEnableVShape) {
+      float vshapescaler = 123;
+      pc.outputs().snapshot(Output{header::gDataOriginTPC, "TPCVShapeScaler"}, vshapescaler);
+    }
   }
 
   void finaliseCCDB(o2::framework::ConcreteDataMatcher& matcher, void* obj) final
@@ -73,19 +85,27 @@ class TPCScalerSpec : public Task
         LOGP(info, "Setting ion drift time to: {}", mIonDriftTimeMS);
         mTPCScaler.setIonDriftTimeMS(mIonDriftTimeMS);
       }
+    } else if (matcher == ConcreteDataMatcher(o2::header::gDataOriginTPC, "VSHAPESCALERCCDB", 0)) {
+      LOGP(info, "Updating V-shape TPC scaler");
+      mVShapeTPCScaler.setFromTree(*((TTree*)obj));
     }
   }
 
  private:
   std::shared_ptr<o2::base::GRPGeomRequest> mCCDBRequest; ///< info for CCDB request
+  const bool mEnableVShape{};                             ///< enable v shape scalers
   float mIonDriftTimeMS{-1};                              ///< ion drift time
   TPCScaler mTPCScaler;                                   ///< tpc scaler
+  TPCVShapeScaler mVShapeTPCScaler;                       ///< TPC V-shape scalers
 };
 
-o2::framework::DataProcessorSpec getTPCScalerSpec()
+o2::framework::DataProcessorSpec getTPCScalerSpec(bool enableVShape)
 {
   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 (enableVShape) {
+    inputs.emplace_back("vshape", o2::header::gDataOriginTPC, "VSHAPESCALERCCDB", 0, Lifetime::Condition, ccdbParamSpec(o2::tpc::CDBTypeMap.at(o2::tpc::CDBType::CalVShape), {}, 1)); // time-dependent
+  }
 
   auto ccdbRequest = std::make_shared<o2::base::GRPGeomRequest>(true,                           // orbitResetTime
                                                                 false,                          // GRPECS=true for nHBF per TF
@@ -102,7 +122,7 @@ o2::framework::DataProcessorSpec getTPCScalerSpec()
     "tpc-scaler",
     inputs,
     outputs,
-    AlgorithmSpec{adaptFromTask<TPCScalerSpec>(ccdbRequest)},
+    AlgorithmSpec{adaptFromTask<TPCScalerSpec>(ccdbRequest, enableVShape)},
     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,10 @@ 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{
+    {"enable-V-shape-correction", VariantType::Bool, false, {"Enable V-shape distortion correction"}},
+    {"configKeyValues", VariantType::String, "", {"Semicolon separated key=value strings"}}};
+
   std::swap(workflowOptions, options);
 }
 
@@ -31,6 +34,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 auto enableVShape = config.options().get<bool>("enable-V-shape-correction");
+  workflow.emplace_back(o2::tpc::getTPCScalerSpec(enableVShape));
   return workflow;
 }