TPC Splines: correct biased voxels; features for debugging

cbmsw · davidrohr · commit 7dd7c39b4cf9 · 2025-12-15T10:50:04.000+01:00
diff --git a/Detectors/TPC/calibration/include/TPCCalibration/TPCFastSpaceChargeCorrectionHelper.h b/Detectors/TPC/calibration/include/TPCCalibration/TPCFastSpaceChargeCorrectionHelper.h
@@ -90,11 +90,16 @@ class TPCFastSpaceChargeCorrectionHelper
   /// \param voxResTreeInverse TTree with inverse voxel residuals
   /// \param useSmoothed if true, use smoothed residuals
   /// \param invertSigns if true, invert the signs of the residuals
+  /// \param fitPointsDirect debug: pointer to the data used for the direct correction
+  /// \param fitPointsInverse debug: pointer to the data used for the inverse correction
   /// \return pointer to the created TPCFastSpaceChargeCorrection object
   /// \note voxel trees wont be changed. They are read as non-const because of the ROOT::TTreeProcessorMT interface
+  ///
   std::unique_ptr<o2::gpu::TPCFastSpaceChargeCorrection> createFromTrackResiduals(
     const o2::tpc::TrackResiduals& trackResiduals, TTree* voxResTree, TTree* voxResTreeInverse, //
-    bool useSmoothed, bool invertSigns);
+    bool useSmoothed, bool invertSigns,                                                         //
+    TPCFastSpaceChargeCorrectionMap* fitPointsDirect = nullptr,
+    TPCFastSpaceChargeCorrectionMap* fitPointsInverse = nullptr);
 
   /// _______________  Utilities   ________________________
 
@@ -116,10 +121,28 @@ class TPCFastSpaceChargeCorrectionHelper
   /// \param additionalCorrections vector of pairs of additional corrections and their scaling factors
   /// \param prn printout flag
   /// \return main correction merged with additional corrections
-  void MergeCorrections(
+  void mergeCorrections(
     o2::gpu::TPCFastSpaceChargeCorrection& mainCorrection, float scale,
     const std::vector<std::pair<const o2::gpu::TPCFastSpaceChargeCorrection*, float>>& additionalCorrections, bool prn);
 
+  /// how far the voxel mean is allowed to be outside of the voxel (1.1 means 10%)
+  void setVoxelMeanValidityRange(double range)
+  {
+    mVoxelMeanValidityRange = range;
+  }
+
+  double getVoxelMeanValidityRange() const { return mVoxelMeanValidityRange; }
+
+  /// debug: if true, use voxel centers instead of the fitted positions for correction
+  void setDebugUseVoxelCenters();
+
+  bool isDebugUseVoxelCenters() const { return mDebugUseVoxelCenters; }
+
+  /// debug: if true, mirror the data from the A side to the C side of the TPC
+  void setDebugMirrorAdata2C();
+
+  bool isDebugMirrorAdata2C() const { return mDebugMirrorAdata2C; }
+
  private:
   /// geometry initialization
   void initGeometry();
@@ -133,6 +156,11 @@ class TPCFastSpaceChargeCorrectionHelper
 
   TPCFastSpaceChargeCorrectionMap mCorrectionMap{0, 0};
 
+  double mVoxelMeanValidityRange{1.1}; ///< debug: how far the voxel mean is allowed to be outside of the voxel (1.1 means 10%)
+
+  bool mDebugUseVoxelCenters{false}; ///< debug: if true, use voxel centers instead of the fitted positions for correction
+  bool mDebugMirrorAdata2C{false};   ///< debug: if true, mirror the data from the A side to the C side of the TPC
+
   ClassDefNV(TPCFastSpaceChargeCorrectionHelper, 0);
 };
 
diff --git a/Detectors/TPC/calibration/src/TPCFastSpaceChargeCorrectionHelper.cxx b/Detectors/TPC/calibration/src/TPCFastSpaceChargeCorrectionHelper.cxx
@@ -33,6 +33,7 @@
 #include "TTreeReaderValue.h"
 #include "ROOT/TTreeProcessorMT.hxx"
 #include <algorithm>
