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"
#ifndef GPUCA_GPUCODE_DEVICE
#include "GPUCommonArray.h" // Would work on GPU, but yields performance regressions
#endif
 
namespace o2::gpu
{
 
class TPCFastSpaceChargeCorrection : public FlatObject
{
  friend class TPCFastTransformPOD;
 
 public:
  // obsolete structure, declared here only for the backward compatibility
  struct SliceInfo {
    ClassDefNV(SliceInfo, 2);
  };
 
  struct GridInfo {
   private:
    float y0{0.f};                 
    float yScale{0.f};             //< scale Y to U-grid coordinate
    float z0{0.f};                 
    float zScale{0.f};             //< scale Z to V-grid coordinate
    float zOut{0.f};               // outer z of the grid;
    float splineScalingWithZ{0.f}; 
 
   public:
    void set(float y0_, float yScale_, float z0_, float zScale_, float zOut_, float zReadout_)
    {
      this->y0 = y0_;
      this->yScale = yScale_;
      this->z0 = z0_;
      this->zScale = zScale_;
      this->zOut = zOut_;
      // no scaling when the distance to the readout is too small
      this->splineScalingWithZ = fabs(zReadout_ - zOut_) > 1.f ? 1.f / (zReadout_ - zOut_) : 0.f;
    }
 
    float getY0() const { return y0; }
    float getYscale() const { return yScale; }
    float getZ0() const { return z0; }
    float getZscale() const { return zScale; }
 
    GPUdi() float getSpineScaleForZ(float z) const
    {
      return 1.f - GPUCommonMath::Clamp((z - zOut) * splineScalingWithZ, 0.f, 1.f);
    }
 
    GPUdi() void convLocalToGridUntruncated(int sector, float y, float z, float& u, float& v, float& s) const
    {
      if (sector >= TPCFastTransformGeo::getNumberOfSectorsA()) {
        z = -z;
      }
      u = (y - y0) * yScale;
      v = (z - z0) * zScale;
      s = getSpineScaleForZ(z);
    }
 
    GPUdi() void convGridToLocal(int sector, float gridU, float gridV, float& y, float& z) const
    {
      y = y0 + gridU / yScale;
      z = z0 + gridV / zScale;
      if (sector >= TPCFastTransformGeo::getNumberOfSectorsA()) {
        z = -z;
      }
    }
    ClassDefNV(GridInfo, 1);
  };
 
  struct RowInfo {
    int32_t splineScenarioID{0};  
    size_t dataOffsetBytes[3]{0}; 
 
    GridInfo gridMeasured; 
    GridInfo gridReal;     
 
    ClassDefNV(RowInfo, 2);
  };
 
  typedef Spline2D<float, 3> SplineTypeXYZ;
  typedef Spline2D<float, 1> SplineTypeInvX;
  typedef Spline2D<float, 2> SplineTypeInvYZ;
 
  typedef SplineTypeXYZ SplineType;
 
  using SlimSplineTypeXYZ = Spline2D<float, 3, NoFlatObject>;
  using SlimSplineTypeInvX = Spline2D<float, 1, NoFlatObject>;
  using SlimSplineTypeInvYZ = Spline2D<float, 2, NoFlatObject>;
  using SlimSplineType = SlimSplineTypeXYZ;
 
 
  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();
 
  GPUdi() void setTimeStamp(int64_t v) { mTimeStamp = v; }
 
  GPUd() const SplineType& getSplineForRow(int32_t row) const;
 
  GPUd() SplineType& getSplineForRow(int32_t row);
 
  GPUd() float* getCorrectionData(int32_t sector, int32_t row, int32_t iSpline = 0);
 
  GPUd() const float* getCorrectionData(int32_t sector, int32_t row, int32_t iSpline = 0) const;
 
  GPUd() const SplineTypeInvX& getSplineInvXforRow(int32_t row) const;
 
  GPUd() SplineTypeInvX& getSplineInvXforRow(int32_t row);
 
  GPUd() float* getCorrectionDataInvX(int32_t sector, int32_t row);
 
  GPUd() const float* getCorrectionDataInvX(int32_t sector, int32_t row) const;
 
  GPUd() const SplineTypeInvYZ& getSplineInvYZforRow(int32_t row) const;
 
  GPUd() SplineTypeInvYZ& getSplineInvYZforRow(int32_t row);
 
  GPUd() float* getCorrectionDataInvYZ(int32_t sector, int32_t row);
 
  GPUd() const float* getCorrectionDataInvYZ(int32_t sector, int32_t row) const;
 
  // GPUd() int32_t getCorrectionInternal(int32_t sector, int32_t row, float u, float v, float& dx, float& du, float& dv) const;
 
  GPUdi() void getCorrectionLocal(int32_t sector, int32_t row, float y, float z, float& dx, float& dy, float& dz) const;
 
  GPUd() float getCorrectionXatRealYZ(int32_t sector, int32_t row, float realY, float realZ) const;
 
  GPUd() void getCorrectionYZatRealYZ(int32_t sector, int32_t row, float realY, float realZ, float& y, float& z) const;
 
 
  GPUd() void convLocalToGrid(int32_t sector, int32_t row, float y, float z, float& u, float& v, float& s) const;
 
