TPC: Implement ad-hoc correction for r and z in old SCD map creation (#15103)

* Implement ad-hoc correction for r and z

* fix warning: usage of abs instead of std::abs

Author: wiechula

Detectors/TPC/calibration/SpacePoints/include/SpacePoints/TrackResiduals.h

@@ -443,6 +443,12 @@ class TrackResiduals
   /// output tree
   TTree* getOutputTree() { return mTreeOut.get(); }
 
+  /// Ad-hoc radial scaling factor A/C-Side
+  void setAdhocScalingFactorX(const std::array<float, 2>& scaling) { mAdhocScalingX = scaling; }
+
+  /// Ad-hoc correction of Z/X
+  void doAdhocCorrectionZ2X(bool corr) { mDoAdhocCorrectionZ2X = corr; }
+
  private:
   std::bitset<SECTORSPERSIDE * SIDES> mInitResultsContainer{};
 
@@ -502,6 +508,8 @@ class TrackResiduals
   std::array<std::vector<VoxRes>, SECTORSPERSIDE * SIDES> mVoxelResults{};                  ///< results per sector and per voxel for 3-D distortions
   VoxRes mVoxelResultsOut{};                                                                ///< the results from mVoxelResults are copied in here to be able to stream them
   VoxRes* mVoxelResultsOutPtr{&mVoxelResultsOut};                                           ///< pointer to set the branch address to for the output
+  std::array<float, 2> mAdhocScalingX{0, 0};                                                ///< Ad-hoc radial scaling factor
+  bool mDoAdhocCorrectionZ2X{false};                                                        ///< If to do ad-hoc correction for Z/X
 
   ClassDefNV(TrackResiduals, 3);
 };

Detectors/TPC/calibration/SpacePoints/src/TrackInterpolation.cxx

@@ -820,7 +820,7 @@ void TrackInterpolation::interpolateTrack(int iSeed)
           float xv = vtx.X() * cs + vtx.Y() * sn, yv = -vtx.X() * sn + vtx.Y() * cs, zv = vtx.Z();
           auto dy = yv - trkWorkITS.getY();
           auto dz = zv - trkWorkITS.getZ();
-          if ((std::abs(dy) < param::MaxResid) && (std::abs(dz) < param::MaxResid) && (std::abs(trkWorkITS.getY()) < param::MaxY) && (std::abs(trkWorkITS.getZ()) < param::MaxZ) && abs(xv) < param::MaxVtxX) {
+          if ((std::abs(dy) < param::MaxResid) && (std::abs(dz) < param::MaxResid) && (std::abs(trkWorkITS.getY()) < param::MaxY) && (std::abs(trkWorkITS.getZ()) < param::MaxZ) && std::abs(xv) < param::MaxVtxX) {
             short compXV = static_cast<short>(xv * 0x7fff / param::MaxVtxX);
             mClRes.emplace_back(dy, dz, alpha / TMath::Pi(), trkWorkITS.getY(), trkWorkITS.getZ(), 190, -1, compXV);
             if (!gidTable[GTrackID::ITSTPC].isIndexSet()) {
@@ -1168,7 +1168,7 @@ void TrackInterpolation::extrapolateTrack(int iSeed)
           float xv = vtx.X() * cs + vtx.Y() * sn, yv = -vtx.X() * sn + vtx.Y() * cs, zv = vtx.Z();
           auto dy = yv - trkWorkITS.getY();
           auto dz = zv - trkWorkITS.getZ();
-          if ((std::abs(dy) < param::MaxResid) && (std::abs(dz) < param::MaxResid) && (std::abs(trkWorkITS.getY()) < param::MaxY) && (std::abs(trkWorkITS.getZ()) < param::MaxZ) && abs(xv) < param::MaxVtxX) {
+          if ((std::abs(dy) < param::MaxResid) && (std::abs(dz) < param::MaxResid) && (std::abs(trkWorkITS.getY()) < param::MaxY) && (std::abs(trkWorkITS.getZ()) < param::MaxZ) && std::abs(xv) < param::MaxVtxX) {
             short compXV = static_cast<short>(xv * 0x7fff / param::MaxVtxX);
             mClRes.emplace_back(dy, dz, alpha / TMath::Pi(), trkWorkITS.getY(), trkWorkITS.getZ(), 190, -1, compXV);
             if (!gidTableFull[GTrackID::ITSTPC].isIndexSet()) {

Detectors/TPC/calibration/SpacePoints/src/TrackResiduals.cxx

@@ -719,8 +719,26 @@ void TrackResiduals::smooth(int iSec)
         if (!(resVox.flags & SmoothDone)) {
           continue;
         }
-        resVox.DS[ResZ] += resVox.stat[VoxZ] * resVox.DS[ResX]; // remove slope*dX contribution from dZ
-        resVox.D[ResZ] += resVox.stat[VoxZ] * resVox.DS[ResX];  // remove slope*dX contribution from dZ
+        // TODO: Usage of Z/X is bug???
+        float z2x = resVox.stat[VoxZ];
+        if (mDoAdhocCorrectionZ2X) {
+          //
+          const float z = z2x * resVox.stat[VoxX] - resVox.DS[ResZ];
+          const float x = resVox.stat[VoxX] - resVox.DS[ResX]; // is subration of DS[ResX] correct?
+          z2x = z / x;
+        }
+        resVox.DS[ResZ] += z2x * resVox.DS[ResX]; // remove slope*dX contribution from dZ
+        resVox.D[ResZ] += z2x * resVox.DS[ResX];  // remove slope*dX contribution from dZ
+                                                  //
+        if (mAdhocScalingX[iSec >= 18] != 0) {
+          const float aDX = resVox.DS[ResX] * mAdhocScalingX[iSec >= 18];
+          resVox.D[ResX] += aDX;
+          resVox.DS[ResX] += aDX;
+          resVox.D[ResY] += aDX * resVox.stat[VoxF];
+          resVox.DS[ResY] += aDX * resVox.stat[VoxF];
+          resVox.D[ResZ] += aDX * z2x;
+          resVox.DS[ResZ] += aDX * z2x;
+        }
       }
     }
   }