Skip to content

Commit 8edc590

Browse files
alcalivaalcaliva
authored
[PWGLF] Added deuteron generator using advanced coalescence (#2019)
1 parent ef82e47 commit 8edc590

File tree

3 files changed

+220
-0
lines changed

3 files changed

+220
-0
lines changed

‎MC/config/PWGLF/ini/GeneratorLF_deuteron_wigner.ini‎

Lines changed: 6 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,6 @@
1+
[GeneratorExternal]
2+
fileName = ${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGLF/pythia8/generator_pythia8_deuteron_wigner.C
3+
funcName = generateAntideuteronsWignerCoalescence(1.2)
4+
5+
[GeneratorPythia8]
6+
config=${O2_ROOT}/share/Generators/egconfig/pythia8_inel.cfg

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

Lines changed: 41 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,41 @@
1+
2+
int External() {
3+
std::string path{"o2sim_Kine.root"};
4+
5+
TFile file(path.c_str(), "READ");
6+
if (file.IsZombie()) {
7+
std::cerr << "Cannot open ROOT file " << path << "\n";
8+
return 1;
9+
}
10+
11+
auto tree = (TTree *)file.Get("o2sim");
12+
if (!tree) {
13+
std::cerr << "Cannot find tree o2sim in file " << path << "\n";
14+
return 1;
15+
}
16+
17+
std::vector<o2::MCTrack>* tracks = nullptr;
18+
tree->SetBranchAddress("MCTrack", &tracks);
19+
20+
Long64_t nEntries = tree->GetEntries();
21+
int nSelected = 0;
22+
23+
for (Long64_t i = 0; i < nEntries; ++i) {
24+
tree->GetEntry(i);
25+
for (const auto& track : *tracks) {
26+
if (std::abs(track.GetPdgCode()) == 1000010020) { // Deuteron
27+
++nSelected;
28+
break; // Found in this event
29+
}
30+
}
31+
}
32+
33+
if (nSelected == 0) {
34+
std::cerr << "No event of interest\n";
35+
return 1;
36+
}
37+
38+
std::cout << "Found " << nSelected << " events with deuterons\n";
39+
return 0;
40+
}
41+

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

Lines changed: 173 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,173 @@
1+
#if !defined(__CLING__) || defined(__ROOTCLING__)
2+
#include "FairGenerator.h"
3+
#include "FairPrimaryGenerator.h"
4+
#include "Generators/GeneratorPythia8.h"
5+
#include "Pythia8/Pythia.h"
6+
#include "TDatabasePDG.h"
7+
#include "TMath.h"
8+
#include "TParticlePDG.h"
9+
#include "TRandom3.h"
10+
#include "TSystem.h"
11+
#include "TH2D.h"
12+
#include "TH1D.h"
13+
#include "TFile.h"
14+
#include "fairlogger/Logger.h"
15+
#include "CCDB/BasicCCDBManager.h"
16+
#include <cmath>
17+
#include <fstream>
18+
#include <string>
19+
#include <vector>
20+
using namespace Pythia8;
21+
#endif
22+
23+
/// Event generator for proton-proton (pp) collisions using Pythia8.
24+
/// (Anti)deuterons are formed via nucleon coalescence modeled using the Wigner density formalism.
25+
26+
class GeneratorPythia8DeuteronWigner : public o2::eventgen::GeneratorPythia8 {
27+
public:
28+
GeneratorPythia8DeuteronWigner(double sourceRadius = 1.2)
29+
: o2::eventgen::GeneratorPythia8(), mSourceRadius(sourceRadius) {
30+
31+
// Connect to CCDB and retrieve coalescence probability two-dimensional table
32+
o2::ccdb::CcdbApi ccdb_api;
33+
ccdb_api.init("https://alice-ccdb.cern.ch");
34+
TFile* file = ccdb_api.retrieveFromTFileAny<TFile>("Users/a/alcaliva/WignerCoalescence/ArgonneProbability.root");
35+
36+
if (!file) {
37+
LOG(fatal) << "Could not retrieve ArgonneProbability.root from CCDB!";
38+
}
39+
40+
mTwoDimCoalProbability = dynamic_cast<TH2D*>(file->FindObject("AddedSDWave"));
41+
if (!mTwoDimCoalProbability) {
42+
LOG(fatal) << "Could not find 'AddedSDWave' histogram in the input file!";
43+
}
44+
}
45+
46+
~GeneratorPythia8DeuteronWigner() override = default;
47+
48+
bool Init() override {
49+
addSubGenerator(0, "Pythia8 with (anti)deuterons formed via coalescence using the Wigner density formalism");
50+
return o2::eventgen::GeneratorPythia8::Init();
51+
}
52+
53+
protected:
54+
bool generateEvent() override {
55+
if (mGeneratedEvents % 100 == 0) {
56+
LOG(info) << ">> Generating event " << mGeneratedEvents;
57+
}
58+
59+
bool genOk = false;
60+
int localCounter = 0;
61+
while (!genOk) {
62+
if (GeneratorPythia8::generateEvent()) {
63+
genOk = selectEvent(mPythia.event);
64+
}
65+
localCounter++;
66+
}
67+
68+
LOG(debug) << ">> Generation of event of interest successful after " << localCounter << " iterations";
69+
notifySubGenerator(0);
70+
++mGeneratedEvents;
71+
return true;
72+
}
73+
74+
bool selectEvent(Pythia8::Event& event) {
75+
76+
const double md = 1.87561294257; // Deuteron mass [GeV]
77+
bool deuteronIsFormed = false;
78+
79+
std::vector<int> proton_ID, neutron_ID;
80+
std::vector<int> proton_status, neutron_status;
81+
82+
for (int iPart = 0; iPart < event.size(); ++iPart) {
83+
if (event[iPart].status() <= 0) {
84+
continue;
85+
}
86+
87+
int absID = std::abs(event[iPart].id());
88+
if (absID == 2212) {
89+
proton_ID.push_back(iPart);
90+
proton_status.push_back(0);
91+
} else if (absID == 2112) {
92+
neutron_ID.push_back(iPart);
93+
neutron_status.push_back(0);
94+
}
95+
}
96+
97+
if (proton_ID.empty() || neutron_ID.empty()) {
98+
return false;
99+
}
100+
101+
int radiusBin = mTwoDimCoalProbability->GetXaxis()->FindBin(mSourceRadius);
102+
TH1D* prob_vs_q = mTwoDimCoalProbability->ProjectionY("prob_vs_q", radiusBin, radiusBin, "E");
103+
prob_vs_q->SetDirectory(nullptr);
104+
105+
for (size_t ip = 0; ip < proton_ID.size(); ip++) {
106+
if (proton_status[ip] < 0) continue;
107+
108+
for (size_t in = 0; in < neutron_ID.size(); in++) {
109+
if (neutron_status[in] < 0) continue;
110+
111+
int protonID = proton_ID[ip];
112+
int neutronID = neutron_ID[in];
113+
int sign_p = event[protonID].id() / std::abs(event[protonID].id());
114+
int sign_n = event[neutronID].id() / std::abs(event[neutronID].id());
115+
if (sign_p != sign_n) continue;
116+
117+
auto p1 = event[protonID].p();
118+
auto p2 = event[neutronID].p();
119+
auto p = p1 + p2;
120+
p1.bstback(p);
121+
p2.bstback(p);
122+
123+
Vec4 deltaPVec = p1 - p2;
124+
double deltaP = 0.5 * deltaPVec.pAbs();
125+
126+
int binQ = prob_vs_q->FindBin(deltaP);
127+
128+
// Skip underflow and overflow bins
129+
if (binQ < 1 || binQ > prob_vs_q->GetNbinsX()) {
130+
continue;
131+
}
132+
double coalProb = prob_vs_q->GetBinContent(prob_vs_q->FindBin(deltaP));
133+
double rnd = gRandom->Uniform(0.0, 1.0);
134+
135+
if (rnd < coalProb) {
136+
double energy = std::hypot(p.pAbs(), md);
137+
p.e(energy);
138+
int deuteronPDG = sign_p * 1000010020;
139+
140+
event.append(deuteronPDG, 121, 0, 0, 0, 0, 0, 0, p.px(), p.py(), p.pz(), p.e(), md);
141+
142+
event[protonID].statusNeg();
143+
event[protonID].daughter1(event.size() - 1);
144+
proton_status[ip] = -1;
145+
146+
event[neutronID].statusNeg();
147+
event[neutronID].daughter1(event.size() - 1);
148+
neutron_status[in] = -1;
149+
150+
deuteronIsFormed = true;
151+
}
152+
}
153+
}
154+
155+
// free allocated memory
156+
delete prob_vs_q;
157+
return deuteronIsFormed;
158+
}
159+
160+
private:
161+
double mSourceRadius = 1.2;
162+
uint64_t mGeneratedEvents = 0;
163+
TH2D* mTwoDimCoalProbability = nullptr;
164+
};
165+
166+
///________________________________________________________________________________________________________
167+
FairGenerator* generateAntideuteronsWignerCoalescence(double sourceRadius = 1.2) {
168+
auto myGenerator = new GeneratorPythia8DeuteronWigner(sourceRadius);
169+
auto seed = gRandom->TRandom::GetSeed() % 900000000;
170+
myGenerator->readString("Random:setSeed on");
171+
myGenerator->readString("Random:seed " + std::to_string(seed));
172+
return myGenerator;
173+
}

0 commit comments

Comments
 (0)