  GPUd() void convGridToLocal(int32_t sector, int32_t row, float u, float v, float& y, float& z) const;
 
  GPUd() void convRealLocalToGrid(int32_t sector, int32_t row, float y, float z, float& u, float& v, float& s) const;
 
  GPUd() void convGridToRealLocal(int32_t sector, int32_t row, float u, float v, float& y, float& z) const;
 
  GPUd() bool isLocalInsideGrid(int32_t sector, int32_t row, float y, float z) const;
  GPUd() bool isRealLocalInsideGrid(int32_t sector, int32_t row, float y, float z) const;
 
  GPUdi() const TPCFastTransformGeo& getGeometry() const
  {
    return mGeo;
  }
 
  int64_t getTimeStamp() const { return mTimeStamp; }
 
  GPUdi() const RowInfo& getRowInfo(int32_t row) const
  {
    return mRowInfos[row];
  }
 
  GPUdi() RowInfo& getRowInfo(int32_t row)
  {
    return mRowInfos[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();
 
  static constexpr float kMaxCorrection = 100.f; 
 
 
 
  SplineType* mConstructionScenarios = nullptr; 
 
 
  TPCFastTransformGeo mGeo; 
 
  int32_t mNumberOfScenarios; 
 
  SplineType* mScenarioPtr; 
 
 
  int64_t mTimeStamp; 
 
  char* mCorrectionData[3]; 
 
  size_t mSectorDataSizeBytes[3]; 
 
  int32_t mClassVersion{3};
 
  RowInfo mRowInfos[TPCFastTransformGeo::getMaxNumberOfRows()]; 
 
  ClassDefNV(TPCFastSpaceChargeCorrection, 4);
};
 
 
GPUdi() const TPCFastSpaceChargeCorrection::SplineType& TPCFastSpaceChargeCorrection::getSplineForRow(int32_t row) const
{
  return mScenarioPtr[getRowInfo(row).splineScenarioID];
}
 
GPUdi() TPCFastSpaceChargeCorrection::SplineType& TPCFastSpaceChargeCorrection::getSplineForRow(int32_t row)
{
  return mScenarioPtr[getRowInfo(row).splineScenarioID];
}
 
GPUdi() float* TPCFastSpaceChargeCorrection::getCorrectionData(int32_t sector, int32_t row, int32_t iSpline)
{
  size_t offset = sector * mSectorDataSizeBytes[iSpline] + getRowInfo(row).dataOffsetBytes[iSpline];
  return reinterpret_cast<float*>(mCorrectionData[iSpline] + offset);
}
 
GPUdi() const float* TPCFastSpaceChargeCorrection::getCorrectionData(int32_t sector, int32_t row, int32_t iSpline) const
{
  size_t offset = sector * mSectorDataSizeBytes[iSpline] + getRowInfo(row).dataOffsetBytes[iSpline];
  return reinterpret_cast<const float*>(mCorrectionData[iSpline] + offset);
}
 
GPUdi() TPCFastSpaceChargeCorrection::SplineTypeInvX& TPCFastSpaceChargeCorrection::getSplineInvXforRow(int32_t row)
{
  return reinterpret_cast<SplineTypeInvX&>(getSplineForRow(row));
}
 
GPUdi() const TPCFastSpaceChargeCorrection::SplineTypeInvX& TPCFastSpaceChargeCorrection::getSplineInvXforRow(int32_t row) const
{
  return reinterpret_cast<const SplineTypeInvX&>(getSplineForRow(row));
}
 
GPUdi() float* TPCFastSpaceChargeCorrection::getCorrectionDataInvX(int32_t sector, int32_t row)
{
  return getCorrectionData(sector, row, 1);
}
 
GPUdi() const float* TPCFastSpaceChargeCorrection::getCorrectionDataInvX(int32_t sector, int32_t row) const
{
  return getCorrectionData(sector, row, 1);
}
 
GPUdi() TPCFastSpaceChargeCorrection::SplineTypeInvYZ& TPCFastSpaceChargeCorrection::getSplineInvYZforRow(int32_t row)
{
  return reinterpret_cast<SplineTypeInvYZ&>(getSplineForRow(row));
}
 
GPUdi() const TPCFastSpaceChargeCorrection::SplineTypeInvYZ& TPCFastSpaceChargeCorrection::getSplineInvYZforRow(int32_t row) const
{
  return reinterpret_cast<const SplineTypeInvYZ&>(getSplineForRow(row));
}
 
GPUdi() float* TPCFastSpaceChargeCorrection::getCorrectionDataInvYZ(int32_t sector, int32_t row)
{
  return getCorrectionData(sector, row, 2);
}
 
GPUdi() const float* TPCFastSpaceChargeCorrection::getCorrectionDataInvYZ(int32_t sector, int32_t row) const
{
  return getCorrectionData(sector, row, 2);
}
 
GPUdi() void TPCFastSpaceChargeCorrection::convLocalToGrid(int32_t sector, int32_t row, float y, float z, float& u, float& v, float& s) const
{
  const SplineType& spline = getSplineForRow(row);
  getRowInfo(row).gridMeasured.convLocalToGridUntruncated(sector, y, z, u, v, s);
  // shrink to the grid
  u = GPUCommonMath::Clamp(u, 0.f, (float)spline.getGridX1().getUmax());
  v = GPUCommonMath::Clamp(v, 0.f, (float)spline.getGridX2().getUmax());
}
 
GPUdi() bool TPCFastSpaceChargeCorrection::isLocalInsideGrid(int32_t sector, int32_t row, float y, float z) const
{
  float u, v, s;
  getRowInfo(row).gridMeasured.convLocalToGridUntruncated(sector, y, z, u, v, s);
  const auto& spline = getSplineForRow(row);
  // shrink to the grid
  if (u < 0.f || u > (float)spline.getGridX1().getUmax() || //
      v < 0.f || v > (float)spline.getGridX2().getUmax()) {
    return false;
  }
  return true;
}
 
GPUdi() bool TPCFastSpaceChargeCorrection::isRealLocalInsideGrid(int32_t sector, int32_t row, float y, float z) const
{
  float u, v, s;
  getRowInfo(row).gridReal.convLocalToGridUntruncated(sector, y, z, u, v, s);
  const auto& spline = getSplineForRow(row);
  // shrink to the grid
  if (u < 0.f || u > (float)spline.getGridX1().getUmax() || //
      v < 0.f || v > (float)spline.getGridX2().getUmax()) {
    return false;
  }
  return true;
}
 
GPUdi() void TPCFastSpaceChargeCorrection::convGridToLocal(int32_t sector, int32_t row, float gridU, float gridV, float& y, float& z) const
{
  getRowInfo(row).gridMeasured.convGridToLocal(sector, gridU, gridV, y, z);
}
 
GPUdi() void TPCFastSpaceChargeCorrection::convRealLocalToGrid(int32_t sector, int32_t row, float y, float z, float& u, float& v, float& s) const
{
  const SplineType& spline = getSplineForRow(row);
  getRowInfo(row).gridReal.convLocalToGridUntruncated(sector, y, z, u, v, s);
  // shrink to the grid
  u = GPUCommonMath::Clamp(u, 0.f, (float)spline.getGridX1().getUmax());
  v = GPUCommonMath::Clamp(v, 0.f, (float)spline.getGridX2().getUmax());
}
 
GPUdi() void TPCFastSpaceChargeCorrection::convGridToRealLocal(int32_t sector, int32_t row, float gridU, float gridV, float& y, float& z) const
{
  getRowInfo(row).gridReal.convGridToLocal(sector, gridU, gridV, y, z);
}
 
GPUdi() void TPCFastSpaceChargeCorrection::getCorrectionLocal(int32_t sector, int32_t row, float y, float z, float& dx, float& dy, float& dz) const
{
  const SplineType& spline = getSplineForRow(row);
  const float* splineData = getCorrectionData(sector, row);
 
  float u, v, s;
  convLocalToGrid(sector, row, y, z, u, v, s);
 
  float dxyz[3];
  spline.interpolateAtU(splineData, u, v, dxyz);
 
  if (CAMath::Abs(dxyz[0]) > kMaxCorrection || CAMath::Abs(dxyz[1]) > kMaxCorrection || CAMath::Abs(dxyz[2]) > kMaxCorrection) {
    s = 0.f; // TODO: DR: Protect from FPEs, fix upstream and remove once guaranteed that it is fixed
  }
 
  dx = s * dxyz[0];
  dy = s * dxyz[1];
  dz = s * dxyz[2];
}
 
GPUdi() float TPCFastSpaceChargeCorrection::getCorrectionXatRealYZ(int32_t sector, int32_t row, float realY, float realZ) const
{
  float u, v, s;
  convRealLocalToGrid(sector, row, realY, realZ, u, v, s);
  float dx = 0;
  getSplineInvXforRow(row).interpolateAtU(getCorrectionDataInvX(sector, row), u, v, &dx);
  if (CAMath::Abs(dx) > kMaxCorrection) {
    s = 0.f; // TODO: DR: Protect from FPEs, fix upstream and remove once guaranteed that it is fixed
  }
  return s * dx;
}
 
GPUdi() void TPCFastSpaceChargeCorrection::getCorrectionYZatRealYZ(int32_t sector, int32_t row, float realY, float realZ, float& y, float& z) const
{
  float u, v, s;
  convRealLocalToGrid(sector, row, realY, realZ, u, v, s);
  float dyz[2];
  getSplineInvYZforRow(row).interpolateAtU(getCorrectionDataInvYZ(sector, row), u, v, dyz);
  if (CAMath::Abs(dyz[0]) > kMaxCorrection || CAMath::Abs(dyz[1]) > kMaxCorrection) {
    s = 0.f; // TODO: DR: Protect from FPEs, fix upstream and remove once guaranteed that it is fixed
  }
  y = s * dyz[0];
  z = s * dyz[1];
}
 
} // namespace o2::gpu
 
#endif