diff --git a/PWGLF/Tasks/Nuspex/hadronnucleicorrelation.cxx b/PWGLF/Tasks/Nuspex/hadronnucleicorrelation.cxx index ee383900901..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" @@ -105,6 +106,7 @@ 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"}; @@ -139,25 +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_antidantip; - std::map> mixbinsMC_dp; + 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>(); @@ -169,18 +162,20 @@ 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; + 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; @@ -191,6 +186,8 @@ struct hadronnucleicorrelation { Service ccdb; o2::ccdb::CcdbApi ccdbApi; + Service pdgDB; + void init(o2::framework::InitContext&) { ccdb->setURL(url.value); @@ -227,32 +224,9 @@ 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), 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(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"); @@ -665,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 ispap) + void mixMCParticles(Type const& particles1, Type const& particles2, int mode) { for (auto const& it1 : particles1) { for (auto const& it2 : particles2) { @@ -735,15 +690,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()); } } } @@ -752,7 +719,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) { @@ -760,7 +727,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; } @@ -769,9 +736,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()); } } } @@ -921,6 +894,7 @@ struct hadronnucleicorrelation { continue; registry.fill(HIST("hNEvents"), 0.5); + registry.fill(HIST("hMult"), collision.mult()); if (selectedtracks_antid.find(collision.globalIndex()) != selectedtracks_antid.end() && selectedtracks_antip.find(collision.globalIndex()) != selectedtracks_antip.end()) { @@ -1212,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) @@ -1574,74 +1548,141 @@ struct hadronnucleicorrelation { track.sign() < 0) registry.fill(HIST("hDenominatorPurity_Deuteron_TPCEl_or_TOF"), track.pt() * -1); } + } // track + } + PROCESS_SWITCH(hadronnucleicorrelation, processMC, "processMC", false); + + Preslice perMCCol = aod::mcparticle::mcCollisionId; - if (!mcCorrelation) { + void processGen(SimCollisions const& mcCollisions, + SimParticles const& mcParticles) + { + for (auto particle : mcParticles) { + + if (std::abs(particle.template mcCollision_as().posZ()) > cutzvertex) continue; + + if (particle.pdgCode() == pdgProton) { + registry.fill(HIST("Generated/hQAProtons"), 0.5); + } + if (particle.pdgCode() == pdgDeuteron) { + registry.fill(HIST("Generated/hQADeuterons"), 0.5); } - if (isDe) { - selectedtracksPIDMC_d[track.singleCollSelId()].push_back(std::make_shared(track)); + if (isPrim && !particle.isPhysicalPrimary()) { + continue; } - if (track.pdgCode() == pdgDeuteron) { - selectedtracksMC_d[track.singleCollSelId()].push_back(std::make_shared(track)); + if (particle.pdgCode() == pdgProton) { + registry.fill(HIST("Generated/hQAProtons"), 1.5); } - if (isAntiDe) { - selectedtracksPIDMC_antid[track.singleCollSelId()].push_back(std::make_shared(track)); + if (particle.pdgCode() == pdgDeuteron) { + registry.fill(HIST("Generated/hQADeuterons"), 1.5); } - if (track.pdgCode() == -pdgDeuteron) { - selectedtracksMC_antid[track.singleCollSelId()].push_back(std::make_shared(track)); + + if (particle.pdgCode() == pdgDeuteron && std::abs(particle.y()) < 0.5) { + registry.fill(HIST("Generated/hDeuteronsVsPt"), particle.pt()); } - if (isPr) { - selectedtracksPIDMC_p[track.singleCollSelId()].push_back(std::make_shared(track)); + if (particle.pdgCode() == -pdgDeuteron && std::abs(particle.y()) < 0.5) { + registry.fill(HIST("Generated/hAntiDeuteronsVsPt"), particle.pt()); } - if (track.pdgCode() == pdgProton) { - selectedtracksMC_p[track.singleCollSelId()].push_back(std::make_shared(track)); + + if (std::abs(particle.eta()) > etacut) { + continue; } - if (isAntiPr) { - selectedtracksPIDMC_antip[track.singleCollSelId()].push_back(std::make_shared(track)); + if (particle.pdgCode() == pdgProton) { + registry.fill(HIST("Generated/hQAProtons"), 2.5); } - if (track.pdgCode() == -pdgProton) { - selectedtracksMC_antip[track.singleCollSelId()].push_back(std::make_shared(track)); + if (particle.pdgCode() == pdgDeuteron) { + registry.fill(HIST("Generated/hQADeuterons"), 2.5); } - } // track - if (!mcCorrelation) { - return; + 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)); + } } - for (auto collision : collisions) { - if (std::abs(collision.posZ()) > cutzvertex) + for (auto collision1 : mcCollisions) { // loop on collisions + + registry.fill(HIST("Generated/hNEventsMC"), 0.5); + + if (std::abs(collision1.posZ()) > cutzvertex) { continue; - registry.fill(HIST("hNEvents"), 0.5); + } - int vertexBinToMix = std::floor((collision.posZ() + cutzvertex) / (2 * cutzvertex / _vertexNbinsToMix)); - int centBinToMix = std::floor(collision.multPerc() / (100.0 / _multNsubBins)); + const auto particlesInCollision = mcParticles.sliceBy(perMCCol, collision1.globalIndex()); - if (selectedtracksMC_antid.find(collision.globalIndex()) != selectedtracksMC_antid.end()) { - mixbins_antid[std::pair{vertexBinToMix, centBinToMix}].push_back(std::make_shared(collision)); + 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++; + } + } + + 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)); + } + + if (selectedparticlesMC_d.find(collision1.globalIndex()) != selectedparticlesMC_d.end()) { + mixbinsMC_d[std::pair{vertexBinToMix, centBinToMix}].push_back(std::make_shared(collision1)); } - if (selectedtracksPIDMC_antid.find(collision.globalIndex()) != selectedtracksPIDMC_antid.end()) { - mixbinsPID_antidantip[std::pair{vertexBinToMix, centBinToMix}].push_back(std::make_shared(collision)); + 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)); } - PairMC->SetMagField1(collision.magField()); - PairMC->SetMagField2(collision.magField()); } // coll - if (!mixbins_antid.empty()) { + if (!mixbinsMC_antip.empty()) { - for (auto i = mixbins_antid.begin(); i != mixbins_antid.end(); i++) { // iterating over all vertex&mult bins + // 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 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 + if (selectedparticlesMC_antip.find(col1->index()) != selectedparticlesMC_antip.end()) { + 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 @@ -1652,27 +1693,25 @@ struct hadronnucleicorrelation { 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 (selectedparticlesMC_antip.find(col2->index()) != selectedparticlesMC_antip.end()) { + mixMCParticlesIdentical<1>(selectedparticlesMC_antip[col1->index()], selectedparticlesMC_antip[col2->index()], 0); // mixing SE } } } } - } - if (!mixbinsPID_antidantip.empty()) { + // antip-p + for (auto i = mixbinsMC_antip.begin(); i != mixbinsMC_antip.end(); i++) { // iterating over all vertex&mult bins - for (auto i = mixbinsPID_antidantip.begin(); i != mixbinsPID_antidantip.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 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 + if (selectedparticlesMC_antip.find(col1->index()) != selectedparticlesMC_antip.end()) { + 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 @@ -1683,247 +1722,170 @@ struct hadronnucleicorrelation { continue; } - if (selectedtracksPIDMC_antip.find(col2->index()) != selectedtracksPIDMC_antip.end()) { - mixTracksMC<1>(selectedtracksPIDMC_antid[col1->index()], selectedtracksPIDMC_antip[col2->index()], 0, 1); // mixing ME + if (selectedparticlesMC_antip.find(col2->index()) != selectedparticlesMC_antip.end()) { + mixMCParticles<1>(selectedparticlesMC_antip[col1->index()], selectedparticlesMC_p[col2->index()], 0); // mixing SE } } } } - } - - // 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(); + } // mixbinsMC_antip - for (auto i = selectedtracksMC_antip.begin(); i != selectedtracksMC_antip.end(); i++) - (i->second).clear(); - selectedtracksMC_antip.clear(); + if (!mixbinsMC_p.empty()) { - for (auto i = selectedtracksMC_p.begin(); i != selectedtracksMC_p.end(); i++) - (i->second).clear(); - selectedtracksMC_p.clear(); + // p-p + for (auto i = mixbinsMC_p.begin(); i != mixbinsMC_p.end(); i++) { // iterating over all vertex&mult bins - for (auto i = selectedtracksPIDMC_antid.begin(); i != selectedtracksPIDMC_antid.end(); i++) - (i->second).clear(); - selectedtracksPIDMC_antid.clear(); + std::vector value = i->second; + int EvPerBin = value.size(); // number of collisions in each vertex&mult bin - for (auto i = selectedtracksPIDMC_d.begin(); i != selectedtracksPIDMC_d.end(); i++) - (i->second).clear(); - selectedtracksPIDMC_d.clear(); + for (int indx1 = 0; indx1 < EvPerBin; indx1++) { // loop over all the events in each vertex&mult bin - for (auto i = selectedtracksPIDMC_antip.begin(); i != selectedtracksPIDMC_antip.end(); i++) - (i->second).clear(); - selectedtracksPIDMC_antip.clear(); + auto col1 = value[indx1]; - for (auto i = selectedtracksPIDMC_p.begin(); i != selectedtracksPIDMC_p.end(); i++) - (i->second).clear(); - selectedtracksPIDMC_p.clear(); + if (selectedparticlesMC_p.find(col1->index()) != selectedparticlesMC_p.end()) { + mixMCParticlesIdentical<0>(selectedparticlesMC_p[col1->index()], selectedparticlesMC_p[col1->index()], 1); // mixing SE + } - for (auto& pair : mixbinsPID_antidantip) { - pair.second.clear(); // clear the vector associated with the key - } - mixbinsPID_antidantip.clear(); // clear the map + for (int indx2 = indx1 + 1; indx2 < EvPerBin; indx2++) { // nested loop for all the combinations of collisions in a chosen mult/vertex bin - 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); + auto col2 = (i->second)[indx2]; - void processGen(SimCollisions const& mcCollisions, - SimParticles const& mcParticles) - { - for (auto particle : mcParticles) { + if (col1 == col2) { + continue; + } - if (particle.pdgCode() == pdgProton) { - registry.fill(HIST("Generated/hQAProtons"), 0.5); - } - if (particle.pdgCode() == pdgDeuteron) { - registry.fill(HIST("Generated/hQADeuterons"), 0.5); + if (selectedparticlesMC_p.find(col2->index()) != selectedparticlesMC_p.end()) { + mixMCParticlesIdentical<1>(selectedparticlesMC_p[col1->index()], selectedparticlesMC_p[col2->index()], 1); // mixing SE + } + } + } } + } // mixbinsMC_p - if (isPrim && !particle.isPhysicalPrimary()) { - continue; - } - if (particle.pdgCode() == pdgProton) { - registry.fill(HIST("Generated/hQAProtons"), 1.5); - } - if (particle.pdgCode() == pdgDeuteron) { - registry.fill(HIST("Generated/hQADeuterons"), 1.5); - } + if (!mixbinsMC_antid.empty()) { - if (particle.pdgCode() == pdgDeuteron && std::abs(particle.y()) < 0.5) { - registry.fill(HIST("Generated/hDeuteronsVsPt"), particle.pt()); - } - if (particle.pdgCode() == -pdgDeuteron && std::abs(particle.y()) < 0.5) { - registry.fill(HIST("Generated/hAntiDeuteronsVsPt"), particle.pt()); - } + // antid-antip + for (auto i = mixbinsMC_antid.begin(); i != mixbinsMC_antid.end(); i++) { // iterating over all vertex&mult bins - if (std::abs(particle.eta()) > etacut) { - continue; - } - if (particle.pdgCode() == pdgProton) { - registry.fill(HIST("Generated/hQAProtons"), 2.5); - } - if (particle.pdgCode() == pdgDeuteron) { - registry.fill(HIST("Generated/hQADeuterons"), 2.5); - } + std::vector value = i->second; + int EvPerBin = value.size(); // number of collisions in each vertex&mult bin - if (particle.pdgCode() == pdgDeuteron) { - selectedparticlesMC_d[particle.mcCollisionId()].push_back(std::make_shared(particle)); - } - if (particle.pdgCode() == -pdgDeuteron) { - selectedparticlesMC_antid[particle.mcCollisionId()].push_back(std::make_shared(particle)); - } - if (particle.pdgCode() == pdgProton) { - selectedparticlesMC_p[particle.mcCollisionId()].push_back(std::make_shared(particle)); - } - if (particle.pdgCode() == -pdgProton) { - selectedparticlesMC_antip[particle.mcCollisionId()].push_back(std::make_shared(particle)); - } - } + for (int indx1 = 0; indx1 < EvPerBin; indx1++) { // loop over all the events in each vertex&mult bin - for (auto collision1 : mcCollisions) { // loop on collisions + auto col1 = value[indx1]; - registry.fill(HIST("Generated/hNEventsMC"), 0.5); + if (selectedparticlesMC_antid.find(col1->index()) != selectedparticlesMC_antid.end()) { + mixMCParticles<0>(selectedparticlesMC_antid[col1->index()], selectedparticlesMC_antip[col1->index()], 1); // mixing SE + } - // 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 - } + for (int indx2 = indx1 + 1; indx2 < EvPerBin; indx2++) { // nested loop for all the combinations of collisions in a chosen mult/vertex bin - int stop1 = 0; + auto col2 = (i->second)[indx2]; - for (auto collision2 : mcCollisions) { // nested loop on collisions + if (col1 == col2) { + continue; + } - if (collision1.globalIndex() == collision2.globalIndex()) { - continue; + if (selectedparticlesMC_antid.find(col2->index()) != selectedparticlesMC_antid.end()) { + mixMCParticles<1>(selectedparticlesMC_antid[col1->index()], selectedparticlesMC_antip[col2->index()], 1); // mixing SE + } } + } + } - if (stop1 > maxmixcollsGen) { - break; - } + // antid-p + for (auto i = mixbinsMC_antid.begin(); i != mixbinsMC_antid.end(); i++) { // iterating over all vertex&mult bins - if (selectedparticlesMC_antip.find(collision2.globalIndex()) != selectedparticlesMC_antip.end()) { - mixMCParticles<1>(selectedparticlesMC_antid[collision1.globalIndex()], selectedparticlesMC_antip[collision2.globalIndex()], 0); // mixing ME - } + std::vector value = i->second; + int EvPerBin = value.size(); // number of collisions in each vertex&mult bin - stop1++; - } - } + for (int indx1 = 0; indx1 < EvPerBin; indx1++) { // loop over all the events in each vertex&mult bin + + auto col1 = value[indx1]; - // 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 + if (selectedparticlesMC_antid.find(col1->index()) != selectedparticlesMC_antid.end()) { + mixMCParticles<0>(selectedparticlesMC_antid[col1->index()], selectedparticlesMC_p[col1->index()], 2); // 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_p.find(collision2.globalIndex()) != selectedparticlesMC_p.end()) { - mixMCParticles<1>(selectedparticlesMC_d[collision1.globalIndex()], selectedparticlesMC_p[collision2.globalIndex()], 0); // mixing ME + if (selectedparticlesMC_antid.find(col2->index()) != selectedparticlesMC_antid.end()) { + mixMCParticles<1>(selectedparticlesMC_antid[col1->index()], selectedparticlesMC_p[col2->index()], 2); // mixing SE } - - stop2++; } } } - // anti-p - anti-p correlation - if (selectedparticlesMC_antip.find(collision1.globalIndex()) != selectedparticlesMC_antip.end()) { + } // mixbinsMC_antid - mixMCParticlesIdentical<0>(selectedparticlesMC_antip[collision1.globalIndex()], selectedparticlesMC_antip[collision1.globalIndex()]); // mixing SE + if (!mixbinsMC_d.empty()) { - int stop3 = 0; + // d-antip + for (auto i = mixbinsMC_d.begin(); i != mixbinsMC_d.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 (stop3 > maxmixcollsGen) { - break; - } + auto col1 = value[indx1]; - if (selectedparticlesMC_antip.find(collision2.globalIndex()) != selectedparticlesMC_antip.end()) { - mixMCParticlesIdentical<1>(selectedparticlesMC_antip[collision1.globalIndex()], selectedparticlesMC_antip[collision2.globalIndex()]); // mixing ME + if (selectedparticlesMC_d.find(col1->index()) != selectedparticlesMC_d.end()) { + mixMCParticles<0>(selectedparticlesMC_d[col1->index()], selectedparticlesMC_antip[col1->index()], 3); // mixing SE } - stop3++; - } - } - // 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 - - int stop4 = 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 (stop4 > maxmixcollsGen) { - break; + if (selectedparticlesMC_d.find(col2->index()) != selectedparticlesMC_d.end()) { + mixMCParticles<1>(selectedparticlesMC_d[col1->index()], selectedparticlesMC_antip[col2->index()], 3); // mixing SE } - - if (selectedparticlesMC_p.find(collision2.globalIndex()) != selectedparticlesMC_p.end()) { - mixMCParticlesIdentical<1>(selectedparticlesMC_p[collision1.globalIndex()], selectedparticlesMC_p[collision2.globalIndex()]); // mixing ME - } - - stop4++; } } } - // 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[col1->index()], selectedparticlesMC_p[col1->index()], 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[col1->index()], selectedparticlesMC_p[col2->index()], 4); // mixing SE + } + } } } - } + + } // mixbinsMC_d // clearing up for (auto i = selectedparticlesMC_antid.begin(); i != selectedparticlesMC_antid.end(); i++) @@ -1941,6 +1903,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); };