-
Notifications
You must be signed in to change notification settings - Fork 0
/
calculate_MG_xsection.py
260 lines (218 loc) · 10.6 KB
/
calculate_MG_xsection.py
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
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
# this script creates a MadGraph configuration file, runs Madgraph and outputs the calculated quantities for a given process
import sys, os
import json
import copy
import subprocess as sp
import string
import numpy as np
import math
from optparse import OptionParser
from ParameterSpace import ParameterSpace
from time import time
class Timer():
def __init__(self):
self.start_time = time()
def elapsed_time(self):
return time() - self.start_time
def restart(self):
self.start_time = time()
madgraph_executable = './bin/mg5_aMC'
def make_folder_if_not_exists(folder):
if not os.path.exists(folder):
try:
os.makedirs(folder)
except:
print "Could not create a folder ", folder
def write_data_to_JSON(data, JSON_output_file, indent = True):
output_file = open(JSON_output_file, 'w')
if indent:
output_file.write(json.dumps(data, indent=4, sort_keys = True))
else:
output_file.write(json.dumps(data, sort_keys = True))
output_file.close()
def create_MG_config(subprocess, parameters, file_name_suffix, auto_width = True):
global workspace_path
name = workspace_path + subprocess + '_' + file_name_suffix
filename = name + '.dat'
file = open(filename,'w')
file.write('#************************************************************\n')
file.write('#* MadGraph 5 *\n')
file.write('#* *\n')
file.write('#* * * *\n')
file.write('#* * * * * *\n')
file.write('#* * * * * 5 * * * * *\n')
file.write('#* * * * * *\n')
file.write('#* * * *\n')
file.write('#* *\n')
file.write('#* *\n')
file.write('#* The MadGraph Development Team - Please visit us at *\n')
file.write('#* https://server06.fynu.ucl.ac.be/projects/madgraph *\n')
file.write('#* *\n')
file.write('#************************************************************\n')
file.write('#* *\n')
file.write('#* This is an auto-generated file for MadGraph 5 *\n')
file.write('#* *\n')
file.write('#* run as ./bin/mg5 filename *\n')
file.write('#* *\n')
file.write('#************************************************************\n')
file.write('set automatic_html_opening False\n')
file.write('import model models/MonotopDMF_UFO -modelname\n')
file.write('define j = g u c d s b u~ c~ d~ s~ b~\n')
file.write('define l+ = e+ mu+ ta+\n')
file.write('define l- = e- mu- ta-\n')
file.write('define vl = ve vm vt\n')
file.write('define vl~ = ve~ vm~ vt~\n')
file.write('\n')
if 'tt_excl' in subprocess:
file.write('# tt exclusive decay\n')
file.write('generate p p > t t, (t > b W+, W+ > l+ vl)\n')
file.write('add process p p > t~ t~, (t~ > b~ W-, W- > l- vl~)\n')
file.write('\n')
elif 'onshell' in subprocess:
file.write('# Visible on-shell V decay\n')
file.write('generate p p > V > t t u~, (t > b W+, W+ > l+ vl)\n')
file.write('add process p p > V > t~ t~ u, (t~ > b~ W-, W- > l- vl~)\n')
file.write('\n')
elif 'offshell' in subprocess:
file.write('# Off-shell V decay\n')
file.write('generate p p > t t u~ $$ V, (t > b W+, W+ > l+ vl)\n')
file.write('add process p p > t~ t~ u $$ V, (t~ > b~ W-, W- > l- vl~)\n')
file.write('\n')
elif 'Monotop' in subprocess or 'monotop' in subprocess:
file.write('# Monotop\n')
file.write('generate p p > t psi psibar, (t > b W+, W+ > l+ vl)\n')
file.write('add process p p > t~ psi psibar, (t~ > b~ W-, W- > l- vl~)\n')
file.write('\n')
file.write('# Output processes to MadEvent directory\n')
file.write('output %s\n' % name)
file.write('\n')
file.write('launch %s\n' % name)
file.write('set Mpsi %e # changing the psi mass\n' % parameters.mDM)
file.write('set ar %e # changing the a_r coupling constant\n' % parameters.a_r)
file.write('set gg %e # changing the gg coupling constant\n' % parameters.g)
file.write('set MV %e # changing the V mass\n' % parameters.mV)
if auto_width:
file.write('set WV Auto # changing the V width\n')
else:
file.write('set WV %e # changing the V width\n' % parameters.G_tot)
file.write('launch %s -i\n' % name)
file.write('print_results --path=%s.txt --format=short\n' % name)
file.close
def run_MG_config(subprocess, parameters, file_name_suffix, auto_width = True):
global workspace_path
name = workspace_path + subprocess + '_' + file_name_suffix
config_filename = name + '.dat'
output_filename = name + '.out'
output_file = open(output_filename,'w')
print 'Executing command: ', madgraph_executable + ' ./' + config_filename
print 'Output file: ', output_filename
run = sp.Popen([madgraph_executable, config_filename], cwd='./', stdout=output_file)
run.communicate()
output_file.close
def read_MG_output(subprocess, parameters, file_name_suffix, auto_width = True, overwrite_with_MG = False):
global workspace_path, output_path
name = workspace_path + subprocess + '_' + file_name_suffix
output_banner_filename = name + '/Events/run_01/run_01_tag_1_banner.txt'
output_banner = open(output_banner_filename,'r')
# new parameters instance for MG-calculated quantities
parameters_MG = ParameterSpace()
# parse the output for numbers
for line in output_banner.readlines():
columns = string.split(line)
if '# ar' in line:
parameters_MG.set_a_r(columns[1])
if '# gg' in line:
parameters_MG.set_g(columns[1])
if '# mv' in line:
parameters_MG.set_mV(columns[1])
if '# mpsi' in line:
parameters_MG.set_mDM(columns[1])
if 'DECAY 32' in line:
parameters_MG.set_G_tot(columns[2])
if '1000023 -1000023' in line:
parameters_MG.set_BR(columns[0])
if 'Integrated weight (pb)' in line:
cross_section = float(columns[-1])
output_banner.close
print '*'*100
print 'MadGraph calculated quantities:'
parameters_MG.print_initial_parameters()
print '='*100
print 'Calculated cross-section [pb] = ', cross_section
if overwrite_with_MG:
print 'Using MG calculated parameters in the output.'
output_parameters = parameters_MG
else:
print 'Using analytically calculated parameters in the output.'
output_parameters = parameters
data_to_write = copy.copy(output_parameters.__dict__)
data_to_write['xsection'] = cross_section
data_to_write['process'] = subprocess
JSON_file_name = output_path + subprocess + '_' + file_name_suffix + '.txt'
write_data_to_JSON(data_to_write, JSON_file_name)
print 'Data written into JSON file ', JSON_file_name
return parameters_MG
if __name__ == '__main__':
timer = Timer()
timer_value = timer.elapsed_time()
parser = OptionParser()
parser.add_option( "-o", "--output_path", dest = "output_path", default = 'output_JSON',
help = "Set the output path where all the JSON files with results are stored (default: output_JSON)" )
parser.add_option( "-w", "--workspace_path", dest = "workspace_path", default = 'workspace',
help = "Set the workspace path where all the MG files are stored (default: workspace)" )
parser.add_option( "-M", "--mV", dest = "mV", default = 2000,
help = "Set the mediator mass (mV), default: 2000 GeV" )
parser.add_option( "-m", "--mDM", dest = "mDM", default = 1,
help = "Set the DM candidate mass (mDM), default: 1 GeV" )
parser.add_option( "-a", "--a_r", dest = "a_r", default = 0.5,
help = "Set the a_r coupling constant (V_tu vertex)" )
parser.add_option( "-g", "--g", dest = "g", default = '',
help = "Set the g coupling constant (V_chi_chi vertex)" )
parser.add_option( "-G", "--total_width", dest = "total_width", default = '',
help = "Set the total width in GeV (default - automatic)" )
parser.add_option( "-B", "--BR", dest = "BR", default = '',
help = "Set the invisible BR (default - automatic)" )
parser.add_option( "-s", "--signal", dest = "signal", default = 'tt_exclusive',
help = "Choose the signal of interest. Default: tt_exclusive (prompt tt production), alternatively ttu_offshellV, ttu_onshellV, monotop" )
( options, args ) = parser.parse_args()
workspace_path = options.workspace_path + '/'
make_folder_if_not_exists(workspace_path)
output_path = options.output_path + '/'
make_folder_if_not_exists(output_path)
current_process = options.signal
if not options.g and not options.total_width and not options.BR:
print 'Insufficient input parameters, please enter the total width (G), branching ratio (BR) or the g coupling constant (g).'
sys.exit(1)
parameters = ParameterSpace()
parameters.set_mV(options.mV)
parameters.set_mDM(options.mDM)
parameters.set_a_r(options.a_r)
auto_width = True
if options.total_width:
auto_width = False
if options.g:
parameters.set_g(options.g)
if options.BR:
parameters.set_BR(options.BR)
if options.total_width:
parameters.set_G_tot(options.total_width)
file_name_suffix = copy.copy(parameters.parameter_space_name())
print 'Input parameters:'
parameters.print_initial_parameters()
parameters.calculate_all()
print '*'*100
print 'Analytically calculated parameters:'
parameters.print_calculated_parameters()
print '.'*100
print 'Creating the new MadGraph config for %s process.' % current_process
create_MG_config(current_process, parameters, file_name_suffix, auto_width)
run_MG_config(current_process, parameters, file_name_suffix, auto_width)
print 'Runtime for process %s, parameter set %s : %.1f min' % (current_process, file_name_suffix, (timer.elapsed_time()-timer_value)/60)
timer_value = timer.elapsed_time()
parameters_copy = copy.copy(parameters)
MG_calculated_parameters = read_MG_output(current_process, parameters, file_name_suffix, auto_width)
if parameters_copy!=MG_calculated_parameters:
print 'Warning: analytically calculated parameters do not exactly match with MadGraph calculated ones.'
print 'Analytical parameters: ', parameters_copy.__dict__
print 'MG-calculated parameters: ', MG_calculated_parameters.__dict__
sys.exit(1)