Skip to content

Commit e9c834b

Browse files
fgrosafgrosa
committed
Add configuration for c-deuteron simulation
1 parent b7d1025 commit e9c834b

3 files changed

Lines changed: 242 additions & 0 deletions

File tree

‎MC/config/PWGHF/external/generator/generator_pythia8_embed_charmnuclei.C‎

Lines changed: 208 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,208 @@
1+
#include "FairGenerator.h"
2+
#include "Generators/GeneratorPythia8.h"
3+
#include "Pythia8/Pythia.h"
4+
#include "TRandom.h"
5+
#include "TF1.h"
6+
#include "TMath.h"
7+
#include <fairlogger/Logger.h>
8+
#include <algorithm>
9+
#include <string>
10+
#include <vector>
11+
12+
R__ADD_INCLUDE_PATH($O2DPG_MC_CONFIG_ROOT)
13+
#include "MC/config/common/external/generator/CoalescencePythia8.h"
14+
15+
using namespace Pythia8;
16+
17+
class GeneratorPythia8HFEmbedCharmNuclei : public o2::eventgen::GeneratorPythia8
18+
{
19+
public:
20+
/// default constructor
21+
GeneratorPythia8HFEmbedCharmNuclei() = default;
22+
23+
/// constructor
24+
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)
25+
{
26+
nNumberOfCharmNucleiPerEvent = nCharmNucleiPerEvent;
27+
mRapidityMinCharmNuclei = yMin;
28+
mRapidityMaxCharmNuclei = yMax;
29+
mTrivialCoal = trivialCoal;
30+
mCoalMomentum = coalMomentum;
31+
mPdgCharmNucleus = pdgCode;
32+
if (std::abs(mPdgCharmNucleus) == 2010010020) {
33+
mMassCharmNucleus = 3.226f;
34+
} else {
35+
LOG(fatal) << "********** [GeneratorPythia8HFEmbedCharmNuclei] Only c-deuteron (pdg=2010010020) currently supported! Exit **********";
36+
}
37+
mLifetimeCharmNucleus = lifetime;
38+
mDecayDistr = TF1("mDecayDistr", Form("exp(-x/%f)", mLifetimeCharmNucleus), 0., 1000.)
39+
40+
Print();
41+
42+
/** switch off process level **/
43+
mPythiaGun.readString("ProcessLevel:all off");
44+
auto& param = o2::eventgen::GeneratorPythia8::Instance();
45+
LOG(info) << "Init \'GeneratorPythia8HFEmbedCharmNuclei\' with following parameters";
46+
LOG(info) << param;
47+
for (int iCfg{0}; iCfg < 8; ++iCfg) {
48+
if (param.config[iCfg].empty()) {
49+
continue;
50+
}
51+
std::string config = gSystem->ExpandPathName(param.config[iCfg].c_str());
52+
LOG(info) << "GeneratorPythia8HFEmbedCharmNuclei Reading configuration from file: " << config;
53+
if (!mPythiaGun.readFile(config, true)) {
54+
LOG(fatal) << "Failed to init \'GeneratorPythia8\': problems with configuration file "
55+
<< config;
56+
return;
57+
}
58+
}
59+
if (!mPythiaGun.init()) {
60+
LOG(fatal) << "Failed to init \'GeneratorPythia8\': init returned with error";
61+
return;
62+
}
63+
}
64+
65+
/// Destructor
66+
~GeneratorPythia8HFEmbedCharmNuclei() = default;
67+
68+
/// Init
69+
bool Init() override
70+
{
71+
return o2::eventgen::GeneratorPythia8::Init();
72+
}
73+
74+
/// Print the input
75+
void Print()
76+
{
77+
LOG(info) << "********** GeneratorPythia8HFEmbedCharmNuclei configuration dump **********";
78+
LOG(info) << Form("* PDG code of charm nuclei to be injected: %d", mPdgCharmNucleus);
79+
LOG(info) << Form("* Mass of charm nuclei to be injected (GeV/c2): %f", mMassCharmNucleus);
80+
LOG(info) << Form("* Lifetime of charm nuclei to be injected (mm): %f", mLifetimeCharmNucleus);
81+
LOG(info) << Form("* Number of charm nuclei injected per event: %d", nNumberOfCharmNucleiPerEvent);
82+
LOG(info) << Form("* Hadron rapidity: %f - %f", mHadRapidityMin, mHadRapidityMax);
83+
LOG(info) << Form("* Trivial coalescence: %d", mTrivialCoal);
84+
LOG(info) << Form("* Coalescence momentum: %f", mCoalMomentum);
85+
LOG(info) << "***********************************************************************";
86+
}
87+
88+
void setHadronRapidity(float yMin, float yMax)
89+
{
90+
mRapidityMinCharmNuclei = yMin;
91+
mRapidityMaxCharmNuclei = yMax;
92+
};
93+
94+
void setUsedSeed(unsigned int seed)
95+
{
96+
mUsedSeed = seed;
97+
};
98+
99+
unsigned int getUsedSeed() const
100+
{
101+
return mUsedSeed;
102+
};
103+
104+
protected:
105+
//__________________________________________________________________
106+
bool generateEvent() override
107+
{
108+
// we start from an empty event
109+
mPythia.event.reset();
110+
111+
// we simulate c-deuteron decays
112+
for (int iCharmNuclei{0}; iCharmNuclei<nNumberOfCharmNucleiPerEvent; ++iCharmNuclei) {
113+
int sign = 1;
114+
// we alternate positive and negative
115+
if (nNumberOfCharmNucleiPerEvent % 2 == 0) {
116+
if (iCharmNuclei % 2 != 0) {
117+
sign = -1;
118+
}
119+
} else {
120+
if (gRandom->Rndm() < 0.5) {
121+
sign = -1;
122+
}
123+
}
124+
125+
auto pt = gRandom->Uniform(0., 50.); // placeholder, to be modified
126+
auto y = gRandom->Uniform(mHadRapidityMin, mHadRapidityMax);
127+
auto phi = gRandom->Uniform(0, TMath::TwoPi());
128+
auto px = pt * TMath::Cos(phi);
129+
auto py = pt * TMath::Sin(phi);
130+
auto mt = TMath::Sqrt(mMassCharmNucleus * mMassCharmNucleus + pt * pt);
131+
auto pz = mt * TMath::SinH(y);
132+
auto p = TMath::Sqrt(pt * pt + pz * pz);
133+
auto e = TMath::Sqrt(mass * mass + p * p);
134+
135+
Particle particle;
136+
particle.id(sign * mPdgCharmNucleus);
137+
particle.status(81);
138+
particle.m(mMassCharmNucleus);
139+
particle.px(px);
140+
particle.py(py);
141+
particle.pz(pz);
142+
particle.e(e);
143+
particle.xProd(0.f);
144+
particle.yProd(0.f);
145+
particle.zProd(0.f);
146+
particle.tau(mDecayDistr->GetRandom());
147+
mPythiaGun.particleData.mayDecay(mPdgCharmNucleus, true); // force decay
148+
149+
bool isCoalSuccess{false};
150+
while(!isCoalSuccess) {
151+
mPythiaGun.event.reset();
152+
mPythiaGun.event.append(particle);
153+
mPythiaGun.moreDecays();
154+
std::array<int, 2> dausToCoal = {-1, -1};
155+
int iDau{-1};
156+
for (int iPart{0}; iPart<mPythiaGun.event.size(); ++iPart) {
157+
auto part = mPythiaGun.event[iPart];
158+
if (std::abs(iPart) == 2212 || std::abs(iPart) == 2112) { // coalescence of protons and deuterons
159+
dausToCoal[++iDau] = iPart;
160+
}
161+
if (iDau == 1) { // we found the proton and the neutron to coalesce
162+
break;
163+
}
164+
}
165+
// we try the coalescence here, if successful we copy particles in the pythia event and we move to the next charm nucleus
166+
isCoalSuccess = CoalescencePythia8(mPythiaGun.event, std::vector<int>{1000010020}, mTrivialCoal, mCoalMomentum, dausToCoal[0], dausToCoal[1]);
167+
if (isCoalSuccess) {
168+
for (int iPart{0}; iPart<mPythiaGun.event.size(); ++iPart) {
169+
mPythia.append(mPythiaGun.event[iPart]);
170+
}
171+
}
172+
mPythiaGun.next();
173+
}
174+
}
175+
176+
mPythia.next();
177+
178+
return true;
179+
}
180+
181+
private:
182+
// Properties of selection
183+
float mMassCharmNucleus;
184+
int mPdgCharmNucleus;
185+
float mLifetimeCharmNucleus;
186+
int nNumberOfCharmNucleiPerEvent;
187+
float mRapidityMinCharmNuclei;
188+
float mRapidityMaxCharmNuclei;
189+
unsigned int mUsedSeed;
190+
191+
bool mTrivialCoal = false; /// if true, the coalescence is done without checking the distance in the phase space of the nucleons
192+
float mCoalMomentum; /// coalescence momentum
193+
194+
Pythia8::Pythia mPythiaGun; // Gun generator with decay support
195+
196+
TF1* mDecayDistr;
197+
};
198+
199+
200+
///___________________________________________________________
201+
FairGenerator *GenerateHFEmbedCDeuteron(float lifetime = 1.f, int nCharmNucleiPerEvent = 10, float yMin = -1.f, float yMax = 1.f, bool trivialCoal = false, float coalMomentum = 0.2)
202+
{
203+
auto myGen = new GeneratorPythia8HFEmbedCharmNuclei(2010010020, lifetime, nCharmNucleiPerEvent, yMin, yMax, trivialCoal, coalMomentum);
204+
auto seed = (gRandom->TRandom::GetSeed() % 900000000);
205+
myGen->readString("Random:setSeed on");
206+
myGen->readString("Random:seed " + std::to_string(seed));
207+
return myGen;
208+
}

