TPCFastSpaceChargeCorrection.h

Go to the documentation of this file.

// 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.
 
 
#ifndef ALICEO2_GPUCOMMON_TPCFASTTRANSFORMATION_TPCFASTSPACECHARGECORRECTION_H
#define ALICEO2_GPUCOMMON_TPCFASTTRANSFORMATION_TPCFASTSPACECHARGECORRECTION_H
 
#include "Spline2D.h"
#include "TPCFastTransformGeo.h"
#include "FlatObject.h"
#include "GPUCommonDef.h"
#include "GPUCommonMath.h"
 
namespace o2
{
namespace gpu
{
 
class TPCFastSpaceChargeCorrection : public FlatObject
{
 public:
  struct RowInfo {
    int32_t splineScenarioID{0}; 
    size_t dataOffsetBytes[3]{0}; 
    ClassDefNV(RowInfo, 1);
  };
 
  struct RowActiveArea {
    float maxDriftLengthCheb[5]{0.f};
    float vMax{0.f};
    float cuMin{0.f};
    float cuMax{0.f};
    float cvMax{0.f};
    ClassDefNV(RowActiveArea, 1);
  };
 
  struct SliceRowInfo {
    float gridV0{0.f};           
    float gridCorrU0{0.f};       
    float gridCorrV0{0.f};       
    float scaleCorrUtoGrid{0.f}; 
    float scaleCorrVtoGrid{0.f}; 
    RowActiveArea activeArea;
    ClassDefNV(SliceRowInfo, 1);
  };
 
  struct SliceInfo {
    float vMax{0.f}; 
    ClassDefNV(SliceInfo, 1);
  };
 
  typedef Spline2D<float, 3> SplineType;
 
 
  TPCFastSpaceChargeCorrection();
 
  TPCFastSpaceChargeCorrection(const TPCFastSpaceChargeCorrection&) = delete;
 
  TPCFastSpaceChargeCorrection& operator=(const TPCFastSpaceChargeCorrection&) = delete;
 
  ~TPCFastSpaceChargeCorrection();
 
 
 
  using FlatObject::getBufferAlignmentBytes;
  using FlatObject::getClassAlignmentBytes;
 
 
  void cloneFromObject(const TPCFastSpaceChargeCorrection& obj, char* newFlatBufferPtr);
  void destroy();
 
 
  using FlatObject::releaseInternalBuffer;
  void moveBufferTo(char* newBufferPtr);
 
 
  void setActualBufferAddress(char* actualFlatBufferPtr);
  void setFutureBufferAddress(char* futureFlatBufferPtr);
 
 
  void startConstruction(const TPCFastTransformGeo& geo, int32_t numberOfSplineScenarios);
 
  void setRowScenarioID(int32_t iRow, int32_t iScenario);
 
  void setSplineScenario(int32_t scenarioIndex, const SplineType& spline);
 
  void finishConstruction();
 
  void constructWithNoCorrection(const TPCFastTransformGeo& geo);
 
 
  GPUd() void setNoCorrection();
 
  GPUd() void setTimeStamp(int64_t v) { mTimeStamp = v; }
 
  GPUd() void setInterpolationSafetyMargin(float val) { fInterpolationSafetyMargin = val; }
 
  GPUd() const SplineType& getSpline(int32_t slice, int32_t row) const;
 
  GPUd() SplineType& getSpline(int32_t slice, int32_t row);
 
  GPUd() float* getSplineData(int32_t slice, int32_t row, int32_t iSpline = 0);
 
  GPUd() const float* getSplineData(int32_t slice, int32_t row, int32_t iSpline = 0) const;
 
  GPUd() int32_t getCorrection(int32_t slice, int32_t row, float u, float v, float& dx, float& du, float& dv) const;
 
  GPUd() void getCorrectionInvCorrectedX(int32_t slice, int32_t row, float corrU, float corrV, float& corrX) const;
 
  GPUd() void getCorrectionInvUV(int32_t slice, int32_t row, float corrU, float corrV, float& nomU, float& nomV) const;
 
  GPUd() float getMaxDriftLength(int32_t slice, int32_t row, float pad) const;
 
  GPUd() float getMaxDriftLength(int32_t slice, int32_t row) const;
 
  GPUd() float getMaxDriftLength(int32_t slice) const;
 
 
  GPUd() void schrinkUV(int32_t slice, int32_t row, float& u, float& v) const;
 
