Skip to content

Commit df01596

Browse files
sgorbunodavidrohr
sgorbuno
authored and
davidrohr
committed
TPC Splines: fast merge of SC corrections
1 parent 5d802af commit df01596

File tree

5 files changed

+401
-122
lines changed

5 files changed

+401
-122
lines changed

‎Detectors/TPC/calibration/include/TPCCalibration/TPCFastSpaceChargeCorrectionHelper.h‎

Lines changed: 9 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -102,7 +102,15 @@ class TPCFastSpaceChargeCorrectionHelper
102102
/// initialise inverse transformation from linear combination of several input corrections
103103
void initInverse(std::vector<o2::gpu::TPCFastSpaceChargeCorrection*>& corrections, const std::vector<float>& scaling, bool prn);
104104

105-
void MergeCorrections(std::vector<o2::gpu::TPCFastSpaceChargeCorrection*>& corrections, const std::vector<float>& scaling, bool prn);
105+
/// merge several corrections
106+
/// \param mainCorrection main correction
107+
/// \param scale scaling factor for the main correction
108+
/// \param additionalCorrections vector of pairs of additional corrections and their scaling factors
109+
/// \param prn printout flag
110+
/// \return main correction merged with additional corrections
111+
void MergeCorrections(
112+
o2::gpu::TPCFastSpaceChargeCorrection& mainCorrection, float scale,
113+
const std::vector<std::pair<const o2::gpu::TPCFastSpaceChargeCorrection*, float>>& additionalCorrections, bool prn);
106114

107115
private:
108116
/// geometry initialization

‎Detectors/TPC/calibration/src/TPCFastSpaceChargeCorrectionHelper.cxx‎

Lines changed: 107 additions & 56 deletions
Original file line numberDiff line numberDiff line change
@@ -176,8 +176,8 @@ void TPCFastSpaceChargeCorrectionHelper::fillSpaceChargeCorrectionFromMap(TPCFas
176176
}
177177

