Skip to content

Commit

Permalink
TPC: Optimising simulation of distortions in MC
Browse files Browse the repository at this point in the history
  • Loading branch information
@matthias-kleiner
matthias-kleiner committed Sep 26, 2024
1 parent c8e9b40 commit 7886676
Show file tree
Hide file tree
Showing 3 changed files with 202 additions and 13 deletions.
2 changes: 1 addition & 1 deletion Detectors/TPC/simulation/src/Digitizer.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -257,7 +257,7 @@ void Digitizer::recalculateDistortions()
mSpaceChargeDer->calcGlobalCorrWithGlobalDistIterative(side, nullptr, 0);

LOGP(info, "Calculating scaled distortions with scaling factor {}", mLumiScaleFactor);
mSpaceCharge->calcGlobalDistWithGlobalCorrIterative(side, mSpaceChargeDer.get(), mLumiScaleFactor);
mSpaceCharge->calcGlobalDistWithGlobalCorrIterativeLinearCartesian(side, mSpaceChargeDer.get(), mLumiScaleFactor);
}
// set new lumi of avg map
mSpaceCharge->setMeanLumi(CorrMapParam::Instance().lumiInst);
Expand Down
8 changes: 8 additions & 0 deletions Detectors/TPC/spacecharge/include/TPCSpaceCharge/SpaceCharge.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -369,6 +369,12 @@ class SpaceCharge
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);

/// step 5: calculate global distortions using the global corrections with scaling with second distortion object, which will be consistent with scaled corrections in cartesian coordinates (as done in the tracking)
/// \param scSCale possible second sc object
/// \param scale scaling for second sc object
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);

/// \return returns number of vertices in z direction
unsigned short getNZVertices() const { return mParamGrid.NZVertices; }

Expand Down Expand Up @@ -1210,6 +1216,7 @@ class SpaceCharge
float getMeanLumi() const { return mMeta.meanLumi; }
void setMeanLumi(float lumi) { mMeta.meanLumi = lumi; }
void initAfterReadingFromFile();
void addCorrections(const SpaceCharge<DataT>& otherSC);

