From cc8ab57560565ac83281baa52e1f71648acd695f Mon Sep 17 00:00:00 2001 From: DanielSamitz Date: Sun, 12 May 2024 14:13:16 +0200 Subject: [PATCH] HF cocktail with forced decays --- .../Generator_pythia8_forcedDecays.C | 151 ++++++++++++++++++ .../PWGEM/ini/Pythia8_Beauty_Cocktail.ini | 12 ++ .../PWGEM/ini/Pythia8_Charm_Cocktail.ini | 12 ++ .../pythia8/decayer/force_semileptonic.cfg | 80 ++++++++++ .../pythia8/generator/pythia8_hf_cocktail.cfg | 48 ++++++ .../pythia8/hooks/pythia8_userhooks_HF.C | 63 -------- 6 files changed, 303 insertions(+), 63 deletions(-) create mode 100644 MC/config/PWGEM/external/generator/Generator_pythia8_forcedDecays.C create mode 100644 MC/config/PWGEM/ini/Pythia8_Beauty_Cocktail.ini create mode 100644 MC/config/PWGEM/ini/Pythia8_Charm_Cocktail.ini create mode 100644 MC/config/PWGEM/pythia8/decayer/force_semileptonic.cfg create mode 100644 MC/config/PWGEM/pythia8/generator/pythia8_hf_cocktail.cfg delete mode 100644 MC/config/PWGEM/pythia8/hooks/pythia8_userhooks_HF.C diff --git a/MC/config/PWGEM/external/generator/Generator_pythia8_forcedDecays.C b/MC/config/PWGEM/external/generator/Generator_pythia8_forcedDecays.C new file mode 100644 index 000000000..02b78f2fb --- /dev/null +++ b/MC/config/PWGEM/external/generator/Generator_pythia8_forcedDecays.C @@ -0,0 +1,151 @@ +#include "Generators/DecayerPythia8.h" +#include "Generators/GeneratorPythia8.h" +#include "SimulationDataFormat/MCGenProperties.h" +#include "SimulationDataFormat/ParticleStatus.h" +#include "TLorentzVector.h" + +namespace o2 { +namespace eventgen { + +class DecayerPythia8ForceDecays : public DecayerPythia8 { +public: + DecayerPythia8ForceDecays() : DecayerPythia8(){}; + ~DecayerPythia8ForceDecays() = default; + + void calculateWeights(std::vector &pdgs) { + TLorentzVector mom = TLorentzVector(0., 0., 0., 9999999.); + for (int pdg : pdgs) { + Decay(pdg, &mom); // do one fake decay to initalize everything correctly + auto particleData = mPythia.particleData.particleDataEntryPtr(pdg); + double weight = 0.; + for (int c = 0; c < particleData->sizeChannels(); c++) { + weight += particleData->channel(c).currentBR(); + } + LOG(info) << "PDG = " << pdg + << ": sum of branching ratios of active decay channels = " + << weight; + mWeights[pdg] = weight; + mPythia.particleData.mayDecay( + pdg, false); // make sure it is not decayed in any other way (should + // already be switched off in the decayer config file) + } + } + + void forceDecays(std::vector &mParticles, int mother_pos) { + TParticle p = mParticles[mother_pos]; + int pdg = p.GetPdgCode(); + TLorentzVector mom = TLorentzVector(p.Px(), p.Py(), p.Pz(), p.Energy()); + Decay(pdg, &mom); + TClonesArray daughters = TClonesArray("TParticle"); + int nParticles = ImportParticles(&daughters); + p.SetStatusCode(o2::mcgenstatus::MCGenStatusEncoding(2, -11).fullEncoding); + p.SetBit(ParticleStatus::kToBeDone, false); + double mother_weight = p.GetWeight(); + TParticle *mother = static_cast(daughters[0]); + int mParticles_size = mParticles.size(); + p.SetFirstDaughter(mother->GetFirstDaughter() + mParticles_size - 1); + p.SetLastDaughter(mother->GetLastDaughter() + mParticles_size - 1); + for (int j = 1; j < nParticles; + j++) { // start loop at 1 to not include mother + TParticle *d = static_cast(daughters[j]); + double decay_weight = mWeights[abs(pdg)]; + if (decay_weight == 0) { + LOG(error) << "Decaying particle (PDG = " << pdg + << ") with decay weight = 0. Did you set the pdg codes for " + "calculating weights correctly?"; + } + d->SetWeight(decay_weight * mother_weight); + if (d->GetStatusCode() == 1) { + d->SetStatusCode( + o2::mcgenstatus::MCGenStatusEncoding(1, 91).fullEncoding); + d->SetBit(ParticleStatus::kToBeDone, true); + } else { + d->SetStatusCode( + o2::mcgenstatus::MCGenStatusEncoding(2, -91).fullEncoding); + d->SetBit(ParticleStatus::kToBeDone, false); + } + int firstmother = d->GetFirstMother(); + int firstdaughter = d->GetFirstDaughter(); + int lastdaughter = d->GetLastDaughter(); + if (firstmother == 0) { + d->SetFirstMother(mother_pos); + d->SetLastMother(mother_pos); + } else { + d->SetFirstMother(firstmother + mParticles_size - 1); + d->SetLastMother(firstmother + mParticles_size - 1); + } + if (firstdaughter == 0) { + d->SetFirstDaughter(-1); + } else { + d->SetFirstDaughter(firstdaughter + mParticles_size - 1); + } + if (lastdaughter == 0) { + d->SetLastDaughter(-1); + } else { + d->SetLastDaughter(lastdaughter + mParticles_size - 1); + } + mParticles.push_back(*d); + } + } + +private: + std::map mWeights; +}; + +class GeneratorPythia8ForcedDecays : public GeneratorPythia8 { + +public: + GeneratorPythia8ForcedDecays() : GeneratorPythia8(){}; + ~GeneratorPythia8ForcedDecays() = default; + + // overriden methods + bool Init() override { return GeneratorPythia8::Init() && InitDecayer(); }; + bool importParticles() override { + return GeneratorPythia8::importParticles() && makeForcedDecays(); + }; + + void setPDGs(TString pdgs) { + TObjArray *obj = pdgs.Tokenize(";"); + for (int i = 0; i < obj->GetEntriesFast(); i++) { + std::string spdg = obj->At(i)->GetName(); + int pdg = std::stoi(spdg); + mPdgCodes.push_back(pdg); + LOG(info) << "Force decay of PDG = " << pdg; + } + } + +protected: + bool InitDecayer() { + mDecayer = new DecayerPythia8ForceDecays(); + mDecayer->Init(); + mDecayer->calculateWeights(mPdgCodes); + return true; + } + + bool makeForcedDecays() { + int mParticles_size = mParticles.size(); + for (int i = 0; i < mParticles_size; i++) { + auto p = mParticles[i]; + int pdg = p.GetPdgCode(); + if (std::find(mPdgCodes.begin(), mPdgCodes.end(), abs(pdg)) != + mPdgCodes.end()) { + mDecayer->forceDecays(mParticles, i); + mParticles_size = mParticles.size(); + } + } + return true; + } + +private: + DecayerPythia8ForceDecays *mDecayer; + std::vector mPdgCodes; +}; + +} // namespace eventgen +} // namespace o2 + +FairGenerator *GeneratePythia8ForcedDecays(TString pdgs) { + auto gen = new o2::eventgen::GeneratorPythia8ForcedDecays(); + gen->setPDGs(pdgs); + return gen; +} \ No newline at end of file diff --git a/MC/config/PWGEM/ini/Pythia8_Beauty_Cocktail.ini b/MC/config/PWGEM/ini/Pythia8_Beauty_Cocktail.ini new file mode 100644 index 000000000..d49648089 --- /dev/null +++ b/MC/config/PWGEM/ini/Pythia8_Beauty_Cocktail.ini @@ -0,0 +1,12 @@ +[GeneratorExternal] +fileName = ${O2DPG_ROOT}/MC/config/PWGEM/external/generator/Generator_pythia8_forcedDecays.C +funcName=GeneratePythia8ForcedDecays("411;421;431;4122;4232;4132;4332;511;521;531;5122;5132;5232;5332") + +[GeneratorPythia8] +config = ${O2DPG_ROOT}/MC/config/PWGEM/pythia8/generator/pythia8_hf_cocktail.cfg +hooksFileName = ${O2DPG_ROOT}/MC/config/PWGHF/pythia8/hooks/pythia8_userhooks_qqbar.C +hooksFuncName = pythia8_userhooks_bbbar(-9999.,9999.) + +[DecayerPythia8] +config[0] = ${O2DPG_ROOT}/MC/config/common/pythia8/decayer/base.cfg +config[1] = ${O2DPG_ROOT}/MC/config/PWGEM/pythia8/decayer/force_semileptonic.cfg diff --git a/MC/config/PWGEM/ini/Pythia8_Charm_Cocktail.ini b/MC/config/PWGEM/ini/Pythia8_Charm_Cocktail.ini new file mode 100644 index 000000000..fc63c5759 --- /dev/null +++ b/MC/config/PWGEM/ini/Pythia8_Charm_Cocktail.ini @@ -0,0 +1,12 @@ +[GeneratorExternal] +fileName = ${O2DPG_ROOT}/MC/config/PWGEM/external/generator/Generator_pythia8_forcedDecays.C +funcName=GeneratePythia8ForcedDecays("411;421;431;4122;4232;4132;4332") + +[GeneratorPythia8] +config = ${O2DPG_ROOT}/MC/config/PWGEM/pythia8/generator/pythia8_hf_cocktail.cfg +hooksFileName = ${O2DPG_ROOT}/MC/config/PWGHF/pythia8/hooks/pythia8_userhooks_qqbar.C +hooksFuncName = pythia8_userhooks_ccbar(-9999.,9999.) + +[DecayerPythia8] +config[0] = ${O2DPG_ROOT}/MC/config/common/pythia8/decayer/base.cfg +config[1] = ${O2DPG_ROOT}/MC/config/PWGEM/pythia8/decayer/force_semileptonic.cfg diff --git a/MC/config/PWGEM/pythia8/decayer/force_semileptonic.cfg b/MC/config/PWGEM/pythia8/decayer/force_semileptonic.cfg new file mode 100644 index 000000000..c15283caa --- /dev/null +++ b/MC/config/PWGEM/pythia8/decayer/force_semileptonic.cfg @@ -0,0 +1,80 @@ +### only semileptonic decays for charm +### D+ +411:oneChannel = 1 0.087 0 -311 -11 12 +411:addChannel = 1 0.040 0 -321 211 -11 12 +411:addChannel = 1 0.037 0 -313 -11 12 +### D0 +421:oneChannel = 1 0.035 0 -321 -11 12 +421:addChannel = 1 0.022 0 -323 -11 12 +421:addChannel = 1 0.016 0 -321 111 -11 12 +421:addChannel = 1 0.014 0 -311 -211 -11 12 +### Ds +431:oneChannel = 1 0.025 0 333 -11 12 +431:addChannel = 1 0.027 0 221 -11 12 +### Lambdac +4122:oneChannel = 1 0.036 0 3122 -11 12 +### chi_{c}^{+} +4232:oneChannel = 1 0.07 0 3322 -11 12 +### chi_{c}^{0} +4132:oneChannel = 1 0.014 0 3312 -11 12 +### Omega_{c} +4332:oneChannel = 1 0.01224 0 3334 -11 12 + +### only semileptonic decays for beauty +### B0 +511:oneChannel = 1 0.0207000 0 12 -11 -411 +511:addChannel = 1 0.0570000 0 12 -11 -413 +511:addChannel = 1 0.0023000 0 12 -11 -415 +511:addChannel = 1 0.0001330 0 12 -11 -211 +511:addChannel = 1 0.0002690 0 12 -11 -213 +511:addChannel = 1 0.0045000 0 12 -11 -10411 +511:addChannel = 1 0.0052000 0 12 -11 -10413 +511:addChannel = 1 0.0083000 0 12 -11 -20413 + +### B+ +521:oneChannel = 1 0.0000720 0 12 -11 111 +521:addChannel = 1 0.0001450 0 12 -11 113 +521:addChannel = 1 0.0000840 0 12 -11 221 +521:addChannel = 1 0.0001450 0 12 -11 223 +521:addChannel = 1 0.0000840 0 12 -11 331 +521:addChannel = 1 0.0224000 0 12 -11 -421 +521:addChannel = 1 0.0617000 0 12 -11 -423 +521:addChannel = 1 0.0030000 0 12 -11 -425 +521:addChannel = 1 0.0049000 0 12 -11 -10421 +521:addChannel = 1 0.0056000 0 12 -11 -10423 +521:addChannel = 1 0.0090000 0 12 -11 -20423 + +### Bs +531:oneChannel = 1 0.0002000 0 12 -11 -321 +531:addChannel = 1 0.0003000 0 12 -11 -323 +531:addChannel = 1 0.0210000 0 12 -11 -431 +531:addChannel = 1 0.0490000 0 12 -11 -433 +531:addChannel = 1 0.0070000 0 12 -11 -435 +531:addChannel = 1 0.0003000 0 12 -11 -10323 +531:addChannel = 1 0.0040000 0 12 -11 -10431 +531:addChannel = 1 0.0070000 0 12 -11 -10433 +531:addChannel = 1 0.0002000 0 12 -11 -20323 +531:addChannel = 1 0.0040000 0 12 -11 -20433 + +### Lambdab +5122:oneChannel = 1 0.0546000 0 -12 11 4122 +5122:addChannel = 1 0.0096000 0 -12 11 4124 +5122:addChannel = 1 0.0128000 0 -12 11 14122 + +### Chi_{b}^{-} +5132:oneChannel = 1 0.1080010 0 -12 11 4 3101 +5132:addChannel = 1 0.0020000 0 -12 11 2 3101 +### Chi_{b}^{0} +5232:oneChannel = 1 0.1080010 0 -12 11 4 3201 +5232:addChannel = 1 0.0020000 0 -12 11 2 3201 +### Omega_{b}^{-} +5332:oneChannel = 1 0.1080010 1 -12 11 4 3303 +5332:oneChannel = 1 0.0020000 1 -12 11 2 3303 + + +# Correct OmegaC decay length (wrong in PYTHIA8 decay table) (mm/c) +4332:tau0 = 0.08000000000 +# Correct Lb decay length (wrong in PYTHIA8 decay table) +5122:tau0 = 4.41000e-01 + + diff --git a/MC/config/PWGEM/pythia8/generator/pythia8_hf_cocktail.cfg b/MC/config/PWGEM/pythia8/generator/pythia8_hf_cocktail.cfg new file mode 100644 index 000000000..2f4c6c14c --- /dev/null +++ b/MC/config/PWGEM/pythia8/generator/pythia8_hf_cocktail.cfg @@ -0,0 +1,48 @@ +### beams +Beams:idA 2212 # proton +Beams:idB 2212 # proton +Beams:eCM 13600. # GeV + +### processes +SoftQCD:inelastic on + +### decays +ParticleDecays:limitTau0 on +ParticleDecays:tau0Max 10. + +### switch on color reconnection in mode 2 (https://arxiv.org/pdf/1505.01681.pdf) +Tune:pp = 14 +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 + +### switch off decays +411:mayDecay off # D+ +421:mayDecay off # D0 +431:mayDecay off # D_s +4122:mayDecay off # Lambda_c +4232:mayDecay off # Xi_c +4132:mayDecay off # Xi_c_0 +4332:mayDecay off # Omega_c +511:mayDecay off # B0 +521:mayDecay off # B+ +531:mayDecay off # B_s0 +541:mayDecay off # B_c+ +5112:mayDecay off # Sigma_b- +5122:mayDecay off # Lambda_b0 +5132:mayDecay off # Xi_b- +5232:mayDecay off # Xi_b0 +5332:mayDecay off # Omega_b- \ No newline at end of file diff --git a/MC/config/PWGEM/pythia8/hooks/pythia8_userhooks_HF.C b/MC/config/PWGEM/pythia8/hooks/pythia8_userhooks_HF.C deleted file mode 100644 index 83a1930cf..000000000 --- a/MC/config/PWGEM/pythia8/hooks/pythia8_userhooks_HF.C +++ /dev/null @@ -1,63 +0,0 @@ -#include "Pythia8/Pythia.h" - -class UserHooks_HF : public Pythia8::UserHooks -{ - - public: - UserHooks_HF() = default; - ~UserHooks_HF() = default; - bool canVetoPartonLevel() override { return true; }; - bool doVetoPartonLevel(const Pythia8::Event& event) override - { - for (int ipa = 0; ipa < event.size(); ++ipa) { - auto daughterList = event[ipa].daughterList(); - bool hasQ1 = false; - bool hasQ1bar = false; - bool hasQ2 = false; - bool hasQ2bar = false; - for (auto ida : daughterList) { - if (event[ida].id() == mPDG1) - hasQ1 = true; - if (event[ida].id() == -mPDG1) - hasQ1bar = true; - if (event[ida].id() == mPDG2) - hasQ2 = true; - if (event[ida].id() == -mPDG2) - hasQ2bar = true; - } - if ( (hasQ1 && hasQ1bar) || (hasQ2 && hasQ2bar) ) - return false; // do not veto event - } - return true; // veto event - }; - - void setPDG(int pdg1, int pdg2) { mPDG1 = pdg1; mPDG2 = pdg2 }; - - private: - int mPDG1 = 4; - int mPDG2 = 5; -}; - -Pythia8::UserHooks* - pythia8_userhooks_ccbar() -{ - auto hooks = new UserHooks_HF(); - hooks->setPDG(4, 4); - return hooks; -} - -Pythia8::UserHooks* - pythia8_userhooks_bbbar() -{ - auto hooks = new UserHooks_HF(); - hooks->setPDG(5, 5); - return hooks; -} - -Pythia8::UserHooks* - pythia8_userhooks_ccbar_or_bbbar() -{ - auto hooks = new UserHooks_HF(); - hooks->setPDG(4, 5); - return hooks; -}