[Rivet] ATLAS analysis to be included into Rivet

[Roman Lysak]

  Dear Rivet authors,

we have prepared an analysis from Atlas experiment which we would like 
to include into Rivet.
We validated that the routine produces consistent results comparing to 
results in the paper.
All necessary information can be found in the attachment.

Please let me know if anything else is required.

Thanks,
   Roman Lysak

BTW, while using rivet tools, I needed to make simple modification to 
two scripts in order to be more useful:
1) aidamerge - I added possibility to sum the .aida files (the original 
version makes an average)
2) root2flat - I modified it to work with TGraphErrors

I'm attaching both in case you would think they could be useful for the 
others.

-------------- next part --------------
A non-text attachment was scrubbed...
Name: ATLAS_2012_I1188891.aida
Type: text/xml
Size: 9778 bytes
Desc: not available
URL: <http://www.hepforge.org/lists-archive/rivet/attachments/20121128/a0c341d6/attachment.xml>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: ATLAS_2012_I1188891.cc
Type: text/x-c++src
Size: 6785 bytes
Desc: not available
URL: <http://www.hepforge.org/lists-archive/rivet/attachments/20121128/a0c341d6/attachment.cc>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: ATLAS_2012_I1188891.info
Type: application/x-info
Size: 1467 bytes
Desc: not available
URL: <http://www.hepforge.org/lists-archive/rivet/attachments/20121128/a0c341d6/attachment.info>
-------------- next part --------------
# BEGIN PLOT /ATLAS_2012_I1188891/d01-x01-y01
Title=
XLabel=Jet $p_T$~[GeV]
YLabel=BB fraction [$\%$]
# END PLOT
# BEGIN PLOT /ATLAS_2012_I1188891/d02-x01-y01
Title=
XLabel=Jet $p_T$~[GeV]
YLabel=BC fraction [$\%$]
# END PLOT
# BEGIN PLOT /ATLAS_2012_I1188891/d03-x01-y01
Title=
XLabel=Jet $p_T$~[GeV]
YLabel=CC fraction [$\%$]
# END PLOT
# BEGIN PLOT /ATLAS_2012_I1188891/d04-x01-y01
Title=
XLabel=Jet $p_T$~[GeV]
YLabel=BU fraction [$\%$]
# END PLOT
# BEGIN PLOT /ATLAS_2012_I1188891/d05-x01-y01
Title=
XLabel=Jet $p_T$~[GeV]
YLabel=CU fraction [$\%$]
# END PLOT
# BEGIN PLOT /ATLAS_2012_I1188891/d06-x01-y01
Title=
XLabel=Jet $p_T$~[GeV]
YLabel=UU fraction [$\%$]
# END PLOT
-------------- next part --------------
#! /usr/bin/env python

import sys, os, copy
from math import sqrt
from subprocess import Popen
import lighthisto

## 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)

from optparse import OptionParser
parser = OptionParser(usage="%prog aidafile [aidafile2 ...]")
parser.add_option("-o", "--outfile", dest="OUTFILE",
                  default="merged.aida", help="file for merged aida output.")
parser.add_option("-s", "--sum",
                  action="store_true", dest="performSum", default=False,
                  help="sum the bin values instead of averaging")
opts, args = parser.parse_args()
headerprefix = ""

if len(args) < 1:
    sys.stderr.write("Must specify at least one AIDA histogram file\n")
    sys.exit(1)


try:
    outaida = open(opts.OUTFILE, "w")
except:
    sys.stderr.write("Couldn't open outfile %s for writing." % opts.OUTFILE)

try:
    outdat = open(opts.OUTFILE.replace(".aida", ".dat"), "w")
except:
    sys.stderr.write("Couldn't open outfile %s for writing." % opts.OUTFILE.replace(".aida", ".dat"))

## Get histos
inhistos = {}
weights = {}
for aidafile in args:
    tree = ET.parse(aidafile)
    for dps in tree.findall("dataPointSet"):
        h = lighthisto.Histo.fromDPS(dps)
        if not inhistos.has_key(h.fullPath()):
            inhistos[h.fullPath()] = {}
        inhistos[h.fullPath()][aidafile] = h

