From 82628a3c68714b8f4a31aad114905e4fd5fa7b34 Mon Sep 17 00:00:00 2001 From: Pravata Panigrahi <158029188+pravatp98@users.noreply.github.com> Date: Fri, 14 Aug 2026 19:32:42 +0530 Subject: [PATCH 01/13] Updated for efficiency correction Updated for efficiency correction and sub-sampling --- .../EbyEFluctuations/Tasks/nchCumulantsId.cxx | 519 ++++++++++++++++-- 1 file changed, 464 insertions(+), 55 deletions(-) diff --git a/PWGCF/EbyEFluctuations/Tasks/nchCumulantsId.cxx b/PWGCF/EbyEFluctuations/Tasks/nchCumulantsId.cxx index 1fd39ebf2fc..3c23ef7558f 100644 --- a/PWGCF/EbyEFluctuations/Tasks/nchCumulantsId.cxx +++ b/PWGCF/EbyEFluctuations/Tasks/nchCumulantsId.cxx @@ -23,6 +23,7 @@ #include #include +#include #include #include #include @@ -36,10 +37,15 @@ #include #include +#include +#include +#include +#include #include #include #include +#include #include #include #include @@ -79,6 +85,49 @@ enum ChargeEnum { kNeg = 1 }; +enum class FCPrefixEnum { + Pr = 0, + APr = 1, + PiPos = 2, + PiNeg = 3, + KaPos = 4, + KaNeg = 5, + Pos = 6, + Neg = 7, + NetPi = 8, + NetKa = 9, + NetPr = 10, + NetCh = 11 +}; + +static constexpr std::string_view FCRecoDir[] = { + "Reco/Pr/", + "Reco/APr/", + "Reco/PiPos/", + "Reco/PiNeg/", + "Reco/KaPos/", + "Reco/KaNeg/", + "Reco/Pos/", + "Reco/Neg/", + "Reco/NetPi/", + "Reco/NetKa/", + "Reco/NetPr/", + "Reco/NetCh/"}; + +static constexpr std::string_view FCGenDir[] = { + "Gen/Pr/", + "Gen/APr/", + "Gen/PiPos/", + "Gen/PiNeg/", + "Gen/KaPos/", + "Gen/KaNeg/", + "Gen/Pos/", + "Gen/Neg/", + "Gen/NetPi/", + "Gen/NetKa/", + "Gen/NetPr/", + "Gen/NetCh/"}; + static constexpr std::string_view PidDire[] = { "Ch/", "Pi/", @@ -105,10 +154,10 @@ std::string getModifiedStr(const std::string& myString) struct NchCumulantsId { HistogramRegistry hist{"hist", {}, OutputObjHandlingPolicy::AnalysisObject}; - HistogramRegistry recoTracks{"recoTracks", {}, OutputObjHandlingPolicy::AnalysisObject}; HistogramRegistry genAnalysis{"genAnalysis", {}, OutputObjHandlingPolicy::AnalysisObject}; HistogramRegistry recoAnalysis{"recoAnalysis", {}, OutputObjHandlingPolicy::AnalysisObject}; HistogramRegistry purityAnalysis{"purityAnalysis", {}, OutputObjHandlingPolicy::AnalysisObject}; + HistogramRegistry registry{"registry", {}, OutputObjHandlingPolicy::AnalysisObject}; // PDG data base Service pdgDB; @@ -130,14 +179,41 @@ struct NchCumulantsId { ConfigurableAxis axisAKaCh{"axisAKaCh", {3010, -0.5, 300.5}, "AKa_charge"}; ConfigurableAxis axisPiCh{"axisPiCh", {3010, -0.5, 300.5}, "Pion_Positive"}; ConfigurableAxis axisAPiCh{"axisAPiCh", {3010, -0.5, 300.5}, "Pion_Negative"}; + // QA axes as ConfigurableAxis + ConfigurableAxis axisEvents{"axisEvents", {1, 0, 1}, "Counts"}; + ConfigurableAxis axisEta{"axisEta", {100, -1., +1.}, "#eta"}; + ConfigurableAxis axisRapidity{"axisRapidity", {200, -5, 5}, "Rapidity (y)"}; + ConfigurableAxis axisPt{"axisPt", {100, 0., 5.}, "p_{T} (GeV/c)"}; + ConfigurableAxis axisP{"axisP", {100, 0., 5.}, "p (GeV/c)"}; + ConfigurableAxis axisTPCInnerParam{"axisTPCInnerParam", {100, 0, 3}, "P_innerParam_Gev"}; + ConfigurableAxis axisdEdx{"axisdEdx", {100, 20, 500}, "#frac{dE}{dx}"}; + ConfigurableAxis axisVtxZ{"axisVtxZ", {80, -20., 20.}, "V_{Z} (cm)"}; + ConfigurableAxis axisDCAz{"axisDCAz", {200, -3., 3.}, "DCA_{Z} (cm)"}; + ConfigurableAxis axisDCAxy{"axisDCAxy", {200, -3., 3.}, "DCA_{XY} (cm)"}; + ConfigurableAxis axisMultFT0{"axisMultFT0", {150, 0, 1500}, "MultFT0"}; + ConfigurableAxis axisCent{"axisCent", {103, -1., 102.}, "FT0C(%)"}; + ConfigurableAxis axisPhi{"axisPhi", {80, -1, 7}, "phi"}; + ConfigurableAxis axisTOFBeta{"axisTOFBeta", {40, -2.0, 2.0}, "tofBeta"}; + ConfigurableAxis axisTPCSignal{"axisTPCSignal", {100, -1, 1000}, "tpcSignal"}; + ConfigurableAxis axisTPCNSigma{"axisTPCNSigma", {200, -10.0, 10.0}, "n#sigma_{TPC}"}; + ConfigurableAxis axisTOFNSigma{"axisTOFNSigma", {200, -10.0, 10.0}, "n#sigma_{TOF}"}; + ConfigurableAxis axisTOFExpMom{"axisTOFExpMom", {200, 0.0f, 10.0f}, "#it{p}_{tofExpMom} (GeV/#it{c})"}; - Configurable checkCollPosZMc{"checkCollPosZMc", false, "checkCollPosZMc"}; - Configurable flagUnusedVariableError{"flagUnusedVariableError", false, "flagUnusedVariableError"}; - Configurable cfgDoRejectionForId{"cfgDoRejectionForId", false, "Apply rejection cut before PID selection (selTrackForId)"}; - - Configurable cfgEvSel01doNoSameBunchPileup{"cfgEvSel01doNoSameBunchPileup", true, "apply kNoSameBunchPileup"}; - Configurable cfgEvSel02doIsGoodZvtxFT0vsPV{"cfgEvSel02doIsGoodZvtxFT0vsPV", true, "apply kIsGoodZvtxFT0vsPV"}; - Configurable cfgEvSel03doIsGoodITSLayersAll{"cfgEvSel03doIsGoodITSLayersAll", true, "apply kIsGoodITSLayersAll"}; + struct : ConfigurableGroup { + Configurable checkCollPosZMc{"checkCollPosZMc", false, "checkCollPosZMc"}; + Configurable flagUnusedVariableError{"flagUnusedVariableError", false, "flagUnusedVariableError"}; + Configurable cfgDoRejectionForId{"cfgDoRejectionForId", false, "Apply rejection cut before PID selection (selTrackForId)"}; + Configurable fillSparseForReco{"fillSparseForReco", false, "Fill sparse for reconstructed tracks"}; + Configurable fillSparseForGen{"fillSparseForGen", false, "Fill sparse for generated tracks"}; + Configurable fillSparseForPurity{"fillSparseForPurity", false, "Fill sparse for purity tracks"}; + Configurable cfgEvSel01doNoSameBunchPileup{"cfgEvSel01doNoSameBunchPileup", true, "apply kNoSameBunchPileup"}; + Configurable cfgEvSel02doIsGoodZvtxFT0vsPV{"cfgEvSel02doIsGoodZvtxFT0vsPV", true, "apply kIsGoodZvtxFT0vsPV"}; + Configurable cfgEvSel03doIsGoodITSLayersAll{"cfgEvSel03doIsGoodITSLayersAll", true, "apply kIsGoodITSLayersAll"}; + Configurable cfgDoSubsampling{"cfgDoSubsampling", true, "do subsampling for error estimation"}; + ConfigurableAxis subSampleAxis{"subSampleAxis", {10, 0., 10.}, "Subsample"}; + TRandom3* fRandom = new TRandom3(0); // Random number generator for subsampling + int currentSubsample = 0; + } cfgEventSelection; // Configurables for particle Identification Configurable cfgId01CheckVetoCut{"cfgId01CheckVetoCut", true, "cfgId01CheckVetoCut"}; @@ -213,6 +289,263 @@ struct NchCumulantsId { TH2F* hPtEtaForBinSearch = nullptr; std::vector> hPtEtaForEffCorrection{kDe + 1, std::array{}}; + // AxisSpec subSampleAxis = {static_cast(cfgNSubsamples), 0., static_cast(cfgNSubsamples), "Subsample"}; + + struct EffPowerSums { + double q1 = 0.; + double q2 = 0.; + double q3 = 0.; + double q4 = 0.; + }; + + inline void fillEffPower(EffPowerSums& p, float weight) + { + if (weight <= 0.f) { + return; + } + + p.q1 += weight; + p.q2 += weight * weight; + p.q3 += weight * weight * weight; + p.q4 += weight * weight * weight * weight; + } + + // NEW TEMPLATE FUNCTION FOR RECO PROFILES for individual species + template + void addFCRecoProfiles() + { + constexpr std::string_view Dir = FCRecoDir[static_cast(Prefix)]; + + registry.add(std::string(Dir) + "Q1", "", kTProfile, {axisCent}); + registry.add(std::string(Dir) + "Q1Sq", "", kTProfile, {axisCent}); + registry.add(std::string(Dir) + "Q2", "", kTProfile, {axisCent}); + registry.add(std::string(Dir) + "Q1Cube", "", kTProfile, {axisCent}); + registry.add(std::string(Dir) + "Q1Q2", "", kTProfile, {axisCent}); + registry.add(std::string(Dir) + "Q3", "", kTProfile, {axisCent}); + registry.add(std::string(Dir) + "Q1Pow4", "", kTProfile, {axisCent}); + registry.add(std::string(Dir) + "Q1SqQ2", "", kTProfile, {axisCent}); + registry.add(std::string(Dir) + "Q2Sq", "", kTProfile, {axisCent}); + registry.add(std::string(Dir) + "Q1Q3", "", kTProfile, {axisCent}); + registry.add(std::string(Dir) + "Q4", "", kTProfile, {axisCent}); + + if (cfgEventSelection.cfgDoSubsampling) { + registry.add(std::string(Dir) + "Q1_subsample", "", kTProfile2D, {axisCent, cfgEventSelection.subSampleAxis}); + registry.add(std::string(Dir) + "Q1Sq_subsample", "", kTProfile2D, {axisCent, cfgEventSelection.subSampleAxis}); + registry.add(std::string(Dir) + "Q2_subsample", "", kTProfile2D, {axisCent, cfgEventSelection.subSampleAxis}); + registry.add(std::string(Dir) + "Q1Cube_subsample", "", kTProfile2D, {axisCent, cfgEventSelection.subSampleAxis}); + registry.add(std::string(Dir) + "Q1Q2_subsample", "", kTProfile2D, {axisCent, cfgEventSelection.subSampleAxis}); + registry.add(std::string(Dir) + "Q3_subsample", "", kTProfile2D, {axisCent, cfgEventSelection.subSampleAxis}); + registry.add(std::string(Dir) + "Q1Pow4_subsample", "", kTProfile2D, {axisCent, cfgEventSelection.subSampleAxis}); + registry.add(std::string(Dir) + "Q1SqQ2_subsample", "", kTProfile2D, {axisCent, cfgEventSelection.subSampleAxis}); + registry.add(std::string(Dir) + "Q2Sq_subsample", "", kTProfile2D, {axisCent, cfgEventSelection.subSampleAxis}); + registry.add(std::string(Dir) + "Q1Q3_subsample", "", kTProfile2D, {axisCent, cfgEventSelection.subSampleAxis}); + registry.add(std::string(Dir) + "Q4_subsample", "", kTProfile2D, {axisCent, cfgEventSelection.subSampleAxis}); + } + } + + template + void addNetQVectorProfileHistograms() + { + constexpr std::string_view Dir = FCRecoDir[static_cast(Prefix)]; + + registry.add(std::string(Dir) + "Q_net_1", "Net Q1", kTProfile, {axisCent}); + registry.add(std::string(Dir) + "Q_net_1Sq", "Net Q1²", kTProfile, {axisCent}); + registry.add(std::string(Dir) + "Q_net_2", "Net Q2", kTProfile, {axisCent}); + registry.add(std::string(Dir) + "Q_net_1Cube", "Net Q1³", kTProfile, {axisCent}); + registry.add(std::string(Dir) + "Q_net_1Q_net_2", "Net Q1·Q2", kTProfile, {axisCent}); + registry.add(std::string(Dir) + "Q_net_3", "Net Q3", kTProfile, {axisCent}); + registry.add(std::string(Dir) + "Q_net_1Pow4", "Net Q1⁴", kTProfile, {axisCent}); + registry.add(std::string(Dir) + "Q_net_1SqQ_net_2", "Net Q1²·Q2", kTProfile, {axisCent}); + registry.add(std::string(Dir) + "Q_net_2Sq", "Net Q2²", kTProfile, {axisCent}); + registry.add(std::string(Dir) + "Q_net_1Q_net_3", "Net Q1·Q3", kTProfile, {axisCent}); + registry.add(std::string(Dir) + "Q_net_4", "Net Q4", kTProfile, {axisCent}); + + // joint pos-neg correction profiles, they should be stored here. F11 needed for F3 and F21 and F12 for F4 correction + registry.add(std::string(Dir) + "JointF11", "Joint ", kTProfile, {axisCent}); + registry.add(std::string(Dir) + "JointF12", "Joint ", kTProfile, {axisCent}); + registry.add(std::string(Dir) + "JointF21", "Joint ", kTProfile, {axisCent}); + + // Factorial moments + if (cfgEventSelection.cfgDoSubsampling) { + registry.add(std::string(Dir) + "Q_net_1_subsample", "Net Q1 Subsample", kTProfile2D, {axisCent, cfgEventSelection.subSampleAxis}); + registry.add(std::string(Dir) + "Q_net_1Sq_subsample", "Net Q1² Subsample", kTProfile2D, {axisCent, cfgEventSelection.subSampleAxis}); + registry.add(std::string(Dir) + "Q_net_2_subsample", "Net Q2 Subsample", kTProfile2D, {axisCent, cfgEventSelection.subSampleAxis}); + registry.add(std::string(Dir) + "Q_net_1Cube_subsample", "Net Q1³ Subsample", kTProfile2D, {axisCent, cfgEventSelection.subSampleAxis}); + registry.add(std::string(Dir) + "Q_net_1Q_net_2_subsample", "Net Q1·Q2 Subsample", kTProfile2D, {axisCent, cfgEventSelection.subSampleAxis}); + registry.add(std::string(Dir) + "Q_net_3_subsample", "Net Q3 Subsample", kTProfile2D, {axisCent, cfgEventSelection.subSampleAxis}); + registry.add(std::string(Dir) + "Q_net_1Pow4_subsample", "Net Q1⁴ Subsample", kTProfile2D, {axisCent, cfgEventSelection.subSampleAxis}); + registry.add(std::string(Dir) + "Q_net_1SqQ_net_2_subsample", "Net Q1²·Q2 Subsample", kTProfile2D, {axisCent, cfgEventSelection.subSampleAxis}); + registry.add(std::string(Dir) + "Q_net_2Sq_subsample", "Net Q2² Subsample", kTProfile2D, {axisCent, cfgEventSelection.subSampleAxis}); + registry.add(std::string(Dir) + "Q_net_1Q_net_3_subsample", "Net Q1·Q3 Subsample", kTProfile2D, {axisCent, cfgEventSelection.subSampleAxis}); + registry.add(std::string(Dir) + "Q_net_4_subsample", "Net Q4 Subsample", kTProfile2D, {axisCent, cfgEventSelection.subSampleAxis}); + + registry.add(std::string(Dir) + "JointF11_subsample", "Joint Subsample", kTProfile2D, {axisCent, cfgEventSelection.subSampleAxis}); + registry.add(std::string(Dir) + "JointF12_subsample", "Joint Subsample", kTProfile2D, {axisCent, cfgEventSelection.subSampleAxis}); + registry.add(std::string(Dir) + "JointF21_subsample", "Joint Subsample", kTProfile2D, {axisCent, cfgEventSelection.subSampleAxis}); + } + } + + // NEW TEMPLATE FUNCTION FOR GEN PROFILES + template + void addFCGenProfiles() + { + constexpr std::string_view Dir = FCGenDir[static_cast(Prefix)]; + + registry.add(std::string(Dir) + "F1", "Gen FactorialMoment1", kTProfile, {axisCent}); + registry.add(std::string(Dir) + "F2", "Gen FactorialMoment2", kTProfile, {axisCent}); + registry.add(std::string(Dir) + "F3", "Gen FactorialMoment3", kTProfile, {axisCent}); + registry.add(std::string(Dir) + "F4", "Gen FactorialMoment4", kTProfile, {axisCent}); + + if (cfgEventSelection.cfgDoSubsampling) { + registry.add(std::string(Dir) + "F1_subsample", "Gen FactorialMoment1 Subsample", kTProfile2D, {axisCent, cfgEventSelection.subSampleAxis}); + registry.add(std::string(Dir) + "F2_subsample", "Gen FactorialMoment2 Subsample", kTProfile2D, {axisCent, cfgEventSelection.subSampleAxis}); + registry.add(std::string(Dir) + "F3_subsample", "Gen FactorialMoment3 Subsample", kTProfile2D, {axisCent, cfgEventSelection.subSampleAxis}); + registry.add(std::string(Dir) + "F4_subsample", "Gen FactorialMoment4 Subsample", kTProfile2D, {axisCent, cfgEventSelection.subSampleAxis}); + } + } + + template + void fillFCBasis(const EffPowerSums& p, float cent, Registry& registry) + { + auto base = HIST(FCRecoDir[static_cast(Prefix)]); + + // Fill Factorial Moments for Reco directly to ensure Reco stores the same type of data as Gen + // Compute Efficiency Corrected Factorial Moments from power sums + // These formulas come from elementary symmetric polynomials and their relation to power sums + registry.fill(base + HIST("Q1"), cent, p.q1); + registry.fill(base + HIST("Q1Sq"), cent, p.q1 * p.q1); + registry.fill(base + HIST("Q2"), cent, p.q2); + registry.fill(base + HIST("Q1Cube"), cent, p.q1 * p.q1 * p.q1); + registry.fill(base + HIST("Q1Q2"), cent, p.q1 * p.q2); + registry.fill(base + HIST("Q3"), cent, p.q3); + registry.fill(base + HIST("Q1Pow4"), cent, p.q1 * p.q1 * p.q1 * p.q1); + registry.fill(base + HIST("Q1SqQ2"), cent, p.q1 * p.q1 * p.q2); + registry.fill(base + HIST("Q2Sq"), cent, p.q2 * p.q2); + registry.fill(base + HIST("Q1Q3"), cent, p.q1 * p.q3); + registry.fill(base + HIST("Q4"), cent, p.q4); + if (cfgEventSelection.cfgDoSubsampling) { + registry.fill(base + HIST("Q1_subsample"), cent, cfgEventSelection.currentSubsample, p.q1); + registry.fill(base + HIST("Q1Sq_subsample"), cent, cfgEventSelection.currentSubsample, p.q1 * p.q1); + registry.fill(base + HIST("Q2_subsample"), cent, cfgEventSelection.currentSubsample, p.q2); + registry.fill(base + HIST("Q1Cube_subsample"), cent, cfgEventSelection.currentSubsample, p.q1 * p.q1 * p.q1); + registry.fill(base + HIST("Q1Q2_subsample"), cent, cfgEventSelection.currentSubsample, p.q1 * p.q2); + registry.fill(base + HIST("Q3_subsample"), cent, cfgEventSelection.currentSubsample, p.q3); + registry.fill(base + HIST("Q1Pow4_subsample"), cent, cfgEventSelection.currentSubsample, p.q1 * p.q1 * p.q1 * p.q1); + registry.fill(base + HIST("Q1SqQ2_subsample"), cent, cfgEventSelection.currentSubsample, p.q1 * p.q1 * p.q2); + registry.fill(base + HIST("Q2Sq_subsample"), cent, cfgEventSelection.currentSubsample, p.q2 * p.q2); + registry.fill(base + HIST("Q1Q3_subsample"), cent, cfgEventSelection.currentSubsample, p.q1 * p.q3); + registry.fill(base + HIST("Q4_subsample"), cent, cfgEventSelection.currentSubsample, p.q4); + } + } + + // fill function for reco net species + template + void fillNetQVectorProfileHistograms(const EffPowerSums& posPow, const EffPowerSums& negPow, float cent, Registry& registry) + { + // Construct base directory path + auto base = HIST(FCRecoDir[static_cast(Prefix)]); + + // ─── Compute NET q-vectors (sign alternation by power) ─── + // Odd powers: alternate sign (Q_pos - Q_neg) + // Even powers: same sign (Q_pos + Q_neg) + + double qNet1 = posPow.q1 - negPow.q1; // Odd + double qNet2 = posPow.q2 + negPow.q2; // Even + double qNet3 = posPow.q3 - negPow.q3; // Odd + double qNet4 = posPow.q4 + negPow.q4; // Even + + // Fill joint pos-neg correction profiles + // capture inline from pos and neg + double jointF11 = posPow.q1 * negPow.q1; + double cf2pos = posPow.q1 * posPow.q1 - posPow.q2; + double cf2neg = negPow.q1 * negPow.q1 - negPow.q2; + double jointF12 = posPow.q1 * cf2neg; + double jointF21 = negPow.q1 * cf2pos; + // ─── Fill raw q-vectors (for cross-verification) + registry.fill(base + HIST("Q_net_1"), cent, qNet1); + registry.fill(base + HIST("Q_net_1Sq"), cent, qNet1 * qNet1); + registry.fill(base + HIST("Q_net_2"), cent, qNet2); + registry.fill(base + HIST("Q_net_1Cube"), cent, qNet1 * qNet1 * qNet1); + registry.fill(base + HIST("Q_net_1Q_net_2"), cent, qNet1 * qNet2); + registry.fill(base + HIST("Q_net_3"), cent, qNet3); + registry.fill(base + HIST("Q_net_1Pow4"), cent, qNet1 * qNet1 * qNet1 * qNet1); + registry.fill(base + HIST("Q_net_1SqQ_net_2"), cent, qNet1 * qNet1 * qNet2); + registry.fill(base + HIST("Q_net_2Sq"), cent, qNet2 * qNet2); + registry.fill(base + HIST("Q_net_1Q_net_3"), cent, qNet1 * qNet3); + registry.fill(base + HIST("Q_net_4"), cent, qNet4); + + // ─── Fill factorial moments ─── + registry.fill(base + HIST("JointF11"), cent, jointF11); + registry.fill(base + HIST("JointF12"), cent, jointF12); + registry.fill(base + HIST("JointF21"), cent, jointF21); + + if (cfgEventSelection.cfgDoSubsampling) { + registry.fill(base + HIST("Q_net_1_subsample"), cent, cfgEventSelection.currentSubsample, qNet1); + registry.fill(base + HIST("Q_net_1Sq_subsample"), cent, cfgEventSelection.currentSubsample, qNet1 * qNet1); + registry.fill(base + HIST("Q_net_2_subsample"), cent, cfgEventSelection.currentSubsample, qNet2); + registry.fill(base + HIST("Q_net_1Cube_subsample"), cent, cfgEventSelection.currentSubsample, qNet1 * qNet1 * qNet1); + registry.fill(base + HIST("Q_net_1Q_net_2_subsample"), cent, cfgEventSelection.currentSubsample, qNet1 * qNet2); + registry.fill(base + HIST("Q_net_3_subsample"), cent, cfgEventSelection.currentSubsample, qNet3); + registry.fill(base + HIST("Q_net_1Pow4_subsample"), cent, cfgEventSelection.currentSubsample, qNet1 * qNet1 * qNet1 * qNet1); + registry.fill(base + HIST("Q_net_1SqQ_net_2_subsample"), cent, cfgEventSelection.currentSubsample, qNet1 * qNet1 * qNet2); + registry.fill(base + HIST("Q_net_2Sq_subsample"), cent, cfgEventSelection.currentSubsample, qNet2 * qNet2); + registry.fill(base + HIST("Q_net_1Q_net_3_subsample"), cent, cfgEventSelection.currentSubsample, qNet1 * qNet3); + registry.fill(base + HIST("Q_net_4_subsample"), cent, cfgEventSelection.currentSubsample, qNet4); + + registry.fill(base + HIST("JointF11_subsample"), cent, cfgEventSelection.currentSubsample, jointF11); + registry.fill(base + HIST("JointF12_subsample"), cent, cfgEventSelection.currentSubsample, jointF12); + registry.fill(base + HIST("JointF21_subsample"), cent, cfgEventSelection.currentSubsample, jointF21); + } + } + + template + void fillGenFactorialMoments(int n, float cent, Registry& registry) + { + auto base = HIST(FCGenDir[static_cast(Prefix)]); + + double f1 = n; + double f2 = n * (n - 1.); + double f3 = n * (n - 1.) * (n - 2.); + double f4 = n * (n - 1.) * (n - 2.) * (n - 3.); + + registry.fill(base + HIST("F1"), cent, f1); + registry.fill(base + HIST("F2"), cent, f2); + registry.fill(base + HIST("F3"), cent, f3); + registry.fill(base + HIST("F4"), cent, f4); + + if (cfgEventSelection.cfgDoSubsampling) { + registry.fill(base + HIST("F1_subsample"), cent, cfgEventSelection.currentSubsample, f1); + registry.fill(base + HIST("F2_subsample"), cent, cfgEventSelection.currentSubsample, f2); + registry.fill(base + HIST("F3_subsample"), cent, cfgEventSelection.currentSubsample, f3); + registry.fill(base + HIST("F4_subsample"), cent, cfgEventSelection.currentSubsample, f4); + } + } + // ─── GEN level (net speciestruth) ─── + template + void fillGenNetFactorialMoments(float nPos, float nNeg, float cent, Registry& registry) + { + auto base = HIST(FCGenDir[static_cast(Prefix)]); + + double nNet = nPos - nNeg; + + double f1 = nNet; + double f2 = nNet * (nNet - 1.0); + double f3 = nNet * (nNet - 1.0) * (nNet - 2.0); + double f4 = nNet * (nNet - 1.0) * (nNet - 2.0) * (nNet - 3.0); + + registry.fill(base + HIST("F1"), cent, f1); + registry.fill(base + HIST("F2"), cent, f2); + registry.fill(base + HIST("F3"), cent, f3); + registry.fill(base + HIST("F4"), cent, f4); + + if (cfgEventSelection.cfgDoSubsampling) { + registry.fill(base + HIST("F1_subsample"), cent, cfgEventSelection.currentSubsample, f1); + registry.fill(base + HIST("F2_subsample"), cent, cfgEventSelection.currentSubsample, f2); + registry.fill(base + HIST("F3_subsample"), cent, cfgEventSelection.currentSubsample, f3); + registry.fill(base + HIST("F4_subsample"), cent, cfgEventSelection.currentSubsample, f4); + } + } + void init(InitContext const&) { auto& mgr = o2::ccdb::BasicCCDBManager::instance(); @@ -244,27 +577,6 @@ struct NchCumulantsId { mgr.setURL("http://alice-ccdb.cern.ch"); // RESET the URL otherwise the other process functions which contains ccdb lookups will fail - // QA check axes - const AxisSpec axisEvents{1, 0, 1, "Counts"}; - const AxisSpec axisEta{100, -1., +1., "#eta"}; - const AxisSpec axisRapidity{200, -5, 5, "Rapidity (y)"}; - const AxisSpec axisPt{100, 0., 5., "p_{T} (GeV/c)"}; - const AxisSpec axisP{100, 0., 5., "p (GeV/c)"}; - const AxisSpec axisTPCInnerParam{100, 0, 3, "P_innerParam_Gev"}; - const AxisSpec axisdEdx(100, 20, 500, {"#frac{dE}{dx}"}); - const AxisSpec axisVtxZ{80, -20., 20., "V_{Z} (cm)"}; - const AxisSpec axisDCAz{200, -3., 3., "DCA_{Z} (cm)"}; - const AxisSpec axisDCAxy{200, -3., 3., "DCA_{XY} (cm)"}; - const AxisSpec axisMultFT0(150, 0, 1500, "MultFT0"); - const AxisSpec axisCent(103, -1., 102., "FT0C(%)"); - const AxisSpec axisPhi(80, -1, 7, "phi"); - - const AxisSpec axisTOFBeta = {40, -2.0, 2.0, "tofBeta"}; - const AxisSpec axisTPCSignal = {100, -1, 1000, "tpcSignal"}; - const AxisSpec axisTPCNSigma = {200, -10.0, 10.0, "n#sigma_{TPC}"}; - const AxisSpec axisTOFNSigma = {200, -10.0, 10.0, "n#sigma_{TOF}"}; - const AxisSpec axisTOFExpMom = {200, 0.0f, 10.0f, "#it{p}_{tofExpMom} (GeV/#it{c})"}; - const AxisSpec axisIdTag = {32, -0.5f, 31.5f, "idTag"}; const AxisSpec axisMcTag = {32, -0.5f, 31.5f, "mcTag"}; @@ -290,7 +602,6 @@ struct NchCumulantsId { HistogramConfigSpec histTpcNSigmaTofNSigma({HistType::kTH2F, {axisTPCNSigma, axisTOFNSigma}}); HistogramConfigSpec histPtMc({HistType::kTH1F, {axisPt}}); - // Register histograms for PID validation // Register TPC spares per species hist.add("PIDValidation/tpcSparse_Pi", "p vs tpcNSigmaPi vs idTag vs mcTag (Pion)", histTPCPIDSparse); @@ -473,6 +784,36 @@ struct NchCumulantsId { // purityAnalysis.addClone("purityAnalysis/Pi/", "purityAnalysis/El/"); // purityAnalysis.addClone("purityAnalysis/Pi/", "purityAnalysis/De/"); + // profiles adding for Reco + addFCRecoProfiles(); + addFCRecoProfiles(); + addFCRecoProfiles(); + addFCRecoProfiles(); + addFCRecoProfiles(); + addFCRecoProfiles(); + addFCRecoProfiles(); + addFCRecoProfiles(); + + addNetQVectorProfileHistograms(); + addNetQVectorProfileHistograms(); + addNetQVectorProfileHistograms(); + addNetQVectorProfileHistograms(); + + // profiles adding for Gen + addFCGenProfiles(); + addFCGenProfiles(); + addFCGenProfiles(); + addFCGenProfiles(); + addFCGenProfiles(); + addFCGenProfiles(); + addFCGenProfiles(); + addFCGenProfiles(); + + addFCGenProfiles(); + addFCGenProfiles(); + addFCGenProfiles(); + addFCGenProfiles(); + } // init ends static constexpr std::string_view HistRegDire2[] = { @@ -1028,7 +1369,7 @@ struct NchCumulantsId { const int& idMethodKa, const bool& trackIsKaon, float& nAKa, float& nKa, const int& idMethodPr, const bool& trackIsProton, float& nPr, float& nAPr, H& recoAnalysis) { - if (flagUnusedVariableError) + if (cfgEventSelection.flagUnusedVariableError) LOG(info) << trackIdTag << idMethodPi << ":" << idMethodKa << ":" << idMethodPr; if (track.sign() > 0) { // fillRecoTrackQA(recoAnalysis, track); @@ -1105,6 +1446,7 @@ struct NchCumulantsId { // recoAnalysis.fill(HIST("recoAnalysis/SelectedTrack_IdentificationTag"), trackIdTag); } + // fill the basis for factorial cumulants using MyAllTracks = soa::Join bool isEventSelected(const CollisionType& col) { - if (cfgEvSel01doNoSameBunchPileup && + if (cfgEventSelection.cfgEvSel01doNoSameBunchPileup && !col.selection_bit(aod::evsel::kNoSameBunchPileup)) return false; - if (cfgEvSel02doIsGoodZvtxFT0vsPV && + if (cfgEventSelection.cfgEvSel02doIsGoodZvtxFT0vsPV && !col.selection_bit(aod::evsel::kIsGoodZvtxFT0vsPV)) return false; - if (cfgEvSel03doIsGoodITSLayersAll && + if (cfgEventSelection.cfgEvSel03doIsGoodITSLayersAll && !col.selection_bit(aod::evsel::kIsGoodITSLayersAll)) return false; return true; @@ -1211,7 +1553,7 @@ struct NchCumulantsId { } // Reject electrons first - if (cfgDoRejectionForId && !selTrackForId(track)) { + if (cfgEventSelection.cfgDoRejectionForId && !selTrackForId(track)) { continue; } @@ -1317,7 +1659,7 @@ struct NchCumulantsId { int mcTag = getMCTag(track); // Reject electrons first - if (cfgDoRejectionForId && !selTrackForId(track)) { + if (cfgEventSelection.cfgDoRejectionForId && !selTrackForId(track)) { continue; } @@ -1514,7 +1856,8 @@ struct NchCumulantsId { void processSim(MyFilteredColsWithMcLabels const& collisions, MyFilteredTracksWithMcLabels const& tracks, aod::McCollisions const& mcCollisions, aod::McParticles const& mcParticles) { - if (flagUnusedVariableError) + + if (cfgEventSelection.flagUnusedVariableError) LOG(info) << mcCollisions.size(); bool trackIsPion = false; bool trackIsKaon = false; @@ -1537,7 +1880,7 @@ struct NchCumulantsId { continue; auto mcCollision = col.mcCollision(); - if (checkCollPosZMc && std::abs(mcCollision.posZ()) > cfgCutPosZ) + if (cfgEventSelection.checkCollPosZMc && std::abs(mcCollision.posZ()) > cfgCutPosZ) continue; // slice reco tracks to this collision @@ -1546,6 +1889,18 @@ struct NchCumulantsId { // slice mc particles to mc collisions const auto mcTracksTablePerMcColl = mcParticles.sliceBy(mcTracksPerMcCollisionPreslice, mcCollision.globalIndex()); + float cent = col.centFT0M(); + + EffPowerSums prPow; + EffPowerSums aprPow; + EffPowerSums pipPow; + EffPowerSums pimPow; + EffPowerSums kapPow; + EffPowerSums kamPow; + EffPowerSums posPow; + EffPowerSums negPow; + + cfgEventSelection.currentSubsample = static_cast(cfgEventSelection.fRandom->Uniform(0, cfgEventSelection.subSampleAxis.value[0])); // Denominator -- Generator level(truth) @@ -1617,12 +1972,28 @@ struct NchCumulantsId { nTGen = nPGen + nMGen; // ── Fill GEN sparse (denominator) ──────────────────────── - hist.fill(HIST("sim/gen/sparse1"), nChGen, nPGen, nMGen, - nPrGen, nAPrGen, nKaGen, nAKaGen, nTGen, - col.centFT0M()); - hist.fill(HIST("sim/gen/sparse2"), nChGen, nPGen, nMGen, - nPiGen, nAPiGen, nKaGen, nAKaGen, nTGen, - col.centFT0M()); + if (cfgEventSelection.fillSparseForGen) { + hist.fill(HIST("sim/gen/sparse1"), nChGen, nPGen, nMGen, + nPrGen, nAPrGen, nKaGen, nAKaGen, nTGen, + cent); + hist.fill(HIST("sim/gen/sparse2"), nChGen, nPGen, nMGen, + nPiGen, nAPiGen, nKaGen, nAKaGen, nTGen, + cent); + } + fillGenFactorialMoments(nPrGen, cent, registry); + fillGenFactorialMoments(nAPrGen, cent, registry); + fillGenFactorialMoments(nPiGen, cent, registry); + fillGenFactorialMoments(nAPiGen, cent, registry); + fillGenFactorialMoments(nKaGen, cent, registry); + fillGenFactorialMoments(nAKaGen, cent, registry); + fillGenFactorialMoments(nPGen, cent, registry); + fillGenFactorialMoments(nMGen, cent, registry); + + fillGenNetFactorialMoments(nPGen, nMGen, cent, registry); + fillGenNetFactorialMoments(nPrGen, nAPrGen, cent, registry); + fillGenNetFactorialMoments(nPiGen, nAPiGen, cent, registry); + fillGenNetFactorialMoments(nKaGen, nAKaGen, cent, registry); + // ── End of GEN level filling // // Numerator - Reconstructed + truth matched // reco->selFunc passed, no pdg @@ -1663,7 +2034,7 @@ struct NchCumulantsId { idMethodPr = kUnidentified; // Reject electrons first - if (cfgDoRejectionForId && !selTrackForId(track)) { + if (cfgEventSelection.cfgDoRejectionForId && !selTrackForId(track)) { continue; } @@ -1716,6 +2087,9 @@ struct NchCumulantsId { recoAnalysis.fill(HIST("recoAnalysis/Charge/Pos/h15_phi"), track.phi()); recoAnalysis.fill(HIST("recoAnalysis/Charge/Pos/h16_rapidity"), track.y()); recoAnalysis.fill(HIST("recoAnalysis/Charge/Pos/h20_pt_eta"), track.pt(), track.eta()); + + float weight = hPtEtaForEffCorrection[kCh][kPos]->GetBinContent(ptEtaBin); + fillEffPower(posPow, weight); } else if (track.sign() < 0) { nMRec += hPtEtaForEffCorrection[kCh][kNeg]->GetBinContent(ptEtaBin); recoAnalysis.fill(HIST("recoAnalysis/Charge/Neg/h12_p"), track.p()); @@ -1724,6 +2098,9 @@ struct NchCumulantsId { recoAnalysis.fill(HIST("recoAnalysis/Charge/Neg/h15_phi"), track.phi()); recoAnalysis.fill(HIST("recoAnalysis/Charge/Neg/h16_rapidity"), track.y()); recoAnalysis.fill(HIST("recoAnalysis/Charge/Neg/h20_pt_eta"), track.pt(), track.eta()); + + float weight = hPtEtaForEffCorrection[kCh][kNeg]->GetBinContent(ptEtaBin); + fillEffPower(negPow, weight); } // species reco — sel passes, PDG not checked (raw reco) @@ -1731,9 +2108,15 @@ struct NchCumulantsId { if (track.sign() > 0) { nPiRec += hPtEtaForEffCorrection[kPi][kPos]->GetBinContent(ptEtaBin); fillRecoTrackQA(recoAnalysis, track); + + float weight = hPtEtaForEffCorrection[kPi][kPos]->GetBinContent(ptEtaBin); + fillEffPower(pipPow, weight); } else if (track.sign() < 0) { nAPiRec += hPtEtaForEffCorrection[kPi][kNeg]->GetBinContent(ptEtaBin); fillRecoTrackQA(recoAnalysis, track); + + float weight = hPtEtaForEffCorrection[kPi][kNeg]->GetBinContent(ptEtaBin); + fillEffPower(pimPow, weight); } // PID band QA for pions if (idMethodPi == kTPCidentified) @@ -1744,9 +2127,13 @@ struct NchCumulantsId { if (track.sign() > 0) { nKaRec += hPtEtaForEffCorrection[kKa][kPos]->GetBinContent(ptEtaBin); fillRecoTrackQA(recoAnalysis, track); + float weight = hPtEtaForEffCorrection[kKa][kPos]->GetBinContent(ptEtaBin); + fillEffPower(kapPow, weight); } else if (track.sign() < 0) { nAKaRec += hPtEtaForEffCorrection[kKa][kNeg]->GetBinContent(ptEtaBin); fillRecoTrackQA(recoAnalysis, track); + float weight = hPtEtaForEffCorrection[kKa][kNeg]->GetBinContent(ptEtaBin); + fillEffPower(kamPow, weight); } // PID band QA for kaons if (idMethodKa == kTPCidentified) @@ -1757,9 +2144,13 @@ struct NchCumulantsId { if (track.sign() > 0) { nPrRec += hPtEtaForEffCorrection[kPr][kPos]->GetBinContent(ptEtaBin); fillRecoTrackQA(recoAnalysis, track); + float weight = hPtEtaForEffCorrection[kPr][kPos]->GetBinContent(ptEtaBin); + fillEffPower(prPow, weight); } else if (track.sign() < 0) { nAPrRec += hPtEtaForEffCorrection[kPr][kNeg]->GetBinContent(ptEtaBin); fillRecoTrackQA(recoAnalysis, track); + float weight = hPtEtaForEffCorrection[kPr][kNeg]->GetBinContent(ptEtaBin); + fillEffPower(aprPow, weight); } // PID band QA for protons if (idMethodPr == kTPCidentified) @@ -1839,24 +2230,42 @@ struct NchCumulantsId { fillPurityTrackQA(purityAnalysis, track); } } - } + } // reconstructed track loop ends nChRec = nPRec - nMRec; nTRec = nPRec + nMRec; nChPur = nPPur - nMPur; nTPur = nPPur + nMPur; // ── fill reco histos ───────────────────────────────────── - hist.fill(HIST("sim/reco/sparse1"), nChRec, nPRec, nMRec, - nPrRec, nAPrRec, nKaRec, nAKaRec, nTRec, col.centFT0M()); - hist.fill(HIST("sim/reco/sparse2"), nChRec, nPRec, nMRec, - nPiRec, nAPiRec, nKaRec, nAKaRec, nTRec, col.centFT0M()); + if (cfgEventSelection.fillSparseForReco) { + hist.fill(HIST("sim/reco/sparse1"), nChRec, nPRec, nMRec, + nPrRec, nAPrRec, nKaRec, nAKaRec, nTRec, cent); + hist.fill(HIST("sim/reco/sparse2"), nChRec, nPRec, nMRec, + nPiRec, nAPiRec, nKaRec, nAKaRec, nTRec, cent); + } // ── fill purity histos ─────────────────────────────────── - hist.fill(HIST("sim/purity/sparse1"), nChPur, nPPur, nMPur, - nPrPur, nAPrPur, nKaPur, nAKaPur, nTPur, col.centFT0M()); - hist.fill(HIST("sim/purity/sparse2"), nChPur, nPPur, nMPur, - nPiPur, nAPiPur, nKaPur, nAKaPur, nTPur, col.centFT0M()); - + if (cfgEventSelection.fillSparseForPurity) { + hist.fill(HIST("sim/purity/sparse1"), nChPur, nPPur, nMPur, + nPrPur, nAPrPur, nKaPur, nAKaPur, nTPur, cent); + hist.fill(HIST("sim/purity/sparse2"), nChPur, nPPur, nMPur, + nPiPur, nAPiPur, nKaPur, nAKaPur, nTPur, cent); + } + fillCollQA(col, nChRec, nTRec); + + fillFCBasis(prPow, cent, registry); + fillFCBasis(aprPow, cent, registry); + fillFCBasis(pipPow, cent, registry); + fillFCBasis(pimPow, cent, registry); + fillFCBasis(kapPow, cent, registry); + fillFCBasis(kamPow, cent, registry); + fillFCBasis(posPow, cent, registry); + fillFCBasis(negPow, cent, registry); + + fillNetQVectorProfileHistograms(posPow, negPow, cent, registry); + fillNetQVectorProfileHistograms(prPow, aprPow, cent, registry); + fillNetQVectorProfileHistograms(pipPow, pimPow, cent, registry); + fillNetQVectorProfileHistograms(kapPow, kamPow, cent, registry); } // common collision loop ends } // process sim ends PROCESS_SWITCH(NchCumulantsId, processSim, "Process Sim: Gen + Reco + Purity", true); From d0814e28006e240e873e7498b20b8f1e1fae2db8 Mon Sep 17 00:00:00 2001 From: Pravata Panigrahi <158029188+pravatp98@users.noreply.github.com> Date: Fri, 14 Aug 2026 22:18:46 +0530 Subject: [PATCH 02/13] Refactor registry parameter to hReg in functions --- .../EbyEFluctuations/Tasks/nchCumulantsId.cxx | 162 +++++++++--------- 1 file changed, 81 insertions(+), 81 deletions(-) diff --git a/PWGCF/EbyEFluctuations/Tasks/nchCumulantsId.cxx b/PWGCF/EbyEFluctuations/Tasks/nchCumulantsId.cxx index 3c23ef7558f..1476f306e64 100644 --- a/PWGCF/EbyEFluctuations/Tasks/nchCumulantsId.cxx +++ b/PWGCF/EbyEFluctuations/Tasks/nchCumulantsId.cxx @@ -405,42 +405,42 @@ struct NchCumulantsId { } template - void fillFCBasis(const EffPowerSums& p, float cent, Registry& registry) + void fillFCBasis(const EffPowerSums& p, float cent, Registry& hReg) { auto base = HIST(FCRecoDir[static_cast(Prefix)]); // Fill Factorial Moments for Reco directly to ensure Reco stores the same type of data as Gen // Compute Efficiency Corrected Factorial Moments from power sums // These formulas come from elementary symmetric polynomials and their relation to power sums - registry.fill(base + HIST("Q1"), cent, p.q1); - registry.fill(base + HIST("Q1Sq"), cent, p.q1 * p.q1); - registry.fill(base + HIST("Q2"), cent, p.q2); - registry.fill(base + HIST("Q1Cube"), cent, p.q1 * p.q1 * p.q1); - registry.fill(base + HIST("Q1Q2"), cent, p.q1 * p.q2); - registry.fill(base + HIST("Q3"), cent, p.q3); - registry.fill(base + HIST("Q1Pow4"), cent, p.q1 * p.q1 * p.q1 * p.q1); - registry.fill(base + HIST("Q1SqQ2"), cent, p.q1 * p.q1 * p.q2); - registry.fill(base + HIST("Q2Sq"), cent, p.q2 * p.q2); - registry.fill(base + HIST("Q1Q3"), cent, p.q1 * p.q3); - registry.fill(base + HIST("Q4"), cent, p.q4); + hReg.fill(base + HIST("Q1"), cent, p.q1); + hReg.fill(base + HIST("Q1Sq"), cent, p.q1 * p.q1); + hReg.fill(base + HIST("Q2"), cent, p.q2); + hReg.fill(base + HIST("Q1Cube"), cent, p.q1 * p.q1 * p.q1); + hReg.fill(base + HIST("Q1Q2"), cent, p.q1 * p.q2); + hReg.fill(base + HIST("Q3"), cent, p.q3); + hReg.fill(base + HIST("Q1Pow4"), cent, p.q1 * p.q1 * p.q1 * p.q1); + hReg.fill(base + HIST("Q1SqQ2"), cent, p.q1 * p.q1 * p.q2); + hReg.fill(base + HIST("Q2Sq"), cent, p.q2 * p.q2); + hReg.fill(base + HIST("Q1Q3"), cent, p.q1 * p.q3); + hReg.fill(base + HIST("Q4"), cent, p.q4); if (cfgEventSelection.cfgDoSubsampling) { - registry.fill(base + HIST("Q1_subsample"), cent, cfgEventSelection.currentSubsample, p.q1); - registry.fill(base + HIST("Q1Sq_subsample"), cent, cfgEventSelection.currentSubsample, p.q1 * p.q1); - registry.fill(base + HIST("Q2_subsample"), cent, cfgEventSelection.currentSubsample, p.q2); - registry.fill(base + HIST("Q1Cube_subsample"), cent, cfgEventSelection.currentSubsample, p.q1 * p.q1 * p.q1); - registry.fill(base + HIST("Q1Q2_subsample"), cent, cfgEventSelection.currentSubsample, p.q1 * p.q2); - registry.fill(base + HIST("Q3_subsample"), cent, cfgEventSelection.currentSubsample, p.q3); - registry.fill(base + HIST("Q1Pow4_subsample"), cent, cfgEventSelection.currentSubsample, p.q1 * p.q1 * p.q1 * p.q1); - registry.fill(base + HIST("Q1SqQ2_subsample"), cent, cfgEventSelection.currentSubsample, p.q1 * p.q1 * p.q2); - registry.fill(base + HIST("Q2Sq_subsample"), cent, cfgEventSelection.currentSubsample, p.q2 * p.q2); - registry.fill(base + HIST("Q1Q3_subsample"), cent, cfgEventSelection.currentSubsample, p.q1 * p.q3); - registry.fill(base + HIST("Q4_subsample"), cent, cfgEventSelection.currentSubsample, p.q4); + hReg.fill(base + HIST("Q1_subsample"), cent, cfgEventSelection.currentSubsample, p.q1); + hReg.fill(base + HIST("Q1Sq_subsample"), cent, cfgEventSelection.currentSubsample, p.q1 * p.q1); + hReg.fill(base + HIST("Q2_subsample"), cent, cfgEventSelection.currentSubsample, p.q2); + hReg.fill(base + HIST("Q1Cube_subsample"), cent, cfgEventSelection.currentSubsample, p.q1 * p.q1 * p.q1); + hReg.fill(base + HIST("Q1Q2_subsample"), cent, cfgEventSelection.currentSubsample, p.q1 * p.q2); + hReg.fill(base + HIST("Q3_subsample"), cent, cfgEventSelection.currentSubsample, p.q3); + hReg.fill(base + HIST("Q1Pow4_subsample"), cent, cfgEventSelection.currentSubsample, p.q1 * p.q1 * p.q1 * p.q1); + hReg.fill(base + HIST("Q1SqQ2_subsample"), cent, cfgEventSelection.currentSubsample, p.q1 * p.q1 * p.q2); + hReg.fill(base + HIST("Q2Sq_subsample"), cent, cfgEventSelection.currentSubsample, p.q2 * p.q2); + hReg.fill(base + HIST("Q1Q3_subsample"), cent, cfgEventSelection.currentSubsample, p.q1 * p.q3); + hReg.fill(base + HIST("Q4_subsample"), cent, cfgEventSelection.currentSubsample, p.q4); } } // fill function for reco net species template - void fillNetQVectorProfileHistograms(const EffPowerSums& posPow, const EffPowerSums& negPow, float cent, Registry& registry) + void fillNetQVectorProfileHistograms(const EffPowerSums& posPow, const EffPowerSums& negPow, float cent, Registry& hReg) { // Construct base directory path auto base = HIST(FCRecoDir[static_cast(Prefix)]); @@ -456,50 +456,50 @@ struct NchCumulantsId { // Fill joint pos-neg correction profiles // capture inline from pos and neg - double jointF11 = posPow.q1 * negPow.q1; - double cf2pos = posPow.q1 * posPow.q1 - posPow.q2; - double cf2neg = negPow.q1 * negPow.q1 - negPow.q2; - double jointF12 = posPow.q1 * cf2neg; - double jointF21 = negPow.q1 * cf2pos; + double jointF11 = posPow.q1 * negPow.q1; + double cf2pos = posPow.q1 * posPow.q1 - posPow.q2; + double cf2neg = negPow.q1 * negPow.q1 - negPow.q2; + double jointF12 = posPow.q1 * cf2neg; + double jointF21 = negPow.q1 * cf2pos; // ─── Fill raw q-vectors (for cross-verification) - registry.fill(base + HIST("Q_net_1"), cent, qNet1); - registry.fill(base + HIST("Q_net_1Sq"), cent, qNet1 * qNet1); - registry.fill(base + HIST("Q_net_2"), cent, qNet2); - registry.fill(base + HIST("Q_net_1Cube"), cent, qNet1 * qNet1 * qNet1); - registry.fill(base + HIST("Q_net_1Q_net_2"), cent, qNet1 * qNet2); - registry.fill(base + HIST("Q_net_3"), cent, qNet3); - registry.fill(base + HIST("Q_net_1Pow4"), cent, qNet1 * qNet1 * qNet1 * qNet1); - registry.fill(base + HIST("Q_net_1SqQ_net_2"), cent, qNet1 * qNet1 * qNet2); - registry.fill(base + HIST("Q_net_2Sq"), cent, qNet2 * qNet2); - registry.fill(base + HIST("Q_net_1Q_net_3"), cent, qNet1 * qNet3); - registry.fill(base + HIST("Q_net_4"), cent, qNet4); + hReg.fill(base + HIST("Q_net_1"), cent, qNet1); + hReg.fill(base + HIST("Q_net_1Sq"), cent, qNet1 * qNet1); + hReg.fill(base + HIST("Q_net_2"), cent, qNet2); + hReg.fill(base + HIST("Q_net_1Cube"), cent, qNet1 * qNet1 * qNet1); + hReg.fill(base + HIST("Q_net_1Q_net_2"), cent, qNet1 * qNet2); + hReg.fill(base + HIST("Q_net_3"), cent, qNet3); + hReg.fill(base + HIST("Q_net_1Pow4"), cent, qNet1 * qNet1 * qNet1 * qNet1); + hReg.fill(base + HIST("Q_net_1SqQ_net_2"), cent, qNet1 * qNet1 * qNet2); + hReg.fill(base + HIST("Q_net_2Sq"), cent, qNet2 * qNet2); + hReg.fill(base + HIST("Q_net_1Q_net_3"), cent, qNet1 * qNet3); + hReg.fill(base + HIST("Q_net_4"), cent, qNet4); // ─── Fill factorial moments ─── - registry.fill(base + HIST("JointF11"), cent, jointF11); - registry.fill(base + HIST("JointF12"), cent, jointF12); - registry.fill(base + HIST("JointF21"), cent, jointF21); + hReg.fill(base + HIST("JointF11"), cent, jointF11); + hReg.fill(base + HIST("JointF12"), cent, jointF12); + hReg.fill(base + HIST("JointF21"), cent, jointF21); if (cfgEventSelection.cfgDoSubsampling) { - registry.fill(base + HIST("Q_net_1_subsample"), cent, cfgEventSelection.currentSubsample, qNet1); - registry.fill(base + HIST("Q_net_1Sq_subsample"), cent, cfgEventSelection.currentSubsample, qNet1 * qNet1); - registry.fill(base + HIST("Q_net_2_subsample"), cent, cfgEventSelection.currentSubsample, qNet2); - registry.fill(base + HIST("Q_net_1Cube_subsample"), cent, cfgEventSelection.currentSubsample, qNet1 * qNet1 * qNet1); - registry.fill(base + HIST("Q_net_1Q_net_2_subsample"), cent, cfgEventSelection.currentSubsample, qNet1 * qNet2); - registry.fill(base + HIST("Q_net_3_subsample"), cent, cfgEventSelection.currentSubsample, qNet3); - registry.fill(base + HIST("Q_net_1Pow4_subsample"), cent, cfgEventSelection.currentSubsample, qNet1 * qNet1 * qNet1 * qNet1); - registry.fill(base + HIST("Q_net_1SqQ_net_2_subsample"), cent, cfgEventSelection.currentSubsample, qNet1 * qNet1 * qNet2); - registry.fill(base + HIST("Q_net_2Sq_subsample"), cent, cfgEventSelection.currentSubsample, qNet2 * qNet2); - registry.fill(base + HIST("Q_net_1Q_net_3_subsample"), cent, cfgEventSelection.currentSubsample, qNet1 * qNet3); - registry.fill(base + HIST("Q_net_4_subsample"), cent, cfgEventSelection.currentSubsample, qNet4); - - registry.fill(base + HIST("JointF11_subsample"), cent, cfgEventSelection.currentSubsample, jointF11); - registry.fill(base + HIST("JointF12_subsample"), cent, cfgEventSelection.currentSubsample, jointF12); - registry.fill(base + HIST("JointF21_subsample"), cent, cfgEventSelection.currentSubsample, jointF21); + hReg.fill(base + HIST("Q_net_1_subsample"), cent, cfgEventSelection.currentSubsample, qNet1); + hReg.fill(base + HIST("Q_net_1Sq_subsample"), cent, cfgEventSelection.currentSubsample, qNet1 * qNet1); + hReg.fill(base + HIST("Q_net_2_subsample"), cent, cfgEventSelection.currentSubsample, qNet2); + hReg.fill(base + HIST("Q_net_1Cube_subsample"), cent, cfgEventSelection.currentSubsample, qNet1 * qNet1 * qNet1); + hReg.fill(base + HIST("Q_net_1Q_net_2_subsample"), cent, cfgEventSelection.currentSubsample, qNet1 * qNet2); + hReg.fill(base + HIST("Q_net_3_subsample"), cent, cfgEventSelection.currentSubsample, qNet3); + hReg.fill(base + HIST("Q_net_1Pow4_subsample"), cent, cfgEventSelection.currentSubsample, qNet1 * qNet1 * qNet1 * qNet1); + hReg.fill(base + HIST("Q_net_1SqQ_net_2_subsample"), cent, cfgEventSelection.currentSubsample, qNet1 * qNet1 * qNet2); + hReg.fill(base + HIST("Q_net_2Sq_subsample"), cent, cfgEventSelection.currentSubsample, qNet2 * qNet2); + hReg.fill(base + HIST("Q_net_1Q_net_3_subsample"), cent, cfgEventSelection.currentSubsample, qNet1 * qNet3); + hReg.fill(base + HIST("Q_net_4_subsample"), cent, cfgEventSelection.currentSubsample, qNet4); + + hReg.fill(base + HIST("JointF11_subsample"), cent, cfgEventSelection.currentSubsample, jointF11); + hReg.fill(base + HIST("JointF12_subsample"), cent, cfgEventSelection.currentSubsample, jointF12); + hReg.fill(base + HIST("JointF21_subsample"), cent, cfgEventSelection.currentSubsample, jointF21); } } template - void fillGenFactorialMoments(int n, float cent, Registry& registry) + void fillGenFactorialMoments(int n, float cent, Registry& hReg) { auto base = HIST(FCGenDir[static_cast(Prefix)]); @@ -508,21 +508,21 @@ struct NchCumulantsId { double f3 = n * (n - 1.) * (n - 2.); double f4 = n * (n - 1.) * (n - 2.) * (n - 3.); - registry.fill(base + HIST("F1"), cent, f1); - registry.fill(base + HIST("F2"), cent, f2); - registry.fill(base + HIST("F3"), cent, f3); - registry.fill(base + HIST("F4"), cent, f4); + hReg.fill(base + HIST("F1"), cent, f1); + hReg.fill(base + HIST("F2"), cent, f2); + hReg.fill(base + HIST("F3"), cent, f3); + hReg.fill(base + HIST("F4"), cent, f4); if (cfgEventSelection.cfgDoSubsampling) { - registry.fill(base + HIST("F1_subsample"), cent, cfgEventSelection.currentSubsample, f1); - registry.fill(base + HIST("F2_subsample"), cent, cfgEventSelection.currentSubsample, f2); - registry.fill(base + HIST("F3_subsample"), cent, cfgEventSelection.currentSubsample, f3); - registry.fill(base + HIST("F4_subsample"), cent, cfgEventSelection.currentSubsample, f4); + hReg.fill(base + HIST("F1_subsample"), cent, cfgEventSelection.currentSubsample, f1); + hReg.fill(base + HIST("F2_subsample"), cent, cfgEventSelection.currentSubsample, f2); + hReg.fill(base + HIST("F3_subsample"), cent, cfgEventSelection.currentSubsample, f3); + hReg.fill(base + HIST("F4_subsample"), cent, cfgEventSelection.currentSubsample, f4); } } // ─── GEN level (net speciestruth) ─── template - void fillGenNetFactorialMoments(float nPos, float nNeg, float cent, Registry& registry) + void fillGenNetFactorialMoments(float nPos, float nNeg, float cent, Registry& hReg) { auto base = HIST(FCGenDir[static_cast(Prefix)]); @@ -533,16 +533,16 @@ struct NchCumulantsId { double f3 = nNet * (nNet - 1.0) * (nNet - 2.0); double f4 = nNet * (nNet - 1.0) * (nNet - 2.0) * (nNet - 3.0); - registry.fill(base + HIST("F1"), cent, f1); - registry.fill(base + HIST("F2"), cent, f2); - registry.fill(base + HIST("F3"), cent, f3); - registry.fill(base + HIST("F4"), cent, f4); + hReg.fill(base + HIST("F1"), cent, f1); + hReg.fill(base + HIST("F2"), cent, f2); + hReg.fill(base + HIST("F3"), cent, f3); + hReg.fill(base + HIST("F4"), cent, f4); if (cfgEventSelection.cfgDoSubsampling) { - registry.fill(base + HIST("F1_subsample"), cent, cfgEventSelection.currentSubsample, f1); - registry.fill(base + HIST("F2_subsample"), cent, cfgEventSelection.currentSubsample, f2); - registry.fill(base + HIST("F3_subsample"), cent, cfgEventSelection.currentSubsample, f3); - registry.fill(base + HIST("F4_subsample"), cent, cfgEventSelection.currentSubsample, f4); + hReg.fill(base + HIST("F1_subsample"), cent, cfgEventSelection.currentSubsample, f1); + hReg.fill(base + HIST("F2_subsample"), cent, cfgEventSelection.currentSubsample, f2); + hReg.fill(base + HIST("F3_subsample"), cent, cfgEventSelection.currentSubsample, f3); + hReg.fill(base + HIST("F4_subsample"), cent, cfgEventSelection.currentSubsample, f4); } } @@ -1257,15 +1257,15 @@ struct NchCumulantsId { } template - void fillCollQA(const T& col, const int& nCh, const int& nT) + void fillCollQA(const T& coll, const int& nCh, const int& nT) { - hist.fill(HIST(HistRegDire[mode]) + HIST("h_VtxZ"), col.posZ()); + hist.fill(HIST(HistRegDire[mode]) + HIST("h_VtxZ"), coll.posZ()); hist.fill(HIST(HistRegDire[mode]) + HIST("h_Counts"), 0.5); - hist.fill(HIST(HistRegDire[mode]) + HIST("multFT0"), col.multFT0C()); - hist.fill(HIST(HistRegDire[mode]) + HIST("centFT0"), col.centFT0M()); + hist.fill(HIST(HistRegDire[mode]) + HIST("multFT0"), coll.multFT0C()); + hist.fill(HIST(HistRegDire[mode]) + HIST("centFT0"), coll.centFT0M()); if (mode == qaEventPostSel) { hist.fill(HIST(HistRegDire[mode]) + HIST("net_charge"), nCh); - hist.fill(HIST(HistRegDire[mode]) + HIST("Nt_centFT"), col.centFT0M(), nT); + hist.fill(HIST(HistRegDire[mode]) + HIST("Nt_centFT"), coll.centFT0M(), nT); } } From 0f85e5abb9e0519ecfa49992d95c4437b351bf95 Mon Sep 17 00:00:00 2001 From: Pravata Panigrahi <158029188+pravatp98@users.noreply.github.com> Date: Fri, 14 Aug 2026 22:37:12 +0530 Subject: [PATCH 03/13] Refactor track analysis to use hReg for histograms --- .../EbyEFluctuations/Tasks/nchCumulantsId.cxx | 82 +++++++++---------- 1 file changed, 41 insertions(+), 41 deletions(-) diff --git a/PWGCF/EbyEFluctuations/Tasks/nchCumulantsId.cxx b/PWGCF/EbyEFluctuations/Tasks/nchCumulantsId.cxx index 1476f306e64..db822fa9d08 100644 --- a/PWGCF/EbyEFluctuations/Tasks/nchCumulantsId.cxx +++ b/PWGCF/EbyEFluctuations/Tasks/nchCumulantsId.cxx @@ -1367,34 +1367,34 @@ struct NchCumulantsId { void executeTrackAnalysisPart(const T& track, const int& trackIdTag, float& nP, float& nM, const int& idMethodPi, const bool& trackIsPion, float& nAPi, float& nPi, const int& idMethodKa, const bool& trackIsKaon, float& nAKa, float& nKa, - const int& idMethodPr, const bool& trackIsProton, float& nPr, float& nAPr, H& recoAnalysis) + const int& idMethodPr, const bool& trackIsProton, float& nPr, float& nAPr, H& hReg) { if (cfgEventSelection.flagUnusedVariableError) LOG(info) << trackIdTag << idMethodPi << ":" << idMethodKa << ":" << idMethodPr; if (track.sign() > 0) { - // fillRecoTrackQA(recoAnalysis, track); + // fillRecoTrackQA(hReg, track); nP++; - recoAnalysis.fill(HIST("recoAnalysis/Charge/Pos/h12_p"), track.p()); - recoAnalysis.fill(HIST("recoAnalysis/Charge/Pos/h13_pt"), track.pt()); - recoAnalysis.fill(HIST("recoAnalysis/Charge/Pos/h14_eta"), track.eta()); - recoAnalysis.fill(HIST("recoAnalysis/Charge/Pos/h15_phi"), track.phi()); - recoAnalysis.fill(HIST("recoAnalysis/Charge/Pos/h16_rapidity"), track.y()); - recoAnalysis.fill(HIST("recoAnalysis/Charge/Pos/h20_pt_eta"), track.pt(), track.eta()); + hReg.fill(HIST("recoAnalysis/Charge/Pos/h12_p"), track.p()); + hReg.fill(HIST("recoAnalysis/Charge/Pos/h13_pt"), track.pt()); + hReg.fill(HIST("recoAnalysis/Charge/Pos/h14_eta"), track.eta()); + hReg.fill(HIST("recoAnalysis/Charge/Pos/h15_phi"), track.phi()); + hReg.fill(HIST("recoAnalysis/Charge/Pos/h16_rapidity"), track.y()); + hReg.fill(HIST("recoAnalysis/Charge/Pos/h20_pt_eta"), track.pt(), track.eta()); } if (track.sign() < 0) { - // fillRecoTrackQA(recoAnalysis, track); + // fillRecoTrackQA(hReg, track); nM++; - recoAnalysis.fill(HIST("recoAnalysis/Charge/Neg/h12_p"), track.p()); - recoAnalysis.fill(HIST("recoAnalysis/Charge/Neg/h13_pt"), track.pt()); - recoAnalysis.fill(HIST("recoAnalysis/Charge/Neg/h14_eta"), track.eta()); - recoAnalysis.fill(HIST("recoAnalysis/Charge/Neg/h15_phi"), track.phi()); - recoAnalysis.fill(HIST("recoAnalysis/Charge/Neg/h16_rapidity"), track.y()); - recoAnalysis.fill(HIST("recoAnalysis/Charge/Neg/h20_pt_eta"), track.pt(), track.eta()); + hReg.fill(HIST("recoAnalysis/Charge/Neg/h12_p"), track.p()); + hReg.fill(HIST("recoAnalysis/Charge/Neg/h13_pt"), track.pt()); + hReg.fill(HIST("recoAnalysis/Charge/Neg/h14_eta"), track.eta()); + hReg.fill(HIST("recoAnalysis/Charge/Neg/h15_phi"), track.phi()); + hReg.fill(HIST("recoAnalysis/Charge/Neg/h16_rapidity"), track.y()); + hReg.fill(HIST("recoAnalysis/Charge/Neg/h20_pt_eta"), track.pt(), track.eta()); } if (trackIsPion) { // if (idMethodPi == kTPCidentified) { - // fillIdentificationQA(hist, track); // set hist as recoAnalysis after tpcId etc add true + // fillIdentificationQA(hist, track); // set hist as hReg after tpcId etc add true // } else if (idMethodPi == kTPCTOFidentified) { // fillIdentificationQA(hist, track); // } else if (idMethodPi == kUnidentified) { @@ -1402,12 +1402,12 @@ struct NchCumulantsId { // } if (track.sign() > 0) { nPi++; - fillRecoTrackQA(recoAnalysis, track); + fillRecoTrackQA(hReg, track); } else if (track.sign() < 0) { nAPi++; - fillRecoTrackQA(recoAnalysis, track); + fillRecoTrackQA(hReg, track); } - // fillRecoTrackQA(recoAnalysis, track); + // fillRecoTrackQA(hReg, track); } if (trackIsKaon) { // if (idMethodKa == kTPCidentified) { @@ -1419,12 +1419,12 @@ struct NchCumulantsId { // } if (track.sign() > 0) { nKa++; - fillRecoTrackQA(recoAnalysis, track); + fillRecoTrackQA(hReg, track); } else if (track.sign() < 0) { nAKa++; - fillRecoTrackQA(recoAnalysis, track); + fillRecoTrackQA(hReg, track); } - // fillRecoTrackQA(recoAnalysis, track); + // fillRecoTrackQA(hReg, track); } if (trackIsProton) { // if (idMethodPr == kTPCidentified) { @@ -1436,15 +1436,15 @@ struct NchCumulantsId { // } if (track.sign() > 0) { nPr++; - fillRecoTrackQA(recoAnalysis, track); + fillRecoTrackQA(hReg, track); } else if (track.sign() < 0) { nAPr++; - fillRecoTrackQA(recoAnalysis, track); + fillRecoTrackQA(hReg, track); } - // fillRecoTrackQA(recoAnalysis, track); + // fillRecoTrackQA(hReg, track); } - // recoAnalysis.fill(HIST("recoAnalysis/SelectedTrack_IdentificationTag"), trackIdTag); + // hReg.fill(HIST("recoAnalysis/SelectedTrack_IdentificationTag"), trackIdTag); } // fill the basis for factorial cumulants @@ -1459,22 +1459,22 @@ struct NchCumulantsId { using MyCollisionsWithMcLabels = soa::Join; template - bool isEventSelected(const CollisionType& col) + bool isEventSelected(const CollisionType& coll) { if (cfgEventSelection.cfgEvSel01doNoSameBunchPileup && - !col.selection_bit(aod::evsel::kNoSameBunchPileup)) + !coll.selection_bit(aod::evsel::kNoSameBunchPileup)) return false; if (cfgEventSelection.cfgEvSel02doIsGoodZvtxFT0vsPV && - !col.selection_bit(aod::evsel::kIsGoodZvtxFT0vsPV)) + !coll.selection_bit(aod::evsel::kIsGoodZvtxFT0vsPV)) return false; if (cfgEventSelection.cfgEvSel03doIsGoodITSLayersAll && - !col.selection_bit(aod::evsel::kIsGoodITSLayersAll)) + !coll.selection_bit(aod::evsel::kIsGoodITSLayersAll)) return false; return true; } // tracks and collision filters - Filter col = aod::evsel::sel8 == true; + Filter colSel8 = aod::evsel::sel8 == true; Filter colFilter = nabs(aod::collision::posZ) < cfgCutPosZ; Filter trackFilter = requireGlobalTrackInFilter(); Filter trackPt = (aod::track::pt > cfgCutPtMin) && (aod::track::pt < cfgCutPtMax); @@ -1620,16 +1620,16 @@ struct NchCumulantsId { } if (!isEventSelected(col)) continue; - float nP = 0; - float nM = 0; - float nCh = 0; - float nT = 0; - float nPr = 0; - float nAPr = 0; - float nKa = 0; - float nAKa = 0; - float nPi = 0; - float nAPi = 0; + nP = 0; + nM = 0; + nCh = 0; + nT = 0; + nPr = 0; + nAPr = 0; + nKa = 0; + nAKa = 0; + nPi = 0; + nAPi = 0; // group tracks manually with corresponding collision using col id; const uint64_t collIdx = col.globalIndex(); const auto tracksTablePerColl = tracks.sliceBy(mctracksPerCollisionPreslice, collIdx); From ef11977571dc5de1cf97b4097d9c968e7debf5a9 Mon Sep 17 00:00:00 2001 From: Pravata Panigrahi <158029188+pravatp98@users.noreply.github.com> Date: Mon, 17 Aug 2026 16:40:06 +0530 Subject: [PATCH 04/13] Refactor identification bits and clean up includes and updated for Efficiency correction Refactor identification bits to use constexpr instead of #define, update includes to comment out unused headers, and improve readability of various functions. --- .../EbyEFluctuations/Tasks/nchCumulantsId.cxx | 284 +++++++++++------- 1 file changed, 171 insertions(+), 113 deletions(-) diff --git a/PWGCF/EbyEFluctuations/Tasks/nchCumulantsId.cxx b/PWGCF/EbyEFluctuations/Tasks/nchCumulantsId.cxx index db822fa9d08..11c5f42df95 100644 --- a/PWGCF/EbyEFluctuations/Tasks/nchCumulantsId.cxx +++ b/PWGCF/EbyEFluctuations/Tasks/nchCumulantsId.cxx @@ -23,7 +23,7 @@ #include #include -#include +// #include #include #include #include @@ -37,11 +37,11 @@ #include #include -#include -#include +// #include +// #include #include -#include +// #include #include #include #include @@ -56,17 +56,17 @@ using namespace o2::framework::expressions; using namespace o2::constants::physics; // for constants using namespace std; -#define ID_BIT_PI 0 // Identificationi bits for PID checks -#define ID_BIT_KA 1 -#define ID_BIT_PR 2 -#define ID_BIT_EL 3 -#define ID_BIT_DE 4 +constexpr int IdBitPI = 0; // Identificationi bits for PID checks +constexpr int IdBitKA = 1; +constexpr int IdBitPR = 2; +constexpr int IdBitEL = 3; +constexpr int IdBitDE = 4; -#define MC_BIT_PI 0 // MC particle identification bits -#define MC_BIT_KA 1 -#define MC_BIT_PR 2 -#define MC_BIT_EL 3 -#define MC_BIT_DE 4 +constexpr int McBitPI = 0; // MC particle identification bits +constexpr int McBitKA = 1; +constexpr int McBitPR = 2; +constexpr int McBitEL = 3; +constexpr int McBitDE = 4; #define BITSET(mask, ithBit) ((mask) |= (1 << (ithBit))) // avoid name bitset as std::bitset is already there #define BITCHECK(mask, ithBit) ((mask) & (1 << (ithBit))) // bit check will return int value, not bool, use BITCHECK != 0 in Analysi @@ -100,7 +100,7 @@ enum class FCPrefixEnum { NetCh = 11 }; -static constexpr std::string_view FCRecoDir[] = { +static constexpr std::array FCRecoDir = { "Reco/Pr/", "Reco/APr/", "Reco/PiPos/", @@ -114,7 +114,7 @@ static constexpr std::string_view FCRecoDir[] = { "Reco/NetPr/", "Reco/NetCh/"}; -static constexpr std::string_view FCGenDir[] = { +static constexpr std::array FCGenDir = { "Gen/Pr/", "Gen/APr/", "Gen/PiPos/", @@ -128,7 +128,7 @@ static constexpr std::string_view FCGenDir[] = { "Gen/NetPr/", "Gen/NetCh/"}; -static constexpr std::string_view PidDire[] = { +static constexpr std::array PidDire = { "Ch/", "Pi/", "Ka/", @@ -136,7 +136,7 @@ static constexpr std::string_view PidDire[] = { "El/", "De/"}; -static constexpr std::string_view ChargeDire[] = { +static constexpr std::array ChargeDire = { "Pos/", "Neg/"}; @@ -146,9 +146,8 @@ std::string getModifiedStr(const std::string& myString) if (pos != std::string::npos) { std::string subString = myString.substr(0, pos); // remove "/" from end of the string return subString; - } else { - return myString; } + return myString; } struct NchCumulantsId { @@ -160,7 +159,7 @@ struct NchCumulantsId { HistogramRegistry registry{"registry", {}, OutputObjHandlingPolicy::AnalysisObject}; // PDG data base - Service pdgDB; + // Service pdgDB; Configurable cfgCutPosZ{"cfgCutPosZ", 10.0, "cut for vertex Z"}; Configurable cfgCutDcaXY{"cfgCutDcaXY", 0.12, "cut for dcaXY"}; @@ -499,7 +498,7 @@ struct NchCumulantsId { } template - void fillGenFactorialMoments(int n, float cent, Registry& hReg) + void fillGenFactorialMoments(float n, float cent, Registry& hReg) { auto base = HIST(FCGenDir[static_cast(Prefix)]); @@ -553,25 +552,29 @@ struct NchCumulantsId { mgr.setCaching(true); auto ccdbObj = mgr.getForTimeStamp(cfgCCDB.cfgCCDB02Path, cfgCCDB.cfgCCDB03SOR); if (!ccdbObj) { - if (cfgDebug.printDebugMessages) + if (cfgDebug.printDebugMessages) { LOG(info) << "DEBUG :: CCDB OBJECT NOT FOUND"; + } } else { - if (cfgDebug.printDebugMessages) + if (cfgDebug.printDebugMessages) { LOG(info) << "DEBUG :: CCDB OBJECT FOUND"; + } } ccdbObj->Print(); - hPtEtaForBinSearch = reinterpret_cast(ccdbObj->FindObject("hPtEta")); - if (cfgDebug.printDebugMessages) + hPtEtaForBinSearch = dynamic_cast(ccdbObj->FindObject("hPtEta")); + if (cfgDebug.printDebugMessages) { LOG(info) << "DEBUG :: Obj Name = " << hPtEtaForBinSearch->GetName() << " :: entries = " << hPtEtaForBinSearch->GetEntries(); - std::string name = ""; + } + std::string name; for (int i = 0; i <= kDe; i++) { for (int j = 0; j < (kNeg + 1); j++) { name = "hPtEta" + getModifiedStr(static_cast(PidDire[i])) + getModifiedStr(static_cast(ChargeDire[j])); - hPtEtaForEffCorrection[i][j] = reinterpret_cast(ccdbObj->FindObject(name.c_str())); - if (cfgDebug.printDebugMessages) + hPtEtaForEffCorrection[i][j] = dynamic_cast(ccdbObj->FindObject(name.c_str())); + if (cfgDebug.printDebugMessages) { LOG(info) << "DEBUG :: Obj Name = " << hPtEtaForEffCorrection[i][j]->GetName() << " :: entries = " << hPtEtaForBinSearch->GetEntries(); + } } } @@ -816,7 +819,7 @@ struct NchCumulantsId { } // init ends - static constexpr std::string_view HistRegDire2[] = { + static constexpr std::array HistRegDire2 = { "v0Table/Full/", "v0Table/postK0sCheck/", "v0Table/postMassCut/", @@ -863,7 +866,7 @@ struct NchCumulantsId { NoId }; - static constexpr std::string_view DetDire[] = { + static constexpr std::array DetDire = { "tpcId/", "tofId/", "tpctofId/", @@ -877,7 +880,7 @@ struct NchCumulantsId { qaTracksIdfd, }; - static constexpr std::string_view HistRegDire[] = { + static constexpr std::array HistRegDire = { "QA/events/preSel/", "QA/events/postSel/", "QA/tracks/preSel/", @@ -899,27 +902,32 @@ struct NchCumulantsId { bool vetoIdOthersTPC(const T& track) { if (pidMode != kPi) { - if (std::fabs(track.tpcNSigmaPi()) < cfgVetoIdCut.cfgVetoId01PiTPC) + if (std::fabs(track.tpcNSigmaPi()) < cfgVetoIdCut.cfgVetoId01PiTPC) { return false; + } } if (pidMode != kKa) { - if (std::fabs(track.tpcNSigmaKa()) < cfgVetoIdCut.cfgVetoId03KaTPC) + if (std::fabs(track.tpcNSigmaKa()) < cfgVetoIdCut.cfgVetoId03KaTPC) { return false; + } } if (pidMode != kPr) { - if (std::fabs(track.tpcNSigmaPr()) < cfgVetoIdCut.cfgVetoId05PrTPC) + if (std::fabs(track.tpcNSigmaPr()) < cfgVetoIdCut.cfgVetoId05PrTPC) { return false; + } } if (cfgId02DoElRejection) { if (pidMode != kEl) { - if (std::fabs(track.tpcNSigmaEl()) < cfgVetoIdCut.cfgVetoId07ElTPC) + if (std::fabs(track.tpcNSigmaEl()) < cfgVetoIdCut.cfgVetoId07ElTPC) { return false; + } } } if (cfgId03DoDeRejection) { if (pidMode != kDe) { - if (std::fabs(track.tpcNSigmaDe()) < cfgVetoIdCut.cfgVetoId09DeTPC) + if (std::fabs(track.tpcNSigmaDe()) < cfgVetoIdCut.cfgVetoId09DeTPC) { return false; + } } } return true; @@ -929,27 +937,32 @@ struct NchCumulantsId { bool vetoIdOthersTOF(const T& track) { if (pidMode != kPi) { - if (std::fabs(track.tofNSigmaPi()) < cfgVetoIdCut.cfgVetoId02PiTOF) + if (std::fabs(track.tofNSigmaPi()) < cfgVetoIdCut.cfgVetoId02PiTOF) { return false; + } } if (pidMode != kKa) { - if (std::fabs(track.tofNSigmaKa()) < cfgVetoIdCut.cfgVetoId04KaTOF) + if (std::fabs(track.tofNSigmaKa()) < cfgVetoIdCut.cfgVetoId04KaTOF) { return false; + } } if (pidMode != kPr) { - if (std::fabs(track.tofNSigmaPr()) < cfgVetoIdCut.cfgVetoId06PrTOF) + if (std::fabs(track.tofNSigmaPr()) < cfgVetoIdCut.cfgVetoId06PrTOF) { return false; + } } if (cfgId02DoElRejection) { if (pidMode != kEl) { - if (std::fabs(track.tofNSigmaEl()) < cfgVetoIdCut.cfgVetoId08ElTOF) + if (std::fabs(track.tofNSigmaEl()) < cfgVetoIdCut.cfgVetoId08ElTOF) { return false; + } } } if (cfgId03DoDeRejection) { if (pidMode != kDe) { - if (std::fabs(track.tofNSigmaDe()) < cfgVetoIdCut.cfgVetoId10DeTOF) + if (std::fabs(track.tofNSigmaDe()) < cfgVetoIdCut.cfgVetoId10DeTOF) { return false; + } } } return true; @@ -959,27 +972,32 @@ struct NchCumulantsId { bool vetoIdOthersTPCTOF(const T& track) { if (pidMode != kPi) { - if (std::fabs(track.tpcNSigmaPi()) < cfgVetoIdCut.cfgVetoId01PiTPC && std::fabs(track.tofNSigmaPi()) < cfgVetoIdCut.cfgVetoId02PiTOF) + if (std::fabs(track.tpcNSigmaPi()) < cfgVetoIdCut.cfgVetoId01PiTPC && std::fabs(track.tofNSigmaPi()) < cfgVetoIdCut.cfgVetoId02PiTOF) { return false; + } } if (pidMode != kKa) { - if (std::fabs(track.tpcNSigmaKa()) < cfgVetoIdCut.cfgVetoId03KaTPC && std::fabs(track.tofNSigmaKa()) < cfgVetoIdCut.cfgVetoId04KaTOF) + if (std::fabs(track.tpcNSigmaKa()) < cfgVetoIdCut.cfgVetoId03KaTPC && std::fabs(track.tofNSigmaKa()) < cfgVetoIdCut.cfgVetoId04KaTOF) { return false; + } } if (pidMode != kPr) { - if (std::fabs(track.tpcNSigmaPr()) < cfgVetoIdCut.cfgVetoId05PrTPC && std::fabs(track.tofNSigmaPr()) < cfgVetoIdCut.cfgVetoId06PrTOF) + if (std::fabs(track.tpcNSigmaPr()) < cfgVetoIdCut.cfgVetoId05PrTPC && std::fabs(track.tofNSigmaPr()) < cfgVetoIdCut.cfgVetoId06PrTOF) { return false; + } } if (cfgId02DoElRejection) { if (pidMode != kEl) { - if (std::fabs(track.tpcNSigmaEl()) < cfgVetoIdCut.cfgVetoId07ElTPC && std::fabs(track.tofNSigmaEl()) < cfgVetoIdCut.cfgVetoId08ElTOF) + if (std::fabs(track.tpcNSigmaEl()) < cfgVetoIdCut.cfgVetoId07ElTPC && std::fabs(track.tofNSigmaEl()) < cfgVetoIdCut.cfgVetoId08ElTOF) { return false; + } } } if (cfgId03DoDeRejection) { if (pidMode != kDe) { - if (std::fabs(track.tpcNSigmaDe()) < cfgVetoIdCut.cfgVetoId09DeTPC && std::fabs(track.tofNSigmaDe()) < cfgVetoIdCut.cfgVetoId10DeTOF) + if (std::fabs(track.tpcNSigmaDe()) < cfgVetoIdCut.cfgVetoId09DeTPC && std::fabs(track.tofNSigmaDe()) < cfgVetoIdCut.cfgVetoId10DeTOF) { return false; + } } } return true; @@ -1019,16 +1037,19 @@ struct NchCumulantsId { { switch (pidMode) { case kPi: - if (std::pow(track.tpcNSigmaPi() / nSigmaTPC, 2) + std::pow(track.tofNSigmaPi() / nSigmaTOF, 2) < 1.0) + if (std::pow(track.tpcNSigmaPi() / nSigmaTPC, 2) + std::pow(track.tofNSigmaPi() / nSigmaTOF, 2) < 1.0) { return true; + } break; case kKa: - if (std::pow(track.tpcNSigmaKa() / nSigmaTPC, 2) + std::pow(track.tofNSigmaKa() / nSigmaTOF, 2) < 1.0) + if (std::pow(track.tpcNSigmaKa() / nSigmaTPC, 2) + std::pow(track.tofNSigmaKa() / nSigmaTOF, 2) < 1.0) { return true; + } break; case kPr: - if (std::pow(track.tpcNSigmaPr() / nSigmaTPC, 2) + std::pow(track.tofNSigmaPr() / nSigmaTOF, 2) < 1.0) + if (std::pow(track.tpcNSigmaPr() / nSigmaTPC, 2) + std::pow(track.tofNSigmaPr() / nSigmaTOF, 2) < 1.0) { return true; + } break; default: return false; @@ -1042,16 +1063,19 @@ struct NchCumulantsId { { switch (pidMode) { case kPi: - if (std::pow(track.tpcNSigmaPi(), 2) + std::pow(track.tofNSigmaPi(), 2) < nSigmaSquaredRad) + if (std::pow(track.tpcNSigmaPi(), 2) + std::pow(track.tofNSigmaPi(), 2) < nSigmaSquaredRad) { return true; + } break; case kKa: - if (std::pow(track.tpcNSigmaKa(), 2) + std::pow(track.tofNSigmaKa(), 2) < nSigmaSquaredRad) + if (std::pow(track.tpcNSigmaKa(), 2) + std::pow(track.tofNSigmaKa(), 2) < nSigmaSquaredRad) { return true; + } break; case kPr: - if (std::pow(track.tpcNSigmaPr(), 2) + std::pow(track.tofNSigmaPr(), 2) < nSigmaSquaredRad) + if (std::pow(track.tpcNSigmaPr(), 2) + std::pow(track.tofNSigmaPr(), 2) < nSigmaSquaredRad) { return true; + } break; default: return false; @@ -1063,29 +1087,30 @@ struct NchCumulantsId { template bool checkReliableTOF(const T& track) { - if (track.hasTOF()) - return true; // which check makes the information of TOF relaiable? should track.beta() be checked? - else - return false; + return track.hasTOF(); // which check makes the information of TOF relaiable? should track.beta() be checked? } template bool idTPC(const T& track, const float& nSigmaTPC) { - if (cfgId01CheckVetoCut && !vetoIdOthersTPC(track)) + if (cfgId01CheckVetoCut && !vetoIdOthersTPC(track)) { return false; + } switch (pidMode) { case kPi: - if (std::fabs(track.tpcNSigmaPi()) < nSigmaTPC) + if (std::fabs(track.tpcNSigmaPi()) < nSigmaTPC) { return true; + } break; case kKa: - if (std::fabs(track.tpcNSigmaKa()) < nSigmaTPC) + if (std::fabs(track.tpcNSigmaKa()) < nSigmaTPC) { return true; + } break; case kPr: - if (std::fabs(track.tpcNSigmaPr()) < nSigmaTPC) + if (std::fabs(track.tpcNSigmaPr()) < nSigmaTPC) { return true; + } break; default: return false; @@ -1144,21 +1169,23 @@ struct NchCumulantsId { int getMCTag(const T& track) { int mcTag = 0; - if (!track.has_mcParticle()) + if (!track.has_mcParticle()) { return mcTag; + } auto mcPart = track.mcParticle(); int pdgCode = std::abs(mcPart.pdgCode()); - if (pdgCode == kPiPlus || pdgCode == kPiMinus) - BITSET(mcTag, MC_BIT_PI); - else if (pdgCode == kKPlus || pdgCode == kKMinus) - BITSET(mcTag, MC_BIT_KA); - else if (pdgCode == kProton || pdgCode == kProtonBar) - BITSET(mcTag, MC_BIT_PR); - else if (pdgCode == kElectron || pdgCode == kPositron) - BITSET(mcTag, MC_BIT_EL); - else if (pdgCode == kDeuteron || pdgCode == -kDeuteron) - BITSET(mcTag, MC_BIT_DE); + if (pdgCode == kPiPlus || pdgCode == kPiMinus) { + BITSET(mcTag, McBitPI); + } else if (pdgCode == kKPlus || pdgCode == kKMinus) { + BITSET(mcTag, McBitKA); + } else if (pdgCode == kProton || pdgCode == kProtonBar) { + BITSET(mcTag, McBitPR); + } else if (pdgCode == kElectron || pdgCode == kPositron) { + BITSET(mcTag, McBitEL); + } else if (pdgCode == kDeuteron || pdgCode == -kDeuteron) { + BITSET(mcTag, McBitDE); + } return mcTag; } @@ -1173,9 +1200,8 @@ struct NchCumulantsId { std::fabs(track.tpcNSigmaKa()) > cfgIdKaRejNSigma && std::fabs(track.tpcNSigmaPr()) > cfgIdPrRejNSigma) { return false; - } else { - return true; } + return true; } // Pion @@ -1462,14 +1488,17 @@ struct NchCumulantsId { bool isEventSelected(const CollisionType& coll) { if (cfgEventSelection.cfgEvSel01doNoSameBunchPileup && - !coll.selection_bit(aod::evsel::kNoSameBunchPileup)) + !coll.selection_bit(aod::evsel::kNoSameBunchPileup)) { return false; + } if (cfgEventSelection.cfgEvSel02doIsGoodZvtxFT0vsPV && - !coll.selection_bit(aod::evsel::kIsGoodZvtxFT0vsPV)) + !coll.selection_bit(aod::evsel::kIsGoodZvtxFT0vsPV)) { return false; + } if (cfgEventSelection.cfgEvSel03doIsGoodITSLayersAll && - !coll.selection_bit(aod::evsel::kIsGoodITSLayersAll)) + !coll.selection_bit(aod::evsel::kIsGoodITSLayersAll)) { return false; + } return true; } @@ -1521,8 +1550,9 @@ struct NchCumulantsId { if constexpr (analysisType == doDataProcessing) { for (const auto& col : collisions) { - if (!isEventSelected(col)) + if (!isEventSelected(col)) { continue; + } nP = 0; nM = 0; nCh = 0; @@ -1557,7 +1587,7 @@ struct NchCumulantsId { continue; } - int idMethod; + int idMethod = -1; // pion if (selPion(track, idMethod)) { if (track.sign() == 1) { @@ -1567,10 +1597,12 @@ struct NchCumulantsId { nAPi += hPtEtaForEffCorrection[kPi][kNeg]->GetBinContent(ptEtaBin); } - if (idMethod == kTPCidentified) + if (idMethod == kTPCidentified) { fillIdentificationQA(hist, track); - if (idMethod == kTPCTOFidentified) + } + if (idMethod == kTPCTOFidentified) { fillIdentificationQA(hist, track); + } } // kaon if (selKaon(track, idMethod)) { @@ -1581,10 +1613,12 @@ struct NchCumulantsId { nAKa += hPtEtaForEffCorrection[kKa][kNeg]->GetBinContent(ptEtaBin); } - if (idMethod == kTPCidentified) + if (idMethod == kTPCidentified) { fillIdentificationQA(hist, track); - if (idMethod == kTPCTOFidentified) + } + if (idMethod == kTPCTOFidentified) { fillIdentificationQA(hist, track); + } } // proton if (selProton(track, idMethod)) { @@ -1595,10 +1629,12 @@ struct NchCumulantsId { nAPr += hPtEtaForEffCorrection[kPr][kNeg]->GetBinContent(ptEtaBin); } - if (idMethod == kTPCidentified) + if (idMethod == kTPCidentified) { fillIdentificationQA(hist, track); - if (idMethod == kTPCTOFidentified) + } + if (idMethod == kTPCTOFidentified) { fillIdentificationQA(hist, track); + } } } // track loop ends nCh = nP - nM; @@ -1618,8 +1654,9 @@ struct NchCumulantsId { LOG(warning) << "No MC collision for this collision, skip..."; continue; } - if (!isEventSelected(col)) + if (!isEventSelected(col)) { continue; + } nP = 0; nM = 0; nCh = 0; @@ -1665,24 +1702,27 @@ struct NchCumulantsId { if (selPion(track, idMethodPi)) { trackIsPion = true; - BITSET(trackIdTag, ID_BIT_PI); + BITSET(trackIdTag, IdBitPI); hist.fill(HIST("PIDValidation/tpcSparse_Pi"), track.p(), track.tpcNSigmaPi(), trackIdTag, mcTag); - if (track.hasTOF()) + if (track.hasTOF()) { hist.fill(HIST("PIDValidation/tofSparse_Pi"), track.p(), track.tofNSigmaPi(), trackIdTag, mcTag); + } } if (selKaon(track, idMethodKa)) { trackIsKaon = true; - BITSET(trackIdTag, ID_BIT_KA); + BITSET(trackIdTag, IdBitKA); hist.fill(HIST("PIDValidation/tpcSparse_Ka"), track.p(), track.tpcNSigmaKa(), trackIdTag, mcTag); - if (track.hasTOF()) + if (track.hasTOF()) { hist.fill(HIST("PIDValidation/tofSparse_Ka"), track.p(), track.tofNSigmaKa(), trackIdTag, mcTag); + } } if (selProton(track, idMethodPr)) { trackIsProton = true; - BITSET(trackIdTag, ID_BIT_PR); + BITSET(trackIdTag, IdBitPR); hist.fill(HIST("PIDValidation/tpcSparse_Pr"), track.p(), track.tpcNSigmaPr(), trackIdTag, mcTag); - if (track.hasTOF()) + if (track.hasTOF()) { hist.fill(HIST("PIDValidation/tofSparse_Pr"), track.p(), track.tofNSigmaPr(), trackIdTag, mcTag); + } } if constexpr (analysisType == doPurityProcessing) { @@ -1763,8 +1803,9 @@ struct NchCumulantsId { LOG(warning) << "No MC collision for this event, skip..."; continue; } - if (!isEventSelected(col)) + if (!isEventSelected(col)) { continue; + } const auto& mcColl = col.mcCollision(); // ---- apply same Vz cut as data/reco ---- @@ -1857,8 +1898,9 @@ struct NchCumulantsId { void processSim(MyFilteredColsWithMcLabels const& collisions, MyFilteredTracksWithMcLabels const& tracks, aod::McCollisions const& mcCollisions, aod::McParticles const& mcParticles) { - if (cfgEventSelection.flagUnusedVariableError) + if (cfgEventSelection.flagUnusedVariableError) { LOG(info) << mcCollisions.size(); + } bool trackIsPion = false; bool trackIsKaon = false; bool trackIsProton = false; @@ -1876,12 +1918,14 @@ struct NchCumulantsId { LOG(warning) << "No MC collision for this collision, skip..."; continue; } - if (!isEventSelected(col)) + if (!isEventSelected(col)) { continue; + } auto mcCollision = col.mcCollision(); - if (cfgEventSelection.checkCollPosZMc && std::abs(mcCollision.posZ()) > cfgCutPosZ) + if (cfgEventSelection.checkCollPosZMc && std::abs(mcCollision.posZ()) > cfgCutPosZ) { continue; + } // slice reco tracks to this collision const uint64_t collIdx = col.globalIndex(); @@ -1910,12 +1954,14 @@ struct NchCumulantsId { float nPiGen = 0, nAPiGen = 0; for (const auto& mcTrack : mcTracksTablePerMcColl) { - if (!mcTrack.isPhysicalPrimary()) + if (!mcTrack.isPhysicalPrimary()) { continue; + } if (mcTrack.pt() <= cfgCutPtMin || mcTrack.pt() >= cfgCutPtMax || - std::abs(mcTrack.eta()) >= cfgCutEta) + std::abs(mcTrack.eta()) >= cfgCutEta) { continue; + } int pdg = mcTrack.pdgCode(); @@ -2041,7 +2087,7 @@ struct NchCumulantsId { // Fill separate spares for each species if it passes the cut if (selPion(track, idMethodPi)) { trackIsPion = true; - BITSET(trackIdTag, ID_BIT_PI); + BITSET(trackIdTag, IdBitPI); // Fill TPC sparse for pion hist.fill(HIST("PIDValidation/tpcSparse_Pi"), track.p(), track.tpcNSigmaPi(), trackIdTag, mcTag); // Fill TOF sparse for pion if has TOF @@ -2051,7 +2097,7 @@ struct NchCumulantsId { } if (selKaon(track, idMethodKa)) { trackIsKaon = true; - BITSET(trackIdTag, ID_BIT_KA); + BITSET(trackIdTag, IdBitKA); // Fill TPC sparse for kaon hist.fill(HIST("PIDValidation/tpcSparse_Ka"), track.p(), track.tpcNSigmaKa(), trackIdTag, mcTag); // Fill TOF sparse for kaon if has TOF @@ -2061,7 +2107,7 @@ struct NchCumulantsId { } if (selProton(track, idMethodPr)) { trackIsProton = true; - BITSET(trackIdTag, ID_BIT_PR); + BITSET(trackIdTag, IdBitPR); // Fill TPC sparse for proton hist.fill(HIST("PIDValidation/tpcSparse_Pr"), track.p(), track.tpcNSigmaPr(), trackIdTag, mcTag); // Fill TOF sparse for proton if has TOF @@ -2119,10 +2165,12 @@ struct NchCumulantsId { fillEffPower(pimPow, weight); } // PID band QA for pions - if (idMethodPi == kTPCidentified) + if (idMethodPi == kTPCidentified) { fillIdentificationQA(hist, track); - if (idMethodPi == kTPCTOFidentified) + } + if (idMethodPi == kTPCTOFidentified) { fillIdentificationQA(hist, track); + } } else if (trackIsKaon) { if (track.sign() > 0) { nKaRec += hPtEtaForEffCorrection[kKa][kPos]->GetBinContent(ptEtaBin); @@ -2136,10 +2184,12 @@ struct NchCumulantsId { fillEffPower(kamPow, weight); } // PID band QA for kaons - if (idMethodKa == kTPCidentified) + if (idMethodKa == kTPCidentified) { fillIdentificationQA(hist, track); - if (idMethodKa == kTPCTOFidentified) + } + if (idMethodKa == kTPCTOFidentified) { fillIdentificationQA(hist, track); + } } else if (trackIsProton) { if (track.sign() > 0) { nPrRec += hPtEtaForEffCorrection[kPr][kPos]->GetBinContent(ptEtaBin); @@ -2153,10 +2203,12 @@ struct NchCumulantsId { fillEffPower(aprPow, weight); } // PID band QA for protons - if (idMethodPr == kTPCidentified) + if (idMethodPr == kTPCidentified) { fillIdentificationQA(hist, track); - if (idMethodPr == kTPCTOFidentified) + } + if (idMethodPr == kTPCTOFidentified) { fillIdentificationQA(hist, track); + } } // purity check - check pdg aginst sign bool purityPion = false; @@ -2164,20 +2216,26 @@ struct NchCumulantsId { bool purityProton = false; if (trackIsPion) { - if (track.sign() > 0 && pdg == kPiPlus) + if (track.sign() > 0 && pdg == kPiPlus) { purityPion = true; - if (track.sign() < 0 && pdg == kPiMinus) + } + if (track.sign() < 0 && pdg == kPiMinus) { purityPion = true; + } } else if (trackIsKaon) { - if (track.sign() > 0 && pdg == kKPlus) + if (track.sign() > 0 && pdg == kKPlus) { purityKaon = true; - if (track.sign() < 0 && pdg == kKMinus) + } + if (track.sign() < 0 && pdg == kKMinus) { purityKaon = true; + } } else if (trackIsProton) { - if (track.sign() > 0 && pdg == kProton) + if (track.sign() > 0 && pdg == kProton) { purityProton = true; - if (track.sign() < 0 && pdg == kProtonBar) + } + if (track.sign() < 0 && pdg == kProtonBar) { purityProton = true; + } } // charge purity — track.sign() + isKnownCharged + PDG sign consistency From 5ceacf7fe0d22e4116cf9cdbf6d1f09c6ea6442b Mon Sep 17 00:00:00 2001 From: Pravata Panigrahi <158029188+pravatp98@users.noreply.github.com> Date: Mon, 17 Aug 2026 16:54:55 +0530 Subject: [PATCH 05/13] Comment out unused IdBitEL and IdBitDE Commented out unused identification bits for clarity. --- PWGCF/EbyEFluctuations/Tasks/nchCumulantsId.cxx | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/PWGCF/EbyEFluctuations/Tasks/nchCumulantsId.cxx b/PWGCF/EbyEFluctuations/Tasks/nchCumulantsId.cxx index 11c5f42df95..05ffe3c2448 100644 --- a/PWGCF/EbyEFluctuations/Tasks/nchCumulantsId.cxx +++ b/PWGCF/EbyEFluctuations/Tasks/nchCumulantsId.cxx @@ -59,8 +59,8 @@ using namespace std; constexpr int IdBitPI = 0; // Identificationi bits for PID checks constexpr int IdBitKA = 1; constexpr int IdBitPR = 2; -constexpr int IdBitEL = 3; -constexpr int IdBitDE = 4; +// constexpr int IdBitEL = 3; +// constexpr int IdBitDE = 4; constexpr int McBitPI = 0; // MC particle identification bits constexpr int McBitKA = 1; From 2917414cf40b574b7a99ede3f643797d415af8da Mon Sep 17 00:00:00 2001 From: Pravata Panigrahi <158029188+pravatp98@users.noreply.github.com> Date: Mon, 17 Aug 2026 18:21:48 +0530 Subject: [PATCH 06/13] Refactor conditional statements and parameter types --- PWGCF/EbyEFluctuations/Tasks/nchCumulantsId.cxx | 15 +++++++++------ 1 file changed, 9 insertions(+), 6 deletions(-) diff --git a/PWGCF/EbyEFluctuations/Tasks/nchCumulantsId.cxx b/PWGCF/EbyEFluctuations/Tasks/nchCumulantsId.cxx index 05ffe3c2448..027f38957c9 100644 --- a/PWGCF/EbyEFluctuations/Tasks/nchCumulantsId.cxx +++ b/PWGCF/EbyEFluctuations/Tasks/nchCumulantsId.cxx @@ -25,13 +25,13 @@ #include // #include #include -#include +// #include #include #include #include #include #include -#include +// #include #include #include @@ -1122,8 +1122,9 @@ struct NchCumulantsId { template bool idTPCTOF(const T& track, const int& pidCutType, const float& nSigmaTPC, const float& nSigmaTOF, const float& nSigmaSquaredRad) { - if (cfgId01CheckVetoCut && !vetoIdOthersTPCTOF(track)) + if (cfgId01CheckVetoCut && !vetoIdOthersTPCTOF(track)) { return false; + } if (pidCutType == kRectangularCut) { return selIdRectangularCut(track, nSigmaTPC, nSigmaTOF); } else if (pidCutType == kCircularCut) { @@ -1283,7 +1284,7 @@ struct NchCumulantsId { } template - void fillCollQA(const T& coll, const int& nCh, const int& nT) + void fillCollQA(const T& coll, const float& nCh, const float& nT) { hist.fill(HIST(HistRegDire[mode]) + HIST("h_VtxZ"), coll.posZ()); hist.fill(HIST(HistRegDire[mode]) + HIST("h_Counts"), 0.5); @@ -1395,8 +1396,9 @@ struct NchCumulantsId { const int& idMethodKa, const bool& trackIsKaon, float& nAKa, float& nKa, const int& idMethodPr, const bool& trackIsProton, float& nPr, float& nAPr, H& hReg) { - if (cfgEventSelection.flagUnusedVariableError) + if (cfgEventSelection.flagUnusedVariableError) { LOG(info) << trackIdTag << idMethodPi << ":" << idMethodKa << ":" << idMethodPr; + } if (track.sign() > 0) { // fillRecoTrackQA(hReg, track); nP++; @@ -2063,8 +2065,9 @@ struct NchCumulantsId { continue; } auto mcPart = track.mcParticle(); - if (!mcPart.isPhysicalPrimary()) + if (!mcPart.isPhysicalPrimary()) { continue; + } int pdg = mcPart.pdgCode(); int mcTag = getMCTag(track); From 8f848642e29783ca3424a7af9d60c78de57c0c9c Mon Sep 17 00:00:00 2001 From: Pravata Panigrahi <158029188+pravatp98@users.noreply.github.com> Date: Tue, 8 Sep 2026 07:38:45 +0530 Subject: [PATCH 07/13] Add configurable pT cutoffs for kaon and proton Added configurable pT cutoffs for kaons and protons, and updated relevant sections to use these cutoffs for efficiency corrections. --- .../EbyEFluctuations/Tasks/nchCumulantsId.cxx | 56 +++++++++++++------ 1 file changed, 39 insertions(+), 17 deletions(-) diff --git a/PWGCF/EbyEFluctuations/Tasks/nchCumulantsId.cxx b/PWGCF/EbyEFluctuations/Tasks/nchCumulantsId.cxx index 027f38957c9..a30b976e07e 100644 --- a/PWGCF/EbyEFluctuations/Tasks/nchCumulantsId.cxx +++ b/PWGCF/EbyEFluctuations/Tasks/nchCumulantsId.cxx @@ -212,6 +212,10 @@ struct NchCumulantsId { ConfigurableAxis subSampleAxis{"subSampleAxis", {10, 0., 10.}, "Subsample"}; TRandom3* fRandom = new TRandom3(0); // Random number generator for subsampling int currentSubsample = 0; + Configurable cfgCutPtMinKa{"cfgCutPtMinKa", 0.30, "min pT for kaon eff-correction inclusion"}; + Configurable cfgCutPtMaxKa{"cfgCutPtMaxKa", 1.20, "max pT for kaon eff-correction inclusion (purity cutoff)"}; + Configurable cfgCutPtMinPr{"cfgCutPtMinPr", 0.40, "min pT for proton eff-correction inclusion"}; + Configurable cfgCutPtMaxPr{"cfgCutPtMaxPr", 1.50, "max pT for proton eff-correction inclusion (purity cutoff)"}; } cfgEventSelection; // Configurables for particle Identification @@ -297,9 +301,11 @@ struct NchCumulantsId { double q4 = 0.; }; + static constexpr float KmaxWeight = 10.0f; // secondary safeguard for pr and kaon eff fix + inline void fillEffPower(EffPowerSums& p, float weight) { - if (weight <= 0.f) { + if (weight <= 0.f || weight > KmaxWeight) { return; } @@ -1999,19 +2005,27 @@ struct NchCumulantsId { fillGenTrackQA(genAnalysis, mcTrack); } else if (pdg == kKPlus) { // fillGenTrackQA(genAnalysis, mcTrack); - nKaGen++; + if (mcTrack.pt() >= cfgEventSelection.cfgCutPtMinKa && mcTrack.pt() <= cfgEventSelection.cfgCutPtMaxKa) { + nKaGen++; + } fillGenTrackQA(genAnalysis, mcTrack); } else if (pdg == kKMinus) { // fillGenTrackQA(genAnalysis, mcTrack); - nAKaGen++; + if (mcTrack.pt() >= cfgEventSelection.cfgCutPtMinKa && mcTrack.pt() <= cfgEventSelection.cfgCutPtMaxKa) { + nAKaGen++; + } fillGenTrackQA(genAnalysis, mcTrack); } else if (pdg == kProton) { // fillGenTrackQA(genAnalysis, mcTrack); - nPrGen++; + if (mcTrack.pt() >= cfgEventSelection.cfgCutPtMinPr && mcTrack.pt() <= cfgEventSelection.cfgCutPtMaxPr) { + nPrGen++; + } fillGenTrackQA(genAnalysis, mcTrack); } else if (pdg == kProtonBar) { // fillGenTrackQA(genAnalysis, mcTrack); - nAPrGen++; + if (mcTrack.pt() >= cfgEventSelection.cfgCutPtMinPr && mcTrack.pt() <= cfgEventSelection.cfgCutPtMaxPr) { + nAPrGen++; + } fillGenTrackQA(genAnalysis, mcTrack); } @@ -2176,15 +2190,19 @@ struct NchCumulantsId { } } else if (trackIsKaon) { if (track.sign() > 0) { - nKaRec += hPtEtaForEffCorrection[kKa][kPos]->GetBinContent(ptEtaBin); fillRecoTrackQA(recoAnalysis, track); - float weight = hPtEtaForEffCorrection[kKa][kPos]->GetBinContent(ptEtaBin); - fillEffPower(kapPow, weight); + if (track.pt() >= cfgEventSelection.cfgCutPtMinKa && track.pt() <= cfgEventSelection.cfgCutPtMaxKa) { + nKaRec += hPtEtaForEffCorrection[kKa][kPos]->GetBinContent(ptEtaBin); + float weight = hPtEtaForEffCorrection[kKa][kPos]->GetBinContent(ptEtaBin); + fillEffPower(kapPow, weight); + } } else if (track.sign() < 0) { - nAKaRec += hPtEtaForEffCorrection[kKa][kNeg]->GetBinContent(ptEtaBin); fillRecoTrackQA(recoAnalysis, track); - float weight = hPtEtaForEffCorrection[kKa][kNeg]->GetBinContent(ptEtaBin); - fillEffPower(kamPow, weight); + if (track.pt() >= cfgEventSelection.cfgCutPtMinKa && track.pt() <= cfgEventSelection.cfgCutPtMaxKa) { + nAKaRec += hPtEtaForEffCorrection[kKa][kNeg]->GetBinContent(ptEtaBin); + float weight = hPtEtaForEffCorrection[kKa][kNeg]->GetBinContent(ptEtaBin); + fillEffPower(kamPow, weight); + } } // PID band QA for kaons if (idMethodKa == kTPCidentified) { @@ -2195,15 +2213,19 @@ struct NchCumulantsId { } } else if (trackIsProton) { if (track.sign() > 0) { - nPrRec += hPtEtaForEffCorrection[kPr][kPos]->GetBinContent(ptEtaBin); fillRecoTrackQA(recoAnalysis, track); - float weight = hPtEtaForEffCorrection[kPr][kPos]->GetBinContent(ptEtaBin); - fillEffPower(prPow, weight); + if (track.pt() >= cfgEventSelection.cfgCutPtMinPr && track.pt() <= cfgEventSelection.cfgCutPtMaxPr) { + nPrRec += hPtEtaForEffCorrection[kPr][kPos]->GetBinContent(ptEtaBin); + float weight = hPtEtaForEffCorrection[kPr][kPos]->GetBinContent(ptEtaBin); + fillEffPower(prPow, weight); + } } else if (track.sign() < 0) { - nAPrRec += hPtEtaForEffCorrection[kPr][kNeg]->GetBinContent(ptEtaBin); fillRecoTrackQA(recoAnalysis, track); - float weight = hPtEtaForEffCorrection[kPr][kNeg]->GetBinContent(ptEtaBin); - fillEffPower(aprPow, weight); + if (track.pt() >= cfgEventSelection.cfgCutPtMinPr && track.pt() <= cfgEventSelection.cfgCutPtMaxPr) { + nAPrRec += hPtEtaForEffCorrection[kPr][kNeg]->GetBinContent(ptEtaBin); + float weight = hPtEtaForEffCorrection[kPr][kNeg]->GetBinContent(ptEtaBin); + fillEffPower(aprPow, weight); + } } // PID band QA for protons if (idMethodPr == kTPCidentified) { From cf08cd338098adffe38ccd554ff7369c2a080007 Mon Sep 17 00:00:00 2001 From: Pravata Panigrahi <158029188+pravatp98@users.noreply.github.com> Date: Tue, 8 Sep 2026 20:45:47 +0530 Subject: [PATCH 08/13] Refactor bit manipulation functions and simplify conditions --- .../EbyEFluctuations/Tasks/nchCumulantsId.cxx | 74 +++++++++---------- 1 file changed, 34 insertions(+), 40 deletions(-) diff --git a/PWGCF/EbyEFluctuations/Tasks/nchCumulantsId.cxx b/PWGCF/EbyEFluctuations/Tasks/nchCumulantsId.cxx index a30b976e07e..d19f975b2e3 100644 --- a/PWGCF/EbyEFluctuations/Tasks/nchCumulantsId.cxx +++ b/PWGCF/EbyEFluctuations/Tasks/nchCumulantsId.cxx @@ -68,8 +68,15 @@ constexpr int McBitPR = 2; constexpr int McBitEL = 3; constexpr int McBitDE = 4; -#define BITSET(mask, ithBit) ((mask) |= (1 << (ithBit))) // avoid name bitset as std::bitset is already there -#define BITCHECK(mask, ithBit) ((mask) & (1 << (ithBit))) // bit check will return int value, not bool, use BITCHECK != 0 in Analysi +constexpr void BITSET(int& mask, int ithBit) +{ + mask |= (1 << ithBit); +} + +constexpr bool BITCHECK(int mask, int ithBit) +{ + return (mask & (1 << ithBit)) != 0; +} enum PidEnum { kCh = 0, @@ -1014,28 +1021,20 @@ struct NchCumulantsId { { switch (pidMode) { case kPi: - if (std::fabs(track.tpcNSigmaPi()) < nSigmaTPC && - std::fabs(track.tofNSigmaPi()) < nSigmaTOF) { - return true; - } - break; + return std::fabs(track.tpcNSigmaPi()) < nSigmaTPC && + std::fabs(track.tofNSigmaPi()) < nSigmaTOF; + case kKa: - if (std::fabs(track.tpcNSigmaKa()) < nSigmaTPC && - std::fabs(track.tofNSigmaKa()) < nSigmaTOF) { - return true; - } - break; + return std::fabs(track.tpcNSigmaKa()) < nSigmaTPC && + std::fabs(track.tofNSigmaKa()) < nSigmaTOF; + case kPr: - if (std::fabs(track.tpcNSigmaPr()) < nSigmaTPC && - std::fabs(track.tofNSigmaPr()) < nSigmaTOF) { - return true; - } - break; + return std::fabs(track.tpcNSigmaPr()) < nSigmaTPC && + std::fabs(track.tofNSigmaPr()) < nSigmaTOF; + default: return false; - break; } - return false; } template @@ -1043,25 +1042,23 @@ struct NchCumulantsId { { switch (pidMode) { case kPi: - if (std::pow(track.tpcNSigmaPi() / nSigmaTPC, 2) + std::pow(track.tofNSigmaPi() / nSigmaTOF, 2) < 1.0) { - return true; - } - break; + return std::pow(track.tpcNSigmaPi() / nSigmaTPC, 2) + + std::pow(track.tofNSigmaPi() / nSigmaTOF, 2) < + 1.0; + case kKa: - if (std::pow(track.tpcNSigmaKa() / nSigmaTPC, 2) + std::pow(track.tofNSigmaKa() / nSigmaTOF, 2) < 1.0) { - return true; - } - break; + return std::pow(track.tpcNSigmaKa() / nSigmaTPC, 2) + + std::pow(track.tofNSigmaKa() / nSigmaTOF, 2) < + 1.0; + case kPr: - if (std::pow(track.tpcNSigmaPr() / nSigmaTPC, 2) + std::pow(track.tofNSigmaPr() / nSigmaTOF, 2) < 1.0) { - return true; - } - break; + return std::pow(track.tpcNSigmaPr() / nSigmaTPC, 2) + + std::pow(track.tofNSigmaPr() / nSigmaTOF, 2) < + 1.0; + default: return false; - break; } - return false; } template @@ -1202,13 +1199,10 @@ struct NchCumulantsId { template bool selTrackForId(const T& track) { - if (cfgIdElRejLowNSigma < track.tpcNSigmaEl() && track.tpcNSigmaEl() < cfgIdElRejHighNSigma && - std::fabs(track.tpcNSigmaPi()) > cfgIdPiRejNSigma && - std::fabs(track.tpcNSigmaKa()) > cfgIdKaRejNSigma && - std::fabs(track.tpcNSigmaPr()) > cfgIdPrRejNSigma) { - return false; - } - return true; + return !(cfgIdElRejLowNSigma < track.tpcNSigmaEl() && track.tpcNSigmaEl() < cfgIdElRejHighNSigma && + std::fabs(track.tpcNSigmaPi()) > cfgIdPiRejNSigma && + std::fabs(track.tpcNSigmaKa()) > cfgIdKaRejNSigma && + std::fabs(track.tpcNSigmaPr()) > cfgIdPrRejNSigma); } // Pion From d8d29ebd5350b2e86af62bc5815909c914960158 Mon Sep 17 00:00:00 2001 From: Pravata Panigrahi <158029188+pravatp98@users.noreply.github.com> Date: Tue, 10 Nov 2026 10:55:27 +0530 Subject: [PATCH 09/13] Updated for Particle Identification and Correction --- .../EbyEFluctuations/Tasks/nchCumulantsId.cxx | 1003 ++++++++++++----- 1 file changed, 711 insertions(+), 292 deletions(-) diff --git a/PWGCF/EbyEFluctuations/Tasks/nchCumulantsId.cxx b/PWGCF/EbyEFluctuations/Tasks/nchCumulantsId.cxx index d19f975b2e3..3ea1b590327 100644 --- a/PWGCF/EbyEFluctuations/Tasks/nchCumulantsId.cxx +++ b/PWGCF/EbyEFluctuations/Tasks/nchCumulantsId.cxx @@ -48,6 +48,7 @@ #include #include #include +#include #include using namespace o2; @@ -135,6 +136,16 @@ static constexpr std::array FCGenDir = { "Gen/NetPr/", "Gen/NetCh/"}; +static constexpr std::array CountDir = { + "CountCheck/Pr/", + "CountCheck/APr/", + "CountCheck/PiPos/", + "CountCheck/PiNeg/", + "CountCheck/KaPos/", + "CountCheck/KaNeg/", + "CountCheck/Pos/", + "CountCheck/Neg/"}; + static constexpr std::array PidDire = { "Ch/", "Pi/", @@ -157,6 +168,425 @@ std::string getModifiedStr(const std::string& myString) return myString; } +enum PidSelectionEnum { + kPidPi = 0, + kPidKa, + kPidPr, + kPidEl, + kPidMu, + kPidDe, + kNPidSelection +}; + +enum CutSettingEnum { + kThrPforTOF = 0, + kIdCutTypeLowP, + kNSigmaTPCLowP, + kNSigmaTOFLowP, + kNSigmaRadLowP, + kIdCutTypeHighP, + kNSigmaTPCHighP, + kNSigmaTOFHighP, + kNSigmaRadHighP, + kDoVetoOthers, + kDoRelativeTPCcheck, + kDoRelativeTOFcheck, + kDoRelativeTPCTOFcheck, + kNCutSettings +}; + +enum VetoSettingEnum { + kDoVetoTPC = 0, + kDoVetoTOF, + kVetoTPC, + kVetoTOF, + kNVetoSettings +}; + +static constexpr std::array, kNPidSelection> DefaultPIDcheckValues = {{ + {0.7, 0, 3.0, 3.0, 9.0, 0, 3.0, 3.0, 9.0, 0, 1, 0, 1}, // Pi + {0.5, 0, 3.0, 3.0, 9.0, 0, 3.0, 3.0, 9.0, 0, 1, 0, 1}, // Ka + {0.8, 0, 3.0, 3.0, 9.0, 0, 3.0, 3.0, 9.0, 0, 1, 0, 1}, // Pr + {0.4, 0, 3.0, 3.0, 9.0, 0, 3.0, 3.0, 9.0, 0, 1, 0, 1}, // El + {0.4, 0, 3.0, 3.0, 9.0, 0, 3.0, 3.0, 9.0, 0, 1, 0, 1}, // Mu + {0.4, 0, 3.0, 3.0, 9.0, 0, 3.0, 3.0, 9.0, 0, 1, 0, 1} // De +}}; + +static constexpr std::array, kNPidSelection> DefaultPidVetoValues = {{ + {1, 1, 3.0, 3.0}, // Pi + {1, 1, 3.0, 3.0}, // Ka + {1, 1, 3.0, 3.0}, // Pr + {0, 0, 3.0, 3.0}, // El + {0, 0, 3.0, 3.0}, // Mu + {0, 0, 3.0, 3.0} // De +}}; + +template +inline T getCfg(const auto& cfgLabelledArray, int32_t row, int32_t col) +{ + double val = cfgLabelledArray->get(row, col); // Always stored as double + + if constexpr (std::is_same_v) { + return val != 0.0; // will return anything nonzero as true + } else if constexpr (std::is_same_v) { + // cppcheck-suppress suspiciousFloatingPointCast + return static_cast(val); + } else if constexpr (std::is_same_v) { + return val; + } else if constexpr (std::is_same_v) { + return static_cast(val); + } else { + static_assert( + std::is_same_v || std::is_same_v || std::is_same_v || std::is_same_v, + "getCfg() only supports T = bool, float, double or int"); + } +} + +template +bool relativeIdOthersTPC(const T& track) +{ + float distSpecies = 1.e6f; + + if constexpr (pidMode == kPidPi) { + distSpecies = std::fabs(track.tpcNSigmaPi()); + } else if constexpr (pidMode == kPidKa) { + distSpecies = std::fabs(track.tpcNSigmaKa()); + } else if constexpr (pidMode == kPidPr) { + distSpecies = std::fabs(track.tpcNSigmaPr()); + } else { + // FIX ME! : Relative check will be updated in future, currently implemented for pion, kaon and proton only + // else if constexpr (pidMode == kPidEl) { + // distSpecies = std::fabs(track.tpcNSigmaEl()); + // } else if constexpr (pidMode == kPidMu) { + // distSpecies = std::fabs(track.tpcNSigmaMu()); + // } else if constexpr (pidMode == kPidDe) { + // distSpecies = std::fabs(track.tpcNSigmaDe()); + // } + + return true; + } + + if constexpr (pidMode != kPidPi) { + if (std::fabs(track.tpcNSigmaPi()) <= distSpecies) { + return false; + } + } + + if constexpr (pidMode != kPidKa) { + if (std::fabs(track.tpcNSigmaKa()) <= distSpecies) { + return false; + } + } + + if constexpr (pidMode != kPidPr) { + if (std::fabs(track.tpcNSigmaPr()) <= distSpecies) { + return false; + } + } + + // FIX ME! : Relative check will be updated in future, currently implemented for pion, kaon and proton only + // if constexpr (pidMode != kPidEl) + // if (std::fabs(track.tpcNSigmaEl()) <= distSpecies) + // return false; + + // if constexpr (pidMode != kPidMu) + // if (std::fabs(track.tpcNSigmaMu()) <= distSpecies) + // return false; + + // if constexpr (pidMode != kPidDe) + // if (std::fabs(track.tpcNSigmaDe()) <= distSpecies) + // return false; + + return true; +} + +template +bool relativeIdOthersTOF(const T& track) +{ + float distSpecies = 1.e6f; + + if constexpr (pidMode == kPidPi) { + distSpecies = std::fabs(track.tofNSigmaPi()); + } else if constexpr (pidMode == kPidKa) { + distSpecies = std::fabs(track.tofNSigmaKa()); + } else if constexpr (pidMode == kPidPr) { + distSpecies = std::fabs(track.tofNSigmaPr()); + + // FIX ME! : Relative check will be updated in future, currently implemented for pion, kaon and proton only + // } else if constexpr (pidMode == kPidEl) { + // distSpecies = std::fabs(track.tofNSigmaEl()); + // } else if constexpr (pidMode == kPidMu) { + // distSpecies = std::fabs(track.tofNSigmaMu()); + // } else if constexpr (pidMode == kPidDe) { + // distSpecies = std::fabs(track.tofNSigmaDe()); + } else { + return true; + } + + if constexpr (pidMode != kPidPi) { + if (std::fabs(track.tofNSigmaPi()) <= distSpecies) { + return false; + } + } + + if constexpr (pidMode != kPidKa) { + if (std::fabs(track.tofNSigmaKa()) <= distSpecies) { + return false; + } + } + + if constexpr (pidMode != kPidPr) { + if (std::fabs(track.tofNSigmaPr()) <= distSpecies) { + return false; + } + } + // FIX ME! : Relative check will be updated in future, currently implemented for pion, kaon and proton only + // if constexpr (pidMode != kPidEl) + // if (std::fabs(track.tofNSigmaEl()) <= distSpecies) + // return false; + + // if constexpr (pidMode != kPidMu) + // if (std::fabs(track.tofNSigmaMu()) <= distSpecies) + // return false; + + // if constexpr (pidMode != kPidDe) + // if (std::fabs(track.tofNSigmaDe()) <= distSpecies) + // return false; + + return true; +} + +template +bool relativeIdOthersTPCTOF(const T& track) +{ + float distSpeciesSq = 1.e6f; + + if constexpr (pidMode == kPidPi) { + distSpeciesSq = track.tpcNSigmaPi() * track.tpcNSigmaPi() + track.tofNSigmaPi() * track.tofNSigmaPi(); + } else if constexpr (pidMode == kPidKa) { + distSpeciesSq = track.tpcNSigmaKa() * track.tpcNSigmaKa() + track.tofNSigmaKa() * track.tofNSigmaKa(); + } else if constexpr (pidMode == kPidPr) { + distSpeciesSq = track.tpcNSigmaPr() * track.tpcNSigmaPr() + track.tofNSigmaPr() * track.tofNSigmaPr(); + // FIX ME! : Relative check will be updated in future, currently implemented for pion, kaon and proton only + // } else if constexpr (pidMode == kPidEl) { + // distSpeciesSq = track.tpcNSigmaEl() * track.tpcNSigmaEl() + track.tofNSigmaEl() * track.tofNSigmaEl(); + // } else if constexpr (pidMode == kPidMu) { + // distSpeciesSq = track.tpcNSigmaMu() * track.tpcNSigmaMu() + track.tofNSigmaMu() * track.tofNSigmaMu(); + // } else if constexpr (pidMode == kPidDe) { + // distSpeciesSq = track.tpcNSigmaDe() * track.tpcNSigmaDe() + track.tofNSigmaDe() * track.tofNSigmaDe(); + } else { + return true; + } + + if constexpr (pidMode != kPidPi) { + const float distSq = track.tpcNSigmaPi() * track.tpcNSigmaPi() + track.tofNSigmaPi() * track.tofNSigmaPi(); + if (distSq <= distSpeciesSq) { + return false; + } + } + + if constexpr (pidMode != kPidKa) { + const float distSq = track.tpcNSigmaKa() * track.tpcNSigmaKa() + track.tofNSigmaKa() * track.tofNSigmaKa(); + if (distSq <= distSpeciesSq) { + return false; + } + } + + if constexpr (pidMode != kPidPr) { + const float distSq = track.tpcNSigmaPr() * track.tpcNSigmaPr() + track.tofNSigmaPr() * track.tofNSigmaPr(); + if (distSq <= distSpeciesSq) { + return false; + } + } + // FIX ME! : Relative check will be updated in future, currently implemented for pion, kaon and proton only + // if constexpr (pidMode != kPidEl) { + // const float distSq = track.tpcNSigmaEl() * track.tpcNSigmaEl() + track.tofNSigmaEl() * track.tofNSigmaEl(); + // if (distSq <= distSpeciesSq) + // return false; + // } + + // if constexpr (pidMode != kPidMu) { + // const float distSq = track.tpcNSigmaMu() * track.tpcNSigmaMu() + track.tofNSigmaMu() * track.tofNSigmaMu(); + // if (distSq <= distSpeciesSq) + // return false; + // } + + // if constexpr (pidMode != kPidDe) { + // const float distSq = track.tpcNSigmaDe() * track.tpcNSigmaDe() + track.tofNSigmaDe() * track.tofNSigmaDe(); + // if (distSq <= distSpeciesSq) + // return false; + // } + + return true; +} + +template +bool selIdRectangularCut(const T& track, const float& nSigmaTPC, const float& nSigmaTOF) +{ + if constexpr (pidMode == kPidPi) { + return std::fabs(track.tpcNSigmaPi()) < nSigmaTPC && std::fabs(track.tofNSigmaPi()) < nSigmaTOF; + } else if constexpr (pidMode == kPidKa) { + return std::fabs(track.tpcNSigmaKa()) < nSigmaTPC && std::fabs(track.tofNSigmaKa()) < nSigmaTOF; + } else if constexpr (pidMode == kPidPr) { + return std::fabs(track.tpcNSigmaPr()) < nSigmaTPC && std::fabs(track.tofNSigmaPr()) < nSigmaTOF; + } else if constexpr (pidMode == kPidEl) { + return std::fabs(track.tpcNSigmaEl()) < nSigmaTPC && std::fabs(track.tofNSigmaEl()) < nSigmaTOF; + } else if constexpr (pidMode == kPidMu) { + return std::fabs(track.tpcNSigmaMu()) < nSigmaTPC && std::fabs(track.tofNSigmaMu()) < nSigmaTOF; + } else if constexpr (pidMode == kPidDe) { + return std::fabs(track.tpcNSigmaDe()) < nSigmaTPC && std::fabs(track.tofNSigmaDe()) < nSigmaTOF; + } else { + return false; + } +} + +template +bool selIdEllipsoidalCut(const T& track, const float& nSigmaTPC, const float& nSigmaTOF) +{ + if constexpr (pidMode == kPidPi) { + float tpc = track.tpcNSigmaPi() / nSigmaTPC; + float tof = track.tofNSigmaPi() / nSigmaTOF; + return (tpc * tpc + tof * tof) < 1.0f; + } else if constexpr (pidMode == kPidKa) { + float tpc = track.tpcNSigmaKa() / nSigmaTPC; + float tof = track.tofNSigmaKa() / nSigmaTOF; + return (tpc * tpc + tof * tof) < 1.0f; + } else if constexpr (pidMode == kPidPr) { + float tpc = track.tpcNSigmaPr() / nSigmaTPC; + float tof = track.tofNSigmaPr() / nSigmaTOF; + return (tpc * tpc + tof * tof) < 1.0f; + } else if constexpr (pidMode == kPidEl) { + float tpc = track.tpcNSigmaEl() / nSigmaTPC; + float tof = track.tofNSigmaEl() / nSigmaTOF; + return (tpc * tpc + tof * tof) < 1.0f; + } else if constexpr (pidMode == kPidMu) { + float tpc = track.tpcNSigmaMu() / nSigmaTPC; + float tof = track.tofNSigmaMu() / nSigmaTOF; + return (tpc * tpc + tof * tof) < 1.0f; + } else if constexpr (pidMode == kPidDe) { + float tpc = track.tpcNSigmaDe() / nSigmaTPC; + float tof = track.tofNSigmaDe() / nSigmaTOF; + return (tpc * tpc + tof * tof) < 1.0f; + } else { + return false; + } +} + +// Circular cut : (tpcNSigma^2 + tofNSigma^2) < nSigmaSquaredRad ; for 3 sigma, nSigmaSquaredRad = 9.0 +template +constexpr bool selIdCircularCut(const T& track, const float& nSigmaSquaredRad) +{ + float tpc = 0.f; + float tof = 0.f; + + if constexpr (pidMode == kPidPi) { + tpc = track.tpcNSigmaPi(); + tof = track.tofNSigmaPi(); + } else if constexpr (pidMode == kPidKa) { + tpc = track.tpcNSigmaKa(); + tof = track.tofNSigmaKa(); + } else if constexpr (pidMode == kPidPr) { + tpc = track.tpcNSigmaPr(); + tof = track.tofNSigmaPr(); + } else if constexpr (pidMode == kPidEl) { + tpc = track.tpcNSigmaEl(); + tof = track.tofNSigmaEl(); + } else if constexpr (pidMode == kPidMu) { + tpc = track.tpcNSigmaMu(); + tof = track.tofNSigmaMu(); + } else if constexpr (pidMode == kPidDe) { + tpc = track.tpcNSigmaDe(); + tof = track.tofNSigmaDe(); + } else { + return false; // unknown pidMode + } + return (tpc * tpc + tof * tof) < nSigmaSquaredRad; +} + +template +bool applyVetoOthersTPC(const T& track, const std::array& doVetoTPC, const std::array& vetoTPC) +{ + if constexpr (pidMode != kPidPi) { + if (doVetoTPC[kPidPi] && std::fabs(track.tpcNSigmaPi()) < vetoTPC[kPidPi]) { + return false; + } + } + + if constexpr (pidMode != kPidKa) { + if (doVetoTPC[kPidKa] && std::fabs(track.tpcNSigmaKa()) < vetoTPC[kPidKa]) { + return false; + } + } + + if constexpr (pidMode != kPidPr) { + if (doVetoTPC[kPidPr] && std::fabs(track.tpcNSigmaPr()) < vetoTPC[kPidPr]) { + return false; + } + } + + if constexpr (pidMode != kPidEl) { + if (doVetoTPC[kPidEl] && std::fabs(track.tpcNSigmaEl()) < vetoTPC[kPidEl]) { + return false; + } + } + + if constexpr (pidMode != kPidMu) { + if (doVetoTPC[kPidMu] && std::fabs(track.tpcNSigmaMu()) < vetoTPC[kPidMu]) { + return false; + } + } + + if constexpr (pidMode != kPidDe) { + if (doVetoTPC[kPidDe] && std::fabs(track.tpcNSigmaDe()) < vetoTPC[kPidDe]) { + return false; + } + } + + return true; +} + +template +bool applyVetoOthersTOF(const T& track, const std::array& doVetoTOF, const std::array& vetoTOF) +{ + if constexpr (pidMode != kPidPi) { + if (doVetoTOF[kPidPi] && std::fabs(track.tofNSigmaPi()) < vetoTOF[kPidPi]) { + return false; + } + } + + if constexpr (pidMode != kPidKa) { + if (doVetoTOF[kPidKa] && std::fabs(track.tofNSigmaKa()) < vetoTOF[kPidKa]) { + return false; + } + } + + if constexpr (pidMode != kPidPr) { + if (doVetoTOF[kPidPr] && std::fabs(track.tofNSigmaPr()) < vetoTOF[kPidPr]) { + return false; + } + } + + if constexpr (pidMode != kPidEl) { + if (doVetoTOF[kPidEl] && std::fabs(track.tofNSigmaEl()) < vetoTOF[kPidEl]) { + return false; + } + } + + if constexpr (pidMode != kPidMu) { + if (doVetoTOF[kPidMu] && std::fabs(track.tofNSigmaMu()) < vetoTOF[kPidMu]) { + return false; + } + } + + if constexpr (pidMode != kPidDe) { + if (doVetoTOF[kPidDe] && std::fabs(track.tofNSigmaDe()) < vetoTOF[kPidDe]) { + return false; + } + } + + return true; +} + struct NchCumulantsId { HistogramRegistry hist{"hist", {}, OutputObjHandlingPolicy::AnalysisObject}; @@ -223,64 +653,25 @@ struct NchCumulantsId { Configurable cfgCutPtMaxKa{"cfgCutPtMaxKa", 1.20, "max pT for kaon eff-correction inclusion (purity cutoff)"}; Configurable cfgCutPtMinPr{"cfgCutPtMinPr", 0.40, "min pT for proton eff-correction inclusion"}; Configurable cfgCutPtMaxPr{"cfgCutPtMaxPr", 1.50, "max pT for proton eff-correction inclusion (purity cutoff)"}; + Configurable cfgMaxWeight{"cfgMaxWeight", 20.0f, "max reco weights, use high values to disable the cap"}; + Configurable cfgGateGenByWeight{"cfgGateGenByWeight", true, "count a gen particle iff that is counted by reco looking at pt eta"}; + Configurable cfgRequireSameMcColl{"cfgRequireSameMcColl", true, "ignore reco track that belongs to a diff collision"}; } cfgEventSelection; - // Configurables for particle Identification - Configurable cfgId01CheckVetoCut{"cfgId01CheckVetoCut", true, "cfgId01CheckVetoCut"}; - Configurable cfgId02DoElRejection{"cfgId02DoElRejection", true, "cfgId02DoElRejection"}; - Configurable cfgId03DoDeRejection{"cfgId03DoDeRejection", false, "cfgId03DoDeRejection"}; - Configurable cfgId04DoPdependentId{"cfgId04DoPdependentId", true, "cfgId04DoPdependentId"}; - Configurable cfgId05DoTpcInnerParamId{"cfgId05DoTpcInnerParamId", false, "cfgId05DoTpcInnerParamId"}; - - Configurable cfgIdPi01ThrPforTOF{"cfgIdPi01ThrPforTOF", 0.7, "cfgIdPi01ThrPforTOF"}; - Configurable cfgIdPi02IdCutTypeLowP{"cfgIdPi02IdCutTypeLowP", 1, "cfgIdPi02IdCutTypeLowP"}; - Configurable cfgIdPi03NSigmaTPCLowP{"cfgIdPi03NSigmaTPCLowP", 2.0, "cfgIdPi03NSigmaTPCLowP"}; - Configurable cfgIdPi04NSigmaTOFLowP{"cfgIdPi04NSigmaTOFLowP", 2.0, "cfgIdPi04NSigmaTOFLowP"}; - Configurable cfgIdPi05NSigmaRadLowP{"cfgIdPi05NSigmaRadLowP", 4.0, "cfgIdPi05NSigmaRadLowP"}; - Configurable cfgIdPi06IdCutTypeHighP{"cfgIdPi06IdCutTypeHighP", 1, "cfgIdPi06IdCutTypeHighP"}; - Configurable cfgIdPi07NSigmaTPCHighP{"cfgIdPi07NSigmaTPCHighP", 2.0, "cfgIdPi07NSigmaTPCHighP"}; - Configurable cfgIdPi08NSigmaTOFHighP{"cfgIdPi08NSigmaTOFHighP", 2.0, "cfgIdPi08NSigmaTOFHighP"}; - Configurable cfgIdPi09NSigmaRadHighP{"cfgIdPi09NSigmaRadHighP", 4.0, "cfgIdPi09NSigmaRadHighP"}; - - Configurable cfgIdKa01ThrPforTOF{"cfgIdKa01ThrPforTOF", 0.8, "cfgIdKa01ThrPforTOF"}; - Configurable cfgIdKa02IdCutTypeLowP{"cfgIdKa02IdCutTypeLowP", 1, "cfgIdKa02IdCutTypeLowP"}; - Configurable cfgIdKa03NSigmaTPCLowP{"cfgIdKa03NSigmaTPCLowP", 2.0, "cfgIdKa03NSigmaTPCLowP"}; - Configurable cfgIdKa04NSigmaTOFLowP{"cfgIdKa04NSigmaTOFLowP", 2.0, "cfgIdKa04NSigmaTOFLowP"}; - Configurable cfgIdKa05NSigmaRadLowP{"cfgIdKa05NSigmaRadLowP", 4.0, "cfgIdKa05NSigmaRadLowP"}; - Configurable cfgIdKa06IdCutTypeHighP{"cfgIdKa06IdCutTypeHighP", 1, "cfgIdKa06IdCutTypeHighP"}; - Configurable cfgIdKa07NSigmaTPCHighP{"cfgIdKa07NSigmaTPCHighP", 2.0, "cfgIdKa07NSigmaTPCHighP"}; - Configurable cfgIdKa08NSigmaTOFHighP{"cfgIdKa08NSigmaTOFHighP", 2.0, "cfgIdKa08NSigmaTOFHighP"}; - Configurable cfgIdKa09NSigmaRadHighP{"cfgIdKa09NSigmaRadHighP", 4.0, "cfgIdKa09NSigmaRadHighP"}; - - Configurable cfgIdPr01ThrPforTOF{"cfgIdPr01ThrPforTOF", 0.8, "cfgIdPr01ThrPforTOF"}; - Configurable cfgIdPr02IdCutTypeLowP{"cfgIdPr02IdCutTypeLowP", 1, "cfgIdPr02IdCutTypeLowP"}; - Configurable cfgIdPr03NSigmaTPCLowP{"cfgIdPr03NSigmaTPCLowP", 2.0, "cfgIdPr03NSigmaTPCLowP"}; - Configurable cfgIdPr04NSigmaTOFLowP{"cfgIdPr04NSigmaTOFLowP", 2.0, "cfgIdPr04NSigmaTOFLowP"}; - Configurable cfgIdPr05NSigmaRadLowP{"cfgIdPr05NSigmaRadLowP", 4.0, "cfgIdPr05NSigmaRadLowP"}; - Configurable cfgIdPr06IdCutTypeHighP{"cfgIdPr06IdCutTypeHighP", 1, "cfgIdPr06IdCutTypeHighP"}; - Configurable cfgIdPr07NSigmaTPCHighP{"cfgIdPr07NSigmaTPCHighP", 2.0, "cfgIdPr07NSigmaTPCHighP"}; - Configurable cfgIdPr08NSigmaTOFHighP{"cfgIdPr08NSigmaTOFHighP", 2.0, "cfgIdPr08NSigmaTOFHighP"}; - Configurable cfgIdPr09NSigmaRadHighP{"cfgIdPr09NSigmaRadHighP", 4.0, "cfgIdPr09NSigmaRadHighP"}; + // PID selection: same defaults and configurables as hParticleCorrelation.cxx. + struct : ConfigurableGroup { + Configurable> c1142pidConfigSetting{"c1142pidConfigSetting", {DefaultPIDcheckValues[0].data(), kNPidSelection, kNCutSettings, {"Pi", "Ka", "Pr", "El", "Mu", "De"}, {"ThrPforTOF", "IdCutTypeLowP", "NSigmaTPCLowP", "NSigmaTOFLowP", "NSigmaRadLowP", "IdCutTypeHighP", "NSigmaTPCHighP", "NSigmaTOFHighP", "NSigmaRadHighP", "doVetoOthers", "doRelativeTPCCheck", "doRelativeTOFcheck", "doRelativeTPCTOFcheck"}}, "Cut values for particle identification"}; + Configurable> c1143pidVetoSetting{"c1143pidVetoSetting", {DefaultPidVetoValues[0].data(), kNPidSelection, kNVetoSettings, {"Pi", "Ka", "Pr", "El", "Mu", "De"}, {"doVetoTPC", "doVetoTOF", "vetoTPC", "vetoTOF"}}, "Veto cuts for particle identification"}; + Configurable c1144cfgId07CheckTofBeta{"c1144cfgId07CheckTofBeta", false, "Require beta > 0 for reliable TOF"}; + } cfgIdCut; + // Optional pre-PID electron rejection already present in NchCumulantsId. Configurable cfgIdElRejLowNSigma{"cfgIdElRejLowNSigma", -3.0, "cfgIdElRejLowNSigma"}; Configurable cfgIdElRejHighNSigma{"cfgIdElRejHighNSigma", 5.0, "cfgIdElRejHighNSigma"}; Configurable cfgIdPiRejNSigma{"cfgIdPiRejNSigma", 3.0, "cfgIdPiRejNSigma"}; Configurable cfgIdKaRejNSigma{"cfgIdKaRejNSigma", 3.0, "cfgIdKaRejNSigma"}; Configurable cfgIdPrRejNSigma{"cfgIdPrRejNSigma", 3.0, "cfgIdPrRejNSigma"}; - struct : ConfigurableGroup { - Configurable cfgVetoId01PiTPC{"cfgVetoId01PiTPC", 3.0, "cfgVetoId01PiTPC"}; - Configurable cfgVetoId02PiTOF{"cfgVetoId02PiTOF", 3.0, "cfgVetoId02PiTOF"}; - Configurable cfgVetoId03KaTPC{"cfgVetoId03KaTPC", 3.0, "cfgVetoId03KaTPC"}; - Configurable cfgVetoId04KaTOF{"cfgVetoId04KaTOF", 3.0, "cfgVetoId04KaTOF"}; - Configurable cfgVetoId05PrTPC{"cfgVetoId05PrTPC", 3.0, "cfgVetoId05PrTPC"}; - Configurable cfgVetoId06PrTOF{"cfgVetoId06PrTOF", 3.0, "cfgVetoId06PrTOF"}; - Configurable cfgVetoId07ElTPC{"cfgVetoId07ElTPC", 3.0, "cfgVetoId07ElTPC"}; - Configurable cfgVetoId08ElTOF{"cfgVetoId08ElTOF", 3.0, "cfgVetoId08ElTOF"}; - Configurable cfgVetoId09DeTPC{"cfgVetoId09DeTPC", 3.0, "cfgVetoId09DeTPC"}; - Configurable cfgVetoId10DeTOF{"cfgVetoId10DeTOF", 3.0, "cfgVetoId10DeTOF"}; - } cfgVetoIdCut; - struct : ConfigurableGroup { Configurable cfgCCDB01URL{"cfgCCDB01URL", "http://ccdb-test.cern.ch:8080", "cfgCCDB01URL"}; Configurable cfgCCDB02Path{"cfgCCDB02Path", "Users/p/ppanigra/NchCumulants/correctionWeights/", "cfgCCDB02Path"}; @@ -306,13 +697,27 @@ struct NchCumulantsId { double q2 = 0.; double q3 = 0.; double q4 = 0.; + // bookeeping for counts check not to be used for cumulants calculation + double n = 0.; // track that entered this species (raw, before correction) + double nDrop = 0.; // track that was dropped due to weight > cfgMaxWeight + double sumAll = 0; // sum of all weights (including dropped) }; - static constexpr float KmaxWeight = 10.0f; // secondary safeguard for pr and kaon eff fix + //static constexpr float KmaxWeight = 10.0f; // secondary safeguard for pr and kaon + + inline bool usableWeight(float weight) const + { + return weight > 0.f && weight <= cfgEventSelection.cfgMaxWeight.value; + } inline void fillEffPower(EffPowerSums& p, float weight) { - if (weight <= 0.f || weight > KmaxWeight) { + p.n += 1; + if (weight > 0.f) { + p.sumAll += weight; + } + if (!usableWeight(weight)) { + p.nDrop += 1; return; } @@ -416,6 +821,30 @@ struct NchCumulantsId { } } + template + void addCountCheckProfiles() + { + constexpr std::string_view Dir = CountDir[static_cast(Prefix)]; + registry.add(std::string(Dir) + "nGenAll", "Generated counted Before Weight Gate", kTProfile, {axisCent}); + registry.add(std::string(Dir) + "nGenUsed", "Generated counted After Weight Gate", kTProfile, {axisCent}); + registry.add(std::string(Dir) + "nRecRaw", "Reco tracks counted Before Weight Gate", kTProfile, {axisCent}); + registry.add(std::string(Dir) + "nRecDropped", "Reco tracks dropped after correction(weight <0 or weight > maxWeight)", kTProfile, {axisCent}); + registry.add(std::string(Dir) + "sumWUsable", "reco corrected: sum of usable weights (= Q1)", kTProfile, {axisCent}); + registry.add(std::string(Dir) + "sumWNoCap", "reco corrected: sum of all positive weights, no cap", kTProfile, {axisCent}); + } + + template + void fillCountCheck(float genAll, float genUsed, const EffPowerSums& p, float cent, Registry& hReg) + { + auto base = HIST(CountDir[static_cast(Prefix)]); + hReg.fill(base + HIST("nGenAll"), cent, genAll); + hReg.fill(base + HIST("nGenUsed"), cent, genUsed); + hReg.fill(base + HIST("nRecRaw"), cent, p.n); + hReg.fill(base + HIST("nRecDropped"), cent, p.nDrop); + hReg.fill(base + HIST("sumWUsable"), cent, p.q1); + hReg.fill(base + HIST("sumWNoCap"), cent, p.sumAll); + } + template void fillFCBasis(const EffPowerSums& p, float cent, Registry& hReg) { @@ -588,8 +1017,12 @@ struct NchCumulantsId { if (cfgDebug.printDebugMessages) { LOG(info) << "DEBUG :: Obj Name = " << hPtEtaForEffCorrection[i][j]->GetName() << " :: entries = " << hPtEtaForBinSearch->GetEntries(); } + if (hPtEtaForEffCorrection[i][j]) { + LOG(info) << "MAPCHECK " << name << " integral=" << hPtEtaForEffCorrection[i][j]->Integral(); + } } } + LOG(info) << "MAPCHECK cfgMaxWeight=" << cfgEventSelection.cfgMaxWeight; mgr.setURL("http://alice-ccdb.cern.ch"); // RESET the URL otherwise the other process functions which contains ccdb lookups will fail @@ -830,6 +1263,16 @@ struct NchCumulantsId { addFCGenProfiles(); addFCGenProfiles(); + // profiles for counts + addCountCheckProfiles(); + addCountCheckProfiles(); + addCountCheckProfiles(); + addCountCheckProfiles(); + addCountCheckProfiles(); + addCountCheckProfiles(); + addCountCheckProfiles(); + addCountCheckProfiles(); + } // init ends static constexpr std::array HistRegDire2 = { @@ -867,8 +1310,8 @@ struct NchCumulantsId { }; enum TpcTofCutType { - kRectangularCut = 0, - kCircularCut, + kCircularCut = 0, + kRectangularCut, kEllipsoidalCut }; @@ -908,265 +1351,218 @@ struct NchCumulantsId { doSimProcessing }; - // particle identifications - // tpc Selections - + // Particle identification: ported from hParticleCorrelation.cxx. + // 0 = circular, 1 = rectangular, 2 = ellipsoidal for combined TPC-TOF cuts. template bool vetoIdOthersTPC(const T& track) { - if (pidMode != kPi) { - if (std::fabs(track.tpcNSigmaPi()) < cfgVetoIdCut.cfgVetoId01PiTPC) { - return false; - } - } - if (pidMode != kKa) { - if (std::fabs(track.tpcNSigmaKa()) < cfgVetoIdCut.cfgVetoId03KaTPC) { - return false; - } - } - if (pidMode != kPr) { - if (std::fabs(track.tpcNSigmaPr()) < cfgVetoIdCut.cfgVetoId05PrTPC) { - return false; - } - } - if (cfgId02DoElRejection) { - if (pidMode != kEl) { - if (std::fabs(track.tpcNSigmaEl()) < cfgVetoIdCut.cfgVetoId07ElTPC) { - return false; - } - } - } - if (cfgId03DoDeRejection) { - if (pidMode != kDe) { - if (std::fabs(track.tpcNSigmaDe()) < cfgVetoIdCut.cfgVetoId09DeTPC) { - return false; - } - } - } - return true; + // Static is only run once, ever. + static const std::array doVetoTPC = { + getCfg(cfgIdCut.c1143pidVetoSetting, kPidPi, kDoVetoTPC), + getCfg(cfgIdCut.c1143pidVetoSetting, kPidKa, kDoVetoTPC), + getCfg(cfgIdCut.c1143pidVetoSetting, kPidPr, kDoVetoTPC), + getCfg(cfgIdCut.c1143pidVetoSetting, kPidEl, kDoVetoTPC), + getCfg(cfgIdCut.c1143pidVetoSetting, kPidMu, kDoVetoTPC), + getCfg(cfgIdCut.c1143pidVetoSetting, kPidDe, kDoVetoTPC)}; + + static const std::array vetoTPC = { + getCfg(cfgIdCut.c1143pidVetoSetting, kPidPi, kVetoTPC), + getCfg(cfgIdCut.c1143pidVetoSetting, kPidKa, kVetoTPC), + getCfg(cfgIdCut.c1143pidVetoSetting, kPidPr, kVetoTPC), + getCfg(cfgIdCut.c1143pidVetoSetting, kPidEl, kVetoTPC), + getCfg(cfgIdCut.c1143pidVetoSetting, kPidMu, kVetoTPC), + getCfg(cfgIdCut.c1143pidVetoSetting, kPidDe, kVetoTPC)}; + + return applyVetoOthersTPC(track, doVetoTPC, vetoTPC); } template bool vetoIdOthersTOF(const T& track) { - if (pidMode != kPi) { - if (std::fabs(track.tofNSigmaPi()) < cfgVetoIdCut.cfgVetoId02PiTOF) { - return false; - } - } - if (pidMode != kKa) { - if (std::fabs(track.tofNSigmaKa()) < cfgVetoIdCut.cfgVetoId04KaTOF) { - return false; - } - } - if (pidMode != kPr) { - if (std::fabs(track.tofNSigmaPr()) < cfgVetoIdCut.cfgVetoId06PrTOF) { - return false; - } - } - if (cfgId02DoElRejection) { - if (pidMode != kEl) { - if (std::fabs(track.tofNSigmaEl()) < cfgVetoIdCut.cfgVetoId08ElTOF) { - return false; - } - } - } - if (cfgId03DoDeRejection) { - if (pidMode != kDe) { - if (std::fabs(track.tofNSigmaDe()) < cfgVetoIdCut.cfgVetoId10DeTOF) { - return false; - } - } - } - return true; + // Only computed once + static const std::array doVetoTOF = { + getCfg(cfgIdCut.c1143pidVetoSetting, kPidPi, kDoVetoTOF), + getCfg(cfgIdCut.c1143pidVetoSetting, kPidKa, kDoVetoTOF), + getCfg(cfgIdCut.c1143pidVetoSetting, kPidPr, kDoVetoTOF), + getCfg(cfgIdCut.c1143pidVetoSetting, kPidEl, kDoVetoTOF), + getCfg(cfgIdCut.c1143pidVetoSetting, kPidMu, kDoVetoTOF), + getCfg(cfgIdCut.c1143pidVetoSetting, kPidDe, kDoVetoTOF)}; + + static const std::array vetoTOF = { + getCfg(cfgIdCut.c1143pidVetoSetting, kPidPi, kVetoTOF), + getCfg(cfgIdCut.c1143pidVetoSetting, kPidKa, kVetoTOF), + getCfg(cfgIdCut.c1143pidVetoSetting, kPidPr, kVetoTOF), + getCfg(cfgIdCut.c1143pidVetoSetting, kPidEl, kVetoTOF), + getCfg(cfgIdCut.c1143pidVetoSetting, kPidMu, kVetoTOF), + getCfg(cfgIdCut.c1143pidVetoSetting, kPidDe, kVetoTOF)}; + + return applyVetoOthersTOF(track, doVetoTOF, vetoTOF); } template bool vetoIdOthersTPCTOF(const T& track) { - if (pidMode != kPi) { - if (std::fabs(track.tpcNSigmaPi()) < cfgVetoIdCut.cfgVetoId01PiTPC && std::fabs(track.tofNSigmaPi()) < cfgVetoIdCut.cfgVetoId02PiTOF) { - return false; - } - } - if (pidMode != kKa) { - if (std::fabs(track.tpcNSigmaKa()) < cfgVetoIdCut.cfgVetoId03KaTPC && std::fabs(track.tofNSigmaKa()) < cfgVetoIdCut.cfgVetoId04KaTOF) { - return false; - } - } - if (pidMode != kPr) { - if (std::fabs(track.tpcNSigmaPr()) < cfgVetoIdCut.cfgVetoId05PrTPC && std::fabs(track.tofNSigmaPr()) < cfgVetoIdCut.cfgVetoId06PrTOF) { - return false; - } - } - if (cfgId02DoElRejection) { - if (pidMode != kEl) { - if (std::fabs(track.tpcNSigmaEl()) < cfgVetoIdCut.cfgVetoId07ElTPC && std::fabs(track.tofNSigmaEl()) < cfgVetoIdCut.cfgVetoId08ElTOF) { - return false; - } - } - } - if (cfgId03DoDeRejection) { - if (pidMode != kDe) { - if (std::fabs(track.tpcNSigmaDe()) < cfgVetoIdCut.cfgVetoId09DeTPC && std::fabs(track.tofNSigmaDe()) < cfgVetoIdCut.cfgVetoId10DeTOF) { - return false; - } - } - } - return true; + // If either veto fails, reject + return vetoIdOthersTPC(track) && vetoIdOthersTOF(track); } - template - bool selIdRectangularCut(const T& track, const float& nSigmaTPC, const float& nSigmaTOF) + // Check if TOF is reliable + template + inline bool checkReliableTOF(const T& track) { - switch (pidMode) { - case kPi: - return std::fabs(track.tpcNSigmaPi()) < nSigmaTPC && - std::fabs(track.tofNSigmaPi()) < nSigmaTOF; - - case kKa: - return std::fabs(track.tpcNSigmaKa()) < nSigmaTPC && - std::fabs(track.tofNSigmaKa()) < nSigmaTOF; - - case kPr: - return std::fabs(track.tpcNSigmaPr()) < nSigmaTPC && - std::fabs(track.tofNSigmaPr()) < nSigmaTOF; - - default: - return false; + // which check makes the information of TOF relaiable? should track.beta() be checked? e.g.: + if (cfgIdCut.c1144cfgId07CheckTofBeta) { + return (track.hasTOF() && track.beta() > 0.0f); } + return track.hasTOF(); } template - bool selIdEllipsoidalCut(const T& track, const float& nSigmaTPC, const float& nSigmaTOF) + bool idTPC(const T& track, const float& nSigmaTPC, float& nSigmaIdDistSq) { - switch (pidMode) { - case kPi: - return std::pow(track.tpcNSigmaPi() / nSigmaTPC, 2) + - std::pow(track.tofNSigmaPi() / nSigmaTOF, 2) < - 1.0; - - case kKa: - return std::pow(track.tpcNSigmaKa() / nSigmaTPC, 2) + - std::pow(track.tofNSigmaKa() / nSigmaTOF, 2) < - 1.0; - - case kPr: - return std::pow(track.tpcNSigmaPr() / nSigmaTPC, 2) + - std::pow(track.tofNSigmaPr() / nSigmaTOF, 2) < - 1.0; + if constexpr (pidMode == kPidPi) { + nSigmaIdDistSq = track.tpcNSigmaPi() * track.tpcNSigmaPi(); + } else if constexpr (pidMode == kPidKa) { + nSigmaIdDistSq = track.tpcNSigmaKa() * track.tpcNSigmaKa(); + } else if constexpr (pidMode == kPidPr) { + nSigmaIdDistSq = track.tpcNSigmaPr() * track.tpcNSigmaPr(); + } else if constexpr (pidMode == kPidEl) { + nSigmaIdDistSq = track.tpcNSigmaEl() * track.tpcNSigmaEl(); + } else if constexpr (pidMode == kPidMu) { + nSigmaIdDistSq = track.tpcNSigmaMu() * track.tpcNSigmaMu(); + } else if constexpr (pidMode == kPidDe) { + nSigmaIdDistSq = track.tpcNSigmaDe() * track.tpcNSigmaDe(); + } else { + nSigmaIdDistSq = 1000000; + } - default: + static const bool doVetoOthers = getCfg(cfgIdCut.c1142pidConfigSetting, pidMode, kDoVetoOthers); + if (doVetoOthers) { + if (!vetoIdOthersTPC(track)) { + // If vetoIdOthers = true; it passed all veto checks + // If vetoIdOthers = false; it failed veto check with some other particle return false; + } } - } - template - bool selIdCircularCut(const T& track, const float& nSigmaSquaredRad) - { - switch (pidMode) { - case kPi: - if (std::pow(track.tpcNSigmaPi(), 2) + std::pow(track.tofNSigmaPi(), 2) < nSigmaSquaredRad) { - return true; - } - break; - case kKa: - if (std::pow(track.tpcNSigmaKa(), 2) + std::pow(track.tofNSigmaKa(), 2) < nSigmaSquaredRad) { - return true; - } - break; - case kPr: - if (std::pow(track.tpcNSigmaPr(), 2) + std::pow(track.tofNSigmaPr(), 2) < nSigmaSquaredRad) { - return true; - } - break; - default: + static const bool doRelativeTPCcheck = getCfg(cfgIdCut.c1142pidConfigSetting, pidMode, kDoRelativeTPCcheck); + if (doRelativeTPCcheck) { + if (!relativeIdOthersTPC(track)) { + // If relativeIdOthersTPC = true; particle has stronger nSigma compared to others + // If relativeIdOthersTPC = false; some particle has stronger nSigma compared to it return false; - break; + } } - return false; - } - template - bool checkReliableTOF(const T& track) - { - return track.hasTOF(); // which check makes the information of TOF relaiable? should track.beta() be checked? + if constexpr (pidMode == kPidPi) { + return std::fabs(track.tpcNSigmaPi()) < nSigmaTPC; + } else if constexpr (pidMode == kPidKa) { + return std::fabs(track.tpcNSigmaKa()) < nSigmaTPC; + } else if constexpr (pidMode == kPidPr) { + return std::fabs(track.tpcNSigmaPr()) < nSigmaTPC; + } else if constexpr (pidMode == kPidEl) { + return std::fabs(track.tpcNSigmaEl()) < nSigmaTPC; + } else if constexpr (pidMode == kPidMu) { + return std::fabs(track.tpcNSigmaMu()) < nSigmaTPC; + } else if constexpr (pidMode == kPidDe) { + return std::fabs(track.tpcNSigmaDe()) < nSigmaTPC; + } else { + return false; + } } template - bool idTPC(const T& track, const float& nSigmaTPC) + bool idTPCTOF(const T& track, const int& pidCutType, const float& nSigmaTPC, const float& nSigmaTOF, const float& nSigmaSquaredRad, float& nSigmaIdDistSq) { - if (cfgId01CheckVetoCut && !vetoIdOthersTPC(track)) { - return false; + if constexpr (pidMode == kPidPi) { + nSigmaIdDistSq = track.tpcNSigmaPi() * track.tpcNSigmaPi() + track.tofNSigmaPi() * track.tofNSigmaPi(); + } else if constexpr (pidMode == kPidKa) { + nSigmaIdDistSq = track.tpcNSigmaKa() * track.tpcNSigmaKa() + track.tofNSigmaKa() * track.tofNSigmaKa(); + } else if constexpr (pidMode == kPidPr) { + nSigmaIdDistSq = track.tpcNSigmaPr() * track.tpcNSigmaPr() + track.tofNSigmaPr() * track.tofNSigmaPr(); + } else if constexpr (pidMode == kPidEl) { + nSigmaIdDistSq = track.tpcNSigmaEl() * track.tpcNSigmaEl() + track.tofNSigmaEl() * track.tofNSigmaEl(); + } else if constexpr (pidMode == kPidMu) { + nSigmaIdDistSq = track.tpcNSigmaMu() * track.tpcNSigmaMu() + track.tofNSigmaMu() * track.tofNSigmaMu(); + } else if constexpr (pidMode == kPidDe) { + nSigmaIdDistSq = track.tpcNSigmaDe() * track.tpcNSigmaDe() + track.tofNSigmaDe() * track.tofNSigmaDe(); + } else { + nSigmaIdDistSq = 1000000; } - switch (pidMode) { - case kPi: - if (std::fabs(track.tpcNSigmaPi()) < nSigmaTPC) { - return true; - } - break; - case kKa: - if (std::fabs(track.tpcNSigmaKa()) < nSigmaTPC) { - return true; - } - break; - case kPr: - if (std::fabs(track.tpcNSigmaPr()) < nSigmaTPC) { - return true; - } - break; - default: + static const bool doVetoOthers = getCfg(cfgIdCut.c1142pidConfigSetting, pidMode, kDoVetoOthers); + if (doVetoOthers) { + if (!vetoIdOthersTPCTOF(track)) { + // If vetoIdOthers = true; it passed all veto checks + // If vetoIdOthers = false; it failed veto check with some other particle return false; - break; + } } - return false; - } - template - bool idTPCTOF(const T& track, const int& pidCutType, const float& nSigmaTPC, const float& nSigmaTOF, const float& nSigmaSquaredRad) - { - if (cfgId01CheckVetoCut && !vetoIdOthersTPCTOF(track)) { - return false; + static const bool doRelativeTOFcheck = getCfg(cfgIdCut.c1142pidConfigSetting, pidMode, kDoRelativeTOFcheck); + if (doRelativeTOFcheck) { + if (!relativeIdOthersTOF(track)) { + // If relativeIdOthersTOF = true; particle has stronger nSigma compared to others + // If relativeIdOthersTOF = false; some particle has stronger nSigma compared to it + return false; + } } - if (pidCutType == kRectangularCut) { - return selIdRectangularCut(track, nSigmaTPC, nSigmaTOF); - } else if (pidCutType == kCircularCut) { - return selIdCircularCut(track, nSigmaSquaredRad); - } else if (pidCutType == kEllipsoidalCut) { - return selIdEllipsoidalCut(track, nSigmaTPC, nSigmaTOF); + + static const bool doRelativeTPCTOFcheck = getCfg(cfgIdCut.c1142pidConfigSetting, pidMode, kDoRelativeTPCTOFcheck); + if (doRelativeTPCTOFcheck) { + if (!relativeIdOthersTPCTOF(track)) { + // If relativeIdOthersTPCTOF = true; particle has stronger nSigma compared to others + // If relativeIdOthersTPCTOF = false; some particle has stronger nSigma compared to it + return false; + } + } + + switch (pidCutType) { + case kCircularCut: + return selIdCircularCut(track, nSigmaSquaredRad); + case kRectangularCut: + return selIdRectangularCut(track, nSigmaTPC, nSigmaTOF); + case kEllipsoidalCut: + return selIdEllipsoidalCut(track, nSigmaTPC, nSigmaTOF); + default: + return false; } - return false; } template - bool selPdependent(const T& track, int& IdMethod, const float& cfgIdThrPforTOF, - const int& idCutTypeLowP, const float& nSigmaTPCLowP, const float& nSigmaTOFLowP, const float& nSigmaRadLowP, - const int& idCutTypeHighP, const float& nSigmaTPCHighP, const float& nSigmaTOFHighP, const float& nSigmaRadHighP) + bool selPdependent(const T& track, int& IdMethod, float& nSigmaIdDistSq) { - if (track.p() < cfgIdThrPforTOF) { + // Static cache inside function - initialized once on first call + static const auto thrPforTOF = getCfg(cfgIdCut.c1142pidConfigSetting, pidMode, kThrPforTOF); + static const auto idCutTypeLowP = getCfg(cfgIdCut.c1142pidConfigSetting, pidMode, kIdCutTypeLowP); + static const auto nSigmaTPCLowP = getCfg(cfgIdCut.c1142pidConfigSetting, pidMode, kNSigmaTPCLowP); + static const auto nSigmaTOFLowP = getCfg(cfgIdCut.c1142pidConfigSetting, pidMode, kNSigmaTOFLowP); + static const auto nSigmaRadLowP = getCfg(cfgIdCut.c1142pidConfigSetting, pidMode, kNSigmaRadLowP); + static const auto idCutTypeHighP = getCfg(cfgIdCut.c1142pidConfigSetting, pidMode, kIdCutTypeHighP); + static const auto nSigmaTPCHighP = getCfg(cfgIdCut.c1142pidConfigSetting, pidMode, kNSigmaTPCHighP); + static const auto nSigmaTOFHighP = getCfg(cfgIdCut.c1142pidConfigSetting, pidMode, kNSigmaTOFHighP); + static const auto nSigmaRadHighP = getCfg(cfgIdCut.c1142pidConfigSetting, pidMode, kNSigmaRadHighP); + + if (track.p() < thrPforTOF) { if (checkReliableTOF(track)) { - if (idTPCTOF(track, idCutTypeLowP, nSigmaTPCLowP, nSigmaTOFLowP, nSigmaRadLowP)) { + if (idTPCTOF(track, idCutTypeLowP, nSigmaTPCLowP, nSigmaTOFLowP, nSigmaRadLowP, nSigmaIdDistSq)) { IdMethod = kTPCTOFidentified; return true; } - return false; } else { - if (idTPC(track, nSigmaTPCLowP)) { + if (idTPC(track, nSigmaTPCLowP, nSigmaIdDistSq)) { IdMethod = kTPCidentified; return true; } - return false; } } else { if (checkReliableTOF(track)) { - if (idTPCTOF(track, idCutTypeHighP, nSigmaTPCHighP, nSigmaTOFHighP, nSigmaRadHighP)) { + if (idTPCTOF(track, idCutTypeHighP, nSigmaTPCHighP, nSigmaTOFHighP, nSigmaRadHighP, nSigmaIdDistSq)) { IdMethod = kTPCTOFidentified; return true; } - return false; } - return false; } + return false; } template @@ -1205,40 +1601,26 @@ struct NchCumulantsId { std::fabs(track.tpcNSigmaPr()) > cfgIdPrRejNSigma); } - // Pion + // Existing NchCumulantsId two-argument PID interface (data / MC / purity). template bool selPion(const T& track, int& IdMethod) { - if (cfgId04DoPdependentId) { - return selPdependent(track, IdMethod, - cfgIdPi01ThrPforTOF, cfgIdPi02IdCutTypeLowP, cfgIdPi03NSigmaTPCLowP, cfgIdPi04NSigmaTOFLowP, cfgIdPi05NSigmaRadLowP, - cfgIdPi06IdCutTypeHighP, cfgIdPi07NSigmaTPCHighP, cfgIdPi08NSigmaTOFHighP, cfgIdPi09NSigmaRadHighP); - } - return false; + float nSigmaIdDistSq = -1.0f; + return selPdependent(track, IdMethod, nSigmaIdDistSq); } - // Kaon template bool selKaon(const T& track, int& IdMethod) { - if (cfgId04DoPdependentId) { - return selPdependent(track, IdMethod, - cfgIdKa01ThrPforTOF, cfgIdKa02IdCutTypeLowP, cfgIdKa03NSigmaTPCLowP, cfgIdKa04NSigmaTOFLowP, cfgIdKa05NSigmaRadLowP, - cfgIdKa06IdCutTypeHighP, cfgIdKa07NSigmaTPCHighP, cfgIdKa08NSigmaTOFHighP, cfgIdKa09NSigmaRadHighP); - } - return false; + float nSigmaIdDistSq = -1.0f; + return selPdependent(track, IdMethod, nSigmaIdDistSq); } - // Proton template bool selProton(const T& track, int& IdMethod) { - if (cfgId04DoPdependentId) { - return selPdependent(track, IdMethod, - cfgIdPr01ThrPforTOF, cfgIdPr02IdCutTypeLowP, cfgIdPr03NSigmaTPCLowP, cfgIdPr04NSigmaTOFLowP, cfgIdPr05NSigmaRadLowP, - cfgIdPr06IdCutTypeHighP, cfgIdPr07NSigmaTPCHighP, cfgIdPr08NSigmaTOFHighP, cfgIdPr09NSigmaRadHighP); - } - return false; + float nSigmaIdDistSq = -1.0f; + return selPdependent(track, IdMethod, nSigmaIdDistSq); } template @@ -1478,12 +1860,12 @@ struct NchCumulantsId { using MyAllTracks = soa::Join; + aod::pidTPCFullMu, aod::pidTPCFullDe, aod::pidTOFFullPi, aod::pidTOFFullKa, aod::pidTOFFullPr, aod::pidTOFFullEl, aod::pidTOFFullMu, aod::pidTOFFullDe>; using MyCollisions = soa::Join; using MyTracksWithMclabels = soa::Join; + aod::pidTPCFullPi, aod::pidTPCFullKa, aod::pidTPCFullPr, aod::pidTPCFullEl, aod::pidTPCFullMu, aod::pidTPCFullDe, + aod::pidTOFFullPi, aod::pidTOFFullKa, aod::pidTOFFullPr, aod::pidTOFFullEl, aod::pidTOFFullMu, aod::pidTOFFullDe, aod::McTrackLabels>; using MyCollisionsWithMcLabels = soa::Join; template @@ -1955,6 +2337,16 @@ struct NchCumulantsId { float nKaGen = 0, nAKaGen = 0; float nPiGen = 0, nAPiGen = 0; + std::array genAll{}; // candidates per species (index = FCPrefixEnum 0..7), before the weight gate + auto okCell = [&](FCPrefixEnum sp, int pid, int sgn, double pt, double eta) { + genAll[static_cast(sp)] += 1.f; + if (!cfgEventSelection.cfgGateGenByWeight) { + return true; + } + const int bin = hPtEtaForBinSearch->FindBin(pt, eta); + return usableWeight(hPtEtaForEffCorrection[pid][sgn]->GetBinContent(bin)); + }; + for (const auto& mcTrack : mcTracksTablePerMcColl) { if (!mcTrack.isPhysicalPrimary()) { continue; @@ -1969,58 +2361,71 @@ struct NchCumulantsId { if (pdg == kPiPlus || pdg == kKPlus || pdg == kProton || pdg == kPositron || pdg == kMuonPlus || pdg == kDeuteron) { // fillGenTrackQA(genAnalysis, mcTrack); - nPGen++; - + // nPGen++; genAnalysis.fill(HIST("genAnalysis/Charge/Pos/h12_p"), mcTrack.p()); genAnalysis.fill(HIST("genAnalysis/Charge/Pos/h13_pt"), mcTrack.pt()); genAnalysis.fill(HIST("genAnalysis/Charge/Pos/h14_eta"), mcTrack.eta()); genAnalysis.fill(HIST("genAnalysis/Charge/Pos/h15_phi"), mcTrack.phi()); genAnalysis.fill(HIST("genAnalysis/Charge/Pos/h16_rapidity"), mcTrack.y()); genAnalysis.fill(HIST("genAnalysis/Charge/Pos/h20_pt_eta"), mcTrack.pt(), mcTrack.eta()); + + if (okCell(FCPrefixEnum::Pos, kCh, kPos, mcTrack.pt(), mcTrack.eta())) { + nPGen++; + } } else if (pdg == kPiMinus || pdg == kKMinus || pdg == kProtonBar || pdg == kElectron || pdg == kMuonMinus || pdg == -kDeuteron) { // fillGenTrackQA(genAnalysis, mcTrack); - nMGen++; - + // nMGen++; genAnalysis.fill(HIST("genAnalysis/Charge/Neg/h12_p"), mcTrack.p()); genAnalysis.fill(HIST("genAnalysis/Charge/Neg/h13_pt"), mcTrack.pt()); genAnalysis.fill(HIST("genAnalysis/Charge/Neg/h14_eta"), mcTrack.eta()); genAnalysis.fill(HIST("genAnalysis/Charge/Neg/h15_phi"), mcTrack.phi()); genAnalysis.fill(HIST("genAnalysis/Charge/Neg/h16_rapidity"), mcTrack.y()); genAnalysis.fill(HIST("genAnalysis/Charge/Neg/h20_pt_eta"), mcTrack.pt(), mcTrack.eta()); + + if (okCell(FCPrefixEnum::Neg, kCh, kNeg, mcTrack.pt(), mcTrack.eta())) { + nMGen++; + } } // ----- Pions ----- if (pdg == kPiPlus) { // fillGenTrackQA(genAnalysis, mcTrack); - nPiGen++; + // nPiGen++; fillGenTrackQA(genAnalysis, mcTrack); + if (okCell(FCPrefixEnum::PiPos, kPi, kPos, mcTrack.pt(), mcTrack.eta())) { + nPiGen++; + } } else if (pdg == kPiMinus) { // fillGenTrackQA(genAnalysis, mcTrack); - nAPiGen++; + // nAPiGen++; fillGenTrackQA(genAnalysis, mcTrack); + if (okCell(FCPrefixEnum::PiNeg, kPi, kNeg, mcTrack.pt(), mcTrack.eta())) { + nAPiGen++; + } } else if (pdg == kKPlus) { // fillGenTrackQA(genAnalysis, mcTrack); - if (mcTrack.pt() >= cfgEventSelection.cfgCutPtMinKa && mcTrack.pt() <= cfgEventSelection.cfgCutPtMaxKa) { + fillGenTrackQA(genAnalysis, mcTrack); + if (mcTrack.pt() >= cfgEventSelection.cfgCutPtMinKa && mcTrack.pt() <= cfgEventSelection.cfgCutPtMaxKa && okCell(FCPrefixEnum::KaPos, kKa, kPos, mcTrack.pt(), mcTrack.eta())) { nKaGen++; } - fillGenTrackQA(genAnalysis, mcTrack); + } else if (pdg == kKMinus) { // fillGenTrackQA(genAnalysis, mcTrack); - if (mcTrack.pt() >= cfgEventSelection.cfgCutPtMinKa && mcTrack.pt() <= cfgEventSelection.cfgCutPtMaxKa) { + fillGenTrackQA(genAnalysis, mcTrack); + if (mcTrack.pt() >= cfgEventSelection.cfgCutPtMinKa && mcTrack.pt() <= cfgEventSelection.cfgCutPtMaxKa && okCell(FCPrefixEnum::KaNeg, kKa, kNeg, mcTrack.pt(), mcTrack.eta())) { nAKaGen++; } - fillGenTrackQA(genAnalysis, mcTrack); } else if (pdg == kProton) { // fillGenTrackQA(genAnalysis, mcTrack); - if (mcTrack.pt() >= cfgEventSelection.cfgCutPtMinPr && mcTrack.pt() <= cfgEventSelection.cfgCutPtMaxPr) { + fillGenTrackQA(genAnalysis, mcTrack); + if (mcTrack.pt() >= cfgEventSelection.cfgCutPtMinPr && mcTrack.pt() <= cfgEventSelection.cfgCutPtMaxPr && okCell(FCPrefixEnum::Pr, kPr, kPos, mcTrack.pt(), mcTrack.eta())) { nPrGen++; } - fillGenTrackQA(genAnalysis, mcTrack); } else if (pdg == kProtonBar) { // fillGenTrackQA(genAnalysis, mcTrack); - if (mcTrack.pt() >= cfgEventSelection.cfgCutPtMinPr && mcTrack.pt() <= cfgEventSelection.cfgCutPtMaxPr) { + fillGenTrackQA(genAnalysis, mcTrack); + if (mcTrack.pt() >= cfgEventSelection.cfgCutPtMinPr && mcTrack.pt() <= cfgEventSelection.cfgCutPtMaxPr && okCell(FCPrefixEnum::APr, kPr, kNeg, mcTrack.pt(), mcTrack.eta())) { nAPrGen++; } - fillGenTrackQA(genAnalysis, mcTrack); } } // particle/track loop for gen ends @@ -2076,6 +2481,11 @@ struct NchCumulantsId { if (!mcPart.isPhysicalPrimary()) { continue; } + + if (cfgEventSelection.cfgRequireSameMcColl && mcPart.mcCollisionId() != mcCollision.globalIndex()) { + continue; + } + int pdg = mcPart.pdgCode(); int mcTag = getMCTag(track); @@ -2330,6 +2740,15 @@ struct NchCumulantsId { } fillCollQA(col, nChRec, nTRec); + fillCountCheck(genAll[0], nPrGen, prPow, cent, registry); + fillCountCheck(genAll[1], nAPrGen, aprPow, cent, registry); + fillCountCheck(genAll[2], nPiGen, pipPow, cent, registry); + fillCountCheck(genAll[3], nAPiGen, pimPow, cent, registry); + fillCountCheck(genAll[4], nKaGen, kapPow, cent, registry); + fillCountCheck(genAll[5], nAKaGen, kamPow, cent, registry); + fillCountCheck(genAll[6], nPGen, posPow, cent, registry); + fillCountCheck(genAll[7], nMGen, negPow, cent, registry); + fillFCBasis(prPow, cent, registry); fillFCBasis(aprPow, cent, registry); fillFCBasis(pipPow, cent, registry); From b440ab7b575acb4ceeac0f3aebd8d1252f0d6805 Mon Sep 17 00:00:00 2001 From: Pravata Panigrahi <158029188+pravatp98@users.noreply.github.com> Date: Tue, 10 Nov 2026 11:17:47 +0530 Subject: [PATCH 10/13] Updated for Particle Indentification and correction --- PWGCF/EbyEFluctuations/Tasks/nchCumulantsId.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/PWGCF/EbyEFluctuations/Tasks/nchCumulantsId.cxx b/PWGCF/EbyEFluctuations/Tasks/nchCumulantsId.cxx index 3ea1b590327..a0879e639b6 100644 --- a/PWGCF/EbyEFluctuations/Tasks/nchCumulantsId.cxx +++ b/PWGCF/EbyEFluctuations/Tasks/nchCumulantsId.cxx @@ -703,7 +703,7 @@ struct NchCumulantsId { double sumAll = 0; // sum of all weights (including dropped) }; - //static constexpr float KmaxWeight = 10.0f; // secondary safeguard for pr and kaon + // static constexpr float KmaxWeight = 10.0f; // secondary safeguard for pr and kaon inline bool usableWeight(float weight) const { From 32e84b7650510f1301076d87aecb56c9a8c8e58d Mon Sep 17 00:00:00 2001 From: Pravata Panigrahi <158029188+pravatp98@users.noreply.github.com> Date: Tue, 10 Nov 2026 20:08:10 +0530 Subject: [PATCH 11/13] Add LabeledArray include to nchCumulantsId.cxx --- PWGCF/EbyEFluctuations/Tasks/nchCumulantsId.cxx | 1 + 1 file changed, 1 insertion(+) diff --git a/PWGCF/EbyEFluctuations/Tasks/nchCumulantsId.cxx b/PWGCF/EbyEFluctuations/Tasks/nchCumulantsId.cxx index a0879e639b6..5181492c937 100644 --- a/PWGCF/EbyEFluctuations/Tasks/nchCumulantsId.cxx +++ b/PWGCF/EbyEFluctuations/Tasks/nchCumulantsId.cxx @@ -31,6 +31,7 @@ #include #include #include +#include // #include #include #include From bbb98e68b6ace1f39b005b6cdbfb15b9eb86eb28 Mon Sep 17 00:00:00 2001 From: Pravata Panigrahi <158029188+pravatp98@users.noreply.github.com> Date: Wed, 11 Nov 2026 08:05:04 +0530 Subject: [PATCH 12/13] Add Array2D include and remove LabeledArray include --- PWGCF/EbyEFluctuations/Tasks/nchCumulantsId.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/PWGCF/EbyEFluctuations/Tasks/nchCumulantsId.cxx b/PWGCF/EbyEFluctuations/Tasks/nchCumulantsId.cxx index 5181492c937..37d4ad9f337 100644 --- a/PWGCF/EbyEFluctuations/Tasks/nchCumulantsId.cxx +++ b/PWGCF/EbyEFluctuations/Tasks/nchCumulantsId.cxx @@ -23,6 +23,7 @@ #include #include +#include // #include #include // #include @@ -31,7 +32,6 @@ #include #include #include -#include // #include #include #include From 5608a42acc1332be988d298a500bd032311356fe Mon Sep 17 00:00:00 2001 From: Pravata Panigrahi <158029188+pravatp98@users.noreply.github.com> Date: Wed, 11 Nov 2026 08:29:51 +0530 Subject: [PATCH 13/13] Updated the Code for Particle Identification and fixed for the missing header files