From ae16d723a96c24c2cfb3ee74761a683f238f4d05 Mon Sep 17 00:00:00 2001 From: Stefano Cannito Date: Fri, 11 Jul 2025 11:46:06 +0200 Subject: [PATCH 1/2] Substituted collision loop with iterator --- .../Tasks/Strangeness/phik0shortanalysis.cxx | 240 +++++++++--------- 1 file changed, 116 insertions(+), 124 deletions(-) diff --git a/PWGLF/Tasks/Strangeness/phik0shortanalysis.cxx b/PWGLF/Tasks/Strangeness/phik0shortanalysis.cxx index e3f1c0dcdae..568f8291675 100644 --- a/PWGLF/Tasks/Strangeness/phik0shortanalysis.cxx +++ b/PWGLF/Tasks/Strangeness/phik0shortanalysis.cxx @@ -2691,168 +2691,160 @@ struct Phik0shortanalysis { PROCESS_SWITCH(Phik0shortanalysis, processPhiK0SPionMCClosure2D, "Process function for Phi-K0S and Phi-Pion Correlations in MCClosure2D", false); - void processAllPartMCReco(SimCollisions const& collisions, FullMCTracks const& fullMCTracks, FullMCV0s const& V0s, V0DauMCTracks const&, MCCollisions const&, aod::McParticles const& mcParticles) + void processAllPartMCReco(SimCollisions::iterator const& collision, FullMCTracks const& fullMCTracks, FullMCV0s const& V0s, V0DauMCTracks const&, MCCollisions const&, aod::McParticles const& mcParticles) { - for (const auto& collision : collisions) { - if (!acceptEventQA(collision, false)) - continue; + if (!acceptEventQA(collision, false)) + return; - if (!collision.has_mcCollision()) - continue; + if (!collision.has_mcCollision()) + return; - const auto& mcCollision = collision.mcCollision_as(); - float genmultiplicity = mcCollision.centFT0M(); + const auto& mcCollision = collision.mcCollision_as(); + float genmultiplicity = mcCollision.centFT0M(); - // Defining positive and negative tracks for phi reconstruction - auto posThisColl = posMCTracks->sliceByCached(aod::track::collisionId, collision.globalIndex(), cache); - auto negThisColl = negMCTracks->sliceByCached(aod::track::collisionId, collision.globalIndex(), cache); + // Defining positive and negative tracks for phi reconstruction + auto posThisColl = posMCTracks->sliceByCached(aod::track::collisionId, collision.globalIndex(), cache); + auto negThisColl = negMCTracks->sliceByCached(aod::track::collisionId, collision.globalIndex(), cache); - for (const auto& track1 : posThisColl) { // loop over all selected tracks - if (!selectionTrackResonance(track1, false) || !selectionPIDKaonpTdependent(track1)) - continue; // topological and PID selection + for (const auto& track1 : posThisColl) { // loop over all selected tracks + if (!selectionTrackResonance(track1, false) || !selectionPIDKaonpTdependent(track1)) + continue; // topological and PID selection - auto track1ID = track1.globalIndex(); + auto track1ID = track1.globalIndex(); - if (!track1.has_mcParticle()) - continue; - auto mcTrack1 = track1.mcParticle_as(); - if (mcTrack1.pdgCode() != PDG_t::kKPlus || !mcTrack1.isPhysicalPrimary()) - continue; + if (!track1.has_mcParticle()) + continue; + auto mcTrack1 = track1.mcParticle_as(); + if (mcTrack1.pdgCode() != PDG_t::kKPlus || !mcTrack1.isPhysicalPrimary()) + continue; - for (const auto& track2 : negThisColl) { - if (!selectionTrackResonance(track2, false) || !selectionPIDKaonpTdependent(track2)) - continue; // topological and PID selection + for (const auto& track2 : negThisColl) { + if (!selectionTrackResonance(track2, false) || !selectionPIDKaonpTdependent(track2)) + continue; // topological and PID selection - auto track2ID = track2.globalIndex(); - if (track2ID == track1ID) - continue; // condition to avoid double counting of pair + auto track2ID = track2.globalIndex(); + if (track2ID == track1ID) + continue; // condition to avoid double counting of pair - if (!track2.has_mcParticle()) - continue; - auto mcTrack2 = track2.mcParticle_as(); - if (mcTrack2.pdgCode() != PDG_t::kKMinus || !mcTrack2.isPhysicalPrimary()) - continue; + if (!track2.has_mcParticle()) + continue; + auto mcTrack2 = track2.mcParticle_as(); + if (mcTrack2.pdgCode() != PDG_t::kKMinus || !mcTrack2.isPhysicalPrimary()) + continue; - float pTMother = -1.0f; - float yMother = -1.0f; - bool isMCMotherPhi = false; - for (const auto& motherOfMcTrack1 : mcTrack1.mothers_as()) { - for (const auto& motherOfMcTrack2 : mcTrack2.mothers_as()) { - if (motherOfMcTrack1.pdgCode() != motherOfMcTrack2.pdgCode()) - continue; - if (motherOfMcTrack1.globalIndex() != motherOfMcTrack2.globalIndex()) - continue; - if (motherOfMcTrack1.pdgCode() != o2::constants::physics::Pdg::kPhi) - continue; + float pTMother = -1.0f; + float yMother = -1.0f; + bool isMCMotherPhi = false; + for (const auto& motherOfMcTrack1 : mcTrack1.mothers_as()) { + for (const auto& motherOfMcTrack2 : mcTrack2.mothers_as()) { + if (motherOfMcTrack1.pdgCode() != motherOfMcTrack2.pdgCode()) + continue; + if (motherOfMcTrack1.globalIndex() != motherOfMcTrack2.globalIndex()) + continue; + if (motherOfMcTrack1.pdgCode() != o2::constants::physics::Pdg::kPhi) + continue; - pTMother = motherOfMcTrack1.pt(); - yMother = motherOfMcTrack1.y(); - isMCMotherPhi = true; - } + pTMother = motherOfMcTrack1.pt(); + yMother = motherOfMcTrack1.y(); + isMCMotherPhi = true; } - - if (!isMCMotherPhi) - continue; - if (pTMother < minPhiPt || std::abs(yMother) > cfgYAcceptance) - continue; - - mcPhiHist.fill(HIST("h3PhiMCRecoNewProc"), genmultiplicity, pTMother, yMother); } - } - // Defining V0s in the collision - auto v0sThisColl = V0s.sliceByCached(aod::track::collisionId, collision.globalIndex(), cache); - - for (const auto& v0 : v0sThisColl) { - if (!v0.has_mcParticle()) + if (!isMCMotherPhi) continue; - - auto v0mcparticle = v0.mcParticle(); - if (v0mcparticle.pdgCode() != PDG_t::kK0Short || !v0mcparticle.isPhysicalPrimary()) + if (pTMother < minPhiPt || std::abs(yMother) > cfgYAcceptance) continue; - const auto& posDaughterTrack = v0.posTrack_as(); - const auto& negDaughterTrack = v0.negTrack_as(); + mcPhiHist.fill(HIST("h3PhiMCRecoNewProc"), genmultiplicity, pTMother, yMother); + } + } - if (!selectionV0(v0, posDaughterTrack, negDaughterTrack)) - continue; - if (v0Configs.cfgFurtherV0Selection && !furtherSelectionV0(v0, collision)) - continue; - if (std::abs(v0mcparticle.y()) > cfgYAcceptance) - continue; + for (const auto& v0 : V0s) { + if (!v0.has_mcParticle()) + continue; - mcK0SHist.fill(HIST("h3K0SMCRecoNewProc"), genmultiplicity, v0mcparticle.pt(), v0mcparticle.y()); - } + auto v0mcparticle = v0.mcParticle(); + if (v0mcparticle.pdgCode() != PDG_t::kK0Short || !v0mcparticle.isPhysicalPrimary()) + continue; - // Defining tracks in the collision - auto mcTracksThisColl = fullMCTracks.sliceByCached(aod::track::collisionId, collision.globalIndex(), cache); + const auto& posDaughterTrack = v0.posTrack_as(); + const auto& negDaughterTrack = v0.negTrack_as(); - for (const auto& track : mcTracksThisColl) { - // Pion selection - if (!selectionPion(track, false)) - continue; + if (!selectionV0(v0, posDaughterTrack, negDaughterTrack)) + continue; + if (v0Configs.cfgFurtherV0Selection && !furtherSelectionV0(v0, collision)) + continue; + if (std::abs(v0mcparticle.y()) > cfgYAcceptance) + continue; - if (!track.has_mcParticle()) - continue; + mcK0SHist.fill(HIST("h3K0SMCRecoNewProc"), genmultiplicity, v0mcparticle.pt(), v0mcparticle.y()); + } - auto mcTrack = track.mcParticle_as(); - if (std::abs(mcTrack.pdgCode()) != PDG_t::kPiPlus) - continue; + for (const auto& track : fullMCTracks) { + // Pion selection + if (!selectionPion(track, false)) + continue; - if (std::abs(mcTrack.y()) > cfgYAcceptance) - continue; + if (!track.has_mcParticle()) + continue; - // Primary pion selection - if (mcTrack.isPhysicalPrimary()) { - mcPionHist.fill(HIST("h3RecMCDCAxyPrimPi"), track.pt(), track.dcaXY()); - } else { - if (mcTrack.getProcess() == 4) { // Selection of secondary pions from weak decay - mcPionHist.fill(HIST("h3RecMCDCAxySecWeakDecayPi"), track.pt(), track.dcaXY()); - } else { // Selection of secondary pions from material interactions - mcPionHist.fill(HIST("h3RecMCDCAxySecMaterialPi"), track.pt(), track.dcaXY()); - } - continue; + auto mcTrack = track.mcParticle_as(); + if (std::abs(mcTrack.pdgCode()) != PDG_t::kPiPlus) + continue; + + if (std::abs(mcTrack.y()) > cfgYAcceptance) + continue; + + // Primary pion selection + if (mcTrack.isPhysicalPrimary()) { + mcPionHist.fill(HIST("h3RecMCDCAxyPrimPi"), track.pt(), track.dcaXY()); + } else { + if (mcTrack.getProcess() == 4) { // Selection of secondary pions from weak decay + mcPionHist.fill(HIST("h3RecMCDCAxySecWeakDecayPi"), track.pt(), track.dcaXY()); + } else { // Selection of secondary pions from material interactions + mcPionHist.fill(HIST("h3RecMCDCAxySecMaterialPi"), track.pt(), track.dcaXY()); } + continue; + } - mcPionHist.fill(HIST("h3PiMCRecoNewProc"), genmultiplicity, mcTrack.pt(), mcTrack.y()); + mcPionHist.fill(HIST("h3PiMCRecoNewProc"), genmultiplicity, mcTrack.pt(), mcTrack.y()); - if (track.pt() >= trackConfigs.pTToUseTOF && !track.hasTOF()) - continue; + if (track.pt() >= trackConfigs.pTToUseTOF && !track.hasTOF()) + continue; - mcPionHist.fill(HIST("h3PiMCReco2NewProc"), genmultiplicity, mcTrack.pt(), mcTrack.y()); - } + mcPionHist.fill(HIST("h3PiMCReco2NewProc"), genmultiplicity, mcTrack.pt(), mcTrack.y()); + } - // Defining McParticles in the collision - auto mcParticlesThisColl = mcParticles.sliceByCached(aod::mcparticle::mcCollisionId, mcCollision.globalIndex(), cache); + // Defining McParticles in the collision + auto mcParticlesThisColl = mcParticles.sliceByCached(aod::mcparticle::mcCollisionId, mcCollision.globalIndex(), cache); - for (const auto& mcParticle : mcParticlesThisColl) { - if (std::abs(mcParticle.y()) > cfgYAcceptance) - continue; + for (const auto& mcParticle : mcParticlesThisColl) { + if (std::abs(mcParticle.y()) > cfgYAcceptance) + continue; - // Phi selection - if (mcParticle.pdgCode() != o2::constants::physics::Pdg::kPhi) - continue; - if (mcParticle.pt() < minPhiPt) - continue; + // Phi selection + if (mcParticle.pdgCode() != o2::constants::physics::Pdg::kPhi) + continue; + if (mcParticle.pt() < minPhiPt) + continue; - mcPhiHist.fill(HIST("h3PhiMCGenAssocRecoCheckNewProc"), genmultiplicity, mcParticle.pt(), mcParticle.y()); + mcPhiHist.fill(HIST("h3PhiMCGenAssocRecoCheckNewProc"), genmultiplicity, mcParticle.pt(), mcParticle.y()); - // K0S selection - if (mcParticle.pdgCode() != PDG_t::kK0Short) - continue; - if (!mcParticle.isPhysicalPrimary() || mcParticle.pt() < v0Configs.v0SettingMinPt) - continue; + // K0S selection + if (mcParticle.pdgCode() != PDG_t::kK0Short) + continue; + if (!mcParticle.isPhysicalPrimary() || mcParticle.pt() < v0Configs.v0SettingMinPt) + continue; - mcK0SHist.fill(HIST("h3K0SMCGenAssocRecoCheckNewProc"), genmultiplicity, mcParticle.pt(), mcParticle.y()); + mcK0SHist.fill(HIST("h3K0SMCGenAssocRecoCheckNewProc"), genmultiplicity, mcParticle.pt(), mcParticle.y()); - // Pion selection - if (std::abs(mcParticle.pdgCode()) != PDG_t::kPiPlus) - continue; - if (!mcParticle.isPhysicalPrimary() || mcParticle.pt() < trackConfigs.cMinPionPtcut) - continue; + // Pion selection + if (std::abs(mcParticle.pdgCode()) != PDG_t::kPiPlus) + continue; + if (!mcParticle.isPhysicalPrimary() || mcParticle.pt() < trackConfigs.cMinPionPtcut) + continue; - mcPionHist.fill(HIST("h3PiMCGenAssocRecoCheckNewProc"), genmultiplicity, mcParticle.pt(), mcParticle.y()); - } + mcPionHist.fill(HIST("h3PiMCGenAssocRecoCheckNewProc"), genmultiplicity, mcParticle.pt(), mcParticle.y()); } } From 87aaa2401dd507c38dc634537386117b6dd85560 Mon Sep 17 00:00:00 2001 From: Stefano Cannito Date: Fri, 11 Jul 2025 14:43:34 +0200 Subject: [PATCH 2/2] Study signal split + histo renaming --- .../Tasks/Strangeness/phik0shortanalysis.cxx | 79 +++++++++++++++---- 1 file changed, 64 insertions(+), 15 deletions(-) diff --git a/PWGLF/Tasks/Strangeness/phik0shortanalysis.cxx b/PWGLF/Tasks/Strangeness/phik0shortanalysis.cxx index 568f8291675..b0fdf4f8389 100644 --- a/PWGLF/Tasks/Strangeness/phik0shortanalysis.cxx +++ b/PWGLF/Tasks/Strangeness/phik0shortanalysis.cxx @@ -339,7 +339,8 @@ struct Phik0shortanalysis { mcEventHist.add("h2RecMCEtaDistribution", "Eta vs multiplicity in MCReco", kTH2F, {binnedmultAxis, etaAxis}); mcEventHist.add("h2GenMCEtaDistribution", "Eta vs multiplicity in MCGen", kTH2F, {binnedmultAxis, etaAxis}); mcEventHist.add("h2GenMCEtaDistributionAssocReco", "Eta vs multiplicity in MCGen Assoc Reco", kTH2F, {binnedmultAxis, etaAxis}); - mcEventHist.add("h2GenMCEtaDistributionAssocReco2", "Eta vs multiplicity in MCGen Assoc Reco", kTH2F, {binnedmultAxis, etaAxis}); + mcEventHist.add("h2GenMCEtaDistributionReco", "Eta vs multiplicity in MCGen Reco", kTH2F, {binnedmultAxis, etaAxis}); + mcEventHist.add("h2GenMCEtaDistributionRecoCheck", "Eta vs multiplicity in MCGen Reco Check", kTH2F, {binnedmultAxis, etaAxis}); // Phi topological/PID cuts dataPhiHist.add("h2DauTracksPhiDCAxyPreCutData", "Dcaxy distribution vs pt before DCAxy cut", kTH2F, {{100, 0.0, 5.0, "#it{p}_{T} (GeV/#it{c})"}, {2000, -0.05, 0.05, "DCA_{xy} (cm)"}}); @@ -592,9 +593,13 @@ struct Phik0shortanalysis { mcK0SHist.add("h3K0SMCGenAssocRecoNewProc", "K0S in MCGen Associated MCReco", kTH3F, {binnedmultAxis, pTK0SAxis, yAxis}); mcPionHist.add("h3PiMCGenAssocRecoNewProc", "Pion in MCGen Associated MCReco", kTH3F, {binnedmultAxis, pTPiAxis, yAxis}); - mcPhiHist.add("h3PhiMCGenAssocRecoCheckNewProc", "Phi in MCGen Associated MCReco Check", kTH3F, {binnedmultAxis, pTPhiAxis, yAxis}); - mcK0SHist.add("h3K0SMCGenAssocRecoCheckNewProc", "K0S in MCGen Associated MCReco Check", kTH3F, {binnedmultAxis, pTK0SAxis, yAxis}); - mcPionHist.add("h3PiMCGenAssocRecoCheckNewProc", "Pion in MCGen Associated MCReco Check", kTH3F, {binnedmultAxis, pTPiAxis, yAxis}); + mcPhiHist.add("h3PhiMCGenRecoNewProc", "Phi in MCGen MCReco", kTH3F, {binnedmultAxis, pTPhiAxis, yAxis}); + mcK0SHist.add("h3K0SMCGenRecoNewProc", "K0S in MCGen MCReco", kTH3F, {binnedmultAxis, pTK0SAxis, yAxis}); + mcPionHist.add("h3PiMCGenRecoNewProc", "Pion in MCGen MCReco", kTH3F, {binnedmultAxis, pTPiAxis, yAxis}); + + mcPhiHist.add("h3PhiMCGenRecoCheckNewProc", "Phi in MCGen MCReco Check", kTH3F, {binnedmultAxis, pTPhiAxis, yAxis}); + mcK0SHist.add("h3K0SMCGenRecoCheckNewProc", "K0S in MCGen MCReco Check", kTH3F, {binnedmultAxis, pTK0SAxis, yAxis}); + mcPionHist.add("h3PiMCGenRecoCheckNewProc", "Pion in MCGen MCReco Check", kTH3F, {binnedmultAxis, pTPiAxis, yAxis}); // Initialize CCDB only if purity or efficiencies are requested in the task if (useCCDB) { @@ -2405,7 +2410,7 @@ struct Phik0shortanalysis { if (pdgTrack->Charge() == trackConfigs.cfgCutCharge) continue; - mcEventHist.fill(HIST("h2GenMCEtaDistributionAssocReco"), genmultiplicity, mcParticle.eta()); + mcEventHist.fill(HIST("h2GenMCEtaDistributionRecoCheck"), genmultiplicity, mcParticle.eta()); } } @@ -2420,17 +2425,32 @@ struct Phik0shortanalysis { if (!eventHasMCPhi(mcParticles)) return; - bool isAssocColl = false; + float genmultiplicity = mcCollision.centFT0M(); + + uint64_t numberAssocColl = 0; for (const auto& collision : collisions) { if (acceptEventQA(collision, false)) { - isAssocColl = true; - break; + mcEventHist.fill(HIST("hGenMCRecoMultiplicityPercent"), genmultiplicity); + + for (const auto& mcParticle : mcParticles) { + if (!mcParticle.isPhysicalPrimary() || std::abs(mcParticle.eta()) > trackConfigs.etaMax) + continue; + + auto pdgTrack = pdgDB->GetParticle(mcParticle.pdgCode()); + if (pdgTrack == nullptr) + continue; + if (pdgTrack->Charge() == trackConfigs.cfgCutCharge) + continue; + + mcEventHist.fill(HIST("h2GenMCEtaDistributionReco"), genmultiplicity, mcParticle.eta()); + } + + numberAssocColl++; } } - float genmultiplicity = mcCollision.centFT0M(); mcEventHist.fill(HIST("hGenMCMultiplicityPercent"), genmultiplicity); - if (isAssocColl) + if (numberAssocColl > 0) mcEventHist.fill(HIST("hGenMCAssocRecoMultiplicityPercent"), genmultiplicity); for (const auto& mcParticle : mcParticles) { @@ -2444,8 +2464,8 @@ struct Phik0shortanalysis { continue; mcEventHist.fill(HIST("h2GenMCEtaDistribution"), genmultiplicity, mcParticle.eta()); - if (isAssocColl) - mcEventHist.fill(HIST("h2GenMCEtaDistributionAssocReco2"), genmultiplicity, mcParticle.eta()); + if (numberAssocColl > 0) + mcEventHist.fill(HIST("h2GenMCEtaDistributionAssocReco"), genmultiplicity, mcParticle.eta()); } } @@ -2828,7 +2848,7 @@ struct Phik0shortanalysis { if (mcParticle.pt() < minPhiPt) continue; - mcPhiHist.fill(HIST("h3PhiMCGenAssocRecoCheckNewProc"), genmultiplicity, mcParticle.pt(), mcParticle.y()); + mcPhiHist.fill(HIST("h3PhiMCGenRecoCheckNewProc"), genmultiplicity, mcParticle.pt(), mcParticle.y()); // K0S selection if (mcParticle.pdgCode() != PDG_t::kK0Short) @@ -2836,7 +2856,7 @@ struct Phik0shortanalysis { if (!mcParticle.isPhysicalPrimary() || mcParticle.pt() < v0Configs.v0SettingMinPt) continue; - mcK0SHist.fill(HIST("h3K0SMCGenAssocRecoCheckNewProc"), genmultiplicity, mcParticle.pt(), mcParticle.y()); + mcK0SHist.fill(HIST("h3K0SMCGenRecoCheckNewProc"), genmultiplicity, mcParticle.pt(), mcParticle.y()); // Pion selection if (std::abs(mcParticle.pdgCode()) != PDG_t::kPiPlus) @@ -2844,7 +2864,7 @@ struct Phik0shortanalysis { if (!mcParticle.isPhysicalPrimary() || mcParticle.pt() < trackConfigs.cMinPionPtcut) continue; - mcPionHist.fill(HIST("h3PiMCGenAssocRecoCheckNewProc"), genmultiplicity, mcParticle.pt(), mcParticle.y()); + mcPionHist.fill(HIST("h3PiMCGenRecoCheckNewProc"), genmultiplicity, mcParticle.pt(), mcParticle.y()); } } @@ -2863,6 +2883,35 @@ struct Phik0shortanalysis { for (const auto& collision : collisions) { if (acceptEventQA(collision, false)) { mcEventHist.fill(HIST("hGenMCRecoMultiplicityPercent"), genmultiplicity); // Event split numerator + + for (const auto& mcParticle : mcParticles) { + // The inclusive number of particles is the signal loss denominator, + // while the number of associated particles is the signal loss numerator + if (std::abs(mcParticle.y()) > cfgYAcceptance) + continue; + + // Phi selection + if (mcParticle.pdgCode() != o2::constants::physics::Pdg::kPhi) + continue; + if (mcParticle.pt() < minPhiPt) + continue; + mcPhiHist.fill(HIST("h3PhiMCGenRecoNewProc"), genmultiplicity, mcParticle.pt(), mcParticle.y()); + + // K0S selection + if (mcParticle.pdgCode() != PDG_t::kK0Short) + continue; + if (!mcParticle.isPhysicalPrimary() || mcParticle.pt() < v0Configs.v0SettingMinPt) + continue; + mcK0SHist.fill(HIST("h3K0SMCGenRecoNewProc"), genmultiplicity, mcParticle.pt(), mcParticle.y()); + + // Pion selection + if (std::abs(mcParticle.pdgCode()) != PDG_t::kPiPlus) + continue; + if (!mcParticle.isPhysicalPrimary() || mcParticle.pt() < trackConfigs.cMinPionPtcut) + continue; + mcPionHist.fill(HIST("h3PiMCGenRecoNewProc"), genmultiplicity, mcParticle.pt(), mcParticle.y()); + } + numberAssocColl++; } }