Skip to content

[PWGHF] Fix matching to MC for B2JPsi analyses to account for photons produced by EVTGEN #18668

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
fgrosa merged 3 commits into from
Nov 10, 2026
Merged

[PWGHF] Fix matching to MC for B2JPsi analyses to account for photons produced by EVTGEN #18668

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
e0b5305
Fix matching to MC for B2JPsi analyses to account for photons produce…
fgrosa Nov 9, 2026
e4ab123
Remove debug logs
fgrosa Nov 9, 2026
42316fd
Dummy commit
fgrosa Nov 10, 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
97 changes: 58 additions & 39 deletions PWGHF/D2H/TableProducer/dataCreatorJpsiHadReduced.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -535,11 +535,11 @@ struct HfDataCreatorJpsiHadReduced {

// Now we check that subtracting gammas from bremsstrahlung, the number of daughters of the J/Psi is 2 (µ+µ-)
// In EVTGEN, bremmstralhung photons are considered as daughters of the J/Psi
std::vector<int> jpsiDauIndices{};
std::vector<int> jpsiDauPdgs{};
for (const auto& dau : jpsiPart.template daughters_as<PParticles>()) {
jpsiDauIndices.push_back(dau.pdgCode());
jpsiDauPdgs.push_back(dau.pdgCode());
}
if (jpsiDauIndices.size() - std::count(jpsiDauIndices.begin(), jpsiDauIndices.end(), kGamma) == nProngsJpsi) {
if (jpsiDauPdgs.size() - std::count(jpsiDauPdgs.begin(), jpsiDauPdgs.end(), kGamma) == nProngsJpsi) {
flag = signB * o2::hf_decay::hf_cand_beauty::BplusToJpsiKToMuMuK;
}
}
Expand Down Expand Up @@ -568,11 +568,11 @@ struct HfDataCreatorJpsiHadReduced {

// Now we check that subtracting gammas from bremsstrahlung, the number of daughters of the J/Psi is 2 (e+e-)
// In EVTGEN, bremmstralhung photons are considered as daughters of the J/Psi
std::vector<int> jpsiDauIndices{};
std::vector<int> jpsiDauPdgs{};
for (const auto& dau : jpsiPart.template daughters_as<PParticles>()) {
jpsiDauIndices.push_back(dau.pdgCode());
jpsiDauPdgs.push_back(dau.pdgCode());
}
if (jpsiDauIndices.size() - std::count(jpsiDauIndices.begin(), jpsiDauIndices.end(), kGamma) == nProngsJpsi) {
if (jpsiDauPdgs.size() - std::count(jpsiDauPdgs.begin(), jpsiDauPdgs.end(), kGamma) == nProngsJpsi) {
flag = signB * o2::hf_decay::hf_cand_beauty::BplusToJpsiKToEEK;
}
}
Expand Down Expand Up @@ -615,16 +615,16 @@ struct HfDataCreatorJpsiHadReduced {
std::vector<int> dauIndices{};
const int8_t signB = sign;
// Since J/Psi is self-conjugate, we need to set acceptAntiParticles to false
if (RecoDecay::isMatchedMCGen<false>(particlesMc, bParticle, signB * Pdg::kB0, std::array{+kJPsi, signB * Pdg::kK0Star892}, false, nullptr, 1, &dauIndices)) {
if (RecoDecay::isMatchedMCGen(particlesMc, bParticle, signB * Pdg::kB0, std::array{+kJPsi, signB * Pdg::kK0Star892}, false, nullptr, 1, &dauIndices)) {
const auto& jpsiPart = particlesMc.rawIteratorAt(dauIndices[0]);

// Now we check that subtracting gammas from bremsstrahlung, the number of daughters of the J/Psi is 2 (e+e-)
// In EVTGEN, bremmstralhung photons are considered as daughters of the J/Psi
std::vector<int> jpsiDauIndices{};
std::vector<int> jpsiDauPdgs{};
for (const auto& dau : jpsiPart.template daughters_as<PParticles>()) {
jpsiDauIndices.push_back(dau.pdgCode());
jpsiDauPdgs.push_back(dau.pdgCode());
}
if (jpsiDauIndices.size() - std::count(jpsiDauIndices.begin(), jpsiDauIndices.end(), kGamma) == nProngsJpsi) {
if (jpsiDauPdgs.size() - std::count(jpsiDauPdgs.begin(), jpsiDauPdgs.end(), kGamma) == nProngsJpsi) {
flag = signB * o2::hf_decay::hf_cand_beauty::B0ToJpsiPiKToMuMuPiK;
indexRec = RecoDecay::getMatchedMCRec<false, false, false, true, true>(particlesMc, std::array{vecDaughtersB[2], vecDaughtersB[3]}, Pdg::kK0Star892, std::array{+kKPlus, -kPiPlus}, true, &sign, 1);
if (indexRec > -1) {
Expand Down Expand Up @@ -652,16 +652,16 @@ struct HfDataCreatorJpsiHadReduced {
const int8_t signB = sign;

// Since J/Psi is self-conjugate, we need to set acceptAntiParticles to false
if (RecoDecay::isMatchedMCGen<false>(particlesMc, bParticle, signB * Pdg::kB0, std::array{+kJPsi, signB * Pdg::kK0Star892}, false, nullptr, 1, &dauIndices)) {
if (RecoDecay::isMatchedMCGen(particlesMc, bParticle, signB * Pdg::kB0, std::array{+kJPsi, signB * Pdg::kK0Star892}, false, nullptr, 1, &dauIndices)) {
const auto& jpsiPart = particlesMc.rawIteratorAt(dauIndices[0]);

// Now we check that subtracting gammas from bremsstrahlung, the number of daughters of the J/Psi is 2 (e+e-)
// In EVTGEN, bremmstralhung photons are considered as daughters of the J/Psi
std::vector<int> jpsiDauIndices{};
std::vector<int> jpsiDauPdgs{};
for (const auto& dau : jpsiPart.template daughters_as<PParticles>()) {
jpsiDauIndices.push_back(dau.pdgCode());
jpsiDauPdgs.push_back(dau.pdgCode());
}
if (jpsiDauIndices.size() - std::count(jpsiDauIndices.begin(), jpsiDauIndices.end(), kGamma) == nProngsJpsi) {
if (jpsiDauPdgs.size() - std::count(jpsiDauPdgs.begin(), jpsiDauPdgs.end(), kGamma) == nProngsJpsi) {
flag = signB * o2::hf_decay::hf_cand_beauty::B0ToJpsiPiKToEEPiK;
indexRec = RecoDecay::getMatchedMCRec<false, false, false, true, true>(particlesMc, std::array{vecDaughtersB[2], vecDaughtersB[3]}, Pdg::kK0Star892, std::array{+kKPlus, -kPiPlus}, true, &sign, 1);
if (indexRec > -1) {
Expand Down Expand Up @@ -715,16 +715,16 @@ struct HfDataCreatorJpsiHadReduced {
std::vector<int> dauIndices{};
const int8_t signB = sign;
// Since J/Psi and phi are self-conjugate, we need to set acceptAntiParticles to false
if (RecoDecay::isMatchedMCGen<false>(particlesMc, bParticle, signB * Pdg::kBS, std::array{+kJPsi, +Pdg::kPhi}, false, nullptr, 1, &dauIndices)) {
if (RecoDecay::isMatchedMCGen(particlesMc, bParticle, signB * Pdg::kBS, std::array{+kJPsi, +Pdg::kPhi}, false, nullptr, 1, &dauIndices)) {
const auto& jpsiPart = particlesMc.rawIteratorAt(dauIndices[0]);

// Now we check that subtracting gammas from bremsstrahlung, the number of daughters of the J/Psi is 2 (e+e-)
// In EVTGEN, bremmstralhung photons are considered as daughters of the J/Psi
std::vector<int> jpsiDauIndices{};
std::vector<int> jpsiDauPdgs{};
for (const auto& dau : jpsiPart.template daughters_as<PParticles>()) {
jpsiDauIndices.push_back(dau.pdgCode());
jpsiDauPdgs.push_back(dau.pdgCode());
}
if (jpsiDauIndices.size() - std::count(jpsiDauIndices.begin(), jpsiDauIndices.end(), kGamma) == nProngsJpsi) {
if (jpsiDauPdgs.size() - std::count(jpsiDauPdgs.begin(), jpsiDauPdgs.end(), kGamma) == nProngsJpsi) {
flag = signB * o2::hf_decay::hf_cand_beauty::BsToJpsiKKToMuMuKK;
indexRec = RecoDecay::getMatchedMCRec<false, false, false, true, true>(particlesMc, std::array{vecDaughtersB[2], vecDaughtersB[3]}, Pdg::kPhi, std::array{+kKPlus, -kKPlus}, true, &sign, 1);
if (indexRec > -1) {
Expand Down Expand Up @@ -752,16 +752,16 @@ struct HfDataCreatorJpsiHadReduced {
const int8_t signB = sign;

// Since J/Psi and phi are self-conjugate, we need to set acceptAntiParticles to false
if (RecoDecay::isMatchedMCGen<false>(particlesMc, bParticle, signB * Pdg::kBS, std::array{+kJPsi, +Pdg::kPhi}, false, nullptr, 1, &dauIndices)) {
if (RecoDecay::isMatchedMCGen(particlesMc, bParticle, signB * Pdg::kBS, std::array{+kJPsi, +Pdg::kPhi}, false, nullptr, 1, &dauIndices)) {
const auto& jpsiPart = particlesMc.rawIteratorAt(dauIndices[0]);

// Now we check that subtracting gammas from bremsstrahlung, the number of daughters of the J/Psi is 2 (e+e-)
// In EVTGEN, bremmstralhung photons are considered as daughters of the J/Psi
std::vector<int> jpsiDauIndices{};
std::vector<int> jpsiDauPdgs{};
for (const auto& dau : jpsiPart.template daughters_as<PParticles>()) {
jpsiDauIndices.push_back(dau.pdgCode());
jpsiDauPdgs.push_back(dau.pdgCode());
}
if (jpsiDauIndices.size() - std::count(jpsiDauIndices.begin(), jpsiDauIndices.end(), kGamma) == nProngsJpsi) {
if (jpsiDauPdgs.size() - std::count(jpsiDauPdgs.begin(), jpsiDauPdgs.end(), kGamma) == nProngsJpsi) {
flag = signB * o2::hf_decay::hf_cand_beauty::BsToJpsiKKToEEKK;
indexRec = RecoDecay::getMatchedMCRec<false, false, false, true, true>(particlesMc, std::array{vecDaughtersB[2], vecDaughtersB[3]}, Pdg::kPhi, std::array{+kKPlus, -kKPlus}, true, &sign, 1);
if (indexRec > -1) {
Expand Down Expand Up @@ -844,14 +844,21 @@ struct HfDataCreatorJpsiHadReduced {

// Match generated particles.
for (const auto& particle : mcParticlesPerMcColl) {
int8_t sign{0}, flag{0}, channel{0};
int8_t sign = particle.pdgCode() > 0 ? 1 : -1;
int8_t flag{0}, channel{0};
if constexpr (DecChannel == DecayChannel::BplusToJpsiK) {
// B+ → J/Psi K+ → (µ+µ-) K+
if (RecoDecay::isMatchedMCGen<false>(particlesMc, particle, Pdg::kBPlus, std::array{static_cast<int>(Pdg::kJPsi), +kKPlus}, true, &sign)) {
// Since J/Psi is self-conjugate, we need to set acceptAntiParticles to false
if (RecoDecay::isMatchedMCGen(particlesMc, particle, sign * Pdg::kBPlus, std::array{static_cast<int>(Pdg::kJPsi), sign * kKPlus}, false)) {
// Match J/Psi -> µ+µ-
auto candJpsiMC = particlesMc.rawIteratorAt(particle.daughtersIds().front());
// Printf("Checking J/Psi -> µ+µ-");
if (RecoDecay::isMatchedMCGen(particlesMc, candJpsiMC, static_cast<int>(Pdg::kJPsi), std::array{-kMuonMinus, +kMuonMinus}, true)) {
const auto& jpsiPart = particlesMc.rawIteratorAt(particle.daughtersIds().front());
// Now we check that subtracting gammas from bremsstrahlung, the number of daughters of the J/Psi is 2 (µ+µ-)
// In EVTGEN, bremmstralhung photons are considered as daughters of the J/Psi
std::vector<int> jpsiDauPdgs{};
for (const auto& dau : jpsiPart.template daughters_as<aod::McParticles>()) {
jpsiDauPdgs.push_back(dau.pdgCode());
}
if (std::count(jpsiDauPdgs.begin(), jpsiDauPdgs.end(), kMuonMinus) == 1 && std::count(jpsiDauPdgs.begin(), jpsiDauPdgs.end(), -kMuonMinus) && jpsiDauPdgs.size() - std::count(jpsiDauPdgs.begin(), jpsiDauPdgs.end(), kGamma) == nProngsJpsi) {
flag = sign * o2::hf_decay::hf_cand_beauty::BplusToJpsiKToMuMuK;
}
}
Expand Down Expand Up @@ -880,16 +887,22 @@ struct HfDataCreatorJpsiHadReduced {
ptProngs[1], yProngs[1], etaProngs[1], hfRejMap, centFT0C, centFT0M);
} else if constexpr (DecChannel == DecayChannel::B0ToJpsiK0Star) {
// B0 → J/Psi K*0 → (µ+µ-) (K+pi-)
if (RecoDecay::isMatchedMCGen<true>(particlesMc, particle, Pdg::kB0, std::array{static_cast<int>(Pdg::kJPsi), +kKPlus, -kPiPlus}, true, &sign, 2)) {
// Since J/Psi is self-conjugate, we need to set acceptAntiParticles to false
if (RecoDecay::isMatchedMCGen(particlesMc, particle, sign * Pdg::kB0, std::array{static_cast<int>(Pdg::kJPsi), sign * kKPlus, -1 * sign * kPiPlus}, false, nullptr, 2)) {
// Match J/Psi -> µ+µ- and K*0 -> K+pi-
auto candJpsiMC = particlesMc.rawIteratorAt(particle.daughtersIds().front());
auto candK0StarMC = particlesMc.rawIteratorAt(particle.daughtersIds().back());
// Printf("Checking J/Psi -> µ+µ- and K*0 -> K+pi-");
if (RecoDecay::isMatchedMCGen(particlesMc, candJpsiMC, static_cast<int>(Pdg::kJPsi), std::array{-kMuonMinus, +kMuonMinus}, true)) {
const auto& jpsiPart = particlesMc.rawIteratorAt(particle.daughtersIds().front());
// Now we check that subtracting gammas from bremsstrahlung, the number of daughters of the J/Psi is 2 (µ+µ-)
// In EVTGEN bremmstralhung photons are considered as daughters of the J/Psi
std::vector<int> jpsiDauPdgs{};
for (const auto& dau : jpsiPart.template daughters_as<aod::McParticles>()) {
jpsiDauPdgs.push_back(dau.pdgCode());
}
if (std::count(jpsiDauPdgs.begin(), jpsiDauPdgs.end(), kMuonMinus) == 1 && std::count(jpsiDauPdgs.begin(), jpsiDauPdgs.end(), -kMuonMinus) && jpsiDauPdgs.size() - std::count(jpsiDauPdgs.begin(), jpsiDauPdgs.end(), kGamma) == nProngsJpsi) {
flag = sign * o2::hf_decay::hf_cand_beauty::B0ToJpsiPiKToMuMuPiK;
}
const auto& kStarPart = particlesMc.rawIteratorAt(particle.daughtersIds().back());
// Check K*0 -> K+pi-
if (RecoDecay::isMatchedMCGen(particlesMc, candK0StarMC, static_cast<int>(Pdg::kK0Star892), std::array{+kKPlus, -kPiPlus}, true)) {
if (RecoDecay::isMatchedMCGen(particlesMc, kStarPart, static_cast<int>(Pdg::kK0Star892), std::array{+kKPlus, -kPiPlus}, true)) {
channel = o2::hf_decay::hf_cand_beauty::B0ToJpsiKstar0;
}
}
Expand Down Expand Up @@ -918,16 +931,22 @@ struct HfDataCreatorJpsiHadReduced {
ptProngs[1], yProngs[1], etaProngs[1], hfRejMap, centFT0C, centFT0M);
} else if constexpr (DecChannel == DecayChannel::BsToJpsiPhi) {
// Bs → J/Psi phi → (µ+µ-) (K+K-)
if (RecoDecay::isMatchedMCGen<true>(particlesMc, particle, Pdg::kBS, std::array{static_cast<int>(Pdg::kJPsi), +kKPlus, -kKPlus}, true, &sign, 2)) {
// Since J/Psi is self-conjugate, we need to set acceptAntiParticles to false
if (RecoDecay::isMatchedMCGen(particlesMc, particle, sign * Pdg::kBS, std::array{static_cast<int>(Pdg::kJPsi), sign * kKPlus, -1 * sign * kKPlus}, false, nullptr, 2)) {
// Match J/Psi -> µ+µ- and phi -> K+K-
auto candJpsiMC = particlesMc.rawIteratorAt(particle.daughtersIds().front());
auto candPhiMC = particlesMc.rawIteratorAt(particle.daughtersIds().back());
// Printf("Checking J/Psi -> µ+µ- and phi -> K+K-");
if (RecoDecay::isMatchedMCGen(particlesMc, candJpsiMC, static_cast<int>(Pdg::kJPsi), std::array{-kMuonMinus, +kMuonMinus}, true)) {
const auto& jpsiPart = particlesMc.rawIteratorAt(particle.daughtersIds().front());
// Now we check that subtracting gammas from bremsstrahlung, the number of daughters of the J/Psi is 2 (µ+µ-)
// In EVTGEN, bremmstralhung photons are considered as daughters of the J/Psi
std::vector<int> jpsiDauPdgs{};
for (const auto& dau : jpsiPart.template daughters_as<aod::McParticles>()) {
jpsiDauPdgs.push_back(dau.pdgCode());
}
if (std::count(jpsiDauPdgs.begin(), jpsiDauPdgs.end(), kMuonMinus) == 1 && std::count(jpsiDauPdgs.begin(), jpsiDauPdgs.end(), -kMuonMinus) && jpsiDauPdgs.size() - std::count(jpsiDauPdgs.begin(), jpsiDauPdgs.end(), kGamma) == nProngsJpsi) {
flag = sign * o2::hf_decay::hf_cand_beauty::BsToJpsiKKToMuMuKK;
}
const auto& phiPart = particlesMc.rawIteratorAt(particle.daughtersIds().back());
// Check phi -> K+K-
if (RecoDecay::isMatchedMCGen(particlesMc, candPhiMC, static_cast<int>(Pdg::kPhi), std::array{-kKPlus, +kKPlus}, true)) {
if (RecoDecay::isMatchedMCGen(particlesMc, phiPart, static_cast<int>(Pdg::kPhi), std::array{-kKPlus, +kKPlus}, true)) {
channel = o2::hf_decay::hf_cand_beauty::BsToJpsiPhi;
}
}
Expand Down
8 changes: 5 additions & 3 deletions PWGHF/D2H/Tasks/taskB0ToJpsiK0StarReduced.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -342,7 +342,7 @@ struct HfTaskB0ToJpsiK0StarReduced {
const AxisSpec axisImpactPar{200, -0.05, 0.05};
const AxisSpec axisPtJpsi{100, 0., 50.};
const AxisSpec axisRapidity{100, -2., 2.};
const AxisSpec axisPtB{std::vector<double>(binsPt), "#it{p}_{T}^{B^{0}} (GeV/#it{c})"};
const AxisSpec axisPtB{500, 0., 50., "#it{p}_{T}^{B^{0}} (GeV/#it{c})"};
const AxisSpec axisPtK0Star{100, 0.f, 10.f};

registry.add("hMass", b0CandTitle + "inv. mass J/#Psi K*^{0} (GeV/#it{c}^{2});" + stringPt, {HistType::kTH2F, {axisMassB0, axisPtB}});
Expand All @@ -353,7 +353,8 @@ struct HfTaskB0ToJpsiK0StarReduced {
if (doprocessMc || doprocessMcWithB0Ml) {
registry.add("hPtJpsiGen", mcParticleMatched + "J/#Psi #it{p}_{T}^{gen} (GeV/#it{c}); B^{0} " + stringPt, {HistType::kTH2F, {axisPtProng, axisPtB}});
registry.add("hPtK0StarGen", mcParticleMatched + "K*^{0} #it{p}_{T}^{gen} (GeV/#it{c}); B^{0} " + stringPt, {HistType::kTH2F, {axisPtK0Star, axisPtB}});
registry.add("hYGenWithProngsInAcceptance", mcParticleMatched + "B^{0} #it{p}_{T}^{gen} (GeV/#it{c}); B^{0} #it{y}", {HistType::kTH2F, {axisPtProng, axisRapidity}});
registry.add("hPtYGenSig", mcParticleMatched + "B^{0} #it{p}_{T}^{gen} (GeV/#it{c}); B^{0} #it{y}", {HistType::kTH2F, {axisPtB, axisRapidity}});
registry.add("hPtYWithProngsInAccepanceGenSig", mcParticleMatched + "B^{0} #it{p}_{T}^{gen} (GeV/#it{c}); B^{0} #it{y}", {HistType::kTH2F, {axisPtB, axisRapidity}});
registry.add("hMassRecSig", b0CandMatch + "inv. mass J/#Psi K*^{0} (GeV/#it{c}^{2}); B^{0} " + stringPt, {HistType::kTH2F, {axisMassB0, axisPtB}});
registry.add("hMassJpsiRecSig", b0CandMatch + "inv. mass #mu^{+}#mu^{#minus} (GeV/#it{c}^{2}); J/#Psi " + stringPt, {HistType::kTH2F, {axisMassJpsi, axisPtJpsi}});
registry.add("hMassK0StarRecSig", b0CandMatch + "inv. mass K^{+}#pi^{#minus} (GeV/#it{c}^{2}); K*^{0} " + stringPt, {HistType::kTH2F, {axisMassK0Star, axisPtK0Star}});
Expand Down Expand Up @@ -679,9 +680,10 @@ struct HfTaskB0ToJpsiK0StarReduced {
registry.fill(HIST("hPtJpsiGen"), ptProngs[0], ptParticle);
registry.fill(HIST("hPtK0StarGen"), ptProngs[1], ptParticle);

registry.fill(HIST("hPtYGenSig"), ptParticle, yParticle);
// generated B0 with daughters in geometrical acceptance
if (prongsInAcc) {
registry.fill(HIST("hYGenWithProngsInAcceptance"), ptParticle, yParticle);
registry.fill(HIST("hPtYWithProngsInAccepanceGenSig"), ptParticle, yParticle);
}
}

Expand Down
Loading
Loading