diff --git a/DataFormats/Detectors/TPC/CMakeLists.txt b/DataFormats/Detectors/TPC/CMakeLists.txt index 0567dd094ab1f..af0ca3386e981 100644 --- a/DataFormats/Detectors/TPC/CMakeLists.txt +++ b/DataFormats/Detectors/TPC/CMakeLists.txt @@ -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) diff --git a/DataFormats/Detectors/TPC/include/DataFormatsTPC/CalibdEdxCorrection.h b/DataFormats/Detectors/TPC/include/DataFormatsTPC/CalibdEdxCorrection.h index c61623b0aa7e0..08017427a5585 100644 --- a/DataFormats/Detectors/TPC/include/DataFormatsTPC/CalibdEdxCorrection.h +++ b/DataFormats/Detectors/TPC/include/DataFormatsTPC/CalibdEdxCorrection.h @@ -16,6 +16,8 @@ #define ALICEO2_TPC_CALIBDEDXCORRECTION_H_ #include "GPUCommonDef.h" +#include "GPUCommonRtypes.h" + #ifndef GPUCA_GPUCODE_DEVICE #include #endif diff --git a/DataFormats/Detectors/TPC/include/DataFormatsTPC/CalibdEdxTrackTopologyPol.h b/DataFormats/Detectors/TPC/include/DataFormatsTPC/CalibdEdxTrackTopologyPol.h index 77b351f312636..8a6c030baae30 100644 --- a/DataFormats/Detectors/TPC/include/DataFormatsTPC/CalibdEdxTrackTopologyPol.h +++ b/DataFormats/Detectors/TPC/include/DataFormatsTPC/CalibdEdxTrackTopologyPol.h @@ -12,33 +12,57 @@ /// \file CalibdEdxTrackTopologyPol.h /// \author Matthias Kleiner -#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 #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 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(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 diff --git a/DataFormats/Detectors/TPC/include/DataFormatsTPC/CalibdEdxTrackTopologySpline.h b/DataFormats/Detectors/TPC/include/DataFormatsTPC/CalibdEdxTrackTopologySpline.h index 591a5e35ac22d..2483d7d4f3501 100644 --- a/DataFormats/Detectors/TPC/include/DataFormatsTPC/CalibdEdxTrackTopologySpline.h +++ b/DataFormats/Detectors/TPC/include/DataFormatsTPC/CalibdEdxTrackTopologySpline.h @@ -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 diff --git a/DataFormats/Detectors/TPC/src/CalibdEdxContainer.cxx b/DataFormats/Detectors/TPC/src/CalibdEdxContainer.cxx deleted file mode 100644 index 3fc11496c5e55..0000000000000 --- a/DataFormats/Detectors/TPC/src/CalibdEdxContainer.cxx +++ /dev/null @@ -1,148 +0,0 @@ -// 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 CalibdEdxContainer.cxx -/// \author Matthias Kleiner - -#include "DataFormatsTPC/CalibdEdxContainer.h" - -#if !defined(GPUCA_GPUCODE) && !defined(GPUCA_STANDALONE) -#include "TFile.h" -#include "GPUCommonLogger.h" -#endif - -using namespace GPUCA_NAMESPACE::gpu; -using namespace o2::tpc; - -void CalibdEdxContainer::cloneFromObject(const CalibdEdxContainer& obj, char* newFlatBufferPtr) -{ - const char* oldFlatBufferPtr = obj.mFlatBufferPtr; - FlatObject::cloneFromObject(obj, newFlatBufferPtr); - mCalibResidualdEdx = obj.mCalibResidualdEdx; - if (obj.mCalibTrackTopologyPol) { - mCalibTrackTopologyPol = FlatObject::relocatePointer(oldFlatBufferPtr, mFlatBufferPtr, obj.mCalibTrackTopologyPol); - } - if (obj.mCalibTrackTopologySpline) { - mCalibTrackTopologySpline = FlatObject::relocatePointer(oldFlatBufferPtr, mFlatBufferPtr, obj.mCalibTrackTopologySpline); - } -} - -void CalibdEdxContainer::moveBufferTo(char* newFlatBufferPtr) -{ - char* oldFlatBufferPtr = mFlatBufferPtr; - FlatObject::moveBufferTo(newFlatBufferPtr); - char* currFlatBufferPtr = mFlatBufferPtr; - mFlatBufferPtr = oldFlatBufferPtr; - setActualBufferAddress(currFlatBufferPtr); -} - -void CalibdEdxContainer::destroy() -{ - if (mCalibTrackTopologySpline) { - mCalibTrackTopologySpline->destroy(); - } - mCalibTrackTopologySpline = nullptr; - mCalibTrackTopologyPol = nullptr; - FlatObject::destroy(); -} - -void CalibdEdxContainer::setActualBufferAddress(char* actualFlatBufferPtr) -{ - FlatObject::setActualBufferAddress(actualFlatBufferPtr); - const size_t buffsize = getFlatBufferSize(); - if (buffsize == 0) { - mCalibTrackTopologyPol = nullptr; - mCalibTrackTopologySpline = nullptr; - } else if (buffsize == sizeOfCalibdEdxTrackTopologyPol()) { - // if the size of the buffer is equal to the size of the pol class - // set the pointer to the new location of the buffer - mCalibTrackTopologyPol = reinterpret_cast(mFlatBufferPtr); - } else { - // set the pointer to the new location of the buffer - mCalibTrackTopologySpline = reinterpret_cast(mFlatBufferPtr); - - // set buffer of the spline container class to the correct position - const std::size_t offset = sizeOfCalibdEdxTrackTopologySpline(); - mCalibTrackTopologySpline->setActualBufferAddress(mFlatBufferPtr + offset); - } -} - -void CalibdEdxContainer::setFutureBufferAddress(char* futureFlatBufferPtr) -{ - const size_t buffsize = getFlatBufferSize(); - if (buffsize == 0) { - mCalibTrackTopologyPol = nullptr; - mCalibTrackTopologySpline = nullptr; - } else if (buffsize == sizeOfCalibdEdxTrackTopologyPol()) { - // set member to correct new flat buffer - mCalibTrackTopologyPol = FlatObject::relocatePointer(mFlatBufferPtr, futureFlatBufferPtr, mCalibTrackTopologyPol); - } else { - // set pointer of the spline container to correct new flat buffer - char* distBuffer = FlatObject::relocatePointer(mFlatBufferPtr, futureFlatBufferPtr, mCalibTrackTopologySpline->getFlatBufferPtr()); - mCalibTrackTopologySpline->setFutureBufferAddress(distBuffer); - - // set member to correct new flat buffer - mCalibTrackTopologySpline = FlatObject::relocatePointer(mFlatBufferPtr, futureFlatBufferPtr, mCalibTrackTopologySpline); - } - FlatObject::setFutureBufferAddress(futureFlatBufferPtr); -} - -#if !defined(GPUCA_GPUCODE) && !defined(GPUCA_STANDALONE) && !defined(GPUCA_ALIROOT_LIB) - -void CalibdEdxContainer::loadPolTopologyCorrectionFromFile(std::string_view fileName) -{ - FlatObject::startConstruction(); - CalibdEdxTrackTopologyPol calibTrackTopologyPolTmp(fileName); - - // create mFlatBuffer with correct size - const std::size_t bufferSize = sizeOfCalibdEdxTrackTopologyPol(); - FlatObject::finishConstruction(bufferSize); - - // CalibdEdxTrackTopologyPol* ptrToBuffer = reinterpret_cast(mFlatBufferPtr); // set pointer to flat buffer - // *ptrToBuffer = cal; // deep copy of CalibdEdxTrackTopologyPol to buffer - - // setting member of CalibdEdxTrackTopologyPol to correct buffer address - mCalibTrackTopologyPol = reinterpret_cast(mFlatBufferPtr); - *mCalibTrackTopologyPol = calibTrackTopologyPolTmp; // deep copy of CalibdEdxTrackTopologyPol to buffer -} - -void CalibdEdxContainer::loadSplineTopologyCorrectionFromFile(std::string_view fileName) -{ - FlatObject::startConstruction(); - - // load and set-up spline container - CalibdEdxTrackTopologySpline calibTrackTopologySplineTmp(fileName.data()); - - // get size of the flat buffer of the splines - const std::size_t flatbufferSize = calibTrackTopologySplineTmp.getFlatBufferSize(); - - // size of the dEdx spline container without taking flat buffer into account - const std::size_t objSize = sizeOfCalibdEdxTrackTopologySpline(); - - // create mFlatBuffer with correct size - const std::size_t totalSize = flatbufferSize + objSize; - FlatObject::finishConstruction(totalSize); - - // CalibdEdxTrackTopologySpline* ptrToBuffer = reinterpret_cast(mFlatBufferPtr); // set pointer to flat buffer - // *ptrToBuffer = calibTrackTopologySplineTmp; // deep copy of CalibdEdxTrackTopologyPol to buffer - // delete cal; - - // setting member of CalibdEdxTrackTopologyPol to correct buffer address - mCalibTrackTopologySpline = reinterpret_cast(mFlatBufferPtr); - - // deep copy of CalibdEdxTrackTopologyPol to buffer without moving the flat buffer to correct address - *mCalibTrackTopologySpline = calibTrackTopologySplineTmp; - - // seting the buffer of the splines to current buffer - mCalibTrackTopologySpline->moveBufferTo(objSize + mFlatBufferPtr); -} - -#endif diff --git a/DataFormats/Detectors/TPC/src/CalibdEdxTrackTopologyPol.cxx b/DataFormats/Detectors/TPC/src/CalibdEdxTrackTopologyPol.cxx index 9a324a8d7ddc3..adbc287cf4ccb 100644 --- a/DataFormats/Detectors/TPC/src/CalibdEdxTrackTopologyPol.cxx +++ b/DataFormats/Detectors/TPC/src/CalibdEdxTrackTopologyPol.cxx @@ -11,36 +11,184 @@ #include "DataFormatsTPC/CalibdEdxTrackTopologyPol.h" -#include -#include +#if !defined(GPUCA_GPUCODE) && !defined(GPUCA_STANDALONE) // code invisible on GPU and in the standalone compilation +#include "TFile.h" #include #include - -// root includes -#include "TFile.h" +#endif using namespace o2::tpc; -void CalibdEdxTrackTopologyPol::clear() +#if !defined(GPUCA_GPUCODE) && !defined(GPUCA_STANDALONE) // code invisible on GPU and in the standalone compilation + +void CalibdEdxTrackTopologyPol::cloneFromObject(const CalibdEdxTrackTopologyPol& obj, char* newFlatBufferPtr) +{ + const char* oldFlatBufferPtr = obj.mFlatBufferPtr; + FlatObject::cloneFromObject(obj, newFlatBufferPtr); + + for (unsigned int i = 0; i < FFits; i++) { + char* buffer = FlatObject::relocatePointer(oldFlatBufferPtr, mFlatBufferPtr, obj.mCalibPolsqTot[i].getFlatBufferPtr()); + mCalibPolsqTot[i].cloneFromObject(obj.mCalibPolsqTot[i], buffer); + } + + for (unsigned int i = 0; i < FFits; i++) { + char* buffer = FlatObject::relocatePointer(oldFlatBufferPtr, mFlatBufferPtr, obj.mCalibPolsqMax[i].getFlatBufferPtr()); + mCalibPolsqMax[i].cloneFromObject(obj.mCalibPolsqMax[i], buffer); + } + + mMaxTanTheta = obj.mMaxTanTheta; + mMaxSinPhi = obj.mMaxSinPhi; + mThresholdMin = obj.mThresholdMin; + mThresholdMax = obj.mThresholdMax; +} + +void CalibdEdxTrackTopologyPol::moveBufferTo(char* newFlatBufferPtr) +{ + char* oldFlatBufferPtr = mFlatBufferPtr; + FlatObject::moveBufferTo(newFlatBufferPtr); + char* currFlatBufferPtr = mFlatBufferPtr; + mFlatBufferPtr = oldFlatBufferPtr; + setActualBufferAddress(currFlatBufferPtr); +} + +#endif + +void CalibdEdxTrackTopologyPol::destroy() +{ + for (unsigned int i = 0; i < FFits; i++) { + mCalibPolsqTot[i].destroy(); + mCalibPolsqMax[i].destroy(); + } + FlatObject::destroy(); +} + +void CalibdEdxTrackTopologyPol::setActualBufferAddress(char* actualFlatBufferPtr) +{ + FlatObject::setActualBufferAddress(actualFlatBufferPtr); + size_t offset = 0; + for (unsigned int i = 0; i < FFits; i++) { + offset = alignSize(offset, mCalibPolsqTot[i].getBufferAlignmentBytes()); + mCalibPolsqTot[i].setActualBufferAddress(mFlatBufferPtr + offset); + offset += mCalibPolsqTot[i].getFlatBufferSize(); + } + for (unsigned int i = 0; i < FFits; i++) { + offset = alignSize(offset, mCalibPolsqMax[i].getBufferAlignmentBytes()); + mCalibPolsqMax[i].setActualBufferAddress(mFlatBufferPtr + offset); + offset += mCalibPolsqMax[i].getFlatBufferSize(); + } +} + +void CalibdEdxTrackTopologyPol::setFutureBufferAddress(char* futureFlatBufferPtr) +{ + for (unsigned int i = 0; i < FFits; i++) { + char* buffer = relocatePointer(mFlatBufferPtr, futureFlatBufferPtr, mCalibPolsqTot[i].getFlatBufferPtr()); + mCalibPolsqTot[i].setFutureBufferAddress(buffer); + } + for (unsigned int i = 0; i < FFits; i++) { + char* buffer = relocatePointer(mFlatBufferPtr, futureFlatBufferPtr, mCalibPolsqMax[i].getFlatBufferPtr()); + mCalibPolsqMax[i].setFutureBufferAddress(buffer); + } + FlatObject::setFutureBufferAddress(futureFlatBufferPtr); +} + +#if !defined(GPUCA_GPUCODE) && !defined(GPUCA_STANDALONE) // code invisible on GPU and in the standalone compilation + +void CalibdEdxTrackTopologyPol::construct() +{ + FlatObject::startConstruction(); + + size_t buffSize = 0; + size_t offsets1[FFits]; + size_t offsets2[FFits]; + + for (int index = 0; index < FFits; ++index) { + buffSize = alignSize(buffSize, mCalibPolsqTot[index].getBufferAlignmentBytes()); + offsets1[index] = buffSize; + buffSize += mCalibPolsqTot[index].getFlatBufferSize(); + } + for (int index = 0; index < FFits; ++index) { + buffSize = alignSize(buffSize, mCalibPolsqMax[index].getBufferAlignmentBytes()); + offsets2[index] = buffSize; + buffSize += mCalibPolsqMax[index].getFlatBufferSize(); + } + + FlatObject::finishConstruction(buffSize); + + for (unsigned int i = 0; i < FFits; i++) { + mCalibPolsqTot[i].moveBufferTo(mFlatBufferPtr + offsets1[i]); + } + for (unsigned int i = 0; i < FFits; i++) { + mCalibPolsqMax[i].moveBufferTo(mFlatBufferPtr + offsets2[i]); + } +} + +void CalibdEdxTrackTopologyPol::writeToFile(TFile& outf, const char* name) const +{ + CalibdEdxTrackTopologyPolContainer cont(mMaxTanTheta, mMaxSinPhi, mThresholdMin, mThresholdMax); + cont.mCalibPols.reserve(FFits); + + for (const auto& par : mCalibPolsqTot) { + cont.mCalibPols.emplace_back(par.getContainer()); + } + + for (const auto& par : mCalibPolsqMax) { + cont.mCalibPols.emplace_back(par.getContainer()); + } + + outf.WriteObject(&cont, name); +} + +void CalibdEdxTrackTopologyPol::setFromContainer(const CalibdEdxTrackTopologyPolContainer& container) { - for (auto& row : mParams) { - for (auto& x : row) { - x = 0.f; - } + if (2 * FFits != container.mCalibPols.size()) { + LOGP(info, fmt::format("wrong number of polynomials stored! this {} container {}", 2 * FFits, container.mCalibPols.size())); + return; + } + + for (int i = 0; i < FFits; ++i) { + mCalibPolsqTot[i].setFromContainer(container.mCalibPols[i]); } + + for (int i = 0; i < FFits; ++i) { + mCalibPolsqMax[i].setFromContainer(container.mCalibPols[FFits + i]); + } + + mMaxTanTheta = container.mMaxTanTheta; + mMaxSinPhi = container.mMaxSinPhi; + mThresholdMin = container.mThresholdMin; + mThresholdMax = container.mThresholdMax; } -void CalibdEdxTrackTopologyPol::saveFile(std::string_view fileName) const +void CalibdEdxTrackTopologyPol::loadFromFile(const char* fileName, const char* name) { - std::unique_ptr file(TFile::Open(fileName.data(), "recreate")); - file->WriteObject(this, "CalibdEdxTrackTopologyPol"); + TFile inpf(fileName, "READ"); + CalibdEdxTrackTopologyPolContainer* polTmp = nullptr; + inpf.GetObject(name, polTmp); + if (polTmp) { + setFromContainer(*polTmp); + construct(); + delete polTmp; + } else { + LOGP(info, fmt::format("couldnt load object {} from input file", name)); + } } -void CalibdEdxTrackTopologyPol::loadFromFile(std::string_view fileName) +void CalibdEdxTrackTopologyPol::setPolynomialsFromFile(TFile& inpf) { - std::unique_ptr file(TFile::Open(fileName.data())); - auto tmp = file->Get("CalibdEdxTrackTopologyPol"); - if (tmp != nullptr) { - *this = *tmp; + for (int ireg = 0; ireg < FFits; ++ireg) { + const auto polnameqTot = getPolyName(ireg, ChargeType::Tot); + mCalibPolsqTot[ireg].loadFromFile(inpf, polnameqTot.data()); + const auto polnameqMax = getPolyName(ireg, ChargeType::Max); + mCalibPolsqMax[ireg].loadFromFile(inpf, polnameqMax.data()); } + construct(); } + +std::string CalibdEdxTrackTopologyPol::getPolyName(const int region, const ChargeType charge) +{ + const std::string typeName[2] = {"qMax", "qTot"}; + const std::string polname = fmt::format("spline_{}_region{}", typeName[charge], region).data(); + return polname; +} + +#endif diff --git a/DataFormats/Detectors/TPC/src/DataFormatsTPCLinkDef.h b/DataFormats/Detectors/TPC/src/DataFormatsTPCLinkDef.h index 5e1430e72d9cc..e2230be8abf48 100644 --- a/DataFormats/Detectors/TPC/src/DataFormatsTPCLinkDef.h +++ b/DataFormats/Detectors/TPC/src/DataFormatsTPCLinkDef.h @@ -58,6 +58,6 @@ #pragma link C++ class o2::tpc::CalibdEdxCorrection + ; #pragma link C++ class o2::tpc::CalibdEdxTrackTopologyPol + ; #pragma link C++ class o2::tpc::CalibdEdxTrackTopologySpline + ; -#pragma link C++ class o2::tpc::CalibdEdxContainer + ; +#pragma link C++ struct o2::tpc::CalibdEdxTrackTopologyPolContainer + ; #endif diff --git a/Detectors/TPC/reconstruction/test/testGPUCATracking.cxx b/Detectors/TPC/reconstruction/test/testGPUCATracking.cxx index 73de8a0be0897..9e297984a0406 100644 --- a/Detectors/TPC/reconstruction/test/testGPUCATracking.cxx +++ b/Detectors/TPC/reconstruction/test/testGPUCATracking.cxx @@ -25,10 +25,10 @@ #include "TPCReconstruction/TPCFastTransformHelperO2.h" #include "TPCFastTransform.h" -#include "DataFormatsTPC/CalibdEdxContainer.h" #include "GPUO2Interface.h" #include "GPUO2InterfaceConfiguration.h" #include "TPCPadGainCalib.h" +#include "CalibdEdxContainer.h" using namespace o2::gpu; @@ -76,7 +76,7 @@ BOOST_AUTO_TEST_CASE(CATracking_test1) std::unique_ptr fastTransform(TPCFastTransformHelperO2::instance()->create(0)); config.configCalib.fastTransform = fastTransform.get(); - auto dEdxCalibContainer = std::make_unique(); + auto dEdxCalibContainer = GPUO2Interface::getCalibdEdxContainerDefault(); config.configCalib.dEdxCalibContainer = dEdxCalibContainer.get(); std::unique_ptr gainCalib = GPUO2Interface::getPadGainCalibDefault(); config.configCalib.tpcPadGain = gainCalib.get(); diff --git a/GPU/GPUTracking/CMakeLists.txt b/GPU/GPUTracking/CMakeLists.txt index 23799349d390e..0348b981dbfd6 100644 --- a/GPU/GPUTracking/CMakeLists.txt +++ b/GPU/GPUTracking/CMakeLists.txt @@ -191,6 +191,7 @@ if(ALIGPU_BUILD_TYPE STREQUAL "O2" OR CONFIG_O2_EXTENSIONS) TPCClusterFinder/GPUTPCCFDecodeZS.cxx TPCClusterFinder/GPUTPCCFGather.cxx DataTypes/TPCPadGainCalib.cxx + DataTypes/CalibdEdxContainer.cxx Refit/GPUTrackingRefit.cxx Refit/GPUTrackingRefitKernel.cxx Merger/GPUTPCGMO2Output.cxx) diff --git a/GPU/GPUTracking/DataTypes/CalibdEdxContainer.cxx b/GPU/GPUTracking/DataTypes/CalibdEdxContainer.cxx new file mode 100644 index 0000000000000..11ff98cde4986 --- /dev/null +++ b/GPU/GPUTracking/DataTypes/CalibdEdxContainer.cxx @@ -0,0 +1,208 @@ +// 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 CalibdEdxContainer.cxx +/// \author Matthias Kleiner + +#include "CalibdEdxContainer.h" + +#if !defined(GPUCA_GPUCODE) && !defined(GPUCA_STANDALONE) +#include "TFile.h" +#include "TPCBase/CalDet.h" +#include "Framework/Logger.h" +#include "clusterFinderDefs.h" +#endif + +using namespace GPUCA_NAMESPACE::gpu; +using namespace o2::tpc; + +#if !defined(GPUCA_GPUCODE) && !defined(GPUCA_STANDALONE) +void CalibdEdxContainer::cloneFromObject(const CalibdEdxContainer& obj, char* newFlatBufferPtr) +{ + const char* oldFlatBufferPtr = obj.mFlatBufferPtr; + FlatObject::cloneFromObject(obj, newFlatBufferPtr); + mCalibResidualdEdx = obj.mCalibResidualdEdx; + mThresholdMap = obj.mThresholdMap; + if (obj.mCalibTrackTopologyPol) { + cloneFromObject(mCalibTrackTopologyPol, obj.mCalibTrackTopologyPol, newFlatBufferPtr, oldFlatBufferPtr); + } + if (obj.mCalibTrackTopologySpline) { + cloneFromObject(mCalibTrackTopologySpline, obj.mCalibTrackTopologySpline, newFlatBufferPtr, oldFlatBufferPtr); + } +} + +template +void CalibdEdxContainer::cloneFromObject(Type*& obj, const Type* objOld, char* newFlatBufferPtr, const char* oldFlatBufferPtr) +{ + obj = FlatObject::relocatePointer(oldFlatBufferPtr, mFlatBufferPtr, objOld); + obj->cloneFromObject(*objOld, newFlatBufferPtr); +} +#endif + +void CalibdEdxContainer::moveBufferTo(char* newFlatBufferPtr) +{ + char* oldFlatBufferPtr = mFlatBufferPtr; + FlatObject::moveBufferTo(newFlatBufferPtr); + char* currFlatBufferPtr = mFlatBufferPtr; + mFlatBufferPtr = oldFlatBufferPtr; + setActualBufferAddress(currFlatBufferPtr); +} + +void CalibdEdxContainer::destroy() +{ + if (mCalibTrackTopologySpline) { + mCalibTrackTopologySpline->destroy(); + } + if (mCalibTrackTopologyPol) { + mCalibTrackTopologyPol->destroy(); + } + mCalibTrackTopologySpline = nullptr; + mCalibTrackTopologyPol = nullptr; + FlatObject::destroy(); +} + +void CalibdEdxContainer::setActualBufferAddress(char* actualFlatBufferPtr) +{ + FlatObject::setActualBufferAddress(actualFlatBufferPtr); + if (mCalibTrackTopologyPol) { + setActualBufferAddress(mCalibTrackTopologyPol); + } else if (mCalibTrackTopologySpline) { + setActualBufferAddress(mCalibTrackTopologySpline); + } else { + mCalibTrackTopologyPol = nullptr; + mCalibTrackTopologySpline = nullptr; + } +} + +template +void CalibdEdxContainer::setActualBufferAddress(Type*& obj) +{ + // set the pointer to the new location of the buffer + obj = reinterpret_cast(mFlatBufferPtr); + + // set buffer of the spline container class to the correct position + obj->setActualBufferAddress(mFlatBufferPtr + sizeOfCalibdEdxTrackTopologyObj()); +} + +void CalibdEdxContainer::setFutureBufferAddress(char* futureFlatBufferPtr) +{ + if (mCalibTrackTopologyPol) { + setFutureBufferAddress(mCalibTrackTopologyPol, futureFlatBufferPtr); + } else if (mCalibTrackTopologySpline) { + setFutureBufferAddress(mCalibTrackTopologySpline, futureFlatBufferPtr); + } else { + mCalibTrackTopologyPol = nullptr; + mCalibTrackTopologySpline = nullptr; + } + FlatObject::setFutureBufferAddress(futureFlatBufferPtr); +} + +template +void CalibdEdxContainer::setFutureBufferAddress(Type*& obj, char* futureFlatBufferPtr) +{ + // set pointer of the polynomial container to correct new flat buffer + char* distBuffer = FlatObject::relocatePointer(mFlatBufferPtr, futureFlatBufferPtr, obj->getFlatBufferPtr()); + obj->setFutureBufferAddress(distBuffer); + + // set member to correct new flat buffer + obj = FlatObject::relocatePointer(mFlatBufferPtr, futureFlatBufferPtr, obj); +} + +#if !defined(GPUCA_GPUCODE) && !defined(GPUCA_STANDALONE) + +float CalibdEdxContainer::getMinZeroSupresssionThreshold() const +{ + if (mCalibTrackTopologyPol) { + return mCalibTrackTopologyPol->getMinThreshold(); + } else { + const float minThr = 0; + LOGP(info, "Topology correction not set! Returning defualt min threshold of: {}", minThr); + return minThr; + } +} + +float CalibdEdxContainer::getMaxZeroSupresssionThreshold() const +{ + if (mCalibTrackTopologyPol) { + return mCalibTrackTopologyPol->getMaxThreshold(); + } else { + const float maxThr = 1; + LOGP(info, "Topology correction not set! Returning defualt max threshold of: {}", maxThr); + return maxThr; + } +} + +void CalibdEdxContainer::loadPolTopologyCorrectionFromFile(std::string_view fileName) +{ + loadTopologyCorrectionFromFile(fileName, mCalibTrackTopologyPol); +} + +void CalibdEdxContainer::loadSplineTopologyCorrectionFromFile(std::string_view fileName) +{ + loadTopologyCorrectionFromFile(fileName, mCalibTrackTopologySpline); +} + +void CalibdEdxContainer::loadZeroSupresssionThresholdFromFile(std::string_view fileName, std::string_view objName, const float minCorrectionFactor, const float maxCorrectionFactor) +{ + TFile fInp(fileName.data(), "READ"); + CalDet* threshold = nullptr; + fInp.GetObject(objName.data(), threshold); + setZeroSupresssionThreshold(*threshold, minCorrectionFactor, maxCorrectionFactor); + delete threshold; +} + +void CalibdEdxContainer::setZeroSupresssionThreshold(const CalDet& thresholdMap, const float minCorrectionFactor, const float maxCorrectionFactor) +{ + o2::gpu::TPCPadGainCalib thresholdMapTmp(thresholdMap, minCorrectionFactor, maxCorrectionFactor, false); + mThresholdMap = thresholdMapTmp; +} + +void CalibdEdxContainer::setDefaultZeroSupresssionThreshold() +{ + const float defaultVal = getMinZeroSupresssionThreshold(); + mThresholdMap.setMinCorrectionFactor(defaultVal - 0.1f); + mThresholdMap.setMaxCorrectionFactor(defaultVal + 0.1f); + for (int sector = 0; sector < o2::tpc::constants::MAXSECTOR; ++sector) { + for (unsigned short globPad = 0; globPad < TPC_PADS_IN_SECTOR; ++globPad) { + mThresholdMap.setGainCorrection(sector, globPad, defaultVal); + } + } +} + +template +void CalibdEdxContainer::loadTopologyCorrectionFromFile(std::string_view fileName, Type*& obj) +{ + FlatObject::startConstruction(); + + // load and set-up container + Type calibTrackTopologyTmp(fileName.data()); + + // get size of the flat buffer of the splines + const std::size_t flatbufferSize = calibTrackTopologyTmp.getFlatBufferSize(); + + // size of the dEdx container without taking flat buffer into account + const std::size_t objSize = sizeOfCalibdEdxTrackTopologyObj(); + + // create mFlatBuffer with correct size + const std::size_t totalSize = flatbufferSize + objSize; + FlatObject::finishConstruction(totalSize); + + // setting member of CalibdEdxTrackTopologyPol to correct buffer address + obj = reinterpret_cast(mFlatBufferPtr); + + // deep copy of CalibdEdxTrackTopologyPol to buffer without moving the flat buffer to correct address + obj->cloneFromObject(calibTrackTopologyTmp, nullptr); + + // seting the buffer of the splines to current buffer + obj->moveBufferTo(objSize + mFlatBufferPtr); +} + +#endif diff --git a/DataFormats/Detectors/TPC/include/DataFormatsTPC/CalibdEdxContainer.h b/GPU/GPUTracking/DataTypes/CalibdEdxContainer.h similarity index 62% rename from DataFormats/Detectors/TPC/include/DataFormatsTPC/CalibdEdxContainer.h rename to GPU/GPUTracking/DataTypes/CalibdEdxContainer.h index 21d56d9fb6edf..b19ce9a15e805 100644 --- a/DataFormats/Detectors/TPC/include/DataFormatsTPC/CalibdEdxContainer.h +++ b/GPU/GPUTracking/DataTypes/CalibdEdxContainer.h @@ -23,6 +23,7 @@ #include "DataFormatsTPC/CalibdEdxTrackTopologyPol.h" #include "DataFormatsTPC/CalibdEdxTrackTopologySpline.h" #include "FlatObject.h" +#include "TPCPadGainCalib.h" #ifndef GPUCA_ALIGPUCODE #include @@ -61,9 +62,18 @@ class CalibdEdxContainer : public o2::gpu::FlatObject /// \param z z position /// \param relPad relative pad position of the cluster /// \param relTime relative time position of the cluster - GPUd() float getTopologyCorrection(const int region, const ChargeType charge, const float tanTheta, const float sinPhi, const float z, const float relPad, const float relTime) const + GPUd() float getTopologyCorrection(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 { - return mCalibTrackTopologyPol ? mCalibTrackTopologyPol->getCorrection(region, charge, tanTheta, sinPhi, z, relPad, relTime) : (mCalibTrackTopologySpline ? mCalibTrackTopologySpline->getCorrection(region, charge, tanTheta, sinPhi, z) : getDefaultTopologyCorrection(tanTheta, sinPhi)); + return mCalibTrackTopologyPol ? mCalibTrackTopologyPol->getCorrection(region, charge, tanTheta, sinPhi, z, relPad, relTime, threshold) : (mCalibTrackTopologySpline ? mCalibTrackTopologySpline->getCorrection(region, charge, tanTheta, sinPhi, z) : getDefaultTopologyCorrection(tanTheta, sinPhi)); + } + + /// \return returns the topology correction for the cluster charge + /// \param region region of the TPC + /// \param charge type of the charge (qMax or qTot) + /// \param x coordinates where the correction is evaluated + GPUd() float getTopologyCorrection(const int region, const ChargeType charge, const float x[]) const + { + return mCalibTrackTopologyPol ? mCalibTrackTopologyPol->getCorrection(region, charge, x) : (mCalibTrackTopologySpline ? mCalibTrackTopologySpline->getCorrection(region, charge, x) : getDefaultTopologyCorrection(x[0], x[1])); } /// \return returns analytical default correction @@ -76,6 +86,19 @@ class CalibdEdxContainer : public o2::gpu::FlatObject /// \return returns maximum sinPhi for which the topology correction is valid GPUd() float getMaxSinPhiTopologyCorrection() const { return mCalibTrackTopologyPol ? mCalibTrackTopologyPol->getMaxSinPhi() : (mCalibTrackTopologySpline ? mCalibTrackTopologySpline->getMaxSinPhi() : 1); } +#if !defined(GPUCA_GPUCODE) + /// \returns the minimum zero supression threshold for which the track topology correction is valid + float getMinZeroSupresssionThreshold() const; + + /// \returns the maximum zero supression threshold for which the track topology correction is valid + float getMaxZeroSupresssionThreshold() const; +#endif + + /// \return returns zero supression threshold + /// \param sector tpc sector + /// \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 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) @@ -105,30 +128,61 @@ class CalibdEdxContainer : public o2::gpu::FlatObject #if !defined(GPUCA_GPUCODE) && !defined(GPUCA_STANDALONE) // loading the polynomial track topology correction from a file - /// \param fileName input file containg the correction + /// \param fileName input file containing the correction void loadPolTopologyCorrectionFromFile(std::string_view fileName); // loading the spline track topology correction from a file - /// \param fileName input file containg the correction + /// \param fileName input file containing the correction void loadSplineTopologyCorrectionFromFile(std::string_view fileName); // loading the residual dE/dx correction from a file - /// \param fileName input file containg the correction + /// \param fileName input file containing the correction void loadResidualCorrectionFromFile(std::string_view fileName) { mCalibResidualdEdx.loadFromFile(fileName); } + + // 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 + void setZeroSupresssionThreshold(const CalDet& thresholdMap) { setZeroSupresssionThreshold(thresholdMap, getMinZeroSupresssionThreshold(), getMaxZeroSupresssionThreshold()); } + + // loading the zero supression threshold map from a file + /// \param fileName input file containing the CalDet map + void setZeroSupresssionThreshold(const CalDet& thresholdMap, const float minCorrectionFactor, const float maxCorrectionFactor); + + /// setting default zero supression threshold map (all values are set to getMinZeroSupresssionThreshold()) + /// \param fileName input file containing the CalDet map + void setDefaultZeroSupresssionThreshold(); #endif // !GPUCA_GPUCODE private: CalibdEdxTrackTopologySpline* mCalibTrackTopologySpline{nullptr}; ///< calibration for the track topology correction (splines) CalibdEdxTrackTopologyPol* mCalibTrackTopologyPol{nullptr}; ///< calibration for the track topology correction (polynomial) - CalibdEdxCorrection mCalibResidualdEdx; ///< calibration for the residual dE/dx correction + o2::gpu::TPCPadGainCalib mThresholdMap{}; ///< calibration object containing the zero supression threshold map + CalibdEdxCorrection mCalibResidualdEdx{}; ///< calibration for the residual dE/dx correction #if !defined(GPUCA_GPUCODE) - /// \return returns size of the CalibdEdxTrackTopologyPol class - std::size_t sizeOfCalibdEdxTrackTopologyPol() const { return alignSize(sizeof(CalibdEdxTrackTopologyPol), FlatObject::getClassAlignmentBytes()); } + template + std::size_t sizeOfCalibdEdxTrackTopologyObj() const + { + return alignSize(sizeof(Type), FlatObject::getClassAlignmentBytes()); + } - /// \return returns size of the CalibdEdxTrackTopologySpline class (without taking the size of the flat buffer into acocunt) - std::size_t sizeOfCalibdEdxTrackTopologySpline() const { return alignSize(sizeof(CalibdEdxTrackTopologyPol), FlatObject::getClassAlignmentBytes()); } -#endif // !GPUCA_GPUCODE + template + void loadTopologyCorrectionFromFile(std::string_view fileName, Type*& obj); +#endif + + template + void setActualBufferAddress(Type*& obj); + + template + void setFutureBufferAddress(Type*& obj, char* futureFlatBufferPtr); + +#if !defined(GPUCA_GPUCODE) + template + void cloneFromObject(Type*& obj, const Type* objOld, char* newFlatBufferPtr, const char* oldFlatBufferPtr); +#endif #ifndef GPUCA_ALIROOT_LIB ClassDefNV(CalibdEdxContainer, 1); diff --git a/GPU/GPUTracking/DataTypes/TPCPadGainCalib.cxx b/GPU/GPUTracking/DataTypes/TPCPadGainCalib.cxx index d3eff4254f3c0..d516e3fac89c3 100644 --- a/GPU/GPUTracking/DataTypes/TPCPadGainCalib.cxx +++ b/GPU/GPUTracking/DataTypes/TPCPadGainCalib.cxx @@ -31,11 +31,23 @@ TPCPadGainCalib::TPCPadGainCalib() } TPCPadGainCalib::TPCPadGainCalib(const o2::tpc::CalDet& gainMap) : TPCPadGainCalib() +{ + setFromMap(gainMap); +} + +TPCPadGainCalib::TPCPadGainCalib(const o2::tpc::CalDet& gainMap, const float minValue, const float maxValue, const bool inv) : TPCPadGainCalib() +{ + setMinCorrectionFactor(minValue); + setMaxCorrectionFactor(maxValue); + setFromMap(gainMap, inv); +} + +void TPCPadGainCalib::setFromMap(const o2::tpc::CalDet& gainMap, const bool inv) { for (int sector = 0; sector < o2::tpc::constants::MAXSECTOR; sector++) { for (int p = 0; p < TPC_PADS_IN_SECTOR; p++) { const float gainVal = gainMap.getValue(sector, p); - mGainCorrection[sector].set(p, (gainVal > 1.e-5f) ? 1.f / gainVal : 1.f); + inv ? mGainCorrection[sector].set(p, (gainVal > 1.e-5f) ? 1.f / gainVal : 1.f) : mGainCorrection[sector].set(p, gainVal); } } } diff --git a/GPU/GPUTracking/DataTypes/TPCPadGainCalib.h b/GPU/GPUTracking/DataTypes/TPCPadGainCalib.h index 5912e042c6c9d..a06a1a0bf71d9 100644 --- a/GPU/GPUTracking/DataTypes/TPCPadGainCalib.h +++ b/GPU/GPUTracking/DataTypes/TPCPadGainCalib.h @@ -46,6 +46,16 @@ struct TPCPadGainCalib { #ifndef GPUCA_GPUCODE TPCPadGainCalib(); TPCPadGainCalib(const o2::tpc::CalDet&); + + /// constructor + /// \param minValue minimum value which will be stored + /// \param maxValue maximum value which will be stored + /// \param inv setting the inverse value + TPCPadGainCalib(const o2::tpc::CalDet&, const float minValue, const float maxValue, const bool inv); + + /// setting the stored values from CalDet + /// \param inv setting the inverse value + void setFromMap(const o2::tpc::CalDet&, const bool inv = true); #endif // Deal with pad gain correction from here on @@ -54,6 +64,11 @@ struct TPCPadGainCalib { mGainCorrection[sector].set(globalPad(row, pad), c); } + GPUdi() void setGainCorrection(int sector, unsigned short globalPad, float c) + { + mGainCorrection[sector].set(globalPad, c); + } + GPUdi() float getGainCorrection(int sector, tpccf::Row row, tpccf::Pad pad) const { return mGainCorrection[sector].get(globalPad(row, pad)); @@ -64,14 +79,28 @@ struct TPCPadGainCalib { return mPadOffsetPerRow[row] + pad; } + GPUdi() void setMinCorrectionFactor(const float minCorrectionFactor) + { + for (int sector = 0; sector < GPUCA_NSLICES; sector++) { + mGainCorrection[sector].mMinCorrectionFactor = minCorrectionFactor; + } + } + + GPUdi() void setMaxCorrectionFactor(const float maxCorrectionFactor) + { + for (int sector = 0; sector < GPUCA_NSLICES; sector++) { + mGainCorrection[sector].mMaxCorrectionFactor = maxCorrectionFactor; + } + } + private: template class SectorPadGainCorrection { public: - constexpr static float MinCorrectionFactor = 0.f; - constexpr static float MaxCorrectionFactor = 2.f; + float mMinCorrectionFactor = 0.f; + float mMaxCorrectionFactor = 2.f; constexpr static int NumOfSteps = TPCPadGainCorrectionStepNum::value; GPUdi() SectorPadGainCorrection() @@ -97,18 +126,18 @@ struct TPCPadGainCalib { } private: - GPUd() static T pack(float f) + GPUd() T pack(float f) const { - f = CAMath::Clamp(f, MinCorrectionFactor, MaxCorrectionFactor); - f -= MinCorrectionFactor; + f = CAMath::Clamp(f, mMinCorrectionFactor, mMaxCorrectionFactor); + f -= mMinCorrectionFactor; f *= float(NumOfSteps); - f /= (MaxCorrectionFactor - MinCorrectionFactor); + f /= (mMaxCorrectionFactor - mMinCorrectionFactor); return CAMath::Nint(f); } - GPUd() static float unpack(T c) + GPUd() float unpack(T c) const { - return MinCorrectionFactor + (MaxCorrectionFactor - MinCorrectionFactor) * float(c) / float(NumOfSteps); + return mMinCorrectionFactor + (mMaxCorrectionFactor - mMinCorrectionFactor) * float(c) / float(NumOfSteps); } T mGainCorrection[TPC_PADS_IN_SECTOR]; diff --git a/GPU/GPUTracking/Definitions/GPUSettingsList.h b/GPU/GPUTracking/Definitions/GPUSettingsList.h index 5d6139674f524..b79f78059a0d8 100644 --- a/GPU/GPUTracking/Definitions/GPUSettingsList.h +++ b/GPU/GPUTracking/Definitions/GPUSettingsList.h @@ -410,6 +410,7 @@ AddOption(dEdxPolTopologyCorrFile, std::string, "", "", 0, "File name of the dE/ AddOption(transformationFile, std::string, "", "", 0, "File name of TPC fast transformation map") AddOption(matLUTFile, std::string, "", "", 0, "File name of material LUT file") AddOption(gainCalibFile, std::string, "", "", 0, "File name of TPC pad gain calibration") +AddOption(thresholdCalibFile, std::string, "", "", 0, "File name of TPC zero supression threshold map") AddOption(allocateOutputOnTheFly, bool, true, "", 0, "Allocate shm output buffers on the fly, instead of using preallocated buffer with upper bound size") AddOption(outputBufferSize, unsigned long, 200000000ul, "", 0, "Size of the output buffers to be allocated") AddOption(mutexMemReg, bool, false, "", 0, "Global mutex to serialize GPU memory registration") diff --git a/GPU/GPUTracking/Global/GPUChainTracking.cxx b/GPU/GPUTracking/Global/GPUChainTracking.cxx index 42eb814fabec9..9ecdf60a2563e 100644 --- a/GPU/GPUTracking/Global/GPUChainTracking.cxx +++ b/GPU/GPUTracking/Global/GPUChainTracking.cxx @@ -45,7 +45,7 @@ #include "GPUHostDataTypes.h" #include "GPUTPCCFChainContext.h" #include "GPUTrackingRefit.h" -#include "DataFormatsTPC/CalibdEdxContainer.h" +#include "CalibdEdxContainer.h" #else #include "GPUO2FakeClasses.h" #endif diff --git a/GPU/GPUTracking/Global/GPUChainTrackingIO.cxx b/GPU/GPUTracking/Global/GPUChainTrackingIO.cxx index abb37d673a889..e9f33335314e5 100644 --- a/GPU/GPUTracking/Global/GPUChainTrackingIO.cxx +++ b/GPU/GPUTracking/Global/GPUChainTrackingIO.cxx @@ -41,7 +41,7 @@ #include "DataFormatsTPC/ZeroSuppression.h" #include "GPUHostDataTypes.h" #include "DataFormatsTPC/Digit.h" -#include "DataFormatsTPC/CalibdEdxContainer.h" +#include "CalibdEdxContainer.h" #else #include "GPUO2FakeClasses.h" #endif diff --git a/GPU/GPUTracking/Interface/GPUO2Interface.cxx b/GPU/GPUTracking/Interface/GPUO2Interface.cxx index 96484b0950cc5..681965cce4324 100644 --- a/GPU/GPUTracking/Interface/GPUO2Interface.cxx +++ b/GPU/GPUTracking/Interface/GPUO2Interface.cxx @@ -22,6 +22,7 @@ #include "GPUQA.h" #include "GPUOutputControl.h" #include "TPCPadGainCalib.h" +#include "CalibdEdxContainer.h" #include #include @@ -170,6 +171,11 @@ std::unique_ptr GPUO2Interface::getPadGainCalib(const o2::tpc:: return std::make_unique(in); } +std::unique_ptr GPUO2Interface::getCalibdEdxContainerDefault() +{ + return std::make_unique(); +} + int GPUO2Interface::UpdateCalibration(const GPUCalibObjectsConst& newCalib) { return 1; diff --git a/GPU/GPUTracking/Interface/GPUO2Interface.h b/GPU/GPUTracking/Interface/GPUO2Interface.h index 19d77b09a28d2..859a1caa5ac9a 100644 --- a/GPU/GPUTracking/Interface/GPUO2Interface.h +++ b/GPU/GPUTracking/Interface/GPUO2Interface.h @@ -68,6 +68,8 @@ class GPUO2Interface static std::unique_ptr getPadGainCalibDefault(); static std::unique_ptr getPadGainCalib(const o2::tpc::CalDet& in); + static std::unique_ptr getCalibdEdxContainerDefault(); + int registerMemoryForGPU(const void* ptr, size_t size); int unregisterMemoryForGPU(const void* ptr); diff --git a/GPU/GPUTracking/Merger/GPUTPCGMTrackParam.cxx b/GPU/GPUTracking/Merger/GPUTPCGMTrackParam.cxx index f1fa07ca93fbc..5c825d9b73757 100644 --- a/GPU/GPUTracking/Merger/GPUTPCGMTrackParam.cxx +++ b/GPU/GPUTracking/Merger/GPUTPCGMTrackParam.cxx @@ -318,21 +318,20 @@ GPUd() bool GPUTPCGMTrackParam::Fit(GPUTPCGMMerger* GPUrestrict() merger, int iT prop.SetTrack(this, prop.GetAlpha()); } if (merger->Param().par.dodEdx && iWay == nWays - 1 && clusters[ihit].leg == clusters[maxN - 1].leg) { - float qtot, qmax, relPad, relTime; + float qtot, qmax, pad, relTime; if (merger->GetConstantMem()->ioPtrs.clustersNative == nullptr) { qtot = clustersXYZ[ihit].amp; qmax = 0; - relPad = 0; + pad = 0; relTime = 0; } else { const ClusterNative& cl = merger->GetConstantMem()->ioPtrs.clustersNative->clustersLinear[clusters[ihit].num]; qtot = cl.qTot; qmax = cl.qMax; - relPad = cl.getPad() - int(cl.getPad() + 0.5f); + pad = cl.getPad(); relTime = cl.getTime() - int(cl.getTime() + 0.5f); - //zPos = std::abs(std::abs(cl.getTime()) * 0.199606f * 2.58 - 250.f); // std::abs(time * eleParam.ZbinWidth * gasParam.DriftV - zMaxTPC); } - dEdx.fillCluster(qtot, qmax, clusters[ihit].row, clusters[ihit].slice, mP[2], mP[3], param, merger->GetConstantMem()->calibObjects, zz, relPad, relTime); + dEdx.fillCluster(qtot, qmax, clusters[ihit].row, clusters[ihit].slice, mP[2], mP[3], param, merger->GetConstantMem()->calibObjects, zz, pad, relTime); } } else if (retVal == 2) { // cluster far away form the track if (allowModification) { diff --git a/GPU/GPUTracking/Standalone/CMakeLists.txt b/GPU/GPUTracking/Standalone/CMakeLists.txt index 1d10e75f75197..08cb8392800bc 100644 --- a/GPU/GPUTracking/Standalone/CMakeLists.txt +++ b/GPU/GPUTracking/Standalone/CMakeLists.txt @@ -191,7 +191,7 @@ if(CONFIG_O2_EXTENSIONS) target_sources(standalone_support PRIVATE ${O2_DIR}/Common/Field/src/MagFieldFast.cxx ${O2_DIR}/DataFormats/Detectors/TPC/src/CompressedClusters.cxx - ${O2_DIR}/DataFormats/Detectors/TPC/src/CalibdEdxContainer.cxx + ${O2_DIR}/GPU/GPUTracking/DataTypes/CalibdEdxContainer.cxx ${O2_DIR}/DataFormats/Detectors/TPC/src/CalibdEdxTrackTopologySpline.cxx ${O2_DIR}/DataFormats/simulation/src/MCCompLabel.cxx ${O2_DIR}/DataFormats/Reconstruction/src/TrackParametrization.cxx diff --git a/GPU/GPUTracking/dEdx/GPUdEdx.h b/GPU/GPUTracking/dEdx/GPUdEdx.h index 5c21d33c10849..3e101e4101329 100644 --- a/GPU/GPUTracking/dEdx/GPUdEdx.h +++ b/GPU/GPUTracking/dEdx/GPUdEdx.h @@ -22,7 +22,7 @@ #include "GPUdEdxInfo.h" #if defined(GPUCA_HAVE_O2HEADERS) && !defined(GPUCA_OPENCL1) #include "DataFormatsTPC/Defs.h" -#include "DataFormatsTPC/CalibdEdxContainer.h" +#include "CalibdEdxContainer.h" #endif namespace GPUCA_NAMESPACE @@ -35,7 +35,7 @@ class GPUdEdx { public: GPUd() void clear() {} - GPUd() void fillCluster(float qtot, float qmax, int padRow, unsigned char slice, float trackSnp, float trackTgl, const GPUParam& param, const GPUCalibObjectsConst& calib, float z, float relPad, float relTime) {} + GPUd() void fillCluster(float qtot, float qmax, int padRow, unsigned char slice, float trackSnp, float trackTgl, const GPUParam& param, const GPUCalibObjectsConst& calib, float z, float pad, float relTime) {} GPUd() void fillSubThreshold(int padRow, const GPUParam& param) {} GPUd() void computedEdx(GPUdEdxInfo& output, const GPUParam& param) {} }; @@ -47,7 +47,7 @@ class GPUdEdx public: // The driver must call clear(), fill clusters row by row outside-in, then run computedEdx() to get the result GPUd() void clear(); - GPUd() void fillCluster(float qtot, float qmax, int padRow, unsigned char slice, float trackSnp, float trackTgl, const GPUParam& param, const GPUCalibObjectsConst& calib, float z, float relPad, float relTime); + GPUd() void fillCluster(float qtot, float qmax, int padRow, unsigned char slice, float trackSnp, float trackTgl, const GPUParam& param, const GPUCalibObjectsConst& calib, float z, float pad, float relTime); GPUd() void fillSubThreshold(int padRow, const GPUParam& param); GPUd() void computedEdx(GPUdEdxInfo& output, const GPUParam& param); @@ -107,7 +107,7 @@ GPUdi() void GPUdEdx::checkSubThresh(int roc) mLastROC = roc; } -GPUdnii() void GPUdEdx::fillCluster(float qtot, float qmax, int padRow, unsigned char slice, float trackSnp, float trackTgl, const GPUParam& GPUrestrict() param, const GPUCalibObjectsConst& calib, float z, float relPad, float relTime) +GPUdnii() void GPUdEdx::fillCluster(float qtot, float qmax, int padRow, unsigned char slice, float trackSnp, float trackTgl, const GPUParam& GPUrestrict() param, const GPUCalibObjectsConst& calib, float z, float pad, float relTime) { if (mCount >= MAX_NCL) { return; @@ -140,11 +140,12 @@ GPUdnii() void GPUdEdx::fillCluster(float qtot, float qmax, int padRow, unsigned } // getting the topology correction - const float absRelPad = CAMath::Abs(relPad); + const float absRelPad = CAMath::Abs(pad - int(pad + 0.5f)); const int region = param.tpcGeometry.GetRegion(padRow); z = CAMath::Abs(z); - const float qMaxTopologyCorr = calibContainer->getTopologyCorrection(region, o2::tpc::ChargeType::Max, tanTheta, snp, z, absRelPad, relTime); - const float qTotTopologyCorr = calibContainer->getTopologyCorrection(region, o2::tpc::ChargeType::Tot, tanTheta, snp, z, absRelPad, relTime); + const float threshold = calibContainer->getZeroSupressionThreshold(slice, padRow, pad); // TODO: Use the mean zero supresion threshold of all pads in the cluster? + const float qMaxTopologyCorr = calibContainer->getTopologyCorrection(region, o2::tpc::ChargeType::Max, tanTheta, snp, z, absRelPad, relTime, threshold); + const float qTotTopologyCorr = calibContainer->getTopologyCorrection(region, o2::tpc::ChargeType::Tot, tanTheta, snp, z, absRelPad, relTime, threshold); qmax /= qMaxTopologyCorr; qtot /= qTotTopologyCorr; diff --git a/GPU/TPCFastTransformation/CMakeLists.txt b/GPU/TPCFastTransformation/CMakeLists.txt index 0ff6d8f39640a..4eab98e5b5a11 100644 --- a/GPU/TPCFastTransformation/CMakeLists.txt +++ b/GPU/TPCFastTransformation/CMakeLists.txt @@ -29,6 +29,8 @@ set(SRCS TPCFastTransformGeo.cxx TPCFastSpaceChargeCorrection.cxx TPCFastTransform.cxx + MultivariatePolynomial.cxx + MultivariatePolynomialHelper.cxx ) string(REPLACE ".cxx" ".h" HDRS_CINT_O2 "${SRCS}") @@ -58,6 +60,14 @@ if(${ALIGPU_BUILD_TYPE} STREQUAL "O2") COMPONENT_NAME GPU LABELS gpu) + o2_add_test(MultivarPolynomials + COMPONENT_NAME GPU + PUBLIC_LINK_LIBRARIES O2::${MODULE} + SOURCES test/testMultivarPolynomials.cxx + ENVIRONMENT O2_ROOT=${CMAKE_BINARY_DIR}/stage + LABELS gpu + CONFIGURATIONS RelWithDebInfo Release MinRelSize) + foreach(m SplineDemo.C fastTransformQA.C diff --git a/GPU/TPCFastTransformation/MultivariatePolynomial.cxx b/GPU/TPCFastTransformation/MultivariatePolynomial.cxx new file mode 100644 index 0000000000000..b1ffe616fb65e --- /dev/null +++ b/GPU/TPCFastTransformation/MultivariatePolynomial.cxx @@ -0,0 +1,15 @@ +// 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 MultivariatePolynomial.cxx +/// \author Matthias Kleiner + +#include "MultivariatePolynomial.h" diff --git a/GPU/TPCFastTransformation/MultivariatePolynomial.h b/GPU/TPCFastTransformation/MultivariatePolynomial.h new file mode 100644 index 0000000000000..814fffbc54c5d --- /dev/null +++ b/GPU/TPCFastTransformation/MultivariatePolynomial.h @@ -0,0 +1,240 @@ +// 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 MultivariatePolynomial.h +/// \author Matthias Kleiner + +#ifndef ALICEO2_TPC_MULTIVARIATEPOLYNOMIAL +#define ALICEO2_TPC_MULTIVARIATEPOLYNOMIAL + +#include "GPUCommonDef.h" +#include "FlatObject.h" +#include "MultivariatePolynomialHelper.h" + +#if !defined(GPUCA_GPUCODE) +#include +#include +#include "TFile.h" +#include "Framework/Logger.h" +#endif + +namespace GPUCA_NAMESPACE::gpu +{ + +/// Class for multivariate polynomials. +/// The parameters of the coefficients have to be provided as input and can be obtained from TLinear fitter or from sklearn (PolynomialFeatures) etc. +/// The evaluation of the polynomials can be speed up by providing the dimensions and degree during compile time! +/// +/// Usage: see example in testMultivarPolynomials.cxx +/// Dim > 0 && Degree > 0 : the number of dimensions and the degree is known at compile time +/// Dim = 0 && Degree = 0 : the number of dimensions and the degree will be set during runtime +template +class MultivariatePolynomial : public FlatObject, public MultivariatePolynomialHelper +{ + public: +#if !defined(GPUCA_GPUCODE) + /// constructor for runtime evaluation of polynomial formula + /// \param nDim number of dimensions + /// \param degree degree of the polynomial + template ::type = 0> + MultivariatePolynomial(const unsigned int nDim, const unsigned int degree) : MultivariatePolynomialHelper{nDim, degree}, mNParams{this->getNParameters(degree, nDim)} + { + construct(); + } + + /// constructor for compile time evaluation of polynomial formula + template ::type = 0> + MultivariatePolynomial() : mNParams{this->getNParameters(Degree, Dim)} + { + construct(); + } +#else + /// default constructor + MultivariatePolynomial() CON_DEFAULT; +#endif + + /// default destructor + ~MultivariatePolynomial() CON_DEFAULT; + + /// Copy constructor + MultivariatePolynomial(const MultivariatePolynomial& obj) { this->cloneFromObject(obj, nullptr); } + + /// ========== FlatObject functionality, see FlatObject class for description ================= +#if !defined(GPUCA_GPUCODE) + /// cloning a container object (use newFlatBufferPtr=nullptr for simple copy) + void cloneFromObject(const MultivariatePolynomial& obj, char* newFlatBufferPtr); + + /// 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); + /// ================================================================================================ + + /// evaluates the polynomial for given coordinates + /// \param x query coordinates + GPUd() float eval(const float x[/*Dim*/]) const { return this->evalPol(mParams, x); } + +#if !defined(GPUCA_GPUCODE) + /// \return returns number of parameters of the polynomials + unsigned int getNParams() const { return mNParams; } + + /// set the parameters for the coefficients of the polynomial + /// \param params parameter for the coefficients + void setParams(const float params[/*mNParams*/]) { std::copy(params, params + mNParams, mParams); } + + /// \return returns the paramaters of the coefficients + const float* getParams() const { return mParams; } + + /// load parameters from input file (which were written using the writeToFile method) + /// \param inpf input file + /// \parma name name of the object in the file + void loadFromFile(TFile& inpf, const char* name); + + /// write parameters to file + /// \param outf output file + /// \param name name of the output object + void writeToFile(TFile& outf, const char* name); + + /// converts the parameters to a container which can be written to a root file + MultivariatePolynomialContainer getContainer() const { return MultivariatePolynomialContainer{this->getDim(), this->getDegree(), mNParams, mParams}; } + + /// set the parameters from MultivariatePolynomialContainer + /// \param container container for the parameters + void setFromContainer(const MultivariatePolynomialContainer& container); +#endif + + private: + using DataTParams = float; ///< data type of the parameters of the polynomials + unsigned int mNParams{}; ///< number of parameters of the polynomial + DataTParams* mParams{nullptr}; ///< parameters of the coefficients of the polynomial + +#if !defined(GPUCA_GPUCODE) + /// \return returns the size of the parameters + std::size_t sizeOfParameters() const { return mNParams * sizeof(DataTParams); } + + // construct the object (flatbuffer) + void construct(); +#endif +}; + +//================================================================================= +//============================ inline implementations ============================= +//================================================================================= + +#if !defined(GPUCA_GPUCODE) +template +void MultivariatePolynomial::loadFromFile(TFile& inpf, const char* name) +{ + MultivariatePolynomialContainer* polTmp = nullptr; + inpf.GetObject(name, polTmp); + if (polTmp) { + setFromContainer(*polTmp); + delete polTmp; + } else { + LOGP(info, fmt::format("couldnt load object {} from input file", name)); + } +} + +template +void MultivariatePolynomial::setFromContainer(const MultivariatePolynomialContainer& container) +{ + if constexpr (Dim > 0 && Degree > 0) { + if (this->getDim() != container.mDim) { + LOGP(info, fmt::format("wrong number of dimensions! this {} container {}", this->getDim(), container.mDim)); + return; + } + if (this->getDegree() != container.mDegree) { + LOGP(info, fmt::format("wrong number of degrees! this {} container {}", this->getDegree(), container.mDegree)); + return; + } + setParams(container.mParams.data()); + } else { + MultivariatePolynomial polTmp(container.mDim, container.mDegree); + polTmp.setParams(container.mParams.data()); + this->cloneFromObject(polTmp, nullptr); + } +} + +template +void MultivariatePolynomial::writeToFile(TFile& outf, const char* name) +{ + const MultivariatePolynomialContainer cont = getContainer(); + outf.WriteObject(&cont, name); +} + +template +void MultivariatePolynomial::cloneFromObject(const MultivariatePolynomial& obj, char* newFlatBufferPtr) +{ + const char* oldFlatBufferPtr = obj.mFlatBufferPtr; + FlatObject::cloneFromObject(obj, newFlatBufferPtr); + mNParams = obj.mNParams; + if constexpr (Dim == 0 && Degree == 0) { + this->mDim = obj.mDim; + this->mDegree = obj.mDegree; + } + if (obj.mParams) { + mParams = FlatObject::relocatePointer(oldFlatBufferPtr, mFlatBufferPtr, obj.mParams); + } +} + +template +void MultivariatePolynomial::moveBufferTo(char* newFlatBufferPtr) +{ + char* oldFlatBufferPtr = mFlatBufferPtr; + FlatObject::moveBufferTo(newFlatBufferPtr); + char* currFlatBufferPtr = mFlatBufferPtr; + mFlatBufferPtr = oldFlatBufferPtr; + setActualBufferAddress(currFlatBufferPtr); +} + +template +void MultivariatePolynomial::construct() +{ + FlatObject::startConstruction(); + const std::size_t flatbufferSize = sizeOfParameters(); + FlatObject::finishConstruction(flatbufferSize); + mParams = reinterpret_cast(mFlatBufferPtr); +} +#endif + +template +void MultivariatePolynomial::destroy() +{ + mParams = nullptr; + FlatObject::destroy(); +} + +template +void MultivariatePolynomial::setActualBufferAddress(char* actualFlatBufferPtr) +{ + FlatObject::setActualBufferAddress(actualFlatBufferPtr); + mParams = reinterpret_cast(mFlatBufferPtr); +} + +template +void MultivariatePolynomial::setFutureBufferAddress(char* futureFlatBufferPtr) +{ + mParams = FlatObject::relocatePointer(mFlatBufferPtr, futureFlatBufferPtr, mParams); + FlatObject::setFutureBufferAddress(futureFlatBufferPtr); +} + +} // namespace GPUCA_NAMESPACE::gpu + +#endif diff --git a/GPU/TPCFastTransformation/MultivariatePolynomialHelper.cxx b/GPU/TPCFastTransformation/MultivariatePolynomialHelper.cxx new file mode 100644 index 0000000000000..db209cfbdcfae --- /dev/null +++ b/GPU/TPCFastTransformation/MultivariatePolynomialHelper.cxx @@ -0,0 +1,45 @@ +// 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 MultivariatePolynomialHelper.cxx +/// \author Matthias Kleiner + +#include "MultivariatePolynomialHelper.h" + +#if !defined(GPUCA_GPUCODE) && !defined(GPUCA_STANDALONE) +#include "Framework/Logger.h" +#endif + +using namespace GPUCA_NAMESPACE::gpu; + +#if !defined(GPUCA_GPUCODE) && !defined(GPUCA_STANDALONE) +void MultivariatePolynomialHelper<0, 0>::print() const +{ + const auto terms = getTerms(); + std::string formula = ""; + for (int i = 0; i < terms.size() - 1; ++i) { + formula += fmt::format("{} + ", terms[i]); + } + formula += terms.back(); + LOGP(info, formula); +} + +std::vector MultivariatePolynomialHelper<0, 0>::getTerms() const +{ + std::vector terms{"par[0]"}; + int indexPar = 1; + for (unsigned int deg = 1; deg <= mDegree; ++deg) { + const auto strTmp = combination_with_repetiton>(deg, mDim, nullptr, indexPar, nullptr); + terms.insert(terms.end(), strTmp.begin(), strTmp.end()); + } + return terms; +} +#endif diff --git a/GPU/TPCFastTransformation/MultivariatePolynomialHelper.h b/GPU/TPCFastTransformation/MultivariatePolynomialHelper.h new file mode 100644 index 0000000000000..59889d3e1ad96 --- /dev/null +++ b/GPU/TPCFastTransformation/MultivariatePolynomialHelper.h @@ -0,0 +1,319 @@ +// 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 MultivariatePolynomialHelper.h +/// \author Matthias Kleiner + +#ifndef ALICEO2_TPC_MULTIVARIATEPOLYNOMIALHELPER +#define ALICEO2_TPC_MULTIVARIATEPOLYNOMIALHELPER + +#include "GPUCommonDef.h" + +#if !defined(GPUCA_GPUCODE) && !defined(GPUCA_STANDALONE) +#include +#include +#include +#include +#endif + +namespace GPUCA_NAMESPACE::gpu +{ + +#if !defined(GPUCA_GPUCODE) + +/// simple struct to enable writing the MultivariatePolynomial to file +struct MultivariatePolynomialContainer { + + /// constructor + /// \param dim number of dimensions of the polynomial + /// \param degree degree of the polynomials + /// \param nParameters number of parameters + /// \param params parmaeters + MultivariatePolynomialContainer(const unsigned int dim, const unsigned int degree, const unsigned int nParameters, const float params[/* nParameters*/]) : mDim{dim}, mDegree{degree}, mParams{params, params + nParameters} {}; + + /// for ROOT I/O + MultivariatePolynomialContainer() = default; + + const unsigned int mDim{}; ///< number of dimensions of the polynomial + const unsigned int mDegree{}; ///< degree of the polynomials + const std::vector mParams{}; ///< parameters of the polynomial +}; +#endif + +/// Helper class for calculating the number of parameters for a multidimensional polynomial +class MultivariatePolynomialParametersHelper +{ + public: + /// \returns number of parameters for given dimension and degree of polynomials + /// calculates the number of parameters for a multivariate polynomial for given degree: nParameters = (n+d-1 d) -> binomial coefficient + /// see: https://mathoverflow.net/questions/225953/number-of-polynomial-terms-for-certain-degree-and-certain-number-of-variables + GPUd() static constexpr unsigned int getNParameters(const unsigned int degree, const unsigned int dim) { return (degree == 0) ? binomialCoeff(dim - 1, 0) : binomialCoeff(dim - 1 + degree, degree) + getNParameters(degree - 1, dim); } + + private: + /// calculate factorial of n + /// \return returns n! + GPUd() static constexpr unsigned int factorial(const unsigned int n) { return (n == 0) || (n == 1) ? 1 : n * factorial(n - 1); } + + /// calculates binomial coefficient + /// \return returns (n k) + GPUd() static constexpr unsigned int binomialCoeff(const unsigned int n, const unsigned int k) { return factorial(n) / (factorial(k) * factorial(n - k)); } +}; + +/// Helper struct for evaluating a multidimensional polynomial using compile time evaluated formula +/// Compile time method to extract the formula is obtained from run time method (check combination_with_repetiton() and evalPol()) +/// by performing all loops during compile time and replacing the array to keep track of the dimensions for given term (pos[FMaxdegree + 1]) +/// to a simple unsigned int called Pos where each digit represents the dimension for a given term e.g. pos = 2234 -> x[2]*x[2]*x[3]*x[4] +/// +template +class MultivariatePolynomialHelper : public MultivariatePolynomialParametersHelper +{ + static constexpr unsigned short FMaxdim = 10; ///< maximum dimensionality of the polynomials (number of different digits: 0,1,2,3....9 ) + static constexpr unsigned short FMaxdegree = 9; ///< maximum degree of the polynomials (maximum number of digits in unsigned integer - 1) + +#if !defined(GPUCA_GPUCODE) + static_assert(Dim <= MultivariatePolynomialHelper::FMaxdim && Degree <= MultivariatePolynomialHelper::FMaxdegree, "Max. number of dimensions or degrees exceeded!"); +#endif + + public: + /// evaluates the polynomial for given parameters and coordinates + /// \param par parameters of the polynomials + /// \param x input coordinates + GPUd() static constexpr float evalPol(GPUgeneric() const float par[/*number of parameters*/], const float x[/*number of dimensions*/]) { return par[0] + loopDegrees<1>(par, x); } + + /// \return returns number of dimensions of the polynomials + GPUd() static constexpr unsigned int getDim() { return Dim; } + + /// \return returns the degree of the polynomials + GPUd() static constexpr unsigned int getDegree() { return Degree; } + + private: + /// computes power of 10 + GPUd() static constexpr unsigned int pow10(const unsigned int n) { return n == 0 ? 1 : 10 * pow10(n - 1); } + + /// helper for modulo to extract the digit in an integer a at position b (can be obtained with pow10(digitposition)): e.g. a=1234 b=pow10(2)=100 -> returns 2 + GPUd() static constexpr unsigned int mod10(const unsigned int a, const unsigned int b) { return (a / b) % 10; } + + /// resetting digits of pos for given position to refDigit + GPUd() static constexpr unsigned int resetIndices(const unsigned int degreePol, const unsigned int pos, const unsigned int leftDigit, const unsigned int iter, const unsigned int refDigit); + + GPUd() static constexpr unsigned int getNewPos(const unsigned int degreePol, const unsigned int pos, const unsigned int digitPos); + + /// calculates term e.g. x^3*y + /// \tparam DegreePol max degree of the polynomials + /// \pos decoded information about the current term e.g. 1233 -> x[1]*x[2]*x[3]*x[3] (otherwise an array could be used) + template + GPUd() static constexpr float prodTerm(const float x[], const unsigned int pos); + + /// calculate sum of the terms for given degree -> summation of the par[]*x^4 + par[]*x^3*y + par[]*x^3*z + par[]*x^2*y^2..... terms + /// \tparam DegreePol max degree of the polynomials + /// \Pos decoded information about the current term e.g. 1233 -> x[1]*x[2]*x[3]*x[3] (otherwise an array could be used) + /// \tparam Index index for accessing the parameters + template + GPUd() static constexpr float sumTerms(GPUgeneric() const float par[], const float x[]); + + /// loop over the degrees of the polynomials (for formula see https://math.stackexchange.com/questions/1234240/equation-that-defines-multi-dimensional-polynomial) + /// \tparam degree iteration of the loop which starts from 1 to the max degree of the polynomial e.g. Iter=4 -> summation of the par[]*x^4 + par[]*x^3*y + par[]*x^3*z + par[]*x^2*y^2..... terms + /// \param par parameters of the pokynomial + /// \param x input coordinates + template + GPUd() static constexpr float loopDegrees(GPUgeneric() const float par[], const float x[]); +}; + +/// Helper struct for evaluating a multidimensional polynomial using run time evaluated formula +template <> +class MultivariatePolynomialHelper<0, 0> : public MultivariatePolynomialParametersHelper +{ + public: +#if !defined(GPUCA_GPUCODE) + /// constructor + /// \param nDim dimensionality of the polynomials + /// \param degree degree of the polynomials + MultivariatePolynomialHelper(const unsigned int nDim, const unsigned int degree) : mDim{nDim}, mDegree{degree} { assert(mDegree <= FMaxdegree); }; +#endif + + /// default constructor + MultivariatePolynomialHelper() CON_DEFAULT; + + /// Destructor + ~MultivariatePolynomialHelper() CON_DEFAULT; + +#if !defined(GPUCA_GPUCODE) + /// printing the formula of the polynomial + void print() const; + + /// \return returns the terms which are used to evaluate the polynomial + std::vector getTerms() const; +#endif + + /// evaluating the polynomial + /// \param par coefficients of the polynomial + /// \param x input coordinates + GPUd() float evalPol(GPUgeneric() const float par[/*number of parameters*/], const float x[/*number of dimensions*/]) const { return evalPol(par, x, mDegree, mDim); } + + /// \return returns number of dimensions of the polynomials + GPUd() unsigned int getDim() const { return mDim; } + + /// \return returns the degree of the polynomials + GPUd() unsigned int getDegree() const { return mDegree; } + + protected: + unsigned int mDim{}; ///< dimensionality of the polynomial + unsigned int mDegree{}; ///< maximum degree of the polynomial + + private: + static constexpr unsigned short FMaxdegree = 9; ///< maximum degree of the polynomials (can be increased if desired: size of array in combination_with_repetiton: pos[FMaxdegree + 1]) + + /// evalutes the polynomial + GPUd() static constexpr float evalPol(GPUgeneric() const float par[], const float x[], const unsigned int degree, const unsigned int dim); + + /// helper function to get all combinations + template + GPUd() static Type combination_with_repetiton(const unsigned int degree, const unsigned int dim, GPUgeneric() const float par[], int& indexPar, const float x[]); +}; + +//================================================================================= +//============================ inline implementations ============================= +//================================================================================= + +template +GPUd() constexpr unsigned int MultivariatePolynomialHelper::resetIndices(const unsigned int degreePol, const unsigned int pos, const unsigned int leftDigit, const unsigned int iter, const unsigned int refDigit) +{ + if (iter <= degreePol) { + const int powTmp = pow10(leftDigit); + const int rightDigit = mod10(pos, powTmp); + const int posTmp = pos - (rightDigit - refDigit) * powTmp; + return resetIndices(degreePol, posTmp, leftDigit - 1, iter + 1, refDigit); + } + return pos; +} + +template +GPUd() constexpr unsigned int MultivariatePolynomialHelper::getNewPos(const unsigned int degreePol, const unsigned int pos, const unsigned int digitPos) +{ + if (degreePol > digitPos) { + // check if digit of current position is at is max position + if (mod10(pos, pow10(digitPos)) == Dim) { + // increase digit of left position + const unsigned int leftDigit = digitPos + 1; + const unsigned int posTmp = pos + pow10(leftDigit); + const unsigned int refDigit = mod10(posTmp, pow10(digitPos + 1)); + + // resetting digits to the right if digit exceeds number of dimensions + const unsigned int posReset = resetIndices(degreePol, posTmp, leftDigit - 1, degreePol - digitPos, refDigit); + + // check next digit + return getNewPos(degreePol, posReset, digitPos + 1); + } + return getNewPos(degreePol, pos, digitPos + 1); + } + return pos; +} + +template +template +GPUd() constexpr float MultivariatePolynomialHelper::prodTerm(const float x[], const unsigned int pos) +{ + if constexpr (DegreePol > 0) { + // extract index of the dimension which is decoded in the digit + const unsigned int index = mod10(pos, pow10(DegreePol - 1)); + return x[index] * prodTerm(x, pos); + } + return 1; +} + +template +template +GPUd() constexpr float MultivariatePolynomialHelper::sumTerms(GPUgeneric() const float par[], const float x[]) +{ + // checking if the current position is reasonable e.g. if the max dimension is x[4]: for Pos=15 -> x[1]*x[5] the position is set to 22 -> x[2]*x[2] + constexpr unsigned int posNew = getNewPos(DegreePol, Pos, 0); + if constexpr (mod10(posNew, pow10(DegreePol)) != 1) { + // sum up the term for corrent term and set posotion for next combination + return par[Index] * prodTerm(x, posNew) + sumTerms(par, x); + } + return 0; +} + +template +template +GPUd() constexpr float MultivariatePolynomialHelper::loopDegrees(GPUgeneric() const float par[], const float x[]) +{ + if constexpr (DegreePol <= Degree) { + constexpr unsigned int index{getNParameters(DegreePol - 1, Dim)}; // offset of the index for accessing the parameters + return sumTerms(par, x) + loopDegrees(par, x); + } + return 0; +} + +template +GPUd() Type MultivariatePolynomialHelper<0, 0>::combination_with_repetiton(const unsigned int degree, const unsigned int dim, GPUgeneric() const float par[], int& indexPar, const float x[]) +{ + { + const unsigned int size = degree + 1; + unsigned int pos[FMaxdegree + 1]{0, 0, 0, 0, 0, 0, 0, 0, 0, 0}; + + // return value is either the sum of all polynomials or a vector of strings containing the formula for each polynomial + Type val(0); + for (;;) { + // starting on the rightmost digit + for (unsigned int i = degree; i > 0; --i) { + // check if digit of current position is at is max position + if (pos[i] == dim) { + // increase digit of left position + ++pos[i - 1]; + // resetting the indices of the digits to the right + for (unsigned int j = i; j <= degree; ++j) { + pos[j] = pos[i - 1]; + } + } + } + + // check if all combinations are processed + if (pos[0] == 1) { + break; + } else { + if constexpr (std::is_same_v) { + float term = par[indexPar++]; + for (size_t i = 1; i < size; ++i) { + term *= x[pos[i]]; + } + val += term; + } else { +#if !defined(GPUCA_GPUCODE) + std::string term{}; + for (size_t i = 1; i < size; ++i) { + term += fmt::format("x[{}] * ", pos[i]); + } + term += fmt::format("par[{}]", indexPar++); + val.emplace_back(term); +#endif + } + } + // increase the rightmost digit + ++pos[degree]; + } + return val; + } +} + +GPUd() constexpr float MultivariatePolynomialHelper<0, 0>::evalPol(GPUgeneric() const float par[], const float x[], const unsigned int degree, const unsigned int dim) +{ + float val = par[0]; + int indexPar = 1; + for (unsigned int deg = 1; deg <= degree; ++deg) { + val += combination_with_repetiton(deg, dim, par, indexPar, x); + } + return val; +} + +} // namespace GPUCA_NAMESPACE::gpu + +#endif diff --git a/GPU/TPCFastTransformation/TPCFastTransformationLinkDef_O2.h b/GPU/TPCFastTransformation/TPCFastTransformationLinkDef_O2.h index c072a5201ad56..211d749d94049 100644 --- a/GPU/TPCFastTransformation/TPCFastTransformationLinkDef_O2.h +++ b/GPU/TPCFastTransformation/TPCFastTransformationLinkDef_O2.h @@ -59,5 +59,6 @@ #pragma link C++ class o2::gpu::TPCFastTransform + ; #pragma link C++ class o2::gpu::TPCFastSpaceChargeCorrection::SliceInfo + ; #pragma link C++ class o2::gpu::TPCFastSpaceChargeCorrection + ; +#pragma link C++ struct o2::gpu::MultivariatePolynomialContainer +; #endif diff --git a/GPU/TPCFastTransformation/test/testMultivarPolynomials.cxx b/GPU/TPCFastTransformation/test/testMultivarPolynomials.cxx new file mode 100644 index 0000000000000..bce935178aa72 --- /dev/null +++ b/GPU/TPCFastTransformation/test/testMultivarPolynomials.cxx @@ -0,0 +1,85 @@ +// 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 testMultivarPolynomials.cxx +/// \author Matthias Kleiner + +#define BOOST_TEST_MODULE Test TPC Fast Transformation +#define BOOST_TEST_MAIN +#define BOOST_TEST_DYN_LINK + +#include +#include "MultivariatePolynomial.h" + +namespace o2::gpu +{ + +// evaluate the polynomial of 4th degree in 5 dimensions for given coordinates and parameters +float evalPol4_5D(const float* x, const float* par) +{ + return par[0] + par[1] * x[0] + par[2] * x[1] + par[3] * x[2] + par[4] * x[3] + par[5] * x[4] + par[6] * x[0] * x[0] + par[7] * x[0] * x[1] + par[8] * x[0] * x[2] + par[9] * x[0] * x[3] + par[10] * x[0] * x[4] + par[11] * x[1] * x[1] + par[12] * x[1] * x[2] + par[13] * x[1] * x[3] + par[14] * x[1] * x[4] + par[15] * x[2] * x[2] + par[16] * x[2] * x[3] + par[17] * x[2] * x[4] + par[18] * x[3] * x[3] + par[19] * x[3] * x[4] + par[20] * x[4] * x[4] + par[21] * x[0] * x[0] * x[0] + par[22] * x[0] * x[0] * x[1] + par[23] * x[0] * x[0] * x[2] + par[24] * x[0] * x[0] * x[3] + par[25] * x[0] * x[0] * x[4] + par[26] * x[0] * x[1] * x[1] + par[27] * x[0] * x[1] * x[2] + par[28] * x[0] * x[1] * x[3] + par[29] * x[0] * x[1] * x[4] + par[30] * x[0] * x[2] * x[2] + par[31] * x[0] * x[2] * x[3] + par[32] * x[0] * x[2] * x[4] + par[33] * x[0] * x[3] * x[3] + par[34] * x[0] * x[3] * x[4] + par[35] * x[0] * x[4] * x[4] + par[36] * x[1] * x[1] * x[1] + par[37] * x[1] * x[1] * x[2] + par[38] * x[1] * x[1] * x[3] + par[39] * x[1] * x[1] * x[4] + par[40] * x[1] * x[2] * x[2] + par[41] * x[1] * x[2] * x[3] + par[42] * x[1] * x[2] * x[4] + par[43] * x[1] * x[3] * x[3] + par[44] * x[1] * x[3] * x[4] + par[45] * x[1] * x[4] * x[4] + par[46] * x[2] * x[2] * x[2] + par[47] * x[2] * x[2] * x[3] + par[48] * x[2] * x[2] * x[4] + par[49] * x[2] * x[3] * x[3] + par[50] * x[2] * x[3] * x[4] + par[51] * x[2] * x[4] * x[4] + par[52] * x[3] * x[3] * x[3] + par[53] * x[3] * x[3] * x[4] + par[54] * x[3] * x[4] * x[4] + par[55] * x[4] * x[4] * x[4] + par[56] * x[0] * x[0] * x[0] * x[0] + par[57] * x[0] * x[0] * x[0] * x[1] + par[58] * x[0] * x[0] * x[0] * x[2] + par[59] * x[0] * x[0] * x[0] * x[3] + par[60] * x[0] * x[0] * x[0] * x[4] + par[61] * x[0] * x[0] * x[1] * x[1] + par[62] * x[0] * x[0] * x[1] * x[2] + par[63] * x[0] * x[0] * x[1] * x[3] + par[64] * x[0] * x[0] * x[1] * x[4] + par[65] * x[0] * x[0] * x[2] * x[2] + par[66] * x[0] * x[0] * x[2] * x[3] + par[67] * x[0] * x[0] * x[2] * x[4] + par[68] * x[0] * x[0] * x[3] * x[3] + par[69] * x[0] * x[0] * x[3] * x[4] + par[70] * x[0] * x[0] * x[4] * x[4] + par[71] * x[0] * x[1] * x[1] * x[1] + par[72] * x[0] * x[1] * x[1] * x[2] + par[73] * x[0] * x[1] * x[1] * x[3] + par[74] * x[0] * x[1] * x[1] * x[4] + par[75] * x[0] * x[1] * x[2] * x[2] + par[76] * x[0] * x[1] * x[2] * x[3] + par[77] * x[0] * x[1] * x[2] * x[4] + par[78] * x[0] * x[1] * x[3] * x[3] + par[79] * x[0] * x[1] * x[3] * x[4] + par[80] * x[0] * x[1] * x[4] * x[4] + par[81] * x[0] * x[2] * x[2] * x[2] + par[82] * x[0] * x[2] * x[2] * x[3] + par[83] * x[0] * x[2] * x[2] * x[4] + par[84] * x[0] * x[2] * x[3] * x[3] + par[85] * x[0] * x[2] * x[3] * x[4] + par[86] * x[0] * x[2] * x[4] * x[4] + par[87] * x[0] * x[3] * x[3] * x[3] + par[88] * x[0] * x[3] * x[3] * x[4] + par[89] * x[0] * x[3] * x[4] * x[4] + par[90] * x[0] * x[4] * x[4] * x[4] + par[91] * x[1] * x[1] * x[1] * x[1] + par[92] * x[1] * x[1] * x[1] * x[2] + par[93] * x[1] * x[1] * x[1] * x[3] + par[94] * x[1] * x[1] * x[1] * x[4] + par[95] * x[1] * x[1] * x[2] * x[2] + par[96] * x[1] * x[1] * x[2] * x[3] + par[97] * x[1] * x[1] * x[2] * x[4] + par[98] * x[1] * x[1] * x[3] * x[3] + par[99] * x[1] * x[1] * x[3] * x[4] + par[100] * x[1] * x[1] * x[4] * x[4] + par[101] * x[1] * x[2] * x[2] * x[2] + par[102] * x[1] * x[2] * x[2] * x[3] + par[103] * x[1] * x[2] * x[2] * x[4] + par[104] * x[1] * x[2] * x[3] * x[3] + par[105] * x[1] * x[2] * x[3] * x[4] + par[106] * x[1] * x[2] * x[4] * x[4] + par[107] * x[1] * x[3] * x[3] * x[3] + par[108] * x[1] * x[3] * x[3] * x[4] + par[109] * x[1] * x[3] * x[4] * x[4] + par[110] * x[1] * x[4] * x[4] * x[4] + par[111] * x[2] * x[2] * x[2] * x[2] + par[112] * x[2] * x[2] * x[2] * x[3] + par[113] * x[2] * x[2] * x[2] * x[4] + par[114] * x[2] * x[2] * x[3] * x[3] + par[115] * x[2] * x[2] * x[3] * x[4] + par[116] * x[2] * x[2] * x[4] * x[4] + par[117] * x[2] * x[3] * x[3] * x[3] + par[118] * x[2] * x[3] * x[3] * x[4] + par[119] * x[2] * x[3] * x[4] * x[4] + par[120] * x[2] * x[4] * x[4] * x[4] + par[121] * x[3] * x[3] * x[3] * x[3] + par[122] * x[3] * x[3] * x[3] * x[4] + par[123] * x[3] * x[3] * x[4] * x[4] + par[124] * x[3] * x[4] * x[4] * x[4] + par[125] * x[4] * x[4] * x[4] * x[4]; +} + +float genRand() +{ + const float minVal = -5; + const float maxVal = 5; + const float val = minVal + static_cast(rand()) / (static_cast(RAND_MAX / (maxVal - minVal))); + return val; +} + +BOOST_AUTO_TEST_CASE(Polynomials5D) +{ + std::srand(std::time(nullptr)); + const int nPar5D4Deg = 126; // number of parameters + const int nDim = 5; // dimensions + const int nDegree = 4; // degree + const float abstolerance = 0.0001f; // abosulte difference between refernce to polynomial class + + MultivariatePolynomial polCT; // compile time polynomial + MultivariatePolynomial<0, 0> polRT(nDim, nDegree); // run time polynomial + + // compare number of parameters + BOOST_CHECK(nPar5D4Deg == polCT.getNParams()); + BOOST_CHECK(nPar5D4Deg == polRT.getNParams()); + + float par[nPar5D4Deg]{20}; + for (int iter = 0; iter < 10; ++iter) { + + // draw random parameters + for (int i = 1; i < nPar5D4Deg; ++i) { + par[i] = genRand(); + } + + polCT.setParams(par); + polRT.setParams(par); + + // compare evaluated polynomials + for (float a = 0; a < 1; a += 0.2f) { + for (float b = 0; b < 1; b += 0.2f) { + for (float c = 0; c < 1; c += 0.2f) { + for (float d = 0; d < 1; d += 0.2f) { + for (float e = 0; e < 1; e += 0.2f) { + const float arr[nDim]{a, b, c, d, e}; + const float valCT = polCT.eval(arr); + const float valRT = polRT.eval(arr); + const float valRef = evalPol4_5D(arr, par); + BOOST_CHECK_SMALL(valCT - valRef, abstolerance); + BOOST_CHECK_SMALL(valRT - valRef, abstolerance); + } + } + } + } + } + } +} + +} // namespace o2::gpu diff --git a/GPU/Workflow/src/GPUWorkflowSpec.cxx b/GPU/Workflow/src/GPUWorkflowSpec.cxx index dffc82d77d46b..ec943e528d776 100644 --- a/GPU/Workflow/src/GPUWorkflowSpec.cxx +++ b/GPU/Workflow/src/GPUWorkflowSpec.cxx @@ -36,7 +36,6 @@ #include "TPCReconstruction/TPCFastTransformHelperO2.h" #include "DataFormatsTPC/Digit.h" #include "TPCFastTransform.h" -#include "DataFormatsTPC/CalibdEdxContainer.h" #include "DPLUtils/DPLRawParser.h" #include "DPLUtils/DPLRawPageSequencer.h" #include "DetectorsBase/MatLayerCylSet.h" @@ -48,6 +47,7 @@ #include "GPUO2InterfaceConfiguration.h" #include "GPUO2InterfaceQA.h" #include "GPUO2Interface.h" +#include "CalibdEdxContainer.h" #include "TPCPadGainCalib.h" #include "GPUDisplayBackend.h" #ifdef GPUCA_BUILD_EVENT_DISPLAY @@ -265,6 +265,18 @@ DataProcessorSpec getGPURecoWorkflowSpec(gpuworkflow::CompletionPolicyData* poli LOGP(info, "Loading dEdx correction from file: {}", confParam.dEdxCorrFile); processAttributes->dEdxCalibContainer->loadResidualCorrectionFromFile(confParam.dEdxCorrFile); } + if (std::filesystem::exists(confParam.thresholdCalibFile)) { + LOG(info) << "Loading tpc zero supression map from file " << confParam.thresholdCalibFile; + const auto* thresholdMap = o2::tpc::utils::readCalPads(confParam.thresholdCalibFile, "ThresholdMap")[0]; + processAttributes->dEdxCalibContainer->setZeroSupresssionThreshold(*thresholdMap); + } else { + if (not confParam.thresholdCalibFile.empty()) { + LOG(warn) << "Couldn't find tpc zero supression file " << confParam.thresholdCalibFile << ". Not setting any zero supression."; + } + LOG(info) << "Setting default zero supression map"; + processAttributes->dEdxCalibContainer->setDefaultZeroSupresssionThreshold(); + } + } else { processAttributes->dEdxCalibContainer.reset(new o2::tpc::CalibdEdxContainer()); }