From e88262234ce7f7000e982c47f42c4005515fc216 Mon Sep 17 00:00:00 2001 From: Omar Vazquez Date: Fri, 25 Jul 2025 22:46:36 -0600 Subject: [PATCH] Added more histos for MC closure --- PWGLF/Tasks/GlobalEventProperties/uccZdc.cxx | 916 ++++++++++--------- 1 file changed, 501 insertions(+), 415 deletions(-) diff --git a/PWGLF/Tasks/GlobalEventProperties/uccZdc.cxx b/PWGLF/Tasks/GlobalEventProperties/uccZdc.cxx index 584c3ecfc1b..ea894e460c3 100644 --- a/PWGLF/Tasks/GlobalEventProperties/uccZdc.cxx +++ b/PWGLF/Tasks/GlobalEventProperties/uccZdc.cxx @@ -69,8 +69,8 @@ using SimTracks = soa::Join, kSizeBootStrapEnsemble> hPoisson{}; std::array, kSizeBootStrapEnsemble> hNch{}; +std::array, kSizeBootStrapEnsemble> hPoisson{}; std::array, kSizeBootStrapEnsemble> pNchVsOneParCorrVsZN{}; std::array, kSizeBootStrapEnsemble> pNchVsTwoParCorrVsZN{}; std::array, kSizeBootStrapEnsemble> pNchVsThreeParCorrVsZN{}; @@ -83,18 +83,24 @@ std::array, kSizeBootStrapEnsemble> pOneParCorrVsV0A{} std::array, kSizeBootStrapEnsemble> pTwoParCorrVsV0A{}; std::array, kSizeBootStrapEnsemble> pThreeParCorrVsV0A{}; +std::array, kSizeBootStrapEnsemble> pOneParCorrVsNch{}; +std::array, kSizeBootStrapEnsemble> pTwoParCorrVsNch{}; +std::array, kSizeBootStrapEnsemble> pThreeParCorrVsNch{}; + std::array, kSizeBootStrapEnsemble> hPoissonMC{}; std::array, kSizeBootStrapEnsemble> hNchGen{}; -std::array, kSizeBootStrapEnsemble> pNchvsOneParCorrGen{}; -std::array, kSizeBootStrapEnsemble> pNchvsTwoParCorrGen{}; -std::array, kSizeBootStrapEnsemble> pNchvsThreeParCorrGen{}; -std::array, kSizeBootStrapEnsemble> pNchvsFourParCorrGen{}; -std::array, kSizeBootStrapEnsemble> hNchRec{}; -std::array, kSizeBootStrapEnsemble> pNchvsOneParCorrRec{}; -std::array, kSizeBootStrapEnsemble> pNchvsTwoParCorrRec{}; -std::array, kSizeBootStrapEnsemble> pNchvsThreeParCorrRec{}; -std::array, kSizeBootStrapEnsemble> pNchvsFourParCorrRec{}; +// std::array, kSizeBootStrapEnsemble> pOneParCorrVsT0MGen{}; +// std::array, kSizeBootStrapEnsemble> pTwoParCorrVsT0MGen{}; +// std::array, kSizeBootStrapEnsemble> pThreeParCorrVsT0MGen{}; +// +// std::array, kSizeBootStrapEnsemble> pOneParCorrVsV0AGen{}; +// std::array, kSizeBootStrapEnsemble> pTwoParCorrVsV0AGen{}; +// std::array, kSizeBootStrapEnsemble> pThreeParCorrVsV0AGen{}; + +std::array, kSizeBootStrapEnsemble> pOneParCorrVsNchGen{}; +std::array, kSizeBootStrapEnsemble> pTwoParCorrVsNchGen{}; +std::array, kSizeBootStrapEnsemble> pThreeParCorrVsNchGen{}; struct UccZdc { @@ -317,34 +323,32 @@ struct UccZdc { xEvtsDiv->SetBinLabel(1, "MC closure"); xEvtsDiv->SetBinLabel(2, "Corrections"); // MC closure - registry.add("NchGen", "MC closure;#it{N}_{ch} (|#eta| < 0.8);Entries;", kTH1F, {{nBinsNch, minNch, maxNch}}); - registry.add("NchvsOneParCorrGen", "MC closure;#it{N}_{ch} (|#eta| < 0.8);#LT[#it{p}_{T}^{(1)}]#GT (GeV/#it{c})", kTProfile, {{nBinsNch, minNch, maxNch}}); - registry.add("NchvsTwoParCorrGen", "MC closure;#it{N}_{ch} (|#eta| < 0.8);#LT[#it{p}_{T}^{(2)}]#GT", kTProfile, {{nBinsNch, minNch, maxNch}}); - registry.add("NchvsThreeParCorrGen", "MC closure;#it{N}_{ch} (|#eta| < 0.8);#LT[#it{p}_{T}^{(3)}]#GT", kTProfile, {{nBinsNch, minNch, maxNch}}); - registry.add("NchvsFourParCorrGen", "MC closure;#it{N}_{ch} (|#eta| < 0.8);#LT[#it{p}_{T}^{(4)}]#GT", kTProfile, {{nBinsNch, minNch, maxNch}}); - registry.add("NchVsTwoParCorr", "MC closure;#it{N}_{ch} (|#eta| < 0.8, Corrected);#LT[#it{p}_{T}^{(2)}]#GT", kTProfile, {{nBinsNch, minNch, maxNch}}); - registry.add("NchVsThreeParCorr", "MC closure;#it{N}_{ch} (|#eta| < 0.8, Corrected);#LT[#it{p}_{T}^{(3)}]#GT", kTProfile, {{nBinsNch, minNch, maxNch}}); - registry.add("NchVsFourParCorr", "MC closure;#it{N}_{ch} (|#eta| < 0.8, Corrected);#LT[#it{p}_{T}^{(4)}]#GT", kTProfile, {{nBinsNch, minNch, maxNch}}); + registry.add("NchGen", Form("MC Closure;%s;Entries", tiNch), kTH1F, {{nBinsNch, minNch, maxNch}}); + registry.add("NchvsOneParCorrGen", Form("MC Closure;%s;%s", tiNch, tiOneParCorr), kTProfile, {{nBinsNch, minNch, maxNch}}); + registry.add("NchvsTwoParCorrGen", Form("MC Closure;%s;%s", tiNch, tiTwoParCorr), kTProfile, {{nBinsNch, minNch, maxNch}}); + registry.add("NchvsThreeParCorrGen", Form("MC Closure;%s;%s", tiNch, tiThreeParCorr), kTProfile, {{nBinsNch, minNch, maxNch}}); + registry.add("NchVsTwoParCorr", Form("MC Closure;%s;%s", tiNch, tiTwoParCorr), kTProfile, {{nBinsNch, minNch, maxNch}}); + registry.add("NchVsThreeParCorr", Form("MC Closure;%s;%s", tiNch, tiThreeParCorr), kTProfile, {{nBinsNch, minNch, maxNch}}); // Corrections registry.add("zPosMC", "Filled at MC closure + Corrections;;Entries;", kTH1F, {axisZpos}); registry.add("hEventCounterMC", "Event counter", kTH1F, {axisEvent}); registry.add("nRecColvsCent", "", kTH2F, {{6, -0.5, 5.5}, {{axisCent}}}); - registry.add("Pt_all_ch", "Corrections;#it{N}_{ch} (|#eta|<0.8);;", kTH2F, {{nBinsNch, minNch, maxNch}, {axisPt}}); - registry.add("Pt_ch", "Corrections;#it{N}_{ch} (|#eta|<0.8);;", kTH2F, {{nBinsNch, minNch, maxNch}, {axisPt}}); - registry.add("Pt_pi", "Corrections;#it{N}_{ch} (|#eta|<0.8);;", kTH2F, {{nBinsNch, minNch, maxNch}, {axisPt}}); - registry.add("Pt_ka", "Corrections;#it{N}_{ch} (|#eta|<0.8);;", kTH2F, {{nBinsNch, minNch, maxNch}, {axisPt}}); - registry.add("Pt_pr", "Corrections;#it{N}_{ch} (|#eta|<0.8);;", kTH2F, {{nBinsNch, minNch, maxNch}, {axisPt}}); - registry.add("Pt_sigpos", "Corrections;;;", kTH2F, {{nBinsNch, minNch, maxNch}, {axisPt}}); - registry.add("Pt_signeg", "Corrections;;;", kTH2F, {{nBinsNch, minNch, maxNch}, {axisPt}}); - registry.add("Pt_re", "Corrections;;;", kTH2F, {{nBinsNch, minNch, maxNch}, {axisPt}}); - registry.add("PtMC_ch", "Corrections;;;", kTH2F, {{nBinsNch, minNch, maxNch}, {axisPt}}); - registry.add("PtMC_pi", "Corrections;;;", kTH2F, {{nBinsNch, minNch, maxNch}, {axisPt}}); - registry.add("PtMC_ka", "Corrections;;;", kTH2F, {{nBinsNch, minNch, maxNch}, {axisPt}}); - registry.add("PtMC_pr", "Corrections;;;", kTH2F, {{nBinsNch, minNch, maxNch}, {axisPt}}); - registry.add("PtMC_sigpos", "Corrections;;;", kTH2F, {{nBinsNch, minNch, maxNch}, {axisPt}}); - registry.add("PtMC_signeg", "Corrections;;;", kTH2F, {{nBinsNch, minNch, maxNch}, {axisPt}}); - registry.add("PtMC_re", "Corrections;;;", kTH2F, {{nBinsNch, minNch, maxNch}, {axisPt}}); - registry.add("McNchVsFT0M", ";T0A+T0C (#times 1/100, -3.3 < #eta < -2.1 and 3.5 < #eta < 4.9);#it{N}_{ch} (|#eta|<0.8);", kTH2F, {{{nBinsAmpFT0, 0., maxAmpFT0}, {nBinsNch, minNch, maxNch}}}); + registry.add("Pt_all_ch", Form("Corrections;%s;%s", tiNch, tiPt), kTH2F, {{nBinsNch, minNch, maxNch}, {axisPt}}); + registry.add("Pt_ch", Form("Corrections;%s;%s", tiNch, tiPt), kTH2F, {{nBinsNch, minNch, maxNch}, {axisPt}}); + registry.add("Pt_pi", Form("Corrections;%s;%s", tiNch, tiPt), kTH2F, {{nBinsNch, minNch, maxNch}, {axisPt}}); + registry.add("Pt_ka", Form("Corrections;%s;%s", tiNch, tiPt), kTH2F, {{nBinsNch, minNch, maxNch}, {axisPt}}); + registry.add("Pt_pr", Form("Corrections;%s;%s", tiNch, tiPt), kTH2F, {{nBinsNch, minNch, maxNch}, {axisPt}}); + registry.add("Pt_sigpos", Form("Corrections;%s;%s", tiNch, tiPt), kTH2F, {{nBinsNch, minNch, maxNch}, {axisPt}}); + registry.add("Pt_signeg", Form("Corrections;%s;%s", tiNch, tiPt), kTH2F, {{nBinsNch, minNch, maxNch}, {axisPt}}); + registry.add("Pt_re", Form("Corrections;%s;%s", tiNch, tiPt), kTH2F, {{nBinsNch, minNch, maxNch}, {axisPt}}); + registry.add("PtMC_ch", Form("Corrections;%s;%s", tiNch, tiPt), kTH2F, {{nBinsNch, minNch, maxNch}, {axisPt}}); + registry.add("PtMC_pi", Form("Corrections;%s;%s", tiNch, tiPt), kTH2F, {{nBinsNch, minNch, maxNch}, {axisPt}}); + registry.add("PtMC_ka", Form("Corrections;%s;%s", tiNch, tiPt), kTH2F, {{nBinsNch, minNch, maxNch}, {axisPt}}); + registry.add("PtMC_pr", Form("Corrections;%s;%s", tiNch, tiPt), kTH2F, {{nBinsNch, minNch, maxNch}, {axisPt}}); + registry.add("PtMC_sigpos", Form("Corrections;%s;%s", tiNch, tiPt), kTH2F, {{nBinsNch, minNch, maxNch}, {axisPt}}); + registry.add("PtMC_signeg", Form("Corrections;%s;%s", tiNch, tiPt), kTH2F, {{nBinsNch, minNch, maxNch}, {axisPt}}); + registry.add("PtMC_re", Form("Corrections;%s;%s", tiNch, tiPt), kTH2F, {{nBinsNch, minNch, maxNch}, {axisPt}}); + registry.add("McNchVsFT0M", Form("Corrections;%s;%s", tiT0M, tiNch), kTH2F, {{{nBinsAmpFT0, 0., maxAmpFT0}, {nBinsNch, minNch, maxNch}}}); auto hECMC = registry.get(HIST("hEventCounterMC")); auto* x = hECMC->GetXaxis(); @@ -352,18 +356,36 @@ struct UccZdc { x->SetBinLabel(13, "VtxZ cut"); for (int i = 0; i < kSizeBootStrapEnsemble; i++) { - hPoissonMC[i] = registry.add(Form("PoissonMC_Replica%d", i), ";#it{k};Entries", kTH1F, {{21, -0.5, 20.5}}); - hNchGen[i] = registry.add(Form("NchGen_Replica%d", i), ";#it{N}_{ch} (|#eta| < 0.8);Entries", kTH1F, {{nBinsNch, minNch, maxNch}}); - pNchvsOneParCorrGen[i] = registry.add(Form("NchvsOneParCorrGen_Replica%d", i), ";#it{N}_{ch}, |#eta| < 0.8; One-particle #it{p}_{T} correlator", kTProfile, {{nBinsNch, minNch, maxNch}}); - pNchvsTwoParCorrGen[i] = registry.add(Form("NchvsTwoParCorrGen_Replica%d", i), ";#it{N}_{ch}, |#eta| < 0.8; Two-particle #it{p}_{T} correlator", kTProfile, {{nBinsNch, minNch, maxNch}}); - pNchvsThreeParCorrGen[i] = registry.add(Form("NchvsThreeParCorrGen_Replica%d", i), ";#it{N}_{ch}, |#eta| < 0.8; Three-particle #it{p}_{T} correlator", kTProfile, {{nBinsNch, minNch, maxNch}}); - pNchvsFourParCorrGen[i] = registry.add(Form("NchvsFourParCorrGen_Replica%d", i), ";#it{N}_{ch}, |#eta| < 0.8; Four-particle #it{p}_{T} correlator", kTProfile, {{nBinsNch, minNch, maxNch}}); - - hNchRec[i] = registry.add(Form("NchRec_Replica%d", i), ";#it{N}_{ch} (|#eta| < 0.8);Entries", kTH1F, {{nBinsNch, minNch, maxNch}}); - pNchvsOneParCorrRec[i] = registry.add(Form("NchvsOneParCorrRec_Replica%d", i), ";#it{N}_{ch}, |#eta| < 0.8; One-particle #it{p}_{T} correlator", kTProfile, {{nBinsNch, minNch, maxNch}}); - pNchvsTwoParCorrRec[i] = registry.add(Form("NchvsTwoParCorrRec_Replica%d", i), ";#it{N}_{ch}, |#eta| < 0.8; Two-particle #it{p}_{T} correlator", kTProfile, {{nBinsNch, minNch, maxNch}}); - pNchvsThreeParCorrRec[i] = registry.add(Form("NchvsThreeParCorrRec_Replica%d", i), ";#it{N}_{ch}, |#eta| < 0.8; Three-particle #it{p}_{T} correlator", kTProfile, {{nBinsNch, minNch, maxNch}}); - pNchvsFourParCorrRec[i] = registry.add(Form("NchvsFourParCorrRec_Replica%d", i), ";#it{N}_{ch}, |#eta| < 0.8; Four-particle #it{p}_{T} correlator", kTProfile, {{nBinsNch, minNch, maxNch}}); + + hPoissonMC[i] = registry.add(Form("PoissonMC_Replica%d", i), ";#it{k};Entries", kTH1F, {{11, -0.5, 10.5}}); + hNchGen[i] = registry.add(Form("NchGen_Replica%d", i), Form(";%s;Entries", tiNch), kTH1F, {{nBinsNch, minNch, maxNch}}); + + pOneParCorrVsNchGen[i] = registry.add(Form("OneParCorrVsNchGen_Replica%d", i), Form(";%s;%s;", tiNch, tiOneParCorr), kTProfile, {{nBinsNch, minNch, maxNch}}); + pTwoParCorrVsNchGen[i] = registry.add(Form("TwoParCorrVsNchGen_Replica%d", i), Form(";%s;%s;", tiNch, tiTwoParCorr), kTProfile, {{nBinsNch, minNch, maxNch}}); + pThreeParCorrVsNchGen[i] = registry.add(Form("ThreeParCorrVsNchGen_Replica%d", i), Form(";%s;%s;", tiNch, tiThreeParCorr), kTProfile, {{nBinsNch, minNch, maxNch}}); + + // pOneParCorrVsT0MGen[i] = registry.add(Form("OneParCorrVsT0MGen_Replica%d",i),Form(";%s;%s;",tiT0M,tiOneParCorr), kTProfile,{{nBinsAmpFT0,0.,maxAmpFT0}}); + // pTwoParCorrVsT0MGen[i] = registry.add(Form("TwoParCorrVsT0MGen_Replica%d",i),Form(";%s;%s;",tiT0M,tiTwoParCorr), kTProfile,{{nBinsAmpFT0,0.,maxAmpFT0}}); + // pThreeParCorrVsT0MGen[i] = registry.add(Form("ThreeParCorrVsT0MGen_Replica%d",i),Form(";%s;%s;",tiT0M,tiThreeParCorr), kTProfile,{{nBinsAmpFT0,0.,maxAmpFT0}}); + + // pOneParCorrVsV0AGen[i] = registry.add(Form("OneParCorrVsV0AGen_Replica%d",i),Form(";%s;%s;",tiT0M,tiOneParCorr), kTProfile,{{nBinsAmpFT0,0.,maxAmpFT0}}); + // pTwoParCorrVsV0AGen[i] = registry.add(Form("TwoParCorrVsV0AGen_Replica%d",i),Form(";%s;%s;",tiT0M,tiTwoParCorr), kTProfile,{{nBinsAmpFT0,0.,maxAmpFT0}}); + // pThreeParCorrVsV0AGen[i] = registry.add(Form("ThreeParCorrVsV0AGen_Replica%d",i),Form(";%s;%s;",tiT0M,tiThreeParCorr), kTProfile,{{nBinsAmpFT0,0.,maxAmpFT0}}); + + hNch[i] = registry.add(Form("Nch_Replica%d", i), Form(";%s;Entries", tiNch), kTH1F, {{nBinsNch, minNch, maxNch}}); + hPoisson[i] = registry.add(Form("Poisson_Replica%d", i), ";#it{k};Entries", kTH1F, {{11, -0.5, 10.5}}); + + pOneParCorrVsNch[i] = registry.add(Form("OneParCorrVsNch_Replica%d", i), Form(";%s;%s;", tiNch, tiOneParCorr), kTProfile, {{nBinsNch, minNch, maxNch}}); + pTwoParCorrVsNch[i] = registry.add(Form("TwoParCorrVsNch_Replica%d", i), Form(";%s;%s;", tiNch, tiTwoParCorr), kTProfile, {{nBinsNch, minNch, maxNch}}); + pThreeParCorrVsNch[i] = registry.add(Form("ThreeParCorrVsNch_Replica%d", i), Form(";%s;%s;", tiNch, tiTwoParCorr), kTProfile, {{nBinsNch, minNch, maxNch}}); + + // pOneParCorrVsT0M[i] = registry.add(Form("OneParCorrVsT0M_Replica%d",i),Form(";%s;%s;",tiT0M,tiOneParCorr),kTProfile,{{nBinsAmpFT0,0.,maxAmpFT0}}); + // pTwoParCorrVsT0M[i] = registry.add(Form("TwoParCorrVsT0M_Replica%d",i),Form(";%s;%s;",tiT0M,tiTwoParCorr),kTProfile,{{nBinsAmpFT0,0.,maxAmpFT0}}); + // pThreeParCorrVsT0M[i] = registry.add(Form("ThreeParCorrVsT0M_Replica%d",i),Form(";%s;%s;",tiT0M,tiThreeParCorr),kTProfile,{{nBinsAmpFT0,0.,maxAmpFT0}}); + // + // pOneParCorrVsV0A[i] = registry.add(Form("OneParCorrVsV0A_Replica%d",i),Form(";%s;%s;",tiV0A,tiOneParCorr),kTProfile,{{nBinsAmpV0A,0.,maxAmpV0A}}); + // pTwoParCorrVsV0A[i] = registry.add(Form("TwoParCorrVsV0A_Replica%d",i),Form(";%s;%s;",tiV0A,tiTwoParCorr),kTProfile,{{nBinsAmpV0A,0.,maxAmpV0A}}); + // pThreeParCorrVsV0A[i] = registry.add(Form("ThreeParCorrVsV0A_Replica%d",i),Form(";%s;%s;",tiV0A,tiThreeParCorr),kTProfile,{{nBinsAmpV0A,0.,maxAmpV0A}}); } } @@ -756,7 +778,7 @@ struct UccZdc { } // Nch-based selection - int glbTracks = 0; + int glbTracks{0}; for (const auto& track : tracks) { // Track Selection if (!track.isGlobalTrack()) { @@ -809,7 +831,7 @@ struct UccZdc { return; } - double nchMult{0.}; + double nchMult{static_cast(glbTracks)}; std::vector pTs; std::vector vecFD; std::vector vecEff; @@ -953,165 +975,9 @@ struct UccZdc { } PROCESS_SWITCH(UccZdc, processZdcCollAss, "Process ZDC W/Coll Ass.", true); - template - void eventSampling(const T& tracks, const U& normV0A, const U& normT0M, const U& sumZNs, const V& timeStamp) - { - TRandom3 rndGen(timeStamp); - std::vector vPoisson; - for (int replica = 0; replica < kSizeBootStrapEnsemble; ++replica) - vPoisson.emplace_back(rndGen.Poisson(1.)); - - for (int replica = 0; replica < kSizeBootStrapEnsemble; ++replica) { - - hPoisson[replica]->Fill(vPoisson.at(replica)); - - for (uint64_t evtRep = 0; evtRep < vPoisson.at(replica); ++evtRep) { - - double nchMult{0.}; - int glbTracks{0}; - std::vector pTs; - std::vector vecFD; - std::vector vecEff; - - // Calculates the uncorrected Nch multiplicity - for (const auto& track : tracks) { - // Track Selection - if (!track.isGlobalTrack()) { - continue; - } - if ((track.pt() < minPt) || (track.pt() > maxPt)) { - continue; - } - if ((track.eta() < minEta) || (track.eta() > maxEta)) { - continue; - } - glbTracks++; - } - - if (glbTracks < minNchSel) { - continue; - } - - // Calculates the Nch multiplicity if corrections are loaded - if (cfg.correctionsLoaded) { - const int foundNchBin{cfg.hEfficiency->GetXaxis()->FindBin(glbTracks)}; - for (const auto& track : tracks) { - // Track Selection - if (!track.isGlobalTrack()) { - continue; - } - if ((track.pt() < minPt) || (track.pt() > maxPt)) { - continue; - } - if ((track.eta() < minEta) || (track.eta() > maxEta)) { - continue; - } - - float pt{track.pt()}; - double fdValue{1.}; - int foundPtBin{cfg.hEfficiency->GetYaxis()->FindBin(pt)}; - double effValue{cfg.hEfficiency->GetBinContent(foundNchBin, foundPtBin)}; - - if (applyFD) - fdValue = cfg.hFeedDown->GetBinContent(foundNchBin, foundPtBin); - if ((effValue > 0.) && (fdValue > 0.)) { - nchMult += (std::pow(effValue, -1.) * fdValue); - } - } - } - - if (!applyEff) - nchMult = static_cast(glbTracks); - if (applyEff && !correctNch) - nchMult = static_cast(glbTracks); - - // Fill vectors for [pT] measurement - if (cfg.correctionsLoaded) { - const int foundNchBin{cfg.hEfficiency->GetXaxis()->FindBin(glbTracks)}; - // Fill vectors for [pT] measurement - for (const auto& track : tracks) { - // Track Selection - if (!track.isGlobalTrack()) { - continue; - } - if ((track.pt() < minPt) || (track.pt() > maxPtSpectra)) { - continue; - } - if ((track.eta() < minEta) || (track.eta() > maxEta)) { - continue; - } - - float pt{track.pt()}; - double fdValue{1.}; - int foundPtBin{cfg.hEfficiency->GetYaxis()->FindBin(pt)}; - double effValue{cfg.hEfficiency->GetBinContent(foundNchBin, foundPtBin)}; - - if (applyFD) - fdValue = cfg.hFeedDown->GetBinContent(foundNchBin, foundPtBin); - - if ((effValue > 0.) && (fdValue > 0.)) { - pTs.emplace_back(pt); - vecEff.emplace_back(effValue); - vecFD.emplace_back(fdValue); - // To calculate event-averaged - registry.fill(HIST("NchVsZNVsPt"), nchMult, sumZNs, pt * (fdValue / effValue)); - } - } - } else { - for (const auto& track : tracks) { - // Track Selection - if (!track.isGlobalTrack()) { - continue; - } - if ((track.pt() < minPt) || (track.pt() > maxPtSpectra)) { - continue; - } - - pTs.emplace_back(track.pt()); - vecEff.emplace_back(1.); - vecFD.emplace_back(1.); - - // To calculate event-averaged - registry.fill(HIST("NchVsZNVsPt"), nchMult, sumZNs, track.pt()); - } - } - - double p1, p2, p3, p4, w1, w2, w3, w4; - p1 = p2 = p3 = p4 = w1 = w2 = w3 = w4 = 0.0; - getPTpowers(pTs, vecEff, vecFD, p1, w1, p2, w2, p3, w3, p4, w4); - - // EbE one-particle pT correlation - const double oneParCorr{p1 / w1}; - - // EbE two-particle pT correlation - const double denTwoParCorr{std::pow(w1, 2.) - w2}; - const double numTwoParCorr{std::pow(p1, 2.) - p2}; - const double twoParCorr{numTwoParCorr / denTwoParCorr}; - - // EbE three-particle pT correlation - const double denThreeParCorr{std::pow(w1, 3.) - 3. * w2 * w1 + 2. * w3}; - const double numThreeParCorr{std::pow(p1, 3.) - 3. * p2 * p1 + 2. * p3}; - const double threeParCorr{numThreeParCorr / denThreeParCorr}; - - hNch[replica]->Fill(nchMult); - pNchVsOneParCorrVsZN[replica]->Fill(nchMult, sumZNs, oneParCorr, w1); - pNchVsTwoParCorrVsZN[replica]->Fill(nchMult, sumZNs, twoParCorr, denTwoParCorr); - pNchVsThreeParCorrVsZN[replica]->Fill(nchMult, sumZNs, threeParCorr, denThreeParCorr); - - pOneParCorrVsV0A[replica]->Fill(normV0A, oneParCorr, w1); - pTwoParCorrVsV0A[replica]->Fill(normV0A, twoParCorr, denTwoParCorr); - pThreeParCorrVsV0A[replica]->Fill(normV0A, threeParCorr, denThreeParCorr); - - pOneParCorrVsT0M[replica]->Fill(normT0M, oneParCorr, w1); - pTwoParCorrVsT0M[replica]->Fill(normT0M, twoParCorr, denTwoParCorr); - pThreeParCorrVsT0M[replica]->Fill(normT0M, threeParCorr, denThreeParCorr); - } // event per replica - } // replica's loop - } - Preslice perCollision = aod::track::collisionId; Service pdg; - void processMCclosure(aod::McCollisions::iterator const& mccollision, soa::SmallGroups const& collisions, o2::aod::BCsRun3 const& /*bcs*/, aod::FT0s const& /*ft0s*/, aod::McParticles const& mcParticles, TheFilteredSimTracks const& simTracks) + void processMCclosure(aod::McCollisions::iterator const& mccollision, soa::SmallGroups const& collisions, o2::aod::BCsRun3 const& /*bcs*/, aod::FT0s const& /*ft0s*/, aod::FV0As const& /*fv0as*/, aod::McParticles const& mcParticles, TheFilteredSimTracks const& simTracks) { for (const auto& collision : collisions) { // Event selection @@ -1136,7 +1002,7 @@ struct UccZdc { const double rndNum{rndGen.Uniform(0.0, 1.0)}; registry.fill(HIST("RandomNumber"), rndNum); - float aT0A = 0., aT0C = 0.; + double aT0A = 0., aT0C = 0.; if (foundBC.has_ft0()) { for (const auto& amplitude : foundBC.ft0().amplitudeA()) { aT0A += amplitude; @@ -1148,11 +1014,17 @@ struct UccZdc { return; } + // double aV0A{-999.}; + // if (foundBC.has_fv0a()) { + // for (const auto& amplitude : foundBC.fv0a().amplitude()) { aV0A += amplitude; } + // } + + const double normT0M{(aT0A + aT0C) / 100.}; + // const double normV0A{aV0A/100.}; + double nchRaw{0.}; double nchMult{0.}; double nchMC{0.}; - double normT0M{0.}; - normT0M = (aT0A + aT0C) / 100.; registry.fill(HIST("zPos"), collision.posZ()); registry.fill(HIST("zPosMC"), mccollision.posZ()); @@ -1174,15 +1046,10 @@ struct UccZdc { registry.fill(HIST("EvtsDivided"), 0); - // To use run-by-run efficiency - auto efficiency = ccdb->getForTimeStamp(paTHEff.value, foundBC.timestamp()); - if (!efficiency) { - return; - } - - auto feedDown = ccdb->getForTimeStamp(paTHFD.value, foundBC.timestamp()); - if (!feedDown) { - return; + // Run-by-run efficiency + loadCorrections(foundBC.timestamp()); + if (!(cfg.hEfficiency && cfg.hFeedDown)) { + continue; } std::vector pTs; @@ -1204,8 +1071,13 @@ struct UccZdc { nchRaw++; } + // Reject event if nchRaw less than a lower cutoff + if (nchRaw < minNchSel) { + return; + } + // Calculates the event weight, W_k - const int foundNchBin{efficiency->GetXaxis()->FindBin(nchRaw)}; + const int foundNchBin{cfg.hEfficiency->GetXaxis()->FindBin(nchRaw)}; for (const auto& track : groupedTracks) { // Track Selection @@ -1240,14 +1112,12 @@ struct UccZdc { // if (!particle.isPhysicalPrimary()) { continue; } const double pt{static_cast(track.pt())}; - const int foundPtBin{efficiency->GetYaxis()->FindBin(pt)}; - double effValue{1.}; + const int foundPtBin{cfg.hEfficiency->GetYaxis()->FindBin(pt)}; + const double effValue{cfg.hEfficiency->GetBinContent(foundNchBin, foundPtBin)}; double fdValue{1.}; - if (applyEff) { - effValue = efficiency->GetBinContent(foundNchBin, foundPtBin); - fdValue = feedDown->GetBinContent(foundNchBin, foundPtBin); - } + if (applyFD) + fdValue = cfg.hFeedDown->GetBinContent(foundNchBin, foundPtBin); if ((effValue > 0.) && (fdValue > 0.)) { pTs.emplace_back(pt); vecEff.emplace_back(effValue); @@ -1256,10 +1126,6 @@ struct UccZdc { } } - if (nchMult < minNchSel) { - return; - } - double p1, p2, p3, p4, w1, w2, w3, w4; p1 = p2 = p3 = p4 = w1 = w2 = w3 = w4 = 0.0; getPTpowers(pTs, vecEff, vecFD, p1, w1, p2, w2, p3, w3, p4, w4); @@ -1268,20 +1134,18 @@ struct UccZdc { const double numTwoParCorr{std::pow(p1, 2.) - p2}; const double denThreeParCorr{std::pow(w1, 3.) - 3. * w2 * w1 + 2. * w3}; const double numThreeParCorr{std::pow(p1, 3.) - 3. * p2 * p1 + 2. * p3}; - const double denFourParCorr{std::pow(w1, 4.) - 6. * w2 * std::pow(w1, 2.) + 3. * std::pow(w2, 2.) + 8 * w3 * w1 - 6. * w4}; - const double numFourParCorr{std::pow(p1, 4.) - 6. * p2 * std::pow(p1, 2.) + 3. * std::pow(p2, 2.) + 8 * p3 * p1 - 6. * p4}; + // const double denFourParCorr{std::pow(w1, 4.) - 6. * w2 * std::pow(w1, 2.) + 3. * std::pow(w2, 2.) + 8 * w3 * w1 - 6. * w4}; + // const double numFourParCorr{std::pow(p1, 4.) - 6. * p2 * std::pow(p1, 2.) + 3. * std::pow(p2, 2.) + 8 * p3 * p1 - 6. * p4}; const double oneParCorr{p1 / w1}; const double twoParCorr{numTwoParCorr / denTwoParCorr}; const double threeParCorr{numThreeParCorr / denThreeParCorr}; - const double fourParCorr{numFourParCorr / denFourParCorr}; registry.fill(HIST("Nch"), nchMult); registry.fill(HIST("NchUncorrected"), nchRaw); registry.fill(HIST("NchVsOneParCorr"), nchMult, oneParCorr, w1); registry.fill(HIST("NchVsTwoParCorr"), nchMult, twoParCorr, denTwoParCorr); registry.fill(HIST("NchVsThreeParCorr"), nchMult, threeParCorr, denThreeParCorr); - registry.fill(HIST("NchVsFourParCorr"), nchMult, fourParCorr, denFourParCorr); //--------------------------- Generated MC --------------------------- std::vector pTsMC; @@ -1335,200 +1199,22 @@ struct UccZdc { const double numTwoParCorrMC{std::pow(p1MC, 2.) - p2MC}; const double denThreeParCorrMC{std::pow(w1MC, 3.) - 3. * w2MC * w1MC + 2. * w3MC}; const double numThreeParCorrMC{std::pow(p1MC, 3.) - 3. * p2MC * p1MC + 2. * p3MC}; - const double denFourParCorrMC{std::pow(w1MC, 4.) - 6. * w2MC * std::pow(w1MC, 2.) + 3. * std::pow(w2MC, 2.) + 8 * w3MC * w1MC - 6. * w4MC}; - const double numFourParCorrMC{std::pow(p1MC, 4.) - 6. * p2MC * std::pow(p1MC, 2.) + 3. * std::pow(p2MC, 2.) + 8 * p3MC * p1MC - 6. * p4MC}; + // const double denFourParCorrMC{std::pow(w1MC, 4.) - 6. * w2MC * std::pow(w1MC, 2.) + 3. * std::pow(w2MC, 2.) + 8 * w3MC * w1MC - 6. * w4MC}; + // const double numFourParCorrMC{std::pow(p1MC, 4.) - 6. * p2MC * std::pow(p1MC, 2.) + 3. * std::pow(p2MC, 2.) + 8 * p3MC * p1MC - 6. * p4MC}; const double oneParCorrMC{p1MC / w1MC}; const double twoParCorrMC{numTwoParCorrMC / denTwoParCorrMC}; const double threeParCorrMC{numThreeParCorrMC / denThreeParCorrMC}; - const double fourParCorrMC{numFourParCorrMC / denFourParCorrMC}; + // const double fourParCorrMC{numFourParCorrMC / denFourParCorrMC}; registry.fill(HIST("NchGen"), nchMC); registry.fill(HIST("NchvsOneParCorrGen"), nchMC, oneParCorrMC, w1MC); registry.fill(HIST("NchvsTwoParCorrGen"), nchMC, twoParCorrMC, denTwoParCorrMC); registry.fill(HIST("NchvsThreeParCorrGen"), nchMC, threeParCorrMC, denThreeParCorrMC); - registry.fill(HIST("NchvsFourParCorrGen"), nchMC, fourParCorrMC, denFourParCorrMC); //------------------ Poisson sampling - std::vector vPoisson; - for (int replica = 0; replica < kSizeBootStrapEnsemble; ++replica) { - vPoisson.emplace_back(rndGen.Poisson(1.)); - } - - for (int replica = 0; replica < kSizeBootStrapEnsemble; ++replica) { - hPoissonMC[replica]->Fill(vPoisson.at(replica)); - - for (uint64_t evtRep = 0; evtRep < vPoisson.at(replica); ++evtRep) { - - double nchRaw{0.0}; - double nchMult{0.0}; - std::vector pTs; - std::vector vecFD; - std::vector vecEff; - - // const auto& groupedTracks{simTracks.sliceBy(perCollision, collision.globalIndex())}; - - // Calculates the event's Nch to evaluate the efficiency - for (const auto& track : groupedTracks) { - // Track Selection - if (track.eta() < minEta || track.eta() > maxEta) { - continue; - } - if (track.pt() < minPt || track.pt() > maxPt) { - continue; - } - if (!track.isGlobalTrack()) { - continue; - } - nchRaw++; - } - - // Calculates the event weight, W_k - const int foundNchBin{efficiency->GetXaxis()->FindBin(nchRaw)}; - - for (const auto& track : groupedTracks) { - // Track Selection - if (track.eta() < minEta || track.eta() > maxEta) { - continue; - } - if (track.pt() < minPt || track.pt() > maxPt) { - continue; - } - if (!track.isGlobalTrack()) { - continue; - } - if (!track.has_mcParticle()) { - continue; - } - const auto& particle{track.mcParticle()}; - - auto charge{0.}; - // Get the MC particle - auto* pdgParticle = pdg->GetParticle(particle.pdgCode()); - if (pdgParticle != nullptr) { - charge = pdgParticle->Charge(); - } else { - continue; - } - - // Is it a charged particle? - if (std::abs(charge) < kMinCharge) { - continue; - } - // Is it a primary particle? - // if (!particle.isPhysicalPrimary()) { continue; } - - const double pt{static_cast(track.pt())}; - const int foundPtBin{efficiency->GetYaxis()->FindBin(pt)}; - double effValue{1.}; - double fdValue{1.}; - - if (applyEff) { - effValue = efficiency->GetBinContent(foundNchBin, foundPtBin); - fdValue = feedDown->GetBinContent(foundNchBin, foundPtBin); - } - if ((effValue > 0.) && (fdValue > 0.)) { - pTs.emplace_back(pt); - vecEff.emplace_back(effValue); - vecFD.emplace_back(fdValue); - nchMult += (std::pow(effValue, -1.0) * fdValue); - } - } - - if (nchMult < minNchSel) { - return; - } - - double p1, p2, p3, p4, w1, w2, w3, w4; - p1 = p2 = p3 = p4 = w1 = w2 = w3 = w4 = 0.0; - getPTpowers(pTs, vecEff, vecFD, p1, w1, p2, w2, p3, w3, p4, w4); - - const double denTwoParCorr{std::pow(w1, 2.) - w2}; - const double numTwoParCorr{std::pow(p1, 2.) - p2}; - const double denThreeParCorr{std::pow(w1, 3.) - 3. * w2 * w1 + 2. * w3}; - const double numThreeParCorr{std::pow(p1, 3.) - 3. * p2 * p1 + 2. * p3}; - const double denFourParCorr{std::pow(w1, 4.) - 6. * w2 * std::pow(w1, 2.) + 3. * std::pow(w2, 2.) + 8 * w3 * w1 - 6. * w4}; - const double numFourParCorr{std::pow(p1, 4.) - 6. * p2 * std::pow(p1, 2.) + 3. * std::pow(p2, 2.) + 8 * p3 * p1 - 6. * p4}; - - const double oneParCorr{p1 / w1}; - const double twoParCorr{numTwoParCorr / denTwoParCorr}; - const double threeParCorr{numThreeParCorr / denThreeParCorr}; - const double fourParCorr{numFourParCorr / denFourParCorr}; - - hNchRec[replica]->Fill(nchMult); - pNchvsOneParCorrRec[replica]->Fill(nchMult, oneParCorr, w1); - pNchvsTwoParCorrRec[replica]->Fill(nchMult, twoParCorr, denTwoParCorr); - pNchvsThreeParCorrRec[replica]->Fill(nchMult, threeParCorr, denThreeParCorr); - pNchvsFourParCorrRec[replica]->Fill(nchMult, fourParCorr, denFourParCorr); - - //--------------------------- Generated MC --------------------------- - double nchMC{0.0}; - std::vector pTsMC; - std::vector vecFullEff; - std::vector vecFDEqualOne; - - // Calculates the event weight, W_k - for (const auto& particle : mcParticles) { - if (particle.eta() < minEta || particle.eta() > maxEta) { - continue; - } - if (particle.pt() < minPt || particle.pt() > maxPt) { - continue; - } - - auto charge{0.}; - // Get the MC particle - auto* pdgParticle = pdg->GetParticle(particle.pdgCode()); - if (pdgParticle != nullptr) { - charge = pdgParticle->Charge(); - } else { - continue; - } - - // Is it a charged particle? - if (std::abs(charge) < kMinCharge) { - continue; - } - // Is it a primary particle? - if (!particle.isPhysicalPrimary()) { - continue; - } - - float pt{particle.pt()}; - pTsMC.emplace_back(pt); - vecFullEff.emplace_back(1.); - vecFDEqualOne.emplace_back(1.); - nchMC++; - } - - if (nchMC < minNchSel) { - continue; - } - // printf("nchMult = %f | nchMC = %f | nchMult/nchMc = %f\n",nchMult,nchMC,nchMult/nchMC); - - double p1MC, p2MC, p3MC, p4MC, w1MC, w2MC, w3MC, w4MC; - p1MC = p2MC = p3MC = p4MC = w1MC = w2MC = w3MC = w4MC = 0.0; - getPTpowers(pTsMC, vecFullEff, vecFDEqualOne, p1MC, w1MC, p2MC, w2MC, p3MC, w3MC, p4MC, w4MC); - - const double denTwoParCorrMC{std::pow(w1MC, 2.) - w2MC}; - const double numTwoParCorrMC{std::pow(p1MC, 2.) - p2MC}; - const double denThreeParCorrMC{std::pow(w1MC, 3.) - 3. * w2MC * w1MC + 2. * w3MC}; - const double numThreeParCorrMC{std::pow(p1MC, 3.) - 3. * p2MC * p1MC + 2. * p3MC}; - const double denFourParCorrMC{std::pow(w1MC, 4.) - 6. * w2MC * std::pow(w1MC, 2.) + 3. * std::pow(w2MC, 2.) + 8 * w3MC * w1MC - 6. * w4MC}; - const double numFourParCorrMC{std::pow(p1MC, 4.) - 6. * p2MC * std::pow(p1MC, 2.) + 3. * std::pow(p2MC, 2.) + 8 * p3MC * p1MC - 6. * p4MC}; - - const double oneParCorrMC{p1MC / w1MC}; - const double twoParCorrMC{numTwoParCorrMC / denTwoParCorrMC}; - const double threeParCorrMC{numThreeParCorrMC / denThreeParCorrMC}; - const double fourParCorrMC{numFourParCorrMC / denFourParCorrMC}; - - hNchGen[replica]->Fill(nchMC); - pNchvsOneParCorrGen[replica]->Fill(nchMC, oneParCorrMC, w1MC); - pNchvsTwoParCorrGen[replica]->Fill(nchMC, twoParCorrMC, w1MC); - pNchvsThreeParCorrGen[replica]->Fill(nchMC, threeParCorrMC, w1MC); - pNchvsFourParCorrGen[replica]->Fill(nchMC, fourParCorrMC, w1MC); - } // events per replica - } // replica's loop + eventSamplingMC(mcParticles, timeStamp); + eventSamplingMCRec(groupedTracks, timeStamp); } else { // Correction with the remaining half of the sample registry.fill(HIST("EvtsDivided"), 1); //----- MC reconstructed -----// @@ -1670,6 +1356,406 @@ struct UccZdc { } } + template + void eventSamplingMC(const T& mcParticles, const V& timeStamp) + { + TRandom3 rndGen(timeStamp); + std::vector vPoisson; + for (int replica = 0; replica < kSizeBootStrapEnsemble; ++replica) + vPoisson.emplace_back(rndGen.Poisson(1.)); + + for (int replica = 0; replica < kSizeBootStrapEnsemble; ++replica) { + + hPoissonMC[replica]->Fill(vPoisson.at(replica)); + + for (uint64_t evtRep = 0; evtRep < vPoisson.at(replica); ++evtRep) { + + double nchMult{0.}; + std::vector pTs; + std::vector vecFD; + std::vector vecEff; + + // Calculates the event weight, W_k + for (const auto& particle : mcParticles) { + if (particle.eta() < minEta || particle.eta() > maxEta) { + continue; + } + if (particle.pt() < minPt || particle.pt() > maxPt) { + continue; + } + + auto charge{0.}; + // Get the MC particle + auto* pdgParticle = pdg->GetParticle(particle.pdgCode()); + if (pdgParticle != nullptr) { + charge = pdgParticle->Charge(); + } else { + continue; + } + + // Is it a charged particle? + if (std::abs(charge) < kMinCharge) { + continue; + } + // Is it a primary particle? + if (!particle.isPhysicalPrimary()) { + continue; + } + + float pt{particle.pt()}; + pTs.emplace_back(pt); + vecEff.emplace_back(1.); + vecFD.emplace_back(1.); + nchMult++; + } + + if (nchMult < minNchSel) { + continue; + } + // printf("nchMult = %f | nchMC = %f | nchMult/nchMc = %f\n",nchMult,nchMC,nchMult/nchMC); + + double p1, p2, p3, p4, w1, w2, w3, w4; + p1 = p2 = p3 = p4 = w1 = w2 = w3 = w4 = 0.0; + getPTpowers(pTs, vecEff, vecFD, p1, w1, p2, w2, p3, w3, p4, w4); + + // EbE one-particle pT correlation + const double oneParCorr{p1 / w1}; + + // EbE two-particle pT correlation + const double denTwoParCorr{std::pow(w1, 2.) - w2}; + const double numTwoParCorr{std::pow(p1, 2.) - p2}; + const double twoParCorr{numTwoParCorr / denTwoParCorr}; + + // EbE three-particle pT correlation + const double denThreeParCorr{std::pow(w1, 3.) - 3. * w2 * w1 + 2. * w3}; + const double numThreeParCorr{std::pow(p1, 3.) - 3. * p2 * p1 + 2. * p3}; + const double threeParCorr{numThreeParCorr / denThreeParCorr}; + + hNchGen[replica]->Fill(nchMult); + + pOneParCorrVsNchGen[replica]->Fill(nchMult, oneParCorr, w1); + pTwoParCorrVsNchGen[replica]->Fill(nchMult, twoParCorr, denTwoParCorr); + pThreeParCorrVsNchGen[replica]->Fill(nchMult, threeParCorr, denThreeParCorr); + + // pOneParCorrVsV0AGen[replica]->Fill(normV0A, oneParCorr, w1); + // pTwoParCorrVsV0AGen[replica]->Fill(normV0A, twoParCorr, denTwoParCorr); + // pThreeParCorrVsV0AGen[replica]->Fill(normV0A, threeParCorr, denThreeParCorr); + + // pOneParCorrVsT0MGen[replica]->Fill(normT0M, oneParCorr, w1); + // pTwoParCorrVsT0MGen[replica]->Fill(normT0M, twoParCorr, denTwoParCorr); + // pThreeParCorrVsT0MGen[replica]->Fill(normT0M, threeParCorr, denThreeParCorr); + } // event per replica + } // replica's loop + } + + template + void eventSamplingMCRec(const T& tracks, const V& timeStamp) + { + TRandom3 rndGen(timeStamp); + std::vector vPoisson; + for (int replica = 0; replica < kSizeBootStrapEnsemble; ++replica) + vPoisson.emplace_back(rndGen.Poisson(1.)); + + for (int replica = 0; replica < kSizeBootStrapEnsemble; ++replica) { + + hPoisson[replica]->Fill(vPoisson.at(replica)); + + for (uint64_t evtRep = 0; evtRep < vPoisson.at(replica); ++evtRep) { + + double nchMult{0.}; + int glbTracks{0}; + std::vector pTs; + std::vector vecFD; + std::vector vecEff; + + // Calculates the uncorrected Nch multiplicity + for (const auto& track : tracks) { + // Track Selection + if (!track.isGlobalTrack()) { + continue; + } + if ((track.pt() < minPt) || (track.pt() > maxPt)) { + continue; + } + if ((track.eta() < minEta) || (track.eta() > maxEta)) { + continue; + } + glbTracks++; + } + + if (glbTracks < minNchSel) { + continue; + } + + // Calculates the Nch multiplicity if corrections are loaded + if (cfg.correctionsLoaded) { + const int foundNchBin{cfg.hEfficiency->GetXaxis()->FindBin(glbTracks)}; + for (const auto& track : tracks) { + // Track Selection + if (!track.isGlobalTrack()) { + continue; + } + if ((track.pt() < minPt) || (track.pt() > maxPt)) { + continue; + } + if ((track.eta() < minEta) || (track.eta() > maxEta)) { + continue; + } + + float pt{track.pt()}; + double fdValue{1.}; + int foundPtBin{cfg.hEfficiency->GetYaxis()->FindBin(pt)}; + double effValue{cfg.hEfficiency->GetBinContent(foundNchBin, foundPtBin)}; + + if (applyFD) + fdValue = cfg.hFeedDown->GetBinContent(foundNchBin, foundPtBin); + if ((effValue > 0.) && (fdValue > 0.)) { + nchMult += (std::pow(effValue, -1.) * fdValue); + } + } + } + + if (!applyEff) + nchMult = static_cast(glbTracks); + if (applyEff && !correctNch) + nchMult = static_cast(glbTracks); + + // Fill vectors for [pT] measurement + if (cfg.correctionsLoaded) { + const int foundNchBin{cfg.hEfficiency->GetXaxis()->FindBin(glbTracks)}; + // Fill vectors for [pT] measurement + for (const auto& track : tracks) { + // Track Selection + if (!track.isGlobalTrack()) { + continue; + } + if ((track.pt() < minPt) || (track.pt() > maxPtSpectra)) { + continue; + } + if ((track.eta() < minEta) || (track.eta() > maxEta)) { + continue; + } + + float pt{track.pt()}; + double fdValue{1.}; + int foundPtBin{cfg.hEfficiency->GetYaxis()->FindBin(pt)}; + double effValue{cfg.hEfficiency->GetBinContent(foundNchBin, foundPtBin)}; + + if (applyFD) + fdValue = cfg.hFeedDown->GetBinContent(foundNchBin, foundPtBin); + + if ((effValue > 0.) && (fdValue > 0.)) { + pTs.emplace_back(pt); + vecEff.emplace_back(effValue); + vecFD.emplace_back(fdValue); + } + } + } else { + for (const auto& track : tracks) { + // Track Selection + if (!track.isGlobalTrack()) { + continue; + } + if ((track.pt() < minPt) || (track.pt() > maxPtSpectra)) { + continue; + } + + pTs.emplace_back(track.pt()); + vecEff.emplace_back(1.); + vecFD.emplace_back(1.); + } + } + + double p1, p2, p3, p4, w1, w2, w3, w4; + p1 = p2 = p3 = p4 = w1 = w2 = w3 = w4 = 0.0; + getPTpowers(pTs, vecEff, vecFD, p1, w1, p2, w2, p3, w3, p4, w4); + + // EbE one-particle pT correlation + const double oneParCorr{p1 / w1}; + + // EbE two-particle pT correlation + const double denTwoParCorr{std::pow(w1, 2.) - w2}; + const double numTwoParCorr{std::pow(p1, 2.) - p2}; + const double twoParCorr{numTwoParCorr / denTwoParCorr}; + + // EbE three-particle pT correlation + const double denThreeParCorr{std::pow(w1, 3.) - 3. * w2 * w1 + 2. * w3}; + const double numThreeParCorr{std::pow(p1, 3.) - 3. * p2 * p1 + 2. * p3}; + const double threeParCorr{numThreeParCorr / denThreeParCorr}; + + hNch[replica]->Fill(nchMult); + + pOneParCorrVsNch[replica]->Fill(nchMult, oneParCorr, w1); + pTwoParCorrVsNch[replica]->Fill(nchMult, twoParCorr, denTwoParCorr); + pThreeParCorrVsNch[replica]->Fill(nchMult, threeParCorr, denThreeParCorr); + + // pOneParCorrVsV0A[replica]->Fill(normV0A, oneParCorr, w1); + // pTwoParCorrVsV0A[replica]->Fill(normV0A, twoParCorr, denTwoParCorr); + // pThreeParCorrVsV0A[replica]->Fill(normV0A, threeParCorr, denThreeParCorr); + // + // pOneParCorrVsT0M[replica]->Fill(normT0M, oneParCorr, w1); + // pTwoParCorrVsT0M[replica]->Fill(normT0M, twoParCorr, denTwoParCorr); + // pThreeParCorrVsT0M[replica]->Fill(normT0M, threeParCorr, denThreeParCorr); + } // event per replica + } // replica's loop + } + + template + void eventSampling(const T& tracks, const U& normV0A, const U& normT0M, const U& sumZNs, const V& timeStamp) + { + TRandom3 rndGen(timeStamp); + std::vector vPoisson; + for (int replica = 0; replica < kSizeBootStrapEnsemble; ++replica) + vPoisson.emplace_back(rndGen.Poisson(1.)); + + for (int replica = 0; replica < kSizeBootStrapEnsemble; ++replica) { + + hPoisson[replica]->Fill(vPoisson.at(replica)); + + for (uint64_t evtRep = 0; evtRep < vPoisson.at(replica); ++evtRep) { + + double nchMult{0.}; + int glbTracks{0}; + std::vector pTs; + std::vector vecFD; + std::vector vecEff; + + // Calculates the uncorrected Nch multiplicity + for (const auto& track : tracks) { + // Track Selection + if (!track.isGlobalTrack()) { + continue; + } + if ((track.pt() < minPt) || (track.pt() > maxPt)) { + continue; + } + if ((track.eta() < minEta) || (track.eta() > maxEta)) { + continue; + } + glbTracks++; + } + + if (glbTracks < minNchSel) { + continue; + } + + // Calculates the Nch multiplicity if corrections are loaded + if (cfg.correctionsLoaded) { + const int foundNchBin{cfg.hEfficiency->GetXaxis()->FindBin(glbTracks)}; + for (const auto& track : tracks) { + // Track Selection + if (!track.isGlobalTrack()) { + continue; + } + if ((track.pt() < minPt) || (track.pt() > maxPt)) { + continue; + } + if ((track.eta() < minEta) || (track.eta() > maxEta)) { + continue; + } + + float pt{track.pt()}; + double fdValue{1.}; + int foundPtBin{cfg.hEfficiency->GetYaxis()->FindBin(pt)}; + double effValue{cfg.hEfficiency->GetBinContent(foundNchBin, foundPtBin)}; + + if (applyFD) + fdValue = cfg.hFeedDown->GetBinContent(foundNchBin, foundPtBin); + if ((effValue > 0.) && (fdValue > 0.)) { + nchMult += (std::pow(effValue, -1.) * fdValue); + } + } + } + + if (!applyEff) + nchMult = static_cast(glbTracks); + if (applyEff && !correctNch) + nchMult = static_cast(glbTracks); + + // Fill vectors for [pT] measurement + if (cfg.correctionsLoaded) { + const int foundNchBin{cfg.hEfficiency->GetXaxis()->FindBin(glbTracks)}; + // Fill vectors for [pT] measurement + for (const auto& track : tracks) { + // Track Selection + if (!track.isGlobalTrack()) { + continue; + } + if ((track.pt() < minPt) || (track.pt() > maxPtSpectra)) { + continue; + } + if ((track.eta() < minEta) || (track.eta() > maxEta)) { + continue; + } + + float pt{track.pt()}; + double fdValue{1.}; + int foundPtBin{cfg.hEfficiency->GetYaxis()->FindBin(pt)}; + double effValue{cfg.hEfficiency->GetBinContent(foundNchBin, foundPtBin)}; + + if (applyFD) + fdValue = cfg.hFeedDown->GetBinContent(foundNchBin, foundPtBin); + + if ((effValue > 0.) && (fdValue > 0.)) { + pTs.emplace_back(pt); + vecEff.emplace_back(effValue); + vecFD.emplace_back(fdValue); + // To calculate event-averaged + registry.fill(HIST("NchVsZNVsPt"), nchMult, sumZNs, pt * (fdValue / effValue)); + } + } + } else { + for (const auto& track : tracks) { + // Track Selection + if (!track.isGlobalTrack()) { + continue; + } + if ((track.pt() < minPt) || (track.pt() > maxPtSpectra)) { + continue; + } + + pTs.emplace_back(track.pt()); + vecEff.emplace_back(1.); + vecFD.emplace_back(1.); + + // To calculate event-averaged + registry.fill(HIST("NchVsZNVsPt"), nchMult, sumZNs, track.pt()); + } + } + + double p1, p2, p3, p4, w1, w2, w3, w4; + p1 = p2 = p3 = p4 = w1 = w2 = w3 = w4 = 0.0; + getPTpowers(pTs, vecEff, vecFD, p1, w1, p2, w2, p3, w3, p4, w4); + + // EbE one-particle pT correlation + const double oneParCorr{p1 / w1}; + + // EbE two-particle pT correlation + const double denTwoParCorr{std::pow(w1, 2.) - w2}; + const double numTwoParCorr{std::pow(p1, 2.) - p2}; + const double twoParCorr{numTwoParCorr / denTwoParCorr}; + + // EbE three-particle pT correlation + const double denThreeParCorr{std::pow(w1, 3.) - 3. * w2 * w1 + 2. * w3}; + const double numThreeParCorr{std::pow(p1, 3.) - 3. * p2 * p1 + 2. * p3}; + const double threeParCorr{numThreeParCorr / denThreeParCorr}; + + hNch[replica]->Fill(nchMult); + pNchVsOneParCorrVsZN[replica]->Fill(nchMult, sumZNs, oneParCorr, w1); + pNchVsTwoParCorrVsZN[replica]->Fill(nchMult, sumZNs, twoParCorr, denTwoParCorr); + pNchVsThreeParCorrVsZN[replica]->Fill(nchMult, sumZNs, threeParCorr, denThreeParCorr); + + pOneParCorrVsV0A[replica]->Fill(normV0A, oneParCorr, w1); + pTwoParCorrVsV0A[replica]->Fill(normV0A, twoParCorr, denTwoParCorr); + pThreeParCorrVsV0A[replica]->Fill(normV0A, threeParCorr, denThreeParCorr); + + pOneParCorrVsT0M[replica]->Fill(normT0M, oneParCorr, w1); + pTwoParCorrVsT0M[replica]->Fill(normT0M, twoParCorr, denTwoParCorr); + pThreeParCorrVsT0M[replica]->Fill(normT0M, threeParCorr, denThreeParCorr); + } // event per replica + } // replica's loop + } + void loadCorrections(uint64_t timeStamp) { // if (cfg.correctionsLoaded) return;