From 6a44babc7b930c13c93980c44806a08c94c42d2e Mon Sep 17 00:00:00 2001 From: Tanu Gahlaut <154991749+TGahlaut1@users.noreply.github.com> Date: Fri, 28 Aug 2026 10:48:56 +0530 Subject: [PATCH 1/3] update correction method and remove extra matrices histograms --- PWGCF/EbyEFluctuations/Tasks/meanPtFlucId.cxx | 1216 ++++++++++------- 1 file changed, 722 insertions(+), 494 deletions(-) diff --git a/PWGCF/EbyEFluctuations/Tasks/meanPtFlucId.cxx b/PWGCF/EbyEFluctuations/Tasks/meanPtFlucId.cxx index 7665f440bdb..b704037695f 100644 --- a/PWGCF/EbyEFluctuations/Tasks/meanPtFlucId.cxx +++ b/PWGCF/EbyEFluctuations/Tasks/meanPtFlucId.cxx @@ -24,6 +24,7 @@ #include "Common/DataModel/PIDResponseTPC.h" #include "Common/DataModel/TrackSelectionTables.h" +#include "CCDB/BasicCCDBManager.h" #include #include #include @@ -37,10 +38,15 @@ #include #include +#include +#include #include +#include #include +#include #include +#include #include #include #include @@ -53,34 +59,37 @@ using namespace std; struct MeanPtFlucId { Configurable nPBins{"nPBins", 300, ""}; - Configurable nPartBins{"nPartBins", 100, ""}; Configurable nPhiBins{"nPhiBins", 100, ""}; + Configurable nPtBins{"nPtBins", 50, "Number of pT bins"}; Configurable cfgCutPtMax{"cfgCutPtMax", 2.0, "maximum pT"}; Configurable cfgCutPtMin{"cfgCutPtMin", 0.2, "minimum pT"}; Configurable cfgCutEta{"cfgCutEta", 0.8, "Eta cut"}; - Configurable cfgCutDcaZ{"cfgCutDcaZ", 0.15, "DCAz cut"}; + Configurable cfgDcaPtDepOffset{"cfgDcaPtDepOffset", 0.004, "Offset of the pT-dependent DCA cut (cm)"}; + Configurable cfgDcaPtDepScale{"cfgDcaPtDepScale", 0.013, "1/pT coefficient of the pT-dependent DCA cut (cm GeV/c)"}; Configurable cfgCutPosZ{"cfgCutPosZ", 7.0, "cut for vertex Z"}; Configurable cfgPosZ{"cfgPosZ", true, "Position Z"}; Configurable cfgCutNSig2{"cfgCutNSig2", 2.0, "nSigma cut: 2"}; Configurable cfgCutNSig3{"cfgCutNSig3", 3.0, "nSigma cut: 3"}; Configurable cfgCutNSig5{"cfgCutNSig5", 5.0, "nSigma cut: 5"}; - Configurable cfgSelCutNSigPi{"cfgSelCutNSigPi", 2.0, "nSigma cut for pion selection"}; - Configurable cfgSelCutNSigKa{"cfgSelCutNSigKa", 3.0, "nSigma cut for kaon selection"}; - Configurable cfgSelCutNSigPr{"cfgSelCutNSigPr", 2.0, "nSigma cut for proton selection"}; + Configurable cfgCutNSigTpcTof{"cfgCutNSigTpcTof", 3.0, "nSigma TPC cut when TPC + TOF "}; + Configurable cfgSelCutNSigTpcPi{"cfgSelCutNSigTpcPi", 2.0, "nSigma TPC cut for pion selection"}; + Configurable cfgSelCutNSigTpcKa{"cfgSelCutNSigTpcKa", 2.0, "nSigma TPC cut for kaon selection"}; + Configurable cfgSelCutNSigTpcPr{"cfgSelCutNSigTpcPr", 2.0, "nSigma TPC cut for proton selection"}; + Configurable cfgSelCutNSigTofPi{"cfgSelCutNSigTofPi", 2.0, "nSigma TOF cut for pion selection"}; + Configurable cfgSelCutNSigTofKa{"cfgSelCutNSigTofKa", 2.0, "nSigma TOF cut for kaon selection"}; + Configurable cfgSelCutNSigTofPr{"cfgSelCutNSigTofPr", 2.0, "nSigma TOF cut for proton selection"}; Configurable cfgCutPiPtMin{"cfgCutPiPtMin", 0.2, "Minimum pion p_{T} cut"}; Configurable cfgCutKaPtMin{"cfgCutKaPtMin", 0.3, "Minimum kaon p_{T} cut"}; Configurable cfgCutPrPtMin{"cfgCutPrPtMin", 0.5, "Minimum proton p_{T} cut"}; - Configurable cfgCutPiThrsldP{"cfgCutPiThrsldP", 0.6, "Threshold p cut pion"}; - Configurable cfgCutKaThrsldP{"cfgCutKaThrsldP", 0.6, "Threshold p cut kaon"}; - Configurable cfgCutPrThrsldP{"cfgCutPrThrsldP", 1.0, "Threshold p cut proton "}; + Configurable cfgCutPiThrsldP{"cfgCutPiThrsldP", 0.7, "Momentum threshold for requiring pion TOF PID"}; + Configurable cfgCutKaThrsldP{"cfgCutKaThrsldP", 0.8, "Momentum threshold for requiring kaon TOF PID"}; + Configurable cfgCutPrThrsldP{"cfgCutPrThrsldP", 0.8, "Momentum threshold for requiring proton TOF PID"}; Configurable cfgP0U{"cfgP0U", 6.36269, "p_{0} for upper cut for N_{TPC} vs N{sim}"}; Configurable cfgP1U{"cfgP1U", 2.22108, "p_{1} for upper cut for N_{TPC} vs N{sim}"}; Configurable cfgP2U{"cfgP2U", -0.0118451, "p_{2} for upper cut for N_{TPC} vs N{sim}"}; Configurable cfgP0L{"cfgP0L", 23.8213, "p_{0} for lower cut for N_{TPC} vs N{sim}"}; Configurable cfgP1L{"cfgP1L", -0.402328, "p_{1} for lower cut for N_{TPC} vs N{sim}"}; Configurable cfgP2L{"cfgP2L", 0.0128423, "p_{2} for lower cut for N_{TPC} vs N{sim}"}; - Configurable cfgP0L1{"cfgP0L1", -10.5149, "p_{0} for lower cut for N_{TPC} vs N{sim}"}; - Configurable cfgP1L1{"cfgP1L1", 0.936596, "p_{1} for lower cut for N_{TPC} vs N{sim}"}; Configurable cfgP0corr{"cfgP0corr", 0.0001, "p_{0} for corrected N_{TPC} "}; Configurable cfgP1corr{"cfgP1corr", 0.0001, "p_{1} for corrected N_{TPC} "}; Configurable cfgP2corr{"cfgP2corr", 0.0001, "p_{2} for corrected N_{TPC} "}; @@ -92,25 +101,32 @@ struct MeanPtFlucId { Configurable cfgIsGoodZvtxFT0vsPV{"cfgIsGoodZvtxFT0vsPV", true, "kIsGoodZvtxFT0vsPV"}; Configurable cfgRejTrk{"cfgRejTrk", true, "Rejected Tracks"}; Configurable cfgNtpcEventCut{"cfgNtpcEventCut", true, "N_{TPC} Event Cut for reco"}; + Configurable cfgApplyEfficiencyCorrection{"cfgApplyEfficiencyCorrection", false, "Apply measured track-by-track efficiency weights"}; + Configurable cfgUseCombinedPurityEfficiencyWeight{"cfgUseCombinedPurityEfficiencyWeight", false, "Use the legacy Ngen/Nselected weight instead of the efficiency-only 1/epsilon weight"}; + Configurable cfgEfficiencyCCDBUrl{"cfgEfficiencyCCDBUrl", "http://ccdb-test.cern.ch:8080", "Efficiency CCDB URL"}; + Configurable cfgEfficiencyCCDBPath{"cfgEfficiencyCCDBPath", "Users/t/tgahlaut/weightCorr", "CCDB path containing ccdb_object"}; + Configurable cfgEfficiencyTimestamp{"cfgEfficiencyTimestamp", -1, "Efficiency CCDB timestamp; -1 selects the latest object"}; ConfigurableAxis multTPCBins{"multTPCBins", {150, 0, 150}, "TPC Multiplicity bins"}; ConfigurableAxis multMCBins{"multMCBins", {300, 0, 300}, "MC Multiplicity bins"}; ConfigurableAxis multCorrBins{"multCorrBins", {100, 0., 150.}, "Corrected TPC Multiplicity bins"}; ConfigurableAxis multFT0MBins{"multFT0MBins", {1000, 0, 5000}, "Forward Multiplicity bins"}; ConfigurableAxis dcaXYBins{"dcaXYBins", {100, -0.15, 0.15}, "dcaXY bins"}; ConfigurableAxis dcaZBins{"dcaZBins", {500, -1.2, 1.2}, "dcaZ bins"}; - ConfigurableAxis centBins{"centBins", {0, 5, 10, 15, 20, 25, 30, 35, 40, 45, 50, 55, 60, 65, 70, 75, 80, 85, 90, 95, 100}, "Centrality bins"}; + ConfigurableAxis centBins{"centBins", {VARIABLE_WIDTH, 0, 5, 10, 15, 20, 25, 30, 35, 40, 45, 50, 55, 60, 65, 70, 75, 80, 85, 90, 95, 100}, "Centrality bins"}; Configurable> ptBins{"ptBins", {0.20, 0.25, 0.30, 0.35, 0.40, 0.45, 0.50, 0.55, 0.60, 0.65, 0.70, 0.75, 0.80, 0.85, 0.90, 0.95, 1.00, 1.05, 1.10, 1.15, 1.20, 1.25, 1.30, 1.35, 1.40, 1.45, 1.50, 1.55, 1.60, 1.65, 1.70, 1.75, 1.80, 1.85, 1.90, 1.95, 2.00}, "p_{T} bins"}; - Configurable> etaBins{"etaBins", {-0.8, -0.75, -0.7, -0.65, -0.6, -0.55, -0.5, -0.45, -0.4, -0.35, -0.3, -0.25, -0.2, -0.15, -0.1, -0.05, 0.0, 0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5, 0.55, 0.6, 0.65, 0.7, 0.75, 0.8}, "#eta bins"}; - - Configurable> ptBinsPi{"ptBinsPi", {0.20, 0.225, 0.25, 0.275, 0.30, 0.325, 0.35, 0.375, 0.40, 0.425, 0.45, 0.475, 0.50, 0.525, 0.55, 0.575, 0.60, 0.70, 0.90, 1.1, 1.3, 1.4, 1.45, 1.5, 1.55, 1.6, 1.65, 1.70, 1.75, 1.80, 1.85, 1.90, 1.95, 2.00}, "p_{T} bins for pions"}; - Configurable> ptBinsKa{"ptBinsKa", {0.30, 0.325, 0.35, 0.375, 0.40, 0.425, 0.45, 0.475, 0.50, 0.55, 0.60, 0.65, 0.70, 0.75, 0.8, 0.85, 0.90, 0.95, 1.0, 1.05, 1.1, 1.15, 1.2, 1.25, 1.3, 1.35, 1.4, 1.45, 1.5, 1.55, 1.6, 1.70, 1.80, 1.90, 2.00}, "p_{T} bins for kaons"}; - Configurable> ptBinsPr{"ptBinsPr", {0.40, 0.425, 0.45, 0.475, 0.5, 0.525, 0.55, 0.575, 0.60, 0.625, 0.65, 0.675, 0.70, 0.725, 0.75, 0.775, 0.80, 0.825, 0.85, 0.875, 0.90, 0.925, 0.95, 1.0, 1.025, 1.05, 1.075, 1.1, 1.15, 1.2, 1.4, 1.6, 1.8, 2.0}, "p_{T} bins for protons"}; + Configurable> etaBins{"etaBins", {-0.8, -0.7, -0.6, -0.5, -0.4, -0.3, -0.2, -0.1, 0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8}, "#eta bins"}; Service pdg{}; + Service ccdb; HistogramRegistry hist{"hist", {}, OutputObjHandlingPolicy::AnalysisObject}; + TH3* efficiencyWeightCh = nullptr; + TH3* efficiencyWeightPi = nullptr; + TH3* efficiencyWeightKa = nullptr; + TH3* efficiencyWeightPr = nullptr; + enum CollisionLabels { kTotCol = 1, kPassSelCol @@ -154,7 +170,7 @@ struct MeanPtFlucId { Gen_Proton }; - static constexpr std::array Dire{ + static constexpr std::array Directories{ "QA/after/", "QA/Pion/", "QA/Kaon/", @@ -167,10 +183,34 @@ struct MeanPtFlucId { "Gen/Pion/", "Gen/Kaon/", "Gen/Proton/"}; + static constexpr auto const& dire = Directories; + + enum CorrectedMode { + Corrected_Charged = 0, + Corrected_Pion, + Corrected_Kaon, + Corrected_Proton + }; + + static constexpr std::array CorrectedDirectories{ + "Corrected/Charged/", + "Corrected/Pion/", + "Corrected/Kaon/", + "Corrected/Proton/"}; + + static constexpr std::array TruthMatchedCorrectedDirectories{ + "CorrectedTruthMatched/Charged/", + "CorrectedTruthMatched/Pion/", + "CorrectedTruthMatched/Kaon/", + "CorrectedTruthMatched/Proton/"}; struct TrackInfo { double pt; - double eta; + }; + + struct WeightedTrackInfo { + double pt; + double weight; }; std::vector chargedTracks; @@ -183,19 +223,85 @@ struct MeanPtFlucId { std::vector kaonTracksGen; std::vector protonTracksGen; + std::vector chargedTracksCorrected; + std::vector pionTracksCorrected; + std::vector kaonTracksCorrected; + std::vector protonTracksCorrected; + std::vector chargedTracksTruthMatchedCorrected; + std::vector pionTracksTruthMatchedCorrected; + std::vector kaonTracksTruthMatchedCorrected; + std::vector protonTracksTruthMatchedCorrected; + + // Reused event-wise storage for the factorial-moment matrices. + std::vector matrixNInBin; + std::vector matrixOccupiedBins; + std::vector matrixWInBin; + std::vector matrixW2InBin; + void init(InitContext const&) { - const int nEta = etaBins->size() - 1; - const int nPtCh = ptBins->size() - 1; - const int nPtPi = ptBinsPi->size() - 1; - const int nPtKa = ptBinsKa->size() - 1; - const int nPtPr = ptBinsPr->size() - 1; + if (nPtBins <= 0 || cfgCutPtMin >= cfgCutPtMax || cfgCutPiPtMin >= cfgCutPtMax || + cfgCutKaPtMin >= cfgCutPtMax || cfgCutPrPtMin >= cfgCutPtMax) { + LOGF(fatal, "Invalid pT matrix binning: nPtBins=%d, max=%g, minima=(%g,%g,%g,%g)", + static_cast(nPtBins), static_cast(cfgCutPtMax), + static_cast(cfgCutPtMin), static_cast(cfgCutPiPtMin), + static_cast(cfgCutKaPtMin), static_cast(cfgCutPrPtMin)); + } - const AxisSpec axisKCh{nEta * nPtCh, 0, static_cast(nEta * nPtCh), "k_{ch}"}; - const AxisSpec axisKPi{nEta * nPtPi, 0, static_cast(nEta * nPtPi), "k_{#pi}"}; - const AxisSpec axisKKa{nEta * nPtKa, 0, static_cast(nEta * nPtKa), "k_{K}"}; - const AxisSpec axisKPr{nEta * nPtPr, 0, static_cast(nEta * nPtPr), "k_{p}"}; + if (cfgApplyEfficiencyCorrection) { + ccdb->setURL(cfgEfficiencyCCDBUrl.value); + ccdb->setCaching(true); + TList* efficiencyList = ccdb->getForTimeStamp(cfgEfficiencyCCDBPath.value, cfgEfficiencyTimestamp.value); + if (!efficiencyList) { + LOGF(fatal, "Cannot load ccdb_object containing the efficiency maps"); + } + + const char* chargedMapName = cfgUseCombinedPurityEfficiencyWeight ? "hWeightPtEtaCent" : "hInverseEfficiencyPtEtaCent"; + const char* pionMapName = cfgUseCombinedPurityEfficiencyWeight ? "hWeightPtEtaCentPi" : "hInverseEfficiencyPtEtaCentPi"; + const char* kaonMapName = cfgUseCombinedPurityEfficiencyWeight ? "hWeightPtEtaCentKa" : "hInverseEfficiencyPtEtaCentKa"; + const char* protonMapName = cfgUseCombinedPurityEfficiencyWeight ? "hWeightPtEtaCentPr" : "hInverseEfficiencyPtEtaCentPr"; + efficiencyWeightCh = dynamic_cast(efficiencyList->FindObject(chargedMapName)); + efficiencyWeightPi = dynamic_cast(efficiencyList->FindObject(pionMapName)); + efficiencyWeightKa = dynamic_cast(efficiencyList->FindObject(kaonMapName)); + efficiencyWeightPr = dynamic_cast(efficiencyList->FindObject(protonMapName)); + if (!efficiencyWeightCh || !efficiencyWeightPi || !efficiencyWeightKa || !efficiencyWeightPr) { + LOGF(fatal, "Missing efficiency map(s): charged=%s, pion=%s, kaon=%s, proton=%s", + chargedMapName, pionMapName, kaonMapName, protonMapName); + } + + auto validateEfficiencyMap = [&](const TH3* map, const char* name, double minPt) { + constexpr double AxisTolerance = 1.e-6; + constexpr double MaximumCentrality = 100.; + const bool coversEta = map->GetXaxis()->GetXmin() <= -cfgCutEta + AxisTolerance && + map->GetXaxis()->GetXmax() >= cfgCutEta - AxisTolerance; + const bool coversPt = map->GetYaxis()->GetXmin() <= minPt + AxisTolerance && + map->GetYaxis()->GetXmax() >= cfgCutPtMax - AxisTolerance; + const bool coversCentrality = map->GetZaxis()->GetXmin() <= AxisTolerance && + map->GetZaxis()->GetXmax() >= MaximumCentrality - AxisTolerance; + if (!coversEta || !coversPt || !coversCentrality) { + LOGF(fatal, + "Efficiency map %s does not cover the analysis phase space: " + "eta=[%g,%g], pT=[%g,%g], centrality=[%g,%g]. " + "Expected eta=[%g,%g], pT=[%g,%g], centrality=[0,100].", + name, + map->GetXaxis()->GetXmin(), map->GetXaxis()->GetXmax(), + map->GetYaxis()->GetXmin(), map->GetYaxis()->GetXmax(), + map->GetZaxis()->GetXmin(), map->GetZaxis()->GetXmax(), + -static_cast(cfgCutEta), static_cast(cfgCutEta), + minPt, static_cast(cfgCutPtMax)); + } + }; + validateEfficiencyMap(efficiencyWeightCh, chargedMapName, cfgCutPtMin); + validateEfficiencyMap(efficiencyWeightPi, pionMapName, cfgCutPiPtMin); + validateEfficiencyMap(efficiencyWeightKa, kaonMapName, cfgCutKaPtMin); + validateEfficiencyMap(efficiencyWeightPr, protonMapName, cfgCutPrPtMin); + } + + const AxisSpec axisPtCh{nPtBins, cfgCutPtMin, cfgCutPtMax, "p_{T} (GeV/c)"}; + const AxisSpec axisPtPi{nPtBins, cfgCutPiPtMin, cfgCutPtMax, "p_{T} (GeV/c)"}; + const AxisSpec axisPtKa{nPtBins, cfgCutKaPtMin, cfgCutPtMax, "p_{T} (GeV/c)"}; + const AxisSpec axisPtPr{nPtBins, cfgCutPrPtMin, cfgCutPtMax, "p_{T} (GeV/c)"}; const AxisSpec axisCol{3, 1, 4, ""}; const AxisSpec axisTrack{5, 1, 6, ""}; @@ -204,13 +310,8 @@ struct MeanPtFlucId { const AxisSpec axisPhi{nPhiBins, 0., +7., "#phi (rad)"}; const AxisSpec axisY{100, -0.6, 0.6, "y"}; const AxisSpec axisPt{ptBins, "p_{T} (GeV/c)"}; - const AxisSpec axisPtPi{ptBinsPi, "p_{T} (GeV/c)"}; - const AxisSpec axisPtKa{ptBinsKa, "p_{T}(GeV/c)"}; - const AxisSpec axisPtPr{ptBinsPr, "p_{T} (GeV/c)"}; const AxisSpec axisP{nPBins, 0., 3., "p (GeV/c)"}; const AxisSpec axisInnerParam{nPBins, 0., 3., "p_{InnerParam} (GeV/c)"}; - const AxisSpec axisPart{nPartBins, 0., 18., " "}; - const AxisSpec axisMult{100, 0, 100, "N_{ch}"}; const AxisSpec axisMultTPC{multTPCBins, "N_{TPC} "}; const AxisSpec axisMultCorr{multCorrBins, "N_{Corr} "}; const AxisSpec axisMultMC{multMCBins, "N_{Gen} "}; @@ -229,8 +330,6 @@ struct MeanPtFlucId { const AxisSpec axisM2{100, 0., 1.4, "#it{m}^{2} (GeV/#it{c}^{2})^{2}"}; const AxisSpec axisPid{300, 0, 3000, "PID"}; - HistogramConfigSpec partHistCent({HistType::kTHnSparseD, {axisCentFT0M, axisPart}}); - HistogramConfigSpec partMCHist({HistType::kTHnSparseD, {axisCentFT0M, axisPart}}); HistogramConfigSpec tofNSigmaHist({HistType::kTH2D, {axisP, axisTOFNsigma}}); HistogramConfigSpec tofSignalHist({HistType::kTH2D, {axisP, axisTOFSignal}}); HistogramConfigSpec tpcNSigmaHist({HistType::kTH2D, {axisP, axisTPCNsigma}}); @@ -238,219 +337,226 @@ struct MeanPtFlucId { HistogramConfigSpec tpcTofHist({HistType::kTH2D, {axisTPCNsigma, axisTOFNsigma}}); HistogramConfigSpec pvsM2Hist({HistType::kTH2D, {axisM2, axisP}}); HistogramConfigSpec tpcSignalHist1({HistType::kTH2D, {axisInnerParam, axisTPCSignal}}); - HistogramConfigSpec pvsM2Hist1({HistType::kTH2D, {axisM2, axisInnerParam}}); - - // QA Plots - hist.add("QA/before/h_Counts", "Counts", kTH1D, {axisEvents}); - hist.add("QA/before/h_VtxZ", "V_{Z}", kTH1D, {axisVtxZ}); - hist.add("QA/before/h_NTPC", "N_{TPC}", kTH1D, {axisMultTPC}); - hist.add("QA/before/h_NFT0M", "FT0M Multiplicity", kTH1D, {axisMultFT0M}); - hist.add("QA/before/h_CentM", "FT0M (%)", kTH1D, {axisCentFT0M}); - - hist.add("QA/before/h2_TPCSignal", "TPC Signal", tpcSignalHist); - hist.add("QA/before/h2_TOFSignal", "TOF Signal", tofSignalHist); - hist.add("QA/before/h2_pvsm2", "p vs m^{2}", pvsM2Hist); - - hist.addClone("QA/before/", "QA/after/"); - - hist.add("QA/before/h_Pt", "p_{T}", kTH1D, {axisPt}); - hist.add("QA/before/h_Eta", "#eta ", kTH1D, {axisEta}); - hist.add("QA/before/h_Phi", "#phi ", kTH1D, {axisPhi}); - hist.add("QA/before/h_DcaZ", "DCA_{Z}", kTH1D, {axisDCAz}); - hist.add("QA/before/h_DcaXY", "DCA_{XY}", kTH1D, {axisDCAxy}); - hist.add("QA/before/h2_DcaZ", "DCA_{Z}", kTH2D, {{axisPt}, {axisDCAz}}); - hist.add("QA/before/h2_DcaXY", "DCA_{XY}", kTH2D, {{axisPt}, {axisDCAxy}}); - - hist.add("QA/after/h_Ncorr", "N_{corr}", kTH1D, {axisMultCorr}); - hist.add("QA/after/h_counts_evSelCuts", "Event selection cuts", kTH1D, {axisEvents}); - hist.add("QA/after/h_TPCChi2perCluster", "TPC #Chi^{2}/Cluster", kTH1D, {axisChi2}); - hist.add("QA/after/h_ITSChi2perCluster", "ITS #Chi^{2}/Cluster", kTH1D, {axisChi2}); - hist.add("QA/after/h_crossedTPC", "Crossed TPC", kTH1D, {axisCrossedTPC}); - hist.add("QA/after/h2_NTPC_CentM", "N_{TPC} vs FT0M(%)", kTH2D, {{axisCentFT0M}, {axisMultTPC}}); - hist.add("QA/after/h2_NTPC_NFT0M", "N_{TPC} vs N_{FT0M}", kTH2D, {{axisMultFT0M}, {axisMultTPC}}); - hist.add("QA/after/p_NTPC_NFT0M", "N_{TPC} vs N_{FT0M} (Profile)", kTProfile, {axisMultFT0M}); - hist.add("QA/after/p_NTPC_CentM", "N_{TPC} vs FT0M(%) (Profile)", kTProfile, {axisCentFT0M}); - hist.add("QA/after/h_DCAxy_primary", "DCA_{XY} (Primary)", kTH1D, {axisDCAxy}); - hist.add("QA/after/h_DCAz_primary", "DCA_{Z} (Primary)", kTH1D, {axisDCAz}); - hist.add("QA/after/h_DCAxy_secondary", "DCA_{XY} (Secondary)", kTH1D, {axisDCAxy}); - hist.add("QA/after/h_DCAz_secondary", "DCA_{Z} (Secondary)", kTH1D, {axisDCAz}); - hist.add("QA/after/innerParam/h2_TPCSignal", "TPC Signal", tpcSignalHist1); - hist.add("QA/after/h2_NSim_NTPC", "Reco vs Truth Multiplicty TPC", kTH2D, {{axisMultMC}, {axisMultTPC}}); - hist.add("QA/after/h2_NTPC_NSim", "Truth vs Reco Multiplicty TPC", kTH2D, {{axisMultTPC}, {axisMultMC}}); - hist.add("QA/after/p_NTPC_NSim", "Truth vs Reco Multiplicty TPC", kTProfile, {axisMultTPC}); - hist.add("QA/after/h2_Ncorr_NSim", "Truth vs Corrected Multiplicty TPC", kTH2D, {{axisMultCorr}, {axisMultMC}}); - hist.add("QA/after/h2_Ncorr_CentM", "N_{corr} vs FT0M(%)", kTH2D, {{axisCentFT0M}, {axisMultCorr}}); - - hist.add("QA/Charged/h_Pt", "p_{T}", kTH1D, {axisPt}); - hist.add("QA/Charged/h_Eta", "#eta ", kTH1D, {axisEta}); - hist.add("QA/Charged/h_Phi", "#phi ", kTH1D, {axisPhi}); - hist.add("QA/Charged/h_DcaZ", "DCA_{Z}", kTH1D, {axisDCAz}); - hist.add("QA/Charged/h_DcaXY", "DCA_{XY}", kTH1D, {axisDCAxy}); - hist.add("QA/Charged/h2_Pt_PtMC", "p_{T} vs p_{T} (MC)", kTH2D, {{axisPt}, {axisPt}}); - hist.add("QA/Charged/h2_DcaZ", "DCA_{Z}", kTH2D, {{axisPt}, {axisDCAz}}); - hist.add("QA/Charged/h2_DcaXY", "DCA_{XY}", kTH2D, {{axisPt}, {axisDCAxy}}); - hist.add("QA/Charged/h2_Pt_centFT0M", "p_{T} in centrality Classes ", kTH2D, {{axisCentFT0M}, {axisPt}}); - hist.add("QA/Charged/h3_Pt_Eta_Phi", "p_{T}, #eta, #phi ", kTHnSparseD, {{axisPt}, {axisEta}, {axisPhi}}); - - hist.addClone("QA/Charged/", "QA/Pion/"); - - hist.add("QA/Pion/before/h2_TPCNsigma", "n #sigma_{TPC}", tpcNSigmaHist); - hist.add("QA/Pion/before/h2_TPCNsigma_nottof", "n #sigma_{TPC}", tpcNSigmaHist); - hist.add("QA/Pion/before/h2_TPCNsigma_tof", "n #sigma_{TPC}", tpcNSigmaHist); - hist.add("QA/Pion/before/h2_TOFNsigma", "n #sigma_{TOF}", tofNSigmaHist); - hist.add("QA/Pion/before/h2_TpcTofNsigma", "n #sigma_{TPC} vs n #sigma_{TOF}", tpcTofHist); - hist.add("QA/Pion/h_Rap", "y ", kTH1D, {axisY}); - hist.add("QA/Pion/h2_TPCNsigma", "n #sigma_{TPC}", tpcNSigmaHist); - hist.add("QA/Pion/h2_TOFNsigma", "n #sigma_{TOF}", tofNSigmaHist); - hist.add("QA/Pion/h2_TpcTofNsigma", "n #sigma_{TPC} vs n #sigma_{TOF}", tpcTofHist); - hist.add("QA/Pion/h2_TPCSignal", "TPC Signal ", tpcSignalHist); - hist.add("QA/Pion/h2_TOFSignal", "TOF Signal", tofSignalHist); - hist.add("QA/Pion/innerParam/h2_TPCSignal", "TPC Signal", tpcSignalHist1); - hist.add("QA/Pion/h2_pvsm2", "p vs m^{2}", pvsM2Hist); - hist.addClone("QA/Pion/", "QA/Kaon/"); - hist.addClone("QA/Pion/", "QA/Proton/"); - - hist.add("QA/Charged/h_PtMC", "p_{T} (MC)", kTH1D, {axisPt}); - hist.add("QA/Charged/h2_Pt_Eta", "p_{T} vs #eta ", kTH2D, {{axisEta}, {axisPt}}); - hist.add("QA/Charged/h3_Pt_Eta_centFT0M", "p_{T} vs #eta in centrality classes", kTH3D, {{axisEta}, {axisPt}, {axisCent}}); - hist.add("QA/Charged/h3_Pt_EtaMC_centFT0M", "p_{T} vs #eta in centrality classes", kTH3D, {{axisEta}, {axisPt}, {axisCent}}); - hist.add("QA/Charged/h2_Pt_EtaMC", "p_{T} vs #eta (MC)", kTH2D, {{axisEta}, {axisPt}}); - - hist.add("QA/Pion/h_Pt_e", "p_{T}", kTH1D, {axisPtPi}); - hist.add("QA/Pion/h_PtMC", "p_{T} (MC)", kTH1D, {axisPtPi}); - hist.add("QA/Pion/h2_Pt_Eta", "p_{T} vs #eta ", kTH2D, {{axisEta}, {axisPtPi}}); - hist.add("QA/Pion/h3_Pt_Eta_centFT0M", "p_{T} vs #eta in centrality classes", kTH3D, {{axisEta}, {axisPtPi}, {axisCent}}); - hist.add("QA/Pion/h3_Pt_EtaMC_centFT0M", "p_{T} vs #eta in centrality classes", kTH3D, {{axisEta}, {axisPtPi}, {axisCent}}); - hist.add("QA/Pion/h2_Pt_EtaMC", "p_{T} vs #eta (MC)", kTH2D, {{axisEta}, {axisPtPi}}); - hist.add("QA/Pion/h_PtPos", "p_{T} (positive) ", kTH1D, {axisPtPi}); - hist.add("QA/Pion/h_PtNeg", "p_{T} (negative) ", kTH1D, {axisPtPi}); - hist.add("QA/Pion/h_PtTruth", "p_{T} (Truth)", kTH1D, {axisPtPi}); - hist.add("QA/Pion/h_PtPosTruth", "p_{T} (positive) (Truth)", kTH1D, {axisPtPi}); - hist.add("QA/Pion/h_PtNegTruth", "p_{T} (negative) (Truth) ", kTH1D, {axisPtPi}); - hist.add("QA/Pion/h_PtTruth_primary", "p_{T} (Truth Primary)", kTH1D, {axisPtPi}); - hist.add("QA/Pion/h_PtTruth_secondary", "p_{T} (Truth Secondary)", kTH1D, {axisPtPi}); - hist.add("QA/Pion/h2_Pt_EtaTruth", "p_{T} vs #eta", kTH2D, {{axisEta}, {axisPtPi}}); - - hist.add("QA/Kaon/h_Pt_e", "p_{T}", kTH1D, {axisPtKa}); - hist.add("QA/Kaon/h_PtMC", "p_{T} (MC)", kTH1D, {axisPtKa}); - hist.add("QA/Kaon/h2_Pt_Eta", "p_{T} vs #eta ", kTH2D, {{axisEta}, {axisPtKa}}); - hist.add("QA/Kaon/h3_Pt_Eta_centFT0M", "p_{T} vs #eta in centrality classes", kTH3D, {{axisEta}, {axisPtKa}, {axisCent}}); - hist.add("QA/Kaon/h3_Pt_EtaMC_centFT0M", "p_{T} vs #eta in centrality classes", kTH3D, {{axisEta}, {axisPtKa}, {axisCent}}); - hist.add("QA/Kaon/h2_Pt_EtaMC", "p_{T} vs #eta (MC)", kTH2D, {{axisEta}, {axisPtKa}}); - hist.add("QA/Kaon/h_PtPos", "p_{T} (positive) ", kTH1D, {axisPtKa}); - hist.add("QA/Kaon/h_PtNeg", "p_{T} (negative) ", kTH1D, {axisPtKa}); - hist.add("QA/Kaon/h_PtTruth", "p_{T} (Truth)", kTH1D, {axisPtKa}); - hist.add("QA/Kaon/h_PtPosTruth", "p_{T} (positive) (Truth)", kTH1D, {axisPtKa}); - hist.add("QA/Kaon/h_PtNegTruth", "p_{T} (negative) (Truth) ", kTH1D, {axisPtKa}); - hist.add("QA/Kaon/h_PtTruth_primary", "p_{T} (Truth Primary)", kTH1D, {axisPtKa}); - hist.add("QA/Kaon/h_PtTruth_secondary", "p_{T} (Truth Secondary)", kTH1D, {axisPtKa}); - hist.add("QA/Kaon/h2_Pt_EtaTruth", "p_{T} vs #eta", kTH2D, {{axisEta}, {axisPtKa}}); - - hist.add("QA/Proton/h_Pt_e", "p_{T}", kTH1D, {axisPtPr}); - hist.add("QA/Proton/h_PtMC", "p_{T} (MC)", kTH1D, {axisPtPr}); - hist.add("QA/Proton/h2_Pt_Eta", "p_{T} vs #eta ", kTH2D, {{axisEta}, {axisPtPr}}); - hist.add("QA/Proton/h3_Pt_Eta_centFT0M", "p_{T} vs #eta in centrality classes", kTH3D, {{axisEta}, {axisPtPr}, {axisCent}}); - hist.add("QA/Proton/h3_Pt_EtaMC_centFT0M", "p_{T} vs #eta in centrality classes", kTH3D, {{axisEta}, {axisPtPr}, {axisCent}}); - hist.add("QA/Proton/h2_Pt_EtaMC", "p_{T} vs #eta (MC)", kTH2D, {{axisEta}, {axisPtPr}}); - hist.add("QA/Proton/h_PtPos", "p_{T} (positive) ", kTH1D, {axisPtPr}); - hist.add("QA/Proton/h_PtNeg", "p_{T} (negative) ", kTH1D, {axisPtPr}); - hist.add("QA/Proton/h_PtTruth", "p_{T} (Truth)", kTH1D, {axisPtPr}); - hist.add("QA/Proton/h_PtPosTruth", "p_{T} (positive) (Truth)", kTH1D, {axisPtPr}); - hist.add("QA/Proton/h_PtNegTruth", "p_{T} (negative) (Truth) ", kTH1D, {axisPtPr}); - hist.add("QA/Proton/h_PtTruth_primary", "p_{T} (Truth Primary)", kTH1D, {axisPtPr}); - hist.add("QA/Proton/h_PtTruth_secondary", "p_{T} (Truth Secondary)", kTH1D, {axisPtPr}); - hist.add("QA/Proton/h2_Pt_EtaTruth", "p_{T} vs #eta", kTH2D, {{axisEta}, {axisPtPr}}); - - // AnalysisPlots - hist.add("Analysis/Charged/h_mult", "Multiplicity", kTH1D, {axisMult}); - hist.add("Analysis/Charged/p_Q1_Mult", " Q1 vs Mult ", kTProfile, {axisMultTPC}); - hist.add("Analysis/Charged/p_Q2_Mult", " Q2 vs Mult ", kTProfile, {axisMultTPC}); - hist.add("Analysis/Charged/p_mean_pT_Mult", " vs Mult ", kTProfile, {axisMultTPC}); - hist.add("Analysis/Charged/p_mean_pT_Cent", " vs CentFT0M ", kTProfile, {axisCentFT0M}); - hist.add("Analysis/Charged/p_mean_pT_Cent_var", " vs CentFT0M ", kTProfile, {axisCentFT0M}); - hist.add("Analysis/Charged/p_twopart_Cent_var", "Twopart vs CentFT0M ", kTProfile, {axisCentFT0M}); - hist.add("Analysis/Charged/h_mean_pT_Cent", " vs CentFT0M ", kTH1D, {axisCentFT0M}); - hist.add("Analysis/Charged/h_mean_pT_Cent_var", " vs CentFT0M ", kTH1D, {axisCentFT0M}); - hist.add("Analysis/Charged/h_twopart_Cent_var", "Twopart vs CentFT0M ", kTH1D, {axisCentFT0M}); - hist.add("Analysis/Charged/p_mean_pT_Num_Cent_var", " Numerator vs CentFT0M ", kTProfile, {axisCentFT0M}); - hist.add("Analysis/Charged/p_mean_pT_Den_Cent_var", " Denominator vs CentFT0M ", kTProfile, {axisCentFT0M}); - hist.add("Analysis/Charged/p_twopart_Num_Cent_var", "Twopart Numerator vs CentFT0M ", kTProfile, {axisCentFT0M}); - hist.add("Analysis/Charged/p_twopart_Den_Cent_var", "Twopart Denominator vs CentFT0M ", kTProfile, {axisCentFT0M}); - hist.add("Analysis/Charged/h_mean_pT_Num_Cent_var", " Numerator vs CentFT0M ", kTH1D, {axisCentFT0M}); - hist.add("Analysis/Charged/h_mean_pT_Den_Cent_var", " Denominator vs CentFT0M ", kTH1D, {axisCentFT0M}); - hist.add("Analysis/Charged/h_twopart_Num_Cent_var", "Twopart Numerator vs CentFT0M ", kTH1D, {axisCentFT0M}); - hist.add("Analysis/Charged/h_twopart_Den_Cent_var", "Twopart Denominator vs CentFT0M ", kTH1D, {axisCentFT0M}); - - hist.addClone("Analysis/Charged/", "Analysis/Pion/"); - hist.addClone("Analysis/Charged/", "Analysis/Kaon/"); - hist.addClone("Analysis/Charged/", "Analysis/Proton/"); - - hist.add("Analysis/Charged/hN1Matrix", "", kTH2D, {{axisCent}, {axisKCh}}); - hist.add("Analysis/Charged/hPt1Matrix", "", kTH2D, {{axisCent}, {axisKCh}}); - hist.add("Analysis/Charged/hN2Matrix", "", kTH3D, {{axisCent}, {axisKCh}, {axisKCh}}); - hist.add("Analysis/Charged/hPtPtMatrix", "", kTH3D, {{axisCent}, {axisKCh}, {axisKCh}}); - - hist.add("Analysis/Pion/hN1Matrix", "", kTH2D, {{axisCent}, {axisKPi}}); - hist.add("Analysis/Pion/hPt1Matrix", "", kTH2D, {{axisCent}, {axisKPi}}); - hist.add("Analysis/Pion/hN2Matrix", "", kTH3D, {{axisCent}, {axisKPi}, {axisKPi}}); - hist.add("Analysis/Pion/hPtPtMatrix", "", kTH3D, {{axisCent}, {axisKPi}, {axisKPi}}); - - hist.add("Analysis/Kaon/hN1Matrix", "", kTH2D, {{axisCent}, {axisKKa}}); - hist.add("Analysis/Kaon/hPt1Matrix", "", kTH2D, {{axisCent}, {axisKKa}}); - hist.add("Analysis/Kaon/hN2Matrix", "", kTH3D, {{axisCent}, {axisKKa}, {axisKKa}}); - hist.add("Analysis/Kaon/hPtPtMatrix", "", kTH3D, {{axisCent}, {axisKKa}, {axisKKa}}); - - hist.add("Analysis/Proton/hN1Matrix", "", kTH2D, {{axisCent}, {axisKPr}}); - hist.add("Analysis/Proton/hPt1Matrix", "", kTH2D, {{axisCent}, {axisKPr}}); - hist.add("Analysis/Proton/hN2Matrix", "", kTH3D, {{axisCent}, {axisKPr}, {axisKPr}}); - hist.add("Analysis/Proton/hPtPtMatrix", "", kTH3D, {{axisCent}, {axisKPr}, {axisKPr}}); - - // MC Generated - hist.add("Gen/h_Counts", "Counts", kTH1D, {axisEvents}); - hist.add("Gen/h_VtxZ", "Vertex Z ", kTH1D, {axisVtxZ}); - hist.add("Gen/h_VtxZ_b", "Vertex Z ", kTH1D, {axisVtxZ}); - hist.add("Gen/h_NSim", "Truth Multiplicity TPC", kTH1D, {axisMultMC}); - hist.add("Gen/h2_NSim_NTPC", "Reco vs Truth Multiplicty TPC", kTH2D, {{axisMultMC}, {axisMultTPC}}); - hist.add("Gen/h2_NTPC_NSim", "Truth vs Reco Multiplicty TPC", kTH2D, {{axisMultTPC}, {axisMultMC}}); - - hist.add("Gen/Charged/h_EtaTruth", "#eta ", kTH1D, {axisEta}); - hist.add("Gen/Charged/h_PhiTruth", "#phi ", kTH1D, {axisPhi}); - hist.add("Gen/Charged/h2_PtTruth_centFT0M", "p_{T} in centrality Classes ", kTH2D, {{axisCentFT0M}, {axisPt}}); - hist.add("Gen/Charged/h3_Pt_Eta_PhiTruth", "p_{T}, #eta, #phi ", kTHnSparseD, {{axisPt}, {axisEta}, {axisPhi}}); - - hist.addClone("Gen/Charged/", "Gen/Pion/"); - - hist.add("Gen/Pion/h_RapTruth", "y", kTH1D, {axisY}); - - hist.addClone("Gen/Pion/", "Gen/Kaon/"); - hist.addClone("Gen/Pion/", "Gen/Proton/"); - - hist.addClone("Analysis/Charged/", "Gen/Charged/"); - hist.addClone("Analysis/Pion/", "Gen/Pion/"); - hist.addClone("Analysis/Kaon/", "Gen/Kaon/"); - hist.addClone("Analysis/Proton/", "Gen/Proton/"); - - hist.add("Gen/Charged/h_PtTruth", "p_{T} ", kTH1D, {axisPt}); - hist.add("Gen/Charged/h2_Pt_EtaTruth", "p_{T} vs #eta", kTH2D, {{axisEta}, {axisPt}}); - hist.add("Gen/Charged/h3_Pt_EtaTruth_centFT0M", "p_{T} vs #eta in centrality classes", kTH3D, {{axisEta}, {axisPt}, {axisCent}}); - - hist.add("Gen/Pion/h_PtTruth", "p_{T} ", kTH1D, {axisPtPi}); - hist.add("Gen/Pion/h_PtPosTruth", "p_{T} (positive) ", kTH1D, {axisPtPi}); - hist.add("Gen/Pion/h_PtNegTruth", "p_{T} (negative) ", kTH1D, {axisPtPi}); - hist.add("Gen/Pion/h2_Pt_EtaTruth", "p_{T} vs #eta", kTH2D, {{axisEta}, {axisPtPi}}); - hist.add("Gen/Pion/h3_Pt_EtaTruth_centFT0M", "p_{T} vs #eta in centrality classes", kTH3D, {{axisEta}, {axisPtPi}, {axisCent}}); - - hist.add("Gen/Kaon/h_PtTruth", "p_{T} ", kTH1D, {axisPtKa}); - hist.add("Gen/Kaon/h_PtPosTruth", "p_{T} (positive) ", kTH1D, {axisPtKa}); - hist.add("Gen/Kaon/h_PtNegTruth", "p_{T} (negative) ", kTH1D, {axisPtKa}); - hist.add("Gen/Kaon/h2_Pt_EtaTruth", "p_{T} vs #eta", kTH2D, {{axisEta}, {axisPtKa}}); - hist.add("Gen/Kaon/h3_Pt_EtaTruth_centFT0M", "p_{T} vs #eta in centrality classes", kTH3D, {{axisEta}, {axisPtKa}, {axisCent}}); - - hist.add("Gen/Proton/h_PtTruth", "p_{T} ", kTH1D, {axisPtPr}); - hist.add("Gen/Proton/h_PtPosTruth", "p_{T} (positive) ", kTH1D, {axisPtPr}); - hist.add("Gen/Proton/h_PtNegTruth", "p_{T} (negative) ", kTH1D, {axisPtPr}); - hist.add("Gen/Proton/h2_Pt_EtaTruth", "p_{T} vs #eta", kTH2D, {{axisEta}, {axisPtPr}}); - hist.add("Gen/Proton/h3_Pt_EtaTruth_centFT0M", "p_{T} vs #eta in centrality classes", kTH3D, {{axisEta}, {axisPtPr}, {axisCent}}); + + auto addEventQAHistograms = [&]() { + // Event QA plots + hist.add("QA/before/h_Counts", "Counts", kTH1D, {axisEvents}); + hist.add("QA/before/h_VtxZ", "V_{Z}", kTH1D, {axisVtxZ}); + hist.add("QA/before/h_NTPC", "N_{TPC}", kTH1D, {axisMultTPC}); + hist.add("QA/before/h_NFT0M", "FT0M Multiplicity", kTH1D, {axisMultFT0M}); + hist.add("QA/before/h_CentM", "FT0M (%)", kTH1D, {axisCentFT0M}); + + hist.add("QA/before/h2_TPCSignal", "TPC Signal", tpcSignalHist); + hist.add("QA/before/h2_TOFSignal", "TOF Signal", tofSignalHist); + hist.add("QA/before/h2_pvsm2", "p vs m^{2}", pvsM2Hist); + + hist.addClone("QA/before/", "QA/after/"); + + hist.add("QA/before/h_Pt", "p_{T}", kTH1D, {axisPt}); + hist.add("QA/before/h_Eta", "#eta ", kTH1D, {axisEta}); + hist.add("QA/before/h_Phi", "#phi ", kTH1D, {axisPhi}); + hist.add("QA/before/h_DcaZ", "DCA_{Z}", kTH1D, {axisDCAz}); + hist.add("QA/before/h_DcaXY", "DCA_{XY}", kTH1D, {axisDCAxy}); + hist.add("QA/before/h2_DcaZ", "DCA_{Z}", kTH2D, {{axisPt}, {axisDCAz}}); + hist.add("QA/before/h2_DcaXY", "DCA_{XY}", kTH2D, {{axisPt}, {axisDCAxy}}); + + hist.add("QA/after/h_Ncorr", "N_{corr}", kTH1D, {axisMultCorr}); + hist.add("QA/after/h_counts_evSelCuts", "Event selection cuts", kTH1D, {axisEvents}); + hist.add("QA/after/h_TPCChi2perCluster", "TPC #Chi^{2}/Cluster", kTH1D, {axisChi2}); + hist.add("QA/after/h_ITSChi2perCluster", "ITS #Chi^{2}/Cluster", kTH1D, {axisChi2}); + hist.add("QA/after/h_crossedTPC", "Crossed TPC", kTH1D, {axisCrossedTPC}); + hist.add("QA/after/h2_NTPC_CentM", "N_{TPC} vs FT0M(%)", kTH2D, {{axisCentFT0M}, {axisMultTPC}}); + hist.add("QA/after/h2_NTPC_NFT0M", "N_{TPC} vs N_{FT0M}", kTH2D, {{axisMultFT0M}, {axisMultTPC}}); + hist.add("QA/after/p_NTPC_NFT0M", "N_{TPC} vs N_{FT0M} (Profile)", kTProfile, {axisMultFT0M}); + hist.add("QA/after/p_NTPC_CentM", "N_{TPC} vs FT0M(%) (Profile)", kTProfile, {axisCentFT0M}); + hist.add("QA/after/h_DCAxy_primary", "DCA_{XY} (Primary)", kTH1D, {axisDCAxy}); + hist.add("QA/after/h_DCAz_primary", "DCA_{Z} (Primary)", kTH1D, {axisDCAz}); + hist.add("QA/after/h_DCAxy_secondary", "DCA_{XY} (Secondary)", kTH1D, {axisDCAxy}); + hist.add("QA/after/h_DCAz_secondary", "DCA_{Z} (Secondary)", kTH1D, {axisDCAz}); + hist.add("QA/after/innerParam/h2_TPCSignal", "TPC Signal", tpcSignalHist1); + hist.add("QA/after/h2_NSim_NTPC", "Reco vs Truth Multiplicty TPC", kTH2D, {{axisMultMC}, {axisMultTPC}}); + hist.add("QA/after/h2_NTPC_NSim", "Truth vs Reco Multiplicty TPC", kTH2D, {{axisMultTPC}, {axisMultMC}}); + hist.add("QA/after/p_NTPC_NSim", "Truth vs Reco Multiplicty TPC", kTProfile, {axisMultTPC}); + hist.add("QA/after/h2_Ncorr_NSim", "Truth vs Corrected Multiplicty TPC", kTH2D, {{axisMultCorr}, {axisMultMC}}); + hist.add("QA/after/h2_Ncorr_CentM", "N_{corr} vs FT0M(%)", kTH2D, {{axisCentFT0M}, {axisMultCorr}}); + }; + addEventQAHistograms(); + + auto addTrackQAHistograms = [&]() { + hist.add("QA/Charged/h_Pt", "p_{T}", kTH1D, {axisPt}); + hist.add("QA/Charged/h_Eta", "#eta ", kTH1D, {axisEta}); + hist.add("QA/Charged/h_Phi", "#phi ", kTH1D, {axisPhi}); + hist.add("QA/Charged/h_DcaZ", "DCA_{Z}", kTH1D, {axisDCAz}); + hist.add("QA/Charged/h_DcaXY", "DCA_{XY}", kTH1D, {axisDCAxy}); + hist.add("QA/Charged/h2_Pt_PtMC", "p_{T} vs p_{T} (MC)", kTH2D, {{axisPt}, {axisPt}}); + hist.add("QA/Charged/h2_DcaZ", "DCA_{Z}", kTH2D, {{axisPt}, {axisDCAz}}); + hist.add("QA/Charged/h2_DcaXY", "DCA_{XY}", kTH2D, {{axisPt}, {axisDCAxy}}); + hist.add("QA/Charged/h2_Pt_centFT0M", "p_{T} in centrality Classes ", kTH2D, {{axisCentFT0M}, {axisPt}}); + hist.add("QA/Charged/h3_Pt_Eta_Phi", "p_{T}, #eta, #phi ", kTHnSparseD, {{axisPt}, {axisEta}, {axisPhi}}); + + hist.addClone("QA/Charged/", "QA/Pion/"); + + hist.add("QA/Pion/before/h2_TPCNsigma", "n #sigma_{TPC}", tpcNSigmaHist); + hist.add("QA/Pion/before/h2_TPCNsigma_nottof", "n #sigma_{TPC}", tpcNSigmaHist); + hist.add("QA/Pion/before/h2_TPCNsigma_tof", "n #sigma_{TPC}", tpcNSigmaHist); + hist.add("QA/Pion/before/h2_TOFNsigma", "n #sigma_{TOF}", tofNSigmaHist); + hist.add("QA/Pion/before/h2_TpcTofNsigma", "n #sigma_{TPC} vs n #sigma_{TOF}", tpcTofHist); + hist.add("QA/Pion/h_Rap", "y ", kTH1D, {axisY}); + hist.add("QA/Pion/h2_TPCNsigma", "n #sigma_{TPC}", tpcNSigmaHist); + hist.add("QA/Pion/h2_TOFNsigma", "n #sigma_{TOF}", tofNSigmaHist); + hist.add("QA/Pion/h2_TpcTofNsigma", "n #sigma_{TPC} vs n #sigma_{TOF}", tpcTofHist); + hist.add("QA/Pion/h2_TPCSignal", "TPC Signal ", tpcSignalHist); + hist.add("QA/Pion/h2_TOFSignal", "TOF Signal", tofSignalHist); + hist.add("QA/Pion/innerParam/h2_TPCSignal", "TPC Signal", tpcSignalHist1); + hist.add("QA/Pion/h2_pvsm2", "p vs m^{2}", pvsM2Hist); + hist.addClone("QA/Pion/", "QA/Kaon/"); + hist.addClone("QA/Pion/", "QA/Proton/"); + + hist.add("QA/Charged/h_Pt_e", "p_{T}", kTH1D, {axisPtCh}); + hist.add("QA/Charged/h_PtMC", "p_{T} (MC)", kTH1D, {axisPtCh}); + hist.add("QA/Charged/h_PtTruth_primary", "p_{T} (Truth Primary)", kTH1D, {axisPtCh}); + hist.add("QA/Charged/h2_Pt_Eta", "p_{T} vs #eta ", kTH2D, {{axisEta}, {axisPtCh}}); + hist.add("QA/Charged/h3_Pt_Eta_centFT0M", "p_{T} vs #eta in centrality classes", kTH3D, {{axisEta}, {axisPtCh}, {axisCent}}); + hist.add("QA/Charged/h3_Pt_EtaMC_centFT0M", "p_{T} vs #eta in centrality classes", kTH3D, {{axisEta}, {axisPtCh}, {axisCent}}); + hist.add("QA/Charged/h3_Pt_EtaMC_primary_centFT0M", "Selected primary MC particles;#eta_{MC};p_{T,MC};FT0M (%)", kTH3D, {{axisEta}, {axisPtCh}, {axisCent}}); + hist.add("QA/Charged/h2_Pt_EtaMC", "p_{T} vs #eta (MC)", kTH2D, {{axisEta}, {axisPtCh}}); + }; + addTrackQAHistograms(); + + auto addEfficiencyHistograms = [&]() { + hist.add("QA/Pion/h_Pt_e", "p_{T}", kTH1D, {axisPtPi}); + hist.add("QA/Pion/h_PtMC", "p_{T} (MC)", kTH1D, {axisPtPi}); + hist.add("QA/Pion/h2_Pt_Eta", "p_{T} vs #eta ", kTH2D, {{axisEta}, {axisPtPi}}); + hist.add("QA/Pion/h3_Pt_Eta_centFT0M", "p_{T} vs #eta in centrality classes", kTH3D, {{axisEta}, {axisPtPi}, {axisCent}}); + hist.add("QA/Pion/h3_Pt_EtaMC_centFT0M", "p_{T} vs #eta in centrality classes", kTH3D, {{axisEta}, {axisPtPi}, {axisCent}}); + hist.add("QA/Pion/h3_Pt_EtaMC_primary_centFT0M", "Correctly selected primary pions;#eta_{MC};p_{T,MC};FT0M (%)", kTH3D, {{axisEta}, {axisPtPi}, {axisCent}}); + hist.add("QA/Pion/h2_Pt_EtaMC", "p_{T} vs #eta (MC)", kTH2D, {{axisEta}, {axisPtPi}}); + hist.add("QA/Pion/h_PtPos", "p_{T} (positive) ", kTH1D, {axisPtPi}); + hist.add("QA/Pion/h_PtNeg", "p_{T} (negative) ", kTH1D, {axisPtPi}); + hist.add("QA/Pion/h_PtTruth", "p_{T} (Truth)", kTH1D, {axisPtPi}); + hist.add("QA/Pion/h_PtPosTruth", "p_{T} (positive) (Truth)", kTH1D, {axisPtPi}); + hist.add("QA/Pion/h_PtNegTruth", "p_{T} (negative) (Truth) ", kTH1D, {axisPtPi}); + hist.add("QA/Pion/h_PtTruth_primary", "p_{T} (Truth Primary)", kTH1D, {axisPtPi}); + hist.add("QA/Pion/h_PtTruth_secondary", "p_{T} (Truth Secondary)", kTH1D, {axisPtPi}); + hist.add("QA/Pion/h2_Pt_EtaTruth", "p_{T} vs #eta", kTH2D, {{axisEta}, {axisPtPi}}); + hist.add("QA/Pion/h2_Pt_EtaTruth_primary", "p_{T} vs #eta (Truth Primary)", kTH2D, {{axisEta}, {axisPtPi}}); + + hist.add("QA/Kaon/h_Pt_e", "p_{T}", kTH1D, {axisPtKa}); + hist.add("QA/Kaon/h_PtMC", "p_{T} (MC)", kTH1D, {axisPtKa}); + hist.add("QA/Kaon/h2_Pt_Eta", "p_{T} vs #eta ", kTH2D, {{axisEta}, {axisPtKa}}); + hist.add("QA/Kaon/h3_Pt_Eta_centFT0M", "p_{T} vs #eta in centrality classes", kTH3D, {{axisEta}, {axisPtKa}, {axisCent}}); + hist.add("QA/Kaon/h3_Pt_EtaMC_centFT0M", "p_{T} vs #eta in centrality classes", kTH3D, {{axisEta}, {axisPtKa}, {axisCent}}); + hist.add("QA/Kaon/h3_Pt_EtaMC_primary_centFT0M", "Correctly selected primary kaons;#eta_{MC};p_{T,MC};FT0M (%)", kTH3D, {{axisEta}, {axisPtKa}, {axisCent}}); + hist.add("QA/Kaon/h2_Pt_EtaMC", "p_{T} vs #eta (MC)", kTH2D, {{axisEta}, {axisPtKa}}); + hist.add("QA/Kaon/h_PtPos", "p_{T} (positive) ", kTH1D, {axisPtKa}); + hist.add("QA/Kaon/h_PtNeg", "p_{T} (negative) ", kTH1D, {axisPtKa}); + hist.add("QA/Kaon/h_PtTruth", "p_{T} (Truth)", kTH1D, {axisPtKa}); + hist.add("QA/Kaon/h_PtPosTruth", "p_{T} (positive) (Truth)", kTH1D, {axisPtKa}); + hist.add("QA/Kaon/h_PtNegTruth", "p_{T} (negative) (Truth) ", kTH1D, {axisPtKa}); + hist.add("QA/Kaon/h_PtTruth_primary", "p_{T} (Truth Primary)", kTH1D, {axisPtKa}); + hist.add("QA/Kaon/h_PtTruth_secondary", "p_{T} (Truth Secondary)", kTH1D, {axisPtKa}); + hist.add("QA/Kaon/h2_Pt_EtaTruth", "p_{T} vs #eta", kTH2D, {{axisEta}, {axisPtKa}}); + hist.add("QA/Kaon/h2_Pt_EtaTruth_primary", "p_{T} vs #eta (Truth Primary)", kTH2D, {{axisEta}, {axisPtKa}}); + + hist.add("QA/Proton/h_Pt_e", "p_{T}", kTH1D, {axisPtPr}); + hist.add("QA/Proton/h_PtMC", "p_{T} (MC)", kTH1D, {axisPtPr}); + hist.add("QA/Proton/h2_Pt_Eta", "p_{T} vs #eta ", kTH2D, {{axisEta}, {axisPtPr}}); + hist.add("QA/Proton/h3_Pt_Eta_centFT0M", "p_{T} vs #eta in centrality classes", kTH3D, {{axisEta}, {axisPtPr}, {axisCent}}); + hist.add("QA/Proton/h3_Pt_EtaMC_centFT0M", "p_{T} vs #eta in centrality classes", kTH3D, {{axisEta}, {axisPtPr}, {axisCent}}); + hist.add("QA/Proton/h3_Pt_EtaMC_primary_centFT0M", "Correctly selected primary protons;#eta_{MC};p_{T,MC};FT0M (%)", kTH3D, {{axisEta}, {axisPtPr}, {axisCent}}); + hist.add("QA/Proton/h2_Pt_EtaMC", "p_{T} vs #eta (MC)", kTH2D, {{axisEta}, {axisPtPr}}); + hist.add("QA/Proton/h_PtPos", "p_{T} (positive) ", kTH1D, {axisPtPr}); + hist.add("QA/Proton/h_PtNeg", "p_{T} (negative) ", kTH1D, {axisPtPr}); + hist.add("QA/Proton/h_PtTruth", "p_{T} (Truth)", kTH1D, {axisPtPr}); + hist.add("QA/Proton/h_PtPosTruth", "p_{T} (positive) (Truth)", kTH1D, {axisPtPr}); + hist.add("QA/Proton/h_PtNegTruth", "p_{T} (negative) (Truth) ", kTH1D, {axisPtPr}); + hist.add("QA/Proton/h_PtTruth_primary", "p_{T} (Truth Primary)", kTH1D, {axisPtPr}); + hist.add("QA/Proton/h_PtTruth_secondary", "p_{T} (Truth Secondary)", kTH1D, {axisPtPr}); + hist.add("QA/Proton/h2_Pt_EtaTruth", "p_{T} vs #eta", kTH2D, {{axisEta}, {axisPtPr}}); + hist.add("QA/Proton/h2_Pt_EtaTruth_primary", "p_{T} vs #eta (Truth Primary)", kTH2D, {{axisEta}, {axisPtPr}}); + }; + addEfficiencyHistograms(); + + auto addMatrixHistograms = [&]() { + hist.add("Analysis/Charged/hN1Matrix", "Single-particle density", kTH2D, {{axisCent}, {axisPtCh}}); + hist.add("Analysis/Charged/hN2Matrix", "Ordered-pair density", kTH3D, {{axisCent}, {axisPtCh}, {axisPtCh}}); + + hist.add("Analysis/Pion/hN1Matrix", "Single-particle density", kTH2D, {{axisCent}, {axisPtPi}}); + hist.add("Analysis/Pion/hN2Matrix", "Ordered-pair density", kTH3D, {{axisCent}, {axisPtPi}, {axisPtPi}}); + + hist.add("Analysis/Kaon/hN1Matrix", "Single-particle density", kTH2D, {{axisCent}, {axisPtKa}}); + hist.add("Analysis/Kaon/hN2Matrix", "Ordered-pair density", kTH3D, {{axisCent}, {axisPtKa}, {axisPtKa}}); + + hist.add("Analysis/Proton/hN1Matrix", "Single-particle density", kTH2D, {{axisCent}, {axisPtPr}}); + hist.add("Analysis/Proton/hN2Matrix", "Ordered-pair density", kTH3D, {{axisCent}, {axisPtPr}, {axisPtPr}}); + }; + addMatrixHistograms(); + + if (cfgApplyEfficiencyCorrection) { + auto addCorrectedHistograms = [&](const std::string& directory, const AxisSpec& axisPtId) { + hist.add(directory + "hN1Matrix", "Efficiency-corrected Single-particle density", kTH2D, {{axisCent}, {axisPtId}}); + hist.add(directory + "hN2Matrix", "Efficiency-corrected Ordered-pair density", kTH3D, {{axisCent}, {axisPtId}, {axisPtId}}); + }; + addCorrectedHistograms("Corrected/Charged/", axisPtCh); + addCorrectedHistograms("Corrected/Pion/", axisPtPi); + addCorrectedHistograms("Corrected/Kaon/", axisPtKa); + addCorrectedHistograms("Corrected/Proton/", axisPtPr); + addCorrectedHistograms("CorrectedTruthMatched/Charged/", axisPtCh); + addCorrectedHistograms("CorrectedTruthMatched/Pion/", axisPtPi); + addCorrectedHistograms("CorrectedTruthMatched/Kaon/", axisPtKa); + addCorrectedHistograms("CorrectedTruthMatched/Proton/", axisPtPr); + } + + auto addGeneratedHistograms = [&]() { + // MC generated + hist.add("Gen/h_Counts", "Counts", kTH1D, {axisEvents}); + hist.add("Gen/h_VtxZ", "Vertex Z ", kTH1D, {axisVtxZ}); + hist.add("Gen/h_VtxZ_b", "Vertex Z ", kTH1D, {axisVtxZ}); + hist.add("Gen/h_NSim", "Truth Multiplicity TPC", kTH1D, {axisMultMC}); + hist.add("Gen/h2_NSim_NTPC", "Reco vs Truth Multiplicty TPC", kTH2D, {{axisMultMC}, {axisMultTPC}}); + hist.add("Gen/h2_NTPC_NSim", "Truth vs Reco Multiplicty TPC", kTH2D, {{axisMultTPC}, {axisMultMC}}); + + hist.add("Gen/Charged/h_PtTruth", "p_{T} ", kTH1D, {axisPt}); + hist.add("Gen/Charged/h_EtaTruth", "#eta ", kTH1D, {axisEta}); + hist.add("Gen/Charged/h_PhiTruth", "#phi ", kTH1D, {axisPhi}); + hist.add("Gen/Charged/h2_PtTruth_centFT0M", "p_{T} in centrality Classes ", kTH2D, {{axisCentFT0M}, {axisPt}}); + hist.add("Gen/Charged/h3_Pt_Eta_PhiTruth", "p_{T}, #eta, #phi ", kTHnSparseD, {{axisPt}, {axisEta}, {axisPhi}}); + + hist.addClone("Gen/Charged/", "Gen/Pion/"); + + hist.add("Gen/Pion/h_RapTruth", "y", kTH1D, {axisY}); + + hist.addClone("Gen/Pion/", "Gen/Kaon/"); + hist.addClone("Gen/Pion/", "Gen/Proton/"); + + hist.addClone("Analysis/Charged/", "Gen/Charged/"); + hist.addClone("Analysis/Pion/", "Gen/Pion/"); + hist.addClone("Analysis/Kaon/", "Gen/Kaon/"); + hist.addClone("Analysis/Proton/", "Gen/Proton/"); + + hist.add("Gen/Charged/h_PtTruth_e", "p_{T} ", kTH1D, {axisPtCh}); + hist.add("Gen/Charged/h2_Pt_EtaTruth", "p_{T} vs #eta", kTH2D, {{axisEta}, {axisPtCh}}); + hist.add("Gen/Charged/h3_Pt_EtaTruth_centFT0M", "p_{T} vs #eta in centrality classes", kTH3D, {{axisEta}, {axisPtCh}, {axisCent}}); + + hist.add("Gen/Pion/h_PtTruth_e", "p_{T} ", kTH1D, {axisPtPi}); + hist.add("Gen/Pion/h_PtPosTruth", "p_{T} (positive) ", kTH1D, {axisPtPi}); + hist.add("Gen/Pion/h_PtNegTruth", "p_{T} (negative) ", kTH1D, {axisPtPi}); + hist.add("Gen/Pion/h2_Pt_EtaTruth", "p_{T} vs #eta", kTH2D, {{axisEta}, {axisPtPi}}); + hist.add("Gen/Pion/h3_Pt_EtaTruth_centFT0M", "p_{T} vs #eta in centrality classes", kTH3D, {{axisEta}, {axisPtPi}, {axisCent}}); + + hist.add("Gen/Kaon/h_PtTruth_e", "p_{T} ", kTH1D, {axisPtKa}); + hist.add("Gen/Kaon/h_PtPosTruth", "p_{T} (positive) ", kTH1D, {axisPtKa}); + hist.add("Gen/Kaon/h_PtNegTruth", "p_{T} (negative) ", kTH1D, {axisPtKa}); + hist.add("Gen/Kaon/h2_Pt_EtaTruth", "p_{T} vs #eta", kTH2D, {{axisEta}, {axisPtKa}}); + hist.add("Gen/Kaon/h3_Pt_EtaTruth_centFT0M", "p_{T} vs #eta in centrality classes", kTH3D, {{axisEta}, {axisPtKa}, {axisCent}}); + + hist.add("Gen/Proton/h_PtTruth_e", "p_{T} ", kTH1D, {axisPtPr}); + hist.add("Gen/Proton/h_PtPosTruth", "p_{T} (positive) ", kTH1D, {axisPtPr}); + hist.add("Gen/Proton/h_PtNegTruth", "p_{T} (negative) ", kTH1D, {axisPtPr}); + hist.add("Gen/Proton/h2_Pt_EtaTruth", "p_{T} vs #eta", kTH2D, {{axisEta}, {axisPtPr}}); + hist.add("Gen/Proton/h3_Pt_EtaTruth_centFT0M", "p_{T} vs #eta in centrality classes", kTH3D, {{axisEta}, {axisPtPr}, {axisCent}}); + }; + addGeneratedHistograms(); hist.add("QA/h_collisions_info", "Collisions info", kTH1D, {axisCol}); hist.add("Gen/h_collisions_info", "Collisions info", kTH1D, {axisCol}); @@ -556,11 +662,12 @@ struct MeanPtFlucId { return false; } - if (std::fabs(track.dcaZ()) > cfgCutDcaZ) { - return false; - } + const float maxDca = cfgDcaPtDepOffset + cfgDcaPtDepScale / track.pt(); + // if (std::abs(track.dcaXY()) >= maxDca) { + // return false; + // } - if (std::fabs(track.dcaZ()) > (0.0105 + 0.035 / std::pow(track.p(), 1.1))) { + if (std::abs(track.dcaZ()) >= maxDca) { return false; } @@ -586,62 +693,50 @@ struct MeanPtFlucId { template bool identifyParticle(T const& track, PIDType species, PIDType reject1, PIDType reject2, float p, float momThreshold) { - std::vector vTpcNSigma = {-999.f, track.tpcNSigmaPi(), track.tpcNSigmaKa(), track.tpcNSigmaPr()}; - std::vector vTofNSigma = {-999.f, track.tofNSigmaPi(), track.tofNSigmaKa(), track.tofNSigmaPr()}; - - const int sp = static_cast(species); - const int sq = static_cast(reject1); - const int sr = static_cast(reject2); + const std::array tpcNSigma = {-999.f, track.tpcNSigmaPi(), track.tpcNSigmaKa(), track.tpcNSigmaPr()}; + const std::array tofNSigma = {-999.f, track.tofNSigmaPi(), track.tofNSigmaKa(), track.tofNSigmaPr()}; - bool isTofPidFlag = false; - bool isTpcPidFlag = false; + const int selSpecies = static_cast(species); + const int rejSpecies1 = static_cast(reject1); + const int rejSpecies2 = static_cast(reject2); - float nSigmaSelCut = 0.f; + float tpcSelectionCut = 0.f; + float tofSelectionCut = 0.f; switch (species) { case PIDType::kPions: - nSigmaSelCut = cfgSelCutNSigPi; + tpcSelectionCut = cfgSelCutNSigTpcPi; + tofSelectionCut = cfgSelCutNSigTofPi; break; - case PIDType::kKaons: - nSigmaSelCut = cfgSelCutNSigKa; + tpcSelectionCut = cfgSelCutNSigTpcKa; + tofSelectionCut = cfgSelCutNSigTofKa; break; - case PIDType::kProtons: - nSigmaSelCut = cfgSelCutNSigPr; + tpcSelectionCut = cfgSelCutNSigTpcPr; + tofSelectionCut = cfgSelCutNSigTofPr; break; - default: return false; } if (track.hasTOF()) { + const float selectedTofNSigma = std::abs(tofNSigma[selSpecies]); + const bool passesTOF = selectedTofNSigma < tofSelectionCut && + std::abs(tofNSigma[rejSpecies1]) > selectedTofNSigma && + std::abs(tofNSigma[rejSpecies2]) > selectedTofNSigma; + const bool passesTPC = std::abs(tpcNSigma[selSpecies]) < cfgCutNSigTpcTof; + return passesTPC && passesTOF; + } - if (std::abs(vTofNSigma[sp]) < nSigmaSelCut && - std::abs(vTofNSigma[sq]) > std::abs(vTofNSigma[sp]) && - std::abs(vTofNSigma[sr]) > std::abs(vTofNSigma[sp])) { - isTofPidFlag = true; - } - - if (std::abs(vTpcNSigma[sp]) < cfgCutNSig2) { - isTpcPidFlag = true; - } - - } else { - - if (p >= momThreshold) { - return false; - } - - if (std::abs(vTpcNSigma[sp]) < nSigmaSelCut && - std::abs(vTpcNSigma[sq]) > std::abs(vTpcNSigma[sp]) && - std::abs(vTpcNSigma[sr]) > std::abs(vTpcNSigma[sp])) { - isTofPidFlag = true; - isTpcPidFlag = true; - } + if (p >= momThreshold) { + return false; } - return (isTofPidFlag && isTpcPidFlag); + const float selTpcNSigma = std::abs(tpcNSigma[selSpecies]); + return selTpcNSigma < tpcSelectionCut && + std::abs(tpcNSigma[rejSpecies1]) > selTpcNSigma && + std::abs(tpcNSigma[rejSpecies2]) > selTpcNSigma; } // Fill QA histograms before selection cuts: @@ -706,82 +801,95 @@ struct MeanPtFlucId { template void fillBeforePIDQAHistos(T const& track) { - hist.fill(HIST("QA/before/h2_TOFSignal"), track.p(), track.beta()); hist.fill(HIST("QA/before/h2_TPCSignal"), track.p(), track.tpcSignal()); - hist.fill(HIST("QA/before/h2_pvsm2"), track.mass() * track.mass(), track.p()); + if (track.hasTOF()) { + hist.fill(HIST("QA/before/h2_TOFSignal"), track.p(), track.beta()); + hist.fill(HIST("QA/before/h2_pvsm2"), track.mass() * track.mass(), track.p()); + } hist.fill(HIST("QA/Pion/before/h2_TPCNsigma"), track.p(), track.tpcNSigmaPi()); - if (!track.hasTOF()) + if (!track.hasTOF()) { hist.fill(HIST("QA/Pion/before/h2_TPCNsigma_nottof"), track.p(), track.tpcNSigmaPi()); - if (track.hasTOF()) + } + if (track.hasTOF()) { hist.fill(HIST("QA/Pion/before/h2_TPCNsigma_tof"), track.p(), track.tpcNSigmaPi()); - hist.fill(HIST("QA/Pion/before/h2_TOFNsigma"), track.p(), track.tofNSigmaPi()); - hist.fill(HIST("QA/Pion/before/h2_TpcTofNsigma"), track.tpcNSigmaPi(), track.tofNSigmaPi()); + hist.fill(HIST("QA/Pion/before/h2_TOFNsigma"), track.p(), track.tofNSigmaPi()); + hist.fill(HIST("QA/Pion/before/h2_TpcTofNsigma"), track.tpcNSigmaPi(), track.tofNSigmaPi()); + } hist.fill(HIST("QA/Kaon/before/h2_TPCNsigma"), track.p(), track.tpcNSigmaKa()); - if (!track.hasTOF()) + if (!track.hasTOF()) { hist.fill(HIST("QA/Kaon/before/h2_TPCNsigma_nottof"), track.p(), track.tpcNSigmaKa()); - if (track.hasTOF()) + } + if (track.hasTOF()) { hist.fill(HIST("QA/Kaon/before/h2_TPCNsigma_tof"), track.p(), track.tpcNSigmaKa()); - hist.fill(HIST("QA/Kaon/before/h2_TOFNsigma"), track.p(), track.tofNSigmaKa()); - hist.fill(HIST("QA/Kaon/before/h2_TpcTofNsigma"), track.tpcNSigmaKa(), track.tofNSigmaKa()); + hist.fill(HIST("QA/Kaon/before/h2_TOFNsigma"), track.p(), track.tofNSigmaKa()); + hist.fill(HIST("QA/Kaon/before/h2_TpcTofNsigma"), track.tpcNSigmaKa(), track.tofNSigmaKa()); + } hist.fill(HIST("QA/Proton/before/h2_TPCNsigma"), track.p(), track.tpcNSigmaPr()); - if (!track.hasTOF()) + if (!track.hasTOF()) { hist.fill(HIST("QA/Proton/before/h2_TPCNsigma_nottof"), track.p(), track.tpcNSigmaPr()); - if (track.hasTOF()) + } + if (track.hasTOF()) { hist.fill(HIST("QA/Proton/before/h2_TPCNsigma_tof"), track.p(), track.tpcNSigmaPr()); - hist.fill(HIST("QA/Proton/before/h2_TOFNsigma"), track.p(), track.tofNSigmaPr()); - hist.fill(HIST("QA/Proton/before/h2_TpcTofNsigma"), track.tpcNSigmaPr(), track.tofNSigmaPr()); + hist.fill(HIST("QA/Proton/before/h2_TOFNsigma"), track.p(), track.tofNSigmaPr()); + hist.fill(HIST("QA/Proton/before/h2_TpcTofNsigma"), track.tpcNSigmaPr(), track.tofNSigmaPr()); + } } // Fill identified QA histograms after PID cuts: template - void fillIdParticleQAHistos(T const& track, float ptMC, float etaMC, float rap, float nSigmaTPC, float nSigmaTOF, std::vector& Tracks) + void fillIdParticleQAHistos(T const& track, float ptMC, float etaMC, float rap, float nSigmaTPC, float nSigmaTOF, std::vector& tracks) { float pt = track.pt(); float eta = track.eta(); float phi = track.phi(); - Tracks.push_back({pt, eta}); + tracks.push_back({pt}); if (track.sign() > 0) { - hist.fill(HIST(Dire[Mode]) + HIST("h_PtPos"), pt); + hist.fill(HIST(dire[Mode]) + HIST("h_PtPos"), pt); } if (track.sign() < 0) { - hist.fill(HIST(Dire[Mode]) + HIST("h_PtNeg"), pt); + hist.fill(HIST(dire[Mode]) + HIST("h_PtNeg"), pt); } - hist.fill(HIST(Dire[Mode]) + HIST("h_Pt"), pt); - hist.fill(HIST(Dire[Mode]) + HIST("h_PtMC"), ptMC); - hist.fill(HIST(Dire[Mode]) + HIST("h2_Pt_PtMC"), ptMC, pt); - hist.fill(HIST(Dire[Mode]) + HIST("h2_Pt_centFT0M"), centFT0M, pt); - hist.fill(HIST(Dire[Mode]) + HIST("h3_Pt_Eta_Phi"), pt, eta, phi); - hist.fill(HIST(Dire[Mode]) + HIST("h2_Pt_Eta"), eta, pt); - hist.fill(HIST(Dire[Mode]) + HIST("h3_Pt_Eta_centFT0M"), eta, pt, centFT0M); - hist.fill(HIST(Dire[Mode]) + HIST("h2_Pt_EtaMC"), etaMC, ptMC); - hist.fill(HIST(Dire[Mode]) + HIST("h3_Pt_EtaMC_centFT0M"), etaMC, ptMC, centFT0M); - hist.fill(HIST(Dire[Mode]) + HIST("h_Eta"), eta); - hist.fill(HIST(Dire[Mode]) + HIST("h_Phi"), phi); - hist.fill(HIST(Dire[Mode]) + HIST("h_Rap"), rap); - hist.fill(HIST(Dire[Mode]) + HIST("h_DcaZ"), track.dcaZ()); - hist.fill(HIST(Dire[Mode]) + HIST("h_DcaXY"), track.dcaXY()); - hist.fill(HIST(Dire[Mode]) + HIST("h2_DcaZ"), pt, track.dcaZ()); - hist.fill(HIST(Dire[Mode]) + HIST("h2_DcaXY"), pt, track.dcaXY()); - - hist.fill(HIST(Dire[Mode]) + HIST("h2_TPCNsigma"), track.p(), nSigmaTPC); - hist.fill(HIST(Dire[Mode]) + HIST("h2_TOFNsigma"), track.p(), nSigmaTOF); - hist.fill(HIST(Dire[Mode]) + HIST("h2_TpcTofNsigma"), nSigmaTPC, nSigmaTOF); - hist.fill(HIST(Dire[Mode]) + HIST("h2_TPCSignal"), track.p(), track.tpcSignal()); - hist.fill(HIST(Dire[Mode]) + HIST("innerParam/h2_TPCSignal"), track.tpcInnerParam(), track.tpcSignal()); - hist.fill(HIST(Dire[Mode]) + HIST("h2_TOFSignal"), track.p(), track.beta()); - hist.fill(HIST(Dire[Mode]) + HIST("h2_pvsm2"), track.mass() * track.mass(), track.p()); + hist.fill(HIST(dire[Mode]) + HIST("h_Pt"), pt); + hist.fill(HIST(dire[Mode]) + HIST("h_Pt_e"), pt); + hist.fill(HIST(dire[Mode]) + HIST("h_PtMC"), ptMC); + hist.fill(HIST(dire[Mode]) + HIST("h2_Pt_PtMC"), ptMC, pt); + hist.fill(HIST(dire[Mode]) + HIST("h2_Pt_centFT0M"), centFT0M, pt); + hist.fill(HIST(dire[Mode]) + HIST("h3_Pt_Eta_Phi"), pt, eta, phi); + hist.fill(HIST(dire[Mode]) + HIST("h2_Pt_Eta"), eta, pt); + hist.fill(HIST(dire[Mode]) + HIST("h3_Pt_Eta_centFT0M"), eta, pt, centFT0M); + hist.fill(HIST(dire[Mode]) + HIST("h2_Pt_EtaMC"), etaMC, ptMC); + hist.fill(HIST(dire[Mode]) + HIST("h3_Pt_EtaMC_centFT0M"), etaMC, ptMC, centFT0M); + hist.fill(HIST(dire[Mode]) + HIST("h_Eta"), eta); + hist.fill(HIST(dire[Mode]) + HIST("h_Phi"), phi); + hist.fill(HIST(dire[Mode]) + HIST("h_Rap"), rap); + hist.fill(HIST(dire[Mode]) + HIST("h_DcaZ"), track.dcaZ()); + hist.fill(HIST(dire[Mode]) + HIST("h_DcaXY"), track.dcaXY()); + hist.fill(HIST(dire[Mode]) + HIST("h2_DcaZ"), pt, track.dcaZ()); + hist.fill(HIST(dire[Mode]) + HIST("h2_DcaXY"), pt, track.dcaXY()); + + hist.fill(HIST(dire[Mode]) + HIST("h2_TPCNsigma"), track.p(), nSigmaTPC); + hist.fill(HIST(dire[Mode]) + HIST("h2_TPCSignal"), track.p(), track.tpcSignal()); + hist.fill(HIST(dire[Mode]) + HIST("innerParam/h2_TPCSignal"), track.tpcInnerParam(), track.tpcSignal()); + if (track.hasTOF()) { + hist.fill(HIST(dire[Mode]) + HIST("h2_TOFNsigma"), track.p(), nSigmaTOF); + hist.fill(HIST(dire[Mode]) + HIST("h2_TpcTofNsigma"), nSigmaTPC, nSigmaTOF); + hist.fill(HIST(dire[Mode]) + HIST("h2_TOFSignal"), track.p(), track.beta()); + hist.fill(HIST(dire[Mode]) + HIST("h2_pvsm2"), track.mass() * track.mass(), track.p()); + } hist.fill(HIST("QA/after/h2_TPCSignal"), track.p(), track.tpcSignal()); hist.fill(HIST("QA/after/innerParam/h2_TPCSignal"), track.tpcInnerParam(), track.tpcSignal()); - hist.fill(HIST("QA/after/h2_TOFSignal"), track.p(), track.beta()); - hist.fill(HIST("QA/after/h2_pvsm2"), track.mass() * track.mass(), track.p()); + if (track.hasTOF()) { + hist.fill(HIST("QA/after/h2_TOFSignal"), track.p(), track.beta()); + hist.fill(HIST("QA/after/h2_pvsm2"), track.mass() * track.mass(), track.p()); + } } // Fill Truth identified particles histograms @@ -789,142 +897,158 @@ struct MeanPtFlucId { void fillPtMCHist(bool cfgGen, float pt, float eta, float rap, float phi, int pid, int pdgCodePos, int pdgCodeNeg) { if (cfgGen) { - hist.fill(HIST(Dire[Mode]) + HIST("h_PtTruth"), pt); - hist.fill(HIST(Dire[Mode]) + HIST("h_EtaTruth"), eta); - hist.fill(HIST(Dire[Mode]) + HIST("h_RapTruth"), rap); - hist.fill(HIST(Dire[Mode]) + HIST("h_PhiTruth"), phi); - hist.fill(HIST(Dire[Mode]) + HIST("h2_PtTruth_centFT0M"), centFT0M, pt); - hist.fill(HIST(Dire[Mode]) + HIST("h3_Pt_EtaTruth_centFT0M"), eta, pt, centFT0M); - hist.fill(HIST(Dire[Mode]) + HIST("h3_Pt_Eta_PhiTruth"), pt, eta, phi); + hist.fill(HIST(dire[Mode]) + HIST("h_PtTruth"), pt); + hist.fill(HIST(dire[Mode]) + HIST("h_PtTruth_e"), pt); + hist.fill(HIST(dire[Mode]) + HIST("h_EtaTruth"), eta); + hist.fill(HIST(dire[Mode]) + HIST("h_RapTruth"), rap); + hist.fill(HIST(dire[Mode]) + HIST("h_PhiTruth"), phi); + hist.fill(HIST(dire[Mode]) + HIST("h2_PtTruth_centFT0M"), centFT0M, pt); + hist.fill(HIST(dire[Mode]) + HIST("h3_Pt_EtaTruth_centFT0M"), eta, pt, centFT0M); + hist.fill(HIST(dire[Mode]) + HIST("h3_Pt_Eta_PhiTruth"), pt, eta, phi); } - hist.fill(HIST(Dire[Mode]) + HIST("h2_Pt_EtaTruth"), eta, pt); + hist.fill(HIST(dire[Mode]) + HIST("h2_Pt_EtaTruth"), eta, pt); if (pid == pdgCodePos) { - hist.fill(HIST(Dire[Mode]) + HIST("h_PtPosTruth"), pt); + hist.fill(HIST(dire[Mode]) + HIST("h_PtPosTruth"), pt); } if (pid == pdgCodeNeg) { - hist.fill(HIST(Dire[Mode]) + HIST("h_PtNegTruth"), pt); - } - } - - int getPtBin(double pt, const std::vector& ptB) - { - for (size_t i = 0; i < ptB.size() - 1; ++i) { - if (pt >= ptB[i] && pt < ptB[i + 1]) { - return i; - } - } - return -1; - } - - int getEtaBin(double eta) - { - for (size_t i = 0; i < etaBins->size() - 1; ++i) { - if (eta >= etaBins->at(i) && eta < etaBins->at(i + 1)) { - return i; - } + hist.fill(HIST(dire[Mode]) + HIST("h_PtNegTruth"), pt); } - return -1; } - // fill differential matrices for two-particle analysis template - void fillDifferentialMatrices(const std::vector& tracks, const std::vector& ptB) + void fillMatrices(const std::vector& tracks, double ptMin) { - const int nPtBins = ptB.size() - 1; - - auto getK = [&](int ipt, int ieta) { - return ieta * nPtBins + ipt; - }; + const int nBins = nPtBins; + const double binWidth = (cfgCutPtMax - ptMin) / nBins; + matrixNInBin.assign(nBins, 0); + matrixOccupiedBins.clear(); + matrixOccupiedBins.reserve(std::min(static_cast(nBins), tracks.size())); - for (const auto& trk : tracks) { - int ipt = getPtBin(trk.pt, ptB); - int ieta = getEtaBin(trk.eta); - - if (ipt < 0 || ieta < 0) + for (const auto& track : tracks) { + const int ptBin = static_cast((track.pt - ptMin) / binWidth); + if (ptBin < 0 || ptBin >= nBins) { continue; - - int k = getK(ipt, ieta); - - hist.fill(HIST(Dire[Mode]) + HIST("hN1Matrix"), centFT0M, k, 1.0); - hist.fill(HIST(Dire[Mode]) + HIST("hPt1Matrix"), centFT0M, k, trk.pt); + } + if (matrixNInBin[ptBin] == 0) { + matrixOccupiedBins.push_back(ptBin); + } + ++matrixNInBin[ptBin]; } - for (size_t i = 0; i < tracks.size(); ++i) { - int ipt1 = getPtBin(tracks[i].pt, ptB); - int ieta1 = getEtaBin(tracks[i].eta); - - if (ipt1 < 0 || ieta1 < 0) - continue; - - int k1 = getK(ipt1, ieta1); - - for (size_t j = i + 1; j < tracks.size(); ++j) { - int ipt2 = getPtBin(tracks[j].pt, ptB); - int ieta2 = getEtaBin(tracks[j].eta); - - if (ipt2 < 0 || ieta2 < 0) - continue; + for (const int& ptBin : matrixOccupiedBins) { + const double ptCenter = ptMin + (ptBin + 0.5) * binWidth; + hist.fill(HIST(dire[Mode]) + HIST("hN1Matrix"), centFT0M, ptCenter, static_cast(matrixNInBin[ptBin])); + } - int k2 = getK(ipt2, ieta2); + for (size_t i = 0; i < matrixOccupiedBins.size(); ++i) { + const int ptBin1 = matrixOccupiedBins[i]; + const double ptCenter1 = ptMin + (ptBin1 + 0.5) * binWidth; + const double nPairDiagonal = static_cast(matrixNInBin[ptBin1]) * (matrixNInBin[ptBin1] - 1); + if (nPairDiagonal > 0.0) { + hist.fill(HIST(dire[Mode]) + HIST("hN2Matrix"), centFT0M, ptCenter1, ptCenter1, nPairDiagonal); + } - hist.fill(HIST(Dire[Mode]) + HIST("hN2Matrix"), centFT0M, k1, k2, 1.0); - hist.fill(HIST(Dire[Mode]) + HIST("hN2Matrix"), centFT0M, k2, k1, 1.0); - hist.fill(HIST(Dire[Mode]) + HIST("hPtPtMatrix"), centFT0M, k1, k2, tracks[i].pt * tracks[j].pt); - hist.fill(HIST(Dire[Mode]) + HIST("hPtPtMatrix"), centFT0M, k2, k1, tracks[j].pt * tracks[i].pt); + for (size_t j = i + 1; j < matrixOccupiedBins.size(); ++j) { + const int ptBin2 = matrixOccupiedBins[j]; + const double ptCenter2 = ptMin + (ptBin2 + 0.5) * binWidth; + const double nPair = static_cast(matrixNInBin[ptBin1]) * matrixNInBin[ptBin2]; + hist.fill(HIST(dire[Mode]) + HIST("hN2Matrix"), centFT0M, ptCenter1, ptCenter2, nPair); + hist.fill(HIST(dire[Mode]) + HIST("hN2Matrix"), centFT0M, ptCenter2, ptCenter1, nPair); } } } - // Calculate moments for two-particle correlations - void moments(double pt, int& n, double& q1, double& q2) - { - n++; - q1 += pt; - q2 += pt * pt; - } - - // Fill terms calculted from moments method for two-particle correlations - template - void fillAnalysis(int n, double q1, double q2) + bool getEfficiencyWeight(PIDType species, double pt, double eta, double& weight) const { + TH3* weightMap = nullptr; + switch (species) { + case PIDType::kNone: + weightMap = efficiencyWeightCh; + break; + case PIDType::kPions: + weightMap = efficiencyWeightPi; + break; + case PIDType::kKaons: + weightMap = efficiencyWeightKa; + break; + case PIDType::kProtons: + weightMap = efficiencyWeightPr; + break; + } + if (!weightMap) { + return false; + } - double npair = n * (n - 1); - - if (n > 0) { - double meanPt = q1 / n; - - hist.fill(HIST(Dire[Mode]) + HIST("h_mult"), n); - - hist.fill(HIST(Dire[Mode]) + HIST("h_mean_pT_Cent"), centFT0M, meanPt); - hist.fill(HIST(Dire[Mode]) + HIST("p_mean_pT_Cent"), centFT0M, meanPt); + const int etaBin = weightMap->GetXaxis()->FindFixBin(eta); + const int ptBin = weightMap->GetYaxis()->FindFixBin(pt); + const int centBin = weightMap->GetZaxis()->FindFixBin(centFT0M); + if (etaBin < 1 || etaBin > weightMap->GetNbinsX() || + ptBin < 1 || ptBin > weightMap->GetNbinsY() || + centBin < 1 || centBin > weightMap->GetNbinsZ()) { + return false; + } - hist.fill(HIST(Dire[Mode]) + HIST("p_mean_pT_Num_Cent_var"), centFT0M, q1); - hist.fill(HIST(Dire[Mode]) + HIST("p_mean_pT_Den_Cent_var"), centFT0M, n); + weight = weightMap->GetBinContent(etaBin, ptBin, centBin); + return std::isfinite(weight) && weight > 0.; + } - hist.fill(HIST(Dire[Mode]) + HIST("h_mean_pT_Num_Cent_var"), centFT0M, q1); - hist.fill(HIST(Dire[Mode]) + HIST("h_mean_pT_Den_Cent_var"), centFT0M, n); + template + void fillCorrectedMatrices(const std::vector& tracks, double ptMin) + { + const int nBins = nPtBins; + const double binWidth = (cfgCutPtMax - ptMin) / nBins; + matrixWInBin.assign(nBins, 0.); + matrixW2InBin.assign(nBins, 0.); + matrixOccupiedBins.clear(); + matrixOccupiedBins.reserve(std::min(static_cast(nBins), tracks.size())); - hist.fill(HIST(Dire[Mode]) + HIST("p_Q1_Mult"), n, q1); - hist.fill(HIST(Dire[Mode]) + HIST("p_Q2_Mult"), n, q2); - hist.fill(HIST(Dire[Mode]) + HIST("p_mean_pT_Mult"), n, meanPt); + for (const auto& track : tracks) { + const int ptBin = static_cast((track.pt - ptMin) / binWidth); + if (ptBin < 0 || ptBin >= nBins) { + continue; + } + if (matrixWInBin[ptBin] == 0.) { + matrixOccupiedBins.push_back(ptBin); + } + matrixWInBin[ptBin] += track.weight; + matrixW2InBin[ptBin] += track.weight * track.weight; } - if (npair != 0) { - double meanPt = q1 / n; - double twopart = (q1 * q1 - q2); - double varPt = twopart / npair; - hist.fill(HIST(Dire[Mode]) + HIST("p_twopart_Num_Cent_var"), centFT0M, twopart); - hist.fill(HIST(Dire[Mode]) + HIST("p_twopart_Den_Cent_var"), centFT0M, npair); + for (const int& ptBin : matrixOccupiedBins) { + const double ptCenter = ptMin + (ptBin + 0.5) * binWidth; + if constexpr (TruthMatched) { + hist.fill(HIST(TruthMatchedCorrectedDirectories[Mode]) + HIST("hN1Matrix"), centFT0M, ptCenter, matrixWInBin[ptBin]); + } else { + hist.fill(HIST(CorrectedDirectories[Mode]) + HIST("hN1Matrix"), centFT0M, ptCenter, matrixWInBin[ptBin]); + } + } - hist.fill(HIST(Dire[Mode]) + HIST("h_twopart_Num_Cent_var"), centFT0M, twopart); - hist.fill(HIST(Dire[Mode]) + HIST("h_twopart_Den_Cent_var"), centFT0M, npair); + for (size_t i = 0; i < matrixOccupiedBins.size(); ++i) { + const int ptBin1 = matrixOccupiedBins[i]; + const double ptCenter1 = ptMin + (ptBin1 + 0.5) * binWidth; + const double diagonalPairDensity = matrixWInBin[ptBin1] * matrixWInBin[ptBin1] - matrixW2InBin[ptBin1]; + if (diagonalPairDensity > 0.) { + if constexpr (TruthMatched) { + hist.fill(HIST(TruthMatchedCorrectedDirectories[Mode]) + HIST("hN2Matrix"), centFT0M, ptCenter1, ptCenter1, diagonalPairDensity); + } else { + hist.fill(HIST(CorrectedDirectories[Mode]) + HIST("hN2Matrix"), centFT0M, ptCenter1, ptCenter1, diagonalPairDensity); + } + } - hist.fill(HIST(Dire[Mode]) + HIST("p_mean_pT_Cent_var"), centFT0M, meanPt); - hist.fill(HIST(Dire[Mode]) + HIST("p_twopart_Cent_var"), centFT0M, varPt); - hist.fill(HIST(Dire[Mode]) + HIST("h_twopart_Cent_var"), centFT0M, varPt); - hist.fill(HIST(Dire[Mode]) + HIST("h_mean_pT_Cent_var"), centFT0M, meanPt); + for (size_t j = i + 1; j < matrixOccupiedBins.size(); ++j) { + const int ptBin2 = matrixOccupiedBins[j]; + const double ptCenter2 = ptMin + (ptBin2 + 0.5) * binWidth; + const double pairDensity = matrixWInBin[ptBin1] * matrixWInBin[ptBin2]; + if constexpr (TruthMatched) { + hist.fill(HIST(TruthMatchedCorrectedDirectories[Mode]) + HIST("hN2Matrix"), centFT0M, ptCenter1, ptCenter2, pairDensity); + hist.fill(HIST(TruthMatchedCorrectedDirectories[Mode]) + HIST("hN2Matrix"), centFT0M, ptCenter2, ptCenter1, pairDensity); + } else { + hist.fill(HIST(CorrectedDirectories[Mode]) + HIST("hN2Matrix"), centFT0M, ptCenter1, ptCenter2, pairDensity); + hist.fill(HIST(CorrectedDirectories[Mode]) + HIST("hN2Matrix"), centFT0M, ptCenter2, ptCenter1, pairDensity); + } + } } } @@ -936,10 +1060,16 @@ struct MeanPtFlucId { pionTracks.clear(); kaonTracks.clear(); protonTracks.clear(); + chargedTracksCorrected.clear(); + pionTracksCorrected.clear(); + kaonTracksCorrected.clear(); + protonTracksCorrected.clear(); + chargedTracksTruthMatchedCorrected.clear(); + pionTracksTruthMatchedCorrected.clear(); + kaonTracksTruthMatchedCorrected.clear(); + protonTracksTruthMatchedCorrected.clear(); - int nCh = 0, nPi = 0, nKa = 0, nPr = 0; float pt = 0., eta = 0., phi = 0., p = 0.; - double q1Ch = 0., q2Ch = 0., q1Pi = 0., q2Pi = 0., q1Ka = 0., q2Ka = 0., q1Pr = 0., q2Pr = 0.; hist.fill(HIST("QA/h_collisions_info"), kTotCol); @@ -949,14 +1079,6 @@ struct MeanPtFlucId { } } - if constexpr (RecoFlag) { - float lower = cfgP2L * nSim * nSim + cfgP1L * nSim - cfgP0L; - float upper = cfgP2U * nSim * nSim + cfgP1U * nSim + cfgP0U; - if (cfgNtpcEventCut && (nTPC < lower || nTPC > upper)) { - return; // Reject event - } - } - hist.fill(HIST("QA/h_collisions_info"), kPassSelCol); fillAfterQAHistos(col); @@ -1002,22 +1124,39 @@ struct MeanPtFlucId { auto mc = track.template mcParticle_as(); ptMC = mc.pt(); etaMC = mc.eta(); - p = mc.p(); if (mc.isPhysicalPrimary()) { hist.fill(HIST("QA/after/h_DCAxy_primary"), track.dcaXY()); hist.fill(HIST("QA/after/h_DCAz_primary"), track.dcaZ()); + hist.fill(HIST("QA/Charged/h_PtTruth_primary"), track.pt()); + const bool mcInChargedAcceptance = ptMC >= cfgCutPtMin && ptMC < cfgCutPtMax && std::abs(etaMC) < cfgCutEta; + if (mcInChargedAcceptance) { + hist.fill(HIST("QA/Charged/h3_Pt_EtaMC_primary_centFT0M"), etaMC, ptMC, centFT0M); + if (cfgApplyEfficiencyCorrection && !cfgUseCombinedPurityEfficiencyWeight) { + double weight = 0.; + if (!getEfficiencyWeight(PIDType::kNone, ptMC, etaMC, weight)) { + LOGF(fatal, "Undefined truth-matched charged efficiency weight at pT=%g, eta=%g, centrality=%g", ptMC, etaMC, centFT0M); + } + chargedTracksTruthMatchedCorrected.push_back({ptMC, weight}); + } + } } else { hist.fill(HIST("QA/after/h_DCAxy_secondary"), track.dcaXY()); hist.fill(HIST("QA/after/h_DCAz_secondary"), track.dcaZ()); } } - nCh++; - chargedTracks.push_back({pt, eta}); - moments(pt, nCh, q1Ch, q2Ch); + if (cfgApplyEfficiencyCorrection) { + double weight = 0.; + if (!getEfficiencyWeight(PIDType::kNone, pt, eta, weight)) { + LOGF(fatal, "Undefined charged-particle efficiency weight at pT=%g, eta=%g, centrality=%g", pt, eta, centFT0M); + } + chargedTracksCorrected.push_back({pt, weight}); + } + chargedTracks.push_back({pt}); hist.fill(HIST("QA/Charged/h_Pt"), pt); + hist.fill(HIST("QA/Charged/h_Pt_e"), pt); hist.fill(HIST("QA/Charged/h_PtMC"), ptMC); hist.fill(HIST("QA/Charged/h2_Pt_PtMC"), ptMC, pt); hist.fill(HIST("QA/Charged/h2_Pt_EtaMC"), etaMC, ptMC); @@ -1035,15 +1174,33 @@ struct MeanPtFlucId { auto selIDProton = identifyParticle(track, PIDType::kProtons, PIDType::kPions, PIDType::kKaons, p, cfgCutPrThrsldP); if (selIDPion && track.pt() >= cfgCutPiPtMin) { - moments(pt, nPi, q1Pi, q2Pi); + if (cfgApplyEfficiencyCorrection) { + double weight = 0.; + if (!getEfficiencyWeight(PIDType::kPions, pt, eta, weight)) { + LOGF(fatal, "Undefined pion efficiency weight at pT=%g, eta=%g, centrality=%g", pt, eta, centFT0M); + } + pionTracksCorrected.push_back({pt, weight}); + } fillIdParticleQAHistos(track, ptMC, etaMC, rapPi, nSigmaTPCPi, nSigmaTOFPi, pionTracks); } if (selIDKaon && track.pt() >= cfgCutKaPtMin) { - moments(pt, nKa, q1Ka, q2Ka); + if (cfgApplyEfficiencyCorrection) { + double weight = 0.; + if (!getEfficiencyWeight(PIDType::kKaons, pt, eta, weight)) { + LOGF(fatal, "Undefined kaon efficiency weight at pT=%g, eta=%g, centrality=%g", pt, eta, centFT0M); + } + kaonTracksCorrected.push_back({pt, weight}); + } fillIdParticleQAHistos(track, ptMC, etaMC, rapKa, nSigmaTPCKa, nSigmaTOFKa, kaonTracks); } if (selIDProton && track.pt() >= cfgCutPrPtMin) { - moments(pt, nPr, q1Pr, q2Pr); + if (cfgApplyEfficiencyCorrection) { + double weight = 0.; + if (!getEfficiencyWeight(PIDType::kProtons, pt, eta, weight)) { + LOGF(fatal, "Undefined proton efficiency weight at pT=%g, eta=%g, centrality=%g", pt, eta, centFT0M); + } + protonTracksCorrected.push_back({pt, weight}); + } fillIdParticleQAHistos(track, ptMC, etaMC, rapPr, nSigmaTPCPr, nSigmaTOFPr, protonTracks); } @@ -1056,6 +1213,18 @@ struct MeanPtFlucId { fillPtMCHist(false, pt, eta, rapPi, phi, pid, kPiPlus, kPiMinus); if (mc.isPhysicalPrimary()) { hist.fill(HIST("QA/Pion/h_PtTruth_primary"), track.pt()); + hist.fill(HIST("QA/Pion/h2_Pt_EtaTruth_primary"), track.eta(), track.pt()); + const bool mcInPionAcceptance = ptMC >= cfgCutPiPtMin && ptMC < cfgCutPtMax && std::abs(etaMC) < cfgCutEta; + if (mcInPionAcceptance) { + hist.fill(HIST("QA/Pion/h3_Pt_EtaMC_primary_centFT0M"), etaMC, ptMC, centFT0M); + if (cfgApplyEfficiencyCorrection && !cfgUseCombinedPurityEfficiencyWeight) { + double weight = 0.; + if (!getEfficiencyWeight(PIDType::kPions, ptMC, etaMC, weight)) { + LOGF(fatal, "Undefined truth-matched pion efficiency weight at pT=%g, eta=%g, centrality=%g", ptMC, etaMC, centFT0M); + } + pionTracksTruthMatchedCorrected.push_back({ptMC, weight}); + } + } } else { hist.fill(HIST("QA/Pion/h_PtTruth_secondary"), track.pt()); } @@ -1067,6 +1236,18 @@ struct MeanPtFlucId { fillPtMCHist(false, pt, eta, rapKa, phi, pid, kKPlus, kKMinus); if (mc.isPhysicalPrimary()) { hist.fill(HIST("QA/Kaon/h_PtTruth_primary"), track.pt()); + hist.fill(HIST("QA/Kaon/h2_Pt_EtaTruth_primary"), track.eta(), track.pt()); + const bool mcInKaonAcceptance = ptMC >= cfgCutKaPtMin && ptMC < cfgCutPtMax && std::abs(etaMC) < cfgCutEta; + if (mcInKaonAcceptance) { + hist.fill(HIST("QA/Kaon/h3_Pt_EtaMC_primary_centFT0M"), etaMC, ptMC, centFT0M); + if (cfgApplyEfficiencyCorrection && !cfgUseCombinedPurityEfficiencyWeight) { + double weight = 0.; + if (!getEfficiencyWeight(PIDType::kKaons, ptMC, etaMC, weight)) { + LOGF(fatal, "Undefined truth-matched kaon efficiency weight at pT=%g, eta=%g, centrality=%g", ptMC, etaMC, centFT0M); + } + kaonTracksTruthMatchedCorrected.push_back({ptMC, weight}); + } + } } else { hist.fill(HIST("QA/Kaon/h_PtTruth_secondary"), track.pt()); } @@ -1078,6 +1259,18 @@ struct MeanPtFlucId { fillPtMCHist(false, pt, eta, rapPr, phi, pid, kProton, kProtonBar); if (mc.isPhysicalPrimary()) { hist.fill(HIST("QA/Proton/h_PtTruth_primary"), track.pt()); + hist.fill(HIST("QA/Proton/h2_Pt_EtaTruth_primary"), track.eta(), track.pt()); + const bool mcInProtonAcceptance = ptMC >= cfgCutPrPtMin && ptMC < cfgCutPtMax && std::abs(etaMC) < cfgCutEta; + if (mcInProtonAcceptance) { + hist.fill(HIST("QA/Proton/h3_Pt_EtaMC_primary_centFT0M"), etaMC, ptMC, centFT0M); + if (cfgApplyEfficiencyCorrection && !cfgUseCombinedPurityEfficiencyWeight) { + double weight = 0.; + if (!getEfficiencyWeight(PIDType::kProtons, ptMC, etaMC, weight)) { + LOGF(fatal, "Undefined truth-matched proton efficiency weight at pT=%g, eta=%g, centrality=%g", ptMC, etaMC, centFT0M); + } + protonTracksTruthMatchedCorrected.push_back({ptMC, weight}); + } + } } else { hist.fill(HIST("QA/Proton/h_PtTruth_secondary"), track.pt()); } @@ -1085,14 +1278,51 @@ struct MeanPtFlucId { } } } - fillDifferentialMatrices(chargedTracks, ptBins.value); - fillDifferentialMatrices(pionTracks, ptBinsPi.value); - fillDifferentialMatrices(kaonTracks, ptBinsKa.value); - fillDifferentialMatrices(protonTracks, ptBinsPr.value); - fillAnalysis(nCh, q1Ch, q2Ch); - fillAnalysis(nPi, q1Pi, q2Pi); - fillAnalysis(nKa, q1Ka, q2Ka); - fillAnalysis(nPr, q1Pr, q2Pr); + fillMatrices(chargedTracks, cfgCutPtMin); + fillMatrices(pionTracks, cfgCutPiPtMin); + fillMatrices(kaonTracks, cfgCutKaPtMin); + fillMatrices(protonTracks, cfgCutPrPtMin); + + if (cfgApplyEfficiencyCorrection) { + fillCorrectedMatrices(chargedTracksCorrected, cfgCutPtMin); + fillCorrectedMatrices(pionTracksCorrected, cfgCutPiPtMin); + fillCorrectedMatrices(kaonTracksCorrected, cfgCutKaPtMin); + fillCorrectedMatrices(protonTracksCorrected, cfgCutPrPtMin); + if constexpr (RecoFlag) { + if (!cfgUseCombinedPurityEfficiencyWeight) { + fillCorrectedMatrices(chargedTracksTruthMatchedCorrected, cfgCutPtMin); + fillCorrectedMatrices(pionTracksTruthMatchedCorrected, cfgCutPiPtMin); + fillCorrectedMatrices(kaonTracksTruthMatchedCorrected, cfgCutKaPtMin); + fillCorrectedMatrices(protonTracksTruthMatchedCorrected, cfgCutPrPtMin); + } + } + } + } + + template + int countSimMultiplicity(M const& mcParticles) + { + int multiplicity = 0; + for (auto const& mcPart : mcParticles) { + if (!mcPart.isPhysicalPrimary() || std::abs(mcPart.eta()) >= cfgCutEta) { + continue; + } + auto* particle = pdg->GetParticle(mcPart.pdgCode()); + if (particle && particle->Charge() != 0.) { + ++multiplicity; + } + } + return multiplicity; + } + + bool passesNtpcEventCut() + { + if (!cfgNtpcEventCut) { + return true; + } + const float lower = cfgP2L * nSim * nSim + cfgP1L * nSim - cfgP0L; + const float upper = cfgP2U * nSim * nSim + cfgP1U * nSim + cfgP0U; + return nTPC >= lower && nTPC <= upper; } template @@ -1102,10 +1332,9 @@ struct MeanPtFlucId { pionTracksGen.clear(); kaonTracksGen.clear(); protonTracksGen.clear(); + nSim = 0; - int nChSim = 0, nPiSim = 0, nKaSim = 0, nPrSim = 0; float pt = 0., eta = 0, phi = 0., rap = 0.; - double q1Ch = 0., q2Ch = 0., q1Pi = 0., q2Pi = 0., q1Ka = 0., q2Ka = 0., q1Pr = 0., q2Pr = 0.; for (auto const& mcPart : mcParticles) { if (!mcPart.isPhysicalPrimary()) { continue; @@ -1113,11 +1342,13 @@ struct MeanPtFlucId { auto* particle = pdg->GetParticle(mcPart.pdgCode()); - if (!particle) + if (!particle) { continue; + } - if (particle->Charge() == 0) + if (particle->Charge() == 0) { continue; + } pt = mcPart.pt(); eta = mcPart.eta(); @@ -1127,15 +1358,10 @@ struct MeanPtFlucId { nSim++; } - auto pid = mcPart.pdgCode(); - if (std::abs(pid) != kElectron && std::abs(pid) != kMuonMinus && std::abs(pid) != kPiPlus && std::abs(pid) != kKPlus && std::abs(pid) != kProton) { - continue; - } - if (pt >= cfgCutPtMin && pt < cfgCutPtMax && std::abs(eta) < cfgCutEta) { - moments(pt, nChSim, q1Ch, q2Ch); - chargedTracksGen.push_back({pt, eta}); + chargedTracksGen.push_back({pt}); hist.fill(HIST("Gen/Charged/h_PtTruth"), pt); + hist.fill(HIST("Gen/Charged/h_PtTruth_e"), pt); hist.fill(HIST("Gen/Charged/h_EtaTruth"), eta); hist.fill(HIST("Gen/Charged/h_PhiTruth"), phi); hist.fill(HIST("Gen/Charged/h2_Pt_EtaTruth"), eta, pt); @@ -1144,32 +1370,26 @@ struct MeanPtFlucId { hist.fill(HIST("Gen/Charged/h3_Pt_Eta_PhiTruth"), pt, eta, phi); rap = mcPart.y(); + auto pid = mcPart.pdgCode(); if (std::abs(pid) == kPiPlus && pt >= cfgCutPiPtMin) { - moments(pt, nPiSim, q1Pi, q2Pi); - pionTracksGen.push_back({pt, eta}); - fillPtMCHist(true, pt, eta, rap, phi, pid, kPiMinus, kPiMinus); + pionTracksGen.push_back({pt}); + fillPtMCHist(true, pt, eta, rap, phi, pid, kPiPlus, kPiMinus); } if (std::abs(pid) == kKPlus && pt >= cfgCutKaPtMin) { - moments(pt, nKaSim, q1Ka, q2Ka); - kaonTracksGen.push_back({pt, eta}); - fillPtMCHist(true, pt, eta, rap, phi, pid, kKMinus, kKMinus); + kaonTracksGen.push_back({pt}); + fillPtMCHist(true, pt, eta, rap, phi, pid, kKPlus, kKMinus); } if (std::abs(pid) == kProton && pt >= cfgCutPrPtMin) { - moments(pt, nPrSim, q1Pr, q2Pr); - protonTracksGen.push_back({pt, eta}); - fillPtMCHist(true, pt, eta, rap, phi, pid, kProtonBar, kProtonBar); + protonTracksGen.push_back({pt}); + fillPtMCHist(true, pt, eta, rap, phi, pid, kProton, kProtonBar); } } } - fillDifferentialMatrices(chargedTracksGen, ptBins.value); - fillDifferentialMatrices(pionTracksGen, ptBinsPi.value); - fillDifferentialMatrices(kaonTracksGen, ptBinsKa.value); - fillDifferentialMatrices(protonTracksGen, ptBinsPr.value); - fillAnalysis(nChSim, q1Ch, q2Ch); - fillAnalysis(nPiSim, q1Pi, q2Pi); - fillAnalysis(nKaSim, q1Ka, q2Ka); - fillAnalysis(nPrSim, q1Pr, q2Pr); + fillMatrices(chargedTracksGen, cfgCutPtMin); + fillMatrices(pionTracksGen, cfgCutPiPtMin); + fillMatrices(kaonTracksGen, cfgCutKaPtMin); + fillMatrices(protonTracksGen, cfgCutPrPtMin); hist.fill(HIST("Gen/h_Counts"), 2); hist.fill(HIST("Gen/h_VtxZ"), mcCol.posZ()); @@ -1200,6 +1420,14 @@ struct MeanPtFlucId { hist.fill(HIST("Gen/h_collisions_info"), kPassSelCol); hist.fill(HIST("Gen/h2_collision_posZ"), mcCol.posZ(), cols.begin().posZ()); + // The generated and reconstructed moments must be evaluated for the same + // event ensemble. Determine nSim and apply the correlated multiplicity cut + // before filling either level. + nSim = countSimMultiplicity(mcParts); + if (!passesNtpcEventCut()) { + return; + } + auto sTracks = tracks.sliceBy(perCollision, cols.begin().globalIndex()); fillGenHistos(mcCol, mcParts); fillRecoHistos(cols.begin(), sTracks); @@ -1209,8 +1437,8 @@ struct MeanPtFlucId { aod::pidTOFFullPi, aod::pidTPCFullPi, aod::pidTOFFullPr, aod::pidTPCFullPr, aod::pidTOFFullKa, aod::pidTPCFullKa, aod::pidTOFFullEl, aod::pidTPCFullEl, aod::pidTOFbeta, aod::pidTOFmass>; - using MyRun3Collisions = soa::Join; - using MyRun3MCCollisions = soa::Join; + using MyRun3Collisions = soa::Join; + using MyRun3MCCollisions = soa::Join; using MyMCTracks = soa::Join Date: Fri, 28 Aug 2026 11:05:02 +0530 Subject: [PATCH 2/3] Fix include directive formatting in meanPtFlucId.cxx --- PWGCF/EbyEFluctuations/Tasks/meanPtFlucId.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/PWGCF/EbyEFluctuations/Tasks/meanPtFlucId.cxx b/PWGCF/EbyEFluctuations/Tasks/meanPtFlucId.cxx index b704037695f..6e258f1575d 100644 --- a/PWGCF/EbyEFluctuations/Tasks/meanPtFlucId.cxx +++ b/PWGCF/EbyEFluctuations/Tasks/meanPtFlucId.cxx @@ -24,7 +24,7 @@ #include "Common/DataModel/PIDResponseTPC.h" #include "Common/DataModel/TrackSelectionTables.h" -#include "CCDB/BasicCCDBManager.h" +#include #include #include #include From 3a0a030c446a72f518478631f061e068cb023d37 Mon Sep 17 00:00:00 2001 From: Tanu Gahlaut <154991749+TGahlaut1@users.noreply.github.com> Date: Fri, 28 Aug 2026 18:15:10 +0530 Subject: [PATCH 3/3] Fix build error --- PWGCF/EbyEFluctuations/Tasks/meanPtFlucId.cxx | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/PWGCF/EbyEFluctuations/Tasks/meanPtFlucId.cxx b/PWGCF/EbyEFluctuations/Tasks/meanPtFlucId.cxx index 6e258f1575d..7f4432ae05e 100644 --- a/PWGCF/EbyEFluctuations/Tasks/meanPtFlucId.cxx +++ b/PWGCF/EbyEFluctuations/Tasks/meanPtFlucId.cxx @@ -118,7 +118,7 @@ struct MeanPtFlucId { Configurable> etaBins{"etaBins", {-0.8, -0.7, -0.6, -0.5, -0.4, -0.3, -0.2, -0.1, 0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8}, "#eta bins"}; Service pdg{}; - Service ccdb; + Service ccdb{}; HistogramRegistry hist{"hist", {}, OutputObjHandlingPolicy::AnalysisObject}; @@ -252,7 +252,7 @@ struct MeanPtFlucId { if (cfgApplyEfficiencyCorrection) { ccdb->setURL(cfgEfficiencyCCDBUrl.value); ccdb->setCaching(true); - TList* efficiencyList = ccdb->getForTimeStamp(cfgEfficiencyCCDBPath.value, cfgEfficiencyTimestamp.value); + auto* efficiencyList = ccdb->getForTimeStamp(cfgEfficiencyCCDBPath.value, cfgEfficiencyTimestamp.value); if (!efficiencyList) { LOGF(fatal, "Cannot load ccdb_object containing the efficiency maps"); } @@ -1460,7 +1460,7 @@ struct MeanPtFlucId { PROCESS_SWITCH(MeanPtFlucId, processMCRecoSimRun3, "process MC Reconstructed & Truth Run-3", true); }; -WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +WorkflowSpec defineDataProcessing(ConfigContext const& context) { - return WorkflowSpec{adaptAnalysisTask(cfgc)}; + return WorkflowSpec{adaptAnalysisTask(context)}; }