Skip to content

Commit 06d6e85

Browse files
alcalivaalcaliva
committed
Add hybrid genberator for antideuterons in jets
1 parent 6ab618e commit 06d6e85

4 files changed

Lines changed: 229 additions & 0 deletions

File tree

‎MC/config/PWGLF/hybrid/hybrid_coalescence.json‎

Lines changed: 24 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,24 @@
1+
2+
{
3+
"mode": "sequential",
4+
"generators": [
5+
{
6+
"name": "pythia8",
7+
"config": {
8+
"config": "${O2_ROOT}/share/Generators/egconfig/pythia8_inel.cfg",
9+
"hooksFileName": "",
10+
"hooksFuncName": "",
11+
"includePartonEvent": false,
12+
"particleFilter": "",
13+
"verbose": 0
14+
}
15+
},
16+
{
17+
"name": "external",
18+
"config": {
19+
"iniFile": "${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGLF/ini/GeneratorAntideuteronCoalescenceJet.ini"
20+
}
21+
}
22+
],
23+
"fractions": [4, 1]
24+
}

‎MC/config/PWGLF/ini/GeneratorAntideuteronCoalescenceJet.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_antideuteron_injection.C
3+
funcName = generateAntideuterons(0.8)
4+
5+
[GeneratorPythia8]
6+
config = ${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGGAJE/pythia8/generator/pythia8_jet.cfg

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

Lines changed: 85 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,85 @@
1+
2+
#if !defined(__CLING__) || defined(__ROOTCLING__)
3+
#include "SimulationDataFormat/MCTrack.h"
4+
#include "TFile.h"
5+
#include "TTree.h"
6+
7+
#include <algorithm>
8+
#include <iostream>
9+
#include <string>
10+
#include <vector>
11+
#endif // !defined(__CLING__) || defined(__ROOTCLING__)
12+
13+
int External()
14+
{
15+
std::string path{"o2sim_Kine.root"};
16+
17+
std::vector<int> injectedPDGs = {
18+
-1000010020 // anti-deuteron
19+
};
20+
21+
TFile file(path.c_str(), "READ");
22+
if (file.IsZombie()) {
23+
std::cerr << "Cannot open ROOT file " << path << "\n";
24+
return 1;
25+
}
26+
27+
auto* tree = (TTree*)file.Get("o2sim");
28+
if (!tree) {
29+
std::cerr << "Cannot find tree o2sim in file " << path << "\n";
30+
return 1;
31+
}
32+
33+
std::vector<o2::MCTrack>* tracks = nullptr;
34+
tree->SetBranchAddress("MCTrack", &tracks);
35+
36+
std::vector<int> nSignal(injectedPDGs.size(), 0);
37+
std::vector<int> nEventsWithSignal(injectedPDGs.size(), 0);
38+
39+
const auto nEvents = tree->GetEntries();
40+
41+
for (Long64_t i = 0; i < nEvents; ++i) {
42+
tree->GetEntry(i);
43+
44+
std::vector<bool> found(injectedPDGs.size(), false);
45+
46+
for (const auto& track : *tracks) {
47+
const int pdg = track.GetPdgCode();
48+
49+
auto it = std::find(
50+
injectedPDGs.begin(),
51+
injectedPDGs.end(),
52+
pdg);
53+
54+
if (it != injectedPDGs.end()) {
55+
const auto index = std::distance(injectedPDGs.begin(), it);
56+
++nSignal[index];
57+
found[index] = true;
58+
}
59+
}
60+
61+
for (size_t j = 0; j < injectedPDGs.size(); ++j) {
62+
if (found[j]) {
63+
++nEventsWithSignal[j];
64+
}
65+
}
66+
}
67+
68+
std::cout << "--------------------------------\n";
69+
std::cout << "# Events: " << nEvents << "\n";
70+
71+
for (size_t i = 0; i < injectedPDGs.size(); ++i) {
72+
std::cout << "PDG: " << injectedPDGs[i] << "\n";
73+
std::cout << "# Generated particles: " << nSignal[i] << "\n";
74+
std::cout << "# Events containing particle: "
75+
<< nEventsWithSignal[i] << "\n";
76+
77+
if (nSignal[i] == 0) {
78+
std::cerr << "No generated particles with PDG "
79+
<< injectedPDGs[i] << "\n";
80+
return 1;
81+
}
82+
}
83+
84+
return 0;
85+
}

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

Lines changed: 114 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,114 @@
1+
#if !defined(__CLING__) || defined(__ROOTCLING__)
2+
3+
#include "FairGenerator.h"
4+
#include "Generators/GeneratorPythia8.h"
5+
#include "Pythia8/Pythia.h"
6+
#include "TRandom.h"
7+
#include "fairlogger/Logger.h"
8+
9+
#include <cmath>
10+
#include <cstdint>
11+
#include <string>
12+
13+
#endif // !defined(__CLING__) || defined(__ROOTCLING__)
14+
15+
/// PYTHIA8 event generator with antideuteron injection.
16+
/// Events are generated until a final-state antiproton within the specified etaMax is found.
17+
/// The antiproton is then replaced by an antideuteron with twice its momentum.
18+
19+
class GeneratorPythia8AntideuteronInjection
20+
: public o2::eventgen::GeneratorPythia8 {
21+
22+
public:
23+
24+
GeneratorPythia8AntideuteronInjection(double etaMax = 0.8)
25+
: o2::eventgen::GeneratorPythia8(),
26+
mEtaMax(etaMax) {}
27+
28+
~GeneratorPythia8AntideuteronInjection() override = default;
29+
30+
bool Init() override {
31+
addSubGenerator(0, "Pythia8 events with antideuterons");
32+
return o2::eventgen::GeneratorPythia8::Init();
33+
}
34+
35+
protected:
36+
37+
38+
bool generateEvent() override {
39+
constexpr int maxAttempts = 10000;
40+
for (int attempt = 0; attempt < maxAttempts; ++attempt) {
41+
if (!GeneratorPythia8::generateEvent()) {
42+
return false;
43+
}
44+
45+
if (EventHasAntideuteron(mPythia.event)) {
46+
LOG(debug) << ">> An antideuteron was formed!";
47+
48+
notifySubGenerator(0);
49+
++mGeneratedEvents;
50+
51+
return true;
52+
}
53+
}
54+
55+
LOG(error) << "Failed to generate an antideuteron after " << maxAttempts << " attempts";
56+
return false;
57+
}
58+
59+
bool EventHasAntideuteron(Pythia8::Event& event) {
60+
61+
const double md = 1.87561294257; // Deuteron mass in GeV/c^2
62+
const int pdgAntid = -1000010020; // Antideuteron PDG code
63+
64+
for (int iPart = 0; iPart < event.size(); ++iPart) {
65+
66+
if (!event[iPart].isFinal()) {
67+
continue;
68+
}
69+
70+
if (std::abs(event[iPart].eta()) > mEtaMax) {
71+
continue;
72+
}
73+
74+
if (event[iPart].id() != -2212) {
75+
continue;
76+
}
77+
78+
double px = 2.0 * event[iPart].px();
79+
double py = 2.0 * event[iPart].py();
80+
double pz = 2.0 * event[iPart].pz();
81+
double energy = std::sqrt(md * md + px * px + py * py + pz * pz);
82+
83+
// Remove antiproton from final state
84+
event[iPart].statusNeg();
85+
86+
// Add antideuteron
87+
const int iAntid = event.size();
88+
event.append(pdgAntid, 121, iPart, 0, 0, 0, 0, 0, px, py, pz, energy, md);
89+
// Update daughter indices
90+
event[iPart].daughter1(iAntid);
91+
event[iPart].daughter2(iAntid);
92+
return true;
93+
}
94+
return false;
95+
}
96+
97+
private:
98+
99+
double mEtaMax = 0.8;
100+
uint64_t mGeneratedEvents = 0;
101+
};
102+
103+
/// Generator factory
104+
105+
FairGenerator* generateAntideuterons(double etaMax = 0.8) {
106+
107+
auto* myGenerator = new GeneratorPythia8AntideuteronInjection(etaMax);
108+
auto seed = gRandom->TRandom::GetSeed() % 900000000;
109+
110+
myGenerator->readString("Random:setSeed = on");
111+
myGenerator->readString("Random:seed = " + std::to_string(seed));
112+
113+
return myGenerator;
114+
}

0 commit comments

Comments
 (0)