From 672f30b5236aa7fd2a7923c937db44099a75233a Mon Sep 17 00:00:00 2001 From: Matthias Kleiner Date: Sat, 13 Jan 2024 23:21:05 +0100 Subject: [PATCH] Calculating spline object in TPCScalerSpec (crashing) --- .../calibration/src/CorrectionMapsLoader.cxx | 43 +++------------ Detectors/TPC/workflow/src/TPCScalerSpec.cxx | 54 +++++++++++++++++-- .../CorrectionMapsHelper.cxx | 9 ++-- GPU/TPCFastTransformation/TPCFastTransform.h | 6 +++ GPU/Utils/FlatObject.h | 9 ++-- 5 files changed, 73 insertions(+), 48 deletions(-) diff --git a/Detectors/TPC/calibration/src/CorrectionMapsLoader.cxx b/Detectors/TPC/calibration/src/CorrectionMapsLoader.cxx index 46eeedf2244fb..f97217889bd26 100644 --- a/Detectors/TPC/calibration/src/CorrectionMapsLoader.cxx +++ b/Detectors/TPC/calibration/src/CorrectionMapsLoader.cxx @@ -69,43 +69,14 @@ void CorrectionMapsLoader::extractCCDBInputs(ProcessingContext& pc) setInstLumi(mInstLumiFactor * tpcScaler, !getUseMShapeCorrection()); } if (getUseMShapeCorrection()) { - auto mshapescaler = pc.inputs().get>("mshape"); - float potential = mshapescaler.first; - float deltaPotPos = mshapescaler.second; - // setMShapeScaler(potential); + auto mapMShape = pc.inputs().get("mshape"); + static int counter = 0; + LOGP(info, "Set MSHape map"); + setCorrMapMShape(std::move(o2::tpc::TPCFastTransformHelperO2::instance()->create(0))); + // setCorrMapMShape(const_cast(mapMShape.get())); + mCorrMapMShape->rectifyAfterReadingFromFile(); + setUpdatedMapMShape(); setMShapeScaler(1); // not needed or as some custom scaling - if (potential != 0) { - LOGP(info, "Recalculating M-shape correction map for z position of delta potential {}", deltaPotPos); - // update analytical map! - const int nthreads = 5; - using namespace o2::tpc; - o2::tpc::SpaceCharge::setNThreads(nthreads); - const Side side = Side::A; - int field = 5; - o2::tpc::SpaceCharge sc = o2::tpc::SpaceCharge(field, 129, 129, 5); - - sc.setIFCChargeUpFallingPot(potential, deltaPotPos, 0, 250, 0, side); - sc.setIFCChargeUpRisingPot(potential, deltaPotPos, 0, 0, 0, side); - - sc.poissonSolver(side); - sc.calcEField(side); - const auto numEFields = sc.getElectricFieldsInterpolator(side); - sc.calcGlobalCorrections(numEFields); - - std::function getCorrections = [&sc = sc](const int roc, double x, double y, double z, double& dx, double& dy, double& dz) { - Side side = roc < 18 ? Side::A : Side::C; - sc.getCorrections(x, y, z, side, dx, dy, dz); - }; - - TPCFastSpaceChargeCorrectionHelper::instance()->setNthreads(nthreads); - const int nKnotsY = 4; - const int nKnotsZ = 4; - std::unique_ptr spCorrection = TPCFastSpaceChargeCorrectionHelper::instance()->createFromGlobalCorrection(getCorrections, nKnotsY, nKnotsZ); - std::unique_ptr fastTransform(TPCFastTransformHelperO2::instance()->create(0, *spCorrection)); - LOGP(info, "Setting M-shape correction from calculated model"); - setCorrMapMShape(std::move(fastTransform)); - setUpdatedMapMShape(); - } } // update inverse in case it is requested if (!mScaleInverse) { diff --git a/Detectors/TPC/workflow/src/TPCScalerSpec.cxx b/Detectors/TPC/workflow/src/TPCScalerSpec.cxx index e776c0284f475..80d730a32f052 100644 --- a/Detectors/TPC/workflow/src/TPCScalerSpec.cxx +++ b/Detectors/TPC/workflow/src/TPCScalerSpec.cxx @@ -25,6 +25,8 @@ #include "TPCCalibration/TPCMShapeScaler.h" #include "CommonDataFormat/Pair.h" #include "TTree.h" +#include "TPCCalibration/TPCFastSpaceChargeCorrectionHelper.h" +#include "TPCSpaceCharge/SpaceCharge.h" using namespace o2::framework; @@ -82,11 +84,53 @@ class TPCScalerSpec : public Task if (pc.services().get().runNumber != mMShapeTPCScaler.getRun()) { LOGP(error, "Run number {} of processed data and run number {} of loaded TPC M-shape scaler doesnt match!", pc.services().get().runNumber, mMShapeTPCScaler.getRun()); } - const auto scalers = mMShapeTPCScaler.getScaler(timestamp); - float mshapescaler = mMShapeScalingFac * scalers.first; - float zPos = scalers.first; - LOGP(info, "Publishing TPC M-shape scaler: {} for timestamp {}, firstTFOrbit: {} with MShapeScalingFac: {}", mshapescaler, long(timestamp), firstTFOrbit, mMShapeScalingFac); - pc.outputs().snapshot(Output{header::gDataOriginTPC, "TPCMShapeScaler"}, dataformats::Pair{mshapescaler, zPos}); + + float potential = 0; + float deltaPotPos = 150; + // setMShapeScaler(potential); + LOGP(info, "Recalculating M-shape correction map for z position of delta potential {}", deltaPotPos); + // update analytical map! + const int nthreads = 5; + using namespace o2::tpc; + o2::tpc::SpaceCharge::setNThreads(nthreads); + const Side side = Side::A; + int field = 5; + o2::tpc::SpaceCharge sc = o2::tpc::SpaceCharge(field, 129, 129, 5); + + if (potential != 0) { + sc.setIFCChargeUpFallingPot(potential, deltaPotPos, 0, 250, 0, side); + sc.setIFCChargeUpRisingPot(potential, deltaPotPos, 0, 0, 0, side); + + sc.poissonSolver(side); + sc.calcEField(side); + const auto numEFields = sc.getElectricFieldsInterpolator(side); + sc.calcGlobalCorrections(numEFields); + } + + std::function getCorrections = [potential, &sc = sc](const int roc, double x, double y, double z, double& dx, double& dy, double& dz) { + Side side = roc < 18 ? Side::A : Side::C; + if (potential == 0 || side == Side::C) { + dx = 0; + dy = 0; + dz = 0; + } else { + sc.getCorrections(x, y, z, side, dx, dy, dz); + } + }; + + TPCFastSpaceChargeCorrectionHelper::instance()->setNthreads(nthreads); + const int nKnotsY = 4; + const int nKnotsZ = 4; + // std::unique_ptr spCorrection = TPCFastSpaceChargeCorrectionHelper::instance()->createFromGlobalCorrection(getCorrections, nKnotsY, nKnotsZ); + // std::unique_ptr fastTransform(TPCFastTransformHelperO2::instance()->create(0, *spCorrection)); + std::unique_ptr fastTransform(TPCFastTransformHelperO2::instance()->create(0)); + pc.outputs().snapshot(Output{header::gDataOriginTPC, "TPCMShapeScaler"}, *fastTransform); + + // const auto scalers = mMShapeTPCScaler.getScaler(timestamp); + // float mshapescaler = mMShapeScalingFac * scalers.first; + // float zPos = scalers.first; + // LOGP(info, "Publishing TPC M-shape scaler: {} for zMax: {} timestamp {}, firstTFOrbit: {} with MShapeScalingFac: {}", mshapescaler, zPos, long(timestamp), firstTFOrbit, mMShapeScalingFac); + // pc.outputs().snapshot(Output{header::gDataOriginTPC, "TPCMShapeScaler"}, dataformats::Pair{mshapescaler, zPos}); } } diff --git a/GPU/TPCFastTransformation/CorrectionMapsHelper.cxx b/GPU/TPCFastTransformation/CorrectionMapsHelper.cxx index a8d2cdef9be8b..06a1bd79dc6d2 100644 --- a/GPU/TPCFastTransformation/CorrectionMapsHelper.cxx +++ b/GPU/TPCFastTransformation/CorrectionMapsHelper.cxx @@ -20,8 +20,9 @@ void CorrectionMapsHelper::clear() if (mOwner) { delete mCorrMap; delete mCorrMapRef; - delete mCorrMapMShape; + // delete mCorrMapMShape; } + delete mCorrMapMShape; mCorrMap = nullptr; mCorrMapRef = nullptr; mCorrMapMShape = nullptr; @@ -89,9 +90,9 @@ void CorrectionMapsHelper::setCorrMapRef(std::unique_ptr&& m) void CorrectionMapsHelper::setCorrMapMShape(std::unique_ptr&& m) { - if (!mOwner) { - throw std::runtime_error("we must not take the ownership from a unique ptr if mOwner is not set"); - } + // if (!mOwner) { + // throw std::runtime_error("we must not take the ownership from a unique ptr if mOwner is not set"); + // } delete mCorrMapMShape; mCorrMapMShape = m.release(); } diff --git a/GPU/TPCFastTransformation/TPCFastTransform.h b/GPU/TPCFastTransformation/TPCFastTransform.h index 6bf4680b3a2c5..2a110cc61aaaf 100644 --- a/GPU/TPCFastTransformation/TPCFastTransform.h +++ b/GPU/TPCFastTransformation/TPCFastTransform.h @@ -109,6 +109,12 @@ class TPCFastTransform : public FlatObject /// Assignment operator: disabled to avoid ambiguity. Use cloneFromObject() instead TPCFastTransform& operator=(const TPCFastTransform&) CON_DELETE; + inline void destroy() + { + mCorrection.destroy(); + FlatObject::destroy(); + } + /// Destructor #if !defined(GPUCA_GPUCODE) && !defined(GPUCA_STANDALONE) && defined(GPUCA_O2_LIB) ~TPCFastTransform() diff --git a/GPU/Utils/FlatObject.h b/GPU/Utils/FlatObject.h index c4123034c57d3..b0179b73d0f29 100644 --- a/GPU/Utils/FlatObject.h +++ b/GPU/Utils/FlatObject.h @@ -28,7 +28,7 @@ #include "GPUCommonRtypes.h" #include "GPUCommonLogger.h" -//#define GPUCA_GPUCODE // uncomment to test "GPU" mode +// #define GPUCA_GPUCODE // uncomment to test "GPU" mode namespace GPUCA_NAMESPACE { @@ -132,7 +132,10 @@ T* resizeArray(T*& ptr, int oldSize, int newSize, T* newPtr = nullptr) newPtr = new T[newSize]; } int mcp = std::min(newSize, oldSize); - std::memmove(newPtr, ptr, mcp * sizeof(T)); + if (mcp) { + assert(ptr); + std::memmove(newPtr, ptr, mcp * sizeof(T)); + } if (newSize > oldSize) { std::memset(newPtr + mcp, 0, (newSize - oldSize) * sizeof(T)); } @@ -564,7 +567,7 @@ inline void FlatObject::setFutureBufferAddress(char* futureFlatBufferPtr) mFlatBufferContainer = nullptr; } -#endif //GPUCA_GPUCODE_DEVICE +#endif // GPUCA_GPUCODE_DEVICE } // namespace gpu } // namespace GPUCA_NAMESPACE