Skip to content

Commit e71700c

Browse files
fgrosafgrosa
committed
Add possibility to simulate also non-prompt c-deuterons
1 parent ab26c14 commit e71700c

4 files changed

Lines changed: 169 additions & 39 deletions

File tree

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

Lines changed: 95 additions & 32 deletions
Original file line numberDiff line numberDiff line change
@@ -19,26 +19,30 @@ class GeneratorPythia8HFEmbedCharmNuclei : public o2::eventgen::GeneratorPythia8
1919
public:
2020

2121
/// constructor
22-
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)
22+
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)
2323
{
2424
nNumberOfCharmNucleiPerEvent = nCharmNucleiPerEvent;
2525
mRapidityMinCharmNuclei = yMin;
2626
mRapidityMaxCharmNuclei = yMax;
27+
mPtMaxCharmNuclei = ptMax;
2728
mTrivialCoal = trivialCoal;
2829
mCoalMomentum = coalMomentum;
30+
mFractionFromBeauty = fracFromB;
2931
mPdgCharmNucleus = pdgCode;
32+
mSign = 1;
3033
if (std::abs(mPdgCharmNucleus) == 2010010020) {
3134
mMassCharmNucleus = 3.226f;
3235
} else {
3336
LOG(fatal) << "********** [GeneratorPythia8HFEmbedCharmNuclei] Only c-deuteron (pdg=2010010020) currently supported! Exit **********";
3437
}
3538
mLifetimeCharmNucleus = lifetime;
3639
mDecayDistr = new TF1("mDecayDistr", Form("exp(-x/%f)", mLifetimeCharmNucleus), 0., 1000.);
40+
mDecayDistrLb = new TF1("mDecayDistrLb", "exp(-x/0.440)", 0., 1000.);
41+
mPtDistrLb = new TF1("mPtDistrLb","[0]*x/TMath::Power((1+TMath::Power(x/[1],[3])),[2])",0.,1000.);
42+
mPtDistrLb->SetParameters(1000., 6.97355, 3.20721, 1.71678);
3743

3844
Print();
3945

40-
/** switch off process level **/
41-
// mPythiaGun.readString("ProcessLevel:all off");
4246
auto& param = o2::eventgen::GeneratorPythia8Param::Instance();
4347
LOG(info) << "Init \'GeneratorPythia8HFEmbedCharmNuclei\' with following parameters";
4448
LOG(info) << param;
@@ -67,9 +71,11 @@ class GeneratorPythia8HFEmbedCharmNuclei : public o2::eventgen::GeneratorPythia8
6771
LOG(info) << Form("* Mass of charm nuclei to be injected (GeV/c2): %f", mMassCharmNucleus);
6872
LOG(info) << Form("* Lifetime of charm nuclei to be injected (mm): %f", mLifetimeCharmNucleus);
6973
LOG(info) << Form("* Number of charm nuclei injected per event: %d", nNumberOfCharmNucleiPerEvent);
70-
LOG(info) << Form("* Hadron rapidity: %f - %f", mRapidityMinCharmNuclei, mRapidityMaxCharmNuclei);
74+
LOG(info) << Form("* Charmed nucleus rapidity: %f - %f", mRapidityMinCharmNuclei, mRapidityMaxCharmNuclei);
75+
LOG(info) << Form("* Charmed nucleus pT max (prompt): %f", mPtMaxCharmNuclei);
7176
LOG(info) << Form("* Trivial coalescence: %d", mTrivialCoal);
7277
LOG(info) << Form("* Coalescence momentum: %f", mCoalMomentum);
78+
LOG(info) << Form("* Fraction from beauty: %f", mFractionFromBeauty);
7379
LOG(info) << "***********************************************************************";
7480
}
7581

@@ -97,44 +103,58 @@ class GeneratorPythia8HFEmbedCharmNuclei : public o2::eventgen::GeneratorPythia8
97103

