diff --git a/Detectors/TPC/spacecharge/include/TPCSpaceCharge/DataContainer3D.h b/Detectors/TPC/spacecharge/include/TPCSpaceCharge/DataContainer3D.h index 9c457141e5cee..79400c3d3d214 100644 --- a/Detectors/TPC/spacecharge/include/TPCSpaceCharge/DataContainer3D.h +++ b/Detectors/TPC/spacecharge/include/TPCSpaceCharge/DataContainer3D.h @@ -189,6 +189,7 @@ struct DataContainer3D { /// operator overload DataContainer3D& operator*=(const DataT value); DataContainer3D& operator+=(const DataContainer3D& other); + DataContainer3D& operator*=(const DataContainer3D& other); DataContainer3D& operator-=(const DataContainer3D& other); private: diff --git a/Detectors/TPC/spacecharge/include/TPCSpaceCharge/SpaceCharge.h b/Detectors/TPC/spacecharge/include/TPCSpaceCharge/SpaceCharge.h index 62109f009e436..014be1f2cc50e 100644 --- a/Detectors/TPC/spacecharge/include/TPCSpaceCharge/SpaceCharge.h +++ b/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& 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& other, const Side side, const float scaling = 1); + + /// setting the boundary potential to 0 for zzMax + 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& 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& otherSC); @@ -320,6 +331,9 @@ class SpaceCharge template > 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 getPotentialCyl(const std::vector& z, const std::vector& r, const std::vector& 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 diff --git a/Detectors/TPC/spacecharge/src/DataContainer3D.cxx b/Detectors/TPC/spacecharge/src/DataContainer3D.cxx index c3ca2cbe75b11..c9a39e940873d 100644 --- a/Detectors/TPC/spacecharge/src/DataContainer3D.cxx +++ b/Detectors/TPC/spacecharge/src/DataContainer3D.cxx @@ -247,6 +247,13 @@ DataContainer3D& DataContainer3D::operator-=(const DataContainer3D return *this; } +template +DataContainer3D& DataContainer3D::operator*=(const DataContainer3D& other) +{ + std::transform(mData.begin(), mData.end(), other.mData.begin(), mData.begin(), std::multiplies<>()); + return *this; +} + template size_t DataContainer3D::getIndexZ(size_t index, const int nz, const int nr, const int nphi) { diff --git a/Detectors/TPC/spacecharge/src/SpaceCharge.cxx b/Detectors/TPC/spacecharge/src/SpaceCharge.cxx index 96e86853deae5..3ed00345159bc 100644 --- a/Detectors/TPC/spacecharge/src/SpaceCharge.cxx +++ b/Detectors/TPC/spacecharge/src/SpaceCharge.cxx @@ -1290,7 +1290,8 @@ void SpaceCharge::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::getPotentialCyl(const DataT z, const DataT r, const Da return mInterpolatorPotential[side](z, r, phi); } +template +std::vector SpaceCharge::getPotentialCyl(const std::vector& z, const std::vector& r, const std::vector& phi, const Side side) const +{ + const auto nPoints = z.size(); + std::vector 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 void SpaceCharge::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::getCorrectionsCyl(const std::vector& 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::initRodAlignmentVoltages(const MisalignmentType misalig } template -void SpaceCharge::setIFCChargeUpLinear(const float deltaPot, const float zMaxDeltaPot, const bool cutOff, const Side side) +void SpaceCharge::addBoundaryPotential(const SpaceCharge& 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 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 +void SpaceCharge::resetBoundaryPotentialToZeroInRangeZ(float zMin, float zMax, const Side side) { - std::function 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 +void SpaceCharge::setIFCChargeUpRisingPot(const float deltaPot, const float zMaxDeltaPot, const int type, const float zStart, const float offs, const Side side) +{ + std::function 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(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(deltaPot / (absZMaxDeltaPot - zStart) * (absZ - zStart) + offs); + } else if (type == 2) { + // flat + return DataT(deltaPot); + } else if (type == 3) { + // linear falling + return static_cast(-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 +void SpaceCharge::setIFCChargeUpFallingPot(const float deltaPot, const float zMaxDeltaPot, const int type, const float zEnd, const float offs, const Side side) +{ + std::function 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(deltaPot / (absZMaxDeltaPot - zEnd) * zPos + offs); + } else { return DataT(0); } - const float zPos = absZ - zMax; - return static_cast(deltaPot / (absZMaxDeltaPot - zMax) * zPos); + } else if (type == 3) { + return DataT(offs); + } else { + return DataT(0); } }; setPotentialBoundaryInnerRadius(chargeUpIFCLinear, side); @@ -3519,6 +3669,31 @@ void SpaceCharge::scaleChargeDensitySector(const float scalingFactor, con } } +template +void SpaceCharge::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;