ATLAS Offline Software
Loading...
Searching...
No Matches
madspin.py
Go to the documentation of this file.
1# Copyright (C) 2002-2022 CERN for the benefit of the ATLAS collaboration
2
3import glob
4import os
5import shutil
6import subprocess
7from AthenaCommon import Logging
8from ...decorators import timed
9from ...utility import LHE, ProcessManager, SingleProcessThread
10
11
12logger = Logging.logging.getLogger("PowhegControl")
13
14
15@timed("MadSpin post-processing")
16def MadSpin(process, powheg_LHE_output):
17 """! Move output to correctly named file.
18
19 @param process External MadSpin process.
20 @param powheg_LHE_output Name of LHE file produced by PowhegBox.
21
22 @author James Robinson <james.robinson@cern.ch>
23 """
24 process.expose()
25 __construct_inputs(powheg_LHE_output, process)
26 __run_executable(process.executable)
27 __prepare_outputs(powheg_LHE_output)
28
29
30def __construct_inputs(input_LHE_events, process):
31 """! Construct MadSpin runcard.
32
33 @param input_LHE_events Input LHE file name.
34 @param process MadSpin process.
35
36 @author James Robinson <james.robinson@cern.ch>
37 """
38 # Find insertion point for MadSpin header
39 logger.info("Constructing MadSpin runcard header")
40 pre_header = "\n".join([elem for elem in [LHE.opening_tag(input_LHE_events), LHE.comment_block(input_LHE_events)] if elem])
41 header_contents = LHE.header_block(input_LHE_events).replace("<header>", "").replace("</header>", "").strip("\n")
42 post_header = LHE.init_block(input_LHE_events)
43 postamble = LHE.postamble(input_LHE_events)
44
45 # Write events to LHE file with MadSpin information
46 with open("madspin_LHE_input.lhe", "w") as f_madspin_LHE:
47 f_madspin_LHE.write("{}\n".format(pre_header))
48 f_madspin_LHE.write("<header>\n")
49 f_madspin_LHE.write("<mgversion>\n")
50 f_madspin_LHE.write("{}\n".format(os.environ["MADPATH"].split("/")[-1].split("v")[-1].split("_p")[0].replace("_", ".")))
51 f_madspin_LHE.write("</mgversion>\n")
52 f_madspin_LHE.write("<mg5proccard>\n")
53 f_madspin_LHE.write("set group_subprocesses Auto\n")
54 f_madspin_LHE.write("set ignore_six_quark_processes False\n")
55 f_madspin_LHE.write("set loop_optimized_output True\n")
56 f_madspin_LHE.write("set gauge unitary\n")
57 f_madspin_LHE.write("set complex_mass_scheme False\n")
58 f_madspin_LHE.write("import model {}\n".format(process.MadSpin_model))
59 if int(process.MadSpin_nFlavours) == 4:
60 f_madspin_LHE.write("define p = g u c d s u~ c~ d~ s~\n")
61 f_madspin_LHE.write("define j = g u c d s u~ c~ d~ s~\n")
62 elif int(process.MadSpin_nFlavours) == 5:
63 f_madspin_LHE.write("define p = g u c d s b u~ c~ d~ s~ b~\n")
64 f_madspin_LHE.write("define j = g u c d s b u~ c~ d~ s~ b~\n")
65 else:
66 raise ValueError("'MadSpin_nFlavours' must be set to 4 or 5")
67 if process.MadSpin_taus_are_leptons:
68 f_madspin_LHE.write("define l+ = e+ mu+ ta+\n")
69 f_madspin_LHE.write("define l- = e- mu- ta-\n")
70 else:
71 f_madspin_LHE.write("define l+ = e+ mu+\n")
72 f_madspin_LHE.write("define l- = e- mu-\n")
73 f_madspin_LHE.write("{}\n".format(process.MadSpin_process))
74 f_madspin_LHE.write("output tchan\n")
75 f_madspin_LHE.write("</mg5proccard>\n")
76 f_madspin_LHE.write("<mgruncard>\n")
77 f_madspin_LHE.write("#0.01 = req_acc_FO ! needed for determining LO/NLO - see AGENE-1459\n")
78 f_madspin_LHE.write("{} = nevents\n".format(LHE.event_counter(input_LHE_events)))
79 f_madspin_LHE.write("1 = lpp1 ! beam 1 type (0 = no PDF)\n")
80 f_madspin_LHE.write("1 = lpp2 ! beam 2 type (0 = no PDF)\n")
81 f_madspin_LHE.write("{} = ebeam1 ! beam 1 energy in GeV\n".format(process.beam_energy))
82 f_madspin_LHE.write("{} = ebeam2 ! beam 2 energy in GeV\n".format(process.beam_energy))
83 f_madspin_LHE.write("{} = bwcutoff\n".format(process.bwcutoff))
84 f_madspin_LHE.write("</mgruncard>\n")
85 f_madspin_LHE.write("<slha>\n")
86 if process.MadSpin_paramcard != "":
87 with open(process.MadSpin_paramcard) as file:
88 for line in file:
89 f_madspin_LHE.write(line)
90 else:
91 f_madspin_LHE.write("######################################################################\n")
92 f_madspin_LHE.write("## PARAM_CARD AUTOMATICALY GENERATED BY MG5 ####\n")
93 f_madspin_LHE.write("######################################################################\n")
94 f_madspin_LHE.write("###################################\n")
95 f_madspin_LHE.write("## INFORMATION FOR LOOP\n")
96 f_madspin_LHE.write("###################################\n")
97 f_madspin_LHE.write("BLOCK LOOP #\n")
98 f_madspin_LHE.write(" 1 {:e} # mu_r\n".format(process.mass_Z))
99 f_madspin_LHE.write("###################################\n")
100 f_madspin_LHE.write("## INFORMATION FOR MASS\n")
101 f_madspin_LHE.write("###################################\n")
102 f_madspin_LHE.write("BLOCK MASS #\n")
103 f_madspin_LHE.write(" 1 0.000000e+00 # d : 0.0\n")
104 f_madspin_LHE.write(" 2 0.000000e+00 # u : 0.0\n")
105 f_madspin_LHE.write(" 3 0.000000e+00 # s : 0.0\n")
106 f_madspin_LHE.write(" 4 0.000000e+00 # c : 0.0\n")
107 f_madspin_LHE.write(" 5 {:e} # mb\n".format(process.mass_b))
108 f_madspin_LHE.write(" 6 {:e} # mt\n".format(process.mass_t))
109 f_madspin_LHE.write(" 11 0.000000e+00 # e- : 0.0\n")
110 f_madspin_LHE.write(" 12 0.000000e+00 # ve : 0.0\n")
111 f_madspin_LHE.write(" 13 0.000000e+00 # mu- : 0.0\n")
112 f_madspin_LHE.write(" 14 0.000000e+00 # vm : 0.0\n")
113 f_madspin_LHE.write(" 15 {:e} # mta\n".format(process.mass_tau))
114 f_madspin_LHE.write(" 16 0.000000e+00 # vt : 0.0\n")
115 f_madspin_LHE.write(" 21 0.000000e+00 # g : 0.0\n")
116 f_madspin_LHE.write(" 22 0.000000e+00 # a : 0.0\n")
117 f_madspin_LHE.write(" 23 {:e} # mz\n".format(process.mass_Z))
118 f_madspin_LHE.write(" 24 {:e} # w+\n".format(process.mass_W))
119 f_madspin_LHE.write(" 25 {:e} # mh\n".format(process.mass_H))
120 f_madspin_LHE.write(" 82 0.000000e+00 # gh : 0.0\n")
121 f_madspin_LHE.write("###################################\n")
122 f_madspin_LHE.write("## INFORMATION FOR SMINPUTS\n")
123 f_madspin_LHE.write("###################################\n")
124 f_madspin_LHE.write("BLOCK SMINPUTS #\n")
125 f_madspin_LHE.write(" 1 {:e} # aewm1\n".format(process.alphaem_inv))
126 f_madspin_LHE.write(" 2 {:e} # gf\n".format(process.G_F))
127 f_madspin_LHE.write(" 3 {:e} # as\n".format(process.alphaqcd))
128 f_madspin_LHE.write("###################################\n")
129 f_madspin_LHE.write("## INFORMATION FOR YUKAWA\n")
130 f_madspin_LHE.write("###################################\n")
131 f_madspin_LHE.write("BLOCK YUKAWA #\n")
132 f_madspin_LHE.write(" 5 {:e} # ymb\n".format(process.mass_b))
133 f_madspin_LHE.write(" 6 {:e} # ymt\n".format(process.mass_t))
134 f_madspin_LHE.write(" 15 {:e} # ymtau\n".format(process.mass_tau))
135 f_madspin_LHE.write("###################################\n")
136 f_madspin_LHE.write("## INFORMATION FOR QNUMBERS 82\n")
137 f_madspin_LHE.write("###################################\n")
138 f_madspin_LHE.write("BLOCK QNUMBERS 82 # gh\n")
139 f_madspin_LHE.write(" 1 0 # 3 times electric charge\n")
140 f_madspin_LHE.write(" 2 1 # number of spin states (2s+1)\n")
141 f_madspin_LHE.write(" 3 8 # colour rep (1: singlet, 3: triplet, 8: octet)\n")
142 f_madspin_LHE.write(" 4 1 # particle/antiparticle distinction (0=own anti)\n")
143 if "ckm" in process.MadSpin_model or "CKM" in process.MadSpin_model:
144 logger.info("Adding Wolfenstein parameters in param_card for model {}".format(process.MadSpin_model))
145 f_madspin_LHE.write("###################################\n")
146 f_madspin_LHE.write("INFORMATION FOR WOLFENSTEIN\n")
147 f_madspin_LHE.write("###################################\n")
148 f_madspin_LHE.write("Block wolfenstein\n")
149 f_madspin_LHE.write(" 1 {:e} # lamWS\n".format(process.Wolfenstein_lambda))
150 f_madspin_LHE.write(" 2 {:e} # AWS\n".format(process.Wolfenstein_A))
151 f_madspin_LHE.write(" 3 {:e} # rhoWS\n".format(process.Wolfenstein_rho))
152 f_madspin_LHE.write(" 4 {:e} # etaWS\n".format(process.Wolfenstein_eta))
153 f_madspin_LHE.write("#\n")
154 f_madspin_LHE.write("#*************************\n")
155 f_madspin_LHE.write("# Decay widths *\n")
156 f_madspin_LHE.write("#*************************\n")
157 f_madspin_LHE.write("#\n")
158 f_madspin_LHE.write("# PDG Width\n")
159 f_madspin_LHE.write("DECAY 1 0.000000e+00\n")
160 f_madspin_LHE.write("#\n")
161 f_madspin_LHE.write("# PDG Width\n")
162 f_madspin_LHE.write("DECAY 2 0.000000e+00\n")
163 f_madspin_LHE.write("#\n")
164 f_madspin_LHE.write("# PDG Width\n")
165 f_madspin_LHE.write("DECAY 3 0.000000e+00\n")
166 f_madspin_LHE.write("#\n")
167 f_madspin_LHE.write("# PDG Width\n")
168 f_madspin_LHE.write("DECAY 4 0.000000e+00\n")
169 f_madspin_LHE.write("#\n")
170 f_madspin_LHE.write("# PDG Width\n")
171 f_madspin_LHE.write("DECAY 5 0.000000e+00\n")
172 f_madspin_LHE.write("#\n")
173 if "ckm" in process.MadSpin_model or "CKM" in process.MadSpin_model:
174 logger.info("Top decay modes not set in param_card for model {} - BRs will be calculated at LO by MadSpin".format(process.MadSpin_model))
175 f_madspin_LHE.write("# PDG Width\n")
176 f_madspin_LHE.write("DECAY 6 {:e}\n".format(process.width_t))
177 f_madspin_LHE.write("#\n")
178 else:
179 # Scale t BRs down so that sum is 1.0
180 BR_t_to_any = process.BR_t_to_Wb + process.BR_t_to_Ws + process.BR_t_to_Wd
181 f_madspin_LHE.write("# PDG Width\n")
182 f_madspin_LHE.write("DECAY 6 {:e}\n".format(process.width_t))
183 f_madspin_LHE.write("# BR NDA ID1 ID2 ...\n")
184 f_madspin_LHE.write(" {:e} 2 5 24 # 1.32\n".format(process.BR_t_to_Wb / BR_t_to_any))
185 f_madspin_LHE.write(" {:e} 2 3 24 # 1.32\n".format(process.BR_t_to_Ws / BR_t_to_any))
186 f_madspin_LHE.write(" {:e} 2 1 24 # 1.32\n".format(process.BR_t_to_Wd / BR_t_to_any))
187 f_madspin_LHE.write("#\n")
188 f_madspin_LHE.write("# PDG Width\n")
189 f_madspin_LHE.write("DECAY 11 0.000000e+00\n")
190 f_madspin_LHE.write("#\n")
191 f_madspin_LHE.write("# PDG Width\n")
192 f_madspin_LHE.write("DECAY 12 0.000000e+00\n")
193 f_madspin_LHE.write("#\n")
194 f_madspin_LHE.write("# PDG Width\n")
195 f_madspin_LHE.write("DECAY 13 0.000000e+00\n")
196 f_madspin_LHE.write("#\n")
197 f_madspin_LHE.write("# PDG Width\n")
198 f_madspin_LHE.write("DECAY 14 0.000000e+00\n")
199 f_madspin_LHE.write("#\n")
200 f_madspin_LHE.write("# PDG Width\n")
201 f_madspin_LHE.write("DECAY 15 0.000000e+00\n")
202 f_madspin_LHE.write("#\n")
203 f_madspin_LHE.write("# PDG Width\n")
204 f_madspin_LHE.write("DECAY 16 0.000000e+00\n")
205 f_madspin_LHE.write("#\n")
206 f_madspin_LHE.write("# PDG Width\n")
207 f_madspin_LHE.write("DECAY 21 0.000000e+00\n")
208 f_madspin_LHE.write("#\n")
209 f_madspin_LHE.write("# PDG Width\n")
210 f_madspin_LHE.write("DECAY 22 0.000000e+00\n")
211 f_madspin_LHE.write("#\n")
212 f_madspin_LHE.write("# PDG Width\n")
213 f_madspin_LHE.write("DECAY 23 {:e}\n".format(process.width_Z))
214 f_madspin_LHE.write("#\n")
215 if "ckm" in process.MadSpin_model or "CKM" in process.MadSpin_model:
216 logger.info("W decay modes not set in param_card for model {} - BRs will be calculated at LO by MadSpin".format(process.MadSpin_model))
217 f_madspin_LHE.write("# PDG Width\n")
218 f_madspin_LHE.write("DECAY 24 {}\n".format(process.width_W))
219 f_madspin_LHE.write("#\n")
220 else:
221 # Scale W BRs down so that sum is 1.0
222 BR_W_to_any = process.BR_W_to_hadrons + process.BR_W_to_leptons
223 f_madspin_LHE.write("# PDG Width\n")
224 f_madspin_LHE.write("DECAY 24 {}\n".format(process.width_W))
225 f_madspin_LHE.write("# BR NDA ID1 ID2 ...\n")
226 f_madspin_LHE.write(" {:e} 2 -1 2\n".format((process.BR_W_to_hadrons / 2.) / BR_W_to_any))
227 f_madspin_LHE.write(" {:e} 2 -3 4\n".format((process.BR_W_to_hadrons / 2.) / BR_W_to_any))
228 f_madspin_LHE.write(" {:e} 2 -11 12\n".format((process.BR_W_to_leptons / 3.) / BR_W_to_any))
229 f_madspin_LHE.write(" {:e} 2 -13 14\n".format((process.BR_W_to_leptons / 3.) / BR_W_to_any))
230 f_madspin_LHE.write(" {:e} 2 -15 16\n".format((process.BR_W_to_leptons / 3.) / BR_W_to_any))
231 f_madspin_LHE.write("#\n")
232 f_madspin_LHE.write("# PDG Width\n")
233 f_madspin_LHE.write("DECAY 25 {:e}\n".format(process.width_H))
234 f_madspin_LHE.write("#\n")
235 f_madspin_LHE.write("# PDG Width\n")
236 f_madspin_LHE.write("DECAY 82 0.000000e+00\n")
237 f_madspin_LHE.write("</slha>\n")
238 f_madspin_LHE.write("{}\n".format(header_contents))
239 f_madspin_LHE.write("</header>\n")
240 f_madspin_LHE.write("{}\n".format(post_header))
241 for event in LHE.event_iterator(input_LHE_events):
242 f_madspin_LHE.write(event)
243 f_madspin_LHE.write(postamble)
244
245 # Rename LHE files
246 shutil.move(input_LHE_events, "{}.undecayed".format(input_LHE_events))
247 shutil.move("madspin_LHE_input.lhe", input_LHE_events)
248
249 # Write MadSpin runcard
250 with open("madspin_runcard.txt", "w") as f_madspin_runcard:
251 f_madspin_runcard.write("import {}\n".format(input_LHE_events))
252 f_madspin_runcard.write("set spinmode {}\n".format(process.MadSpin_mode))
253 if process.MadSpin_max_weight_ps_point>0:
254 f_madspin_runcard.write("set max_weight_ps_point {}\n".format(int(process.MadSpin_max_weight_ps_point)))
255 if process.MadSpin_Nevents_for_max_weight>0:
256 f_madspin_runcard.write("set Nevents_for_max_weight {}\n".format(int(process.MadSpin_Nevents_for_max_weight)))
257 for decay in process.MadSpin_decays:
258 f_madspin_runcard.write("{0}\n".format(decay))
259 f_madspin_runcard.write("launch\n")
260 f_madspin_runcard.write("quit\n")
261
262
263@timed("MadSpin executable")
264def __run_executable(executable):
265 """! Run MadSpin executable.
266
267 @author James Robinson <james.robinson@cern.ch>
268 """
269 if not os.path.isfile(executable):
270 raise OSError("MadSpin executable {} not found!".format(executable))
271 logger.info("MadSpin executable: {}".format(executable))
272 with open("madspin_runcard.txt", "r") as runcard_input:
273 processes = [SingleProcessThread(['python', executable], stdin=runcard_input, ignore_output=["INFO:","MadSpin>"],
274 error_output=["Command \"launch\" interrupted with error:","MadSpinError"],
275 info_output=["generating the production square matrix element"],
276 warning_output=["stty: standard input: Invalid argument",
277 "stty: standard input: Inappropriate ioctl for device",
278 "no version information available",
279 "CRITICAL: Branching ratio larger than one for"])]
280 manager = ProcessManager(processes)
281 while manager.monitor():
282 pass
283
284
285def __prepare_outputs(input_LHE_events):
286 """! Prepare MadSpin output.
287
288 @author James Robinson <james.robinson@cern.ch>
289 """
290 logger.info("Preparing MadSpin output")
291
292 # Unzip MadSpin events
293 subprocess.call("gunzip pwgevents_decayed.lhe.gz 2> /dev/null", shell=True)
294 shutil.move("pwgevents_decayed.lhe", "{}.decayed".format(input_LHE_events))
295
296 # Get appropriate header sections from Powheg and MadSpin LHE
297 powheg_LHE, madspin_LHE = "{}.undecayed".format(input_LHE_events), "{}.decayed".format(input_LHE_events)
298 pre_header = "\n".join([elem for elem in [LHE.opening_tag(powheg_LHE), LHE.comment_block(powheg_LHE)] if elem])
299 header = LHE.header_block(madspin_LHE)
300 post_header = LHE.init_block(powheg_LHE)
301 postamble = LHE.postamble(powheg_LHE)
302
303 # Write events to LHE file with MadSpin information
304 with open("pwgevents_with_MadSpin.lhe", "w") as f_combined_LHE:
305 f_combined_LHE.write("{}\n".format(pre_header))
306 f_combined_LHE.write("{}\n".format(header))
307 f_combined_LHE.write("{}\n".format(post_header))
308 for event in LHE.event_iterator("{}.decayed".format(input_LHE_events)):
309 f_combined_LHE.write(event)
310 f_combined_LHE.write(postamble)
311
312 # Rename output file
313 if os.path.isfile(input_LHE_events):
314 shutil.move(input_LHE_events, "{}.undecayed_ready_for_madspin".format(input_LHE_events))
315 shutil.move("pwgevents_with_MadSpin.lhe", input_LHE_events)
316
317 for LHE_tarball in list(set(glob.glob("*lhe*.gz") + glob.glob("*events*.gz"))):
318 logger.info("Cleaning up MadSpin tarball: {}".format(LHE_tarball))
319 os.remove(LHE_tarball)
python.utility.process_handling.ProcessManager
Wrapper to handle multiple Powheg subprocesses.
Definition process_handling.py:11
python.utility.process_handling.SingleProcessThread
Single executable running in a subprocess (usually PowhegBox).
Definition process_handling.py:40
set
STL class.
replace
std::string replace(std::string s, const std::string &s2, const std::string &s3)
Definition hcg.cxx:310
split
std::vector< std::string > split(const std::string &s, const std::string &t=":")
Definition hcg.cxx:177
python.algorithms.postprocessors.madspin.__prepare_outputs
__prepare_outputs(input_LHE_events)
Prepare MadSpin output.
Definition madspin.py:285
python.algorithms.postprocessors.madspin.__construct_inputs
__construct_inputs(input_LHE_events, process)
Construct MadSpin runcard.
Definition madspin.py:30
python.algorithms.postprocessors.madspin.__run_executable
__run_executable(executable)
Run MadSpin executable.
Definition madspin.py:264
python.algorithms.postprocessors.madspin.MadSpin
MadSpin(process, powheg_LHE_output)
Move output to correctly named file.
Definition madspin.py:16
Generated on for ATLAS Offline Software by doxygen 1.14.0