|
[Rivet-svn] r3261 - in trunk: . binblackhole at projects.hepforge.org blackhole at projects.hepforge.orgFri Jul 29 11:23:31 BST 2011
Author: buckley Date: Fri Jul 29 11:23:30 2011 New Revision: 3261 Log: Moving James' new aida2root2D to be the new aida2root. Deleted: trunk/bin/aida2root2D Modified: trunk/ChangeLog trunk/bin/aida2root Modified: trunk/ChangeLog ============================================================================== --- trunk/ChangeLog Thu Jul 28 18:43:40 2011 (r3260) +++ trunk/ChangeLog Fri Jul 29 11:23:30 2011 (r3261) @@ -1,3 +1,7 @@ +2011-07-29 Andy Buckley <andy at insectnation.org> + + * New version of aida2root from James Monk. + 2011-07-27 David Mallows <dave.mallows at gmail.com> * Updated MC_TTBAR.plot to reflect updated analysis. Modified: trunk/bin/aida2root ============================================================================== --- trunk/bin/aida2root Thu Jul 28 18:43:40 2011 (r3260) +++ trunk/bin/aida2root Fri Jul 29 11:23:30 2011 (r3261) @@ -5,11 +5,11 @@ Verify in the ROOT user manual what needs to be setup for use of ROOT with python -E.g. do the following additional setup steps: +For example, additional setup steps such as this may be required: setup setenv PYTHONDIR /usr setenv PATH $ROOTSYS/bin:$PYTHONDIR/bin:$PATH - setenv LD_LIBRARY_PATH $ROOTSYS/lib:$PYTHONDIR/lib/python2.3:$LD_LIBRARY_PATH - setenv PYTHONPATH $ROOTSYS/lib:$PYTHONDIR/lib/python2.3 + setenv LD_LIBRARY_PATH $ROOTSYS/lib:$PYTHONDIR/lib/python2.6:$LD_LIBRARY_PATH + setenv PYTHONPATH $ROOTSYS/lib:$PYTHONDIR/lib/python2.6 """ import sys @@ -18,11 +18,10 @@ sys.exit(1) +import os, array -import os -from array import array try: - from ROOT import TGraphAsymmErrors, TFile + from ROOT import ROOT, TGraphAsymmErrors, TGraph2DErrors, TH1F, TH2F, TFile except: sys.stderr.write("Couldn't find PyROOT module. If ROOT is installed, PyROOT probably is too, but it's not currently accessible\n") foundpyroot = False @@ -40,160 +39,6 @@ sys.stderr.write("You can test if PyROOT is properly set up by calling 'import ROOT' at the interactive Python shell\n") sys.exit(1) - -## TODO: replace with lighthisto -class Histo: - def __init__(self): - self._bins = [] - self.path = None - self.name = None - self.title = None - - def __cmp__(self, other): - """Sort by $path/$name string""" - return self.fullPath() > other.fullPath() - - def __str__(self): - out = "Histogram '%s' with %d bins\n" % (self.fullPath(), self.numBins()) - out += "Title: %s\n" % self.title - out += "\n".join([str(b) for b in self.getBins()]) - return out - - def fullPath(self): - return os.path.join(self.path, self.name) - - def asFlat(self): - out = "### %s\n" % self.fullPath() - out += "## Title: %s\n" % self.title - out += "## Area: %s\n" % self.area() - out += "## Num bins: %d\n" % self.numBins() - out += "## xlow xhigh yval yerrminus yerrplus\n" - out += "\n".join([b.asFlat() for b in self.getBins()]) - return out - - def asRoot(self): - xerrminus = array('f',self.numBins()*[0.]) - xerrplus = array('f',self.numBins()*[0.]) - xval = array('f',self.numBins()*[0.]) - yval = array('f',self.numBins()*[0.]) - yerrminus = array('f',self.numBins()*[0.]) - yerrplus= array('f',self.numBins()*[0.]) - ibin = 0 - for b in self.getBins(): - xval[ibin] = b.getXval() - yval[ibin] = b.getYval() - xerrminus[ibin], xerrplus[ibin] = b.getXErrors() - yerrminus[ibin], yerrplus[ibin] = b.getYErrors() - #print 'xerrminus=',xerrminus[ibin],' xerrplus=',xerrplus[ibin],' xval=',xval[ibin],' yval=',yval[ibin],' yerrplus=',yerrplus[ibin],' yerrminus=',yerrminus[ibin] - ibin = ibin +1 - tg = TGraphAsymmErrors(self.numBins(), xval, yval, xerrminus, xerrplus, yerrminus, yerrplus); - tg.SetTitle(self.title); - tg.SetName(self.name.replace("-", "_")); - return tg - - def numBins(self): - return len(self._bins) - - def getBins(self): - return sorted(self._bins) - - def setBins(self, bins): - self._bins = bins - return self - - def addBin(self, bin): - self._bins.append(bin) - return self - - def getBin(self, index): - self._bins.sort() - return self.getBins()[index] - - bins = property(getBins, setBins) - - def area(self): - return sum([bin.area() for bin in self.bins]) - - def __iter__(self): - return iter(self.getBins()) - - def __len__(self): - return len(self._bins) - - def __getitem__(self, index): - return self.getBin(index) - - -## TODO: replace with lighthisto -class Bin: - """A simple container for a binned value with an error.""" - def __init__(self, xval=0, xerrminus=None, xerrplus=None, yval=0, yerrminus=0, yerrplus=0, focus=None): - self.xval = xval - self.xerrplus = xerrplus - self.xerrminus= xerrminus - self.yval = yval - self.yerrplus = yerrplus - self.yerrminus = yerrminus - self.focus= focus - - def __str__(self): - out = "%e to %e: %e +- %e" % (self.xval-self.xerrminus, self.xval+self.xerrplus, self._yval, self._yerr) - return out - - def asFlat(self): - out = "%e %e %e %e %e" % (self.xval-self.xerrminus, self.xval+self.xerrplus, self.yval, self.yerrminus, self.yerrplus) - return out - - def __cmp__(self, other): - """Sort by mean x value (yeah, I know...)""" - return (self.xval) > (other.xval) - - def getXRange(self): - return (self.xval-self.xerrminus, self.xval+self.xerrplus) - - def getXErrors(self): - return (self.xerrminus, self.xerrplus) - - def setXRange(self, xlow, xhigh): - self.xlow = xlow - self.xhigh = xhigh - return self - - def getYErrors(self): - return (self.yerrminus, self.yerrplus) - - def getBinCenter(self): - """Geometric middle of the bin range.""" - return self.xlow + .5*(self.xhigh - self.xlow) - - def getFocus(self): - """Mean x-value of the bin.""" - if self.focus is not None: - return (self.xlow + self.xhigh)/2.0 - else: - return self.focus - - def getXval(self): - return self.xval - - def getYval(self): - return self.yval - - def area(self): - return self.yval * (self.xerrplus - self.xerrminus) - - def getYErr(self): - """Get mean of +ve and -ve y-errors.""" - return (self.yerrplus + self.yerrminus)/2.0 - - def setYErr(self, yerr): - """Set both +ve and -ve y-errors simultaneously.""" - self.yerrplus = yerr - self.yerrminus = yerr - return self - - - ## Try to load faster but non-standard cElementTree module try: import xml.etree.cElementTree as ET @@ -208,33 +53,285 @@ sys.exit(1) +class Histo: + + def __init__(self, nDim): + self._points = [] + self.name = "" + self.title = "" + self._nDim = nDim + + def addPoint(self, dp): + if dp.dimensionality() != self._nDim: + er = "Tried to add a datapoint of dimensionality " + str(dp.dimensionality()) + " to a histogram of dimensionality " + str(self._nDim) + sys.stderr.write(er) + sys.exit(1) + self._points.append(dp) + + def numPts(self): + return len(self._points) + + def asTGraph(self): + tg = TGraph() + tg.SetName(self.name) + tg.SetTitle(self.title) + return tg + + def asHisto(self): + tg = self.asTGraph() + histo = tg.Histogram().Clone() + return histo + + @staticmethod + def equalFloats(left, right, precision=1.e-6): + + try: + test = abs((left - right) / (left + right)) + return test < precision + except ZeroDivisionError: + if left * right < 0.: + return False + else: + return True + + return False + + +class Histo2D(Histo): + + def __init__(self): + Histo.__init__(self,3) + + def asTGraph(self): + xs = array.array("d", []) + ex = array.array("d", []) + ys = array.array("d", []) + ey = array.array("d", []) + zs = array.array("d", []) + ez = array.array("d", []) + + for pt in self._points: + x = pt.mean(0) + erx = pt.er(0) + y = pt.mean(1) + ery = pt.er(1) + z = pt.mean(2) + erz = pt.er(2) + + xs.append(x) + ex.append(erx) + ys.append(y) + ey.append(ery) + zs.append(z) + ez.append(erz) + + if self.numPts() == 0: + tg = TGraph2DErrors() + er = "Tried to create TGraph2DErrors called " + self.name + " with zero datapoints" + else: + tg = TGraph2DErrors(self.numPts(), xs, ys, zs, ex, ey, ez) + tg.SetTitle(self.title) + tg.SetName(self.name.replace("-", "_")) + return tg + + def asTHisto(self): + + if self.numPts() == 0: + histo = TH2F() + histo.SetName(self.name) + return histo + + tmpXEdges = [] + tmpYEdges = [] + + for pt in self._points: + tmpXEdges.append(pt.lowEdge(0)) + tmpXEdges.append(pt.highEdge(0)) + tmpYEdges.append(pt.lowEdge(1)) + tmpYEdges.append(pt.highEdge(1)) + + sortedX = sorted(tmpXEdges) + sortedY = sorted(tmpYEdges) + + xBinEdges = array.array("d", [sortedX[0]]) + yBinEdges = array.array("d", [sortedY[0]]) + + for edge in sortedX: + if not Histo.equalFloats(edge, xBinEdges[-1]): + xBinEdges.append(edge) + + for edge in sortedY: + if not Histo.equalFloats(edge, yBinEdges[-1]): + yBinEdges.append(edge) + + histo = TH2F(self.name, self.title, len(xBinEdges)-1, xBinEdges, len(yBinEdges)-1, yBinEdges) + histo.Sumw2() + + for pt in self._points: + bin = histo.FindBin(pt.value(0), pt.value(1)) + histo.SetBinContent(bin, pt.value(2)) + histo.SetBinError(bin, pt.er(2)) + + return histo + + +class Histo1D(Histo): + def __init__(self): + Histo.__init__(self,2) + + def asTGraph(self): + xerrminus = array.array("d", []) + xerrplus = array.array("d", []) + xval = array.array("d", []) + yval = array.array("d", []) + yerrminus = array.array("d", []) + yerrplus = array.array("d", []) + + for pt in self._points: + x = pt.value(0) + xplus = pt.erUp(0) + xminus = pt.erDn(0) + + y = pt.value(1) + yplus = pt.erUp(1) + yminus = pt.erDn(1) + + xval.append(x) + xerrminus.append(xminus) + xerrplus.append(xplus) + yval.append(y) + yerrminus.append(yminus) + yerrplus.append(yplus) + + tg = TGraphAsymmErrors(self.numPts(), xval, yval, xerrminus, xerrplus, yerrminus, yerrplus) + tg.SetTitle(self.title) + tg.SetName(self.name.replace("-", "_")) + return tg + + def asTHisto(self): + + if self.numPts() == 0: + histo = TH1F() + histo.SetName(self.name) + return histo + + binEdges = array.array("d", []) + binEdges.append(self._points[0].lowEdge(0)) + + bin = 0 + binNumbers = [] + + for pt in self._points: + lowEdge = pt.lowEdge(0) + highEdge = pt.highEdge(0) + if not Histo1D.equalFloats(lowEdge, binEdges[-1]): + binEdges.append(lowEdge) + bin = bin + 1 + + bin = bin + 1 + binEdges.append(highEdge) + binNumbers.append(bin) + + histo = TH1F(self.name, self.title, self.numPts(), binEdges) + histo.Sumw2() + + for i, pt in enumerate(self._points): + histo.SetBinContent(binNumbers[i], pt.value(1)) + histo.SetBinError(binNumbers[i], pt.er(1)) + + return histo + + +class DataPoint: + + def __init__(self): + self._dims = 0 + self._coords = [] + self._erUps = [] + self._erDns = [] + + def setCoord(self, val, up, down): + self._dims = self._dims + 1 + self._coords.append(val) + self._erUps.append(up) + self._erDns.append(down) + + def dimensionality(self): + return self._dims + + def th(self, dim): + th = "th" + if dim == 1: + th = "st" + elif dim == 2: + th = "nd" + elif dim == 3: + th = "rd" + return th + + def checkDimensionality(self, dim): + if dim >= self.dimensionality(): + er = "Tried to obtain the " + str(dim) + self.th(dim) + " dimension of a " + str(self.dimensionality()) + " dimension DataPoint" + sys.stderr.write(er) + sys.exit(1) + + def value(self, dim): + self.checkDimensionality(dim) + return self._coords[dim] + + def erUp(self, dim): + self.checkDimensionality(dim) + return self._erUps[dim] + + def erDn(self, dim): + self.checkDimensionality(dim) + return self._erDns[dim] + + def mean(self, dim): + val = self.value(dim) + 0.5 * (self.erUp(dim) - self.erDn(dim)) + return val + + def er(self, dim): + ee = 0.5 * (self.erUp(dim) + self.erDn(dim)) + return ee + + def lowEdge(self, dim): + return self.value(dim) - self.erDn(dim) + + def highEdge(self, dim): + return self.value(dim) + self.erUp(dim) + def mkHistoFromDPS(dps): - """Make a mini histo representation from an AIDA dataPointSet tag.""" - myhist = Histo() - myhist.name = dps.get("name") - myhist.title = dps.get("title") - myhist.path = dps.get("path") - points = dps.findall("dataPoint") - numbins = len(points) - for binnum, point in enumerate(points): - bin = Bin() - for d, m in enumerate(point.findall("measurement")): - val = float(m.get("value")) - down = float(m.get("errorMinus")) - up = float(m.get("errorPlus")) - if d == 0: - low = val - down - high = val + up - bin.setXRange(low, high) - bin.xval = val - bin.xerrplus = up - bin.xerrminus = down - elif d == 1: - bin.yval = val - bin.yerrplus = up - bin.yerrminus = down - myhist.addBin(bin) - return myhist + + dim = dps.get("dimension") + + is3D = False + if dim == "3": + myhist = Histo2D() + is3D = True + else: + myhist = Histo1D() + + myhist.name = dps.get("name") + myhist.title = dps.get("title") + myhist.path = dps.get("path") + + points = dps.findall("dataPoint") + numbins = len(points) + + for ptNum, point in enumerate(points): + dp = DataPoint() + for d, m in enumerate(point.findall("measurement")): + val = float(m.get("value")) + down = float(m.get("errorMinus")) + up = float(m.get("errorPlus")) + dp.setCoord(val, up, down) + myhist.addPoint(dp) + + return myhist + + +################################## from optparse import OptionParser @@ -245,6 +342,13 @@ parser.add_option("-m", "--match", action="append", help="Only write out histograms whose $path/$name string matches these regexes", dest="PATHPATTERNS") +parser.add_option("-g", "--tgraph", action="store_true", default=True, + help="Store output as ROOT TGraphAsymmErrors or TGraph2DErrors", + dest="TGRAPH") +parser.add_option("-t", "--thisto", action="store_false", + help="Store output as ROOT TH1 or TH2", + dest="TGRAPH") + opts, args = parser.parse_args() if opts.PATHPATTERNS is None: opts.PATHPATTERNS = [] @@ -253,11 +357,11 @@ sys.stderr.write("Must specify at least one AIDA histogram file\n") sys.exit(1) -import re for aidafile in args: - out = sys.stdout + tree = ET.parse(aidafile) histos = [] + for dps in tree.findall("dataPointSet"): useThisDps = True if len(opts.PATHPATTERNS) > 0: @@ -268,13 +372,17 @@ useThisDps = True break if useThisDps: - histos.append(mkHistoFromDPS(dps)) + histos.append(mkHistoFromDPS(dps)) if len(histos) > 0: if opts.SMARTOUTPUT: outfile = os.path.basename(aidafile).replace(".aida", ".root") out = TFile(outfile,"RECREATE") for h in sorted(histos): - h.asRoot().Write() + if opts.TGRAPH: + h.asTGraph().Write() + else: + h.asTHisto().Write() + out.Close() else: sys.stderr.write("ROOT objects must be written to a file")
More information about the Rivet-svn mailing list |