From 52114c98178fdb91746b0fe0f60891cb2bafde05 Mon Sep 17 00:00:00 2001 From: Himanshu Sharma Date: Fri, 26 Apr 2024 09:30:14 +0200 Subject: [PATCH] Common: TrackTuner with Pt smearing in MC (#5823) --- Common/TableProducer/trackPropagation.cxx | 71 ++- Common/Tools/TrackTuner.h | 692 +++++++++++++--------- 2 files changed, 457 insertions(+), 306 deletions(-) diff --git a/Common/TableProducer/trackPropagation.cxx b/Common/TableProducer/trackPropagation.cxx index c20790016c1..ac8a4e60e0e 100644 --- a/Common/TableProducer/trackPropagation.cxx +++ b/Common/TableProducer/trackPropagation.cxx @@ -38,6 +38,8 @@ struct TrackPropagation { Produces tracksDCA; Produces tracksDCACov; + Produces tunertable; + Service ccdb; bool fillTracksDCA = false; @@ -58,12 +60,19 @@ struct TrackPropagation { Configurable mVtxPath{"mVtxPath", "GLO/Calib/MeanVertex", "Path of the mean vertex file"}; Configurable minPropagationRadius{"minPropagationDistance", o2::constants::geom::XTPCInnerRef + 0.1, "Only tracks which are at a smaller radius will be propagated, defaults to TPC inner wall"}; // for TrackTuner only (MC smearing) - Configurable useTrackTuner{"useTrackTuner", false, "Apply Improver/DCA corrections to MC"}; - Configurable trackTunerParams{"trackTunerParams", "debugInfo=0|updateTrackCovMat=1|updateCurvature=0|updatePulls=0|isInputFileFromCCDB=1|pathInputFile=Users/m/mfaggin/test/inputsTrackTuner/PbPb2022|nameInputFile=trackTuner_DataLHC22sPass5_McLHC22l1b2_run529397.root|usePvRefitCorrections=0|oneOverPtCurrent=0|oneOverPtUpgr=0", "TrackTuner parameter initialization (format: =|=)"}; - OutputObj trackTunedTracks{TH1D("trackTunedTracks", "", 4, 0.5, 4.5), OutputObjHandlingPolicy::AnalysisObject}; + Configurable useTrackTuner{"useTrackTuner", false, "Apply track tuner corrections to MC"}; + Configurable fillTrackTunerTable{"fillTrackTunerTable", false, "flag to fill track tuner table"}; + Configurable trackTunerParams{"trackTunerParams", "debugInfo=0|updateTrackDCAs=1|updateTrackCovMat=1|updateCurvature=0|updateCurvatureIU=0|updatePulls=0|isInputFileFromCCDB=1|pathInputFile=Users/m/mfaggin/test/inputsTrackTuner/PbPb2022|nameInputFile=trackTuner_DataLHC22sPass5_McLHC22l1b2_run529397.root|pathFileQoverPt=Users/h/hsharma/qOverPtGraphs|nameFileQoverPt=D0sigma_Data_removal_itstps_MC_LHC22b1b.root|usePvRefitCorrections=0|qOverPtMC=-1.|qOverPtData=-1.", "TrackTuner parameter initialization (format: =|=)"}; + ConfigurableAxis axisPtQA{"axisPtQA", {VARIABLE_WIDTH, 0.0f, 0.1f, 0.2f, 0.3f, 0.4f, 0.5f, 0.6f, 0.7f, 0.8f, 0.9f, 1.0f, 1.1f, 1.2f, 1.3f, 1.4f, 1.5f, 1.6f, 1.7f, 1.8f, 1.9f, 2.0f, 2.2f, 2.4f, 2.6f, 2.8f, 3.0f, 3.2f, 3.4f, 3.6f, 3.8f, 4.0f, 4.4f, 4.8f, 5.2f, 5.6f, 6.0f, 6.5f, 7.0f, 7.5f, 8.0f, 9.0f, 10.0f, 11.0f, 12.0f, 13.0f, 14.0f, 15.0f, 17.0f, 19.0f, 21.0f, 23.0f, 25.0f, 30.0f, 35.0f, 40.0f, 50.0f}, "pt axis for QA histograms"}; + OutputObj trackTunedTracks{TH1D("trackTunedTracks", "", 1, 0.5, 1.5), OutputObjHandlingPolicy::AnalysisObject}; + + // OutputObj hDCAxyVsPtRec{TH2F("hDCAxyVsPtRec", ";DCAxy;PtRec", 600, -0.15, 0.15, axisPtQA)}; + // OutputObj hDCAxyVsPtMC{TH2F("hDCAxyVsPtMC", ";DCAxy;PtMC", 600, -0.15, 0.15, axisPtQA)}; using TracksIUWithMc = soa::Join; + HistogramRegistry registry{"registry"}; + void init(o2::framework::InitContext& initContext) { int nEnabledProcesses = 0; @@ -101,6 +110,13 @@ struct TrackPropagation { ccdb->setLocalObjectValidityChecking(); lut = o2::base::MatLayerCylSet::rectifyPtrFromFile(ccdb->get(lutPath)); + // Histograms for track tuner + AxisSpec axisBinsDCA = {600, -0.15f, 0.15f, "#it{dca}_{xy} (cm)"}; + registry.add("hDCAxyVsPtRec", "hDCAxyVsPtRec", kTH2F, {axisBinsDCA, axisPtQA}); + registry.add("hDCAxyVsPtMC", "hDCAxyVsPtMC", kTH2F, {axisBinsDCA, axisPtQA}); + registry.add("hDCAzVsPtRec", "hDCAzVsPtRec", kTH2F, {axisBinsDCA, axisPtQA}); + registry.add("hDCAzVsPtMC", "hDCAzVsPtMC", kTH2F, {axisBinsDCA, axisPtQA}); + /// TrackTuner initialization if (useTrackTuner) { std::string outputStringParams = trackTunerObj.configParams(trackTunerParams); @@ -179,45 +195,65 @@ struct TrackPropagation { } // auto trackParCov = getTrackParCov(track); aod::track::TrackTypeEnum trackType = (aod::track::TrackTypeEnum)track.trackType(); + // std::array trackPxPyPz; + // std::array trackPxPyPzTuned = {0.0, 0.0, 0.0}; + double q2OverPtNew = -9999.; // Only propagate tracks which have passed the innermost wall of the TPC (e.g. skipping loopers etc). Others fill unpropagated. if (track.trackType() == aod::track::TrackIU && track.x() < minPropagationRadius) { - if constexpr (isMc && fillCovMat) { /// track tuner ok only if cov. matrix is used + if constexpr (isMc && fillCovMat) { // checking MC and fillCovMat block begins + // bool hasMcParticle = track.has_mcParticle(); if (useTrackTuner) { trackTunedTracks->Fill(1); // all tracks - // call track propagator - // this function reads many many things - // - reads track params bool hasMcParticle = track.has_mcParticle(); if (hasMcParticle) { - // LOG(info) << " MC particle exists... "; - // LOG(info) << "Inside trackPropagation: before calling tuneTrackParams trackParCov.getY(): " << trackParCov.getY(); auto mcParticle = track.mcParticle(); trackTunerObj.tuneTrackParams(mcParticle, mTrackParCov, matCorr, &mDcaInfoCov, trackTunedTracks); - // LOG(info) << "Inside trackPropagation: after calling tuneTrackParams trackParCov.getY(): " << trackParCov.getY(); - // trackTunedTracks->Fill(1); + q2OverPtNew = mTrackParCov.getQ2Pt(); } } - } + } // MC and fillCovMat block ends + bool isPropagationOK = true; + if (track.has_collision()) { auto const& collision = track.collision(); if constexpr (fillCovMat) { mVtx.setPos({collision.posX(), collision.posY(), collision.posZ()}); mVtx.setCov(collision.covXX(), collision.covXY(), collision.covYY(), collision.covXZ(), collision.covYZ(), collision.covZZ()); - o2::base::Propagator::Instance()->propagateToDCABxByBz(mVtx, mTrackParCov, 2.f, matCorr, &mDcaInfoCov); + isPropagationOK = o2::base::Propagator::Instance()->propagateToDCABxByBz(mVtx, mTrackParCov, 2.f, matCorr, &mDcaInfoCov); } else { - o2::base::Propagator::Instance()->propagateToDCABxByBz({collision.posX(), collision.posY(), collision.posZ()}, mTrackPar, 2.f, matCorr, &mDcaInfo); + isPropagationOK = o2::base::Propagator::Instance()->propagateToDCABxByBz({collision.posX(), collision.posY(), collision.posZ()}, mTrackPar, 2.f, matCorr, &mDcaInfo); } } else { if constexpr (fillCovMat) { mVtx.setPos({mMeanVtx->getX(), mMeanVtx->getY(), mMeanVtx->getZ()}); mVtx.setCov(mMeanVtx->getSigmaX() * mMeanVtx->getSigmaX(), 0.0f, mMeanVtx->getSigmaY() * mMeanVtx->getSigmaY(), 0.0f, 0.0f, mMeanVtx->getSigmaZ() * mMeanVtx->getSigmaZ()); - o2::base::Propagator::Instance()->propagateToDCABxByBz(mVtx, mTrackParCov, 2.f, matCorr, &mDcaInfoCov); + isPropagationOK = o2::base::Propagator::Instance()->propagateToDCABxByBz(mVtx, mTrackParCov, 2.f, matCorr, &mDcaInfoCov); } else { - o2::base::Propagator::Instance()->propagateToDCABxByBz({mMeanVtx->getX(), mMeanVtx->getY(), mMeanVtx->getZ()}, mTrackPar, 2.f, matCorr, &mDcaInfo); + isPropagationOK = o2::base::Propagator::Instance()->propagateToDCABxByBz({mMeanVtx->getX(), mMeanVtx->getY(), mMeanVtx->getZ()}, mTrackPar, 2.f, matCorr, &mDcaInfo); } } - trackType = aod::track::Track; + if (isPropagationOK) { + trackType = aod::track::Track; + } + // filling some QA histograms for track tuner test purpose + if constexpr (isMc && fillCovMat) { // checking MC and fillCovMat block begins + if (track.has_mcParticle() && isPropagationOK) { + auto mcParticle1 = track.mcParticle(); + // && abs(mcParticle1.pdgCode())==211 + if (mcParticle1.isPhysicalPrimary()) { + registry.fill(HIST("hDCAxyVsPtRec"), mDcaInfoCov.getY(), mTrackParCov.getPt()); + registry.fill(HIST("hDCAxyVsPtMC"), mDcaInfoCov.getY(), mcParticle1.pt()); + registry.fill(HIST("hDCAzVsPtRec"), mDcaInfoCov.getZ(), mTrackParCov.getPt()); + registry.fill(HIST("hDCAzVsPtMC"), mDcaInfoCov.getZ(), mcParticle1.pt()); + } + } + } // MC and fillCovMat block ends + } + // Filling modified Q/Pt values at IU/production point by track tuner in track tuner table + if (useTrackTuner && fillTrackTunerTable) { + tunertable(q2OverPtNew); } + // LOG(info) << " trackPropagation (this value filled in tuner table)--> " << q2OverPtNew; if constexpr (fillCovMat) { tracksParPropagated(track.collisionId(), trackType, mTrackParCov.getX(), mTrackParCov.getAlpha(), mTrackParCov.getY(), mTrackParCov.getZ(), mTrackParCov.getSnp(), mTrackParCov.getTgl(), mTrackParCov.getQ2Pt()); tracksParExtensionPropagated(mTrackParCov.getPt(), mTrackParCov.getP(), mTrackParCov.getEta(), mTrackParCov.getPhi()); @@ -259,6 +295,7 @@ struct TrackPropagation { // ----------------------- void processCovarianceMc(TracksIUWithMc const& tracks, aod::McParticles const& mcParticles, aod::Collisions const& collisions, aod::BCsWithTimestamps const& bcs) { + // auto table_extension = soa::Extend(tracks); fillTrackTables(tracks, mcParticles, collisions, bcs); } PROCESS_SWITCH(TrackPropagation, processCovarianceMc, "Process with covariance on MC", false); diff --git a/Common/Tools/TrackTuner.h b/Common/Tools/TrackTuner.h index d3e5861962e..0d5d8452c68 100644 --- a/Common/Tools/TrackTuner.h +++ b/Common/Tools/TrackTuner.h @@ -23,6 +23,7 @@ #include #include #include +#include #include "CCDB/BasicCCDBManager.h" #include "CCDB/CcdbApi.h" @@ -41,55 +42,67 @@ #include "Framework/RunningWorkflowInfo.h" #include "ReconstructionDataFormats/DCA.h" #include "ReconstructionDataFormats/Track.h" - #include -struct TrackTuner { +namespace o2::aod +{ +namespace track_tuner +{ +DECLARE_SOA_COLUMN(TunedQOverPt, tunedQOverPt, float); +} // namespace track_tuner + +DECLARE_SOA_TABLE(TrackTunerTable, "AOD", "TRACKTUNERTABLE", //! + track_tuner::TunedQOverPt); +} // namespace o2::aod +struct TrackTuner { /////////////////////////////// /// parameters to be configured bool debugInfo = false; + bool updateTrackDCAs = false; // To update the track DCAs; bool updateTrackCovMat = false; bool updateCurvature = false; + bool updateCurvatureIU = false; // To update the track parameter Q/Pt in trackIU table, particularly used for V0 mass width dependence on Q/Pt bool updatePulls = false; bool isInputFileFromCCDB = false; // query input file from CCDB or local folder - std::string pathInputFile = ""; // Path to file containing DCAxy, DCAz graphs from data (upgr) and MC (current) + std::string pathInputFile = ""; // Path to file containing DCAxy, DCAz graphs from data and MC std::string nameInputFile = ""; // Common Name of different files containing graphs, found in the above paths + std::string pathFileQoverPt = ""; // Path to file containing D0 sigma graphs from data and MC + std::string nameFileQoverPt = ""; // file name containing Q/Pt correction graphs from data and MC bool usePvRefitCorrections = false; // establish whether to use corrections obtained with or w/o PV refit - float oneOverPtCurrent = 0.; // 1/pt old - float oneOverPtUpgr = 0.; // 1/pt new + float qOverPtMC = -1.; // 1/pt old + float qOverPtData = -1.; // 1/pt new /////////////////////////////// o2::ccdb::CcdbApi ccdbApi; std::map metadata; - std::unique_ptr grDcaXYResVsPtPionCurrent; - std::unique_ptr grDcaXYResVsPtPionUpgr; + std::unique_ptr grDcaXYResVsPtPionMC; + std::unique_ptr grDcaXYResVsPtPionData; - std::unique_ptr grDcaZResVsPtPionCurrent; - std::unique_ptr grDcaZResVsPtPionUpgr; + std::unique_ptr grDcaZResVsPtPionMC; + std::unique_ptr grDcaZResVsPtPionData; - std::unique_ptr grDcaXYMeanVsPtPionCurrent; - std::unique_ptr grDcaXYMeanVsPtPionUpgr; + std::unique_ptr grDcaXYMeanVsPtPionMC; + std::unique_ptr grDcaXYMeanVsPtPionData; - std::unique_ptr grDcaZMeanVsPtPionCurrent; - std::unique_ptr grDcaZMeanVsPtPionUpgr; + std::unique_ptr grDcaZMeanVsPtPionMC; + std::unique_ptr grDcaZMeanVsPtPionData; - std::unique_ptr grOneOverPtPionCurrent; - std::unique_ptr grOneOverPtPionUpgr; + std::unique_ptr grOneOverPtPionMC; // MC + std::unique_ptr grOneOverPtPionData; // Data - std::unique_ptr grDcaXYPullVsPtPionCurrent; - std::unique_ptr grDcaXYPullVsPtPionUpgr; + std::unique_ptr grDcaXYPullVsPtPionMC; + std::unique_ptr grDcaXYPullVsPtPionData; - std::unique_ptr grDcaZPullVsPtPionCurrent; - std::unique_ptr grDcaZPullVsPtPionUpgr; + std::unique_ptr grDcaZPullVsPtPionMC; + std::unique_ptr grDcaZPullVsPtPionData; /// @brief Function to configure the TrackTuner parameters /// @param inputString Input string with all parameter configuration. Format: =|= /// @return String with the values of all parameters after configurations are listed, to cross check that everything worked well std::string configParams(std::string inputString) { - std::string delimiter = "|"; std::string assignmentSymbol = "="; @@ -122,39 +135,51 @@ struct TrackTuner { /// +++ to be manually updated every time one adds a new parameter to the TrackTuner.h +++ enum Pars : uint8_t { DebugInfo = 0, UpdateTrackCovMat, + UpdateTrackDCAs, UpdateCurvature, + UpdateCurvatureIU, UpdatePulls, - PathInputFile, IsInputFileFromCCDB, + PathInputFile, NameInputFile, + PathFileQoverPt, + NameFileQoverPt, UsePvRefitCorrections, - OneOverPtCurrent, - OneOverPtUpgr, + QOverPtMC, + QOverPtData, NPars }; std::map mapParNames = { std::make_pair(DebugInfo, "debugInfo"), + std::make_pair(UpdateTrackDCAs, "updateTrackDCAs"), std::make_pair(UpdateTrackCovMat, "updateTrackCovMat"), std::make_pair(UpdateCurvature, "updateCurvature"), + std::make_pair(UpdateCurvatureIU, "updateCurvatureIU"), std::make_pair(UpdatePulls, "updatePulls"), std::make_pair(IsInputFileFromCCDB, "isInputFileFromCCDB"), std::make_pair(PathInputFile, "pathInputFile"), + std::make_pair(PathFileQoverPt, "pathFileQoverPt"), std::make_pair(NameInputFile, "nameInputFile"), + std::make_pair(NameFileQoverPt, "nameFileQoverPt"), std::make_pair(UsePvRefitCorrections, "usePvRefitCorrections"), - std::make_pair(OneOverPtCurrent, "oneOverPtCurrent"), - std::make_pair(OneOverPtUpgr, "oneOverPtUpgr")}; + std::make_pair(QOverPtMC, "qOverPtMC"), + std::make_pair(QOverPtData, "qOverPtData")}; /////////////////////////////////////////////////////////////////////////////////// LOG(info) << "[TrackTuner]"; LOG(info) << "[TrackTuner] >>> Parameters before the custom settings"; LOG(info) << "[TrackTuner] debugInfo = " << debugInfo; + LOG(info) << "[TrackTuner] updateTrackDCAs = " << UpdateTrackDCAs; LOG(info) << "[TrackTuner] updateTrackCovMat = " << updateTrackCovMat; LOG(info) << "[TrackTuner] updateCurvature = " << updateCurvature; + LOG(info) << "[TrackTuner] updateCurvatureIU = " << updateCurvatureIU; LOG(info) << "[TrackTuner] updatePulls = " << updatePulls; LOG(info) << "[TrackTuner] isInputFileFromCCDB = " << isInputFileFromCCDB; LOG(info) << "[TrackTuner] pathInputFile = " << pathInputFile; LOG(info) << "[TrackTuner] nameInputFile = " << nameInputFile; + LOG(info) << "[TrackTuner] pathFileQoverPt = " << pathFileQoverPt; + LOG(info) << "[TrackTuner] nameFileQoverPt = " << nameFileQoverPt; LOG(info) << "[TrackTuner] usePvRefitCorrections = " << usePvRefitCorrections; - LOG(info) << "[TrackTuner] oneOverPtCurrent = " << oneOverPtCurrent; - LOG(info) << "[TrackTuner] oneOverPtUpgr = " << oneOverPtUpgr; + LOG(info) << "[TrackTuner] qOverPtMC = " << qOverPtMC; + LOG(info) << "[TrackTuner] qOverPtData = " << qOverPtData; // ############################################################################################## // ######## split the original string, separating substrings delimited by "|" symbol ######## @@ -216,6 +241,10 @@ struct TrackTuner { setBoolFromString(debugInfo, getValueString(DebugInfo)); LOG(info) << "[TrackTuner] debugInfo = " << debugInfo; outputString += "debugInfo=" + std::to_string(debugInfo); + // Configure updateTrackDCAs + setBoolFromString(updateTrackDCAs, getValueString(UpdateTrackDCAs)); + LOG(info) << "[TrackTuner] updateTrackDCAs = " << updateTrackDCAs; + outputString += ", updateTrackDCAs=" + std::to_string(updateTrackDCAs); // Configure updateTrackCovMat setBoolFromString(updateTrackCovMat, getValueString(UpdateTrackCovMat)); LOG(info) << "[TrackTuner] updateTrackCovMat = " << updateTrackCovMat; @@ -224,6 +253,10 @@ struct TrackTuner { setBoolFromString(updateCurvature, getValueString(UpdateCurvature)); LOG(info) << "[TrackTuner] updateCurvature = " << updateCurvature; outputString += ", updateCurvature=" + std::to_string(updateCurvature); + // Configure updateCurvatureIU + setBoolFromString(updateCurvatureIU, getValueString(UpdateCurvatureIU)); + LOG(info) << "[TrackTuner] updateCurvatureIU = " << updateCurvatureIU; + outputString += ", updateCurvatureIU=" + std::to_string(updateCurvatureIU); // Configure updatePulls setBoolFromString(updatePulls, getValueString(UpdatePulls)); LOG(info) << "[TrackTuner] updatePulls = " << updatePulls; @@ -236,22 +269,34 @@ struct TrackTuner { pathInputFile = getValueString(PathInputFile); outputString += ", pathInputFile=" + pathInputFile; LOG(info) << "[TrackTuner] pathInputFile = " << pathInputFile; + // Configure pathInputFile + pathFileQoverPt = getValueString(PathFileQoverPt); + outputString += ", pathFileQoverPt=" + pathFileQoverPt; + LOG(info) << "[TrackTuner] pathFileQoverPt = " << pathFileQoverPt; // Configure nameInputFile nameInputFile = getValueString(NameInputFile); outputString += ", nameInputFile=" + nameInputFile; LOG(info) << "[TrackTuner] nameInputFile = " << nameInputFile; + // Configure nameFileQoverPt + nameFileQoverPt = getValueString(NameFileQoverPt); + outputString += ", nameFileQoverPt=" + nameFileQoverPt; + LOG(info) << "[TrackTuner] nameFileQoverPt = " << nameFileQoverPt; // Configure usePvRefitCorrections setBoolFromString(usePvRefitCorrections, getValueString(UsePvRefitCorrections)); - outputString += ", usePvRefitCorrections=" + usePvRefitCorrections; + outputString += ", usePvRefitCorrections=" + std::to_string(usePvRefitCorrections); LOG(info) << "[TrackTuner] usePvRefitCorrections = " << usePvRefitCorrections; - // Configure oneOverPtCurr - oneOverPtCurrent = std::stof(getValueString(OneOverPtCurrent)); - outputString += ", oneOverPtCurrent=" + std::to_string(oneOverPtCurrent); - LOG(info) << "[TrackTuner] oneOverPtCurrent = " << oneOverPtCurrent; - // Configure oneOverPtUpgr - oneOverPtUpgr = std::stof(getValueString(OneOverPtUpgr)); - outputString += ", oneOverPtUpgr=" + std::to_string(oneOverPtUpgr); - LOG(info) << "[TrackTuner] oneOverPtUpgr = " << oneOverPtUpgr; + // Configure qOverPtMC + qOverPtMC = std::stof(getValueString(QOverPtMC)); + outputString += ", qOverPtMC=" + std::to_string(qOverPtMC); + LOG(info) << "[TrackTuner] qOverPtMC = " << qOverPtMC; + // Configure qOverPtData + qOverPtData = std::stof(getValueString(QOverPtData)); + outputString += ", qOverPtData=" + std::to_string(qOverPtData); + LOG(info) << "[TrackTuner] qOverPtData = " << qOverPtData; + + if ((updateCurvatureIU) && (updateCurvature)) { + LOG(fatal) << " [ updateCurvatureIU==kTRUE and updateCurvature==kTRUE ] -> Only one of them can be set to kTRUE at once! Please refer to the trackTuner documentation."; + } return outputString; } @@ -259,6 +304,7 @@ struct TrackTuner { void getDcaGraphs() { std::string fullNameInputFile = ""; + std::string fullNameFileQoverPt = ""; if (isInputFileFromCCDB) { /// use input correction file from CCDB @@ -267,23 +313,32 @@ struct TrackTuner { std::string tmpDir = "."; ccdbApi.init("http://alice-ccdb.cern.ch"); - // get the file from CCDB + // get the DCA correction file from CCDB if (!ccdbApi.retrieveBlob(pathInputFile.data(), tmpDir, metadata, 0, false, nameInputFile.data())) { - LOG(fatal) << "[TrackTuner] input file not found on CCDB, please check the pathInputFile and nameInputFile!"; + LOG(fatal) << "[TrackTuner] input file for DCA corrections not found on CCDB, please check the pathInputFile and nameInputFile!"; } + // get the Q/Pt correction file from CCDB + if (!ccdbApi.retrieveBlob(pathFileQoverPt.data(), tmpDir, metadata, 0, false, nameFileQoverPt.data())) { + LOG(fatal) << "[TrackTuner] input file for Q/Pt corrections not found on CCDB, please check the pathFileQoverPt and nameFileQoverPt!"; + } // point to the file in the tmp local folder fullNameInputFile = tmpDir + std::string("/") + nameInputFile; + fullNameFileQoverPt = tmpDir + std::string("/") + nameFileQoverPt; } else { /// use input correction file from local filesystem fullNameInputFile = pathInputFile + std::string("/") + nameInputFile; + fullNameFileQoverPt = pathFileQoverPt + std::string("/") + nameFileQoverPt; } - /// open the input correction file std::unique_ptr inputFile(TFile::Open(fullNameInputFile.c_str(), "READ")); if (!inputFile.get()) { LOG(fatal) << "Something wrong with the input file" << fullNameInputFile << " for dca correction. Fix it!"; } + std::unique_ptr inputFileQoverPt(TFile::Open(fullNameFileQoverPt.c_str(), "READ")); + if (!inputFileQoverPt.get() && (updateCurvature || updateCurvatureIU)) { + LOG(fatal) << "Something wrong with the Q/Pt input file" << fullNameFileQoverPt << " for Q/Pt correction. Fix it!"; + } // choose wheter to use corrections w/ PV refit or w/o it, and retrieve the proper TDirectory std::string dir = "woPvRefit"; @@ -295,330 +350,389 @@ struct TrackTuner { LOG(fatal) << "TDirectory " << td << " not found in input file" << inputFile->GetName() << ". Fix it!"; } - std::string grDcaXYResNameCurr = "resCurrentDcaXY"; - std::string grDcaXYMeanNameCurr = "meanCurrentDcaXY"; - std::string grDcaXYPullNameCurr = "pullsCurrentDcaXY"; - std::string grDcaXYResNameUpgr = "resUpgrDcaXY"; - std::string grDcaXYMeanNameUpgr = "meanUpgrDcaXY"; - std::string grDcaXYPullNameUpgr = "pullsUpgrDcaXY"; - - grDcaXYResVsPtPionCurrent.reset(dynamic_cast(td->Get(grDcaXYResNameCurr.c_str()))); - grDcaXYResVsPtPionUpgr.reset(dynamic_cast(td->Get(grDcaXYResNameUpgr.c_str()))); - grDcaXYMeanVsPtPionCurrent.reset(dynamic_cast(td->Get(grDcaXYMeanNameCurr.c_str()))); - grDcaXYMeanVsPtPionUpgr.reset(dynamic_cast(td->Get(grDcaXYMeanNameUpgr.c_str()))); - grDcaXYPullVsPtPionCurrent.reset(dynamic_cast(td->Get(grDcaXYPullNameCurr.c_str()))); - grDcaXYPullVsPtPionUpgr.reset(dynamic_cast(td->Get(grDcaXYPullNameUpgr.c_str()))); - if (!grDcaXYResVsPtPionCurrent.get() || !grDcaXYResVsPtPionUpgr.get() || !grDcaXYMeanVsPtPionCurrent.get() || !grDcaXYMeanVsPtPionUpgr.get() || !grDcaXYPullVsPtPionCurrent.get() || !grDcaXYPullVsPtPionUpgr.get()) { + std::string grDcaXYResNameMC = "resCurrentDcaXY"; + std::string grDcaXYMeanNameMC = "meanCurrentDcaXY"; + std::string grDcaXYPullNameMC = "pullsCurrentDcaXY"; + std::string grDcaXYResNameData = "resUpgrDcaXY"; + std::string grDcaXYMeanNameData = "meanUpgrDcaXY"; + std::string grDcaXYPullNameData = "pullsUpgrDcaXY"; + + grDcaXYResVsPtPionMC.reset(dynamic_cast(td->Get(grDcaXYResNameMC.c_str()))); + grDcaXYResVsPtPionData.reset(dynamic_cast(td->Get(grDcaXYResNameData.c_str()))); + grDcaXYMeanVsPtPionMC.reset(dynamic_cast(td->Get(grDcaXYMeanNameMC.c_str()))); + grDcaXYMeanVsPtPionData.reset(dynamic_cast(td->Get(grDcaXYMeanNameData.c_str()))); + grDcaXYPullVsPtPionMC.reset(dynamic_cast(td->Get(grDcaXYPullNameMC.c_str()))); + grDcaXYPullVsPtPionData.reset(dynamic_cast(td->Get(grDcaXYPullNameData.c_str()))); + if (!grDcaXYResVsPtPionMC.get() || !grDcaXYResVsPtPionData.get() || !grDcaXYMeanVsPtPionMC.get() || !grDcaXYMeanVsPtPionData.get() || !grDcaXYPullVsPtPionMC.get() || !grDcaXYPullVsPtPionData.get()) { LOG(fatal) << "Something wrong with the names of the correction graphs for dcaXY. Fix it!"; } - std::string grDcaZResNameCurr = "resCurrentDcaZ"; - std::string grDcaZMeanNameCurr = "meanCurrentDcaZ"; - std::string grDcaZPullNameCurr = "pullsCurrentDcaZ"; - std::string grDcaZResNameUpgr = "resUpgrDcaZ"; - std::string grDcaZMeanNameUpgr = "meanUpgrDcaZ"; - std::string grDcaZPullNameUpgr = "pullsUpgrDcaZ"; - - grDcaZResVsPtPionCurrent.reset(dynamic_cast(td->Get(grDcaZResNameCurr.c_str()))); - grDcaZResVsPtPionUpgr.reset(dynamic_cast(td->Get(grDcaZResNameUpgr.c_str()))); - grDcaZMeanVsPtPionCurrent.reset(dynamic_cast(td->Get(grDcaZMeanNameCurr.c_str()))); - grDcaZMeanVsPtPionUpgr.reset(dynamic_cast(td->Get(grDcaZMeanNameUpgr.c_str()))); - grDcaZPullVsPtPionCurrent.reset(dynamic_cast(td->Get(grDcaZPullNameCurr.c_str()))); - grDcaZPullVsPtPionUpgr.reset(dynamic_cast(td->Get(grDcaZPullNameUpgr.c_str()))); - if (!grDcaZResVsPtPionCurrent.get() || !grDcaZResVsPtPionUpgr.get() || !grDcaZMeanVsPtPionCurrent.get() || !grDcaZMeanVsPtPionUpgr.get() || !grDcaZPullVsPtPionCurrent.get() || !grDcaZPullVsPtPionUpgr.get()) { + std::string grDcaZResNameMC = "resCurrentDcaZ"; + std::string grDcaZMeanNameMC = "meanCurrentDcaZ"; + std::string grDcaZPullNameMC = "pullsCurrentDcaZ"; + std::string grDcaZResNameData = "resUpgrDcaZ"; + std::string grDcaZMeanNameData = "meanUpgrDcaZ"; + std::string grDcaZPullNameData = "pullsUpgrDcaZ"; + + grDcaZResVsPtPionMC.reset(dynamic_cast(td->Get(grDcaZResNameMC.c_str()))); + grDcaZResVsPtPionData.reset(dynamic_cast(td->Get(grDcaZResNameData.c_str()))); + grDcaZMeanVsPtPionMC.reset(dynamic_cast(td->Get(grDcaZMeanNameMC.c_str()))); + grDcaZMeanVsPtPionData.reset(dynamic_cast(td->Get(grDcaZMeanNameData.c_str()))); + grDcaZPullVsPtPionMC.reset(dynamic_cast(td->Get(grDcaZPullNameMC.c_str()))); + grDcaZPullVsPtPionData.reset(dynamic_cast(td->Get(grDcaZPullNameData.c_str()))); + if (!grDcaZResVsPtPionMC.get() || !grDcaZResVsPtPionData.get() || !grDcaZMeanVsPtPionMC.get() || !grDcaZMeanVsPtPionData.get() || !grDcaZPullVsPtPionMC.get() || !grDcaZPullVsPtPionData.get()) { LOG(fatal) << "Something wrong with the names of the correction graphs for dcaZ. Fix it!"; } - } + + std::string grOneOverPtPionNameMC = "sigmaVsPtMc"; + std::string grOneOverPtPionNameData = "sigmaVsPtData"; + + if (updateCurvature || updateCurvatureIU) { + grOneOverPtPionMC.reset(dynamic_cast(inputFileQoverPt->Get(grOneOverPtPionNameMC.c_str()))); + grOneOverPtPionData.reset(dynamic_cast(inputFileQoverPt->Get(grOneOverPtPionNameData.c_str()))); + } + } // getDcaGraphs() ends here template void tuneTrackParams(T1 const& mcparticle, T2& trackParCov, T3 const& matCorr, T4 dcaInfoCov, H hQA) { - double ptMC = mcparticle.pt(); - - double dcaXYResCurrent = 0.0; // sd0rpo=0.; - double dcaZResCurrent = 0.0; // sd0zo =0.; - - double dcaXYResUpgr = 0.0; // sd0rpn=0.; - double dcaZResUpgr = 0.0; // sd0zn =0.; - - // double OneOverPtCurrent = 0.0; // spt1o =0.; - // double OneOverPtUpgr = 0.0; // spt1n =0.; - - double dcaXYMeanCurrent = 0.0; // sd0mrpo=0.; - double dcaXYMeanUpgr = 0.0; // sd0mrpn=0.; - - double dcaXYPullCurrent = 1.0; - double dcaXYPullUpgr = 1.0; - - double dcaZPullCurrent = 1.0; - double dcaZPullUpgr = 1.0; - - dcaXYResCurrent = evalGraph(ptMC, grDcaXYResVsPtPionCurrent.get()); - dcaXYResUpgr = evalGraph(ptMC, grDcaXYResVsPtPionUpgr.get()); - - // dcaXYResCurrent = 1.0; - // dcaXYResUpgr = 1.5; - - dcaZResCurrent = evalGraph(ptMC, grDcaZResVsPtPionCurrent.get()); - dcaZResUpgr = evalGraph(ptMC, grDcaZResVsPtPionUpgr.get()); - - // dcaZResCurrent = 1.0; - // dcaZResUpgr = 1.0; - - // OneOverPtCurrent = evalGraph(ptMC, grOneOverPtPionCurrent.get() ); - // OneOverPtUpgr = evalGraph(ptMC, grOneOverPtPionUpgr.get() ); - - // OneOverPtCurrent = 1.0; - // OneOverPtUpgr = 2.0; - - dcaXYMeanCurrent = evalGraph(ptMC, grDcaXYMeanVsPtPionCurrent.get()); - dcaXYMeanUpgr = evalGraph(ptMC, grDcaXYMeanVsPtPionUpgr.get()); - - // dcaXYMeanCurrent = 0.0; - // dcaXYMeanUpgr = 0.0; - - dcaXYPullCurrent = evalGraph(ptMC, grDcaXYPullVsPtPionCurrent.get()); - dcaXYPullUpgr = evalGraph(ptMC, grDcaXYPullVsPtPionUpgr.get()); - - dcaZPullCurrent = evalGraph(ptMC, grDcaZPullVsPtPionCurrent.get()); - dcaZPullUpgr = evalGraph(ptMC, grDcaZPullVsPtPionUpgr.get()); - + double dcaXYResMC = 0.0; // sd0rpo=0.; + double dcaZResMC = 0.0; // sd0zo =0.; + + double dcaXYResData = 0.0; // sd0rpn=0.; + double dcaZResData = 0.0; // sd0zn =0.; + + double dcaXYMeanMC = 0.0; // sd0mrpo=0.; + double dcaXYMeanData = 0.0; // sd0mrpn=0.; + + double dcaXYPullMC = 1.0; + double dcaXYPullData = 1.0; + + double dcaZPullMC = 1.0; + double dcaZPullData = 1.0; + + dcaXYResMC = evalGraph(ptMC, grDcaXYResVsPtPionMC.get()); + dcaXYResData = evalGraph(ptMC, grDcaXYResVsPtPionData.get()); + + dcaZResMC = evalGraph(ptMC, grDcaZResVsPtPionMC.get()); + dcaZResData = evalGraph(ptMC, grDcaZResVsPtPionData.get()); + + // For Q/Pt corrections, files on CCDB will be used if both qOverPtMC and qOverPtData are null + if (updateCurvature || updateCurvatureIU) { + if ((qOverPtMC < 0) || (qOverPtData < 0)) { + if (debugInfo) { + LOG(info) << "### q/pt smearing: qOverPtMC=" << qOverPtMC << ", qOverPtData=" << qOverPtData << ". One of them is negative. Retrieving then values from graphs from input .root file"; + } + /// check that input graphs for q/pt smearing are correctly retrieved + if (!grOneOverPtPionData.get() || !grOneOverPtPionMC.get()) { + LOG(fatal) << "### q/pt smearing: input graphs not correctly retrieved. Aborting."; + } + qOverPtMC = std::max(0.0, evalGraph(ptMC, grOneOverPtPionMC.get())); + qOverPtData = std::max(0.0, evalGraph(ptMC, grOneOverPtPionData.get())); + } // qOverPtMC, qOverPtData block ends here + } // updateCurvature, updateCurvatureIU block ends here + + if (updateTrackDCAs) { + dcaXYMeanMC = evalGraph(ptMC, grDcaXYMeanVsPtPionMC.get()); + dcaXYMeanData = evalGraph(ptMC, grDcaXYMeanVsPtPionData.get()); + + dcaXYPullMC = evalGraph(ptMC, grDcaXYPullVsPtPionMC.get()); + dcaXYPullData = evalGraph(ptMC, grDcaXYPullVsPtPionData.get()); + + dcaZPullMC = evalGraph(ptMC, grDcaZPullVsPtPionMC.get()); + dcaZPullData = evalGraph(ptMC, grDcaZPullVsPtPionData.get()); + } // Unit conversion, is it required ?? - dcaXYResCurrent *= 1.e-4; - dcaZResCurrent *= 1.e-4; + dcaXYResMC *= 1.e-4; + dcaZResMC *= 1.e-4; - dcaXYResUpgr *= 1.e-4; - dcaZResUpgr *= 1.e-4; + dcaXYResData *= 1.e-4; + dcaZResData *= 1.e-4; - dcaXYMeanCurrent *= 1.e-4; - dcaXYMeanUpgr *= 1.e-4; + dcaXYMeanMC *= 1.e-4; + dcaXYMeanData *= 1.e-4; // Apply the smearing // --------------------------------------------- // double pt1o =param [4]; - double trackParOneOverPtCurrent = trackParCov.getQ2Pt(); + double trackParQPtMCRec = trackParCov.getQ2Pt(); int sign = trackParCov.getQ2Pt() / std::abs(trackParCov.getQ2Pt()); // double pt1mc =parammc[4]; - double trackParOneOverPtMC = sign / mcparticle.pt(); + double trackParQPtMC = sign / mcparticle.pt(); o2::dataformats::VertexBase vtxMC; vtxMC.setPos({mcparticle.vx(), mcparticle.vy(), mcparticle.vz()}); vtxMC.setCov(0, 0, 0, 0, 0, 0); // ??? or All ZEROs // == 1 cm2? wrt prop point if (debugInfo) { - // LOG(info) << " sign " << sign; - // LOG(info) << " trackParCov.getQ2Pt() " << trackParOneOverPtCurrent << " " << trackParCov.getQ2Pt(); - // LOG(info) << " sign/mcparticle.pt() " << trackParOneOverPtMC; - // LOG(info) << " (curvReco-curvMC)/curvMC " << (trackParOneOverPtCurrent - trackParOneOverPtMC)/trackParOneOverPtMC * 100.0 << "%"; - // LOG(info) << " trackParCov.getPtInv() " << trackParCov.getPtInv() << std::endl; - // LOG(info) << " 1/trackParCov.getPtInv() " << 1./trackParCov.getPtInv() << " & mcparticle.pt() " << mcparticle.pt(); - LOG(info) << mcparticle.pt() << " " << 1 / trackParCov.getPtInv() << " " << (trackParOneOverPtCurrent - trackParOneOverPtMC) / trackParOneOverPtMC * 100.0; - // LOG(info) << "Before Propagation to Production Point -> alpha: " << trackParCov.getAlpha() << ", DCAxy: " << trackParCov.getY() << ", DCAz: " << trackParCov.getZ(); - } - // propagate to DCA with respect to the Production point - o2::base::Propagator::Instance()->propagateToDCABxByBz(vtxMC, trackParCov, 2.f, matCorr, dcaInfoCov); - if (debugInfo) { - // LOG(info) << "After Propagation to Production Point -> alpha: " << trackParCov.getAlpha() << ", DCAxy: " << trackParCov.getY() << ", DCAz: " << trackParCov.getZ(); - LOG(info) << "track.y(): " << trackParCov.getY(); + LOG(info) << mcparticle.pt() << " " << 1 / trackParCov.getPtInv() << " " << (trackParQPtMCRec - trackParQPtMC) / trackParQPtMC * 100.0; } - ////////////////////////////// DCAs modifications Start ///////////////////////////////// - - // double d0zo =param [1]; - double trackParDcaZCurrent = trackParCov.getZ(); - - // double d0rpo =param [0]; - double trackParDcaXYCurrent = trackParCov.getY(); + // for updating the Q/pT from the tracksIU. + // This is placed before track propagation to production point, so that Q/pT can be updated before propagation + // Q/Pt is modified and set before track propagation - float mcVxRotated = mcparticle.vx() * std::cos(trackParCov.getAlpha()) + mcparticle.vy() * std::sin(trackParCov.getAlpha()); // invert - float mcVyRotated = mcparticle.vy() * std::cos(trackParCov.getAlpha()) - mcparticle.vx() * std::sin(trackParCov.getAlpha()); - - if (debugInfo) { - // LOG(info) << "mcVy " << mcparticle.vy() << std::endl; - LOG(info) << "mcVxRotated " << mcVxRotated; - LOG(info) << "mcVyRotated " << mcVyRotated; - } + double deltaQpt = 0.0; + double deltaQptTuned = 0.0; + double trackParQPtTuned = 0.0; - // std::array arrayXYZ = { mcVxRotated, mcVyRotated , mcparticle.vz()}; - // std::array arrayPxPyPz = {mcparticle.px(), mcparticle.py(), mcparticle.pz()}; + // variables for track cov matrix elements update + double sigmaY2 = 0.0; + double sigmaZY = 0.0; + double sigmaZ2 = 0.0; + double sigmaSnpY = 0.0; + double sigmaSnpZ = 0.0; + double sigmaTglY = 0.0; + double sigmaTglZ = 0.0; + double sigma1PtY = 0.0; + double sigma1PtZ = 0.0; + double sigma1PtSnp = 0.0; + double sigma1PtTgl = 0.0; + double sigma1Pt2 = 0.0; - // const int matSize = 21; - // std::array arrayCovMatMC = {0.,0.,0.,0.,0., 0.,0.,0.,0.,0., 0.,0.,0.,0.,0., 0.,0.,0.,0.,0.,0. }; - // o2::track::TrackParametrizationWithError trackParCovMC{arrayXYZ, arrayPxPyPz, arrayCovMatMC, trackParCov.getSign()}; + double sigmaY2orig = trackParCov.getSigmaY2(); + double sigmaZYorig = trackParCov.getSigmaZY(); + double sigmaZ2orig = trackParCov.getSigmaZ2(); - // double d0rpmc=parammc[0]; - // double trackParDcaXYMC = trackParCovMC.getY(); // here + if (updateCurvatureIU) { + // double dpt1o =pt1o-pt1mc; + deltaQpt = trackParQPtMCRec - trackParQPtMC; + // double dpt1n =dpt1o *(spt1o >0. ? (spt1n /spt1o ) : 1.); + deltaQptTuned = deltaQpt * (qOverPtMC > 0. ? (qOverPtData / qOverPtMC) : 1.); + // double pt1n = pt1mc+dpt1n; + trackParQPtTuned = trackParQPtMC + deltaQptTuned; + trackParCov.setQ2Pt(trackParQPtTuned); - double trackParDcaXYMC = mcVyRotated; // here - // LOG(info) << "trackParCovMC.getY() " << trackParCovMC.getY() << std::endl; + // updating track cov matrix elements for 1/Pt at innermost update point + // if(sd0rpo>0. && spt1o>0.)covar[10]*=(sd0rpn/sd0rpo)*(spt1n/spt1o);//ypt + sigma1PtY = trackParCov.getSigma1PtY(); + if (dcaXYResMC > 0. && qOverPtMC > 0.) { + sigma1PtY *= ((dcaXYResData / dcaXYResMC) * (qOverPtData / qOverPtMC)); + trackParCov.setCov(sigma1PtY, 10); + } - // double d0zmc =parammc[1]; - double trackParDcaZMC = mcparticle.vz(); + // if(sd0zo>0. && spt1o>0.) covar[11]*=(sd0zn/sd0zo)*(spt1n/spt1o);//zpt + sigma1PtZ = trackParCov.getSigma1PtZ(); + if (dcaZResMC > 0. && qOverPtMC > 0.) { + sigma1PtZ *= ((dcaZResData / dcaZResMC) * (qOverPtData / qOverPtMC)); + trackParCov.setCov(sigma1PtZ, 11); + } - // double dd0zo =d0zo-d0zmc; - double diffDcaZFromMCCurent = trackParDcaZCurrent - trackParDcaZMC; + // if(spt1o>0.) covar[12]*=(spt1n/spt1o);//sinPhipt + sigma1PtSnp = trackParCov.getSigma1PtSnp(); + if (qOverPtMC > 0.) { + sigma1PtSnp *= (qOverPtData / qOverPtMC); + trackParCov.setCov(sigma1PtSnp, 12); + } - // double dd0zn =dd0zo *(sd0zo >0. ? (sd0zn /sd0zo ) : 1.); - double diffDcaZFromMCUpgr = diffDcaZFromMCCurent * (dcaZResCurrent > 0. ? (dcaZResUpgr / dcaZResCurrent) : 1.); + // if(spt1o>0.) covar[13]*=(spt1n/spt1o);//tanTpt + sigma1PtTgl = trackParCov.getSigma1PtTgl(); + if (qOverPtMC > 0.) { + sigma1PtTgl *= (qOverPtData / qOverPtMC); + trackParCov.setCov(sigma1PtTgl, 13); + } - // double d0zn =d0zmc+dd0zn; - double trackParDcaZUpgr = trackParDcaZMC + diffDcaZFromMCUpgr; + // if(spt1o>0.) covar[14]*=(spt1n/spt1o)*(spt1n/spt1o);//ptpt + sigma1Pt2 = trackParCov.getSigma1Pt2(); + if (qOverPtMC > 0.) { + sigma1Pt2 *= (qOverPtData / qOverPtMC); + trackParCov.setCov(sigma1Pt2, 14); + } + } // updateCurvatureIU block ends here + // propagate to DCA with respect to the Production point + // if (!updateCurvatureIU) { + // o2::base::Propagator::Instance()->propagateToDCABxByBz(vtxMC, trackParCov, 2.f, matCorr, dcaInfoCov); + // } - // double dd0rpo=d0rpo-d0rpmc; - double diffDcaXYFromMCCurent = trackParDcaXYCurrent - trackParDcaXYMC; + double trackParDcaXYoriginal = trackParCov.getY(); + double trackParDcaZoriginal = trackParCov.getZ(); - // double dd0rpn=dd0rpo*(sd0rpo>0. ? (sd0rpn/sd0rpo) : 1.); - double diffDcaXYFromMCUpgr = diffDcaXYFromMCCurent * (dcaXYResCurrent > 0. ? (dcaXYResUpgr / dcaXYResCurrent) : 1.); + if (updateTrackDCAs) { + // propagate to DCA with respect to the Production point + o2::base::Propagator::Instance()->propagateToDCABxByBz(vtxMC, trackParCov, 2.f, matCorr, dcaInfoCov); + // double d0zo =param [1]; + double trackParDcaZRec = trackParCov.getZ(); + // double d0rpo =param [0]; + double trackParDcaXYRec = trackParCov.getY(); + float mcVxRotated = mcparticle.vx() * std::cos(trackParCov.getAlpha()) + mcparticle.vy() * std::sin(trackParCov.getAlpha()); // invert + float mcVyRotated = mcparticle.vy() * std::cos(trackParCov.getAlpha()) - mcparticle.vx() * std::sin(trackParCov.getAlpha()); + if (debugInfo) { + // LOG(info) << "mcVy " << mcparticle.vy() << std::endl; + LOG(info) << "mcVxRotated " << mcVxRotated; + LOG(info) << "mcVyRotated " << mcVyRotated; + } + // std::array arrayXYZ = { mcVxRotated, mcVyRotated , mcparticle.vz()}; + // std::array arrayPxPyPz = {mcparticle.px(), mcparticle.py(), mcparticle.pz()}; - // double dd0mrpn=std::abs(sd0mrpn)-std::abs(sd0mrpo); - // double diffDcaXYMeanUpgMinusCur = std::abs(dcaXYMeanUpgr) - std::abs(dcaXYMeanCurrent) ; - double diffDcaXYMeanUpgMinusCur = dcaXYMeanUpgr - dcaXYMeanCurrent; + // const int matSize = 21; + // std::array arrayCovMatMC = {0.,0.,0.,0.,0., 0.,0.,0.,0.,0., 0.,0.,0.,0.,0., 0.,0.,0.,0.,0.,0. }; + // o2::track::TrackParametrizationWithError trackParCovMC{arrayXYZ, arrayPxPyPz, arrayCovMatMC, trackParCov.getSign()}; - // double d0rpn =d0rpmc+dd0rpn-dd0mrpn; - double trackParDcaXYUpgr = trackParDcaXYMC + diffDcaXYFromMCUpgr - diffDcaXYMeanUpgMinusCur; + // double d0rpmc=parammc[0]; + // double trackParDcaXYMC = trackParCovMC.getY(); // here - if (debugInfo) { - LOG(info) << dcaZResCurrent << ", " << dcaZResUpgr << ", diff(DcaZ - DcaZMC): " << diffDcaZFromMCCurent << ", diff upgraded: " << diffDcaZFromMCUpgr << ", DcaZ Upgr : " << trackParDcaZUpgr; - LOG(info) << dcaXYResCurrent << ", " << dcaXYResUpgr << ", diff(DcaY - DcaYMC): " << diffDcaXYFromMCCurent << ", diff upgraded: " << diffDcaXYFromMCUpgr << ", DcaY Upgr :" << trackParDcaXYUpgr; - } + double trackParDcaXYMC = mcVyRotated; // here - // option mimic data - // ---------------------- - // if(fMimicData){ - // // dd0mrpn=sd0mrpn-sd0mrpo; - // diffDcaXYMeanUpgMinusCur = dcaXYMeanUpgr - dcaXYMeanCurrent; - // // d0rpn = d0rpmc+dd0rpn+dd0mrpn; - // trackParDcaXYUpgr = diffDcaXYFromMCCurent + diffDcaXYFromMCUpgr + diffDcaXYMeanUpgMinusCur; - // } + // double d0zmc =parammc[1]; + double trackParDcaZMC = mcparticle.vz(); - // setting updated track parameters - // -------------------------------- - // param[0]=d0rpn; - // double oldDCAxyValue = trackParCov.getY(); - double trackParDcaXYoriginal = trackParCov.getY(); - trackParCov.setY(trackParDcaXYUpgr); - // trackParCov.setY(oldDCAxyValue); - // param[1]=d0zn ; - double trackParDcaZoriginal = trackParCov.getZ(); - trackParCov.setZ(trackParDcaZUpgr); + // double dd0zo =d0zo-d0zmc; + double deltaDcaZ = trackParDcaZRec - trackParDcaZMC; - if (updateCurvature) { - // -------------------------------------- - // double dpt1o =pt1o-pt1mc; - double diffOneOverPtFromMCCurent = trackParOneOverPtCurrent - trackParOneOverPtMC; + // double dd0zn =dd0zo *(sd0zo >0. ? (sd0zn /sd0zo ) : 1.); + double deltaDcaZTuned = deltaDcaZ * (dcaZResMC > 0. ? (dcaZResData / dcaZResMC) : 1.); - // double dpt1n =dpt1o *(spt1o >0. ? (spt1n /spt1o ) : 1.); - double diffOneOverPtFromMCUpgr = diffOneOverPtFromMCCurent * (oneOverPtCurrent > 0. ? (oneOverPtUpgr / oneOverPtCurrent) : 1.); + // double d0zn =d0zmc+dd0zn; + double trackParDcaZTuned = trackParDcaZMC + deltaDcaZTuned; - // double pt1n = pt1mc+dpt1n; - double trackParOneOverPtUpgr = trackParOneOverPtMC + diffOneOverPtFromMCUpgr; + // double dd0rpo=d0rpo-d0rpmc; + double deltaDcaXY = trackParDcaXYRec - trackParDcaXYMC; - // param[4]=pt1n ; - trackParCov.setQ2Pt(trackParOneOverPtUpgr); - } - // if(debugInfo){ - // LOG(info) << "Inside tuneTrackParams() before modifying trackParCov.getY(): " << trackParCov.getY() << " trackParOneOverPtMC = " << trackParOneOverPtMC << " diffOneOverPtFromMCUpgr = " << diffOneOverPtFromMCUpgr << std::endl; - //} + // double dd0rpn=dd0rpo*(sd0rpo>0. ? (sd0rpn/sd0rpo) : 1.); + double deltaDcaXYTuned = deltaDcaXY * (dcaXYResMC > 0. ? (dcaXYResData / dcaXYResMC) : 1.); - // Updating Single Track Covariance matrices + // double dd0mrpn=std::abs(sd0mrpn)-std::abs(sd0mrpo); + // double deltaDcaXYmean = std::abs(dcaXYMeanData) - std::abs(dcaXYMeanMC) ; + double deltaDcaXYmean = dcaXYMeanData - dcaXYMeanMC; - double sigmaY2 = 0.0; - double sigmaZY = 0.0; - double sigmaZ2 = 0.0; - double sigmaSnpY = 0.0; - double sigmaSnpZ = 0.0; - double sigmaTglY = 0.0; - double sigmaTglZ = 0.0; - double sigma1PtY = 0.0; - double sigma1PtZ = 0.0; - double sigma1PtSnp = 0.0; - double sigma1PtTgl = 0.0; - double sigma1Pt2 = 0.0; + // double d0rpn =d0rpmc+dd0rpn-dd0mrpn; + double trackParDcaXYTuned = trackParDcaXYMC + deltaDcaXYTuned - deltaDcaXYmean; - double sigmaY2orig = trackParCov.getSigmaY2(); - double sigmaZYorig = trackParCov.getSigmaZY(); - double sigmaZ2orig = trackParCov.getSigmaZ2(); + if (debugInfo) { + LOG(info) << dcaZResMC << ", " << dcaZResData << ", diff(DcaZ - DcaZMC): " << deltaDcaZ << ", diff upgraded: " << deltaDcaZTuned << ", DcaZ Data : " << trackParDcaZTuned; + LOG(info) << dcaXYResMC << ", " << dcaXYResData << ", diff(DcaY - DcaYMC): " << deltaDcaXY << ", diff upgraded: " << deltaDcaXYTuned << ", DcaY Data :" << trackParDcaXYTuned; + } + // option mimic data + // ---------------------- + // if(fMimicData){ + // // dd0mrpn=sd0mrpn-sd0mrpo; + // deltaDcaXYmean = dcaXYMeanData - dcaXYMeanMC; + // // d0rpn = d0rpmc+dd0rpn+dd0mrpn; + // trackParDcaXYTuned = deltaDcaXY + deltaDcaXYTuned + deltaDcaXYmean; + // } + + // setting updated track parameters + // -------------------------------- + trackParDcaXYoriginal = trackParCov.getY(); + trackParCov.setY(trackParDcaXYTuned); + trackParDcaZoriginal = trackParCov.getZ(); + trackParCov.setZ(trackParDcaZTuned); + } // ----> updateTrackDCAs block ends here + + if ((updateCurvature) && (!updateCurvatureIU)) { // ...block begins here + if (!updateTrackDCAs) { + /// propagation to production point not done yet, doing it now + o2::base::Propagator::Instance()->propagateToDCABxByBz(vtxMC, trackParCov, 2.f, matCorr, dcaInfoCov); + } + deltaQpt = trackParQPtMCRec - trackParQPtMC; + // double dpt1n =dpt1o *(spt1o >0. ? (spt1n /spt1o ) : 1.); + deltaQptTuned = deltaQpt * (qOverPtMC > 0. ? (qOverPtData / qOverPtMC) : 1.); + // double pt1n = pt1mc+dpt1n; + trackParQPtTuned = trackParQPtMC + deltaQptTuned; + trackParCov.setQ2Pt(trackParQPtTuned); + } // ...block ends here if (updateTrackCovMat) { // if(sd0rpo>0.) covar[0]*=(sd0rpn/sd0rpo)*(sd0rpn/sd0rpo);//yy sigmaY2 = trackParCov.getSigmaY2(); - if (dcaXYResCurrent > 0.) - sigmaY2 *= ((dcaXYResUpgr / dcaXYResCurrent) * (dcaXYResUpgr / dcaXYResCurrent)); - trackParCov.setCov(sigmaY2, 0); + if (dcaXYResMC > 0.) { + sigmaY2 *= ((dcaXYResData / dcaXYResMC) * (dcaXYResData / dcaXYResMC)); + trackParCov.setCov(sigmaY2, 0); + } // if(sd0zo>0. && sd0rpo>0.)covar[1]*=(sd0rpn/sd0rpo)*(sd0zn/sd0zo);//yz sigmaZY = trackParCov.getSigmaZY(); - if (dcaZResCurrent > 0. && dcaXYResCurrent > 0.) - sigmaZY *= ((dcaXYResUpgr / dcaXYResCurrent) * (dcaZResUpgr / dcaZResCurrent)); - trackParCov.setCov(sigmaZY, 1); + if (dcaZResMC > 0. && dcaXYResMC > 0.) { + // sigmaZY *= ((dcaXYResData / dcaXYResMC) * (dcaZResData / dcaZResMC)); + sigmaZY *= ((dcaXYResData / dcaXYResMC) * (dcaZResData / dcaZResMC)); + trackParCov.setCov(sigmaZY, 1); + } // if(sd0zo>0.) covar[2]*=(sd0zn/sd0zo)*(sd0zn/sd0zo);//zz sigmaZ2 = trackParCov.getSigmaZ2(); - if (dcaZResCurrent > 0.) - sigmaZ2 *= ((dcaZResUpgr / dcaZResCurrent) * (dcaZResUpgr / dcaZResCurrent)); - trackParCov.setCov(sigmaZ2, 2); + if (dcaZResMC > 0.) { + sigmaZ2 *= ((dcaZResData / dcaZResMC) * (dcaZResData / dcaZResMC)); + trackParCov.setCov(sigmaZ2, 2); + } // if(sd0rpo>0.) covar[3]*=(sd0rpn/sd0rpo);//yl sigmaSnpY = trackParCov.getSigmaSnpY(); - if (dcaXYResCurrent > 0.) - sigmaSnpY *= ((dcaXYResUpgr / dcaXYResCurrent)); - trackParCov.setCov(sigmaSnpY, 3); + if (dcaXYResMC > 0.) { + sigmaSnpY *= ((dcaXYResData / dcaXYResMC)); + trackParCov.setCov(sigmaSnpY, 3); + } // if(sd0zo>0.) covar[4]*=(sd0zn/sd0zo);//zl sigmaSnpZ = trackParCov.getSigmaSnpZ(); - if (dcaZResCurrent > 0.) - sigmaSnpZ *= ((dcaZResUpgr / dcaZResCurrent)); - trackParCov.setCov(sigmaSnpZ, 4); + if (dcaZResMC > 0.) { + sigmaSnpZ *= ((dcaZResData / dcaZResMC)); + trackParCov.setCov(sigmaSnpZ, 4); + } // if(sd0rpo>0.) covar[6]*=(sd0rpn/sd0rpo);//ysenT sigmaTglY = trackParCov.getSigmaTglY(); - if (dcaXYResCurrent > 0.) - sigmaTglY *= ((dcaXYResUpgr / dcaXYResCurrent)); - trackParCov.setCov(sigmaTglY, 6); + if (dcaXYResMC > 0.) { + sigmaTglY *= ((dcaXYResData / dcaXYResMC)); + trackParCov.setCov(sigmaTglY, 6); + } // if(sd0zo>0.) covar[7]*=(sd0zn/sd0zo);//zsenT sigmaTglZ = trackParCov.getSigmaTglZ(); - if (dcaZResCurrent > 0.) - sigmaTglZ *= ((dcaZResUpgr / dcaZResCurrent)); - trackParCov.setCov(sigmaTglZ, 7); - - // if(sd0rpo>0. && spt1o>0.)covar[10]*=(sd0rpn/sd0rpo)*(spt1n/spt1o);//ypt - sigma1PtY = trackParCov.getSigma1PtY(); - if (dcaXYResCurrent > 0. && oneOverPtCurrent > 0.) - sigma1PtY *= ((dcaXYResUpgr / dcaXYResCurrent) * (oneOverPtUpgr / oneOverPtCurrent)); - trackParCov.setCov(sigma1PtY, 10); - - // if(sd0zo>0. && spt1o>0.) covar[11]*=(sd0zn/sd0zo)*(spt1n/spt1o);//zpt - sigma1PtZ = trackParCov.getSigma1PtZ(); - if (dcaZResCurrent > 0. && oneOverPtCurrent > 0.) - sigma1PtZ *= ((dcaZResUpgr / dcaZResCurrent) * (oneOverPtUpgr / oneOverPtCurrent)); - trackParCov.setCov(sigma1PtZ, 11); - - // if(spt1o>0.) covar[12]*=(spt1n/spt1o);//sinPhipt - sigma1PtSnp = trackParCov.getSigma1PtSnp(); - if (oneOverPtCurrent > 0.) - sigma1PtSnp *= (oneOverPtUpgr / oneOverPtCurrent); - trackParCov.setCov(sigma1PtSnp, 12); - - // if(spt1o>0.) covar[13]*=(spt1n/spt1o);//tanTpt - sigma1PtTgl = trackParCov.getSigma1PtTgl(); - if (oneOverPtCurrent > 0.) - sigma1PtTgl *= (oneOverPtUpgr / oneOverPtCurrent); - trackParCov.setCov(sigma1PtTgl, 13); + if (dcaZResMC > 0.) { + sigmaTglZ *= ((dcaZResData / dcaZResMC)); + trackParCov.setCov(sigmaTglZ, 7); + } - // if(spt1o>0.) covar[14]*=(spt1n/spt1o)*(spt1n/spt1o);//ptpt - sigma1Pt2 = trackParCov.getSigma1Pt2(); - if (oneOverPtCurrent > 0.) - sigma1Pt2 *= (oneOverPtUpgr / oneOverPtCurrent); - trackParCov.setCov(sigma1Pt2, 14); - } + // checking and updating track cov matrix elements for 1/Pt begins + if ((updateCurvature) && (!updateCurvatureIU)) { + // if(sd0rpo>0. && spt1o>0.)covar[10]*=(sd0rpn/sd0rpo)*(spt1n/spt1o);//ypt + sigma1PtY = trackParCov.getSigma1PtY(); + if (dcaXYResMC > 0. && qOverPtMC > 0.) { + sigma1PtY *= ((dcaXYResData / dcaXYResMC) * (qOverPtData / qOverPtMC)); + trackParCov.setCov(sigma1PtY, 10); + } + + // if(sd0zo>0. && spt1o>0.) covar[11]*=(sd0zn/sd0zo)*(spt1n/spt1o);//zpt + sigma1PtZ = trackParCov.getSigma1PtZ(); + if (dcaZResMC > 0. && qOverPtMC > 0.) { + sigma1PtZ *= ((dcaZResData / dcaZResMC) * (qOverPtData / qOverPtMC)); + trackParCov.setCov(sigma1PtZ, 11); + } + + // if(spt1o>0.) covar[12]*=(spt1n/spt1o);//sinPhipt + sigma1PtSnp = trackParCov.getSigma1PtSnp(); + if (qOverPtMC > 0.) { + sigma1PtSnp *= (qOverPtData / qOverPtMC); + trackParCov.setCov(sigma1PtSnp, 12); + } + + // if(spt1o>0.) covar[13]*=(spt1n/spt1o);//tanTpt + sigma1PtTgl = trackParCov.getSigma1PtTgl(); + if (qOverPtMC > 0.) { + sigma1PtTgl *= (qOverPtData / qOverPtMC); + trackParCov.setCov(sigma1PtTgl, 13); + } + + // if(spt1o>0.) covar[14]*=(spt1n/spt1o)*(spt1n/spt1o);//ptpt + sigma1Pt2 = trackParCov.getSigma1Pt2(); + if (qOverPtMC > 0.) { + sigma1Pt2 *= (qOverPtData / qOverPtMC); + trackParCov.setCov(sigma1Pt2, 14); + } + } // ---> track cov matrix elements for 1/Pt ends here + } // ---> updateTrackCovMat block ends here if (updatePulls) { - double ratioDCAxyPulls = dcaXYPullCurrent / dcaXYPullUpgr; - double ratioDCAzPulls = dcaZPullCurrent / dcaZPullUpgr; + double ratioDCAxyPulls = 1.0; + double ratioDCAzPulls = 1.0; + if (dcaZPullData > 0.0) { + ratioDCAzPulls = dcaZPullMC / dcaZPullData; + } + if (dcaXYPullData > 0.0) { + ratioDCAxyPulls = dcaXYPullMC / dcaXYPullData; + } // covar[0]*=pullcorr*pullcorr;//yy - sigmaY2 *= (ratioDCAxyPulls * ratioDCAxyPulls); trackParCov.setCov(sigmaY2, 0); @@ -649,7 +763,7 @@ struct TrackTuner { sigma1PtZ *= ratioDCAzPulls; trackParCov.setCov(sigma1PtZ, 11); - } + } // ---> updatePulls block ends here /// sanity check for track covariance matrix element /// see https://github.com/AliceO2Group/AliceO2/blob/66de30958153cd7badf522150e8554f9fcf975ff/Common/DCAFitter/include/DCAFitter/DCAFitterN.h#L38-L54 @@ -689,7 +803,7 @@ struct TrackTuner { } else { hQA->Fill(2); } - } + } // tuneTrackParams() ends here // to be declared // ---------------