  GPUd() void schrinkCorrectedUV(int32_t slice, int32_t row, float& corrU, float& corrV) const;
 
  GPUd() void convUVtoGrid(int32_t slice, int32_t row, float u, float v, float& gridU, float& gridV) const;
 
  GPUd() void convGridToUV(int32_t slice, int32_t row, float gridU, float gridV, float& u, float& v) const;
 
  GPUd() void convCorrectedUVtoGrid(int32_t slice, int32_t row, float cu, float cv, float& gridU, float& gridV) const;
 
  GPUd() const TPCFastTransformGeo& getGeometry() const
  {
    return mGeo;
  }
 
  int64_t getTimeStamp() const { return mTimeStamp; }
 
  GPUd() float getInterpolationSafetyMargin() const { return fInterpolationSafetyMargin; }
 
  GPUd() const RowInfo& getRowInfo(int32_t row) const { return mRowInfoPtr[row]; }
 
  GPUd() const SliceInfo& getSliceInfo(int32_t slice) const
  {
    return mSliceInfo[slice];
  }
 
  GPUd() SliceInfo& getSliceInfo(int32_t slice)
  {
    return mSliceInfo[slice];
  }
 
  GPUd() const SliceRowInfo& getSliceRowInfo(int32_t slice, int32_t row) const
  {
    return mSliceRowInfoPtr[mGeo.getNumberOfRows() * slice + row];
  }
 
  GPUd() SliceRowInfo& getSliceRowInfo(int32_t slice, int32_t row)
  {
    return mSliceRowInfoPtr[mGeo.getNumberOfRows() * slice + row];
  }
 
#if !defined(GPUCA_GPUCODE)
  void print() const;
  GPUh() double testInverse(bool prn = 0);
#endif
 
 private:
  void relocateBufferPointers(const char* oldBuffer, char* newBuffer);
  void releaseConstructionMemory();
 
  GPUd() int32_t getCorrectionOld(int32_t slice, int32_t row, float u, float v, float& dx, float& du, float& dv) const;
 
 
 
  RowInfo* mConstructionRowInfos = nullptr;     
  SplineType* mConstructionScenarios = nullptr; 
 
 
  TPCFastTransformGeo mGeo; 
 
  int32_t mNumberOfScenarios; 
 
  SliceInfo mSliceInfo[TPCFastTransformGeo::getNumberOfSlices()]; 
 
  SplineType* mScenarioPtr;       
  RowInfo* mRowInfoPtr;           
  SliceRowInfo* mSliceRowInfoPtr; 
 
 
  int64_t mTimeStamp; 
 
  char* mSplineData[3]; 
 
  size_t mSliceDataSizeBytes[3]; 
 
  float fInterpolationSafetyMargin{0.1f}; // 10% area around the TPC row. Outside of this area the interpolation returns the boundary values.
 
