ATLAS Offline Software
Loading...
Searching...
No Matches
SuperChicUtils.py
Go to the documentation of this file.
1# Copyright (C) 2002-2023 CERN for the benefit of the ATLAS collaboration
2
3
4import subprocess, os, shlex, re, shutil
5from packaging.version import Version
6
7from AthenaCommon import Logging
8
9
10logger = Logging.logging.getLogger("Superchic_i")
11
12
13class SuperChicConfig:
14
15 def __init__(self, runArgs):
16 self.superchicpath = os.environ['SUPERCHIC_SOURCE_PATH']
17
18 #SuperChic specific variables for input.DAT, see writeInputDAT function for more elaboration
19 self.rts = 13000. #collision energy (GeV)
20 if hasattr(runArgs,"ecmEnergy"):
21 self.rts = runArgs.ecmEnergy
22
23 self.isurv = 4 #Model of soft survival (from 1 -> 4, corresponding to arXiv:1306.2149)
24 self.intag = 'in' #tag for input files
25 self.PDFname = 'MMHT2014lo68cl'
26 self.PDFmember = 0
27 self.proc = 5 #Please consult Superchic Manual https://superchic.hepforge.org/
28 self.beam = 'prot'
29 self.outtg = 'out'
30 self.sfaci = True
31 self.diff = 'el' #interaction: elastic ('el'), single ('sd','sda','sdb') and double ('dd') dissociation.
32 self.an = 208
33 self.az = 82
34 self.rz = 6.68
35 self.dz = 0.447
36 self.rn = 6.7
37 self.dn = 0.55
38 self.ionqcd = 'coh'
39 self.ionbreakup = True
40 self.fAA = '00'
41 self.fracsigX = 1
42 if (Version(os.environ['SUPERCHICVER']) >= Version("5.7")):
43 self.veto = False
44 self.veto_ecent = 53.9655
45 self.veto_eden = 263.313
46 self.rshad = 0.4
47 if (Version(os.environ['SUPERCHICVER']) >= Version("5.5.1")):
48 self.wrho = False
49 self.yrho = 2.5
50 self.accrho = 1.0
51 if (Version(os.environ['SUPERCHICVER']) >= Version("5.6.1")):
52 self.ion_incoh = False
53 self.ion_incoh_type = 'inel'
54 self.ncall = 10000
55 self.itmx = 10
56 self.prec = 0.5
57 self.ncall1 = 10000
58 self.inccall = 10000
59 self.itend = 1000
60
61 self.iseed = 34
62 if hasattr(runArgs,"randomSeed"):
63 self.iseed = runArgs.randomSeed
64
65 self.genunw = True
66
67 self.nev = "500"
68 if hasattr(runArgs,"maxEvents"):
69 self.nev = runArgs.maxEvents
70
71 self.erec = 'lhe' #output file type
72 self.readwt = False
73 self.wtmax = 0
74 self.ymin = -2.4 # Minimum object rapidity Y_X
75 self.ymax = 2.4 # Maximum object rapidity Y_X
76 self.mmin = 6. # Minimum object mass M_X
77 self.mmax = 500. # Maximum object mass M_X
78 self.gencuts = True # Generate cuts below
79 self.scorr = True
80 self.fwidth = True
81 self.ptxmax = 100.
82 self.ptamin = 3.0 # Minimum pT of outgoing object a (gamma)
83 self.ptbmin = 3.0 # Minimum pT of outgoing object b (gamma)
84 self.etaamin = -2.4 # Minimum eta of outgoing object a
85 self.etaamax = 2.4 # Maximum eta of outgoing object a
86 self.etabmin = -2.4 # Minimum eta of outgoing object b
87 self.etabmax = 2.4 # Maximum eta of outgoing object b
88 self.acoabmax = 100.
89 self.ptcmin = 0
90 self.etacmin = -2.5
91 self.etacmax = 2.5
92 self.ptdmin = 0
93 self.etadmin = -2.5
94 self.etadmax = 2.5
95 self.ptemin = 0
96 self.etaemin = -2.5
97 self.etaemax = 2.5
98 self.ptfmin = 0
99 self.etafmin = -2.5
100 self.etafmax = 2.5
101 self.rjet = 0.6
102 self.jalg = 'antikt'
103 self.m2b = 0.133
104 self.pdgid1 = 211
105 self.pdgid2 = -211
106 self.malp = 1000
107 self.gax = 0.001
108 self.alpt = 'ps'
109 self.mpol = 500
110 self.mmon = 933
111 self.gamm = 10
112 self.mcharg = 100
113 self.mneut = 80
114 self.wlp = 'el'
115 self.wlm = 'el'
116 self.tau = 0.04
117 self.mxs = 100
118 self.atau = 0.0
119 self.dtau = 0.0
120 self.calc_tau_coeff = False
121 self.tau_mom = 'atau'
122 self.tau_coeff = 1
123
124 def toFortran(self):
125
126 def fortDouble(x):
127 return str(x)+"d0"
128 def fortInt(x):
129 return str(x)
130 def fortBool(x):
131 return '.true.' if x else '.false.'
132 def fortStr(x):
133 return "'{}'".format(x)
134
135 conf = ""
136 conf+="***********************************************************************************\n"
137 conf+="****** Initialize afain if FIRST FIVE PARAMETERS ARE CHANGED (and beam = 'prot'):******\n"
138 conf+="***********************************************************************************\n"
139 conf+=fortDouble(self.rts) + " ! [rts] : CMS collision energy (GeV) \n"
140 conf+=fortInt(self.isurv) + " ! [isurv] : Model of soft survival (from 1 -> 4)\n"
141 conf+=fortStr(self.intag) + " ! [intag] for input files \n"
142 conf+="***********************************************************************************\n"
143 conf+="***********************************************************************************\n"
144 conf+=fortStr(self.PDFname) + " ! [PDFname] : PDF set \n"
145 conf+=fortInt(self.PDFmember) + " ! [PDFmember] : PDF member \n"
146 conf+="***********************************************************************************\n"
147 conf+=fortInt(self.proc) + " ! [proc] : Process number \n"
148 conf+=fortStr(self.beam) + " ! [beam] : Beam type ('prot', 'ion', 'ionp' or 'el') \n"
149 conf+=fortStr(self.outtg) + " ! [outtg] : for output file \n"
150 conf+=fortBool(self.sfaci) + " ! [sfaci] : Include soft survival effects \n"
151 conf+="***********************************************************************************\n"
152 conf+="***********************************************************************************\n"
153 conf+="***********************************************************************************\n"
154 conf+=fortStr(self.diff) + " ! [diff] : dissociation flag \n"
155 conf+="***********************************************************************************\n"
156 conf+=fortDouble(self.an) + " ! [an] : Ion mass number \n"
157 conf+=fortDouble(self.az) + " ! [az] : Ion atomic number \n"
158 conf+=fortDouble(self.rz) + " ! [rz] : Ion proton density - radius \n"
159 conf+=fortDouble(self.dz) + " ! [dz] : Ion proton density - skin thickness \n"
160 conf+=fortDouble(self.rn) + " ! [rn] : Ion neutron density - radius \n"
161 conf+=fortDouble(self.dn) + " ! [dn] : Ion neutron density - skin thickness \n"
162 conf+=fortStr(self.ionqcd) + " ! [ionqcd] : Coherent ('coh') or incoherent ('incoh') for QCD-induced processes \n"
163 conf+=fortBool(self.ionbreakup) + " ! [ionbreakup] \n"
164 conf+=fortStr(self.fAA) + " ! [fAA] \n"
165 conf+=fortDouble(self.fracsigX) + " ! [fracsigX] : multiply sig_(gamA) by this factor (1d0 - default) \n"
166 if (Version(os.environ['SUPERCHICVER']) >= Version("5.7")):
167 conf+=fortBool(self.veto) + " ! [veto] : impose exclusivity veto in ion-ion case due to EMD \n"
168 conf+=fortDouble(self.veto_ecent) + " ! [veto_ecent] \n"
169 conf+=fortDouble(self.veto_eden) + " ! [veto_eden] \n"
170 conf+=fortDouble(self.rshad) + " ! [rshad] \n"
171 if (Version(os.environ['SUPERCHICVER']) >= Version("5.5.1")):
172 conf+=fortBool(self.wrho) + " ! [wrho] : generate concident rho production in AA collisions \n"
173 conf+=fortDouble(self.yrho) + " ! [yrho] : maximum rho |y| \n"
174 conf+=fortDouble(self.accrho) + " ! [accrho] : additional correction factor to account for rho->pipi acceptance, non-resonant dipions and different beams \n"
175 if (Version(os.environ['SUPERCHICVER']) >= Version("5.6.1")):
176 conf+=fortBool(self.ion_incoh) + " ! [ion_incoh] \n"
177 conf+=fortStr(self.ion_incoh_type) + " ! [ion_incoh_type] \n"
178 conf+="***********************************************************************************\n"
179 conf+="*************Integration parameters************************************************\n"
180 conf+="***********************************************************************************\n"
181 conf+=fortInt(self.ncall) + " ! [ncall] : Number of calls for preconditioning \n"
182 conf+=fortInt(self.itmx) + " ! [itmx] : Number of iterations for preconditioning \n"
183 conf+=fortDouble(self.prec) + " ! [prec] : Relative accuracy (in %) in main run \n"
184 conf+=fortInt(self.ncall1) + " ! [ncall1] : Number of calls in first iteration \n"
185 conf+=fortInt(self.inccall) + " ! [inccall] : Number of increase calls per iteration \n"
186 conf+=fortInt(self.itend) + " ! [itend] : Maximum number of iterations \n"
187 conf+=fortInt(self.iseed) + " ! [iseed] : Random number seed (integer > 0) \n"
188 conf+="***********************************************************************************\n"
189 conf+="********************Unweighted events**********************************************\n"
190 conf+="***********************************************************************************\n"
191 conf+=fortBool(self.genunw) + " ! [genunw] : Generate unweighted events \n"
192 conf+=fortInt(int(self.nev*1.01)) + " ! [nev] : Number of events (preferably controlled by maxEvents option in Gen_tf command) \n"
193 conf+=fortStr(self.erec) + " ! [erec] : Event record format ('hepmc','lhe','hepevt') \n"
194 conf+=fortBool(self.readwt) + " ! [readwt] : Set to true to read in pre-calculated maxium weight below \n"
195 conf+=fortDouble(self.wtmax) + " ! [wtmax] : Maximum weight \n"
196 conf+="***********************************************************************************\n"
197 conf+="***********************************************************************************\n"
198 conf+="******************* general cuts ************************************************\n"
199 conf+=fortDouble(self.ymin) + " ! [ymin] : Minimum object rapidity Y_X \n"
200 conf+=fortDouble(self.ymax) + " ! [ymax] : Maximum object rapidity Y_X \n"
201 conf+=fortDouble(self.mmin) + " ! [mmin] : Minimum object mass M_X (redundant for resonance production) \n"
202 conf+=fortDouble(self.mmax) + " ! [mmax] : Maximum object mass M_X (redundant for resonance production) \n"
203 conf+=fortBool(self.gencuts) + " ! [gencuts] : Generate cuts below \n"
204 conf+=fortBool(self.scorr) + " ! [scorr] : Include spin correlations (for chi_c/b decays) \n"
205 conf+=fortBool(self.fwidth) + " ! [fwidth] : Include finite width (for chi_c decays) \n"
206 conf+="***********************************************************************************\n"
207 conf+="************ See manual for momentum assignments***********************************\n"
208 conf+="***********************************************************************************\n"
209 conf+="******************* Proton Cuts ***************************************************\n"
210 conf+="***********************************************************************************\n"
211 conf+=fortDouble(self.ptxmax) + " ! [ptxmax] : max pT of the system \n"
212 conf+="***********************************************************************************\n"
213 conf+="**********2 body final states : p(a) + p(b) ***************************************\n"
214 conf+="***********************************************************************************\n"
215 conf+=fortDouble(self.ptamin) + " ! [ptamin] \n"
216 conf+=fortDouble(self.ptbmin) + " ! [ptbmin] \n"
217 conf+=fortDouble(self.etaamin) + " ! [etaamin] \n"
218 conf+=fortDouble(self.etaamax) + " ! [etaamax] \n"
219 conf+=fortDouble(self.etabmin) + " ! [etabmin] \n"
220 conf+=fortDouble(self.etabmax) + " ! [etabmax] \n"
221 conf+=fortDouble(self.acoabmax) + " ! [acoabmax] \n"
222 conf+="***********************************************************************************\n"
223 conf+="****** 3 body final states : p(a) + p(b) + p(c) ***********************************\n"
224 conf+="***********************************************************************************\n"
225 conf+=fortDouble(self.ptamin) + " ! [ptamin] \n"
226 conf+=fortDouble(self.ptbmin) + " ! [ptbmin] \n"
227 conf+=fortDouble(self.ptcmin) + " ! [ptcmin] \n"
228 conf+=fortDouble(self.etaamin) + " ! [etaamin] \n"
229 conf+=fortDouble(self.etaamax) + " ! [etaamax] \n"
230 conf+=fortDouble(self.etabmin) + " ! [etabmin] \n"
231 conf+=fortDouble(self.etabmax) + " ! [etabmax] \n"
232 conf+=fortDouble(self.etacmin) + " ! [etacmin] \n"
233 conf+=fortDouble(self.etacmax) + " ! [etacmax] \n"
234 conf+="***********************************************************************************\n"
235 conf+="****** 4 body final states : p(a) + p(b) + p(c) + p(d) ****************************\n"
236 conf+="***********************************************************************************\n"
237 conf+=fortDouble(self.ptamin) + " ! [ptamin] \n"
238 conf+=fortDouble(self.ptbmin) + " ! [ptbmin] \n"
239 conf+=fortDouble(self.ptcmin) + " ! [ptcmin] \n"
240 conf+=fortDouble(self.ptdmin) + " ! [ptdmin] \n"
241 conf+=fortDouble(self.etaamin) + " ! [etaamin] \n"
242 conf+=fortDouble(self.etaamax) + " ! [etaamax] \n"
243 conf+=fortDouble(self.etabmin) + " ! [etabmin] \n"
244 conf+=fortDouble(self.etabmax) + " ! [etabmax] \n"
245 conf+=fortDouble(self.etacmin) + " ! [etacmin] \n"
246 conf+=fortDouble(self.etacmax) + " ! [etacmax] \n"
247 conf+=fortDouble(self.etadmin) + " ! [etadmin] \n"
248 conf+=fortDouble(self.etadmax) + " ! [etadmax] \n"
249 conf+="***********************************************************************************\n"
250 conf+="****** 6 body final states : p(a) + p(b) + p(c) + p(d) + p(e) + p(f) **************\n"
251 conf+="***********************************************************************************\n"
252 conf+=fortDouble(self.ptamin) + " ! [ptamin] \n"
253 conf+=fortDouble(self.ptbmin) + " ! [ptbmin] \n"
254 conf+=fortDouble(self.ptcmin) + " ! [ptcmin] \n"
255 conf+=fortDouble(self.ptdmin) + " ! [ptdmin] \n"
256 conf+=fortDouble(self.ptemin) + " ! [ptemin] \n"
257 conf+=fortDouble(self.ptfmin) + " ! [ptfmin] \n"
258 conf+=fortDouble(self.etaamin) + " ! [etaamin] \n"
259 conf+=fortDouble(self.etaamax) + " ! [etaamax] \n"
260 conf+=fortDouble(self.etabmin) + " ! [etabmin] \n"
261 conf+=fortDouble(self.etabmax) + " ! [etabmax] \n"
262 conf+=fortDouble(self.etacmin) + " ! [etacmin] \n"
263 conf+=fortDouble(self.etacmax) + " ! [etacmax] \n"
264 conf+=fortDouble(self.etadmin) + " ! [etadmin] \n"
265 conf+=fortDouble(self.etadmax) + " ! [etadmax] \n"
266 conf+=fortDouble(self.etaemin) + " ! [etaemin] \n"
267 conf+=fortDouble(self.etaemax) + " ! [etaemax] \n"
268 conf+=fortDouble(self.etaemin) + " ! [etafmin] \n"
269 conf+=fortDouble(self.etafmax) + " ! [etafmax] \n"
270 conf+="***********************************************************************************\n"
271 conf+="******************* jet algorithm parameters **************************************\n"
272 conf+="***********************************************************************************\n"
273 conf+=fortDouble(self.rjet) + " ! [rjet] : Jet Radius \n"
274 conf+=fortStr(self.jalg) + " ! [jalg] : Jet algorithm ('antikt','kt') \n"
275 conf+="***********************************************************************************\n"
276 conf+="***** chi_c/b two-body decays *****************************************************\n"
277 conf+="***********************************************************************************\n"
278 conf+=fortDouble(self.m2b) + " ! [m2b] : mass of decay particles \n"
279 conf+=fortInt(self.pdgid1) + " ! [pdgid1] : PDG number particle 1 \n"
280 conf+=fortInt(self.pdgid2) + " ! [pdgid2] : PDG number particle 2 \n"
281 conf+="***********************************************************************************\n"
282 conf+="******* ALP Parameters ***********************************************************\n"
283 conf+="***********************************************************************************\n"
284 conf+=fortDouble(self.malp) + " ! [malp] : ALP mass (GeV) \n"
285 conf+=fortDouble(self.gax) + " ! [gax] : ALP coupling (GeV^-1) \n"
286 conf+=fortStr(self.alpt) + " ! [alpt] : AlP type (scalar - 'sc', pseudoscalar - 'ps') \n"
287 conf+="***********************************************************************************\n"
288 conf+="**** Monopole Parameters **********************************************************\n"
289 conf+="***********************************************************************************\n"
290 conf+=fortDouble(self.mpol) + " ! [mpol] : Monopole mass \n"
291 conf+=fortDouble(self.mmon) + " ! [mmon] : Monopolium mass \n"
292 conf+=fortDouble(self.gamm) + " ! [gamm] : Monopolium width \n"
293 conf+="***********************************************************************************\n"
294 conf+="**** SUSY Parameters ************************************************************\n"
295 conf+="***********************************************************************************\n"
296 conf+=fortDouble(self.mcharg) + " ! [mcharg] : Chargino/Slepton mass \n"
297 conf+=fortDouble(self.mneut) + " ! [mneut] : Neutralino mass \n"
298 conf+="***********************************************************************************\n"
299 conf+="***********************************************************************************\n"
300 conf+="***********************************************************************************\n"
301 conf+=fortStr(self.wlp) + " ! [wlp] : leptonic decay (either 'mu' or 'el') for Wplus \n"
302 conf+=fortStr(self.wlm) + " ! [wlm] : leptonic decay (either 'mu' or 'el') for Wminus \n"
303
304 conf+="****************************************************************************************\n"
305 conf+="****** V+X simplified model \n"
306 conf+="****************************************************************************************\n"
307 conf+=fortDouble(self.tau) + " ! [tau] : mass distribution decay constant (GeV^-1) \n"
308 conf+=fortDouble(self.mxs) + " ! [mxs] : mass of MX \n"
309
310 conf+="****************************************************************************************\n"
311 conf+="****** tau anomalous moments\n"
312 conf+="****************************************************************************************\n"
313 conf+=fortDouble(self.atau) + " ! [atau] : magnetic dipole moment\n"
314 conf+=fortDouble(self.dtau) + " ! [dtau] : electric dipole moment [e cm]\n"
315 conf+="****** flags for calculating individual coeffs - SEE MANUAL for explanation\n"
316 conf+=fortBool(self.calc_tau_coeff) + " ! [calc_tau_coeff] : if true calculate O(a_tau^n) or O(d_tau^n) coefficients\n"
317 conf+=fortStr(self.atau) + " ! [tau_mom] : 'atau','dtau' - coeffecients for magnetic/electric dipole moments\n"
318 conf+=fortInt(self.tau_coeff) + " ! [tau_coeff] : order 'n' in coefficient (0-4)\n"
319
320 return conf
321
322 def outputLHEFile(self):
323 return "evrecs/evrec"+self.outtg+".dat"
324
325
326def writeInputDAT(Init):
327
328 with open("input.DAT", "w") as outF:
329 outF.write(Init.toFortran())
330
331 return
332
333
334def run_command(command, stdin = None):
335 """
336 Run a command and print output continuously
337 """
338 process = subprocess.Popen(shlex.split(command), stdout=subprocess.PIPE, stdin=stdin)
339 while True:
340 output = process.stdout.readline().decode("utf-8")
341 if output == '' and process.poll() is not None:
342 break
343 if output:
344 # remove ANSI escape formatting characters
345 reaesc = re.compile(r'(\x9B|\x1B\‍[)[0-?]*[ -\/]*[@-~]')
346 text = reaesc.sub('', output.strip())
347 logger.info(text)
348
349 rc = process.poll()
350 return rc
351
352
353def SuperChicInitialize(Init, stdin=None):
354
355 logger.info("Starting SuperChic Initialization")
356
357 os.makedirs('inputs', exist_ok=True)
358 os.makedirs('evrecs', exist_ok=True)
359 os.makedirs('outputs', exist_ok=True)
360 if not os.path.exists("param_card.dat"):
361 if os.path.exists(Init.superchicpath+"/share/doc/SuperChic/Cards/param_card.dat"):
362 shutil.copyfile(Init.superchicpath+"/share/doc/SuperChic/Cards/param_card.dat","param_card.dat")
363 else:
364 raise Exception('Unexpected error in superchic init execution: probably absent param_card.dat')
365 if not os.path.exists("ident_card.dat"):
366 if os.path.exists(Init.superchicpath+"/share/doc/SuperChic/Cards/ident_card.dat"):
367 shutil.copyfile(Init.superchicpath+"/share/doc/SuperChic/Cards/ident_card.dat","ident_card.dat")
368 else:
369 raise Exception('Unexpected error in superchic init execution: probably absent ident_card.dat')
370
371 try:
372 inputDAT = open('input.DAT')
373
374 except IOError:
375 raise Exception("problem with file IO; potentially input.DAT not created correctly")
376 else:
377
378 try:
379 rc = run_command(Init.superchicpath+"/bin/init", inputDAT)
380
381 except OSError:
382 raise Exception("init executable or file not found")
383
384 except Exception:
385 raise Exception("Non-OSError or IOError in init execution block")
386
387 if rc:
388 raise Exception('Unexpected error in superchic init execution')
389
390 return
391
392
393def SuperChicExecute(Init):
394
395 logger.info("Starting SuperChic Itself")
396 os.makedirs('inputs', exist_ok=True)
397 os.makedirs('evrecs', exist_ok=True)
398 os.makedirs('outputs', exist_ok=True)
399
400 try:
401 inputDAT = open('input.DAT')
402 except IOError:
403 raise Exception ("problem with IO; potentially input.DAT not created correctly")
404 else:
405
406 try:
407 rc = run_command(Init.superchicpath+'/bin/superchic', stdin=inputDAT)
408
409 except OSError:
410 raise Exception("Superchic executable or file not found")
411
412 except Exception:
413 raise Exception("Non-OSError or IOError in Superchic execution block")
414
415 if rc:
416 raise Exception('Unexpected error in superchic execution')
417
418 return
419
420
421
422def SuperChicRun(Init, genSeq):
423
424 # dump the job configuration for fortran code
425 print(Init.toFortran())
426
427 # attach config to genSequence for later usin in showering
428 genSeq.SuperChicConfig = Init
429
430 writeInputDAT(Init)
431 SuperChicInitialize(Init)
432 SuperChicExecute(Init)
433
434 return
print
void print(char *figname, TCanvas *c1)
Definition TRTCalib_StrawStatusPlots.cxx:26
SuperChicUtils.SuperChicInitialize
SuperChicInitialize(Init, stdin=None)
Definition SuperChicUtils.py:353
SuperChicUtils.SuperChicRun
SuperChicRun(Init, genSeq)
Definition SuperChicUtils.py:422
SuperChicUtils.run_command
run_command(command, stdin=None)
Definition SuperChicUtils.py:334
Generated on for ATLAS Offline Software by doxygen 1.16.1