Skip to content

[PWGEM] Differentiate between b2c2e and b2cc2e, add function hasNCharmHadrons… #14941

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
feisenhu merged 4 commits into from
Feb 13, 2026
+91 −11
Merged

[PWGEM] Differentiate between b2c2e and b2cc2e, add function hasNCharmHadrons… #14941

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
5ac6cfd
Differentiate between b2c2e and b2cc2e, add function hasNCharmHadrons…
feisenhu Feb 12, 2026
2a3aaca
Apply changes for b2cc2l for CheckMCPairTemplates
feisenhu Feb 12, 2026
6fd19f4
Merge pull request #44 from feisenhu/b2CC2ee_Signal_CheckMCPairTemplates
feisenhu Feb 12, 2026
8a448e4
add if constexpr doMoreDifferentially
dsekihat Feb 12, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
30 changes: 21 additions & 9 deletions PWGEM/Dilepton/Tasks/checkMCPairTemplate.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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/"};

Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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++) {
Expand Down Expand Up @@ -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<true>(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);
Expand Down Expand Up @@ -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<true>(t1mc, t2mc, mcparticles);
if (mother_id < 0 && hfee_type < 0) {
return false;
}
Expand Down Expand Up @@ -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<true>(t1, t2, mcparticles);
if (mother_id < 0 && hfee_type < 0) {
return false;
}
Expand Down Expand Up @@ -2326,6 +2332,12 @@ struct checkMCPairTemplate {
case static_cast<int>(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<int>(EM_HFeeType::kBCCe_BCe):
fillGenHistograms<20>(sign1, sign2, mp1.pdgCode(), mp2.pdgCode(), v12.M(), v12.Pt(), weight); // b2cc2l_b2c2l
break;
case static_cast<int>(EM_HFeeType::kBCCe_BCCe):
fillGenHistograms<21>(sign1, sign2, mp1.pdgCode(), mp2.pdgCode(), v12.M(), v12.Pt(), weight); // b2cc2l_b2cc2l
break;
default:
break;
}
Expand Down Expand Up @@ -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<true>(t1mc, t2mc, mcparticles);
if (mother_id < 0 && hfee_type < 0) {
return false;
}
Expand Down
72 changes: 70 additions & 2 deletions PWGEM/Dilepton/Utils/MCUtilities.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
};

//_______________________________________________________________________
Expand Down Expand Up @@ -461,6 +463,56 @@ int IsFromCharm(TMCParticle const& p, TMCParticles const& mcparticles)
return -999;
}
//_______________________________________________________________________
template <typename T, typename U>
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 <typename TMCParticle>
bool isFlavorOscillationB(TMCParticle const& mcParticle)
{
Expand Down Expand Up @@ -551,7 +603,7 @@ int find1stHadron(TMCParticle const& mcParticle, TMCParticles const& mcParticles
return hadronId;
}
//_______________________________________________________________________
template <typename TMCParticle1, typename TMCParticle2, typename TMCParticles>
template <bool doMoreDifferentially = false, typename TMCParticle1, typename TMCParticle2, typename TMCParticles>
int IsHF(TMCParticle1 const& p1, TMCParticle2 const& p2, TMCParticles const& mcparticles)
{
if (!p1.has_mothers() || !p2.has_mothers()) {
Expand Down Expand Up @@ -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<int>(EM_HFeeType::kBCe_BCe); // b->c->e and b->c->e, decay type = 1

if constexpr (!doMoreDifferentially) {
return static_cast<int>(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<int>(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<int>(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<int>(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<int>(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)) {
Expand Down
Loading