SpaceCharge.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_SPACECHARGE_H_
#define ALICEO2_TPC_SPACECHARGE_H_
 
#include "TPCSpaceCharge/TriCubic.h"
#include "TPCSpaceCharge/SpaceChargeHelpers.h"
#include "TPCSpaceCharge/PoissonSolverHelpers.h"
#include "TPCSpaceCharge/RegularGrid3D.h"
#include "TPCSpaceCharge/DataContainer3D.h"
#include "TPCSpaceCharge/SpaceChargeParameter.h"
#include "DataFormatsTPC/Defs.h"
 
class TH3;
class TH3D;
class TH3F;
class TH2F;
 
namespace o2
{
 
namespace parameters
{
class GRPMagField;
}
 
namespace utils
{
class TreeStreamRedirector;
}
 
namespace tpc
{
 
class Sector;
 
template <class T>
class CalDet;
 
enum class SCDistortionType : int {
  SCDistortionsConstant = 0, // space-charge distortions constant over time
  SCDistortionsRealistic = 1 // realistic evolution of space-charge distortions over time
};
 
 
template <typename DataT = double>
class SpaceCharge
{
  using RegularGrid = RegularGrid3D<DataT>;
  using DataContainer = DataContainer3D<DataT>;
  using GridProp = GridProperties<DataT>;
  using TriCubic = TriCubicInterpolator<DataT>;
  using TH3DataT = std::conditional_t<std::is_same<DataT, double>::value, TH3D, TH3F>; // datatype for TH3 (TH3F for DataT==float and TH3D for DataT==double)
 
 public:
  SpaceCharge(const int bfield, const unsigned short nZVertices = 129, const unsigned short nRVertices = 129, const unsigned short nPhiVertices = 360, const bool initBuffers = false);
 
  SpaceCharge();
 
  SpaceCharge(SpaceCharge&&) = default;
 
  SpaceCharge& operator=(SpaceCharge&&) = default;
 
  enum class IntegrationStrategy { Trapezoidal = 0,     
                                   Simpson = 1,         
                                   Root = 2,            
                                   SimpsonIterative = 3 
  };
 
  enum class Type {
    Distortions = 0, 
    Corrections = 1  
  };
 
  enum class GlobalDistType {
    Standard = 0, 
    Fast = 1,     
    None = 2      
  };
 
  enum class GlobalDistCorrMethod {
    LocalDistCorr,  
    ElectricalField 
  };
 
  enum class MisalignmentType {
    ShiftedClip, 
    RotatedClip, 
    RodShift     
  };
 
  enum class FCType {
    IFC, 
    OFC  
  };
 
  void fillChargeDensityFromHisto(const TH3& hisSCDensity3D);
 
  void fillChargeDensityFromHisto(const TH3& hisSCDensity3D_A, const TH3& hisSCDensity3D_C);
 
  void fillChargeDensityFromHisto(const char* file, const char* nameA, const char* nameC);
 
  void fillChargeDensityFromCalDet(const std::vector<CalDet<float>>& calSCDensity3D);
 
  static void convertIDCsToCharge(std::vector<CalDet<float>>& idcZero, const CalDet<float>& mapIBF, const float ionDriftTimeMS = 200, const bool normToPadArea = true);
 
  void fillChargeFromIDCs(std::vector<CalDet<float>>& idcZero, const CalDet<float>& mapIBF, const float ionDriftTimeMS = 200, const bool normToPadArea = true);
 
  void fillChargeFromCalDet(const std::vector<CalDet<float>>& calCharge3D);
 
  void fillChargeDensityFromFile(TFile& fInp, const char* name);
 
  void calculateDistortionsCorrections(const o2::tpc::Side side, const bool calcVectors = false);
 
  void setChargeDensityFromFormula(const AnalyticalFields<DataT>& formulaStruct);
 
  void setPotentialBoundaryFromFormula(const AnalyticalFields<DataT>& formulaStruct);
 
  void addBoundaryPotential(const SpaceCharge<DataT>& other, const Side side, const float scaling = 1);
 
  void resetBoundaryPotentialToZeroInRangeZ(float zMin, float zMax, const Side side);
 
  void setPotentialFromFormula(const AnalyticalFields<DataT>& formulaStruct);
 
  void mirrorPotential(const Side sideRef, const Side sideMirrored);
 
  void setSimOneSector();
 
  static void unsetSimOneSector();
 
  void setDefaultStaticDistortionsGEMFrameChargeUp(const Side side, const DataT deltaPotential = 1000);
 
  void setPotentialBoundaryGEMFrameAlongR(const std::function<DataT(DataT)>& potentialFunc, const Side side);
 
  void setPotentialBoundaryGEMFrameIROCBottomAlongPhi(const std::function<DataT(DataT)>& potentialFunc, const Side side) { setPotentialBoundaryGEMFrameAlongPhi(potentialFunc, GEMstack::IROCgem, true, side); }
 
  void setPotentialBoundaryGEMFrameIROCTopAlongPhi(const std::function<DataT(DataT)>& potentialFunc, const Side side) { setPotentialBoundaryGEMFrameAlongPhi(potentialFunc, GEMstack::IROCgem, false, side); }
 
  void setPotentialBoundaryGEMFrameOROC1BottomAlongPhi(const std::function<DataT(DataT)>& potentialFunc, const Side side) { setPotentialBoundaryGEMFrameAlongPhi(potentialFunc, GEMstack::OROC1gem, true, side); }
 
  void setPotentialBoundaryGEMFrameOROC1TopAlongPhi(const std::function<DataT(DataT)>& potentialFunc, const Side side) { setPotentialBoundaryGEMFrameAlongPhi(potentialFunc, GEMstack::OROC1gem, false, side); }
 
  void setPotentialBoundaryGEMFrameOROC2BottomAlongPhi(const std::function<DataT(DataT)>& potentialFunc, const Side side) { setPotentialBoundaryGEMFrameAlongPhi(potentialFunc, GEMstack::OROC2gem, true, side); }
 
