AliceO2Group · davidrohr · Mar 17, 2022 · Mar 4, 2022
@@ -49,6 +49,8 @@ struct CalibdEdxTrackTopologyPolContainer {
   float mQTotMax{200};                                            ///< maximum qTot for which the correction is valid
   std::vector<float> mScalingFactorsqTot{};                       ///< value which is used to scale the result of the polynomial for qTot (can be used for normalization)
   std::vector<float> mScalingFactorsqMax{};                       ///< value which is used to scale the result of the polynomial for qMax (can be used for normalization)
+
+  ClassDefNV(CalibdEdxTrackTopologyPolContainer, 1);
 };
 #endif
 
@@ -206,6 +208,9 @@ class CalibdEdxTrackTopologyPol : public o2::gpu::FlatObject
   /// sets the polynomials from an input file. The names of the objects have to be the same as in the getPolyName() function
   /// \param inpf file where the polynomials are stored
   void setPolynomialsFromFile(TFile& inpf);
+
+  /// setting a default topology correction which just returns 1
+  void setDefaultPolynomials();
 #endif
 
   /// \return returns the name of the polynomial object which can be read in with the setPolynomialsFromFile() function

@@ -131,6 +131,15 @@ void CalibdEdxTrackTopologyPol::construct()
   }
 }
 
+void CalibdEdxTrackTopologyPol::setDefaultPolynomials()
+{
+  for (int i = 0; i < FFits; ++i) {
+    mCalibPolsqTot[i].setParam(0, 1);
+    mCalibPolsqMax[i].setParam(0, 1);
+  }
+  construct();
+}
+
 void CalibdEdxTrackTopologyPol::writeToFile(TFile& outf, const char* name) const
 {
   CalibdEdxTrackTopologyPolContainer cont(mMaxTanTheta, mMaxSinPhi, mThresholdMin, mThresholdMax, mQTotMin, mQTotMax);
@@ -187,6 +196,7 @@ void CalibdEdxTrackTopologyPol::setFromContainer(const CalibdEdxTrackTopologyPol
     mScalingFactorsqTot[i] = container.mScalingFactorsqTot[i];
     mScalingFactorsqMax[i] = container.mScalingFactorsqMax[i];
   }
+  construct();
 }
 
 void CalibdEdxTrackTopologyPol::loadFromFile(const char* fileName, const char* name)
