RegularGrid3D.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_TPC_REGULARGRID3D_H_
#define ALICEO2_TPC_REGULARGRID3D_H_
 
#include "TPCSpaceCharge/Vector.h"
#include "Rtypes.h" // for ClassDefNV
#include "TPCSpaceCharge/SpaceChargeParameter.h"
#include <vector>
 
namespace o2
{
namespace tpc
{
 
template <typename DataT = double>
struct RegularGridHelper {
  DataT zmin;
  DataT rmin;
  DataT phimin;
  DataT spacingZ;
  DataT spacingR;
  DataT spacingPhi;
  ParamSpaceCharge params;
  ClassDefNV(RegularGridHelper, 1)
};
 
template <typename DataT = double>
struct RegularGrid3D {
 
 public:
  RegularGrid3D(const DataT zmin, const DataT rmin, const DataT phimin, const DataT spacingZ, const DataT spacingR, const DataT spacingPhi, const ParamSpaceCharge& params) : mParamGrid{params}, mMin{{zmin, rmin, phimin}}, mMax{{zmin + (mParamGrid.NZVertices - 1) * spacingZ, rmin + (mParamGrid.NRVertices - 1) * spacingR, phimin + (mParamGrid.NPhiVertices - 1) * spacingPhi}}, mSpacing{{spacingZ, spacingR, spacingPhi}}, mInvSpacing{{static_cast<DataT>(1 / spacingZ), static_cast<DataT>(1 / spacingR), static_cast<DataT>(1 / spacingPhi)}}
  {
    mRVertices.reserve(mParamGrid.NRVertices);
    mPhiVertices.reserve(mParamGrid.NPhiVertices);
    mZVertices.reserve(mParamGrid.NZVertices);
    initLists();
  }
 
  int getDeltaXDataIndex(const int deltaX) const { return deltaX; }
 
  int getDeltaYDataIndex(const int deltaY) const { return mParamGrid.NZVertices * deltaY; }
 
  int getDeltaZDataIndex(const int deltaZ) const { return deltaZ * mParamGrid.NRVertices * mParamGrid.NZVertices; }
 
  // same as above
  int getDeltaDataIndex(const int delta, const int dim) const;
 
  // check if the specified index for given dimension lies in the grid
  bool isIndexInGrid(const int index, const unsigned int dim) const { return index < 0 ? false : (index > (sNdim[dim] - 1) ? false : true); }
 
  size_t getN(unsigned int dim) const { return sNdim[dim]; }
  size_t getNZ() const { return sNdim[FZ]; }
  size_t getNR() const { return sNdim[FR]; }
  size_t getNPhi() const { return sNdim[FPHI]; }
 
  static constexpr unsigned int getDim() { return FDIM; }  
  static constexpr unsigned int getFZ() { return FZ; }     
  static constexpr unsigned int getFR() { return FR; }     
  static constexpr unsigned int getFPhi() { return FPHI; } 
 
  const Vector<DataT, 3>& getGridMin() const { return mMin; } 
  DataT getGridMinZ() const { return mMin[FZ]; }              
  DataT getGridMinR() const { return mMin[FR]; }              
  DataT getGridMinPhi() const { return mMin[FPHI]; }          
 
  DataT getGridMaxZ() const { return mMax[FZ]; }
  DataT getGridMaxR() const { return mMax[FR]; }
  DataT getGridMaxPhi() const { return mMax[FPHI]; }
 
  const Vector<DataT, 3>& getInvSpacing() const { return mInvSpacing; }
  DataT getInvSpacingZ() const { return mInvSpacing[FZ]; }
  DataT getInvSpacingR() const { return mInvSpacing[FR]; }
  DataT getInvSpacingPhi() const { return mInvSpacing[FPHI]; }
 
  DataT getSpacingZ() const { return mSpacing[FZ]; }
  DataT getSpacingR() const { return mSpacing[FR]; }
  DataT getSpacingPhi() const { return mSpacing[FPHI]; }
 
  // clamp coordinates to the grid (not circular)
  DataT clampToGrid(const DataT pos, const unsigned int dim) const;
 
  // clamp coordinates to the grid (not circular)
  DataT clampToGridRel(const DataT pos, const unsigned int dim) const;
 
  // clamp coordinates to the grid circular
  DataT clampToGridCircular(DataT pos, const unsigned int dim) const;
 
  // clamp coordinates to the grid circular
  DataT clampToGridCircularRel(DataT pos, const unsigned int dim) const;
 
  DataT getZVertex(const size_t vertexX) const { return mZVertices[vertexX]; }
 
  DataT getRVertex(const size_t vertexY) const { return mRVertices[vertexY]; }
 
  DataT getPhiVertex(const size_t vertexZ) const { return mPhiVertices[vertexZ]; }
 
  const Vector<DataT, 3>& getMaxIndices() const { return sMaxIndex; } 
 
