Skip to content

Commit e882201

Browse files
fgrosafgrosa
authored
Add generator and config for charmed nuclei (c-deuteron) (#2368)
* Add configuration for c-deuteron simulation * Fix coalescence and definition of event containing decayed charmed nuclei * Add possibility to simulate also non-prompt c-deuterons * Add test macro and change lifetime from 1 mm to 0.5 mm
1 parent 216e56d commit e882201

5 files changed

Lines changed: 546 additions & 2 deletions

File tree

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

Lines changed: 305 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,305 @@
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+
21+
/// constructor
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)
23+
{
24+
nNumberOfCharmNucleiPerEvent = nCharmNucleiPerEvent;
25+
mRapidityMinCharmNuclei = yMin;
26+
mRapidityMaxCharmNuclei = yMax;
27+
mPtMaxCharmNuclei = ptMax;
28+
mTrivialCoal = trivialCoal;
29+
mCoalMomentum = coalMomentum;
30+
mFractionFromBeauty = fracFromB;
31+
mPdgCharmNucleus = pdgCode;
32+
mSign = 1;
33+
if (std::abs(mPdgCharmNucleus) == 2010010020) {
34+
mMassCharmNucleus = 3.226f;
35+
} else {
36+
LOG(fatal) << "********** [GeneratorPythia8HFEmbedCharmNuclei] Only c-deuteron (pdg=2010010020) currently supported! Exit **********";
37+
}
38+
mLifetimeCharmNucleus = lifetime;
39+
mDecayDistr = new TF1("mDecayDistr", "TMath::Exp(-x * 1./[0])", 0., mLifetimeCharmNucleus * 100);
40+
mDecayDistr->SetNpx(10000);
41+
mDecayDistr->SetParameter(0, mLifetimeCharmNucleus);
42+
mDecayDistrLb = new TF1("mDecayDistrLb", "TMath::Exp(-x * 1./[0])", 0., 44.f);
43+
mDecayDistrLb->SetParameter(0, 0.440f); // lifetime of Lambda_b in mm/c
44+
mDecayDistrLb->SetNpx(10000);
45+
mPtDistrLb = new TF1("mPtDistrLb","[0]*x/TMath::Power((1+TMath::Power(x/[1],[3])),[2])",0.,100.);
46+
mPtDistrLb->SetParameters(1000., 6.97355, 3.20721, 1.71678);
47+
mPtDistrLb->SetNpx(10000);
48+
49+
Print();
50+
51+
auto& param = o2::eventgen::GeneratorPythia8Param::Instance();
52+
LOG(info) << "Init \'GeneratorPythia8HFEmbedCharmNuclei\' with following parameters";
53+
LOG(info) << param;
54+
if (param.config.empty()) {
55+
LOG(fatal) << "Failed to init \'GeneratorPythia8\': problems with configuration file ";
56+
}
57+
std::string cfg = gSystem->ExpandPathName(param.config.c_str());
58+
LOG(info) << "GeneratorPythia8HFEmbedCharmNuclei Reading configuration from file: " << cfg;
59+
if (!mPythiaGun.readFile(cfg, true)) {
60+
LOG(fatal) << "Failed to init \'GeneratorPythia8\': problems with configuration file " << cfg;
61+
}
62+
63+
if (!mPythiaGun.init()) {
64+
LOG(fatal) << "Failed to init \'GeneratorPythia8\': init returned with error";
65+
}
66+
}
67+
68+
/// Destructor
69+
~GeneratorPythia8HFEmbedCharmNuclei() = default;
70+
71+
/// Print the input
72+
void Print()
73+
{
74+
LOG(info) << "********** GeneratorPythia8HFEmbedCharmNuclei configuration dump **********";
75+
LOG(info) << Form("* PDG code of charm nuclei to be injected: %d", mPdgCharmNucleus);
76+
LOG(info) << Form("* Mass of charm nuclei to be injected (GeV/c2): %f", mMassCharmNucleus);
77+
LOG(info) << Form("* Lifetime of charm nuclei to be injected (mm): %f", mLifetimeCharmNucleus);
78+
LOG(info) << Form("* Number of charm nuclei injected per event: %d", nNumberOfCharmNucleiPerEvent);
79+
LOG(info) << Form("* Charmed nucleus rapidity: %f - %f", mRapidityMinCharmNuclei, mRapidityMaxCharmNuclei);
80+
LOG(info) << Form("* Charmed nucleus pT max (prompt): %f", mPtMaxCharmNuclei);
81+
LOG(info) << Form("* Trivial coalescence: %d", mTrivialCoal);
82+
LOG(info) << Form("* Coalescence momentum: %f", mCoalMomentum);
83+
LOG(info) << Form("* Fraction from beauty: %f", mFractionFromBeauty);
84+
LOG(info) << "***********************************************************************";
85+
}
86+
87+
void setHadronRapidity(float yMin, float yMax)
88+
{
89+
mRapidityMinCharmNuclei = yMin;
90+
mRapidityMaxCharmNuclei = yMax;
91+
};
92+
93+
void setUsedSeed(unsigned int seed)
94+
{
95+
mUsedSeed = seed;
96+
};
97+
98+
unsigned int getUsedSeed() const
99+
{
100+
return mUsedSeed;
101+
};
102+
103+
//__________________________________________________________________
104+
bool generateEvent() override
105+
{
106+
// we start from an empty event
107+
mPythia.event.reset();
108+
109+
// we simulate c-deuteron decays
110+
for (int iCharmNuclei{0}; iCharmNuclei<nNumberOfCharmNucleiPerEvent; ++iCharmNuclei) {
111+
112+
// 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
113+
(mSign > 0) ? mSign = -1 : mSign = 1;
114+
if (nNumberOfCharmNucleiPerEvent % 2 != 0 && iCharmNuclei == nNumberOfCharmNucleiPerEvent - 1) {
115+
if (gRandom->Rndm() < 0.5) {
116+
mSign = 1;
117+
}
118+
}
119+
120+
int pdgToGen = mPdgCharmNucleus;
121+
float massToGen = mMassCharmNucleus;
122+
float lifetimeToGen = 0.f;
123+
float minRapToGen = mRapidityMinCharmNuclei;
124+
float maxRapToGen = mRapidityMaxCharmNuclei;
125+
bool isFromB = gRandom->Rndm() < mFractionFromBeauty;
126+
// we determine if it's prompt or non-prompt
127+
if (isFromB) {
128+
pdgToGen = 5122; // we generate a Lambda_b and we let it decay into the charmed nucleus, no other beauty hadrons are considered
129+
massToGen = 5.61940f; // mass of Lambda_b (GeV/c2)
130+
lifetimeToGen = mDecayDistrLb->GetRandom();
131+
minRapToGen *= 2;
132+
maxRapToGen *= 2;
133+
} else {
134+
lifetimeToGen = mDecayDistr->GetRandom();
135+
}
136+
137+
auto pt = (!isFromB) ? gRandom->Uniform(0., mPtMaxCharmNuclei) : mPtDistrLb->GetRandom();
138+
auto y = gRandom->Uniform(minRapToGen, maxRapToGen);
139+
auto phi = gRandom->Uniform(0, TMath::TwoPi());
140+
auto px = pt * TMath::Cos(phi);
141+
auto py = pt * TMath::Sin(phi);
142+
auto mt = TMath::Sqrt(massToGen * massToGen + pt * pt);
143+
auto pz = mt * TMath::SinH(y);
144+
auto p = TMath::Sqrt(pt * pt + pz * pz);
145+
auto e = TMath::Sqrt(massToGen * massToGen + p * p);
146+
147+
Particle particle;
148+
particle.id(mSign * pdgToGen);
149+
particle.status(83);
150+
particle.m(massToGen);
151+
particle.px(px);
152+
particle.py(py);
153+
particle.pz(pz);
154+
particle.e(e);
155+
particle.xProd(0.f);
156+
particle.yProd(0.f);
157+
particle.zProd(0.f);
158+
particle.tau(lifetimeToGen);
159+
mPythiaGun.particleData.mayDecay(5122, true); // force decay
160+
mPythiaGun.particleData.mayDecay(mPdgCharmNucleus, true); // force decay
161+
162+
bool isCoalSuccess{false};
163+
int nTrials{0};
164+
while(!isCoalSuccess || nTrials > 1e4) {
165+
mPythiaGun.event.reset();
166+
mPythiaGun.event.append(particle);
167+
mPythiaGun.moreDecays();
168+
std::array<int, 2> dausToCoal = {-1, -1};
169+
std::vector<int> pdgShortLivedResos = {313, 2224, 102134};
170+
bool isResoFound{false};
171+
int idxCharmNucleus{-1};
172+
for (int iPart{0}; iPart<mPythiaGun.event.size(); ++iPart) {
173+
auto part = mPythiaGun.event[iPart];
174+
auto absPdg = std::abs(part.id());
175+
if (absPdg == mPdgCharmNucleus) {
176+
idxCharmNucleus = iPart;
177+
}
178+
auto mother = part.mother1();
179+
// if we find a resonance, we remove it, otherwise we prevent the coalescence of daughters from resonances and daughters from charmed nucleus directly
180+
if (std::find(pdgShortLivedResos.begin(), pdgShortLivedResos.end(), absPdg) != pdgShortLivedResos.end() && mother >= idxCharmNucleus) {
181+
// we need to change the indices of the daughter particles to point to the charmed nucleus
182+
auto dauList = part.daughterList();
183+
for (auto const& dau : dauList) {
184+
mPythiaGun.event[dau].mother1(idxCharmNucleus);
185+
}
186+
mPythiaGun.event.remove(iPart, iPart, true);
187+
isResoFound=true;
188+
}
189+
}
190+
if (isResoFound) { // we have to reset all the particles as daughters of the charm nucleus
191+
std::vector<int> idxDausCharmNucleus{};
192+
for (int iPart{0}; iPart<mPythiaGun.event.size(); ++iPart) {
193+
auto mother = mPythiaGun.event[iPart].mother1();
194+
if (mother == idxCharmNucleus) {
195+
idxDausCharmNucleus.push_back(iPart);
196+
}
197+
}
198+
if (idxDausCharmNucleus.size() == idxDausCharmNucleus.back() - idxDausCharmNucleus[0] + 1) {
199+
mPythiaGun.event[idxCharmNucleus].daughter1(idxDausCharmNucleus[0]);
200+
mPythiaGun.event[idxCharmNucleus].daughter2(idxDausCharmNucleus.back());
201+
} else { // the history is broken, we need to restore it
202+
std::vector<Particle> newPartList{};
203+
std::vector<int> idxToRemove{};
204+
for (int iPart{idxDausCharmNucleus[0]}; iPart<mPythiaGun.event.size(); ++iPart) {
205+
if (std::find(idxDausCharmNucleus.begin(), idxDausCharmNucleus.end(), iPart) == idxDausCharmNucleus.end()) {
206+
newPartList.push_back(mPythiaGun.event[iPart]);
207+
idxToRemove.push_back(iPart);
208+
}
209+
}
210+
int removed{0};
211+
for (auto const& idx : idxToRemove) {
212+
mPythiaGun.event.remove(idx - removed, idx - removed, false);
213+
++removed;
214+
}
215+
std::vector<int> updatedMothers{};
216+
for (int iPart{0}; iPart<newPartList.size(); ++iPart) {
217+
mPythiaGun.event.append(newPartList[iPart]);
218+
auto mother = newPartList[iPart].mother1();
219+
if (std::find(updatedMothers.begin(), updatedMothers.end(), mother) == updatedMothers.end()) {
220+
updatedMothers.push_back(mother);
221+
int delta = mPythiaGun.event.size() - idxToRemove[iPart] - 1;
222+
mPythiaGun.event[mother].daughters(mPythiaGun.event[mother].daughter1() + delta, mPythiaGun.event[mother].daughter2() + delta);
223+
}
224+
}
225+
mPythiaGun.event[idxCharmNucleus].daughter1(idxDausCharmNucleus[0]);
226+
mPythiaGun.event[idxCharmNucleus].daughter2(idxDausCharmNucleus[0] + idxDausCharmNucleus.size() - 1);
227+
}
228+
}
229+
230+
int iDau{-1};
231+
for (int iPart{0}; iPart<mPythiaGun.event.size(); ++iPart) {
232+
auto absPdg = std::abs(mPythiaGun.event[iPart].id());
233+
auto mother = mPythiaGun.event[iPart].mother1();
234+
235+
if ((absPdg == 2212 || absPdg == 2112) && mother == idxCharmNucleus) { // coalescence of protons and deuterons
236+
dausToCoal[++iDau] = iPart;
237+
}
238+
}
239+
240+
// we try the coalescence here, if successful we copy particles in the pythia event and we move to the next charm nucleus
241+
isCoalSuccess = CoalescencePythia8(mPythiaGun.event, std::vector<unsigned int>{1000010020}, mTrivialCoal, mCoalMomentum, dausToCoal[0], dausToCoal[1], 10.);
242+
if (isCoalSuccess) {
243+
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
244+
for (int iPart{0}; iPart<mPythiaGun.event.size(); ++iPart) {
245+
auto part = mPythiaGun.event[iPart];
246+
if (part.id() == 90) {
247+
continue;
248+
}
249+
auto mother1 = part.mother1();
250+
auto mother2 = part.mother2();
251+
auto daughter1 = part.daughter1();
252+
auto daughter2 = part.daughter2();
253+
if (mother1 > 0) {
254+
part.mother1(mother1 + offset - 1);
255+
}
256+
if (mother2 > 0) {
257+
part.mother2(mother2 + offset - 1);
258+
}
259+
if (daughter1 > 0) {
260+
part.daughter1(daughter1 + offset - 1);
261+
}
262+
if (daughter2 > 0) {
263+
part.daughter2(daughter2 + offset - 1);
264+
}
265+
mPythia.event.append(part);
266+
}
267+
}
268+
nTrials++;
269+
}
270+
}
271+
272+
return true;
273+
}
274+
275+
276+
private:
277+
// Properties of selection
278+
float mMassCharmNucleus; /// mass of the charmed nucleus
279+
int mPdgCharmNucleus; /// pdg code of the charmed nucleus
280+
float mLifetimeCharmNucleus; /// lifetime of the charmed nucleus
281+
int nNumberOfCharmNucleiPerEvent; /// number of charmed nuclei injected per event
282+
float mRapidityMinCharmNuclei; /// rapidity min of the generated charmed nuclei
283+
float mRapidityMaxCharmNuclei; /// rapidity max of the generated charmed nuclei
284+
float mPtMaxCharmNuclei; /// pT max of the generated charmed nuclei
285+
unsigned int mUsedSeed; /// seed
286+
bool mTrivialCoal; /// if true, the coalescence is done without checking the distance in the phase space of the nucleons
287+
float mCoalMomentum; /// coalescence momentum
288+
Pythia8::Pythia mPythiaGun; /// Gun generator with decay support
289+
TF1* mDecayDistr; /// Lifetime distribution
290+
TF1* mDecayDistrLb; /// Lifetime distribution for Lb
291+
TF1* mPtDistrLb; /// pt distribution for Lb (power-law fit to FONLL)
292+
float mFractionFromBeauty; /// fraction of charmed nuclei coming from beauty hadrons
293+
int mSign; /// sign of the charmed nuclei to be generated, if 0 they are generated with 50% of probability as particle or antiparticle
294+
};
295+
296+
297+
///___________________________________________________________
298+
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)
299+
{
300+
auto myGen = new GeneratorPythia8HFEmbedCharmNuclei(2010010020, lifetime, nCharmNucleiPerEvent, yMin, yMax, ptMax, trivialCoal, coalMomentum, fracFromB);
301+
auto seed = (gRandom->TRandom::GetSeed() % 900000000);
302+
myGen->readString("Random:setSeed on");
303+
myGen->readString("Random:seed " + std::to_string(seed));
304+
return myGen;
305+
}

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

Lines changed: 10 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,10 @@
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(0.5, 10, -1., 1., 25., 0, 0.4, 0.25)
6+
7+
### includePartonEvent is needed to keep the c-deuteron in the event record, even if there are no partons in the event
8+
[GeneratorPythia8]
9+
config=${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGHF/pythia8/generator/pythia8_cdeuteron_to_dkpi.cfg
10+
includePartonEvent=true

0 commit comments

Comments
 (0)