-
Notifications
You must be signed in to change notification settings - Fork 0
/
pythia8.C
110 lines (86 loc) · 2.93 KB
/
pythia8.C
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
#include "TTree.h"
#include "TFile.h"
#include <vector>
#include <chrono>
void pythia8(Int_t nev = 100)
{
auto start = std::chrono::high_resolution_clock::now();
// Load libraries
//gSystem->Load("pythia/lib/libpythia8.so");
//gSystem->Load("pythia/examples/main92.so");
//gSystem->Load("libEG");
//gSystem->Load("libEGPythia8");
// Create pythia8 object
Pythia8::Pythia *pythia8 = new Pythia8::Pythia();
// File
TFile* file = TFile::Open("pytree.root","recreate");
// Prepare Tree
Int_t np;
std::vector<Int_t> id;
std::vector<Int_t> status;
std::vector<Float_t> e;
std::vector<Float_t> m;
std::vector<Float_t> pt;
std::vector<Float_t> phi;
std::vector<Float_t> eta;
TTree* tree = new TTree("tree","tree");
auto np_branch = tree->Branch("np",&np);
auto id_branch = tree->Branch("id",&id);
auto status_branch = tree->Branch("status",&status);
auto e_branch = tree->Branch("e",&e);
auto m_branch = tree->Branch("m",&m);
auto pt_branch = tree->Branch("pt",&pt);
auto phi_branch = tree->Branch("phi",&phi);
auto eta_branch = tree->Branch("eta",&eta);
// Configure
//pythia8->readString("HardQCD:all = on");
//pythia8->readString("HiggsSM:gg2H = on");
// pythia8->readString("25:onMode = off");
// pythia8->readString("25:onIfAll= 11 -11");
// pythia8->readString("Random:setSeed = on");
// use a reproducible seed: always the same results for the tutorial.
// pythia8->readString("Random:seed = 42");
pythia8->readString("Beams:eCM = 13000.");
pythia8->readString("Top:gg2ttbar = on");
pythia8->readString("Top:qqbar2ttbar = on");
// pythia8->readString("PartonLevel:FSR = off");
// pythia8->readString("PartonLevel:ISR = off");
// pythia8->readString("PartonLevel:MPI = off");
// pythia8->readString("HadronLevel:all = off");
// pythia8->readString("PartonLevel:all = off");
// Initialize
//pythia8->init(2212 /* p */, 2212 /* p */, 14000. /* TeV */);
pythia8->init();
// Event loop
for (Int_t iev = 0; iev < nev; iev++) {
status.clear();
id.clear();
e.clear();
m.clear();
pt.clear();
phi.clear();
eta.clear();
pythia8->next();
np = pythia8->event.size();
// std::cout << "Event Nr.: " << iev << std::endl;
// std::cout << "Nr. of Particles: " << np << std::endl;
//Particle loop
for (Int_t ip = 0; ip < np; ip++) {
status.push_back(pythia8->event[ip].status());
id.push_back(pythia8->event[ip].id());
e.push_back(pythia8->event[ip].e());
m.push_back(pythia8->event[ip].m());
pt.push_back(pythia8->event[ip].pT());
phi.push_back(pythia8->event[ip].phi());
eta.push_back(pythia8->event[ip].eta());
}
tree->Fill();
}
pythia8->stat();
tree->Print();
tree->Write();
delete file;
auto stop = std::chrono::high_resolution_clock::now();
auto duration = std::chrono::duration_cast<std::chrono::seconds>(stop - start);
std::cout << duration.count() << std::endl;
}