@@ -196,7 +206,6 @@ void CalibdEdxTrackTopologyPol::loadFromFile(const char* fileName, const char* n
   inpf.GetObject(name, polTmp);
   if (polTmp) {
     setFromContainer(*polTmp);
-    construct();
     delete polTmp;
   } else {
     LOGP(info, fmt::format("couldnt load object {} from input file", name));

@@ -49,6 +49,7 @@ enum class CDBType {
   CalGas,             ///< DCS gas measurements
   CalTemperature,     ///< DCS temperature measurements
   CalHV,              ///< DCS HV measurements
+  CalTopologyGain,    ///< Q cluster topology correction
                       ///
   ConfigFEEPad,       ///< FEE pad-by-pad configuration map
                       ///
@@ -78,6 +79,7 @@ const std::unordered_map<CDBType, const std::string> CDBTypeMap{
   {CDBType::CalGas, "TPC/Calib/Gas"},
   {CDBType::CalTemperature, "TPC/Calib/Temperature"},
   {CDBType::CalHV, "TPC/Calib/HV"},
+  {CDBType::CalTopologyGain, "TPC/Calib/TopologyGain"},
   //
   {CDBType::ConfigFEEPad, "TPC/Config/FEEPad"},
   //

@@ -150,6 +150,23 @@ void CalibdEdxContainer::loadSplineTopologyCorrectionFromFile(std::string_view f
   loadTopologyCorrectionFromFile(fileName, mCalibTrackTopologySpline);
 }
 
+void CalibdEdxContainer::setPolTopologyCorrection(const CalibdEdxTrackTopologyPol& calibTrackTopology)
+{
+  setTopologyCorrection(calibTrackTopology, mCalibTrackTopologyPol);
+}
+
+void CalibdEdxContainer::setDefaultPolTopologyCorrection()
+{
+  CalibdEdxTrackTopologyPol calibTrackTopology;
+  calibTrackTopology.setDefaultPolynomials();
+  setTopologyCorrection(calibTrackTopology, mCalibTrackTopologyPol);
+}
+
+void CalibdEdxContainer::setSplineTopologyCorrection(const CalibdEdxTrackTopologySpline& calibTrackTopology)
+{
+  setTopologyCorrection(calibTrackTopology, mCalibTrackTopologySpline);
+}
+
 void CalibdEdxContainer::loadZeroSupresssionThresholdFromFile(std::string_view fileName, std::string_view objName, const float minCorrectionFactor, const float maxCorrectionFactor)
 {
   TFile fInp(fileName.data(), "READ");
@@ -161,10 +178,75 @@ void CalibdEdxContainer::loadZeroSupresssionThresholdFromFile(std::string_view f
 
 void CalibdEdxContainer::setZeroSupresssionThreshold(const CalDet<float>& thresholdMap, const float minCorrectionFactor, const float maxCorrectionFactor)
 {
-  o2::gpu::TPCPadGainCalib thresholdMapTmp(thresholdMap, minCorrectionFactor, maxCorrectionFactor, false);
+  const auto thresholdMapProcessed = processThresholdMap(thresholdMap, maxCorrectionFactor);
+  o2::gpu::TPCPadGainCalib thresholdMapTmp(thresholdMapProcessed, minCorrectionFactor, maxCorrectionFactor, false);
   mThresholdMap = thresholdMapTmp;
 }
 
+CalDet<float> CalibdEdxContainer::processThresholdMap(const CalDet<float>& thresholdMap, const float maxThreshold, const int nPadsInRowCl, const int nPadsInPadCl) const
+{
+  CalDet<float> thresholdMapProcessed(thresholdMap);
+
+  for (int sector = 0; sector < Mapper::NSECTORS; ++sector) {
+    for (int region = 0; region < Mapper::NREGIONS; ++region) {
+      const int maxRow = Mapper::ROWSPERREGION[region] - 1;
+      for (int lrow = 0; lrow <= maxRow; ++lrow) {
+        // find first row of the cluster
+        const int rowStart = std::clamp(lrow - nPadsInRowCl, 0, maxRow);
+        const int rowEnd = std::clamp(lrow + nPadsInRowCl, 0, maxRow);
+        const int addPadsStart = Mapper::ADDITIONALPADSPERROW[region][lrow];
+
+        for (int pad = 0; pad < Mapper::PADSPERROW[region][lrow]; ++pad) {
+          float sumThr = 0;
+          int countThr = 0;
+          // loop ove the rows from the cluster
+          for (int rowCl = rowStart; rowCl <= rowEnd; ++rowCl) {
+            // shift local pad in row in case current row from the cluster has more pads in the row
+            const int addPadsCl = Mapper::ADDITIONALPADSPERROW[region][rowCl];
+            const int diffAddPads = addPadsCl - addPadsStart;
+            const int padClCentre = pad + diffAddPads;
+
+            const int maxPad = Mapper::PADSPERROW[region][rowCl] - 1;
+            const int padStart = std::clamp(padClCentre - nPadsInPadCl, 0, maxPad);
+            const int padEnd = std::clamp(padClCentre + nPadsInPadCl, 0, maxPad);
+            for (int padCl = padStart; padCl <= padEnd; ++padCl) {
+              const int globalPad = Mapper::getGlobalPadNumber(rowCl, padCl, region);
+              // skip for current cluster position as the charge there is not effected from the thresold
+              if (padCl == pad && rowCl == lrow) {
+                continue;
+              }
+
+              float threshold = thresholdMap.getValue(sector, globalPad);
+              if (threshold > maxThreshold) {
+                threshold = maxThreshold;
+              }
+
+              sumThr += threshold;
+              ++countThr;
+            }
+          }
+          const float meanThresold = sumThr / countThr;
+          const int globalPad = Mapper::getGlobalPadNumber(lrow, pad, region);
+          thresholdMapProcessed.setValue(sector, globalPad, meanThresold);
+        }
+      }
+    }
+  }
+  return thresholdMapProcessed;
+}
+
+void CalibdEdxContainer::setGainMap(const CalDet<float>& gainMap, const float minGain, const float maxGain)
+{
+  o2::gpu::TPCPadGainCalib gainMapTmp(gainMap, minGain, maxGain, false);
+  mGainMap = gainMapTmp;
+}
+
+void CalibdEdxContainer::setGainMapResidual(const CalDet<float>& gainMapResidual, const float minResidualGain, const float maxResidualGain)
+{
+  o2::gpu::TPCPadGainCalib gainMapResTmp(gainMapResidual, minResidualGain, maxResidualGain, false);
+  mGainMapResidual = gainMapResTmp;
+}
+
 void CalibdEdxContainer::setDefaultZeroSupresssionThreshold()
 {
   const float defaultVal = getMinZeroSupresssionThreshold() + (getMaxZeroSupresssionThreshold() - getMinZeroSupresssionThreshold()) / 2;
@@ -180,10 +262,15 @@ void CalibdEdxContainer::setDefaultZeroSupresssionThreshold()
 template <class Type>
 void CalibdEdxContainer::loadTopologyCorrectionFromFile(std::string_view fileName, Type*& obj)
 {
-  FlatObject::startConstruction();
-
   // load and set-up container
   Type calibTrackTopologyTmp(fileName.data());
+  setTopologyCorrection(calibTrackTopologyTmp, obj);
+}
+
+template <class Type>
+void CalibdEdxContainer::setTopologyCorrection(const Type& calibTrackTopologyTmp, Type*& obj)
+{
+  FlatObject::startConstruction();
 
   // get size of the flat buffer of the splines
   const std::size_t flatbufferSize = calibTrackTopologyTmp.getFlatBufferSize();

@@ -32,6 +32,22 @@
 namespace o2::tpc
 {
 
+/// flags to set which corrections will be loaded from the CCDB
+enum class CalibsdEdx : unsigned short {
+  CalTopologySpline = 1 << 0,  ///< flag for a topology correction using splines
+  CalTopologyPol = 1 << 1,     ///< flag for a topology correction using polynomials
+  CalThresholdMap = 1 << 2,    ///< flag for using threshold map
+  CalGainMap = 1 << 3,         ///< flag for using the gain map to get the correct cluster charge
+  CalResidualGainMap = 1 << 4, ///< flag for applying residual gain map
+  CalTimeGain = 1 << 5,        ///< flag for residual dE/dx time dependent gain correction
+};
+
+inline CalibsdEdx operator|(CalibsdEdx a, CalibsdEdx b) { return static_cast<CalibsdEdx>(static_cast<int>(a) | static_cast<int>(b)); }
+
+inline CalibsdEdx operator&(CalibsdEdx a, CalibsdEdx b) { return static_cast<CalibsdEdx>(static_cast<int>(a) & static_cast<int>(b)); }
+
+inline CalibsdEdx operator~(CalibsdEdx a) { return static_cast<CalibsdEdx>(~static_cast<int>(a)); }
+
 ///
 /// This container class contains all necessary corrections for the dE/dx
 /// Currently it holds the residual dE/dx correction and the track topology correction, which can be either provided by 2D-splines or by 5D-polynomials.
@@ -107,13 +123,26 @@ class CalibdEdxContainer : public o2::gpu::FlatObject
   /// \param row global pad row
   GPUd() float getZeroSupressionThreshold(const int sector, const gpu::tpccf::Row row, const gpu::tpccf::Pad pad) const { return mThresholdMap.getGainCorrection(sector, row, pad); }
 
+  /// \return returns gain from the full gain map
+  /// \param sector tpc sector
+  /// \param row global pad row
+  GPUd() float getGain(const int sector, const gpu::tpccf::Row row, const gpu::tpccf::Pad pad) const { return mGainMap.getGainCorrection(sector, row, pad); }
+
+  /// \return returns gain from residual gain map
+  /// \param sector tpc sector
+  /// \param row global pad row
+  GPUd() float getResidualGain(const int sector, const gpu::tpccf::Row row, const gpu::tpccf::Pad pad) const { return mGainMapResidual.getGainCorrection(sector, row, pad); }
+
   /// \return returns the residual dE/dx correction for the cluster charge
   /// \param stack ID of the GEM stack
   /// \param charge type of the charge (qMax or qTot)
   /// \param z z position
   /// \param tgl tracking parameter tgl
   GPUd() float getResidualCorrection(const StackID& stack, const ChargeType charge, const float z = 0, const float tgl = 0) const { return mCalibResidualdEdx.getCorrection(stack, charge, z, tgl); }
 
+  /// \return returns if the full gain map will be used during the calculation of the dE/dx to correct the cluster charge
+  GPUd() bool isUsageOfFullGainMap() const { return mApplyFullGainMap; }
+
   /// ========== FlatObject functionality, see FlatObject class for description  =================
 #if !defined(GPUCA_GPUCODE)
   /// cloning a container object (use newFlatBufferPtr=nullptr for simple copy)
@@ -143,32 +172,78 @@ class CalibdEdxContainer : public o2::gpu::FlatObject
   /// \param fileName input file containing the correction
   void loadSplineTopologyCorrectionFromFile(std::string_view fileName);
 
+  /// set the polynomial track topology
+  /// \param calibTrackTopology polynomial track topology correction
+  void setPolTopologyCorrection(const CalibdEdxTrackTopologyPol& calibTrackTopology);
+
+  /// setting a default topology correction which just returns 1
+  void setDefaultPolTopologyCorrection();
+
+  /// set the spline track topology
+  /// \param calibTrackTopology spline track topology correction
+  void setSplineTopologyCorrection(const CalibdEdxTrackTopologySpline& calibTrackTopology);
+
   // loading the residual dE/dx correction from a file
   /// \param fileName input file containing the correction
   void loadResidualCorrectionFromFile(std::string_view fileName) { mCalibResidualdEdx.loadFromFile(fileName); }
 
+  /// setting the residual dEdx correction
+  /// \param residualCorr residual gain calibration object
+  void setResidualCorrection(const CalibdEdxCorrection& residualCorr) { mCalibResidualdEdx = residualCorr; }
+
   // loading the zero supression threshold map from a file
   /// \param fileName input file containing the CalDet map
   void loadZeroSupresssionThresholdFromFile(std::string_view fileName, std::string_view objName, const float minCorrectionFactor, const float maxCorrectionFactor);
 
-  // loading the zero supression threshold map from a file
-  /// \param fileName input file containing the CalDet map
+  /// setting the zero supression threshold map from a CalDet
+  /// \param thresholdMap CalDet containing the zero supression threshold
   void setZeroSupresssionThreshold(const CalDet<float>& thresholdMap) { setZeroSupresssionThreshold(thresholdMap, getMinZeroSupresssionThreshold(), getMaxZeroSupresssionThreshold()); }
 
-  // loading the zero supression threshold map from a file
-  /// \param fileName input file containing the CalDet map
+  /// setting the zero supression threshold map from a CalDet
+  /// \param thresholdMap CalDet containing the zero supression threshold
   void setZeroSupresssionThreshold(const CalDet<float>& thresholdMap, const float minCorrectionFactor, const float maxCorrectionFactor);
 
+  /// setting the gain map map from a CalDet
+  void setGainMap(const CalDet<float>& gainMap, const float minGain, const float maxGain);
+
+  /// setting the gain map map from a CalDet
+  void setGainMapResidual(const CalDet<float>& gainMapResidual, const float minResidualGain = 0.7f, const float maxResidualGain = 1.3f);
+
   /// setting default zero supression threshold map (all values are set to getMinZeroSupresssionThreshold())
   /// \param fileName input file containing the CalDet map
   void setDefaultZeroSupresssionThreshold();
+
+  /// returns status if the spline correction is set
+  bool isTopologyCorrectionSplinesSet() const { return mCalibTrackTopologySpline ? true : false; }
+
+  /// returns status if the polynomials correction is set
+  bool isTopologyCorrectionPolynomialsSet() const { return mCalibTrackTopologyPol ? true : false; }
+
+  /// set loading of a correction from CCDB
+  /// \param calib calibration which will be loaded from CCDB
+  void setCorrectionCCDB(const CalibsdEdx calib) { mCalibsLoad = calib | mCalibsLoad; }
+
+  /// disable loading of a correction from CCDB
+  /// \param calib calibration which will not be loaded from CCDB
+  void disableCorrectionCCDB(const CalibsdEdx calib) { mCalibsLoad = ~calib & mCalibsLoad; }
+
+  /// check if a correction will be loaded from CCDB
+  /// \param calib calibration which will be loaded from CCDB
+  bool isCorrectionCCDB(const CalibsdEdx calib) const { return ((mCalibsLoad & calib) == calib) ? true : false; }
+
+  /// \param applyFullGainMap if set to true the cluster charge will be corrected with the full gain map
+  void setUsageOfFullGainMap(const bool applyFullGainMap) { mApplyFullGainMap = applyFullGainMap; }
 #endif // !GPUCA_GPUCODE
 
  private:
-  CalibdEdxTrackTopologySpline* mCalibTrackTopologySpline{nullptr}; ///< calibration for the track topology correction (splines)
-  CalibdEdxTrackTopologyPol* mCalibTrackTopologyPol{nullptr};       ///< calibration for the track topology correction (polynomial)
-  o2::gpu::TPCPadGainCalib mThresholdMap{};                         ///< calibration object containing the zero supression threshold map
-  CalibdEdxCorrection mCalibResidualdEdx{};                         ///< calibration for the residual dE/dx correction
+  CalibdEdxTrackTopologySpline* mCalibTrackTopologySpline{nullptr};                                                                                                     ///< calibration for the track topology correction (splines)
+  CalibdEdxTrackTopologyPol* mCalibTrackTopologyPol{nullptr};                                                                                                           ///< calibration for the track topology correction (polynomial)
+  o2::gpu::TPCPadGainCalib mThresholdMap{};                                                                                                                             ///< calibration object containing the zero supression threshold map
+  o2::gpu::TPCPadGainCalib mGainMap{};                                                                                                                                  ///< calibration object containing the gain map
+  o2::gpu::TPCPadGainCalib mGainMapResidual{};                                                                                                                          ///< calibration object containing the residual gain map
+  CalibdEdxCorrection mCalibResidualdEdx{};                                                                                                                             ///< calibration for the residual dE/dx correction
+  bool mApplyFullGainMap{false};                                                                                                                                        ///< if set to true the cluster charge will be corrected with the full gain map (when the gain map was not applied during the clusterizer)
+  CalibsdEdx mCalibsLoad{CalibsdEdx::CalTopologyPol | CalibsdEdx::CalThresholdMap | CalibsdEdx::CalGainMap | CalibsdEdx::CalResidualGainMap | CalibsdEdx::CalTimeGain}; ///< flags to set which corrections will be loaded from the CCDB and used during calculation of the dE/dx
 
 #if !defined(GPUCA_GPUCODE)
   template <class Type>
@@ -179,6 +254,9 @@ class CalibdEdxContainer : public o2::gpu::FlatObject
 
   template <class Type>
   void loadTopologyCorrectionFromFile(std::string_view fileName, Type*& obj);
+
+  template <class Type>
+  void setTopologyCorrection(const Type& calibTrackTopologyTmp, Type*& obj);
 #endif
 
   template <class Type>
@@ -190,6 +268,14 @@ class CalibdEdxContainer : public o2::gpu::FlatObject
 #if !defined(GPUCA_GPUCODE)
   template <class Type>
   void cloneFromObject(Type*& obj, const Type* objOld, char* newFlatBufferPtr, const char* oldFlatBufferPtr);
+
+  /// this functions 'smoothes' a CalDet by calculating for each value in a pad the average value using the neighbouring pads, but do not take into account the current pad
+  /// \return returns 'smoothed' CalDet object
+  /// \param thresholdMap zero supression threshold map which will be 'smoothed'
+  /// \param maxThreshold max threshold value which will be considered for averaging
+  /// \param nPadsInRowCl number of pads in row direction which will be taken into account (+- nPadsInRowCl)
+  /// \param nPadsInPadCl number of pads in pad direction which will be taken into account (+- nPadsInPadCl)
+  CalDet<float> processThresholdMap(const CalDet<float>& thresholdMap, const float maxThreshold, const int nPadsInRowCl = 2, const int nPadsInPadCl = 2) const;
 #endif
 
 #ifndef GPUCA_ALIROOT_LIB