4 import subprocess, os, shlex, re
6 from AthenaCommon
import Logging
9 logger = Logging.logging.getLogger(
"Superchic_i")
19 if hasattr(runArgs,
"ecmEnergy"):
20 self.
rts = runArgs.ecmEnergy
49 if hasattr(runArgs,
"randomSeed"):
50 self.
iseed = runArgs.randomSeed
55 if hasattr(runArgs,
"maxEvents"):
56 self.
nev = runArgs.maxEvents
113 return '.true.' if x
else '.false.'
118 conf+=
"***********************************************************************************\n"
119 conf+=
"****** Initialize afain if FIRST FIVE PARAMETERS ARE CHANGED (and beam = 'prot'):******\n"
120 conf+=
"***********************************************************************************\n"
121 conf+=fortDouble(self.
rts) +
" ! [rts] : CMS collision energy (GeV) \n"
122 conf+=fortInt(self.
isurv) +
" ! [isurv] : Model of soft survival (from 1 -> 4)\n"
123 conf+=fortStr(self.
intag) +
" ! [intag] for input files \n"
124 conf+=
"***********************************************************************************\n"
125 conf+=
"***********************************************************************************\n"
126 conf+=fortStr(self.
PDFname) +
" ! [PDFname] : PDF set \n"
127 conf+=fortInt(self.
PDFmember) +
" ! [PDFmember] : PDF member \n"
128 conf+=
"***********************************************************************************\n"
129 conf+=fortInt(self.
proc) +
" ! [proc] : Process number \n"
130 conf+=fortStr(self.
beam) +
" ! [beam] : Beam type ('prot', 'ion', 'ionp' or 'el') \n"
131 conf+=fortStr(self.
outtg) +
" ! [outtg] : for output file \n"
132 conf+=fortBool(self.
sfaci) +
" ! [sfaci] : Include soft survival effects \n"
133 conf+=
"***********************************************************************************\n"
134 conf+=
"***********************************************************************************\n"
135 conf+=
"***********************************************************************************\n"
136 conf+=fortStr(self.
diff) +
" ! [diff] : dissociation flag \n"
137 conf+=
"***********************************************************************************\n"
138 conf+=fortDouble(self.
an) +
" ! [an] : Ion mass number \n"
139 conf+=fortDouble(self.
az) +
" ! [az] : Ion atomic number \n"
140 conf+=fortDouble(self.
rz) +
" ! [rz] : Ion proton density - radius \n"
141 conf+=fortDouble(self.
dz) +
" ! [dz] : Ion proton density - skin thickness \n"
142 conf+=fortDouble(self.
rn) +
" ! [rn] : Ion neutron density - radius \n"
143 conf+=fortDouble(self.
dn) +
" ! [dn] : Ion neutron density - skin thickness \n"
144 conf+=fortStr(self.
ionqcd) +
" ! [ionqcd] : Coherent ('coh') or incoherent ('incoh') for QCD-induced processes \n"
145 conf+=fortBool(self.
ionbreakup) +
" ! [ionbreakup] \n"
146 conf+=fortStr(self.
fAA) +
" ! [fAA] \n"
147 conf+=fortDouble(self.
fracsigX) +
" ! [fracsigX] : multiply sig_(gamA) by this factor (1d0 - default) \n"
148 conf+=
"***********************************************************************************\n"
149 conf+=
"*************Integration parameters************************************************\n"
150 conf+=
"***********************************************************************************\n"
151 conf+=fortInt(self.
ncall) +
" ! [ncall] : Number of calls for preconditioning \n"
152 conf+=fortInt(self.
itmx) +
" ! [itmx] : Number of iterations for preconditioning \n"
153 conf+=fortDouble(self.
prec) +
" ! [prec] : Relative accuracy (in %) in main run \n"
154 conf+=fortInt(self.
ncall1) +
" ! [ncall1] : Number of calls in first iteration \n"
155 conf+=fortInt(self.
inccall) +
" ! [inccall] : Number of increase calls per iteration \n"
156 conf+=fortInt(self.
itend) +
" ! [itend] : Maximum number of iterations \n"
157 conf+=fortInt(self.
iseed) +
" ! [iseed] : Random number seed (integer > 0) \n"
158 conf+=
"***********************************************************************************\n"
159 conf+=
"********************Unweighted events**********************************************\n"
160 conf+=
"***********************************************************************************\n"
161 conf+=fortBool(self.
genunw) +
" ! [genunw] : Generate unweighted events \n"
162 conf+=fortInt(
int(self.
nev*1.01)) +
" ! [nev] : Number of events (preferably controlled by maxEvents option in Gen_tf command) \n"
163 conf+=fortStr(self.
erec) +
" ! [erec] : Event record format ('hepmc','lhe','hepevt') \n"
164 conf+=fortBool(self.
readwt) +
" ! [readwt] : Set to true to read in pre-calculated maxium weight below \n"
165 conf+=fortDouble(self.
wtmax) +
" ! [wtmax] : Maximum weight \n"
166 conf+=
"***********************************************************************************\n"
167 conf+=
"***********************************************************************************\n"
168 conf+=
"******************* general cuts ************************************************\n"
169 conf+=fortDouble(self.
ymin) +
" ! [ymin] : Minimum object rapidity Y_X \n"
170 conf+=fortDouble(self.
ymax) +
" ! [ymax] : Maximum object rapidity Y_X \n"
171 conf+=fortDouble(self.
mmin) +
" ! [mmin] : Minimum object mass M_X (redundant for resonance production) \n"
172 conf+=fortDouble(self.
mmax) +
" ! [mmax] : Maximum object mass M_X (redundant for resonance production) \n"
173 conf+=fortBool(self.
gencuts) +
" ! [gencuts] : Generate cuts below \n"
174 conf+=fortBool(self.
scorr) +
" ! [scorr] : Include spin correlations (for chi_c/b decays) \n"
175 conf+=fortBool(self.
fwidth) +
" ! [fwidth] : Include finite width (for chi_c decays) \n"
176 conf+=
"***********************************************************************************\n"
177 conf+=
"************ See manual for momentum assignments***********************************\n"
178 conf+=
"***********************************************************************************\n"
179 conf+=
"******************* Proton Cuts ***************************************************\n"
180 conf+=
"***********************************************************************************\n"
181 conf+=fortDouble(self.
ptxmax) +
" ! [ptxmax] : max pT of the system \n"
182 conf+=
"***********************************************************************************\n"
183 conf+=
"**********2 body final states : p(a) + p(b) ***************************************\n"
184 conf+=
"***********************************************************************************\n"
185 conf+=fortDouble(self.
ptamin) +
" ! [ptamin] \n"
186 conf+=fortDouble(self.
ptbmin) +
" ! [ptbmin] \n"
187 conf+=fortDouble(self.
etaamin) +
" ! [etaamin] \n"
188 conf+=fortDouble(self.
etaamax) +
" ! [etaamax] \n"
189 conf+=fortDouble(self.
etabmin) +
" ! [etabmin] \n"
190 conf+=fortDouble(self.
etabmax) +
" ! [etabmax] \n"
191 conf+=fortDouble(self.
acoabmax) +
" ! [acoabmax] \n"
192 conf+=
"***********************************************************************************\n"
193 conf+=
"****** 3 body final states : p(a) + p(b) + p(c) ***********************************\n"
194 conf+=
"***********************************************************************************\n"
195 conf+=fortDouble(self.
ptamin) +
" ! [ptamin] \n"
196 conf+=fortDouble(self.
ptbmin) +
" ! [ptbmin] \n"
197 conf+=fortDouble(self.
ptcmin) +
" ! [ptcmin] \n"
198 conf+=fortDouble(self.
etaamin) +
" ! [etaamin] \n"
199 conf+=fortDouble(self.
etaamax) +
" ! [etaamax] \n"
200 conf+=fortDouble(self.
etabmin) +
" ! [etabmin] \n"
201 conf+=fortDouble(self.
etabmax) +
" ! [etabmax] \n"
202 conf+=fortDouble(self.
etacmin) +
" ! [etacmin] \n"
203 conf+=fortDouble(self.
etacmax) +
" ! [etacmax] \n"
204 conf+=
"***********************************************************************************\n"
205 conf+=
"****** 4 body final states : p(a) + p(b) + p(c) + p(d) ****************************\n"
206 conf+=
"***********************************************************************************\n"
207 conf+=fortDouble(self.
ptamin) +
" ! [ptamin] \n"
208 conf+=fortDouble(self.
ptbmin) +
" ! [ptbmin] \n"
209 conf+=fortDouble(self.
ptcmin) +
" ! [ptcmin] \n"
210 conf+=fortDouble(self.
ptdmin) +
" ! [ptdmin] \n"
211 conf+=fortDouble(self.
etaamin) +
" ! [etaamin] \n"
212 conf+=fortDouble(self.
etaamax) +
" ! [etaamax] \n"
213 conf+=fortDouble(self.
etabmin) +
" ! [etabmin] \n"
214 conf+=fortDouble(self.
etabmax) +
" ! [etabmax] \n"
215 conf+=fortDouble(self.
etacmin) +
" ! [etacmin] \n"
216 conf+=fortDouble(self.
etacmax) +
" ! [etacmax] \n"
217 conf+=fortDouble(self.
etadmin) +
" ! [etadmin] \n"
218 conf+=fortDouble(self.
etadmax) +
" ! [etadmax] \n"
219 conf+=
"***********************************************************************************\n"
220 conf+=
"****** 6 body final states : p(a) + p(b) + p(c) + p(d) + p(e) + p(f) **************\n"
221 conf+=
"***********************************************************************************\n"
222 conf+=fortDouble(self.
ptamin) +
" ! [ptamin] \n"
223 conf+=fortDouble(self.
ptbmin) +
" ! [ptbmin] \n"
224 conf+=fortDouble(self.
ptcmin) +
" ! [ptcmin] \n"
225 conf+=fortDouble(self.
ptdmin) +
" ! [ptdmin] \n"
226 conf+=fortDouble(self.
ptemin) +
" ! [ptemin] \n"
227 conf+=fortDouble(self.
ptfmin) +
" ! [ptfmin] \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+=fortDouble(self.
etadmin) +
" ! [etadmin] \n"
235 conf+=fortDouble(self.
etadmax) +
" ! [etadmax] \n"
236 conf+=fortDouble(self.
etaemin) +
" ! [etaemin] \n"
237 conf+=fortDouble(self.
etaemax) +
" ! [etaemax] \n"
238 conf+=fortDouble(self.
etaemin) +
" ! [etafmin] \n"
239 conf+=fortDouble(self.
etafmax) +
" ! [etafmax] \n"
240 conf+=
"***********************************************************************************\n"
241 conf+=
"******************* jet algorithm parameters **************************************\n"
242 conf+=
"***********************************************************************************\n"
243 conf+=fortDouble(self.
rjet) +
" ! [rjet] : Jet Radius \n"
244 conf+=fortStr(self.
jalg) +
" ! [jalg] : Jet algorithm ('antikt','kt') \n"
245 conf+=
"***********************************************************************************\n"
246 conf+=
"***** chi_c/b two-body decays *****************************************************\n"
247 conf+=
"***********************************************************************************\n"
248 conf+=fortDouble(self.
m2b) +
" ! [m2b] : mass of decay particles \n"
249 conf+=fortInt(self.
pdgid1) +
" ! [pdgid1] : PDG number particle 1 \n"
250 conf+=fortInt(self.
pdgid2) +
" ! [pdgid2] : PDG number particle 2 \n"
251 conf+=
"***********************************************************************************\n"
252 conf+=
"******* ALP Parameters ***********************************************************\n"
253 conf+=
"***********************************************************************************\n"
254 conf+=fortDouble(self.
malp) +
" ! [malp] : ALP mass (GeV) \n"
255 conf+=fortDouble(self.
gax) +
" ! [gax] : ALP coupling (GeV^-1) \n"
256 conf+=fortStr(self.
alpt) +
" ! [alpt] : AlP type (scalar - 'sc', pseudoscalar - 'ps') \n"
257 conf+=
"***********************************************************************************\n"
258 conf+=
"**** Monopole Parameters **********************************************************\n"
259 conf+=
"***********************************************************************************\n"
260 conf+=fortDouble(self.
mpol) +
" ! [mpol] : Monopole mass \n"
261 conf+=fortDouble(self.
mmon) +
" ! [mmon] : Monopolium mass \n"
262 conf+=fortDouble(self.
gamm) +
" ! [gamm] : Monopolium width \n"
263 conf+=
"***********************************************************************************\n"
264 conf+=
"**** SUSY Parameters ************************************************************\n"
265 conf+=
"***********************************************************************************\n"
266 conf+=fortDouble(self.
mcharg) +
" ! [mcharg] : Chargino/Slepton mass \n"
267 conf+=fortDouble(self.
mneut) +
" ! [mneut] : Neutralino mass \n"
268 conf+=
"***********************************************************************************\n"
269 conf+=
"***********************************************************************************\n"
270 conf+=
"***********************************************************************************\n"
271 conf+=fortStr(self.
wlp) +
" ! [wlp] : leptonic decay (either 'mu' or 'el') for Wplus \n"
272 conf+=fortStr(self.
wlm) +
" ! [wlm] : leptonic decay (either 'mu' or 'el') for Wminus \n"
274 conf+=
"****************************************************************************************\n"
275 conf+=
"****** V+X simplified model \n"
276 conf+=
"****************************************************************************************\n"
277 conf+=fortDouble(self.
tau) +
" ! [tau] : mass distribution decay constant (GeV^-1) \n"
278 conf+=fortDouble(self.
mxs) +
" ! [mxs] : mass of MX \n"
284 return "evrecs/evrec"+self.
outtg+
".dat"
289 with open(
"input.DAT",
"w")
as outF:
290 outF.write(Init.toFortran())
297 Run a command and print output continuously
299 process = subprocess.Popen(shlex.split(command), stdout=subprocess.PIPE, stdin=stdin)
301 output = process.stdout.readline().
decode(
"utf-8")
302 if output ==
'' and process.poll()
is not None:
306 reaesc = re.compile(
r'(\x9B|\x1B\[)[0-?]*[ -\/]*[@-~]')
307 text = reaesc.sub(
'', output.strip())
316 logger.info(
"Starting SuperChic Initialization")
318 if not os.path.exists(
'inputs'):
319 os.makedirs(
'inputs')
320 if not os.path.exists(
'evrecs'):
321 os.makedirs(
'evrecs')
322 if not os.path.exists(
'outputs'):
323 os.makedirs(
'outputs')
327 inputDAT =
open(
'input.DAT')
330 raise Exception(
"problem with file IO; potentially input.DAT not created correctly")
334 rc =
run_command(Init.superchicpath+
"/bin/init", inputDAT)
337 raise Exception(
"init executable or file not found")
340 raise Exception(
"Non-OSError or IOError in init execution block")
343 raise Exception(
'Unexpected error in superchic init execution')
350 logger.info(
"Starting SuperChic Itself")
352 if not os.path.exists(
'inputs'):
353 os.makedirs(
'inputs')
354 if not os.path.exists(
'evrecs'):
355 os.makedirs(
'evrecs')
356 if not os.path.exists(
'outputs'):
357 os.makedirs(
'outputs')
361 inputDAT =
open(
'input.DAT')
363 raise Exception (
"problem with IO; potentially input.DAT not created correctly")
367 rc =
run_command(Init.superchicpath+
'/bin/superchic', stdin=inputDAT)
370 raise Exception(
"Superchic executable or file not found")
373 raise Exception(
"Non-OSError or IOError in Superchic execution block")
376 raise Exception(
'Unexpected error in superchic execution')
385 print(Init.toFortran())
388 genSeq.SuperChicConfig = Init