From 542ed5ba88ee8cfd6f2d3eb39f49f830db078b7f Mon Sep 17 00:00:00 2001 From: Francesca Ercolessi Date: Thu, 30 Jan 2025 15:17:21 +0100 Subject: [PATCH] [PWGLF] nuclei-proton correlations: Add process Gen for simulation studies (#9622) --- .../Tasks/Nuspex/hadronnucleicorrelation.cxx | 207 ++++++++++++++++-- 1 file changed, 194 insertions(+), 13 deletions(-) diff --git a/PWGLF/Tasks/Nuspex/hadronnucleicorrelation.cxx b/PWGLF/Tasks/Nuspex/hadronnucleicorrelation.cxx index e3e426fcad1..8d0e9cbe231 100644 --- a/PWGLF/Tasks/Nuspex/hadronnucleicorrelation.cxx +++ b/PWGLF/Tasks/Nuspex/hadronnucleicorrelation.cxx @@ -59,6 +59,8 @@ struct hadronnucleicorrelation { Configurable doQA{"doQA", true, "save QA histograms"}; Configurable doMCQA{"doMCQA", false, "save MC QA histograms"}; Configurable isMC{"isMC", false, "is MC"}; + Configurable isMCGen{"isMCGen", false, "is isMCGen"}; + Configurable isPrim{"isPrim", true, "is isPrim"}; Configurable mcCorrelation{"mcCorrelation", false, "true: build the correlation function only for SE"}; Configurable docorrection{"docorrection", false, "do efficiency correction"}; @@ -101,18 +103,19 @@ struct hadronnucleicorrelation { ConfigurableAxis AxisNSigma{"AxisNSigma", {35, -7.f, 7.f}, "n#sigma"}; using FilteredCollisions = soa::Filtered; - using FilteredTracks = soa::Filtered>; // old tables (v2) - using FilteredTracksMC = soa::Filtered>; // old tables (v2) - - // using FilteredTracks = soa::Filtered>; // new tables (v3) - // using FilteredTracksMC = soa::Filtered>; // new tables (v3) + using SimCollisions = aod::McCollisions; + using SimParticles = aod::McParticles; + using FilteredTracks = soa::Filtered>; // new tables (v3) + using FilteredTracksMC = soa::Filtered>; // new tables (v3) HistogramRegistry registry{"registry"}; HistogramRegistry QA{"QA"}; typedef std::shared_ptr trkType; typedef std::shared_ptr trkTypeMC; + typedef std::shared_ptr partTypeMC; typedef std::shared_ptr colType; + typedef std::shared_ptr MCcolType; // key: int64_t - value: vector of trkType objects std::map> selectedtracks_p; @@ -120,6 +123,12 @@ struct hadronnucleicorrelation { std::map> selectedtracks_antid; std::map> selectedtracks_antip; + // key: int64_t - value: vector of trkType objects + std::map> selectedparticlesMC_d; + std::map> selectedparticlesMC_p; + std::map> selectedparticlesMC_antid; + std::map> selectedparticlesMC_antip; + // key: int64_t - value: vector of trkType objects std::map> selectedtracksMC_d; std::map> selectedtracksMC_p; @@ -134,6 +143,8 @@ struct hadronnucleicorrelation { // for each key I have a vector of collisions std::map, std::vector> mixbins_antidantip; std::map, std::vector> mixbinsPID_antidantip; + std::map> mixbinsMC_antidantip; + std::map> mixbinsMC_dp; std::vector> hEtaPhi_AntiDeAntiPr_SE; std::vector> hEtaPhi_AntiDeAntiPr_ME; @@ -197,16 +208,10 @@ struct hadronnucleicorrelation { 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("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("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("hEtaPhiGen_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("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("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_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("hEtaPhi_AntiDeAntiPr_SE_pt%02.0f%02.0f", pTBins.value.at(i) * 10, pTBins.value.at(i + 1) * 10), Form("Raw #Delta#eta#Delta#phi (%.1f(HIST("Generated/hNEventsMC"))->GetXaxis()->SetBinLabel(1, "All"); + registry.get(HIST("Generated/hNEventsMC"))->GetXaxis()->SetBinLabel(2, "|z_{vtx}| < 10 cm"); + } } // Filters @@ -349,8 +371,8 @@ struct hadronnucleicorrelation { o2::aod::singletrackselector::unPack(o2::aod::singletrackselector::storedTpcChi2NCl) <= max_chi2_TPC && o2::aod::singletrackselector::unPack(o2::aod::singletrackselector::storedTpcCrossedRowsOverFindableCls) >= min_TPC_nCrossedRowsOverFindableCls && o2::aod::singletrackselector::unPack(o2::aod::singletrackselector::storedItsChi2NCl) <= max_chi2_ITS && - nabs(o2::aod::singletrackselector::unPack(o2::aod::singletrackselector::storedDcaXY)) <= max_dcaxy && // For now no filtering on the DCAxy or DCAz (casting not supported) - nabs(o2::aod::singletrackselector::unPack(o2::aod::singletrackselector::storedDcaXY)) <= max_dcaz && // For now no filtering on the DCAxy or DCAz (casting not supported) + nabs(o2::aod::singletrackselector::unPack(o2::aod::singletrackselector::storedDcaXY)) <= max_dcaxy && + nabs(o2::aod::singletrackselector::unPack(o2::aod::singletrackselector::storedDcaXY)) <= max_dcaz && nabs(o2::aod::singletrackselector::eta) <= etacut; template @@ -536,6 +558,32 @@ struct hadronnucleicorrelation { } // tracks 1 } + template + void mixMCParticles(Type const& particles1, Type const& particles2) + { // last value: 0 -- SE; 1 -- ME + for (auto it1 : particles1) { + for (auto it2 : particles2) { + + // Calculate Delta-eta Delta-phi (gen) + float deltaEtaGen = it2->eta() - it1->eta(); + float deltaPhiGen = it2->phi() - it1->phi(); + deltaPhiGen = getDeltaPhi(deltaPhiGen); + + for (int k = 0; k < nBinspT; k++) { + + if (it1->pt() >= pTBins.value.at(k) && it1->pt() < pTBins.value.at(k + 1)) { + + if (ME) { + hEtaPhiGen_AntiDeAntiPr_ME[k]->Fill(deltaEtaGen, deltaPhiGen, it2->pt()); + } else { + hEtaPhiGen_AntiDeAntiPr_SE[k]->Fill(deltaEtaGen, deltaPhiGen, it2->pt()); + } // SE + } // pT condition + } // nBinspT loop + } // particles 2 + } // particles 1 + } + float getDeltaPhi(float deltaPhi) { if (deltaPhi < -TMath::Pi() / 2) { @@ -1236,6 +1284,139 @@ struct hadronnucleicorrelation { mixbins_antidantip.clear(); // clear the map } PROCESS_SWITCH(hadronnucleicorrelation, processMC, "processMC", false); + + void processGen(SimCollisions const& mcCollisions, + SimParticles const& mcParticles) + { + for (auto particle : mcParticles) { + + if (isPrim && !particle.isPhysicalPrimary()) { + continue; + } + + 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 (auto collision : mcCollisions) { + + registry.fill(HIST("Generated/hNEventsMC"), 0.5); + if (std::abs(collision.posZ()) < 10.f) { + // return; + registry.fill(HIST("Generated/hNEventsMC"), 1.5); + } + + int vertexBinToMix = std::floor((collision.posZ() + cutzvertex) / (2 * cutzvertex / _vertexNbinsToMix)); + // int centBinToMix = std::floor(collision.centFT0M() / (100.0 / _multNsubBins)); + + if (selectedparticlesMC_antid.find(collision.globalIndex()) != selectedparticlesMC_antid.end()) { + mixbinsMC_antidantip[vertexBinToMix].push_back(std::make_shared(collision)); + } + if (selectedparticlesMC_d.find(collision.globalIndex()) != selectedparticlesMC_d.end()) { + mixbinsMC_dp[vertexBinToMix].push_back(std::make_shared(collision)); + } + } // coll + + if (!mixbinsMC_antidantip.empty()) { + + for (auto i = mixbinsMC_antidantip.begin(); i != mixbinsMC_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 (selectedparticlesMC_antip.find(col1->index()) != selectedparticlesMC_antip.end()) { + mixMCParticles<0>(selectedparticlesMC_antid[col1->index()], selectedparticlesMC_antip[col1->index()]); // 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 (selectedparticlesMC_antip.find(col2->index()) != selectedparticlesMC_antip.end()) { + mixMCParticles<1>(selectedparticlesMC_antid[col1->index()], selectedparticlesMC_antip[col2->index()]); // mixing ME + } + } + } + } + } + + if (!mixbinsMC_dp.empty()) { + + for (auto i = mixbinsMC_dp.begin(); i != mixbinsMC_dp.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()) { + mixMCParticles<0>(selectedparticlesMC_d[col1->index()], selectedparticlesMC_p[col1->index()]); // 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 (selectedparticlesMC_p.find(col2->index()) != selectedparticlesMC_p.end()) { + mixMCParticles<1>(selectedparticlesMC_d[col1->index()], selectedparticlesMC_p[col2->index()]); // mixing ME + } + } + } + } + } + + // clearing up + for (auto i = selectedparticlesMC_antid.begin(); i != selectedparticlesMC_antid.end(); i++) + (i->second).clear(); + selectedparticlesMC_antid.clear(); + + for (auto i = selectedparticlesMC_d.begin(); i != selectedparticlesMC_d.end(); i++) + (i->second).clear(); + selectedparticlesMC_d.clear(); + + for (auto i = selectedparticlesMC_antip.begin(); i != selectedparticlesMC_antip.end(); i++) + (i->second).clear(); + selectedparticlesMC_antip.clear(); + + for (auto i = selectedparticlesMC_p.begin(); i != selectedparticlesMC_p.end(); i++) + (i->second).clear(); + selectedparticlesMC_p.clear(); + + for (auto& pair : mixbinsMC_antidantip) { + pair.second.clear(); // clear the vector associated with the key + } + mixbinsMC_antidantip.clear(); // clear the map + + for (auto& pair : mixbinsMC_dp) { + pair.second.clear(); // clear the vector associated with the key + } + mixbinsMC_dp.clear(); // clear the map + } + PROCESS_SWITCH(hadronnucleicorrelation, processGen, "processGen", false); }; WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)