diff --git a/PWGLF/Tasks/Strangeness/phik0shortanalysis.cxx b/PWGLF/Tasks/Strangeness/phik0shortanalysis.cxx index f7b80611267..8c65cca608d 100644 --- a/PWGLF/Tasks/Strangeness/phik0shortanalysis.cxx +++ b/PWGLF/Tasks/Strangeness/phik0shortanalysis.cxx @@ -142,7 +142,9 @@ struct Phik0shortanalysis { Configurable cfgiskNoITSROFrameBorder{"cfgiskNoITSROFrameBorder", false, "kNoITSROFrameBorder request on RecMC collisions"}; // Configurable for MC closure - Configurable cfgisGenMCForClosure{"cfgisGenMCForClosure", false, "isGenMCForClosure"}; + Configurable cfgisRecMCWPDGForClosure1{"cfgisRecMCWPDGForClosure1", false, "RecoMC with PDG Codes for Closure only for Associated particles"}; + Configurable cfgisRecMCWPDGForClosure2{"cfgisRecMCWPDGForClosure2", false, "RecoMC with PDG Codes for Closure"}; + Configurable cfgisGenMCForClosure{"cfgisGenMCForClosure", false, "GenMC for Closure"}; // Constants double massKa = o2::constants::physics::MassKPlus; @@ -1369,7 +1371,7 @@ struct Phik0shortanalysis { PROCESS_SWITCH(Phik0shortanalysis, processRecMCClosurePhiQA, "Process for ReCMCQA and Phi in RecMCClosure", false); - void processRecMCClosurePhiK0S(SimCollisions::iterator const& collision, FullMCTracks const&, FullV0s const& V0s, V0DauMCTracks const&, MCCollisions const&) + void processRecMCClosurePhiK0S(SimCollisions::iterator const& collision, FullMCTracks const&, FullMCV0s const& V0s, V0DauMCTracks const&, MCCollisions const&) { if (!acceptEventQA(collision, false)) return; @@ -1386,6 +1388,14 @@ struct Phik0shortanalysis { // V0 already reconstructed by the builder for (const auto& v0 : V0s) { + if (cfgisRecMCWPDGForClosure1) { + if (!v0.has_mcParticle()) + continue; + auto v0mcparticle = v0.mcParticle(); + if (v0mcparticle.pdgCode() != 310 || !v0mcparticle.isPhysicalPrimary()) + continue; + } + const auto& posDaughterTrack = v0.posTrack_as(); const auto& negDaughterTrack = v0.negTrack_as(); @@ -1415,6 +1425,35 @@ struct Phik0shortanalysis { if (track2ID == track1ID) continue; // condition to avoid double counting of pair + if (cfgisRecMCWPDGForClosure2) { + if (!track1.has_mcParticle()) + continue; + auto mcTrack1 = track1.mcParticle_as(); + if (mcTrack1.pdgCode() != 321 || !mcTrack1.isPhysicalPrimary()) + continue; + + if (!track2.has_mcParticle()) + continue; + auto mcTrack2 = track2.mcParticle_as(); + if (mcTrack2.pdgCode() != -321 || !mcTrack2.isPhysicalPrimary()) + continue; + + 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() != 333) + continue; + isMCMotherPhi = true; + } + } + if (!isMCMotherPhi) + continue; + } + TLorentzVector recPhi = recMother(track1, track2, massKa, massKa); if (recPhi.M() < lowMPhi || recPhi.M() > upMPhi) @@ -1461,6 +1500,13 @@ struct Phik0shortanalysis { // Loop over all primary pion candidates for (const auto& track : fullMCTracks) { + if (cfgisRecMCWPDGForClosure1) { + if (!track.has_mcParticle()) + continue; + auto mcTrack = track.mcParticle_as(); + if (std::abs(mcTrack.pdgCode()) != 211 || !mcTrack.isPhysicalPrimary()) + continue; + } // Pion selection if (!selectionPion(track)) @@ -1487,6 +1533,35 @@ struct Phik0shortanalysis { if (track2ID == track1ID) continue; // condition to avoid double counting of pair + if (cfgisRecMCWPDGForClosure2) { + if (!track1.has_mcParticle()) + continue; + auto mcTrack1 = track1.mcParticle_as(); + if (mcTrack1.pdgCode() != 321 || !mcTrack1.isPhysicalPrimary()) + continue; + + if (!track2.has_mcParticle()) + continue; + auto mcTrack2 = track2.mcParticle_as(); + if (mcTrack2.pdgCode() != -321 || !mcTrack2.isPhysicalPrimary()) + continue; + + 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() != 333) + continue; + isMCMotherPhi = true; + } + } + if (!isMCMotherPhi) + continue; + } + TLorentzVector recPhi = recMother(track1, track2, massKa, massKa); if (recPhi.M() < lowMPhi || recPhi.M() > upMPhi)