diff --git a/PWGEM/Dilepton/Tasks/checkMCPairTemplate.cxx b/PWGEM/Dilepton/Tasks/checkMCPairTemplate.cxx index 436a0aa5777..d13a75eb9b8 100644 --- a/PWGEM/Dilepton/Tasks/checkMCPairTemplate.cxx +++ b/PWGEM/Dilepton/Tasks/checkMCPairTemplate.cxx @@ -309,7 +309,7 @@ struct checkMCPairTemplate { HistogramRegistry fRegistry{"output", {}, OutputObjHandlingPolicy::AnalysisObject, false, false}; static constexpr std::string_view event_cut_types[2] = {"before/", "after/"}; static constexpr std::string_view pair_sign_types[3] = {"uls/", "lspp/", "lsmm/"}; - static constexpr std::string_view dilepton_source_types[20] = { + static constexpr std::string_view dilepton_source_types[22] = { "sm/Photon/", // 0 "sm/PromptPi0/", // 1 "sm/NonPromptPi0/", // 2 @@ -329,7 +329,9 @@ struct checkMCPairTemplate { "bbbar/b2l_b2l/", // 16 "bbbar/b2c2l_b2c2l/", // 17 "bbbar/b2c2l_b2l_sameb/", // 18 - "bbbar/b2c2l_b2l_diffb/" // 19 + "bbbar/b2c2l_b2l_diffb/", // 19 + "bbbar/b2cc2l_b2c2l/", // 20 + "bbbar/b2cc2l_b2cc2l/", // 21 }; // unordered_map is better, but cannot be constexpr. static constexpr std::string_view unfolding_dilepton_source_types[3] = {"sm/", "ccbar/", "bbbar/"}; @@ -424,6 +426,8 @@ struct checkMCPairTemplate { fRegistry.addClone("Generated/ccbar/c2l_c2l/", "Generated/bbbar/b2c2l_b2c2l/"); fRegistry.addClone("Generated/ccbar/c2l_c2l/", "Generated/bbbar/b2c2l_b2l_sameb/"); fRegistry.addClone("Generated/ccbar/c2l_c2l/", "Generated/bbbar/b2c2l_b2l_diffb/"); // LS + fRegistry.addClone("Generated/ccbar/c2l_c2l/", "Generated/bbbar/b2cc2l_b2c2l/"); + fRegistry.addClone("Generated/ccbar/c2l_c2l/", "Generated/bbbar/b2cc2l_b2cc2l/"); // for charmed hadrons // create 28 combinations static constexpr std::string_view charmed_mesons[] = {"Dplus", "D0", "Dsplus"}; // 411, 421, 431 @@ -562,6 +566,8 @@ struct checkMCPairTemplate { fRegistry.addClone("Pair/ccbar/c2l_c2l/", "Pair/bbbar/b2c2l_b2c2l/"); fRegistry.addClone("Pair/ccbar/c2l_c2l/", "Pair/bbbar/b2c2l_b2l_sameb/"); fRegistry.addClone("Pair/ccbar/c2l_c2l/", "Pair/bbbar/b2c2l_b2l_diffb/"); // LS + fRegistry.addClone("Pair/ccbar/c2l_c2l/", "Pair/bbbar/b2cc2l_b2c2l/"); + fRegistry.addClone("Pair/ccbar/c2l_c2l/", "Pair/bbbar/b2cc2l_b2cc2l/"); if (cfgFillSeparateCharmHadronPairs) { for (int im = 0; im < nm_c; im++) { @@ -1949,10 +1955,10 @@ struct checkMCPairTemplate { // o2::math_utils::bringToPMPi(phiPol); // float quadmom = (3.f * std::pow(cos_thetaPol, 2) - 1.f) / 2.f; - if ((FindCommonMotherFrom2ProngsWithoutPDG(t1mc, t2mc) > 0 || IsHF(t1mc, t2mc, mcparticles) > 0) && is_pair_from_same_mcevent) { // for bkg study - if (std::abs(t1mc.pdgCode()) != pdg_lepton || std::abs(t2mc.pdgCode()) != pdg_lepton) { // hh or lh correlated bkg - if (std::abs(t1mc.pdgCode()) != pdg_lepton && std::abs(t2mc.pdgCode()) != pdg_lepton) { // hh correlated bkg - if (t1.sign() * t2.sign() < 0) { // ULS + if ((FindCommonMotherFrom2ProngsWithoutPDG(t1mc, t2mc) > 0 || IsHF(t1mc, t2mc, mcparticles) > 0) && is_pair_from_same_mcevent) { // for bkg study + if (std::abs(t1mc.pdgCode()) != pdg_lepton || std::abs(t2mc.pdgCode()) != pdg_lepton) { // hh or lh correlated bkg + if (std::abs(t1mc.pdgCode()) != pdg_lepton && std::abs(t2mc.pdgCode()) != pdg_lepton) { // hh correlated bkg + if (t1.sign() * t2.sign() < 0) { // ULS fRegistry.fill(HIST("Pair/corr_bkg_hh/uls/hs"), v12.M(), v12.Pt(), pair_dca, weight); } else if (t1.sign() > 0 && t2.sign() > 0) { // LS++ fRegistry.fill(HIST("Pair/corr_bkg_hh/lspp/hs"), v12.M(), v12.Pt(), pair_dca, weight); @@ -1990,7 +1996,7 @@ struct checkMCPairTemplate { return false; } int mother_id = std::max({FindSMULS(t1mc, t2mc, mcparticles), FindSMULS(t2mc, t1mc, mcparticles), FindSMLSPP(t1mc, t2mc, mcparticles), FindSMLSMM(t1mc, t2mc, mcparticles)}); - int hfee_type = IsHF(t1mc, t2mc, mcparticles); + int hfee_type = IsHF(t1mc, t2mc, mcparticles); if (mother_id < 0 && hfee_type < 0) { return false; } @@ -2159,7 +2165,7 @@ struct checkMCPairTemplate { } int mother_id = std::max({FindSMULS(t1, t2, mcparticles), FindSMULS(t2, t1, mcparticles), FindSMLSPP(t1, t2, mcparticles), FindSMLSMM(t1, t2, mcparticles)}); - int hfee_type = IsHF(t1, t2, mcparticles); + int hfee_type = IsHF(t1, t2, mcparticles); if (mother_id < 0 && hfee_type < 0) { return false; } @@ -2326,6 +2332,12 @@ struct checkMCPairTemplate { case static_cast(EM_HFeeType::kBCe_Be_DiffB): fillGenHistograms<19>(sign1, sign2, mp1.pdgCode(), mp2.pdgCode(), v12.M(), v12.Pt(), weight); // b2c2l_b2l_diffb break; + case static_cast(EM_HFeeType::kBCCe_BCe): + fillGenHistograms<20>(sign1, sign2, mp1.pdgCode(), mp2.pdgCode(), v12.M(), v12.Pt(), weight); // b2cc2l_b2c2l + break; + case static_cast(EM_HFeeType::kBCCe_BCCe): + fillGenHistograms<21>(sign1, sign2, mp1.pdgCode(), mp2.pdgCode(), v12.M(), v12.Pt(), weight); // b2cc2l_b2cc2l + break; default: break; } @@ -2791,7 +2803,7 @@ struct checkMCPairTemplate { return false; } int mother_id = std::max({FindSMULS(t1mc, t2mc, mcparticles), FindSMULS(t2mc, t1mc, mcparticles), FindSMLSPP(t1mc, t2mc, mcparticles), FindSMLSMM(t1mc, t2mc, mcparticles)}); - int hfee_type = IsHF(t1mc, t2mc, mcparticles); + int hfee_type = IsHF(t1mc, t2mc, mcparticles); if (mother_id < 0 && hfee_type < 0) { return false; } diff --git a/PWGEM/Dilepton/Utils/MCUtilities.h b/PWGEM/Dilepton/Utils/MCUtilities.h index 1026576a805..ae4f7682d7a 100644 --- a/PWGEM/Dilepton/Utils/MCUtilities.h +++ b/PWGEM/Dilepton/Utils/MCUtilities.h @@ -32,6 +32,8 @@ enum class EM_HFeeType : int { kBCe_BCe = 2, // ULS kBCe_Be_SameB = 3, // ULS kBCe_Be_DiffB = 4, // LS + kBCCe_BCe = 5, // ULS + kBCCe_BCCe = 6, // ULS }; //_______________________________________________________________________ @@ -461,6 +463,56 @@ int IsFromCharm(TMCParticle const& p, TMCParticles const& mcparticles) return -999; } //_______________________________________________________________________ +template +int hasNCharmHadronsInBeautyDecay(T const& mcParticle, U const& mcParticles) +{ + // require that the direct mother is beauty hadron via semileptonice decay. e.g. hb->e, not hb->X->pi0->eegamma + if (!mcParticle.has_mothers()) { + return -999; + } + if (!IsFromBeauty(mcParticle, mcParticles)) { + return -999; + } + + auto mp = mcParticles.iteratorAt(mcParticle.mothersIds()[0]); + if (!isCharmMeson(mp) && !isCharmBaryon(mp)) { + return -999; + } + + int motherid = mcParticle.mothersIds()[0]; // first mother index + while (motherid > -1) { + if (motherid < mcParticles.size()) { // protect against bad mother indices. why is this needed? + mp = mcParticles.iteratorAt(motherid); + if (std::abs(mp.pdgCode()) < 1e+9 && (std::to_string(std::abs(mp.pdgCode()))[std::to_string(std::abs(mp.pdgCode())).length() - 3] == '5' || std::to_string(std::abs(mp.pdgCode()))[std::to_string(std::abs(mp.pdgCode())).length() - 4] == '5')) { + // check if mp has two charm hadrons as daughters + if (mp.has_daughters()) { + const auto& daughtersIds = mp.daughtersIds(); + int count_charm_hadron = 0; + for (const auto& daughterId : daughtersIds) { + if (daughterId >= 0 && daughterId < mcParticles.size()) { + auto daughter = mcParticles.iteratorAt(daughterId); + if (isCharmMeson(daughter) || isCharmBaryon(daughter)) { + count_charm_hadron++; + } + } + } + return count_charm_hadron; + } else { + LOGF(debug, "Something went wrong: Did not find any daughter for the current mother! Can't be a mother if there are no daughters\n"); + } + } + if (mp.has_mothers()) { + motherid = mp.mothersIds()[0]; + } else { + return -999; + } + } else { + LOGF(info, "Mother label(%d) exceeds the McParticles size(%d)", motherid, mcParticles.size()); + } + } + return -999; +} +//_______________________________________________________________________ template bool isFlavorOscillationB(TMCParticle const& mcParticle) { @@ -551,7 +603,7 @@ int find1stHadron(TMCParticle const& mcParticle, TMCParticles const& mcParticles return hadronId; } //_______________________________________________________________________ -template +template int IsHF(TMCParticle1 const& p1, TMCParticle2 const& p2, TMCParticles const& mcparticles) { if (!p1.has_mothers() || !p2.has_mothers()) { @@ -661,7 +713,23 @@ int IsHF(TMCParticle1 const& p1, TMCParticle2 const& p2, TMCParticles const& mcp mothers_pdg1.shrink_to_fit(); mothers_id2.shrink_to_fit(); mothers_pdg2.shrink_to_fit(); - return static_cast(EM_HFeeType::kBCe_BCe); // b->c->e and b->c->e, decay type = 1 + + if constexpr (!doMoreDifferentially) { + return static_cast(EM_HFeeType::kBCe_BCe); // default to b->c->e and b->c->e, decay type = 1 + } else { + int n_c_from_b1 = hasNCharmHadronsInBeautyDecay(p1, mcparticles); + int n_c_from_b2 = hasNCharmHadronsInBeautyDecay(p2, mcparticles); + if (n_c_from_b1 == 1 && n_c_from_b2 == 1) { + return static_cast(EM_HFeeType::kBCe_BCe); // b->c->e and b->c->e, decay type = 1 + } else if (n_c_from_b1 == 2 && n_c_from_b2 == 2) { + return static_cast(EM_HFeeType::kBCCe_BCCe); // b->cc->e and b->cc->e, decay type = 6 + } else if ((n_c_from_b1 == 1 && n_c_from_b2 == 2) || (n_c_from_b1 == 2 && n_c_from_b2 == 1)) { + return static_cast(EM_HFeeType::kBCCe_BCe); // b->cc->e and b->c->e, decay type = 5 + } else { + LOGF(debug, "Unexpected number of charm hadrons from beauty decay: n_c_from_b1 = %d, n_c_from_b2 = %d. Return kBCe_BCe as default.", n_c_from_b1, n_c_from_b2); + return static_cast(EM_HFeeType::kBCe_BCe); // default to b->c->e and b->c->e, decay type = 1 + } + } } if ((is_direct_from_b1 && is_c_from_b2) || (is_direct_from_b2 && is_c_from_b1)) {