################################################################################
#
# Copyright (c) 2011 The MadGraph5_aMC@NLO Development team and Contributors
#
# This file is a part of the MadGraph5_aMC@NLO project, an application which
# automatically generates Feynman diagrams and matrix elements for arbitrary
# high-energy processes in the Standard Model and beyond.
#
# It is subject to the MadGraph5_aMC@NLO license which should accompany this
# distribution.
#
# For more information, visit madgraph.phys.ucl.ac.be and amcatnlo.web.cern.ch
#
################################################################################
"""A user friendly command line interface to access MadGraph5_aMC@NLO features.
Uses the cmd package for command interpretation and tab completion.
"""
from __future__ import division
import atexit
import glob
import logging
import math
import optparse
import os
import pydoc
import random
import re
import shutil
import subprocess
import sys
import traceback
import time
import signal
import tarfile
import copy
import datetime
import tarfile
import traceback
import StringIO
try:
import cpickle as pickle
except:
import pickle
try:
import readline
GNU_SPLITTING = ('GNU' in readline.__doc__)
except:
GNU_SPLITTING = True
root_path = os.path.split(os.path.dirname(os.path.realpath( __file__ )))[0]
root_path = os.path.split(root_path)[0]
sys.path.insert(0, os.path.join(root_path,'bin'))
# usefull shortcut
pjoin = os.path.join
# Special logger for the Cmd Interface
logger = logging.getLogger('madgraph.stdout') # -> stdout
logger_stderr = logging.getLogger('madgraph.stderr') # ->stderr
try:
import madgraph
except ImportError:
aMCatNLO = True
import internal.extended_cmd as cmd
import internal.common_run_interface as common_run
import internal.banner as banner_mod
import internal.misc as misc
from internal import InvalidCmd, MadGraph5Error
import internal.files as files
import internal.cluster as cluster
import internal.save_load_object as save_load_object
import internal.gen_crossxhtml as gen_crossxhtml
import internal.sum_html as sum_html
import internal.shower_card as shower_card
import internal.FO_analyse_card as analyse_card
import internal.lhe_parser as lhe_parser
else:
# import from madgraph directory
aMCatNLO = False
import madgraph.interface.extended_cmd as cmd
import madgraph.interface.common_run_interface as common_run
import madgraph.iolibs.files as files
import madgraph.iolibs.save_load_object as save_load_object
import madgraph.madevent.gen_crossxhtml as gen_crossxhtml
import madgraph.madevent.sum_html as sum_html
import madgraph.various.banner as banner_mod
import madgraph.various.cluster as cluster
import madgraph.various.misc as misc
import madgraph.various.shower_card as shower_card
import madgraph.various.FO_analyse_card as analyse_card
import madgraph.various.lhe_parser as lhe_parser
from madgraph import InvalidCmd, aMCatNLOError, MadGraph5Error,MG5DIR
class aMCatNLOError(Exception):
pass
def compile_dir(*arguments):
"""compile the direcory p_dir
arguments is the tuple (me_dir, p_dir, mode, options, tests, exe, run_mode)
this function needs not to be a class method in order to do
the compilation on multicore"""
if len(arguments) == 1:
(me_dir, p_dir, mode, options, tests, exe, run_mode) = arguments[0]
elif len(arguments)==7:
(me_dir, p_dir, mode, options, tests, exe, run_mode) = arguments
else:
raise aMCatNLOError, 'not correct number of argument'
logger.info(' Compiling %s...' % p_dir)
this_dir = pjoin(me_dir, 'SubProcesses', p_dir)
try:
#compile everything
# compile and run tests
for test in tests:
# skip check_poles for LOonly dirs
if test == 'check_poles' and os.path.exists(pjoin(this_dir, 'parton_lum_0.f')):
continue
misc.compile([test], cwd = this_dir, job_specs = False)
input = pjoin(me_dir, '%s_input.txt' % test)
#this can be improved/better written to handle the output
misc.call(['./%s' % (test)], cwd=this_dir,
stdin = open(input), stdout=open(pjoin(this_dir, '%s.log' % test), 'w'),
close_fds=True)
if test == 'check_poles' and os.path.exists(pjoin(this_dir,'MadLoop5_resources')) :
tf=tarfile.open(pjoin(this_dir,'MadLoop5_resources.tar.gz'),'w:gz',
dereference=True)
tf.add(pjoin(this_dir,'MadLoop5_resources'),arcname='MadLoop5_resources')
tf.close()
if not options['reweightonly']:
misc.compile(['gensym'], cwd=this_dir, job_specs = False)
open(pjoin(this_dir, 'gensym_input.txt'), 'w').write('%s\n' % run_mode)
misc.call(['./gensym'],cwd= this_dir,
stdin=open(pjoin(this_dir, 'gensym_input.txt')),
stdout=open(pjoin(this_dir, 'gensym.log'), 'w'),
close_fds=True)
#compile madevent_mintMC/mintFO
misc.compile([exe], cwd=this_dir, job_specs = False)
if mode in ['aMC@NLO', 'aMC@LO', 'noshower', 'noshowerLO']:
misc.compile(['reweight_xsec_events'], cwd=this_dir, job_specs = False)
logger.info(' %s done.' % p_dir)
return 0
except MadGraph5Error, msg:
return msg
def check_compiler(options, block=False):
"""check that the current fortran compiler is gfortran 4.6 or later.
If block, stops the execution, otherwise just print a warning"""
msg = 'In order to be able to run at NLO MadGraph5_aMC@NLO, you need to have ' + \
'gfortran 4.6 or later installed.\n%s has been detected\n'+\
'Note that You can still run all MadEvent run without any problem!'
#first check that gfortran is installed
if options['fortran_compiler']:
compiler = options['fortran_compiler']
elif misc.which('gfortran'):
compiler = 'gfortran'
else:
compiler = ''
if 'gfortran' not in compiler:
if block:
raise aMCatNLOError(msg % compiler)
else:
logger.warning(msg % compiler)
else:
curr_version = misc.get_gfortran_version(compiler)
if not ''.join(curr_version.split('.')) >= '46':
if block:
raise aMCatNLOError(msg % (compiler + ' ' + curr_version))
else:
logger.warning(msg % (compiler + ' ' + curr_version))
#===============================================================================
# CmdExtended
#===============================================================================
class CmdExtended(common_run.CommonRunCmd):
"""Particularisation of the cmd command for aMCatNLO"""
#suggested list of command
next_possibility = {
'start': [],
}
debug_output = 'ME5_debug'
error_debug = 'Please report this bug on https://bugs.launchpad.net/mg5amcnlo\n'
error_debug += 'More information is found in \'%(debug)s\'.\n'
error_debug += 'Please attach this file to your report.'
config_debug = 'If you need help with this issue please contact us on https://answers.launchpad.net/mg5amcnlo\n'
keyboard_stop_msg = """stopping all operation
in order to quit MadGraph5_aMC@NLO please enter exit"""
# Define the Error
InvalidCmd = InvalidCmd
ConfigurationError = aMCatNLOError
def __init__(self, me_dir, options, *arg, **opt):
"""Init history and line continuation"""
# Tag allowing/forbiding question
self.force = False
# If possible, build an info line with current version number
# and date, from the VERSION text file
info = misc.get_pkg_info()
info_line = ""
if info and info.has_key('version') and info.has_key('date'):
len_version = len(info['version'])
len_date = len(info['date'])
if len_version + len_date < 30:
info_line = "#* VERSION %s %s %s *\n" % \
(info['version'],
(30 - len_version - len_date) * ' ',
info['date'])
else:
version = open(pjoin(root_path,'MGMEVersion.txt')).readline().strip()
info_line = "#* VERSION %s %s *\n" % \
(version, (24 - len(version)) * ' ')
# Create a header for the history file.
# Remember to fill in time at writeout time!
self.history_header = \
'#************************************************************\n' + \
'#* MadGraph5_aMC@NLO *\n' + \
'#* *\n' + \
"#* * * *\n" + \
"#* * * * * *\n" + \
"#* * * * * 5 * * * * *\n" + \
"#* * * * * *\n" + \
"#* * * *\n" + \
"#* *\n" + \
"#* *\n" + \
info_line + \
"#* *\n" + \
"#* The MadGraph5_aMC@NLO Development Team - Find us at *\n" + \
"#* https://server06.fynu.ucl.ac.be/projects/madgraph *\n" + \
"#* and *\n" + \
"#* http://amcatnlo.cern.ch *\n" + \
'#* *\n' + \
'#************************************************************\n' + \
'#* *\n' + \
'#* Command File for aMCatNLO *\n' + \
'#* *\n' + \
'#* run as ./bin/aMCatNLO.py filename *\n' + \
'#* *\n' + \
'#************************************************************\n'
if info_line:
info_line = info_line[1:]
logger.info(\
"************************************************************\n" + \
"* *\n" + \
"* W E L C O M E to M A D G R A P H 5 *\n" + \
"* a M C @ N L O *\n" + \
"* *\n" + \
"* * * *\n" + \
"* * * * * *\n" + \
"* * * * * 5 * * * * *\n" + \
"* * * * * *\n" + \
"* * * *\n" + \
"* *\n" + \
info_line + \
"* *\n" + \
"* The MadGraph5_aMC@NLO Development Team - Find us at *\n" + \
"* http://amcatnlo.cern.ch *\n" + \
"* *\n" + \
"* Type 'help' for in-line help. *\n" + \
"* *\n" + \
"************************************************************")
super(CmdExtended, self).__init__(me_dir, options, *arg, **opt)
def get_history_header(self):
"""return the history header"""
return self.history_header % misc.get_time_info()
def stop_on_keyboard_stop(self):
"""action to perform to close nicely on a keyboard interupt"""
try:
if hasattr(self, 'cluster'):
logger.info('rm jobs on queue')
self.cluster.remove()
if hasattr(self, 'results'):
self.update_status('Stop by the user', level=None, makehtml=True, error=True)
self.add_error_log_in_html(KeyboardInterrupt)
except:
pass
def postcmd(self, stop, line):
""" Update the status of the run for finishing interactive command """
# relaxing the tag forbidding question
self.force = False
if not self.use_rawinput:
return stop
arg = line.split()
if len(arg) == 0:
return stop
elif str(arg[0]) in ['exit','quit','EOF']:
return stop
try:
self.update_status('Command \'%s\' done. Waiting for instruction.' % arg[0],
level=None, error=True)
except Exception:
misc.sprint('self.update_status fails', log=logger)
pass
def nice_user_error(self, error, line):
"""If a ME run is currently running add a link in the html output"""
self.add_error_log_in_html()
cmd.Cmd.nice_user_error(self, error, line)
def nice_config_error(self, error, line):
"""If a ME run is currently running add a link in the html output"""
self.add_error_log_in_html()
cmd.Cmd.nice_config_error(self, error, line)
def nice_error_handling(self, error, line):
"""If a ME run is currently running add a link in the html output"""
self.add_error_log_in_html()
cmd.Cmd.nice_error_handling(self, error, line)
#===============================================================================
# HelpToCmd
#===============================================================================
class HelpToCmd(object):
""" The Series of help routine for the aMCatNLOCmd"""
def help_launch(self):
"""help for launch command"""
_launch_parser.print_help()
def help_banner_run(self):
logger.info("syntax: banner_run Path|RUN [--run_options]")
logger.info("-- Reproduce a run following a given banner")
logger.info(" One of the following argument is require:")
logger.info(" Path should be the path of a valid banner.")
logger.info(" RUN should be the name of a run of the current directory")
self.run_options_help([('-f','answer all question by default'),
('--name=X', 'Define the name associated with the new run')])
def help_compile(self):
"""help for compile command"""
_compile_parser.print_help()
def help_generate_events(self):
"""help for generate_events commandi
just call help_launch"""
_generate_events_parser.print_help()
def help_calculate_xsect(self):
"""help for generate_events command"""
_calculate_xsect_parser.print_help()
def help_shower(self):
"""help for shower command"""
_shower_parser.print_help()
def help_open(self):
logger.info("syntax: open FILE ")
logger.info("-- open a file with the appropriate editor.")
logger.info(' If FILE belongs to index.html, param_card.dat, run_card.dat')
logger.info(' the path to the last created/used directory is used')
def run_options_help(self, data):
if data:
logger.info('-- local options:')
for name, info in data:
logger.info(' %s : %s' % (name, info))
logger.info("-- session options:")
logger.info(" Note that those options will be kept for the current session")
logger.info(" --cluster : Submit to the cluster. Current cluster: %s" % self.options['cluster_type'])
logger.info(" --multicore : Run in multi-core configuration")
logger.info(" --nb_core=X : limit the number of core to use to X.")
#===============================================================================
# CheckValidForCmd
#===============================================================================
class CheckValidForCmd(object):
""" The Series of check routine for the aMCatNLOCmd"""
def check_shower(self, args, options):
"""Check the validity of the line. args[0] is the run_directory"""
if options['force']:
self.force = True
if len(args) == 0:
self.help_shower()
raise self.InvalidCmd, 'Invalid syntax, please specify the run name'
if not os.path.isdir(pjoin(self.me_dir, 'Events', args[0])):
raise self.InvalidCmd, 'Directory %s does not exists' % \
pjoin(os.getcwd(), 'Events', args[0])
self.set_run_name(args[0], level= 'shower')
args[0] = pjoin(self.me_dir, 'Events', args[0])
def check_plot(self, args):
"""Check the argument for the plot command
plot run_name modes"""
madir = self.options['madanalysis_path']
td = self.options['td_path']
if not madir or not td:
logger.info('Retry to read configuration file to find madanalysis/td')
self.set_configuration()
madir = self.options['madanalysis_path']
td = self.options['td_path']
if not madir:
error_msg = 'No Madanalysis path correctly set.'
error_msg += 'Please use the set command to define the path and retry.'
error_msg += 'You can also define it in the configuration file.'
raise self.InvalidCmd(error_msg)
if not td:
error_msg = 'No path to td directory correctly set.'
error_msg += 'Please use the set command to define the path and retry.'
error_msg += 'You can also define it in the configuration file.'
raise self.InvalidCmd(error_msg)
if len(args) == 0:
if not hasattr(self, 'run_name') or not self.run_name:
self.help_plot()
raise self.InvalidCmd('No run name currently define. Please add this information.')
args.append('all')
return
if args[0] not in self._plot_mode:
self.set_run_name(args[0], level='plot')
del args[0]
if len(args) == 0:
args.append('all')
elif not self.run_name:
self.help_plot()
raise self.InvalidCmd('No run name currently define. Please add this information.')
for arg in args:
if arg not in self._plot_mode and arg != self.run_name:
self.help_plot()
raise self.InvalidCmd('unknown options %s' % arg)
def check_pgs(self, arg):
"""Check the argument for pythia command
syntax: pgs [NAME]
Note that other option are already remove at this point
"""
# If not pythia-pgs path
if not self.options['pythia-pgs_path']:
logger.info('Retry to read configuration file to find pythia-pgs path')
self.set_configuration()
if not self.options['pythia-pgs_path'] or not \
os.path.exists(pjoin(self.options['pythia-pgs_path'],'src')):
error_msg = 'No pythia-pgs path correctly set.'
error_msg += 'Please use the set command to define the path and retry.'
error_msg += 'You can also define it in the configuration file.'
raise self.InvalidCmd(error_msg)
tag = [a for a in arg if a.startswith('--tag=')]
if tag:
arg.remove(tag[0])
tag = tag[0][6:]
if len(arg) == 0 and not self.run_name:
if self.results.lastrun:
arg.insert(0, self.results.lastrun)
else:
raise self.InvalidCmd('No run name currently define. Please add this information.')
if len(arg) == 1 and self.run_name == arg[0]:
arg.pop(0)
if not len(arg) and \
not os.path.exists(pjoin(self.me_dir,'Events','pythia_events.hep')):
self.help_pgs()
raise self.InvalidCmd('''No file file pythia_events.hep currently available
Please specify a valid run_name''')
lock = None
if len(arg) == 1:
prev_tag = self.set_run_name(arg[0], tag, 'pgs')
filenames = misc.glob('events_*.hep.gz', pjoin(self.me_dir, 'Events', self.run_name))
if not filenames:
raise self.InvalidCmd('No events file corresponding to %s run with tag %s. '% (self.run_name, prev_tag))
else:
input_file = filenames[0]
output_file = pjoin(self.me_dir, 'Events', 'pythia_events.hep')
lock = cluster.asyncrone_launch('gunzip',stdout=open(output_file,'w'),
argument=['-c', input_file],
close_fds=True)
else:
if tag:
self.run_card['run_tag'] = tag
self.set_run_name(self.run_name, tag, 'pgs')
return lock
def check_delphes(self, arg):
"""Check the argument for pythia command
syntax: delphes [NAME]
Note that other option are already remove at this point
"""
# If not pythia-pgs path
if not self.options['delphes_path']:
logger.info('Retry to read configuration file to find delphes path')
self.set_configuration()
if not self.options['delphes_path']:
error_msg = 'No delphes path correctly set.'
error_msg += 'Please use the set command to define the path and retry.'
error_msg += 'You can also define it in the configuration file.'
raise self.InvalidCmd(error_msg)
tag = [a for a in arg if a.startswith('--tag=')]
if tag:
arg.remove(tag[0])
tag = tag[0][6:]
if len(arg) == 0 and not self.run_name:
if self.results.lastrun:
arg.insert(0, self.results.lastrun)
else:
raise self.InvalidCmd('No run name currently define. Please add this information.')
if len(arg) == 1 and self.run_name == arg[0]:
arg.pop(0)
if not len(arg) and \
not os.path.exists(pjoin(self.me_dir,'Events','pythia_events.hep')):
self.help_pgs()
raise self.InvalidCmd('''No file file pythia_events.hep currently available
Please specify a valid run_name''')
if len(arg) == 1:
prev_tag = self.set_run_name(arg[0], tag, 'delphes')
filenames = misc.glob('events_*.hep.gz', pjoin(self.me_dir, 'Events'))
if not filenames:
raise self.InvalidCmd('No events file corresponding to %s run with tag %s.:%s '\
% (self.run_name, prev_tag,
pjoin(self.me_dir,'Events',self.run_name, '%s_pythia_events.hep.gz' % prev_tag)))
else:
input_file = filenames[0]
output_file = pjoin(self.me_dir, 'Events', 'pythia_events.hep')
lock = cluster.asyncrone_launch('gunzip',stdout=open(output_file,'w'),
argument=['-c', input_file],
close_fds=True)
else:
if tag:
self.run_card['run_tag'] = tag
self.set_run_name(self.run_name, tag, 'delphes')
def check_calculate_xsect(self, args, options):
"""check the validity of the line. args is ORDER,
ORDER being LO or NLO. If no mode is passed, NLO is used"""
# modify args in order to be DIR
# mode being either standalone or madevent
if options['force']:
self.force = True
if not args:
args.append('NLO')
return
if len(args) > 1:
self.help_calculate_xsect()
raise self.InvalidCmd, 'Invalid Syntax: Too many argument'
elif len(args) == 1:
if not args[0] in ['NLO', 'LO']:
raise self.InvalidCmd, '%s is not a valid mode, please use "LO" or "NLO"' % args[1]
mode = args[0]
# check for incompatible options/modes
if options['multicore'] and options['cluster']:
raise self.InvalidCmd, 'options -m (--multicore) and -c (--cluster)' + \
' are not compatible. Please choose one.'
def check_generate_events(self, args, options):
"""check the validity of the line. args is ORDER,
ORDER being LO or NLO. If no mode is passed, NLO is used"""
# modify args in order to be DIR
# mode being either standalone or madevent
if not args:
args.append('NLO')
return
if len(args) > 1:
self.help_generate_events()
raise self.InvalidCmd, 'Invalid Syntax: Too many argument'
elif len(args) == 1:
if not args[0] in ['NLO', 'LO']:
raise self.InvalidCmd, '%s is not a valid mode, please use "LO" or "NLO"' % args[1]
mode = args[0]
# check for incompatible options/modes
if options['multicore'] and options['cluster']:
raise self.InvalidCmd, 'options -m (--multicore) and -c (--cluster)' + \
' are not compatible. Please choose one.'
def check_banner_run(self, args):
"""check the validity of line"""
if len(args) == 0:
self.help_banner_run()
raise self.InvalidCmd('banner_run requires at least one argument.')
tag = [a[6:] for a in args if a.startswith('--tag=')]
if os.path.exists(args[0]):
type ='banner'
format = self.detect_card_type(args[0])
if format != 'banner':
raise self.InvalidCmd('The file is not a valid banner.')
elif tag:
args[0] = pjoin(self.me_dir,'Events', args[0], '%s_%s_banner.txt' % \
(args[0], tag))
if not os.path.exists(args[0]):
raise self.InvalidCmd('No banner associates to this name and tag.')
else:
name = args[0]
type = 'run'
banners = misc.glob('*_banner.txt', pjoin(self.me_dir,'Events', args[0]))
if not banners:
raise self.InvalidCmd('No banner associates to this name.')
elif len(banners) == 1:
args[0] = banners[0]
else:
#list the tag and propose those to the user
tags = [os.path.basename(p)[len(args[0])+1:-11] for p in banners]
tag = self.ask('which tag do you want to use?', tags[0], tags)
args[0] = pjoin(self.me_dir,'Events', args[0], '%s_%s_banner.txt' % \
(args[0], tag))
run_name = [arg[7:] for arg in args if arg.startswith('--name=')]
if run_name:
try:
self.exec_cmd('remove %s all banner -f' % run_name)
except Exception:
pass
self.set_run_name(args[0], tag=None, level='parton', reload_card=True)
elif type == 'banner':
self.set_run_name(self.find_available_run_name(self.me_dir))
elif type == 'run':
if not self.results[name].is_empty():
run_name = self.find_available_run_name(self.me_dir)
logger.info('Run %s is not empty so will use run_name: %s' % \
(name, run_name))
self.set_run_name(run_name)
else:
try:
self.exec_cmd('remove %s all banner -f' % run_name)
except Exception:
pass
self.set_run_name(name)
def check_launch(self, args, options):
"""check the validity of the line. args is MODE
MODE being LO, NLO, aMC@NLO or aMC@LO. If no mode is passed, auto is used"""
# modify args in order to be DIR
# mode being either standalone or madevent
if options['force']:
self.force = True
if not args:
args.append('auto')
return
if len(args) > 1:
self.help_launch()
raise self.InvalidCmd, 'Invalid Syntax: Too many argument'
elif len(args) == 1:
if not args[0] in ['LO', 'NLO', 'aMC@NLO', 'aMC@LO','auto']:
raise self.InvalidCmd, '%s is not a valid mode, please use "LO", "NLO", "aMC@NLO" or "aMC@LO"' % args[0]
mode = args[0]
# check for incompatible options/modes
if options['multicore'] and options['cluster']:
raise self.InvalidCmd, 'options -m (--multicore) and -c (--cluster)' + \
' are not compatible. Please choose one.'
if mode == 'NLO' and options['reweightonly']:
raise self.InvalidCmd, 'option -r (--reweightonly) needs mode "aMC@NLO" or "aMC@LO"'
def check_compile(self, args, options):
"""check the validity of the line. args is MODE
MODE being FO or MC. If no mode is passed, MC is used"""
# modify args in order to be DIR
# mode being either standalone or madevent
if options['force']:
self.force = True
if not args:
args.append('MC')
return
if len(args) > 1:
self.help_compile()
raise self.InvalidCmd, 'Invalid Syntax: Too many argument'
elif len(args) == 1:
if not args[0] in ['MC', 'FO']:
raise self.InvalidCmd, '%s is not a valid mode, please use "FO" or "MC"' % args[0]
mode = args[0]
# check for incompatible options/modes
#===============================================================================
# CompleteForCmd
#===============================================================================
class CompleteForCmd(CheckValidForCmd):
""" The Series of help routine for the MadGraphCmd"""
def complete_launch(self, text, line, begidx, endidx):
"""auto-completion for launch command"""
args = self.split_arg(line[0:begidx])
if len(args) == 1:
#return mode
return self.list_completion(text,['LO','NLO','aMC@NLO','aMC@LO'],line)
elif len(args) == 2 and line[begidx-1] == '@':
return self.list_completion(text,['LO','NLO'],line)
else:
opts = []
for opt in _launch_parser.option_list:
opts += opt._long_opts + opt._short_opts
return self.list_completion(text, opts, line)
def complete_banner_run(self, text, line, begidx, endidx, formatting=True):
"Complete the banner run command"
try:
args = self.split_arg(line[0:begidx], error=False)
if args[-1].endswith(os.path.sep):
return self.path_completion(text,
os.path.join('.',*[a for a in args \
if a.endswith(os.path.sep)]))
if len(args) > 1:
# only options are possible
tags = misc.glob('%s_*_banner.txt' % args[1],pjoin(self.me_dir, 'Events' , args[1]))
tags = ['%s' % os.path.basename(t)[len(args[1])+1:-11] for t in tags]
if args[-1] != '--tag=':
tags = ['--tag=%s' % t for t in tags]
else:
return self.list_completion(text, tags)
return self.list_completion(text, tags +['--name=','-f'], line)
# First argument
possibilites = {}
comp = self.path_completion(text, os.path.join('.',*[a for a in args \
if a.endswith(os.path.sep)]))
if os.path.sep in line:
return comp
else:
possibilites['Path from ./'] = comp
run_list = misc.glob(pjoin('*','*_banner.txt'), pjoin(self.me_dir, 'Events'))
run_list = [n.rsplit('/',2)[1] for n in run_list]
possibilites['RUN Name'] = self.list_completion(text, run_list)
return self.deal_multiple_categories(possibilites, formatting)
except Exception, error:
print error
def complete_compile(self, text, line, begidx, endidx):
"""auto-completion for launch command"""
args = self.split_arg(line[0:begidx])
if len(args) == 1:
#return mode
return self.list_completion(text,['FO','MC'],line)
else:
opts = []
for opt in _compile_parser.option_list:
opts += opt._long_opts + opt._short_opts
return self.list_completion(text, opts, line)
def complete_calculate_xsect(self, text, line, begidx, endidx):
"""auto-completion for launch command"""
args = self.split_arg(line[0:begidx])
if len(args) == 1:
#return mode
return self.list_completion(text,['LO','NLO'],line)
else:
opts = []
for opt in _calculate_xsect_parser.option_list:
opts += opt._long_opts + opt._short_opts
return self.list_completion(text, opts, line)
def complete_generate_events(self, text, line, begidx, endidx):
"""auto-completion for generate_events command
call the compeltion for launch"""
self.complete_launch(text, line, begidx, endidx)
def complete_shower(self, text, line, begidx, endidx):
args = self.split_arg(line[0:begidx])
if len(args) == 1:
#return valid run_name
data = misc.glob(pjoin('*','events.lhe.gz', pjoin(self.me_dir, 'Events')))
data = [n.rsplit('/',2)[1] for n in data]
tmp1 = self.list_completion(text, data)
if not self.run_name:
return tmp1
def complete_plot(self, text, line, begidx, endidx):
""" Complete the plot command """
args = self.split_arg(line[0:begidx], error=False)
if len(args) == 1:
#return valid run_name
data = misc.glob(pjoin('*','events.lhe*', pjoin(self.me_dir, 'Events')))
data = [n.rsplit('/',2)[1] for n in data]
tmp1 = self.list_completion(text, data)
if not self.run_name:
return tmp1
if len(args) > 1:
return self.list_completion(text, self._plot_mode)
def complete_pgs(self,text, line, begidx, endidx):
"Complete the pgs command"
args = self.split_arg(line[0:begidx], error=False)
if len(args) == 1:
#return valid run_name
data = misc.glob(pjoin('*', 'events_*.hep.gz'),
pjoin(self.me_dir, 'Events'))
data = [n.rsplit('/',2)[1] for n in data]
tmp1 = self.list_completion(text, data)
if not self.run_name:
return tmp1
else:
tmp2 = self.list_completion(text, self._run_options + ['-f',
'--tag=' ,'--no_default'], line)
return tmp1 + tmp2
else:
return self.list_completion(text, self._run_options + ['-f',
'--tag=','--no_default'], line)
complete_delphes = complete_pgs
class aMCatNLOAlreadyRunning(InvalidCmd):
pass
#===============================================================================
# aMCatNLOCmd
#===============================================================================
class aMCatNLOCmd(CmdExtended, HelpToCmd, CompleteForCmd, common_run.CommonRunCmd):
"""The command line processor of MadGraph"""
# Truth values
true = ['T','.true.',True,'true']
# Options and formats available
_run_options = ['--cluster','--multicore','--nb_core=','--nb_core=2', '-c', '-m']
_generate_options = ['-f', '--laststep=parton', '--laststep=pythia', '--laststep=pgs', '--laststep=delphes']
_calculate_decay_options = ['-f', '--accuracy=0.']
_set_options = ['stdout_level','fortran_compiler','cpp_compiler','timeout']
_plot_mode = ['all', 'parton','shower','pgs','delphes']
_clean_mode = _plot_mode + ['channel', 'banner']
_display_opts = ['run_name', 'options', 'variable']
# survey options, dict from name to type, default value, and help text
# Variables to store object information
web = False
cluster_mode = 0
queue = 'madgraph'
nb_core = None
make_opts_var = {}
next_possibility = {
'start': ['generate_events [OPTIONS]', 'calculate_crossx [OPTIONS]', 'launch [OPTIONS]',
'help generate_events'],
'generate_events': ['generate_events [OPTIONS]', 'shower'],
'launch': ['launch [OPTIONS]', 'shower'],
'shower' : ['generate_events [OPTIONS]']
}
############################################################################
def __init__(self, me_dir = None, options = {}, *completekey, **stdin):
""" add information to the cmd """
self.start_time = 0
CmdExtended.__init__(self, me_dir, options, *completekey, **stdin)
#common_run.CommonRunCmd.__init__(self, me_dir, options)
self.mode = 'aMCatNLO'
self.nb_core = 0
self.prompt = "%s>"%os.path.basename(pjoin(self.me_dir))
self.load_results_db()
self.results.def_web_mode(self.web)
# check that compiler is gfortran 4.6 or later if virtuals have been exported
proc_card = open(pjoin(self.me_dir, 'Cards', 'proc_card_mg5.dat')).read()
if not '[real=QCD]' in proc_card:
check_compiler(self.options, block=True)
############################################################################
def do_shower(self, line):
""" run the shower on a given parton level file """
argss = self.split_arg(line)
(options, argss) = _launch_parser.parse_args(argss)
# check argument validity and normalise argument
options = options.__dict__
options['reweightonly'] = False
self.check_shower(argss, options)
evt_file = pjoin(os.getcwd(), argss[0], 'events.lhe')
self.ask_run_configuration('onlyshower', options)
self.run_mcatnlo(evt_file, options)
self.update_status('', level='all', update_results=True)
################################################################################
def do_plot(self, line):
"""Create the plot for a given run"""
# Since in principle, all plot are already done automaticaly
args = self.split_arg(line)
# Check argument's validity
self.check_plot(args)
logger.info('plot for run %s' % self.run_name)
if not self.force:
self.ask_edit_cards([], args, plot=True)
if any([arg in ['parton'] for arg in args]):
filename = pjoin(self.me_dir, 'Events', self.run_name, 'events.lhe')
if os.path.exists(filename+'.gz'):
misc.gunzip(filename)
if os.path.exists(filename):
logger.info('Found events.lhe file for run %s' % self.run_name)
shutil.move(filename, pjoin(self.me_dir, 'Events', 'unweighted_events.lhe'))
self.create_plot('parton')
shutil.move(pjoin(self.me_dir, 'Events', 'unweighted_events.lhe'), filename)
misc.gzip(filename)
if any([arg in ['all','parton'] for arg in args]):
filename = pjoin(self.me_dir, 'Events', self.run_name, 'MADatNLO.top')
if os.path.exists(filename):
logger.info('Found MADatNLO.top file for run %s' % \
self.run_name)
output = pjoin(self.me_dir, 'HTML',self.run_name, 'plots_parton.html')
plot_dir = pjoin(self.me_dir, 'HTML', self.run_name, 'plots_parton')
if not os.path.isdir(plot_dir):
os.makedirs(plot_dir)
top_file = pjoin(plot_dir, 'plots.top')
files.cp(filename, top_file)
madir = self.options['madanalysis_path']
tag = self.run_card['run_tag']
td = self.options['td_path']
misc.call(['%s/plot' % self.dirbin, madir, td],
stdout = open(pjoin(plot_dir, 'plot.log'),'a'),
stderr = subprocess.STDOUT,
cwd=plot_dir)
misc.call(['%s/plot_page-pl' % self.dirbin,
os.path.basename(plot_dir),
'parton'],
stdout = open(pjoin(plot_dir, 'plot.log'),'a'),
stderr = subprocess.STDOUT,
cwd=pjoin(self.me_dir, 'HTML', self.run_name))
shutil.move(pjoin(self.me_dir, 'HTML',self.run_name ,'plots.html'),
output)
os.remove(pjoin(self.me_dir, 'Events', 'plots.top'))
if any([arg in ['all','shower'] for arg in args]):
filenames = misc.glob('events_*.lhe.gz', pjoin(self.me_dir, 'Events', self.run_name))
if len(filenames) != 1:
filenames = misc.glob('events_*.hep.gz', pjoin(self.me_dir, 'Events', self.run_name))
if len(filenames) != 1:
logger.info('No shower level file found for run %s' % \
self.run_name)
return
filename = filenames[0]
misc.gunzip(filename, keep=True, stdout=pjoin(self.me_dir, 'Events','pythia_events.hep'))
if not os.path.exists(pjoin(self.me_dir, 'Cards', 'pythia_card.dat')):
if aMCatNLO and not self.options['mg5_path']:
raise "plotting NLO HEP file needs MG5 utilities"
files.cp(pjoin(self.options['mg5_path'], 'Template','LO', 'Cards', 'pythia_card_default.dat'),
pjoin(self.me_dir, 'Cards', 'pythia_card.dat'))
self.run_hep2lhe()
else:
filename = filenames[0]
misc.gunzip(filename, keep=True, stdout=pjoin(self.me_dir, 'Events','pythia_events.hep'))
self.create_plot('shower')
lhe_file_name = filename.replace('.hep.gz', '.lhe')
shutil.move(pjoin(self.me_dir, 'Events','pythia_events.lhe'),
lhe_file_name)
misc.gzip(lhe_file_name)
if any([arg in ['all','pgs'] for arg in args]):
filename = pjoin(self.me_dir, 'Events', self.run_name,
'%s_pgs_events.lhco' % self.run_tag)
if os.path.exists(filename+'.gz'):
misc.gunzip(filename)
if os.path.exists(filename):
self.create_plot('PGS')
misc.gzip(filename)
else:
logger.info('No valid files for pgs plot')
if any([arg in ['all','delphes'] for arg in args]):
filename = pjoin(self.me_dir, 'Events', self.run_name,
'%s_delphes_events.lhco' % self.run_tag)
if os.path.exists(filename+'.gz'):
misc.gunzip(filename)
if os.path.exists(filename):
#shutil.move(filename, pjoin(self.me_dir, 'Events','delphes_events.lhco'))
self.create_plot('Delphes')
#shutil.move(pjoin(self.me_dir, 'Events','delphes_events.lhco'), filename)
misc.gzip(filename)
else:
logger.info('No valid files for delphes plot')
############################################################################
def do_calculate_xsect(self, line):
"""Main commands: calculates LO/NLO cross-section, using madevent_mintFO
this function wraps the do_launch one"""
self.start_time = time.time()
argss = self.split_arg(line)
# check argument validity and normalise argument
(options, argss) = _calculate_xsect_parser.parse_args(argss)
options = options.__dict__
options['reweightonly'] = False
options['parton'] = True
self.check_calculate_xsect(argss, options)
self.do_launch(line, options, argss)
############################################################################
def do_banner_run(self, line):
"""Make a run from the banner file"""
args = self.split_arg(line)
#check the validity of the arguments
self.check_banner_run(args)
# Remove previous cards
for name in ['shower_card.dat', 'madspin_card.dat']:
try:
os.remove(pjoin(self.me_dir, 'Cards', name))
except Exception:
pass
banner_mod.split_banner(args[0], self.me_dir, proc_card=False)
# Check if we want to modify the run
if not self.force:
ans = self.ask('Do you want to modify the Cards/Run Type?', 'n', ['y','n'])
if ans == 'n':
self.force = True
# Compute run mode:
if self.force:
mode_status = {'order': 'NLO', 'fixed_order': False, 'madspin':False, 'shower':True}
banner = banner_mod.Banner(args[0])
for line in banner['run_settings']:
if '=' in line:
mode, value = [t.strip() for t in line.split('=')]
mode_status[mode] = value
else:
mode_status = {}
# Call Generate events
self.do_launch('-n %s %s' % (self.run_name, '-f' if self.force else ''),
switch=mode_status)
############################################################################
def do_generate_events(self, line):
"""Main commands: generate events
this function just wraps the do_launch one"""
self.do_launch(line)
############################################################################
def do_treatcards(self, line, amcatnlo=True,mode=''):
"""Advanced commands: this is for creating the correct run_card.inc from the nlo format"""
#check if no 'Auto' are present in the file
self.check_param_card(pjoin(self.me_dir, 'Cards','param_card.dat'))
# propagate the FO_card entry FO_LHE_weight_ratio to the run_card.
# this variable is system only in the run_card
# can not be done in EditCard since this parameter is not written in the
# run_card directly.
if mode in ['LO', 'NLO']:
name = 'fo_lhe_weight_ratio'
FO_card = analyse_card.FOAnalyseCard(pjoin(self.me_dir,'Cards', 'FO_analyse_card.dat'))
if name in FO_card:
self.run_card.set(name, FO_card[name], user=False)
name = 'fo_lhe_postprocessing'
if name in FO_card:
self.run_card.set(name, FO_card[name], user=False)
return super(aMCatNLOCmd,self).do_treatcards(line, amcatnlo)
############################################################################
def set_configuration(self, amcatnlo=True, **opt):
"""assign all configuration variable from file
loop over the different config file if config_file not define """
return super(aMCatNLOCmd,self).set_configuration(amcatnlo=amcatnlo, **opt)
############################################################################
def do_launch(self, line, options={}, argss=[], switch={}):
"""Main commands: launch the full chain
options and args are relevant if the function is called from other
functions, such as generate_events or calculate_xsect
mode gives the list of switch needed for the computation (usefull for banner_run)
"""
if not argss and not options:
self.start_time = time.time()
argss = self.split_arg(line)
# check argument validity and normalise argument
(options, argss) = _launch_parser.parse_args(argss)
options = options.__dict__
self.check_launch(argss, options)
if 'run_name' in options.keys() and options['run_name']:
self.run_name = options['run_name']
# if a dir with the given run_name already exists
# remove it and warn the user
if os.path.isdir(pjoin(self.me_dir, 'Events', self.run_name)):
logger.warning('Removing old run information in \n'+
pjoin(self.me_dir, 'Events', self.run_name))
files.rm(pjoin(self.me_dir, 'Events', self.run_name))
self.results.delete_run(self.run_name)
else:
self.run_name = '' # will be set later
if options['multicore']:
self.cluster_mode = 2
elif options['cluster']:
self.cluster_mode = 1
if not switch:
mode = argss[0]
if mode in ['LO', 'NLO']:
options['parton'] = True
mode = self.ask_run_configuration(mode, options)
else:
mode = self.ask_run_configuration('auto', options, switch)
self.results.add_detail('run_mode', mode)
self.update_status('Starting run', level=None, update_results=True)
if self.options['automatic_html_opening']:
misc.open_file(os.path.join(self.me_dir, 'crossx.html'))
self.options['automatic_html_opening'] = False
if '+' in mode:
mode = mode.split('+')[0]
self.compile(mode, options)
evt_file = self.run(mode, options)
if self.run_card['nevents'] == 0 and not mode in ['LO', 'NLO']:
logger.info('No event file generated: grids have been set-up with a '\
'relative precision of %s' % self.run_card['req_acc'])
return
if not mode in ['LO', 'NLO']:
assert evt_file == pjoin(self.me_dir,'Events', self.run_name, 'events.lhe'), '%s != %s' %(evt_file, pjoin(self.me_dir,'Events', self.run_name, 'events.lhe.gz'))
if self.run_card['systematics_program'] == 'systematics':
self.exec_cmd('systematics %s %s ' % (self.run_name, ' '.join(self.run_card['systematics_arguments'])))
self.exec_cmd('reweight -from_cards', postcmd=False)
self.exec_cmd('decay_events -from_cards', postcmd=False)
evt_file = pjoin(self.me_dir,'Events', self.run_name, 'events.lhe')
if not mode in ['LO', 'NLO', 'noshower', 'noshowerLO'] \
and not options['parton']:
self.run_mcatnlo(evt_file, options)
self.exec_cmd('madanalysis5_hadron --no_default', postcmd=False, printcmd=False)
elif mode == 'noshower':
logger.warning("""You have chosen not to run a parton shower. NLO events without showering are NOT physical.
Please, shower the Les Houches events before using them for physics analyses.""")
self.update_status('', level='all', update_results=True)
if self.run_card['ickkw'] == 3 and \
(mode in ['noshower'] or \
(('PYTHIA8' not in self.run_card['parton_shower'].upper()) and (mode in ['aMC@NLO']))):
logger.warning("""You are running with FxFx merging enabled.
To be able to merge samples of various multiplicities without double counting,
you have to remove some events after showering 'by hand'.
Please read http://amcatnlo.cern.ch/FxFx_merging.htm for more details.""")
self.store_result()
#check if the param_card defines a scan.
if self.param_card_iterator:
param_card_iterator = self.param_card_iterator
self.param_card_iterator = [] #avoid to next generate go trough here
param_card_iterator.store_entry(self.run_name, self.results.current['cross'])
orig_name = self.run_name
#go trough the scal
with misc.TMP_variable(self, 'allow_notification_center', False):
for i,card in enumerate(param_card_iterator):
card.write(pjoin(self.me_dir,'Cards','param_card.dat'))
self.check_param_card(pjoin(self.me_dir,'Cards','param_card.dat'), dependent=True)
if not options['force']:
options['force'] = True
if options['run_name']:
options['run_name'] = '%s_%s' % (orig_name, i+1)
if not argss:
argss = [mode, "-f"]
elif argss[0] == "auto":
argss[0] = mode
self.do_launch("", options=options, argss=argss, switch=switch)
#self.exec_cmd("launch -f ",precmd=True, postcmd=True,errorhandling=False)
param_card_iterator.store_entry(self.run_name, self.results.current['cross'])
#restore original param_card
param_card_iterator.write(pjoin(self.me_dir,'Cards','param_card.dat'))
name = misc.get_scan_name(orig_name, self.run_name)
path = pjoin(self.me_dir, 'Events','scan_%s.txt' % name)
logger.info("write all cross-section results in %s" % path, '$MG:color:BLACK')
param_card_iterator.write_summary(path)
if self.allow_notification_center:
misc.apple_notify('Run %s finished' % os.path.basename(self.me_dir),
'%s: %s +- %s ' % (self.results.current['run_name'],
self.results.current['cross'],
self.results.current['error']))
############################################################################
def do_compile(self, line):
"""Advanced commands: just compile the executables """
argss = self.split_arg(line)
# check argument validity and normalise argument
(options, argss) = _compile_parser.parse_args(argss)
options = options.__dict__
options['reweightonly'] = False
options['nocompile'] = False
self.check_compile(argss, options)
mode = {'FO': 'NLO', 'MC': 'aMC@NLO'}[argss[0]]
self.ask_run_configuration(mode, options)
self.compile(mode, options)
self.update_status('', level='all', update_results=True)
def update_random_seed(self):
"""Update random number seed with the value from the run_card.
If this is 0, update the number according to a fresh one"""
iseed = self.run_card['iseed']
if iseed == 0:
randinit = open(pjoin(self.me_dir, 'SubProcesses', 'randinit'))
iseed = int(randinit.read()[2:]) + 1
randinit.close()
randinit = open(pjoin(self.me_dir, 'SubProcesses', 'randinit'), 'w')
randinit.write('r=%d' % iseed)
randinit.close()
def run(self, mode, options):
"""runs aMC@NLO. Returns the name of the event file created"""
logger.info('Starting run')
if not 'only_generation' in options.keys():
options['only_generation'] = False
# for second step in applgrid mode, do only the event generation step
if mode in ['LO', 'NLO'] and self.run_card['iappl'] == 2 and not options['only_generation']:
options['only_generation'] = True
self.get_characteristics(pjoin(self.me_dir, 'SubProcesses', 'proc_characteristics'))
self.setup_cluster_or_multicore()
self.update_random_seed()
#find and keep track of all the jobs
folder_names = {'LO': ['born_G*'], 'NLO': ['all_G*'],
'aMC@LO': ['GB*'], 'aMC@NLO': ['GF*']}
folder_names['noshower'] = folder_names['aMC@NLO']
folder_names['noshowerLO'] = folder_names['aMC@LO']
p_dirs = [d for d in \
open(pjoin(self.me_dir, 'SubProcesses', 'subproc.mg')).read().split('\n') if d]
#Clean previous results
self.clean_previous_results(options,p_dirs,folder_names[mode])
mcatnlo_status = ['Setting up grids', 'Computing upper envelope', 'Generating events']
if options['reweightonly']:
event_norm=self.run_card['event_norm']
nevents=self.run_card['nevents']
return self.reweight_and_collect_events(options, mode, nevents, event_norm)
if mode in ['LO', 'NLO']:
# this is for fixed order runs
mode_dict = {'NLO': 'all', 'LO': 'born'}
logger.info('Doing fixed order %s' % mode)
req_acc = self.run_card['req_acc_FO']
# Re-distribute the grids for the 2nd step of the applgrid
# running
if self.run_card['iappl'] == 2:
self.applgrid_distribute(options,mode_dict[mode],p_dirs)
# create a list of dictionaries "jobs_to_run" with all the
# jobs that need to be run
integration_step=-1
jobs_to_run,jobs_to_collect,integration_step = self.create_jobs_to_run(options,p_dirs, \
req_acc,mode_dict[mode],integration_step,mode,fixed_order=True)
self.prepare_directories(jobs_to_run,mode)
# loop over the integration steps. After every step, check
# if we have the required accuracy. If this is the case,
# stop running, else do another step.
while True:
integration_step=integration_step+1
self.run_all_jobs(jobs_to_run,integration_step)
self.collect_log_files(jobs_to_run,integration_step)
jobs_to_run,jobs_to_collect=self.collect_the_results(options,req_acc,jobs_to_run, \
jobs_to_collect,integration_step,mode,mode_dict[mode])
if not jobs_to_run:
# there are no more jobs to run (jobs_to_run is empty)
break
# We are done.
self.finalise_run_FO(folder_names[mode],jobs_to_collect)
self.update_status('Run complete', level='parton', update_results=True)
return
elif mode in ['aMC@NLO','aMC@LO','noshower','noshowerLO']:
if self.ninitial == 1:
raise aMCatNLOError('Decay processes can only be run at fixed order.')
mode_dict = {'aMC@NLO': 'all', 'aMC@LO': 'born',\
'noshower': 'all', 'noshowerLO': 'born'}
shower = self.run_card['parton_shower'].upper()
nevents = self.run_card['nevents']
req_acc = self.run_card['req_acc']
if nevents == 0 and req_acc < 0 :
raise aMCatNLOError('Cannot determine the required accuracy from the number '\
'of events, because 0 events requested. Please set '\
'the "req_acc" parameter in the run_card to a value '\
'between 0 and 1')
elif req_acc >1 or req_acc == 0 :
raise aMCatNLOError('Required accuracy ("req_acc" in the run_card) should '\
'be between larger than 0 and smaller than 1, '\
'or set to -1 for automatic determination. Current '\
'value is %f' % req_acc)
# For more than 1M events, set req_acc to 0.001 (except when it was explicitly set in the run_card)
elif req_acc < 0 and nevents > 1000000 :
req_acc=0.001
shower_list = ['HERWIG6', 'HERWIGPP', 'PYTHIA6Q', 'PYTHIA6PT', 'PYTHIA8']
if not shower in shower_list:
raise aMCatNLOError('%s is not a valid parton shower. '\
'Please use one of the following: %s' \
% (shower, ', '.join(shower_list)))
# check that PYTHIA6PT is not used for processes with FSR
if shower == 'PYTHIA6PT' and self.proc_characteristics['has_fsr']:
raise aMCatNLOError('PYTHIA6PT does not support processes with FSR')
if mode in ['aMC@NLO', 'aMC@LO']:
logger.info('Doing %s matched to parton shower' % mode[4:])
elif mode in ['noshower','noshowerLO']:
logger.info('Generating events without running the shower.')
elif options['only_generation']:
logger.info('Generating events starting from existing results')
jobs_to_run,jobs_to_collect,integration_step = self.create_jobs_to_run(options,p_dirs, \
req_acc,mode_dict[mode],1,mode,fixed_order=False)
# Make sure to update all the jobs to be ready for the event generation step
if options['only_generation']:
jobs_to_run,jobs_to_collect=self.collect_the_results(options,req_acc,jobs_to_run, \
jobs_to_collect,1,mode,mode_dict[mode],fixed_order=False)
else:
self.prepare_directories(jobs_to_run,mode,fixed_order=False)
# Main loop over the three MINT generation steps:
for mint_step, status in enumerate(mcatnlo_status):
if options['only_generation'] and mint_step < 2:
continue
self.update_status(status, level='parton')
self.run_all_jobs(jobs_to_run,mint_step,fixed_order=False)
self.collect_log_files(jobs_to_run,mint_step)
if mint_step+1==2 and nevents==0:
self.print_summary(options,2,mode)
return
jobs_to_run,jobs_to_collect=self.collect_the_results(options,req_acc,jobs_to_run, \
jobs_to_collect,mint_step,mode,mode_dict[mode],fixed_order=False)
# Sanity check on the event files. If error the jobs are resubmitted
self.check_event_files(jobs_to_collect)
if self.cluster_mode == 1:
#if cluster run, wait 10 sec so that event files are transferred back
self.update_status(
'Waiting while files are transferred back from the cluster nodes',
level='parton')
time.sleep(10)
event_norm=self.run_card['event_norm']
return self.reweight_and_collect_events(options, mode, nevents, event_norm)
def create_jobs_to_run(self,options,p_dirs,req_acc,run_mode,\
integration_step,mode,fixed_order=True):
"""Creates a list of dictionaries with all the jobs to be run"""
jobs_to_run=[]
if not options['only_generation']:
# Fresh, new run. Check all the P*/channels.txt files
# (created by the 'gensym' executable) to set-up all the
# jobs using the default inputs.
npoints = self.run_card['npoints_FO_grid']
niters = self.run_card['niters_FO_grid']
for p_dir in p_dirs:
try:
with open(pjoin(self.me_dir,'SubProcesses',p_dir,'channels.txt')) as chan_file:
channels=chan_file.readline().split()
except IOError:
logger.warning('No integration channels found for contribution %s' % p_dir)
continue
if fixed_order:
lch=len(channels)
maxchannels=20 # combine up to 20 channels in a single job
if self.run_card['iappl'] != 0: maxchannels=1
njobs=(int(lch/maxchannels)+1 if lch%maxchannels!= 0 \
else int(lch/maxchannels))
for nj in range(1,njobs+1):
job={}
job['p_dir']=p_dir
job['channel']=str(nj)
job['nchans']=(int(lch/njobs)+1 if nj <= lch%njobs else int(lch/njobs))
job['configs']=' '.join(channels[:job['nchans']])
del channels[:job['nchans']]
job['split']=0
if req_acc == -1:
job['accuracy']=0
job['niters']=niters
job['npoints']=npoints
elif req_acc > 0:
job['accuracy']=0.05
job['niters']=6
job['npoints']=-1
else:
raise aMCatNLOError('No consistent "req_acc_FO" set. Use a value '+
'between 0 and 1 or set it equal to -1.')
job['mint_mode']=0
job['run_mode']=run_mode
job['wgt_frac']=1.0
job['wgt_mult']=1.0
jobs_to_run.append(job)
if channels:
raise aMCatNLOError('channels is not empty %s' % channels)
else:
for channel in channels:
job={}
job['p_dir']=p_dir
job['channel']=channel
job['split']=0
job['accuracy']=0.03
job['niters']=12
job['npoints']=-1
job['mint_mode']=0
job['run_mode']=run_mode
job['wgt_frac']=1.0
jobs_to_run.append(job)
jobs_to_collect=copy.copy(jobs_to_run) # These are all jobs
else:
# if options['only_generation'] is true, just read the current jobs from file
try:
with open(pjoin(self.me_dir,"SubProcesses","job_status.pkl"),'rb') as f:
jobs_to_collect=pickle.load(f)
jobs_to_run=copy.copy(jobs_to_collect)
except:
raise aMCatNLOError('Cannot reconstruct saved job status in %s' % \
pjoin(self.me_dir,'SubProcesses','job_status.pkl'))
# Update cross sections and determine which jobs to run next
if fixed_order:
jobs_to_run,jobs_to_collect=self.collect_the_results(options,req_acc,jobs_to_run,
jobs_to_collect,integration_step,mode,run_mode)
# Update the integration_step to make sure that nothing will be overwritten
integration_step=1
for job in jobs_to_run:
while os.path.exists(pjoin(job['dirname'],'res_%s.dat' % integration_step)):
integration_step=integration_step+1
integration_step=integration_step-1
else:
self.append_the_results(jobs_to_collect,integration_step)
return jobs_to_run,jobs_to_collect,integration_step
def prepare_directories(self,jobs_to_run,mode,fixed_order=True):
"""Set-up the G* directories for running"""
name_suffix={'born' :'B' , 'all':'F'}
for job in jobs_to_run:
if job['split'] == 0:
if fixed_order :
dirname=pjoin(self.me_dir,'SubProcesses',job['p_dir'],
job['run_mode']+'_G'+job['channel'])
else:
dirname=pjoin(self.me_dir,'SubProcesses',job['p_dir'],
'G'+name_suffix[job['run_mode']]+job['channel'])
else:
if fixed_order :
dirname=pjoin(self.me_dir,'SubProcesses',job['p_dir'],
job['run_mode']+'_G'+job['channel']+'_'+str(job['split']))
else:
dirname=pjoin(self.me_dir,'SubProcesses',job['p_dir'],
'G'+name_suffix[job['run_mode']]+job['channel']+'_'+str(job['split']))
job['dirname']=dirname
if not os.path.isdir(dirname):
os.makedirs(dirname)
self.write_input_file(job,fixed_order)
# link or copy the grids from the base directory to the split directory:
if not fixed_order:
if job['split'] != 0:
for f in ['grid.MC_integer','mint_grids','res_1']:
if not os.path.isfile(pjoin(job['dirname'],f)):
files.ln(pjoin(job['dirname'].rsplit("_",1)[0],f),job['dirname'])
else:
if job['split'] != 0:
for f in ['grid.MC_integer','mint_grids']:
files.cp(pjoin(job['dirname'].rsplit("_",1)[0],f),job['dirname'])
def write_input_file(self,job,fixed_order):
"""write the input file for the madevent_mint* executable in the appropriate directory"""
if fixed_order:
content= \
"""NPOINTS = %(npoints)s
NITERATIONS = %(niters)s
ACCURACY = %(accuracy)s
ADAPT_GRID = 2
MULTICHANNEL = 1
SUM_HELICITY = 1
NCHANS = %(nchans)s
CHANNEL = %(configs)s
SPLIT = %(split)s
WGT_MULT= %(wgt_mult)s
RUN_MODE = %(run_mode)s
RESTART = %(mint_mode)s
""" \
% job
else:
content = \
"""-1 12 ! points, iterations
%(accuracy)s ! desired fractional accuracy
1 -0.1 ! alpha, beta for Gsoft
1 -0.1 ! alpha, beta for Gazi
1 ! Suppress amplitude (0 no, 1 yes)?
1 ! Exact helicity sum (0 yes, n = number/event)?
%(channel)s ! Enter Configuration Number:
%(mint_mode)s ! MINT imode: 0 to set-up grids, 1 to perform integral, 2 generate events
1 1 1 ! if imode is 1: Folding parameters for xi_i, phi_i and y_ij
%(run_mode)s ! all, born, real, virt
""" \
% job
with open(pjoin(job['dirname'], 'input_app.txt'), 'w') as input_file:
input_file.write(content)
def run_all_jobs(self,jobs_to_run,integration_step,fixed_order=True):
"""Loops over the jobs_to_run and executes them using the function 'run_exe'"""
if fixed_order:
if integration_step == 0:
self.update_status('Setting up grids', level=None)
else:
self.update_status('Refining results, step %i' % integration_step, level=None)
self.ijob = 0
name_suffix={'born' :'B', 'all':'F'}
if fixed_order:
run_type="Fixed order integration step %s" % integration_step
else:
run_type="MINT step %s" % integration_step
self.njobs=len(jobs_to_run)
for job in jobs_to_run:
executable='ajob1'
if fixed_order:
arguments=[job['channel'],job['run_mode'], \
str(job['split']),str(integration_step)]
else:
arguments=[job['channel'],name_suffix[job['run_mode']], \
str(job['split']),str(integration_step)]
self.run_exe(executable,arguments,run_type,
cwd=pjoin(self.me_dir,'SubProcesses',job['p_dir']))
if self.cluster_mode == 2:
time.sleep(1) # security to allow all jobs to be launched
self.wait_for_complete(run_type)
def collect_the_results(self,options,req_acc,jobs_to_run,jobs_to_collect,\
integration_step,mode,run_mode,fixed_order=True):
"""Collect the results, make HTML pages, print the summary and
determine if there are more jobs to run. Returns the list
of the jobs that still need to be run, as well as the
complete list of jobs that need to be collected to get the
final answer.
"""
# Get the results of the current integration/MINT step
self.append_the_results(jobs_to_run,integration_step)
self.cross_sect_dict = self.write_res_txt_file(jobs_to_collect,integration_step)
# Update HTML pages
if fixed_order:
cross, error = sum_html.make_all_html_results(self, jobs=jobs_to_collect)
else:
name_suffix={'born' :'B' , 'all':'F'}
cross, error = sum_html.make_all_html_results(self, folder_names=['G%s*' % name_suffix[run_mode]])
self.results.add_detail('cross', cross)
self.results.add_detail('error', error)
# Combine grids from split fixed order jobs
if fixed_order:
jobs_to_run=self.combine_split_order_run(jobs_to_run)
# Set-up jobs for the next iteration/MINT step
jobs_to_run_new=self.update_jobs_to_run(req_acc,integration_step,jobs_to_run,fixed_order)
# IF THERE ARE NO MORE JOBS, WE ARE DONE!!!
if fixed_order:
# Write the jobs_to_collect directory to file so that we
# can restart them later (with only-generation option)
with open(pjoin(self.me_dir,"SubProcesses","job_status.pkl"),'wb') as f:
pickle.dump(jobs_to_collect,f)
# Print summary
if (not jobs_to_run_new) and fixed_order:
# print final summary of results (for fixed order)
scale_pdf_info=self.collect_scale_pdf_info(options,jobs_to_collect)
self.print_summary(options,integration_step,mode,scale_pdf_info,done=True)
return jobs_to_run_new,jobs_to_collect
elif jobs_to_run_new:
# print intermediate summary of results
scale_pdf_info=[]
self.print_summary(options,integration_step,mode,scale_pdf_info,done=False)
else:
# When we are done for (N)LO+PS runs, do not print
# anything yet. This will be done after the reweighting
# and collection of the events
scale_pdf_info=[]
# Prepare for the next integration/MINT step
if (not fixed_order) and integration_step+1 == 2 :
# Write the jobs_to_collect directory to file so that we
# can restart them later (with only-generation option)
with open(pjoin(self.me_dir,"SubProcesses","job_status.pkl"),'wb') as f:
pickle.dump(jobs_to_collect,f)
# next step is event generation (mint_step 2)
jobs_to_run_new,jobs_to_collect_new= \
self.check_the_need_to_split(jobs_to_run_new,jobs_to_collect)
self.prepare_directories(jobs_to_run_new,mode,fixed_order)
self.write_nevents_unweighted_file(jobs_to_collect_new,jobs_to_collect)
self.write_nevts_files(jobs_to_run_new)
else:
if fixed_order and self.run_card['iappl'] == 0 \
and self.run_card['req_acc_FO'] > 0:
jobs_to_run_new,jobs_to_collect= \
self.split_jobs_fixed_order(jobs_to_run_new,jobs_to_collect)
self.prepare_directories(jobs_to_run_new,mode,fixed_order)
jobs_to_collect_new=jobs_to_collect
return jobs_to_run_new,jobs_to_collect_new
def write_nevents_unweighted_file(self,jobs,jobs0events):
"""writes the nevents_unweighted file in the SubProcesses directory.
We also need to write the jobs that will generate 0 events,
because that makes sure that the cross section from those channels
is taken into account in the event weights (by collect_events.f).
"""
content=[]
for job in jobs:
path=pjoin(job['dirname'].split('/')[-2],job['dirname'].split('/')[-1])
lhefile=pjoin(path,'events.lhe')
content.append(' %s %d %9e %9e' % \
(lhefile.ljust(40),job['nevents'],job['resultABS']*job['wgt_frac'],job['wgt_frac']))
for job in jobs0events:
if job['nevents']==0:
path=pjoin(job['dirname'].split('/')[-2],job['dirname'].split('/')[-1])
lhefile=pjoin(path,'events.lhe')
content.append(' %s %d %9e %9e' % \
(lhefile.ljust(40),job['nevents'],job['resultABS'],1.))
with open(pjoin(self.me_dir,'SubProcesses',"nevents_unweighted"),'w') as f:
f.write('\n'.join(content)+'\n')
def write_nevts_files(self,jobs):
"""write the nevts files in the SubProcesses/P*/G*/ directories"""
for job in jobs:
with open(pjoin(job['dirname'],'nevts'),'w') as f:
f.write('%i\n' % job['nevents'])
def combine_split_order_run(self,jobs_to_run):
"""Combines jobs and grids from split jobs that have been run"""
# combine the jobs that need to be combined in job
# groups. Simply combine the ones that have the same p_dir and
# same channel.
jobgroups_to_combine=[]
jobs_to_run_new=[]
for job in jobs_to_run:
if job['split'] == 0:
job['combined']=1
jobs_to_run_new.append(job) # this jobs wasn't split
elif job['split'] == 1:
jobgroups_to_combine.append(filter(lambda j: j['p_dir'] == job['p_dir'] and \
j['channel'] == job['channel'], jobs_to_run))
else:
continue
for job_group in jobgroups_to_combine:
# Combine the grids (mint-grids & MC-integer grids) first
self.combine_split_order_grids(job_group)
jobs_to_run_new.append(self.combine_split_order_jobs(job_group))
return jobs_to_run_new
def combine_split_order_jobs(self,job_group):
"""combine the jobs in job_group and return a single summed job"""
# first copy one of the jobs in 'jobs'
sum_job=copy.copy(job_group[0])
# update the information to have a 'non-split' job:
sum_job['dirname']=pjoin(sum_job['dirname'].rsplit('_',1)[0])
sum_job['split']=0
sum_job['wgt_mult']=1.0
sum_job['combined']=len(job_group)
# information to be summed:
keys=['niters_done','npoints_done','niters','npoints',\
'result','resultABS','time_spend']
keys2=['error','errorABS']
# information to be summed in quadrature:
for key in keys2:
sum_job[key]=math.pow(sum_job[key],2)
# Loop over the jobs and sum the information
for i,job in enumerate(job_group):
if i==0 : continue # skip the first
for key in keys:
sum_job[key]+=job[key]
for key in keys2:
sum_job[key]+=math.pow(job[key],2)
for key in keys2:
sum_job[key]=math.sqrt(sum_job[key])
sum_job['err_percABS'] = sum_job['errorABS']/sum_job['resultABS']*100.
sum_job['err_perc'] = sum_job['error']/sum_job['result']*100.
sum_job['niters']=int(sum_job['niters_done']/len(job_group))
sum_job['niters_done']=int(sum_job['niters_done']/len(job_group))
return sum_job
def combine_split_order_grids(self,job_group):
"""Combines the mint_grids and MC-integer grids from the split order
jobs (fixed order only).
"""
files_mint_grids=[]
files_MC_integer=[]
location=None
for job in job_group:
files_mint_grids.append(open(pjoin(job['dirname'],'mint_grids'),'r+'))
files_MC_integer.append(open(pjoin(job['dirname'],'grid.MC_integer'),'r+'))
if not location:
location=pjoin(job['dirname'].rsplit('_',1)[0])
else:
if location != pjoin(job['dirname'].rsplit('_',1)[0]) :
raise aMCatNLOError('Not all jobs have the same location. '\
+'Cannot combine them.')
# Needed to average the grids (both xgrids, ave_virt and
# MC_integer grids), but sum the cross section info. The
# latter is only the only line that contains integers.
for j,fs in enumerate([files_mint_grids,files_MC_integer]):
linesoffiles=[f.readlines() for f in fs]
to_write=[]
for rowgrp in zip(*linesoffiles):
try:
# check that last element on the line is an
# integer (will raise ValueError if not the
# case). If integer, this is the line that
# contains information that needs to be
# summed. All other lines can be averaged.
is_integer = [[int(row.strip().split()[-1])] for row in rowgrp]
floatsbyfile = [[float(a) for a in row.strip().split()] for row in rowgrp]
floatgrps = zip(*floatsbyfile)
special=[]
for i,floatgrp in enumerate(floatgrps):
if i==0: # sum X-sec
special.append(sum(floatgrp))
elif i==1: # sum unc in quadrature
special.append(math.sqrt(sum([err**2 for err in floatgrp])))
elif i==2: # average number of PS per iteration
special.append(int(sum(floatgrp)/len(floatgrp)))
elif i==3: # sum the number of iterations
special.append(int(sum(floatgrp)))
elif i==4: # average the nhits_in_grids
special.append(int(sum(floatgrp)/len(floatgrp)))
else:
raise aMCatNLOError('"mint_grids" files not in correct format. '+\
'Cannot combine them.')
to_write.append(" ".join(str(s) for s in special) + "\n")
except ValueError:
# just average all
floatsbyfile = [[float(a) for a in row.strip().split()] for row in rowgrp]
floatgrps = zip(*floatsbyfile)
averages = [sum(floatgrp)/len(floatgrp) for floatgrp in floatgrps]
to_write.append(" ".join(str(a) for a in averages) + "\n")
# close the files
for f in fs:
f.close
# write the data over the master location
if j==0:
with open(pjoin(location,'mint_grids'),'w') as f:
f.writelines(to_write)
elif j==1:
with open(pjoin(location,'grid.MC_integer'),'w') as f:
f.writelines(to_write)
def split_jobs_fixed_order(self,jobs_to_run,jobs_to_collect):
"""Looks in the jobs_to_run to see if there is the need to split the
jobs, depending on the expected time they take. Updates
jobs_to_run and jobs_to_collect to replace the split-job by
its splits.
"""
# determine the number jobs we should have (this is per p_dir)
if self.options['run_mode'] ==2:
nb_submit = int(self.options['nb_core'])
elif self.options['run_mode'] ==1:
nb_submit = int(self.options['cluster_size'])
else:
nb_submit =1
# total expected aggregated running time
time_expected=0
for job in jobs_to_run:
time_expected+=job['time_spend']*(job['niters']*job['npoints'])/ \
(job['niters_done']*job['npoints_done'])
# this means that we must expect the following per job (in
# ideal conditions)
time_per_job=time_expected/(nb_submit*(1+len(jobs_to_run)/2))
jobs_to_run_new=[]
jobs_to_collect_new=copy.copy(jobs_to_collect)
for job in jobs_to_run:
# remove current job from jobs_to_collect. Make sure
# to remove all the split ones in case the original
# job had been a split one (before it was re-combined)
for j in filter(lambda j: j['p_dir'] == job['p_dir'] and \
j['channel'] == job['channel'], jobs_to_collect_new):
jobs_to_collect_new.remove(j)
time_expected=job['time_spend']*(job['niters']*job['npoints'])/ \
(job['niters_done']*job['npoints_done'])
# if the time expected for this job is (much) larger than
# the time spend in the previous iteration, and larger
# than the expected time per job, split it
if time_expected > max(2*job['time_spend']/job['combined'],time_per_job):
# determine the number of splits needed
nsplit=min(max(int(time_expected/max(2*job['time_spend']/job['combined'],time_per_job)),2),nb_submit)
for i in range(1,nsplit+1):
job_new=copy.copy(job)
job_new['split']=i
job_new['wgt_mult']=1./float(nsplit)
job_new['dirname']=job['dirname']+'_%i' % job_new['split']
job_new['accuracy']=min(job['accuracy']*math.sqrt(float(nsplit)),0.1)
if nsplit >= job['niters']:
job_new['npoints']=int(job['npoints']*job['niters']/nsplit)
job_new['niters']=1
else:
job_new['npoints']=int(job['npoints']/nsplit)
jobs_to_collect_new.append(job_new)
jobs_to_run_new.append(job_new)
else:
jobs_to_collect_new.append(job)
jobs_to_run_new.append(job)
return jobs_to_run_new,jobs_to_collect_new
def check_the_need_to_split(self,jobs_to_run,jobs_to_collect):
"""Looks in the jobs_to_run to see if there is the need to split the
event generation step. Updates jobs_to_run and
jobs_to_collect to replace the split-job by its
splits. Also removes jobs that do not need any events.
"""
nevt_job=self.run_card['nevt_job']
if nevt_job > 0:
jobs_to_collect_new=copy.copy(jobs_to_collect)
for job in jobs_to_run:
nevents=job['nevents']
if nevents == 0:
jobs_to_collect_new.remove(job)
elif nevents > nevt_job:
jobs_to_collect_new.remove(job)
if nevents % nevt_job != 0 :
nsplit=int(nevents/nevt_job)+1
else:
nsplit=int(nevents/nevt_job)
for i in range(1,nsplit+1):
job_new=copy.copy(job)
left_over=nevents % nsplit
if i <= left_over:
job_new['nevents']=int(nevents/nsplit)+1
job_new['wgt_frac']=float(job_new['nevents'])/float(nevents)
else:
job_new['nevents']=int(nevents/nsplit)
job_new['wgt_frac']=float(job_new['nevents'])/float(nevents)
job_new['split']=i
job_new['dirname']=job['dirname']+'_%i' % job_new['split']
jobs_to_collect_new.append(job_new)
jobs_to_run_new=copy.copy(jobs_to_collect_new)
else:
jobs_to_run_new=copy.copy(jobs_to_collect)
for job in jobs_to_collect:
if job['nevents'] == 0:
jobs_to_run_new.remove(job)
jobs_to_collect_new=copy.copy(jobs_to_run_new)
return jobs_to_run_new,jobs_to_collect_new
def update_jobs_to_run(self,req_acc,step,jobs,fixed_order=True):
"""
For (N)LO+PS: determines the number of events and/or the required
accuracy per job.
For fixed order: determines which jobs need higher precision and
returns those with the newly requested precision.
"""
err=self.cross_sect_dict['errt']
tot=self.cross_sect_dict['xsect']
errABS=self.cross_sect_dict['erra']
totABS=self.cross_sect_dict['xseca']
jobs_new=[]
if fixed_order:
if req_acc == -1:
if step+1 == 1:
npoints = self.run_card['npoints_FO']
niters = self.run_card['niters_FO']
for job in jobs:
job['mint_mode']=-1
job['niters']=niters
job['npoints']=npoints
jobs_new.append(job)
elif step+1 == 2:
pass
elif step+1 > 2:
raise aMCatNLOError('Cannot determine number of iterations and PS points '+
'for integration step %i' % step )
elif ( req_acc > 0 and err/abs(tot) > req_acc*1.2 ) or step <= 0:
req_accABS=req_acc*abs(tot)/totABS # overal relative required accuracy on ABS Xsec.
for job in jobs:
job['mint_mode']=-1
# Determine relative required accuracy on the ABS for this job
job['accuracy']=req_accABS*math.sqrt(totABS/job['resultABS'])
# If already accurate enough, skip the job (except when doing the first
# step for the iappl=2 run: we need to fill all the applgrid grids!)
if (job['accuracy'] > job['errorABS']/job['resultABS'] and step != 0) \
and not (step==-1 and self.run_card['iappl'] == 2):
continue
# Update the number of PS points based on errorABS, ncall and accuracy
itmax_fl=job['niters_done']*math.pow(job['errorABS']/
(job['accuracy']*job['resultABS']),2)
if itmax_fl <= 4.0 :
job['niters']=max(int(round(itmax_fl)),2)
job['npoints']=job['npoints_done']*2
elif itmax_fl > 4.0 and itmax_fl <= 16.0 :
job['niters']=4
job['npoints']=int(round(job['npoints_done']*itmax_fl/4.0))*2
else:
if itmax_fl > 100.0 : itmax_fl=50.0
job['niters']=int(round(math.sqrt(itmax_fl)))
job['npoints']=int(round(job['npoints_done']*itmax_fl/
round(math.sqrt(itmax_fl))))*2
# Add the job to the list of jobs that need to be run
jobs_new.append(job)
return jobs_new
elif step+1 <= 2:
nevents=self.run_card['nevents']
# Total required accuracy for the upper bounding envelope
if req_acc<0:
req_acc2_inv=nevents
else:
req_acc2_inv=1/(req_acc*req_acc)
if step+1 == 1 or step+1 == 2 :
# determine the req. accuracy for each of the jobs for Mint-step = 1
for job in jobs:
accuracy=min(math.sqrt(totABS/(req_acc2_inv*job['resultABS'])),0.2)
job['accuracy']=accuracy
if step+1 == 2:
# Randomly (based on the relative ABS Xsec of the job) determine the
# number of events each job needs to generate for MINT-step = 2.
r=self.get_randinit_seed()
random.seed(r)
totevts=nevents
for job in jobs:
job['nevents'] = 0
while totevts :
target = random.random() * totABS
crosssum = 0.
i = 0
while i= 0 :
with open(pjoin(job['dirname'],'res_%s.dat' % integration_step)) as res_file:
results=res_file.readline().split()
else:
# should only be here when doing fixed order with the 'only_generation'
# option equal to True. Take the results from the final run done.
with open(pjoin(job['dirname'],'res.dat')) as res_file:
results=res_file.readline().split()
except IOError:
if not error_found:
error_found=True
error_log=[]
error_log.append(pjoin(job['dirname'],'log.txt'))
continue
job['resultABS']=float(results[0])
job['errorABS']=float(results[1])
job['result']=float(results[2])
job['error']=float(results[3])
job['niters_done']=int(results[4])
job['npoints_done']=int(results[5])
job['time_spend']=float(results[6])
job['err_percABS'] = job['errorABS']/job['resultABS']*100.
job['err_perc'] = job['error']/job['result']*100.
if error_found:
raise aMCatNLOError('An error occurred during the collection of results.\n' +
'Please check the .log files inside the directories which failed:\n' +
'\n'.join(error_log)+'\n')
def write_res_txt_file(self,jobs,integration_step):
"""writes the res.txt files in the SubProcess dir"""
jobs.sort(key = lambda job: -job['errorABS'])
content=[]
content.append('\n\nCross section per integration channel:')
for job in jobs:
content.append('%(p_dir)20s %(channel)15s %(result)10.8e %(error)6.4e %(err_perc)6.4f%% ' % job)
content.append('\n\nABS cross section per integration channel:')
for job in jobs:
content.append('%(p_dir)20s %(channel)15s %(resultABS)10.8e %(errorABS)6.4e %(err_percABS)6.4f%% ' % job)
totABS=0
errABS=0
tot=0
err=0
for job in jobs:
totABS+= job['resultABS']*job['wgt_frac']
errABS+= math.pow(job['errorABS'],2)*job['wgt_frac']
tot+= job['result']*job['wgt_frac']
err+= math.pow(job['error'],2)*job['wgt_frac']
if jobs:
content.append('\nTotal ABS and \nTotal: \n %10.8e +- %6.4e (%6.4e%%)\n %10.8e +- %6.4e (%6.4e%%) \n' %\
(totABS, math.sqrt(errABS), math.sqrt(errABS)/totABS *100.,\
tot, math.sqrt(err), math.sqrt(err)/tot *100.))
with open(pjoin(self.me_dir,'SubProcesses','res_%s.txt' % integration_step),'w') as res_file:
res_file.write('\n'.join(content))
randinit=self.get_randinit_seed()
return {'xsect':tot,'xseca':totABS,'errt':math.sqrt(err),\
'erra':math.sqrt(errABS),'randinit':randinit}
def collect_scale_pdf_info(self,options,jobs):
"""read the scale_pdf_dependence.dat files and collects there results"""
scale_pdf_info=[]
if any(self.run_card['reweight_scale']) or any(self.run_card['reweight_PDF']) or \
len(self.run_card['dynamical_scale_choice']) > 1 or len(self.run_card['lhaid']) > 1:
evt_files=[]
evt_wghts=[]
for job in jobs:
evt_files.append(pjoin(job['dirname'],'scale_pdf_dependence.dat'))
evt_wghts.append(job['wgt_frac'])
scale_pdf_info = self.pdf_scale_from_reweighting(evt_files,evt_wghts)
return scale_pdf_info
def combine_plots_FO(self,folder_name,jobs):
"""combines the plots and puts then in the Events/run* directory"""
devnull = open(os.devnull, 'w')
if self.analyse_card['fo_analysis_format'].lower() == 'topdrawer':
misc.call(['./combine_plots_FO.sh'] + folder_name, \
stdout=devnull,
cwd=pjoin(self.me_dir, 'SubProcesses'))
files.cp(pjoin(self.me_dir, 'SubProcesses', 'MADatNLO.top'),
pjoin(self.me_dir, 'Events', self.run_name))
logger.info('The results of this run and the TopDrawer file with the plots' + \
' have been saved in %s' % pjoin(self.me_dir, 'Events', self.run_name))
elif self.analyse_card['fo_analysis_format'].lower() == 'hwu':
out=pjoin(self.me_dir,'Events',self.run_name,'MADatNLO')
self.combine_plots_HwU(jobs,out)
try:
misc.call(['gnuplot','MADatNLO.gnuplot'],\
stdout=devnull,stderr=devnull,\
cwd=pjoin(self.me_dir, 'Events', self.run_name))
except Exception:
pass
logger.info('The results of this run and the HwU and GnuPlot files with the plots' + \
' have been saved in %s' % pjoin(self.me_dir, 'Events', self.run_name))
elif self.analyse_card['fo_analysis_format'].lower() == 'root':
misc.call(['./combine_root.sh'] + folder_name, \
stdout=devnull,
cwd=pjoin(self.me_dir, 'SubProcesses'))
files.cp(pjoin(self.me_dir, 'SubProcesses', 'MADatNLO.root'),
pjoin(self.me_dir, 'Events', self.run_name))
logger.info('The results of this run and the ROOT file with the plots' + \
' have been saved in %s' % pjoin(self.me_dir, 'Events', self.run_name))
elif self.analyse_card['fo_analysis_format'].lower() == 'lhe':
self.combine_FO_lhe(jobs)
logger.info('The results of this run and the LHE File (to be used for plotting only)' + \
' have been saved in %s' % pjoin(self.me_dir, 'Events', self.run_name))
else:
logger.info('The results of this run' + \
' have been saved in %s' % pjoin(self.me_dir, 'Events', self.run_name))
def combine_FO_lhe(self,jobs):
"""combine the various lhe file generated in each directory.
They are two steps:
1) banner
2) reweight each sample by the factor written at the end of each file
3) concatenate each of the new files (gzip those).
"""
logger.info('Combining lhe events for plotting analysis')
start = time.time()
self.run_card['fo_lhe_postprocessing'] = [i.lower() for i in self.run_card['fo_lhe_postprocessing']]
output = pjoin(self.me_dir, 'Events', self.run_name, 'events.lhe.gz')
if os.path.exists(output):
os.remove(output)
# 1. write the banner
text = open(pjoin(jobs[0]['dirname'],'header.txt'),'r').read()
i1, i2 = text.find(''),text.find('')
self.banner['initrwgt'] = text[10+i1:i2]
#
#
# 2212 2212 6.500000e+03 6.500000e+03 0 0 247000 247000 -4 1
# 8.430000e+02 2.132160e+00 8.430000e+02 1
# please cite 1405.0301
#
cross = sum(j['result'] for j in jobs)
error = math.sqrt(sum(j['error'] for j in jobs))
self.banner['init'] = "0 0 0e0 0e0 0 0 0 0 -4 1\n %s %s %s 1" % (cross, error, cross)
self.banner.write(output[:-3], close_tag=False)
misc.gzip(output[:-3])
fsock = lhe_parser.EventFile(output,'a')
if 'nogrouping' in self.run_card['fo_lhe_postprocessing']:
fsock.eventgroup = False
else:
fsock.eventgroup = True
if 'norandom' in self.run_card['fo_lhe_postprocessing']:
for job in jobs:
dirname = job['dirname']
#read last line
lastline = misc.BackRead(pjoin(dirname,'events.lhe')).readline()
nb_event, sumwgt, cross = [float(i) for i in lastline.split()]
# get normalisation ratio
ratio = cross/sumwgt
lhe = lhe_parser.EventFile(pjoin(dirname,'events.lhe'))
lhe.eventgroup = True # read the events by eventgroup
for eventsgroup in lhe:
neweventsgroup = []
for i,event in enumerate(eventsgroup):
event.rescale_weights(ratio)
if i>0 and 'noidentification' not in self.run_card['fo_lhe_postprocessing'] \
and event == neweventsgroup[-1]:
neweventsgroup[-1].wgt += event.wgt
for key in event.reweight_data:
neweventsgroup[-1].reweight_data[key] += event.reweight_data[key]
else:
neweventsgroup.append(event)
fsock.write_events(neweventsgroup)
lhe.close()
os.remove(pjoin(dirname,'events.lhe'))
else:
lhe = []
lenlhe = []
misc.sprint('need to combine %s event file' % len(jobs))
globallhe = lhe_parser.MultiEventFile()
globallhe.eventgroup = True
for job in jobs:
dirname = job['dirname']
lastline = misc.BackRead(pjoin(dirname,'events.lhe')).readline()
nb_event, sumwgt, cross = [float(i) for i in lastline.split()]
lastlhe = globallhe.add(pjoin(dirname,'events.lhe'),cross, 0, cross,
nb_event=int(nb_event), scale=cross/sumwgt)
for eventsgroup in globallhe:
neweventsgroup = []
for i,event in enumerate(eventsgroup):
event.rescale_weights(event.sample_scale)
if i>0 and 'noidentification' not in self.run_card['fo_lhe_postprocessing'] \
and event == neweventsgroup[-1]:
neweventsgroup[-1].wgt += event.wgt
for key in event.reweight_data:
neweventsgroup[-1].reweight_data[key] += event.reweight_data[key]
else:
neweventsgroup.append(event)
fsock.write_events(neweventsgroup)
globallhe.close()
fsock.write('\n')
fsock.close()
misc.sprint('combining lhe file done in ', time.time()-start)
for job in jobs:
dirname = job['dirname']
os.remove(pjoin(dirname,'events.lhe'))
misc.sprint('combining lhe file done in ', time.time()-start)
def combine_plots_HwU(self,jobs,out,normalisation=None):
"""Sums all the plots in the HwU format."""
logger.debug('Combining HwU plots.')
command = []
command.append(pjoin(self.me_dir, 'bin', 'internal','histograms.py'))
for job in jobs:
if job['dirname'].endswith('.HwU'):
command.append(job['dirname'])
else:
command.append(pjoin(job['dirname'],'MADatNLO.HwU'))
command.append("--out="+out)
command.append("--gnuplot")
command.append("--band=[]")
command.append("--lhapdf-config="+self.options['lhapdf'])
if normalisation:
command.append("--multiply="+(','.join([str(n) for n in normalisation])))
command.append("--sum")
command.append("--keep_all_weights")
command.append("--no_open")
p = misc.Popen(command, stdout = subprocess.PIPE, stderr = subprocess.STDOUT, cwd=self.me_dir)
while p.poll() is None:
line = p.stdout.readline()
if any(t in line for t in ['INFO:','WARNING:','CRITICAL:','ERROR:','KEEP:']):
print line[:-1]
elif __debug__ and line:
logger.debug(line[:-1])
def applgrid_combine(self,cross,error,jobs):
"""Combines the APPLgrids in all the SubProcess/P*/all_G*/ directories"""
logger.debug('Combining APPLgrids \n')
applcomb=pjoin(self.options['applgrid'].rstrip('applgrid-config'),
'applgrid-combine')
all_jobs=[]
for job in jobs:
all_jobs.append(job['dirname'])
ngrids=len(all_jobs)
nobs =len([name for name in os.listdir(all_jobs[0]) if name.endswith("_out.root")])
for obs in range(0,nobs):
gdir = [pjoin(job,"grid_obs_"+str(obs)+"_out.root") for job in all_jobs]
# combine APPLgrids from different channels for observable 'obs'
if self.run_card["iappl"] == 1:
misc.call([applcomb,'-o', pjoin(self.me_dir,"Events",self.run_name,
"aMCfast_obs_"+str(obs)+"_starting_grid.root"), '--optimise']+ gdir)
elif self.run_card["iappl"] == 2:
unc2_inv=pow(cross/error,2)
unc2_inv_ngrids=pow(cross/error,2)*ngrids
misc.call([applcomb,'-o', pjoin(self.me_dir,"Events",
self.run_name,"aMCfast_obs_"+str(obs)+".root"),'-s',
str(unc2_inv),'--weight',str(unc2_inv)]+ gdir)
for job in all_jobs:
os.remove(pjoin(job,"grid_obs_"+str(obs)+"_in.root"))
else:
raise aMCatNLOError('iappl parameter can only be 0, 1 or 2')
# after combining, delete the original grids
for ggdir in gdir:
os.remove(ggdir)
def applgrid_distribute(self,options,mode,p_dirs):
"""Distributes the APPLgrids ready to be filled by a second run of the code"""
# if no appl_start_grid argument given, guess it from the time stamps
# of the starting grid files
if not('appl_start_grid' in options.keys() and options['appl_start_grid']):
gfiles = misc.glob(pjoin('*', 'aMCfast_obs_0_starting_grid.root'),
pjoin(self.me_dir,'Events'))
time_stamps={}
for root_file in gfiles:
time_stamps[root_file]=os.path.getmtime(root_file)
options['appl_start_grid']= \
max(time_stamps.iterkeys(), key=(lambda key:
time_stamps[key])).split('/')[-2]
logger.info('No --appl_start_grid option given. '+\
'Guessing that start grid from run "%s" should be used.' \
% options['appl_start_grid'])
if 'appl_start_grid' in options.keys() and options['appl_start_grid']:
self.appl_start_grid = options['appl_start_grid']
start_grid_dir=pjoin(self.me_dir, 'Events', self.appl_start_grid)
# check that this dir exists and at least one grid file is there
if not os.path.exists(pjoin(start_grid_dir,
'aMCfast_obs_0_starting_grid.root')):
raise self.InvalidCmd('APPLgrid file not found: %s' % \
pjoin(start_grid_dir,'aMCfast_obs_0_starting_grid.root'))
else:
all_grids=[pjoin(start_grid_dir,name) for name in os.listdir( \
start_grid_dir) if name.endswith("_starting_grid.root")]
nobs =len(all_grids)
gstring=" ".join(all_grids)
if not hasattr(self, 'appl_start_grid') or not self.appl_start_grid:
raise self.InvalidCmd('No APPLgrid name currently defined.'+
'Please provide this information.')
#copy the grid to all relevant directories
for pdir in p_dirs:
g_dirs = [file for file in os.listdir(pjoin(self.me_dir,
"SubProcesses",pdir)) if file.startswith(mode+'_G') and
os.path.isdir(pjoin(self.me_dir,"SubProcesses",pdir, file))]
for g_dir in g_dirs:
for grid in all_grids:
obs=grid.split('_')[-3]
files.cp(grid,pjoin(self.me_dir,"SubProcesses",pdir,g_dir,
'grid_obs_'+obs+'_in.root'))
def collect_log_files(self, jobs, integration_step):
"""collect the log files and put them in a single, html-friendly file
inside the Events/run_.../ directory"""
log_file = pjoin(self.me_dir, 'Events', self.run_name,
'alllogs_%d.html' % integration_step)
outfile = open(log_file, 'w')
content = ''
content += '\n'
for job in jobs:
# put an anchor
log=pjoin(job['dirname'],'log_MINT%s.txt' % integration_step)
content += '\n' % (os.path.dirname(log).replace(
pjoin(self.me_dir,'SubProcesses'),''))
# and put some nice header
content += '\n'
content += ' LOG file for integration channel %s, %s ' % \
(os.path.dirname(log).replace(pjoin(self.me_dir,
'SubProcesses'), ''),
integration_step)
content += '\n'
#then just flush the content of the small log inside the big log
#the PRE tag prints everything verbatim
with open(log) as l:
content += '
\n' + l.read() + '\n
'
content +=' \n'
outfile.write(content)
content=''
outfile.write('\n\n')
outfile.close()
def finalise_run_FO(self,folder_name,jobs):
"""Combine the plots and put the res*.txt files in the Events/run.../ folder."""
# Copy the res_*.txt files to the Events/run* folder
res_files = misc.glob('res_*.txt', pjoin(self.me_dir, 'SubProcesses'))
for res_file in res_files:
files.mv(res_file,pjoin(self.me_dir, 'Events', self.run_name))
# Collect the plots and put them in the Events/run* folder
self.combine_plots_FO(folder_name,jobs)
# If doing the applgrid-stuff, also combine those grids
# and put those in the Events/run* folder
if self.run_card['iappl'] != 0:
cross=self.cross_sect_dict['xsect']
error=self.cross_sect_dict['errt']
self.applgrid_combine(cross,error,jobs)
def setup_cluster_or_multicore(self):
"""setup the number of cores for multicore, and the cluster-type for cluster runs"""
if self.cluster_mode == 1:
cluster_name = self.options['cluster_type']
try:
self.cluster = cluster.from_name[cluster_name](**self.options)
except KeyError:
if aMCatNLO and ('mg5_path' not in self.options or not self.options['mg5_path']):
if not self.plugin_path:
raise self.InvalidCmd('%s not native cluster type and no plugin directory available.' % cluster_name)
elif aMCatNLO:
mg5dir = self.options['mg5_path']
if mg5dir not in sys.path:
sys.path.append(mg5dir)
if pjoin(mg5dir, 'PLUGIN') not in self.plugin_path:
self.plugin_path.append(pjoin(mg5dir))
else:
mg5dir = MG5DIR
# Check if a plugin define this type of cluster
# check for PLUGIN format
for plugpath in self.plugin_path:
plugindirname = os.path.basename(plugpath)
for plug in os.listdir(plugpath):
if os.path.exists(pjoin(plugpath, plug, '__init__.py')):
try:
__import__('%s.%s' % (plugindirname, plug))
except Exception, error:
logger.critical('plugin directory %s/%s fail to be loaded. Please check it',plugindirname, plug)
continue
plugin = sys.modules['%s.%s' % (plugindirname,plug)]
if not hasattr(plugin, 'new_cluster'):
continue
if not misc.is_plugin_supported(plugin):
continue
if cluster_name in plugin.new_cluster:
logger.info("cluster handling will be done with PLUGIN: %s" % plug,'$MG:color:BLACK')
self.cluster = plugin.new_cluster[cluster_name](**self.options)
break
if self.cluster_mode == 2:
try:
import multiprocessing
if not self.nb_core:
try:
self.nb_core = int(self.options['nb_core'])
except TypeError:
self.nb_core = multiprocessing.cpu_count()
logger.info('Using %d cores' % self.nb_core)
except ImportError:
self.nb_core = 1
logger.warning('Impossible to detect the number of cores => Using One.\n'+
'Use set nb_core X in order to set this number and be able to'+
'run in multicore.')
self.cluster = cluster.MultiCore(**self.options)
def clean_previous_results(self,options,p_dirs,folder_name):
"""Clean previous results.
o. If doing only the reweighting step, do not delete anything and return directlty.
o. Always remove all the G*_* files (from split event generation).
o. Remove the G* (or born_G* or all_G*) only when NOT doing only_generation or reweight_only."""
if options['reweightonly']:
return
if not options['only_generation']:
self.update_status('Cleaning previous results', level=None)
for dir in p_dirs:
#find old folders to be removed
for obj in folder_name:
# list all the G* (or all_G* or born_G*) directories
to_rm = [file for file in \
os.listdir(pjoin(self.me_dir, 'SubProcesses', dir)) \
if file.startswith(obj[:-1]) and \
(os.path.isdir(pjoin(self.me_dir, 'SubProcesses', dir, file)) or \
os.path.exists(pjoin(self.me_dir, 'SubProcesses', dir, file)))]
# list all the G*_* directories (from split event generation)
to_always_rm = [file for file in \
os.listdir(pjoin(self.me_dir, 'SubProcesses', dir)) \
if file.startswith(obj[:-1]) and
'_' in file and not '_G' in file and \
(os.path.isdir(pjoin(self.me_dir, 'SubProcesses', dir, file)) or \
os.path.exists(pjoin(self.me_dir, 'SubProcesses', dir, file)))]
if not options['only_generation']:
to_always_rm.extend(to_rm)
if os.path.exists(pjoin(self.me_dir, 'SubProcesses', dir,'MadLoop5_resources.tar.gz')):
to_always_rm.append(pjoin(self.me_dir, 'SubProcesses', dir,'MadLoop5_resources.tar.gz'))
files.rm([pjoin(self.me_dir, 'SubProcesses', dir, d) for d in to_always_rm])
return
def print_summary(self, options, step, mode, scale_pdf_info=[], done=True):
"""print a summary of the results contained in self.cross_sect_dict.
step corresponds to the mintMC step, if =2 (i.e. after event generation)
some additional infos are printed"""
# find process name
proc_card_lines = open(pjoin(self.me_dir, 'Cards', 'proc_card_mg5.dat')).read().split('\n')
process = ''
for line in proc_card_lines:
if line.startswith('generate') or line.startswith('add process'):
process = process+(line.replace('generate ', '')).replace('add process ','')+' ; '
lpp = {0:'l', 1:'p', -1:'pbar'}
if self.ninitial == 1:
proc_info = '\n Process %s' % process[:-3]
else:
proc_info = '\n Process %s\n Run at %s-%s collider (%s + %s GeV)' % \
(process[:-3], lpp[self.run_card['lpp1']], lpp[self.run_card['lpp2']],
self.run_card['ebeam1'], self.run_card['ebeam2'])
if self.ninitial == 1:
self.cross_sect_dict['unit']='GeV'
self.cross_sect_dict['xsec_string']='(Partial) decay width'
self.cross_sect_dict['axsec_string']='(Partial) abs(decay width)'
else:
self.cross_sect_dict['unit']='pb'
self.cross_sect_dict['xsec_string']='Total cross section'
self.cross_sect_dict['axsec_string']='Total abs(cross section)'
if mode in ['aMC@NLO', 'aMC@LO', 'noshower', 'noshowerLO']:
status = ['Determining the number of unweighted events per channel',
'Updating the number of unweighted events per channel',
'Summary:']
computed='(computed from LHE events)'
elif mode in ['NLO', 'LO']:
status = ['Results after grid setup:','Current results:',
'Final results and run summary:']
computed='(computed from histogram information)'
if step != 2 and mode in ['aMC@NLO', 'aMC@LO', 'noshower', 'noshowerLO']:
message = status[step] + '\n\n Intermediate results:' + \
('\n Random seed: %(randinit)d' + \
'\n %(xsec_string)s: %(xsect)8.3e +- %(errt)6.1e %(unit)s' + \
'\n %(axsec_string)s: %(xseca)8.3e +- %(erra)6.1e %(unit)s \n') \
% self.cross_sect_dict
elif mode in ['NLO','LO'] and not done:
if step == 0:
message = '\n ' + status[0] + \
'\n %(xsec_string)s: %(xsect)8.3e +- %(errt)6.1e %(unit)s' % \
self.cross_sect_dict
else:
message = '\n ' + status[1] + \
'\n %(xsec_string)s: %(xsect)8.3e +- %(errt)6.1e %(unit)s' % \
self.cross_sect_dict
else:
message = '\n --------------------------------------------------------------'
message = message + \
'\n ' + status[2] + proc_info
if mode not in ['LO', 'NLO']:
message = message + \
'\n Number of events generated: %s' % self.run_card['nevents']
message = message + \
'\n %(xsec_string)s: %(xsect)8.3e +- %(errt)6.1e %(unit)s' % \
self.cross_sect_dict
message = message + \
'\n --------------------------------------------------------------'
if scale_pdf_info and (self.run_card['nevents']>=10000 or mode in ['NLO', 'LO']):
if scale_pdf_info[0]:
# scale uncertainties
message = message + '\n Scale variation %s:' % computed
for s in scale_pdf_info[0]:
if s['unc']:
if self.run_card['ickkw'] != -1:
message = message + \
('\n Dynamical_scale_choice %(label)i (envelope of %(size)s values): '\
'\n %(cen)8.3e pb +%(max)0.1f%% -%(min)0.1f%%') % s
else:
message = message + \
('\n Soft and hard scale dependence (added in quadrature): '\
'\n %(cen)8.3e pb +%(max_q)0.1f%% -%(min_q)0.1f%%') % s
else:
message = message + \
('\n Dynamical_scale_choice %(label)i: '\
'\n %(cen)8.3e pb') % s
if scale_pdf_info[1]:
message = message + '\n PDF variation %s:' % computed
for p in scale_pdf_info[1]:
if p['unc']=='none':
message = message + \
('\n %(name)s (central value only): '\
'\n %(cen)8.3e pb') % p
elif p['unc']=='unknown':
message = message + \
('\n %(name)s (%(size)s members; combination method unknown): '\
'\n %(cen)8.3e pb') % p
else:
message = message + \
('\n %(name)s (%(size)s members; using %(unc)s method): '\
'\n %(cen)8.3e pb +%(max)0.1f%% -%(min)0.1f%%') % p
# pdf uncertainties
message = message + \
'\n --------------------------------------------------------------'
if (mode in ['NLO', 'LO'] and not done) or \
(mode in ['aMC@NLO', 'aMC@LO', 'noshower', 'noshowerLO'] and step!=2):
logger.info(message+'\n')
return
# Some advanced general statistics are shown in the debug message at the
# end of the run
# Make sure it never stops a run
# Gather some basic statistics for the run and extracted from the log files.
if mode in ['aMC@NLO', 'aMC@LO', 'noshower', 'noshowerLO']:
log_GV_files = misc.glob(pjoin('P*','G*','log_MINT*.txt'),
pjoin(self.me_dir, 'SubProcesses'))
all_log_files = log_GV_files
elif mode == 'NLO':
log_GV_files = misc.glob(pjoin('P*','all_G*','log_MINT*.txt'),
pjoin(self.me_dir, 'SubProcesses'))
all_log_files = log_GV_files
elif mode == 'LO':
log_GV_files = ''
all_log_files = misc.glob(pjoin('P*','born_G*','log_MINT*.txt'),
pjoin(self.me_dir, 'SubProcesses'))
else:
raise aMCatNLOError, 'Running mode %s not supported.'%mode
try:
message, debug_msg = \
self.compile_advanced_stats(log_GV_files, all_log_files, message)
except Exception as e:
debug_msg = 'Advanced statistics collection failed with error "%s"\n'%str(e)
err_string = StringIO.StringIO()
traceback.print_exc(limit=4, file=err_string)
debug_msg += 'Please report this backtrace to a MadGraph developer:\n%s'\
%err_string.getvalue()
logger.debug(debug_msg+'\n')
logger.info(message+'\n')
# Now copy relevant information in the Events/Run_ directory
evt_path = pjoin(self.me_dir, 'Events', self.run_name)
open(pjoin(evt_path, 'summary.txt'),'w').write(message+'\n')
open(pjoin(evt_path, '.full_summary.txt'),
'w').write(message+'\n\n'+debug_msg+'\n')
self.archive_files(evt_path,mode)
def archive_files(self, evt_path, mode):
""" Copies in the Events/Run_ directory relevant files characterizing
the run."""
files_to_arxiv = [pjoin('Cards','param_card.dat'),
pjoin('Cards','MadLoopParams.dat'),
pjoin('Cards','FKS_params.dat'),
pjoin('Cards','run_card.dat'),
pjoin('Subprocesses','setscales.f'),
pjoin('Subprocesses','cuts.f')]
if mode in ['NLO', 'LO']:
files_to_arxiv.append(pjoin('Cards','FO_analyse_card.dat'))
if not os.path.exists(pjoin(evt_path,'RunMaterial')):
os.mkdir(pjoin(evt_path,'RunMaterial'))
for path in files_to_arxiv:
if os.path.isfile(pjoin(self.me_dir,path)):
files.cp(pjoin(self.me_dir,path),pjoin(evt_path,'RunMaterial'))
misc.call(['tar','-czpf','RunMaterial.tar.gz','RunMaterial'],cwd=evt_path)
shutil.rmtree(pjoin(evt_path,'RunMaterial'))
def compile_advanced_stats(self,log_GV_files,all_log_files,message):
""" This functions goes through the log files given in arguments and
compiles statistics about MadLoop stability, virtual integration
optimization and detection of potential error messages into a nice
debug message to printed at the end of the run """
def safe_float(str_float):
try:
return float(str_float)
except ValueError:
logger.debug('Could not convert the following float during'+
' advanced statistics printout: %s'%str(str_float))
return -1.0
# > UPS is a dictionary of tuples with this format {channel:[nPS,nUPS]}
# > Errors is a list of tuples with this format (log_file,nErrors)
stats = {'UPS':{}, 'Errors':[], 'virt_stats':{}, 'timings':{}}
mint_search = re.compile(r"MINT(?P\d*).txt")
# ==================================
# == MadLoop stability statistics ==
# ==================================
# Recuperate the fraction of unstable PS points found in the runs for
# the virtuals
UPS_stat_finder = re.compile(
r"Satistics from MadLoop:.*"+\
r"Total points tried\:\s+(?P\d+).*"+\
r"Stability unknown\:\s+(?P\d+).*"+\
r"Stable PS point\:\s+(?P\d+).*"+\
r"Unstable PS point \(and rescued\)\:\s+(?P\d+).*"+\
r"Exceptional PS point \(unstable and not rescued\)\:\s+(?P\d+).*"+\
r"Double precision used\:\s+(?P\d+).*"+\
r"Quadruple precision used\:\s+(?P\d+).*"+\
r"Initialization phase\-space points\:\s+(?P\d+).*"+\
r"Unknown return code \(100\)\:\s+(?P\d+).*"+\
r"Unknown return code \(10\)\:\s+(?P\d+).*",re.DOTALL)
unit_code_meaning = { 0 : 'Not identified (CTModeRun != -1)',
1 : 'CutTools (double precision)',
2 : 'PJFry++',
3 : 'IREGI',
4 : 'Golem95',
5 : 'Samurai',
6 : 'Ninja (double precision)',
7 : 'COLLIER',
8 : 'Ninja (quadruple precision)',
9 : 'CutTools (quadruple precision)'}
RetUnit_finder =re.compile(
r"#Unit\s*(?P\d+)\s*=\s*(?P\d+)")
#Unit
for gv_log in log_GV_files:
channel_name = '/'.join(gv_log.split('/')[-5:-1])
log=open(gv_log,'r').read()
UPS_stats = re.search(UPS_stat_finder,log)
for retunit_stats in re.finditer(RetUnit_finder, log):
if channel_name not in stats['UPS'].keys():
stats['UPS'][channel_name] = [0]*10+[[0]*10]
stats['UPS'][channel_name][10][int(retunit_stats.group('unit'))] \
+= int(retunit_stats.group('n_occurences'))
if not UPS_stats is None:
try:
stats['UPS'][channel_name][0] += int(UPS_stats.group('ntot'))
stats['UPS'][channel_name][1] += int(UPS_stats.group('nsun'))
stats['UPS'][channel_name][2] += int(UPS_stats.group('nsps'))
stats['UPS'][channel_name][3] += int(UPS_stats.group('nups'))
stats['UPS'][channel_name][4] += int(UPS_stats.group('neps'))
stats['UPS'][channel_name][5] += int(UPS_stats.group('nddp'))
stats['UPS'][channel_name][6] += int(UPS_stats.group('nqdp'))
stats['UPS'][channel_name][7] += int(UPS_stats.group('nini'))
stats['UPS'][channel_name][8] += int(UPS_stats.group('n100'))
stats['UPS'][channel_name][9] += int(UPS_stats.group('n10'))
except KeyError:
stats['UPS'][channel_name] = [int(UPS_stats.group('ntot')),
int(UPS_stats.group('nsun')),int(UPS_stats.group('nsps')),
int(UPS_stats.group('nups')),int(UPS_stats.group('neps')),
int(UPS_stats.group('nddp')),int(UPS_stats.group('nqdp')),
int(UPS_stats.group('nini')),int(UPS_stats.group('n100')),
int(UPS_stats.group('n10')),[0]*10]
debug_msg = ""
if len(stats['UPS'].keys())>0:
nTotPS = sum([chan[0] for chan in stats['UPS'].values()],0)
nTotsun = sum([chan[1] for chan in stats['UPS'].values()],0)
nTotsps = sum([chan[2] for chan in stats['UPS'].values()],0)
nTotups = sum([chan[3] for chan in stats['UPS'].values()],0)
nToteps = sum([chan[4] for chan in stats['UPS'].values()],0)
nTotddp = sum([chan[5] for chan in stats['UPS'].values()],0)
nTotqdp = sum([chan[6] for chan in stats['UPS'].values()],0)
nTotini = sum([chan[7] for chan in stats['UPS'].values()],0)
nTot100 = sum([chan[8] for chan in stats['UPS'].values()],0)
nTot10 = sum([chan[9] for chan in stats['UPS'].values()],0)
nTot1 = [sum([chan[10][i] for chan in stats['UPS'].values()],0) \
for i in range(10)]
UPSfracs = [(chan[0] , 0.0 if chan[1][0]==0 else \
safe_float(chan[1][4]*100)/chan[1][0]) for chan in stats['UPS'].items()]
maxUPS = max(UPSfracs, key = lambda w: w[1])
tmpStr = ""
tmpStr += '\n Number of loop ME evaluations (by MadLoop): %d'%nTotPS
tmpStr += '\n Stability unknown: %d'%nTotsun
tmpStr += '\n Stable PS point: %d'%nTotsps
tmpStr += '\n Unstable PS point (and rescued): %d'%nTotups
tmpStr += '\n Unstable PS point (and not rescued): %d'%nToteps
tmpStr += '\n Only double precision used: %d'%nTotddp
tmpStr += '\n Quadruple precision used: %d'%nTotqdp
tmpStr += '\n Initialization phase-space points: %d'%nTotini
tmpStr += '\n Reduction methods used:'
red_methods = [(unit_code_meaning[i],nTot1[i]) for i in \
unit_code_meaning.keys() if nTot1[i]>0]
for method, n in sorted(red_methods, key= lambda l: l[1], reverse=True):
tmpStr += '\n > %s%s%s'%(method,' '*(33-len(method)),n)
if nTot100 != 0:
debug_msg += '\n Unknown return code (100): %d'%nTot100
if nTot10 != 0:
debug_msg += '\n Unknown return code (10): %d'%nTot10
nUnknownUnit = sum(nTot1[u] for u in range(10) if u \
not in unit_code_meaning.keys())
if nUnknownUnit != 0:
debug_msg += '\n Unknown return code (1): %d'\
%nUnknownUnit
if maxUPS[1]>0.001:
message += tmpStr
message += '\n Total number of unstable PS point detected:'+\
' %d (%4.2f%%)'%(nToteps,safe_float(100*nToteps)/nTotPS)
message += '\n Maximum fraction of UPS points in '+\
'channel %s (%4.2f%%)'%maxUPS
message += '\n Please report this to the authors while '+\
'providing the file'
message += '\n %s'%str(pjoin(os.path.dirname(self.me_dir),
maxUPS[0],'UPS.log'))
else:
debug_msg += tmpStr
# ====================================================
# == aMC@NLO virtual integration optimization stats ==
# ====================================================
virt_tricks_finder = re.compile(
r"accumulated results Virtual ratio\s*=\s*-?(?P[\d\+-Eed\.]*)"+\
r"\s*\+/-\s*-?[\d\+-Eed\.]*\s*\(\s*-?(?P[\d\+-Eed\.]*)\s*\%\)\s*\n"+\
r"accumulated results ABS virtual\s*=\s*-?(?P[\d\+-Eed\.]*)"+\
r"\s*\+/-\s*-?[\d\+-Eed\.]*\s*\(\s*-?(?P[\d\+-Eed\.]*)\s*\%\)")
virt_frac_finder = re.compile(r"update virtual fraction to\s*:\s*"+\
"-?(?P[\d\+-Eed\.]*)\s*-?(?P[\d\+-Eed\.]*)")
channel_contr_finder = re.compile(r"Final result \[ABS\]\s*:\s*-?(?P[\d\+-Eed\.]*)")
channel_contr_list = {}
for gv_log in log_GV_files:
logfile=open(gv_log,'r')
log = logfile.read()
logfile.close()
channel_name = '/'.join(gv_log.split('/')[-3:-1])
vf_stats = None
for vf_stats in re.finditer(virt_frac_finder, log):
pass
if not vf_stats is None:
v_frac = safe_float(vf_stats.group('v_frac'))
v_average = safe_float(vf_stats.group('v_average'))
try:
if v_frac < stats['virt_stats']['v_frac_min'][0]:
stats['virt_stats']['v_frac_min']=(v_frac,channel_name)
if v_frac > stats['virt_stats']['v_frac_max'][0]:
stats['virt_stats']['v_frac_max']=(v_frac,channel_name)
stats['virt_stats']['v_frac_avg'][0] += v_frac
stats['virt_stats']['v_frac_avg'][1] += 1
except KeyError:
stats['virt_stats']['v_frac_min']=[v_frac,channel_name]
stats['virt_stats']['v_frac_max']=[v_frac,channel_name]
stats['virt_stats']['v_frac_avg']=[v_frac,1]
ccontr_stats = None
for ccontr_stats in re.finditer(channel_contr_finder, log):
pass
if not ccontr_stats is None:
contrib = safe_float(ccontr_stats.group('v_contr'))
try:
if contrib>channel_contr_list[channel_name]:
channel_contr_list[channel_name]=contrib
except KeyError:
channel_contr_list[channel_name]=contrib
# Now build the list of relevant virt log files to look for the maxima
# of virt fractions and such.
average_contrib = 0.0
for value in channel_contr_list.values():
average_contrib += value
if len(channel_contr_list.values()) !=0:
average_contrib = average_contrib / len(channel_contr_list.values())
relevant_log_GV_files = []
excluded_channels = set([])
all_channels = set([])
for log_file in log_GV_files:
channel_name = '/'.join(log_file.split('/')[-3:-1])
all_channels.add(channel_name)
try:
if channel_contr_list[channel_name] > (0.1*average_contrib):
relevant_log_GV_files.append(log_file)
else:
excluded_channels.add(channel_name)
except KeyError:
relevant_log_GV_files.append(log_file)
# Now we want to use the latest occurence of accumulated result in the log file
for gv_log in relevant_log_GV_files:
logfile=open(gv_log,'r')
log = logfile.read()
logfile.close()
channel_name = '/'.join(gv_log.split('/')[-3:-1])
vt_stats = None
for vt_stats in re.finditer(virt_tricks_finder, log):
pass
if not vt_stats is None:
vt_stats_group = vt_stats.groupdict()
v_ratio = safe_float(vt_stats.group('v_ratio'))
v_ratio_err = safe_float(vt_stats.group('v_ratio_err'))
v_contr = safe_float(vt_stats.group('v_abs_contr'))
v_contr_err = safe_float(vt_stats.group('v_abs_contr_err'))
try:
if v_ratio < stats['virt_stats']['v_ratio_min'][0]:
stats['virt_stats']['v_ratio_min']=(v_ratio,channel_name)
if v_ratio > stats['virt_stats']['v_ratio_max'][0]:
stats['virt_stats']['v_ratio_max']=(v_ratio,channel_name)
if v_ratio < stats['virt_stats']['v_ratio_err_min'][0]:
stats['virt_stats']['v_ratio_err_min']=(v_ratio_err,channel_name)
if v_ratio > stats['virt_stats']['v_ratio_err_max'][0]:
stats['virt_stats']['v_ratio_err_max']=(v_ratio_err,channel_name)
if v_contr < stats['virt_stats']['v_contr_min'][0]:
stats['virt_stats']['v_contr_min']=(v_contr,channel_name)
if v_contr > stats['virt_stats']['v_contr_max'][0]:
stats['virt_stats']['v_contr_max']=(v_contr,channel_name)
if v_contr_err < stats['virt_stats']['v_contr_err_min'][0]:
stats['virt_stats']['v_contr_err_min']=(v_contr_err,channel_name)
if v_contr_err > stats['virt_stats']['v_contr_err_max'][0]:
stats['virt_stats']['v_contr_err_max']=(v_contr_err,channel_name)
except KeyError:
stats['virt_stats']['v_ratio_min']=[v_ratio,channel_name]
stats['virt_stats']['v_ratio_max']=[v_ratio,channel_name]
stats['virt_stats']['v_ratio_err_min']=[v_ratio_err,channel_name]
stats['virt_stats']['v_ratio_err_max']=[v_ratio_err,channel_name]
stats['virt_stats']['v_contr_min']=[v_contr,channel_name]
stats['virt_stats']['v_contr_max']=[v_contr,channel_name]
stats['virt_stats']['v_contr_err_min']=[v_contr_err,channel_name]
stats['virt_stats']['v_contr_err_max']=[v_contr_err,channel_name]
vf_stats = None
for vf_stats in re.finditer(virt_frac_finder, log):
pass
if not vf_stats is None:
v_frac = safe_float(vf_stats.group('v_frac'))
v_average = safe_float(vf_stats.group('v_average'))
try:
if v_average < stats['virt_stats']['v_average_min'][0]:
stats['virt_stats']['v_average_min']=(v_average,channel_name)
if v_average > stats['virt_stats']['v_average_max'][0]:
stats['virt_stats']['v_average_max']=(v_average,channel_name)
stats['virt_stats']['v_average_avg'][0] += v_average
stats['virt_stats']['v_average_avg'][1] += 1
except KeyError:
stats['virt_stats']['v_average_min']=[v_average,channel_name]
stats['virt_stats']['v_average_max']=[v_average,channel_name]
stats['virt_stats']['v_average_avg']=[v_average,1]
try:
debug_msg += '\n\n Statistics on virtual integration optimization : '
debug_msg += '\n Maximum virt fraction computed %.3f (%s)'\
%tuple(stats['virt_stats']['v_frac_max'])
debug_msg += '\n Minimum virt fraction computed %.3f (%s)'\
%tuple(stats['virt_stats']['v_frac_min'])
debug_msg += '\n Average virt fraction computed %.3f'\
%safe_float(stats['virt_stats']['v_frac_avg'][0]/safe_float(stats['virt_stats']['v_frac_avg'][1]))
debug_msg += '\n Stats below exclude negligible channels (%d excluded out of %d)'%\
(len(excluded_channels),len(all_channels))
debug_msg += '\n Maximum virt ratio used %.2f (%s)'\
%tuple(stats['virt_stats']['v_average_max'])
debug_msg += '\n Maximum virt ratio found from grids %.2f (%s)'\
%tuple(stats['virt_stats']['v_ratio_max'])
tmpStr = '\n Max. MC err. on virt ratio from grids %.1f %% (%s)'\
%tuple(stats['virt_stats']['v_ratio_err_max'])
debug_msg += tmpStr
# After all it was decided that it is better not to alarm the user unecessarily
# with such printout of the statistics.
# if stats['virt_stats']['v_ratio_err_max'][0]>100.0 or \
# stats['virt_stats']['v_ratio_err_max'][0]>100.0:
# message += "\n Suspiciously large MC error in :"
# if stats['virt_stats']['v_ratio_err_max'][0]>100.0:
# message += tmpStr
tmpStr = '\n Maximum MC error on abs virt %.1f %% (%s)'\
%tuple(stats['virt_stats']['v_contr_err_max'])
debug_msg += tmpStr
# if stats['virt_stats']['v_contr_err_max'][0]>100.0:
# message += tmpStr
except KeyError:
debug_msg += '\n Could not find statistics on the integration optimization. '
# =======================================
# == aMC@NLO timing profile statistics ==
# =======================================
timing_stat_finder = re.compile(r"\s*Time spent in\s*(?P\w*)\s*:\s*"+\
"(?P