|
|
[Rivet-svn] r2196 - trunk/bin
blackhole at projects.hepforge.org
blackhole at projects.hepforge.org
Thu Dec 17 14:36:00 GMT 2009
Author: holsch
Date: Thu Dec 17 14:36:00 2009
New Revision: 2196
Log:
Changing name of normalisation script
Added:
trunk/bin/rivet-normalise-histos (contents, props changed)
Deleted:
trunk/bin/normalise-to-refdata
Modified:
trunk/bin/Makefile.am
Modified: trunk/bin/Makefile.am
==============================================================================
--- trunk/bin/Makefile.am Thu Dec 17 14:34:10 2009 (r2195)
+++ trunk/bin/Makefile.am Thu Dec 17 14:36:00 2009 (r2196)
@@ -1,5 +1,5 @@
bin_SCRIPTS = rivet-config rivet-buildplugin
-dist_bin_SCRIPTS = aida2flat aida2root compare-histos make-html make-plots rivet-mkanalysis chop_bins normalise-to-refdata
+dist_bin_SCRIPTS = aida2flat aida2root compare-histos make-html make-plots rivet-mkanalysis chop_bins rivet-normalise-histos
EXTRA_DIST = flat2aida
if ENABLE_PYEXT
Added: trunk/bin/rivet-normalise-histos
==============================================================================
--- /dev/null 00:00:00 1970 (empty, because file is newly added)
+++ trunk/bin/rivet-normalise-histos Thu Dec 17 14:36:00 2009 (r2196)
@@ -0,0 +1,275 @@
+#!/usr/bin/env python
+
+"""%prog -r <REFDATAPATH> -O <observable-file> <AIDAFILE>
+
+Normalise histos in observable-file of AIDAFILE to the area of the
+corresponding histos in REFAIDAPATH. REFAIDAPATH can either be
+a single AIDA-file or a directory containing AIDA-files.
+
+Or, if no REFAIDAPATH is given, normalise to new area given in
+observable-file "observables_and_areas", e.g.:
+CDF_2000_S4155203_d01-x01-y01 1.0
+
+You can also define bin to chop out in the observable-file in the
+same way as in chop_bins:
+CDF_2000_S4155203_d01-x01-y01:0:35 1.0
+This will chop the bins with Z-pT > 35 GeV and afterwards perform
+the renormalisation.
+
+
+Example:
+ %prog -r /usr/share/Rivet/ -O observables out.aida
+This will return the Z-boson pT-distribution, normalised to the histo-area
+measured by CDF.
+ %prog -O observables_and_areas out.aida
+This will return the Z-boson pT-distribution, normalised to 1.0
+
+"""
+
+import os, re, sys, logging
+from lighthisto import Histo
+## Make "sorted" a builtin function on Python < 2.4
+if not 'sorted' in dir(__builtins__):
+ def sorted(iterable, cmp=None, key=None, reverse=None):
+ rtn = iterable
+ rtn.sort(cmp)
+ return rtn
+
+## Add logging.log if needed
+if not 'log' in dir(logging):
+ def _logit(level, msg):
+ l = logging.getLogger()
+ l.log(level, msg)
+ logging.log = _logit
+
+try:
+ from IPython.Shell import IPShellEmbed
+ ipshell = IPShellEmbed([])
+except:
+ logging.info("Ipython shell not available.")
+
+
+## Try to load faster but non-standard cElementTree module
+try:
+ import xml.etree.cElementTree as ET
+except ImportError:
+ try:
+ import cElementTree as ET
+ except ImportError:
+ try:
+ import xml.etree.ElementTree as ET
+ except:
+ sys.stderr.write("Can't load the ElementTree XML parser: please install it!\n")
+ sys.exit(1)
+
+
+def getHistosFromAIDA(aidafile):
+ '''Get a dictionary of histograms indexed by name.'''
+ if not re.match(r'.*\.aida$', aidafile):
+ logging.error("Error: input file '%s' is not an AIDA file" % aidafile)
+ sys.exit(2)
+ aidafilepath = os.path.abspath(aidafile)
+ if not os.access(aidafilepath, os.R_OK):
+ logging.error("Error: cannot read from %s" % aidafile)
+ sys.exit(2)
+
+ histos = {}
+ tree = ET.parse(aidafilepath)
+ for dps in tree.findall("dataPointSet"):
+ ## Get this histogram's path name
+ dpsname = os.path.join(dps.get("path"), dps.get("name"))
+ ## Is it a data histo?
+ h = Histo.fromDPS(dps)
+ h.isdata = dpsname.upper().startswith("/REF")
+ if h.isdata:
+ dpsname = dpsname.replace("/REF", "")
+ histos[dpsname] = h
+ return histos
+
+def getRefHistos(refpath):
+ """ Return dictionary of reference histos {name: histo}.
+ Refpath can either be a single file or a directory.
+ """
+ refhistos = {}
+ if not refpath is None:
+ try:
+ refhistos=getHistosFromAIDA(refpath)
+ except:
+ for aida in os.listdir(refpath):
+ if aida.endswith(".aida"):
+ temp = getHistosFromAIDA(os.path.join(refpath, aida))
+ for k, v in temp.iteritems():
+ if not k in refhistos.keys():
+ refhistos[k]=v
+ return refhistos
+
+def readObservableFile(obsfile):
+ """ Read observables to normalise from file obsfile.
+ Return-values are a list of the histo-names to normalise and a
+ dictionary with name:newarea entries.
+ """
+ obslist = []
+ obsnorms = {}
+ bindefs = {}
+
+ if obsfile is not None:
+ try:
+ f = open(obsfile, 'r')
+ except:
+ logging.error("Cannot open histo list file %s" % opts.obsfile)
+ sys.exit(2)
+ for line in f:
+ stripped = line.strip()
+ # Skip empty or commented lines
+ if len(stripped) == 0 or stripped.startswith("#"):
+ continue
+
+ # Split the line to find out whether newarea is given in obsfile
+ splitline = stripped.split()
+
+ # Try to read bin specs for chopping
+ try:
+ path, low, high = splitline[0].split(":")
+ except:
+ path = splitline[0].split(":")[0]
+ low = ""
+ high = ""
+ logging.warning("No bin definition given for %s" % (line))
+ #sys.exit(1)
+ if low == "":
+ low = None
+ else:
+ low = float(low)
+ if high == "":
+ high = None
+ else:
+ high = float(high)
+ bindefs[path] = (low, high)
+
+ # Try to read areas to normalise to from obsfile
+ if len(splitline) == 2:
+ obsnorms[path] = float(splitline[1])
+ obslist.append(path)
+ f.close()
+ return obslist, obsnorms, bindefs
+
+
+if __name__ == "__main__":
+ from optparse import OptionParser, OptionGroup
+ parser = OptionParser(usage=__doc__)
+
+ parser.add_option("-O", "--obsfile", default=None,
+ help="Specify a file with histograms (and bin ranges) " +
+ " that are to be normalised.")
+ parser.add_option("-r", "--refdir", default=None,
+ help="File of folder with reference histos")
+ parser.add_option("-o", "--outdir",
+ dest="outdir", default="renormalised",
+ help="output directory (default: %default)")
+ parser.add_option("-a", dest="AIDA", default=True, action="store_true",
+ help="Produce AIDA output rather than FLAT")
+ parser.add_option("-f", dest="AIDA", action="store_false",
+ help="Produce FLAT output rather than AIDA")
+ verbgroup = OptionGroup(parser, "Verbosity control")
+ verbgroup.add_option("-V", "--verbose", action="store_const",
+ const=logging.DEBUG, dest="LOGLEVEL",
+ help="print debug (very verbose) messages")
+ verbgroup.add_option("-Q", "--quiet", action="store_const",
+ const=logging.WARNING, dest="LOGLEVEL",
+ help="be very quiet")
+ parser.set_defaults(bins=[],
+ outdir=".",
+ LOGLEVEL=logging.INFO)
+ opts, args = parser.parse_args()
+
+
+ # Configure logging
+ try:
+ logging.basicConfig(level=opts.LOGLEVEL, format="%(message)s")
+ except:
+ pass
+ h = logging.StreamHandler()
+ h.setFormatter(logging.Formatter("%(message)s"))
+ logging.getLogger().setLevel(opts.LOGLEVEL)
+ if logging.getLogger().handlers:
+ logging.getLogger().handlers[0] = h
+ else:
+ logging.getLogger().addHandler(h)
+
+
+
+ # Read in histogram files to normalise
+ histos = getHistosFromAIDA(args[0])
+
+ # Read in reference histos to get reference areas to normalise to
+ refhistos = getRefHistos(opts.refdir)
+
+ # Read in observables, if no observable file is given or observable file
+ # is empty all observables will be renormalised if possible
+ obslist, obsnorms, bindefs = readObservableFile(opts.obsfile)
+ if len(obslist) == 0:
+ logging.warning("No observable file given or empty, will try to ",
+ "normalise all histos.")
+ obslist = histos.keys()
+
+
+ # Create output filename
+ base = args[0].split("/")[-1].split(".aida")[0]
+ if opts.AIDA:
+ outfile = os.path.join(opts.outdir, base + "-renormed.aida")
+ else:
+ outfile = os.path.join(opts.outdir, base + "-renormed.dat")
+ if not os.path.exists(opts.outdir):
+ logging.info("Creating outdir %s"%opts.outdir)
+ os.mkdir(opts.outdir)
+
+ aidaheader = """<?xml version="1.0" encoding="ISO-8859-1" ?>
+<!DOCTYPE aida SYSTEM "http://aida.freehep.org/schemas/3.3/aida.dtd">
+<aida version="3.3">
+ <implementation version="1.1" package="FreeHEP"/>
+ """
+ # Open file for output
+ f = open(outfile, "w")
+ if opts.AIDA:
+ f.write(aidaheader)
+
+
+ # Iterate over all histos
+ for name, histo in histos.iteritems():
+ # Find out whether ref-histo is given
+ if name in refhistos.keys():
+ tempref = refhistos[name]
+ else:
+ tempref = histos[name]
+
+ # Try to chop bins
+ if name in bindefs.keys():
+ tempref = tempref.chop(bindefs[name])
+ tempold = histos[name].chop(bindefs[name])
+ else:
+ #tempref = refhistos[name]
+ tempold = histos[name]
+
+ # Get old and new histogram areas
+ oldarea = tempold.getArea()
+ if name in obsnorms.keys():
+ newarea = obsnorms[name]
+ else:
+ newarea = tempref.getArea()
+
+ # Renormalise histo
+ logging.info("Renormalising %s from %f to %f"%(name, oldarea,
+ newarea))
+ renormed = tempold.renormalise(newarea)
+
+ # Write to file
+ if opts.AIDA:
+ f.write(renormed.asAIDA())
+ else:
+ f.write(renormed.asFlat())
+
+
+ if opts.AIDA:
+ f.write("</aida>")
+
+ logging.info("Done. Written output to %s."%outfile)
More information about the Rivet-svn
mailing list
|