TPC SpaceCharge: add scaling of sc density per stack and IFC models

- Adding function to add the boundary potential from other space-charge objects to be able to get the combined model from single models
- Adding function to scale the space-charge density for single stacks. This can be used to emulate the impact on the distortions of IBF variations across chambers
- Adding some functions which can be used with the Python interface
- Adding functions for the inner field cage charging up effect
- Fixing some small bugs

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

@@ -189,6 +189,7 @@ struct DataContainer3D {
   /// operator overload
   DataContainer3D<DataT>& operator*=(const DataT value);
   DataContainer3D<DataT>& operator+=(const DataContainer3D<DataT>& other);
+  DataContainer3D<DataT>& operator*=(const DataContainer3D<DataT>& other);
   DataContainer3D<DataT>& operator-=(const DataContainer3D<DataT>& other);
 
  private:

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

@@ -180,6 +180,14 @@ class SpaceCharge
   /// \param formulaStruct struct containing a method to evaluate the potential
   void setPotentialBoundaryFromFormula(const AnalyticalFields<DataT>& formulaStruct);
 
+  /// adding the boundary potential from other sc object which same number of vertices!
+  /// \param other other SC object which boundary potential witll be added
+  /// \param scaling sclaing factor which used to scale the others boundary potential
+  void addBoundaryPotential(const SpaceCharge<DataT>& other, const Side side, const float scaling = 1);
+
+  /// setting the boundary potential to 0 for z<zMin and z>zMax
+  void resetBoundaryPotentialToZeroInRangeZ(float zMin, float zMax, const Side side);
+
   /// step 0: this function fills the potential using an analytical formula
   /// \param formulaStruct struct containing a method to evaluate the potential
   void setPotentialFromFormula(const AnalyticalFields<DataT>& formulaStruct);
@@ -292,6 +300,9 @@ class SpaceCharge
   /// \param scalingFactor scaling factor for the space-charge density
   void scaleChargeDensitySector(const float scalingFactor, const Sector sector);
 
+  /// scaling the space-charge density for given stack
+  void scaleChargeDensityStack(const float scalingFactor, const Sector sector, const GEMstack stack);
+
   /// add space charge density from other object (this.mDensity = this.mDensity + other.mDensity)
   /// \param otherSC other space-charge object, which charge will be added to current object
   void addChargeDensity(const SpaceCharge<DataT>& otherSC);
@@ -320,6 +331,9 @@ class SpaceCharge
   template <typename Fields = AnalyticalFields<DataT>>
   void calcGlobalCorrections(const Fields& formulaStruct, const int type = 3);
 
+  /// calculate the global corrections using the electric fields (interface for python)
+  void calcGlobalCorrectionsEField(const Side side, const int type = 3) { calcGlobalCorrections(getElectricFieldsInterpolator(side), type); }
+
   /// step 5: calculate global distortions by using the electric field or the local distortions (SLOW)
   /// \param formulaStruct struct containing a method to evaluate the electric field Er, Ez, Ephi or the local distortions
   /// \param maxIterations maximum steps which are are performed to reach the central electrode (in general this is not necessary, but in case of problems this value aborts the calculation)
@@ -372,6 +386,9 @@ class SpaceCharge
   /// \param phi global phi coordinate
   DataT getPotentialCyl(const DataT z, const DataT r, const DataT phi, const Side side) const;
 
+  /// get the potential for list of given coordinate
+  std::vector<float> getPotentialCyl(const std::vector<DataT>& z, const std::vector<DataT>& r, const std::vector<DataT>& phi, const Side side) const;
+
   /// get the electric field for given coordinate
   /// \param z global z coordinate
   /// \param r global r coordinate
@@ -1112,11 +1129,19 @@ class SpaceCharge
   /// \param deltaPot delta potential which will be set at the copper rod
   void setDeltaVoltageRotatedClipOFC(const float deltaPot, const Side side, const int sector) { initRodAlignmentVoltages(MisalignmentType::RotatedClip, FCType::OFC, sector, side, deltaPot); }
 
+  /// IFC charge up: set a linear rising delta potential from the CE to given z
+  /// \param deltaPot maximum value of the delta potential in V
+  /// \param zMaxDeltaPot z position where the maximum delta potential of deltaPot will be set
+  /// \param type functional form of the delta potential: 0=linear, 1= 1/x, 2=flat, 3=linear falling, 4=flat no z dependence
+  /// \param zStart usually at 0 to start the rising of the potential at the IFC
+  void setIFCChargeUpRisingPot(const float deltaPot, const float zMaxDeltaPot, const int type, const float zStart, const float offs, const Side side);
+
   /// IFC charge up: set a linear rising delta potential from the CE to given z position which falls linear down to 0 at the readout
   /// \param deltaPot maximum value of the delta potential in V
-  /// \param cutOff set the delta potential to 0 for z>zMaxDeltaPot
   /// \param zMaxDeltaPot z position where the maximum delta potential of deltaPot will be set
-  void setIFCChargeUpLinear(const float deltaPot, const float zMaxDeltaPot, const bool cutOff, const Side side);
+  /// \param type function which is used to set falling potential: 0=linear falling off, 1=1/x falling off, 2=1/x steeply falling, 3=linear with offset
+  /// \param offs if offs != 0 the potential doesnt fall to 0. E.g. deltaPot=400V and offs=-10V -> Potential falls from 400V at zMaxDeltaPot to -10V at z=250cm
+  void setIFCChargeUpFallingPot(const float deltaPot, const float zMaxDeltaPot, const int type, const float zEnd, const float offs, const Side side);
 
   /// setting the global corrections directly from input function provided in global coordinates
   /// \param gCorr function returning global corrections for given global coordinate

Detectors/TPC/spacecharge/src/DataContainer3D.cxx

@@ -247,6 +247,13 @@ DataContainer3D<DataT>& DataContainer3D<DataT>::operator-=(const DataContainer3D
   return *this;
 }
 
+template <typename DataT>
+DataContainer3D<DataT>& DataContainer3D<DataT>::operator*=(const DataContainer3D<DataT>& other)
+{
+  std::transform(mData.begin(), mData.end(), other.mData.begin(), mData.begin(), std::multiplies<>());
+  return *this;
+}
+
 template <typename DataT>
 size_t DataContainer3D<DataT>::getIndexZ(size_t index, const int nz, const int nr, const int nphi)
 {

Detectors/TPC/spacecharge/src/SpaceCharge.cxx

@@ -1290,7 +1290,8 @@ void SpaceCharge<DataT>::calcGlobalDistortions(const Fields& formulaStruct, cons
           }
           const DataT z0Tmp = z0 + dzDist + iter * stepSize; // starting z position
 
-          if (getSide(z0Tmp) != side) {
+          // do not do check for first iteration
+          if ((getSide(z0Tmp) != side) && iter) {
             LOGP(error, "Aborting calculation of distortions for iZ: {}, iR: {}, iPhi: {} due to change in the sides!", iZ, iR, iPhi);
             break;
           }
@@ -1561,6 +1562,18 @@ DataT SpaceCharge<DataT>::getPotentialCyl(const DataT z, const DataT r, const Da
   return mInterpolatorPotential[side](z, r, phi);
 }
 
+template <typename DataT>
+std::vector<float> SpaceCharge<DataT>::getPotentialCyl(const std::vector<DataT>& z, const std::vector<DataT>& r, const std::vector<DataT>& phi, const Side side) const
+{
+  const auto nPoints = z.size();
+  std::vector<float> potential(nPoints);
+#pragma omp parallel for num_threads(sNThreads)
+  for (size_t i = 0; i < nPoints; ++i) {
+    potential[i] = getPotentialCyl(z[i], r[i], phi[i], side);
+  }
+  return potential;
+}
+
 template <typename DataT>
 void SpaceCharge<DataT>::getElectricFieldsCyl(const DataT z, const DataT r, const DataT phi, const Side side, DataT& eZ, DataT& eR, DataT& ePhi) const
 {
@@ -1607,7 +1620,7 @@ void SpaceCharge<DataT>::getCorrectionsCyl(const std::vector<DataT>& z, const st
   corrRPhi.resize(nPoints);
 #pragma omp parallel for num_threads(sNThreads)
   for (size_t i = 0; i < nPoints; ++i) {
-    getLocalCorrectionsCyl(z[i], r[i], phi[i], side, corrZ[i], corrR[i], corrRPhi[i]);
+    getCorrectionsCyl(z[i], r[i], phi[i], side, corrZ[i], corrR[i], corrRPhi[i]);
   }
 }
 
@@ -3296,19 +3309,156 @@ void SpaceCharge<DataT>::initRodAlignmentVoltages(const MisalignmentType misalig
 }
 
 template <typename DataT>
-void SpaceCharge<DataT>::setIFCChargeUpLinear(const float deltaPot, const float zMaxDeltaPot, const bool cutOff, const Side side)
+void SpaceCharge<DataT>::addBoundaryPotential(const SpaceCharge<DataT>& other, const Side side, const float scaling)
+{
+  if (other.mPotential[side].getData().empty()) {
+    LOGP(info, "Other space-charge object is empty!");
+    return;
+  }
+
+  if ((mParamGrid.NRVertices != other.mParamGrid.NRVertices) || (mParamGrid.NZVertices != other.mParamGrid.NZVertices) || (mParamGrid.NPhiVertices != other.mParamGrid.NPhiVertices)) {
+    LOGP(info, "Different number of vertices found in input file. Initializing new space charge object with nR {} nZ {} nPhi {} vertices", other.mParamGrid.NRVertices, other.mParamGrid.NZVertices, other.mParamGrid.NPhiVertices);
+    SpaceCharge<DataT> scTmp(mBField.getBField(), other.mParamGrid.NZVertices, other.mParamGrid.NRVertices, other.mParamGrid.NPhiVertices, false);
+    scTmp.mC0 = mC0;
+    scTmp.mC1 = mC1;
+    scTmp.mC2 = mC2;
+    *this = std::move(scTmp);
+  }
+
+  initContainer(mPotential[side], true);
+
+  for (size_t iPhi = 0; iPhi < mParamGrid.NPhiVertices; ++iPhi) {
+    for (size_t iZ = 1; iZ < mParamGrid.NZVertices; ++iZ) {
+      const size_t iRFirst = 0;
+      mPotential[side](iZ, iRFirst, iPhi) += scaling * other.mPotential[side](iZ, iRFirst, iPhi);
+
+      const size_t iRLast = mParamGrid.NRVertices - 1;
+      mPotential[side](iZ, iRLast, iPhi) += scaling * other.mPotential[side](iZ, iRLast, iPhi);
+    }
+  }
+
+  for (size_t iPhi = 0; iPhi < mParamGrid.NPhiVertices; ++iPhi) {
+    for (size_t iR = 0; iR < mParamGrid.NRVertices; ++iR) {
+      const size_t iZFirst = 0;
+      mPotential[side](iZFirst, iR, iPhi) += scaling * other.mPotential[side](iZFirst, iR, iPhi);
+
+      const size_t iZLast = mParamGrid.NZVertices - 1;
+      mPotential[side](iZLast, iR, iPhi) += scaling * other.mPotential[side](iZLast, iR, iPhi);
+    }
+  }
+}
+
+template <typename DataT>
+void SpaceCharge<DataT>::resetBoundaryPotentialToZeroInRangeZ(float zMin, float zMax, const Side side)
 {
-  std::function<DataT(DataT)> chargeUpIFCLinear = [cutOff, zMax = getZMax(Side::A), deltaPot, zMaxDeltaPot](const DataT z) {
+  const float zMaxAbs = std::abs(zMax);
+  for (size_t iPhi = 0; iPhi < mParamGrid.NPhiVertices; ++iPhi) {
+    for (size_t iZ = 1; iZ < mParamGrid.NZVertices; ++iZ) {
+      const DataT z = std::abs(getZVertex(iZ, side));
+      if ((z < zMin) || (z > zMax)) {
+        const size_t iRFirst = 0;
+        mPotential[side](iZ, iRFirst, iPhi) = 0;
+
+        const size_t iRLast = mParamGrid.NRVertices - 1;
+        mPotential[side](iZ, iRLast, iPhi) = 0;
+      }
+    }
+  }
+
+  for (size_t iPhi = 0; iPhi < mParamGrid.NPhiVertices; ++iPhi) {
+    for (size_t iR = 0; iR < mParamGrid.NRVertices; ++iR) {
+      const size_t iZFirst = 0;
+      const float zFirst = std::abs(getZVertex(iZFirst, side));
+      if ((zFirst < zMin) || (zFirst > zMax)) {
+        mPotential[side](iZFirst, iR, iPhi) = 0;
+      }
+
+      const size_t iZLast = mParamGrid.NZVertices - 1;
+      const float zLast = std::abs(getZVertex(iZLast, side));
+      if ((zLast < zMin) || (zLast > zMax)) {
+        mPotential[side](iZLast, iR, iPhi) = 0;
+      }
+    }
+  }
+}
+
+template <typename DataT>
+void SpaceCharge<DataT>::setIFCChargeUpRisingPot(const float deltaPot, const float zMaxDeltaPot, const int type, const float zStart, const float offs, const Side side)
+{
+  std::function<DataT(DataT)> chargeUpIFCLinear = [zStart, type, offs, deltaPot, zMaxDeltaPot](const DataT z) {
     const float absZ = std::abs(z);
     const float absZMaxDeltaPot = std::abs(zMaxDeltaPot);
-    if (absZ <= absZMaxDeltaPot) {
-      return static_cast<DataT>(deltaPot / absZMaxDeltaPot * absZ);
+    if ((absZ <= absZMaxDeltaPot) && (absZ >= zStart)) {
+      // 1/x
+      if (type == 1) {
+        const float offsZ = 1;
+        const float zMaxDeltaPotTmp = zMaxDeltaPot - zStart + offsZ;
+        const float p1 = deltaPot / (1 / offsZ - 1 / zMaxDeltaPotTmp);
+        const float p2 = -p1 / zMaxDeltaPotTmp;
+        const float absZShifted = zMaxDeltaPotTmp - (absZ - zStart);
+        DataT pot = p2 + p1 / absZShifted;
+        return pot;
+      } else if (type == 0 || type == 4) {
+        // linearly rising potential
+        return static_cast<DataT>(deltaPot / (absZMaxDeltaPot - zStart) * (absZ - zStart) + offs);
+      } else if (type == 2) {
+        // flat
+        return DataT(deltaPot);
+      } else if (type == 3) {
+        // linear falling
+        return static_cast<DataT>(-deltaPot / (absZMaxDeltaPot - zStart) * (absZ - zStart) + deltaPot);
+      } else {
+        return DataT(0);
+      }
+    } else if (type == 4) {
+      // flat no z dependence
+      return DataT(offs);
     } else {
-      if (cutOff) {
+      return DataT(0);
+    }
+  };
+  setPotentialBoundaryInnerRadius(chargeUpIFCLinear, side);
+}
+
+template <typename DataT>
+void SpaceCharge<DataT>::setIFCChargeUpFallingPot(const float deltaPot, const float zMaxDeltaPot, const int type, const float zEnd, const float offs, const Side side)
+{
+  std::function<DataT(DataT)> chargeUpIFCLinear = [zEnd, type, offs, zMax = getZMax(Side::A), deltaPot, zMaxDeltaPot](const DataT z) {
+    const float absZ = std::abs(z);
+    const float absZMaxDeltaPot = std::abs(zMaxDeltaPot);
+
+    bool check = (absZ >= absZMaxDeltaPot);
+    if (type == 0 || type == 3) {
+      check = (absZ >= absZMaxDeltaPot);
+    }
+
+    if (check && (absZ <= zEnd)) {
+      // 1/x dependency
+      if (type == 1) {
+        const float p1 = (deltaPot - offs) / (1 / zMaxDeltaPot - 1 / zEnd);
+        const float p2 = offs - p1 / zEnd;
+        DataT pot = p2 + p1 / absZ;
+        return pot;
+      } else if (type == 2) {
+        // 1/x dependency steep fall off!
+        const float offsZ = 1;
+        const float zEndTmp = zEnd - zMaxDeltaPot + offsZ;
+        const float p1 = deltaPot / (1 / offsZ - 1 / zEndTmp);
+        const float p2 = -p1 / zEndTmp;
+        const float absZShifted = absZ - zMaxDeltaPot;
+        DataT pot = p2 + p1 / absZShifted;
+        return pot;
+      } else if (type == 0 || type == 3) {
+        // linearly falling potential
+        const float zPos = absZ - zEnd;
+        return static_cast<DataT>(deltaPot / (absZMaxDeltaPot - zEnd) * zPos + offs);
+      } else {
         return DataT(0);
       }
-      const float zPos = absZ - zMax;
-      return static_cast<DataT>(deltaPot / (absZMaxDeltaPot - zMax) * zPos);
+    } else if (type == 3) {
+      return DataT(offs);
+    } else {
+      return DataT(0);
     }
   };
   setPotentialBoundaryInnerRadius(chargeUpIFCLinear, side);
@@ -3519,6 +3669,31 @@ void SpaceCharge<DataT>::scaleChargeDensitySector(const float scalingFactor, con
   }
 }
 
+template <typename DataT>
+void SpaceCharge<DataT>::scaleChargeDensityStack(const float scalingFactor, const Sector sector, const GEMstack stack)
+{
+  const Side side = sector.side();
+  initContainer(mDensity[side], true);
+  const int verticesPerSector = mParamGrid.NPhiVertices / SECTORSPERSIDE;
+  const int sectorInSide = sector % SECTORSPERSIDE;
+  const int iPhiFirst = sectorInSide * verticesPerSector;
+  const int iPhiLast = iPhiFirst + verticesPerSector;
+  for (unsigned int iR = 0; iR < mParamGrid.NRVertices; ++iR) {
+    const DataT radius = getRVertex(iR, side);
+    for (unsigned int iPhi = iPhiFirst; iPhi < iPhiLast; ++iPhi) {
+      const DataT phi = getPhiVertex(iR, side);
+      const GlobalPosition3D pos(getXFromPolar(radius, phi), getYFromPolar(radius, phi), ((side == Side::A) ? 10 : -10));
+      const auto& mapper = o2::tpc::Mapper::instance();
+      const o2::tpc::DigitPos digiPadPos = mapper.findDigitPosFromGlobalPosition(pos);
+      if (digiPadPos.isValid() && digiPadPos.getCRU().gemStack() == stack) {
+        for (unsigned int iZ = 0; iZ < mParamGrid.NZVertices; ++iZ) {
+          mDensity[side](iZ, iR, iPhi) *= scalingFactor;
+        }
+      }
+    }
+  }
+}
+
 using DataTD = double;
 template class o2::tpc::SpaceCharge<DataTD>;