include of check for gamma gamma contamination in Dalitz decays

Laura Gansbartl · Laura Gansbartl · commit 8e432b0e2ef8 · 2026-05-20T19:11:03.000+02:00
diff --git a/PWGEM/PhotonMeson/Core/Pi0EtaToGammaGammaMC.h b/PWGEM/PhotonMeson/Core/Pi0EtaToGammaGammaMC.h
@@ -56,6 +56,7 @@
 #include <Math/Vector4Dfwd.h>
 #include <TF1.h>
 #include <TH2.h>
+#include <TPDGCode.h>
 #include <TString.h>
 
 #include <cmath>
@@ -74,24 +75,25 @@ enum AlphaMesonCutOption {
 
 template <o2::aod::pwgem::photonmeson::photonpair::PairType pairtype, o2::soa::is_table... Types>
 struct Pi0EtaToGammaGammaMC {
-  o2::framework::Configurable<std::string> ccdburl{"ccdb-url", "http://alice-ccdb.cern.ch", "url of the ccdb repository"};
+  o2::framework::Configurable<std::string> ccdbUrl{"ccdbUrl", "http://alice-ccdb.cern.ch", "url of the ccdb repository"};
   o2::framework::Configurable<std::string> grpPath{"grpPath", "GLO/GRP/GRP", "Path of the grp file"};
   o2::framework::Configurable<std::string> grpmagPath{"grpmagPath", "GLO/Config/GRPMagField", "CCDB path of the GRPMagField object"};
   o2::framework::Configurable<bool> skipGRPOquery{"skipGRPOquery", true, "skip grpo query"};
-  o2::framework::Configurable<float> d_bz_input{"d_bz_input", -999, "bz field in kG, -999 is automatic"};
+  o2::framework::Configurable<float> dBzInput{"dBzInput", -999, "bz field in kG, -999 is automatic"};
 
   o2::framework::Configurable<int> cfgQvecEstimator{"cfgQvecEstimator", 0, "FT0M:0, FT0A:1, FT0C:2"};
   o2::framework::Configurable<int> cfgCentEstimator{"cfgCentEstimator", 2, "FT0M:0, FT0A:1, FT0C:2"};
   o2::framework::Configurable<float> cfgCentMin{"cfgCentMin", -1, "min. centrality"};
   o2::framework::Configurable<float> cfgCentMax{"cfgCentMax", 999, "max. centrality"};
-  o2::framework::Configurable<float> maxY_rec{"maxY_rec", 0.9, "maximum rapidity for reconstructed particles"};
-  o2::framework::Configurable<std::string> fd_k0s_to_pi0{"fd_k0s_pi0", "1.0", "feed down correction to pi0"};
+  o2::framework::Configurable<float> maxYRec{"maxYRec", 0.9, "maximum rapidity for reconstructed particles"};
+  o2::framework::Configurable<std::string> fdK0sToPi0{"fdK0sPi0", "1.0", "feed down correction to pi0"};
   o2::framework::Configurable<bool> cfgRequireTrueAssociation{"cfgRequireTrueAssociation", false, "flag to require true mc collision association"};
 
   o2::framework::Configurable<int> cfgAlphaMesonCut{"cfgAlphaMesonCut", 0, "flag for photon energy asymmetry distribution cut: 0: no cut, 1: cut specific value, 2: cut depending on pT"};
   o2::framework::Configurable<float> cfgAlphaMeson{"cfgAlphaMeson", 0.65, "photon energy asymmetry distribution parameter for specific value cut"};
   o2::framework::Configurable<float> cfgAlphaMesonA{"cfgAlphaMesonA", 0.65, "photon energy asymmetry distribution parameter A for pT dependent cut (A * tanh(B*pT))"};
   o2::framework::Configurable<float> cfgAlphaMesonB{"cfgAlphaMesonB", 1.2, "photon energy asymmetry distribution parameter B for pT dependent cut (A * tanh(B*pT))"};
+  o2::framework::Configurable<bool> cfgGGContaCheck{"cfgGGContaCheck", false, "check gamma gamma contamination of dalitz"};
 
   EMPhotonEventCut fEMEventCut;
   struct : o2::framework::ConfigurableGroup {
@@ -224,7 +226,7 @@ struct Pi0EtaToGammaGammaMC {
     o2::framework::Configurable<float> cfg_min_Ecluster{"cfg_min_Ecluster", 0.3, "Minimum cluster energy for PHOS in GeV"};
   } phoscuts;
 
-  TF1* f1fd_k0s_to_pi0;
+  TF1* f1fdK0sToPi0;
   o2::framework::HistogramRegistry fRegistry{"output", {}, o2::framework::OutputObjHandlingPolicy::AnalysisObject, false, false};
   // static constexpr std::string_view event_types[2] = {"before/", "after/"};
   // static constexpr std::string_view event_pair_types[2] = {"same/", "mix/"};
@@ -249,14 +251,20 @@ struct Pi0EtaToGammaGammaMC {
     DefineEMCCut();
     DefinePHOSCut();
 
-    f1fd_k0s_to_pi0 = new TF1("f1fd_k0s_to_pi0", TString(fd_k0s_to_pi0), 0.f, 100.f);
+    f1fdK0sToPi0 = new TF1("f1fdK0sToPi0", TString(fdK0sToPi0), 0.f, 100.f);
 
     fRegistry.add("Event/hNrecPerMCCollision", "Nrec per mc collision;N_{rec} collisions per MC collision", o2::framework::HistType::kTH1F, {{21, -0.5f, 20.5f}}, false);
+    if (cfgGGContaCheck) {
+      fRegistry.add("Event/hNGGContamEta", "Number of Eta from etaToGammaGamma; p_{T, #eta} (GeV/#it{c}); N", o2::framework::HistType::kTH1F, {{40, -0.5f, 20.5f}}, false);
+      fRegistry.add("Event/hNGGContamPion", "Number of Pion from etaToGammaGamma; p_{T, #pi} (GeV/#it{c}); N", o2::framework::HistType::kTH1F, {{40, -0.5f, 20.5f}}, false);
+    }
+    fRegistry.add("Event/hNDalitzEtaPt", "Number of DalitzEta; p_{T, #eta} (GeV/#it{c}); N", o2::framework::HistType::kTH1F, {{40, -0.5f, 20.5f}}, false);
+    fRegistry.add("Event/hNDalitzPionPt", "Number of DalitzPion; p_{T, #pi} (GeV/#it{c}) ; N", o2::framework::HistType::kTH1F, {{40, -0.5f, 20.5f}}, false);
 
     mRunNumber = 0;
     d_bz = 0;
 
-    ccdb->setURL(ccdburl);
+    ccdb->setURL(ccdbUrl);
     ccdb->setCaching(true);
     ccdb->setLocalObjectValidityChecking();
     ccdb->setFatalWhenNull(false);
@@ -270,10 +278,12 @@ struct Pi0EtaToGammaGammaMC {
     }
 
     // In case override, don't proceed, please - no CCDB access required
-    if (d_bz_input > -990) {
-      d_bz = d_bz_input;
+    constexpr float bzInput = -990.0f;
+    if (dBzInput > bzInput) {
+      d_bz = dBzInput;
       o2::parameters::GRPMagField grpmag;
-      if (std::fabs(d_bz) > 1e-5) {
+      float bzThreshold = 1e-5;
+      if (std::fabs(d_bz) > bzThreshold) {
         grpmag.setL3Current(30000.f / (d_bz / 5.0f));
       }
       mRunNumber = collision.runNumber();
@@ -304,8 +314,8 @@ struct Pi0EtaToGammaGammaMC {
 
   ~Pi0EtaToGammaGammaMC()
   {
-    delete f1fd_k0s_to_pi0;
-    f1fd_k0s_to_pi0 = 0x0;
+    delete f1fdK0sToPi0;
+    f1fdK0sToPi0 = 0x0;
   }
 
   void DefineEMEventCut()
@@ -359,7 +369,7 @@ struct Pi0EtaToGammaGammaMC {
     fV0PhotonCut.SetLoadMlModelsFromCCDB(pcmcuts.cfg_load_ml_models_from_ccdb);
     fV0PhotonCut.SetNClassesMl(pcmcuts.cfg_nclasses_ml);
     fV0PhotonCut.SetMlTimestampCCDB(pcmcuts.cfg_timestamp_ccdb);
-    fV0PhotonCut.SetCcdbUrl(ccdburl);
+    fV0PhotonCut.SetCcdbUrl(ccdbUrl);
     CentType mCentralityTypeMlEnum;
     mCentralityTypeMlEnum = static_cast<CentType>(cfgCentEstimator.value);
     fV0PhotonCut.SetCentralityTypeMl(mCentralityTypeMlEnum);
@@ -572,18 +582,19 @@ struct Pi0EtaToGammaGammaMC {
                       TMCCollisions const& mccollisions, TMCParticles const& mcparticles,
                       TLegs const& /*legs*/ = nullptr, TMatchedTracks const& matchedTracks = nullptr, TMatchedSecondaries const& matchedSecondaries = nullptr)
   {
-    for (auto& collision : collisions) {
+    for (auto const& collision : collisions) {
       initCCDB(collision);
       if ((pairtype == o2::aod::pwgem::photonmeson::photonpair::PairType::kPHOSPHOS || pairtype == o2::aod::pwgem::photonmeson::photonpair::PairType::kPCMPHOS) && !collision.alias_bit(triggerAliases::kTVXinPHOS)) {
         continue;
       }
 
       float weight = 1.f;
+      float weightThreshold = 1e-10;
       if constexpr (std::is_same_v<std::decay_t<TCollisions>, o2::soa::Filtered<o2::soa::Join<o2::soa::Join<o2::aod::PMEvents, o2::aod::EMEventsAlias, o2::aod::EMEventsMult_000, o2::aod::EMEventsCent_000, o2::aod::EMMCEventLabels>, o2::aod::EMEventsWeight>>>) {
         weight = collision.weight();
       }
 
-      if (eventcuts.onlyKeepWeightedEvents && std::fabs(weight - 1.0) < 1e-10) {
+      if (eventcuts.onlyKeepWeightedEvents && std::fabs(weight - 1.0) < weightThreshold) {
         continue;
       }
 
@@ -680,7 +691,7 @@ struct Pi0EtaToGammaGammaMC {
           ROOT::Math::PtEtaPhiMVector v1(g1.pt(), g1.eta(), g1.phi(), 0.);
           ROOT::Math::PtEtaPhiMVector v2(g2.pt(), g2.eta(), g2.phi(), 0.);
           ROOT::Math::PtEtaPhiMVector v12 = v1 + v2;
-          if (std::fabs(v12.Rapidity()) > maxY_rec) {
+          if (std::fabs(v12.Rapidity()) > maxYRec) {
             continue;
           }
 
@@ -720,9 +731,9 @@ struct Pi0EtaToGammaGammaMC {
           }
 
           if (g1mc.globalIndex() == g2mc.globalIndex()) {
-            if (o2::aod::pwgem::dilepton::utils::mcutil::getMotherPDGCode(g1mc, mcparticles) == 111)
+            if (o2::aod::pwgem::dilepton::utils::mcutil::getMotherPDGCode(g1mc, mcparticles) == PDG_t::kPi0)
               fRegistry.fill(HIST("Pair/Pi0/hs_FromSameGamma"), v12.M(), v12.Pt(), wpair);
-            else if (o2::aod::pwgem::dilepton::utils::mcutil::getMotherPDGCode(g1mc, mcparticles) == 221)
+            else if (o2::aod::pwgem::dilepton::utils::mcutil::getMotherPDGCode(g1mc, mcparticles) == o2::constants::physics::Pdg::kEta)
               fRegistry.fill(HIST("Pair/Eta/hs_FromSameGamma"), v12.M(), v12.Pt(), wpair);
             continue;
           }
@@ -732,13 +743,13 @@ struct Pi0EtaToGammaGammaMC {
             if (cfgRequireTrueAssociation && (pi0mc.emmceventId() != collision.emmceventId())) {
               continue;
             }
-            o2::aod::pwgem::photonmeson::utils::nmhistogram::fillTruePairInfo(&fRegistry, v12, pi0mc, mcparticles, mccollisions, f1fd_k0s_to_pi0, wpair);
+            o2::aod::pwgem::photonmeson::utils::nmhistogram::fillTruePairInfo(&fRegistry, v12, pi0mc, mcparticles, mccollisions, f1fdK0sToPi0, wpair);
           } else if (etaid > 0) {
             auto etamc = mcparticles.iteratorAt(etaid);
             if (cfgRequireTrueAssociation && (etamc.emmceventId() != collision.emmceventId())) {
               continue;
             }
-            o2::aod::pwgem::photonmeson::utils::nmhistogram::fillTruePairInfo(&fRegistry, v12, etamc, mcparticles, mccollisions, f1fd_k0s_to_pi0, wpair);
+            o2::aod::pwgem::photonmeson::utils::nmhistogram::fillTruePairInfo(&fRegistry, v12, etamc, mcparticles, mccollisions, f1fdK0sToPi0, wpair);
           }
         } // end of pairing loop
       } else if constexpr (pairtype == o2::aod::pwgem::photonmeson::photonpair::PairType::kPCMDalitzEE) {
@@ -792,6 +803,25 @@ struct Pi0EtaToGammaGammaMC {
 
             auto pos2mc = mcparticles.iteratorAt(pos2.emmcparticleId());
             auto ele2mc = mcparticles.iteratorAt(ele2.emmcparticleId());
+            if (cfgGGContaCheck) {
+              photonid2 = o2::aod::pwgem::dilepton::utils::mcutil::FindCommonMotherFrom2Prongs(pos2mc, ele2mc, -11, 11, 22, mcparticles); // check possible contamination
+              if (photonid2 > 0) {
+                auto photon2 = mcparticles.iteratorAt(photonid2);
+                int photon2pdg = photon2.pdgCode();
+                int photon2mothid = photon2.mothersIds()[0];
+                auto photon2moth = mcparticles.iteratorAt(photon2mothid);
+                if (photon2pdg == PDG_t::kGamma && (o2::aod::pwgem::photonmeson::utils::mcutil::isGammaGammaDecay(photon2moth, mcparticles))) {
+                  int mothID = o2::aod::pwgem::dilepton::utils::mcutil::getMotherPDGCode(photon2, mcparticles);
+                  if (mothID == o2::constants::physics::Pdg::kEta) {
+                    fRegistry.fill(HIST("Event/hNGGContamEta"), photon2moth.pt());
+                  }
+                  if (mothID == PDG_t::kPi0) {
+                    fRegistry.fill(HIST("Event/hNGGContamPion"), photon2moth.pt());
+                  }
+                }
+              }
+            }
+
             pi0id = o2::aod::pwgem::dilepton::utils::mcutil::FindCommonMotherFrom3Prongs(g1mc, pos2mc, ele2mc, 22, -11, 11, 111, mcparticles);
             etaid = o2::aod::pwgem::dilepton::utils::mcutil::FindCommonMotherFrom3Prongs(g1mc, pos2mc, ele2mc, 22, -11, 11, 221, mcparticles);
             if (pi0id < 0 && etaid < 0) {
@@ -800,21 +830,23 @@ struct Pi0EtaToGammaGammaMC {
             ROOT::Math::PtEtaPhiMVector v_pos(pos2.pt(), pos2.eta(), pos2.phi(), o2::constants::physics::MassElectron);
             ROOT::Math::PtEtaPhiMVector v_ele(ele2.pt(), ele2.eta(), ele2.phi(), o2::constants::physics::MassElectron);
             ROOT::Math::PtEtaPhiMVector veeg = v_gamma + v_pos + v_ele;
-            if (std::fabs(veeg.Rapidity()) > maxY_rec) {
+            if (std::fabs(veeg.Rapidity()) > maxYRec) {
               continue;
             }
             if (pi0id > 0) {
               auto pi0mc = mcparticles.iteratorAt(pi0id);
+              fRegistry.fill(HIST("Event/hNDalitzPionPt"), pi0mc.pt());
               if (cfgRequireTrueAssociation && (pi0mc.emmceventId() != collision.emmceventId())) {
                 continue;
               }
-              o2::aod::pwgem::photonmeson::utils::nmhistogram::fillTruePairInfo(&fRegistry, veeg, pi0mc, mcparticles, mccollisions, f1fd_k0s_to_pi0, weight);
+              o2::aod::pwgem::photonmeson::utils::nmhistogram::fillTruePairInfo(&fRegistry, veeg, pi0mc, mcparticles, mccollisions, f1fdK0sToPi0, weight);
             } else if (etaid > 0) {
               auto etamc = mcparticles.iteratorAt(etaid);
+              fRegistry.fill(HIST("Event/hNDalitzEtaPt"), etamc.pt());
               if (cfgRequireTrueAssociation && (etamc.emmceventId() != collision.emmceventId())) {
                 continue;
               }
-              o2::aod::pwgem::photonmeson::utils::nmhistogram::fillTruePairInfo(&fRegistry, veeg, etamc, mcparticles, mccollisions, f1fd_k0s_to_pi0, weight);
+              o2::aod::pwgem::photonmeson::utils::nmhistogram::fillTruePairInfo(&fRegistry, veeg, etamc, mcparticles, mccollisions, f1fdK0sToPi0, weight);
             }
           } // end of dielectron loop
         } // end of pcm loop
@@ -853,15 +885,15 @@ struct Pi0EtaToGammaGammaMC {
           ROOT::Math::PtEtaPhiMVector v1(g1.pt(), g1.eta(), g1.phi(), 0.);
           ROOT::Math::PtEtaPhiMVector v2(g2.pt(), g2.eta(), g2.phi(), 0.);
           ROOT::Math::PtEtaPhiMVector v12 = v1 + v2;
-          if (std::fabs(v12.Rapidity()) > maxY_rec) {
+          if (std::fabs(v12.Rapidity()) > maxYRec) {
             continue;
           }
           // if (pi0id > 0) {
           //   auto pi0mc = mcparticles.iteratorAt(pi0id);
-          //   o2::aod::pwgem::photonmeson::utils::nmhistogram::fillTruePairInfo(&fRegistry, v12, pi0mc, mcparticles, mccollisions, f1fd_k0s_to_pi0, weight);
+          //   o2::aod::pwgem::photonmeson::utils::nmhistogram::fillTruePairInfo(&fRegistry, v12, pi0mc, mcparticles, mccollisions, f1fdK0sToPi0, weight);
           // } else if (etaid > 0) {
           //   auto etamc = mcparticles.iteratorAt(etaid);
-          //   o2::aod::pwgem::photonmeson::utils::nmhistogram::fillTruePairInfo(&fRegistry, v12, etamc, mcparticles, mccollisions, f1fd_k0s_to_pi0, weight);
+          //   o2::aod::pwgem::photonmeson::utils::nmhistogram::fillTruePairInfo(&fRegistry, v12, etamc, mcparticles, mccollisions, f1fdK0sToPi0, weight);
           // }
         } // end of pairing loop
       } // end of pairing in same event
@@ -900,23 +932,24 @@ struct Pi0EtaToGammaGammaMC {
     // loop over mc stack and fill histograms for pure MC truth signals
     // all MC tracks which belong to the MC event corresponding to the current reconstructed event
 
-    for (auto& mccollision : mccollisions) {
+    for (auto const& mccollision : mccollisions) {
       auto collision_per_mccoll = collisions.sliceBy(rec_perMcCollision, mccollision.globalIndex());
       int nrec_per_mc = collision_per_mccoll.size();
       fRegistry.fill(HIST("Event/hNrecPerMCCollision"), nrec_per_mc);
     }
 
-    for (auto& collision : collisions) {
+    for (auto const& collision : collisions) {
       if ((pairtype == o2::aod::pwgem::photonmeson::photonpair::kPHOSPHOS || pairtype == o2::aod::pwgem::photonmeson::photonpair::kPCMPHOS) && !collision.alias_bit(triggerAliases::kTVXinPHOS)) {
         continue; // I don't know why this is necessary in simulation.
       }
 
       float weight = 1.f;
+      float weightTreshhold = 1e-10;
       if constexpr (std::is_same_v<std::decay_t<TCollisions>, o2::soa::Filtered<o2::soa::Join<o2::soa::Join<o2::aod::PMEvents, o2::aod::EMEventsAlias, o2::aod::EMEventsMult_000, o2::aod::EMEventsCent_000, o2::aod::EMMCEventLabels>, o2::aod::EMEventsWeight>>>) {
         weight = collision.weight();
       }
 
-      if (eventcuts.onlyKeepWeightedEvents && std::fabs(weight - 1.0) < 1e-10) {
+      if (eventcuts.onlyKeepWeightedEvents && std::fabs(weight - 1.0) < weightTreshhold) {
         continue;
       }
 