  void setPotentialBoundaryGEMFrameOROC2TopAlongPhi(const std::function<DataT(DataT)>& potentialFunc, const Side side) { setPotentialBoundaryGEMFrameAlongPhi(potentialFunc, GEMstack::OROC2gem, false, side); }
 
  void setPotentialBoundaryGEMFrameOROC3BottomAlongPhi(const std::function<DataT(DataT)>& potentialFunc, const Side side) { setPotentialBoundaryGEMFrameAlongPhi(potentialFunc, GEMstack::OROC3gem, true, side); }
 
  void setPotentialBoundaryGEMFrameOROC3TopAlongPhi(const std::function<DataT(DataT)>& potentialFunc, const Side side) { setPotentialBoundaryGEMFrameAlongPhi(potentialFunc, GEMstack::OROC3gem, false, side); }
 
  void setPotentialBoundaryGEMFrameOROC3ToOFCPhi(const std::function<DataT(DataT)>& potentialFunc, const Side side) { setPotentialBoundaryGEMFrameAlongPhi(potentialFunc, GEMstack::OROC3gem, false, side, true); }
 
  void setPotentialBoundaryGEMFrameIROCToIFCPhi(const std::function<DataT(DataT)>& potentialFunc, const Side side) { setPotentialBoundaryGEMFrameAlongPhi(potentialFunc, GEMstack::IROCgem, true, side, true); }
 
  void setPotentialBoundaryInnerRadius(const std::function<DataT(DataT)>& potentialFunc, const Side side);
 
  void setPotentialBoundaryOuterRadius(const std::function<DataT(DataT)>& potentialFunc, const Side side);
 
  void poissonSolver(const Side side, const DataT stoppingConvergence = 1e-6, const int symmetry = 0);
 
  void poissonSolver(const DataT stoppingConvergence = 1e-6, const int symmetry = 0);
 
  void calcEField(const Side side);
 
  void setEFieldFromFormula(const AnalyticalFields<DataT>& formulaStruct);
 
  void scaleChargeDensity(const DataT scalingFactor, const Side side) { mDensity[side] *= scalingFactor; }
 
  void scaleChargeDensitySector(const float scalingFactor, const Sector sector);
 
  void scaleChargeDensityStack(const float scalingFactor, const Sector sector, const GEMstack stack);
 
  void addChargeDensity(const SpaceCharge<DataT>& otherSC);
 
  template <typename ElectricFields = AnalyticalFields<DataT>>
  void calcLocalDistortionsCorrections(const Type type, const ElectricFields& formulaStruct);
 
  template <typename ElectricFields = AnalyticalFields<DataT>>
  void calcLocalDistortionCorrectionVector(const ElectricFields& formulaStruct);
 
  template <typename ElectricFields = AnalyticalFields<DataT>>
  void calcLocalDistortionsCorrectionsRK4(const Type type, const Side side);
 
  template <typename Fields = AnalyticalFields<DataT>>
  void calcGlobalCorrections(const Fields& formulaStruct, const int type = 3);
 
  void calcGlobalCorrectionsEField(const Side side, const int type = 3) { calcGlobalCorrections(getElectricFieldsInterpolator(side), type); }
 
  template <typename Fields = AnalyticalFields<DataT>>
  void calcGlobalDistortions(const Fields& formulaStruct, const int maxIterations = 3 * sSteps * 129);
 
  void init();
 
  void calcGlobalDistWithGlobalCorrIterative(const DistCorrInterpolator<DataT>& globCorr, const int maxIter = 100, const DataT approachZ = 1, const DataT approachR = 1, const DataT approachPhi = 1, const DataT diffCorr = 50e-6, const SpaceCharge<DataT>* scSCale = nullptr, float scale = 0);
 
  void calcGlobalDistWithGlobalCorrIterative(const Side side, const SpaceCharge<DataT>* scSCale = nullptr, float scale = 0, const int maxIter = 100, const DataT approachZ = 1, const DataT approachR = 1, const DataT approachPhi = 1, const DataT diffCorr = 50e-6);
  void calcGlobalDistWithGlobalCorrIterative(const SpaceCharge<DataT>* scSCale = nullptr, float scale = 0, const int maxIter = 100, const DataT approachZ = 1, const DataT approachR = 1, const DataT approachPhi = 1, const DataT diffCorr = 50e-6);
 
  void calcGlobalCorrWithGlobalDistIterative(const Side side, const SpaceCharge<DataT>* scSCale = nullptr, float scale = 0, const int maxIter = 100, const DataT approachZ = 1, const DataT approachR = 1, const DataT approachPhi = 1, const DataT diffCorr = 50e-6);
  void calcGlobalCorrWithGlobalDistIterative(const SpaceCharge<DataT>* scSCale = nullptr, float scale = 0, const int maxIter = 100, const DataT approachZ = 1, const DataT approachR = 1, const DataT approachPhi = 1, const DataT diffCorr = 50e-6);
 
  void calcGlobalDistWithGlobalCorrIterativeLinearCartesian(const Side side, const SpaceCharge<DataT>* scSCale = nullptr, float scale = 0, const int maxIter = 100, const DataT approachX = 1, const DataT approachY = 1, const DataT approachZ = 1, const DataT diffCorr = 50e-6);
  void calcGlobalDistWithGlobalCorrIterativeLinearCartesian(const SpaceCharge<DataT>* scSCale = nullptr, float scale = 0, const int maxIter = 100, const DataT approachX = 1, const DataT approachY = 1, const DataT approachZ = 1, const DataT diffCorr = 50e-6);
 
  void calcGlobalCorrWithGlobalDistIterativeLinearCartesian(const Side side, const SpaceCharge<DataT>* scSCale = nullptr, float scale = 0, const int maxIter = 100, const DataT approachX = 1, const DataT approachY = 1, const DataT approachZ = 1, const DataT diffCorr = 50e-6);
  void calcGlobalCorrWithGlobalDistIterativeLinearCartesian(const SpaceCharge<DataT>* scSCale = nullptr, float scale = 0, const int maxIter = 100, const DataT approachX = 1, const DataT approachY = 1, const DataT approachZ = 1, const DataT diffCorr = 50e-6);
 
