From dec1e31e22b045ac9539b363760ffc637ac7746d Mon Sep 17 00:00:00 2001 From: Sergey Gorbunov Date: Thu, 20 Jun 2024 15:15:41 +0000 Subject: [PATCH] TPC Splines: add limits for SP correction values per TPC row --- .../TPCFastSpaceChargeCorrectionHelper.cxx | 18 +++- .../TPCFastSpaceChargeCorrection.h | 85 ++++++++++++++++--- .../macro/TPCFastTransformInit.C | 15 ++-- 3 files changed, 97 insertions(+), 21 deletions(-) diff --git a/Detectors/TPC/calibration/src/TPCFastSpaceChargeCorrectionHelper.cxx b/Detectors/TPC/calibration/src/TPCFastSpaceChargeCorrectionHelper.cxx index ce0954120281f..e71340a555227 100644 --- a/Detectors/TPC/calibration/src/TPCFastSpaceChargeCorrectionHelper.cxx +++ b/Detectors/TPC/calibration/src/TPCFastSpaceChargeCorrectionHelper.cxx @@ -144,6 +144,9 @@ void TPCFastSpaceChargeCorrectionHelper::fillSpaceChargeCorrectionFromMap(TPCFas float* splineParameters = correction.getSplineData(slice, row); const std::vector& data = mCorrectionMap.getPoints(slice, row); int nDataPoints = data.size(); + auto& info = correction.getSliceRowInfo(slice, row); + info.resetMaxValues(); + info.resetMaxValuesInv(); if (nDataPoints >= 4) { std::vector pointSU(nDataPoints); std::vector pointSV(nDataPoints); @@ -156,6 +159,8 @@ void TPCFastSpaceChargeCorrectionHelper::fillSpaceChargeCorrectionFromMap(TPCFas pointCorr[3 * i + 0] = dx; pointCorr[3 * i + 1] = du; pointCorr[3 * i + 2] = dv; + info.updateMaxValues(2. * dx, 2. * du, 2. * dv); + info.updateMaxValuesInv(-2. * dx, -2. * du, -2. * dv); } helper.approximateDataPoints(spline, splineParameters, 0., spline.getGridX1().getUmax(), 0., spline.getGridX2().getUmax(), &pointSU[0], &pointSV[0], &pointCorr[0], nDataPoints); @@ -767,9 +772,20 @@ std::unique_ptr TPCFastSpaceChargeCorrect double yStep = (yLast - yFirst) / 2; + double zFirst = z - dz / 2.; + double zLast = z + dz / 2.; + double zStep = (zLast - zFirst) / 2.; + + if (0) { // no smoothing + yFirst = y; + yLast = y; + zFirst = z; + zLast = z; + } + for (double py = yFirst; py <= yLast + yStep / 2.; py += yStep) { - for (double pz = z - dz / 2.; pz <= z + dz / 2. + 1.e-4; pz += dz / 2.) { + for (double pz = zFirst; pz <= zLast + zStep / 2.; pz += zStep) { map.addCorrectionPoint(iRoc, iRow, py, pz, correctionX, correctionY, correctionZ); } diff --git a/GPU/TPCFastTransformation/TPCFastSpaceChargeCorrection.h b/GPU/TPCFastTransformation/TPCFastSpaceChargeCorrection.h index ba9494f24a8ee..8a9a470e81888 100644 --- a/GPU/TPCFastTransformation/TPCFastSpaceChargeCorrection.h +++ b/GPU/TPCFastTransformation/TPCFastSpaceChargeCorrection.h @@ -21,6 +21,7 @@ #include "TPCFastTransformGeo.h" #include "FlatObject.h" #include "GPUCommonDef.h" +#include "GPUCommonMath.h" namespace GPUCA_NAMESPACE { @@ -61,15 +62,64 @@ class TPCFastSpaceChargeCorrection : public FlatObject }; struct SliceRowInfo { - float gridU0{0.f}; //< U coordinate of the U-grid start - float scaleUtoGrid{0.f}; //< scale U to U-grid coordinate - float gridV0{0.f}; ///< V coordinate of the V-grid start - float scaleVtoGrid{0.f}; //< scale V to V-grid coordinate - float gridCorrU0{0.f}; ///< U coordinate of the U-grid start for corrected U - float scaleCorrUtoGrid{0.f}; ///< scale corrected U to U-grid coordinate - float gridCorrV0{0.f}; ///< V coordinate of the V-grid start for corrected V - float scaleCorrVtoGrid{0.f}; ///< scale corrected V to V-grid coordinate + float gridU0{0.f}; //< U coordinate of the U-grid start + float scaleUtoGrid{0.f}; //< scale U to U-grid coordinate + float gridV0{0.f}; ///< V coordinate of the V-grid start + float scaleVtoGrid{0.f}; //< scale V to V-grid coordinate + float gridCorrU0{0.f}; ///< U coordinate of the U-grid start for corrected U + float scaleCorrUtoGrid{0.f}; ///< scale corrected U to U-grid coordinate + float gridCorrV0{0.f}; ///< V coordinate of the V-grid start for corrected V + float scaleCorrVtoGrid{0.f}; ///< scale corrected V to V-grid coordinate + float maxCorr[3]{10.f, 10.f, 10.f}; ///< max correction for dX, dU, dV + float minCorr[3]{-10.f, -10.f, -10.f}; ///< min correction for dX, dU, dV + float maxInvCorr[3]{10.f, 10.f, 10.f}; ///< max inverse correction for dX, dU, dV + float minInvCorr[3]{-10.f, -10.f, -10.f}; ///< min inverse correction for dX, dU, dV RowActiveArea activeArea; + + void resetMaxValues() + { + maxCorr[0] = 1.f; + minCorr[0] = -1.f; + maxCorr[1] = 1.f; + minCorr[1] = -1.f; + maxCorr[2] = 1.f; + minCorr[2] = -1.f; + } + + void updateMaxValues(float dx, float du, float dv) + { + maxCorr[0] = GPUCommonMath::Max(maxCorr[0], dx); + minCorr[0] = GPUCommonMath::Min(minCorr[0], dx); + + maxCorr[1] = GPUCommonMath::Max(maxCorr[1], du); + minCorr[1] = GPUCommonMath::Min(minCorr[1], du); + + maxCorr[2] = GPUCommonMath::Max(maxCorr[2], dv); + minCorr[2] = GPUCommonMath::Min(minCorr[2], dv); + } + + void resetMaxValuesInv() + { + maxInvCorr[0] = 1.f; + minInvCorr[0] = -1.f; + maxInvCorr[1] = 1.f; + minInvCorr[1] = -1.f; + maxInvCorr[2] = 1.f; + minInvCorr[2] = -1.f; + } + + void updateMaxValuesInv(float dx, float du, float dv) + { + maxInvCorr[0] = GPUCommonMath::Max(maxInvCorr[0], dx); + minInvCorr[0] = GPUCommonMath::Min(minInvCorr[0], dx); + + maxInvCorr[1] = GPUCommonMath::Max(maxInvCorr[1], du); + minInvCorr[1] = GPUCommonMath::Min(minInvCorr[1], du); + + maxInvCorr[2] = GPUCommonMath::Max(maxInvCorr[2], dv); + minInvCorr[2] = GPUCommonMath::Min(minInvCorr[2], dv); + } + #ifndef GPUCA_ALIROOT_LIB ClassDefNV(SliceRowInfo, 2); #endif @@ -406,9 +456,10 @@ GPUdi() int TPCFastSpaceChargeCorrection::getCorrection(int slice, int row, floa convUVtoGrid(slice, row, u, v, gridU, gridV); float dxuv[3]; spline.interpolateU(splineData, gridU, gridV, dxuv); - dx = dxuv[0]; - du = dxuv[1]; - dv = dxuv[2]; + const auto& info = getSliceRowInfo(slice, row); + dx = GPUCommonMath::Max(info.minCorr[0], GPUCommonMath::Min(info.maxCorr[0], dxuv[0])); + du = GPUCommonMath::Max(info.minCorr[1], GPUCommonMath::Min(info.maxCorr[1], dxuv[1])); + dv = GPUCommonMath::Max(info.minCorr[2], GPUCommonMath::Min(info.maxCorr[2], dxuv[2])); return 0; } @@ -420,9 +471,10 @@ GPUdi() int TPCFastSpaceChargeCorrection::getCorrectionOld(int slice, int row, f convUVtoGrid(slice, row, u, v, gridU, gridV); float dxuv[3]; spline.interpolateUold(splineData, gridU, gridV, dxuv); - dx = dxuv[0]; - du = dxuv[1]; - dv = dxuv[2]; + const auto& info = getSliceRowInfo(slice, row); + dx = GPUCommonMath::Max(info.minCorr[0], GPUCommonMath::Min(info.maxCorr[0], dxuv[0])); + du = GPUCommonMath::Max(info.minCorr[1], GPUCommonMath::Min(info.maxCorr[1], dxuv[1])); + dv = GPUCommonMath::Max(info.minCorr[2], GPUCommonMath::Min(info.maxCorr[2], dxuv[2])); return 0; } @@ -436,6 +488,8 @@ GPUdi() void TPCFastSpaceChargeCorrection::getCorrectionInvCorrectedX( const float* splineData = getSplineData(slice, row, 1); float dx = 0; spline.interpolateU(splineData, gridU, gridV, &dx); + const auto& info = getSliceRowInfo(slice, row); + dx = GPUCommonMath::Max(info.minInvCorr[0], GPUCommonMath::Min(info.maxInvCorr[0], dx)); x = mGeo.getRowInfo(row).x + dx; } @@ -450,6 +504,9 @@ GPUdi() void TPCFastSpaceChargeCorrection::getCorrectionInvUV( float duv[2]; spline.interpolateU(splineData, gridU, gridV, duv); + const auto& info = getSliceRowInfo(slice, row); + duv[0] = GPUCommonMath::Max(info.minInvCorr[1], GPUCommonMath::Min(info.maxInvCorr[1], duv[0])); + duv[1] = GPUCommonMath::Max(info.minInvCorr[2], GPUCommonMath::Min(info.maxInvCorr[2], duv[1])); nomU = corrU - duv[0]; nomV = corrV - duv[1]; } diff --git a/GPU/TPCFastTransformation/macro/TPCFastTransformInit.C b/GPU/TPCFastTransformation/macro/TPCFastTransformInit.C index db2a06d5b9a73..cd6b51d631be8 100644 --- a/GPU/TPCFastTransformation/macro/TPCFastTransformInit.C +++ b/GPU/TPCFastTransformation/macro/TPCFastTransformInit.C @@ -21,8 +21,6 @@ /// root -l TPCFastTransformInit.C'("debugVoxRes.root")' /// -#include "Algorithm/RangeTokenizer.h" - #if !defined(__CLING__) || defined(__ROOTCLING__) #include @@ -41,6 +39,8 @@ #include "TPCCalibration/TPCFastSpaceChargeCorrectionHelper.h" #endif +#include "Algorithm/RangeTokenizer.h" + using namespace o2::tpc; using namespace o2::gpu; @@ -99,8 +99,9 @@ void TPCFastTransformInit(const char* fileName = "debugVoxRes.root", trackResiduals.setZ2XBinning(z2xBins); trackResiduals.init(); - { - std::cout << "input track residuals: " << std::endl; + { // debug output + + std::cout << " ===== input track residuals ==== " << std::endl; std::cout << "voxel tree y2xBins: " << y2xBins.size() << std::endl; for (auto y2x : y2xBins) { @@ -127,6 +128,7 @@ void TPCFastTransformInit(const char* fileName = "debugVoxRes.root", for (int i = 0; i < nZ2Xbins; i++) { std::cout << "getZ2X(bin) : " << trackResiduals.getZ2X(i) << std::endl; } + std::cout << " ==================================== " << std::endl; } std::cout << "create fast transformation ... " << std::endl; @@ -310,6 +312,7 @@ void TPCFastTransformInit(const char* fileName = "debugVoxRes.root", geo.convUVtoLocal(iRoc, u + cu, v + cv, cy, cz); cy -= y; cz -= z; + double d[3] = {cx - correctionX, cy - correctionY, cz - correctionZ}; if (voxEntries >= 1.) { for (int i = 0; i < 3; i++) { @@ -317,8 +320,8 @@ void TPCFastTransformInit(const char* fileName = "debugVoxRes.root", maxDiff[i] = d[i]; maxDiffRoc[i] = iRoc; maxDiffRow[i] = iRow; - std::cout << " roc " << iRoc << " row " << iRow << " xyz " << i - << " diff " << d[i] << " entries " << voxEntries << " y " << y2xBin << " z " << z2xBin << std::endl; + // std::cout << " roc " << iRoc << " row " << iRow << " xyz " << i + // << " diff " << d[i] << " entries " << voxEntries << " y " << y2xBin << " z " << z2xBin << std::endl; } sumDiff[i] += d[i] * d[i]; }