178178
if (processingInverseCorrection) {
179-
float* splineX = correction.getSplineData(sector, row, 1);
180-
float* splineYZ = correction.getSplineData(sector, row, 2);
179+
float* splineX = correction.getSplineDataInvX(sector, row);
180+
float* splineYZ = correction.getSplineDataInvYZ(sector, row);
181181
for (int i = 0; i < spline.getNumberOfParameters() / 3; i++) {
182182
splineX[i] = splineParameters[3 * i + 0];
183183
splineYZ[2 * i + 0] = splineParameters[3 * i + 1];
@@ -940,8 +940,8 @@ void TPCFastSpaceChargeCorrectionHelper::initInverse(std::vector<o2::gpu::TPCFas
940940
dataPointGridU.data(), dataPointGridV.data(),
941941
dataPointF.data(), nDataPoints);
942942

943-
float* splineX = correction.getSplineData(sector, row, 1);
944-
float* splineUV = correction.getSplineData(sector, row, 2);
943+
float* splineX = correction.getSplineDataInvX(sector, row);
944+
float* splineUV = correction.getSplineDataInvYZ(sector, row);
945945
for (int i = 0; i < spline.getNumberOfParameters() / 3; i++) {
946946
splineX[i] = splineParameters[3 * i + 0];
947947
splineUV[2 * i + 0] = splineParameters[3 * i + 1];
@@ -967,77 +967,129 @@ void TPCFastSpaceChargeCorrectionHelper::initInverse(std::vector<o2::gpu::TPCFas
967967
LOGP(info, "Inverse tooks: {}s", duration);
968968
}
969969

970-
void TPCFastSpaceChargeCorrectionHelper::MergeCorrections(std::vector<o2::gpu::TPCFastSpaceChargeCorrection*>& corrections, const std::vector<float>& scaling, bool prn)
970+
void TPCFastSpaceChargeCorrectionHelper::MergeCorrections(
971+
o2::gpu::TPCFastSpaceChargeCorrection& mainCorrection, float mainScale,
972+
const std::vector<std::pair<const o2::gpu::TPCFastSpaceChargeCorrection*, float>>& additionalCorrections, bool /*prn*/)
971973
{
972974
/// merge several corrections
973-
/*
975+
974976
TStopwatch watch;
975977
LOG(info) << "fast space charge correction helper: Merge corrections";
976978

977-
if (corrections.size() != scaling.size()) {
978-
LOGP(error, "Input corrections and scaling values have different size");
979-
return;
980-
}
981-
982-
auto& correction = *(corrections.front());
979+
const auto& geo = mainCorrection.getGeometry();
983980

984-
for (int sector = 0; sector < mGeo.getNumberOfSectors(); sector++) {
981+
for (int sector = 0; sector < geo.getNumberOfSectors(); sector++) {
985982

986983
auto myThread = [&](int iThread) {
987-
for (int row = iThread; row < mGeo.getNumberOfRows(); row += mNthreads) {
988-
TPCFastSpaceChargeCorrection::SplineType spline = correction.getSpline(sector, row);
984+
for (int row = iThread; row < geo.getNumberOfRows(); row += mNthreads) {
985+
const auto& spline = mainCorrection.getSpline(sector, row);
989986

990-
std::vector<float> splineParameters(spline.getNumberOfParameters());
991-
std::vector<float> splineParametersInvX(spline.getNumberOfParameters());
992-
std::vector<float> splineParametersInvYZ(spline.getNumberOfParameters());
987+
float* splineParameters = mainCorrection.getSplineData(sector, row);
988+
float* splineParametersInvX = mainCorrection.getSplineDataInvX(sector, row);
989+
float* splineParametersInvYZ = mainCorrection.getSplineDataInvYZ(sector, row);
993990

994-
const auto& gridU = spline.getGridX1();
995-
const auto& gridV = spline.getGridX2();
996-
997-
for (int iu = 0; iu < gridU.getNumberOfKnots(); iu++) {
998-
double u = gridU.getKnot(iu).u;
999-
for (int iv = 0; iv < gridV.getNumberOfKnots(); iv++) {
1000-
int knotIndex = spline.getKnotIndex(iu, iv);
991+
auto& secRowInfo = mainCorrection.getSectorRowInfo(sector, row);
1001992

1002-
double v = gridV.getKnot(iu).u;
1003-
auto [y, z] = correction.convGridToLocal(sector, row, u, v);
1004-
constexpr int nKnotPar1d = 4;
1005-
constexpr int nKnotPar2d = nKnotPar1d * 2;
1006-
constexpr int nKnotPar3d = nKnotPar1d * 3;
993+
constexpr int nKnotPar1d = 4;
994+
constexpr int nKnotPar2d = nKnotPar1d * 2;
995+
constexpr int nKnotPar3d = nKnotPar1d * 3;
1007996

1008-
for (int i = 0; i < corrections.size(); ++i) {
1009-
double s = scaling[i];
1010-
auto p = corrections[i]->getCorrectionParameters(sector, row, y, z);
1011-
for (int j = 0; j < nKnotPar3d; ++j) {
1012-
splineParameters[knotIndex * nKnotPar3d + j] += s * p[j];
997+
{ // scale the main correction
998+
for (int i = 0; i < 3; i++) {
999+
secRowInfo.maxCorr[i] *= mainScale;
1000+
secRowInfo.minCorr[i] *= mainScale;
1001+
}
1002+
double parscale[4] = {mainScale, mainScale, mainScale, mainScale * mainScale};
1003+
for (int iknot = 0, ind = 0; iknot < spline.getNumberOfKnots(); iknot++) {
1004+
for (int ipar = 0; ipar < nKnotPar1d; ++ipar) {
1005+
for (int idim = 0; idim < 3; idim++, ind++) {
1006+
splineParameters[ind] *= parscale[ipar];
10131007
}
1014-
auto pInvX = corrections[i]->getCorrectionParametersInvX(sector, row, y, z);
1015-
for (int j = 0; j < nKnotPar1d; ++j) {
1016-
splineParametersInvX[knotIndex * nKnotPar1d + j] += s * pInvX[j];
1008+
}
1009+
}
1010+
for (int iknot = 0, ind = 0; iknot < spline.getNumberOfKnots(); iknot++) {
1011+
for (int ipar = 0; ipar < nKnotPar1d; ++ipar) {
1012+
for (int idim = 0; idim < 1; idim++, ind++) {
1013+
splineParametersInvX[ind] *= parscale[ipar];
10171014
}
1018-
auto pInvYZ = corrections[i]->getCorrectionParametersInvYZ(sector, row, y, z);
1019-
for (int j = 0; j < nKnotPar2d; ++j) {
1020-
splineParametersInvYZ[knotIndex * nKnotPar2d + j] += s * pInvYZ[j];
1015+
}
1016+
}
1017+
for (int iknot = 0, ind = 0; iknot < spline.getNumberOfKnots(); iknot++) {
1018+
for (int ipar = 0; ipar < nKnotPar1d; ++ipar) {
1019+
for (int idim = 0; idim < 2; idim++, ind++) {
1020+
splineParametersInvYZ[ind] *= parscale[ipar];
10211021
}
10221022
}
1023-
} // iv
1024-
} // iu
1023+
}
1024+
}
10251025

1026-
float* splineXYZ = correction.getSplineData(sector, row, 0);
1027-
float* splineInvX = correction.getSplineData(sector, row, 1);
1028-
float* splineInvYZ = correction.getSplineData(sector, row, 2);
1026+
// add the other corrections
10291027

1030-
for (int i = 0; i < spline.getNumberOfParameters(); i++) {
1031-
splineXYZ[i] = splineParameters[i];
1032-
}
1033-
for (int i = 0; i < spline.getNumberOfParameters() / 3; i++) {
1034-
splineX[i] = splineParametersInvX[i];
1035-
splineYZ[2 * i + 0] = splineParametersInvYZ[2 * i + 0];
1036-
splineYZ[2 * i + 1] = splineParametersInvYZ[2 * i + 1];
1037-
}
1028+
const auto& gridU = spline.getGridX1();
1029+
const auto& gridV = spline.getGridX2();
1030+
1031+
for (int icorr = 0; icorr < additionalCorrections.size(); ++icorr) {
1032+
const auto& corr = *(additionalCorrections[icorr].first);
1033+
double scale = additionalCorrections[icorr].second;
1034+
auto& linfo = corr.getSectorRowInfo(sector, row);
1035+
secRowInfo.updateMaxValues(linfo.getMaxValues(), scale);
1036+
secRowInfo.updateMaxValues(linfo.getMinValues(), scale);
1037+
1038+
double scaleU = secRowInfo.scaleUtoGrid / linfo.scaleUtoGrid;
1039+
double scaleV = secRowInfo.scaleVtoGrid / linfo.scaleVtoGrid;
1040+
1041+
for (int iu = 0; iu < gridU.getNumberOfKnots(); iu++) {
1042+
double u = gridU.getKnot(iu).u;
1043+
for (int iv = 0; iv < gridV.getNumberOfKnots(); iv++) {
1044+
double v = gridV.getKnot(iu).u;
1045+
int knotIndex = spline.getKnotIndex(iu, iv);
1046+
float P[nKnotPar3d];
1047+
1048+
{ // direct correction
1049+
auto [y, z] = mainCorrection.convGridToLocal(sector, row, u, v);
1050+
// return values: u, v, scaling factor
1051+
auto [lu, lv, ls] = corr.convLocalToGrid(sector, row, y, z);
1052+
ls *= scale;
1053+
double parscale[4] = {ls, ls * scaleU, ls * scaleV, ls * ls * scaleU * scaleV};
1054+
const auto& spl = corr.getSpline(sector, row);
1055+
spl.interpolateParametersAtU(corr.getSplineData(sector, row), lu, lv, P);
1056+
for (int ipar = 0, ind = 0; ipar < nKnotPar1d; ++ipar) {
1057+
for (int idim = 0; idim < 3; idim++, ind++) {
1058+
splineParameters[knotIndex * nKnotPar3d + ind] += parscale[ipar] * P[ind];
1059+
}
1060+
}
1061+
}
1062+
1063+
auto [y, z] = mainCorrection.convGridToCorrectedLocal(sector, row, u, v);
1064+
// return values: u, v, scaling factor
1065+
auto [lu, lv, ls] = corr.convCorrectedLocalToGrid(sector, row, y, z);
1066+
ls *= scale;
1067+
double parscale[4] = {ls, ls * scaleU, ls * scaleV, ls * ls * scaleU * scaleV};
1068+
1069+
{ // inverse X correction
1070+
corr.getSplineInvX(sector, row).interpolateParametersAtU(corr.getSplineDataInvX(sector, row), lu, lv, P);
1071+
for (int ipar = 0, ind = 0; ipar < nKnotPar1d; ++ipar) {
1072+
for (int idim = 0; idim < 1; idim++, ind++) {
1073+
splineParametersInvX[knotIndex * nKnotPar1d + ind] += parscale[ipar] * P[ind];
1074+
}
1075+
}
1076+
}
1077+
1078+
{ // inverse YZ correction
1079+
corr.getSplineInvYZ(sector, row).interpolateParametersAtU(corr.getSplineDataInvYZ(sector, row), lu, lv, P);
1080+
for (int ipar = 0, ind = 0; ipar < nKnotPar1d; ++ipar) {
1081+
for (int idim = 0; idim < 2; idim++, ind++) {
1082+
splineParametersInvYZ[knotIndex * nKnotPar2d + ind] += parscale[ipar] * P[ind];
1083+
}
1084+
}
1085+
}
1086+
1087+
} // iv
1088+
} // iu
1089+
} // corrections
10381090

10391091
} // row
1040-
}; // thread
1092+
}; // thread
10411093

10421094
std::vector<std::thread> threads(mNthreads);
10431095

@@ -1054,7 +1106,6 @@ void TPCFastSpaceChargeCorrectionHelper::MergeCorrections(std::vector<o2::gpu::T
10541106
} // sector
10551107
float duration = watch.RealTime();
10561108
LOGP(info, "Merge of corrections tooks: {}s", duration);
1057-
*/
10581109
}
10591110

10601111
} // namespace tpc

‎GPU/TPCFastTransformation/Spline1DSpec.h‎

Lines changed: 39 additions & 49 deletions
Original file line numberDiff line numberDiff line change
@@ -318,7 +318,10 @@ class Spline1DSpec<DataT, YdimT, 0> : public Spline1DContainer<DataT>
318318
for (int32_t dim = 0; dim < nYdim; ++dim) {
319319
S[dim] = dSdSr * Sr[dim] + dSdSl * Sl[dim] + dSdDl * Dl[dim] + dSdDr * Dr[dim];
320320
}
321+
321322
/*
323+
another way to calculate f(u):
324+
322325
if (u < (DataT)0) {
323326
u = (DataT)0;
324327
}
@@ -336,18 +339,6 @@ class Spline1DSpec<DataT, YdimT, 0> : public Spline1DContainer<DataT>
336339
S[dim] = ((a * v + b) * v + Dl[dim]) * uu + Sl[dim];
337340
}
338341
*/
339-
/*
340-
another way to calculate f(u):
341-
T uu = T(u - knotL.u);
342-
T v = uu * T(knotL.Li); // scaled u
343-
T vm1 = v-1;
344-
T v2 = v * v;
345-
float cSr = v2*(3-2*v);
346-
float cSl = 1-cSr;
347-
float cDl = v*vm1*vm1*knotL.L;
348-
float cDr = v2*vm1*knotL.L;
349-
return cSl*Sl + cSr*Sr + cDl*Dl + cDr*Dr;
350-
*/
351342
}
352343

353344
template <typename T>
@@ -365,51 +356,50 @@ class Spline1DSpec<DataT, YdimT, 0> : public Spline1DContainer<DataT>
365356

366357
u = u - knotL.u;
367358
T v = u * T(knotL.Li); // scaled u
368-
T vm1 = v - 1.;
359+
T vm1 = v - T(1.);
369360
T a = u * vm1;
370361
T v2 = v * v;
371-
T dSdSr = v2 * (3. - 2 * v);
372-
T dSdSl = 1. - dSdSr;
362+
T dSdSr = v2 * (T(3.) - v - v);
363+
T dSdSl = T(1.) - dSdSr;
373364
T dSdDl = vm1 * a;
374365
T dSdDr = v * a;
375366
// S(u) = dSdSl * Sl + dSdSr * Sr + dSdDl * Dl + dSdDr * Dr;
376367
return std::make_tuple(dSdSl, dSdDl, dSdSr, dSdDr);
377368
}
378-
/*
379-
template <typename T>
380-
GPUd() void getUsecondDerivatives(const Knot& knotL, DataT u,
381-
T& dSl, T& dDl, T& dSr, T& dDr,
382-
T& dSl2, T& dDl2, T& dSr2, T& dDr2) const
383-
{
384-
/// Get derivatives of the interpolated value {S(u): 1D -> nYdim} at the segment [knotL, next knotR]
385-
/// over the spline values Sl, Sr and the slopes Dl, Dr
386-
387-
if (u < (DataT)0) {
388-
u = (DataT)0;
389-
}
390-
if (u > (DataT)TBase::getUmax()) {
391-
u = (DataT)TBase::getUmax();
392-
}
393-
394-
u = u - knotL.u;
395-
T v = u * T(knotL.Li); // scaled u
396-
T vm1 = v - 1.;
397-
T a = u * vm1;
398-
T v2 = v * v;
399-
dSr = v2 * (3. - 2 * v);
400-
dSl = 1. - dSr;
401-
dDl = vm1 * a;
402-
dDr = v * a;
403-
T dv = T(knotL.Li);
404-
dSr2 = 6. * v * (1. - v) * dv;
405-
dSl2 = -dSr2;
406-
dDl2 = (v - 1) * (3 * v - 1);
407-
dDr = u * (v * v - v);
408-
dDr2 = 3.f * v * v - 2.f * v;
409-
// F(u) = dSl * Sl + dSr * Sr + dDl * Dl + dDr * Dr;
410-
// dF(u)/du = dSl2 * Sl + dSr2 * Sr + dDl2 * Dl + dDr2 * Dr;
369+
370+
template <typename T>
371+
GPUd() std::tuple<T, T, T, T, T, T, T, T> getSDderivativesOverParsAtU(const Knot& knotL, DataT u) const
372+
{
373+
/// Get derivatives of the interpolated value {S(u): 1D -> nYdim} at the segment [knotL, next knotR]
374+
/// over the spline values Sl, Sr and the slopes Dl, Dr
375+
376+
if (u < (DataT)0) {
377+
u = (DataT)0;
411378
}
412-
*/
379+
if (u > (DataT)TBase::getUmax()) {
380+
u = (DataT)TBase::getUmax();
381+
}
382+
383+
u = u - knotL.u;
384+
T v = u * T(knotL.Li); // scaled u
385+
T vm1 = v - T(1.);
386+
T a = u * vm1;
387+
T v2 = v * v;
388+
T dSdSr = v2 * (T(3.) - v - v);
389+
T dSdSl = T(1.) - dSdSr;
390+
T dSdDl = vm1 * a;
391+
T dSdDr = v * a;
392+
393+
T dv = T(knotL.Li);
394+
T dDdSr = 6. * v * (T(1.) - v) * dv;
395+
T dDdSl = -dDdSr;
396+
T dDdDl = vm1 * (v + v + vm1);
397+
T dDdDr = v * (v + vm1 + vm1);
398+
// S(u) = dSdSl * Sl + dSdSr * Sr + dSdDl * Dl + dSdDr * Dr;
399+
// D(u) = dS(u)/du = dDdSl * Sl + dDdSr * Sr + dDdDl * Dl + dDdDr * Dr;
400+
return std::make_tuple(dSdSl, dSdDl, dSdSr, dSdDr, dDdSl, dDdDl, dDdSr, dDdDr);
401+
}
402+
413403
using TBase::convXtoU;
414404
using TBase::getKnot;
415405
using TBase::getKnots;

0 commit comments

Comments
 (0)