From 3e4bc95032433b9fe1740ddf5b471bdf4cad53c4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Nicol=C3=B2=20Jacazio?= Date: Mon, 25 Nov 2024 16:05:57 +0100 Subject: [PATCH] [Common] [TOF] introduce generic resolution parametrization (#8616) --- Common/Core/PID/PIDTOF.h | 76 +++++++++++--- Common/TableProducer/PID/pidTOFMerge.cxx | 122 +++++++++++++++++------ 2 files changed, 155 insertions(+), 43 deletions(-) diff --git a/Common/Core/PID/PIDTOF.h b/Common/Core/PID/PIDTOF.h index 5ecd32cbbd0..3c507223feb 100644 --- a/Common/Core/PID/PIDTOF.h +++ b/Common/Core/PID/PIDTOF.h @@ -28,17 +28,19 @@ #include "TMath.h" #include "TGraph.h" #include "TFile.h" +#include "TF2.h" // O2 includes #include "DataFormatsTOF/ParameterContainers.h" #include "Framework/Logger.h" #include "ReconstructionDataFormats/PID.h" #include "Framework/DataTypes.h" +#include "CommonConstants/PhysicsConstants.h" namespace o2::pid::tof { -// Utility values +// Utility values (to remove!) static constexpr float kCSPEED = TMath::C() * 1.0e2f * 1.0e-12f; /// Speed of light in TOF units (cm/ps) static constexpr float kCSPEDDInv = 1.f / kCSPEED; /// Inverse of the Speed of light in TOF units (ps/cm) static constexpr float defaultReturnValue = -999.f; /// Default return value in case TOF measurement is not available @@ -55,7 +57,7 @@ class Beta /// \param length Length in cm of the track /// \param tofSignal TOF signal in ps for the track /// \param collisionTime collision time in ps for the event of the track - static float GetBeta(const float length, const float tofSignal, const float collisionTime) { return length / (tofSignal - collisionTime) * kCSPEDDInv; } + static float GetBeta(const float length, const float tofSignal, const float collisionTime) { return length / (tofSignal - collisionTime) * o2::constants::physics::invLightSpeedCm2PS; } /// Gets the beta for the track of interest /// \param track Track of interest @@ -139,6 +141,11 @@ class TOFResoParamsV2 : public o2::tof::Parameters<13> ~TOFResoParamsV2() = default; + template + float getResolution(const float, const float) const + { + return -1.f; + } // Momentum shift for charge calibration void setMomentumChargeShiftParameters(std::unordered_map const& pars) { @@ -254,14 +261,12 @@ class TOFResoParamsV2 : public o2::tof::Parameters<13> class TOFResoParamsV3 : public o2::tof::Parameters<13> { public: - TOFResoParamsV3() : Parameters(std::array{"TrkRes.Pi.P0", "TrkRes.Pi.P1", "TrkRes.Pi.P2", "TrkRes.Pi.P3", "time_resolution", - "TrkRes.Ka.P0", "TrkRes.Ka.P1", "TrkRes.Ka.P2", "TrkRes.Ka.P3", - "TrkRes.Pr.P0", "TrkRes.Pr.P1", "TrkRes.Pr.P2", "TrkRes.Pr.P3"}, + TOFResoParamsV3() : Parameters(std::array{"time_resolution", "time_resolution", "time_resolution", "time_resolution", "time_resolution", + "time_resolution", "time_resolution", "time_resolution", "time_resolution", + "time_resolution", "time_resolution", "time_resolution", "time_resolution"}, "TOFResoParamsV3") { - setParameters(std::array{0.008, 0.008, 0.002, 40.0, 60.0, - 0.008, 0.008, 0.002, 40.0, - 0.008, 0.008, 0.002, 40.0}); + setParameters(std::array{60.0}); } // Default constructor with default parameters ~TOFResoParamsV3() = default; @@ -313,7 +318,7 @@ class TOFResoParamsV3 : public o2::tof::Parameters<13> } } - // Time shift for post calibration + // Time shift for post calibration to realign as a function of eta void setTimeShiftParameters(std::unordered_map const& pars, bool positive) { std::string baseOpt = positive ? "TimeShift.Pos." : "TimeShift.Neg."; @@ -367,6 +372,38 @@ class TOFResoParamsV3 : public o2::tof::Parameters<13> return gNegEtaTimeCorr->Eval(eta); } + void setResolutionParametrization(std::unordered_map const& pars) + { + static constexpr std::array particleNames = {"El", "Mu", "Pi", "Ka", "Pr", "De", "Tr", "He", "Al"}; + for (int i = 0; i < 9; ++i) { + const std::string baseOpt = Form("tofResTrack.%s_", particleNames[i]); + // Check if a key begins with a string + for (const auto& [key, value] : pars) { + if (key.find(baseOpt) == 0) { + // Remove from the key the baseOpt + const std::string fun = key.substr(baseOpt.size()); + mResolution[i] = new TF2(baseOpt.c_str(), fun.c_str(), 0., 20, -1, 1.); + LOG(info) << "Set the resolution function for " << particleNames[i] << " with formula " << mResolution[i]->GetFormula()->GetExpFormula(); + break; + } + } + } + // Print a summary + for (int i = 0; i < 9; ++i) { + if (!mResolution[i]) { + LOG(info) << "Resolution function for " << particleNames[i] << " not provided, using default " << mDefaultResoParams[i]; + mResolution[i] = new TF2(Form("tofResTrack.%s_Default", particleNames[i]), mDefaultResoParams[i], 0., 20, -1, 1.); + } + LOG(info) << "Resolution function for " << particleNames[i] << " is " << mResolution[i]->GetName() << " with formula " << mResolution[i]->GetFormula()->GetExpFormula(); + } + } + + template + float getResolution(const float p, const float eta) const + { + return mResolution[pid]->Eval(p, eta); + } + private: // Charge calibration int mEtaN = 0; // Number of eta bins, 0 means no correction @@ -374,6 +411,16 @@ class TOFResoParamsV3 : public o2::tof::Parameters<13> float mEtaStop = 0.f; float mInvEtaWidth = 9999.f; std::vector mContent; + std::array mResolution{nullptr}; + static constexpr std::array mDefaultResoParams{"14.3*TMath::Power((TMath::Max(x-0.319,0.1))*(1-0.4235*y*y),-0.8467)", + "14.3*TMath::Power((TMath::Max(x-0.319,0.1))*(1-0.4235*y*y),-0.8467)", + "14.3*TMath::Power((TMath::Max(x-0.319,0.1))*(1-0.4235*y*y),-0.8467)", + "42.66*TMath::Power((TMath::Max(x-0.417,0.1))*(1-0.4235*y*y),-0.7145)", + "99.46*TMath::Power((TMath::Max(x-0.447,0.1))*(1-0.4235*y*y),-0.8094)", + "216*TMath::Power((TMath::Max(x-0.647,0.1))*(1-0.4235*y*y),-0.76)", + "315*TMath::Power((TMath::Max(x-0.811,0.1))*(1-0.4235*y*y),-0.783)", + "157*TMath::Power((TMath::Max(x-0.556,0.1))*(1-0.4235*y*y),-0.783)", + "216*TMath::Power((TMath::Max(x-0.647,0.1))*(1-0.4235*y*y),-0.76)"}; // Time shift for post calibration TGraph* gPosEtaTimeCorr = nullptr; /// Time shift correction for positive tracks @@ -391,7 +438,7 @@ class ExpTimes static constexpr float mMassZSqared = mMassZ * mMassZ; /// (M/z)^2 /// Computes the expected time of a track, given it TOF expected momentum - static float ComputeExpectedTime(const float tofExpMom, const float length) { return length * sqrt((mMassZSqared) + (tofExpMom * tofExpMom)) / (kCSPEED * tofExpMom); } + static float ComputeExpectedTime(const float tofExpMom, const float length) { return length * sqrt((mMassZSqared) + (tofExpMom * tofExpMom)) / (o2::constants::physics::LightSpeedCm2PS * tofExpMom); } /// Gets the expected signal of the track of interest under the PID assumption /// \param track Track of interest @@ -401,7 +448,7 @@ class ExpTimes return defaultReturnValue; } if (track.trackType() == o2::aod::track::Run2Track) { - return ComputeExpectedTime(track.tofExpMom() * kCSPEDDInv, track.length()); + return ComputeExpectedTime(track.tofExpMom() * o2::constants::physics::invLightSpeedCm2PS, track.length()); } return ComputeExpectedTime(track.tofExpMom(), track.length()); } @@ -416,7 +463,7 @@ class ExpTimes return defaultReturnValue; } if (track.trackType() == o2::aod::track::Run2Track) { - return ComputeExpectedTime(track.tofExpMom() * kCSPEDDInv / (1.f + track.sign() * parameters.getMomentumChargeShift(track.eta())), track.length()); + return ComputeExpectedTime(track.tofExpMom() * o2::constants::physics::invLightSpeedCm2PS / (1.f + track.sign() * parameters.getMomentumChargeShift(track.eta())), track.length()); } LOG(debug) << "TOF exp. mom. " << track.tofExpMom() << " shifted = " << track.tofExpMom() / (1.f + track.sign() * parameters.getMomentumChargeShift(track.eta())); return ComputeExpectedTime(track.tofExpMom() / (1.f + track.sign() * parameters.getMomentumChargeShift(track.eta())), track.length()) + parameters.getTimeShift(track.eta(), track.sign()); @@ -432,9 +479,14 @@ class ExpTimes static float GetExpectedSigma(const ParamType& parameters, const TrackType& track, const float tofSignal, const float collisionTimeRes) { const float& mom = track.p(); + const float& eta = track.eta(); if (mom <= 0) { return -999.f; } + const float reso = parameters.template getResolution(mom, eta); + if (reso > 0) { + return std::sqrt(reso * reso + parameters[4] * parameters[4] + collisionTimeRes * collisionTimeRes); + } if constexpr (id <= o2::track::PID::Pion) { LOG(debug) << "Using parameters for the pion hypothesis and ID " << id; const float dpp = parameters[0] + parameters[1] * mom + parameters[2] * mMassZ / mom; // mean relative pt resolution; diff --git a/Common/TableProducer/PID/pidTOFMerge.cxx b/Common/TableProducer/PID/pidTOFMerge.cxx index cfd3c566c2b..047dc041cd5 100644 --- a/Common/TableProducer/PID/pidTOFMerge.cxx +++ b/Common/TableProducer/PID/pidTOFMerge.cxx @@ -19,6 +19,8 @@ #include #include #include +#include +#include // O2 includes #include "Framework/runDataProcessing.h" @@ -71,6 +73,7 @@ struct TOFCalibConfig { mParamFileName = opt.cfgParamFileName.value; mParametrizationPath = opt.cfgParametrizationPath.value; mReconstructionPass = opt.cfgReconstructionPass.value; + mReconstructionPassDefault = opt.cfgReconstructionPassDefault.value; mLoadResponseFromCCDB = opt.cfgLoadResponseFromCCDB.value; mFatalOnPassNotAvailable = opt.cfgFatalOnPassNotAvailable.value; mEnableTimeDependentResponse = opt.cfgEnableTimeDependentResponse.value; @@ -81,7 +84,7 @@ struct TOFCalibConfig { template void getCfg(o2::framework::InitContext& initContext, const std::string name, VType& v, const std::string task) { - if (!getTaskOptionValue(initContext, task, name, v, true)) { + if (!getTaskOptionValue(initContext, task, name, v, false)) { LOG(fatal) << "Could not get " << name << " from " << task << " task"; } } @@ -93,10 +96,11 @@ struct TOFCalibConfig { getCfg(initContext, "ccdb-path-grplhcif", mPathGrpLhcIf, task); getCfg(initContext, "ccdb-timestamp", mTimestamp, task); getCfg(initContext, "timeShiftCCDBPathPos", mTimeShiftCCDBPathPos, task); - getCfg(initContext, "timeShiftCCDBPathNeg", mTimeShiftCCDBPathPos, task); + getCfg(initContext, "timeShiftCCDBPathNeg", mTimeShiftCCDBPathNeg, task); getCfg(initContext, "paramFileName", mParamFileName, task); getCfg(initContext, "parametrizationPath", mParametrizationPath, task); getCfg(initContext, "reconstructionPass", mReconstructionPass, task); + getCfg(initContext, "reconstructionPassDefault", mReconstructionPassDefault, task); getCfg(initContext, "loadResponseFromCCDB", mLoadResponseFromCCDB, task); getCfg(initContext, "fatalOnPassNotAvailable", mFatalOnPassNotAvailable, task); getCfg(initContext, "enableTimeDependentResponse", mEnableTimeDependentResponse, task); @@ -129,49 +133,63 @@ struct TOFCalibConfig { } LOG(info) << "Using parameter collection, starting from pass '" << mReconstructionPass << "'"; - const std::string fname = mParamFileName; - if (!fname.empty()) { // Loading the parametrization from file - LOG(info) << "Loading exp. sigma parametrization from file " << fname << ", using param: " << mParametrizationPath; - if (1) { - o2::tof::ParameterCollection paramCollection; - paramCollection.loadParamFromFile(fname, mParametrizationPath); - LOG(info) << "+++ Loaded parameter collection from file +++"; - if (!paramCollection.retrieveParameters(mRespParamsV3, mReconstructionPass)) { - if (mFatalOnPassNotAvailable) { - LOGF(fatal, "Pass '%s' not available in the retrieved CCDB object", mReconstructionPass.data()); + if (!mParamFileName.empty()) { // Loading the parametrization from file + LOG(info) << "Loading exp. sigma parametrization from file " << mParamFileName << ", using param: " << mParametrizationPath << " and pass " << mReconstructionPass; + o2::tof::ParameterCollection paramCollection; + paramCollection.loadParamFromFile(mParamFileName, mParametrizationPath); + LOG(info) << "+++ Loaded parameter collection from file +++"; + if (!paramCollection.retrieveParameters(mRespParamsV3, mReconstructionPass)) { + if (mFatalOnPassNotAvailable) { + LOG(fatal) << "Pass '" << mReconstructionPass << "' not available in the retrieved object from file"; + } else { + LOG(warning) << "Pass '" << mReconstructionPass << "' not available in the retrieved object from file, fetching '" << mReconstructionPassDefault << "'"; + if (!paramCollection.retrieveParameters(mRespParamsV3, mReconstructionPassDefault)) { + paramCollection.print(); + LOG(fatal) << "Cannot get default pass for calibration " << mReconstructionPassDefault; } else { - LOGF(warning, "Pass '%s' not available in the retrieved CCDB object", mReconstructionPass.data()); + mRespParamsV3.setResolutionParametrization(paramCollection.getPars(mReconstructionPassDefault)); + mRespParamsV3.setMomentumChargeShiftParameters(paramCollection.getPars(mReconstructionPassDefault)); + mRespParamsV3.printMomentumChargeShiftParameters(); } - } else { - mRespParamsV3.setMomentumChargeShiftParameters(paramCollection.getPars(mReconstructionPass)); - mRespParamsV3.printMomentumChargeShiftParameters(); } - } else { - mRespParamsV3.loadParamFromFile(fname.data(), mParametrizationPath); + } else { // Pass is available, load non standard parameters + mRespParamsV3.setResolutionParametrization(paramCollection.getPars(mReconstructionPass)); + mRespParamsV3.setMomentumChargeShiftParameters(paramCollection.getPars(mReconstructionPass)); + mRespParamsV3.printMomentumChargeShiftParameters(); } - } else if (mLoadResponseFromCCDB) { // Loading it from CCDB - LOG(info) << "Loading exp. sigma parametrization from CCDB, using path: " << mParametrizationPath << " for timestamp " << mTimestamp; + } else if (mLoadResponseFromCCDB && !mEnableTimeDependentResponse) { // Loading it from CCDB + LOG(info) << "Loading initial exp. sigma parametrization from CCDB, using path: " << mParametrizationPath << " for timestamp " << mTimestamp; o2::tof::ParameterCollection* paramCollection = ccdb->template getForTimeStamp(mParametrizationPath, mTimestamp); paramCollection->print(); if (!paramCollection->retrieveParameters(mRespParamsV3, mReconstructionPass)) { // Attempt at loading the parameters with the pass defined if (mFatalOnPassNotAvailable) { - LOGF(fatal, "Pass '%s' not available in the retrieved CCDB object", mReconstructionPass.data()); + LOG(fatal) << "Pass '" << mReconstructionPass << "' not available in the retrieved CCDB object"; } else { - LOGF(warning, "Pass '%s' not available in the retrieved CCDB object", mReconstructionPass.data()); + LOG(warning) << "Pass '" << mReconstructionPass << "' not available in the retrieved CCDB object, fetching '" << mReconstructionPassDefault << "'"; + if (!paramCollection->retrieveParameters(mRespParamsV3, mReconstructionPassDefault)) { + paramCollection->print(); + LOG(fatal) << "Cannot get default pass for calibration " << mReconstructionPassDefault; + } else { + mRespParamsV3.setResolutionParametrization(paramCollection->getPars(mReconstructionPassDefault)); + mRespParamsV3.setMomentumChargeShiftParameters(paramCollection->getPars(mReconstructionPassDefault)); + mRespParamsV3.printMomentumChargeShiftParameters(); + } } } else { // Pass is available, load non standard parameters + mRespParamsV3.setResolutionParametrization(paramCollection->getPars(mReconstructionPass)); mRespParamsV3.setMomentumChargeShiftParameters(paramCollection->getPars(mReconstructionPass)); mRespParamsV3.printMomentumChargeShiftParameters(); } + } else { + std::unordered_map m; + mRespParamsV3.setResolutionParametrization(m); } - // Calibration object is defined - mRespParamsV3.print(); // Loading additional calibration objects if (mTimeShiftCCDBPathPos != "") { if (mTimeShiftCCDBPathPos.find(".root") != std::string::npos) { mRespParamsV3.setTimeShiftParameters(mTimeShiftCCDBPathPos, "ccdb_object", true); - } else { + } else if (!mEnableTimeDependentResponse) { if (mReconstructionPass == "") { mRespParamsV3.setTimeShiftParameters(ccdb->template getForTimeStamp(mTimeShiftCCDBPathPos, mTimestamp), true); } else { @@ -184,7 +202,7 @@ struct TOFCalibConfig { if (mTimeShiftCCDBPathNeg != "") { if (mTimeShiftCCDBPathNeg.find(".root") != std::string::npos) { mRespParamsV3.setTimeShiftParameters(mTimeShiftCCDBPathNeg, "ccdb_object", false); - } else { + } else if (!mEnableTimeDependentResponse) { if (mReconstructionPass == "") { mRespParamsV3.setTimeShiftParameters(ccdb->template getForTimeStamp(mTimeShiftCCDBPathNeg, mTimestamp), false); } else { @@ -194,6 +212,10 @@ struct TOFCalibConfig { } } } + + // Calibration object is defined + LOG(info) << "Parametrization at init time:"; + mRespParamsV3.print(); } template @@ -221,14 +243,50 @@ struct TOFCalibConfig { if (!mEnableTimeDependentResponse) { return; } - LOG(debug) << "Updating parametrization from path '" << mParametrizationPath << "' and timestamp " << mTimestamp; - if (!ccdb->template getForTimeStamp(mParametrizationPath, mTimestamp)->retrieveParameters(mRespParamsV3, mReconstructionPass)) { - if (mFatalOnPassNotAvailable) { - LOGF(fatal, "Pass '%s' not available in the retrieved CCDB object", mReconstructionPass.data()); + LOG(info) << "Updating parametrization from path '" << mParametrizationPath << "' and timestamp " << mTimestamp << " and reconstruction pass '" << mReconstructionPass << "'"; + if (mParamFileName.empty()) { // Not loading if parametrization from file + if (!ccdb->template getForTimeStamp(mParametrizationPath, mTimestamp)->retrieveParameters(mRespParamsV3, mReconstructionPass)) { + if (mFatalOnPassNotAvailable) { + LOGF(fatal, "Pass '%s' not available in the retrieved CCDB object", mReconstructionPass.data()); + } else { + LOGF(warning, "Pass '%s' not available in the retrieved CCDB object, fetching '%s'", mReconstructionPass.data(), mReconstructionPassDefault.data()); + if (!ccdb->template getForTimeStamp(mParametrizationPath, mTimestamp)->retrieveParameters(mRespParamsV3, mReconstructionPassDefault)) { + ccdb->template getForTimeStamp(mParametrizationPath, mTimestamp)->print(); + LOG(fatal) << "Cannot get default pass for calibration " << mReconstructionPassDefault; + } + } + } + } + + // Loading additional calibration objects + if (mTimeShiftCCDBPathPos != "") { + if (mTimeShiftCCDBPathPos.find(".root") != std::string::npos) { + mRespParamsV3.setTimeShiftParameters(mTimeShiftCCDBPathPos, "ccdb_object", true); + } else { + if (mReconstructionPass == "") { + mRespParamsV3.setTimeShiftParameters(ccdb->template getForTimeStamp(mTimeShiftCCDBPathPos, mTimestamp), true); + } else { + std::map metadata; + metadata["RecoPassName"] = mReconstructionPass; + mRespParamsV3.setTimeShiftParameters(ccdb->template getSpecific(mTimeShiftCCDBPathPos, mTimestamp, metadata), true); + } + } + } + if (mTimeShiftCCDBPathNeg != "") { + if (mTimeShiftCCDBPathNeg.find(".root") != std::string::npos) { + mRespParamsV3.setTimeShiftParameters(mTimeShiftCCDBPathNeg, "ccdb_object", false); } else { - LOGF(warning, "Pass '%s' not available in the retrieved CCDB object", mReconstructionPass.data()); + if (mReconstructionPass == "") { + mRespParamsV3.setTimeShiftParameters(ccdb->template getForTimeStamp(mTimeShiftCCDBPathNeg, mTimestamp), false); + } else { + std::map metadata; + metadata["RecoPassName"] = mReconstructionPass; + mRespParamsV3.setTimeShiftParameters(ccdb->template getSpecific(mTimeShiftCCDBPathNeg, mTimestamp, metadata), false); + } } + LOG(info) << " test getTimeShift neg: " << mRespParamsV3.getTimeShift(0, false); } + return; } @@ -248,6 +306,7 @@ struct TOFCalibConfig { std::string mParamFileName; std::string mParametrizationPath; std::string mReconstructionPass; + std::string mReconstructionPassDefault; bool mLoadResponseFromCCDB; bool mFatalOnPassNotAvailable; bool mEnableTimeDependentResponse; @@ -289,6 +348,7 @@ struct tofSignal { Configurable cfgParamFileName{"paramFileName", "", "Path to the parametrization object. If empty the parametrization is not taken from file"}; Configurable cfgParametrizationPath{"parametrizationPath", "TOF/Calib/Params", "Path of the TOF parametrization on the CCDB or in the file, if the paramFileName is not empty"}; Configurable cfgReconstructionPass{"reconstructionPass", "", {"Apass to use when fetching the calibration tables. Empty (default) does not check for any pass. Use `metadata` to fetch it from the AO2D metadata. Otherwise it will override the metadata."}}; + Configurable cfgReconstructionPassDefault{"reconstructionPassDefault", "unanchored", {"Default pass to get if the standard one is not found"}}; Configurable cfgLoadResponseFromCCDB{"loadResponseFromCCDB", false, "Flag to load the response from the CCDB"}; Configurable cfgFatalOnPassNotAvailable{"fatalOnPassNotAvailable", true, "Flag to throw a fatal if the pass is not available in the retrieved CCDB object"}; Configurable cfgEnableTimeDependentResponse{"enableTimeDependentResponse", false, "Flag to use the collision timestamp to fetch the PID Response"};