From fb5df5b2297a72e7808d7a51f3a4c2865765e2a2 Mon Sep 17 00:00:00 2001 From: Matthias Kleiner Date: Thu, 12 Nov 2026 15:06:18 +0100 Subject: [PATCH] TPC SCD: explicit spline knots for createFromTrackResiduals The spline knots of the FastTransform built from a voxel map were fixed to a knot at every second voxel bin, counted from both ends. With a non-uniform y/x binning (e.g. bins refined at the sector edges and at the sector centre) this skips knots in the refined bins, where the map is steepest, and with an even number of bins the two central bins get no knot at all (48 y/x bins refined at the centre: the FastTransform deviates from the map by up to ~0.4 mm there). TPCFastSpaceChargeCorrectionHelper::setY2XKnots / setZ2XKnots take the knot positions explicitly, in the coordinate of the voxel binning (y/x scaled to -1..1, z/x to 0..1, as the bin boundaries given to TrackResiduals::setY2XBinning / setZ2XBinning). A row's spline spans its first to its last voxel bin centre, so these two are always knots; knots outside them are dropped with a warning, as are knots closer than half the narrowest voxel bin to their neighbouring knot (the knots are put on an integer grid with 1/10 of the smallest knot distance as step and the spline stores one table entry per grid step, so nearly coinciding knots would make the correction object arbitrarily large). A knot interval without any voxel centre is reported. Without explicit knots (default) the previous rule is used unchanged. TPCFastTransformInit.C / TPCFastTransformInitCPM.C: optional arguments y2xKnots / z2xKnots (comma-separated); if empty, the user info "y2xKnots" / "z2xKnots" of the input voxel tree is used, so a map can carry its knots next to its binning. Co-Authored-By: Claude Opus 5.5 --- .../TPCFastSpaceChargeCorrectionHelper.h | 20 +++++ .../TPCFastSpaceChargeCorrectionHelper.cxx | 85 ++++++++++++++++--- .../macro/TPCFastTransformInit.C | 15 +++- .../macro/TPCFastTransformInitCPM.C | 17 +++- 4 files changed, 124 insertions(+), 13 deletions(-) diff --git a/Detectors/TPC/calibration/include/TPCCalibration/TPCFastSpaceChargeCorrectionHelper.h b/Detectors/TPC/calibration/include/TPCCalibration/TPCFastSpaceChargeCorrectionHelper.h index bd667954fad9a..df1ca6d67b4ad 100644 --- a/Detectors/TPC/calibration/include/TPCCalibration/TPCFastSpaceChargeCorrectionHelper.h +++ b/Detectors/TPC/calibration/include/TPCCalibration/TPCFastSpaceChargeCorrectionHelper.h @@ -73,6 +73,24 @@ class TPCFastSpaceChargeCorrectionHelper /// get the number of threads int getNthreads() const { return mNthreads; } + /// sets the y/x spline knots of createFromTrackResiduals, in the coordinate of the voxel binning (y/x scaled to -1..1, as the + /// bin boundaries of TrackResiduals::setY2XBinning). The spline of a row spans its first to its last voxel bin centre: those two are + /// always knots, knots outside them are dropped, as are knots closer than half the narrowest voxel bin to their neighbouring knot. + /// Empty (default): a knot at every second voxel bin, counted from both ends. + void setY2XKnots(const std::vector& knots) { mY2XKnots = knots; } + + /// get the y/x spline knots of createFromTrackResiduals (empty: the default knots) + const std::vector& getY2XKnots() const { return mY2XKnots; } + + /// sets the z/x spline knots of createFromTrackResiduals, in the coordinate of the voxel binning (z/x scaled to 0..1, i.e. z/x divided by + /// SpacePointsCalibConfParam::maxZ2X, as the bin boundaries of TrackResiduals::setZ2XBinning). As for y/x, the first and the last voxel + /// bin centre are always knots, knots outside them are dropped, as are knots closer than half the narrowest voxel bin to their + /// neighbouring knot. Empty (default): a knot at every second voxel bin, counted from both ends. + void setZ2XKnots(const std::vector& knots) { mZ2XKnots = knots; } + + /// get the z/x spline knots of createFromTrackResiduals (empty: the default knots) + const std::vector& getZ2XKnots() const { return mZ2XKnots; } + /// _______________ Main functionality ________________________ /// creates TPCFastSpaceChargeCorrection object from a continious space charge correction in local coordinates @@ -183,6 +201,8 @@ class TPCFastSpaceChargeCorrectionHelper static TPCFastSpaceChargeCorrectionHelper* sInstance; ///< singleton instance bool mIsInitialized = 0; ///< initialization flag int mNthreads = 1; ///< n of threads to use + std::vector mY2XKnots; ///< y/x spline knots of createFromTrackResiduals (empty: the default) + std::vector mZ2XKnots; ///< z/x spline knots of createFromTrackResiduals (empty: the default) TPCFastTransformGeo mGeo; ///< geometry parameters TPCFastSpaceChargeCorrectionMap mCorrectionMap{0, 0}; diff --git a/Detectors/TPC/calibration/src/TPCFastSpaceChargeCorrectionHelper.cxx b/Detectors/TPC/calibration/src/TPCFastSpaceChargeCorrectionHelper.cxx index 0f55aee690b4f..494d8059707ea 100644 --- a/Detectors/TPC/calibration/src/TPCFastSpaceChargeCorrectionHelper.cxx +++ b/Detectors/TPC/calibration/src/TPCFastSpaceChargeCorrectionHelper.cxx @@ -21,6 +21,7 @@ #include "TPCBase/ParameterGas.h" #include "TPCBase/Sector.h" #include "SpacePoints/TrackResiduals.h" +#include "SpacePoints/SpacePointsCalibConfParam.h" #include "DataFormatsTPC/Defs.h" #include "TPCFastTransform.h" #include "ChebyshevFit1D.h" @@ -513,21 +514,83 @@ std::unique_ptr TPCFastSpaceChargeCorrect knotsDouble[0].reserve(nY2Xbins); knotsDouble[1].reserve(nZ2Xbins); - // to get enouth measurements, make a spline knot at every second bin. Boundary bins are always included. - - for (int i = 0, j = nY2Xbins - 1; i <= j; i += 2, j -= 2) { - knotsDouble[0].push_back(trackResiduals.getY2X(0, i)); - if (j >= i + 1) { - knotsDouble[0].push_back(trackResiduals.getY2X(0, j)); - } + // the spline knots, in the coordinate of the knot calculation below (y/x of row 0, z/x): the bin centres and widths of the voxels + std::vector binCentres[2]; + std::vector binWidths[2]; + for (int i = 0; i < nY2Xbins; i++) { + binCentres[0].push_back(trackResiduals.getY2X(0, i)); + binWidths[0].push_back(1. / trackResiduals.getDY2XI(0, i)); } + for (int i = 0; i < nZ2Xbins; i++) { + binCentres[1].push_back(trackResiduals.getZ2X(i)); + binWidths[1].push_back(trackResiduals.getDZ2X(i)); + } + // explicit knots (setY2XKnots / setZ2XKnots) are given in the coordinate of the voxel binning: y/x scaled to -1..1, z/x scaled to 0..1 + const double knotScale[2] = {trackResiduals.getMaxY2X(0), o2::tpc::SpacePointsCalibConfParam::Instance().maxZ2X}; + const std::vector* explicitKnots[2] = {&mY2XKnots, &mZ2XKnots}; - for (int i = 0, j = nZ2Xbins - 1; i <= j; i += 2, j -= 2) { - knotsDouble[1].push_back(trackResiduals.getZ2X(i)); - if (j >= i + 1) { - knotsDouble[1].push_back(trackResiduals.getZ2X(j)); + for (int dim = 0; dim < 2; dim++) { + const auto& centres = binCentres[dim]; + auto& knots = knotsDouble[dim]; + if (explicitKnots[dim]->empty()) { + // default: to get enough measurements, a spline knot at every second bin, counted from both ends. Boundary bins are always included. + for (int i = 0, j = int(centres.size()) - 1; i <= j; i += 2, j -= 2) { + knots.push_back(centres[i]); + if (j >= i + 1) { + knots.push_back(centres[j]); + } + } + continue; + } + // explicit knots: the spline of a row spans its first to its last bin centre, so those are knots and the others must lie between them. + // A knot closer than half the narrowest voxel bin to its neighbouring knot is dropped: the knots are put on an integer grid with 1/10 of + // the smallest knot distance as step, and the spline stores a table with one entry per grid step, so nearly coinciding knots would make + // the correction object arbitrarily large + const double first = centres.front(); + const double last = centres.back(); + const double eps = 1.e-6 * std::abs(last - first); + const double minKnotDistance = 0.5 * *std::min_element(binWidths[dim].begin(), binWidths[dim].end()); + std::vector innerKnots; + int nOutside = 0; + for (const double k : *explicitKnots[dim]) { + const double u = k * knotScale[dim]; + if (u < first - eps || u > last + eps) { + ++nOutside; + } else { + innerKnots.push_back(u); + } + } + std::sort(innerKnots.begin(), innerKnots.end()); + int nTooClose = 0; + knots.push_back(first); + for (const double u : innerKnots) { + if (std::abs(u - knots.back()) <= eps || std::abs(u - last) <= eps) { + continue; // a repeated knot, or the first / the last bin centre given explicitly + } + if (u - knots.back() < minKnotDistance || last - u < minKnotDistance) { + ++nTooClose; + continue; + } + knots.push_back(u); + } + knots.push_back(last); + if (nOutside > 0) { + LOG(warning) << "fast space charge correction helper: " << nOutside << (dim == 0 ? " y/x" : " z/x") << " knots outside the first and the last bin centre dropped"; + } + if (nTooClose > 0) { + LOG(warning) << "fast space charge correction helper: " << nTooClose << (dim == 0 ? " y/x" : " z/x") << " knots closer than half the narrowest voxel bin (" + << minKnotDistance / knotScale[dim] << ") to the neighbouring knot dropped"; + } + // a knot interval without a voxel centre (endpoints included) leaves the spline there without data + for (size_t i = 1; i < knots.size(); i++) { + const bool covered = std::any_of(centres.begin(), centres.end(), [&](double c) { return c >= knots[i - 1] - eps && c <= knots[i] + eps; }); + if (!covered) { + LOG(warning) << "fast space charge correction helper: no " << (dim == 0 ? "y/x" : "z/x") << " voxel centre between the knots " << knots[i - 1] / knotScale[dim] << " and " << knots[i] / knotScale[dim]; + } } } + LOG(info) << "fast space charge correction helper: " << knotsDouble[0].size() << " y/x knots for " << nY2Xbins << " bins (" << (mY2XKnots.empty() ? "default" : "explicit") << "), " + << knotsDouble[1].size() << " z/x knots for " << nZ2Xbins << " bins (" << (mZ2XKnots.empty() ? "default" : "explicit") << ")"; std::vector knotsInt[2]; diff --git a/GPU/TPCFastTransformation/macro/TPCFastTransformInit.C b/GPU/TPCFastTransformation/macro/TPCFastTransformInit.C index dace124617cac..0c2ff1903ca7d 100644 --- a/GPU/TPCFastTransformation/macro/TPCFastTransformInit.C +++ b/GPU/TPCFastTransformation/macro/TPCFastTransformInit.C @@ -45,7 +45,8 @@ using namespace o2::tpc; using namespace o2::gpu; void TPCFastTransformInit(const char* fileName = "debugVoxRes.root", const char* fileNameInv = "debugVoxResInv.root", - const char* outFileName = "TPCFastTransform_VoxRes.root", bool useSmoothed = false, bool invertSigns = false, bool doDebug = true) + const char* outFileName = "TPCFastTransform_VoxRes.root", bool useSmoothed = false, bool invertSigns = false, bool doDebug = true, + const char* y2xKnots = "", const char* z2xKnots = "") { // Initialise TPCFastTransform object from "voxRes" tree of @@ -119,6 +120,18 @@ void TPCFastTransformInit(const char* fileName = "debugVoxRes.root", const char* o2::tpc::TPCFastSpaceChargeCorrectionHelper* corrHelper = o2::tpc::TPCFastSpaceChargeCorrectionHelper::instance(); corrHelper->setNthreadsToMaximum(); + // the spline knots (setY2XKnots / setZ2XKnots, comma-separated, in the coordinate of the binning): the argument, else the user info + // "y2xKnots" / "z2xKnots" of the map, else the helper's default (every second voxel bin) + auto getKnots = [userInfo](const char* argument, const char* key) { + std::string knots = argument ? argument : ""; + if (knots.empty() && userInfo->FindObject(key)) { + knots = userInfo->FindObject(key)->GetTitle(); + } + std::cout << key << ": " << (knots.empty() ? std::string("default (every second voxel bin)") : knots) << std::endl; + return knots.empty() ? std::vector{} : o2::RangeTokenizer::tokenize(knots); + }; + corrHelper->setY2XKnots(getKnots(y2xKnots, "y2xKnots")); + corrHelper->setZ2XKnots(getKnots(z2xKnots, "z2xKnots")); // corrHelper->setNthreads(1); if (debugMirrorAdata2C) { diff --git a/GPU/TPCFastTransformation/macro/TPCFastTransformInitCPM.C b/GPU/TPCFastTransformation/macro/TPCFastTransformInitCPM.C index b5de036ad4b4f..01f0037af2d2a 100644 --- a/GPU/TPCFastTransformation/macro/TPCFastTransformInitCPM.C +++ b/GPU/TPCFastTransformation/macro/TPCFastTransformInitCPM.C @@ -72,7 +72,9 @@ void TPCFastTransformInitCPM(const char* fileName = "debugVoxRes.root", float meanIDC = 0.f, float meanCTP = 0.f, int debug = 0, - int nThreads = 8) + int nThreads = 8, + const char* y2xKnots = "", + const char* z2xKnots = "") { // Initialise TPCFastTransform object from "voxRes" tree of @@ -156,6 +158,19 @@ void TPCFastTransformInitCPM(const char* fileName = "debugVoxRes.root", trackResiduals.setZ2XBinning(z2xBins); trackResiduals.init(); + // the spline knots of the FT (TPCFastSpaceChargeCorrectionHelper::setY2XKnots / setZ2XKnots, comma-separated, in the coordinate of the + // binning): the argument, else the user info "y2xKnots" / "z2xKnots" of the map, else the helper's default (every second voxel bin) + auto getKnots = [userInfo](const char* argument, const char* key) { + std::string knots = argument ? argument : ""; + if (knots.empty() && userInfo->FindObject(key)) { + knots = userInfo->FindObject(key)->GetTitle(); + } + LOGP(info, "{}: {}", key, knots.empty() ? std::string("default (every second voxel bin)") : knots); + return knots.empty() ? std::vector{} : o2::RangeTokenizer::tokenize(knots); + }; + o2::tpc::TPCFastSpaceChargeCorrectionHelper::instance()->setY2XKnots(getKnots(y2xKnots, "y2xKnots")); + o2::tpc::TPCFastSpaceChargeCorrectionHelper::instance()->setZ2XKnots(getKnots(z2xKnots, "z2xKnots")); + auto getFromUserInfo = [userInfo](std::string value, float& valueF) { if (valueF != 0) { LOGP(info, "{} set to {} via command line, not reading it from userInfo", value, valueF);