PWGEM/PhotonMeson: Include true pi0 for EMC with converted photons

- Add function that checks the decay chain of a cluster to see, whether it originated from a gamma from a pi0 or eta decay - Include true pi0s in output whose photons converted in the detector material (mostly TOF)
nstrangm · Apr 4, 2024 · 4b21b7a · 4b21b7a
1 parent 58cf71b
commit 4b21b7a
Show file tree

Hide file tree

Showing 2 changed files with 34 additions and 15 deletions.
diff --git 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<MyMCV0Legs>();
                 auto ele1 = g1.template negTrack_as<MyMCV0Legs>();
@@ -445,24 +447,25 @@ struct Pi0EtaToGammaGammaMC {
                 auto ele2mc = ele2.template emmcparticle_as<aod::EMMCParticles>();
                 // 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<int>{111, 221});
+                photonid2 = FindMotherInChain(cluster2mcparticle, mcparticles, std::vector<int>{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
@@ -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 <typename T, typename TMCs, typename TTargetPDGs>
+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 <typename T, typename TMCs>
 int IsEleFromPC(T const& mctrack, TMCs const& mcTracks)
 {