diff --git a/PWGEM/PhotonMeson/Tasks/Pi0EtaToGammaGammaMC.cxx b/PWGEM/PhotonMeson/Tasks/Pi0EtaToGammaGammaMC.cxx index d76543291e8..d41ed946660 100644 --- a/PWGEM/PhotonMeson/Tasks/Pi0EtaToGammaGammaMC.cxx +++ b/PWGEM/PhotonMeson/Tasks/Pi0EtaToGammaGammaMC.cxx @@ -433,6 +433,8 @@ struct Pi0EtaToGammaGammaMC { continue; } + int photonid1 = -1; + int photonid2 = -1; if constexpr (pairtype == PairType::kPCMPCM) { // check 2 legs auto pos1 = g1.template posTrack_as(); auto ele1 = g1.template negTrack_as(); @@ -445,24 +447,25 @@ struct Pi0EtaToGammaGammaMC { auto ele2mc = ele2.template emmcparticle_as(); // LOGF(info,"pos1mc.globalIndex() = %d , ele1mc.globalIndex() = %d , pos2mc.globalIndex() = %d , ele2mc.globalIndex() = %d", pos1mc.globalIndex(), ele1mc.globalIndex(), pos2mc.globalIndex(), ele2mc.globalIndex()); - int photonid1 = FindCommonMotherFrom2Prongs(pos1mc, ele1mc, -11, 11, 22, mcparticles); - int photonid2 = FindCommonMotherFrom2Prongs(pos2mc, ele2mc, -11, 11, 22, mcparticles); + photonid1 = FindCommonMotherFrom2Prongs(pos1mc, ele1mc, -11, 11, 22, mcparticles); + photonid2 = FindCommonMotherFrom2Prongs(pos2mc, ele2mc, -11, 11, 22, mcparticles); + } else if constexpr (pairtype == PairType::kEMCEMC) { + auto cluster1mcparticle = mcparticles.iteratorAt(g1.emmcparticleId()); + auto cluster2mcparticle = mcparticles.iteratorAt(g2.emmcparticleId()); - // LOGF(info,"photonid1 = %d , photonid2 = %d", photonid1, photonid2); - if (photonid1 < 0 || photonid2 < 0) { - continue; - } + photonid1 = FindMotherInChain(cluster1mcparticle, mcparticles, std::vector{111, 221}); + photonid2 = FindMotherInChain(cluster2mcparticle, mcparticles, std::vector{111, 221}); + } - auto g1mc = mcparticles.iteratorAt(photonid1); - auto g2mc = mcparticles.iteratorAt(photonid2); - pi0id = FindCommonMotherFrom2Prongs(g1mc, g2mc, 22, 22, 111, mcparticles); - etaid = FindCommonMotherFrom2Prongs(g1mc, g2mc, 22, 22, 221, mcparticles); - } else if constexpr (pairtype == PairType::kEMCEMC) { - auto g1mc = mcparticles.iteratorAt(g1.emmcparticleId()); - auto g2mc = mcparticles.iteratorAt(g2.emmcparticleId()); - pi0id = FindCommonMotherFrom2Prongs(g1mc, g2mc, 22, 22, 111, mcparticles); - etaid = FindCommonMotherFrom2Prongs(g1mc, g2mc, 22, 22, 221, mcparticles); + if (photonid1 < 0 || photonid2 < 0) { + continue; } + + auto g1mc = mcparticles.iteratorAt(photonid1); + auto g2mc = mcparticles.iteratorAt(photonid2); + pi0id = FindCommonMotherFrom2Prongs(g1mc, g2mc, 22, 22, 111, mcparticles); + etaid = FindCommonMotherFrom2Prongs(g1mc, g2mc, 22, 22, 221, mcparticles); + if (pi0id < 0 && etaid < 0) { continue; } diff --git a/PWGEM/PhotonMeson/Utils/MCUtilities.h b/PWGEM/PhotonMeson/Utils/MCUtilities.h index 2956ef02d85..e3f8b45df95 100644 --- a/PWGEM/PhotonMeson/Utils/MCUtilities.h +++ b/PWGEM/PhotonMeson/Utils/MCUtilities.h @@ -85,6 +85,22 @@ int IsXFromY(T const& mctrack, TMCs const& mcTracks, const int pdgX, const int p return -1; } //_______________________________________________________________________ +// Go up the decay chain of a mcparticle looking for a mother with the given pdg codes, if found return this mothers daughter +// E.g. Find the gamma that was created in a pi0 or eta decay +template +int FindMotherInChain(T const& mcparticle, TMCs const& mcparticles, TTargetPDGs const& motherpdgs, const int Depth = 50) +{ + if (!mcparticle.has_mothers() || Depth < 1) + return -1; + + int motherid = mcparticle.mothersIds()[0]; + auto mother = mcparticles.iteratorAt(motherid); + if (std::find(motherpdgs.begin(), motherpdgs.end(), mother.pdgCode()) != motherpdgs.end()) + return mcparticle.globalIndex(); // The mother has the required pdg code, so return its daughters global mc particle code. + else + return FindMotherInChain(mother, mcparticles, motherpdgs, Depth - 1); +} +//_______________________________________________________________________ template int IsEleFromPC(T const& mctrack, TMCs const& mcTracks) {