  unsigned short getNZVertices() const { return mParamGrid.NZVertices; }
 
  unsigned short getNRVertices() const { return mParamGrid.NRVertices; }
 
  unsigned short getNPhiVertices() const { return mParamGrid.NPhiVertices; }
 
  auto getC0() const { return mC0; }
 
  auto getC1() const { return mC1; }
 
  auto getC2() const { return mC2; }
 
  int getBField() const { return mBField.getBField(); }
 
  const auto& getPotential(const Side side) const& { return mPotential[side]; }
 
  void setPotential(int iz, int ir, int iphi, Side side, float val);
 
  DataT getDensityCyl(const DataT z, const DataT r, const DataT phi, const Side side) const;
 
  DataT getPotentialCyl(const DataT z, const DataT r, const DataT phi, const Side side) const;
 
  std::vector<float> getPotentialCyl(const std::vector<DataT>& z, const std::vector<DataT>& r, const std::vector<DataT>& phi, const Side side) const;
 
  void getElectricFieldsCyl(const DataT z, const DataT r, const DataT phi, const Side side, DataT& eZ, DataT& eR, DataT& ePhi) const;
 
  void getLocalCorrectionsCyl(const DataT z, const DataT r, const DataT phi, const Side side, DataT& lcorrZ, DataT& lcorrR, DataT& lcorrRPhi) const;
 
  void getLocalCorrectionsCyl(const std::vector<DataT>& z, const std::vector<DataT>& r, const std::vector<DataT>& phi, const Side side, std::vector<DataT>& lcorrZ, std::vector<DataT>& lcorrR, std::vector<DataT>& lcorrRPhi) const;
 
  void getCorrectionsCyl(const DataT z, const DataT r, const DataT phi, const Side side, DataT& corrZ, DataT& corrR, DataT& corrRPhi) const;
 
  void getCorrectionsCyl(const std::vector<DataT>& z, const std::vector<DataT>& r, const std::vector<DataT>& phi, const Side side, std::vector<DataT>& corrZ, std::vector<DataT>& corrR, std::vector<DataT>& corrRPhi) const;
 
  void getCorrections(const DataT x, const DataT y, const DataT z, const Side side, DataT& corrX, DataT& corrY, DataT& corrZ) const;
 
  void getCorrectionsAnalytical(const DataT x, const DataT y, const DataT z, const Side side, DataT& corrX, DataT& corrY, DataT& corrZ) const { getDistortionsCorrectionsAnalytical(x, y, z, side, corrX, corrY, corrZ, false); }
 
  void getLocalDistortionsCyl(const DataT z, const DataT r, const DataT phi, const Side side, DataT& ldistZ, DataT& ldistR, DataT& ldistRPhi) const;
 
  void getLocalDistortionsCyl(const std::vector<DataT>& z, const std::vector<DataT>& r, const std::vector<DataT>& phi, const Side side, std::vector<DataT>& ldistZ, std::vector<DataT>& ldistR, std::vector<DataT>& ldistRPhi) const;
 
  void getLocalDistortionVectorCyl(const DataT z, const DataT r, const DataT phi, const Side side, DataT& lvecdistZ, DataT& lvecdistR, DataT& lvecdistRPhi) const;
 
  void getLocalDistortionVectorCyl(const std::vector<DataT>& z, const std::vector<DataT>& r, const std::vector<DataT>& phi, const Side side, std::vector<DataT>& lvecdistZ, std::vector<DataT>& lvecdistR, std::vector<DataT>& lvecdistRPhi) const;
 
  void getLocalCorrectionVectorCyl(const DataT z, const DataT r, const DataT phi, const Side side, DataT& lveccorrZ, DataT& lveccorrR, DataT& lveccorrRPhi) const;
 
  void getLocalCorrectionVectorCyl(const std::vector<DataT>& z, const std::vector<DataT>& r, const std::vector<DataT>& phi, const Side side, std::vector<DataT>& lveccorrZ, std::vector<DataT>& lveccorrR, std::vector<DataT>& lveccorrRPhi) const;
 
  void getDistortionsCyl(const DataT z, const DataT r, const DataT phi, const Side side, DataT& distZ, DataT& distR, DataT& distRPhi) const;
 
  void getDistortionsCyl(const std::vector<DataT>& z, const std::vector<DataT>& r, const std::vector<DataT>& phi, const Side side, std::vector<DataT>& distZ, std::vector<DataT>& distR, std::vector<DataT>& distRPhi) const;
 
  void getDistortions(const DataT x, const DataT y, const DataT z, const Side side, DataT& distX, DataT& distY, DataT& distZ) const;
 
  void getDistortionsAnalytical(const DataT x, const DataT y, const DataT z, const Side side, DataT& distX, DataT& distY, DataT& distZ) const { getDistortionsCorrectionsAnalytical(x, y, z, side, distX, distY, distZ, true); }
 
  void setDistortionsCorrectionsAnalytical(const AnalyticalDistCorr<DataT>& formula) { mAnaDistCorr = formula; }
 
  const auto& getDistortionsCorrectionsAnalytical() const& { return mAnaDistCorr; }
 
  void setUseAnalyticalDistCorr(const bool useAna) { mUseAnaDistCorr = useAna; }
 
  bool getUseAnalyticalDistCorr() const { return mUseAnaDistCorr; }
 
  static DataT getRadiusFromCartesian(const DataT x, const DataT y) { return std::sqrt(x * x + y * y); }
 
  static DataT getPhiFromCartesian(const DataT x, const DataT y) { return std::atan2(y, x); }
 
  static DataT getXFromPolar(const DataT r, const DataT phi) { return r * std::cos(phi); }
 
  static DataT getYFromPolar(const DataT r, const DataT phi) { return r * std::sin(phi); }
 
  void correctElectron(GlobalPosition3D& point);
 
  void distortElectron(GlobalPosition3D& point, const SpaceCharge<DataT>* scSCale = nullptr, float scale = 0) const;
 
