diff --git a/MC/config/PWGHF/ini/GeneratorHF_HFe_ccbar_and_bbar_gap5_Mode2.ini b/MC/config/PWGHF/ini/GeneratorHF_HFe_ccbar_and_bbar_gap5_Mode2.ini new file mode 100644 index 000000000..957daf1e8 --- /dev/null +++ b/MC/config/PWGHF/ini/GeneratorHF_HFe_ccbar_and_bbar_gap5_Mode2.ini @@ -0,0 +1,8 @@ +### The external generator derives from GeneratorPythia8. +[GeneratorExternal] +fileName=${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGHF/external/generator/generator_pythia8_gaptriggered_hf.C +funcName=GeneratorPythia8GapTriggeredCharmAndBeauty(5, -1.5, 1.5) + +[GeneratorPythia8] +config=${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGHF/pythia8/generator/pythia8_HFe_Mode2.cfg +includePartonEvent=true diff --git a/MC/config/PWGHF/ini/tests/GeneratorHF_HFe_ccbar_and_bbar_gap5_Mode2.C b/MC/config/PWGHF/ini/tests/GeneratorHF_HFe_ccbar_and_bbar_gap5_Mode2.C new file mode 100644 index 000000000..fb9df800e --- /dev/null +++ b/MC/config/PWGHF/ini/tests/GeneratorHF_HFe_ccbar_and_bbar_gap5_Mode2.C @@ -0,0 +1,128 @@ +int External() { + std::string path{"o2sim_Kine.root"}; + + int checkPdgDecayElectron = 11; + int checkPdgQuarkOne = 4; + int checkPdgQuarkTwo = 5; + float ratioTrigger = 1. / 5; // one event triggered out of 5 + + TFile file(path.c_str(), "READ"); + if (file.IsZombie()) { + std::cerr << "Cannot open ROOT file" << path << "\n"; + return 1; + } + + auto tree = (TTree *)file.Get("o2sim"); + if (!tree) { + std::cerr << "Cannot find tree o2sim in file" << path << "\n"; + return 1; + } + + std::vector *tracks{}; + tree->SetBranchAddress("MCTrack", &tracks); + o2::dataformats::MCEventHeader *eventHeader = nullptr; + tree->SetBranchAddress("MCEventHeader.", &eventHeader); + + int nEventsMB{}, nEventsInjOne{}, nEventsInjTwo{}; + int nQuarksOne{}, nQuarksTwo{}; + int nElectrons{}; + auto nEvents = tree->GetEntries(); + + for (int i = 0; i < nEvents; i++) { + tree->GetEntry(i); + + // check subgenerator information + if (eventHeader->hasInfo(o2::mcgenid::GeneratorProperty::SUBGENERATORID)) { + bool isValid = false; + int subGeneratorId = eventHeader->getInfo( + o2::mcgenid::GeneratorProperty::SUBGENERATORID, isValid); + if (subGeneratorId == 0) { + nEventsMB++; + } else if (subGeneratorId == checkPdgQuarkOne) { + nEventsInjOne++; + } else if (subGeneratorId == checkPdgQuarkTwo) { + nEventsInjTwo++; + } + } // if event header + + int nelectronsev = 0; + + for (auto &track : *tracks) { + auto pdg = track.GetPdgCode(); + if (std::abs(pdg) == checkPdgQuarkOne) { + nQuarksOne++; + continue; + } + if (std::abs(pdg) == checkPdgQuarkTwo) { + nQuarksTwo++; + continue; + } + + auto y = track.GetRapidity(); + if (std::abs(pdg) == checkPdgDecayElectron) { + int igmother = track.getMotherTrackId(); + auto gmTrack = (*tracks)[igmother]; + int gmpdg = gmTrack.GetPdgCode(); + if (int(std::abs(gmpdg) / 100.) == 4 || + int(std::abs(gmpdg) / 1000.) == 4 || + int(std::abs(gmpdg) / 100.) == 5 || + int(std::abs(gmpdg) / 1000.) == 5) { + nElectrons++; + nelectronsev++; + } // gmpdg + } // pdgdecay + } // loop track + // std::cout << "#electrons per event: " << nelectronsev << "\n"; + } + + std::cout << "--------------------------------\n"; + std::cout << "# Events: " << nEvents << "\n"; + std::cout << "# MB events: " << nEventsMB << "\n"; + std::cout << Form("# events injected with %d quark pair: ", checkPdgQuarkOne) + << nEventsInjOne << "\n"; + std::cout << Form("# events injected with %d quark pair: ", checkPdgQuarkTwo) + << nEventsInjTwo << "\n"; + std::cout << Form("# %d (anti)quarks: ", checkPdgQuarkOne) << nQuarksOne + << "\n"; + std::cout << Form("# %d (anti)quarks: ", checkPdgQuarkTwo) << nQuarksTwo + << "\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 + 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) { + 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) { + std::cerr << "Number of generated events injected with " << checkPdgQuarkTwo + << " different than expected\n"; + return 1; + } + if (nQuarksOne < + nEvents * + ratioTrigger) { // we expect anyway more because the same quark is + // repeated several time, after each gluon radiation + std::cerr << "Number of generated (anti)quarks " << checkPdgQuarkOne + << " lower than expected\n"; + return 1; + } + if (nQuarksTwo < + nEvents * + ratioTrigger) { // we expect anyway more because the same quark is + // repeated several time, after each gluon radiation + std::cerr << "Number of generated (anti)quarks " << checkPdgQuarkTwo + << " lower than expected\n"; + return 1; + } + std::cout << "#electrons: " << nElectrons << "\n"; + + return 0; +} // external diff --git a/MC/config/PWGHF/pythia8/generator/pythia8_HFe_Mode2.cfg b/MC/config/PWGHF/pythia8/generator/pythia8_HFe_Mode2.cfg new file mode 100644 index 000000000..fd44ab1e8 --- /dev/null +++ b/MC/config/PWGHF/pythia8/generator/pythia8_HFe_Mode2.cfg @@ -0,0 +1,104 @@ +### authors: Fabrizio Grosa (fabrizio.grosa@cern.ch) +### Grazia Luparello (Grazia.Luparello@cern.ch) +### Antonio Palasciano (antonio.palasciano@cern.ch) +### Jonghan Park (jonghan@cern.ch) +### electrons from heavy-flavour hadrons and from light neutral mesons +### last update: August 2025 + +### beams +Beams:idA 2212 # proton +Beams:idB 2212 # proton +Beams:eCM 13600. # GeV + +### processes +SoftQCD:inelastic on # all inelastic processes + +### decays +ParticleDecays:limitTau0 on +ParticleDecays:tau0Max 10. + +### switching on Pythia Mode2 +ColourReconnection:mode 1 +ColourReconnection:allowDoubleJunRem off +ColourReconnection:m0 0.3 +ColourReconnection:allowJunctions on +ColourReconnection:junctionCorrection 1.20 +ColourReconnection:timeDilationMode 2 +ColourReconnection:timeDilationPar 0.18 +StringPT:sigma 0.335 +StringZ:aLund 0.36 +StringZ:bLund 0.56 +StringFlav:probQQtoQ 0.078 +StringFlav:ProbStoUD 0.2 +StringFlav:probQQ1toQQ0join 0.0275,0.0275,0.0275,0.0275 +MultiPartonInteractions:pT0Ref 2.15 +BeamRemnants:remnantMode 1 +BeamRemnants:saturation 5 + +# Correct decay lengths (wrong in PYTHIA8 decay table) +# Lb +5122:tau0 = 0.4390 +# Xic0 +4132:tau0 = 0.0455 +# OmegaC +4332:tau0 = 0.0803 + + +### switch off all decay channels +111:onMode = off +221:onMode = off + +411:onMode = off +421:onMode = off +431:onMode = off +4122:onMode = off +4132:onMode = off +4232:onMode = off +4332:onMode = off + +511:onMode = off +521:onMode = off +531:onMode = off +5122:onMode = off +5132:onMode = off +5232:onMode = off +5332:onMode = off + +###Semimuonic decays of charm + +### pi0 -> e +111:onIfAny = 11 -11 +111:onIfAny = 22 11 -11 +### eta -> e +221:onIfAny = 11 -11 +221:onIfAny = 22 11 -11 + +### D+/- -> e + X +411:onIfAny = -11 12 +### D0 -> e + X +421:onIfAny = -11 12 +### D_s -> e + X +431:onIfAny = -11 12 +### Lambda_c -> e + X +4122:onIfAny = -11 12 +### Xsi0_c -> e + X +4132:onIfAny = -11 12 +### Xsi+_c -> e + X +4232:onIfAny = -11 12 +### Omega_c -> e + X +4332:onIfAny = -11 12 + +### B0 -> e + X +511:onIfAny = -11 12 +### B+/- -> e + X +521:onIfAny = -11 12 +### B_s -> e + X +531:onIfAny = -11 12 +### Lambda_b -> e + X +5122:onIfAny = -11 12 +### Xsi_b -> e + X +5132:onIfAny = -11 12 +### Xsi0_b -> e + X +5232:onIfAny = -11 12 +### Omega_b -> e + X +5332:onIfAny = -11 12