AliceO2Group · davidrohr · Mar 18, 2022 · Mar 15, 2022
@@ -70,6 +70,15 @@ struct TPCParameters {
   inline static DataT vg1t{-3260};                              ///< GEM 1 Top voltage. (setting with reduced ET1,2,4 = 3.5kV/cm)
 };
 
+template <typename DataT = double>
+struct GEMFrameParameters {
+  static constexpr DataT WIDTHFRAME{1};                             ///< width of the frame 1 cm
+  static constexpr DataT LENGTHFRAMEIROCBOTTOM{29.195};             ///< length of the GEM frame on the bottom side of the IROC
+  static constexpr DataT LENGTHFRAMEOROC3TOP{87.048};               ///< length of the GEM frame on the top side of the OROC3
+  static constexpr DataT POSBOTTOM[]{83.65, 133.5, 169.75, 207.85}; ///< local x position of the GEM frame on the bottom side per stack
+  static constexpr DataT POSTOP[]{133.3, 169.75, 207.85, 247.7};    ///< local x position of the GEM frame on the top side per stack
+};
+
 template <typename DataT = double>
 struct GridProperties {
   static constexpr DataT RMIN{TPCParameters<DataT>::IFCRADIUS};                  ///< min radius

@@ -133,12 +133,66 @@ class SpaceCharge
   /// \param formulaStruct struct containing a method to evaluate the potential
   void setPotentialFromFormula(const AnalyticalFields<DataT>& formulaStruct);
 
+  /// setting the boundary potential of the GEM stack along the radius
+  /// \param potentialFunc potential funtion as a function of the radius
+  /// \side Side of the TPC
+  void setPotentialBoundaryGEMFrameAlongR(const std::function<DataT(DataT)>& potentialFunc, const Side side);
+
+  /// setting the boundary potential of the IROC on the bottom along phi
+  /// \param potentialFunc potential funtion as a function of global phi
+  /// \side Side of the TPC
+  void setPotentialBoundaryGEMFrameIROCBottomAlongPhi(const std::function<DataT(DataT)>& potentialFunc, const Side side) { setPotentialBoundaryGEMFrameAlongPhi(potentialFunc, GEMstack::IROCgem, true, side); }
+
+  /// setting the boundary potential of the IROC on the top along phi
+  /// \param potentialFunc potential funtion as a function of global phi
+  /// \side Side of the TPC
+  void setPotentialBoundaryGEMFrameIROCTopAlongPhi(const std::function<DataT(DataT)>& potentialFunc, const Side side) { setPotentialBoundaryGEMFrameAlongPhi(potentialFunc, GEMstack::IROCgem, false, side); }
+
+  /// setting the boundary potential of the OROC1 on the bottom along phi
+  /// \param potentialFunc potential funtion as a function of global phi
+  /// \side Side of the TPC
+  void setPotentialBoundaryGEMFrameOROC1BottomAlongPhi(const std::function<DataT(DataT)>& potentialFunc, const Side side) { setPotentialBoundaryGEMFrameAlongPhi(potentialFunc, GEMstack::OROC1gem, true, side); }
+
+  /// setting the boundary potential of the OROC1 on the top along phi
+  /// \param potentialFunc potential funtion as a function of global phi
+  /// \side Side of the TPC
+  void setPotentialBoundaryGEMFrameOROC1TopAlongPhi(const std::function<DataT(DataT)>& potentialFunc, const Side side) { setPotentialBoundaryGEMFrameAlongPhi(potentialFunc, GEMstack::OROC1gem, false, side); }
+
+  /// setting the boundary potential of the OROC2 on the bottom along phi
+  /// \param potentialFunc potential funtion as a function of global phi
+  /// \side Side of the TPC
+  void setPotentialBoundaryGEMFrameOROC2BottomAlongPhi(const std::function<DataT(DataT)>& potentialFunc, const Side side) { setPotentialBoundaryGEMFrameAlongPhi(potentialFunc, GEMstack::OROC2gem, true, side); }
+
+  /// setting the boundary potential of the OROC2 on the top along phi
+  /// \param potentialFunc potential funtion as a function of global phi
+  /// \side Side of the TPC
+  void setPotentialBoundaryGEMFrameOROC2TopAlongPhi(const std::function<DataT(DataT)>& potentialFunc, const Side side) { setPotentialBoundaryGEMFrameAlongPhi(potentialFunc, GEMstack::OROC2gem, false, side); }
+
+  /// setting the boundary potential of the OROC3 on the bottom along phi
+  /// \param potentialFunc potential funtion as a function of global phi
+  /// \side Side of the TPC
+  void setPotentialBoundaryGEMFrameOROC3BottomAlongPhi(const std::function<DataT(DataT)>& potentialFunc, const Side side) { setPotentialBoundaryGEMFrameAlongPhi(potentialFunc, GEMstack::OROC3gem, true, side); }
+
+  /// setting the boundary potential of the OROC3 on the top along phi
+  /// \param potentialFunc potential funtion as a function of global phi
+  /// \side Side of the TPC
+  void setPotentialBoundaryGEMFrameOROC3TopAlongPhi(const std::function<DataT(DataT)>& potentialFunc, const Side side) { setPotentialBoundaryGEMFrameAlongPhi(potentialFunc, GEMstack::OROC3gem, false, side); }
+
+  /// setting the boundary potential for the inner TPC radius along r
+  /// \param potentialFunc potential funtion as a function of z
+  /// \side Side of the TPC
+  void setPotentialBoundaryInnerRadius(const std::function<DataT(DataT)>& potentialFunc, const Side side);
+
+  /// setting the boundary potential for the outer TPC radius along r
+  /// \param potentialFunc potential funtion as a function of z
+  /// \side Side of the TPC
+  void setPotentialBoundaryOuterRadius(const std::function<DataT(DataT)>& potentialFunc, const Side side);
+
   /// step 1: use the O2TPCPoissonSolver class to numerically calculate the potential with set space charge density and boundary conditions from potential
   /// \param side side of the TPC
-  /// \param maxIteration maximum number of iterations used in the poisson solver
   /// \param stoppingConvergence stopping criterion used in the poisson solver
   /// \param symmetry use symmetry or not in the poisson solver
-  void poissonSolver(const Side side, const int maxIteration = 300, const DataT stoppingConvergence = 1e-6, const int symmetry = 0);
+  void poissonSolver(const Side side, const DataT stoppingConvergence = 1e-6, const int symmetry = 0);
 
   /// step 2: calculate numerically the electric field from the potential
   /// \param side side of the TPC
@@ -1000,6 +1054,18 @@ class SpaceCharge
 
   /// dump the created electron tracks with calculateElectronDriftPath function to a tree
   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;
+
+  /// \return returns nearest phi vertex for given phi position
+  size_t getNearestPhiVertex(const DataT phi, const Side side) const { return std::round(phi / getGridSpacingPhi(side)); }
+
+  /// \return returns nearest r vertex for given radius position
+  size_t getNearestRVertex(const DataT r, const Side side) const { return std::round((r - getRMin(side)) / getGridSpacingR(side)); }
+
+  /// \return returns number of bins in phi direction for the gap between sectors and for the GEM frame
+  std::pair<size_t, size_t> getPhiBinsGapFrame(const Side side) const;
+
+  /// \return setting the boundary potential for given GEM stack
+  void setPotentialBoundaryGEMFrameAlongPhi(const std::function<DataT(DataT)>& potentialFunc, const GEMstack stack, const bool bottom, const Side side);
 };
 
 ///