  void setDistortionLookupTables(const DataContainer& distdZ, const DataContainer& distdR, const DataContainer& distdRPhi, const Side side);
 
  void setFromFile(std::string_view file, const Side side);
 
  void setFromFile(std::string_view file);
 
  DataT getGridSpacingR(const Side side) const { return mGrid3D[side].getSpacingR(); }
 
  DataT getGridSpacingZ(const Side side) const { return mGrid3D[side].getSpacingZ(); }
 
  DataT getGridSpacingPhi(const Side side) const { return mGrid3D[side].getSpacingPhi(); }
 
  static constexpr DataT getEzField(const Side side) { return getSign(side) * (TPCParameters<DataT>::cathodev - TPCParameters<DataT>::vg1t) / TPCParameters<DataT>::TPCZ0; }
 
  DataT getRMin(const Side side) const { return mGrid3D[side].getGridMinR(); }
 
  DataT getZMin(const Side side) const { return mGrid3D[side].getGridMinZ(); }
 
  DataT getPhiMin(const Side side) const { return mGrid3D[side].getGridMinPhi(); }
 
  DataT getRMax(const Side side) const { return mGrid3D[side].getGridMaxR(); };
 
  DataT getZMax(const Side side) const { return mGrid3D[side].getGridMaxZ(); }
 
  DataT getPhiMax(const Side side) const { return mGrid3D[side].getGridMaxPhi(); }
 
  // get side of TPC for z coordinate TODO rewrite this
  static Side getSide(const DataT z) { return ((z >= 0) ? Side::A : Side::C); }
 
  const RegularGrid& getGrid3D(const Side side) const { return mGrid3D[side]; }
 
  NumericalFields<DataT> getElectricFieldsInterpolator(const Side side) const;
 
  DistCorrInterpolator<DataT> getLocalDistInterpolator(const Side side) const;
 
  DistCorrInterpolator<DataT> getLocalCorrInterpolator(const Side side) const;
 
  DistCorrInterpolator<DataT> getGlobalDistInterpolator(const Side side) const;
 
  DistCorrInterpolator<DataT> getGlobalCorrInterpolator(const Side side) const;
 
  DataT getLocalDistR(const size_t iz, const size_t ir, const size_t iphi, const Side side) const { return mLocalDistdR[side](iz, ir, iphi); }
 
  DataT getLocalDistZ(const size_t iz, const size_t ir, const size_t iphi, const Side side) const { return mLocalDistdZ[side](iz, ir, iphi); }
 
  DataT getLocalDistRPhi(const size_t iz, const size_t ir, const size_t iphi, const Side side) const { return mLocalDistdRPhi[side](iz, ir, iphi); }
 
  DataT getLocalVecDistR(const size_t iz, const size_t ir, const size_t iphi, const Side side) const { return mLocalVecDistdR[side](iz, ir, iphi); }
 
  DataT getLocalVecDistZ(const size_t iz, const size_t ir, const size_t iphi, const Side side) const { return mLocalVecDistdZ[side](iz, ir, iphi); }
 
  DataT getLocalVecDistRPhi(const size_t iz, const size_t ir, const size_t iphi, const Side side) const { return mLocalVecDistdRPhi[side](iz, ir, iphi); }
 
  DataT getLocalCorrR(const size_t iz, const size_t ir, const size_t iphi, const Side side) const { return mLocalCorrdR[side](iz, ir, iphi); }
 
  DataT getLocalCorrZ(const size_t iz, const size_t ir, const size_t iphi, const Side side) const { return mLocalCorrdZ[side](iz, ir, iphi); }
 
  DataT getLocalCorrRPhi(const size_t iz, const size_t ir, const size_t iphi, const Side side) const { return mLocalCorrdRPhi[side](iz, ir, iphi); }
 
  DataT getLocalVecCorrR(const size_t iz, const size_t ir, const size_t iphi, const Side side) const { return -mLocalVecDistdR[side](iz, ir, iphi); }
 
  DataT getLocalVecCorrZ(const size_t iz, const size_t ir, const size_t iphi, const Side side) const { return -mLocalVecDistdZ[side](iz, ir, iphi); }
 
  DataT getLocalVecCorrRPhi(const size_t iz, const size_t ir, const size_t iphi, const Side side) const { return -mLocalVecDistdRPhi[side](iz, ir, iphi); }
 
  DataT getGlobalDistR(const size_t iz, const size_t ir, const size_t iphi, const Side side) const { return mGlobalDistdR[side](iz, ir, iphi); }
 
  DataT getGlobalDistZ(const size_t iz, const size_t ir, const size_t iphi, const Side side) const { return mGlobalDistdZ[side](iz, ir, iphi); }
 
  DataT getGlobalDistRPhi(const size_t iz, const size_t ir, const size_t iphi, const Side side) const { return mGlobalDistdRPhi[side](iz, ir, iphi); }
 
  DataT getGlobalCorrR(const size_t iz, const size_t ir, const size_t iphi, const Side side) const { return mGlobalCorrdR[side](iz, ir, iphi); }
 
  DataT getGlobalCorrZ(const size_t iz, const size_t ir, const size_t iphi, const Side side) const { return mGlobalCorrdZ[side](iz, ir, iphi); }
 
  DataT getGlobalCorrRPhi(const size_t iz, const size_t ir, const size_t iphi, const Side side) const { return mGlobalCorrdRPhi[side](iz, ir, iphi); }
 
  DataT getEr(const size_t iz, const size_t ir, const size_t iphi, const Side side) const { return mElectricFieldEr[side](iz, ir, iphi); }
 
  DataT getEz(const size_t iz, const size_t ir, const size_t iphi, const Side side) const { return mElectricFieldEz[side](iz, ir, iphi); }
 
  DataT getEphi(const size_t iz, const size_t ir, const size_t iphi, const Side side) const { return mElectricFieldEphi[side](iz, ir, iphi); }
 
  DataT getDensity(const size_t iz, const size_t ir, const size_t iphi, const Side side) const { return mDensity[side](iz, ir, iphi); }
 
