No OneTemporary
Actions

Size

69 KB

Subscribers

None

View Options

	diff --git a/bin/rivet b/bin/rivet
	--- a/bin/rivet
	+++ b/bin/rivet
	@@ -1,675 +1,678 @@
	#! /usr/bin/env python

	"""\
	Run Rivet analyses on HepMC events read from a file or Unix pipe

	Examples:

	%prog [options] <hepmcfile> [<hepmcfile2> ...]

	or

	mkfifo fifo.hepmc
	my_generator -o fifo.hepmc &
	%prog [options] fifo.hepmc


	ENVIRONMENT:
	* RIVET_ANALYSIS_PATH: list of paths to be searched for plugin analysis libraries at runtime
	* RIVET_DATA_PATH: list of paths to be searched for data files (defaults to use analysis path)
	* RIVET_WEIGHT_INDEX: the numerical weight-vector index to use in this run
	* See the documentation for more environment variables.
	"""

	from __future__ import print_function
	import os, sys

	## Load the rivet module
	try:
	import rivet
	except:
	## If rivet loading failed, try to bootstrap the Python path!
	try:
	# TODO: Is this a good idea? Maybe just notify the user that their PYTHONPATH is wrong?
	import commands
	modname = sys.modules[__name__].__file__
	binpath = os.path.dirname(modname)
	rivetconfigpath = os.path.join(binpath, "rivet-config")
	rivetpypath = commands.getoutput(rivetconfigpath + " --pythonpath")
	sys.path.append(rivetpypath)
	import rivet
	except:
	sys.stderr.write("The rivet Python module could not be loaded: is your PYTHONPATH set correctly?\n")
	sys.exit(5)

	rivet.util.check_python_version()
	rivet.util.set_process_name("rivet")

	import time, datetime, logging, signal

	## Parse command line options
	from optparse import OptionParser, OptionGroup
	parser = OptionParser(usage=__doc__, version="rivet v%s" % rivet.version())

	anagroup = OptionGroup(parser, "Analysis handling")
	anagroup.add_option("-a", "--analysis", "--analyses", dest="ANALYSES", action="append",
	default=[], metavar="ANA",
	help="add an analysis (or comma-separated list of analyses) to the processing list.")
	anagroup.add_option("--list-analyses", "--list", dest="LIST_ANALYSES", action="store_true",
	default=False, help="show the list of available analyses' names. With -v, it shows the descriptions, too")
	anagroup.add_option("--list-keywords", "--keywords", dest="LIST_KEYWORDS", action="store_true",
	default=False, help="show the list of available keywords.")
	anagroup.add_option("--list-used-analyses", action="store_true", dest="LIST_USED_ANALYSES",
	default=False, help="list the analyses used by this command (after subtraction of inappropriate ones)")
	anagroup.add_option("--show-analysis", "--show-analyses", "--show", dest="SHOW_ANALYSES", action="append",
	default=[], help="show the details of an analysis")
	anagroup.add_option("--show-bibtex", dest="SHOW_BIBTEX", action="store_true",
	default=False, help="show BibTeX entries for all used analyses")
	anagroup.add_option("--analysis-path", dest="ANALYSIS_PATH", metavar="PATH", default=None,
	help="specify the analysis search path (cf. $RIVET_ANALYSIS_PATH).")
	# TODO: remove/deprecate the append?
	anagroup.add_option("--analysis-path-append", dest="ANALYSIS_PATH_APPEND", metavar="PATH", default=None,
	help="append to the analysis search path (cf. $RIVET_ANALYSIS_PATH).")
	anagroup.add_option("--pwd", dest="ANALYSIS_PATH_PWD", action="store_true", default=False,
	help="append the current directory (pwd) to the analysis/data search paths (cf. $RIVET_ANALYSIS_PATH).")
	# TODO: add control for more paths?
	parser.add_option_group(anagroup)


	extragroup = OptionGroup(parser, "Extra run settings")
	extragroup.add_option("-o", "-H", "--histo-file", dest="HISTOFILE",
	default="Rivet.yoda", help="specify the output histo file path (default = %default)")
	extragroup.add_option("-p", "--preload-file", dest="PRELOADFILE",
	default=None, help="specify and old yoda file to initialize (default = %default)")
	extragroup.add_option("--no-histo-file", dest="WRITE_DATA", action="store_false", default=True,
	help="don't write out any histogram file at the end of the run (default = write)")
	extragroup.add_option("-x", "--cross-section", dest="CROSS_SECTION",
	default=None, metavar="XS",
	help="specify the signal process cross-section in pb")
	extragroup.add_option("-n", "--nevts", dest="MAXEVTNUM", type="int",
	default=None, metavar="NUM",
	help="restrict the max number of events to process")
	extragroup.add_option("--nskip", dest="EVTSKIPNUM", type="int",
	default=0, metavar="NUM",
	help="skip NUM events read from input before beginning processing")
	+extragroup.add_option("--skip-weights", dest="SKIP_WEIGHTS", action="store_true",
	+ default=False, help="only run on the nominal weight")
	extragroup.add_option("--runname", dest="RUN_NAME", default=None, metavar="NAME",
	help="give an optional run name, to be prepended as a 'top level directory' in histo paths")
	extragroup.add_option("--ignore-beams", dest="IGNORE_BEAMS", action="store_true", default=False,
	help="ignore input event beams when checking analysis compatibility. "
	"WARNING: analyses may not work correctly, or at all, with inappropriate beams")
	parser.add_option_group(extragroup)


	timinggroup = OptionGroup(parser, "Timeouts and periodic operations")
	timinggroup.add_option("--event-timeout", dest="EVENT_TIMEOUT", type="int",
	default=21600, metavar="NSECS",
	help="max time in whole seconds to wait for an event to be generated from the specified source (default = %default)")
	timinggroup.add_option("--run-timeout", dest="RUN_TIMEOUT", type="int",
	default=None, metavar="NSECS",
	help="max time in whole seconds to wait for the run to finish. This can be useful on batch systems such "
	"as the LCG Grid where tokens expire on a fixed wall-clock and can render long Rivet runs unable to write "
	"out the final histogram file (default = unlimited)")
	timinggroup.add_option("--histo-interval", dest="HISTO_WRITE_INTERVAL", type=int,
	default=1000, help="specify the number of events between histogram file updates, default = %default. "
	"Set to 0 to only write out at the end of the run. Note that intermediate histograms will be those "
	"from the analyze step only: analysis finalizing is currently not executed until the end of the run.")
	parser.add_option_group(timinggroup)


	verbgroup = OptionGroup(parser, "Verbosity control")
	parser.add_option("-l", dest="NATIVE_LOG_STRS", action="append",
	default=[], help="set a log level in the Rivet library")
	verbgroup.add_option("-v", "--verbose", action="store_const", const=logging.DEBUG, dest="LOGLEVEL",
	default=logging.INFO, help="print debug (very verbose) messages")
	verbgroup.add_option("-q", "--quiet", action="store_const", const=logging.WARNING, dest="LOGLEVEL",
	default=logging.INFO, help="be very quiet")
	parser.add_option_group(verbgroup)
	opts, args = parser.parse_args()


	## Override/modify analysis search path
	if opts.ANALYSIS_PATH:
	rivet.setAnalysisLibPaths(opts.ANALYSIS_PATH.split(":"))
	rivet.setAnalysisDataPaths(opts.ANALYSIS_PATH.split(":"))
	if opts.ANALYSIS_PATH_APPEND:
	for ap in opts.ANALYSIS_PATH_APPEND.split(":"):
	rivet.addAnalysisLibPath(ap)
	rivet.addAnalysisDataPath(ap)
	if opts.ANALYSIS_PATH_PWD:
	rivet.addAnalysisLibPath(os.path.abspath("."))
	rivet.addAnalysisDataPath(os.path.abspath("."))


	## Configure logging
	logging.basicConfig(level=opts.LOGLEVEL, format="%(message)s")
	for l in opts.NATIVE_LOG_STRS:
	name, level = None, None
	try:
	name, level = l.split("=")
	except:
	name = "Rivet"
	level = l
	## Fix name
	if name != "Rivet" and not name.startswith("Rivet."):
	name = "Rivet." + name
	try:
	## Get right error type
	level = rivet.LEVELS.get(level.upper(), None)
	logging.debug("Setting log level: %s %d" % (name, level))
	rivet.setLogLevel(name, level)
	except:
	logging.warning("Couldn't process logging string '%s'" % l)



	############################
	## Listing available analyses/keywords


	def getAnalysesByKeyword(alist, kstring):
	add, veto, ret = [], [], []
	bits = [i for i in kstring.replace("^@", "@^").split("@") if len(i) > 0]
	for b in bits:
	if b.startswith("^"):
	veto.append(b.strip("^"))
	else:
	add.append(b)

	add = set(add)
	veto = set(veto)

	for a in alist:
	kwds = set([i.lower() for i in rivet.AnalysisLoader.getAnalysis(a).keywords()])
	if kwds.intersection(veto) and len(kwds.intersection(add)) == len(list(add)):
	ret.append(a)
	return ret


	## List of analyses
	all_analyses = rivet.AnalysisLoader.analysisNames()
	if opts.LIST_ANALYSES:
	## Treat args as case-insensitive regexes if present
	regexes = None
	if args:
	import re
	regexes = [re.compile(arg, re.I) for arg in args]
	# import tempfile, subprocess
	# tf, tfpath = tempfile.mkstemp(prefix="rivet-list.")
	names = []
	msg = []
	for aname in all_analyses:
	if not regexes:
	toshow = True
	else:
	toshow = False
	for regex in regexes:
	if regex.search(aname):
	toshow = True
	break
	if toshow:
	names.append(aname)
	msg.append('')
	if opts.LOGLEVEL <= logging.INFO:
	a = rivet.AnalysisLoader.getAnalysis(aname)
	st = "" if a.status() == "VALIDATED" else ("[%s] " % a.status())
	# detex will very likely introduce some non-ASCII chars from
	# greek names in analysis titles.
	# The u"" prefix and explicit print encoding are necessary for
	# py2 to handle this properly
	msg[-1] = u"%s%s" % (st, a.summary())
	if opts.LOGLEVEL < logging.INFO:
	if a.keywords():
	msg[-1] += u" [%s]" % " ".join(a.keywords())
	if a.luminosityfb():
	msg[-1] += u" [ \int L = %s fb^{-1} ]" % a.luminosityfb()
	msg = rivet.util.detex(msg)
	retlist = '\n'.join([ u"%-25s %s" % (a,m) for a,m in zip(names,msg) ])
	if type(u'') is not str:
	retlist = retlist.encode('utf-8')
	print(retlist)
	sys.exit(0)


	def getKeywords(alist):
	all_keywords = []
	for a in alist:
	all_keywords.extend(rivet.AnalysisLoader.getAnalysis(a).keywords())
	all_keywords = [i.lower() for i in all_keywords]
	return sorted(list(set(all_keywords)))


	## List keywords
	if opts.LIST_KEYWORDS:
	# a = rivet.AnalysisLoader.getAnalysis(aname)
	for k in getKeywords(all_analyses):
	print(k)
	sys.exit(0)


	## Show analyses' details
	if len(opts.SHOW_ANALYSES) > 0:
	toshow = []
	for i, a in enumerate(opts.SHOW_ANALYSES):
	a_up = a.upper()
	if a_up in all_analyses and a_up not in toshow:
	toshow.append(a_up)
	else:
	## Treat as a case-insensitive regex
	import re
	regex = re.compile(a, re.I)
	for ana in all_analyses:
	if regex.search(ana) and a_up not in toshow:
	toshow.append(ana)

	msgs = []
	for i, name in enumerate(sorted(toshow)):
	import textwrap
	ana = rivet.AnalysisLoader.getAnalysis(name)

	msg = ""
	msg += name + "\n"
	msg += (len(name) * "=") + "\n\n"
	msg += rivet.util.detex(ana.summary()) + "\n\n"
	msg += "Status: " + ana.status() + "\n\n"

	# TODO: reduce to only show Inspire in v3
	if ana.inspireId():
	msg += "Inspire ID: " + ana.inspireId() + "\n"
	msg += "Inspire URL: http://inspire-hep.net/record/" + ana.inspireId() + "\n"
	msg += "HepData URL: http://hepdata.net/record/ins" + ana.inspireId() + "\n"
	elif ana.spiresId():
	msg += "Spires ID: " + ana.spiresId() + "\n"
	msg += "Inspire URL: http://inspire-hep.net/search?p=find+key+" + ana.spiresId() + "\n"
	msg += "HepData URL: http://hepdata.cedar.ac.uk/view/irn" + ana.spiresId() + "\n"

	if ana.year():
	msg += "Year of publication: " + ana.year() + "\n"
	if ana.bibKey():
	msg += "BibTeX key: " + ana.bibKey() + "\n"

	msg += "Authors:\n"
	for a in ana.authors():
	msg += " " + a + "\n"
	msg += "\n"

	msg += "Description:\n"
	twrap = textwrap.TextWrapper(width=75, initial_indent=2" ", subsequent_indent=2" ")
	msg += twrap.fill(rivet.util.detex(ana.description())) + "\n\n"

	if ana.experiment():
	msg += "Experiment: " + ana.experiment()
	if ana.collider():
	msg += "(%s)" % ana.collider()
	msg += "\n"
	# TODO: move this formatting into Analysis or a helper function?
	if ana.requiredBeams():
	def pid_to_str(pid):
	if pid == 11:
	return "e-"
	elif pid == -11:
	return "e+"
	elif pid == 2212:
	return "p+"
	elif pid == -2212:
	return "p-"
	elif pid == 10000:
	return "*"
	else:
	return str(pid)
	beamstrs = []
	for bp in ana.requiredBeams():
	beamstrs.append(pid_to_str(bp[0]) + " " + pid_to_str(bp[1]))
	msg += "Beams:" + ", ".join(beamstrs) + "\n"

	if ana.requiredEnergies():
	msg += "Beam energies:" + "; ".join(["(%0.1f, %0.1f) GeV\n" % (epair[0], epair[1]) for epair in ana.requiredEnergies()])
	else:
	msg += "Beam energies: ANY\n"

	if ana.runInfo():
	msg += "Run details:\n"
	twrap = textwrap.TextWrapper(width=75, initial_indent=2" ", subsequent_indent=4" ")
	for l in ana.runInfo().split("\n"):
	msg += twrap.fill(l) + "\n"

	if ana.luminosityfb():
	msg+= "\nIntegrated data luminosity = %s inverse fb.\n"%ana.luminosityfb()

	if ana.keywords():
	msg += "\nAnalysis keywords:"
	for k in ana.keywords():
	msg += " %s"%k
	msg+= "\n\n"

	if ana.references():
	msg += "\n" + "References:\n"
	for r in ana.references():
	url = None
	if r.startswith("arXiv:"):
	code = r.split()[0].replace("arXiv:", "")
	url = "http://arxiv.org/abs/" + code
	elif r.startswith("doi:"):
	code = r.replace("doi:", "")
	url = "http://dx.doi.org/" + code
	if url is not None:
	r += " - " + url
	msg += " " + r + "\n"

	## Add to the output
	msgs.append(msg)

	## Write the combined messages to a temporary file and page it
	if msgs:
	try:
	import tempfile, subprocess
	tffd, tfpath = tempfile.mkstemp(prefix="rivet-show.")
	os.write(tffd, u"\n\n".join(msgs).encode('utf-8'))
	if sys.stdout.isatty():
	pager = subprocess.Popen(["less", "-FX", tfpath]) #, stdin=subprocess.PIPE)
	pager.communicate()
	else:
	f = open(tfpath)
	print(f.read())
	f.close()
	finally:
	os.unlink(tfpath) #< always clean up
	sys.exit(0)



	############################
	## Actual analysis runs



	## We allow comma-separated lists of analysis names -- normalise the list here
	newanas = []
	for a in opts.ANALYSES:
	if "," in a:
	newanas += a.split(",")
	elif "@" in a: #< NB. this bans combination of ana lists and keywords in a single arg
	temp = getAnalysesByKeyword(all_analyses, a)
	for i in temp:
	newanas.append(i)
	else:
	newanas.append(a)
	opts.ANALYSES = newanas


	## Parse supplied cross-section
	if opts.CROSS_SECTION is not None:
	xsstr = opts.CROSS_SECTION
	try:
	opts.CROSS_SECTION = float(xsstr)
	except:
	import re
	suffmatch = re.search(r"[^\d.]", xsstr)
	if not suffmatch:
	raise ValueError("Bad cross-section string: %s" % xsstr)
	factor = base = None
	suffstart = suffmatch.start()
	if suffstart != -1:
	base = xsstr[:suffstart]
	suffix = xsstr[suffstart:].lower()
	if suffix == "mb":
	factor = 1e+9
	elif suffix == "mub":
	factor = 1e+6
	elif suffix == "nb":
	factor = 1e+3
	elif suffix == "pb":
	factor = 1
	elif suffix == "fb":
	factor = 1e-3
	elif suffix == "ab":
	factor = 1e-6
	if factor is None or base is None:
	raise ValueError("Bad cross-section string: %s" % xsstr)
	xs = float(base) * factor
	opts.CROSS_SECTION = xs


	## Print the available CLI options!
	#if opts.LIST_OPTIONS:
	# for o in parser.option_list:
	# print(o.get_opt_string())
	# sys.exit(0)


	## Set up signal handling
	RECVD_KILL_SIGNAL = None
	def handleKillSignal(signum, frame):
	"Declare us as having been signalled, and return to default handling behaviour"
	global RECVD_KILL_SIGNAL
	logging.critical("Signal handler called with signal " + str(signum))
	RECVD_KILL_SIGNAL = signum
	signal.signal(signum, signal.SIG_DFL)
	## Signals to handle
	signal.signal(signal.SIGTERM, handleKillSignal);
	signal.signal(signal.SIGHUP, handleKillSignal);
	signal.signal(signal.SIGINT, handleKillSignal);
	signal.signal(signal.SIGUSR1, handleKillSignal);
	signal.signal(signal.SIGUSR2, handleKillSignal);
	try:
	signal.signal(signal.SIGXCPU, handleKillSignal);
	except:
	pass



	## Identify HepMC files/streams
	## TODO: check readability, deal with stdin
	if len(args) > 0:
	HEPMCFILES = args
	else:
	HEPMCFILES = ["-"]


	## Event number logging
	def logNEvt(n, starttime, maxevtnum):
	if n % 10000 == 0:
	nevtloglevel = logging.CRITICAL
	elif n % 1000 == 0:
	nevtloglevel = logging.WARNING
	elif n % 100 == 0:
	nevtloglevel = logging.INFO
	else:
	nevtloglevel = logging.DEBUG
	currenttime = datetime.datetime.now().replace(microsecond=0)
	elapsedtime = currenttime - starttime
	logging.log(nevtloglevel, "Event %d (%s elapsed)" % (n, str(elapsedtime)))
	# if maxevtnum is None:
	# logging.log(nevtloglevel, "Event %d (%s elapsed)" % (n, str(elapsedtime)))
	# else:
	# remainingtime = (maxevtnum-n) * elapsedtime.total_seconds() / float(n)
	# eta = time.strftime("%a %b %d %H:%M", datetime.localtime(currenttime + remainingtime))
	# logging.log(nevtloglevel, "Event %d (%d s elapsed / %d s left) -> ETA: %s" %
	# (n, elapsedtime, remainingtime, eta))


	## Do some checks on output histo file, before we stat the event loop
	histo_parentdir = os.path.dirname(os.path.abspath(opts.HISTOFILE))
	if not os.path.exists(histo_parentdir):
	logging.error('Parent path of output histogram file does not exist: %s\nExiting.' % histo_parentdir)
	sys.exit(4)
	if not os.access(histo_parentdir,os.W_OK):
	logging.error('Insufficient permissions to write output histogram file to directory %s\nExiting.' % histo_parentdir)
	sys.exit(4)


	## Set up analysis handler
	RUNNAME = opts.RUN_NAME or ""
	ah = rivet.AnalysisHandler(RUNNAME)
	ah.setIgnoreBeams(opts.IGNORE_BEAMS)
	+ah.skipMultiWeights(opts.SKIP_WEIGHTS)

	for a in opts.ANALYSES:
	## Print warning message and exit if not a valid analysis name
	if not a in all_analyses:
	logging.warning("'%s' is not a known Rivet analysis! Do you need to set RIVET_ANALYSIS_PATH or use the --pwd switch?\n" % a)
	# TODO: lay out more neatly, or even try for a "did you mean XXXX?" heuristic?
	logging.warning("There are %d currently available analyses:\n" % len(all_analyses) + ", ".join(all_analyses))
	sys.exit(1)
	logging.debug("Adding analysis '%s'" % a)
	ah.addAnalysis(a)

	if opts.PRELOADFILE is not None:
	ah.readData(opts.PRELOADFILE)


	if opts.SHOW_BIBTEX:
	bibs = []
	for aname in sorted(ah.analysisNames()):
	ana = rivet.AnalysisLoader.getAnalysis(aname)
	bibs.append("% " + aname + "\n" + ana.bibTeX())
	if bibs:
	print("\nBibTeX for used Rivet analyses:\n")
	print("% --------------------------\n")
	print("\n\n".join(bibs) + "\n")
	print("% --------------------------\n")


	## Read and process events
	run = rivet.Run(ah)
	if opts.CROSS_SECTION is not None:
	logging.info("User-supplied cross-section = %e pb" % opts.CROSS_SECTION)
	run.setCrossSection(opts.CROSS_SECTION)
	if opts.LIST_USED_ANALYSES is not None:
	run.setListAnalyses(opts.LIST_USED_ANALYSES)

	## Print platform type
	import platform
	starttime = datetime.datetime.now().replace(microsecond=0)
	logging.info("Rivet %s running on machine %s (%s) at %s" % \
	(rivet.version(), platform.node(), platform.machine(), str(starttime)))


	def min_nonnull(a, b):
	"A version of min which considers None to always be greater than a real number"
	if a is None: return b
	if b is None: return a
	return min(a, b)

	## Set up an event timeout handler
	class TimeoutException(Exception):
	pass
	if opts.EVENT_TIMEOUT or opts.RUN_TIMEOUT:
	def evttimeouthandler(signum, frame):
	logging.warn("It has taken more than %d secs to get an event! Is the input event stream working?" %
	min_nonnull(opts.EVENT_TIMEOUT, opts.RUN_TIMEOUT))
	raise TimeoutException("Event timeout")
	signal.signal(signal.SIGALRM, evttimeouthandler)


	## Init run based on one event
	hepmcfile = HEPMCFILES[0]
	## Apply a file-level weight derived from the filename
	hepmcfileweight = 1.0
	if ":" in hepmcfile:
	hepmcfile, hepmcfileweight = hepmcfile.rsplit(":", 1)
	hepmcfileweight = float(hepmcfileweight)
	try:
	if opts.EVENT_TIMEOUT or opts.RUN_TIMEOUT:
	signal.alarm(min_nonnull(opts.EVENT_TIMEOUT, opts.RUN_TIMEOUT))
	init_ok = run.init(hepmcfile, hepmcfileweight)
	signal.alarm(0)
	if not init_ok:
	logging.error("Failed to initialise using event file '%s'... exiting" % hepmcfile)
	sys.exit(2)
	except TimeoutException as te:
	logging.error("Timeout in initialisation from event file '%s'... exiting" % hepmcfile)
	sys.exit(3)
	except Exception as ex:
	logging.warning("Could not read from '%s' (error=%s)" % (hepmcfile, str(ex)))
	sys.exit(3)

	## Event loop
	evtnum = 0
	for fileidx, hepmcfile in enumerate(HEPMCFILES):
	## Apply a file-level weight derived from the filename
	hepmcfileweight = 1.0
	if ":" in hepmcfile:
	hepmcfile, hepmcfileweight = hepmcfile.rsplit(":", 1)
	hepmcfileweight = float(hepmcfileweight)

	## Open next HepMC file (NB. this doesn't apply to the first file: it was already used for the run init)
	if fileidx > 0:
	try:
	run.openFile(hepmcfile, hepmcfileweight)
	except Exception as ex:
	logging.warning("Could not read from '%s' (error=%s)" % (hepmcfile, ex))
	continue
	if not run.readEvent():
	logging.warning("Could not read events from '%s'" % hepmcfile)
	continue

	## Announce new file
	msg = "Reading events from '%s'" % hepmcfile
	if hepmcfileweight != 1.0:
	msg += " (file weight = %e)" % hepmcfileweight
	logging.info(msg)

	## The event loop
	while opts.MAXEVTNUM is None or evtnum-opts.EVTSKIPNUM < opts.MAXEVTNUM:
	evtnum += 1

	## Optional event skipping
	if evtnum <= opts.EVTSKIPNUM:
	logging.debug("Skipping event #%i" % evtnum)
	run.skipEvent();
	continue

	## Only log the event number once we're actually processing
	logNEvt(evtnum, starttime, opts.MAXEVTNUM)

	## Process this event
	processed_ok = run.processEvent()
	if not processed_ok:
	logging.warn("Event processing failed for evt #%i!" % evtnum)
	break

	## Set flag to exit event loop if run timeout exceeded
	if opts.RUN_TIMEOUT and (time.time() - starttime) > opts.RUN_TIMEOUT:
	logging.warning("Run timeout of %d secs exceeded... exiting gracefully" % opts.RUN_TIMEOUT)
	RECVD_KILL_SIGNAL = True

	## Exit the loop if signalled
	if RECVD_KILL_SIGNAL is not None:
	break

	## Read next event (with timeout handling if requested)
	try:
	if opts.EVENT_TIMEOUT:
	signal.alarm(opts.EVENT_TIMEOUT)
	read_ok = run.readEvent()
	signal.alarm(0)
	if not read_ok:
	break
	except TimeoutException as te:
	logging.error("Timeout in reading event from '%s'... exiting" % hepmcfile)
	sys.exit(3)

	## Write a histo file snapshot if appropriate
	if opts.HISTO_WRITE_INTERVAL is not None and opts.HISTO_WRITE_INTERVAL > 0:
	if evtnum % opts.HISTO_WRITE_INTERVAL == 0:
	ah.writeData(opts.HISTOFILE)


	## Print end-of-loop messages
	loopendtime = datetime.datetime.now().replace(microsecond=0)
	logging.info("Finished event loop at %s" % str(loopendtime))
	logging.info("Generator provided cross-section = %e pb" % ah.nominalCrossSection())


	## Finalize and write out data file
	run.finalize()

	logging.info("Applied cross-section = %e pb" % ah.nominalCrossSection())
	print()

	if opts.WRITE_DATA:
	ah.writeData(opts.HISTOFILE)
	print()
	endtime = datetime.datetime.now().replace(microsecond=0)
	logging.info("Rivet run completed at %s, time elapsed = %s" % (str(endtime), str(endtime-starttime)))
	print()
	logging.info("Histograms written to %s" % os.path.abspath(opts.HISTOFILE))
	diff --git a/include/Rivet/AnalysisHandler.hh b/include/Rivet/AnalysisHandler.hh
	--- a/include/Rivet/AnalysisHandler.hh
	+++ b/include/Rivet/AnalysisHandler.hh
	@@ -1,279 +1,285 @@
	// -- C++ --
	#ifndef RIVET_RivetHandler_HH
	#define RIVET_RivetHandler_HH

	#include "Rivet/Config/RivetCommon.hh"
	#include "Rivet/Particle.hh"
	#include "Rivet/AnalysisLoader.hh"
	#include "Rivet/Tools/RivetYODA.hh"

	namespace Rivet {


	// Forward declaration and smart pointer for Analysis
	class Analysis;
	typedef std::shared_ptr<Analysis> AnaHandle;


	// Needed to make smart pointers compare equivalent in the STL set
	struct CmpAnaHandle {
	bool operator() (const AnaHandle& a, const AnaHandle& b) const {
	return a.get() < b.get();
	}
	};


	/// A class which handles a number of analysis objects to be applied to
	/// generated events. An {@link Analysis}' AnalysisHandler is also responsible
	/// for handling the final writing-out of histograms.
	class AnalysisHandler {
	public:

	/// @name Constructors and destructors. */
	//@{

	/// Preferred constructor, with optional run name.
	AnalysisHandler(const string& runname="");

	/// @brief Destructor
	/// The destructor is not virtual, as this class should not be inherited from.
	~AnalysisHandler();

	//@}


	private:

	/// Get a logger object.
	Log& getLog() const;


	public:

	/// @name Run properties
	//@{

	/// Get the name of this run.
	string runName() const;


	/// Get the number of events seen. Should only really be used by external
	/// steering code or analyses in the finalize phase.
	size_t numEvents() const;

	/// Get the sum of the event weights seen - the weighted equivalent of the
	/// number of events. Should only really be used by external steering code
	/// or analyses in the finalize phase.
	double sumOfWeights() const {
	return _eventCounter->sumW();
	}

	const vector<string>& weightNames() const {
	return _weightNames;
	}

	size_t numWeights() const {
	return _weightNames.size();
	}

	/// Set the cross-section for the process being generated.
	AnalysisHandler& setCrossSection(double xs, double xserr);

	/// Get the cross-section known to the handler.
	Scatter1DPtr crossSection() const {
	return _xs;
	}

	double nominalCrossSection() const {
	_xs.get()->setActiveWeightIdx(_defaultWeightIdx);
	const YODA::Scatter1D::Points& ps = _xs->points();
	if (ps.size() != 1) {
	string errMsg = "cross section missing when requesting nominal cross section";
	throw Error(errMsg);
	}
	double xs = ps[0].x();
	_xs.get()->unsetActiveWeight();
	return xs;
	}


	/// Set the beam particles for this run
	AnalysisHandler& setRunBeams(const ParticlePair& beams) {
	_beams = beams;
	MSG_DEBUG("Setting run beams = " << beams << " @ " << sqrtS()/GeV << " GeV");
	return *this;
	}

	/// Get the beam particles for this run, usually determined from the first event.
	const ParticlePair& beams() const { return _beams; }

	/// Get beam IDs for this run, usually determined from the first event.
	/// @deprecated Use standalone beamIds(ah.beams()), to clean AH interface
	PdgIdPair beamIds() const;

	/// Get energy for this run, usually determined from the first event.
	/// @deprecated Use standalone sqrtS(ah.beams()), to clean AH interface
	double sqrtS() const;

	/// Setter for _ignoreBeams
	void setIgnoreBeams(bool ignore=true);

	+ /// Setter for _skipWeights
	+ void skipMultiWeights(bool ignore=false);
	+
	//@}


	/// @name Handle analyses
	//@{

	/// Get a list of the currently registered analyses' names.
	std::vector<std::string> analysisNames() const;

	/// Get the collection of currently registered analyses.
	const std::set<AnaHandle, CmpAnaHandle>& analyses() const {
	return _analyses;
	}

	/// Get a registered analysis by name.
	const AnaHandle analysis(const std::string& analysisname) const;


	/// Add an analysis to the run list by object
	AnalysisHandler& addAnalysis(Analysis* analysis);

	/// @brief Add an analysis to the run list using its name.
	///
	/// The actual Analysis to be used will be obtained via
	/// AnalysisLoader::getAnalysis(string). If no matching analysis is found,
	/// no analysis is added (i.e. the null pointer is checked and discarded.
	AnalysisHandler& addAnalysis(const std::string& analysisname);

	/// @brief Add analyses to the run list using their names.
	///
	/// The actual {@link Analysis}' to be used will be obtained via
	/// AnalysisHandler::addAnalysis(string), which in turn uses
	/// AnalysisLoader::getAnalysis(string). If no matching analysis is found
	/// for a given name, no analysis is added, but also no error is thrown.
	AnalysisHandler& addAnalyses(const std::vector<std::string>& analysisnames);


	/// Remove an analysis from the run list using its name.
	AnalysisHandler& removeAnalysis(const std::string& analysisname);

	/// Remove analyses from the run list using their names.
	AnalysisHandler& removeAnalyses(const std::vector<std::string>& analysisnames);

	//@}


	/// @name Main init/execute/finalise
	//@{

	/// Initialize a run, with the run beams taken from the example event.
	void init(const GenEvent& event);

	/// @brief Analyze the given \a event by reference.
	///
	/// This function will call the AnalysisBase::analyze() function of all
	/// included analysis objects.
	void analyze(const GenEvent& event);

	/// @brief Analyze the given \a event by pointer.
	///
	/// This function will call the AnalysisBase::analyze() function of all
	/// included analysis objects, after checking the event pointer validity.
	void analyze(const GenEvent* event);

	/// Finalize a run. This function calls the AnalysisBase::finalize()
	/// functions of all included analysis objects.
	void finalize();

	//@}


	/// @name Histogram / data object access
	//@{

	/// Add a vector of analysis objects to the current state.
	void addData(const std::vector<YODA::AnalysisObjectPtr>& aos);

	/// Read analysis plots into the histo collection (via addData) from the named file.
	void readData(const std::string& filename);

	/// Get all analyses' plots as a vector of analysis objects.
	std::vector<YODA::AnalysisObjectPtr> getData() const;
	std::vector<MultiweightAOPtr> getRivetAOs() const;
	std::vector<YODA::AnalysisObjectPtr> getYodaAOs() const;

	/// Get all analyses' plots as a vector of analysis objects.
	void setWeightNames(const GenEvent& ge);

	/// Do we have named weights?
	bool haveNamedWeights();


	/// Write all analyses' plots (via getData) to the named file.
	void writeData(const std::string& filename) const;

	//@}

	/// Indicate which Rivet stage we're in.
	/// At the moment, only INIT is used to enable booking.
	enum class Stage { OTHER, INIT };

	/// Which stage are we in?
	Stage stage() const { return _stage; }

	private:
	/// Current handler stage
	Stage _stage = Stage::OTHER;

	/// The collection of Analysis objects to be used.
	set<AnaHandle, CmpAnaHandle> _analyses;

	/// @name Run properties
	//@{

	/// Weight names
	std::vector<std::string> _weightNames;
	std::vector<std::valarray<double> > _subEventWeights;
	size_t _numWeightTypes; // always == WeightVector.size()

	/// Run name
	std::string _runname;

	mutable CounterPtr _eventCounter;

	/// Cross-section known to AH
	Scatter1DPtr _xs;

	/// Beams used by this run.
	ParticlePair _beams;

	/// Flag to check if init has been called
	bool _initialised;

	/// Flag whether input event beams should be ignored in compatibility check
	bool _ignoreBeams;

	+ /// Flag to check if multiweights should be included
	+ bool _skipWeights;
	+
	/// Current event number
	int _eventNumber;

	size_t _defaultWeightIdx;
	//@}


	private:

	/// The assignment operator is private and must never be called.
	/// In fact, it should not even be implemented.
	AnalysisHandler& operator=(const AnalysisHandler&);

	/// The copy constructor is private and must never be called. In
	/// fact, it should not even be implemented.
	AnalysisHandler(const AnalysisHandler&);

	};


	}

	#endif
	diff --git a/pyext/rivet/core.pyx b/pyext/rivet/core.pyx
	--- a/pyext/rivet/core.pyx
	+++ b/pyext/rivet/core.pyx
	@@ -1,232 +1,235 @@
	# distutils: language = c++

	cimport rivet as c
	from cython.operator cimport dereference as deref
	# Need to be careful with memory management -- perhaps use the base object that
	# we used in YODA?

	cdef extern from "<utility>" namespace "std" nogil:
	cdef c.unique_ptr[c.Analysis] move(c.unique_ptr[c.Analysis])

	cdef class AnalysisHandler:
	cdef c.AnalysisHandler *_ptr

	def __cinit__(self):
	self._ptr = new c.AnalysisHandler()

	def __del__(self):
	del self._ptr

	def setIgnoreBeams(self, ignore=True):
	self._ptr.setIgnoreBeams(ignore)

	+ def skipMultiWeights(self, ignore=False):
	+ self._ptr.skipMultiWeights(ignore)
	+
	def addAnalysis(self, name):
	self._ptr.addAnalysis(name.encode('utf-8'))
	return self

	def analysisNames(self):
	anames = self._ptr.analysisNames()
	return [ a.decode('utf-8') for a in anames ]

	# def analysis(self, aname):
	# cdef c.Analysis* ptr = self._ptr.analysis(aname)
	# cdef Analysis pyobj = Analysis.__new__(Analysis)
	# if not ptr:
	# return None
	# pyobj._ptr = ptr
	# return pyobj

	def readData(self, name):
	self._ptr.readData(name.encode('utf-8'))

	def writeData(self, name):
	self._ptr.writeData(name.encode('utf-8'))

	def nominalCrossSection(self):
	return self._ptr.nominalCrossSection()

	def finalize(self):
	self._ptr.finalize()


	cdef class Run:
	cdef c.Run *_ptr

	def __cinit__(self, AnalysisHandler h):
	self._ptr = new c.Run(h._ptr[0])

	def __del__(self):
	del self._ptr

	def setCrossSection(self, double x):
	self._ptr.setCrossSection(x)
	return self

	def setListAnalyses(self, choice):
	self._ptr.setListAnalyses(choice)
	return self

	def init(self, name, weight=1.0):
	return self._ptr.init(name.encode('utf-8'), weight)

	def openFile(self, name, weight=1.0):
	return self._ptr.openFile(name.encode('utf-8'), weight)

	def readEvent(self):
	return self._ptr.readEvent()

	def skipEvent(self):
	return self._ptr.skipEvent()

	def processEvent(self):
	return self._ptr.processEvent()

	def finalize(self):
	return self._ptr.finalize()


	cdef class Analysis:
	cdef c.unique_ptr[c.Analysis] _ptr

	def __init__(self):
	raise RuntimeError('This class cannot be instantiated')

	def requiredBeams(self):
	return deref(self._ptr).requiredBeams()

	def requiredEnergies(self):
	return deref(self._ptr).requiredEnergies()

	def keywords(self):
	kws = deref(self._ptr).keywords()
	return [ k.decode('utf-8') for k in kws ]

	def authors(self):
	auths = deref(self._ptr).authors()
	return [ a.decode('utf-8') for a in auths ]

	def bibKey(self):
	return deref(self._ptr).bibKey().decode('utf-8')

	def name(self):
	return deref(self._ptr).name().decode('utf-8')

	def bibTeX(self):
	return deref(self._ptr).bibTeX().decode('utf-8')

	def references(self):
	refs = deref(self._ptr).references()
	return [ r.decode('utf-8') for r in refs ]

	def collider(self):
	return deref(self._ptr).collider().decode('utf-8')

	def description(self):
	return deref(self._ptr).description().decode('utf-8')

	def experiment(self):
	return deref(self._ptr).experiment().decode('utf-8')

	def inspireId(self):
	return deref(self._ptr).inspireId().decode('utf-8')

	def spiresId(self):
	return deref(self._ptr).spiresId().decode('utf-8')

	def runInfo(self):
	return deref(self._ptr).runInfo().decode('utf-8')

	def status(self):
	return deref(self._ptr).status().decode('utf-8')

	def summary(self):
	return deref(self._ptr).summary().decode('utf-8')

	def year(self):
	return deref(self._ptr).year().decode('utf-8')

	def luminosityfb(self):
	return deref(self._ptr).luminosityfb().decode('utf-8')

	#cdef object
	LEVELS = dict(TRACE = 0, DEBUG = 10, INFO = 20,
	WARN = 30, WARNING = 30, ERROR = 40,
	CRITICAL = 50, ALWAYS = 50)


	cdef class AnalysisLoader:
	@staticmethod
	def analysisNames():
	names = c.AnalysisLoader_analysisNames()
	return [ n.decode('utf-8') for n in names ]


	@staticmethod
	def getAnalysis(name):
	name = name.encode('utf-8')
	cdef c.unique_ptr[c.Analysis] ptr = c.AnalysisLoader_getAnalysis(name)
	cdef Analysis pyobj = Analysis.__new__(Analysis)
	if not ptr:
	return None
	pyobj._ptr = move(ptr)
	# Create python object
	return pyobj


	def getAnalysisLibPaths():
	ps = c.getAnalysisLibPaths()
	return [ p.decode('utf-8') for p in ps ]

	def setAnalysisLibPaths(xs):
	bs = [ x.encode('utf-8') for x in xs ]
	c.setAnalysisLibPaths(bs)

	def addAnalysisLibPath(path):
	c.addAnalysisLibPath(path.encode('utf-8'))


	def setAnalysisDataPaths(xs):
	bs = [ x.encode('utf-8') for x in xs ]
	c.setAnalysisDataPaths(bs)

	def addAnalysisDataPath(path):
	c.addAnalysisDataPath(path.encode('utf-8'))

	def getAnalysisDataPaths():
	ps = c.getAnalysisDataPaths()
	return [ p.decode('utf-8') for p in ps ]

	def findAnalysisDataFile(q):
	f = c.findAnalysisDataFile(q.encode('utf-8'))
	return f.decode('utf-8')

	def getAnalysisRefPaths():
	ps = c.getAnalysisRefPaths()
	return [ p.decode('utf-8') for p in ps ]

	def findAnalysisRefFile(q):
	f = c.findAnalysisRefFile(q.encode('utf-8'))
	return f.decode('utf-8')


	def getAnalysisInfoPaths():
	ps = c.getAnalysisInfoPaths()
	return [ p.decode('utf-8') for p in ps ]

	def findAnalysisInfoFile(q):
	f = c.findAnalysisInfoFile(q.encode('utf-8'))
	return f.decode('utf-8')

	def getAnalysisPlotPaths():
	ps = c.getAnalysisPlotPaths()
	return [ p.decode('utf-8') for p in ps ]

	def findAnalysisPlotFile(q):
	f = c.findAnalysisPlotFile(q.encode('utf-8'))
	return f.decode('utf-8')

	def version():
	return c.version().decode('utf-8')

	def setLogLevel(name, level):
	c.setLogLevel(name.encode('utf-8'), level)
	diff --git a/pyext/rivet/rivet.pxd b/pyext/rivet/rivet.pxd
	--- a/pyext/rivet/rivet.pxd
	+++ b/pyext/rivet/rivet.pxd
	@@ -1,92 +1,93 @@
	from libcpp.map cimport map
	from libcpp.pair cimport pair
	from libcpp.vector cimport vector
	from libcpp cimport bool
	from libcpp.string cimport string
	from libcpp.memory cimport unique_ptr

	ctypedef int PdgId
	ctypedef pair[PdgId,PdgId] PdgIdPair

	cdef extern from "Rivet/AnalysisHandler.hh" namespace "Rivet":
	cdef cppclass AnalysisHandler:
	void setIgnoreBeams(bool)
	+ void skipMultiWeights(bool)
	AnalysisHandler& addAnalysis(string)
	vector[string] analysisNames() const
	# Analysis* analysis(string)
	void writeData(string&)
	void readData(string&)
	double nominalCrossSection()
	void finalize()

	cdef extern from "Rivet/Run.hh" namespace "Rivet":
	cdef cppclass Run:
	Run(AnalysisHandler)
	Run& setCrossSection(double) # For chaining?
	Run& setListAnalyses(bool)
	bool init(string, double) except + # $2=1.0
	bool openFile(string, double) except + # $2=1.0
	bool readEvent() except +
	bool skipEvent() except +
	bool processEvent() except +
	bool finalize() except +

	cdef extern from "Rivet/Analysis.hh" namespace "Rivet":
	cdef cppclass Analysis:
	vector[PdgIdPair]& requiredBeams()
	vector[pair[double, double]] requiredEnergies()
	vector[string] authors()
	vector[string] references()
	vector[string] keywords()
	string name()
	string bibTeX()
	string bibKey()
	string collider()
	string description()
	string experiment()
	string inspireId()
	string spiresId()
	string runInfo()
	string status()
	string summary()
	string year()
	string luminosityfb()

	# Might need to translate the following errors, although I believe 'what' is now
	# preserved. But often, we need the exception class name.
	#Error
	#RangeError
	#LogicError
	#PidError
	#InfoError
	#WeightError
	#UserError

	cdef extern from "Rivet/AnalysisLoader.hh":
	vector[string] AnalysisLoader_analysisNames "Rivet::AnalysisLoader::analysisNames" ()
	unique_ptr[Analysis] AnalysisLoader_getAnalysis "Rivet::AnalysisLoader::getAnalysis" (string)

	cdef extern from "Rivet/Tools/RivetPaths.hh" namespace "Rivet":
	vector[string] getAnalysisLibPaths()
	void setAnalysisLibPaths(vector[string])
	void addAnalysisLibPath(string)

	vector[string] getAnalysisDataPaths()
	void setAnalysisDataPaths(vector[string])
	void addAnalysisDataPath(string)
	string findAnalysisDataFile(string)

	vector[string] getAnalysisRefPaths()
	string findAnalysisRefFile(string)

	vector[string] getAnalysisInfoPaths()
	string findAnalysisInfoFile(string)

	vector[string] getAnalysisPlotPaths()
	string findAnalysisPlotFile(string)

	cdef extern from "Rivet/Rivet.hh" namespace "Rivet":
	string version()

	cdef extern from "Rivet/Tools/Logging.hh":
	void setLogLevel "Rivet::Log::setLevel" (string, int)
	diff --git a/src/Core/AnalysisHandler.cc b/src/Core/AnalysisHandler.cc
	--- a/src/Core/AnalysisHandler.cc
	+++ b/src/Core/AnalysisHandler.cc
	@@ -1,520 +1,536 @@
	// -- C++ --
	#include "Rivet/Config/RivetCommon.hh"
	#include "Rivet/AnalysisHandler.hh"
	#include "Rivet/Analysis.hh"
	#include "Rivet/Tools/ParticleName.hh"
	#include "Rivet/Tools/BeamConstraint.hh"
	#include "Rivet/Tools/Logging.hh"
	#include "Rivet/Projections/Beam.hh"
	#include "YODA/ReaderYODA.h"
	#include "YODA/WriterYODA.h"
	#include <regex>
	#include <iostream>

	using std::cout;
	using std::cerr;

	namespace {
	inline std::vector<std::string> split(const std::string& input, const std::string& regex) {
	// passing -1 as the submatch index parameter performs splitting
	std::regex re(regex);
	std::sregex_token_iterator
	first{input.begin(), input.end(), re, -1},
	last;
	return {first, last};
	}
	}


	namespace Rivet {


	AnalysisHandler::AnalysisHandler(const string& runname)
	: _runname(runname),
	_eventCounter(vector<string>(), Counter()), _xs(vector<string>(), Scatter1D()),
	- _initialised(false), _ignoreBeams(false),
	+ _initialised(false), _ignoreBeams(false), _skipWeights(false),
	_defaultWeightIdx(0)
	{}


	AnalysisHandler::~AnalysisHandler()
	{}


	Log& AnalysisHandler::getLog() const {
	return Log::getLog("Rivet.Analysis.Handler");
	}


	/// http://stackoverflow.com/questions/4654636/how-to-determine-if-a-string-is-a-number-with-c
	bool is_number(const std::string& s)
	{
	std::string::const_iterator it = s.begin();
	while (it != s.end() && std::isdigit(*it)) ++it;
	return !s.empty() && it == s.end();
	}

	/// Check if any of the weightnames is not a number
	bool AnalysisHandler::haveNamedWeights() {
	bool dec=false;
	for (unsigned int i=0;i<_weightNames.size();++i) {
	string s = _weightNames[i];
	if (!is_number(s)) {
	dec=true;
	break;
	}
	}
	return dec;
	}


	void AnalysisHandler::init(const GenEvent& ge) {
	if (_initialised)
	throw UserError("AnalysisHandler::init has already been called: cannot re-initialize!");

	// TODO TODO
	// should the rivet analysis objects know about weight names?

	setRunBeams(Rivet::beams(ge));
	MSG_DEBUG("Initialising the analysis handler");
	_eventNumber = ge.event_number();

	setWeightNames(ge);
	if (haveNamedWeights())
	MSG_INFO("Using named weights");
	else
	MSG_INFO("NOT using named weights. Using first weight as nominal weight");

	_eventCounter = CounterPtr(weightNames(), Counter("_EVTCOUNT"));

	// set the cross section based on what is reported by this event.
	// if no cross section
	if (ge.cross_section()) {
	MSG_TRACE("getting cross section.");
	double xs = ge.cross_section()->cross_section();
	+ MSG_INFO("get xserr");
	double xserr = ge.cross_section()->cross_section_error();
	+ MSG_INFO("set xsec");
	setCrossSection(xs, xserr);
	+ MSG_INFO("set xsec done");
	}

	+ MSG_INFO("check anas size");
	// Check that analyses are beam-compatible, and remove those that aren't
	const size_t num_anas_requested = analysisNames().size();
	vector<string> anamestodelete;
	for (const AnaHandle a : _analyses) {
	if (!_ignoreBeams && !a->isCompatible(beams())) {
	//MSG_DEBUG(a->name() << " requires beams " << a->requiredBeams() << " @ " << a->requiredEnergies() << " GeV");
	anamestodelete.push_back(a->name());
	}
	}
	for (const string& aname : anamestodelete) {
	MSG_WARNING("Analysis '" << aname << "' is incompatible with the provided beams: removing");
	removeAnalysis(aname);
	}
	if (num_anas_requested > 0 && analysisNames().empty()) {
	cerr << "All analyses were incompatible with the first event's beams\n"
	<< "Exiting, since this probably wasn't intentional!" << '\n';
	exit(1);
	}

	// Warn if any analysis' status is not unblemished
	for (const AnaHandle a : analyses()) {
	if (toUpper(a->status()) == "PRELIMINARY") {
	MSG_WARNING("Analysis '" << a->name() << "' is preliminary: be careful, it may change and/or be renamed!");
	} else if (toUpper(a->status()) == "OBSOLETE") {
	MSG_WARNING("Analysis '" << a->name() << "' is obsolete: please update!");
	} else if (toUpper(a->status()).find("UNVALIDATED") != string::npos) {
	MSG_WARNING("Analysis '" << a->name() << "' is unvalidated: be careful, it may be broken!");
	}
	}

	// Initialize the remaining analyses
	_stage = Stage::INIT;
	for (AnaHandle a : _analyses) {
	MSG_DEBUG("Initialising analysis: " << a->name());
	try {
	// Allow projection registration in the init phase onwards
	a->_allowProjReg = true;
	a->init();
	//MSG_DEBUG("Checking consistency of analysis: " << a->name());
	//a->checkConsistency();
	} catch (const Error& err) {
	cerr << "Error in " << a->name() << "::init method: " << err.what() << '\n';
	exit(1);
	}
	MSG_DEBUG("Done initialising analysis: " << a->name());
	}
	_stage = Stage::OTHER;
	_initialised = true;
	MSG_DEBUG("Analysis handler initialised");
	}

	void AnalysisHandler::setWeightNames(const GenEvent& ge) {
	/// reroute the print output to a std::stringstream and process
	/// The iteration is done over a map in hepmc2 so this is safe
	std::ostringstream stream;
	ge.weights().print(stream); // Super lame, I know
	string str = stream.str();

	std::regex re("(([^()]+))"); // Regex for stuff enclosed by parentheses ()
	size_t idx = 0;
	for(std::sregex_iterator i = std::sregex_iterator(str.begin(), str.end(), re);
	i != std::sregex_iterator(); ++i ) {
	std::smatch m = *i;
	vector<string> temp = ::split(m.str(), "[,]");
	if (temp.size() ==2) {
	MSG_DEBUG("Name of weight #" << _weightNames.size() << ": " << temp[0]);

	// store the default weight based on weight names
	if (temp[0] == "Weight" \|\| temp[0] == "0" \|\| temp[0] == "Default") {
	MSG_DEBUG(_weightNames.size() << " is being used as the nominal.");
	_weightNames.push_back("");
	- _defaultWeightIdx = idx;
	- } else
	+ _defaultWeightIdx = _skipWeights? 0 : idx;
	+ } else if (!_skipWeights) {
	_weightNames.push_back(temp[0]);
	+ }


	idx++;
	}
	}
	}


	void AnalysisHandler::analyze(const GenEvent& ge) {
	// Call init with event as template if not already initialised
	if (!_initialised) init(ge);
	assert(_initialised);

	// Ensure that beam details match those from the first event (if we're checking beams)
	if ( !_ignoreBeams ) {
	const PdgIdPair beams = Rivet::beamIds(ge);
	const double sqrts = Rivet::sqrtS(ge);
	if (!compatible(beams, _beams) \|\| !fuzzyEquals(sqrts, sqrtS())) {
	cerr << "Event beams mismatch: "
	<< PID::toBeamsString(beams) << " @ " << sqrts/GeV << " GeV" << " vs. first beams "
	<< this->beams() << " @ " << this->sqrtS()/GeV << " GeV" << '\n';
	exit(1);
	}
	}


	// Create the Rivet event wrapper
	/// @todo Filter/normalize the event here
	Event event(ge);

	// set the cross section based on what is reported by this event.
	// if no cross section
	MSG_TRACE("getting cross section.");
	if (ge.cross_section()) {
	MSG_TRACE("getting cross section from GenEvent.");
	double xs = ge.cross_section()->cross_section();
	double xserr = ge.cross_section()->cross_section_error();
	setCrossSection(xs, xserr);
	}

	// won't happen for first event because _eventNumber is set in
	// init()
	if (_eventNumber != ge.event_number()) {
	/// @todo
	/// can we get away with not passing a matrix?

	MSG_TRACE("AnalysisHandler::analyze(): Pushing _eventCounter to persistent.");
	_eventCounter.get()->pushToPersistent(_subEventWeights);
	// if this is indeed a new event, push the temporary
	// histograms and reset
	for (const AnaHandle& a : _analyses) {
	for (auto ao : a->analysisObjects()) {
	MSG_TRACE("AnalysisHandler::analyze(): Pushing " << a->name() << "'s " << ao->name() << " to persistent.");
	ao.get()->pushToPersistent(_subEventWeights);
	}
	MSG_TRACE("AnalysisHandler::analyze(): finished pushing " << a->name() << "'s objects to persistent.");
	}

	_eventNumber = ge.event_number();

	MSG_DEBUG("nominal event # " << _eventCounter.get()->_persistent[0]->numEntries());
	MSG_DEBUG("nominal sum of weights: " << _eventCounter.get()->_persistent[0]->sumW());
	MSG_DEBUG("Event has " << _subEventWeights.size() << " sub events.");
	_subEventWeights.clear();
	}


	MSG_TRACE("starting new sub event");
	_eventCounter.get()->newSubEvent();

	for (const AnaHandle& a : _analyses) {
	for (auto ao : a->analysisObjects()) {
	ao.get()->newSubEvent();
	}
	}

	_subEventWeights.push_back(event.weights());
	MSG_DEBUG("Analyzing subevent #" << _subEventWeights.size() - 1 << ".");

	_eventCounter->fill();
	// Run the analyses
	for (AnaHandle a : _analyses) {
	MSG_TRACE("About to run analysis " << a->name());
	try {
	a->analyze(event);
	} catch (const Error& err) {
	cerr << "Error in " << a->name() << "::analyze method: " << err.what() << '\n';
	exit(1);
	}
	MSG_TRACE("Finished running analysis " << a->name());
	}
	}


	void AnalysisHandler::analyze(const GenEvent* ge) {
	if (ge == NULL) {
	MSG_ERROR("AnalysisHandler received null pointer to GenEvent");
	//throw Error("AnalysisHandler received null pointer to GenEvent");
	}
	analyze(*ge);
	}


	void AnalysisHandler::finalize() {
	if (!_initialised) return;
	MSG_INFO("Finalising analyses");

	MSG_TRACE("AnalysisHandler::finalize(): Pushing analysis objects to persistent.");
	_eventCounter.get()->pushToPersistent(_subEventWeights);
	for (const AnaHandle& a : _analyses) {
	for (auto ao : a->analysisObjects())
	ao.get()->pushToPersistent(_subEventWeights);
	}

	for (const AnaHandle& a : _analyses) {
	for (size_t iW = 0; iW < numWeights(); iW++) {
	_eventCounter.get()->setActiveWeightIdx(iW);
	_xs.get()->setActiveWeightIdx(iW);
	for (auto ao : a->analysisObjects())
	ao.get()->setActiveWeightIdx(iW);

	MSG_TRACE("running " << a->name() << "::finalize() for weight " << iW << ".");

	try {
	a->finalize();
	} catch (const Error& err) {
	cerr << "Error in " << a->name() << "::finalize method: " << err.what() << '\n';
	exit(1);
	}
	}
	}

	// Print out number of events processed
	MSG_INFO("Processed " << numEvents() << " event" << (numEvents() == 1 ? "" : "s"));

	// Print out MCnet boilerplate
	cout << '\n';
	cout << "The MCnet usage guidelines apply to Rivet: see http://www.montecarlonet.org/GUIDELINES" << '\n';
	cout << "Please acknowledge plots made with Rivet analyses, and cite arXiv:1003.0694 (http://arxiv.org/abs/1003.0694)" << '\n';
	}


	AnalysisHandler& AnalysisHandler::addAnalysis(const string& analysisname) {
	// Check for a duplicate analysis
	/// @todo Might we want to be able to run an analysis twice, with different params?
	/// Requires avoiding histo tree clashes, i.e. storing the histos on the analysis objects.
	for (const AnaHandle& a : _analyses) {
	if (a->name() == analysisname) {
	MSG_WARNING("Analysis '" << analysisname << "' already registered: skipping duplicate");
	return *this;
	}
	}
	AnaHandle analysis( AnalysisLoader::getAnalysis(analysisname) );
	if (analysis.get() != 0) { // < Check for null analysis.
	MSG_DEBUG("Adding analysis '" << analysisname << "'");
	analysis->_analysishandler = this;
	_analyses.insert(analysis);
	} else {
	MSG_WARNING("Analysis '" << analysisname << "' not found.");
	}
	// MSG_WARNING(_analyses.size());
	// for (const AnaHandle& a : _analyses) MSG_WARNING(a->name());
	return *this;
	}


	AnalysisHandler& AnalysisHandler::removeAnalysis(const string& analysisname) {
	std::shared_ptr<Analysis> toremove;
	for (const AnaHandle a : _analyses) {
	if (a->name() == analysisname) {
	toremove = a;
	break;
	}
	}
	if (toremove.get() != 0) {
	MSG_DEBUG("Removing analysis '" << analysisname << "'");
	_analyses.erase(toremove);
	}
	return *this;
	}


	void AnalysisHandler::addData(const std::vector<YODA::AnalysisObjectPtr>& aos) {
	for (const YODA::AnalysisObjectPtr ao : aos) {
	const string path = ao->path();
	if (path.size() > 1) { // path > "/"
	try {
	const string ananame = ::split(path, "/")[0];
	AnaHandle a = analysis(ananame);
	//MultiweightAOPtr mao = ????; /// @todo generate right Multiweight object from ao
	//a->addAnalysisObject(mao); /// @todo Need to statistically merge...
	} catch (const Error& e) {
	MSG_WARNING(e.what());
	}
	}
	}
	}

	void AnalysisHandler::readData(const string& filename) {
	vector<YODA::AnalysisObjectPtr> aos;
	try {
	/// @todo Use new YODA SFINAE to fill the smart ptr vector directly
	vector<YODA::AnalysisObject*> aos_raw;
	YODA::ReaderYODA::read(filename, aos_raw);
	for (YODA::AnalysisObject* aor : aos_raw) aos.push_back(YODA::AnalysisObjectPtr(aor));
	} catch (const YODA::ReadError & e) {
	throw UserError("Unexpected error in reading file: " + filename);
	}
	if (!aos.empty()) addData(aos);
	}


	vector<MultiweightAOPtr> AnalysisHandler::getRivetAOs() const {
	vector<MultiweightAOPtr> rtn;

	for (AnaHandle a : _analyses) {
	for (const auto & ao : a->analysisObjects()) {
	rtn.push_back(ao);
	}
	}

	rtn.push_back(_eventCounter);
	rtn.push_back(_xs);

	return rtn;
	}

	vector<YODA::AnalysisObjectPtr> AnalysisHandler::getYodaAOs() const {
	vector<YODA::AnalysisObjectPtr> rtn;

	for (auto rao : getRivetAOs()) {
	// need to set the index
	// before we can search the PATH
	rao.get()->setActiveWeightIdx(_defaultWeightIdx);
	if (rao->path().find("/TMP/") != string::npos)
	continue;

	for (size_t iW = 0; iW < numWeights(); iW++) {
	rao.get()->setActiveWeightIdx(iW);
	rtn.push_back(rao.get()->activeYODAPtr());
	}
	}

	sort(rtn.begin(), rtn.end(),
	[](YODA::AnalysisObjectPtr a, YODA::AnalysisObjectPtr b) {
	return a->path() < b->path();
	}
	);

	return rtn;
	}

	vector<YODA::AnalysisObjectPtr> AnalysisHandler::getData() const {

	return getYodaAOs();
	}


	void AnalysisHandler::writeData(const string& filename) const {
	const vector<YODA::AnalysisObjectPtr> aos = getData();
	try {
	YODA::WriterYODA::write(filename, aos.begin(), aos.end());
	} catch ( YODA::WriteError ) {
	throw UserError("Unexpected error in writing file: " + filename);
	}
	}


	string AnalysisHandler::runName() const { return _runname; }
	size_t AnalysisHandler::numEvents() const { return _eventCounter->numEntries(); }


	std::vector<std::string> AnalysisHandler::analysisNames() const {
	std::vector<std::string> rtn;
	for (AnaHandle a : _analyses) {
	rtn.push_back(a->name());
	}
	return rtn;
	}


	const AnaHandle AnalysisHandler::analysis(const std::string& analysisname) const {
	for (const AnaHandle a : analyses())
	if (a->name() == analysisname) return a;
	throw Error("No analysis named '" + analysisname + "' registered in AnalysisHandler");
	}


	AnalysisHandler& AnalysisHandler::addAnalyses(const std::vector<std::string>& analysisnames) {
	for (const string& aname : analysisnames) {
	//MSG_DEBUG("Adding analysis '" << aname << "'");
	addAnalysis(aname);
	}
	return *this;
	}


	AnalysisHandler& AnalysisHandler::removeAnalyses(const std::vector<std::string>& analysisnames) {
	for (const string& aname : analysisnames) {
	removeAnalysis(aname);
	}
	return *this;
	}


	AnalysisHandler& AnalysisHandler::setCrossSection(double xs, double xserr) {
	+ MSG_TRACE("init _xs");
	_xs = Scatter1DPtr(weightNames(), Scatter1D("_XSEC"));
	+ MSG_TRACE("init _evtC");
	_eventCounter.get()->setActiveWeightIdx(_defaultWeightIdx);
	+ MSG_TRACE("get sow");
	double nomwgt = sumOfWeights();

	// The cross section of each weight variation is the nominal cross section
	// times the sumOfWeights(variation) / sumOfWeights(nominal).
	// This way the cross section will work correctly
	+ MSG_TRACE("loop over weights");
	for (size_t iW = 0; iW < numWeights(); iW++) {
	+ MSG_TRACE("... weight " << iW);
	_eventCounter.get()->setActiveWeightIdx(iW);
	double s = sumOfWeights() / nomwgt;
	_xs.get()->setActiveWeightIdx(iW);
	_xs->addPoint(xss, xserrs);
	}

	+ MSG_TRACE("unset");
	_eventCounter.get()->unsetActiveWeight();
	_xs.get()->unsetActiveWeight();
	+ MSG_TRACE("done... return");
	return *this;
	}


	AnalysisHandler& AnalysisHandler::addAnalysis(Analysis* analysis) {
	analysis->_analysishandler = this;
	_analyses.insert(AnaHandle(analysis));
	return *this;
	}


	PdgIdPair AnalysisHandler::beamIds() const {
	return Rivet::beamIds(beams());
	}


	double AnalysisHandler::sqrtS() const {
	return Rivet::sqrtS(beams());
	}

	void AnalysisHandler::setIgnoreBeams(bool ignore) {
	_ignoreBeams=ignore;
	}

	+ void AnalysisHandler::skipMultiWeights(bool ignore) {
	+ _skipWeights=ignore;
	+ }
	+

	}
	diff --git a/src/Projections/InvMassFinalState.cc b/src/Projections/InvMassFinalState.cc
	--- a/src/Projections/InvMassFinalState.cc
	+++ b/src/Projections/InvMassFinalState.cc
	@@ -1,185 +1,185 @@
	// -- C++ --
	#include "Rivet/Projections/InvMassFinalState.hh"

	namespace Rivet {


	InvMassFinalState::InvMassFinalState(const FinalState& fsp,
	const pair<PdgId, PdgId>& idpair, // pair of decay products
	double minmass, // min inv mass
	double maxmass, // max inv mass
	double masstarget)
	: _minmass(minmass), _maxmass(maxmass), _masstarget(masstarget), _useTransverseMass(false)
	{
	setName("InvMassFinalState");
	declare(fsp, "FS");
	_decayids.push_back(idpair);
	}


	InvMassFinalState::InvMassFinalState(const FinalState& fsp,
	const vector<pair<PdgId, PdgId> >& idpairs, // vector of pairs of decay products
	double minmass, // min inv mass
	double maxmass, // max inv mass
	double masstarget)
	: _decayids(idpairs), _minmass(minmass), _maxmass(maxmass), _masstarget(masstarget), _useTransverseMass(false)
	{
	setName("InvMassFinalState");
	declare(fsp, "FS");
	}


	InvMassFinalState::InvMassFinalState(const pair<PdgId, PdgId>& idpair, // pair of decay products
	double minmass, // min inv mass
	double maxmass, // max inv mass
	double masstarget)
	: _minmass(minmass), _maxmass(maxmass), _masstarget(masstarget), _useTransverseMass(false)
	{
	setName("InvMassFinalState");
	_decayids.push_back(idpair);
	}


	InvMassFinalState::InvMassFinalState(const vector<pair<PdgId, PdgId> >& idpairs, // vector of pairs of decay products
	double minmass, // min inv mass
	double maxmass, // max inv mass
	double masstarget)
	: _decayids(idpairs), _minmass(minmass), _maxmass(maxmass), _masstarget(masstarget), _useTransverseMass(false)
	{
	setName("InvMassFinalState");
	}


	CmpState InvMassFinalState::compare(const Projection& p) const {
	// First compare the final states we are running on
	CmpState fscmp = mkNamedPCmp(p, "FS");
	if (fscmp != CmpState::EQ) return fscmp;

	// Then compare the two as final states
	const InvMassFinalState& other = dynamic_cast<const InvMassFinalState&>(p);
	fscmp = FinalState::compare(other);
	if (fscmp != CmpState::EQ) return fscmp;

	// Compare the mass limits
	CmpState masstypecmp = cmp(_useTransverseMass, other._useTransverseMass);
	if (masstypecmp != CmpState::EQ) return masstypecmp;
	CmpState massllimcmp = cmp(_minmass, other._minmass);
	if (massllimcmp != CmpState::EQ) return massllimcmp;
	CmpState masshlimcmp = cmp(_maxmass, other._maxmass);
	if (masshlimcmp != CmpState::EQ) return masshlimcmp;

	// Compare the decay species
	CmpState decaycmp = cmp(_decayids, other._decayids);
	if (decaycmp != CmpState::EQ) return decaycmp;

	// Finally compare them as final states
	return FinalState::compare(other);
	}


	void InvMassFinalState::project(const Event& e) {
	const FinalState& fs = applyProjection<FinalState>(e, "FS");
	calc(fs.particles());
	}


	void InvMassFinalState::calc(const Particles& inparticles) {
	_theParticles.clear();
	_particlePairs.clear();

	// Containers for the particles of type specified in the pair
	vector<const Particle*> type1, type2;
	// Get all the particles of the type specified in the pair from the particle list
	for (const Particle& ipart : inparticles) {
	// Loop around possible particle pairs
	for (const PdgIdPair& ipair : _decayids) {
	if (ipart.pid() == ipair.first) {
	if (accept(ipart)) type1 += &ipart;
	} else if (ipart.pid() == ipair.second) {
	if (accept(ipart)) type2 += &ipart;
	}
	}
	}
	if (type1.empty() \|\| type2.empty()) return;

	// Temporary container of selected particles iterators
	// Useful to compare iterators and avoid double occurrences of the same
	// particle in case it matches with more than another particle
	vector<const Particle*> tmp;

	// Now calculate the inv mass
	pair<double, pair<Particle, Particle> > closestPair;
	closestPair.first = 1e30;
	for (const Particle* i1 : type1) {
	for (const Particle* i2 : type2) {
	// Check this is actually a pair
	// (if more than one pair in vector particles can be unrelated)
	bool found = false;
	for (const PdgIdPair& ipair : _decayids) {
	if (i1->pid() == ipair.first && i2->pid() == ipair.second) {
	found = true;
	break;
	}
	}
	if (!found) continue;

	FourMomentum v4 = i1->momentum() + i2->momentum();
	if (v4.mass2() < 0) {
	MSG_DEBUG("Constructed negative inv mass2: skipping!");
	continue;
	}
	bool passedMassCut = false;
	if (_useTransverseMass) {
	passedMassCut = inRange(mT(i1->momentum(), i2->momentum()), _minmass, _maxmass);
	} else {
	passedMassCut = inRange(v4.mass(), _minmass, _maxmass);
	}

	if (passedMassCut) {
	MSG_DEBUG("Selecting particles with IDs " << i1->pid() << " & " << i2->pid()
	<< " and mass = " << v4.mass()/GeV << " GeV");
	// Store accepted particles, avoiding duplicates
	if (find(tmp.begin(), tmp.end(), i1) == tmp.end()) {
	tmp.push_back(i1);
	_theParticles += *i1;
	}
	if (find(tmp.begin(), tmp.end(), i2) == tmp.end()) {
	tmp.push_back(i2);
	_theParticles += *i2;
	}
	// Store accepted particle pairs
	_particlePairs += make_pair(i1, i2);
	if (_masstarget>0.0) {
	double diff=fabs(v4.mass()-_masstarget);
	if (diff<closestPair.first) {
	closestPair.first=diff;
	closestPair.second=make_pair(i1, i2);
	}
	}
	}
	}
	}
	if (_masstarget > 0.0 && closestPair.first < 1e30) {
	_theParticles.clear();
	_particlePairs.clear();
	_theParticles += closestPair.second.first;
	_theParticles += closestPair.second.second;
	_particlePairs += closestPair.second;
	}

	MSG_DEBUG("Selected " << _theParticles.size() << " particles " << "(" << _particlePairs.size() << " pairs)");
	if (getLog().isActive(Log::TRACE)) {
	for (const Particle& p : _theParticles) {
	- MSG_TRACE("ID: " << p.pid() << ", barcode: " << p.genParticle()->barcode());
	+ MSG_TRACE("ID: " << p.pid());
	}
	}
	}


	/// Constituent pairs
	const std::vector<std::pair<Particle, Particle> >& InvMassFinalState::particlePairs() const {
	return _particlePairs;
	}


	}

File Metadata

Mime Type: text/x-diff
Expires: Sat, Dec 21, 1:26 PM (19 h, 27 m)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 4022954
Default Alt Text: (69 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions