TPC SCD: explicit spline knots for createFromTrackResiduals

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<double>& knots) { mY2XKnots = knots; }
+
+  /// get the y/x spline knots of createFromTrackResiduals (empty: the default knots)
+  const std::vector<double>& 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<double>& knots) { mZ2XKnots = knots; }
+
+  /// get the z/x spline knots of createFromTrackResiduals (empty: the default knots)
+  const std::vector<double>& 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<double> mY2XKnots;                        ///< y/x spline knots of createFromTrackResiduals (empty: the default)
+  std::vector<double> mZ2XKnots;                        ///< z/x spline knots of createFromTrackResiduals (empty: the default)
   TPCFastTransformGeo mGeo;                             ///< geometry parameters
 
   TPCFastSpaceChargeCorrectionMap mCorrectionMap{0, 0};

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<o2::gpu::TPCFastSpaceChargeCorrection> 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<double> binCentres[2];
+  std::vector<double> 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<double>* 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<double> 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<int> knotsInt[2];
 

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<double>{} : o2::RangeTokenizer::tokenize<double>(knots);
+  };
+  corrHelper->setY2XKnots(getKnots(y2xKnots, "y2xKnots"));
+  corrHelper->setZ2XKnots(getKnots(z2xKnots, "z2xKnots"));
   // corrHelper->setNthreads(1);
 
   if (debugMirrorAdata2C) {

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<double>{} : o2::RangeTokenizer::tokenize<double>(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);