POD version of TPCFastSpaceChargeCorrection

GPU/TPCFastTransformation/CMakeLists.txt

@@ -25,6 +25,7 @@ set(SRCS
     TPCFastSpaceChargeCorrection.cxx
     TPCFastSpaceChargeCorrectionMap.cxx
     TPCFastTransform.cxx
+    TPCFastTransformPOD.cxx
     CorrectionMapsHelper.cxx
 )
 

GPU/TPCFastTransformation/TPCFastSpaceChargeCorrection.h

@@ -37,6 +37,8 @@ namespace gpu
 ///
 class TPCFastSpaceChargeCorrection : public FlatObject
 {
+  friend class TPCFastTransformPOD;
+
  public:
   ///
   /// \brief The struct contains necessary info for TPC padrow

GPU/TPCFastTransformation/TPCFastTransform.cxx

@@ -37,7 +37,7 @@
 using namespace GPUCA_NAMESPACE::gpu;
 
 TPCFastTransform::TPCFastTransform()
-  : FlatObject(), mTimeStamp(0), mCorrection(), mApplyCorrection(1), mT0(0.f), mVdrift(0.f), mVdriftCorrY(0.f), mLdriftCorr(0.f), mTOFcorr(0.f), mPrimVtxZ(0.f), mLumi(0.f), mLumiError(0.f), mLumiScaleFactor(1.0f)
+  : FlatObject(), mTimeStamp(0), mCorrection(), mApplyCorrection(1), mT0(0.f), mVdrift(0.f), mVdriftCorrY(0.f), mLdriftCorr(0.f), mTOFcorr(0.f), mPrimVtxZ(0.f), mLumiIDC(0.f), mLumi(0.f), mLumiError(0.f), mLumiScaleFactor(1.0f)
 {
   // Default Constructor: creates an empty uninitialized object
 }
@@ -58,6 +58,7 @@ void TPCFastTransform::cloneFromObject(const TPCFastTransform& obj, char* newFla
   mLdriftCorr = obj.mLdriftCorr;
   mTOFcorr = obj.mTOFcorr;
   mPrimVtxZ = obj.mPrimVtxZ;
+  mLumiIDC = obj.mLumiIDC;
   mLumi = obj.mLumi;
   mLumiError = obj.mLumiError;
   mLumiScaleFactor = obj.mLumiScaleFactor;
@@ -108,6 +109,7 @@ void TPCFastTransform::startConstruction(const TPCFastSpaceChargeCorrection& cor
   mLdriftCorr = 0.f;
   mTOFcorr = 0.f;
   mPrimVtxZ = 0.f;
+  mLumiIDC = 0.f;
   mLumi = 0.f;
   mLumiError = 0.f;
   mLumiScaleFactor = 1.f;
@@ -158,6 +160,7 @@ void TPCFastTransform::print() const
   LOG(info) << "mLdriftCorr = " << mLdriftCorr;
   LOG(info) << "mTOFcorr = " << mTOFcorr;
   LOG(info) << "mPrimVtxZ = " << mPrimVtxZ;
+  LOG(info) << "mLumiIDC = " << mLumiIDC;
   LOG(info) << "mLumi = " << mLumi;
   LOG(info) << "mLumiError = " << mLumiError;
   LOG(info) << "mLumiScaleFactor = " << mLumiScaleFactor;

GPU/TPCFastTransformation/TPCFastTransform.h

@@ -162,6 +162,7 @@ class TPCFastTransform : public FlatObject
 
   /// Set Lumi info
   void setLumi(float l) { mLumi = l; }
+  void setLumiIDC(float l) { mLumiIDC = l; }
   void setLumiError(float e) { mLumiError = e; }
   void setLumiScaleFactor(float s) { mLumiScaleFactor = s; }
 
@@ -250,9 +251,15 @@ class TPCFastTransform : public FlatObject
   /// Return TOF correction (vdrift / C)
   GPUd() float getTOFCorr() const { return mLdriftCorr; }
 
+  /// Return nominal PV Z position correction (vdrift / C)
+  GPUd() float getPrimVtxZ() const { return mPrimVtxZ; }
+
   /// Return map lumi
   GPUd() float getLumi() const { return mLumi; }
 
+  /// Return map lumiIDC
+  GPUd() float getLumiIDC() const { return mLumiIDC; }
+
   /// Return map lumi error
   GPUd() float getLumiError() const { return mLumiError; }
 
@@ -332,6 +339,7 @@ class TPCFastTransform : public FlatObject
 
   float mPrimVtxZ; ///< Z of the primary vertex, needed for the Time-Of-Flight correction
 
+  float mLumiIDC;         ///< luminosity estimator in IDC units
   float mLumi;            ///< luminosity estimator
   float mLumiError;       ///< error on luminosity
   float mLumiScaleFactor; ///< user correction factor for lumi (e.g. normalization, efficiency correction etc.)
@@ -342,7 +350,7 @@ class TPCFastTransform : public FlatObject
   GPUd() void TransformInternal(int slice, int row, float& u, float& v, float& x, const TPCFastTransform* ref, const TPCFastTransform* ref2, float scale, float scale2, int scaleMode) const;
 
 #ifndef GPUCA_ALIROOT_LIB
-  ClassDefNV(TPCFastTransform, 3);
+  ClassDefNV(TPCFastTransform, 4);
 #endif
 };
 

GPU/TPCFastTransformation/TPCFastTransformPOD.cxx

@@ -0,0 +1,243 @@
+// 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  TPCFastTransformPOD.cxx
+/// \brief Implementation of POD correction map
+///
+/// \author  [email protected]
+
+#include "TPCFastTransformPOD.h"
+#include "GPUDebugStreamer.h"
+#if !defined(GPUCA_GPUCODE)
+#include <TRandom.h>
+#endif
+
+namespace GPUCA_NAMESPACE
+{
+namespace gpu
+{
+
+#if !defined(GPUCA_GPUCODE)
+
+size_t TPCFastTransformPOD::estimateSize(const TPCFastSpaceChargeCorrection& origCorr)
+{
+  // estimate size of own buffer
+  const size_t selfSizeFix = sizeof(TPCFastTransformPOD);
+  size_t nextDynOffs = alignOffset(selfSizeFix);
+  nextDynOffs = alignOffset(nextDynOffs + origCorr.mNumberOfScenarios * sizeof(size_t)); // spline scenarios start here
+  // space for splines
+  for (int isc = 0; isc < origCorr.mNumberOfScenarios; isc++) {
+    const auto& spline = origCorr.mScenarioPtr[isc];
+    nextDynOffs = alignOffset(nextDynOffs + sizeof(spline));
+  }
+  // space for splines data
+  for (int is = 0; is < 3; is++) {
+    for (int slice = 0; slice < origCorr.mGeo.getNumberOfSlices(); slice++) {
+      for (int row = 0; row < NROWS; row++) {
+        const auto& spline = origCorr.getSpline(slice, row);
+        int nPar = spline.getNumberOfParameters();
+        if (is == 1) {
+          nPar = nPar / 3;
+        }
+        if (is == 2) {
+          nPar = nPar * 2 / 3;
+        }
+        nextDynOffs += nPar * sizeof(float);
+      }
+    }
+  }
+  nextDynOffs = alignOffset(nextDynOffs);
+  return nextDynOffs;
+}
+
+TPCFastTransformPOD& TPCFastTransformPOD::create(char* buff, size_t buffSize, const TPCFastSpaceChargeCorrection& origCorr)
+{
+  // instantiate object to already created buffer of the right size
+  assert(buffSize > sizeof(TPCFastTransformPOD));
+  auto& podMap = getNonConst(buff);
+  podMap.mApplyCorrections = true; // by default always apply corrections
+
+  // copy fixed size data --- start
+  podMap.mNumberOfScenarios = origCorr.mNumberOfScenarios;
+  std::memcpy(&podMap.mGeo, &origCorr.mGeo, sizeof(TPCFastTransformGeo)); // copy geometry (fixed size)
+  for (int row = 0; row < NROWS; row++) {
+    podMap.mRowInfo[row] = origCorr.getRowInfo(row); // dataOffsetBytes will be modified later
+  }
+  for (int slice = 0; slice < TPCFastTransformGeo::getNumberOfSlices(); slice++) {
+    podMap.mSliceInfo[slice] = origCorr.getSliceInfo(slice);
+    for (int row = 0; row < NROWS; row++) {
+      podMap.mSliceRowInfo[NROWS * slice + row] = origCorr.getSliceRowInfo(slice, row);
+    }
+  }
+  podMap.mInterpolationSafetyMargin = origCorr.fInterpolationSafetyMargin;
+  podMap.mTimeStamp = origCorr.mTimeStamp;
+  //
+  // init data members coming from the TPCFastTrasform
+  podMap.mVdrift = 0.;
+  podMap.mT0 = 0.;
+  podMap.mVdriftCorrY = 0.;
+  podMap.mLdriftCorr = 0.;
+  podMap.mTOFcorr = 0.;
+  podMap.mPrimVtxZ = 0.;
+  // copy fixed size data --- end
+
+  size_t nextDynOffs = alignOffset(sizeof(TPCFastTransformPOD));
+  // copy slice scenarios
+  podMap.mOffsScenariosOffsets = nextDynOffs; // spline scenarios offsets start here
+  LOGP(debug, "Set mOffsScenariosOffsets = {}", podMap.mOffsScenariosOffsets);
+  nextDynOffs = alignOffset(nextDynOffs + podMap.mNumberOfScenarios * sizeof(size_t)); // spline scenarios start here
+
+  // copy spline objects
+  size_t* scenOffs = reinterpret_cast<size_t*>(buff + podMap.mOffsScenariosOffsets);
+  for (int isc = 0; isc < origCorr.mNumberOfScenarios; isc++) {
+    scenOffs[isc] = nextDynOffs;
+    const auto& spline = origCorr.mScenarioPtr[isc];
+    if (buffSize < nextDynOffs + sizeof(spline)) {
+      throw std::runtime_error(fmt::format("attempt to copy {} bytes for spline for scenario {} to {}, overflowing the buffer of size {}", sizeof(spline), isc, nextDynOffs + sizeof(spline), buffSize));
+    }
+    std::memcpy(buff + scenOffs[isc], &spline, sizeof(spline));
+    nextDynOffs = alignOffset(nextDynOffs + sizeof(spline));
+    LOGP(debug, "Copy {} bytes for spline scenario {} (ptr:{}) to offsset {}", sizeof(spline), isc, (void*)&spline, scenOffs[isc]);
+  }
+
+  // copy splines data
+  for (int is = 0; is < 3; is++) {
+    float* data = reinterpret_cast<float*>(buff + nextDynOffs);
+    LOGP(debug, "splinID={} start offset {} -> {}", is, nextDynOffs, (void*)data);
+    for (int slice = 0; slice < origCorr.mGeo.getNumberOfSlices(); slice++) {
+      podMap.mSplineDataOffsets[slice][is] = nextDynOffs;
+      size_t rowDataOffs = 0;
+      for (int row = 0; row < NROWS; row++) {
+        const auto& spline = origCorr.getSpline(slice, row);
+        const float* dataOr = origCorr.getSplineData(slice, row, is);
+        int nPar = spline.getNumberOfParameters();
+        if (is == 1) {
+          nPar = nPar / 3;
+        }
+        if (is == 2) {
+          nPar = nPar * 2 / 3;
+        }
+        LOGP(debug, "Copying {} floats for spline{} of slice:{} row:{} to offset {}", nPar, is, slice, row, nextDynOffs);
+        size_t nbcopy = nPar * sizeof(float);
+        if (buffSize < nextDynOffs + nbcopy) {
+          throw std::runtime_error(fmt::format("attempt to copy {} bytes of data for spline{} of slice{}/row{} to {}, overflowing the buffer of size {}", nbcopy, is, slice, row, nextDynOffs, buffSize));
+        }
+        std::memcpy(data, dataOr, nbcopy);
+        podMap.getRowInfo(row).dataOffsetBytes[is] = rowDataOffs;
+        rowDataOffs += nbcopy;
+        data += nPar;
+        nextDynOffs += nbcopy;
+      }
+    }
+  }
+  podMap.mTotalSize = alignOffset(nextDynOffs);
+  if (buffSize != podMap.mTotalSize) {
+    throw std::runtime_error(fmt::format("Estimated buffer size {} differs from filled one {}", buffSize, podMap.mTotalSize));
+  }
+  return podMap;
+}
+
+TPCFastTransformPOD& TPCFastTransformPOD::create(char* buff, size_t buffSize, const TPCFastTransform& src)
+{
+  // instantiate objec to already created buffer of the right size
+  auto& podMap = create(buff, buffSize, src.getCorrection());
+  // set data members of TPCFastTransform
+  podMap.mVdrift = src.getVDrift();
+  podMap.mT0 = src.getT0();
+  podMap.mVdriftCorrY = src.getVdriftCorrY();
+  podMap.mLdriftCorr = src.getLdriftCorr();
+  podMap.mTOFcorr = src.getTOFCorr();
+  podMap.mPrimVtxZ = src.getPrimVtxZ();
+  // copy fixed size data --- end
+  return podMap;
+}
+
+bool TPCFastTransformPOD::test(const TPCFastSpaceChargeCorrection& origCorr, int npoints) const
+{
+  if (npoints < 1) {
+    return false;
+  }
+  std::vector<unsigned char> slice, row;
+  std::vector<float> u, v, dxO, duO, dvO, dxP, duP, dvP, corrXO, corrXP, nomUO, nomVO, nomUP, nomVP;
+  slice.reserve(npoints);
+  row.reserve(npoints);
+  u.reserve(npoints);
+  v.reserve(npoints);
+  dxO.resize(npoints);
+  duO.resize(npoints);
+  dvO.resize(npoints);
+  corrXO.resize(npoints);
+  nomUO.resize(npoints);
+  nomVO.resize(npoints);
+  dxP.resize(npoints);
+  duP.resize(npoints);
+  dvP.resize(npoints);
+  corrXP.resize(npoints);
+  nomUP.resize(npoints);
+  nomVP.resize(npoints);
+
+  for (int i = 0; i < npoints; i++) {
+    slice.push_back(gRandom->Integer(NSLICES));
+    row.push_back(gRandom->Integer(NROWS));
+    u.push_back(gRandom->Rndm() * 15);
+    v.push_back(gRandom->Rndm() * 200);
+  }
+  long origStart[3], origEnd[3], thisStart[3], thisEnd[3];
+  origStart[0] = std::chrono::time_point_cast<std::chrono::microseconds>(std::chrono::system_clock::now()).time_since_epoch().count();
+  for (int i = 0; i < npoints; i++) {
+    origCorr.getCorrection(slice[i], row[i], u[i], v[i], dxO[i], duO[i], dvO[i]);
+  }
+  origEnd[0] = origStart[1] = std::chrono::time_point_cast<std::chrono::microseconds>(std::chrono::system_clock::now()).time_since_epoch().count();
+  for (int i = 0; i < npoints; i++) {
+    origCorr.getCorrectionInvCorrectedX(slice[i], row[i], u[i], v[i], corrXO[i]);
+  }
+  origEnd[1] = origStart[2] = std::chrono::time_point_cast<std::chrono::microseconds>(std::chrono::system_clock::now()).time_since_epoch().count();
+  for (int i = 0; i < npoints; i++) {
+    origCorr.getCorrectionInvUV(slice[i], row[i], u[i], v[i], nomUO[i], nomVO[i]);
+  }
+  origEnd[2] = thisStart[0] = std::chrono::time_point_cast<std::chrono::microseconds>(std::chrono::system_clock::now()).time_since_epoch().count();
+  for (int i = 0; i < npoints; i++) {
+    this->getCorrection(slice[i], row[i], u[i], v[i], dxP[i], duP[i], dvP[i]);
+  }
+  thisEnd[0] = thisStart[1] = std::chrono::time_point_cast<std::chrono::microseconds>(std::chrono::system_clock::now()).time_since_epoch().count();
+  for (int i = 0; i < npoints; i++) {
+    this->getCorrectionInvCorrectedX(slice[i], row[i], u[i], v[i], corrXP[i]);
+  }
+  thisEnd[1] = thisStart[2] = std::chrono::time_point_cast<std::chrono::microseconds>(std::chrono::system_clock::now()).time_since_epoch().count();
+  for (int i = 0; i < npoints; i++) {
+    this->getCorrectionInvUV(slice[i], row[i], u[i], v[i], nomUP[i], nomVP[i]);
+  }
+  thisEnd[2] = std::chrono::time_point_cast<std::chrono::microseconds>(std::chrono::system_clock::now()).time_since_epoch().count();
+  //
+  size_t ndiff[3] = {};
+  for (int i = 0; i < npoints; i++) {
+    if (dxO[i] != dxP[i] || duO[i] != duP[i] || dvO[i] != dvP[i]) {
+      ndiff[0]++;
+    }
+    if (corrXO[i] != corrXP[i]) {
+      ndiff[1]++;
+    }
+    if (nomUO[i] != nomUP[i] || nomVO[i] != nomVP[i]) {
+      ndiff[2]++;
+    }
+  }
+  //
+  LOGP(info, " (ns per call)              original        this     Nmissmatch");
+  LOGP(info, "getCorrection               {:.3e}    {:.3e}   {}", double(origEnd[0] - origStart[0]) / npoints * 1000., double(thisEnd[0] - thisStart[0]) / npoints * 1000., ndiff[0]);
+  LOGP(info, "getCorrectionInvCorrectedX  {:.3e}    {:.3e}   {}", double(origEnd[1] - origStart[1]) / npoints * 1000., double(thisEnd[1] - thisStart[1]) / npoints * 1000., ndiff[1]);
+  LOGP(info, "getCorrectionInvUV          {:.3e}    {:.3e}   {}", double(origEnd[2] - origStart[2]) / npoints * 1000., double(thisEnd[2] - thisStart[2]) / npoints * 1000., ndiff[2]);
+  return ndiff[0] == 0 && ndiff[1] == 0 && ndiff[2] == 0;
+}
+
+#endif
+
+} // namespace gpu
+} // namespace GPUCA_NAMESPACE