  ClassDefNV(TPCFastSpaceChargeCorrection, 3);
};
 
 
GPUdi() const TPCFastSpaceChargeCorrection::SplineType& TPCFastSpaceChargeCorrection::getSpline(int32_t slice, int32_t row) const
{
  const RowInfo& rowInfo = mRowInfoPtr[row];
  return mScenarioPtr[rowInfo.splineScenarioID];
}
 
GPUdi() TPCFastSpaceChargeCorrection::SplineType& TPCFastSpaceChargeCorrection::getSpline(int32_t slice, int32_t row)
{
  const RowInfo& rowInfo = mRowInfoPtr[row];
  return mScenarioPtr[rowInfo.splineScenarioID];
}
 
GPUdi() float* TPCFastSpaceChargeCorrection::getSplineData(int32_t slice, int32_t row, int32_t iSpline)
{
  const RowInfo& rowInfo = mRowInfoPtr[row];
  return reinterpret_cast<float*>(mSplineData[iSpline] + mSliceDataSizeBytes[iSpline] * slice + rowInfo.dataOffsetBytes[iSpline]);
}
 
GPUdi() const float* TPCFastSpaceChargeCorrection::getSplineData(int32_t slice, int32_t row, int32_t iSpline) const
{
  const RowInfo& rowInfo = mRowInfoPtr[row];
  return reinterpret_cast<float*>(mSplineData[iSpline] + mSliceDataSizeBytes[iSpline] * slice + rowInfo.dataOffsetBytes[iSpline]);
}
 
GPUdi() void TPCFastSpaceChargeCorrection::schrinkUV(int32_t slice, int32_t row, float& u, float& v) const
{
 
  const TPCFastTransformGeo::RowInfo& rowInfo = mGeo.getRowInfo(row);
 
  float uWidth05 = rowInfo.getUwidth() * (0.5f + fInterpolationSafetyMargin);
  float vWidth = mGeo.getTPCzLength(slice);
 
  if (u < -uWidth05) {
    u = -uWidth05;
  }
  if (u > uWidth05) {
    u = uWidth05;
  }
  if (v < -0.1f * vWidth) {
    v = -0.1f * vWidth;
  }
  if (v > 1.1f * vWidth) {
    v = 1.1f * vWidth;
  }
}
 
GPUdi() void TPCFastSpaceChargeCorrection::schrinkCorrectedUV(int32_t slice, int32_t row, float& corrU, float& corrV) const
{
 
  const TPCFastTransformGeo::RowInfo& rowInfo = mGeo.getRowInfo(row);
  const SliceRowInfo& sliceRowInfo = getSliceRowInfo(slice, row);
 
  float uMargin = fInterpolationSafetyMargin * rowInfo.getUwidth();
  float vMargin = fInterpolationSafetyMargin * mGeo.getTPCzLength(slice);
 
  if (corrU < sliceRowInfo.activeArea.cuMin - uMargin) {
    corrU = sliceRowInfo.activeArea.cuMin - uMargin;
  }
 
  if (corrU > sliceRowInfo.activeArea.cuMax + uMargin) {
    corrU = sliceRowInfo.activeArea.cuMax + uMargin;
  }
 
  if (corrV < 0.f - vMargin) {
    corrV = 0.f - vMargin;
  }
 
  if (corrV > sliceRowInfo.activeArea.cvMax + vMargin) {
    corrV = sliceRowInfo.activeArea.cvMax + vMargin;
  }
}
 
GPUdi() void TPCFastSpaceChargeCorrection::convUVtoGrid(int32_t slice, int32_t row, float u, float v, float& gu, float& gv) const
{
  // TODO optimise !!!
  gu = 0.f;
  gv = 0.f;
 
  schrinkUV(slice, row, u, v);
 
  const SliceRowInfo& info = getSliceRowInfo(slice, row);
  const SplineType& spline = getSpline(slice, row);
 
  float su0 = 0.f, sv0 = 0.f;
  mGeo.convUVtoScaledUV(slice, row, u, info.gridV0, su0, sv0);
  mGeo.convUVtoScaledUV(slice, row, u, v, gu, gv);
 
  gv = (gv - sv0) / (1.f - sv0);
  gu *= spline.getGridX1().getUmax();
  gv *= spline.getGridX2().getUmax();
}
 
GPUdi() void TPCFastSpaceChargeCorrection::convGridToUV(int32_t slice, int32_t row, float gridU, float gridV, float& u, float& v) const
{
  // TODO optimise
  float su0 = 0.f, sv0 = 0.f;
  const SliceRowInfo& info = getSliceRowInfo(slice, row);
  const SplineType& spline = getSpline(slice, row);
  mGeo.convUVtoScaledUV(slice, row, 0.f, info.gridV0, su0, sv0);
  float su = gridU / spline.getGridX1().getUmax();
  float sv = sv0 + gridV / spline.getGridX2().getUmax() * (1.f - sv0);
  mGeo.convScaledUVtoUV(slice, row, su, sv, u, v);
}
 
GPUdi() void TPCFastSpaceChargeCorrection::convCorrectedUVtoGrid(int32_t slice, int32_t row, float corrU, float corrV, float& gridU, float& gridV) const
{
  schrinkCorrectedUV(slice, row, corrU, corrV);
 
  const SliceRowInfo& sliceRowInfo = getSliceRowInfo(slice, row);
 
  gridU = (corrU - sliceRowInfo.gridCorrU0) * sliceRowInfo.scaleCorrUtoGrid;
  gridV = (corrV - sliceRowInfo.gridCorrV0) * sliceRowInfo.scaleCorrVtoGrid;
}
 
GPUdi() int32_t TPCFastSpaceChargeCorrection::getCorrection(int32_t slice, int32_t row, float u, float v, float& dx, float& du, float& dv) const
{
  const SplineType& spline = getSpline(slice, row);
  const float* splineData = getSplineData(slice, row);
  float gridU = 0, gridV = 0;
  convUVtoGrid(slice, row, u, v, gridU, gridV);
  float dxuv[3];
  spline.interpolateU(splineData, gridU, gridV, dxuv);
  if (CAMath::Abs(dxuv[0]) > 100 || CAMath::Abs(dxuv[1]) > 100 || CAMath::Abs(dxuv[2]) > 100) {
    dxuv[0] = dxuv[1] = dxuv[2] = 0;
  }
  dx = dxuv[0];
  du = dxuv[1];
  dv = dxuv[2];
  return 0;
}
 
GPUdi() int32_t TPCFastSpaceChargeCorrection::getCorrectionOld(int32_t slice, int32_t row, float u, float v, float& dx, float& du, float& dv) const
{
  const SplineType& spline = getSpline(slice, row);
  const float* splineData = getSplineData(slice, row);
  float gridU = 0, gridV = 0;
  convUVtoGrid(slice, row, u, v, gridU, gridV);
  float dxuv[3];
  spline.interpolateUold(splineData, gridU, gridV, dxuv);
  if (CAMath::Abs(dxuv[0]) > 100 || CAMath::Abs(dxuv[1]) > 100 || CAMath::Abs(dxuv[2]) > 100) {
    dxuv[0] = dxuv[1] = dxuv[2] = 0;
  }
  dx = dxuv[0];
  du = dxuv[1];
  dv = dxuv[2];
  return 0;
}
 
GPUdi() void TPCFastSpaceChargeCorrection::getCorrectionInvCorrectedX(
  int32_t slice, int32_t row, float corrU, float corrV, float& x) const
{
  float gridU, gridV;
  convCorrectedUVtoGrid(slice, row, corrU, corrV, gridU, gridV);
 
  const Spline2D<float, 1>& spline = reinterpret_cast<const Spline2D<float, 1>&>(getSpline(slice, row));
  const float* splineData = getSplineData(slice, row, 1);
  float dx = 0;
  spline.interpolateU(splineData, gridU, gridV, &dx);
  if (CAMath::Abs(dx) > 100) {
    dx = 0;
  }
  x = mGeo.getRowInfo(row).x + dx;
}
 
GPUdi() void TPCFastSpaceChargeCorrection::getCorrectionInvUV(
  int32_t slice, int32_t row, float corrU, float corrV, float& nomU, float& nomV) const
{
  float gridU, gridV;
  convCorrectedUVtoGrid(slice, row, corrU, corrV, gridU, gridV);
 
  const Spline2D<float, 2>& spline = reinterpret_cast<const Spline2D<float, 2>&>(getSpline(slice, row));
  const float* splineData = getSplineData(slice, row, 2);
 
  float duv[2];
  spline.interpolateU(splineData, gridU, gridV, duv);
  if (CAMath::Abs(duv[0]) > 100 || CAMath::Abs(duv[1]) > 100) {
    duv[0] = duv[1] = 0;
  }
  nomU = corrU - duv[0];
  nomV = corrV - duv[1];
}
 
GPUdi() float TPCFastSpaceChargeCorrection::getMaxDriftLength(int32_t slice, int32_t row, float pad) const
{
  const RowActiveArea& area = getSliceRowInfo(slice, row).activeArea;
  const float* c = area.maxDriftLengthCheb;
  float x = -1.f + 2.f * pad / mGeo.getRowInfo(row).maxPad;
  float y = c[0] + c[1] * x;
  float f0 = 1.f;
  float f1 = x;
  x *= 2.f;
  for (int32_t i = 2; i < 5; i++) {
    double f = x * f1 - f0;
    y += c[i] * f;
    f0 = f1;
    f1 = f;
  }
  return y;
}
 
GPUdi() float TPCFastSpaceChargeCorrection::getMaxDriftLength(int32_t slice, int32_t row) const
{
  return getSliceRowInfo(slice, row).activeArea.vMax;
}
 
GPUdi() float TPCFastSpaceChargeCorrection::getMaxDriftLength(int32_t slice) const
{
  return getSliceInfo(slice).vMax;
}
 
} // namespace gpu
} // namespace o2
 
#endif