[SC] Adding possibility to use float for calculations (#6870)

* [SC] Adding possibility to use float for calculations

- setting default values for omega tau
- adding function to calculate the drift path for an electron
- adding simple debug function
- poisson solver: cleaning + using more than 1000 phi bins

* fixing function name getChargeCyl -> getDensityCyl

* fixing typo

* fixing error during compiling of macro

Author: matthias-kleiner

Detectors/TPC/reconstruction/macro/createTPCSpaceChargeCorrection.C

@@ -222,7 +222,7 @@ void debugInterpolation(utils::TreeStreamRedirector& pcstream,
           float gxyzf[3] = {gx, gy, gz};
           float pointCyl[3] = {gz, r, phi};
           double efield[3] = {0.0, 0.0, 0.0};
-          double charge = spaceCharge->getChargeCyl(pointCyl[0], pointCyl[1], pointCyl[2], side);
+          double charge = spaceCharge->getDensityCyl(pointCyl[0], pointCyl[1], pointCyl[2], side);
           double potential = spaceCharge->getPotentialCyl(pointCyl[0], pointCyl[1], pointCyl[2], side);
           spaceCharge->getElectricFieldsCyl(pointCyl[0], pointCyl[1], pointCyl[2], side, efield[0], efield[1], efield[2]);
           spaceCharge->getLocalDistortionsCyl(pointCyl[0], pointCyl[1], pointCyl[2], side, ldD[0], ldD[1], ldD[2]);
@@ -379,7 +379,7 @@ void debugGridpoints(utils::TreeStreamRedirector& pcstream, const o2::gpu::TPCFa
           double efield[3] = {0.0, 0.0, 0.0};
           double distLocal[3] = {0.0, 0.0, 0.0};
           double dist[3] = {0.0, 0.0, 0.0};
-          double charge = spaceCharge->getChargeCyl(pointCyl[0], pointCyl[1], pointCyl[2], side);
+          double charge = spaceCharge->getDensityCyl(pointCyl[0], pointCyl[1], pointCyl[2], side);
           double potential = spaceCharge->getPotentialCyl(pointCyl[0], pointCyl[1], pointCyl[2], side);
           spaceCharge->getElectricFieldsCyl(pointCyl[0], pointCyl[1], pointCyl[2], side, efield[0], efield[1], efield[2]);
           spaceCharge->getLocalDistortionsCyl(pointCyl[0], pointCyl[1], pointCyl[2], side, distLocal[0], distLocal[1], distLocal[2]);

Detectors/TPC/spacecharge/include/TPCSpaceCharge/RegularGrid3D.h

@@ -137,19 +137,19 @@ struct RegularGrid3D {
 
  private:
   inline static auto& mParamGrid = ParameterSpaceCharge::Instance();
-  static constexpr unsigned int FDIM = 3;                                                                                                                                                   ///< dimensions of the grid (only 3 supported)
-  static constexpr unsigned int FZ = 0;                                                                                                                                                     ///< index for x coordinate
-  static constexpr unsigned int FR = 1;                                                                                                                                                     ///< index for y coordinate
-  static constexpr unsigned int FPHI = 2;                                                                                                                                                   ///< index for z coordinate
-  const Vector<DataT, FDIM> mMin{};                                                                                                                                                         ///< min vertices positions of the grid
-  const Vector<DataT, FDIM> mMax{};                                                                                                                                                         ///< max vertices positions of the grid
-  const Vector<DataT, FDIM> mSpacing{};                                                                                                                                                     ///<  spacing of the grid
-  const Vector<DataT, FDIM> mInvSpacing{};                                                                                                                                                  ///< inverse spacing of grid
-  const Vector<DataT, FDIM> sMaxIndex{{static_cast<double>(mParamGrid.NZVertices - 1.), static_cast<double>(mParamGrid.NRVertices - 1), static_cast<double>(mParamGrid.NPhiVertices - 1)}}; ///< max index which is on the grid in all dimensions
-  const Vector<int, FDIM> sNdim{{static_cast<int>(mParamGrid.NZVertices), static_cast<int>(mParamGrid.NRVertices), static_cast<int>(mParamGrid.NPhiVertices)}};                             ///< number of vertices for each dimension
-  std::vector<DataT> mZVertices{};                                                                                                                                                          ///< positions of vertices in x direction
-  std::vector<DataT> mRVertices{};                                                                                                                                                          ///< positions of vertices in y direction
-  std::vector<DataT> mPhiVertices{};                                                                                                                                                        ///< positions of vertices in z direction
+  static constexpr unsigned int FDIM = 3;                                                                                                                                                ///< dimensions of the grid (only 3 supported)
+  static constexpr unsigned int FZ = 0;                                                                                                                                                  ///< index for x coordinate
+  static constexpr unsigned int FR = 1;                                                                                                                                                  ///< index for y coordinate
+  static constexpr unsigned int FPHI = 2;                                                                                                                                                ///< index for z coordinate
+  const Vector<DataT, FDIM> mMin{};                                                                                                                                                      ///< min vertices positions of the grid
+  const Vector<DataT, FDIM> mMax{};                                                                                                                                                      ///< max vertices positions of the grid
+  const Vector<DataT, FDIM> mSpacing{};                                                                                                                                                  ///<  spacing of the grid
+  const Vector<DataT, FDIM> mInvSpacing{};                                                                                                                                               ///< inverse spacing of grid
+  const Vector<DataT, FDIM> sMaxIndex{{static_cast<DataT>(mParamGrid.NZVertices - 1.), static_cast<DataT>(mParamGrid.NRVertices - 1), static_cast<DataT>(mParamGrid.NPhiVertices - 1)}}; ///< max index which is on the grid in all dimensions
+  const Vector<int, FDIM> sNdim{{static_cast<int>(mParamGrid.NZVertices), static_cast<int>(mParamGrid.NRVertices), static_cast<int>(mParamGrid.NPhiVertices)}};                          ///< number of vertices for each dimension
+  std::vector<DataT> mZVertices{};                                                                                                                                                       ///< positions of vertices in x direction
+  std::vector<DataT> mRVertices{};                                                                                                                                                       ///< positions of vertices in y direction
+  std::vector<DataT> mPhiVertices{};                                                                                                                                                     ///< positions of vertices in z direction
 
   void initLists();
 

Detectors/TPC/spacecharge/include/TPCSpaceCharge/SpaceCharge.h

@@ -28,6 +28,7 @@
 
 class TH3;
 class TH3D;
+class TH3F;
 
 namespace o2
 {
@@ -50,11 +51,15 @@ class SpaceCharge
   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:
   /// default constructor
   /// grid granularity has to set before constructing an object using the static function setGrid(nZVertices, nRVertices, nPhiVertices)!
-  SpaceCharge() = default;
+  /// \param omegaTau \omega \tau value
+  /// \param t1 value for t1
+  /// \param t2 value for t2
+  SpaceCharge(const DataT omegaTau = -0.35f, const DataT t1 = 1, const DataT t2 = 1) { setOmegaTauT1T2(omegaTau, t1, t2); };
 
   /// \param nZVertices number of vertices of the grid in z direction
   /// \param nRVertices number of vertices of the grid in z direction
@@ -182,7 +187,7 @@ class SpaceCharge
   /// \param z global z coordinate
   /// \param r global r coordinate
   /// \param phi global phi coordinate
-  DataT getChargeCyl(const DataT z, const DataT r, const DataT phi, const Side side) const;
+  DataT getDensityCyl(const DataT z, const DataT r, const DataT phi, const Side side) const;
 
   /// get the potential for given coordinate
   /// \param z global z coordinate
@@ -538,7 +543,7 @@ class SpaceCharge
   /// \param omegaTau \omega \tau value
   /// \param t1 value for t1 see: ???
   /// \param t2 value for t2 see: ???
-  void setOmegaTauT1T2(const DataT omegaTau, const DataT t1, const DataT t2)
+  void setOmegaTauT1T2(const DataT omegaTau = -0.35f, const DataT t1 = 1, const DataT t2 = 1)
   {
     const DataT wt0 = t2 * omegaTau;
     mC0 = 1 / (1 + wt0 * wt0);
@@ -567,7 +572,7 @@ class SpaceCharge
   static void setNThreads(const int nThreads)
   {
     sNThreads = nThreads;
-    o2::tpc::TriCubicInterpolator<DataT>::setNThreads(nThreads);
+    TriCubic::setNThreads(nThreads);
   }
 
   /// set which kind of numerical integration is used for calcution of the integrals int Er/Ez dz, int Ephi/Ez dz, int Ez dz
@@ -597,6 +602,14 @@ class SpaceCharge
   /// \side side of the TPC
   void dumpToFile(TFile& outf, const Side side) const;
 
+  /// dump sc density, potential, electric fields, global distortions/corrections to tree
+  /// \param outFileName name of the output file
+  /// \side side of the TPC
+  /// \param nZPoints number of vertices of the output in z
+  /// \param nRPoints number of vertices of the output in r
+  /// \param nPhiPoints number of vertices of the output in phi
+  void dumpToTree(const char* outFileName, const Side side, const int nZPoints = 50, const int nRPoints = 150, const int nPhiPoints = 180) const;
+
   /// write electric fields to root file
   /// \param outf output file where the electrical fields will be written to
   /// \side side of the TPC
@@ -716,7 +729,13 @@ class SpaceCharge
   /// \param nBinsZNew number of z bins of rebinned histogram
   /// \param nBinsRNew number of r bins of rebinned histogram
   /// \param nBinsPhiNew number of phi bins of rebinned histogram
-  static TH3D rebinDensityHisto(const TH3& hOrig, const unsigned short nBinsZNew, const unsigned short nBinsRNew, const unsigned short nBinsPhiNew);
+  static TH3DataT rebinDensityHisto(const TH3& hOrig, const unsigned short nBinsZNew, const unsigned short nBinsRNew, const unsigned short nBinsPhiNew);
+
+  /// function to calculate the drift paths of the electron whose starting position is delivered. Electric fields must be set!
+  /// \param elePos global position of the start position of the electron
+  /// \param nSamplingPoints number of output points of the electron drift path
+  /// \return returns the input electron a vector of 3D-points describing the drift path of the electron
+  std::vector<GlobalPosition3D> calculateElectronDriftPath(const GlobalPosition3D& elePos, const int nSamplingPoints) const;
 
  private:
   inline static auto& mParamGrid = ParameterSpaceCharge::Instance(); ///< parameters of the grid on which the calculations are performed

Detectors/TPC/spacecharge/include/TPCSpaceCharge/Vector.h

@@ -17,7 +17,8 @@
 #ifndef ALICEO2_TPC_VECTOR_H_
 #define ALICEO2_TPC_VECTOR_H_
 
-#include <Vc/Vc>
+#include <Vc/vector>
+#include <Vc/Memory>
 
 namespace o2
 {

Detectors/TPC/spacecharge/macro/calculateDistortionsCorrections.C

@@ -6,11 +6,15 @@
 #include "CommonUtils/TreeStreamRedirector.h"
 
 template <typename DataT>
-void calculateDistortionsAnalytical(const int globalEFieldTypeAna = 1, const int globalDistTypeAna = 1, const int eFieldTypeAna = 1, const int usePoissonSolverAna = 1, const int nSteps = 1, const int simpsonIterations = 3, const int nThreads = -1);
+void calculateDistortionsAnalytical(const int sides = 0, const int globalEFieldTypeAna = 1, const int globalDistTypeAna = 1, const int eFieldTypeAna = 1, const int usePoissonSolverAna = 1, const int nSteps = 1, const int simpsonIterations = 3, const int nThreads = -1);
 
 template <typename DataT>
 void calculateDistortionsFromHist(const char* path, const char* histoName, const int sides, const int globalEFieldType = 1, const int globalDistType = 1, const int nSteps = 1, const int simpsonIterations = 3, const int nThreads = -1);
 
+int getSideStart(const int sides);
+int getSideEnd(const int sides);
+
+/// \param sides set which sides will be simulated. sides=0: A- and C-Side, sides=1: A-Side only, sides=2: C-Side only
 /// \param nZVertices number of vertices of the grid in z direction
 /// \param nRVertices number of vertices of the grid in z direction
 /// \param nPhiVertices number of vertices of the grid in z direction
@@ -21,10 +25,11 @@ void calculateDistortionsFromHist(const char* path, const char* histoName, const
 /// \param nSteps number of which are used for calculation of distortions/corrections per z-bin
 /// \param simpsonIterations number of iterations used in the simpson intergration
 /// \param nThreads number of threads which are used (if the value is -1 all threads should be used)
-void calcDistAna(const unsigned short nZVertices = 129, const unsigned short nRVertices = 129, const unsigned short nPhiVertices = 180, const int globalEFieldTypeAna = 1, const int globalDistTypeAna = 1, const int eFieldTypeAna = 1, const int usePoissonSolverAna = 1, const int nSteps = 1, const int simpsonIterations = 3, const int nThreads = -1)
+template <typename DataT>
+void calcDistAna(const int sides = 0, const unsigned short nZVertices = 129, const unsigned short nRVertices = 129, const unsigned short nPhiVertices = 180, const int globalEFieldTypeAna = 1, const int globalDistTypeAna = 1, const int eFieldTypeAna = 1, const int usePoissonSolverAna = 1, const int nSteps = 1, const int simpsonIterations = 3, const int nThreads = -1)
 {
-  o2::tpc::SpaceCharge<double>::setGrid(nZVertices, nRVertices, nPhiVertices);
-  calculateDistortionsAnalytical<double>(globalEFieldTypeAna, globalDistTypeAna, eFieldTypeAna, usePoissonSolverAna, nSteps, simpsonIterations, nThreads);
+  o2::tpc::SpaceCharge<DataT>::setGrid(nZVertices, nRVertices, nPhiVertices);
+  calculateDistortionsAnalytical<DataT>(sides, globalEFieldTypeAna, globalDistTypeAna, eFieldTypeAna, usePoissonSolverAna, nSteps, simpsonIterations, nThreads);
 }
 
 /// \param path path to the root file containing the 3D density histogram
@@ -271,7 +276,7 @@ void writeToTree(o2::tpc::SpaceCharge<DataT>& spaceCharge3D, o2::utils::TreeStre
 /// \param nSteps number of which are used for calculation of distortions/corrections per z-bin
 /// \param simpsonIterations number of iterations used in the simpson intergration
 template <typename DataT = double>
-void calculateDistortionsAnalytical(const int globalEFieldTypeAna, const int globalDistTypeAna, const int eFieldTypeAna, const int usePoissonSolverAna, const int nSteps, const int simpsonIterations, const int nThreads)
+void calculateDistortionsAnalytical(const int sides, const int globalEFieldTypeAna, const int globalDistTypeAna, const int eFieldTypeAna, const int usePoissonSolverAna, const int nSteps, const int simpsonIterations, const int nThreads)
 {
   const auto integrationStrategy = o2::tpc::SpaceCharge<DataT>::IntegrationStrategy::SimpsonIterative;
   o2::tpc::SpaceCharge<DataT> spaceCharge3D;
@@ -286,7 +291,7 @@ void calculateDistortionsAnalytical(const int globalEFieldTypeAna, const int glo
   // write to root file
   o2::utils::TreeStreamRedirector pcstream(TString::Format("distortions_ana_nR%hu_nZ%hu_nPhi%hu_SimpsonsIter%i.root", spaceCharge3D.getNRVertices(), spaceCharge3D.getNZVertices(), spaceCharge3D.getNPhiVertices(), simpsonIterations).Data(), "RECREATE");
 
-  for (int iside = 0; iside < 2; ++iside) {
+  for (int iside = getSideStart(sides); iside < getSideEnd(sides); ++iside) {
     std::cout << "side: " << iside << std::endl;
     o2::tpc::Side side = (iside == 0) ? o2::tpc::Side::A : o2::tpc::Side::C;
     o2::tpc::AnalyticalFields<DataT> anaFields(side);
@@ -327,16 +332,7 @@ void calculateDistortionsFromHist(const char* path, const char* histoName, const
   fileOutSCDensity.Close();
 
   o2::utils::TreeStreamRedirector pcstream(TString::Format("distortions_real_nR%hu_nZ%hu_nPhi%hu_SimpsonsIter%i.root", spaceCharge3D.getNRVertices(), spaceCharge3D.getNZVertices(), spaceCharge3D.getNPhiVertices(), simpsonIterations).Data(), "RECREATE");
-  int iSideStart = 0;
-  int iSideEnd = 2;
-  if (sides == 1) {
-    // a side only
-    iSideEnd = 1;
-  } else if (sides == 2) {
-    // c side only
-    iSideStart = 1;
-  }
-  for (int iside = iSideStart; iside < iSideEnd; ++iside) {
+  for (int iside = getSideStart(sides); iside < getSideEnd(sides); ++iside) {
     std::cout << "side: " << iside << std::endl;
     o2::tpc::Side side = (iside == 0) ? o2::tpc::Side::A : o2::tpc::Side::C;
     const auto distType = globalDistType == 0 ? SC::GlobalDistType::Standard : SC ::GlobalDistType::Fast;
@@ -364,3 +360,23 @@ void calculateDistortionsFromHist(const char* path, const char* histoName, const
   }
   fOut.Close();
 }
+
+/// helper function to set the loop over the sides for the tpc
+/// \param sides set for which sides the distortions/corrections will be calculated. sides=0: A- and C-Side, sides=1: A-Side only, sides=2: C-Side only
+int getSideStart(const int sides)
+{
+  if (sides == 2) {
+    return 1;
+  }
+  return 0;
+}
+
+/// helper function to set the loop over the sides for the tpc
+/// \param sides set for which sides the distortions/corrections will be calculated. sides=0: A- and C-Side, sides=1: A-Side only, sides=2: C-Side only
+int getSideEnd(const int sides)
+{
+  if (sides == 1) {
+    return 1;
+  }
+  return 2;
+}