From 3e86b83164651038646289cd245b81676ae46e1f Mon Sep 17 00:00:00 2001 From: Luca610 Date: Mon, 22 Jul 2024 17:41:07 +0200 Subject: [PATCH 1/4] Simulation workflow for D resonances --- ...ratorHF_D2H_ccbar_and_bbbar_gap5_DReso.ini | 8 + ...neratorHF_D2H_ccbar_and_bbbar_gap5_DReso.C | 140 +++++++++++ ...ythia8_charmhadronic_with_decays_DReso.cfg | 236 ++++++++++++++++++ 3 files changed, 384 insertions(+) create mode 100644 MC/config/PWGHF/ini/GeneratorHF_D2H_ccbar_and_bbbar_gap5_DReso.ini create mode 100644 MC/config/PWGHF/ini/tests/GeneratorHF_D2H_ccbar_and_bbbar_gap5_DReso.C create mode 100644 MC/config/PWGHF/pythia8/generator/pythia8_charmhadronic_with_decays_DReso.cfg diff --git a/MC/config/PWGHF/ini/GeneratorHF_D2H_ccbar_and_bbbar_gap5_DReso.ini b/MC/config/PWGHF/ini/GeneratorHF_D2H_ccbar_and_bbbar_gap5_DReso.ini new file mode 100644 index 000000000..9f3037a4b --- /dev/null +++ b/MC/config/PWGHF/ini/GeneratorHF_D2H_ccbar_and_bbbar_gap5_DReso.ini @@ -0,0 +1,8 @@ +### The external generator derives from GeneratorPythia8. +[GeneratorExternal] +fileName=${O2DPG_ROOT}/MC/config/PWGHF/external/generator/generator_pythia8_gaptriggered_hf.C +funcName=GeneratorPythia8GapTriggeredCharmAndBeauty(5, -1.5, 1.5) + +[GeneratorPythia8] +config=${O2DPG_ROOT}/MC/config/PWGHF/pythia8/generator/pythia8_charmhadronic_with_decays_DReso.cfg +includePartonEvent=true diff --git a/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_ccbar_and_bbbar_gap5_DReso.C b/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_ccbar_and_bbbar_gap5_DReso.C new file mode 100644 index 000000000..c16934408 --- /dev/null +++ b/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_ccbar_and_bbbar_gap5_DReso.C @@ -0,0 +1,140 @@ +int External() { + std::string path{"/home/luca/alice/Ds_reso/MCReso/07_22_16-59-12/tf1/genevents_Kine.root"}; + + int checkPdgQuarkOne{4}; + int checkPdgQuarkTwo{5}; + float ratioTrigger = 1./5.; // one event triggered out of 5 + + std::vector checkPdgHadron{411, 415, 421, 425, 431, 435, 511, 521, 531, 4122, 10411, 10421, 10433, 20423, 20433}; + std::map>> checkHadronDecays{ // sorted pdg of daughters + {411, {{-321, 211, 211}, {-313, 211}, {211, 311}, {211, 333}}}, // D+ + {415, {{211, 421}}}, // D2*(2460)+ + {421, {{-321, 211}, {-321, 111, 211}}}, // D0 + {425, {{-211, 413},{-211, 411}}}, // D2*(2460)0 + {431, {{211, 333}, {-313, 321}}}, // Ds+ + {425, {{311, 413}, {311, 411}, {321,421}}}, // Ds2*(2573) + {511, {{-415, -11, 12}, {-10411, -11, 12}, {-415, -13, 14}, {-10411, -13, 14}, {-415, -15, 16}, {-10411, -15, 16}, + {-10411, 211}, {-10421, 211}, {-415, 433}, { -415, 431}, {-415, 211}, {-415, 213} }}, // B0 + {521, {{-20423, -11, 12}, {-425, -11, 12}, {-10421, -11, 12}, {-20423, -13, 14}, {-425, -13, 14}, {-10421, -13, 14}, + {-20423, -15, 16}, {-425, -15, 16}, {-10421, -15, 16}, {-20423, 211}, {-20423, 213}, {-20423, 431}, {-20423, 433}, + {-425, 211}, {-425, 213}, {-425, 431}, {-425, 433}}}, // B+ + {531, {{-435, -11, 12}, {-10433, -11, 12}, {-435, -13, 14}, {-10433, -13, 14}, {-435, -15, 16}, {-10433, -15, 16}, + {-435, 211}, {-20433, 211}, {-20433, 213}}},// Bs0 + {4122, {{-313, 2212}, {-321, 2224}, {211, 3124}, {-321, 211, 2212}, {311, 2212}}}, // Lc+ + {10411, {{211, 421}}}, // D0*+ + {10421, {{-211, 411}}}, // D0*0 + {10433, {{311, 413}}}, // Ds1(2536) + {20423, {{-211, 413}}}, // D1(2430)0 + {20433, {{22, 431}, {-211, 211, 431} }} // Ds1 (2460) + }; + + 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"); + std::vector *tracks{}; + tree->SetBranchAddress("MCTrack", &tracks); + o2::dataformats::MCEventHeader *eventHeader = nullptr; + tree->SetBranchAddress("MCEventHeader.", &eventHeader); + + int nEventsMB{}, nEventsInjOne{}, nEventsInjTwo{}; + int nQuarksOne{}, nQuarksTwo{}, nSignals{}, nSignalGoodDecay{}; + 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++; + } + } + + for (auto &track : *tracks) { + auto pdg = track.GetPdgCode(); + if (std::abs(pdg) == checkPdgQuarkOne) { + nQuarksOne++; + continue; + } + if (std::abs(pdg) == checkPdgQuarkTwo) { + nQuarksTwo++; + continue; + } + if (std::find(checkPdgHadron.begin(), checkPdgHadron.end(), std::abs(pdg)) != checkPdgHadron.end()) { // found signal + nSignals++; // count signal PDG + + std::vector pdgsDecay{}; + std::vector pdgsDecayAntiPart{}; + for (int j{track.getFirstDaughterTrackId()}; j <= track.getLastDaughterTrackId(); ++j) { + auto pdgDau = tracks->at(j).GetPdgCode(); + pdgsDecay.push_back(pdgDau); + if (pdgDau != 333) { // phi is antiparticle of itself + pdgsDecayAntiPart.push_back(-pdgDau); + } else { + pdgsDecayAntiPart.push_back(pdgDau); + } + } + + std::sort(pdgsDecay.begin(), pdgsDecay.end()); + std::sort(pdgsDecayAntiPart.begin(), pdgsDecayAntiPart.end()); + + for (auto &decay : checkHadronDecays[std::abs(pdg)]) { + if (pdgsDecay == decay || pdgsDecayAntiPart == decay) { + nSignalGoodDecay++; + break; + } + } + } + } + } + + 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"; + 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 + 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; + } + + float fracForcedDecays = float(nSignalGoodDecay) / nSignals; + if (fracForcedDecays < 0.9) { // we put some tolerance (e.g. due to oscillations which might change the final state) + std::cerr << "Fraction of signals decaying into the correct channel " << fracForcedDecays << " lower than expected\n"; + return 1; + } + + return 0; +} diff --git a/MC/config/PWGHF/pythia8/generator/pythia8_charmhadronic_with_decays_DReso.cfg b/MC/config/PWGHF/pythia8/generator/pythia8_charmhadronic_with_decays_DReso.cfg new file mode 100644 index 000000000..a82bd5a39 --- /dev/null +++ b/MC/config/PWGHF/pythia8/generator/pythia8_charmhadronic_with_decays_DReso.cfg @@ -0,0 +1,236 @@ +### author: Stefano Politano (stefano.politano@cern.ch) +### since: July 2024 + +### beams +Beams:idA 2212 # proton +Beams:idB 2212 # proton +Beams:eCM 13500. # GeV + +### processes +SoftQCD:inelastic on # all inelastic processes + +### decays +ParticleDecays:limitTau0 on +ParticleDecays:tau0Max 10. + +### switching on Pythia Mode2 +StringFlav:mesonCL1S0J1=0.14 # from Run 2 +StringFlav:mesonCL1S1J2=0.35 # from Run 2 +StringFlav:mesonCL1S1J0=0.35 + + +### turn off all charm resonances decays +10433:onMode = off # Ds1(2536) +435:onMode = off # Ds2*(2573) +425:onMode = off # D2*(2460)0 +415:onMode = off # D2*(2460)+ +10411:onMode = off # D0*+ +10421:onMode = off # D0*0 +20433:onMode = off # Ds1 (2460) +20423:onMode = off # D1(2430)0 + +### turn off all beauty hadron decays +531:onMode = off +511:onMode = off +521:onMode = off + +### add D1(2430)0 +20423:oneChannel = 1 1 0 413 -211 +20423:onIfMatch = 413 -211 + +### add Ds1(2536) +10433:oneChannel = 1 1 0 413 311 +10433:onIfMatch = 413 311 + +### Ds1 (2460) +20433:oneChannel = 1 0.5 0 431 22 +20433:addChannel = 1 0.5 0 431 211 -211 +20433:onIfMatch = 431 22 +20433:onIfMatch = 431 211 211 + +### add Ds2*(2573) +435:oneChannel = 1 0.0500000 0 413 311 +435:addChannel = 1 0.4500000 0 411 311 +435:addChannel = 1 0.4500000 0 421 321 +435:onIfMatch = 413 311 +435:onIfMatch = 411 311 +435:onIfMatch = 421 321 + +### add D2*(2460)0 +425:oneChannel = 1 0.5 0 413 -211 +425:addChannel = 1 0.5 0 411 -211 +425:onIfMatch = 413 211 +425:onIfMatch = 411 211 + +### add D2*(2460)+ +415:oneChannel = 1 1 0 421 211 +415:onIfMatch = 421 211 + + +### add D0*+ +10411:oneChannel = 1 1 0 421 211 +10411:onIfMatch = 421 211 + +### add D0*0 +10421:oneChannel = 1 1 0 411 -211 +10421:onIfMatch = 411 211 + + +### add Bs0 +531:oneChannel = 1 0.0070000 0 12 -11 -435 +531:addChannel = 1 0.0070000 0 12 -11 -10433 +531:addChannel = 1 0.0070000 0 14 -13 -435 +531:addChannel = 1 0.0070000 0 14 -13 -10433 +531:addChannel = 1 0.0040000 0 14 -13 -20433 +531:addChannel = 1 0.0160000 0 16 -15 -433 +531:addChannel = 1 0.0028000 0 16 -15 -435 +531:addChannel = 1 0.0028000 0 16 -15 -10433 +531:addChannel = 1 0.0013000 0 -435 211 +531:addChannel = 1 0.0008000 0 -20433 211 +531:addChannel = 1 0.0021000 0 -20433 213 + +531:onIfMatch = 12 11 435 +531:onIfMatch = 12 11 10433 +531:onIfMatch = 14 13 435 +531:onIfMatch = 14 13 10433 +531:onIfMatch = 14 13 20433 +531:onIfMatch = 16 15 433 +531:onIfMatch = 16 15 435 +531:onIfMatch = 16 15 10433 +531:onIfMatch = 435 211 +531:onIfMatch = 20433 211 +531:onIfMatch = 20433 213 + +### add B0 +511:oneChannel = 1 0.0023000 0 12 -11 -415 +511:addChannel = 1 0.0045000 0 12 -11 -10411 +511:addChannel = 1 0.0023000 0 14 -13 -415 +511:addChannel = 1 0.0045000 0 14 -13 -10411 +511:addChannel = 1 0.0020000 0 16 -15 -415 +511:addChannel = 1 0.0013000 0 16 -15 -10411 +511:addChannel = 1 0.0002000 0 -10411 211 +511:addChannel = 1 0.0009100 0 -10421 211 +511:addChannel = 1 0.0040000 0 433 -415 +511:addChannel = 1 0.0042000 0 431 -415 +511:addChannel = 1 0.0009000 0 -415 211 +511:addChannel = 1 0.0022000 0 -415 213 + +511:onIfMatch = 12 11 415 +511:onIfMatch = 12 11 10411 +511:onIfMatch = 14 13 415 +511:onIfMatch = 14 13 10411 +511:onIfMatch = 16 15 415 +511:onIfMatch = 16 15 10411 +511:onIfMatch = 10411 211 +511:onIfMatch = 10421 211 +511:onIfMatch = 433 415 +511:onIfMatch = 431 415 +511:onIfMatch = 415 211 +511:onIfMatch = 415 213 + +### add B+ +521:oneChannel = 1 0.0090000 0 12 -11 -20423 +521:addChannel = 1 0.0090000 0 14 -13 -20423 +521:addChannel = 1 0.0020000 0 16 -15 -20423 +521:addChannel = 1 0.0030000 0 12 -11 -425 +521:addChannel = 1 0.0030000 0 14 -13 -425 +521:addChannel = 1 0.0020000 0 16 -15 -425 +521:addChannel = 1 0.0049000 0 12 -11 -10421 +521:addChannel = 1 0.0049000 0 14 -13 -10421 +521:addChannel = 1 0.0013000 0 16 -15 -10421 +521:addChannel = 1 0.0007500 0 -20423 211 +521:addChannel = 1 0.0022000 0 -20423 213 +521:addChannel = 1 0.0006000 0 -20423 431 +521:addChannel = 1 0.0012000 0 -20423 433 +521:addChannel = 1 0.0008000 0 -425 211 +521:addChannel = 1 0.0038000 0 -425 213 +521:addChannel = 1 0.0042000 0 431 -425 +521:addChannel = 1 0.0040000 0 433 -425 + +521:onIfMatch = 12 11 20423 +521:onIfMatch = 14 13 20423 +521:onIfMatch = 16 15 20423 +521:onIfMatch = 12 11 425 +521:onIfMatch = 14 13 425 +521:onIfMatch = 16 15 425 +521:onIfMatch = 12 11 10421 +521:onIfMatch = 14 13 10421 +521:onIfMatch = 16 15 10421 +521:onIfMatch = 20423 211 +521:onIfMatch = 20423 213 +521:onIfMatch = 20423 431 +521:onIfMatch = 20423 433 +521:onIfMatch = 425 211 +521:onIfMatch = 425 213 +521:onIfMatch = 431 425 +521:onIfMatch = 433 425 + +# Correct decay lengths (wrong in PYTHIA8 decay table) +# Lb +5122:tau0 = 0.4390 +# Xic0 +4132:tau0 = 0.0455 +# OmegaC +4332:tau0 = 0.0803 + +### Force golden charm hadrons decay modes for D2H studies +### add D+ decays absent in PYTHIA8 decay table and set BRs from PDG for other +411:oneChannel = 1 0.0752 0 -321 211 211 +411:addChannel = 1 0.0104 0 -313 211 +411:addChannel = 1 0.0156 0 311 211 +411:addChannel = 1 0.0752 0 333 211 # to have the same amount of D+->KKpi and D+->Kpipi +## add Lc decays absent in PYTHIA8 decay table and set BRs from PDG for other +4122:oneChannel = 1 0.0196 100 2212 -313 +4122:addChannel = 1 0.0108 100 2224 -321 +4122:addChannel = 1 0.022 100 102134 211 +4122:addChannel = 1 0.035 0 2212 -321 211 +4122:addChannel = 1 0.0159 0 2212 311 + +### K* -> K pi +313:onMode = off +313:onIfAll = 321 211 +### for Ds -> Phi pi+ +333:onMode = off +333:onIfAll = 321 321 +### for D0 -> rho0 pi+ k- +113:onMode = off +113:onIfAll = 211 211 +### for Lambda_c -> Delta++ K- +2224:onMode = off +2224:onIfAll = 2212 211 +### for Lambda_c -> Lambda(1520) K- +102134:onMode = off +102134:onIfAll = 2212 321 + +### switch off all decay channels +411:onMode = off +421:onMode = off +431:onMode = off +4122:onMode = off + +### D0 -> K pi +421:onIfMatch = 321 211 + +### D+/- -> K pi pi +411:onIfMatch = 321 211 211 +### D+/- -> K* pi +411:onIfMatch = 313 211 +### D+/- -> phi pi +411:onIfMatch = 333 211 + +### D_s -> K K* +431:onIfMatch = 321 313 +### D_s -> Phi pi +431:onIfMatch = 333 211 + +### Lambda_c -> p K* +4122:onIfMatch = 2212 313 +### Lambda_c -> Delta K +4122:onIfMatch = 2224 321 +### Lambda_c -> Lambda(1520) pi +4122:onIfMatch = 102134 211 +### Lambda_c -> p K pi +4122:onIfMatch = 2212 321 211 +### Lambda_c -> pK0s +4122:onIfMatch = 2212 311 + From a51cb3db97e8e1d8f6cf09b0ac75d20f1a8ee645 Mon Sep 17 00:00:00 2001 From: Luca610 Date: Mon, 22 Jul 2024 18:04:00 +0200 Subject: [PATCH 2/4] fixed typo --- .../ini/tests/GeneratorHF_D2H_ccbar_and_bbbar_gap5_DReso.C | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_ccbar_and_bbbar_gap5_DReso.C b/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_ccbar_and_bbbar_gap5_DReso.C index c16934408..83cc1871d 100644 --- a/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_ccbar_and_bbbar_gap5_DReso.C +++ b/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_ccbar_and_bbbar_gap5_DReso.C @@ -1,5 +1,5 @@ int External() { - std::string path{"/home/luca/alice/Ds_reso/MCReso/07_22_16-59-12/tf1/genevents_Kine.root"}; + std::string path{"o2sim_Kine.root"}; int checkPdgQuarkOne{4}; int checkPdgQuarkTwo{5}; From e0138e7632db0903ec4b6692e14f2d7ef832f4fe Mon Sep 17 00:00:00 2001 From: Luca610 Date: Mon, 22 Jul 2024 18:10:11 +0200 Subject: [PATCH 3/4] fixed typo --- .../pythia8_charmhadronic_with_decays_DReso.cfg | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/MC/config/PWGHF/pythia8/generator/pythia8_charmhadronic_with_decays_DReso.cfg b/MC/config/PWGHF/pythia8/generator/pythia8_charmhadronic_with_decays_DReso.cfg index a82bd5a39..166db783b 100644 --- a/MC/config/PWGHF/pythia8/generator/pythia8_charmhadronic_with_decays_DReso.cfg +++ b/MC/config/PWGHF/pythia8/generator/pythia8_charmhadronic_with_decays_DReso.cfg @@ -1,10 +1,10 @@ -### author: Stefano Politano (stefano.politano@cern.ch) +### author: Stefano Politano, Luca Aglietta (stefano.politano@cern.ch, luca.aglietta@cern.ch) ### since: July 2024 ### beams Beams:idA 2212 # proton Beams:idB 2212 # proton -Beams:eCM 13500. # GeV +Beams:eCM 13600. # GeV ### processes SoftQCD:inelastic on # all inelastic processes @@ -13,10 +13,11 @@ SoftQCD:inelastic on # all inelastic processes ParticleDecays:limitTau0 on ParticleDecays:tau0Max 10. -### switching on Pythia Mode2 -StringFlav:mesonCL1S0J1=0.14 # from Run 2 -StringFlav:mesonCL1S1J2=0.35 # from Run 2 -StringFlav:mesonCL1S1J0=0.35 +# parameters to boost resonances production +StringFlav:mesonCL1S0J1= 3 +StringFlav:mesonCL1S1J2= 3.2 +StringFlav:mesonCL1S1J0= 0.75 +StringFlav:mesonCL1S1J1= 0.5 ### turn off all charm resonances decays From eb01fe118029ceaacf8377cc16a00a11e5714955 Mon Sep 17 00:00:00 2001 From: Luca Aglietta <75362880+Luca610@users.noreply.github.com> Date: Mon, 22 Jul 2024 18:30:42 +0200 Subject: [PATCH 4/4] Update pythia8_charmhadronic_with_decays_DReso.cfg, change CL1S1J1 parameter --- .../generator/pythia8_charmhadronic_with_decays_DReso.cfg | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/MC/config/PWGHF/pythia8/generator/pythia8_charmhadronic_with_decays_DReso.cfg b/MC/config/PWGHF/pythia8/generator/pythia8_charmhadronic_with_decays_DReso.cfg index 166db783b..7713a5518 100644 --- a/MC/config/PWGHF/pythia8/generator/pythia8_charmhadronic_with_decays_DReso.cfg +++ b/MC/config/PWGHF/pythia8/generator/pythia8_charmhadronic_with_decays_DReso.cfg @@ -17,7 +17,7 @@ ParticleDecays:tau0Max 10. StringFlav:mesonCL1S0J1= 3 StringFlav:mesonCL1S1J2= 3.2 StringFlav:mesonCL1S1J0= 0.75 -StringFlav:mesonCL1S1J1= 0.5 +StringFlav:mesonCL1S1J1= 1. ### turn off all charm resonances decays