AliceO2Group · davidrohr · Feb 11, 2022 · Jan 20, 2022 · Jan 21, 2022 · Jan 24, 2022
@@ -29,7 +29,6 @@ o2_add_library(
           src/TrackCuts.cxx
           src/CalibdEdxCorrection.cxx
           src/CalibdEdxTrackTopologyPol.cxx
-          src/CalibdEdxContainer.cxx
           src/CalibdEdxTrackTopologySpline.cxx
   PUBLIC_LINK_LIBRARIES O2::GPUCommon
                         O2::TPCFastTransformation
@@ -61,7 +60,6 @@ o2_target_root_dictionary(
           include/DataFormatsTPC/LtrCalibData.h
           include/DataFormatsTPC/CalibdEdxCorrection.h
           include/DataFormatsTPC/CalibdEdxTrackTopologyPol.h
-          include/DataFormatsTPC/CalibdEdxContainer.h
           include/DataFormatsTPC/CalibdEdxTrackTopologySpline.h
           include/DataFormatsTPC/BetheBlochAleph.h)
 

@@ -16,6 +16,8 @@
 #define ALICEO2_TPC_CALIBDEDXCORRECTION_H_
 
 #include "GPUCommonDef.h"
+#include "GPUCommonRtypes.h"
+
 #ifndef GPUCA_GPUCODE_DEVICE
 #include <string_view>
 #endif

@@ -12,33 +12,57 @@
 /// \file CalibdEdxTrackTopologyPol.h
 /// \author Matthias Kleiner <mkleiner@ikf.uni-frankfurt.de>
 
-#ifndef ALICEO2_TPC_CalibdEdxTrackTopologyPol_H_
-#define ALICEO2_TPC_CalibdEdxTrackTopologyPol_H_
+#ifndef ALICEO2_TPC_CALIBDEDXTRACKTOPOLOGYPOL_H_
+#define ALICEO2_TPC_CALIBDEDXTRACKTOPOLOGYPOL_H_
 
 #include "GPUCommonRtypes.h"
+#include "MultivariatePolynomial.h"
 #include "GPUCommonDef.h"
+#include "FlatObject.h"
+#include "DataFormatsTPC/Defs.h"
 #ifndef GPUCA_ALIGPUCODE
 #include <string_view>
 #endif
 
-// o2 includes
-#include "DataFormatsTPC/Defs.h"
-
 namespace o2::tpc
 {
 
-class CalibdEdxTrackTopologyPol
+#if !defined(GPUCA_GPUCODE)
+/// simple struct to enable writing the MultivariatePolynomialCT to file
+struct CalibdEdxTrackTopologyPolContainer {
+  /// constructor
+  /// \param maxTheta maximum tanTheta for which the polynomials are valid
+  /// \param maxSinPhi maximum sinPhi for which the polynomials are valid
+  /// \param thresholdMin minimum zero supression threshold for which the polynomials are valid
+  /// \param thresholdMax maximum zero supression threshold for which the polynomials are valid
+  CalibdEdxTrackTopologyPolContainer(const float maxTheta, const float maxSinPhi, const float thresholdMin, const float thresholdMax) : mMaxTanTheta{maxTheta}, mMaxSinPhi{maxSinPhi}, mThresholdMin{thresholdMin}, mThresholdMax{thresholdMax} {};
+
+  /// for ROOT I/O
+  CalibdEdxTrackTopologyPolContainer() = default;
+
+  std::vector<gpu::MultivariatePolynomialContainer> mCalibPols{}; ///< parameters of the polynomial
+  float mMaxTanTheta{2.f};                                        ///< max tanTheta for which the correction is stored
+  float mMaxSinPhi{0.99f};                                        ///< max snp for which the correction is stored
+  float mThresholdMin{2.5f};                                      ///< min zero supression for which the correction is stored
+  float mThresholdMax{5};                                         ///< max zero supression for which the correction is stored
+};
+#endif
+
+/// calibration class for the track topology correction of the dE/dx using multvariate polynomials
+class CalibdEdxTrackTopologyPol : public o2::gpu::FlatObject
 {
  public:
 #if !defined(GPUCA_GPUCODE)
-  CalibdEdxTrackTopologyPol()
-  {
-    clear();
-  }
-  CalibdEdxTrackTopologyPol(std::string_view fileName) { loadFromFile(fileName); }
-#else
-  CalibdEdxTrackTopologyPol() CON_DEFAULT;
+  /// constructor constructs an object Initialized from file
+  /// \param fileName name of the input file containing the object
+  /// \parma name name of the object
+  CalibdEdxTrackTopologyPol(std::string_view fileName, std::string_view name = "CalibdEdxTrackTopologyPol") { loadFromFile(fileName.data(), name.data()); };
 #endif
+
+  /// Default constructor: creates an empty uninitialized object
+  CalibdEdxTrackTopologyPol() CON_DEFAULT;
+
+  /// destructor
   ~CalibdEdxTrackTopologyPol() CON_DEFAULT;
 
   /// \return returns the track topology correction
@@ -49,71 +73,114 @@ class CalibdEdxTrackTopologyPol
   /// \param z z position of the cluster
   /// \param relPad absolute relative pad position of the track
   /// \param relTime relative time position of the track
-  GPUd() float getCorrection(const int region, const ChargeType charge, const float tanTheta, const float sinPhi, const float z, const float relPad, const float relTime) const
+  GPUd() float getCorrection(const int region, const ChargeType charge, const float tanTheta, const float sinPhi, const float z, const float relPad, const float relTime, const float threshold = 0) const
   {
-    const auto& param = mParams[regionIndex(region, charge)];
-    const float x[FXDim]{tanTheta, sinPhi, z, relPad, relTime};
-    const float corr = evalPol4_5D(x, param);
+    const float x[]{tanTheta, sinPhi, z, relPad, relTime, threshold};
+    const float corr = (charge == ChargeType::Tot) ? mCalibPolsqTot[region].eval(x) : mCalibPolsqMax[region].eval(x);
     return corr;
   }
 
-  /// returns the maximum tanTheta for which the splines are valid
+  /// \return returns the track topology correction
+  /// \param region region of the TPC
+  /// \param charge correction for maximum or total charge
+  /// \param x coordinates where the correction is evaluated
+  GPUd() float getCorrection(const int region, const ChargeType charge, const float x[/*inpXdim*/]) const { return (charge == ChargeType::Tot) ? mCalibPolsqTot[region].eval(x) : mCalibPolsqMax[region].eval(x); }
+
+  /// returns the maximum tanTheta for which the polynomials are valid
   GPUd() float getMaxTanTheta() const { return mMaxTanTheta; };
 
-  /// returns the maximum sinPhi for which the splines are valid
+  /// returns the maximum sinPhi for which the polynomials are valid
   GPUd() float getMaxSinPhi() const { return mMaxSinPhi; };
 
-#if !defined(GPUCA_GPUCODE)
-  /// \return returns number of dimensions of the polynomial
-  int getDims() const { return FXDim; }
+  /// returns the minimum zero supression threshold for which the polynomials are valid
+  GPUd() float getMinThreshold() const { return mThresholdMin; };
 
-  /// set the parameters for the polynomials
-  /// \param region region of the TPC
-  /// \param charge correction for maximum or total charge
-  /// \param params parameter for the coefficients
-  void setParams(const int region, const ChargeType charge, const float* params) { std::copy(params, params + FParams, mParams[regionIndex(region, charge)]); }
+  /// returns the maximum zero supression threshold for which the polynomials are valid
+  GPUd() float getMaxThreshold() const { return mThresholdMax; };
 
-  /// \return returns the paramaters of the coefficients
+#if !defined(GPUCA_GPUCODE)
+  /// \return returns polynomial for qTot
   /// \param region region of the TPC
-  /// \param charge correction for maximum or total charge
-  const float* getParams(const int region, const ChargeType charge) const { return mParams[regionIndex(region, charge)]; }
+  const auto& getPolyqTot(const int region) const { return mCalibPolsqTot[region]; }
 
-  /// resetting the parameter
-  void clear();
+  /// \return returns polynomial for qMax
+  /// \param region region of the TPC
+  const auto& getPolyqMax(const int region) const { return mCalibPolsqMax[region]; }
 
-  /// set the maximum tanTheta for which the splines are valid
+  /// set the maximum tanTheta for which the polynomials are valid
   /// \param maxTanTheta maximum tanTheta
   void setMaxTanTheta(const float maxTanTheta) { mMaxTanTheta = maxTanTheta; };
 
-  /// set the maximum sinPhi for which the splines are valid
+  /// set the maximum sinPhi for which the polynomials are valid
   /// \param maxSinPhi maximum sinPhi
   void setMaxSinPhi(const float maxSinPhi) { mMaxSinPhi = maxSinPhi; };
 
-  /// dump the object to a file
-  /// \param fileName name of the output file
-  void saveFile(std::string_view fileName) const;
+  /// set the the minimum zero supression threshold for which the polynomials are valid
+  /// \param thresholdMin minimum threshold
+  void setMinThreshold(const float thresholdMin) { mThresholdMin = thresholdMin; };
+
+  /// set the the maximum zero supression threshold for which the polynomials are valid
+  /// \param thresholdMax maximum threshold
+  void setMaxThreshold(const float thresholdMax) { mThresholdMax = thresholdMax; };
+
+  /// write a class object to the file
+  /// \param outf file where the object will be written to
+  /// \param name name of the object in the output file
+  void writeToFile(TFile& outf, const char* name) const;
+
+  /// init parameters from CalibdEdxTrackTopologyPolContainer
+  /// \param container container for the members
+  void setFromContainer(const CalibdEdxTrackTopologyPolContainer& container);
+
+  /// load members from a file
+  /// \param fileName file where the object will be read from
+  /// \param name name of the object in the output file
+  void loadFromFile(const char* fileName, const char* name);
+
+  /// 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);
+
+  /// \return returns the name of the polynomial object which can be read in with the setPolynomialsFromFile() function
+  /// \param region region of the TPC
+  /// \param charge correction for maximum or total charge
+  static std::string getPolyName(const int region, const ChargeType charge);
+#endif
+
+/// ========== FlatObject functionality, see FlatObject class for description  =================
+#if !defined(GPUCA_GPUCODE)
+  /// cloning a container object (use newFlatBufferPtr=nullptr for simple copy)
+  void cloneFromObject(const CalibdEdxTrackTopologyPol& obj, char* newFlatBufferPtr);
 
-  /// load an object from a file
-  /// \param fileName name of the file
-  void loadFromFile(std::string_view fileName);
+  /// move flat buffer to new location
+  /// \param newBufferPtr new buffer location
+  void moveBufferTo(char* newBufferPtr);
 #endif
 
+  /// destroy the object (release internal flat buffer)
+  void destroy();
+
+  /// set location of external flat buffer
+  void setActualBufferAddress(char* actualFlatBufferPtr);
+
+  /// set future location of the flat buffer
+  void setFutureBufferAddress(char* futureFlatBufferPtr);
+  /// ================================================================================================
+
  private:
-  /// \return returns the index for the stored parameters
-  GPUd() static size_t regionIndex(const int region, const ChargeType charge) { return static_cast<size_t>(region + charge * 10); }
+  constexpr static int FFits{10};                              ///< total number of fits: 10 regions * 2 charge types
+  o2::gpu::MultivariatePolynomial<5, 4> mCalibPolsqTot[FFits]; ///< polynomial objects storage for the polynomials for qTot
+  o2::gpu::MultivariatePolynomial<5, 4> mCalibPolsqMax[FFits]; ///< polynomial objects storage for the polynomials for qMax
+  float mMaxTanTheta{2.f};                                     ///< max tanTheta for which the correction is stored
+  float mMaxSinPhi{0.99f};                                     ///< max snp for which the correction is stored
+  float mThresholdMin{2.5f};                                   ///< min zero supression for which the correction is stored
+  float mThresholdMax{5};                                      ///< max zero supression for which the correction is stored
 
-  /// evaluate the polynyomial for given coordinates and parameters
-  GPUd() static constexpr float evalPol4_5D(const float* x, const float* param)
-  {
-    return param[0] * 1 + param[1] * x[0] + param[2] * x[1] + param[3] * x[2] + param[4] * x[3] + param[5] * x[4] + param[6] * x[0] * x[0] + param[7] * x[0] * x[1] + param[8] * x[0] * x[2] + param[9] * x[0] * x[3] + param[10] * x[0] * x[4] + param[11] * x[1] * x[1] + param[12] * x[1] * x[2] + param[13] * x[1] * x[3] + param[14] * x[1] * x[4] + param[15] * x[2] * x[2] + param[16] * x[2] * x[3] + param[17] * x[2] * x[4] + param[18] * x[3] * x[3] + param[19] * x[3] * x[4] + param[20] * x[4] * x[4] + param[21] * x[0] * x[0] * x[0] + param[22] * x[0] * x[0] * x[1] + param[23] * x[0] * x[0] * x[2] + param[24] * x[0] * x[0] * x[3] + param[25] * x[0] * x[0] * x[4] + param[26] * x[0] * x[1] * x[1] + param[27] * x[0] * x[1] * x[2] + param[28] * x[0] * x[1] * x[3] + param[29] * x[0] * x[1] * x[4] + param[30] * x[0] * x[2] * x[2] + param[31] * x[0] * x[2] * x[3] + param[32] * x[0] * x[2] * x[4] + param[33] * x[0] * x[3] * x[3] + param[34] * x[0] * x[3] * x[4] + param[35] * x[0] * x[4] * x[4] + param[36] * x[1] * x[1] * x[1] + param[37] * x[1] * x[1] * x[2] + param[38] * x[1] * x[1] * x[3] + param[39] * x[1] * x[1] * x[4] + param[40] * x[1] * x[2] * x[2] + param[41] * x[1] * x[2] * x[3] + param[42] * x[1] * x[2] * x[4] + param[43] * x[1] * x[3] * x[3] + param[44] * x[1] * x[3] * x[4] + param[45] * x[1] * x[4] * x[4] + param[46] * x[2] * x[2] * x[2] + param[47] * x[2] * x[2] * x[3] + param[48] * x[2] * x[2] * x[4] + param[49] * x[2] * x[3] * x[3] + param[50] * x[2] * x[3] * x[4] + param[51] * x[2] * x[4] * x[4] + param[52] * x[3] * x[3] * x[3] + param[53] * x[3] * x[3] * x[4] + param[54] * x[3] * x[4] * x[4] + param[55] * x[4] * x[4] * x[4] + param[56] * x[0] * x[0] * x[0] * x[0] + param[57] * x[0] * x[0] * x[0] * x[1] + param[58] * x[0] * x[0] * x[0] * x[2] + param[59] * x[0] * x[0] * x[0] * x[3] + param[60] * x[0] * x[0] * x[0] * x[4] + param[61] * x[0] * x[0] * x[1] * x[1] + param[62] * x[0] * x[0] * x[1] * x[2] + param[63] * x[0] * x[0] * x[1] * x[3] + param[64] * x[0] * x[0] * x[1] * x[4] + param[65] * x[0] * x[0] * x[2] * x[2] + param[66] * x[0] * x[0] * x[2] * x[3] + param[67] * x[0] * x[0] * x[2] * x[4] + param[68] * x[0] * x[0] * x[3] * x[3] + param[69] * x[0] * x[0] * x[3] * x[4] + param[70] * x[0] * x[0] * x[4] * x[4] + param[71] * x[0] * x[1] * x[1] * x[1] + param[72] * x[0] * x[1] * x[1] * x[2] + param[73] * x[0] * x[1] * x[1] * x[3] + param[74] * x[0] * x[1] * x[1] * x[4] + param[75] * x[0] * x[1] * x[2] * x[2] + param[76] * x[0] * x[1] * x[2] * x[3] + param[77] * x[0] * x[1] * x[2] * x[4] + param[78] * x[0] * x[1] * x[3] * x[3] + param[79] * x[0] * x[1] * x[3] * x[4] + param[80] * x[0] * x[1] * x[4] * x[4] + param[81] * x[0] * x[2] * x[2] * x[2] + param[82] * x[0] * x[2] * x[2] * x[3] + param[83] * x[0] * x[2] * x[2] * x[4] + param[84] * x[0] * x[2] * x[3] * x[3] + param[85] * x[0] * x[2] * x[3] * x[4] + param[86] * x[0] * x[2] * x[4] * x[4] + param[87] * x[0] * x[3] * x[3] * x[3] + param[88] * x[0] * x[3] * x[3] * x[4] + param[89] * x[0] * x[3] * x[4] * x[4] + param[90] * x[0] * x[4] * x[4] * x[4] + param[91] * x[1] * x[1] * x[1] * x[1] + param[92] * x[1] * x[1] * x[1] * x[2] + param[93] * x[1] * x[1] * x[1] * x[3] + param[94] * x[1] * x[1] * x[1] * x[4] + param[95] * x[1] * x[1] * x[2] * x[2] + param[96] * x[1] * x[1] * x[2] * x[3] + param[97] * x[1] * x[1] * x[2] * x[4] + param[98] * x[1] * x[1] * x[3] * x[3] + param[99] * x[1] * x[1] * x[3] * x[4] + param[100] * x[1] * x[1] * x[4] * x[4] + param[101] * x[1] * x[2] * x[2] * x[2] + param[102] * x[1] * x[2] * x[2] * x[3] + param[103] * x[1] * x[2] * x[2] * x[4] + param[104] * x[1] * x[2] * x[3] * x[3] + param[105] * x[1] * x[2] * x[3] * x[4] + param[106] * x[1] * x[2] * x[4] * x[4] + param[107] * x[1] * x[3] * x[3] * x[3] + param[108] * x[1] * x[3] * x[3] * x[4] + param[109] * x[1] * x[3] * x[4] * x[4] + param[110] * x[1] * x[4] * x[4] * x[4] + param[111] * x[2] * x[2] * x[2] * x[2] + param[112] * x[2] * x[2] * x[2] * x[3] + param[113] * x[2] * x[2] * x[2] * x[4] + param[114] * x[2] * x[2] * x[3] * x[3] + param[115] * x[2] * x[2] * x[3] * x[4] + param[116] * x[2] * x[2] * x[4] * x[4] + param[117] * x[2] * x[3] * x[3] * x[3] + param[118] * x[2] * x[3] * x[3] * x[4] + param[119] * x[2] * x[3] * x[4] * x[4] + param[120] * x[2] * x[4] * x[4] * x[4] + param[121] * x[3] * x[3] * x[3] * x[3] + param[122] * x[3] * x[3] * x[3] * x[4] + param[123] * x[3] * x[3] * x[4] * x[4] + param[124] * x[3] * x[4] * x[4] * x[4] + param[125] * x[4] * x[4] * x[4] * x[4];
-  }
+#if !defined(GPUCA_GPUCODE)
+  void construct();
+#endif
 
-  static constexpr unsigned short FXDim{5}; ///< number of dimensionality of the polynomial
-  constexpr static int FParams{126};        ///< number of parameters per polynomial
-  constexpr static int FFits{20};           ///< total number of fits: 10 regions * 2 charge types
-  float mParams[FFits][FParams];            ///< paramters of the polynomial
-  float mMaxTanTheta{2.f};                  ///< max tanTheta for which the correction is stored
-  float mMaxSinPhi{0.99f};                  ///< max snp for which the correction is stored
+  ClassDefNV(CalibdEdxTrackTopologyPol, 1);
 };
 
 } // namespace o2::tpc

@@ -157,10 +157,19 @@ class CalibdEdxTrackTopologySpline : public o2::gpu::FlatObject
     return mCalibSplinesqTot[region].interpolate(x);
   };
 
