-
Notifications
You must be signed in to change notification settings - Fork 0
/
main_for_HH_decay.cc
205 lines (170 loc) · 6.55 KB
/
main_for_HH_decay.cc
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
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
#include <fstream>
#include "Pythia.h"
using namespace Pythia8;
int main() {
//unsigned const nfiles=4;
string path= "/home/xanda/Documents/";
/*
for(unsigned ifile=0; ifile<11; ifile++){
//unsigned ifile=0;
string namefile_in=path + "MR_";
if(ifile==0) namefile_in += "test-MR610.lhe";
if(ifile==1) namefile_in += "500_0.lhe";
if(ifile==2) namefile_in += "600_0.lhe";
if(ifile==3) namefile_in += "700_0.lhe";
if(ifile==4) namefile_in += "800_0.lhe";
if(ifile==5) namefile_in += "900_0.lhe";
if(ifile==6) namefile_in += "1000_0.lhe";
if(ifile==7) namefile_in += "1500_0.lhe";
if(ifile==8) namefile_in += "2000_0.lhe";
if(ifile==9) namefile_in += "2500_0.lhe";
if(ifile==10) namefile_in += "3000_0.lhe";
string namefile_out=namefile_in + ".pythia";
*/
string namefile_in=path + "test-MR610.lhe";
string namefile_out=namefile_in + ".pythia";
//namefile_in += "test-MR610.lhe";
cout<<"\n namefile_in = "<<namefile_in<<endl;
cout<<"\n namefile_out = "<<namefile_out<<endl;
// output file
// we want to store the list of all final state particles
ofstream out_pythia;
// Highest precision required for jet clustering
out_pythia.precision(15);
// Generator. We here stick with default values, but changes
// could be inserted with readString or readFile.
Pythia pythia;
// Initialize Les Houches Event File run. List initialization information.
pythia.readString("Beams:frameType = 4");
// the analysis program
string sfile = "Beams:LHEF ="+namefile_in;
pythia.readString(sfile.c_str());
out_pythia.open(namefile_out.c_str());
// Allow for possibility of a few faulty events.
int nAbort = 10;
int iAbort = 0;
// pythia.settings.listAll();
// Settings
/* turn of hadronization settings - for testing
pythia.readString("PartonLevel:MI = off"); // Off multiple interactions
pythia.readString("PartonLevel:ISR = off"); // Shower on
pythia.readString("PartonLevel:FSR = off"); // Shower on
pythia.readString("HadronLevel:all = off"); // Of hadronization
pythia.readString("25:onMode = off");
// pythia.readString("25:onIfMatch = 5 -5"); // b
// pythia.readString("25:onIfMatch = 4 -4"); // c
// pythia.readString("25:onIfMatch = 21 21"); // gluon
pythia.readString("25:onIfMatch = 24 -24"); // W
pythia.readString("24:onMode = off");
pythia.readString("-24:onMode = off");
pythia.readString("24:onIfMatch = 12 11"); // e ve
pythia.readString("24:onIfMatch = 14 13"); // mu numu
pythia.readString("-24:onIfMatch = 12 -11"); // e ve
pythia.readString("-24:onIfMatch = 14 -13"); // mu numu
pythia.readString("21:mayDecay = no");
pythia.readString("4:mayDecay = no");
pythia.readString("-4:mayDecay = no");
pythia.readString("5:mayDecay = no");
pythia.readString("-5:mayDecay = no");
pythia.readString("5:mayDecay = no");
pythia.readString("-5:mayDecay = no");
pythia.readString("11:mayDecay = no");
pythia.readString("-11:mayDecay = no");
pythia.readString("13:mayDecay = no");
pythia.readString("-13:mayDecay = no");
*/
// turn decays
pythia.readString("HadronLevel:all = on"); // generic hadronization settings
// hadronic Higgs decay
pythia.readString("25:onMode = off");
pythia.readString("25:onIfMatch = 5 -5"); // b
pythia.readString("25:onIfMatch = 4 -4"); // c
pythia.readString("25:onIfMatch = 21 21"); // gluon
// WW "higgs" decay
pythia.readString("35:onMode = off");
pythia.readString("35:onIfMatch = 24 -24"); // W
// W decay
pythia.readString("24:onMode = off");
pythia.readString("24:onIfMatch = 12 11"); // e ve
pythia.readString("24:onIfMatch = 14 13"); // mu numu
pythia.readString("-24:onMode = off");
pythia.readString("-24:onIfMatch = 12 -11"); // e ve
pythia.readString("-24:onIfMatch = 14 -13"); // mu numu
// Initialization
pythia.init();
// Begin event loop; generate until none left in input file.
for (int iEvent = 0; ; ++iEvent) {
cout<<"\n ievent = "<<iEvent<<"\n"<<endl;
// Generate events, and check whether generation failed.
if (!pythia.next()) {
// If failure because reached end of file then exit event loop.
if (pythia.info.atEndOfFile()) break;
// First few failures write off as "acceptable" errors, then quit.
if (++iAbort < nAbort) continue;
break;
}
// Acess event record
cout<<"Number of particles = "<<pythia.event.size()<<endl;
vector<int> pID;
vector<double> px;
vector<double> py;
vector<double> pz;
vector<double> E;
vector<int> mother;
vector<int> code;
// Some checks on the event record
// Check for example that at least we have two bs and two bbars
for (int i = 0; i < pythia.event.size(); i++){
int particle_id = pythia.event[i].id();
int particle_status = pythia.event[i].status();
int particle_mother = pythia.event[i].mother1();
// save only final state particles
if(particle_status>0){
cout<<i<<" "<<particle_id<<" "<<particle_mother<<endl;
double ppx= pythia.event[i].px();
double ppy= pythia.event[i].py();
double ppz= pythia.event[i].pz();
double EE= pythia.event[i].e();
//cout<<px<<" "<<py<<" "<<pz<<" "<<E<<endl;
pID.push_back(particle_id);
px.push_back(ppx);
py.push_back(ppy);
pz.push_back(ppz);
E.push_back(EE);
mother.push_back(particle_mother);
code.push_back(particle_id);
}
}
// Save into file
out_pythia<<"#"<<endl;
cout<<"Number of final state particles = "<<E.size()<<"\n"<<endl;
out_pythia<<E.size()<<endl;
/*
for(unsigned i=0;i<E.size();i++){
// First of all write the b quarks
if(code.at(i)==5 && code.at(i+1)==-5 && mother.at(i)==mother.at(i+1) ){
out_pythia<<px.at(i)<<" "<<py.at(i)<<" "<<pz.at(i)<<" "<<E.at(i)<<" "<<endl;
out_pythia<<px.at(i+1)<<" "<<py.at(i+1)<<" "<<pz.at(i+1)<<" "<<E.at(i+1)<<" "<<endl;
}
}
// Now all the other particles
for(unsigned i=0;i<E.size();i++){
// First of all write the b quarks
if(fabs(code.at(i))<5 || code.at(i)==21){
out_pythia<<px.at(i)<<" "<<py.at(i)<<" "<<pz.at(i)<<" "<<E.at(i)<<" "<<endl;
}
}
*/
for(unsigned i=0;i<E.size();i++){
out_pythia<<pID.at(i)<<" "<<px.at(i)<<" "<<py.at(i)<<" "<<pz.at(i)<<" "<<E.at(i)<<" "<<endl;
}
// End of event loop.
}
out_pythia.close();
// Give statistics. Print histogram.
pythia.stat();
//delete pythia;
// } //for unsigned
// Done.
return 0;
}