Skip to content

Commit 94daa12

Browse files
mpucciompuccio
committed
Generalise coalescence generator in view of event pools
1 parent e40ee48 commit 94daa12

File tree

3 files changed

+65
-48
lines changed

3 files changed

+65
-48
lines changed

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

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

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

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

Lines changed: 4 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -20,10 +20,12 @@ int External()
2020

2121
auto nEvents = tree->GetEntries();
2222
auto nInjected = tree->Scan("MCTrack.GetPdgCode()", "MCTrack.GetPdgCode() == -1000020030"); /// don't check matter, too many secondaries
23-
if (nEvents / nInjected != 10)
23+
nInjected += tree->Scan("MCTrack.GetPdgCode()", "MCTrack.GetPdgCode() == -1000010030"); /// don't check matter, too many secondaries
24+
nInjected += tree->Scan("MCTrack.GetPdgCode()", "TMath::Abs(MCTrack.GetPdgCode()) == 1010010030"); /// don't check matter, too many secondaries
25+
if (nInjected == 0)
2426
{
2527
std::cerr << "Unexpected ratio of events to injected nuclei\n";
2628
return 1;
2729
}
2830
return 0;
29-
}
31+
}

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

Lines changed: 60 additions & 45 deletions
Original file line numberDiff line numberDiff line change
@@ -79,64 +79,79 @@ protected:
7979

8080
bool selectEvent(Pythia8::Event &event)
8181
{
82-
static int sign = -1; // start with antimatter
83-
std::vector<int> protons, neutrons, lambdas;
82+
std::vector<int> protons[2], neutrons[2], lambdas[2];
8483
for (auto iPart{0}; iPart < event.size(); ++iPart)
8584
{
8685
if (std::abs(event[iPart].y()) > 1.) // skip particles with y > 1
8786
{
8887
continue;
8988
}
90-
if (event[iPart].id() == 2212 * sign)
89+
switch (std::abs(event[iPart].id()))
9190
{
92-
protons.push_back(iPart);
93-
}
94-
else if (event[iPart].id() == 2112 * sign)
95-
{
96-
neutrons.push_back(iPart);
97-
}
98-
else if (event[iPart].id() == 3122 * sign)
99-
{
100-
lambdas.push_back(iPart);
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;
101102
}
102103
}
103-
double coalescenceRadius{0.5 * 1.122462 * mCoalMomentum}; /// 1.122462 [2^(1/6)] from PRL 126, 101101 (2021), only for 3 body coalescence
104-
if (protons.size() < 2 || neutrons.size() < 1) // at least 2 protons and 1 neutron
105-
{
106-
return false;
107-
}
108-
for (uint32_t i{0}; i < protons.size(); ++i)
109-
{
110-
for (uint32_t j{i + 1}; j < protons.size(); ++j)
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+
auto p1 = event[iD1].p();
108+
auto p2 = event[iD2].p();
109+
auto p3 = event[iD3].p();
110+
auto p = p1 + p2 + p3;
111+
p1.bstback(p);
112+
p2.bstback(p);
113+
p3.bstback(p);
114+
115+
if (p1.pAbs() <= coalescenceRadius && p2.pAbs() <= coalescenceRadius && p3.pAbs() <= coalescenceRadius)
111116
{
112-
for (uint32_t k{0}; k < neutrons.size(); ++k)
113-
{
114-
auto p1 = event[protons[i]].p();
115-
auto p2 = event[protons[j]].p();
116-
auto p3 = event[neutrons[k]].p();
117-
auto p = p1 + p2 + p3;
118-
p1.bstback(p);
119-
p2.bstback(p);
120-
p3.bstback(p);
117+
p.e(std::hypot(p.pAbs(), mass));
118+
/// 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)
119+
event.append((iC * 2 - 1) * pdgCode, 94, 0, 0, 0, 0, 0, 0, p.px(), p.py(), p.pz(), p.e(), mass);
120+
event[iD1].statusNeg();
121+
event[iD1].daughter1(event.size() - 1);
122+
event[iD2].statusNeg();
123+
event[iD2].daughter1(event.size() - 1);
124+
event[iD3].statusNeg();
125+
event[iD3].daughter1(event.size() - 1);
121126

122-
if (p1.pAbs() <= coalescenceRadius && p2.pAbs() <= coalescenceRadius && p3.pAbs() <= coalescenceRadius)
123-
{
124-
p.e(std::hypot(p.pAbs(), 2.80839160743));
125-
/// 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)
126-
event.append(sign * 1000020030, 94, 0, 0, 0, 0, 0, 0, p.px(), p.py(), p.pz(), p.e(), 2.80839160743);
127-
event[protons[i]].statusNeg();
128-
event[protons[i]].daughter1(event.size() - 1);
129-
event[protons[j]].statusNeg();
130-
event[protons[j]].daughter1(event.size() - 1);
131-
event[neutrons[k]].statusNeg();
132-
event[neutrons[k]].daughter1(event.size() - 1);
127+
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());
133128

134-
fmt::printf(">> Adding a He3 with p = %f, %f, %f, E = %f\n", p.px(), p.py(), p.pz(), p.e());
135-
std::cout << std::endl;
136-
std::cout << std::endl;
129+
return true;
130+
}
131+
return false;
132+
};
137133

138-
sign *= -1;
139-
return true;
134+
for (int iC{0}; iC < 2; ++iC)
135+
{
136+
for (int iP{0}; iP < protons[iC].size(); ++iP) {
137+
for (int iN{0}; iN < neutrons[iC].size(); ++iN) {
138+
/// H3L loop
139+
for (int iL{0}; iL < lambdas[iC].size(); ++iL) {
140+
if (coalescence(iC, 1010010030, 2.991134, protons[iC][iP], neutrons[iC][iN], lambdas[iC][iL])) {
141+
return true;
142+
}
143+
}
144+
/// H3 loop
145+
for (int iN2{iN + 1}; iN2 < neutrons[iC].size(); ++iN2) {
146+
if (coalescence(iC, 1000010030, 2.80892113298, protons[iC][iP], neutrons[iC][iN], neutrons[iC][iN2])) {
147+
return true;
148+
}
149+
}
150+
/// He3 loop
151+
for (int iP2{iP + 1}; iP2 < protons[iC].size(); ++iP2) {
152+
if (coalescence(iC, 1000020030, 2.808391, protons[iC][iP], protons[iC][iP2], neutrons[iC][iN])) {
153+
return true;
154+
}
140155
}
141156
}
142157
}

0 commit comments

Comments
 (0)