[PWGEM] Adding histograms for accepted pi0 and eta into gamma gamma decay cha… (#15465)

Author: amaringarcia

PWGEM/PhotonMeson/TableProducer/associateMCinfoPhoton.cxx

@@ -65,6 +65,7 @@ struct AssociateMCInfoPhoton {
   Configurable<float> margin_z_gen{"margin_z_gen", 15.f, "margin for Z of true photon conversion point to store generated information"};
   Configurable<float> max_rxy_gen{"max_rxy_gen", 100, "max rxy to store generated information"};
   Configurable<bool> requireGammaGammaDecay{"requireGammaGammaDecay", false, "require gamma gamma decay for generated pi0 and eta meson"};
+  Configurable<float> cfg_max_eta_photon{"cfg_max_eta_gamma", 0.8, "max eta for photons at PV"};
 
   HistogramRegistry registry{"EMMCEvent"};
 
@@ -98,6 +99,7 @@ struct AssociateMCInfoPhoton {
 
     for (uint i = 0; i < NParticleNames; i++) {
       registry.add<TH2>(Form("Generated/h2PtY_%s", ParticleNames[i].data()), Form("Generated %s", ParticleNames[i].data()), kTH2F, {axisPt, axisRapidity}, true);
+      registry.add<TH2>(Form("Generated/h2PtY_Accepted_%s", ParticleNames[i].data()), Form("Generated %s", ParticleNames[i].data()), kTH2F, {axisPt, axisRapidity}, true);
     }
 
     // reserve space for generated vectors if that process enabled
@@ -170,6 +172,15 @@ struct AssociateMCInfoPhoton {
       for (const auto& mcParticle : groupedMcParticles) { // store necessary information for denominator of efficiency
         if ((mcParticle.isPhysicalPrimary() || mcParticle.producedByGenerator()) && std::fabs(mcParticle.y()) < 0.9f && mcParticle.pt() < 20.f) {
           auto binNumber = hBinFinder->FindBin(mcParticle.pt(), std::fabs(mcParticle.y())); // caution: pack
+
+          bool isMesonAccepted = false;
+          if (mcParticle.has_daughters()) {
+            auto lDaughters = mcParticle.daughters_as<aod::McParticles>();
+
+            if (areTwoPhotonDaughtersAccepted(lDaughters, cfg_max_eta_photon)) {
+              isMesonAccepted = true;
+            }
+          }
           switch (std::abs(mcParticle.pdgCode())) {
             case PDG_t::kGamma:
               registry.fill(HIST("Generated/h2PtY_Gamma"), mcParticle.pt(), std::fabs(mcParticle.y()));
@@ -180,12 +191,18 @@ struct AssociateMCInfoPhoton {
                 continue;
               registry.fill(HIST("Generated/h2PtY_Pi0"), mcParticle.pt(), std::fabs(mcParticle.y()));
               genPi0[binNumber]++;
+              if (isMesonAccepted) {
+                registry.fill(HIST("Generated/h2PtY_Accepted_Pi0"), mcParticle.pt(), std::fabs(mcParticle.y()));
+              }
               break;
             case Pdg::kEta:
               if (requireGammaGammaDecay && !isGammaGammaDecay(mcParticle, mcParticles))
                 continue;
               registry.fill(HIST("Generated/h2PtY_Eta"), mcParticle.pt(), std::fabs(mcParticle.y()));
               genEta[binNumber]++;
+              if (isMesonAccepted) {
+                registry.fill(HIST("Generated/h2PtY_Accepted_Eta"), mcParticle.pt(), std::fabs(mcParticle.y()));
+              }
               break;
             default:
               break;
@@ -492,6 +509,21 @@ struct AssociateMCInfoPhoton {
     fCounters[1] = 0;
   } // end of skimmingMC
 
+  template <typename Daughters>
+  inline bool areTwoPhotonDaughtersAccepted(const Daughters& lDaughters, float maxEta)
+  {
+    if (lDaughters.size() != 2) {
+      return false;
+    }
+
+    auto it0 = lDaughters.begin();
+    auto it1 = it0;
+    ++it1;
+
+    return (std::fabs(it0.eta()) < maxEta &&
+            std::fabs(it1.eta()) < maxEta);
+  }
+
   void processMC_PCM(MyCollisionsMC const& collisions, aod::BCs const& bcs, aod::McCollisions const& mccollisions, aod::McParticles const& mcParticles, TracksMC const& o2tracks, aod::V0PhotonsKF const& v0photons, aod::V0Legs const& v0legs)
   {
     skimmingMC<kPCM>(collisions, bcs, mccollisions, mcParticles, o2tracks, nullptr, v0photons, v0legs, nullptr, nullptr, nullptr);