Skip to content

Fix test macro for Sc bkg, checking for D*+. #2209

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 2 commits into from
Dec 10, 2025
Merged

Fix test macro for Sc bkg, checking for D*+. #2209

Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
df8e47c
Fix test macro for Sc bkg, checking for D*+.
Dec 9, 2025
76ff153
Fix path.
Dec 9, 2025
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
51 changes: 37 additions & 14 deletions MC/config/PWGHF/ini/tests/GeneratorHF_D2H_ccbar_and_bbbar_gap5_Mode2_corrBkgSigmaC.C
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -5,7 +5,7 @@
#include <iostream>

int External() {
std::string path{"/home/mattia/Documenti/cernbox/Documents/PostDoc/D2H/MC/corrBkgSigmaC/tf1/genevents_Kine.root"};
std::string path{"o2sim_Kine.root"};

int checkPdgQuarkOne{4};
int checkPdgQuarkTwo{5};
Expand All @@ -15,12 +15,13 @@ int External() {
std::array<float, 2> freqRepl = {0.5, 0.5};
std::map<int, int> sumOrigReplacedParticles = {{413, 0}};

std::array<int, 2> checkPdgHadron{14122, 4124};
std::array<int, 3> checkPdgHadron{413, 14122, 4124};
std::map<int, std::vector<std::vector<int>>> checkHadronDecays{ // sorted (!) pdg of daughters
//{14122, {{4222, -211}, {4112, 211}, {4122, 211, -211}}}, // Lc(2595)+
//{4124, {{4222, -211}, {4112, 211}, {4122, 211, -211}}} // Lc(2625)+
{14122, {{-211, 4222}, {211, 4112}, {-211, 211, 4122}}}, // Lc(2595)+
{4124, {{-211, 4222}, {211, 4112}, {-211, 211, 4122}}} // Lc(2625)+
{4124, {{-211, 4222}, {211, 4112}, {-211, 211, 4122}}}, // Lc(2625)+
{413, {{22, 411}, {111, 411}, {211, 421}}} // D*+ as in PYTHIA
};

TFile file(path.c_str(), "READ");
Expand Down Expand Up @@ -64,7 +65,7 @@ int External() {
auto absPdg = std::abs(pdg);
if (std::find(checkPdgHadron.begin(), checkPdgHadron.end(), absPdg) != checkPdgHadron.end()) { // found signal
nSignals++; // count signal PDG
std::cout << "==> signal " << absPdg << " found!" << std::endl;
std::cout << std::endl << "==> signal " << absPdg << " found!" << std::endl;

if (subGeneratorId == checkPdgQuarkOne) { // replacement only for prompt ---> BUT ALSO NON-PROMPT D* SEEM TO BE REPLACED
for (int iRepl{0}; iRepl<2; ++iRepl) {
Expand All @@ -76,8 +77,9 @@ int External() {
sumOrigReplacedParticles[pdgReplParticles[iRepl][0]]++;
}
}
} else if (subGeneratorId == checkPdgQuarkTwo) {
std::cout << " NB: we have a " << absPdg << " also in event with quark " << checkPdgQuarkTwo << std::endl;
}
else if (subGeneratorId == checkPdgQuarkTwo && (absPdg == pdgReplParticles[0][1] || absPdg == pdgReplParticles[1][1])) {
std::cout << " NB: we have a " << absPdg << " also in event with quark " << checkPdgQuarkTwo << ", i.e in an event with c-cbar pair present, tagged with a b-bbar (e.g. double-parton scattering)" << std::endl;
std::cout << " ### mother indices: ";
int idFirstMother = track.getMotherTrackId();
int idSecondMother = track.getSecondMotherTrackId();
Expand All @@ -86,13 +88,29 @@ int External() {
std::cout << i << " ";
motherIds.push_back(i);
}
std::cout << std::endl;

/// To establish if the partonic event is switched on or not, check if the motherIds are all -1 for the current hadron
/// This is ok for Lc excited states deriving from the replacement of a prompt D*+
/// Instead, Lc excited states coming from Lb decays (i.e. not coming from a D*+ replacement))have at least one mother that has idx !=-1 (*)
bool partonicEventOn = false;
if(motherIds != std::vector<int>{-1, -1}) {
bool motherIdsAllMinus1 = true;
for(int idx : motherIds) {
if(idx != -1) {
/// one mother id different from
motherIdsAllMinus1 = false;
break;
}
}
//if(motherIds != std::vector<int>{-1, -1}) {
if(!motherIdsAllMinus1) {
std::cout << "The " << absPdg << " particle has mothers. This should mean that it comes directly from parton hadronization, and that the partonic event was kept in the MC production " << std::endl;
partonicEventOn = true;
}

std::cout << " ### mother PDG codes: ";
std::vector<int> motherPdgCodes = {};
bool updateCounters = true;
if(partonicEventOn) {
for(int i=idFirstMother; i<=idSecondMother; i++) {
motherPdgCodes.push_back(tracks->at(i).GetPdgCode());
Expand All @@ -103,16 +121,21 @@ int External() {
/// This means that the charm hadron comes from the c-quark hadronization, where the c/cbar quark
/// comes from a c-cbar pair present in the current event, tagged with a b-bbar (e.g. double-parton scattering)
if(std::find(motherPdgCodes.begin(), motherPdgCodes.end(), 4) == motherPdgCodes.end() && std::find(motherPdgCodes.begin(), motherPdgCodes.end(), -4) == motherPdgCodes.end()) {
/// if the partinc event is not really saved and we arrive here, it means that motherIds != {-1, -1} because
/// the hadron comes from the decay of a beauty hadron. This can happen if and only if this is not a replaced one (i.e. native from Lambdab0 decay)
/// if we arrive here, it means that the hadron comes from the decay of a beauty hadron.
/// This can happen if and only if this is not a replaced one (i.e. native from Lambdab0 decay)
if (std::find(motherPdgCodes.begin(), motherPdgCodes.end(), 5122) == motherPdgCodes.end() && std::find(motherPdgCodes.begin(), motherPdgCodes.end(), -5122) == motherPdgCodes.end()) {
std::cerr << "The particle " << absPdg << " does not originate neither from a c/c-bar quark (replaced) nor from a Lambda_b0 decay. There is something wrong, aborting..." << std::endl;
return 1;
}
/// since this is a native Lc excited state from Lb0 decay, we must not update the replacement counters
updateCounters = false;
}
}
std::cout << std::endl;

/// update the counters only if the Lc excited state comes from a replaced prompt D*+, i.e. not from a Lb0 decay
if (!updateCounters) continue;

/// only if we arrive here it means that everything is ok, and we can safely update the counters for the final statistics
for (int iRepl{0}; iRepl<2; ++iRepl) {
if (absPdg == pdgReplParticles[iRepl][0]) {
Expand All @@ -133,7 +156,7 @@ int External() {
auto pdgDau = tracks->at(j).GetPdgCode();
pdgsDecay.push_back(pdgDau);
std::cout << " -- daughter " << j << ": " << pdgDau << std::endl;
if (pdgDau != 333) { // phi is antiparticle of itself
if (pdgDau != 333 && pdgDau != 111) { // phi and pi0 are antiparticles of themselves
pdgsDecayAntiPart.push_back(-pdgDau);
} else {
pdgsDecayAntiPart.push_back(pdgDau);
Expand Down Expand Up @@ -163,15 +186,15 @@ int External() {
std::cout <<"# signal hadrons: " << nSignals << "\n";
std::cout <<"# signal hadrons decaying in the correct channel: " << nSignalGoodDecay << "\n";

if (nEventsMB < nEvents * (1 - ratioTrigger) * 0.95 || nEventsMB > nEvents * (1 - ratioTrigger) * 1.05) { // we put some tolerance since the number of generated events is small
if (nEventsMB < nEvents * (1 - ratioTrigger) * 0.90 || nEventsMB > nEvents * (1 - ratioTrigger) * 1.10) { // we put some tolerance since the number of generated events is small
std::cerr << "Number of generated MB events different than expected\n";
return 1;
}
if (nEventsInjOne < nEvents * ratioTrigger * 0.5 * 0.95 || nEventsInjOne > nEvents * ratioTrigger * 0.5 * 1.05) {
if (nEventsInjOne < nEvents * ratioTrigger * 0.5 * 0.90 || nEventsInjOne > nEvents * ratioTrigger * 0.5 * 1.10) {
std::cerr << "Number of generated events injected with " << checkPdgQuarkOne << " different than expected\n";
return 1;
}
if (nEventsInjTwo < nEvents * ratioTrigger * 0.5 * 0.95 || nEventsInjTwo > nEvents * ratioTrigger * 0.5 * 1.05) {
if (nEventsInjTwo < nEvents * ratioTrigger * 0.5 * 0.90 || nEventsInjTwo > nEvents * ratioTrigger * 0.5 * 1.10) {
std::cerr << "Number of generated events injected with " << checkPdgQuarkTwo << " different than expected\n";
return 1;
}
Expand All @@ -184,7 +207,7 @@ int External() {

for (int iRepl{0}; iRepl<2; ++iRepl) {

std::cout << " --- pdgReplPartCounters[" << iRepl << "][1] = " << pdgReplPartCounters[iRepl][1] << ", freqRepl[" << iRepl <<"] = " << freqRepl[iRepl] << ", sumOrigReplacedParticles[pdgReplParticles[" << iRepl << "][0]] =" << sumOrigReplacedParticles[pdgReplParticles[iRepl][0]] << std::endl;
std::cout << " --- pdgReplPartCounters[" << iRepl << "][1] = " << pdgReplPartCounters[iRepl][1] << ", freqRepl[" << iRepl <<"] = " << freqRepl[iRepl] << ", sumOrigReplacedParticles[pdgReplParticles[" << iRepl << "][0]] = " << sumOrigReplacedParticles[pdgReplParticles[iRepl][0]] << std::endl;

if (std::abs(pdgReplPartCounters[iRepl][1] - freqRepl[iRepl] * sumOrigReplacedParticles[pdgReplParticles[iRepl][0]]) > 2 * std::sqrt(freqRepl[iRepl] * sumOrigReplacedParticles[pdgReplParticles[iRepl][0]])) { // 2 sigma compatibility
float fracMeas = 0.;
Expand Down