diff --git a/PWGCF/EbyEFluctuations/Tasks/nchCumulantsId.cxx b/PWGCF/EbyEFluctuations/Tasks/nchCumulantsId.cxx index 1fd39ebf2fc..db822fa9d08 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& 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 + 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) { + 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& hReg) + { + // 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) + 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 ─── + hReg.fill(base + HIST("JointF11"), cent, jointF11); + hReg.fill(base + HIST("JointF12"), cent, jointF12); + hReg.fill(base + HIST("JointF21"), cent, jointF21); + + if (cfgEventSelection.cfgDoSubsampling) { + 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& hReg) + { + 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.); + + 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) { + 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& hReg) + { + 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); + + 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) { + 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); + } + } + 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[] = { @@ -916,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); } } @@ -1026,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 (flagUnusedVariableError) + 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) { @@ -1061,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) { @@ -1078,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) { @@ -1095,16 +1436,17 @@ 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 using MyAllTracks = soa::Join; template - bool isEventSelected(const CollisionType& col) + bool isEventSelected(const CollisionType& coll) { - if (cfgEvSel01doNoSameBunchPileup && - !col.selection_bit(aod::evsel::kNoSameBunchPileup)) + if (cfgEventSelection.cfgEvSel01doNoSameBunchPileup && + !coll.selection_bit(aod::evsel::kNoSameBunchPileup)) return false; - if (cfgEvSel02doIsGoodZvtxFT0vsPV && - !col.selection_bit(aod::evsel::kIsGoodZvtxFT0vsPV)) + if (cfgEventSelection.cfgEvSel02doIsGoodZvtxFT0vsPV && + !coll.selection_bit(aod::evsel::kIsGoodZvtxFT0vsPV)) return false; - if (cfgEvSel03doIsGoodITSLayersAll && - !col.selection_bit(aod::evsel::kIsGoodITSLayersAll)) + if (cfgEventSelection.cfgEvSel03doIsGoodITSLayersAll && + !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); @@ -1211,7 +1553,7 @@ struct NchCumulantsId { } // Reject electrons first - if (cfgDoRejectionForId && !selTrackForId(track)) { + if (cfgEventSelection.cfgDoRejectionForId && !selTrackForId(track)) { continue; } @@ -1278,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); @@ -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);