  DataT getPotential(const size_t iz, const size_t ir, const size_t iphi, const Side side) const { return mPotential[side](iz, ir, iphi); }
 
  static int getStepWidth() { return 1 / sSteps; }
 
  DataT getPhiVertex(const size_t indexPhi, const Side side) const { return mGrid3D[side].getPhiVertex(indexPhi); }
 
  DataT getRVertex(const size_t indexR, const Side side) const { return mGrid3D[side].getRVertex(indexR); }
 
  DataT getZVertex(const size_t indexZ, const Side side) const { return mGrid3D[side].getZVertex(indexZ); }
 
  void setOmegaTauT1T2(const DataT omegaTau = 0.32f, const DataT t1 = 1, const DataT t2 = 1);
 
  void setC0C1(const DataT c0, const DataT c1)
  {
    mC0 = c0;
    mC1 = c1;
  }
 
  void setC2(const DataT c2) { mC0 = c2; }
 
  void initBField(const int field);
 
  void setSimExBMisalignment(const bool simExBMisalignment) { mSimExBMisalignment = simExBMisalignment; }
 
  bool getSimExBMisalignment() const { return mSimExBMisalignment; }
 
  void setSimEDistortions(const bool simEDistortions) { mSimEDistortions = simEDistortions; }
 
  bool getSimEDistortions() const { return mSimEDistortions; }
 
  static void setNStep(const int nSteps) { sSteps = nSteps; }
 
  static int getNStep() { return sSteps; }
 
  static int getNThreads() { return sNThreads; }
 
  static void setNThreads(const int nThreads) { sNThreads = nThreads; }
 
  static void setNumericalIntegrationStrategy(const IntegrationStrategy strategy) { sNumericalIntegrationStrategy = strategy; }
  static IntegrationStrategy getNumericalIntegrationStrategy() { return sNumericalIntegrationStrategy; }
 
  static void setGlobalDistType(const GlobalDistType globalDistType) { sGlobalDistType = globalDistType; }
  static GlobalDistType getGlobalDistType() { return sGlobalDistType; }
 
  static void setGlobalDistCorrMethod(const GlobalDistCorrMethod globalDistCorrMethod) { sGlobalDistCorrCalcMethod = globalDistCorrMethod; }
  static GlobalDistCorrMethod getGlobalDistCorrMethod() { return sGlobalDistCorrCalcMethod; }
 
  static void setSimpsonNIteratives(const int nIter) { sSimpsonNIteratives = nIter; }
  static int getSimpsonNIteratives() { return sSimpsonNIteratives; }
 
  static void setSCDistortionType(SCDistortionType distortionType) { sSCDistortionType = distortionType; }
  static SCDistortionType getSCDistortionType() { return sSCDistortionType; }
 
  void setUseInitialSCDensity(const bool useInitialSCDensity) { mUseInitialSCDensity = useInitialSCDensity; }
 
  void dumpToFile(std::string_view file, const Side side, std::string_view option) const;
 
  void dumpToFile(std::string_view file) const;
 
  void dumpMetaData(std::string_view file, std::string_view option, const bool overwriteExisting = false) const;
 
  void readMetaData(std::string_view file);
 
  void dumpToTree(const char* outFileName, const Side side, const int nZPoints = 50, const int nRPoints = 150, const int nPhiPoints = 180, const bool randomize = false) const;
 
  void dumpToTree(const char* outFileName, const o2::tpc::Sector& sector, const int nZPoints = 50) const;
 
  int dumpElectricFields(std::string_view file, const Side side, std::string_view option) const;
 
  void setElectricFieldsFromFile(std::string_view file, const Side side);
 
  int dumpPotential(std::string_view file, const Side side, std::string_view option) const;
 
  void setPotentialFromFile(std::string_view file, const Side side);
 
  void fillPotential(const DataT potential, const size_t iz, const size_t ir, const size_t iphi, const Side side) { mPotential[side](iz, ir, iphi) = potential; }
 
  int dumpDensity(std::string_view file, const Side side, std::string_view option) const;
 
  void setDensityFromFile(std::string_view file, const Side side);
 
  void fillDensity(const DataT density, const size_t iz, const size_t ir, const size_t iphi, const Side side) { mDensity[side](iz, ir, iphi) = density; }
 
  void averageDensityPerSector(const Side side);
 
  int dumpGlobalDistortions(std::string_view file, const Side side, std::string_view option) const;
 
  int dumpAnalyticalCorrectionsDistortions(TFile& outf) const;
 
  void setAnalyticalCorrectionsDistortionsFromFile(std::string_view inpf);
 
  void setGlobalDistortionsFromFile(std::string_view file, const Side side);
 
  template <typename DataTIn = DataT>
  void setGlobalDistortionsFromFile(TFile& inpf, const Side side);
 
  int dumpGlobalCorrections(std::string_view file, const Side side, std::string_view option) const;
 
  int dumpGlobalCorrections(TFile& outf, const Side side) const;
 
  void setGlobalCorrectionsFromFile(std::string_view file, const Side side);
 
  template <typename DataTIn = DataT>
  void setGlobalCorrectionsFromFile(TFile& inpf, const Side side);
 
  int dumpLocalCorrections(std::string_view file, const Side side, std::string_view option) const;
 
  void setLocalCorrectionsFromFile(std::string_view file, const Side side);
 
  int dumpLocalDistortions(std::string_view file, const Side side, std::string_view option) const;
 
  int dumpLocalDistCorrVectors(std::string_view file, const Side side, std::string_view option) const;
 
  void setLocalDistortionsFromFile(std::string_view file, const Side side);
 
  void setLocalDistCorrVectorsFromFile(std::string_view file, const Side side);
 
  DataT regulateZ(const DataT posZ, const Side side) const { return mGrid3D[side].clampToGrid(posZ, 0); }
 
  DataT regulateR(const DataT posR, const Side side) const;
 
  DataT regulatePhi(const DataT posPhi, const Side side) const { return mGrid3D[side].clampToGridCircular(posPhi, 2); }
 
