Skip to content

Commit f53e8a4

Browse files
fmazzascFrancesco Mazzaschi
fmazzasc
and
Francesco Mazzaschi
authored
Add common coalescence utility + HF -> nuclei generator (#1978)
* Add common coalescence utility + HF -> nuclei generator * Fix typo * Move LOG(info) --> LOG(debug) * Fix tests --------- Co-authored-by: Francesco Mazzaschi <fmazzasc@alipap1.cern.ch>
1 parent 1e4bb87 commit f53e8a4

File tree

4 files changed

+370
-90
lines changed

4 files changed

+370
-90
lines changed

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

Lines changed: 174 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,174 @@
1+
#include "FairGenerator.h"
2+
#include "Generators/GeneratorPythia8.h"
3+
#include "Pythia8/Pythia.h"
4+
#include "TRandom.h"
5+
#include <fairlogger/Logger.h>
6+
#include <string>
7+
#include <vector>
8+
9+
R__ADD_INCLUDE_PATH($O2DPG_MC_CONFIG_ROOT)
10+
#include "MC/config/common/external/generator/CoalescencePythia8.h"
11+
12+
using namespace Pythia8;
13+
14+
class GeneratorPythia8HFHadToNuclei : public o2::eventgen::GeneratorPythia8
15+
{
16+
public:
17+
/// default constructor
18+
GeneratorPythia8HFHadToNuclei() = default;
19+
20+
/// constructor
21+
GeneratorPythia8HFHadToNuclei(int inputTriggerRatio = 5, std::vector<int> hfHadronPdgList = {}, std::vector<unsigned int> nucleiPdgList = {}, bool trivialCoal = false, float coalMomentum = 0.4)
22+
{
23+
24+
mGeneratedEvents = 0;
25+
mInverseTriggerRatio = inputTriggerRatio;
26+
mHadRapidityMin = -1.5;
27+
mHadRapidityMax = 1.5;
28+
mHadronPdg = 0;
29+
mHFHadronPdgList = hfHadronPdgList;
30+
mNucleiPdgList = nucleiPdgList;
31+
mTrivialCoal = trivialCoal;
32+
mCoalMomentum = coalMomentum;
33+
Print();
34+
}
35+
36+
/// Destructor
37+
~GeneratorPythia8HFHadToNuclei() = default;
38+
39+
/// Print the input
40+
void Print()
41+
{
42+
LOG(info) << "********** GeneratorPythia8HFHadToNuclei configuration dump **********";
43+
LOG(info) << Form("* Trigger ratio: %d", mInverseTriggerRatio);
44+
LOG(info) << Form("* Hadron rapidity: %f - %f", mHadRapidityMin, mHadRapidityMax);
45+
LOG(info) << Form("* Hadron pdg list: ");
46+
for (auto pdg : mHFHadronPdgList) {
47+
LOG(info) << Form("* %d ", pdg);
48+
}
49+
LOG(info) << Form("* Trivial coalescence: %d", mTrivialCoal);
50+
LOG(info) << Form("* Coalescence momentum: %f", mCoalMomentum);
51+
LOG(info) << Form("* Nuclei pdg list: ");
52+
for (auto pdg : mNucleiPdgList) {
53+
LOG(info) << Form("* %d ", pdg);
54+
}
55+
LOG(info) << "***********************************************************************";
56+
}
57+
58+
bool Init() override
59+
{
60+
addSubGenerator(0, "Minimum bias");
61+
addSubGenerator(1, "HF + Coalescence");
62+
return o2::eventgen::GeneratorPythia8::Init();
63+
}
64+
65+
void setHadronRapidity(float yMin, float yMax)
66+
{
67+
mHadRapidityMin = yMin;
68+
mHadRapidityMax = yMax;
69+
};
70+
void setUsedSeed(unsigned int seed)
71+
{
72+
mUsedSeed = seed;
73+
};
74+
unsigned int getUsedSeed() const
75+
{
76+
return mUsedSeed;
77+
};
78+
79+
protected:
80+
//__________________________________________________________________
81+
bool generateEvent() override
82+
{
83+
// Simple straightforward check to alternate generators
84+
if (mGeneratedEvents % mInverseTriggerRatio == 0) {
85+
int nInjectedEvents = mGeneratedEvents / mInverseTriggerRatio;
86+
// Alternate hadrons if enabled (with the same ratio)
87+
if (mHFHadronPdgList.size() >= 1) {
88+
int iHadron = nInjectedEvents % mHFHadronPdgList.size();
89+
mHadronPdg = mHFHadronPdgList[iHadron];
90+
LOG(info) << "Selected hadron: " << mHFHadronPdgList[iHadron];
91+
}
92+
93+
// Generate event of interest
94+
bool genOk = false;
95+
while (!genOk) {
96+
if (GeneratorPythia8::generateEvent()) {
97+
genOk = selectEvent(mPythia.event);
98+
}
99+
}
100+
notifySubGenerator(1);
101+
} else {
102+
// Generate minimum-bias event
103+
bool genOk = false;
104+
while (!genOk) {
105+
genOk = GeneratorPythia8::generateEvent();
106+
}
107+
notifySubGenerator(0);
108+
}
109+
110+
mGeneratedEvents++;
111+
112+
return true;
113+
}
114+
115+
bool selectEvent(Pythia8::Event& event)
116+
{
117+
for (auto iPart{0}; iPart < event.size(); ++iPart) {
118+
// search for hadron in rapidity window
119+
int id = std::abs(event[iPart].id());
120+
float rap = event[iPart].y();
121+
if (id == mHadronPdg && rap > mHadRapidityMin && rap < mHadRapidityMax) {
122+
LOG(debug) << "-----------------------------------------------------";
123+
LOG(debug) << "Found hadron " << event[iPart].id() << " with rapidity " << rap << " and daughters " << event[iPart].daughter1() << " " << event[iPart].daughter2();
124+
// print pdg code of daughters
125+
LOG(debug) << "Daughters: ";
126+
for (int iDau = event[iPart].daughter1(); iDau <= event[iPart].daughter2(); ++iDau) {
127+
LOG(debug) << "Daughter " << iDau << ": " << event[iDau].id();
128+
}
129+
bool isCoalDone = CoalescencePythia8(event, mNucleiPdgList, mTrivialCoal, mCoalMomentum, event[iPart].daughter1(), event[iPart].daughter2());
130+
if (isCoalDone) {
131+
LOG(debug) << "Coalescence process found for hadron " << event[iPart].id() << " with daughters " << event[iPart].daughter1() << " " << event[iPart].daughter2();
132+
LOG(debug) << "Check updated daughters: ";
133+
for (int iDau = event[iPart].daughter1(); iDau <= event[iPart].daughter2(); ++iDau) {
134+
LOG(debug) << "Daughter " << iDau << ": " << event[iDau].id();
135+
}
136+
return true;
137+
}
138+
}
139+
}
140+
return false;
141+
};
142+
143+
private:
144+
// Interface to override import particles
145+
Pythia8::Event mOutputEvent;
146+
147+
// Properties of selection
148+
int mHadronPdg;
149+
float mHadRapidityMin;
150+
float mHadRapidityMax;
151+
unsigned int mUsedSeed;
152+
153+
// Control gap-triggering
154+
unsigned long long mGeneratedEvents;
155+
int mInverseTriggerRatio;
156+
157+
// Control alternate trigger on different hadrons
158+
std::vector<int> mHFHadronPdgList = {};
159+
std::vector<unsigned int> mNucleiPdgList = {};
160+
161+
bool mTrivialCoal = false; /// if true, the coalescence is done without checking the distance in the phase space of the nucleons
162+
float mCoalMomentum; /// coalescence momentum
163+
};
164+
165+
166+
///___________________________________________________________
167+
FairGenerator *generateHFHadToNuclei(int input_trigger_ratio = 5, std::vector<int> hf_hadron_pdg_list = {}, std::vector<unsigned int> nuclei_pdg_list = {}, bool trivial_coal = false, float coal_momentum = 0.4)
168+
{
169+
auto myGen = new GeneratorPythia8HFHadToNuclei(input_trigger_ratio, hf_hadron_pdg_list, nuclei_pdg_list, trivial_coal, coal_momentum);
170+
auto seed = (gRandom->TRandom::GetSeed() % 900000000);
171+
myGen->readString("Random:setSeed on");
172+
myGen->readString("Random:seed " + std::to_string(seed));
173+
return myGen;
174+
}

‎MC/config/PWGLF/ini/GeneratorLF_Coalescence.ini‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,7 @@
11
[GeneratorExternal]
22
fileName = ${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGLF/pythia8/generator_pythia8_coalescence.C
3-
funcName = generateCoalescence(1, 0.239)
3+
funcName = generateCoalescence({1000010030, 1000020030, 1010010030}, 1, 0.239)
44

55
[GeneratorPythia8]
66
config=${O2_ROOT}/share/Generators/egconfig/pythia8_inel.cfg
7+

‎MC/config/PWGLF/pythia8/generator_pythia8_coalescence.C‎

Lines changed: 12 additions & 89 deletions
Original file line numberDiff line numberDiff line change
@@ -16,19 +16,21 @@
1616
using namespace Pythia8;
1717
#endif
1818

19+
R__ADD_INCLUDE_PATH($O2DPG_MC_CONFIG_ROOT)
20+
#include "MC/config/common/external/generator/CoalescencePythia8.h"
1921
/// First version of the simple coalescence generator based PYTHIA8
20-
/// TODO: extend to other nuclei (only He3 is implemented now)
2122

2223
class GeneratorPythia8Coalescence : public o2::eventgen::GeneratorPythia8
2324
{
2425
public:
2526
/// Constructor
26-
GeneratorPythia8Coalescence(int input_trigger_ratio = 1, double coal_momentum = 0.4)
27+
GeneratorPythia8Coalescence(std::vector<unsigned int> pdgList, int input_trigger_ratio = 1, double coal_momentum = 0.4)
2728
: o2::eventgen::GeneratorPythia8()
2829
{
2930
fmt::printf(">> Coalescence generator %d\n", input_trigger_ratio);
3031
mInverseTriggerRatio = input_trigger_ratio;
3132
mCoalMomentum = coal_momentum;
33+
mPdgList = pdgList;
3234
}
3335
/// Destructor
3436
~GeneratorPythia8Coalescence() = default;
@@ -47,13 +49,14 @@ protected:
4749
// Simple straightforward check to alternate generators
4850
if (mGeneratedEvents % mInverseTriggerRatio == 0)
4951
{
52+
fmt::printf(">> Generating coalescence event %d\n", mGeneratedEvents);
5053
bool genOk = false;
5154
int localCounter{0};
5255
while (!genOk)
5356
{
5457
if (GeneratorPythia8::generateEvent())
5558
{
56-
genOk = selectEvent(mPythia.event);
59+
genOk = CoalescencePythia8(mPythia.event, mPdgList, mCoalMomentum);
5760
}
5861
localCounter++;
5962
}
@@ -63,6 +66,7 @@ protected:
6366
}
6467
else
6568
{
69+
fmt::printf(">> Generating minimum-bias event %d\n", mGeneratedEvents);
6670
// Generate minimum-bias event
6771
bool genOk = false;
6872
while (!genOk)
@@ -71,106 +75,25 @@ protected:
7175
}
7276
notifySubGenerator(0);
7377
}
74-
7578
mGeneratedEvents++;
76-
7779
return true;
7880
}
7981

80-
bool selectEvent(Pythia8::Event &event)
81-
{
82-
std::vector<int> protons[2], neutrons[2], lambdas[2];
83-
for (auto iPart{0}; iPart < event.size(); ++iPart)
84-
{
85-
if (std::abs(event[iPart].y()) > 1.) // skip particles with y > 1
86-
{
87-
continue;
88-
}
89-
switch (std::abs(event[iPart].id()))
90-
{
91-
case 2212:
92-
protons[event[iPart].id() > 0].push_back(iPart);
93-
break;
94-
case 2112:
95-
neutrons[event[iPart].id() > 0].push_back(iPart);
96-
break;
97-
case 3122:
98-
lambdas[event[iPart].id() > 0].push_back(iPart);
99-
break;
100-
default:
101-
break;
102-
}
103-
}
104-
const double coalescenceRadius{0.5 * 1.122462 * mCoalMomentum}; /// 1.122462 [2^(1/6)] from PRL 126, 101101 (2021), only for 3 body coalescence
105-
106-
auto coalescence = [&](int iC, int pdgCode, float mass, int iD1, int iD2, int iD3) {
107-
if (event[iD1].status() < 0 || event[iD2].status() < 0 || event[iD3].status() < 0)
108-
{
109-
return false;
110-
}
111-
auto p1 = event[iD1].p();
112-
auto p2 = event[iD2].p();
113-
auto p3 = event[iD3].p();
114-
auto p = p1 + p2 + p3;
115-
p1.bstback(p);
116-
p2.bstback(p);
117-
p3.bstback(p);
118-
119-
if (p1.pAbs() <= coalescenceRadius && p2.pAbs() <= coalescenceRadius && p3.pAbs() <= coalescenceRadius)
120-
{
121-
p.e(std::hypot(p.pAbs(), mass));
122-
/// In order to avoid the transport of the mother particles, but to still keep them in the stack, we set the status to negative and we mark the nucleus status as 94 (decay product)
123-
event.append((iC * 2 - 1) * pdgCode, 94, 0, 0, 0, 0, 0, 0, p.px(), p.py(), p.pz(), p.e(), mass);
124-
event[iD1].statusNeg();
125-
event[iD1].daughter1(event.size() - 1);
126-
event[iD2].statusNeg();
127-
event[iD2].daughter1(event.size() - 1);
128-
event[iD3].statusNeg();
129-
event[iD3].daughter1(event.size() - 1);
130-
131-
fmt::printf(">> Adding a %i with p = %f, %f, %f, E = %f\n", (iC * 2 - 1) * pdgCode, p.px(), p.py(), p.pz(), p.e());
132-
133-
return true;
134-
}
135-
return false;
136-
};
137-
138-
bool coalHappened = false;
139-
for (int iC{0}; iC < 2; ++iC)
140-
{
141-
for (int iP{0}; iP < protons[iC].size(); ++iP) {
142-
for (int iN{0}; iN < neutrons[iC].size(); ++iN) {
143-
/// H3L loop
144-
for (int iL{0}; iL < lambdas[iC].size(); ++iL) {
145-
coalHappened |= coalescence(iC, 1010010030, 2.991134, protons[iC][iP], neutrons[iC][iN], lambdas[iC][iL]);
146-
}
147-
/// H3 loop
148-
for (int iN2{iN + 1}; iN2 < neutrons[iC].size(); ++iN2) {
149-
coalHappened |= coalescence(iC, 1000010030, 2.80892113298, protons[iC][iP], neutrons[iC][iN], neutrons[iC][iN2]);
150-
}
151-
/// He3 loop
152-
for (int iP2{iP + 1}; iP2 < protons[iC].size(); ++iP2) {
153-
coalHappened |= coalescence(iC, 1000020030, 2.808391, protons[iC][iP], protons[iC][iP2], neutrons[iC][iN]);
154-
}
155-
}
156-
}
157-
}
158-
return coalHappened;
159-
}
16082

16183
private:
162-
// Control gap-triggering
163-
float mCoalMomentum = 0.4;
84+
std::vector<unsigned int> mPdgList; /// list of pdg codes to be generated
85+
float mCoalMomentum = 0.4; /// coalescence momentum
16486
uint64_t mGeneratedEvents = 0; /// number of events generated so far
16587
int mInverseTriggerRatio = 1; /// injection gap
16688
};
16789

16890
///___________________________________________________________
169-
FairGenerator *generateCoalescence(int input_trigger_ratio, double coal_momentum = 0.4)
91+
FairGenerator *generateCoalescence(std::vector<unsigned int> pdgList, int input_trigger_ratio, double coal_momentum = 0.4)
17092
{
171-
auto myGen = new GeneratorPythia8Coalescence(input_trigger_ratio, coal_momentum);
93+
auto myGen = new GeneratorPythia8Coalescence(pdgList, input_trigger_ratio, coal_momentum);
17294
auto seed = (gRandom->TRandom::GetSeed() % 900000000);
17395
myGen->readString("Random:setSeed on");
17496
myGen->readString("Random:seed " + std::to_string(seed));
17597
return myGen;
17698
}
99+
///___________________________________________________________

0 commit comments

Comments
 (0)