diff --git a/AnalysisTools/contur/bin/contur-BL3plot b/AnalysisTools/contur/bin/contur-BL3plot
--- a/AnalysisTools/contur/bin/contur-BL3plot
+++ b/AnalysisTools/contur/bin/contur-BL3plot
@@ -1,165 +1,209 @@
#! /usr/bin/env python
from optparse import OptionParser
import os, sys
import contur.Utils as util
import contur as ct
from contur.Plotting import *
import pickle
import os, sys
import numpy as np
import matplotlib.pyplot as plt
import scipy.stats as sp
from os.path import join
import errno
from matplotlib.ticker import *
+from collections import defaultdict
parser = OptionParser(usage=__doc__)
parser.add_option("-o", "--outputdir", dest="OUTPUTDIR",
default="plots", help="Specify output directory for output plots.")
-#parser.add_option("-a", "--analysisdir", dest="ANALYSISDIR",
-# default="ANALYSIS", help="Output directory for analysis cards.")
parser.add_option("-v", "--version",
action="store_true", dest="printVersion",
default=False, help="print version number and exit.")
-
opts, mapfiles = parser.parse_args()
if opts.printVersion:
util.writeBanner()
sys.exit(0)
if not mapfiles:
sys.stderr.write("Error: You need to specify some contur.map files to be analysed!\n")
sys.exit(1)
-# decode MH from the x axis variable
+# function to decode MH from the x axis variable
def mh2(x,dummy=0):
return int(x)
-# decode g' from the y axis variable
+# function to decode sin(alpha) from the y axis variable
def salpha(y,dummy=0):
y = np.sin(y*np.pi/40.0)
return float('%3.1f' % y )
- #return float('%4.2f' % np.power(10,-y/4.0))
- #return 1.0*np.power(10,-y/4.0)
-
-# inverse of mzprime
-def xfromMzp(mh2):
- return mh2
if __name__ == "__main__":
util.mkoutdir(opts.OUTPUTDIR)
contourXaxis=[]
contourYaxis=[]
+# Axis labels (should get these from the conturDepot eventually?)
+ xAxisLabel = r"$M_{H2}$ [GeV]"
+ yAxisLabel = r"$\sin \alpha$"
+
for m in mapfiles:
with open(m, 'r+b') as f:
x = pickle.load(f)
- print m
- print len(x)
- # this is maximum possible size (ie too big!)
+
+ n_pools = len(x[0].sortedPoints)
+ print "Loaded file", m, " which has ", n_pools, " pools"
+
+ # array to store the total CL values this is maximum possible size (ie too big!)
confLim = np.zeros((len(x),len(x)))
+ # Build a dictionary to store the CL values for individual pools
+ confLims = defaultdict(list)
+ #bestPlots = defaultdict(list)
+ for ctp in x[0].sortedPoints:
+ pool = ctp.pools
+ confLims[pool] = np.zeros((len(x),len(x)))
+ #bestPlots[pool] = np.array((len(x),len(x)),dtype='string')
+
# Sort them so the parameter values are in order.
x.sort(key=lambda i: (float(i.ModelParam1), float(i.ModelParam2)))
# first loop to find the axes ranges.
for cDepot in x:
- print "New conturDepot ----"
- print cDepot.ModelParam1, cDepot.ModelParam2
+ #print "New conturDepot ----"
+ #print cDepot.ModelParam1, cDepot.ModelParam2
if not cDepot.ModelParam1 in contourXaxis:
contourXaxis.append(cDepot.ModelParam1)
if not cDepot.ModelParam2 in contourYaxis:
contourYaxis.append(cDepot.ModelParam2)
ctp = cDepot.conturPoint
- print "adding", ctp.CLs, contourXaxis.index(cDepot.ModelParam1), contourYaxis.index(cDepot.ModelParam2)
+ #print "adding", ctp.CLs
confLim[contourXaxis.index(cDepot.ModelParam1)][contourYaxis.index(cDepot.ModelParam2)]=ctp.CLs
- #for ctp in cDepot.sortedPoints:
- # print ctp.tags
- # print ctp.pools
- # print ctp.CLs
-
+ for ctp in cDepot.sortedPoints:
+ #print ctp.tags
+ #print ctp.pools
+ #print ctp.CLs
+ confLims[ctp.pools][contourXaxis.index(cDepot.ModelParam1)][contourYaxis.index(cDepot.ModelParam2)] = ctp.CLs
+ #bestPlots[ctp.pools][contourXaxis.index(cDepot.ModelParam1)][contourYaxis.index(cDepot.ModelParam2)] = ctp.tags
+# set up for the plots here.
+
Xaxis = np.array(map(float, contourXaxis))
Yaxis = np.array(map(float, contourYaxis))
-#find the grid spacings:
+ # find the grid spacings:
dx= ( min(filter(lambda x: x> min(Xaxis), (Xaxis))) - min(Xaxis))/2.0
dy= ( min(filter(lambda x: x> min(Yaxis), (Yaxis))) - min(Yaxis))/2.0
print dx, dy
-# xx = []
-# yy = []
-# xx.append(Xaxis[0]-dx)
-
-# for i in range(1,len(Xaxis),1):
-# xx.append(Xaxis[i]+dx)
-
-# yy.append(Yaxis[0]-dy)
-# for j in range(1,len(Yaxis),1):
-# yy.append(Xaxis[i]+dy)
-
yy,xx = np.mgrid[min(Yaxis)-dy:max(Yaxis)+2*dy:2*dy,min(Xaxis)-dx:max(Xaxis)+2*dx:2*dx]
- yy=yy/10.
-
cl_values = confLim[:len(contourXaxis), :len(contourYaxis)]
-# for sina in contourYaxis:
-# sindex = contourYaxis.index(sina)
-# for mass in contourXaxis:
-# massindex = contourXaxis.index(mass)
-# print mass, float(sina)/10.0, cl_values[massindex,sindex]
+ # translate from the parameters of the to something more readable
+ fmt = plt.FuncFormatter(mh2)
+ fmt2 = plt.FuncFormatter(salpha)
+
+# make an html index file for browsing them
+
+ index = open("./plots/index.html", "w")
+
+ index.write('<html>\n<head>\n<title>Heatmaps and contours</title>\n</head>\n<body>')
+
+# make the overall heatmap
fig=plt.figure(figsize=fig_dims)
+ ax = fig.add_subplot(1,1,1)
+ ax.xaxis.set_major_formatter(fmt)
+ ax.yaxis.set_major_formatter(fmt2)
+
plt.pcolormesh(xx,yy,cl_values.T,cmap=plt.cm.magma, vmin=0, vmax=1)
-# plt.pcolormesh(cl_values,cmap=plt.cm.magma, vmin=0, vmax=1)
-# plt.axis([min(Xaxis), max(Xaxis), min(Yaxis)/10., max(Yaxis)/10.])
-##axis labels
- plt.xlabel(r"$M_{H2}$ [GeV]")
- plt.ylabel(r"$\sin \alpha$")
+ # axis labels
+ plt.xlabel(xAxisLabel)
+ plt.ylabel(yAxisLabel)
-##save the fig and pad it for better layout
+ # save the fig and pad it for better layout
fig.tight_layout(pad=0.1)
- plt.savefig("./plots/combinedCL.pdf")
- plt.savefig("./plots/combinedCL.png")
+ pngfile = "./plots/combinedCL.png"
+ pdffile = "./plots/combinedCL.pdf"
+ plt.savefig(pdffile)
+ plt.savefig(pngfile)
+ index.write('<h3>Combined Heatmap</h3>')
+ index.write('<img src="%s">\n' % os.path.basename(pngfile))
+
+ # now the heatmaps for the different analysis pools
+ for pool in confLims:
+
+ cl_values_pool = confLims[pool][:len(contourXaxis), :len(contourYaxis)]
+
+ fig=plt.figure(figsize=fig_dims)
+ ax = fig.add_subplot(1,1,1)
+ ax.xaxis.set_major_formatter(fmt)
+ ax.yaxis.set_major_formatter(fmt2)
+
+ plt.pcolormesh(xx,yy,cl_values_pool.T,cmap=plt.cm.magma, vmin=0, vmax=1)
+
+ # axis labels
+ plt.xlabel(xAxisLabel)
+ plt.ylabel(yAxisLabel)
+
+ # save the fig and pad it for better layout
+ fig.tight_layout(pad=0.1)
+ pngfile = "./plots/combinedCL_"+pool+".png"
+ plt.savefig(pngfile)
+ index.write('<h4>%s</h4>' % pool)
+ index.write('<img src="%s">\n' % os.path.basename(pngfile))
+
+
+
+# Now the overall contour plot ---------------------------------
fig=plt.figure(figsize=fig_dims)
-##draw a filled contour region for the CL excl
+ ax = fig.add_subplot(1,1,1)
+ ax.xaxis.set_major_formatter(fmt)
+ ax.yaxis.set_major_formatter(fmt2)
+
+ # draw a filled contour region for the CL excl
CS=plt.contourf(Xaxis,Yaxis,cl_values.T,levels=[0.95,1.0],label="CL",cmap=plt.cm.magma, alpha =0.8)
-##and a black outline
+ # and a black outline
CS2=plt.contour(CS, colors = 'black')
-##axis labels
- plt.xlabel(r"$M_{H2}$ [GeV]")
- plt.ylabel(r"$\sin \alpha$")
+ # axis labels
+ plt.xlabel(xAxisLabel)
+ plt.ylabel(yAxisLabel)
fig.tight_layout(pad=0.1)
plt.savefig("./plots/contur.pdf")
plt.savefig("./plots/contur.png")
-
-# print the colour bar key
+# Now the colour bar key ----------------------------------------
fig=plt.figure(figsize=[fig_dims[0]*2,0.5])
ax = fig.add_subplot(1,1,1)
import matplotlib as mpl
norm = mpl.colors.Normalize(vmin=0, vmax=1)
cb = mpl.colorbar.ColorbarBase(ax, cmap=plt.cm.magma, orientation='horizontal', norm=norm)
cb.set_label("CL of exclusion")
fig.tight_layout(pad=0.1)
plt.savefig('./plots/colorbarkey.pdf')
+ plt.savefig('./plots/colorbarkey.png')
+
+
+# close the html file
+ index.write("\n </body> \n")
+ index.close()
print "done"
diff --git a/AnalysisTools/contur/bin/contur-mkhtml b/AnalysisTools/contur/bin/contur-mkhtml
--- a/AnalysisTools/contur/bin/contur-mkhtml
+++ b/AnalysisTools/contur/bin/contur-mkhtml
@@ -1,491 +1,493 @@
#! /usr/bin/env python
"""\
%prog [options] <yodafile1> [<yodafile2> <yodafile3>...] [PLOT:Key1=Val1:...]
Make web pages from histogram files written out by Rivet and data files from Contur.
(Adapted from rivet-mkhtml.)
You can specify multiple Monte Carlo YODA files to be compared in the same syntax as for
rivet-cmphistos, i.e. including plotting options.
Reference data, analysis metadata, and plot style information should be found
automatically (if not, set the RIVET_ANALYSIS_PATH or similar variables
appropriately).
Contur data is assumed to be in the directory ./plots and ./ANALYSIS.
Output willbe written to ./contur-plots
Any existing output directory will be overwritten.
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
"""
import rivet, sys, os
rivet.util.check_python_version()
rivet.util.set_process_name(os.path.basename(__file__))
import glob, shutil
from subprocess import Popen,PIPE
from optparse import OptionParser, OptionGroup
parser = OptionParser(usage=__doc__)
parser.add_option("-o", "--outputdir", dest="OUTPUTDIR",
default="./contur-plots", help="directory for Web page output")
parser.add_option("-c", "--config", dest="CONFIGFILES", action="append", default=["~/.make-plots"],
help="plot config file(s) to be used with rivet-cmphistos")
parser.add_option("-n", "--num-threads", metavar="NUMTHREADS", dest="NUMTHREADS", type=int,
default=None, help="request make-plots to use a specific number of threads")
parser.add_option("--ignore-missing", dest="IGNORE_MISSING", action="store_true",
default=False, help="ignore missing YODA files")
parser.add_option("-i", "--ignore-unvalidated", dest="IGNORE_UNVALIDATED", action="store_true",
default=False, help="ignore unvalidated analyses")
# parser.add_option("--ref", "--refid", dest="REF_ID",
# default=None, help="ID of reference data set (file path for non-REF data)")
parser.add_option("--dry-run", help="don't actually do any plotting or HTML building", dest="DRY_RUN",
action="store_true", default=False)
parser.add_option("--no-cleanup", dest="NO_CLEANUP", action="store_true", default=False,
help="keep plotting temporary directory")
parser.add_option("--no-subproc", dest="NO_SUBPROC", action="store_true", default=False,
help="don't use subprocesses to render the plots in parallel -- useful for debugging")
parser.add_option("--pwd", dest="PATH_PWD", action="store_true", default=False,
help="append the current directory (pwd) to the analysis/data search paths (cf. $RIVET_ANALYSIS_PATH)")
stygroup = OptionGroup(parser, "Style options")
stygroup.add_option("-t", "--title", dest="TITLE",
default="Contur Plots: Constraints On New Theories Using Rivet", help="title to be displayed on the main web page")
stygroup.add_option("--reftitle", dest="REFTITLE",
default="Data", help="legend entry for reference data")
stygroup.add_option("--no-plottitle", dest="NOPLOTTITLE", action="store_true",
default=False, help="don't show the plot title on the plot "
"(useful when the plot description should only be given in a caption)")
stygroup.add_option("-s", "--single", dest="SINGLE", action="store_true",
default=False, help="display plots on single webpage.")
stygroup.add_option("--no-ratio", dest="SHOW_RATIO", action="store_false",
default=True, help="don't draw a ratio plot under each main plot.")
stygroup.add_option("--errs", "--mc-errs", dest="MC_ERRS", action="store_true",
default=False, help="plot error bars.")
stygroup.add_option("--offline", dest="OFFLINE", action="store_true",
default=False, help="generate HTML that does not use external URLs.")
stygroup.add_option("--pdf", dest="VECTORFORMAT", action="store_const", const="PDF",
default="PDF", help="use PDF as the vector plot format.")
stygroup.add_option("--ps", dest="VECTORFORMAT", action="store_const", const="PS",
default="PDF", help="use PostScript as the vector plot format. DEPRECATED")
stygroup.add_option("--booklet", dest="BOOKLET", action="store_true",
default=False, help="create booklet (currently only available for PDF with pdftk or pdfmerge).")
stygroup.add_option("--font", dest="OUTPUT_FONT", choices="palatino,cm,times,helvetica,minion".split(","),
default="palatino", help="choose the font to be used in the plots")
stygroup.add_option("--palatino", dest="OUTPUT_FONT", action="store_const", const="palatino", default="palatino",
help="use Palatino as font (default). DEPRECATED: Use --font")
stygroup.add_option("--cm", dest="OUTPUT_FONT", action="store_const", const="cm", default="palatino",
help="use Computer Modern as font. DEPRECATED: Use --font")
stygroup.add_option("--times", dest="OUTPUT_FONT", action="store_const", const="times", default="palatino",
help="use Times as font. DEPRECATED: Use --font")
stygroup.add_option("--helvetica", dest="OUTPUT_FONT", action="store_const", const="helvetica", default="palatino",
help="use Helvetica as font. DEPRECATED: Use --font")
stygroup.add_option("--minion", dest="OUTPUT_FONT", action="store_const", const="minion", default="palatino",
help="use Adobe Minion Pro as font. DEPRECATED: Use --font")
parser.add_option_group(stygroup)
selgroup = OptionGroup(parser, "Selective plotting")
selgroup.add_option("-m", "--match", action="append", dest="PATHPATTERNS", default=[],
help="only write out histograms whose $path/$name string matches any of these regexes")
selgroup.add_option("-M", "--unmatch", action="append", dest="PATHUNPATTERNS", default=[],
help="exclude histograms whose $path/$name string matches any of these regexes")
selgroup.add_option("-a", "--ana-match", action="append", dest="ANAPATTERNS", default=[],
help="only write out histograms from analyses whose name matches any of these regexes")
selgroup.add_option("-A", "--ana-unmatch", action="append", dest="ANAUNPATTERNS", default=[],
help="exclude histograms from analyses whose name matches any of these regexes")
parser.add_option_group(selgroup)
vrbgroup = OptionGroup(parser, "Verbosity control")
vrbgroup.add_option("-v", "--verbose", help="add extra debug messages", dest="VERBOSE",
action="store_true", default=False)
parser.add_option_group(vrbgroup)
opts, yodafiles = parser.parse_args()
## Add pwd to search paths
if opts.PATH_PWD:
rivet.addAnalysisLibPath(os.path.abspath("."))
rivet.addAnalysisDataPath(os.path.abspath("."))
## Check that there are some arguments!
if not yodafiles:
print "Error: You need to specify some YODA files to be plotted!"
sys.exit(1)
## Make output directory
if not opts.DRY_RUN:
if os.path.exists(opts.OUTPUTDIR) and not os.path.realpath(opts.OUTPUTDIR)==os.getcwd():
import shutil
shutil.rmtree(opts.OUTPUTDIR)
try:
os.makedirs(opts.OUTPUTDIR)
except:
print "Error: failed to make new directory '%s'" % opts.OUTPUTDIR
sys.exit(1)
## Get set of analyses involved in the runs
plotarg = None
analyses = set()
blocked_analyses = set()
import yoda
for yodafile in yodafiles:
if yodafile.startswith("PLOT:"):
plotarg = yodafile
continue
yodafilepath = os.path.abspath(yodafile.split(":")[0])
if not os.access(yodafilepath, os.R_OK):
print "Error: cannot read from %s" % yodafilepath
if opts.IGNORE_MISSING:
continue
else:
sys.exit(2)
try:
## Note: we use -m/-M flags here as well as when calling rivet-cmphistos, to potentially speed this initial loading
analysisobjects = yoda.read(yodafilepath, patterns=opts.PATHPATTERNS, unpatterns=opts.PATHUNPATTERNS)
except IOError, e:
print "File reading error: ", e.strerror
sys.exit(1)
for path, ao in analysisobjects.iteritems():
## Make a path object and ensure the path is in standard form
try:
aop = rivet.AOPath(path)
except Exception, e:
#print e
print "Found analysis object with non-standard path structure:", path, "... skipping"
continue
## We don't plot data objects with path components hidden by an underscore prefix
if aop.istmp():
continue
## Identify analysis/histo name parts
analysis = "ANALYSIS"
if aop.basepathparts(keepref=False):
analysis = aop.basepathparts(keepref=False)[0] #< TODO: for compatibility with rivet-cmphistos... generalise?
#analysis = "_".join(aop.dirnameparts(keepref=False)[:-1]) #< TODO: would this be nicer? Currently incompatible with rivet-cmphistos
## Optionally veto on analysis name pattern matching
if analysis in analyses.union(blocked_analyses):
continue
import re
matched = True
if opts.ANAPATTERNS:
matched = False
for patt in opts.ANAPATTERNS:
if re.match(patt, analysis) is not None:
matched = True
break
if matched and opts.ANAUNPATTERNS:
for patt in opts.ANAUNPATTERNS:
if re.match(patt, analysis) is not None:
matched = False
break
if matched:
analyses.add(analysis)
else:
blocked_analyses.add(analysis)
## Sort analyses: group ascending by analysis name (could specialise grouping by collider), then
## descending by year, and finally descending by bibliographic archive ID code (INSPIRE first).
def anasort(name):
rtn = (1, name)
if name.startswith("MC"):
rtn = (99999999, name)
else:
stdparts = name.split("_")
try:
year = int(stdparts[1])
rtn = (0, stdparts[0], -year, 0)
idcode = (0 if stdparts[2][0] == "I" else 1e10) - int(stdparts[2][1:])
rtn = (0, stdparts[0], -year, idcode)
if len(stdparts) > 3:
rtn += stdparts[3:]
except:
pass
return rtn
analyses = sorted(analyses, key=anasort)
## Uncomment to test analysis ordering on index page
# print analyses
# sys.exit(0)
## Write web page containing all (matched) plots
## Make web pages first so that we can load it locally in
## a browser to view the output before all plots are made
if not opts.DRY_RUN:
style = '''<style>
html { font-family: sans-serif; }
img { border: 0; }
a { text-decoration: none; font-weight: bold; }
</style>
'''
## Include MathJax configuration
script = ''
if not opts.OFFLINE:
script = '''\
<script type="text/x-mathjax-config">
MathJax.Hub.Config({
tex2jax: {inlineMath: [["$","$"]]}
});
</script>
<script type="text/javascript"
src="https://cdn.mathjax.org/mathjax/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML">
</script>
'''
## A helper function for metadata LaTeX -> HTML conversion
from rivet.util import htmlify
## A timestamp HTML fragment to be used on each page:
import datetime
timestamp = '<p>Generated at %s</p>\n' % datetime.datetime.now().strftime("%A, %d. %B %Y %I:%M%p")
# Open the top-level index file
index = open(os.path.join(opts.OUTPUTDIR, "index.html"), "w")
# Write title
index.write('<html>\n<head>\n<title>%s</title>\n%s</head>\n<body>' % (opts.TITLE, style + script))
if opts.BOOKLET and opts.VECTORFORMAT == "PDF":
index.write('<h2><a href="booklet.pdf">%s</a></h2>\n\n' % opts.TITLE)
else:
index.write('<h2>%s</h2>\n\n' % opts.TITLE)
try:
summary = open("./ANALYSIS/Summary.txt","r")
index.write("<p>")
heading=True
while heading==True:
text = summary.readline()
if "pools" in text:
heading=False
continue
index.write("%s\n" % text)
index.write("</br>")
index.write("</p>")
index.write("<h2>In each analysis pool, these plots contributed:</h2>\n")
for line in summary:
list = line.split(",")
if list[0][0] == "{":
# This is a line containing extra info which we don't yet use.
continue
for name in list:
parts = name.split("/")
if len(parts) == 1:
# This is a line identifying the analysis pool
# TODO would be good to add a human-readable anapool name to the DB.
index.write('<div style="float:none; overflow:auto; width:100%">\n')
index.write("<h3>Analysis Pool:"+str(parts[0])+"</h3>")
else:
# This is a line identifying histograms
obsname = str(parts[1])+"/"+str(parts[1])+"_"+str(parts[2])
pngfile = obsname+".png"
vecfile = obsname+"."+opts.VECTORFORMAT.lower()
srcfile = obsname+".dat"
index.write(' <div style="float:left; font-size:smaller; font-weight:bold;">\n')
index.write(' <a href="#%s-%s">⚓</a><a href="%s">⌘</a> %s:<br/>\n' %
(analysis, obsname, srcfile, os.path.splitext(vecfile)[0]) )
index.write(' <a name="%s-%s"><a href="%s">\n' % (analysis, obsname, vecfile) )
index.write(' <img src="%s">\n' % pngfile )
index.write(' </a></a>\n')
index.write(' </div>\n')
summary.close()
index.write('</div>\n')
except IOError, e:
print "No Contur summary file found: Summary info omitted"
index.write("<hr/><h3>List of all the analyses which were tried:</h3>\n")
if opts.SINGLE:
## Write table of contents
index.write('<ul>\n')
for analysis in analyses:
summary = analysis
ana = rivet.AnalysisLoader.getAnalysis(analysis)
if ana:
summary = "%s (%s)" % (ana.summary(), analysis)
if opts.IGNORE_UNVALIDATED and ana.status() != "VALIDATED":
continue
index.write('<li><a href="#%s">%s</a>\n' % (analysis, htmlify(summary)) )
index.write('</ul>\n')
for analysis in analyses:
references = []
summary = htmlify(analysis)
description, inspireid, spiresid = None, None, None
if analysis.find("_I") > 0:
inspireid = analysis[analysis.rfind('_I')+2:len(analysis)]
elif analysis.find("_S") > 0:
spiresid = analysis[analysis.rfind('_S')+2:len(analysis)]
ana = rivet.AnalysisLoader.getAnalysis(analysis)
if ana:
if ana.summary():
summary = htmlify("%s (%s)" % (ana.summary(), analysis))
references = ana.references()
description = htmlify(ana.description())
spiresid = ana.spiresId()
if opts.IGNORE_UNVALIDATED and ana.status().upper() != "VALIDATED":
continue
if opts.SINGLE:
index.write('\n<h3 style="clear:left; padding-top:2em;"><a name="%s">%s</a></h3>\n' % (analysis, summary))
else:
index.write('\n<h3><a href="%s/index.html" style="text-decoration:none;">%s</a></h3>\n' % (analysis, summary))
reflist = []
if inspireid:
reflist.append('<a href="http://inspire-hep.net/record/%s">Inspire</a>' % inspireid)
reflist.append('<a href="http://hepdata.cedar.ac.uk/view/ins%s">HepData</a>' % inspireid)
elif spiresid:
# elif ana.spiresId():
reflist.append('<a href="http://inspire-hep.net/search?p=find+key+%s">Inspire</a>' % spiresid)
reflist.append('<a href="http://hepdata.cedar.ac.uk/view/irn%s">HepData</a>' % spiresid)
reflist += references
index.write('<p>%s</p>\n' % " | ".join(reflist))
if description:
#print inspireid
#print description
index.write('<p style="font-size:smaller;">%s</p>\n' % description)
anapath = os.path.join(opts.OUTPUTDIR, analysis)
if not opts.SINGLE:
if not os.path.exists(anapath):
os.makedirs(anapath)
anaindex = open(os.path.join(anapath, "index.html"), 'w')
anaindex.write('<html>\n<head>\n<title>%s – %s</title>\n%s</head>\n<body>\n' %
(htmlify(opts.TITLE), analysis, style + script))
anaindex.write('<h3>%s</h3>\n' % htmlify(analysis))
anaindex.write('<p><a href="../index.html">Back to index</a></p>\n')
if description:
anaindex.write('<p>\n %s\n</p>\n' % description)
else:
anaindex = index
datfiles = glob.glob("./plots/%s*.dat" % analysis)
anaindex.write('<div style="float:none; overflow:auto; width:100%">\n')
for datfile in sorted(datfiles):
obsname = os.path.basename(datfile).replace(".dat", "")
pngfile = obsname+".png"
vecfile = obsname+"."+opts.VECTORFORMAT.lower()
srcfile = obsname+".dat"
if opts.SINGLE:
pngfile = os.path.join(analysis, pngfile)
vecfile = os.path.join(analysis, vecfile)
srcfile = os.path.join(analysis, srcfile)
anaindex.write(' <div style="float:left; font-size:smaller; font-weight:bold;">\n')
anaindex.write(' <a href="#%s-%s">⚓</a><a href="%s">⌘</a> %s:<br/>\n' %
(analysis, obsname, srcfile, os.path.splitext(vecfile)[0]) )
anaindex.write(' <a name="%s-%s"><a href="%s">\n' % (analysis, obsname, vecfile) )
anaindex.write(' <img src="%s">\n' % pngfile )
anaindex.write(' </a></a>\n')
anaindex.write(' </div>\n')
anaindex.write('</div>\n')
if not opts.SINGLE:
anaindex.write('<div style="float:none">%s</body>\n</html></div>\n' % timestamp)
anaindex.close()
index.write('<br>%s</body>\n</html>' % timestamp)
index.close()
# http://stackoverflow.com/questions/377017/test-if-executable-exists-in-python
def which(program):
import os
def is_exe(fpath):
return os.path.isfile(fpath) and os.access(fpath, os.X_OK)
fpath, fname = os.path.split(program)
if fpath:
if is_exe(program):
return program
else:
for path in os.environ["PATH"].split(os.pathsep):
path = path.strip('"')
exe_file = os.path.join(path, program)
if is_exe(exe_file):
return exe_file
return None
## Run make-plots on all generated .dat files
# sys.exit(0)
mp_cmd = ["make-plots"]
if opts.NUMTHREADS:
mp_cmd.append("--num-threads=%d" % opts.NUMTHREADS)
if opts.NO_CLEANUP:
mp_cmd.append("--no-cleanup")
if opts.NO_SUBPROC:
mp_cmd.append("--no-subproc")
if opts.VECTORFORMAT == "PDF":
mp_cmd.append("--pdfpng")
elif opts.VECTORFORMAT == "PS":
mp_cmd.append("--pspng")
if opts.OUTPUT_FONT:
mp_cmd.append("--font=%s" % opts.OUTPUT_FONT)
+if opts.OUTPUT_FONT:
+ mp_cmd.append("--font=%s" % opts.OUTPUT_FONT)
datfiles = []
# Copy the dat files from the plots directory into the contur-plot output directory
for analysis in analyses:
anapath = os.path.join(opts.OUTPUTDIR, analysis)
anadatfiles = glob.glob("./plots/%s_d*.dat" % analysis)
for anadatfile in anadatfiles:
shutil.copyfile(anadatfile,"./contur-plots/%s/%s" % (analysis, os.path.basename(anadatfile)))
anadatfiles = glob.glob("./contur-plots/%s/%s*.dat" % (analysis, analysis))
datfiles += sorted(anadatfiles)
# Make the PDF and PNG files
if datfiles:
mp_cmd += datfiles
if opts.VERBOSE:
mp_cmd.append("--verbose")
print "Calling make-plots with the following options:"
print " ".join(mp_cmd)
if not opts.DRY_RUN:
Popen(mp_cmd).wait()
if opts.BOOKLET and opts.VECTORFORMAT=="PDF":
if which("pdftk") is not None:
bookletcmd = ["pdftk"]
for analysis in analyses:
anapath = os.path.join(opts.OUTPUTDIR, analysis)
bookletcmd += sorted(glob.glob("%s/*.pdf" % anapath))
bookletcmd += ["cat", "output", "%s/booklet.pdf" % opts.OUTPUTDIR]
print bookletcmd
Popen(bookletcmd).wait()
elif which("pdfmerge") is not None:
bookletcmd = ["pdfmerge"]
for analysis in analyses:
anapath = os.path.join(opts.OUTPUTDIR, analysis)
bookletcmd += sorted(glob.glob("%s/*.pdf" % anapath))
bookletcmd += ["%s/booklet.pdf" % opts.OUTPUTDIR]
print bookletcmd
Popen(bookletcmd).wait()
else:
print "Neither pdftk nor pdfmerge available --- not booklet output possible"
diff --git a/AnalysisTools/contur/contur/Utils/Utils.py b/AnalysisTools/contur/contur/Utils/Utils.py
--- a/AnalysisTools/contur/contur/Utils/Utils.py
+++ b/AnalysisTools/contur/contur/Utils/Utils.py
@@ -1,184 +1,186 @@
import os
import yoda
import rivet
import plotinfo
from contur import TestingFunctions as ctr
#set the version here globally
global version
version = " Beta pre-release"
def writeBanner():
"""Write Ascii banner"""
print "Running Contur ",version, "\n"
print "Mail to developers: contur (at) projects.hepforge.org \n"
def mkoutdir(outdir):
"Function to make output directories"
if not os.path.exists(outdir):
try:
os.makedirs(outdir)
except:
msg = "Can't make output directory '%s'" % outdir
raise Exception(msg)
if not os.access(outdir, os.W_OK):
msg = "Can't write to output directory '%s'" % outdir
raise Exception(msg)
def writeOutput(output, h):
mkoutdir("ANALYSIS")
f = open("./ANALYSIS/"+h, 'w')
for item in output:
f.write(str(item) + "\n")
f.close()
def getHistos(filelist):
"""Loop over all input files. Only use the first occurrence of any REF-histogram
and the first occurrence in each MC file for every MC-histogram."""
# return reference histograms,
# mchistos = mc plots,
# xsec = generated xsec and its uncertainty,
# Nev = sum of weights, sum of squared weight, number of events
# (derived from rivet-cmphistos)
refhistos = {}
mchistos = {}
xsec = {}
Nev = {}
# for infile in filelist:
mchistos.setdefault(filelist, {})
analysisobjects = yoda.read(filelist)
print len(analysisobjects), "analysisobjects in", filelist
for path, ao in analysisobjects.iteritems():
if path.startswith('/_EVTCOUNT'):
Nev = ao
if path.startswith('/_XSEC'):
xsec = ao
# Conventionally don't plot data objects whose names start with an
# underscore
if os.path.basename(path).startswith("_"):
continue
if path.startswith('/REF/'):
if path not in refhistos:
refhistos[path] = ao
else:
if path not in mchistos[filelist]:
mchistos[filelist][path] = ao
return refhistos, mchistos, xsec, Nev
def writeHistoDat(mcpath, plotparser, outdir, nostack, histo):
"""Write a .dat file for the histogram in the output directory, for later display."""
anaobjects = []
drawonly = []
mcpath = "/"+mcpath
## Check if we have reference data for the histogram
ratioreference = None
if histo.ref:
refdata = histo.refplot
background = histo.bgplot
if nostack:
sigback = histo.sigplot
else:
sigback = histo.stack
h = sigback.path
sigback.setAnnotation('Path', mcpath+h)
refdata.setAnnotation('ErrorBars', '1')
refdata.setAnnotation('PolyMarker', '*')
refdata.setAnnotation('ConnectBins', '0')
refdata.setAnnotation('Title', 'Data')
background.setAnnotation('LineColor', 'green')
anaobjects.append(background)
anaobjects.append(refdata)
drawonly.append('/REF' + h)
drawonly.append(mcpath + h)
drawonly.append('/THY' + h)
# write the bin number of the most significant bin, and the bin number for the plot legend
if histo.maxbin > 0:
sigback.title='[%s] %5.2f' % ( histo.maxbin, histo.conturPoints[histo.maxcl].CLs )
sigback.setAnnotation('LineColor', 'red')
+ sigback.setAnnotation('ErrorBars', '1')
+
anaobjects.append(sigback)
plot = plotinfo.Plot()
plot['DrawOnly'] = ' '.join(drawonly).strip()
plot['Legend'] = '1'
plot['MainPlot'] = '1'
plot['RatioPlotYMin'] = '1'
plot['LogY'] = '1'
plot['RatioPlot'] = '1'
for key, val in plotparser.getHeaders(h).iteritems():
# Get any updated attributes from plotinfo files
plot[key] = val
ratioreference = '/REF'+h
plot['RatioPlotReference'] = ratioreference
output = ''
output += str(plot)
from cStringIO import StringIO
sio = StringIO()
yoda.writeFLAT(anaobjects, sio)
output += sio.getvalue()
hparts = h.strip("/").split("/")
outfile = '%s.dat' % "_".join(hparts)
mkoutdir(outdir)
outfilepath = os.path.join(outdir, outfile)
f = open(outfilepath, 'w')
f.write(output)
f.close()
def mkScatter2D(s1):
""" Make a Scatter2D from a Scatter1D by treating the points as y values and adding dummy x bins."""
rtn = yoda.Scatter2D()
x = 0.5
for a in s1.annotations:
rtn.setAnnotation(a, s1.annotation(a))
rtn.setAnnotation("Type", "Scatter2D");
for point in s1.points:
ex_m = x-0.5;
ex_p = x+0.5;
y = point.x;
ey_p = point.xMax;
ey_m = point.xMin;
pt = yoda.Point2D(x, y, 0.5, ey_m);
rtn.addPoint(pt);
x = x + 1.0
return rtn;
def walklevel(some_dir, level=1):
"""Like os.walk but can specify a level to walk to
useful for managing directories a bit better
https://stackoverflow.com/questions/229186/os-walk-without-digging-into-directories-below
"""
some_dir = some_dir.rstrip(os.path.sep)
assert os.path.isdir(some_dir)
num_sep = some_dir.count(os.path.sep)
for root, dirs, files in os.walk(some_dir):
yield root, dirs, files
num_sep_this = root.count(os.path.sep)
if num_sep + level <= num_sep_this:
del dirs[:]