  std::vector<std::pair<std::vector<o2::math_utils::Point3D<float>>, std::array<DataT, 3>>> calculateElectronDriftPath(const std::vector<GlobalPosition3D>& elePos, const int nSamplingPoints, const std::string_view outFile = "electron_tracks.root") const;
 
  static void makeElectronDriftPathGif(const char* inpFile, TH2F& hBorder, const int type = 0, const int gifSpeed = 2, const int maxsamplingpoints = 100, const char* outName = "electron_drift_path");
 
  static void normalizeHistoQVEps0(TH3& histoIonsPhiRZ);
 
  static int getOMPMaxThreads();
 
  bool checkGridFromFile(std::string_view file, std::string_view tree);
 
  void setDeltaVoltageCopperRodShiftIFC(const float deltaPot, const Side side, const int sector) { initRodAlignmentVoltages(MisalignmentType::RodShift, FCType::IFC, sector, side, deltaPot); }
 
  void setDeltaVoltageCopperRodShiftOFC(const float deltaPot, const Side side, const int sector) { initRodAlignmentVoltages(MisalignmentType::RodShift, FCType::OFC, sector, side, deltaPot); }
 
  void setDeltaVoltageStripsShiftIFC(const float deltaPot, const Side side, const int sector) { initRodAlignmentVoltages(MisalignmentType::ShiftedClip, FCType::IFC, sector, side, deltaPot); }
 
  void setDeltaVoltageStripsShiftOFC(const float deltaPot, const Side side, const int sector) { initRodAlignmentVoltages(MisalignmentType::ShiftedClip, FCType::OFC, sector, side, deltaPot); }
 
  void setDeltaVoltageRotatedClipIFC(const float deltaPot, const Side side, const int sector) { initRodAlignmentVoltages(MisalignmentType::RotatedClip, FCType::IFC, sector, side, deltaPot); }
 
  void setDeltaVoltageRotatedClipOFC(const float deltaPot, const Side side, const int sector) { initRodAlignmentVoltages(MisalignmentType::RotatedClip, FCType::OFC, sector, side, deltaPot); }
 
  void setIFCChargeUpRisingPot(const float deltaPot, const float zMaxDeltaPot, const int type, const float zStart, const float offs, const Side side);
 
  void setIFCChargeUpFallingPot(const float deltaPot, const float zMaxDeltaPot, const int type, const float zEnd, const float offs, const Side side);
 
  void setGlobalCorrections(const std::function<void(int sector, DataT gx, DataT gy, DataT gz, DataT& gCx, DataT& gCy, DataT& gCz)>& gCorr, const Side side);
 
  void setROCMisalignmentShiftZ(const int sector, const int type, const float potential);
 
  void setROCMisalignmentRotationAlongX(const int sector, const int type, const float potentialMin, const float potentialMax);
 
  void setROCMisalignmentRotationAlongY(const int sector, const int type, const float potentialMin, const float potentialMax);
 
  void subtractGlobalCorrections(const SpaceCharge<DataT>& otherSC, const Side side);
 
  void subtractGlobalDistortions(const SpaceCharge<DataT>& otherSC, const Side side);
 
  void scaleCorrections(const float scaleFac, const Side side);
 
  void setMetaData(const SCMetaData& meta) { mMeta = meta; }
  const auto& getMetaData() const { return mMeta; }
  void printMetaData() const { mMeta.print(); }
  float getMeanLumi() const { return mMeta.meanLumi; }
  void setMeanLumi(float lumi) { mMeta.meanLumi = lumi; }
  void initAfterReadingFromFile();
 
  float getDCAr(float tgl, const int nPoints, const float phi, o2::utils::TreeStreamRedirector* pcstream = nullptr) const;
 
 private:
  ParamSpaceCharge mParamGrid{};                                                                          
  inline static int sNThreads{getOMPMaxThreads()};                                                        
  inline static IntegrationStrategy sNumericalIntegrationStrategy{IntegrationStrategy::SimpsonIterative}; 
  inline static int sSimpsonNIteratives{3};                                                               
  inline static int sSteps{1};                                                                            
  inline static GlobalDistType sGlobalDistType{GlobalDistType::Fast};                                     
  inline static GlobalDistCorrMethod sGlobalDistCorrCalcMethod{GlobalDistCorrMethod::LocalDistCorr};      
  inline static SCDistortionType sSCDistortionType{SCDistortionType::SCDistortionsConstant};              
 
  DataT mC0 = 0;                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                          
  DataT mC1 = 0;                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                          
  DataT mC2 = 0;                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                          
  static constexpr int FNSIDES = SIDES;                                                                                                                                                                                                                                                                                                                                                                                                                                                                                   
  bool mUseInitialSCDensity{false};                                                                                                                                                                                                                                                                                                                                                                                                                                                                                       
  bool mInitLookUpTables{false};                                                                                                                                                                                                                                                                                                                                                                                                                                                                                          
  bool mSimExBMisalignment{false};                                                                                                                                                                                                                                                                                                                                                                                                                                                                                        
  bool mSimEDistortions{true};                                                                                                                                                                                                                                                                                                                                                                                                                                                                                            
  bool mReadMetaData{false};                                                                                                                                                                                                                                                                                                                                                                                                                                                                                              
  RegularGrid mGrid3D[FNSIDES]{{GridProp::ZMIN, GridProp::RMIN, GridProp::PHIMIN, getSign(Side::A) * GridProp::getGridSpacingZ(mParamGrid.NZVertices), GridProp::getGridSpacingR(mParamGrid.NRVertices), GridProp::getGridSpacingPhi(mParamGrid.NPhiVertices), mParamGrid}, {GridProp::ZMIN, GridProp::RMIN, GridProp::PHIMIN, getSign(Side::C) * GridProp::getGridSpacingZ(mParamGrid.NZVertices), GridProp::getGridSpacingR(mParamGrid.NRVertices), GridProp::getGridSpacingPhi(mParamGrid.NPhiVertices), mParamGrid}}; 
 