## Merge histos
outhistos = {}
for path, hs in inhistos.iteritems():
    #print path, hs
    outhistos[path] = copy.deepcopy(hs.values()[0])
    for i, b in enumerate(outhistos[path].getBins()):
        sum_val = 0.0
        sum_err2 = 0.0
        n = 0
        for infile, h in hs.iteritems():
            sum_val += h.getBin(i).val
            try:
                sum_err2 += h.getBin(i).getErr()**2
            except OverflowError:
                # in case the **2 produces overflow errors
                # set sum to 'inf'
                sum_err2 = float('inf')
            n += 1
        if (opts.performSum): outhistos[path].getBin(i).val = sum_val
        else:                 outhistos[path].getBin(i).val = sum_val / n

        try:
            if (opts.performSum): outhistos[path].getBin(i).setErr(sum_err2**0.5 )
            else:                 outhistos[path].getBin(i).setErr(sum_err2**0.5 / n)
        except OverflowError:
            # to get back to numerics, replace an eventual 'inf' 
            # in sum_err2 with max float available ~1.e+308
            outhistos[path].getBin(i).setErr(1.e308)

## Write out merged histos
#print sorted(outhistos.values())
outdat.write("\n\n".join([h.asFlat() for h in sorted(outhistos.values())]))
outdat.write("\n")
outdat.close()

Popen(["flat2aida", opts.OUTFILE.replace(".aida", ".dat")], stdout=outaida).wait()

os.unlink(opts.OUTFILE.replace(".aida", ".dat"))
-------------- next part --------------
#! /usr/bin/env python

"""%prog

Read in .root files and write out the histograms in FLAT
format suitable for make-plots.

Use e.g. 'root2flat file1.root -o flathistos' to produce the flat histogram
files which need to be converted to AIDA files using Rivet's flat2aida tool

Usage:
    %prog [options] file1.root

Use --help to get information for all options.
"""

import sys
if sys.version_info[:3] < (2,4,0):
    print "rivet scripts require Python version >= 2.4.0... exiting"
    sys.exit(1)


import os, optparse, logging, ROOT
# try:
#     from IPython.Shell import IPShellEmbed
#     ipshell = IPShellEmbed([])
# except:
#     print "Ipython shell not available."


## Parse options
parser = optparse.OptionParser(usage = __doc__)
parser.add_option("-o", dest="OUTDIR", default=".",
                  help = "specify directory in which to write out converted files")
(opts, args) = parser.parse_args()


def readROOT(rootfile):
    """ This is the main function that opens a ROOT file, browses its contents
        for histograms and tries to write the converted histos to files.
    """
    # Open the ROOT file
    f = ROOT.TFile(rootfile)

    # Initial browse to see the structure
    subdirs, histonames, tgraphnames = browse(f)

    # Keep browding if there are subdirectories TODO: Make this work for
    # arbitrarily deep directory structures
    if len(subdirs) > 0:
        for sd in subdirs:
            t_s, t_h = browse(f, sd)
            histonames.extend(t_h)

    #primary_keys = f.GetListOfKeys()
    #iter_primaries = ROOT.TListIter(primary_keys)

    # This will convert and write the histos
    for histoname in histonames:
        writeHisto(histoname, f.Get(histoname))
    for tgraphname in tgraphnames:
        writeHisto(tgraphname, f.Get(tgraphname), True)



def browse(f, branch=None):
    """ This function browses a file/branch, trying to find objects that
        either inherit from TH1 or from TProfile.
    """
    # Prepare return values
    histos  = []
    tgraphs = []
    subdirs = []

    # Get Iterator
    if branch:
        f.Cd(branch)

    primary_keys = ROOT.gDirectory.GetListOfKeys()
    iter_primaries = ROOT.TListIter(primary_keys)

    if branch:
        f.Cd('..')

    # Iterate over iterator
    for i in xrange(len(primary_keys)):
        if branch:
            t_n = branch + '/' + iter_primaries.Next().GetName()
        else:
            t_n = iter_primaries.Next().GetName()
        # Make sure we don't have a NoneType object here
        if f.Get(t_n):
            # Check if the curent hing is a directory
            if type(f.Get(t_n)) == ROOT.TDirectoryFile:
                subdirs.append(t_n)
            # Check if the curent hing is a histogram
            elif f.Get(t_n).InheritsFrom("TH1") or f.Get(t_n).InheritsFrom("TProfile"):
                histos.append(t_n)
            elif f.Get(t_n).InheritsFrom("TGraphAsymmErrors") or f.Get(t_n).InheritsFrom("TGraphErrors"):
                tgraphs.append(t_n)

    return subdirs, histos, tgraphs


