Skip to content

Lightweight metadata for the splines #140

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
davidrohr merged 8 commits into from
Apr 18, 2026
Merged

Lightweight metadata for the splines #140

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
8 commits
Select commit Hold shift + click to select a range
cf1e60c
TPCFastTransform: adjust the test macro to the new code
cbmsw Apr 14, 2026
7557023
TPCFastTransform: remove row-wise max correction values
cbmsw Apr 14, 2026
9971148
TPCFastTransform: mirror Z-´of the spline grid for the C side
cbmsw Apr 15, 2026
180df09
TPCFastTransform: prepare for reduced metadata
cbmsw Apr 15, 2026
50cec39
TPCFastTransform: reduce the metadata: the same spline setup for all …
cbmsw Apr 16, 2026
1753e3a
TPCFastTransform: renaming
cbmsw Apr 16, 2026
e77f201
TPCFastTransform: remove unnecessary memory alignment in the spline data
cbmsw Apr 16, 2026
455aaeb
TPCFastTransform: cleanup
cbmsw Apr 16, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
134 changes: 53 additions & 81 deletions Detectors/TPC/calibration/src/TPCFastSpaceChargeCorrectionHelper.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -137,19 +137,13 @@ void TPCFastSpaceChargeCorrectionHelper::fillSpaceChargeCorrectionFromMap(TPCFas
auto myThread = [&](int iThread) {
for (int row = iThread; row < correction.getGeometry().getNumberOfRows(); row += mNthreads) {

TPCFastSpaceChargeCorrection::SplineType& spline = correction.getSpline(sector, row);
TPCFastSpaceChargeCorrection::SplineType& spline = correction.getSplineForRow(row);
Spline2DHelper<float> helper;
std::vector<float> splineParameters;
splineParameters.resize(spline.getNumberOfParameters());

const std::vector<o2::gpu::TPCFastSpaceChargeCorrectionMap::CorrectionPoint>& data = mCorrectionMap.getPoints(sector, row);
int nDataPoints = data.size();
auto& info = correction.getSectorRowInfo(sector, row);
if (!processingInverseCorrection) {
info.resetMaxValues();
}
info.updateMaxValues(1., 1., 1.);
info.updateMaxValues(-1., -1., -1.);

if (nDataPoints >= 4) {
std::vector<double> pointGU(nDataPoints);
Expand All @@ -170,7 +164,6 @@ void TPCFastSpaceChargeCorrectionHelper::fillSpaceChargeCorrectionFromMap(TPCFas
pointCorr[3 * i + 0] = p.mDx;
pointCorr[3 * i + 1] = p.mDy;
pointCorr[3 * i + 2] = p.mDz;
info.updateMaxValues(5. * p.mDx, 5. * p.mDy, 5. * p.mDz);
}
helper.approximateDataPoints(spline, splineParameters.data(), 0., spline.getGridX1().getUmax(), 0., spline.getGridX2().getUmax(), pointGU.data(),
pointGV.data(), pointCorr.data(), pointWeight.data(), nDataPoints);
Expand Down Expand Up @@ -263,14 +256,12 @@ std::unique_ptr<TPCFastSpaceChargeCorrection> TPCFastSpaceChargeCorrectionHelper
correction.startConstruction(mGeo, nCorrectionScenarios);

// assign spline type for TPC rows
for (int sector = 0; sector < mGeo.getNumberOfSectors(); sector++) {
for (int row = 0; row < mGeo.getNumberOfRows(); row++) {
int scenario = row / 10;
if (scenario >= nCorrectionScenarios) {
scenario = nCorrectionScenarios - 1;
}
correction.setRowScenarioID(sector, row, scenario);
for (int row = 0; row < mGeo.getNumberOfRows(); row++) {
int scenario = row / 10;
if (scenario >= nCorrectionScenarios) {
scenario = nCorrectionScenarios - 1;
}
correction.setRowScenarioID(row, scenario);
}

for (int scenario = 0; scenario < nCorrectionScenarios; scenario++) {
Expand Down Expand Up @@ -408,11 +399,10 @@ std::unique_ptr<o2::gpu::TPCFastSpaceChargeCorrection> TPCFastSpaceChargeCorrect
int nY2Xbins = trackResiduals.getNY2XBins();
int nZ2Xbins = trackResiduals.getNZ2XBins();

std::vector<double> knotsDouble[3];
std::vector<double> knotsDouble[2];

knotsDouble[0].reserve(nY2Xbins);
knotsDouble[1].reserve(nZ2Xbins);
knotsDouble[2].reserve(nZ2Xbins);

// to get enouth measurements, make a spline knot at every second bin. Boundary bins are always included.

Expand All @@ -425,16 +415,14 @@ std::unique_ptr<o2::gpu::TPCFastSpaceChargeCorrection> TPCFastSpaceChargeCorrect

for (int i = 0, j = nZ2Xbins - 1; i <= j; i += 2, j -= 2) {
knotsDouble[1].push_back(trackResiduals.getZ2X(i));
knotsDouble[2].push_back(-trackResiduals.getZ2X(i));
if (j >= i + 1) {
knotsDouble[1].push_back(trackResiduals.getZ2X(j));
knotsDouble[2].push_back(-trackResiduals.getZ2X(j));
}
}

std::vector<int> knotsInt[3];
std::vector<int> knotsInt[2];

for (int dim = 0; dim < 3; dim++) {
for (int dim = 0; dim < 2; dim++) {
auto& knotsD = knotsDouble[dim];
std::sort(knotsD.begin(), knotsD.end());

Expand Down Expand Up @@ -470,12 +458,10 @@ std::unique_ptr<o2::gpu::TPCFastSpaceChargeCorrection> TPCFastSpaceChargeCorrect
}

auto& yKnotsInt = knotsInt[0];
auto& zKnotsIntA = knotsInt[1];
auto& zKnotsIntC = knotsInt[2];
auto& zKnotsInt = knotsInt[1];

int nKnotsY = yKnotsInt.size();
int nKnotsZA = zKnotsIntA.size();
int nKnotsZC = zKnotsIntC.size();
int nKnotsZ = zKnotsInt.size();

// std::cout << "n knots Y: " << nKnotsY << std::endl;
// std::cout << "n knots Z: " << nKnotsZA << ", " << nKnotsZC << std::endl;
Expand All @@ -485,53 +471,41 @@ std::unique_ptr<o2::gpu::TPCFastSpaceChargeCorrection> TPCFastSpaceChargeCorrect

{ // create the correction object

const int nCorrectionScenarios = 2; // different grids for TPC A and TPC C sides
const int nCorrectionScenarios = 1;

correction.startConstruction(geo, nCorrectionScenarios);

// init rows
for (int iSector = 0; iSector < nSectors; iSector++) {
int id = iSector < geo.getNumberOfSectorsA() ? 0 : 1;
for (int row = 0; row < geo.getNumberOfRows(); row++) {
correction.setRowScenarioID(iSector, row, id);
}
for (int row = 0; row < geo.getNumberOfRows(); row++) {
correction.setRowScenarioID(row, 0);
}

{ // init spline scenario
TPCFastSpaceChargeCorrection::SplineType spline;
spline.recreate(nKnotsY, &yKnotsInt[0], nKnotsZA, &zKnotsIntA[0]);
spline.recreate(nKnotsY, &yKnotsInt[0], nKnotsZ, &zKnotsInt[0]);
correction.setSplineScenario(0, spline);
spline.recreate(nKnotsY, &yKnotsInt[0], nKnotsZC, &zKnotsIntC[0]);
correction.setSplineScenario(1, spline);
}
correction.finishConstruction();
} // .. create the correction object

// set the grid borders
for (int iSector = 0; iSector < geo.getNumberOfSectors(); iSector++) {
for (int iRow = 0; iRow < geo.getNumberOfRows(); iRow++) {
auto& info = correction.getSectorRowInfo(iSector, iRow);
const auto& spline = correction.getSpline(iSector, iRow);
double rowX = geo.getRowInfo(iRow).x;
double yMin = rowX * trackResiduals.getY2X(iRow, 0);
double yMax = rowX * trackResiduals.getY2X(iRow, trackResiduals.getNY2XBins() - 1);
double zMin = rowX * trackResiduals.getZ2X(0);
double zMax = rowX * trackResiduals.getZ2X(trackResiduals.getNZ2XBins() - 1);
double zOut = zMax;
if (iSector >= geo.getNumberOfSectorsA()) {
// TPC C side
zOut = -zOut;
zMax = -zMin;
zMin = zOut;
}
info.gridMeasured.set(yMin, spline.getGridX1().getUmax() / (yMax - yMin), // y
zMin, spline.getGridX2().getUmax() / (zMax - zMin), // z
zOut, geo.getZreadout(iSector)); // correction scaling region

info.gridReal = info.gridMeasured;

// std::cout << " iSector " << iSector << " iRow " << iRow << " uMin: " << uMin << " uMax: " << uMax << " vMin: " << vMin << " vMax: " << vMax
//<< " grid scale u "<< info.scaleUtoGrid << " grid scale v "<< info.scaleVtoGrid<< std::endl;
}
for (int iRow = 0; iRow < geo.getNumberOfRows(); iRow++) {
auto& info = correction.getRowInfo(iRow);
const auto& spline = correction.getSplineForRow(iRow);
double rowX = geo.getRowInfo(iRow).x;
double yMin = rowX * trackResiduals.getY2X(iRow, 0);
double yMax = rowX * trackResiduals.getY2X(iRow, trackResiduals.getNY2XBins() - 1);
double zMin = rowX * trackResiduals.getZ2X(0);
double zMax = rowX * trackResiduals.getZ2X(trackResiduals.getNZ2XBins() - 1);
double zOut = zMax;
info.gridMeasured.set(yMin, spline.getGridX1().getUmax() / (yMax - yMin), // y
zMin, spline.getGridX2().getUmax() / (zMax - zMin), // z
zOut, geo.getTPCzLength()); // correction scaling region

info.gridReal = info.gridMeasured;

// std::cout << " iSector " << iSector << " iRow " << iRow << " uMin: " << uMin << " uMax: " << uMax << " vMin: " << vMin << " vMax: " << vMax
//<< " grid scale u "<< info.scaleUtoGrid << " grid scale v "<< info.scaleVtoGrid<< std::endl;
}

LOG(info) << "fast space charge correction helper: preparation took " << watch1.RealTime() << "s";
Expand Down Expand Up @@ -783,8 +757,8 @@ std::unique_ptr<o2::gpu::TPCFastSpaceChargeCorrection> TPCFastSpaceChargeCorrect

// feed the row data to the helper

auto& info = correction.getSectorRowInfo(iSector, iRow);
const auto& spline = correction.getSpline(iSector, iRow);
auto& info = correction.getRowInfo(iRow);
const auto& spline = correction.getSplineForRow(iRow);

auto addVoxel = [&](int iy, int iz, double weight) {
auto& vox = vRowVoxels[iy * nZ2Xbins + iz];
Expand Down Expand Up @@ -928,6 +902,11 @@ void TPCFastSpaceChargeCorrectionHelper::initInverse(std::vector<o2::gpu::TPCFas
tpcR2max = tpcR2max / cos(2 * M_PI / mGeo.getNumberOfSectorsA() / 2) + 1.;
tpcR2max = tpcR2max * tpcR2max;

for (int row = 0; row < mGeo.getNumberOfRows(); row++) {
auto& rowInfo = correction.getRowInfo(row);
rowInfo.gridReal = rowInfo.gridMeasured;
}

for (int sector = 0; sector < mGeo.getNumberOfSectors(); sector++) {
// LOG(info) << "inverse transform for sector " << sector ;

Expand All @@ -936,10 +915,9 @@ void TPCFastSpaceChargeCorrectionHelper::initInverse(std::vector<o2::gpu::TPCFas
std::vector<float> splineParameters;

for (int row = iThread; row < mGeo.getNumberOfRows(); row += mNthreads) {
auto& sectorRowInfo = correction.getSectorRowInfo(sector, row);
sectorRowInfo.gridReal = sectorRowInfo.gridMeasured;
auto& rowInfo = correction.getRowInfo(row);

TPCFastSpaceChargeCorrection::SplineType spline = correction.getSpline(sector, row);
TPCFastSpaceChargeCorrection::SplineType spline = correction.getSplineForRow(row);
helper.setSpline(spline, 10, 10);

std::vector<double> gridU;
Expand Down Expand Up @@ -1050,23 +1028,19 @@ void TPCFastSpaceChargeCorrectionHelper::mergeCorrections(

auto myThread = [&](int iThread) {
for (int row = iThread; row < geo.getNumberOfRows(); row += mNthreads) {
const auto& spline = mainCorrection.getSpline(sector, row);
auto& rowInfo = mainCorrection.getRowInfo(row);
const auto& spline = mainCorrection.getSplineForRow(row);

float* splineParameters = mainCorrection.getCorrectionData(sector, row);
float* splineParametersInvX = mainCorrection.getCorrectionDataInvX(sector, row);
float* splineParametersInvYZ = mainCorrection.getCorrectionDataInvYZ(sector, row);

auto& secRowInfo = mainCorrection.getSectorRowInfo(sector, row);

constexpr int nKnotPar1d = 4;
constexpr int nKnotPar2d = nKnotPar1d * 2;
constexpr int nKnotPar3d = nKnotPar1d * 3;

{ // scale the main correction
for (int i = 0; i < 3; i++) {
secRowInfo.maxCorr[i] *= mainScale;
secRowInfo.minCorr[i] *= mainScale;
}

double parscale[4] = {mainScale, mainScale, mainScale, mainScale * mainScale};
for (int iknot = 0, ind = 0; iknot < spline.getNumberOfKnots(); iknot++) {
for (int ipar = 0; ipar < nKnotPar1d; ++ipar) {
Expand Down Expand Up @@ -1099,14 +1073,12 @@ void TPCFastSpaceChargeCorrectionHelper::mergeCorrections(
for (int icorr = 0; icorr < additionalCorrections.size(); ++icorr) {
const auto& corr = *(additionalCorrections[icorr].first);
double scale = additionalCorrections[icorr].second;
auto& linfo = corr.getSectorRowInfo(sector, row);
secRowInfo.updateMaxValues(linfo.getMaxValues(), scale);
secRowInfo.updateMaxValues(linfo.getMinValues(), scale);
auto& linfo = corr.getRowInfo(row);

double scaleU = secRowInfo.gridMeasured.getYscale() / linfo.gridMeasured.getYscale();
double scaleV = secRowInfo.gridMeasured.getZscale() / linfo.gridMeasured.getZscale();
double scaleRealU = secRowInfo.gridReal.getYscale() / linfo.gridReal.getYscale();
double scaleRealV = secRowInfo.gridReal.getZscale() / linfo.gridReal.getZscale();
double scaleU = rowInfo.gridMeasured.getYscale() / linfo.gridMeasured.getYscale();
double scaleV = rowInfo.gridMeasured.getZscale() / linfo.gridMeasured.getZscale();
double scaleRealU = rowInfo.gridReal.getYscale() / linfo.gridReal.getYscale();
double scaleRealV = rowInfo.gridReal.getZscale() / linfo.gridReal.getZscale();

for (int iu = 0; iu < gridU.getNumberOfKnots(); iu++) {
double u = gridU.getKnot(iu).u;
Expand All @@ -1123,7 +1095,7 @@ void TPCFastSpaceChargeCorrectionHelper::mergeCorrections(
corr.convLocalToGrid(sector, row, y, z, lu, lv, ls);
ls *= scale;
double parscale[4] = {ls, ls * scaleU, ls * scaleV, ls * ls * scaleU * scaleV};
const auto& spl = corr.getSpline(sector, row);
const auto& spl = corr.getSplineForRow(row);
spl.interpolateParametersAtU(corr.getCorrectionData(sector, row), lu, lv, P);
for (int ipar = 0, ind = 0; ipar < nKnotPar1d; ++ipar) {
for (int idim = 0; idim < 3; idim++, ind++) {
Expand All @@ -1141,7 +1113,7 @@ void TPCFastSpaceChargeCorrectionHelper::mergeCorrections(
double parscale[4] = {ls, ls * scaleRealU, ls * scaleRealV, ls * ls * scaleRealU * scaleRealV};

{ // inverse X correction
corr.getSplineInvX(sector, row).interpolateParametersAtU(corr.getCorrectionDataInvX(sector, row), lu, lv, P);
corr.getSplineInvXforRow(row).interpolateParametersAtU(corr.getCorrectionDataInvX(sector, row), lu, lv, P);
for (int ipar = 0, ind = 0; ipar < nKnotPar1d; ++ipar) {
for (int idim = 0; idim < 1; idim++, ind++) {
splineParametersInvX[knotIndex * nKnotPar1d + ind] += parscale[ipar] * P[ind];
Expand All @@ -1150,7 +1122,7 @@ void TPCFastSpaceChargeCorrectionHelper::mergeCorrections(
}

{ // inverse YZ correction
corr.getSplineInvYZ(sector, row).interpolateParametersAtU(corr.getCorrectionDataInvYZ(sector, row), lu, lv, P);
corr.getSplineInvYZforRow(row).interpolateParametersAtU(corr.getCorrectionDataInvYZ(sector, row), lu, lv, P);
for (int ipar = 0, ind = 0; ipar < nKnotPar1d; ++ipar) {
for (int idim = 0; idim < 2; idim++, ind++) {
splineParametersInvYZ[knotIndex * nKnotPar2d + ind] += parscale[ipar] * P[ind];
Expand Down
4 changes: 2 additions & 2 deletions Detectors/TPC/reconstruction/src/TPCFastTransformHelperO2.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -96,7 +96,7 @@ std::unique_ptr<TPCFastTransform> TPCFastTransformHelperO2::create(int64_t TimeS
const float t0 = 0.;
const float vDrift = 0.f;
const long int initTimeStamp = -1;
fastTransform.setCalibration1(initTimeStamp, t0, vDrift);
fastTransform.setCalibration(initTimeStamp, t0, vDrift);

fastTransform.finishConstruction();
}
Expand Down Expand Up @@ -153,7 +153,7 @@ int TPCFastTransformHelperO2::updateCalibrationImpl(T& fastTransform, int64_t Ti

const double t0 = (driftTimeOffset + elParam.getAverageShapingTime()) / elParam.ZbinWidth;

fastTransform.setCalibration1(TimeStamp, t0, vDrift);
fastTransform.setCalibration(TimeStamp, t0, vDrift);

return 0;
}
Expand Down
2 changes: 1 addition & 1 deletion GPU/TPCFastTransformation/CorrectionMapsHelper.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -122,7 +122,7 @@ class CorrectionMapsHelper
{
if (mCorrMapMShape) {
// just check for the first spline the number of knots which are 4 in case of default spline object
return mCorrMapMShape->getCorrection().getSpline(0, 0).getNumberOfKnots() == 4;
return mCorrMapMShape->getCorrection().getSplineForRow(0).getNumberOfKnots() == 4;
}
return true;
}
Expand Down
Loading