From d869a663170f9de8cac089ed2fa1c2c2bfd6c843 Mon Sep 17 00:00:00 2001 From: aimeric-landou Date: Sun, 5 Nov 2023 19:16:21 +0100 Subject: [PATCH] HF Jet MC - fix weighting, using new HF jet generator --- .../ini/GeneratorHFJETrigger_ccbar.ini | 8 ++ .../ini/tests/GeneratorHFJETrigger_ccbar.C | 125 ++++++++++++++++++ .../pythia8_jet_charmtriggers_with_decays.cfg | 117 ++++++++++++++++ MC/run/PWGGAJE/run_jets_HF_ccbar.sh | 37 ++---- 4 files changed, 258 insertions(+), 29 deletions(-) create mode 100644 MC/config/PWGGAJE/ini/GeneratorHFJETrigger_ccbar.ini create mode 100644 MC/config/PWGGAJE/ini/tests/GeneratorHFJETrigger_ccbar.C create mode 100644 MC/config/PWGGAJE/pythia8/generator/pythia8_jet_charmtriggers_with_decays.cfg diff --git a/MC/config/PWGGAJE/ini/GeneratorHFJETrigger_ccbar.ini b/MC/config/PWGGAJE/ini/GeneratorHFJETrigger_ccbar.ini new file mode 100644 index 000000000..72b484652 --- /dev/null +++ b/MC/config/PWGGAJE/ini/GeneratorHFJETrigger_ccbar.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=GeneratorPythia8GapTriggeredCharm(3, -5, 5) + +[GeneratorPythia8] +config=${O2DPG_ROOT}/MC/config/PWGGAJE/pythia8/generator/pythia8_jet_charmtriggers_with_decays.cfg + diff --git a/MC/config/PWGGAJE/ini/tests/GeneratorHFJETrigger_ccbar.C b/MC/config/PWGGAJE/ini/tests/GeneratorHFJETrigger_ccbar.C new file mode 100644 index 000000000..ef24e0d3e --- /dev/null +++ b/MC/config/PWGGAJE/ini/tests/GeneratorHFJETrigger_ccbar.C @@ -0,0 +1,125 @@ +int External() { + std::string path{"o2sim_Kine.root"}; + int checkPdgQuark{4}; + float ratioTrigger = 1. / 3; // one event triggered out of 3 + std::vector checkPdgHadron{411, 421, 431, 443, 4122, 4132, 4232, 4332}; + std::map>> checkHadronDecays{ + // sorted pdg of daughters + {411, {{-321, 211, 211}, {-313, 211}, {211, 311}, {211, 333}}}, // D+ + {421, {{-321, 211}, {-321, 111, 211}}}, // D0 + {431, {{211, 333}}}, // Ds+ + {443, {{-11, 11}}}, // Jpsi + {4122, + {{-313, 2212}, {-321, 2224}, {211, 3124}, {-321, 211, 2212}}}, // Lc+ + {4132, {{211, 3312}}}, // Xic0 + {4232, + {{-313, 2212}, + {-321, 3324}, + {211, 211, 3312}, + {-321, 211, 2212}}}, // Xic+ + {4332, {{211, 3334}}} // Omegac+ + }; + + 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); + + int nQuarks{}, nSignals{}, nSignalGoodDecay{}; + auto nEvents = tree->GetEntries(); + + // Setting up event header to access event weight info + o2::dataformats::MCEventHeader *eventHeader = nullptr; + ; + tree->SetBranchAddress("MCEventHeader.", &eventHeader); + int weight1Counter = 0; + bool isvalid; + + for (int i = 0; i < nEvents; i++) { + tree->GetEntry(i); + + float eventWeight = eventHeader->getInfo("weight", isvalid); + if (!isvalid) { + std::cerr << "Could not retrieve event weight from MCEventHeader\n"; + return 1; + } + if (abs(eventWeight - 1) < + 1E-5) { // checks that the event weights are not unitary + weight1Counter++; + } + + for (auto &track : *tracks) { + auto pdg = track.GetPdgCode(); + if (std::abs(pdg) == checkPdgQuark) { + nQuarks++; + continue; + } + if (std::find(checkPdgHadron.begin(), checkPdgHadron.end(), + std::abs(pdg)) != checkPdgHadron.end()) // found signal + { + // count signal PDG + nSignals++; + + 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 << Form("# %d (anti)quarks: ", checkPdgQuark) << nQuarks << "\n"; + std::cout << "# signal hadrons: " << nSignals << "\n"; + std::cout << "# signal hadrons decaying in the correct channel: " + << nSignalGoodDecay << "\n"; + + if (nQuarks < + 2 * 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 " << checkPdgQuark + << " 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; + } + + if (weight1Counter == nEvents) { + std::cerr << "All events have a unitary weight\n"; + return 1; + } + + return 0; +} diff --git a/MC/config/PWGGAJE/pythia8/generator/pythia8_jet_charmtriggers_with_decays.cfg b/MC/config/PWGGAJE/pythia8/generator/pythia8_jet_charmtriggers_with_decays.cfg new file mode 100644 index 000000000..9c0089c1c --- /dev/null +++ b/MC/config/PWGGAJE/pythia8/generator/pythia8_jet_charmtriggers_with_decays.cfg @@ -0,0 +1,117 @@ +### author: Aimeric Landou (aimeric.landou@cern.ch) +### based on Fabrizio Grosa's HF template (fabrizio.grosa@cern.ch) +### since: January 2023 + +### 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. + +### Force golden charm hadrons decay modes for trigger 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.00276 0 333 211 +## 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 3124 211 +4122:addChannel = 1 0.035 0 2212 -321 211 +### add Xic+ decays absent in PYTHIA8 decay table +4232:addChannel = 1 0.2 0 2212 -313 +4232:addChannel = 1 0.2 0 2212 -321 211 +4232:addChannel = 1 0.2 0 3324 211 +4232:addChannel = 1 0.2 0 3312 211 211 +### add Xic0 decays absent in PYTHIA8 decay table +4132:addChannel = 1 0.2 0 3312 211 + +### 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 +4112:onMode = off +4122:onMode = off +4232:onMode = off +4132:onMode = off +443:onMode = off +4332: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 -> 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 = 3124 211 +### Lambda_c -> p K pi +4122:onIfMatch = 2212 321 211 + +### Xic+ -> pK*0 +4232:onIfMatch = 2212 313 +### Xic+ -> p K- pi+ +4232:onIfMatch = 2212 321 211 +### Xic+ -> Xi*0 pi+, Xi*->Xi- pi+ +4232:onIfMatch = 3324 211 +### Xic+ -> Xi- pi+ pi+ +4232:onIfMatch = 3312 211 211 + +### Xic0 -> Xi- pi+ +4132:onIfMatch = 3312 211 + +### Omega_c -> Omega pi +4332:onIfMatch = 3334 211 + +### Jpsi -> ee +443:onIfMatch = 11 11 + +# Correct Lb decay length (wrong in PYTHIA8 decay table) +5122:tau0 = 4.41000e-01 + +# Correct OmegaC decay length (wrong in PYTHIA8 decay table) +4332:tau0 = 0.08000000000 + + +### Jets options +PhaseSpace:pTHatMin = 5.0 +PhaseSpace:pTHatMax = 300.0 +PhaseSpace:bias2Selection = on +PhaseSpace:bias2SelectionPow = 6.0 +SoftQCD:inelastic off # all inelastic processes +HardQCD:all on \ No newline at end of file diff --git a/MC/run/PWGGAJE/run_jets_HF_ccbar.sh b/MC/run/PWGGAJE/run_jets_HF_ccbar.sh index 3149f5092..9b155ba2b 100644 --- a/MC/run/PWGGAJE/run_jets_HF_ccbar.sh +++ b/MC/run/PWGGAJE/run_jets_HF_ccbar.sh @@ -15,42 +15,21 @@ RNDSEED=${RNDSEED:-0} # [default = 0] time-based random seed -NSIGEVENTS=${NSIGEVENTS:-10} +NSIGEVENTS=${NSIGEVENTS:-5} NTIMEFRAMES=${NTIMEFRAMES:-1} NWORKERS=${NWORKERS:-8} MODULES="--skipModules ZDC" #"PIPE ITS TPC EMCAL" CONFIG_ENERGY=${CONFIG_ENERGY:-13600.0} SIMENGINE=${SIMENGINE:-TGeant4} -WEIGHTPOW=${WEIGHTPOW:-6.0} [[ ${SPLITID} != "" ]] && SEED="-seed ${SPLITID}" || SEED="" -# Default for weighted productions -PTHATMIN=${PTHATMIN:-5.0} -PTHATMAX=${PTHATMAX:-300.0} - -# Define the pt hat bin arrays -pthatbin_loweredges=(0 5 7 9 12 16 21 28 36 45 57 70 85 99 115 132 150 169 190 212 235) -pthatbin_higheredges=( 5 7 9 12 16 21 28 36 45 57 70 85 99 115 132 150 169 190 212 235 -1) - -# Recover environmental vars for pt binning -#PTHATBIN=${PTHATBIN:-1} - -if [ -z "$PTHATBIN" ]; then - echo "Open Pt-hat range set" -else - PTHATMIN=${pthatbin_loweredges[$PTHATBIN]} - PTHATMAX=${pthatbin_higheredges[$PTHATBIN]} -fi - - -#ccbar filter -${O2DPG_ROOT}/MC/bin/o2dpg_sim_workflow.py -eCM ${CONFIG_ENERGY} -col pp -gen external -proc "jets" \ - -ptHatMin ${PTHATMIN} -ptHatMax ${PTHATMAX} \ - -tf ${NTIMEFRAMES} -ns ${NSIGEVENTS} -e ${SIMENGINE} \ - -j ${NWORKERS} -mod "--skipModules ZDC" \ - -interactionRate 500000 -confKey "Diamond.width[2]=6." ${SEED} \ - -ini $O2DPG_ROOT/MC/config/PWGHF/ini/GeneratorHFTrigger_ccbar.ini \ - -weightPow ${WEIGHTPOW} +#ccbar filter and bias2SelectionPow and PtHat settings are in the ini file given below +${O2DPG_ROOT}/MC/bin/o2dpg_sim_workflow.py -eCM ${CONFIG_ENERGY} -col pp -gen external -proc "jets" \ + -tf ${NTIMEFRAMES} -ns ${NSIGEVENTS} -e ${SIMENGINE} \ + -j ${NWORKERS} -mod "--skipModules ZDC" \ + -interactionRate 500000 -confKey "Diamond.width[2]=6." ${SEED} \ + -ini $O2DPG_ROOT/MC/config/PWGGAJE/ini/GeneratorHFJETrigger_ccbar.ini + # run workflow # allow increased timeframe parallelism with --cpu-limit 32