def convertHisto(R_histo):
    """ This function reads a single ROOT histo and converts it into the
        FLAT format suitable for plotting with make-plots.
    """
    title = R_histo.GetTitle().replace("#","\\")
    xtitle= R_histo.GetXaxis().GetTitle().replace("#","\\")
    ytitle= R_histo.GetYaxis().GetTitle().replace("#","\\")
    bins = getBinsFromTH1F(R_histo)
    for bin in bins:
        try:
            binstr = ""
            for key in ["xlow", "xhigh", "y", "y_err_low", "y_err_high"]:
                binstr += "%s: %e " % (key, bin[key])
            logging.info(binstr)
        except:
            pass
    return title, xtitle, ytitle, bins


def convertTGraph(TGraph):
    title = TGraph.GetTitle().replace("#","\\")
    xtitle= TGraph.GetXaxis().GetTitle().replace("#","\\")
    ytitle= TGraph.GetYaxis().GetTitle().replace("#","\\")
    bins = getBinsFromTGraph(TGraph)
    return title, xtitle, ytitle, bins


def getBinsFromTH1F(R_histo):
    """ A little helper function that returns a list of bin-dictionaries).
    """
    allbins=[]
    for ii in xrange(R_histo.GetNbinsX()):
        i = ii + 1
        xlow  = R_histo.GetBinLowEdge(i)
        xhigh = xlow + R_histo.GetBinWidth(i)
        y     = R_histo.GetBinContent(i)
        y_err = R_histo.GetBinError(i)
        bin = {"xlow":xlow, "xhigh":xhigh, "y":y, "y_err_low":y_err, "y_err_high":y_err}
        allbins.append(bin)
    return allbins


def getBinsFromTGraph(TGraph):
    allbins=[]
    X = TGraph.GetX()
    Y = TGraph.GetY()
    for ii in xrange(TGraph.GetN()):
        i = ii + 1
        xlow  = X[i] - TGraph.GetErrorXlow(i)
        xhigh = X[i] + TGraph.GetErrorXhigh(i)
        y     = Y[i]
        y_err_low = TGraph.GetErrorYlow(i)
        y_err_high = TGraph.GetErrorYhigh(i)
        bin = {"xlow":xlow, "xhigh":xhigh, "y":y, "y_err_low":y_err_low, "y_err_high":y_err_high}
        allbins.append(bin)
    return allbins


def writeHisto(name, R_histo, tgraph=False):
    """ This writes the histogram into a single file, ready to plot with
    make-plots.
    """
    if tgraph:
        title, xlabel, ylabel, bins = convertTGraph(R_histo)
    else:
        title, xlabel, ylabel, bins = convertHisto(R_histo)
    head = "# BEGIN PLOT\nTitle=%s\nLegend=1\nLogY=1\nDrawOnly=%s\n" % (title, name)
    head += "XLabel=%s\nYLabel=%s\n# END PLOT\n" % (xlabel, ylabel)

    histo = getFlatHisto(bins, name, title)

    flatname = name.replace("/","_") + ".dat"
    if flatname.startswith("_"):
        flatname = flatname[1:]
    flatfile = os.path.join(opts.OUTDIR, flatname)
    f = open(flatfile, "w")
    f.write(head)
    f.write("\n")
    f.write(histo)
    f.close()


def getFlatHisto(bins, name, title):
    """ This returns a histo in the FLAT format. """
    histo= "# BEGIN HISTOGRAM %s\n" % name
    histo += "LineColor=black\n"
    histo += "ErrorBars=1\n"
    histo += "PolyMarker=*\n"
    histo += "Title=%s\n" % title
    for bin in bins:
        histo += "%.8e\t%.8e\t%.8e\t%.8e\t%.8e\n" % (bin["xlow"], bin["xhigh"],
                                                     bin["y"], bin["y_err_low"], bin["y_err_high"])
    histo += "# END HISTOGRAM\n"
    return histo


if __name__ == "__main__":
    for infile in args:
        if not os.path.exists(opts.OUTDIR):
            os.mkdir(opts.OUTDIR)
        readROOT(infile)
    print "Done. Written all plot files to %s" % opts.OUTDIR