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