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