From e9c834b56c2157e32c5efda5b2966e903e72c845 Mon Sep 17 00:00:00 2001 From: Fabrizio Grosa Date: Thu, 28 May 2026 18:33:17 +0200 Subject: [PATCH 1/4] Add configuration for c-deuteron simulation --- .../generator_pythia8_embed_charmnuclei.C | 208 ++++++++++++++++++ .../ini/GeneratorHF_CDeuteron_Injected.ini | 9 + .../generator/pythia8_cdeuteron_to_dkpi.cfg | 25 +++ 3 files changed, 242 insertions(+) create mode 100644 MC/config/PWGHF/external/generator/generator_pythia8_embed_charmnuclei.C create mode 100644 MC/config/PWGHF/ini/GeneratorHF_CDeuteron_Injected.ini create mode 100644 MC/config/PWGHF/pythia8/generator/pythia8_cdeuteron_to_dkpi.cfg diff --git a/MC/config/PWGHF/external/generator/generator_pythia8_embed_charmnuclei.C b/MC/config/PWGHF/external/generator/generator_pythia8_embed_charmnuclei.C new file mode 100644 index 000000000..92315af80 --- /dev/null +++ b/MC/config/PWGHF/external/generator/generator_pythia8_embed_charmnuclei.C @@ -0,0 +1,208 @@ +#include "FairGenerator.h" +#include "Generators/GeneratorPythia8.h" +#include "Pythia8/Pythia.h" +#include "TRandom.h" +#include "TF1.h" +#include "TMath.h" +#include +#include +#include +#include + +R__ADD_INCLUDE_PATH($O2DPG_MC_CONFIG_ROOT) +#include "MC/config/common/external/generator/CoalescencePythia8.h" + +using namespace Pythia8; + +class GeneratorPythia8HFEmbedCharmNuclei : public o2::eventgen::GeneratorPythia8 +{ + public: + /// default constructor + GeneratorPythia8HFEmbedCharmNuclei() = default; + + /// constructor + GeneratorPythia8HFEmbedCharmNuclei(int pdgCode = 2010010020, float lifetime = 1.f, int nCharmNucleiPerEvent = 10, float yMin = -1.f, float yMax = 1.f, bool trivialCoal = false, float coalMomentum = 0.2) + { + nNumberOfCharmNucleiPerEvent = nCharmNucleiPerEvent; + mRapidityMinCharmNuclei = yMin; + mRapidityMaxCharmNuclei = yMax; + mTrivialCoal = trivialCoal; + mCoalMomentum = coalMomentum; + mPdgCharmNucleus = pdgCode; + if (std::abs(mPdgCharmNucleus) == 2010010020) { + mMassCharmNucleus = 3.226f; + } else { + LOG(fatal) << "********** [GeneratorPythia8HFEmbedCharmNuclei] Only c-deuteron (pdg=2010010020) currently supported! Exit **********"; + } + mLifetimeCharmNucleus = lifetime; + mDecayDistr = TF1("mDecayDistr", Form("exp(-x/%f)", mLifetimeCharmNucleus), 0., 1000.) + + Print(); + + /** switch off process level **/ + mPythiaGun.readString("ProcessLevel:all off"); + auto& param = o2::eventgen::GeneratorPythia8::Instance(); + LOG(info) << "Init \'GeneratorPythia8HFEmbedCharmNuclei\' with following parameters"; + LOG(info) << param; + for (int iCfg{0}; iCfg < 8; ++iCfg) { + if (param.config[iCfg].empty()) { + continue; + } + std::string config = gSystem->ExpandPathName(param.config[iCfg].c_str()); + LOG(info) << "GeneratorPythia8HFEmbedCharmNuclei Reading configuration from file: " << config; + if (!mPythiaGun.readFile(config, true)) { + LOG(fatal) << "Failed to init \'GeneratorPythia8\': problems with configuration file " + << config; + return; + } + } + if (!mPythiaGun.init()) { + LOG(fatal) << "Failed to init \'GeneratorPythia8\': init returned with error"; + return; + } + } + + /// Destructor + ~GeneratorPythia8HFEmbedCharmNuclei() = default; + + /// Init + bool Init() override + { + return o2::eventgen::GeneratorPythia8::Init(); + } + + /// Print the input + void Print() + { + LOG(info) << "********** GeneratorPythia8HFEmbedCharmNuclei configuration dump **********"; + LOG(info) << Form("* PDG code of charm nuclei to be injected: %d", mPdgCharmNucleus); + LOG(info) << Form("* Mass of charm nuclei to be injected (GeV/c2): %f", mMassCharmNucleus); + LOG(info) << Form("* Lifetime of charm nuclei to be injected (mm): %f", mLifetimeCharmNucleus); + LOG(info) << Form("* Number of charm nuclei injected per event: %d", nNumberOfCharmNucleiPerEvent); + LOG(info) << Form("* Hadron rapidity: %f - %f", mHadRapidityMin, mHadRapidityMax); + LOG(info) << Form("* Trivial coalescence: %d", mTrivialCoal); + LOG(info) << Form("* Coalescence momentum: %f", mCoalMomentum); + LOG(info) << "***********************************************************************"; + } + + void setHadronRapidity(float yMin, float yMax) + { + mRapidityMinCharmNuclei = yMin; + mRapidityMaxCharmNuclei = yMax; + }; + + void setUsedSeed(unsigned int seed) + { + mUsedSeed = seed; + }; + + unsigned int getUsedSeed() const + { + return mUsedSeed; + }; + + protected: + //__________________________________________________________________ + bool generateEvent() override + { + // we start from an empty event + mPythia.event.reset(); + + // we simulate c-deuteron decays + for (int iCharmNuclei{0}; iCharmNucleiRndm() < 0.5) { + sign = -1; + } + } + + auto pt = gRandom->Uniform(0., 50.); // placeholder, to be modified + auto y = gRandom->Uniform(mHadRapidityMin, mHadRapidityMax); + auto phi = gRandom->Uniform(0, TMath::TwoPi()); + auto px = pt * TMath::Cos(phi); + auto py = pt * TMath::Sin(phi); + auto mt = TMath::Sqrt(mMassCharmNucleus * mMassCharmNucleus + pt * pt); + auto pz = mt * TMath::SinH(y); + auto p = TMath::Sqrt(pt * pt + pz * pz); + auto e = TMath::Sqrt(mass * mass + p * p); + + Particle particle; + particle.id(sign * mPdgCharmNucleus); + particle.status(81); + particle.m(mMassCharmNucleus); + particle.px(px); + particle.py(py); + particle.pz(pz); + particle.e(e); + particle.xProd(0.f); + particle.yProd(0.f); + particle.zProd(0.f); + particle.tau(mDecayDistr->GetRandom()); + mPythiaGun.particleData.mayDecay(mPdgCharmNucleus, true); // force decay + + bool isCoalSuccess{false}; + while(!isCoalSuccess) { + mPythiaGun.event.reset(); + mPythiaGun.event.append(particle); + mPythiaGun.moreDecays(); + std::array dausToCoal = {-1, -1}; + int iDau{-1}; + for (int iPart{0}; iPart{1000010020}, mTrivialCoal, mCoalMomentum, dausToCoal[0], dausToCoal[1]); + if (isCoalSuccess) { + for (int iPart{0}; iPartTRandom::GetSeed() % 900000000); + myGen->readString("Random:setSeed on"); + myGen->readString("Random:seed " + std::to_string(seed)); + return myGen; +} diff --git a/MC/config/PWGHF/ini/GeneratorHF_CDeuteron_Injected.ini b/MC/config/PWGHF/ini/GeneratorHF_CDeuteron_Injected.ini new file mode 100644 index 000000000..3840c9076 --- /dev/null +++ b/MC/config/PWGHF/ini/GeneratorHF_CDeuteron_Injected.ini @@ -0,0 +1,9 @@ +#NEV_TEST> 100 +### The external generator derives from GeneratorPythia8. +[GeneratorExternal] +fileName=${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGHF/external/generator/generator_pythia8_embed_charmnuclei.C +funcName=GenerateHFEmbedCDeuteron(1.,10,-1.,1.,0,0.2) + +[GeneratorPythia8] +config=${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGHF/pythia8/generator/pythia8_cdeuteron_to_dkpi.cfg +includePartonEvent=false diff --git a/MC/config/PWGHF/pythia8/generator/pythia8_cdeuteron_to_dkpi.cfg b/MC/config/PWGHF/pythia8/generator/pythia8_cdeuteron_to_dkpi.cfg new file mode 100644 index 000000000..35ff72709 --- /dev/null +++ b/MC/config/PWGHF/pythia8/generator/pythia8_cdeuteron_to_dkpi.cfg @@ -0,0 +1,25 @@ +### author: Fabrizio Grosa (fabrizio.grosa@cern.ch) +### since: May 2026 + +### beams (not relevant) +Beams:idA 2212 # proton +Beams:idB 2212 # proton +Beams:eCM 13600. # GeV + +### processes +ProcessLevel:all off # turning off all the processes, we do not have an underlying event + +### decays +ParticleDecays:limitTau0 on +ParticleDecays:tau0Max 10. + +### add the c-deuteron +### id:all = name antiName spinType chargeType colType m0 mWidth mMin mMax tau0 +2010010020:all = CDeuteron+ AntiCDeuteron- 1 3 0 3.226 0 3.226 3.226 1. +2010010020:mayDecay = on + +# same as Λc+ -> p K- π+ including resonances + a neutron +2010010020:oneChannel = 1 0.14400 0 2112 2212 -321 211 ### Λc+ -> p K- π+ (non-resonant) 3.5% +2010010020:addChannel = 1 0.08100 100 2112 2212 -313 ### Λc+ -> p K*0(892) 1.96% +2010010020:addChannel = 1 0.04500 100 2112 2224 -321 ### Λc+ -> Delta++ K- 1.08% +2010010020:addChannel = 1 0.09000 100 2112 102134 211 ### Λc+ -> Lambda(1520) K- 2.20e-3 From ab26c148b4fb104b34fa1c8b8f8a82362582552b Mon Sep 17 00:00:00 2001 From: Fabrizio Grosa Date: Fri, 29 May 2026 23:52:23 +0200 Subject: [PATCH 2/4] Fix coalescence and definition of event containing decayed charmed nuclei --- .../generator_pythia8_embed_charmnuclei.C | 131 +++++++++++------- .../ini/GeneratorHF_CDeuteron_Injected.ini | 5 +- .../generator/pythia8_cdeuteron_to_dkpi.cfg | 14 +- 3 files changed, 94 insertions(+), 56 deletions(-) diff --git a/MC/config/PWGHF/external/generator/generator_pythia8_embed_charmnuclei.C b/MC/config/PWGHF/external/generator/generator_pythia8_embed_charmnuclei.C index 92315af80..0d4d8e21a 100644 --- a/MC/config/PWGHF/external/generator/generator_pythia8_embed_charmnuclei.C +++ b/MC/config/PWGHF/external/generator/generator_pythia8_embed_charmnuclei.C @@ -17,8 +17,6 @@ using namespace Pythia8; class GeneratorPythia8HFEmbedCharmNuclei : public o2::eventgen::GeneratorPythia8 { public: - /// default constructor - GeneratorPythia8HFEmbedCharmNuclei() = default; /// constructor GeneratorPythia8HFEmbedCharmNuclei(int pdgCode = 2010010020, float lifetime = 1.f, int nCharmNucleiPerEvent = 10, float yMin = -1.f, float yMax = 1.f, bool trivialCoal = false, float coalMomentum = 0.2) @@ -35,42 +33,32 @@ class GeneratorPythia8HFEmbedCharmNuclei : public o2::eventgen::GeneratorPythia8 LOG(fatal) << "********** [GeneratorPythia8HFEmbedCharmNuclei] Only c-deuteron (pdg=2010010020) currently supported! Exit **********"; } mLifetimeCharmNucleus = lifetime; - mDecayDistr = TF1("mDecayDistr", Form("exp(-x/%f)", mLifetimeCharmNucleus), 0., 1000.) + mDecayDistr = new TF1("mDecayDistr", Form("exp(-x/%f)", mLifetimeCharmNucleus), 0., 1000.); Print(); /** switch off process level **/ - mPythiaGun.readString("ProcessLevel:all off"); - auto& param = o2::eventgen::GeneratorPythia8::Instance(); + // mPythiaGun.readString("ProcessLevel:all off"); + auto& param = o2::eventgen::GeneratorPythia8Param::Instance(); LOG(info) << "Init \'GeneratorPythia8HFEmbedCharmNuclei\' with following parameters"; LOG(info) << param; - for (int iCfg{0}; iCfg < 8; ++iCfg) { - if (param.config[iCfg].empty()) { - continue; - } - std::string config = gSystem->ExpandPathName(param.config[iCfg].c_str()); - LOG(info) << "GeneratorPythia8HFEmbedCharmNuclei Reading configuration from file: " << config; - if (!mPythiaGun.readFile(config, true)) { - LOG(fatal) << "Failed to init \'GeneratorPythia8\': problems with configuration file " - << config; - return; - } + if (param.config.empty()) { + LOG(fatal) << "Failed to init \'GeneratorPythia8\': problems with configuration file "; + } + std::string cfg = gSystem->ExpandPathName(param.config.c_str()); + LOG(info) << "GeneratorPythia8HFEmbedCharmNuclei Reading configuration from file: " << cfg; + if (!mPythiaGun.readFile(cfg, true)) { + LOG(fatal) << "Failed to init \'GeneratorPythia8\': problems with configuration file " << cfg; } + if (!mPythiaGun.init()) { LOG(fatal) << "Failed to init \'GeneratorPythia8\': init returned with error"; - return; } } /// Destructor ~GeneratorPythia8HFEmbedCharmNuclei() = default; - /// Init - bool Init() override - { - return o2::eventgen::GeneratorPythia8::Init(); - } - /// Print the input void Print() { @@ -79,7 +67,7 @@ class GeneratorPythia8HFEmbedCharmNuclei : public o2::eventgen::GeneratorPythia8 LOG(info) << Form("* Mass of charm nuclei to be injected (GeV/c2): %f", mMassCharmNucleus); LOG(info) << Form("* Lifetime of charm nuclei to be injected (mm): %f", mLifetimeCharmNucleus); LOG(info) << Form("* Number of charm nuclei injected per event: %d", nNumberOfCharmNucleiPerEvent); - LOG(info) << Form("* Hadron rapidity: %f - %f", mHadRapidityMin, mHadRapidityMax); + LOG(info) << Form("* Hadron rapidity: %f - %f", mRapidityMinCharmNuclei, mRapidityMaxCharmNuclei); LOG(info) << Form("* Trivial coalescence: %d", mTrivialCoal); LOG(info) << Form("* Coalescence momentum: %f", mCoalMomentum); LOG(info) << "***********************************************************************"; @@ -101,7 +89,6 @@ class GeneratorPythia8HFEmbedCharmNuclei : public o2::eventgen::GeneratorPythia8 return mUsedSeed; }; - protected: //__________________________________________________________________ bool generateEvent() override { @@ -123,18 +110,18 @@ class GeneratorPythia8HFEmbedCharmNuclei : public o2::eventgen::GeneratorPythia8 } auto pt = gRandom->Uniform(0., 50.); // placeholder, to be modified - auto y = gRandom->Uniform(mHadRapidityMin, mHadRapidityMax); + auto y = gRandom->Uniform(mRapidityMinCharmNuclei, mRapidityMaxCharmNuclei); auto phi = gRandom->Uniform(0, TMath::TwoPi()); auto px = pt * TMath::Cos(phi); auto py = pt * TMath::Sin(phi); auto mt = TMath::Sqrt(mMassCharmNucleus * mMassCharmNucleus + pt * pt); auto pz = mt * TMath::SinH(y); auto p = TMath::Sqrt(pt * pt + pz * pz); - auto e = TMath::Sqrt(mass * mass + p * p); + auto e = TMath::Sqrt(mMassCharmNucleus * mMassCharmNucleus + p * p); Particle particle; particle.id(sign * mPdgCharmNucleus); - particle.status(81); + particle.status(83); particle.m(mMassCharmNucleus); particle.px(px); particle.py(py); @@ -143,7 +130,7 @@ class GeneratorPythia8HFEmbedCharmNuclei : public o2::eventgen::GeneratorPythia8 particle.xProd(0.f); particle.yProd(0.f); particle.zProd(0.f); - particle.tau(mDecayDistr->GetRandom()); + particle.tau(0.f); //mDecayDistr->GetRandom()); mPythiaGun.particleData.mayDecay(mPdgCharmNucleus, true); // force decay bool isCoalSuccess{false}; @@ -152,48 +139,88 @@ class GeneratorPythia8HFEmbedCharmNuclei : public o2::eventgen::GeneratorPythia8 mPythiaGun.event.append(particle); mPythiaGun.moreDecays(); std::array dausToCoal = {-1, -1}; + std::vector pdgShortLivedResos = {313, 2224, 102134}; + bool isResoFound{false}; + int idxCharmNucleus{-1}; + for (int iPart{0}; iPart{1000010020}, mTrivialCoal, mCoalMomentum, dausToCoal[0], dausToCoal[1]); + isCoalSuccess = CoalescencePythia8(mPythiaGun.event, std::vector{1000010020}, mTrivialCoal, mCoalMomentum, dausToCoal[0], dausToCoal[1]); if (isCoalSuccess) { + int offset = mPythia.event.size(); // we need to rescale the indices of mothers and daughters, accounting for the particles that are already appended to the event for (int iPart{0}; iPart 0) { + part.mother1(mother1 + offset); + } + if (mother2 > 0) { + part.mother2(mother2 + offset); + } + if (daughter1 > 0) { + part.daughter1(daughter1 + offset); + } + if (daughter2 > 0) { + part.daughter2(daughter2 + offset); + } + mPythia.event.append(part); } } - mPythiaGun.next(); } } - mPythia.next(); - return true; } private: // Properties of selection - float mMassCharmNucleus; - int mPdgCharmNucleus; - float mLifetimeCharmNucleus; - int nNumberOfCharmNucleiPerEvent; - float mRapidityMinCharmNuclei; - float mRapidityMaxCharmNuclei; - unsigned int mUsedSeed; - - bool mTrivialCoal = false; /// if true, the coalescence is done without checking the distance in the phase space of the nucleons - float mCoalMomentum; /// coalescence momentum - - Pythia8::Pythia mPythiaGun; // Gun generator with decay support - - TF1* mDecayDistr; + float mMassCharmNucleus; /// mass of the charmed nucleus + int mPdgCharmNucleus; /// pdg code of the charmed nucleus + float mLifetimeCharmNucleus; /// lifetime of the charmed nucleus + int nNumberOfCharmNucleiPerEvent; /// number of charmed nuclei injected per event + float mRapidityMinCharmNuclei; /// rapidity min + float mRapidityMaxCharmNuclei; /// rapidity max + unsigned int mUsedSeed; /// seed + bool mTrivialCoal; /// if true, the coalescence is done without checking the distance in the phase space of the nucleons + float mCoalMomentum; /// coalescence momentum + Pythia8::Pythia mPythiaGun; /// Gun generator with decay support + TF1* mDecayDistr; /// Lifetime distribution }; diff --git a/MC/config/PWGHF/ini/GeneratorHF_CDeuteron_Injected.ini b/MC/config/PWGHF/ini/GeneratorHF_CDeuteron_Injected.ini index 3840c9076..1da458dc8 100644 --- a/MC/config/PWGHF/ini/GeneratorHF_CDeuteron_Injected.ini +++ b/MC/config/PWGHF/ini/GeneratorHF_CDeuteron_Injected.ini @@ -2,8 +2,9 @@ ### The external generator derives from GeneratorPythia8. [GeneratorExternal] fileName=${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGHF/external/generator/generator_pythia8_embed_charmnuclei.C -funcName=GenerateHFEmbedCDeuteron(1.,10,-1.,1.,0,0.2) +funcName=GenerateHFEmbedCDeuteron(1., 10, -1., 1., 0, 0.4) +### includePartonEvent is needed to keep the c-deuteron in the event record, even if there are no partons in the event [GeneratorPythia8] config=${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGHF/pythia8/generator/pythia8_cdeuteron_to_dkpi.cfg -includePartonEvent=false +includePartonEvent=true diff --git a/MC/config/PWGHF/pythia8/generator/pythia8_cdeuteron_to_dkpi.cfg b/MC/config/PWGHF/pythia8/generator/pythia8_cdeuteron_to_dkpi.cfg index 35ff72709..ddcf62126 100644 --- a/MC/config/PWGHF/pythia8/generator/pythia8_cdeuteron_to_dkpi.cfg +++ b/MC/config/PWGHF/pythia8/generator/pythia8_cdeuteron_to_dkpi.cfg @@ -7,7 +7,7 @@ Beams:idB 2212 # proton Beams:eCM 13600. # GeV ### processes -ProcessLevel:all off # turning off all the processes, we do not have an underlying event +SoftQCD:inelastic on # all inelastic processes ### decays ParticleDecays:limitTau0 on @@ -15,7 +15,7 @@ ParticleDecays:tau0Max 10. ### add the c-deuteron ### id:all = name antiName spinType chargeType colType m0 mWidth mMin mMax tau0 -2010010020:all = CDeuteron+ AntiCDeuteron- 1 3 0 3.226 0 3.226 3.226 1. +2010010020:all = CDeuteron AntiCDeuteron 1 3 0 3.226 0 3.226 3.226 1. 2010010020:mayDecay = on # same as Λc+ -> p K- π+ including resonances + a neutron @@ -23,3 +23,13 @@ ParticleDecays:tau0Max 10. 2010010020:addChannel = 1 0.08100 100 2112 2212 -313 ### Λc+ -> p K*0(892) 1.96% 2010010020:addChannel = 1 0.04500 100 2112 2224 -321 ### Λc+ -> Delta++ K- 1.08% 2010010020:addChannel = 1 0.09000 100 2112 102134 211 ### Λc+ -> Lambda(1520) K- 2.20e-3 + +### K*0(892) -> K- π+ +313:onMode = off +313:onIfAll = 321 211 +### for Λc -> Delta++ K- +2224:onMode = off +2224:onIfAll = 2212 211 +### for Λc -> Lambda(1520) K- +102134:onMode = off +102134:onIfAll = 2212 321 From e71700c71e28e77233683579da4308b922f4cdc4 Mon Sep 17 00:00:00 2001 From: Fabrizio Grosa Date: Sun, 31 May 2026 22:01:55 +0200 Subject: [PATCH 3/4] Add possibility to simulate also non-prompt c-deuterons --- .../generator_pythia8_embed_charmnuclei.C | 127 +++++++++++++----- .../ini/GeneratorHF_CDeuteron_Injected.ini | 2 +- .../generator/pythia8_cdeuteron_to_dkpi.cfg | 75 ++++++++++- .../external/generator/CoalescencePythia8.h | 4 +- 4 files changed, 169 insertions(+), 39 deletions(-) diff --git a/MC/config/PWGHF/external/generator/generator_pythia8_embed_charmnuclei.C b/MC/config/PWGHF/external/generator/generator_pythia8_embed_charmnuclei.C index 0d4d8e21a..da4224489 100644 --- a/MC/config/PWGHF/external/generator/generator_pythia8_embed_charmnuclei.C +++ b/MC/config/PWGHF/external/generator/generator_pythia8_embed_charmnuclei.C @@ -19,14 +19,17 @@ class GeneratorPythia8HFEmbedCharmNuclei : public o2::eventgen::GeneratorPythia8 public: /// constructor - GeneratorPythia8HFEmbedCharmNuclei(int pdgCode = 2010010020, float lifetime = 1.f, int nCharmNucleiPerEvent = 10, float yMin = -1.f, float yMax = 1.f, bool trivialCoal = false, float coalMomentum = 0.2) + GeneratorPythia8HFEmbedCharmNuclei(int pdgCode = 2010010020, float lifetime = 1.f, int nCharmNucleiPerEvent = 10, float yMin = -1.f, float yMax = 1.f, float ptMax = 25.f, bool trivialCoal = false, float coalMomentum = 0.2, float fracFromB = 0.f) { nNumberOfCharmNucleiPerEvent = nCharmNucleiPerEvent; mRapidityMinCharmNuclei = yMin; mRapidityMaxCharmNuclei = yMax; + mPtMaxCharmNuclei = ptMax; mTrivialCoal = trivialCoal; mCoalMomentum = coalMomentum; + mFractionFromBeauty = fracFromB; mPdgCharmNucleus = pdgCode; + mSign = 1; if (std::abs(mPdgCharmNucleus) == 2010010020) { mMassCharmNucleus = 3.226f; } else { @@ -34,11 +37,12 @@ class GeneratorPythia8HFEmbedCharmNuclei : public o2::eventgen::GeneratorPythia8 } mLifetimeCharmNucleus = lifetime; mDecayDistr = new TF1("mDecayDistr", Form("exp(-x/%f)", mLifetimeCharmNucleus), 0., 1000.); + mDecayDistrLb = new TF1("mDecayDistrLb", "exp(-x/0.440)", 0., 1000.); + mPtDistrLb = new TF1("mPtDistrLb","[0]*x/TMath::Power((1+TMath::Power(x/[1],[3])),[2])",0.,1000.); + mPtDistrLb->SetParameters(1000., 6.97355, 3.20721, 1.71678); Print(); - /** switch off process level **/ - // mPythiaGun.readString("ProcessLevel:all off"); auto& param = o2::eventgen::GeneratorPythia8Param::Instance(); LOG(info) << "Init \'GeneratorPythia8HFEmbedCharmNuclei\' with following parameters"; LOG(info) << param; @@ -67,9 +71,11 @@ class GeneratorPythia8HFEmbedCharmNuclei : public o2::eventgen::GeneratorPythia8 LOG(info) << Form("* Mass of charm nuclei to be injected (GeV/c2): %f", mMassCharmNucleus); LOG(info) << Form("* Lifetime of charm nuclei to be injected (mm): %f", mLifetimeCharmNucleus); LOG(info) << Form("* Number of charm nuclei injected per event: %d", nNumberOfCharmNucleiPerEvent); - LOG(info) << Form("* Hadron rapidity: %f - %f", mRapidityMinCharmNuclei, mRapidityMaxCharmNuclei); + LOG(info) << Form("* Charmed nucleus rapidity: %f - %f", mRapidityMinCharmNuclei, mRapidityMaxCharmNuclei); + LOG(info) << Form("* Charmed nucleus pT max (prompt): %f", mPtMaxCharmNuclei); LOG(info) << Form("* Trivial coalescence: %d", mTrivialCoal); LOG(info) << Form("* Coalescence momentum: %f", mCoalMomentum); + LOG(info) << Form("* Fraction from beauty: %f", mFractionFromBeauty); LOG(info) << "***********************************************************************"; } @@ -97,32 +103,44 @@ class GeneratorPythia8HFEmbedCharmNuclei : public o2::eventgen::GeneratorPythia8 // we simulate c-deuteron decays for (int iCharmNuclei{0}; iCharmNuclei 0) ? mSign = -1 : mSign = 1; + if (nNumberOfCharmNucleiPerEvent % 2 != 0 && iCharmNuclei == nNumberOfCharmNucleiPerEvent - 1) { if (gRandom->Rndm() < 0.5) { - sign = -1; + mSign = 1; } } - auto pt = gRandom->Uniform(0., 50.); // placeholder, to be modified - auto y = gRandom->Uniform(mRapidityMinCharmNuclei, mRapidityMaxCharmNuclei); + int pdgToGen = mPdgCharmNucleus; + float massToGen = mMassCharmNucleus; + float lifetimeToGen = mDecayDistr->GetRandom(); + float minRapToGen = mRapidityMinCharmNuclei; + float maxRapToGen = mRapidityMaxCharmNuclei; + bool isFromB = gRandom->Rndm() < mFractionFromBeauty; + // we determine if it's prompt or non-prompt + if (isFromB) { + pdgToGen = 5122; // we generate a Lambda_b and we let it decay into the charmed nucleus, no other beauty hadrons are considered + massToGen = 5.61940f; // mass of Lambda_b (GeV/c2) + lifetimeToGen = mDecayDistrLb->GetRandom(); + minRapToGen *= 2; + maxRapToGen *= 2; + } + + auto pt = (!isFromB) ? gRandom->Uniform(0., mPtMaxCharmNuclei) : mPtDistrLb->GetRandom(); + auto y = gRandom->Uniform(minRapToGen, maxRapToGen); auto phi = gRandom->Uniform(0, TMath::TwoPi()); auto px = pt * TMath::Cos(phi); auto py = pt * TMath::Sin(phi); - auto mt = TMath::Sqrt(mMassCharmNucleus * mMassCharmNucleus + pt * pt); + auto mt = TMath::Sqrt(massToGen * massToGen + pt * pt); auto pz = mt * TMath::SinH(y); auto p = TMath::Sqrt(pt * pt + pz * pz); - auto e = TMath::Sqrt(mMassCharmNucleus * mMassCharmNucleus + p * p); + auto e = TMath::Sqrt(massToGen * massToGen + p * p); Particle particle; - particle.id(sign * mPdgCharmNucleus); + particle.id(mSign * pdgToGen); particle.status(83); - particle.m(mMassCharmNucleus); + particle.m(massToGen); particle.px(px); particle.py(py); particle.pz(pz); @@ -130,11 +148,13 @@ class GeneratorPythia8HFEmbedCharmNuclei : public o2::eventgen::GeneratorPythia8 particle.xProd(0.f); particle.yProd(0.f); particle.zProd(0.f); - particle.tau(0.f); //mDecayDistr->GetRandom()); + particle.tau(lifetimeToGen); + mPythiaGun.particleData.mayDecay(5122, true); // force decay mPythiaGun.particleData.mayDecay(mPdgCharmNucleus, true); // force decay bool isCoalSuccess{false}; - while(!isCoalSuccess) { + int nTrials{0}; + while(!isCoalSuccess || nTrials > 1e4) { mPythiaGun.event.reset(); mPythiaGun.event.append(particle); mPythiaGun.moreDecays(); @@ -148,8 +168,9 @@ class GeneratorPythia8HFEmbedCharmNuclei : public o2::eventgen::GeneratorPythia8 if (absPdg == mPdgCharmNucleus) { idxCharmNucleus = iPart; } + auto mother = part.mother1(); // if we find a resonance, we remove it, otherwise we prevent the coalescence of daughters from resonances and daughters from charmed nucleus directly - if (std::find(pdgShortLivedResos.begin(), pdgShortLivedResos.end(), absPdg) != pdgShortLivedResos.end()) { + if (std::find(pdgShortLivedResos.begin(), pdgShortLivedResos.end(), absPdg) != pdgShortLivedResos.end() && mother >= idxCharmNucleus) { // we need to change the indices of the daughter particles to point to the charmed nucleus auto dauList = part.daughterList(); for (auto const& dau : dauList) { @@ -160,22 +181,57 @@ class GeneratorPythia8HFEmbedCharmNuclei : public o2::eventgen::GeneratorPythia8 } } if (isResoFound) { // we have to reset all the particles as daughters of the charm nucleus - mPythiaGun.event[idxCharmNucleus].daughter1(idxCharmNucleus + 1); - mPythiaGun.event[idxCharmNucleus].daughter2(mPythiaGun.event.size() - 1); + std::vector idxDausCharmNucleus{}; + for (int iPart{0}; iPart newPartList{}; + std::vector idxToRemove{}; + for (int iPart{idxDausCharmNucleus[0]}; iPart updatedMothers{}; + for (int iPart{0}; iPart{1000010020}, mTrivialCoal, mCoalMomentum, dausToCoal[0], dausToCoal[1]); + isCoalSuccess = CoalescencePythia8(mPythiaGun.event, std::vector{1000010020}, mTrivialCoal, mCoalMomentum, dausToCoal[0], dausToCoal[1], 10.); if (isCoalSuccess) { int offset = mPythia.event.size(); // we need to rescale the indices of mothers and daughters, accounting for the particles that are already appended to the event for (int iPart{0}; iPartTRandom::GetSeed() % 900000000); myGen->readString("Random:setSeed on"); myGen->readString("Random:seed " + std::to_string(seed)); diff --git a/MC/config/PWGHF/ini/GeneratorHF_CDeuteron_Injected.ini b/MC/config/PWGHF/ini/GeneratorHF_CDeuteron_Injected.ini index 1da458dc8..5a97ca2a4 100644 --- a/MC/config/PWGHF/ini/GeneratorHF_CDeuteron_Injected.ini +++ b/MC/config/PWGHF/ini/GeneratorHF_CDeuteron_Injected.ini @@ -2,7 +2,7 @@ ### The external generator derives from GeneratorPythia8. [GeneratorExternal] fileName=${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGHF/external/generator/generator_pythia8_embed_charmnuclei.C -funcName=GenerateHFEmbedCDeuteron(1., 10, -1., 1., 0, 0.4) +funcName=GenerateHFEmbedCDeuteron(1., 10, -1., 1., 25., 0, 0.4, 0.25) ### includePartonEvent is needed to keep the c-deuteron in the event record, even if there are no partons in the event [GeneratorPythia8] diff --git a/MC/config/PWGHF/pythia8/generator/pythia8_cdeuteron_to_dkpi.cfg b/MC/config/PWGHF/pythia8/generator/pythia8_cdeuteron_to_dkpi.cfg index ddcf62126..9fda88564 100644 --- a/MC/config/PWGHF/pythia8/generator/pythia8_cdeuteron_to_dkpi.cfg +++ b/MC/config/PWGHF/pythia8/generator/pythia8_cdeuteron_to_dkpi.cfg @@ -19,10 +19,10 @@ ParticleDecays:tau0Max 10. 2010010020:mayDecay = on # same as Λc+ -> p K- π+ including resonances + a neutron -2010010020:oneChannel = 1 0.14400 0 2112 2212 -321 211 ### Λc+ -> p K- π+ (non-resonant) 3.5% -2010010020:addChannel = 1 0.08100 100 2112 2212 -313 ### Λc+ -> p K*0(892) 1.96% -2010010020:addChannel = 1 0.04500 100 2112 2224 -321 ### Λc+ -> Delta++ K- 1.08% -2010010020:addChannel = 1 0.09000 100 2112 102134 211 ### Λc+ -> Lambda(1520) K- 2.20e-3 +2010010020:oneChannel = 1 0.14400 0 2112 2212 -321 211 ### cd+ -> n p K- π+ (non-resonant) 3.5% +2010010020:addChannel = 1 0.08100 100 2112 2212 -313 ### cd+ -> n p K*0(892) 1.96% +2010010020:addChannel = 1 0.04500 100 2112 2224 -321 ### cd+ -> n Delta++ K- 1.08% +2010010020:addChannel = 1 0.09000 100 2112 102134 211 ### cd+ -> n Lambda(1520) K- 2.20e-3 ### K*0(892) -> K- π+ 313:onMode = off @@ -33,3 +33,70 @@ ParticleDecays:tau0Max 10. ### for Λc -> Lambda(1520) K- 102134:onMode = off 102134:onIfAll = 2212 321 + +# Λb0 -> cd+ decays taken from PYTHIA +# NB: cd+ must always be the last daughter, not to screw up replacements in the generator +# 15% 2 body +5122:oneChannel = 1 0.15 0 2010010020 -2212 ### Λb0 -> cd+ anti-p +# 27% 3 body +5122:addChannel = 1 0.054 100 -2212 111 2010010020 ### Λb0 -> cd+ anti-p pi0 +5122:addChannel = 1 0.027 100 -2212 113 2010010020 ### Λb0 -> cd+ anti-p rho0 +5122:addChannel = 1 0.027 100 -2212 221 2010010020 ### Λb0 -> cd+ anti-p eta +5122:addChannel = 1 0.027 100 -2212 223 2010010020 ### Λb0 -> cd+ anti-p omega +5122:addChannel = 1 0.060 100 -2112 -211 2010010020 ### Λb0 -> cd+ anti-n0 pi- +5122:addChannel = 1 0.030 100 -2112 -213 2010010020 ### Λb0 -> cd+ anti-n0 rho- +5122:addChannel = 1 0.010 100 -2214 111 2010010020 ### Λb0 -> cd+ DeltaBar- pi0 +5122:addChannel = 1 0.005 100 -2214 113 2010010020 ### Λb0 -> cd+ DeltaBar- rho0 +5122:addChannel = 1 0.005 100 -2214 221 2010010020 ### Λb0 -> cd+ DeltaBar- eta +5122:addChannel = 1 0.005 100 -2214 223 2010010020 ### Λb0 -> cd+ DeltaBar- omega +5122:addChannel = 1 0.007 100 -2114 -211 2010010020 ### Λb0 -> cd+ DeltaBar0 pi- +5122:addChannel = 1 0.003 100 -2114 -213 2010010020 ### Λb0 -> cd+ DeltaBar0 rho- +5122:addChannel = 1 0.007 100 -2224 211 2010010020 ### Λb0 -> cd+ DeltaBar-- pi+ +5122:addChannel = 1 0.003 100 -2224 213 2010010020 ### Λb0 -> cd+ DeltaBar-- rho+ +5122:addChannel = 1 0.010 0 -3122 -321 2010010020 ### Λb0 -> cd+ Lambdabar0 K- +# 31% 4 body +5122:addChannel = 1 0.054 100 -2212 111 111 2010010020 ### Λb0 -> cd+ anti-p pi0 pi0 +5122:addChannel = 1 0.027 100 -2212 211 -211 2010010020 ### Λb0 -> cd+ anti-p pi+ pi- +5122:addChannel = 1 0.013 100 -2212 213 -211 2010010020 ### Λb0 -> cd+ anti-p rho+ pi- +5122:addChannel = 1 0.018 100 -2212 113 111 2010010020 ### Λb0 -> cd+ anti-p rho0 pi0 +5122:addChannel = 1 0.009 100 -2212 113 113 2010010020 ### Λb0 -> cd+ anti-p rho0 rho0 +5122:addChannel = 1 0.018 100 -2212 221 111 2010010020 ### Λb0 -> cd+ anti-p eta pi0 +5122:addChannel = 1 0.009 100 -2212 221 113 2010010020 ### Λb0 -> cd+ anti-p eta rho0 +5122:addChannel = 1 0.018 100 -2212 223 111 2010010020 ### Λb0 -> cd+ anti-p omega pi0 +5122:addChannel = 1 0.009 100 -2212 223 113 2010010020 ### Λb0 -> cd+ anti-p omega rho0 +5122:addChannel = 1 0.040 100 -2112 -211 111 2010010020 ### Λb0 -> cd+ anti-n0 pi- pi0 +5122:addChannel = 1 0.020 100 -2112 -211 113 2010010020 ### Λb0 -> cd+ anti-n0 pi- rho0 +5122:addChannel = 1 0.020 100 -2112 -213 111 2010010020 ### Λb0 -> cd+ anti-n0 rho- pi0 +5122:addChannel = 1 0.010 100 -2112 -213 113 2010010020 ### Λb0 -> cd+ anti-n0 rho- rho0 +5122:addChannel = 1 0.010 100 -2214 111 111 2010010020 ### Λb0 -> cd+ DeltaBar- pi0 pi0 +5122:addChannel = 1 0.005 100 -2214 113 111 2010010020 ### Λb0 -> cd+ DeltaBar- rho0 pi0 +5122:addChannel = 1 0.005 100 -2214 221 111 2010010020 ### Λb0 -> cd+ DeltaBar- eta pi0 +5122:addChannel = 1 0.005 100 -2214 223 111 2010010020 ### Λb0 -> cd+ DeltaBar- omega pi0 +5122:addChannel = 1 0.007 100 -2114 -211 111 2010010020 ### Λb0 -> cd+ DeltaBar0 pi- pi0 +5122:addChannel = 1 0.003 100 -2114 -213 111 2010010020 ### Λb0 -> cd+ DeltaBar0 rho- pi0 +5122:addChannel = 1 0.007 100 -2224 211 111 2010010020 ### Λb0 -> cd+ DeltaBar-- pi+ pi0 +5122:addChannel = 1 0.003 100 -2224 213 111 2010010020 ### Λb0 -> cd+ DeltaBar-- rho+ pi0 +5122:addChannel = 1 0.010 0 -3122 -321 111 2010010020 ### Λb0 -> cd+ Lambdabar0 K- pi0 +# 19% 5 body +5122:addChannel = 1 0.033 100 -2212 111 111 111 2010010020 ### Λb0 -> cd+ anti-p pi0 pi0 pi0 +5122:addChannel = 1 0.016 100 -2212 211 -211 111 2010010020 ### Λb0 -> cd+ anti-p pi+ pi- pi0 +5122:addChannel = 1 0.008 100 -2212 213 -211 111 2010010020 ### Λb0 -> cd+ anti-p rho+ pi- pi0 +5122:addChannel = 1 0.012 100 -2212 113 111 111 2010010020 ### Λb0 -> cd+ anti-p rho0 pi0 pi0 +5122:addChannel = 1 0.006 100 -2212 113 113 111 2010010020 ### Λb0 -> cd+ anti-p rho0 rho0 pi0 +5122:addChannel = 1 0.012 100 -2212 221 111 111 2010010020 ### Λb0 -> cd+ anti-p eta pi0 pi0 +5122:addChannel = 1 0.006 100 -2212 221 113 111 2010010020 ### Λb0 -> cd+ anti-p eta rho0 pi0 +5122:addChannel = 1 0.012 100 -2212 223 111 111 2010010020 ### Λb0 -> cd+ anti-p omega pi0 pi0 +5122:addChannel = 1 0.006 100 -2212 223 113 111 2010010020 ### Λb0 -> cd+ anti-p omega rho0 pi0 +5122:addChannel = 1 0.030 100 -2112 -211 111 111 2010010020 ### Λb0 -> cd+ anti-n0 pi- pi0 pi0 +5122:addChannel = 1 0.011 100 -2112 -211 113 111 2010010020 ### Λb0 -> cd+ anti-n0 pi- rho0 pi0 +5122:addChannel = 1 0.011 100 -2112 -213 111 111 2010010020 ### Λb0 -> cd+ anti-n0 rho- pi0 pi0 +5122:addChannel = 1 0.006 100 -2112 -213 113 111 2010010020 ### Λb0 -> cd+ anti-n0 rho- rho0 pi0 +5122:addChannel = 1 0.006 100 -2214 111 111 111 2010010020 ### Λb0 -> cd+ DeltaBar- pi0 pi0 pi0 +5122:addChannel = 1 0.002 100 -2214 113 111 111 2010010020 ### Λb0 -> cd+ DeltaBar- rho0 pi0 pi0 +5122:addChannel = 1 0.002 100 -2214 221 111 111 2010010020 ### Λb0 -> cd+ DeltaBar- eta pi0 pi0 +5122:addChannel = 1 0.002 100 -2214 223 111 111 2010010020 ### Λb0 -> cd+ DeltaBar- omega pi0 pi0 +5122:addChannel = 1 0.005 100 -2114 -211 111 111 2010010020 ### Λb0 -> cd+ DeltaBar0 pi- pi0 pi0 +5122:addChannel = 1 0.001 100 -2114 -213 111 111 2010010020 ### Λb0 -> cd+ DeltaBar0 rho- pi0 pi0 +5122:addChannel = 1 0.005 100 -2224 211 111 111 2010010020 ### Λb0 -> cd+ DeltaBar-- pi+ pi0 pi0 +5122:addChannel = 1 0.001 100 -2224 213 111 111 2010010020 ### Λb0 -> cd+ DeltaBar-- rho+ pi0 pi0 +5122:addChannel = 1 0.006 0 -3122 -321 111 111 2010010020 ### Λb0 -> cd+ Lambdabar0 K- pi0 pi0 diff --git a/MC/config/common/external/generator/CoalescencePythia8.h b/MC/config/common/external/generator/CoalescencePythia8.h index b1f392ba4..7e4abc53e 100644 --- a/MC/config/common/external/generator/CoalescencePythia8.h +++ b/MC/config/common/external/generator/CoalescencePythia8.h @@ -104,7 +104,7 @@ bool doCoal(Pythia8::Event& event, int charge, int pdgCode, float mass, bool tri return true; } -bool CoalescencePythia8(Pythia8::Event& event, std::vector inputPdgList = {}, bool trivialCoal = false, double coalMomentum = 0.4, int firstDauID = -1, int lastDauId = -1) +bool CoalescencePythia8(Pythia8::Event& event, std::vector inputPdgList = {}, bool trivialCoal = false, double coalMomentum = 0.4, int firstDauID = -1, int lastDauId = -1, float maxRapidity = 1.) { const double coalescenceRadius{0.5 * 1.122462 * coalMomentum}; // if coalescence from a heavy hadron, loop only between firstDauID and lastDauID @@ -131,7 +131,7 @@ bool CoalescencePythia8(Pythia8::Event& event, std::vector inputPd // fill nucleon pools std::vector protons[2], neutrons[2], lambdas[2]; for (auto iPart{loopStart}; iPart <= loopEnd; ++iPart) { - if (std::abs(event[iPart].y()) > 1.) // skip particles with y > 1 + if (std::abs(event[iPart].y()) > maxRapidity) // skip particles with y > ymax { continue; } From 70179ce8524032b3a5fc37571bff0e7a0b1f092f Mon Sep 17 00:00:00 2001 From: Fabrizio Grosa Date: Mon, 1 Jun 2026 13:36:28 +0200 Subject: [PATCH 4/4] Add test macro and change lifetime from 1 mm to 0.5 mm --- .../generator_pythia8_embed_charmnuclei.C | 23 ++-- .../ini/GeneratorHF_CDeuteron_Injected.ini | 2 +- .../tests/GeneratorHF_CDeuteron_Injected.C | 127 ++++++++++++++++++ .../generator/pythia8_cdeuteron_to_dkpi.cfg | 2 +- 4 files changed, 144 insertions(+), 10 deletions(-) create mode 100644 MC/config/PWGHF/ini/tests/GeneratorHF_CDeuteron_Injected.C diff --git a/MC/config/PWGHF/external/generator/generator_pythia8_embed_charmnuclei.C b/MC/config/PWGHF/external/generator/generator_pythia8_embed_charmnuclei.C index da4224489..4b6e9a3f7 100644 --- a/MC/config/PWGHF/external/generator/generator_pythia8_embed_charmnuclei.C +++ b/MC/config/PWGHF/external/generator/generator_pythia8_embed_charmnuclei.C @@ -36,10 +36,15 @@ class GeneratorPythia8HFEmbedCharmNuclei : public o2::eventgen::GeneratorPythia8 LOG(fatal) << "********** [GeneratorPythia8HFEmbedCharmNuclei] Only c-deuteron (pdg=2010010020) currently supported! Exit **********"; } mLifetimeCharmNucleus = lifetime; - mDecayDistr = new TF1("mDecayDistr", Form("exp(-x/%f)", mLifetimeCharmNucleus), 0., 1000.); - mDecayDistrLb = new TF1("mDecayDistrLb", "exp(-x/0.440)", 0., 1000.); - mPtDistrLb = new TF1("mPtDistrLb","[0]*x/TMath::Power((1+TMath::Power(x/[1],[3])),[2])",0.,1000.); + mDecayDistr = new TF1("mDecayDistr", "TMath::Exp(-x * 1./[0])", 0., mLifetimeCharmNucleus * 100); + mDecayDistr->SetNpx(10000); + mDecayDistr->SetParameter(0, mLifetimeCharmNucleus); + mDecayDistrLb = new TF1("mDecayDistrLb", "TMath::Exp(-x * 1./[0])", 0., 44.f); + mDecayDistrLb->SetParameter(0, 0.440f); // lifetime of Lambda_b in mm/c + mDecayDistrLb->SetNpx(10000); + mPtDistrLb = new TF1("mPtDistrLb","[0]*x/TMath::Power((1+TMath::Power(x/[1],[3])),[2])",0.,100.); mPtDistrLb->SetParameters(1000., 6.97355, 3.20721, 1.71678); + mPtDistrLb->SetNpx(10000); Print(); @@ -114,7 +119,7 @@ class GeneratorPythia8HFEmbedCharmNuclei : public o2::eventgen::GeneratorPythia8 int pdgToGen = mPdgCharmNucleus; float massToGen = mMassCharmNucleus; - float lifetimeToGen = mDecayDistr->GetRandom(); + float lifetimeToGen = 0.f; float minRapToGen = mRapidityMinCharmNuclei; float maxRapToGen = mRapidityMaxCharmNuclei; bool isFromB = gRandom->Rndm() < mFractionFromBeauty; @@ -125,6 +130,8 @@ class GeneratorPythia8HFEmbedCharmNuclei : public o2::eventgen::GeneratorPythia8 lifetimeToGen = mDecayDistrLb->GetRandom(); minRapToGen *= 2; maxRapToGen *= 2; + } else { + lifetimeToGen = mDecayDistr->GetRandom(); } auto pt = (!isFromB) ? gRandom->Uniform(0., mPtMaxCharmNuclei) : mPtDistrLb->GetRandom(); @@ -244,16 +251,16 @@ class GeneratorPythia8HFEmbedCharmNuclei : public o2::eventgen::GeneratorPythia8 auto daughter1 = part.daughter1(); auto daughter2 = part.daughter2(); if (mother1 > 0) { - part.mother1(mother1 + offset); + part.mother1(mother1 + offset - 1); } if (mother2 > 0) { - part.mother2(mother2 + offset); + part.mother2(mother2 + offset - 1); } if (daughter1 > 0) { - part.daughter1(daughter1 + offset); + part.daughter1(daughter1 + offset - 1); } if (daughter2 > 0) { - part.daughter2(daughter2 + offset); + part.daughter2(daughter2 + offset - 1); } mPythia.event.append(part); } diff --git a/MC/config/PWGHF/ini/GeneratorHF_CDeuteron_Injected.ini b/MC/config/PWGHF/ini/GeneratorHF_CDeuteron_Injected.ini index 5a97ca2a4..2a574f340 100644 --- a/MC/config/PWGHF/ini/GeneratorHF_CDeuteron_Injected.ini +++ b/MC/config/PWGHF/ini/GeneratorHF_CDeuteron_Injected.ini @@ -2,7 +2,7 @@ ### The external generator derives from GeneratorPythia8. [GeneratorExternal] fileName=${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGHF/external/generator/generator_pythia8_embed_charmnuclei.C -funcName=GenerateHFEmbedCDeuteron(1., 10, -1., 1., 25., 0, 0.4, 0.25) +funcName=GenerateHFEmbedCDeuteron(0.5, 10, -1., 1., 25., 0, 0.4, 0.25) ### includePartonEvent is needed to keep the c-deuteron in the event record, even if there are no partons in the event [GeneratorPythia8] diff --git a/MC/config/PWGHF/ini/tests/GeneratorHF_CDeuteron_Injected.C b/MC/config/PWGHF/ini/tests/GeneratorHF_CDeuteron_Injected.C new file mode 100644 index 000000000..56d22c72a --- /dev/null +++ b/MC/config/PWGHF/ini/tests/GeneratorHF_CDeuteron_Injected.C @@ -0,0 +1,127 @@ +int External() { + std::string path{"o2sim_Kine.root"}; + + int pdgCDeuteron{2010010020}; + int checkNumberOfCDeuteronPerEvent{10}; + float checkLifetimeCDeuteron{0.05f}; + std::map>> checkDecays{ + {2010010020, {{-321, 211, 1000010020}}} // c-deuteron -> K- + pi+ + deuteron + }; + float checkFracCDeuteronFromBeauty{0.25}; + + 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 nCDeuteron{}, nCDeuteronGoodDecay{}, nCDeuteronFromBeauty{}; + std::array averageLifetimeCDeuteron{0.f, 0.f}; // prompt and non-prompt + float massCDeuteron = 3.226f; + auto nEvents = tree->GetEntries(); + + for (int i = 0; i < nEvents; i++) { + tree->GetEntry(i); + + for (auto &track : *tracks) { + auto pdg = track.GetPdgCode(); + auto absPdg = std::abs(pdg); + // std::cout << "Event " << i << ": found particle with PDG " << pdg << std::endl; + + if (absPdg == pdgCDeuteron) { // found signal + nCDeuteron++; // count signal PDG + + // std::cout << "Event " << i << ": found c-deuteron with PDG " << pdg << std::endl; + + std::vector pdgsDecay{}; + std::vector pdgsDecayAntiPart{}; + if (track.getFirstDaughterTrackId() >= 0 && track.getLastDaughterTrackId() >= 0) { + for (int j{track.getFirstDaughterTrackId()}; j <= track.getLastDaughterTrackId(); ++j) { + // std::cout << "Fetching daughter track with ID " << j << std::endl; + auto pdgDau = tracks->at(j).GetPdgCode(); + // std::cout << "PDG of daughter track " << j << ": " << pdgDau << std::endl; + pdgsDecay.push_back(pdgDau); + if (pdgDau != 333 && pdgDau != 111) { // phi and pi0 are antiparticles of themselves + pdgsDecayAntiPart.push_back(-pdgDau); + } else { + pdgsDecayAntiPart.push_back(pdgDau); + } + } + } + // std::cout << "Daughters fetched" << std::endl; + + auto mother = track.getMotherTrackId(); + bool isFromBeauty{false}; + if (mother >= 0 && std::abs(tracks->at(mother).GetPdgCode()) == 5122) { // check if c-deuteron comes from Lb + nCDeuteronFromBeauty++; + isFromBeauty = true; + } + + auto dauTrack = tracks->at(track.getFirstDaughterTrackId()); + float decayLength = std::sqrt((track.GetStartVertexCoordinatesX() - dauTrack.GetStartVertexCoordinatesX()) * (track.GetStartVertexCoordinatesX() - dauTrack.GetStartVertexCoordinatesX()) + (track.GetStartVertexCoordinatesY() - dauTrack.GetStartVertexCoordinatesY()) * (track.GetStartVertexCoordinatesY() - dauTrack.GetStartVertexCoordinatesY()) + (track.GetStartVertexCoordinatesZ() - dauTrack.GetStartVertexCoordinatesZ()) * (track.GetStartVertexCoordinatesZ() - dauTrack.GetStartVertexCoordinatesZ())); + if (!isFromBeauty) { + averageLifetimeCDeuteron[0] += decayLength * massCDeuteron / track.GetP(); + } else { + averageLifetimeCDeuteron[1] += decayLength * massCDeuteron / track.GetP(); + } + + std::sort(pdgsDecay.begin(), pdgsDecay.end()); + std::sort(pdgsDecayAntiPart.begin(), pdgsDecayAntiPart.end()); + + for (auto &decay : checkDecays[std::abs(pdg)]) { + if (pdgsDecay == decay || pdgsDecayAntiPart == decay) { + nCDeuteronGoodDecay++; + break; + } + } + // std::cout << "Daughters checked " << std::endl; + } + } + } + + averageLifetimeCDeuteron[0] /= nCDeuteron - nCDeuteronFromBeauty; + averageLifetimeCDeuteron[1] /= nCDeuteronFromBeauty; + + std::cout << "--------------------------------\n"; + std::cout << "# Events: " << nEvents << "\n"; + std::cout <<"# signal c-deuteron: " << nCDeuteron << "\n"; + std::cout <<"# signal c-deuteron decaying in the correct channel: " << nCDeuteronGoodDecay << "\n"; + std::cout <<"# signal c-deuteron from beauty: " << nCDeuteronFromBeauty << "\n"; + std::cout <<"Average lifetime of c-deuteron (prompt): " << averageLifetimeCDeuteron[0] << " (cm) \n"; + std::cout <<"Average lifetime of c-deuteron (non-prompt): " << averageLifetimeCDeuteron[1] << " (cm) \n"; + + float numberOfCDeuteronPerEvent = float(nCDeuteron) / nEvents; + float fracCDeuteronGoodDecay = float(nCDeuteronGoodDecay) / nCDeuteron; + float fracCDeuteronFromBeauty = float(nCDeuteronFromBeauty) / nCDeuteron; + + if (std::abs(numberOfCDeuteronPerEvent - checkNumberOfCDeuteronPerEvent) / numberOfCDeuteronPerEvent > 0.05) { // we put some tolerance since the number of generated events is small + std::cerr << "Number of C-deuterons per event " << numberOfCDeuteronPerEvent << " different than expected " << checkNumberOfCDeuteronPerEvent << "\n"; + return 1; + } + + if (fracCDeuteronGoodDecay < 0.95) { // we put some tolerance since the number of generated events is small + std::cerr << "Fraction of signals decaying into the correct channel " << fracCDeuteronGoodDecay << " lower than expected\n"; + return 1; + } + + if (std::abs(fracCDeuteronFromBeauty - checkFracCDeuteronFromBeauty) / checkFracCDeuteronFromBeauty > 0.10) { // we put some tolerance since the number of generated events is small + std::cerr << "Fraction of signals from beauty " << fracCDeuteronFromBeauty << " different than expected " << checkFracCDeuteronFromBeauty << "\n"; + return 1; + } + + if (std::abs(averageLifetimeCDeuteron[0] - checkLifetimeCDeuteron) / checkLifetimeCDeuteron > 0.10) { // we put some tolerance since the number of generated events is small + std::cerr << "Lifetime for prompt c-deuteron " << averageLifetimeCDeuteron[0] << " different than expected " << checkLifetimeCDeuteron << "\n"; + return 1; + } + + if (std::abs(averageLifetimeCDeuteron[1] - checkLifetimeCDeuteron) / checkLifetimeCDeuteron > 0.10) { // we put some tolerance since the number of generated events is small + std::cerr << "Lifetime for non-prompt c-deuteron " << averageLifetimeCDeuteron[1] << " different than expected " << checkLifetimeCDeuteron << "\n"; + return 1; + } + + return 0; +} diff --git a/MC/config/PWGHF/pythia8/generator/pythia8_cdeuteron_to_dkpi.cfg b/MC/config/PWGHF/pythia8/generator/pythia8_cdeuteron_to_dkpi.cfg index 9fda88564..788e56fe3 100644 --- a/MC/config/PWGHF/pythia8/generator/pythia8_cdeuteron_to_dkpi.cfg +++ b/MC/config/PWGHF/pythia8/generator/pythia8_cdeuteron_to_dkpi.cfg @@ -15,7 +15,7 @@ ParticleDecays:tau0Max 10. ### add the c-deuteron ### id:all = name antiName spinType chargeType colType m0 mWidth mMin mMax tau0 -2010010020:all = CDeuteron AntiCDeuteron 1 3 0 3.226 0 3.226 3.226 1. +2010010020:all = CDeuteron AntiCDeuteron 1 3 0 3.226 0 3.226 3.226 0.5 2010010020:mayDecay = on # same as Λc+ -> p K- π+ including resonances + a neutron