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);