‎MC/config/PWGHF/ini/GeneratorHF_CDeuteron_Injected.ini‎

Lines changed: 9 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,9 @@
1+
#NEV_TEST> 100
2+
### The external generator derives from GeneratorPythia8.
3+
[GeneratorExternal]
4+
fileName=${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGHF/external/generator/generator_pythia8_embed_charmnuclei.C
5+
funcName=GenerateHFEmbedCDeuteron(1.,10,-1.,1.,0,0.2)
6+
7+
[GeneratorPythia8]
8+
config=${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGHF/pythia8/generator/pythia8_cdeuteron_to_dkpi.cfg
9+
includePartonEvent=false

‎MC/config/PWGHF/pythia8/generator/pythia8_cdeuteron_to_dkpi.cfg‎

Lines changed: 25 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,25 @@
1+
### author: Fabrizio Grosa (fabrizio.grosa@cern.ch)
2+
### since: May 2026
3+
4+
### beams (not relevant)
5+
Beams:idA 2212 # proton
6+
Beams:idB 2212 # proton
7+
Beams:eCM 13600. # GeV
8+
9+
### processes
10+
ProcessLevel:all off # turning off all the processes, we do not have an underlying event
11+
12+
### decays
13+
ParticleDecays:limitTau0 on
14+
ParticleDecays:tau0Max 10.
15+
16+
### add the c-deuteron
17+
### id:all = name antiName spinType chargeType colType m0 mWidth mMin mMax tau0
18+
2010010020:all = CDeuteron+ AntiCDeuteron- 1 3 0 3.226 0 3.226 3.226 1.
19+
2010010020:mayDecay = on
20+
21+
# same as Λc+ -> p K- π+ including resonances + a neutron
22+
2010010020:oneChannel = 1 0.14400 0 2112 2212 -321 211 ### Λc+ -> p K- π+ (non-resonant) 3.5%
23+
2010010020:addChannel = 1 0.08100 100 2112 2212 -313 ### Λc+ -> p K*0(892) 1.96%
24+
2010010020:addChannel = 1 0.04500 100 2112 2224 -321 ### Λc+ -> Delta++ K- 1.08%
25+
2010010020:addChannel = 1 0.09000 100 2112 102134 211 ### Λc+ -> Lambda(1520) K- 2.20e-3

0 commit comments

Comments
 (0)