-  GPUd() float getCorrection(const int region, const ChargeType charge, const float tanTheta, const float sinPhi, const float z) const
-  {
-    return charge == ChargeType::Max ? interpolateqMax(region, tanTheta, sinPhi, z) : interpolateqTot(region, tanTheta, sinPhi, z);
-  }
+  /// \return returns the track topology correction
+  /// \param region region of the TPC
+  /// \param charge correction for maximum or total charge
+  /// \param tanTheta local dip angle: z angle - dz/dx
+  /// \param sinPhi track parameter sinphi
+  /// \param z drift length
+  GPUd() float getCorrection(const int region, const ChargeType charge, const float tanTheta, const float sinPhi, const float z) const { return charge == ChargeType::Max ? interpolateqMax(region, tanTheta, sinPhi, z) : interpolateqTot(region, tanTheta, sinPhi, z); }
+
+  /// \return returns the track topology correction
+  /// \param region region of the TPC
+  /// \param charge correction for maximum or total charge
+  /// \param x coordinates where the correction is evaluated
+  GPUd() float getCorrection(const int region, const ChargeType charge, const float x[/*inpXdim*/]) const { return charge == ChargeType::Tot ? mCalibSplinesqTot[region].interpolate(x) : mCalibSplinesqMax[region].interpolate(x); }
 
   /// \param region index of the spline (region)
   /// \return returns the spline for qMax