/// get DCA in RPhi for high pt track
/// \param tgl tgl of the track
Expand Down Expand Up @@ -1373,6 +1380,7 @@ class SpaceCharge
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);
};
Expand Down
205 changes: 193 additions & 12 deletions Detectors/TPC/spacecharge/src/SpaceCharge.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -758,13 +758,13 @@ void SpaceCharge<DataT>::calcGlobalDistCorrIterative(const DistCorrInterpolator<
for (unsigned int iZ = 1; iZ < mParamGrid.NZVertices; ++iZ) {
const DataT z = getZVertex(iZ, side);

unsigned int nearestiZ = iZ;
unsigned int nearestiR = iR;
unsigned int nearestiPhi = iPhi;
// unsigned int nearestiZ = iZ;
// unsigned int nearestiR = iR;
// unsigned int nearestiPhi = iPhi;

DataT nearestZ = getZVertex(nearestiZ, side);
DataT nearestR = getRVertex(nearestiR, side);
DataT nearestPhi = getPhiVertex(nearestiPhi, side);
// DataT nearestZ = getZVertex(nearestiZ, side);
// DataT nearestR = getRVertex(nearestiR, side);
// DataT nearestPhi = getPhiVertex(nearestiPhi, side);

//
//==========================================================================================
Expand All @@ -774,9 +774,9 @@ void SpaceCharge<DataT>::calcGlobalDistCorrIterative(const DistCorrInterpolator<
// 1. calculate difference from nearest point to query point with stepwidth factor x
// and approach the new point
//
DataT stepR = (radius - nearestR) * approachR;
DataT stepZ = (z - nearestZ) * approachZ;
DataT stepPhi = (phi - nearestPhi) * approachPhi;
DataT stepR = 0;
DataT stepZ = 0;
DataT stepPhi = 0;

// needed to check for convergence
DataT lastCorrdR = std::numeric_limits<DataT>::max();
Expand All @@ -790,9 +790,9 @@ void SpaceCharge<DataT>::calcGlobalDistCorrIterative(const DistCorrInterpolator<

for (int iter = 0; iter < maxIter; ++iter) {
// 2. get new point coordinates
const DataT rCurrPos = getRVertex(nearestiR, side) + stepR;
const DataT zCurrPos = getZVertex(nearestiZ, side) + stepZ;
const DataT phiCurrPos = getPhiVertex(nearestiPhi, side) + stepPhi;
const DataT rCurrPos = radius + stepR;
const DataT zCurrPos = z + stepZ;
const DataT phiCurrPos = phi + stepPhi;

// abort calculation of drift path if electron reached inner/outer field cage or central electrode
if (rCurrPos <= getRMinSim(side) || rCurrPos >= getRMaxSim(side) || getSide(zCurrPos) != side) {
Expand Down Expand Up @@ -867,6 +867,169 @@ void SpaceCharge<DataT>::calcGlobalDistCorrIterative(const DistCorrInterpolator<
}
}

template <typename DataT>
void SpaceCharge<DataT>::calcGlobalDistWithGlobalCorrIterativeLinearCartesian(const Side side, const SpaceCharge<DataT>* scSCale, float scale, const int maxIter, const DataT approachX, const DataT approachY, const DataT approachZ, const DataT diffCorr)
{
calcGlobalDistCorrIterativeLinearCartesian(getGlobalCorrInterpolator(side), maxIter, approachX, approachY, approachZ, diffCorr, scSCale, scale, Type::Distortions);
}

template <typename DataT>
void SpaceCharge<DataT>::calcGlobalDistWithGlobalCorrIterativeLinearCartesian(const SpaceCharge<DataT>* scSCale, float scale, const int maxIter, const DataT approachX, const DataT approachY, const DataT approachZ, const DataT diffCorr)
{
#pragma omp parallel for num_threads(sNThreads)
for (int iside = 0; iside < FNSIDES; ++iside) {
const o2::tpc::Side side = (iside == 0) ? Side::A : Side::C;
calcGlobalDistWithGlobalCorrIterativeLinearCartesian(side, scSCale, scale, maxIter, approachX, approachY, approachZ, diffCorr);
}
}

template <typename DataT>
void SpaceCharge<DataT>::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)
{
const Side side = globCorr.getSide();
if (type == Type::Distortions) {
initContainer(mGlobalDistdR[side], true);
initContainer(mGlobalDistdZ[side], true);
initContainer(mGlobalDistdRPhi[side], true);
} else {
initContainer(mGlobalCorrdR[side], true);
initContainer(mGlobalCorrdZ[side], true);
initContainer(mGlobalCorrdRPhi[side], true);
}

const auto& scSCaleInterpolator = (type == Type::Distortions) ? scSCale->mInterpolatorGlobalCorr[side] : scSCale->mInterpolatorGlobalDist[side];

#pragma omp parallel for num_threads(sNThreads)
for (unsigned int iPhi = 0; iPhi < mParamGrid.NPhiVertices; ++iPhi) {
const DataT phi = getPhiVertex(iPhi, side);
for (unsigned int iR = 0; iR < mParamGrid.NRVertices; ++iR) {
const DataT radius = getRVertex(iR, side);
const DataT x = getXFromPolar(radius, phi);
const DataT y = getYFromPolar(radius, phi);

for (unsigned int iZ = 1; iZ < mParamGrid.NZVertices; ++iZ) {
const DataT z = getZVertex(iZ, side);

//
//==========================================================================================
//==== start algorithm: use tricubic upsampling to numerically approach the query point ====
//==========================================================================================
//
// 1. calculate difference from nearest point to query point with stepwidth factor x
// and approach the new point
//
DataT stepX = 0;
DataT stepY = 0;
DataT stepZ = 0;

// needed to check for convergence
DataT lastCorrX = std::numeric_limits<DataT>::max();
DataT lastCorrY = std::numeric_limits<DataT>::max();
DataT lastCorrZ = std::numeric_limits<DataT>::max();
DataT lastX = std::numeric_limits<DataT>::max();
DataT lastY = std::numeric_limits<DataT>::max();
DataT lastZ = std::numeric_limits<DataT>::max();

for (int iter = 0; iter < maxIter; ++iter) {
// 2. get new point coordinates
const DataT xCurrPos = x + stepX;
const DataT yCurrPos = y + stepY;
const DataT zCurrPos = z + stepZ;

// abort calculation of drift path if electron reached inner/outer field cage or central electrode
const DataT rCurrPos = getRadiusFromCartesian(xCurrPos, yCurrPos);
if (rCurrPos <= getRMinSim(side) || rCurrPos >= getRMaxSim(side) || getSide(zCurrPos) != side) {
break;
}

// interpolate global correction at new point and calculate position of global correction
DataT corrX = 0;
DataT corrY = 0;
DataT corrZ = 0;
getCorrections(xCurrPos, yCurrPos, zCurrPos, side, corrX, corrY, corrZ);

if (scSCale && scale != 0) {
DataT corrXScale = 0;
DataT corrYScale = 0;
DataT corrZScale = 0;
scSCale->getCorrections(xCurrPos, yCurrPos, zCurrPos, side, corrXScale, corrYScale, corrZScale);
corrX += scale * corrXScale;
corrY += scale * corrYScale;
corrZ += scale * corrZScale;
}

const DataT xNewPos = xCurrPos + corrX;
const DataT yNewPos = yCurrPos + corrY;
const DataT zNewPos = zCurrPos + corrZ;

// approach desired coordinate
stepX += (x - xNewPos) * approachX;
stepY += (y - yNewPos) * approachY;
stepZ += (z - zNewPos) * approachZ;

// check for convergence
const DataT diffCorrX = std::abs(corrX - lastCorrX);
const DataT diffCorrY = std::abs(corrY - lastCorrY);
const DataT diffCorrZ = std::abs(corrZ - lastCorrZ);

lastCorrX = corrX;
lastCorrY = corrY;
lastCorrZ = corrZ;
lastX = xCurrPos;
lastY = yCurrPos;
lastZ = zCurrPos;

// stop algorithm if converged
if ((diffCorrX < diffCorr) && (diffCorrY < diffCorr) && (diffCorrZ < diffCorr)) {
break;
}
}
const DataT xNew = lastX + lastCorrX;
const DataT yNew = lastY + lastCorrY;
const DataT radiusNew = getRadiusFromCartesian(xNew, yNew);
const DataT corrdR = -(radiusNew - getRadiusFromCartesian(lastX, lastY));

float phiNew = getPhiFromCartesian(xNew, yNew);
o2::math_utils::bringTo02PiGen(phiNew);

float phiLast = getPhiFromCartesian(lastX, lastY);
o2::math_utils::bringTo02PiGen(phiLast);

DataT deltaPhi = (phiNew - phiLast);
// handle edge cases
if (deltaPhi > PI) {
deltaPhi -= 2 * PI;
} else if (deltaPhi < -PI) {
deltaPhi += 2 * PI;
}
const DataT corrdRPhi = -deltaPhi * radiusNew;

// set global distortions if algorithm converged or iterations exceed max numbers of iterations
if (type == Type::Distortions) {
mGlobalDistdR[side](iZ, iR, iPhi) = corrdR;
mGlobalDistdRPhi[side](iZ, iR, iPhi) = corrdRPhi;
mGlobalDistdZ[side](iZ, iR, iPhi) = -lastCorrZ;
} else {
mGlobalCorrdR[side](iZ, iR, iPhi) = corrdR;
mGlobalCorrdRPhi[side](iZ, iR, iPhi) = corrdRPhi;
mGlobalCorrdZ[side](iZ, iR, iPhi) = -lastCorrZ;
}
}
}
for (unsigned int iR = 0; iR < mParamGrid.NRVertices; ++iR) {
if (type == Type::Distortions) {
mGlobalDistdR[side](0, iR, iPhi) = 3 * (mGlobalDistdR[side](1, iR, iPhi) - mGlobalDistdR[side](2, iR, iPhi)) + mGlobalDistdR[side](3, iR, iPhi);
mGlobalDistdRPhi[side](0, iR, iPhi) = 3 * (mGlobalDistdRPhi[side](1, iR, iPhi) - mGlobalDistdRPhi[side](2, iR, iPhi)) + mGlobalDistdRPhi[side](3, iR, iPhi);
mGlobalDistdZ[side](0, iR, iPhi) = 3 * (mGlobalDistdZ[side](1, iR, iPhi) - mGlobalDistdZ[side](2, iR, iPhi)) + mGlobalDistdZ[side](3, iR, iPhi);
} else {
mGlobalCorrdR[side](0, iR, iPhi) = 3 * (mGlobalCorrdR[side](1, iR, iPhi) - mGlobalCorrdR[side](2, iR, iPhi)) + mGlobalCorrdR[side](3, iR, iPhi);
mGlobalCorrdRPhi[side](0, iR, iPhi) = 3 * (mGlobalCorrdRPhi[side](1, iR, iPhi) - mGlobalCorrdRPhi[side](2, iR, iPhi)) + mGlobalCorrdRPhi[side](3, iR, iPhi);
mGlobalCorrdZ[side](0, iR, iPhi) = 3 * (mGlobalCorrdZ[side](1, iR, iPhi) - mGlobalCorrdZ[side](2, iR, iPhi)) + mGlobalCorrdZ[side](3, iR, iPhi);
}
}
}
}

template <typename DataT>
NumericalFields<DataT> SpaceCharge<DataT>::getElectricFieldsInterpolator(const Side side) const
{
Expand Down Expand Up @@ -3262,6 +3425,24 @@ void SpaceCharge<DataT>::addChargeDensity(const SpaceCharge<DataT>& otherSC)
mDensity[Side::C] += otherSC.mDensity[Side::C];
}

template <typename DataT>
void SpaceCharge<DataT>::addCorrections(const SpaceCharge<DataT>& otherSC)
{
const bool sameGrid = (getNPhiVertices() == otherSC.getNPhiVertices()) && (getNRVertices() == otherSC.getNRVertices()) && (getNZVertices() == otherSC.getNZVertices());
if (!sameGrid) {
LOGP(warning, "Space charge objects have different grid definition");
return;
}

mGlobalCorrdR[Side::A] += otherSC.mGlobalCorrdR[Side::A];
mGlobalCorrdRPhi[Side::A] += otherSC.mGlobalCorrdRPhi[Side::A];
mGlobalCorrdZ[Side::A] += otherSC.mGlobalCorrdZ[Side::A];

mGlobalCorrdR[Side::C] += otherSC.mGlobalCorrdR[Side::C];
mGlobalCorrdRPhi[Side::C] += otherSC.mGlobalCorrdRPhi[Side::C];
mGlobalCorrdZ[Side::C] += otherSC.mGlobalCorrdZ[Side::C];
}

template <typename DataT>
void SpaceCharge<DataT>::fillChargeDensityFromHisto(const char* file, const char* nameA, const char* nameC)
{
Expand Down

0 comments on commit 7886676

Please sign in to comment.