  DataContainer mLocalDistdR[FNSIDES]{};       
  DataContainer mLocalDistdZ[FNSIDES]{};       
  DataContainer mLocalDistdRPhi[FNSIDES]{};    
  DataContainer mLocalVecDistdR[FNSIDES]{};    
  DataContainer mLocalVecDistdZ[FNSIDES]{};    
  DataContainer mLocalVecDistdRPhi[FNSIDES]{}; 
  DataContainer mLocalCorrdR[FNSIDES]{};       
  DataContainer mLocalCorrdZ[FNSIDES]{};       
  DataContainer mLocalCorrdRPhi[FNSIDES]{};    
  DataContainer mGlobalDistdR[FNSIDES]{};      
  DataContainer mGlobalDistdZ[FNSIDES]{};      
  DataContainer mGlobalDistdRPhi[FNSIDES]{};   
  DataContainer mGlobalCorrdR[FNSIDES]{};      
  DataContainer mGlobalCorrdZ[FNSIDES]{};      
  DataContainer mGlobalCorrdRPhi[FNSIDES]{};   
  DataContainer mDensity[FNSIDES]{};           
  DataContainer mPotential[FNSIDES]{};         
  DataContainer mElectricFieldEr[FNSIDES]{};   
  DataContainer mElectricFieldEz[FNSIDES]{};   
  DataContainer mElectricFieldEphi[FNSIDES]{}; 
 
  TriCubic mInterpolatorPotential[FNSIDES]{{mPotential[Side::A], mGrid3D[Side::A]}, {mPotential[Side::C], mGrid3D[Side::C]}};                                                                                                                                                                 
  TriCubic mInterpolatorDensity[FNSIDES]{{mDensity[Side::A], mGrid3D[Side::A]}, {mDensity[Side::C], mGrid3D[Side::C]}};                                                                                                                                                                       
  DistCorrInterpolator<DataT> mInterpolatorGlobalCorr[FNSIDES]{{mGlobalCorrdR[Side::A], mGlobalCorrdZ[Side::A], mGlobalCorrdRPhi[Side::A], mGrid3D[Side::A], Side::A}, {mGlobalCorrdR[Side::C], mGlobalCorrdZ[Side::C], mGlobalCorrdRPhi[Side::C], mGrid3D[Side::C], Side::C}};               
  DistCorrInterpolator<DataT> mInterpolatorLocalCorr[FNSIDES]{{mLocalCorrdR[Side::A], mLocalCorrdZ[Side::A], mLocalCorrdRPhi[Side::A], mGrid3D[Side::A], Side::A}, {mLocalCorrdR[Side::C], mLocalCorrdZ[Side::C], mLocalCorrdRPhi[Side::C], mGrid3D[Side::C], Side::C}};                      
  DistCorrInterpolator<DataT> mInterpolatorGlobalDist[FNSIDES]{{mGlobalDistdR[Side::A], mGlobalDistdZ[Side::A], mGlobalDistdRPhi[Side::A], mGrid3D[Side::A], Side::A}, {mGlobalDistdR[Side::C], mGlobalDistdZ[Side::C], mGlobalDistdRPhi[Side::C], mGrid3D[Side::C], Side::C}};               
  DistCorrInterpolator<DataT> mInterpolatorLocalDist[FNSIDES]{{mLocalDistdR[Side::A], mLocalDistdZ[Side::A], mLocalDistdRPhi[Side::A], mGrid3D[Side::A], Side::A}, {mLocalDistdR[Side::C], mLocalDistdZ[Side::C], mLocalDistdRPhi[Side::C], mGrid3D[Side::C], Side::C}};                      
  DistCorrInterpolator<DataT> mInterpolatorLocalVecDist[FNSIDES]{{mLocalVecDistdR[Side::A], mLocalVecDistdZ[Side::A], mLocalVecDistdRPhi[Side::A], mGrid3D[Side::A], Side::A}, {mLocalVecDistdR[Side::C], mLocalVecDistdZ[Side::C], mLocalVecDistdRPhi[Side::C], mGrid3D[Side::C], Side::C}}; 
  NumericalFields<DataT> mInterpolatorEField[FNSIDES]{{mElectricFieldEr[Side::A], mElectricFieldEz[Side::A], mElectricFieldEphi[Side::A], mGrid3D[Side::A], Side::A}, {mElectricFieldEr[Side::C], mElectricFieldEz[Side::C], mElectricFieldEphi[Side::C], mGrid3D[Side::C], Side::C}};        
  AnalyticalDistCorr<DataT> mAnaDistCorr;                                                                                                                                                                                                                                                     
  bool mUseAnaDistCorr{false};                                                                                                                                                                                                                                                                
  BField mBField{};                                                                                                                                                                                                                                                                           
  SCMetaData mMeta{};                                                                                                                                                                                                                                                                         
 
  bool isCloseToZero(const DataT valA, const DataT valB) const { return std::abs(valA + valB) < static_cast<DataT>(0.01); }
 
  static int getSign(const Side side) { return side == Side::C ? -1 : 1; }
 
  DataT getInvSpacingZ(const Side side) const { return mGrid3D[side].getInvSpacingZ(); }
 
  DataT getInvSpacingR(const Side side) const { return mGrid3D[side].getInvSpacingR(); }
 
  DataT getInvSpacingPhi(const Side side) const { return mGrid3D[side].getInvSpacingPhi(); }
 
  DataT getRMinSim(const Side side) const { return getRMin(side) - 4 * getGridSpacingR(side); }
 
  DataT getRMaxSim(const Side side) const { return getRMax(side) + 2 * getGridSpacingR(side); }
 
  std::string getSideName(const Side side) const { return side == Side::A ? "A" : "C"; }
 
  template <typename Fields = AnalyticalFields<DataT>>
  void calcDistCorr(const DataT p1r, const DataT p1phi, const DataT p1z, const DataT p2z, DataT& ddR, DataT& ddPhi, DataT& ddZ, const Fields& formulaStruct, const bool localDistCorr, const Side side) const;
 