+#include <sstream>
 
 using namespace o2::gpu;
 
@@ -381,7 +382,9 @@ void TPCFastSpaceChargeCorrectionHelper::testGeometry(const TPCFastTransformGeo&
 }
 
 std::unique_ptr<o2::gpu::TPCFastSpaceChargeCorrection> TPCFastSpaceChargeCorrectionHelper::createFromTrackResiduals(
-  const o2::tpc::TrackResiduals& trackResiduals, TTree* voxResTree, TTree* voxResTreeInverse, bool useSmoothed, bool invertSigns)
+  const o2::tpc::TrackResiduals& trackResiduals, TTree* voxResTree, TTree* voxResTreeInverse, bool useSmoothed, bool invertSigns,
+  TPCFastSpaceChargeCorrectionMap* fitPointsDirect,
+  TPCFastSpaceChargeCorrectionMap* fitPointsInverse)
 {
   // create o2::gpu::TPCFastSpaceChargeCorrection  from o2::tpc::TrackResiduals::VoxRes voxel tree
 
@@ -603,6 +606,24 @@ std::unique_ptr<o2::gpu::TPCFastSpaceChargeCorrection> TPCFastSpaceChargeCorrect
       processor.Process(myThread);
     }
 
+    // debug: mirror the data for TPC C side
+
+    if (mDebugMirrorAdata2C) {
+      for (int iSector = 0; iSector < geo.getNumberOfSectorsA(); iSector++) {
+        for (int iRow = 0; iRow < nRows; iRow++) {
+          for (int iy = 0; iy < nY2Xbins; iy++) {
+            for (int iz = 0; iz < nZ2Xbins; iz++) {
+              auto& dataA = vSectorData[iSector * nRows + iRow][iy * nZ2Xbins + iz];
+              auto& dataC = vSectorData[(iSector + geo.getNumberOfSectorsA()) * nRows + iRow][iy * nZ2Xbins + iz];
+              dataC = dataA;          // copy the data
+              dataC.mZ = -dataC.mZ;   // mirror the Z coordinate
+              dataC.mCz = -dataC.mCz; // mirror the Z correction
+            }
+          }
+        }
+      }
+    }
+
     for (int iSector = 0; iSector < nSectors; iSector++) {
 
       // now process the data row-by-row
@@ -623,6 +644,7 @@ std::unique_ptr<o2::gpu::TPCFastSpaceChargeCorrection> TPCFastSpaceChargeCorrect
           {
             int xBin = iRow;
             double x = trackResiduals.getX(xBin); // radius of the pad row
+            double dx = 1. / trackResiduals.getDXI(xBin);
             bool isDataFound = false;
             for (int iy = 0; iy < nY2Xbins; iy++) {
               for (int iz = 0; iz < nZ2Xbins; iz++) {
@@ -642,13 +664,29 @@ std::unique_ptr<o2::gpu::TPCFastSpaceChargeCorrection> TPCFastSpaceChargeCorrect
                 if (data.mNentries > 0) { // voxel contains data
                   vox.mSmoothingStep = 0; // take original data
                   isDataFound = true;
-                  if (fabs(x - data.mX) > 1. || fabs(vox.mY - data.mY) > 5. || fabs(vox.mZ - data.mZ) > 5.) {
-                    std::cout << directionName << ": fitted voxel is too far from the nominal position: "
-                              << " sector " << iSector << " row " << iRow
-                              << " center x " << x << " y " << vox.mY << " z " << vox.mZ
-                              << " fitted x " << data.mX << " y " << data.mY << " z " << data.mZ
-                              << std::endl;
+
+                  // correct the mean position if it is outside the voxel
+                  std::stringstream msg;
+                  if (fabs(x - data.mX) > mVoxelMeanValidityRange * dx / 2.) {
+                    msg << "\n     x: center " << x << " dx " << data.mX - x << " half bin size: " << dx / 2;
+                  }
+
+                  if (fabs(vox.mY - data.mY) > mVoxelMeanValidityRange * vox.mDy / 2.) {
+                    msg << "\n     y: center " << vox.mY << " dy " << data.mY - vox.mY << " half bin size: " << vox.mDy / 2;
+                    data.mY = vox.mY;
+                  }
+
+                  if (fabs(vox.mZ - data.mZ) > mVoxelMeanValidityRange * vox.mDz / 2.) {
+                    msg << "\n     z: center " << vox.mZ << " dz " << data.mZ - vox.mZ << " half bin size: " << vox.mDz / 2;
+                    data.mZ = vox.mZ;
                   }
+
+                  if (!msg.str().empty()) {
+                    LOG(warning) << directionName << " correction: fitted voxel position is outside the voxel: "
+                                 << " sector " << iSector << " row " << iRow << " bin: " << iy << " " << iz
+                                 << msg.str();
+                  }
+
                 } else { // no data, take voxel center position
                   data.mCx = 0.;
                   data.mCy = 0.;
@@ -658,7 +696,7 @@ std::unique_ptr<o2::gpu::TPCFastSpaceChargeCorrection> TPCFastSpaceChargeCorrect
                   data.mZ = vox.mZ;
                   vox.mSmoothingStep = 100; // fill this data point with smoothed values from the neighbours
                 }
-                if (0) { // debug: always use voxel center instead of the mean position
+                if (mDebugUseVoxelCenters) { // debug: always use voxel center instead of the mean position
                   data.mY = vox.mY;
                   data.mZ = vox.mZ;
                 }
@@ -809,6 +847,13 @@ std::unique_ptr<o2::gpu::TPCFastSpaceChargeCorrection> TPCFastSpaceChargeCorrect
 
     TStopwatch watch4;
 
+    if (!processingInverseCorrection && fitPointsDirect) {
+      *fitPointsDirect = helper->getCorrectionMap();
+    }
+    if (processingInverseCorrection && fitPointsInverse) {
+      *fitPointsInverse = helper->getCorrectionMap();
+    }
+
     helper->fillSpaceChargeCorrectionFromMap(correction, processingInverseCorrection);
 
     LOG(info) << "fast space charge correction helper: creation from the data map took " << watch4.RealTime() << "s";
@@ -956,7 +1001,7 @@ void TPCFastSpaceChargeCorrectionHelper::initInverse(std::vector<o2::gpu::TPCFas
   LOGP(info, "Inverse tooks: {}s", duration);
 }
 
-void TPCFastSpaceChargeCorrectionHelper::MergeCorrections(
+void TPCFastSpaceChargeCorrectionHelper::mergeCorrections(
   o2::gpu::TPCFastSpaceChargeCorrection& mainCorrection, float mainScale,
   const std::vector<std::pair<const o2::gpu::TPCFastSpaceChargeCorrection*, float>>& additionalCorrections, bool /*prn*/)
 {
@@ -1099,5 +1144,17 @@ void TPCFastSpaceChargeCorrectionHelper::MergeCorrections(
   LOGP(info, "Merge of corrections tooks: {}s", duration);
 }
 
+void TPCFastSpaceChargeCorrectionHelper::setDebugUseVoxelCenters()
+{
+  LOG(info) << "fast space charge correction helper: use voxel centers for correction";
+  mDebugUseVoxelCenters = true;
+}
+
+void TPCFastSpaceChargeCorrectionHelper::setDebugMirrorAdata2C()
+{
+  LOG(info) << "fast space charge correction helper: mirror A data to C data";
+  mDebugMirrorAdata2C = true;
+}
+
 } // namespace tpc
 } // namespace o2
diff --git a/GPU/TPCFastTransformation/macro/TPCFastTransformInit.C b/GPU/TPCFastTransformation/macro/TPCFastTransformInit.C
