From 8a7610dc90a000c5fc77e9f1956d86159f75a7e4 Mon Sep 17 00:00:00 2001 From: Francesca Ercolessi Date: Tue, 7 Oct 2025 17:03:52 +0200 Subject: [PATCH 1/6] add generated denominator and z vtx cut in gen collisins --- PWGLF/Tasks/Nuspex/hadronnucleicorrelation.cxx | 18 ++++++++++++++++++ 1 file changed, 18 insertions(+) diff --git a/PWGLF/Tasks/Nuspex/hadronnucleicorrelation.cxx b/PWGLF/Tasks/Nuspex/hadronnucleicorrelation.cxx index ee383900901..424f61b8583 100644 --- a/PWGLF/Tasks/Nuspex/hadronnucleicorrelation.cxx +++ b/PWGLF/Tasks/Nuspex/hadronnucleicorrelation.cxx @@ -158,6 +158,9 @@ struct hadronnucleicorrelation { std::map, std::vector> mixbinsPID_antidantip; std::map> mixbinsMC_antidantip; std::map> mixbinsMC_dp; + std::map> mixbinsMC_antipantip; + std::map> mixbinsMC_pp; + std::map> mixbinsMC_antipp; std::unique_ptr> Pair = std::make_unique>(); std::unique_ptr> PairMC = std::make_unique>(); @@ -376,6 +379,9 @@ struct hadronnucleicorrelation { registry.add("hReco_Pt_Proton", "Reco (anti)protons in reco collisions", {HistType::kTH1F, {pTAxis_small}}); registry.add("hReco_Pt_Deuteron", "Reco (anti)deuterons in reco collisions", {HistType::kTH1F, {pTAxis_small}}); + registry.add("hGen_EtaPhiPt_Proton", "Gen (anti)protons in gen collisions", {HistType::kTH3F, {etaAxis, phiAxis, pTAxis_small}}); + registry.add("hGen_EtaPhiPt_Deuteron", "Gen (anti)deuteron in gen collisions", {HistType::kTH3F, {etaAxis, phiAxis, pTAxis_small}}); + registry.add("hSec_EtaPhiPt_Proton", "Secondary (anti)protons", {HistType::kTH3F, {etaAxis, phiAxis, pTAxis_small}}); registry.add("hPrimSec_EtaPhiPt_Proton", "Primary + Secondary (anti)protons", {HistType::kTH3F, {etaAxis, phiAxis, pTAxis_small}}); @@ -1739,8 +1745,12 @@ struct hadronnucleicorrelation { void processGen(SimCollisions const& mcCollisions, SimParticles const& mcParticles) { + for (auto particle : mcParticles) { + if (std::abs(particle.template singleCollSel_as().posZ()) > cutzvertex) + continue; + if (particle.pdgCode() == pdgProton) { registry.fill(HIST("Generated/hQAProtons"), 0.5); } @@ -1776,15 +1786,19 @@ struct hadronnucleicorrelation { } if (particle.pdgCode() == pdgDeuteron) { + registry.fill(HIST("hGen_EtaPhiPt_Deuteron"), particle.eta(), particle.phi(), particle.pt()); selectedparticlesMC_d[particle.mcCollisionId()].push_back(std::make_shared(particle)); } if (particle.pdgCode() == -pdgDeuteron) { + registry.fill(HIST("hGen_EtaPhiPt_Deuteron"), particle.eta(), particle.phi(), -1. * particle.pt()); selectedparticlesMC_antid[particle.mcCollisionId()].push_back(std::make_shared(particle)); } if (particle.pdgCode() == pdgProton) { + registry.fill(HIST("hGen_EtaPhiPt_Proton"), particle.eta(), particle.phi(), particle.pt()); selectedparticlesMC_p[particle.mcCollisionId()].push_back(std::make_shared(particle)); } if (particle.pdgCode() == -pdgProton) { + registry.fill(HIST("hGen_EtaPhiPt_Proton"), particle.eta(), particle.phi(), -1. * particle.pt()); selectedparticlesMC_antip[particle.mcCollisionId()].push_back(std::make_shared(particle)); } } @@ -1793,6 +1807,10 @@ struct hadronnucleicorrelation { registry.fill(HIST("Generated/hNEventsMC"), 0.5); + if (std::abs(collision1.posZ()) > cutzvertex) { + continue; + } + // anti-d - anti-p correlation if (selectedparticlesMC_antid.find(collision1.globalIndex()) != selectedparticlesMC_antid.end()) { if (selectedparticlesMC_antip.find(collision1.globalIndex()) != selectedparticlesMC_antip.end()) { From 407a9a49a31371e717d356825cc05efeb78955fd Mon Sep 17 00:00:00 2001 From: Francesca Ercolessi Date: Tue, 7 Oct 2025 17:10:20 +0200 Subject: [PATCH 2/6] fix --- PWGLF/Tasks/Nuspex/hadronnucleicorrelation.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/PWGLF/Tasks/Nuspex/hadronnucleicorrelation.cxx b/PWGLF/Tasks/Nuspex/hadronnucleicorrelation.cxx index 424f61b8583..52dd7825f09 100644 --- a/PWGLF/Tasks/Nuspex/hadronnucleicorrelation.cxx +++ b/PWGLF/Tasks/Nuspex/hadronnucleicorrelation.cxx @@ -1748,7 +1748,7 @@ struct hadronnucleicorrelation { for (auto particle : mcParticles) { - if (std::abs(particle.template singleCollSel_as().posZ()) > cutzvertex) + if (std::abs(particle.template mcCollision_as().posZ()) > cutzvertex) continue; if (particle.pdgCode() == pdgProton) { From 1700fac2ea373a23a0fb4e80bfe347c01a29a1aa Mon Sep 17 00:00:00 2001 From: Francesca Ercolessi Date: Fri, 10 Oct 2025 13:25:28 +0200 Subject: [PATCH 3/6] add multiplicity histogram --- PWGLF/Tasks/Nuspex/hadronnucleicorrelation.cxx | 13 ++++++++----- 1 file changed, 8 insertions(+), 5 deletions(-) diff --git a/PWGLF/Tasks/Nuspex/hadronnucleicorrelation.cxx b/PWGLF/Tasks/Nuspex/hadronnucleicorrelation.cxx index 52dd7825f09..76fb5151b85 100644 --- a/PWGLF/Tasks/Nuspex/hadronnucleicorrelation.cxx +++ b/PWGLF/Tasks/Nuspex/hadronnucleicorrelation.cxx @@ -156,11 +156,10 @@ struct hadronnucleicorrelation { std::map, std::vector> mixbins_antip; std::map, std::vector> mixbins_p; std::map, std::vector> mixbinsPID_antidantip; - std::map> mixbinsMC_antidantip; - std::map> mixbinsMC_dp; - std::map> mixbinsMC_antipantip; - std::map> mixbinsMC_pp; - std::map> mixbinsMC_antipp; + std::map> mixbinsMC_antid; + std::map> mixbinsMC_d; + std::map> mixbinsMC_antip; + std::map> mixbinsMC_p; std::unique_ptr> Pair = std::make_unique>(); std::unique_ptr> PairMC = std::make_unique>(); @@ -323,6 +322,7 @@ struct hadronnucleicorrelation { if (doQA) { // Track QA + QA.add("QA/hMult", "multiplicity", {HistType::kTH1D, {{150, 0.f, 150.f, "N_{ch}"}}}); QA.add("QA/hVtxZ_trk", "#it{z}_{vtx}", {HistType::kTH1D, {{150, -15.f, 15.f, "#it{z}_{vtx} (cm)"}}}); QA.add("QA/hTPCnClusters", "N TPC Clusters; N TPC Clusters", {HistType::kTH1D, {{200, 0.f, 200.f}}}); QA.add("QA/hTPCSharedClusters", "N TPC Shared Clusters; N TPC SharedClusters", {HistType::kTH1D, {{100, 0.f, 1.f}}}); @@ -927,6 +927,7 @@ struct hadronnucleicorrelation { continue; registry.fill(HIST("hNEvents"), 0.5); + QA.fill(HIST("QA/hMult"), collision.mult()); if (selectedtracks_antid.find(collision.globalIndex()) != selectedtracks_antid.end() && selectedtracks_antip.find(collision.globalIndex()) != selectedtracks_antip.end()) { @@ -1619,6 +1620,7 @@ struct hadronnucleicorrelation { if (std::abs(collision.posZ()) > cutzvertex) continue; registry.fill(HIST("hNEvents"), 0.5); + QA.fill(HIST("QA/hMult"), collision.mult()); int vertexBinToMix = std::floor((collision.posZ() + cutzvertex) / (2 * cutzvertex / _vertexNbinsToMix)); int centBinToMix = std::floor(collision.multPerc() / (100.0 / _multNsubBins)); @@ -1806,6 +1808,7 @@ struct hadronnucleicorrelation { for (auto collision1 : mcCollisions) { // loop on collisions registry.fill(HIST("Generated/hNEventsMC"), 0.5); + QA.fill(HIST("QA/hMult"), collision1.mult()); if (std::abs(collision1.posZ()) > cutzvertex) { continue; From 3cd64255620a55e93b65864bb3a20e2f341b41ed Mon Sep 17 00:00:00 2001 From: Francesca Ercolessi Date: Fri, 10 Oct 2025 13:31:45 +0200 Subject: [PATCH 4/6] add mixing bins for gen analysis --- .../Tasks/Nuspex/hadronnucleicorrelation.cxx | 24 +++++++++++++++++++ 1 file changed, 24 insertions(+) diff --git a/PWGLF/Tasks/Nuspex/hadronnucleicorrelation.cxx b/PWGLF/Tasks/Nuspex/hadronnucleicorrelation.cxx index 76fb5151b85..5abea2d3a96 100644 --- a/PWGLF/Tasks/Nuspex/hadronnucleicorrelation.cxx +++ b/PWGLF/Tasks/Nuspex/hadronnucleicorrelation.cxx @@ -105,6 +105,8 @@ struct hadronnucleicorrelation { // Mixing parameters Configurable _vertexNbinsToMix{"vertexNbinsToMix", 10, "Number of vertexZ bins for the mixing"}; Configurable _multNsubBins{"multSubBins", 10, "number of sub-bins to perform the mixing within"}; + Configurable maxmultmix{"maxmultmix", 20, "maximum multiplicity to mix"}; + // pT/A bins Configurable> pTBins{"pTBins", {0.6f, 1.0f, 1.2f, 2.f}, "p_{T} bins"}; @@ -1814,6 +1816,28 @@ struct hadronnucleicorrelation { continue; } + int vertexBinToMix = std::floor((collision.posZ() + cutzvertex) / (2 * cutzvertex / _vertexNbinsToMix)); + int centBinToMix = std::floor(collision.mult() / (maxmultmix / _multNsubBins)); + + if (collision.mult()>maxmultmix) centBinToMix=_multNsubBins-1; // to avoid overflow in centrality bin + if (centBinToMix<0) centBinToMix=0; // to avoid underflow in centrality bin + + if (selectedparticlesMC_antid.find(collision1.globalIndex()) != selectedparticlesMC_antid.end()) { + mixbinsMC_antid[std::pair{vertexBinToMix, centBinToMix}].push_back(std::make_shared(collision1)); + } + + if (selectedparticlesMC_d.find(collision1.globalIndex()) != selectedparticlesMC_d.end()) { + mixbinsMC_d[std::pair{vertexBinToMix, centBinToMix}].push_back(std::make_shared(collision1)); + } + + if (selectedparticlesMC_antip.find(collision1.globalIndex()) != selectedparticlesMC_antip.end()) { + mixbinsMC_antip[std::pair{vertexBinToMix, centBinToMix}].push_back(std::make_shared(collision1)); + } + + if (selectedparticlesMC_p.find(collision1.globalIndex()) != selectedparticlesMC_p.end()) { + mixbinsMC_p[std::pair{vertexBinToMix, centBinToMix}].push_back(std::make_shared(collision1)); + } + // anti-d - anti-p correlation if (selectedparticlesMC_antid.find(collision1.globalIndex()) != selectedparticlesMC_antid.end()) { if (selectedparticlesMC_antip.find(collision1.globalIndex()) != selectedparticlesMC_antip.end()) { From 701bde0856ef64deb2b7e918b47dff1f24ee7d8c Mon Sep 17 00:00:00 2001 From: Francesca Ercolessi Date: Sat, 11 Oct 2025 15:04:08 +0200 Subject: [PATCH 5/6] Add correlations in processGen --- .../Tasks/Nuspex/hadronnucleicorrelation.cxx | 336 +++++++++++++----- 1 file changed, 253 insertions(+), 83 deletions(-) diff --git a/PWGLF/Tasks/Nuspex/hadronnucleicorrelation.cxx b/PWGLF/Tasks/Nuspex/hadronnucleicorrelation.cxx index 5abea2d3a96..c5a28e72a30 100644 --- a/PWGLF/Tasks/Nuspex/hadronnucleicorrelation.cxx +++ b/PWGLF/Tasks/Nuspex/hadronnucleicorrelation.cxx @@ -183,8 +183,16 @@ struct hadronnucleicorrelation { std::vector> hPIDEtaPhiGen_AntiDeAntiPr_ME; std::vector> hEtaPhiGen_AntiPrAntiPr_SE; std::vector> hEtaPhiGen_AntiPrAntiPr_ME; + std::vector> hEtaPhiGen_PrPr_SE; + std::vector> hEtaPhiGen_PrPr_ME; std::vector> hEtaPhiGen_AntiPrPr_SE; std::vector> hEtaPhiGen_AntiPrPr_ME; + std::vector> hEtaPhiGen_AntiDePr_SE; + std::vector> hEtaPhiGen_AntiDePr_ME; + std::vector> hEtaPhiGen_DeAntiPr_SE; + std::vector> hEtaPhiGen_DeAntiPr_ME; + std::vector> hEtaPhiGen_DePr_SE; + std::vector> hEtaPhiGen_DePr_ME; int nBinspT; TH2F* hEffpTEta_proton; @@ -257,6 +265,7 @@ struct hadronnucleicorrelation { if (isMCGen || mcCorrelation) { for (int i = 0; i < nBinspT; i++) { + // antid-antip auto htempSEGen_AntiDeAntiPr = registry.add(Form("hEtaPhiGen_AntiDeAntiPr_SE_pt%02.0f%02.0f", pTBins.value.at(i) * 10, pTBins.value.at(i + 1) * 10), Form("Gen #Delta#eta#Delta#phi (%.1f(Form("hEtaPhiGen_AntiPrAntiPr_SE_pt%02.0f%02.0f", pTBins.value.at(i) * 10, pTBins.value.at(i + 1) * 10), Form("Gen #Delta#eta#Delta#phi (%.1f(Form("hEtaPhiGen_PrPr_SE_pt%02.0f%02.0f", pTBins.value.at(i) * 10, pTBins.value.at(i + 1) * 10), + Form("Gen #Delta#eta#Delta#phi (%.1f(Form("hEtaPhiGen_PrPr_ME_pt%02.0f%02.0f", pTBins.value.at(i) * 10, pTBins.value.at(i + 1) * 10), + Form("Gen #Delta#eta#Delta#phi (%.1f(Form("hEtaPhiGen_AntiPrPr_SE_pt%02.0f%02.0f", pTBins.value.at(i) * 10, pTBins.value.at(i + 1) * 10), - Form("Gen #Delta#eta#Delta#phi (%.1f(Form("hEtaPhiGen_AntiPrPr_ME_pt%02.0f%02.0f", pTBins.value.at(i) * 10, pTBins.value.at(i + 1) * 10), - Form("Gen #Delta#eta#Delta#phi (%.1f(Form("hEtaPhiGen_DePr_SE_pt%02.0f%02.0f", pTBins.value.at(i) * 10, pTBins.value.at(i + 1) * 10), + Form("Gen #Delta#eta#Delta#phi (%.1f(Form("hEtaPhiGen_DePr_ME_pt%02.0f%02.0f", pTBins.value.at(i) * 10, pTBins.value.at(i + 1) * 10), + Form("Gen #Delta#eta#Delta#phi (%.1f(Form("hEtaPhiGen_AntiDePr_SE_pt%02.0f%02.0f", pTBins.value.at(i) * 10, pTBins.value.at(i + 1) * 10), + Form("Gen #Delta#eta#Delta#phi (%.1f(Form("hEtaPhiGen_AntiDePr_ME_pt%02.0f%02.0f", pTBins.value.at(i) * 10, pTBins.value.at(i + 1) * 10), + Form("Gen #Delta#eta#Delta#phi (%.1f(Form("hEtaPhiGen_DeAntiPr_SE_pt%02.0f%02.0f", pTBins.value.at(i) * 10, pTBins.value.at(i + 1) * 10), + Form("Gen #Delta#eta#Delta#phi (%.1f(Form("hEtaPhiGen_DeAntiPr_ME_pt%02.0f%02.0f", pTBins.value.at(i) * 10, pTBins.value.at(i + 1) * 10), + Form("Gen #Delta#eta#Delta#phi (%.1f - void mixMCParticles(Type const& particles1, Type const& particles2, bool ispap) + void mixMCParticles(Type const& particles1, Type const& particles2, bool mode) { for (auto const& it1 : particles1) { for (auto const& it2 : particles2) { @@ -743,15 +787,27 @@ struct hadronnucleicorrelation { if (it1->pt() >= pTBins.value.at(k) && it1->pt() < pTBins.value.at(k + 1)) { // Use correct histogram based on ME flag if constexpr (ME) { - if (ispap) + if (mode==0) hEtaPhiGen_AntiPrPr_ME[k]->Fill(deltaEtaGen, deltaPhiGen, it2->pt()); - else + else if (mode==1) hEtaPhiGen_AntiDeAntiPr_ME[k]->Fill(deltaEtaGen, deltaPhiGen, it2->pt()); + else if (mode==2) + hEtaPhiGen_AntiDePr_ME[k]->Fill(deltaEtaGen, deltaPhiGen, it2->pt()); + else if (mode==3) + hEtaPhiGen_DeAntiPr_ME[k]->Fill(deltaEtaGen, deltaPhiGen, it2->pt()); + else if (mode==4) + hEtaPhiGen_DePr_ME[k]->Fill(deltaEtaGen, deltaPhiGen, it2->pt()); } else { - if (ispap) + if (mode==0) hEtaPhiGen_AntiPrPr_SE[k]->Fill(deltaEtaGen, deltaPhiGen, it2->pt()); - else + else if (mode==1) hEtaPhiGen_AntiDeAntiPr_SE[k]->Fill(deltaEtaGen, deltaPhiGen, it2->pt()); + else if (mode==2) + hEtaPhiGen_AntiDePr_SE[k]->Fill(deltaEtaGen, deltaPhiGen, it2->pt()); + else if (mode==3) + hEtaPhiGen_DeAntiPr_SE[k]->Fill(deltaEtaGen, deltaPhiGen, it2->pt()); + else if (mode==4) + hEtaPhiGen_DePr_SE[k]->Fill(deltaEtaGen, deltaPhiGen, it2->pt()); } } } @@ -760,7 +816,7 @@ struct hadronnucleicorrelation { } template - void mixMCParticlesIdentical(Type const& particles1, Type const& particles2) + void mixMCParticlesIdentical(Type const& particles1, Type const& particles2, bool ismatter) { for (auto const& it1 : particles1) { for (auto const& it2 : particles2) { @@ -768,7 +824,7 @@ struct hadronnucleicorrelation { float deltaEtaGen = it2->eta() - it1->eta(); float deltaPhiGen = RecoDecay::constrainAngle(it2->phi() - it1->phi(), -1 * o2::constants::math::PIHalf); - if (!ME && std::abs(deltaPhiGen) < 0.001 && std::abs(deltaEtaGen) < 0.001) { + if (!ME && std::abs(deltaPhiGen) < 0.0001 && std::abs(deltaEtaGen) < 0.0001) { continue; } @@ -777,9 +833,15 @@ struct hadronnucleicorrelation { if (it1->pt() >= pTBins.value.at(k) && it1->pt() < pTBins.value.at(k + 1)) { // Use correct histogram based on ME flag if constexpr (ME) { - hEtaPhiGen_AntiPrAntiPr_ME[k]->Fill(deltaEtaGen, deltaPhiGen, it2->pt()); + if (ismatter) + hEtaPhiGen_PrPr_ME[k]->Fill(deltaEtaGen, deltaPhiGen, it2->pt()); + else + hEtaPhiGen_AntiPrAntiPr_ME[k]->Fill(deltaEtaGen, deltaPhiGen, it2->pt()); } else { - hEtaPhiGen_AntiPrAntiPr_SE[k]->Fill(deltaEtaGen, deltaPhiGen, it2->pt()); + if (ismatter) + hEtaPhiGen_PrPr_SE[k]->Fill(deltaEtaGen, deltaPhiGen, it2->pt()); + else + hEtaPhiGen_AntiPrAntiPr_SE[k]->Fill(deltaEtaGen, deltaPhiGen, it2->pt()); } } } @@ -1837,138 +1899,227 @@ struct hadronnucleicorrelation { if (selectedparticlesMC_p.find(collision1.globalIndex()) != selectedparticlesMC_p.end()) { mixbinsMC_p[std::pair{vertexBinToMix, centBinToMix}].push_back(std::make_shared(collision1)); } + } // coll - // anti-d - anti-p correlation - if (selectedparticlesMC_antid.find(collision1.globalIndex()) != selectedparticlesMC_antid.end()) { - if (selectedparticlesMC_antip.find(collision1.globalIndex()) != selectedparticlesMC_antip.end()) { - mixMCParticles<0>(selectedparticlesMC_antid[collision1.globalIndex()], selectedparticlesMC_antip[collision1.globalIndex()], 0); // mixing SE - } + if (!mixbinsMC_antip.empty()) { - int stop1 = 0; + //antip-antip + for (auto i = mixbinsMC_antip.begin(); i != mixbinsMC_antip.end(); i++) { // iterating over all vertex&mult bins - for (auto collision2 : mcCollisions) { // nested loop on collisions + std::vector value = i->second; + int EvPerBin = value.size(); // number of collisions in each vertex&mult bin - if (collision1.globalIndex() == collision2.globalIndex()) { - continue; - } + for (int indx1 = 0; indx1 < EvPerBin; indx1++) { // loop over all the events in each vertex&mult bin - if (stop1 > maxmixcollsGen) { - break; - } + auto col1 = value[indx1]; - if (selectedparticlesMC_antip.find(collision2.globalIndex()) != selectedparticlesMC_antip.end()) { - mixMCParticles<1>(selectedparticlesMC_antid[collision1.globalIndex()], selectedparticlesMC_antip[collision2.globalIndex()], 0); // mixing ME + if (selectedparticlesMC_antip.find(col1->index()) != selectedparticlesMC_antip.end()) { + mixMCParticlesIdentical<0>(selectedparticlesMC_antip[collision1.globalIndex()], selectedparticlesMC_antip[collision1.globalIndex()], 0); // mixing SE } - stop1++; + for (int indx2 = indx1 + 1; indx2 < EvPerBin; indx2++) { // nested loop for all the combinations of collisions in a chosen mult/vertex bin + + auto col2 = (i->second)[indx2]; + + if (col1 == col2) { + continue; + } + + if (selectedparticlesMC_antip.find(col2->index()) != selectedparticlesMC_antip.end()) { + mixMCParticlesIdentical<1>(selectedparticlesMC_antip[collision1.globalIndex()], selectedparticlesMC_antip[collision1.globalIndex()], 0); // mixing SE + } + } } } - // d - p correlation - if (domatterGen) { - if (selectedparticlesMC_d.find(collision1.globalIndex()) != selectedparticlesMC_d.end()) { - if (selectedparticlesMC_p.find(collision1.globalIndex()) != selectedparticlesMC_p.end()) { - mixMCParticles<0>(selectedparticlesMC_d[collision1.globalIndex()], selectedparticlesMC_p[collision1.globalIndex()], 0); // mixing SE + //antip-p + for (auto i = mixbinsMC_antip.begin(); i != mixbinsMC_antip.end(); i++) { // iterating over all vertex&mult bins + + std::vector value = i->second; + int EvPerBin = value.size(); // number of collisions in each vertex&mult bin + + for (int indx1 = 0; indx1 < EvPerBin; indx1++) { // loop over all the events in each vertex&mult bin + + auto col1 = value[indx1]; + + if (selectedparticlesMC_antip.find(col1->index()) != selectedparticlesMC_antip.end()) { + mixMCParticles<0>(selectedparticlesMC_antip[collision1.globalIndex()], selectedparticlesMC_p[collision1.globalIndex()], 0); // mixing SE } - int stop2 = 0; + for (int indx2 = indx1 + 1; indx2 < EvPerBin; indx2++) { // nested loop for all the combinations of collisions in a chosen mult/vertex bin - for (auto collision2 : mcCollisions) { // nested loop on collisions + auto col2 = (i->second)[indx2]; - if (collision1.globalIndex() == collision2.globalIndex()) { + if (col1 == col2) { continue; } - if (stop2 > maxmixcollsGen) { - break; + if (selectedparticlesMC_antip.find(col2->index()) != selectedparticlesMC_antip.end()) { + mixMCParticles<1>(selectedparticlesMC_antip[collision1.globalIndex()], selectedparticlesMC_p[collision1.globalIndex()], 0); // mixing SE } + } + } + } + + } // mixbinsMC_antip + + if (!mixbinsMC_p.empty()) { + + //p-p + for (auto i = mixbinsMC_p.begin(); i != mixbinsMC_p.end(); i++) { // iterating over all vertex&mult bins + + std::vector value = i->second; + int EvPerBin = value.size(); // number of collisions in each vertex&mult bin + + for (int indx1 = 0; indx1 < EvPerBin; indx1++) { // loop over all the events in each vertex&mult bin + + auto col1 = value[indx1]; + + if (selectedparticlesMC_p.find(col1->index()) != selectedparticlesMC_p.end()) { + mixMCParticlesIdentical<0>(selectedparticlesMC_p[collision1.globalIndex()], selectedparticlesMC_p[collision1.globalIndex()], 1); // mixing SE + } + + for (int indx2 = indx1 + 1; indx2 < EvPerBin; indx2++) { // nested loop for all the combinations of collisions in a chosen mult/vertex bin - if (selectedparticlesMC_p.find(collision2.globalIndex()) != selectedparticlesMC_p.end()) { - mixMCParticles<1>(selectedparticlesMC_d[collision1.globalIndex()], selectedparticlesMC_p[collision2.globalIndex()], 0); // mixing ME + auto col2 = (i->second)[indx2]; + + if (col1 == col2) { + continue; } - stop2++; + if (selectedparticlesMC_p.find(col2->index()) != selectedparticlesMC_p.end()) { + mixMCParticlesIdentical<1>(selectedparticlesMC_p[collision1.globalIndex()], selectedparticlesMC_p[collision1.globalIndex()], 1); // mixing SE + } } } } + } // mixbinsMC_p - // anti-p - anti-p correlation - if (selectedparticlesMC_antip.find(collision1.globalIndex()) != selectedparticlesMC_antip.end()) { + if (!mixbinsMC_antid.empty()) { - mixMCParticlesIdentical<0>(selectedparticlesMC_antip[collision1.globalIndex()], selectedparticlesMC_antip[collision1.globalIndex()]); // mixing SE + //antid-antip + for (auto i = mixbinsMC_antid.begin(); i != mixbinsMC_antid.end(); i++) { // iterating over all vertex&mult bins - int stop3 = 0; + std::vector value = i->second; + int EvPerBin = value.size(); // number of collisions in each vertex&mult bin - for (auto collision2 : mcCollisions) { // nested loop on collisions + for (int indx1 = 0; indx1 < EvPerBin; indx1++) { // loop over all the events in each vertex&mult bin - if (collision1.globalIndex() == collision2.globalIndex()) { - continue; - } + auto col1 = value[indx1]; - if (stop3 > maxmixcollsGen) { - break; + if (selectedparticlesMC_antid.find(col1->index()) != selectedparticlesMC_antid.end()) { + mixMCParticles<0>(selectedparticlesMC_antid[collision1.globalIndex()], selectedparticlesMC_antip[collision1.globalIndex()], 1); // mixing SE } - if (selectedparticlesMC_antip.find(collision2.globalIndex()) != selectedparticlesMC_antip.end()) { - mixMCParticlesIdentical<1>(selectedparticlesMC_antip[collision1.globalIndex()], selectedparticlesMC_antip[collision2.globalIndex()]); // mixing ME - } + for (int indx2 = indx1 + 1; indx2 < EvPerBin; indx2++) { // nested loop for all the combinations of collisions in a chosen mult/vertex bin - stop3++; + auto col2 = (i->second)[indx2]; + + if (col1 == col2) { + continue; + } + + if (selectedparticlesMC_antid.find(col2->index()) != selectedparticlesMC_antid.end()) { + mixMCParticles<1>(selectedparticlesMC_antid[collision1.globalIndex()], selectedparticlesMC_antip[collision1.globalIndex()], 1); // mixing SE + } + } } } - // p - p correlation - if (domatterGen) { - if (selectedparticlesMC_p.find(collision1.globalIndex()) != selectedparticlesMC_p.end()) { - mixMCParticlesIdentical<0>(selectedparticlesMC_p[collision1.globalIndex()], selectedparticlesMC_p[collision1.globalIndex()]); // mixing SE + //antid-p + for (auto i = mixbinsMC_antid.begin(); i != mixbinsMC_antid.end(); i++) { // iterating over all vertex&mult bins + + std::vector value = i->second; + int EvPerBin = value.size(); // number of collisions in each vertex&mult bin + + for (int indx1 = 0; indx1 < EvPerBin; indx1++) { // loop over all the events in each vertex&mult bin + + auto col1 = value[indx1]; - int stop4 = 0; + if (selectedparticlesMC_antid.find(col1->index()) != selectedparticlesMC_antid.end()) { + mixMCParticles<0>(selectedparticlesMC_antid[collision1.globalIndex()], selectedparticlesMC_p[collision1.globalIndex()], 2); // mixing SE + } - for (auto collision2 : mcCollisions) { // nested loop on collisions + for (int indx2 = indx1 + 1; indx2 < EvPerBin; indx2++) { // nested loop for all the combinations of collisions in a chosen mult/vertex bin - if (collision1.globalIndex() == collision2.globalIndex()) { + auto col2 = (i->second)[indx2]; + + if (col1 == col2) { continue; } - if (stop4 > maxmixcollsGen) { - break; + if (selectedparticlesMC_antid.find(col2->index()) != selectedparticlesMC_antid.end()) { + mixMCParticles<1>(selectedparticlesMC_antid[collision1.globalIndex()], selectedparticlesMC_p[collision1.globalIndex()], 2); // mixing SE } + } + } + } + + } // mixbinsMC_antid + + if (!mixbinsMC_d.empty()) { + + //d-antip + for (auto i = mixbinsMC_d.begin(); i != mixbinsMC_d.end(); i++) { // iterating over all vertex&mult bins + + std::vector value = i->second; + int EvPerBin = value.size(); // number of collisions in each vertex&mult bin + + for (int indx1 = 0; indx1 < EvPerBin; indx1++) { // loop over all the events in each vertex&mult bin + + auto col1 = value[indx1]; + + if (selectedparticlesMC_d.find(col1->index()) != selectedparticlesMC_d.end()) { + mixMCParticles<0>(selectedparticlesMC_d[collision1.globalIndex()], selectedparticlesMC_antip[collision1.globalIndex()], 3); // mixing SE + } + + for (int indx2 = indx1 + 1; indx2 < EvPerBin; indx2++) { // nested loop for all the combinations of collisions in a chosen mult/vertex bin + + auto col2 = (i->second)[indx2]; - if (selectedparticlesMC_p.find(collision2.globalIndex()) != selectedparticlesMC_p.end()) { - mixMCParticlesIdentical<1>(selectedparticlesMC_p[collision1.globalIndex()], selectedparticlesMC_p[collision2.globalIndex()]); // mixing ME + if (col1 == col2) { + continue; } - stop4++; + if (selectedparticlesMC_d.find(col2->index()) != selectedparticlesMC_d.end()) { + mixMCParticles<1>(selectedparticlesMC_d[collision1.globalIndex()], selectedparticlesMC_antip[collision1.globalIndex()], 3); // mixing SE + } } } } - // p-antip correlation - if (selectedparticlesMC_p.find(collision1.globalIndex()) != selectedparticlesMC_p.end()) { - if (selectedparticlesMC_antip.find(collision1.globalIndex()) != selectedparticlesMC_antip.end()) { - mixMCParticles<0>(selectedparticlesMC_p[collision1.globalIndex()], selectedparticlesMC_antip[collision1.globalIndex()], 1); // mixing SE - } + //d-p + for (auto i = mixbinsMC_d.begin(); i != mixbinsMC_d.end(); i++) { // iterating over all vertex&mult bins - int stop5 = 0; + std::vector value = i->second; + int EvPerBin = value.size(); // number of collisions in each vertex&mult bin - for (auto collision2 : mcCollisions) { // nested loop on collisions + for (int indx1 = 0; indx1 < EvPerBin; indx1++) { // loop over all the events in each vertex&mult bin - if (collision1.globalIndex() == collision2.globalIndex()) { - continue; - } + auto col1 = value[indx1]; - if (stop5 > maxmixcollsGen) { - break; + if (selectedparticlesMC_d.find(col1->index()) != selectedparticlesMC_d.end()) { + mixMCParticles<0>(selectedparticlesMC_d[collision1.globalIndex()], selectedparticlesMC_p[collision1.globalIndex()], 4); // mixing SE } - if (selectedparticlesMC_antip.find(collision2.globalIndex()) != selectedparticlesMC_antip.end()) { - mixMCParticles<1>(selectedparticlesMC_p[collision1.globalIndex()], selectedparticlesMC_antip[collision2.globalIndex()], 1); // mixing ME - } + for (int indx2 = indx1 + 1; indx2 < EvPerBin; indx2++) { // nested loop for all the combinations of collisions in a chosen mult/vertex bin + + auto col2 = (i->second)[indx2]; + + if (col1 == col2) { + continue; + } - stop5++; + if (selectedparticlesMC_d.find(col2->index()) != selectedparticlesMC_d.end()) { + mixMCParticles<1>(selectedparticlesMC_d[collision1.globalIndex()], selectedparticlesMC_p[collision1.globalIndex()], 4); // mixing SE + } + } } } - } + + } // mixbinsMC_d + + // clearing up for (auto i = selectedparticlesMC_antid.begin(); i != selectedparticlesMC_antid.end(); i++) @@ -1986,6 +2137,25 @@ struct hadronnucleicorrelation { for (auto i = selectedparticlesMC_p.begin(); i != selectedparticlesMC_p.end(); i++) (i->second).clear(); selectedparticlesMC_p.clear(); + + for (auto& pair : mixbinsMC_antip) { + pair.second.clear(); // clear the vector associated with the key + } + mixbinsMC_antip.clear(); // clear the map + + for (auto& pair : mixbinsMC_p) { + pair.second.clear(); // clear the vector associated with the key + } + mixbinsMC_p.clear(); // clear the map + for (auto& pair : mixbinsMC_antid) { + pair.second.clear(); // clear the vector associated with the key + } + mixbinsMC_antid.clear(); // clear the map + + for (auto& pair : mixbinsMC_d) { + pair.second.clear(); // clear the vector associated with the key + } + mixbinsMC_d.clear(); // clear the map } PROCESS_SWITCH(hadronnucleicorrelation, processGen, "processGen", false); }; From 5f398f74807b97e07d354a1e41e3d1e43833fd83 Mon Sep 17 00:00:00 2001 From: Francesca Ercolessi Date: Sat, 11 Oct 2025 18:02:05 +0200 Subject: [PATCH 6/6] remove unused histograms --- .../Tasks/Nuspex/hadronnucleicorrelation.cxx | 408 ++++-------------- 1 file changed, 87 insertions(+), 321 deletions(-) diff --git a/PWGLF/Tasks/Nuspex/hadronnucleicorrelation.cxx b/PWGLF/Tasks/Nuspex/hadronnucleicorrelation.cxx index c5a28e72a30..8c3ef349645 100644 --- a/PWGLF/Tasks/Nuspex/hadronnucleicorrelation.cxx +++ b/PWGLF/Tasks/Nuspex/hadronnucleicorrelation.cxx @@ -27,6 +27,7 @@ #include "Framework/DataTypes.h" #include "Framework/Expressions.h" #include "Framework/HistogramRegistry.h" +#include "Framework/O2DatabasePDGPlugin.h" #include "Framework/StaticFor.h" #include "Framework/runDataProcessing.h" #include "MathUtils/Utils.h" @@ -107,7 +108,6 @@ struct hadronnucleicorrelation { Configurable _multNsubBins{"multSubBins", 10, "number of sub-bins to perform the mixing within"}; Configurable maxmultmix{"maxmultmix", 20, "maximum multiplicity to mix"}; - // pT/A bins Configurable> pTBins{"pTBins", {0.6f, 1.0f, 1.2f, 2.f}, "p_{T} bins"}; @@ -141,27 +141,16 @@ struct hadronnucleicorrelation { std::map> selectedparticlesMC_antid; std::map> selectedparticlesMC_antip; - // key: int64_t - value: vector of trkType objects - std::map> selectedtracksMC_d; - std::map> selectedtracksMC_p; - std::map> selectedtracksMC_antid; - std::map> selectedtracksMC_antip; - std::map> selectedtracksPIDMC_d; - std::map> selectedtracksPIDMC_p; - std::map> selectedtracksPIDMC_antid; - std::map> selectedtracksPIDMC_antip; - // key: pair of an integer and a float - value: vector of colType objects // for each key I have a vector of collisions std::map, std::vector> mixbins_antid; std::map, std::vector> mixbins_d; std::map, std::vector> mixbins_antip; std::map, std::vector> mixbins_p; - std::map, std::vector> mixbinsPID_antidantip; - std::map> mixbinsMC_antid; - std::map> mixbinsMC_d; - std::map> mixbinsMC_antip; - std::map> mixbinsMC_p; + std::map, std::vector> mixbinsMC_antid; + std::map, std::vector> mixbinsMC_d; + std::map, std::vector> mixbinsMC_antip; + std::map, std::vector> mixbinsMC_p; std::unique_ptr> Pair = std::make_unique>(); std::unique_ptr> PairMC = std::make_unique>(); @@ -173,14 +162,8 @@ struct hadronnucleicorrelation { std::vector> hCorrEtaPhi_ME; // MC histograms - std::vector> hEtaPhiRec_AntiDeAntiPr_SE; std::vector> hEtaPhiGen_AntiDeAntiPr_SE; - std::vector> hEtaPhiRec_AntiDeAntiPr_ME; std::vector> hEtaPhiGen_AntiDeAntiPr_ME; - std::vector> hPIDEtaPhiRec_AntiDeAntiPr_SE; - std::vector> hPIDEtaPhiGen_AntiDeAntiPr_SE; - std::vector> hPIDEtaPhiRec_AntiDeAntiPr_ME; - std::vector> hPIDEtaPhiGen_AntiDeAntiPr_ME; std::vector> hEtaPhiGen_AntiPrAntiPr_SE; std::vector> hEtaPhiGen_AntiPrAntiPr_ME; std::vector> hEtaPhiGen_PrPr_SE; @@ -203,6 +186,8 @@ struct hadronnucleicorrelation { Service ccdb; o2::ccdb::CcdbApi ccdbApi; + Service pdgDB; + void init(o2::framework::InitContext&) { ccdb->setURL(url.value); @@ -239,31 +224,7 @@ struct hadronnucleicorrelation { nBinspT = pTBins.value.size() - 1; - if (mcCorrelation) { - for (int i = 0; i < nBinspT; i++) { - auto htempSERec_AntiDeAntiPr = registry.add(Form("hEtaPhiRec_AntiDeAntiPr_SE_pt%02.0f%02.0f", pTBins.value.at(i) * 10, pTBins.value.at(i + 1) * 10), - Form("Rec #Delta#eta#Delta#phi (%.1f(Form("hEtaPhiRec_AntiDeAntiPr_ME_pt%02.0f%02.0f", pTBins.value.at(i) * 10, pTBins.value.at(i + 1) * 10), - Form("Rec #Delta#eta#Delta#phi (%.1f(Form("hPIDEtaPhiRec_AntiDeAntiPr_SE_pt%02.0f%02.0f", pTBins.value.at(i) * 10, pTBins.value.at(i + 1) * 10), - Form("Rec #Delta#eta#Delta#phi (%.1f(Form("hPIDEtaPhiGen_AntiDeAntiPr_SE_pt%02.0f%02.0f", pTBins.value.at(i) * 10, pTBins.value.at(i + 1) * 10), - Form("Gen #Delta#eta#Delta#phi (%.1f(Form("hPIDEtaPhiRec_AntiDeAntiPr_ME_pt%02.0f%02.0f", pTBins.value.at(i) * 10, pTBins.value.at(i + 1) * 10), - Form("Rec #Delta#eta#Delta#phi (%.1f(Form("hPIDEtaPhiGen_AntiDeAntiPr_ME_pt%02.0f%02.0f", pTBins.value.at(i) * 10, pTBins.value.at(i + 1) * 10), - Form("Gen #Delta#eta#Delta#phi (%.1f(Form("hEtaPhiGen_AntiDeAntiPr_SE_pt%02.0f%02.0f", pTBins.value.at(i) * 10, pTBins.value.at(i + 1) * 10), @@ -283,10 +244,10 @@ struct hadronnucleicorrelation { // p-p auto htempSEGen_PrPr = registry.add(Form("hEtaPhiGen_PrPr_SE_pt%02.0f%02.0f", pTBins.value.at(i) * 10, pTBins.value.at(i + 1) * 10), - Form("Gen #Delta#eta#Delta#phi (%.1f(Form("hEtaPhiGen_PrPr_ME_pt%02.0f%02.0f", pTBins.value.at(i) * 10, pTBins.value.at(i + 1) * 10), - Form("Gen #Delta#eta#Delta#phi (%.1f(Form("hEtaPhiGen_DePr_SE_pt%02.0f%02.0f", pTBins.value.at(i) * 10, pTBins.value.at(i + 1) * 10), - Form("Gen #Delta#eta#Delta#phi (%.1f(Form("hEtaPhiGen_DePr_ME_pt%02.0f%02.0f", pTBins.value.at(i) * 10, pTBins.value.at(i + 1) * 10), - Form("Gen #Delta#eta#Delta#phi (%.1f(Form("hEtaPhiGen_DeAntiPr_ME_pt%02.0f%02.0f", pTBins.value.at(i) * 10, pTBins.value.at(i + 1) * 10), Form("Gen #Delta#eta#Delta#phi (%.1f(HIST("Generated/hNEventsMC"))->GetXaxis()->SetBinLabel(1, "All"); + registry.add("hGen_EtaPhiPt_Proton", "Gen (anti)protons in gen collisions", {HistType::kTH3F, {etaAxis, phiAxis, pTAxis_small}}); + registry.add("hGen_EtaPhiPt_Deuteron", "Gen (anti)deuteron in gen collisions", {HistType::kTH3F, {etaAxis, phiAxis, pTAxis_small}}); + registry.add("Generated/hQAProtons", "hQAProtons", {HistType::kTH1D, {{5, 0.f, 5.f}}}); registry.get(HIST("Generated/hQAProtons"))->GetXaxis()->SetBinLabel(1, "All"); registry.get(HIST("Generated/hQAProtons"))->GetXaxis()->SetBinLabel(2, "PhysicalPrimary"); @@ -717,64 +677,7 @@ struct hadronnucleicorrelation { } template - void mixTracksMC(Type const& tracks1, Type const& tracks2, bool isIdentical, bool isMCPID) - { // last value: 0 -- SE; 1 -- ME - for (auto const& it1 : tracks1) { - for (auto const& it2 : tracks2) { - - PairMC->SetPair(it1, it2); - PairMC->SetIdentical(isIdentical); - - // if Identical (pp and antip-antip) - if (isIdentical && PairMC->IsClosePair(deta, dphi, radiusTPC)) { - QA.fill(HIST("QA/hdetadphistar"), PairMC->GetPhiStarDiff(radiusTPC), PairMC->GetEtaDiff()); - continue; - } - - // Calculate Delta-eta Delta-phi (reco) - float deltaEta = it2->eta() - it1->eta(); - float deltaPhi = it2->phi() - it1->phi(); - deltaPhi = RecoDecay::constrainAngle(deltaPhi, -1 * o2::constants::math::PIHalf); - - // Calculate Delta-eta Delta-phi (gen) - float deltaEtaGen = -999.; - float deltaPhiGen = -999.; - deltaEtaGen = it2->eta_MC() - it1->eta_MC(); - deltaPhiGen = it2->phi_MC() - it1->phi_MC(); - deltaPhiGen = RecoDecay::constrainAngle(deltaPhiGen, -1 * o2::constants::math::PIHalf); - - for (int k = 0; k < nBinspT; k++) { - - if (it1->pt() >= pTBins.value.at(k) && it1->pt() < pTBins.value.at(k + 1)) { - - if (ME) { - if (isMCPID) { - hPIDEtaPhiRec_AntiDeAntiPr_ME[k]->Fill(deltaEta, deltaPhi, it2->pt()); - hPIDEtaPhiGen_AntiDeAntiPr_ME[k]->Fill(deltaEtaGen, deltaPhiGen, it2->pt()); - } else { - hEtaPhiRec_AntiDeAntiPr_ME[k]->Fill(deltaEta, deltaPhi, it2->pt()); - hEtaPhiGen_AntiDeAntiPr_ME[k]->Fill(deltaEtaGen, deltaPhiGen, it2->pt()); - } - } else { - if (isMCPID) { - hPIDEtaPhiRec_AntiDeAntiPr_SE[k]->Fill(deltaEta, deltaPhi, it2->pt()); - hPIDEtaPhiGen_AntiDeAntiPr_SE[k]->Fill(deltaEtaGen, deltaPhiGen, it2->pt()); - } else { - hEtaPhiRec_AntiDeAntiPr_SE[k]->Fill(deltaEta, deltaPhi, it2->pt()); - hEtaPhiGen_AntiDeAntiPr_SE[k]->Fill(deltaEtaGen, deltaPhiGen, it2->pt()); - } - } // SE - } // pT condition - } // nBinspT loop - - Pair->ResetPair(); - - } // tracks 2 - } // tracks 1 - } - - template - void mixMCParticles(Type const& particles1, Type const& particles2, bool mode) + void mixMCParticles(Type const& particles1, Type const& particles2, int mode) { for (auto const& it1 : particles1) { for (auto const& it2 : particles2) { @@ -787,26 +690,26 @@ struct hadronnucleicorrelation { if (it1->pt() >= pTBins.value.at(k) && it1->pt() < pTBins.value.at(k + 1)) { // Use correct histogram based on ME flag if constexpr (ME) { - if (mode==0) + if (mode == 0) hEtaPhiGen_AntiPrPr_ME[k]->Fill(deltaEtaGen, deltaPhiGen, it2->pt()); - else if (mode==1) + else if (mode == 1) hEtaPhiGen_AntiDeAntiPr_ME[k]->Fill(deltaEtaGen, deltaPhiGen, it2->pt()); - else if (mode==2) + else if (mode == 2) hEtaPhiGen_AntiDePr_ME[k]->Fill(deltaEtaGen, deltaPhiGen, it2->pt()); - else if (mode==3) + else if (mode == 3) hEtaPhiGen_DeAntiPr_ME[k]->Fill(deltaEtaGen, deltaPhiGen, it2->pt()); - else if (mode==4) + else if (mode == 4) hEtaPhiGen_DePr_ME[k]->Fill(deltaEtaGen, deltaPhiGen, it2->pt()); } else { - if (mode==0) + if (mode == 0) hEtaPhiGen_AntiPrPr_SE[k]->Fill(deltaEtaGen, deltaPhiGen, it2->pt()); - else if (mode==1) + else if (mode == 1) hEtaPhiGen_AntiDeAntiPr_SE[k]->Fill(deltaEtaGen, deltaPhiGen, it2->pt()); - else if (mode==2) + else if (mode == 2) hEtaPhiGen_AntiDePr_SE[k]->Fill(deltaEtaGen, deltaPhiGen, it2->pt()); - else if (mode==3) + else if (mode == 3) hEtaPhiGen_DeAntiPr_SE[k]->Fill(deltaEtaGen, deltaPhiGen, it2->pt()); - else if (mode==4) + else if (mode == 4) hEtaPhiGen_DePr_SE[k]->Fill(deltaEtaGen, deltaPhiGen, it2->pt()); } } @@ -991,7 +894,7 @@ struct hadronnucleicorrelation { continue; registry.fill(HIST("hNEvents"), 0.5); - QA.fill(HIST("QA/hMult"), collision.mult()); + registry.fill(HIST("hMult"), collision.mult()); if (selectedtracks_antid.find(collision.globalIndex()) != selectedtracks_antid.end() && selectedtracks_antip.find(collision.globalIndex()) != selectedtracks_antip.end()) { @@ -1283,7 +1186,7 @@ struct hadronnucleicorrelation { } PROCESS_SWITCH(hadronnucleicorrelation, processData, "processData", true); - void processMC(FilteredCollisions const& collisions, FilteredTracksMC const& tracks) + void processMC(FilteredCollisions const&, FilteredTracksMC const& tracks) { for (auto track : tracks) { if (std::abs(track.template singleCollSel_as().posZ()) > cutzvertex) @@ -1645,173 +1548,15 @@ struct hadronnucleicorrelation { track.sign() < 0) registry.fill(HIST("hDenominatorPurity_Deuteron_TPCEl_or_TOF"), track.pt() * -1); } - - if (!mcCorrelation) { - continue; - } - - if (isDe) { - selectedtracksPIDMC_d[track.singleCollSelId()].push_back(std::make_shared(track)); - } - if (track.pdgCode() == pdgDeuteron) { - selectedtracksMC_d[track.singleCollSelId()].push_back(std::make_shared(track)); - } - if (isAntiDe) { - selectedtracksPIDMC_antid[track.singleCollSelId()].push_back(std::make_shared(track)); - } - if (track.pdgCode() == -pdgDeuteron) { - selectedtracksMC_antid[track.singleCollSelId()].push_back(std::make_shared(track)); - } - if (isPr) { - selectedtracksPIDMC_p[track.singleCollSelId()].push_back(std::make_shared(track)); - } - if (track.pdgCode() == pdgProton) { - selectedtracksMC_p[track.singleCollSelId()].push_back(std::make_shared(track)); - } - if (isAntiPr) { - selectedtracksPIDMC_antip[track.singleCollSelId()].push_back(std::make_shared(track)); - } - if (track.pdgCode() == -pdgProton) { - selectedtracksMC_antip[track.singleCollSelId()].push_back(std::make_shared(track)); - } } // track - - if (!mcCorrelation) { - return; - } - - for (auto collision : collisions) { - if (std::abs(collision.posZ()) > cutzvertex) - continue; - registry.fill(HIST("hNEvents"), 0.5); - QA.fill(HIST("QA/hMult"), collision.mult()); - - int vertexBinToMix = std::floor((collision.posZ() + cutzvertex) / (2 * cutzvertex / _vertexNbinsToMix)); - int centBinToMix = std::floor(collision.multPerc() / (100.0 / _multNsubBins)); - - if (selectedtracksMC_antid.find(collision.globalIndex()) != selectedtracksMC_antid.end()) { - mixbins_antid[std::pair{vertexBinToMix, centBinToMix}].push_back(std::make_shared(collision)); - } - - if (selectedtracksPIDMC_antid.find(collision.globalIndex()) != selectedtracksPIDMC_antid.end()) { - mixbinsPID_antidantip[std::pair{vertexBinToMix, centBinToMix}].push_back(std::make_shared(collision)); - } - - PairMC->SetMagField1(collision.magField()); - PairMC->SetMagField2(collision.magField()); - } // coll - - if (!mixbins_antid.empty()) { - - for (auto i = mixbins_antid.begin(); i != mixbins_antid.end(); i++) { // iterating over all vertex&mult bins - - std::vector value = i->second; - int EvPerBin = value.size(); // number of collisions in each vertex&mult bin - - for (int indx1 = 0; indx1 < EvPerBin; indx1++) { // loop over all the events in each vertex&mult bin - - auto col1 = value[indx1]; - - if (selectedtracksMC_antip.find(col1->index()) != selectedtracksMC_antip.end()) { - mixTracksMC<0>(selectedtracksMC_antid[col1->index()], selectedtracksMC_antip[col1->index()], 0, 0); // mixing SE - } - - for (int indx2 = indx1 + 1; indx2 < EvPerBin; indx2++) { // nested loop for all the combinations of collisions in a chosen mult/vertex bin - - auto col2 = (i->second)[indx2]; - - if (col1 == col2) { - continue; - } - - if (selectedtracksMC_antip.find(col2->index()) != selectedtracksMC_antip.end()) { - mixTracksMC<1>(selectedtracksMC_antid[col1->index()], selectedtracksMC_antip[col2->index()], 0, 0); // mixing ME - } - } - } - } - } - - if (!mixbinsPID_antidantip.empty()) { - - for (auto i = mixbinsPID_antidantip.begin(); i != mixbinsPID_antidantip.end(); i++) { // iterating over all vertex&mult bins - - std::vector value = i->second; - int EvPerBin = value.size(); // number of collisions in each vertex&mult bin - - for (int indx1 = 0; indx1 < EvPerBin; indx1++) { // loop over all the events in each vertex&mult bin - - auto col1 = value[indx1]; - - if (selectedtracksPIDMC_antip.find(col1->index()) != selectedtracksPIDMC_antip.end()) { - mixTracksMC<0>(selectedtracksPIDMC_antid[col1->index()], selectedtracksPIDMC_antip[col1->index()], 0, 1); // mixing SE - } - - for (int indx2 = indx1 + 1; indx2 < EvPerBin; indx2++) { // nested loop for all the combinations of collisions in a chosen mult/vertex bin - - auto col2 = (i->second)[indx2]; - - if (col1 == col2) { - continue; - } - - if (selectedtracksPIDMC_antip.find(col2->index()) != selectedtracksPIDMC_antip.end()) { - mixTracksMC<1>(selectedtracksPIDMC_antid[col1->index()], selectedtracksPIDMC_antip[col2->index()], 0, 1); // mixing ME - } - } - } - } - } - - // clearing up - for (auto i = selectedtracksMC_antid.begin(); i != selectedtracksMC_antid.end(); i++) - (i->second).clear(); - selectedtracksMC_antid.clear(); - - for (auto i = selectedtracksMC_d.begin(); i != selectedtracksMC_d.end(); i++) - (i->second).clear(); - selectedtracksMC_d.clear(); - - for (auto i = selectedtracksMC_antip.begin(); i != selectedtracksMC_antip.end(); i++) - (i->second).clear(); - selectedtracksMC_antip.clear(); - - for (auto i = selectedtracksMC_p.begin(); i != selectedtracksMC_p.end(); i++) - (i->second).clear(); - selectedtracksMC_p.clear(); - - for (auto i = selectedtracksPIDMC_antid.begin(); i != selectedtracksPIDMC_antid.end(); i++) - (i->second).clear(); - selectedtracksPIDMC_antid.clear(); - - for (auto i = selectedtracksPIDMC_d.begin(); i != selectedtracksPIDMC_d.end(); i++) - (i->second).clear(); - selectedtracksPIDMC_d.clear(); - - for (auto i = selectedtracksPIDMC_antip.begin(); i != selectedtracksPIDMC_antip.end(); i++) - (i->second).clear(); - selectedtracksPIDMC_antip.clear(); - - for (auto i = selectedtracksPIDMC_p.begin(); i != selectedtracksPIDMC_p.end(); i++) - (i->second).clear(); - selectedtracksPIDMC_p.clear(); - - for (auto& pair : mixbinsPID_antidantip) { - pair.second.clear(); // clear the vector associated with the key - } - mixbinsPID_antidantip.clear(); // clear the map - - for (auto& pair : mixbins_antid) { - pair.second.clear(); // clear the vector associated with the key - } - mixbins_antid.clear(); // clear the map } PROCESS_SWITCH(hadronnucleicorrelation, processMC, "processMC", false); + Preslice perMCCol = aod::mcparticle::mcCollisionId; + void processGen(SimCollisions const& mcCollisions, SimParticles const& mcParticles) { - for (auto particle : mcParticles) { if (std::abs(particle.template mcCollision_as().posZ()) > cutzvertex) @@ -1872,17 +1617,39 @@ struct hadronnucleicorrelation { for (auto collision1 : mcCollisions) { // loop on collisions registry.fill(HIST("Generated/hNEventsMC"), 0.5); - QA.fill(HIST("QA/hMult"), collision1.mult()); if (std::abs(collision1.posZ()) > cutzvertex) { continue; } - int vertexBinToMix = std::floor((collision.posZ() + cutzvertex) / (2 * cutzvertex / _vertexNbinsToMix)); - int centBinToMix = std::floor(collision.mult() / (maxmultmix / _multNsubBins)); + const auto particlesInCollision = mcParticles.sliceBy(perMCCol, collision1.globalIndex()); + + float Ncharged = 0.; + for (auto& mcParticle : particlesInCollision) { + + if (!mcParticle.isPhysicalPrimary()) { + continue; + } + + if (std::abs(mcParticle.eta()) > 0.5f) { + continue; + } + + TParticlePDG* p = pdgDB->GetParticle(mcParticle.pdgCode()); + if (std::abs(p->Charge()) > 1E-3) { + Ncharged++; + } + } - if (collision.mult()>maxmultmix) centBinToMix=_multNsubBins-1; // to avoid overflow in centrality bin - if (centBinToMix<0) centBinToMix=0; // to avoid underflow in centrality bin + registry.fill(HIST("hMult"), Ncharged); + + int vertexBinToMix = std::floor((collision1.posZ() + cutzvertex) / (2 * cutzvertex / _vertexNbinsToMix)); + int centBinToMix = std::floor(Ncharged / (maxmultmix / _multNsubBins)); + + if (Ncharged > maxmultmix) + centBinToMix = _multNsubBins - 1; // to avoid overflow in centrality bin + if (centBinToMix < 0) + centBinToMix = 0; // to avoid underflow in centrality bin if (selectedparticlesMC_antid.find(collision1.globalIndex()) != selectedparticlesMC_antid.end()) { mixbinsMC_antid[std::pair{vertexBinToMix, centBinToMix}].push_back(std::make_shared(collision1)); @@ -1899,14 +1666,15 @@ struct hadronnucleicorrelation { if (selectedparticlesMC_p.find(collision1.globalIndex()) != selectedparticlesMC_p.end()) { mixbinsMC_p[std::pair{vertexBinToMix, centBinToMix}].push_back(std::make_shared(collision1)); } + } // coll if (!mixbinsMC_antip.empty()) { - //antip-antip + // antip-antip for (auto i = mixbinsMC_antip.begin(); i != mixbinsMC_antip.end(); i++) { // iterating over all vertex&mult bins - std::vector value = i->second; + std::vector value = i->second; int EvPerBin = value.size(); // number of collisions in each vertex&mult bin for (int indx1 = 0; indx1 < EvPerBin; indx1++) { // loop over all the events in each vertex&mult bin @@ -1914,7 +1682,7 @@ struct hadronnucleicorrelation { auto col1 = value[indx1]; if (selectedparticlesMC_antip.find(col1->index()) != selectedparticlesMC_antip.end()) { - mixMCParticlesIdentical<0>(selectedparticlesMC_antip[collision1.globalIndex()], selectedparticlesMC_antip[collision1.globalIndex()], 0); // mixing SE + mixMCParticlesIdentical<0>(selectedparticlesMC_antip[col1->index()], selectedparticlesMC_antip[col1->index()], 0); // mixing SE } for (int indx2 = indx1 + 1; indx2 < EvPerBin; indx2++) { // nested loop for all the combinations of collisions in a chosen mult/vertex bin @@ -1926,16 +1694,16 @@ struct hadronnucleicorrelation { } if (selectedparticlesMC_antip.find(col2->index()) != selectedparticlesMC_antip.end()) { - mixMCParticlesIdentical<1>(selectedparticlesMC_antip[collision1.globalIndex()], selectedparticlesMC_antip[collision1.globalIndex()], 0); // mixing SE + mixMCParticlesIdentical<1>(selectedparticlesMC_antip[col1->index()], selectedparticlesMC_antip[col2->index()], 0); // mixing SE } } } } - //antip-p + // antip-p for (auto i = mixbinsMC_antip.begin(); i != mixbinsMC_antip.end(); i++) { // iterating over all vertex&mult bins - std::vector value = i->second; + std::vector value = i->second; int EvPerBin = value.size(); // number of collisions in each vertex&mult bin for (int indx1 = 0; indx1 < EvPerBin; indx1++) { // loop over all the events in each vertex&mult bin @@ -1943,7 +1711,7 @@ struct hadronnucleicorrelation { auto col1 = value[indx1]; if (selectedparticlesMC_antip.find(col1->index()) != selectedparticlesMC_antip.end()) { - mixMCParticles<0>(selectedparticlesMC_antip[collision1.globalIndex()], selectedparticlesMC_p[collision1.globalIndex()], 0); // mixing SE + mixMCParticles<0>(selectedparticlesMC_antip[col1->index()], selectedparticlesMC_p[col1->index()], 0); // mixing SE } for (int indx2 = indx1 + 1; indx2 < EvPerBin; indx2++) { // nested loop for all the combinations of collisions in a chosen mult/vertex bin @@ -1955,7 +1723,7 @@ struct hadronnucleicorrelation { } if (selectedparticlesMC_antip.find(col2->index()) != selectedparticlesMC_antip.end()) { - mixMCParticles<1>(selectedparticlesMC_antip[collision1.globalIndex()], selectedparticlesMC_p[collision1.globalIndex()], 0); // mixing SE + mixMCParticles<1>(selectedparticlesMC_antip[col1->index()], selectedparticlesMC_p[col2->index()], 0); // mixing SE } } } @@ -1965,10 +1733,10 @@ struct hadronnucleicorrelation { if (!mixbinsMC_p.empty()) { - //p-p + // p-p for (auto i = mixbinsMC_p.begin(); i != mixbinsMC_p.end(); i++) { // iterating over all vertex&mult bins - std::vector value = i->second; + std::vector value = i->second; int EvPerBin = value.size(); // number of collisions in each vertex&mult bin for (int indx1 = 0; indx1 < EvPerBin; indx1++) { // loop over all the events in each vertex&mult bin @@ -1976,7 +1744,7 @@ struct hadronnucleicorrelation { auto col1 = value[indx1]; if (selectedparticlesMC_p.find(col1->index()) != selectedparticlesMC_p.end()) { - mixMCParticlesIdentical<0>(selectedparticlesMC_p[collision1.globalIndex()], selectedparticlesMC_p[collision1.globalIndex()], 1); // mixing SE + mixMCParticlesIdentical<0>(selectedparticlesMC_p[col1->index()], selectedparticlesMC_p[col1->index()], 1); // mixing SE } for (int indx2 = indx1 + 1; indx2 < EvPerBin; indx2++) { // nested loop for all the combinations of collisions in a chosen mult/vertex bin @@ -1988,7 +1756,7 @@ struct hadronnucleicorrelation { } if (selectedparticlesMC_p.find(col2->index()) != selectedparticlesMC_p.end()) { - mixMCParticlesIdentical<1>(selectedparticlesMC_p[collision1.globalIndex()], selectedparticlesMC_p[collision1.globalIndex()], 1); // mixing SE + mixMCParticlesIdentical<1>(selectedparticlesMC_p[col1->index()], selectedparticlesMC_p[col2->index()], 1); // mixing SE } } } @@ -1997,10 +1765,10 @@ struct hadronnucleicorrelation { if (!mixbinsMC_antid.empty()) { - //antid-antip + // antid-antip for (auto i = mixbinsMC_antid.begin(); i != mixbinsMC_antid.end(); i++) { // iterating over all vertex&mult bins - std::vector value = i->second; + std::vector value = i->second; int EvPerBin = value.size(); // number of collisions in each vertex&mult bin for (int indx1 = 0; indx1 < EvPerBin; indx1++) { // loop over all the events in each vertex&mult bin @@ -2008,7 +1776,7 @@ struct hadronnucleicorrelation { auto col1 = value[indx1]; if (selectedparticlesMC_antid.find(col1->index()) != selectedparticlesMC_antid.end()) { - mixMCParticles<0>(selectedparticlesMC_antid[collision1.globalIndex()], selectedparticlesMC_antip[collision1.globalIndex()], 1); // mixing SE + mixMCParticles<0>(selectedparticlesMC_antid[col1->index()], selectedparticlesMC_antip[col1->index()], 1); // mixing SE } for (int indx2 = indx1 + 1; indx2 < EvPerBin; indx2++) { // nested loop for all the combinations of collisions in a chosen mult/vertex bin @@ -2020,16 +1788,16 @@ struct hadronnucleicorrelation { } if (selectedparticlesMC_antid.find(col2->index()) != selectedparticlesMC_antid.end()) { - mixMCParticles<1>(selectedparticlesMC_antid[collision1.globalIndex()], selectedparticlesMC_antip[collision1.globalIndex()], 1); // mixing SE + mixMCParticles<1>(selectedparticlesMC_antid[col1->index()], selectedparticlesMC_antip[col2->index()], 1); // mixing SE } } } } - //antid-p + // antid-p for (auto i = mixbinsMC_antid.begin(); i != mixbinsMC_antid.end(); i++) { // iterating over all vertex&mult bins - std::vector value = i->second; + std::vector value = i->second; int EvPerBin = value.size(); // number of collisions in each vertex&mult bin for (int indx1 = 0; indx1 < EvPerBin; indx1++) { // loop over all the events in each vertex&mult bin @@ -2037,7 +1805,7 @@ struct hadronnucleicorrelation { auto col1 = value[indx1]; if (selectedparticlesMC_antid.find(col1->index()) != selectedparticlesMC_antid.end()) { - mixMCParticles<0>(selectedparticlesMC_antid[collision1.globalIndex()], selectedparticlesMC_p[collision1.globalIndex()], 2); // mixing SE + mixMCParticles<0>(selectedparticlesMC_antid[col1->index()], selectedparticlesMC_p[col1->index()], 2); // mixing SE } for (int indx2 = indx1 + 1; indx2 < EvPerBin; indx2++) { // nested loop for all the combinations of collisions in a chosen mult/vertex bin @@ -2049,7 +1817,7 @@ struct hadronnucleicorrelation { } if (selectedparticlesMC_antid.find(col2->index()) != selectedparticlesMC_antid.end()) { - mixMCParticles<1>(selectedparticlesMC_antid[collision1.globalIndex()], selectedparticlesMC_p[collision1.globalIndex()], 2); // mixing SE + mixMCParticles<1>(selectedparticlesMC_antid[col1->index()], selectedparticlesMC_p[col2->index()], 2); // mixing SE } } } @@ -2059,10 +1827,10 @@ struct hadronnucleicorrelation { if (!mixbinsMC_d.empty()) { - //d-antip + // d-antip for (auto i = mixbinsMC_d.begin(); i != mixbinsMC_d.end(); i++) { // iterating over all vertex&mult bins - std::vector value = i->second; + std::vector value = i->second; int EvPerBin = value.size(); // number of collisions in each vertex&mult bin for (int indx1 = 0; indx1 < EvPerBin; indx1++) { // loop over all the events in each vertex&mult bin @@ -2070,7 +1838,7 @@ struct hadronnucleicorrelation { auto col1 = value[indx1]; if (selectedparticlesMC_d.find(col1->index()) != selectedparticlesMC_d.end()) { - mixMCParticles<0>(selectedparticlesMC_d[collision1.globalIndex()], selectedparticlesMC_antip[collision1.globalIndex()], 3); // mixing SE + mixMCParticles<0>(selectedparticlesMC_d[col1->index()], selectedparticlesMC_antip[col1->index()], 3); // mixing SE } for (int indx2 = indx1 + 1; indx2 < EvPerBin; indx2++) { // nested loop for all the combinations of collisions in a chosen mult/vertex bin @@ -2082,16 +1850,16 @@ struct hadronnucleicorrelation { } if (selectedparticlesMC_d.find(col2->index()) != selectedparticlesMC_d.end()) { - mixMCParticles<1>(selectedparticlesMC_d[collision1.globalIndex()], selectedparticlesMC_antip[collision1.globalIndex()], 3); // mixing SE + mixMCParticles<1>(selectedparticlesMC_d[col1->index()], selectedparticlesMC_antip[col2->index()], 3); // mixing SE } } } } - //d-p + // d-p for (auto i = mixbinsMC_d.begin(); i != mixbinsMC_d.end(); i++) { // iterating over all vertex&mult bins - std::vector value = i->second; + std::vector value = i->second; int EvPerBin = value.size(); // number of collisions in each vertex&mult bin for (int indx1 = 0; indx1 < EvPerBin; indx1++) { // loop over all the events in each vertex&mult bin @@ -2099,7 +1867,7 @@ struct hadronnucleicorrelation { auto col1 = value[indx1]; if (selectedparticlesMC_d.find(col1->index()) != selectedparticlesMC_d.end()) { - mixMCParticles<0>(selectedparticlesMC_d[collision1.globalIndex()], selectedparticlesMC_p[collision1.globalIndex()], 4); // mixing SE + mixMCParticles<0>(selectedparticlesMC_d[col1->index()], selectedparticlesMC_p[col1->index()], 4); // mixing SE } for (int indx2 = indx1 + 1; indx2 < EvPerBin; indx2++) { // nested loop for all the combinations of collisions in a chosen mult/vertex bin @@ -2111,7 +1879,7 @@ struct hadronnucleicorrelation { } if (selectedparticlesMC_d.find(col2->index()) != selectedparticlesMC_d.end()) { - mixMCParticles<1>(selectedparticlesMC_d[collision1.globalIndex()], selectedparticlesMC_p[collision1.globalIndex()], 4); // mixing SE + mixMCParticles<1>(selectedparticlesMC_d[col1->index()], selectedparticlesMC_p[col2->index()], 4); // mixing SE } } } @@ -2119,8 +1887,6 @@ struct hadronnucleicorrelation { } // mixbinsMC_d - - // clearing up for (auto i = selectedparticlesMC_antid.begin(); i != selectedparticlesMC_antid.end(); i++) (i->second).clear();