  DataT getMaxIndexZ() const { return sMaxIndex[0]; }   
  DataT getMaxIndexR() const { return sMaxIndex[1]; }   
  DataT getMaxIndexPhi() const { return sMaxIndex[2]; } 
 
  const auto& getParamSC() const { return mParamGrid; }
 
  unsigned short getRVertices() const { return mParamGrid.NRVertices; }
  unsigned short getZVertices() const { return mParamGrid.NZVertices; }
  unsigned short getPhiVertices() const { return mParamGrid.NPhiVertices; }
 
  RegularGridHelper<double> getHelper() const { return RegularGridHelper<double>{getGridMinZ(), getGridMinR(), getGridMinPhi(), getSpacingZ(), getSpacingR(), getSpacingPhi(), mParamGrid}; }
 
 private:
  ParamSpaceCharge mParamGrid{};
  static constexpr unsigned int FDIM = 3;                                                                                                                                          
  static constexpr unsigned int FZ = 0;                                                                                                                                            
  static constexpr unsigned int FR = 1;                                                                                                                                            
  static constexpr unsigned int FPHI = 2;                                                                                                                                          
  Vector<DataT, FDIM> mMin{};                                                                                                                                                      
  Vector<DataT, FDIM> mMax{};                                                                                                                                                      
  Vector<DataT, FDIM> mSpacing{};                                                                                                                                                  
  Vector<DataT, FDIM> mInvSpacing{};                                                                                                                                               
  Vector<DataT, FDIM> sMaxIndex{{static_cast<DataT>(mParamGrid.NZVertices - 1.), static_cast<DataT>(mParamGrid.NRVertices - 1), static_cast<DataT>(mParamGrid.NPhiVertices - 1)}}; 
  Vector<int, FDIM> sNdim{{static_cast<int>(mParamGrid.NZVertices), static_cast<int>(mParamGrid.NRVertices), static_cast<int>(mParamGrid.NPhiVertices)}};                          
  std::vector<DataT> mZVertices{};                                                                                                                                                 
  std::vector<DataT> mRVertices{};                                                                                                                                                 
  std::vector<DataT> mPhiVertices{};                                                                                                                                               
 
  void initLists();
 
  ClassDefNV(RegularGrid3D, 3)
};
 
 
template <typename DataT>
DataT RegularGrid3D<DataT>::clampToGrid(const DataT pos, const unsigned int dim) const
{
  if (mMin[dim] < mMax[dim]) {
    if (pos < mMin[dim]) {
      return mMin[dim];
    } else if (pos > mMax[dim]) {
      return mMax[dim];
    }
  } else {
    if (pos > mMin[dim]) {
      return mMin[dim];
    } else if (pos < mMax[dim]) {
      return mMax[dim];
    }
  }
  return pos;
}
 
template <typename DataT>
DataT RegularGrid3D<DataT>::clampToGridRel(const DataT pos, const unsigned int dim) const
{
  if (pos < 0) {
    return 0;
  } else if (pos >= sMaxIndex[dim]) {
    return sMaxIndex[dim] - 1; // -1 return second last index. otherwise two additional points have to be extrapolated for tricubic interpolation
  }
  return pos;
}
 
template <typename DataT>
DataT RegularGrid3D<DataT>::clampToGridCircular(DataT pos, const unsigned int dim) const
{
  while (pos < mMin[dim]) {
    pos += mMax[dim] - mMin[dim] + mSpacing[dim];
  }
  while (pos >= mMax[dim] + mSpacing[dim]) {
    pos -= mMax[dim] + mSpacing[dim] - mMin[dim];
  }
  return pos;
}
 
template <typename DataT>
DataT RegularGrid3D<DataT>::clampToGridCircularRel(DataT pos, const unsigned int dim) const
{
  while (pos < 0) {
    pos += sNdim[dim];
  }
  while (pos > sNdim[dim]) {
    pos -= sNdim[dim];
  }
  if (pos == sNdim[dim]) {
    pos = 0;
  }
  return pos;
}
 
template <typename DataT>
void RegularGrid3D<DataT>::initLists()
{
  for (size_t i = 0; i < mParamGrid.NZVertices; ++i) {
    mZVertices.emplace_back(mMin[FZ] + i * mSpacing[FZ]);
  }
  for (size_t i = 0; i < mParamGrid.NRVertices; ++i) {
    mRVertices.emplace_back(mMin[FR] + i * mSpacing[FR]);
  }
  for (size_t i = 0; i < mParamGrid.NPhiVertices; ++i) {
    mPhiVertices.emplace_back(mMin[FPHI] + i * mSpacing[FPHI]);
  }
}
 
template <typename DataT>
int RegularGrid3D<DataT>::getDeltaDataIndex(const int delta, const int dim) const
{
  const int offset[FDIM]{1, mParamGrid.NZVertices, mParamGrid.NRVertices * mParamGrid.NZVertices};
  const int deltaIndex = delta * offset[dim];
  return deltaIndex;
}
 
} // namespace tpc
} // namespace o2
 
#endif