From 535d626c369f30602f3c1695f01139d27415e33f Mon Sep 17 00:00:00 2001 From: Matthias Kleiner Date: Wed, 2 Oct 2024 10:49:08 +0200 Subject: [PATCH] TPC: adding dEdx space-charge correction - adding check for shared clusters - adding cluster occupancy --- Detectors/TPC/calibration/CMakeLists.txt | 4 +- .../include/TPCCalibration/CalculatedEdx.h | 20 +- .../TPCCalibration/CorrectdEdxDistortions.h | 101 ++++++++++ .../TPC/calibration/src/CalculatedEdx.cxx | 63 ++++++- .../src/CorrectdEdxDistortions.cxx | 173 ++++++++++++++++++ 5 files changed, 357 insertions(+), 4 deletions(-) create mode 100644 Detectors/TPC/calibration/include/TPCCalibration/CorrectdEdxDistortions.h create mode 100644 Detectors/TPC/calibration/src/CorrectdEdxDistortions.cxx diff --git a/Detectors/TPC/calibration/CMakeLists.txt b/Detectors/TPC/calibration/CMakeLists.txt index 4e569774c86a5..0ec62e5f323b3 100644 --- a/Detectors/TPC/calibration/CMakeLists.txt +++ b/Detectors/TPC/calibration/CMakeLists.txt @@ -56,6 +56,7 @@ o2_add_library(TPCCalibration src/CorrMapParam.cxx src/TPCMShapeCorrection.cxx src/DigitAdd.cxx + src/CorrectdEdxDistortions.cxx PUBLIC_LINK_LIBRARIES O2::DataFormatsTPC O2::TPCBase O2::TPCReconstruction ROOT::Minuit Microsoft.GSL::GSL @@ -111,7 +112,8 @@ o2_target_root_dictionary(TPCCalibration include/TPCCalibration/TPCScaler.h include/TPCCalibration/CorrMapParam.h include/TPCCalibration/TPCMShapeCorrection.h - include/TPCCalibration/DigitAdd.h) + include/TPCCalibration/DigitAdd.h + include/TPCCalibration/CorrectdEdxDistortions.h) o2_add_test_root_macro(macro/comparePedestalsAndNoise.C PUBLIC_LINK_LIBRARIES O2::TPCBase diff --git a/Detectors/TPC/calibration/include/TPCCalibration/CalculatedEdx.h b/Detectors/TPC/calibration/include/TPCCalibration/CalculatedEdx.h index e5889cc302b38..701c840e060eb 100644 --- a/Detectors/TPC/calibration/include/TPCCalibration/CalculatedEdx.h +++ b/Detectors/TPC/calibration/include/TPCCalibration/CalculatedEdx.h @@ -24,6 +24,7 @@ #include "CalibdEdxContainer.h" #include "CorrectionMapsHelper.h" #include "CommonUtils/TreeStreamRedirector.h" +#include "TPCCalibration/CorrectdEdxDistortions.h" #include @@ -55,6 +56,7 @@ enum class CorrectionFlags : unsigned short { GainFull = 1 << 2, ///< flag for full gain map from calibration container GainResidual = 1 << 3, ///< flag for residuals gain map from calibration container dEdxResidual = 1 << 4, ///< flag for residual dEdx correction + dEdxSC = 1 << 5, ///< flag for space-charge dEdx correction }; enum class ClusterFlags : unsigned short { @@ -64,6 +66,7 @@ enum class ClusterFlags : unsigned short { ExcludeEdgeCl = 1 << 2, ///< flag to exclude sector edge clusters in dEdx calculation ExcludeSubthresholdCl = 1 << 3, ///< flag to exclude subthreshold clusters in dEdx calculation ExcludeSectorBoundaries = 1 << 4, ///< flag to exclude sector boundary clusters in subthreshold cluster treatment + ExcludeSharedCl = 1 << 5, ///< flag to exclude clusters shared between tracks }; inline CorrectionFlags operator&(CorrectionFlags a, CorrectionFlags b) { return static_cast(static_cast(a) & static_cast(b)); } @@ -111,6 +114,9 @@ class CalculatedEdx /// set the debug streamer void setStreamer(const char* debugRootFile) { mStreamer = std::make_unique(debugRootFile, "recreate"); }; + /// set the debug streamer of the space-charge dedx correction + void setSCStreamer(const char* debugRootFile = "debug_sc_corrections.root") { mSCdEdxCorrection.setStreamer(debugRootFile); } + /// \return returns magnetic field in kG float getFieldNominalGPUBz() { return mFieldNominalGPUBz; } @@ -160,7 +166,8 @@ class CalculatedEdx /// load calibration objects from CCDB /// \param runNumberOrTimeStamp run number or time stamp - void loadCalibsFromCCDB(long runNumberOrTimeStamp); + /// \param isMC set if dEdx space-charge corrections will be loaded for MC or real data + void loadCalibsFromCCDB(long runNumberOrTimeStamp, const bool isMC = false); /// load calibration objects from local CCDB folder /// \param localCCDBFolder local CCDB folder @@ -208,6 +215,15 @@ class CalculatedEdx /// \param object name of the object to load void setPropagatorFromFile(const char* folder, const char* file, const char* object); + /// \param lumi set luminosity for space-charge correction map scaling + void setLumi(const float lumi) { mSCdEdxCorrection.setLumi(lumi); } + + /// \return returns space-charge dedx correctin + auto& getSCCorrection() { return mSCdEdxCorrection; } + + /// \return returns cluster occupancy for given cluster + unsigned int getOccupancy(const o2::tpc::ClusterNative& cl) const; + private: std::vector* mTracks{nullptr}; ///< vector containing the tpc tracks which will be processed std::vector* mTPCTrackClIdxVecInput{nullptr}; ///< input vector with TPC tracks cluster indicies @@ -228,6 +244,8 @@ class CalculatedEdx std::array, 5> mChargeTotROC; std::array, 5> mChargeMaxROC; + + CorrectdEdxDistortions mSCdEdxCorrection; ///< for space-charge correction of dE/dx }; } // namespace o2::tpc diff --git a/Detectors/TPC/calibration/include/TPCCalibration/CorrectdEdxDistortions.h b/Detectors/TPC/calibration/include/TPCCalibration/CorrectdEdxDistortions.h new file mode 100644 index 0000000000000..74e4e58fb5bf8 --- /dev/null +++ b/Detectors/TPC/calibration/include/TPCCalibration/CorrectdEdxDistortions.h @@ -0,0 +1,101 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +/// +/// @file CorrectdEdxDistortions.h +/// @author Matthias Kleiner, mkleiner@ikf.uni-frankfurt.de +/// + +#ifndef ALICEO2_TPC_CORRECTDEDXDISTORTIONS_H +#define ALICEO2_TPC_CORRECTDEDXDISTORTIONS_H + +#include "GPUTPCGeometry.h" +#include + +// forward declare classes +namespace o2::gpu +{ +class TPCFastTransform; +} + +namespace o2::utils +{ +class TreeStreamRedirector; +} + +namespace o2::tpc +{ + +class CorrectdEdxDistortions +{ + public: + /// Default constructor + CorrectdEdxDistortions(); + + /// Default destructor + ~CorrectdEdxDistortions(); + + /// setting the space-charge corrections from local file + /// \param scAvgFile average space-charge correction map + /// \param scDerFile derivative space-charge correction map + /// \param lumi luminosity which is used for scaling the correction map + void setSCCorrFromFile(const char* scAvgFile, const char* scDerFile, const float lumi = -1); + + /// setting the space-charge corrections + /// \param avg average space-charge correction map + /// \param der derivative space-charge correction map + /// \param lumi luminosity which is used for scaling the correction map + void setCorrectionMaps(o2::gpu::TPCFastTransform* avg, o2::gpu::TPCFastTransform* der, const float lumi); + + /// setting the space-charge corrections + /// \param avg average space-charge correction map + /// \param der derivative space-charge correction map + void setCorrectionMaps(o2::gpu::TPCFastTransform* avg, o2::gpu::TPCFastTransform* der); + + /// \param lumi luminosity which is used for scaling the correction map + void setLumi(float lumi); + + /// enable the debug streamer + void setStreamer(const char* debugRootFile = "debug_sc_corrections.root"); + + /// \return getting the correction value for the cluster charge + /// \param time true time information of the cluster + /// \param sector TPC sector of the cluster + /// \param padrow global padrow + /// \param pad pad number in the padrow + float getCorrection(const float time, unsigned char sector, unsigned char padrow, int pad) const; + + /// set minimum allowed lx correction at lower pad + void setMinLX0(const float minLX) { mLX0Min = minLX; } + + /// set minimum allowed lx correction upper lower pad + void setMinLX1(const float minLX) { mLX1Min = minLX; } + + /// set minimum allowed correction value + void setMinCorr(const float minCorr) { mScCorrMin = minCorr; } + + /// set maximum allowed correction value + void setMaxCorr(const float maxCorr) { mScCorrMax = maxCorr; } + + private: + float mScaleDer = 0; + std::unique_ptr mCorrAvg{nullptr}; ///< correction map for average distortions + std::unique_ptr mCorrDer{nullptr}; ///< correction map for distortions + o2::gpu::GPUTPCGeometry mTPCGeometry; ///< geometry information of TPC + std::unique_ptr mStreamer{nullptr}; ///< debug streamer + float mLX0Min{82}; ///< minimum allowed lx position after correction at position of the bottom pad + float mLX1Min{82}; ///< minimum allowed lx position after correction at position of the top pad + float mScCorrMin{0.7}; ///< minimum allowed correction values + float mScCorrMax{1.6}; ///< maximum allowed correction values +}; +} // namespace o2::tpc + +#endif diff --git a/Detectors/TPC/calibration/src/CalculatedEdx.cxx b/Detectors/TPC/calibration/src/CalculatedEdx.cxx index 00c8a4c8aa743..9809cc454e94f 100644 --- a/Detectors/TPC/calibration/src/CalculatedEdx.cxx +++ b/Detectors/TPC/calibration/src/CalculatedEdx.cxx @@ -26,6 +26,7 @@ #include "CalibdEdxTrackTopologyPol.h" #include "DataFormatsParameters/GRPMagField.h" #include "GPUO2InterfaceUtils.h" +#include "GPUTPCGMMergedTrackHit.h" using namespace o2::tpc; @@ -111,9 +112,12 @@ void CalculatedEdx::calculatedEdx(o2::tpc::TrackTPC& track, dEdxInfo& output, fl std::vector gainResidualVector; std::vector residualCorrTotVector; std::vector residualCorrMaxVector; + std::vector scCorrVector; std::vector trackVector; std::vector clVector; + std::vector occupancyVector; + std::vector isClusterShared; if (mDebug) { excludeClVector.reserve(nClusters); @@ -134,6 +138,9 @@ void CalculatedEdx::calculatedEdx(o2::tpc::TrackTPC& track, dEdxInfo& output, fl residualCorrMaxVector.reserve(nClusters); trackVector.reserve(nClusters); clVector.reserve(nClusters); + scCorrVector.reserve(nClusters); + occupancyVector.reserve(nClusters); + isClusterShared.reserve(nClusters); } // for missing clusters @@ -154,6 +161,10 @@ void CalculatedEdx::calculatedEdx(o2::tpc::TrackTPC& track, dEdxInfo& output, fl // set sectorIndex, rowIndex, clusterIndexNumb track.getClusterReference(*mTPCTrackClIdxVecInput, iCl, sectorIndex, rowIndex, clusterIndexNumb); + // check if the cluster is shared + const unsigned int absoluteIndex = mClusterIndex->clusterOffset[sectorIndex][rowIndex] + clusterIndexNumb; + const bool isShared = mRefit ? (mTPCRefitterShMap[absoluteIndex] & GPUCA_NAMESPACE::gpu::GPUTPCGMMergedTrackHit::flagShared) : 0; + // get region, pad, stack and stack ID const int region = Mapper::REGION[rowIndex]; const unsigned char pad = std::clamp(static_cast(cl.getPad() + 0.5f), static_cast(0), Mapper::PADSPERROW[region][Mapper::getLocalRowFromGlobalRow(rowIndex)] - 1); // the left side of the pad is defined at e.g. 3.5 and the right side at 4.5 @@ -179,6 +190,9 @@ void CalculatedEdx::calculatedEdx(o2::tpc::TrackTPC& track, dEdxInfo& output, fl if (((clusterMask & ClusterFlags::ExcludeEdgeCl) == ClusterFlags::ExcludeEdgeCl) && ((flagsCl & ClusterNative::flagEdge) == ClusterNative::flagEdge)) { excludeCl += 0b100; // 4 for edge cluster } + if (((clusterMask & ClusterFlags::ExcludeSharedCl) == ClusterFlags::ExcludeSharedCl) && isShared) { + excludeCl += 0b10000; // for shared cluster + } // get the x position of the track const float xPosition = Mapper::instance().getPadCentre(PadPos(rowIndex, 0)).X(); @@ -214,6 +228,8 @@ void CalculatedEdx::calculatedEdx(o2::tpc::TrackTPC& track, dEdxInfo& output, fl offsPadVector.emplace_back(offsPad); trackVector.emplace_back(track); clVector.emplace_back(cl); + occupancyVector.emplace_back(getOccupancy(cl)); + isClusterShared.emplace_back(isShared); topologyCorrVector.emplace_back(-999.f); topologyCorrTotVector.emplace_back(-999.f); @@ -222,6 +238,7 @@ void CalculatedEdx::calculatedEdx(o2::tpc::TrackTPC& track, dEdxInfo& output, fl gainResidualVector.emplace_back(-999.f); residualCorrTotVector.emplace_back(-999.f); residualCorrMaxVector.emplace_back(-999.f); + scCorrVector.emplace_back(-999.f); } // to avoid counting the skipped cluster as a subthreshold cluster rowIndexOld = rowIndex; @@ -325,6 +342,19 @@ void CalculatedEdx::calculatedEdx(o2::tpc::TrackTPC& track, dEdxInfo& output, fl minChargeMax = chargeMax; }; + // space-charge dEdx corrections + const float time = cl.getTime() - track.getTime0(); // ToDo: get correct time from ITS-TPC track if possible + float scCorr = 1.0f; + if ((correctionMask & CorrectionFlags::dEdxSC) == CorrectionFlags::dEdxSC) { + scCorr = mSCdEdxCorrection.getCorrection(time, sectorIndex, rowIndex, pad); + if (scCorr > 0) { + chargeTot /= scCorr; + }; + if (corrMax > 0) { + chargeMax /= scCorr; + }; + } + if (stack == GEMstack::IROCgem) { mChargeTotROC[0].emplace_back(chargeTot); mChargeMaxROC[0].emplace_back(chargeMax); @@ -359,6 +389,8 @@ void CalculatedEdx::calculatedEdx(o2::tpc::TrackTPC& track, dEdxInfo& output, fl offsPadVector.emplace_back(offsPad); trackVector.emplace_back(track); clVector.emplace_back(cl); + occupancyVector.emplace_back(getOccupancy(cl)); + isClusterShared.emplace_back(isShared); topologyCorrVector.emplace_back(effectiveLength); topologyCorrTotVector.emplace_back(effectiveLengthTot); @@ -367,6 +399,7 @@ void CalculatedEdx::calculatedEdx(o2::tpc::TrackTPC& track, dEdxInfo& output, fl gainResidualVector.emplace_back(gainResidual); residualCorrTotVector.emplace_back(corrTot); residualCorrMaxVector.emplace_back(corrMax); + scCorrVector.emplace_back(scCorr); }; } @@ -394,6 +427,10 @@ void CalculatedEdx::calculatedEdx(o2::tpc::TrackTPC& track, dEdxInfo& output, fl output.NHitsOROC3 = nClsROC[3]; } + // copy corrected cluster charges + auto chargeTotVector = mChargeTotROC[4]; + auto chargeMaxVector = mChargeMaxROC[4]; + // calculate dEdx output.dEdxTotIROC = getTruncMean(mChargeTotROC[0], low, high); output.dEdxTotOROC1 = getTruncMean(mChargeTotROC[1], low, high); @@ -428,6 +465,7 @@ void CalculatedEdx::calculatedEdx(o2::tpc::TrackTPC& track, dEdxInfo& output, fl << "gainResidualVector=" << gainResidualVector << "residualCorrTotVector=" << residualCorrTotVector << "residualCorrMaxVector=" << residualCorrMaxVector + << "scCorrVector=" << scCorrVector << "localXVector=" << localXVector << "localYVector=" << localYVector << "offsPadVector=" << offsPadVector @@ -436,6 +474,10 @@ void CalculatedEdx::calculatedEdx(o2::tpc::TrackTPC& track, dEdxInfo& output, fl << "minChargeTot=" << minChargeTot << "minChargeMax=" << minChargeMax << "output=" << output + << "occupancy=" << occupancyVector + << "chargeTotVector=" << chargeTotVector + << "chargeMaxVector=" << chargeMaxVector + << "isClusterShared=" << isClusterShared << "\n"; } } @@ -492,7 +534,7 @@ float CalculatedEdx::getTrackTopologyCorrectionPol(const o2::tpc::TrackTPC& trac return effectiveLength; } -void CalculatedEdx::loadCalibsFromCCDB(long runNumberOrTimeStamp) +void CalculatedEdx::loadCalibsFromCCDB(long runNumberOrTimeStamp, const bool isMC) { // setup CCDB manager auto& cm = o2::ccdb::BasicCCDBManager::instance(); @@ -542,6 +584,15 @@ void CalculatedEdx::loadCalibsFromCCDB(long runNumberOrTimeStamp) auto propagator = o2::base::Propagator::Instance(); const o2::base::MatLayerCylSet* matLut = o2::base::MatLayerCylSet::rectifyPtrFromFile(cm.get("GLO/Param/MatLUT")); propagator->setMatLUT(matLut); + + // load sc correction maps + auto avgMap = isMC ? cm.getForTimeStamp(o2::tpc::CDBTypeMap.at(o2::tpc::CDBType::CalCorrMapMC), tRun) : cm.getForTimeStamp(o2::tpc::CDBTypeMap.at(o2::tpc::CDBType::CalCorrMap), tRun); + avgMap->rectifyAfterReadingFromFile(); + + auto derMap = isMC ? cm.getForTimeStamp(o2::tpc::CDBTypeMap.at(o2::tpc::CDBType::CalCorrDerivMapMC), tRun) : cm.getForTimeStamp(o2::tpc::CDBTypeMap.at(o2::tpc::CDBType::CalCorrDerivMap), tRun); + derMap->rectifyAfterReadingFromFile(); + + mSCdEdxCorrection.setCorrectionMaps(avgMap, derMap); } void CalculatedEdx::loadCalibsFromLocalCCDBFolder(const char* localCCDBFolder) @@ -624,4 +675,12 @@ void CalculatedEdx::setPropagatorFromFile(const char* folder, const char* file, o2::base::MatLayerCylSet* matLut = o2::base::MatLayerCylSet::rectifyPtrFromFile((o2::base::MatLayerCylSet*)matLutFile->Get(object)); propagator->setMatLUT(matLut); } -} \ No newline at end of file +} + +unsigned int CalculatedEdx::getOccupancy(const o2::tpc::ClusterNative& cl) const +{ + const int nTimeBinsPerOccupBin = 16; + const int iBinOcc = cl.getTime() / nTimeBinsPerOccupBin + 2; + const unsigned int occupancy = mTPCRefitterOccMap.empty() ? -1 : mTPCRefitterOccMap[iBinOcc]; + return occupancy; +} diff --git a/Detectors/TPC/calibration/src/CorrectdEdxDistortions.cxx b/Detectors/TPC/calibration/src/CorrectdEdxDistortions.cxx new file mode 100644 index 0000000000000..f6e4fe9d64664 --- /dev/null +++ b/Detectors/TPC/calibration/src/CorrectdEdxDistortions.cxx @@ -0,0 +1,173 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +/// +/// @file CorrectdEdxDistortions.cxx +/// @author Matthias Kleiner, mkleiner@ikf.uni-frankfurt.de +/// + +#include "TPCCalibration/CorrectdEdxDistortions.h" +#include "TPCFastTransform.h" + +o2::tpc::CorrectdEdxDistortions::~CorrectdEdxDistortions() = default; + +o2::tpc::CorrectdEdxDistortions::CorrectdEdxDistortions() = default; + +void o2::tpc::CorrectdEdxDistortions::setStreamer(const char* debugRootFile) +{ + mStreamer = std::make_unique(debugRootFile, "recreate"); +}; + +void o2::tpc::CorrectdEdxDistortions::setSCCorrFromFile(const char* scAvgFile, const char* scDerFile, const float lumi) +{ + auto avg = o2::gpu::TPCFastTransform::loadFromFile(scAvgFile, "ccdb_object"); + auto der = o2::gpu::TPCFastTransform::loadFromFile(scDerFile, "ccdb_object"); + if (!avg || !der) { + LOGP(warn, "Couldnt load all sc correction objects from local file"); + return; + } + avg->rectifyAfterReadingFromFile(); + der->rectifyAfterReadingFromFile(); + setCorrectionMaps(avg, der, lumi); +} + +void o2::tpc::CorrectdEdxDistortions::setCorrectionMaps(o2::gpu::TPCFastTransform* avg, o2::gpu::TPCFastTransform* der) +{ + if (!avg || !der) { + LOGP(warn, "Nullptr detected in setting the correction maps"); + return; + } + mCorrAvg = std::unique_ptr(new o2::gpu::TPCFastTransform); + mCorrAvg->cloneFromObject(*avg, nullptr); + mCorrAvg->rectifyAfterReadingFromFile(); + + mCorrDer = std::unique_ptr(new o2::gpu::TPCFastTransform); + mCorrDer->cloneFromObject(*der, nullptr); + mCorrDer->rectifyAfterReadingFromFile(); + + mCorrAvg->setApplyCorrectionOn(); + mCorrDer->setApplyCorrectionOn(); +} + +void o2::tpc::CorrectdEdxDistortions::setCorrectionMaps(o2::gpu::TPCFastTransform* avg, o2::gpu::TPCFastTransform* der, const float lumi) +{ + setCorrectionMaps(avg, der); + setLumi(lumi); +} + +void o2::tpc::CorrectdEdxDistortions::setLumi(float lumi) +{ + if (!mCorrAvg || !mCorrDer) { + LOGP(warn, "Nullptr detected in accessing the correction maps"); + return; + } + const float lumiAvg = mCorrAvg->getLumi(); + const float lumiDer = mCorrDer->getLumi(); + mScaleDer = (lumi - lumiAvg) / lumiDer; + LOGP(info, "Setting mScaleDer: {} for inst lumi: {} avg lumi: {} deriv. lumi: {}", mScaleDer, lumi, lumiAvg, lumiDer); +} + +float o2::tpc::CorrectdEdxDistortions::getCorrection(const float time, unsigned char sector, unsigned char padrow, int pad) const +{ + // + // Get the corrections at the previous and next padrow and interpolate them to the start and end position of current pad + // Calculate from corrected position the radial distortions and compare the effective length of distorted electrons with pad length + // + + // localY of current pad + const float ly = mTPCGeometry.LinearPad2Y(sector, padrow, pad); + + // get correction at "pad + 0.5*padlength" pos1 and dont extrapolate/interpolate across GEM gaps + const int row1 = ((padrow == mTPCGeometry.EndIROC() - 1) || (padrow == mTPCGeometry.EndOROC1() - 1) || (padrow == mTPCGeometry.EndOROC2() - 1)) ? padrow : std::clamp(padrow + 1, 0, GPUCA_ROW_COUNT - 1); + + float lxT_1 = 0; + float lyT_1 = 0; + float lzT_1 = 0; + mCorrAvg->Transform(sector, row1, pad, time, lxT_1, lyT_1, lzT_1, 0, mCorrDer.get(), nullptr, mScaleDer, 0, 1); + + // correct for different localY position of pads + lyT_1 += ly - mTPCGeometry.LinearPad2Y(sector, row1, pad); + + // get radius of upper pad + const float r_1_f = std::sqrt(lxT_1 * lxT_1 + lyT_1 * lyT_1); + + // get correction at "pad - 0.5*padlength" pos0 and dont extrapolate/interpolate across GEM gaps + const int row0 = ((padrow == mTPCGeometry.EndIROC()) || (padrow == mTPCGeometry.EndOROC1()) || (padrow == mTPCGeometry.EndOROC2())) ? padrow : std::clamp(padrow - 1, 0, GPUCA_ROW_COUNT - 1); + + // check if previous pad row has enough pads + const unsigned char pad0 = std::clamp(static_cast(pad), 0, mTPCGeometry.NPads(row0) - 1); + float lxT_0 = 0; + float lyT_0 = 0; + float lzT_0 = 0; + mCorrAvg->Transform(sector, row0, pad0, time, lxT_0, lyT_0, lzT_0, 0, mCorrDer.get(), nullptr, mScaleDer, 0, 1); + + // correct for different localY position of pads + lyT_0 += ly - mTPCGeometry.LinearPad2Y(sector, row0, pad0); + + // get radius of lower pad + const float r_0_f = std::sqrt(lxT_0 * lxT_0 + lyT_0 * lyT_0); + + // effective radial length of electrons + const float dr_f = r_1_f - r_0_f; + + // position of upper and lower pad edge + const float x_0 = padrow - 0.5; + const float x_1 = padrow + 0.5; + + // interpolate corrections to upper and lower pad edge + const int deltaRow = (row1 - row0); + const float d_StartPad = (r_0_f * (row1 - x_0) + r_1_f * (x_0 - row0)) / deltaRow; + const float d_EndPad = (r_0_f * (row1 - x_1) + r_1_f * (x_1 - row0)) / deltaRow; + const float scCorr = (d_EndPad - d_StartPad) / mTPCGeometry.PadHeight(padrow); + + // check if corrected position is still reasonable + const bool isOk = ((lxT_1 < mLX0Min) || (lxT_0 < mLX1Min) || (scCorr < mScCorrMin) || (scCorr > mScCorrMax)) ? false : true; + + // store debug informations + if (mStreamer) { + const float lx = mTPCGeometry.Row2X(padrow); + + // original correction + float lxT = 0; + float lyT = 0; + float lzT = 0; + mCorrAvg->Transform(sector, padrow, pad, time, lxT, lyT, lzT, 0, mCorrDer.get(), nullptr, mScaleDer, 0, 1); + + (*mStreamer) << "tree" + << "sector=" << sector + << "padrow=" << padrow + << "row0=" << row0 + << "row1=" << row1 + << "pad0=" << pad0 + << "pad=" << pad + << "time=" << time + << "lx=" << lx + << "lxT=" << lxT + << "lyT=" << lyT + << "lzT=" << lzT + << "lxT_0=" << lxT_0 + << "lxT_1=" << lxT_1 + << "ly=" << ly + << "lyT_0=" << lyT_0 + << "lyT_1=" << lyT_1 + << "lzT_0=" << lzT_0 + << "lzT_1=" << lzT_1 + << "d_StartPad=" << d_StartPad + << "d_EndPad=" << d_EndPad + << "r_0_f=" << r_0_f + << "r_1_f=" << r_1_f + << "scCorr=" << scCorr + << "isOk=" << isOk + << "\n"; + } + + return isOk ? scCorr : 1; +}