diff --git a/PWGEM/PhotonMeson/Core/Pi0EtaToGammaGammaMC.h b/PWGEM/PhotonMeson/Core/Pi0EtaToGammaGammaMC.h index 2464ec82ec2..6badb8c6be2 100644 --- a/PWGEM/PhotonMeson/Core/Pi0EtaToGammaGammaMC.h +++ b/PWGEM/PhotonMeson/Core/Pi0EtaToGammaGammaMC.h @@ -56,6 +56,7 @@ #include #include #include +#include #include #include @@ -74,24 +75,25 @@ enum AlphaMesonCutOption { template struct Pi0EtaToGammaGammaMC { - o2::framework::Configurable ccdburl{"ccdb-url", "http://alice-ccdb.cern.ch", "url of the ccdb repository"}; + o2::framework::Configurable ccdbUrl{"ccdbUrl", "http://alice-ccdb.cern.ch", "url of the ccdb repository"}; o2::framework::Configurable grpPath{"grpPath", "GLO/GRP/GRP", "Path of the grp file"}; o2::framework::Configurable grpmagPath{"grpmagPath", "GLO/Config/GRPMagField", "CCDB path of the GRPMagField object"}; o2::framework::Configurable skipGRPOquery{"skipGRPOquery", true, "skip grpo query"}; - o2::framework::Configurable d_bz_input{"d_bz_input", -999, "bz field in kG, -999 is automatic"}; + o2::framework::Configurable dBzInput{"dBzInput", -999, "bz field in kG, -999 is automatic"}; o2::framework::Configurable cfgQvecEstimator{"cfgQvecEstimator", 0, "FT0M:0, FT0A:1, FT0C:2"}; o2::framework::Configurable cfgCentEstimator{"cfgCentEstimator", 2, "FT0M:0, FT0A:1, FT0C:2"}; o2::framework::Configurable cfgCentMin{"cfgCentMin", -1, "min. centrality"}; o2::framework::Configurable cfgCentMax{"cfgCentMax", 999, "max. centrality"}; - o2::framework::Configurable maxY_rec{"maxY_rec", 0.9, "maximum rapidity for reconstructed particles"}; - o2::framework::Configurable fd_k0s_to_pi0{"fd_k0s_pi0", "1.0", "feed down correction to pi0"}; + o2::framework::Configurable maxYRec{"maxYRec", 0.9, "maximum rapidity for reconstructed particles"}; + o2::framework::Configurable fdK0sToPi0{"fdK0sPi0", "1.0", "feed down correction to pi0"}; o2::framework::Configurable cfgRequireTrueAssociation{"cfgRequireTrueAssociation", false, "flag to require true mc collision association"}; o2::framework::Configurable 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 cfgAlphaMeson{"cfgAlphaMeson", 0.65, "photon energy asymmetry distribution parameter for specific value cut"}; o2::framework::Configurable cfgAlphaMesonA{"cfgAlphaMesonA", 0.65, "photon energy asymmetry distribution parameter A for pT dependent cut (A * tanh(B*pT))"}; o2::framework::Configurable cfgAlphaMesonB{"cfgAlphaMesonB", 1.2, "photon energy asymmetry distribution parameter B for pT dependent cut (A * tanh(B*pT))"}; + o2::framework::Configurable cfgGGContaCheck{"cfgGGContaCheck", false, "check gamma gamma contamination of dalitz"}; EMPhotonEventCut fEMEventCut; struct : o2::framework::ConfigurableGroup { @@ -224,7 +226,7 @@ struct Pi0EtaToGammaGammaMC { o2::framework::Configurable 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(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, o2::soa::Filtered, 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, o2::soa::Filtered, 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; }