  void langevinCylindricalE(DataT& ddR, DataT& ddPhi, DataT& ddZ, const DataT radius, const DataT localIntErOverEz, const DataT localIntEPhiOverEz, const DataT localIntDeltaEz) const;
 
  void langevinCylindricalB(DataT& ddR, DataT& ddPhi, const DataT radius, const DataT localIntBrOverBz, const DataT localIntBPhiOverBz) const;
 
  template <typename Fields = AnalyticalFields<DataT>>
  void integrateEFieldsRoot(const DataT p1r, const DataT p1phi, const DataT p1z, const DataT p2z, DataT& localIntErOverEz, DataT& localIntEPhiOverEz, DataT& localIntDeltaEz, const Fields& formulaStruct, const DataT ezField, const Side side) const;
 
  template <typename Fields = AnalyticalFields<DataT>>
  void integrateEFieldsTrapezoidal(const DataT p1r, const DataT p1phi, const DataT p1z, const DataT p2z, DataT& localIntErOverEz, DataT& localIntEPhiOverEz, DataT& localIntDeltaEz, const Fields& formulaStruct, const DataT ezField, const Side side) const;
 
  template <typename Fields = AnalyticalFields<DataT>>
  void integrateEFieldsSimpson(const DataT p1r, const DataT p1phi, const DataT p1z, const DataT p2z, DataT& localIntErOverEz, DataT& localIntEPhiOverEz, DataT& localIntDeltaEz, const Fields& formulaStruct, const DataT ezField, const Side side) const;
 
  template <typename Fields = AnalyticalFields<DataT>>
  void integrateEFieldsSimpsonIterative(const DataT p1r, const DataT p2r, const DataT p1phi, const DataT p2phi, const DataT p1z, const DataT p2z, DataT& localIntErOverEz, DataT& localIntEPhiOverEz, DataT& localIntDeltaEz, const Fields& formulaStruct, const DataT ezField, const Side side) const;
 
  void processGlobalDistCorr(const DataT radius, const DataT phi, const DataT z0Tmp, const DataT z1Tmp, DataT& ddR, DataT& ddPhi, DataT& ddZ, const AnalyticalFields<DataT>& formulaStruct) const { calcDistCorr(radius, phi, z0Tmp, z1Tmp, ddR, ddPhi, ddZ, formulaStruct, false, formulaStruct.getSide()); }
 
  void processGlobalDistCorr(const DataT radius, const DataT phi, const DataT z0Tmp, const DataT z1Tmp, DataT& ddR, DataT& ddPhi, DataT& ddZ, const NumericalFields<DataT>& formulaStruct) const { calcDistCorr(radius, phi, z0Tmp, z1Tmp, ddR, ddPhi, ddZ, formulaStruct, false, formulaStruct.getSide()); }
 
  void processGlobalDistCorr(const DataT radius, const DataT phi, const DataT z0Tmp, [[maybe_unused]] const DataT z1Tmp, DataT& ddR, DataT& ddPhi, DataT& ddZ, const DistCorrInterpolator<DataT>& localDistCorr) const
  {
    ddR = localDistCorr.evaldR(z0Tmp, radius, phi);
    ddZ = localDistCorr.evaldZ(z0Tmp, radius, phi);
    ddPhi = localDistCorr.evaldRPhi(z0Tmp, radius, phi) / radius;
  }
 
  void dumpElectronTracksToTree(const std::vector<std::pair<std::vector<o2::math_utils::Point3D<float>>, std::array<DataT, 3>>>& electronTracks, const int nSamplingPoints, const char* outFile) const;
 
  size_t getNearestPhiVertex(const DataT phi, const Side side) const { return std::round(phi / getGridSpacingPhi(side)); }
 
  size_t getNearestRVertex(const DataT r, const Side side) const { return std::round((r - getRMin(side)) / getGridSpacingR(side) + 0.5); }
 
  size_t getPhiBinsGapFrame(const Side side) const;
 
  void setPotentialBoundaryGEMFrameAlongPhi(const std::function<DataT(DataT)>& potentialFunc, const GEMstack stack, const bool bottom, const Side side, const bool outerFrame = false);
 
  void getDistortionsCorrectionsAnalytical(const DataT x, const DataT y, const DataT z, const Side side, DataT& distX, DataT& distY, DataT& distZ, const bool dist) const;
 
  void setBFields(o2::parameters::GRPMagField& magField);
 
  void initContainer(DataContainer& data, const bool initMem = true);
 
  void initAllBuffers();
 
  void setBoundaryFromIndices(const std::function<DataT(DataT)>& potentialFunc, const std::vector<size_t>& indices, const Side side);
 
  std::vector<size_t> getPotentialBoundaryGEMFrameAlongRIndices(const Side side) const;
 
  std::vector<size_t> getPotentialBoundaryGEMFrameAlongPhiIndices(const GEMstack stack, const bool bottom, const Side side, const bool outerFrame, const bool noGap = false) const;
 
  void setROCMisalignment(int stackType, int misalignmentType, int sector, const float potMin, const float potMax);
  void fillROCMisalignment(const std::vector<size_t>& indicesTop, const std::vector<size_t>& indicesBottom, int sector, int misalignmentType, const std::pair<float, float>& deltaPotPar);
 
  void initRodAlignmentVoltages(const MisalignmentType misalignmentType, const FCType fcType, const int sector, const Side side, const float deltaPot);
 
  void calcGlobalDistCorrIterative(const DistCorrInterpolator<DataT>& globCorr, const int maxIter, const DataT approachZ, const DataT approachR, const DataT approachPhi, const DataT diffCorr, const SpaceCharge<DataT>* scSCale, float scale, const Type type);
  void calcGlobalDistCorrIterativeLinearCartesian(const DistCorrInterpolator<DataT>& globCorr, const int maxIter, const DataT approachX, const DataT approachY, const DataT approachZ, const DataT diffCorr, const SpaceCharge<DataT>* scSCale, float scale, const Type type);
 
  ClassDefNV(SpaceCharge, 6);
};
 
} // namespace tpc
} // namespace o2
 
#endif