@@ -179,6 +179,135 @@ void SpaceCharge<DataT>::setChargeDensityFromFormula(const AnalyticalFields<Data
   mIsChargeSet[side] = true;
 }
 
+template <typename DataT>
+std::pair<size_t, size_t> SpaceCharge<DataT>::getPhiBinsGapFrame(const Side side) const
+{
+  const auto& regInf = Mapper::instance().getPadRegionInfo(0);
+  const float localYEdgeIROC = regInf.getPadsInRowRegion(0) / 2 * regInf.getPadWidth();
+  const auto globalPosGap = Mapper::LocalToGlobal(LocalPosition2D(regInf.getRadiusFirstRow(), -(localYEdgeIROC + GEMFrameParameters<DataT>::WIDTHFRAME)), Sector(0));
+  const auto phiGap = std::atan(globalPosGap.Y() / globalPosGap.X());
+
+  auto nBinsPhiGap = getNearestPhiVertex(phiGap, side);
+  if (nBinsPhiGap == 0) {
+    nBinsPhiGap = 1;
+  }
+
+  const auto globalPosEdgeIROC = Mapper::LocalToGlobal(LocalPosition2D(regInf.getRadiusFirstRow(), -localYEdgeIROC), Sector(0));
+  const auto phiEdge = std::atan(globalPosEdgeIROC.Y() / globalPosEdgeIROC.X());
+  auto nBinsPhiGEMFrame = getNearestPhiVertex(phiEdge - phiGap, side);
+  if (nBinsPhiGEMFrame == 0) {
+    nBinsPhiGEMFrame = 1;
+  }
+
+  return std::pair<size_t, size_t>(nBinsPhiGap, nBinsPhiGEMFrame);
+}
+
+template <typename DataT>
+void SpaceCharge<DataT>::setPotentialBoundaryGEMFrameAlongR(const std::function<DataT(DataT)>& potentialFunc, const Side side)
+{
+  const auto radiusStart = std::sqrt(std::pow(GEMFrameParameters<DataT>::LENGTHFRAMEIROCBOTTOM / 2, 2) + std::pow(GEMFrameParameters<DataT>::POSBOTTOM[0], 2));
+  const auto rStart = getNearestRVertex(radiusStart, side);
+  const auto radiusEnd = std::sqrt(std::pow(GEMFrameParameters<DataT>::LENGTHFRAMEOROC3TOP / 2, 2) + std::pow(GEMFrameParameters<DataT>::POSTOP[GEMstack::OROC3gem], 2));
+  const auto rEnd = getNearestRVertex(radiusEnd, side);
+
+  const auto nBinsPhi = getPhiBinsGapFrame(side);
+  const int verticesPerSector = mParamGrid.NPhiVertices / SECTORSPERSIDE;
+  for (int sector = 0; sector < SECTORSPERSIDE; ++sector) {
+    for (size_t iPhiTmp = 0; iPhiTmp < nBinsPhi.second; ++iPhiTmp) {
+      const size_t iPhiLeft = sector * verticesPerSector + nBinsPhi.first + iPhiTmp;
+      const size_t iPhiRight = (sector + 1) * verticesPerSector - nBinsPhi.first - iPhiTmp;
+      for (size_t iR = rStart; iR < rEnd; ++iR) {
+        const DataT radius = getRVertex(iR, side);
+        const size_t iZ = mParamGrid.NZVertices - 1;
+        mPotential[side](iZ, iR, iPhiLeft) = potentialFunc(radius);
+        mPotential[side](iZ, iR, iPhiRight) = potentialFunc(radius);
+      }
+    }
+  }
+}
+
+template <typename DataT>
+void SpaceCharge<DataT>::setPotentialBoundaryGEMFrameAlongPhi(const std::function<DataT(DataT)>& potentialFunc, const GEMstack stack, const bool bottom, const Side side)
+{
+  int region = 0;
+  if (bottom) {
+    if (stack == GEMstack::IROCgem) {
+      region = 0;
+    } else if (stack == GEMstack::OROC1gem) {
+      region = 4;
+    } else if (stack == GEMstack::OROC2gem) {
+      region = 6;
+    } else if (stack == GEMstack::OROC3gem) {
+      region = 8;
+    }
+  } else {
+    if (stack == GEMstack::IROCgem) {
+      region = 3;
+    } else if (stack == GEMstack::OROC1gem) {
+      region = 5;
+    } else if (stack == GEMstack::OROC2gem) {
+      region = 7;
+    } else if (stack == GEMstack::OROC3gem) {
+      region = 9;
+    }
+  }
+
+  const auto& regInf = Mapper::instance().getPadRegionInfo(region);
+  const auto radiusFirstRow = regInf.getRadiusFirstRow();
+  const DataT radiusStart = bottom ? GEMFrameParameters<DataT>::POSBOTTOM[stack] : radiusFirstRow + regInf.getPadHeight() * regInf.getNumberOfPadRows();
+  const auto radiusMax = bottom ? radiusFirstRow : GEMFrameParameters<DataT>::POSTOP[stack];
+  auto nVerticesR = std::round((radiusMax - radiusStart) / getGridSpacingR(side));
+  if (nVerticesR == 0) {
+    nVerticesR = 1;
+  }
+
+  const int verticesPerSector = mParamGrid.NPhiVertices / SECTORSPERSIDE;
+  const auto nBinsPhi = getPhiBinsGapFrame(side);
+  for (int sector = 0; sector < SECTORSPERSIDE; ++sector) {
+    const auto offsetPhi = sector * verticesPerSector + verticesPerSector / 2;
+    for (size_t iPhiLocal = 0; iPhiLocal <= verticesPerSector / 2 - nBinsPhi.first; ++iPhiLocal) {
+      const auto iPhiLeft = offsetPhi + iPhiLocal;
+      const auto iPhiRight = offsetPhi - iPhiLocal;
+      const DataT phiLeft = getPhiVertex(iPhiLeft, side);
+      const DataT phiRight = getPhiVertex(iPhiRight, side);
+      const DataT localphi = getPhiVertex(iPhiLocal, side);
+      const DataT radiusBottom = radiusStart / std::cos(localphi);
+      const auto rStart = getNearestRVertex(radiusBottom, side);
+      for (size_t iR = rStart; iR < rStart + nVerticesR; ++iR) {
+        const size_t iZ = mParamGrid.NZVertices - 1;
+        mPotential[side](iZ, iR, iPhiLeft) = potentialFunc(phiLeft);
+        mPotential[side](iZ, iR, iPhiRight) = potentialFunc(phiRight);
+      }
+    }
+  }
+}
+
+template <typename DataT>
+void SpaceCharge<DataT>::setPotentialBoundaryInnerRadius(const std::function<DataT(DataT)>& potentialFunc, const Side side)
+{
+  for (size_t iZ = 0; iZ < mParamGrid.NZVertices; ++iZ) {
+    const DataT z = getZVertex(iZ, side);
+    const auto pot = potentialFunc(z);
+    for (size_t iPhi = 0; iPhi < mParamGrid.NPhiVertices; ++iPhi) {
+      const size_t iR = 0;
+      mPotential[side](iZ, iR, iPhi) = pot;
+    }
+  }
+}
+
+template <typename DataT>
+void SpaceCharge<DataT>::setPotentialBoundaryOuterRadius(const std::function<DataT(DataT)>& potentialFunc, const Side side)
+{
+  for (size_t iZ = 0; iZ < mParamGrid.NZVertices; ++iZ) {
+    const DataT z = getZVertex(iZ, side);
+    const auto pot = potentialFunc(z);
+    for (size_t iPhi = 0; iPhi < mParamGrid.NPhiVertices; ++iPhi) {
+      const size_t iR = mParamGrid.NRVertices - 1;
+      mPotential[side](iZ, iR, iPhi) = pot;
+    }
+  }
+}
+
 template <typename DataT>
 void SpaceCharge<DataT>::setPotentialFromFormula(const AnalyticalFields<DataT>& formulaStruct)
 {
@@ -241,7 +370,7 @@ void SpaceCharge<DataT>::setPotentialBoundaryFromFormula(const AnalyticalFields<
 }
 
 template <typename DataT>
-void SpaceCharge<DataT>::poissonSolver(const Side side, const int maxIteration, const DataT stoppingConvergence, const int symmetry)
+void SpaceCharge<DataT>::poissonSolver(const Side side, const DataT stoppingConvergence, const int symmetry)
 {
   PoissonSolver<DataT>::setConvergenceError(stoppingConvergence);
   PoissonSolver<DataT> poissonSolver(mGrid3D[0]);