-
Notifications
You must be signed in to change notification settings - Fork 144
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
PWGHF: add config for Pb-Pb MCs without pt hard bins for signal (#1743)
* PWGHF: add config for Pb-Pb MCs without pt hard bins for signal * Add missing config
- Loading branch information
Showing
3 changed files
with
163 additions
and
18 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
8 changes: 8 additions & 0 deletions
8
MC/config/PWGHF/ini/GeneratorHF_D2H_ccbar_and_bbbar_PbPb_ptHardBins.ini
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,8 @@ | ||
| ### The external generator derives from GeneratorPythia8. | ||
| [GeneratorExternal] | ||
| fileName=${O2DPG_ROOT}/MC/config/PWGHF/external/generator/generator_pythia8_embed_hf.C | ||
| funcName=GeneratorPythia8EmbedHFCharmAndBeauty(true) | ||
|
|
||
| [GeneratorPythia8] | ||
| config=${O2DPG_ROOT}/MC/config/PWGHF/pythia8/generator/pythia8_charmhadronic_with_decays_Mode2_hardQCD_5TeV.cfg | ||
| includePartonEvent=true |
134 changes: 134 additions & 0 deletions
134
MC/config/PWGHF/ini/tests/GeneratorHF_D2H_ccbar_and_bbbar_PbPb_ptHardBins.C
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,134 @@ | ||
| int External() { | ||
| std::string path{"o2sim_Kine.root"}; | ||
| //std::string path{"tf1/sgn_1_Kine.root"}; | ||
|
|
||
| int checkPdgQuarkOne{4}; | ||
| int checkPdgQuarkTwo{5}; | ||
| float ratioTrigger = 1.; // one event triggered out of 1 | ||
| float averagePt = 0.; | ||
|
|
||
| std::vector<int> checkPdgHadron{411, 421, 431, 4122, 4132, 4232, 4332}; | ||
| std::map<int, std::vector<std::vector<int>>> 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}, {-313, 321}}}, // Ds+ | ||
| {4122, {{-313, 2212}, {-321, 2224}, {211, 102134}, {-321, 211, 2212}, {311, 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<o2::MCTrack> *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<int>(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 | ||
| averagePt += track.GetPt(); | ||
|
|
||
| std::vector<int> pdgsDecay{}; | ||
| std::vector<int> 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; | ||
| } | ||
| } | ||
| } | ||
| } | ||
| } | ||
|
|
||
| averagePt /= nSignals; | ||
|
|
||
| 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"; | ||
| std::cout <<"average pT of signal hadrons: " << averagePt << "\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; | ||
| } |