From b82f8ebbb601548a88ee76a82859b9c43342800c Mon Sep 17 00:00:00 2001 From: Matthias Kleiner Date: Thu, 23 Nov 2023 19:13:03 +0100 Subject: [PATCH] Implement scaling of distortions in MC - Scaling of the distortions can be enabled with `o2-sim-digitizer-workflow --tpc-distortion-type 2`, where the lumi is taken from `TPCCorrMap.lumiInst`. - space charge: adding meta data for scaling - In reconstruction the --corrmap-lumi-mode 2 should be used for MC and luminosity set via TPCCorrMap.lumiInst - Fixing inverse transform for scalMode 1 and 2! - other small changes --- Detectors/TPC/base/include/TPCBase/CDBTypes.h | 11 ++- .../calibration/src/CorrectionMapsLoader.cxx | 9 ++- Detectors/TPC/simulation/CMakeLists.txt | 2 +- .../include/TPCSimulation/Digitizer.h | 32 +++++--- Detectors/TPC/simulation/src/Digitizer.cxx | 34 +++++++-- .../include/TPCSpaceCharge/SpaceCharge.h | 21 +++++- .../TPCSpaceCharge/SpaceChargeHelpers.h | 22 ++++++ Detectors/TPC/spacecharge/src/SpaceCharge.cxx | 74 ++++++++++++++++--- .../spacecharge/src/TPCSpacechargeLinkDef.h | 1 + .../CorrectionMapsHelper.h | 2 +- GPU/TPCFastTransformation/TPCFastTransform.h | 30 ++++---- .../src/SimpleDigitizerWorkflow.cxx | 6 +- .../src/TPCDigitizerSpec.cxx | 47 +++++++++--- .../DigitizerWorkflow/src/TPCDigitizerSpec.h | 4 +- 14 files changed, 230 insertions(+), 65 deletions(-) diff --git a/Detectors/TPC/base/include/TPCBase/CDBTypes.h b/Detectors/TPC/base/include/TPCBase/CDBTypes.h index 34851092ee285..5bab9ffaada87 100644 --- a/Detectors/TPC/base/include/TPCBase/CDBTypes.h +++ b/Detectors/TPC/base/include/TPCBase/CDBTypes.h @@ -71,6 +71,8 @@ enum class CDBType { /// CalCorrMap, ///< Cluster correction map (high IR rate distortions) CalCorrMapRef, ///< Cluster correction reference map (static distortions) + CalCorrMapMC, ///< Cluster correction map (high IR rate distortions) for MC + CalCorrMapRefMC, ///< Cluster correction reference map (static distortions) for MC /// CalCorrDerivMap, ///< Cluster correction map (derivative map) /// @@ -79,7 +81,8 @@ enum class CDBType { /// CorrMapParam, ///< parameters for CorrectionMapsLoader configuration /// - DistortionMap, ///< distortions for MC used in the digitizer + DistortionMapMC, ///< full distortions (static + IR dependant) for MC used in the digitizer + DistortionMapDerivMC ///< derivative distortions for MC used in the digitizer for scaling }; /// Storage name in CCDB for each calibration and parameter type @@ -131,6 +134,9 @@ const std::unordered_map CDBTypeMap{ // correction maps {CDBType::CalCorrMap, "TPC/Calib/CorrectionMapV2"}, {CDBType::CalCorrMapRef, "TPC/Calib/CorrectionMapRefV2"}, + // correction maps for MC + {CDBType::CalCorrMapMC, "TPC/Calib/CorrectionMapMCV2"}, + {CDBType::CalCorrMapRefMC, "TPC/Calib/CorrectionMapRefMCV2"}, // derivative map correction {CDBType::CalCorrDerivMap, "TPC/Calib/CorrectionMapDerivativeV2"}, // time series @@ -139,7 +145,8 @@ const std::unordered_map CDBTypeMap{ // correction maps loader params {CDBType::CorrMapParam, "TPC/Calib/CorrMapParam"}, // distortion maps - {CDBType::DistortionMap, "TPC/Calib/DistortionMap"}, + {CDBType::DistortionMapMC, "TPC/Calib/DistortionMapMC"}, + {CDBType::DistortionMapDerivMC, "TPC/Calib/DistortionMapDerivativeMC"}, }; } // namespace o2::tpc diff --git a/Detectors/TPC/calibration/src/CorrectionMapsLoader.cxx b/Detectors/TPC/calibration/src/CorrectionMapsLoader.cxx index 9a2af80934a7e..7c517f33b05e7 100644 --- a/Detectors/TPC/calibration/src/CorrectionMapsLoader.cxx +++ b/Detectors/TPC/calibration/src/CorrectionMapsLoader.cxx @@ -67,11 +67,16 @@ void CorrectionMapsLoader::extractCCDBInputs(ProcessingContext& pc) //________________________________________________________ void CorrectionMapsLoader::requestCCDBInputs(std::vector& inputs, std::vector& options, const CorrectionMapsLoaderGloOpts& gloOpts) { - addInput(inputs, {"tpcCorrMap", "TPC", "CorrMap", 0, Lifetime::Condition, ccdbParamSpec(CDBTypeMap.at(CDBType::CalCorrMap), {}, 1)}); // time-dependent if (gloOpts.lumiMode == 0) { + addInput(inputs, {"tpcCorrMap", "TPC", "CorrMap", 0, Lifetime::Condition, ccdbParamSpec(CDBTypeMap.at(CDBType::CalCorrMap), {}, 1)}); // time-dependent addInput(inputs, {"tpcCorrMapRef", "TPC", "CorrMapRef", 0, Lifetime::Condition, ccdbParamSpec(CDBTypeMap.at(CDBType::CalCorrMapRef), {}, 0)}); // load once } else if (gloOpts.lumiMode == 1) { + addInput(inputs, {"tpcCorrMap", "TPC", "CorrMap", 0, Lifetime::Condition, ccdbParamSpec(CDBTypeMap.at(CDBType::CalCorrMap), {}, 1)}); // time-dependent addInput(inputs, {"tpcCorrMapRef", "TPC", "CorrMapRef", 0, Lifetime::Condition, ccdbParamSpec(CDBTypeMap.at(CDBType::CalCorrDerivMap), {}, 1)}); // time-dependent + } else if (gloOpts.lumiMode == 2) { + // for MC corrections + addInput(inputs, {"tpcCorrMap", "TPC", "CorrMap", 0, Lifetime::Condition, ccdbParamSpec(CDBTypeMap.at(CDBType::CalCorrMapMC), {}, 1)}); // time-dependent + addInput(inputs, {"tpcCorrMapRef", "TPC", "CorrMapRef", 0, Lifetime::Condition, ccdbParamSpec(CDBTypeMap.at(CDBType::CalCorrMapRefMC), {}, 1)}); // time-dependent } else { LOG(fatal) << "Correction mode unknown! Choose either 0 (default) or 1 (derivative map) for flag corrmap-lumi-mode."; } @@ -99,7 +104,7 @@ void CorrectionMapsLoader::addGlobalOptions(std::vector& option { // these are options which should be added at the workflow level, since they modify the inputs of the devices addOption(options, ConfigParamSpec{"lumi-type", o2::framework::VariantType::Int, 0, {"1 = require CTP lumi for TPC correction scaling, 2 = require TPC scalers for TPC correction scaling"}}); - addOption(options, ConfigParamSpec{"corrmap-lumi-mode", o2::framework::VariantType::Int, 0, {"scaling mode: (default) 0 = static + scale * full; 1 = full + scale * derivative"}}); + addOption(options, ConfigParamSpec{"corrmap-lumi-mode", o2::framework::VariantType::Int, 0, {"scaling mode: (default) 0 = static + scale * full; 1 = full + scale * derivative; 2 = full + scale * derivative (for MC)"}}); } //________________________________________________________ diff --git a/Detectors/TPC/simulation/CMakeLists.txt b/Detectors/TPC/simulation/CMakeLists.txt index 1f433469ab769..510d42d2d6d85 100644 --- a/Detectors/TPC/simulation/CMakeLists.txt +++ b/Detectors/TPC/simulation/CMakeLists.txt @@ -22,7 +22,7 @@ o2_add_library(TPCSimulation src/SAMPAProcessing.cxx src/IDCSim.cxx PUBLIC_LINK_LIBRARIES O2::DetectorsBase O2::SimulationDataFormat - O2::TPCBase O2::TPCSpaceCharge + O2::TPCBase O2::TPCSpaceCharge O2::TPCCalibration ROOT::Physics) o2_target_root_dictionary(TPCSimulation diff --git a/Detectors/TPC/simulation/include/TPCSimulation/Digitizer.h b/Detectors/TPC/simulation/include/TPCSimulation/Digitizer.h index 4080968c6018c..e2c7e9ec7d175 100644 --- a/Detectors/TPC/simulation/include/TPCSimulation/Digitizer.h +++ b/Detectors/TPC/simulation/include/TPCSimulation/Digitizer.h @@ -121,6 +121,9 @@ class Digitizer /// \param spaceCharge unique pointer to spaceCharge object void setUseSCDistortions(SC* spaceCharge); + /// \param spaceCharge unique pointer to spaceCharge object + void setSCDistortionsDerivative(SC* spaceCharge); + /// Enable the use of space-charge distortions by providing global distortions and global corrections stored in a ROOT file /// The storage of the values should be done by the methods provided in the SpaceCharge class /// \param file containing distortions @@ -129,17 +132,26 @@ class Digitizer void setVDrift(float v) { mVDrift = v; } void setTDriftOffset(float t) { mTDriftOffset = t; } + void setDistortionScaleType(int distortionScaleType) { mDistortionScaleType = distortionScaleType; } + int getDistortionScaleType() const { return mDistortionScaleType; } + void setLumiScaleFactor(); + void setMeanLumiDistortions(float meanLumi); + void setMeanLumiDistortionsDerivative(float meanLumi); + private: - DigitContainer mDigitContainer; ///< Container for the Digits - std::unique_ptr mSpaceCharge; ///< Handler of space-charge distortions - Sector mSector = -1; ///< ID of the currently processed sector - double mEventTime = 0.f; ///< Time of the currently processed event - double mOutputDigitTimeOffset = 0; ///< Time of the first IR sampled in the digitizer - float mVDrift = 0; ///< VDrift for current timestamp - float mTDriftOffset = 0; ///< drift time additive offset in \mus - bool mIsContinuous; ///< Switch for continuous readout - bool mUseSCDistortions = false; ///< Flag to switch on the use of space-charge distortions - ClassDefNV(Digitizer, 1); + DigitContainer mDigitContainer; ///< Container for the Digits + std::unique_ptr mSpaceCharge; ///< Handler of full distortions (static + IR dependant) + std::unique_ptr mSpaceChargeDer; ///< Handler of reference static distortions + Sector mSector = -1; ///< ID of the currently processed sector + double mEventTime = 0.f; ///< Time of the currently processed event + double mOutputDigitTimeOffset = 0; ///< Time of the first IR sampled in the digitizer + float mVDrift = 0; ///< VDrift for current timestamp + float mTDriftOffset = 0; ///< drift time additive offset in \mus + bool mIsContinuous; ///< Switch for continuous readout + bool mUseSCDistortions = false; ///< Flag to switch on the use of space-charge distortions + int mDistortionScaleType = 0; ///< type=0: no scaling of distortions, type=1 distortions without any scaling, type=2 distortions scaling with lumi + float mLumiScaleFactor = 0; ///< value used to scale the derivative map + ClassDefNV(Digitizer, 2); }; } // namespace tpc } // namespace o2 diff --git a/Detectors/TPC/simulation/src/Digitizer.cxx b/Detectors/TPC/simulation/src/Digitizer.cxx index 7bac9e6ebb3e0..38968b3b38bb9 100644 --- a/Detectors/TPC/simulation/src/Digitizer.cxx +++ b/Detectors/TPC/simulation/src/Digitizer.cxx @@ -27,6 +27,7 @@ #include "TPCBase/CDBInterface.h" #include "TPCSpaceCharge/SpaceCharge.h" #include "TPCBase/Mapper.h" +#include "TPCCalibration/CorrMapParam.h" #include @@ -40,10 +41,6 @@ Digitizer::Digitizer() = default; void Digitizer::init() { - // Calculate distortion lookup tables if initial space-charge density is provided - if (mUseSCDistortions) { - mSpaceCharge->init(); - } auto& gemAmplification = GEMAmplification::instance(); gemAmplification.updateParameters(); auto& electronTransport = ElectronTransport::instance(); @@ -83,8 +80,10 @@ void Digitizer::process(const std::vector& hits, GlobalPosition3D posEle(eh.GetX(), eh.GetY(), eh.GetZ()); // Distort the electron position in case space-charge distortions are used - if (mUseSCDistortions) { + if (mDistortionScaleType == 1) { mSpaceCharge->distortElectron(posEle); + } else if (mDistortionScaleType == 2) { + mSpaceCharge->distortElectron(posEle, mSpaceChargeDer.get(), mLumiScaleFactor); } /// Remove electrons that end up more than three sigma of the hit's average diffusion away from the current sector @@ -190,6 +189,15 @@ void Digitizer::setUseSCDistortions(SC* spaceCharge) { mUseSCDistortions = true; mSpaceCharge.reset(spaceCharge); + mSpaceCharge->initAfterReadingFromFile(); + mSpaceCharge->printMetaData(); +} + +void Digitizer::setSCDistortionsDerivative(SC* spaceCharge) +{ + mSpaceChargeDer.reset(spaceCharge); + mSpaceChargeDer->initAfterReadingFromFile(); + mSpaceChargeDer->printMetaData(); } void Digitizer::setUseSCDistortions(std::string_view finp) @@ -213,3 +221,19 @@ void Digitizer::setStartTime(double time) sampaProcessing.updateParameters(mVDrift); mDigitContainer.setStartTime(sampaProcessing.getTimeBinFromTime(time - mOutputDigitTimeOffset)); } + +void Digitizer::setLumiScaleFactor() +{ + mLumiScaleFactor = (CorrMapParam::Instance().lumiInst - mSpaceCharge->getMeanLumi()) / mSpaceChargeDer->getMeanLumi(); + LOGP(info, "Setting Lumi scale factor: lumiInst: {} lumi mean: {} lumi mean derivative: {} lumi scale factor: {}", CorrMapParam::Instance().lumiInst, mSpaceCharge->getMeanLumi(), mSpaceChargeDer->getMeanLumi(), mLumiScaleFactor); +} + +void Digitizer::setMeanLumiDistortions(float meanLumi) +{ + mSpaceCharge->setMeanLumi(meanLumi); +} + +void Digitizer::setMeanLumiDistortionsDerivative(float meanLumi) +{ + mSpaceChargeDer->setMeanLumi(meanLumi); +} diff --git a/Detectors/TPC/spacecharge/include/TPCSpaceCharge/SpaceCharge.h b/Detectors/TPC/spacecharge/include/TPCSpaceCharge/SpaceCharge.h index 014be1f2cc50e..dc6bb1be3aba9 100644 --- a/Detectors/TPC/spacecharge/include/TPCSpaceCharge/SpaceCharge.h +++ b/Detectors/TPC/spacecharge/include/TPCSpaceCharge/SpaceCharge.h @@ -572,7 +572,9 @@ class SpaceCharge /// Distort electron position using distortion lookup tables /// \param point 3D coordinates of the electron - void distortElectron(GlobalPosition3D& point) const; + /// \param scSCale other sc object which is used for scaling of the distortions + /// \param scale scaling value + void distortElectron(GlobalPosition3D& point, const SpaceCharge* scSCale = nullptr, float scale = 0) const; /// set the distortions directly from a look up table /// \param distdZ distortions in z direction @@ -1170,12 +1172,24 @@ class SpaceCharge /// substract global corrections from other sc object (global corrections -= other.global corrections) /// can be used to calculate the derivative: (this - other)/normalization /// for normalization see scaleCorrections() - void substractGlobalCorrections(const SpaceCharge& otherSC, const Side side); + void subtractGlobalCorrections(const SpaceCharge& otherSC, const Side side); + + /// substract global distortions from other sc object (global distortions -= other.global distortions) + /// can be used to calculate the derivative: (this - other)/normalization + void subtractGlobalDistortions(const SpaceCharge& otherSC, const Side side); /// scale corrections by factor /// \param scaleFac global corrections are multiplied by this factor void scaleCorrections(const float scaleFac, const Side side); + /// setting meta data for this object + void setMetaData(const SCMetaData& meta) { mMeta = meta; } + const auto& getMetaData() const { return mMeta; } + void printMetaData() const { mMeta.print(); } + float getMeanLumi() const { return mMeta.meanLumi; } + void setMeanLumi(float lumi) { mMeta.meanLumi = lumi; } + void initAfterReadingFromFile(); + private: ParamSpaceCharge mParamGrid{}; ///< parameters of the grid on which the calculations are performed inline static int sNThreads{getOMPMaxThreads()}; /// mAnaDistCorr; ///< analytical distortions and corrections bool mUseAnaDistCorr{false}; ///< flag if analytical distortions will be used in the distortElectron() and getCorrections() function BField mBField{}; /// collisionTypes{"PP", "Pb-Pb"}; + if (collisionType < collisionTypes.size()) { + LOGP(info, "meanLumi: {}, IR: {}kHz, run: {}, collisionType {}", meanLumi, ir, run, collisionTypes[collisionType]); + } else { + LOGP(info, "Specified collision type {} not allowed", collisionType); + } + } + + float meanLumi = 0; ///< mean lumi the sc object corresponds to + float ir = 0; ///< IR + int run = 0; ///< run number this object anchored to to + int collisionType = 0; ///< 0=PP, 1-Pb-Pb + + private: + ClassDefNV(SCMetaData, 1); +}; + /// /// this class contains an analytical description of the space charge, potential and the electric fields. /// The analytical functions can be used to test the poisson solver and the caluclation of distortions/corrections. diff --git a/Detectors/TPC/spacecharge/src/SpaceCharge.cxx b/Detectors/TPC/spacecharge/src/SpaceCharge.cxx index 3ed00345159bc..d3fa62b94a1c9 100644 --- a/Detectors/TPC/spacecharge/src/SpaceCharge.cxx +++ b/Detectors/TPC/spacecharge/src/SpaceCharge.cxx @@ -19,7 +19,6 @@ #include "fmt/core.h" #include "Framework/Logger.h" #include "TPCSpaceCharge/PoissonSolver.h" -#include "Framework/Logger.h" #include "TGeoGlobalMagField.h" #include "TPCBase/ParameterGas.h" #include "TPCBase/ParameterElectronics.h" @@ -709,7 +708,7 @@ void SpaceCharge::calcGlobalDistWithGlobalCorrIterative(const DistCorrInt const DataT phi = getPhiVertex(iPhi, side); for (unsigned int iR = 0; iR < mParamGrid.NRVertices; ++iR) { const DataT radius = getRVertex(iR, side); - for (unsigned int iZ = 0; iZ < mParamGrid.NZVertices; ++iZ) { + for (unsigned int iZ = 1; iZ < mParamGrid.NZVertices; ++iZ) { const DataT z = getZVertex(iZ, side); unsigned int nearestiZ = iZ; @@ -789,6 +788,11 @@ void SpaceCharge::calcGlobalDistWithGlobalCorrIterative(const DistCorrInt mGlobalDistdZ[side](iZ, iR, iPhi) = -corrdZ; } } + for (unsigned int iR = 0; iR < mParamGrid.NRVertices; ++iR) { + mGlobalDistdR[side](0, iR, iPhi) = 3 * (mGlobalDistdR[side](1, iR, iPhi) - mGlobalDistdR[side](2, iR, iPhi)) + mGlobalDistdR[side](3, iR, iPhi); + mGlobalDistdRPhi[side](0, iR, iPhi) = 3 * (mGlobalDistdRPhi[side](1, iR, iPhi) - mGlobalDistdRPhi[side](2, iR, iPhi)) + mGlobalDistdRPhi[side](3, iR, iPhi); + mGlobalDistdZ[side](0, iR, iPhi) = 3 * (mGlobalDistdZ[side](1, iR, iPhi) - mGlobalDistdZ[side](2, iR, iPhi)) + mGlobalDistdZ[side](3, iR, iPhi); + } } } @@ -1516,7 +1520,7 @@ void SpaceCharge::correctElectron(GlobalPosition3D& point) } template -void SpaceCharge::distortElectron(GlobalPosition3D& point) const +void SpaceCharge::distortElectron(GlobalPosition3D& point, const SpaceCharge* scSCale, float scale) const { DataT distX{}; DataT distY{}; @@ -1525,6 +1529,18 @@ void SpaceCharge::distortElectron(GlobalPosition3D& point) const // get the distortions for input coordinate getDistortions(point.X(), point.Y(), point.Z(), side, distX, distY, distZ); + DataT distXTmp{}; + DataT distYTmp{}; + DataT distZTmp{}; + + // scale distortions if requested + if (scSCale && scale != 0) { + scSCale->getDistortions(point.X(), point.Y(), point.Z(), side, distXTmp, distYTmp, distZTmp); + distX += distXTmp * scale; + distY += distYTmp * scale; + distZ += distZTmp * scale; + } + GPUCA_DEBUG_STREAMER_CHECK(if (o2::utils::DebugStreamer::checkStream(o2::utils::StreamFlags::streamDistortionsSC)) { GlobalPosition3D pos(point); float phi = std::atan2(pos.Y(), pos.X()); @@ -1543,6 +1559,10 @@ void SpaceCharge::distortElectron(GlobalPosition3D& point) const << "distX=" << distX << "distY=" << distY << "distZ=" << distZ + << "distXDer=" << distXTmp + << "distYDer=" << distYTmp + << "distZDer=" << distZTmp + << "scale=" << scale << "\n"; }) @@ -1734,8 +1754,10 @@ void SpaceCharge::getDistortionsCyl(const std::vector& z, const st template void SpaceCharge::getDistortions(const DataT x, const DataT y, const DataT z, const Side side, DataT& distX, DataT& distY, DataT& distZ) const { + DataT zClamped = regulateZ(z, side); + if (mUseAnaDistCorr) { - getDistortionsAnalytical(x, y, z, side, distX, distY, distZ); + getDistortionsAnalytical(x, y, zClamped, side, distX, distY, distZ); } else { // convert cartesian to polar const DataT radius = getRadiusFromCartesian(x, y); @@ -1743,11 +1765,12 @@ void SpaceCharge::getDistortions(const DataT x, const DataT y, const Data DataT distR{}; DataT distRPhi{}; - getDistortionsCyl(z, radius, phi, side, distZ, distR, distRPhi); + DataT rClamped = regulateR(radius, side); + getDistortionsCyl(zClamped, rClamped, phi, side, distZ, distR, distRPhi); // Calculate distorted position - const DataT radiusDist = radius + distR; - const DataT phiDist = phi + distRPhi / radius; + const DataT radiusDist = rClamped + distR; + const DataT phiDist = phi + distRPhi / rClamped; distX = getXFromPolar(radiusDist, phiDist) - x; // difference between distorted and original x coordinate distY = getYFromPolar(radiusDist, phiDist) - y; // difference between distorted and original y coordinate @@ -2268,7 +2291,13 @@ void SpaceCharge::dumpToTree(const char* outFileName, const Side side, co DataT corrZ{}; DataT corrR{}; DataT corrRPhi{}; - getCorrectionsCyl(zPos, rPos, phiPos, side, corrZ, corrR, corrRPhi); + // getCorrectionsCyl(zPos, rPos, phiPos, side, corrZ, corrR, corrRPhi); + + const DataT zDistorted = zPos + distZ; + const DataT radiusDistorted = rPos + distR; + const DataT phiDistorted = regulatePhi(phiPos + distRPhi / rPos, side); + getCorrectionsCyl(zDistorted, radiusDistorted, phiDistorted, side, corrZ, corrR, corrRPhi); + corrRPhi *= rPos / radiusDistorted; DataT lcorrZ{}; DataT lcorrR{}; @@ -3046,12 +3075,13 @@ void SpaceCharge::dumpMetaData(std::string_view file, std::string_view op dfStore = dfStore.DefineSlotEntry("grid_A", [&helperA = helperA](unsigned int, ULong64_t entry) { return helperA; }); dfStore = dfStore.DefineSlotEntry("grid_C", [&helperC = helperC](unsigned int, ULong64_t entry) { return helperC; }); dfStore = dfStore.DefineSlotEntry("BField", [field = mBField.getBField()](unsigned int, ULong64_t entry) { return field; }); + dfStore = dfStore.DefineSlotEntry("metaInf", [meta = mMeta](unsigned int, ULong64_t entry) { return meta; }); // write to TTree ROOT::RDF::RSnapshotOptions opt; opt.fMode = option; opt.fOverwriteIfExists = true; // overwrite if already exists - dfStore.Snapshot("meta", file, {"paramsC", "grid_A", "grid_C", "BField"}, opt); + dfStore.Snapshot("meta", file, {"paramsC", "grid_A", "grid_C", "BField", "metaInf"}, opt); } template @@ -3079,6 +3109,15 @@ void SpaceCharge::readMetaData(std::string_view file) ROOT::RDataFrame dFrame("meta", file); dFrame.Foreach(readMeta, {"paramsC", "grid_A", "grid_C", "BField"}); + + const auto& cols = dFrame.GetColumnNames(); + if (std::find(cols.begin(), cols.end(), "metaInf") != cols.end()) { + auto readMetaInf = [&mMeta = mMeta](const SCMetaData& meta) { + mMeta = meta; + }; + dFrame.Foreach(readMetaInf, {"metaInf"}); + } + LOGP(info, "Setting meta data: mC0={} mC1={} mC2={}", mC0, mC1, mC2); mReadMetaData = true; } @@ -3608,13 +3647,21 @@ void SpaceCharge::fillROCMisalignment(const std::vector& indicesT } template -void SpaceCharge::substractGlobalCorrections(const SpaceCharge& otherSC, const Side side) +void SpaceCharge::subtractGlobalCorrections(const SpaceCharge& otherSC, const Side side) { mGlobalCorrdR[side] -= otherSC.mGlobalCorrdR[side]; mGlobalCorrdZ[side] -= otherSC.mGlobalCorrdZ[side]; mGlobalCorrdRPhi[side] -= otherSC.mGlobalCorrdRPhi[side]; } +template +void SpaceCharge::subtractGlobalDistortions(const SpaceCharge& otherSC, const Side side) +{ + mGlobalDistdR[side] -= otherSC.mGlobalDistdR[side]; + mGlobalDistdZ[side] -= otherSC.mGlobalDistdZ[side]; + mGlobalDistdRPhi[side] -= otherSC.mGlobalDistdRPhi[side]; +} + template void SpaceCharge::scaleCorrections(const float val, const Side side) { @@ -3694,6 +3741,13 @@ void SpaceCharge::scaleChargeDensityStack(const float scalingFactor, cons } } +template +void SpaceCharge::initAfterReadingFromFile() +{ + mGrid3D[Side::A] = RegularGrid(GridProp::ZMIN, GridProp::RMIN, GridProp::PHIMIN, getSign(Side::A) * GridProp::getGridSpacingZ(mParamGrid.NZVertices), GridProp::getGridSpacingR(mParamGrid.NRVertices), GridProp::getGridSpacingPhi(mParamGrid.NPhiVertices), mParamGrid); + mGrid3D[Side::C] = RegularGrid(GridProp::ZMIN, GridProp::RMIN, GridProp::PHIMIN, getSign(Side::C) * GridProp::getGridSpacingZ(mParamGrid.NZVertices), GridProp::getGridSpacingR(mParamGrid.NRVertices), GridProp::getGridSpacingPhi(mParamGrid.NPhiVertices), mParamGrid); +} + using DataTD = double; template class o2::tpc::SpaceCharge; diff --git a/Detectors/TPC/spacecharge/src/TPCSpacechargeLinkDef.h b/Detectors/TPC/spacecharge/src/TPCSpacechargeLinkDef.h index 7574b6896c315..b1bb33dc1aeb7 100644 --- a/Detectors/TPC/spacecharge/src/TPCSpacechargeLinkDef.h +++ b/Detectors/TPC/spacecharge/src/TPCSpacechargeLinkDef.h @@ -47,4 +47,5 @@ #pragma link C++ class o2::tpc::AnalyticalDistCorr < double> + ; #pragma link C++ class o2::tpc::AnalyticalDistCorr < float> + ; #pragma link C++ struct o2::tpc::ParamSpaceCharge + ; +#pragma link C++ struct o2::tpc::SCMetaData + ; #endif diff --git a/GPU/TPCFastTransformation/CorrectionMapsHelper.h b/GPU/TPCFastTransformation/CorrectionMapsHelper.h index 75d194bef63cb..970f9b4c88a26 100644 --- a/GPU/TPCFastTransformation/CorrectionMapsHelper.h +++ b/GPU/TPCFastTransformation/CorrectionMapsHelper.h @@ -99,7 +99,7 @@ class CorrectionMapsHelper { if (mMeanLumi < 0.f || mInstLumi < 0.f) { mLumiScale = -1.f; - } else if (mLumiScaleMode == 1) { + } else if ((mLumiScaleMode == 1) || (mLumiScaleMode == 2)) { mLumiScale = mMeanLumiRef ? (mInstLumi - mMeanLumi) / mMeanLumiRef : 0.f; LOGP(debug, "mInstLumi: {} mMeanLumi: {} mMeanLumiRef: {}", mInstLumi, mMeanLumi, mMeanLumiRef); } else { diff --git a/GPU/TPCFastTransformation/TPCFastTransform.h b/GPU/TPCFastTransformation/TPCFastTransform.h index 956ae13646ff6..1936ae4f19a99 100644 --- a/GPU/TPCFastTransformation/TPCFastTransform.h +++ b/GPU/TPCFastTransformation/TPCFastTransform.h @@ -442,7 +442,7 @@ GPUdi() void TPCFastTransform::TransformInternal(int slice, int row, float& u, f { if (mApplyCorrection) { float dx = 0.f, du = 0.f, dv = 0.f; - if ((scale >= 0.f) || (scaleMode == 1)) { + if ((scale >= 0.f) || (scaleMode == 1) || (scaleMode == 2)) { #ifndef GPUCA_GPUCODE if (mCorrectionSlow) { float ly, lz; @@ -470,7 +470,7 @@ GPUdi() void TPCFastTransform::TransformInternal(int slice, int row, float& u, f dx = (dx - dxRef) * scale + dxRef; du = (du - duRef) * scale + duRef; dv = (dv - dvRef) * scale + dvRef; - } else if ((scale != 0.f) && (scaleMode == 1)) { + } else if ((scale != 0.f) && ((scaleMode == 1) || (scaleMode == 2))) { float dxRef, duRef, dvRef; ref->mCorrection.getCorrection(slice, row, u, v, dxRef, duRef, dvRef); dx = dxRef * scale + dx; @@ -494,11 +494,11 @@ GPUdi() void TPCFastTransform::TransformInternal(int slice, int row, float& u, f getGeometry().convUVtoLocal(slice, uCorr, vCorr, lyT, lzT); float invYZtoX; - InverseTransformYZtoX(slice, row, ly, lz, invYZtoX); + InverseTransformYZtoX(slice, row, lyT, lzT, invYZtoX, ref, scale, scaleMode); float YZtoNominalY; float YZtoNominalZ; - InverseTransformYZtoNominalYZ(slice, row, ly, lz, YZtoNominalY, YZtoNominalZ); + InverseTransformYZtoNominalYZ(slice, row, lyT, lzT, YZtoNominalY, YZtoNominalZ, ref, scale, scaleMode); float dxRef, duRef, dvRef; ref->mCorrection.getCorrection(slice, row, u, v, dxRef, duRef, dvRef); @@ -745,17 +745,17 @@ GPUdi() void TPCFastTransform::InverseTransformYZtoX(int slice, int row, float y /// Transformation y,z -> x float u = 0, v = 0; getGeometry().convLocalToUV(slice, y, z, u, v); - if (scale >= 0.f) { + if ((scale >= 0.f) || (scaleMode == 1) || (scaleMode == 2)) { mCorrection.getCorrectionInvCorrectedX(slice, row, u, v, x); - if (ref && scale > 0.f) { // scaling was requested - if (scaleMode == 0) { + if (ref) { // scaling was requested + if (scaleMode == 0 && scale > 0.f) { float xr; ref->mCorrection.getCorrectionInvCorrectedX(slice, row, u, v, xr); x = (x - xr) * scale + xr; - } else if (scaleMode == 1) { + } else if ((scale != 0) && ((scaleMode == 1) || (scaleMode == 2))) { float xr; ref->mCorrection.getCorrectionInvCorrectedX(slice, row, u, v, xr); - x = xr * scale + x; + x = (xr - getGeometry().getRowInfo(row).x) * scale + x; // xr=mGeo.getRowInfo(row).x + dx; } } } else { @@ -780,19 +780,19 @@ GPUdi() void TPCFastTransform::InverseTransformYZtoNominalYZ(int slice, int row, /// Transformation y,z -> x float u = 0, v = 0, un = 0, vn = 0; getGeometry().convLocalToUV(slice, y, z, u, v); - if (scale >= 0.f) { + if ((scale >= 0.f) || (scaleMode == 1) || (scaleMode == 2)) { mCorrection.getCorrectionInvUV(slice, row, u, v, un, vn); - if (ref && scale > 0.f) { // scaling was requested - if (scaleMode == 0) { + if (ref) { // scaling was requested + if (scaleMode == 0 && scale > 0.f) { float unr = 0, vnr = 0; ref->mCorrection.getCorrectionInvUV(slice, row, u, v, unr, vnr); un = (un - unr) * scale + unr; vn = (vn - vnr) * scale + vnr; - } else if (scaleMode == 1) { + } else if ((scale != 0) && ((scaleMode == 1) || (scaleMode == 2))) { float unr = 0, vnr = 0; ref->mCorrection.getCorrectionInvUV(slice, row, u, v, unr, vnr); - un = unr * scale + un; - vn = vnr * scale + vn; + un = (unr - u) * scale + un; // unr = u - duv[0]; + vn = (vnr - v) * scale + vn; } } } else { diff --git a/Steer/DigitizerWorkflow/src/SimpleDigitizerWorkflow.cxx b/Steer/DigitizerWorkflow/src/SimpleDigitizerWorkflow.cxx index 5b796ec793272..f1e3bdff9a0c6 100644 --- a/Steer/DigitizerWorkflow/src/SimpleDigitizerWorkflow.cxx +++ b/Steer/DigitizerWorkflow/src/SimpleDigitizerWorkflow.cxx @@ -169,7 +169,7 @@ void customize(std::vector& workflowOptions) // Option to write TPC digits internaly, without forwarding to a special writer instance. // This is useful in GRID productions with small available memory. workflowOptions.push_back(ConfigParamSpec{"tpc-chunked-writer", o2::framework::VariantType::Bool, false, {"Write independent TPC digit chunks as soon as they can be flushed."}}); - workflowOptions.push_back(ConfigParamSpec{"tpc-simulate-distortions", o2::framework::VariantType::Bool, false, {"Simulate distortions in the TPC."}}); + workflowOptions.push_back(ConfigParamSpec{"tpc-distortion-type", o2::framework::VariantType::Int, 0, {"Simulate distortions in the TPC (0=no distortions, 1=distortions without scaling, 2=distortions with CTP scaling)"}}); std::string simhelp("Comma separated list of simulation prefixes (for background, signal productions)"); workflowOptions.push_back( @@ -569,8 +569,8 @@ WorkflowSpec defineDataProcessing(ConfigContext const& configcontext) detList.emplace_back(o2::detectors::DetID::TPC); auto internalwrite = configcontext.options().get("tpc-chunked-writer"); - auto simDis = configcontext.options().get("tpc-simulate-distortions"); - WorkflowSpec tpcPipelines = o2::tpc::getTPCDigitizerSpec(lanes, tpcsectors, mctruth, internalwrite, simDis); + auto distortionType = configcontext.options().get("tpc-distortion-type"); + WorkflowSpec tpcPipelines = o2::tpc::getTPCDigitizerSpec(lanes, tpcsectors, mctruth, internalwrite, distortionType); specs.insert(specs.end(), tpcPipelines.begin(), tpcPipelines.end()); if (configcontext.options().get("tpc-reco-type").empty() == false) { diff --git a/Steer/DigitizerWorkflow/src/TPCDigitizerSpec.cxx b/Steer/DigitizerWorkflow/src/TPCDigitizerSpec.cxx index 92c7b54acd4ba..59f8589a56449 100644 --- a/Steer/DigitizerWorkflow/src/TPCDigitizerSpec.cxx +++ b/Steer/DigitizerWorkflow/src/TPCDigitizerSpec.cxx @@ -108,7 +108,7 @@ using namespace o2::base; class TPCDPLDigitizerTask : public BaseDPLDigitizer { public: - TPCDPLDigitizerTask(bool internalwriter, bool simulateDistortions) : mInternalWriter(internalwriter), BaseDPLDigitizer(InitServices::FIELD | InitServices::GEOM), mSimulateDistortions(simulateDistortions) + TPCDPLDigitizerTask(bool internalwriter, int distortionType) : mInternalWriter(internalwriter), BaseDPLDigitizer(InitServices::FIELD | InitServices::GEOM), mDistortionType(distortionType) { } @@ -119,12 +119,15 @@ class TPCDPLDigitizerTask : public BaseDPLDigitizer mLaneId = ic.services().get().rank; mWithMCTruth = o2::conf::DigiParams::Instance().mctruth; - auto useDistortions = ic.options().get("distortionType"); auto triggeredMode = ic.options().get("TPCtriggered"); mUseCalibrationsFromCCDB = ic.options().get("TPCuseCCDB"); + mMeanLumiDistortions = ic.options().get("meanLumiDistortions"); + mMeanLumiDistortionsDerivative = ic.options().get("meanLumiDistortionsDerivative"); + LOG(info) << "TPC calibrations from CCDB: " << mUseCalibrationsFromCCDB; mDigitizer.setContinuousReadout(!triggeredMode); + mDigitizer.setDistortionScaleType(mDistortionType); // we send the GRP data once if the corresponding output channel is available // and set the flag to false after @@ -190,6 +193,16 @@ class TPCDPLDigitizerTask : public BaseDPLDigitizer if (matcher == ConcreteDataMatcher(o2::header::gDataOriginTPC, "TPCDIST", 0)) { LOGP(info, "Updating distortion map"); mDigitizer.setUseSCDistortions(static_cast(obj)); + if (mMeanLumiDistortions >= 0) { + mDigitizer.setMeanLumiDistortions(mMeanLumiDistortions); + } + } + if (matcher == ConcreteDataMatcher(o2::header::gDataOriginTPC, "TPCDISTDERIV", 0)) { + LOGP(info, "Updating reference distortion map"); + mDigitizer.setSCDistortionsDerivative(static_cast(obj)); + if (mMeanLumiDistortionsDerivative >= 0) { + mDigitizer.setMeanLumiDistortionsDerivative(mMeanLumiDistortionsDerivative); + } } } @@ -202,9 +215,14 @@ class TPCDPLDigitizerTask : public BaseDPLDigitizer cdb.setUseDefaults(!mUseCalibrationsFromCCDB); // whatever are global settings for CCDB usage, we have to extract the TPC vdrift from CCDB for anchored simulations mTPCVDriftHelper.extractCCDBInputs(pc); - if (mSimulateDistortions) { + if (mDistortionType) { pc.inputs().get("tpcdistortions"); + if (mDistortionType == 2) { + pc.inputs().get("tpcdistortionsderiv"); + mDigitizer.setLumiScaleFactor(); + } } + if (mTPCVDriftHelper.isUpdated()) { const auto& vd = mTPCVDriftHelper.getVDriftObject(); LOGP(info, "Updating TPC fast transform map with new VDrift factor of {} wrt reference {} and DriftTimeOffset correction {} wrt {} from source {}", @@ -456,10 +474,12 @@ class TPCDPLDigitizerTask : public BaseDPLDigitizer bool mWithMCTruth = true; bool mInternalWriter = false; bool mUseCalibrationsFromCCDB = false; - bool mSimulateDistortions = false; + int mDistortionType = 0; + float mMeanLumiDistortions = -1; + float mMeanLumiDistortionsDerivative = -1; }; -o2::framework::DataProcessorSpec getTPCDigitizerSpec(int channel, bool writeGRP, bool mctruth, bool internalwriter, bool simulateDistortions) +o2::framework::DataProcessorSpec getTPCDigitizerSpec(int channel, bool writeGRP, bool mctruth, bool internalwriter, int distortionType) { // create the full data processor spec using // a name identifier @@ -489,21 +509,22 @@ o2::framework::DataProcessorSpec getTPCDigitizerSpec(int channel, bool writeGRP, id.str().c_str(), inputs, outputs, - AlgorithmSpec{adaptFromTask(internalwriter, simulateDistortions)}, + AlgorithmSpec{adaptFromTask(internalwriter, distortionType)}, Options{ - {"distortionType", VariantType::Int, 0, {"Distortion type to be used. 0 = no distortions (default), 1 = realistic distortions (not implemented yet), 2 = constant distortions"}}, {"TPCtriggered", VariantType::Bool, false, {"Impose triggered RO mode (default: continuous)"}}, {"TPCuseCCDB", VariantType::Bool, false, {"true: load calibrations from CCDB; false: use random calibratoins"}}, + {"meanLumiDistortions", VariantType::Float, -1.f, {"override lumi of distortion object if >=0"}}, + {"meanLumiDistortionsDerivative", VariantType::Float, -1.f, {"override lumi of derivative distortion object if >=0"}}, }}; } -o2::framework::WorkflowSpec getTPCDigitizerSpec(int nLanes, std::vector const& sectors, bool mctruth, bool internalwriter, bool simulateDistortions) +o2::framework::WorkflowSpec getTPCDigitizerSpec(int nLanes, std::vector const& sectors, bool mctruth, bool internalwriter, int distortionType) { // channel parameter is deprecated in the TPCDigitizer processor, all descendants // are initialized not to publish GRP mode, but the channel will be added to the first // processor after the pipelines have been created. The processor will decide upon // the index in the ParallelContext whether to publish - WorkflowSpec pipelineTemplate{getTPCDigitizerSpec(0, false, mctruth, internalwriter, simulateDistortions)}; + WorkflowSpec pipelineTemplate{getTPCDigitizerSpec(0, false, mctruth, internalwriter, distortionType)}; // override the predefined name, index will be added by parallelPipeline method pipelineTemplate[0].name = "TPCDigitizer"; WorkflowSpec pipelines = parallelPipeline( @@ -511,8 +532,12 @@ o2::framework::WorkflowSpec getTPCDigitizerSpec(int nLanes, std::vector con // add the channel for the GRP information to the first processor for (auto& spec : pipelines) { o2::tpc::VDriftHelper::requestCCDBInputs(spec.inputs); // add the same CCDB request to each pipeline - if (simulateDistortions) { - spec.inputs.emplace_back("tpcdistortions", o2::header::gDataOriginTPC, "TPCDIST", 0, Lifetime::Condition, ccdbParamSpec(o2::tpc::CDBTypeMap.at(o2::tpc::CDBType::DistortionMap), {}, 1)); // time-dependent + if (distortionType) { + spec.inputs.emplace_back("tpcdistortions", o2::header::gDataOriginTPC, "TPCDIST", 0, Lifetime::Condition, ccdbParamSpec(o2::tpc::CDBTypeMap.at(o2::tpc::CDBType::DistortionMapMC), {}, 1)); // time-dependent + // load derivative map in case scaling was requested + if (distortionType == 2) { + spec.inputs.emplace_back("tpcdistortionsderiv", o2::header::gDataOriginTPC, "TPCDISTDERIV", 0, Lifetime::Condition, ccdbParamSpec(o2::tpc::CDBTypeMap.at(o2::tpc::CDBType::DistortionMapDerivMC), {}, 1)); // time-dependent + } } } pipelines[0].outputs.emplace_back("TPC", "ROMode", 0, Lifetime::Timeframe); diff --git a/Steer/DigitizerWorkflow/src/TPCDigitizerSpec.h b/Steer/DigitizerWorkflow/src/TPCDigitizerSpec.h index 5b29c9d917749..65837f7974ba1 100644 --- a/Steer/DigitizerWorkflow/src/TPCDigitizerSpec.h +++ b/Steer/DigitizerWorkflow/src/TPCDigitizerSpec.h @@ -20,9 +20,9 @@ namespace o2 namespace tpc { -o2::framework::DataProcessorSpec getTPCDigitizerSpec(int channel, bool writeGRP, bool mctruth, bool internalwriter, bool simulateDistortions); +o2::framework::DataProcessorSpec getTPCDigitizerSpec(int channel, bool writeGRP, bool mctruth, bool internalwriter, int distortionType); -o2::framework::WorkflowSpec getTPCDigitizerSpec(int nLanes, std::vector const& sectors, bool mctruth, bool internalwriter, bool simulateDistortions); +o2::framework::WorkflowSpec getTPCDigitizerSpec(int nLanes, std::vector const& sectors, bool mctruth, bool internalwriter, int distortionType); } // end namespace tpc } // end namespace o2