98104
// we simulate c-deuteron decays
99105
for (int iCharmNuclei{0}; iCharmNuclei<nNumberOfCharmNucleiPerEvent; ++iCharmNuclei) {
100-
int sign = 1;
101-
// we alternate positive and negative
102-
if (nNumberOfCharmNucleiPerEvent % 2 == 0) {
103-
if (iCharmNuclei % 2 != 0) {
104-
sign = -1;
105-
}
106-
} else {
106+
107+
// we alternate the sign of the generated charmed nuclei, if mSign is set to 0, they are generated with 50% of probability as particle or antiparticle
108+
(mSign > 0) ? mSign = -1 : mSign = 1;
109+
if (nNumberOfCharmNucleiPerEvent % 2 != 0 && iCharmNuclei == nNumberOfCharmNucleiPerEvent - 1) {
107110
if (gRandom->Rndm() < 0.5) {
108-
sign = -1;
111+
mSign = 1;
109112
}
110113
}
111114

112-
auto pt = gRandom->Uniform(0., 50.); // placeholder, to be modified
113-
auto y = gRandom->Uniform(mRapidityMinCharmNuclei, mRapidityMaxCharmNuclei);
115+
int pdgToGen = mPdgCharmNucleus;
116+
float massToGen = mMassCharmNucleus;
117+
float lifetimeToGen = mDecayDistr->GetRandom();
118+
float minRapToGen = mRapidityMinCharmNuclei;
119+
float maxRapToGen = mRapidityMaxCharmNuclei;
120+
bool isFromB = gRandom->Rndm() < mFractionFromBeauty;
121+
// we determine if it's prompt or non-prompt
122+
if (isFromB) {
123+
pdgToGen = 5122; // we generate a Lambda_b and we let it decay into the charmed nucleus, no other beauty hadrons are considered
124+
massToGen = 5.61940f; // mass of Lambda_b (GeV/c2)
125+
lifetimeToGen = mDecayDistrLb->GetRandom();
126+
minRapToGen *= 2;
127+
maxRapToGen *= 2;
128+
}
129+
130+
auto pt = (!isFromB) ? gRandom->Uniform(0., mPtMaxCharmNuclei) : mPtDistrLb->GetRandom();
131+
auto y = gRandom->Uniform(minRapToGen, maxRapToGen);
114132
auto phi = gRandom->Uniform(0, TMath::TwoPi());
115133
auto px = pt * TMath::Cos(phi);
116134
auto py = pt * TMath::Sin(phi);
117-
auto mt = TMath::Sqrt(mMassCharmNucleus * mMassCharmNucleus + pt * pt);
135+
auto mt = TMath::Sqrt(massToGen * massToGen + pt * pt);
118136
auto pz = mt * TMath::SinH(y);
119137
auto p = TMath::Sqrt(pt * pt + pz * pz);
120-
auto e = TMath::Sqrt(mMassCharmNucleus * mMassCharmNucleus + p * p);
138+
auto e = TMath::Sqrt(massToGen * massToGen + p * p);
121139

122140
Particle particle;
123-
particle.id(sign * mPdgCharmNucleus);
141+
particle.id(mSign * pdgToGen);
124142
particle.status(83);
125-
particle.m(mMassCharmNucleus);
143+
particle.m(massToGen);
126144
particle.px(px);
127145
particle.py(py);
128146
particle.pz(pz);
129147
particle.e(e);
130148
particle.xProd(0.f);
131149
particle.yProd(0.f);
132150
particle.zProd(0.f);
133-
particle.tau(0.f); //mDecayDistr->GetRandom());
151+
particle.tau(lifetimeToGen);
152+
mPythiaGun.particleData.mayDecay(5122, true); // force decay
134153
mPythiaGun.particleData.mayDecay(mPdgCharmNucleus, true); // force decay
135154

136155
bool isCoalSuccess{false};
137-
while(!isCoalSuccess) {
156+
int nTrials{0};
157+
while(!isCoalSuccess || nTrials > 1e4) {
138158
mPythiaGun.event.reset();
139159
mPythiaGun.event.append(particle);
140160
mPythiaGun.moreDecays();
@@ -148,8 +168,9 @@ class GeneratorPythia8HFEmbedCharmNuclei : public o2::eventgen::GeneratorPythia8
148168
if (absPdg == mPdgCharmNucleus) {
149169
idxCharmNucleus = iPart;
150170
}
171+
auto mother = part.mother1();
151172
// if we find a resonance, we remove it, otherwise we prevent the coalescence of daughters from resonances and daughters from charmed nucleus directly
152-
if (std::find(pdgShortLivedResos.begin(), pdgShortLivedResos.end(), absPdg) != pdgShortLivedResos.end()) {
173+
if (std::find(pdgShortLivedResos.begin(), pdgShortLivedResos.end(), absPdg) != pdgShortLivedResos.end() && mother >= idxCharmNucleus) {
153174
// we need to change the indices of the daughter particles to point to the charmed nucleus
154175
auto dauList = part.daughterList();
155176
for (auto const& dau : dauList) {
@@ -160,22 +181,57 @@ class GeneratorPythia8HFEmbedCharmNuclei : public o2::eventgen::GeneratorPythia8
160181
}
161182
}
162183
if (isResoFound) { // we have to reset all the particles as daughters of the charm nucleus
163-
mPythiaGun.event[idxCharmNucleus].daughter1(idxCharmNucleus + 1);
164-
mPythiaGun.event[idxCharmNucleus].daughter2(mPythiaGun.event.size() - 1);
184+
std::vector<int> idxDausCharmNucleus{};
185+
for (int iPart{0}; iPart<mPythiaGun.event.size(); ++iPart) {
186+
auto mother = mPythiaGun.event[iPart].mother1();
187+
if (mother == idxCharmNucleus) {
188+
idxDausCharmNucleus.push_back(iPart);
189+
}
190+
}
191+
if (idxDausCharmNucleus.size() == idxDausCharmNucleus.back() - idxDausCharmNucleus[0] + 1) {
192+
mPythiaGun.event[idxCharmNucleus].daughter1(idxDausCharmNucleus[0]);
193+
mPythiaGun.event[idxCharmNucleus].daughter2(idxDausCharmNucleus.back());
194+
} else { // the history is broken, we need to restore it
195+
std::vector<Particle> newPartList{};
196+
std::vector<int> idxToRemove{};
197+
for (int iPart{idxDausCharmNucleus[0]}; iPart<mPythiaGun.event.size(); ++iPart) {
198+
if (std::find(idxDausCharmNucleus.begin(), idxDausCharmNucleus.end(), iPart) == idxDausCharmNucleus.end()) {
199+
newPartList.push_back(mPythiaGun.event[iPart]);
200+
idxToRemove.push_back(iPart);
201+
}
202+
}
203+
int removed{0};
204+
for (auto const& idx : idxToRemove) {
205+
mPythiaGun.event.remove(idx - removed, idx - removed, false);
206+
++removed;
207+
}
208+
std::vector<int> updatedMothers{};
209+
for (int iPart{0}; iPart<newPartList.size(); ++iPart) {
210+
mPythiaGun.event.append(newPartList[iPart]);
211+
auto mother = newPartList[iPart].mother1();
212+
if (std::find(updatedMothers.begin(), updatedMothers.end(), mother) == updatedMothers.end()) {
213+
updatedMothers.push_back(mother);
214+
int delta = mPythiaGun.event.size() - idxToRemove[iPart] - 1;
215+
mPythiaGun.event[mother].daughters(mPythiaGun.event[mother].daughter1() + delta, mPythiaGun.event[mother].daughter2() + delta);
216+
}
217+
}
218+
mPythiaGun.event[idxCharmNucleus].daughter1(idxDausCharmNucleus[0]);
219+
mPythiaGun.event[idxCharmNucleus].daughter2(idxDausCharmNucleus[0] + idxDausCharmNucleus.size() - 1);
220+
}
165221
}
166222

167223
int iDau{-1};
168224
for (int iPart{0}; iPart<mPythiaGun.event.size(); ++iPart) {
169-
auto part = mPythiaGun.event[iPart];
170-
auto absPdg = std::abs(part.id());
225+
auto absPdg = std::abs(mPythiaGun.event[iPart].id());
226+
auto mother = mPythiaGun.event[iPart].mother1();
171227

172-
if (absPdg == 2212 || absPdg == 2112) { // coalescence of protons and deuterons
228+
if ((absPdg == 2212 || absPdg == 2112) && mother == idxCharmNucleus) { // coalescence of protons and deuterons
173229
dausToCoal[++iDau] = iPart;
174230
}
175231
}
176232

177233
// we try the coalescence here, if successful we copy particles in the pythia event and we move to the next charm nucleus
178-
isCoalSuccess = CoalescencePythia8(mPythiaGun.event, std::vector<unsigned int>{1000010020}, mTrivialCoal, mCoalMomentum, dausToCoal[0], dausToCoal[1]);
234+
isCoalSuccess = CoalescencePythia8(mPythiaGun.event, std::vector<unsigned int>{1000010020}, mTrivialCoal, mCoalMomentum, dausToCoal[0], dausToCoal[1], 10.);
179235
if (isCoalSuccess) {
180236
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
181237
for (int iPart{0}; iPart<mPythiaGun.event.size(); ++iPart) {
@@ -202,32 +258,39 @@ class GeneratorPythia8HFEmbedCharmNuclei : public o2::eventgen::GeneratorPythia8
202258
mPythia.event.append(part);
203259
}
204260
}
261+
nTrials++;
205262
}
206263
}
207264

208265
return true;
209266
}
210267

211-
private:
268+
269+
private:
212270
// Properties of selection
213271
float mMassCharmNucleus; /// mass of the charmed nucleus
214272
int mPdgCharmNucleus; /// pdg code of the charmed nucleus
215273
float mLifetimeCharmNucleus; /// lifetime of the charmed nucleus
216274
int nNumberOfCharmNucleiPerEvent; /// number of charmed nuclei injected per event
217-
float mRapidityMinCharmNuclei; /// rapidity min
218-
float mRapidityMaxCharmNuclei; /// rapidity max
275+
float mRapidityMinCharmNuclei; /// rapidity min of the generated charmed nuclei
276+
float mRapidityMaxCharmNuclei; /// rapidity max of the generated charmed nuclei
277+
float mPtMaxCharmNuclei; /// pT max of the generated charmed nuclei
219278
unsigned int mUsedSeed; /// seed
220279
bool mTrivialCoal; /// if true, the coalescence is done without checking the distance in the phase space of the nucleons
221280
float mCoalMomentum; /// coalescence momentum
222281
Pythia8::Pythia mPythiaGun; /// Gun generator with decay support
223282
TF1* mDecayDistr; /// Lifetime distribution
283+
TF1* mDecayDistrLb; /// Lifetime distribution for Lb
284+
TF1* mPtDistrLb; /// pt distribution for Lb (power-law fit to FONLL)
285+
float mFractionFromBeauty; /// fraction of charmed nuclei coming from beauty hadrons
286+
int mSign; /// sign of the charmed nuclei to be generated, if 0 they are generated with 50% of probability as particle or antiparticle
224287
};
225288

226289

227290
///___________________________________________________________
228-
FairGenerator *GenerateHFEmbedCDeuteron(float lifetime = 1.f, int nCharmNucleiPerEvent = 10, float yMin = -1.f, float yMax = 1.f, bool trivialCoal = false, float coalMomentum = 0.2)
291+
FairGenerator *GenerateHFEmbedCDeuteron(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.2f, float fracFromB = 0.25f)
229292
{
230-
auto myGen = new GeneratorPythia8HFEmbedCharmNuclei(2010010020, lifetime, nCharmNucleiPerEvent, yMin, yMax, trivialCoal, coalMomentum);
293+
auto myGen = new GeneratorPythia8HFEmbedCharmNuclei(2010010020, lifetime, nCharmNucleiPerEvent, yMin, yMax, ptMax, trivialCoal, coalMomentum, fracFromB);
231294
auto seed = (gRandom->TRandom::GetSeed() % 900000000);
232295
myGen->readString("Random:setSeed on");
233296
myGen->readString("Random:seed " + std::to_string(seed));

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

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -2,7 +2,7 @@
22
### The external generator derives from GeneratorPythia8.
33
[GeneratorExternal]
44
fileName=${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGHF/external/generator/generator_pythia8_embed_charmnuclei.C
5-
funcName=GenerateHFEmbedCDeuteron(1., 10, -1., 1., 0, 0.4)
5+
funcName=GenerateHFEmbedCDeuteron(1., 10, -1., 1., 25., 0, 0.4, 0.25)
66

77
### includePartonEvent is needed to keep the c-deuteron in the event record, even if there are no partons in the event
88
[GeneratorPythia8]

0 commit comments

Comments
 (0)