adds python tools for examining the output of a filval program

.flake8

@@ -0,0 +1,4 @@
+
+[flake8]
+ignore = E248, E241, E226
+max_line_length=120

.gitignore

@@ -0,0 +1,5 @@
+figures/
+env/
+
+__pycache__/
+.ipynb_checkpoints/

canvas_wrapper.py

@@ -0,0 +1,66 @@
+"""
+Helper module for displaying ROOT canvases in ipython notebooks
+
+Usage example:
+# Save this file as rootnotes.py to your working directory.
+import rootnotes
+c1 = rootnotes.canvas()
+fun1 = TF1( 'fun1', 'abs(sin(x)/x)', 0, 10)
+c1.SetGridx()
+c1.SetGridy()
+fun1.Draw()
+c1
+
+More examples: http://mazurov.github.io/webfest2013/
+
+@author alexander.mazurov@cern.ch
+@author andrey.ustyuzhanin@cern.ch
+@date 2013-08-09
+"""
+
+import ROOT
+ROOT.gROOT.SetBatch()
+
+import tempfile
+from IPython.core import display
+
+
+def canvas(name="icanvas", size=(800, 600)):
+    """Helper method for creating canvas"""
+
+    assert len(size) == 2
+    if len(size) != 2:
+        raise ValueError("Size must be a 2 element tuple")
+    # Check if icanvas already exists
+    canvas = ROOT.gROOT.FindObject(name)
+    if canvas:
+        return canvas
+    else:
+        return ROOT.TCanvas(name, name, size[0], size[1])
+
+
+def _display_canvas(canvas):
+    file = tempfile.NamedTemporaryFile(suffix=".png")
+    canvas.SaveAs(file.name)
+    ip_img = display.Image(filename=file.name, format='png', embed=True)
+    return ip_img._repr_png_()
+
+
+def _display_any(obj):
+    file = tempfile.NamedTemporaryFile(suffix=".png")
+    obj.Draw()
+    ROOT.gPad.SaveAs(file.name)
+    ip_img = display.Image(filename=file.name, format='png', embed=True)
+    return ip_img._repr_png_()
+
+
+# register display function with PNG formatter:
+png_formatter = get_ipython().display_formatter.formatters['image/png'] # noqa
+
+# Register ROOT types in ipython
+#
+# In [1]: canvas = canvas_wrapper.canvas()
+# In [2]: canvas
+# Out [2]: [image will be here]
+png_formatter.for_type(ROOT.TCanvas, _display_canvas)
+png_formatter.for_type(ROOT.TF1, _display_any)

graph_vals.py

@@ -0,0 +1,101 @@
+import pydotplus.graphviz as pdp
+import sys
+import re
+
+
+def parse(str_in, alias=None):
+    str_in = "("+str_in+")"
+
+    functions = []
+    ends = {}
+    nests = {}
+    names = {}
+    styles = {}
+    parens = []
+    name = ""
+    name_start = 0
+    name_end = 0
+    for i, char in enumerate(str_in):
+        if char == "(":
+            nests[name_start] = []
+            if parens:
+                nests[parens[-1]].append(name_start)
+            names[name_start] = name  # save function name
+            styles[name_start] = {"shape": "ellipse"}
+            name = ""
+            parens.append(name_start)
+        elif char == ")":
+            if name:
+                ends[name_start] = name_end
+                names[name_start] = name
+                styles[name_start] = {"shape": "rectangle"}
+                nests[parens[-1]].append(name_start)
+                name = ""
+            ends[parens.pop()] = i
+        elif char in ",:":
+            if name:
+                ends[name_start] = name_end
+                names[name_start] = name
+                if char == ",":
+                    styles[name_start] = {"shape": "rectangle"}
+                else:
+                    styles[name_start] = {"shape": "invhouse"}
+                    functions.append(name)
+                nests[parens[-1]].append(name_start)
+                name = ""
+        else:
+            if not name:
+                name_start = i
+            name += char
+            name_end = i
+
+    # clean up duplicate sub-trees
+    text = {}
+    for start, end in ends.items():
+        s = str_in[start:end+1]
+        if s in text:
+            dup_id = text[s]
+            names.pop(start)
+            if start in nests:
+                nests.pop(start)
+            for l in nests.values():
+                for i in range(len(l)):
+                    if l[i] == start:
+                        l[i] = dup_id
+        else:
+            text[s] = start
+
+    names.pop(0)
+    nests.pop(0)
+
+    g = pdp.Dot()
+
+    for id_, name in names.items():
+        g.add_node(pdp.Node(str(id_), label=name, **styles[id_]))
+    for group_id, children in nests.items():
+        for child_id in children:
+            g.add_edge(pdp.Edge(str(group_id), str(child_id)))
+    if alias:
+        g.add_node(pdp.Node(alias, shape="plain", pos="0,0!"))
+    return g, functions
+
+
+if __name__ == '__main__':
+    aliases = {}
+    ali_re = re.compile(r"ALIAS::\"([^\"]*)\" referring to \"([^\"]*)\"")
+
+    with open(sys.argv[1]) as f:
+        for line in f.readlines():
+            res = ali_re.findall(line)
+            if res:
+                aliases[res[0][1]] = res[0][0]
+                continue
+    for name, alias in aliases.items():
+        graph, _ = parse(name, alias)
+        fname = "val_graph_{}.gif".format(alias)
+        with open(fname, "wb") as f:
+            try:
+                f.write(graph.create_gif())
+            except Exception as e:
+                print(e)
+                print(graph.to_string())

plotter.py

@@ -0,0 +1,248 @@
+#!/usr/bin/env python3
+import math
+import matplotlib as mpl
+mpl.rc('font', **{'family': 'sans-serif', 'sans-serif': ['Helvetica']})
+mpl.rc('font', **{'family': 'serif', 'serif': ['Palatino']})
+mpl.rc('text', usetex=True)
+mpl.rc('savefig', dpi=120)
+
+
+def add_decorations(axes, luminosity, energy):
+    cms_prelim = r'{\raggedright{}\textsf{\textbf{CMS}}\\ \emph{Preliminary}}'
+    axes.text(0.01, 0.98, cms_prelim,
+              horizontalalignment='left',
+              verticalalignment='top',
+              transform=axes.transAxes)
+
+    lumi = ""
+    energy_str = ""
+    if luminosity is not None:
+        lumi = r'${} \mathrm{{fb}}^{{-1}}$'.format(luminosity)
+    if energy is not None:
+        energy_str = r'({} TeV)'.format(energy)
+
+    axes.text(1, 1, ' '.join([lumi, energy_str]),
+              horizontalalignment='right',
+              verticalalignment='bottom',
+              transform=axes.transAxes)
+
+
+class StackHist:
+
+    def __init__(self, title=""):
+        self.title = title
+        self.xlabel = ""
+        self.ylabel = ""
+        self.xlim = (None, None)
+        self.ylim = (None, None)
+        self.logx = False
+        self.logy = False
+        self.backgrounds = []
+        self.signal = None
+        self.signal_stack = True
+        self.data = None
+
+    @staticmethod
+    def to_bin_list(th1):
+        bins = []
+        for i in range(th1.GetNbinsX()):
+            center = th1.GetBinCenter(i + 1)
+            width = th1.GetBinWidth(i + 1)
+            content = th1.GetBinContent(i + 1)
+            bins.append((center-width/2, center+width/2, content))
+        return bins
+
+    def add_mc_background(self, th1, label, lumi=None, plot_color=''):
+        self.backgrounds.append((label, lumi, self.to_bin_list(th1), plot_color))
+
+    def set_mc_signal(self, th1, label, lumi=None, stack=True, scale=1, plot_color=''):
+        self.signal = (label, lumi, self.to_bin_list(th1), plot_color)
+        self.signal_stack = stack
+        self.signal_scale = scale
+
+    def set_data(self, th1, lumi=None, plot_color=''):
+        self.data = ('data', lumi, self.to_bin_list(th1), plot_color)
+        self.luminosity = lumi
+
+    def _verify_binning_match(self):
+        bins_count = [len(bins) for _, _, bins, _ in self.backgrounds]
+        if self.signal is not None:
+            bins_count.append(len(self.signal[2]))
+        if self.data is not None:
+            bins_count.append(len(self.data[2]))
+        n_bins = bins_count[0]
+        if any(bin_count != n_bins for bin_count in bins_count):
+            raise ValueError("all histograms must have the same number of bins")
+        self.n_bins = n_bins
+
+    def save(self, filename, **kwargs):
+        import matplotlib.pyplot as plt
+        plt.ioff()
+        fig = plt.figure()
+        ax = fig.gca()
+        self.do_draw(ax, **kwargs)
+        fig.savefig("figures/"+filename, transparent=True)
+        plt.close(fig)
+        plt.ion()
+
+    def do_draw(self, axes):
+        self.axeses = [axes]
+        self._verify_binning_match()
+        bottoms = [0]*self.n_bins
+
+        if self.logx:
+            axes.set_xscale('log')
+        if self.logy:
+            axes.set_yscale('log')
+
+        def draw_bar(label, lumi, bins, plot_color, scale=1, stack=True, **kwargs):
+            if stack:
+                lefts = []
+                widths = []
+                heights = []
+                for left, right, content in bins:
+                    lefts.append(left)
+                    widths.append(right-left)
+                    if lumi is not None:
+                        content *= self.luminosity/lumi
+                    content *= scale
+                    heights.append(content)
+
+                axes.bar(lefts, heights, widths, bottoms, label=label, color=plot_color, **kwargs)
+                for i, (_, _, content) in enumerate(bins):
+                    if lumi is not None:
+                        content *= self.luminosity/lumi
+                    content *= scale
+                    bottoms[i] += content
+            else:
+                xs = [bins[0][0] - (bins[0][1]-bins[0][0])/2]
+                ys = [0]
+                for left, right, content in bins:
+                    width2 = (right-left)/2
+                    if lumi is not None:
+                        content *= self.luminosity/lumi
+                    content *= scale
+                    xs.append(left-width2)
+                    ys.append(content)
+                    xs.append(right-width2)
+                    ys.append(content)
+                xs.append(bins[-1][0] + (bins[-1][1]-bins[-1][0])/2)
+                ys.append(0)
+                axes.plot(xs, ys, label=label, color=plot_color, **kwargs)
+
+        if self.signal is not None and self.signal_stack:
+            label, lumi, bins, plot_color = self.signal
+            if self.signal_scale != 1:
+                label = r"{}$\times{:d}$".format(label, self.signal_scale)
+            draw_bar(label, lumi, bins, plot_color, scale=self.signal_scale, hatch='/')
+
+        for background in self.backgrounds:
+            draw_bar(*background)
+
+        if self.signal is not None and not self.signal_stack:
+            # draw_bar(*self.signal, stack=False, color='k')
+            label, lumi, bins, plot_color = self.signal
+            if self.signal_scale != 1:
+                label = r"{}$\times{:d}$".format(label, self.signal_scale)
+            draw_bar(label, lumi, bins, plot_color, scale=self.signal_scale, stack=False)
+
+        axes.set_title(self.title)
+        axes.set_xlabel(self.xlabel)
+        axes.set_ylabel(self.ylabel)
+        axes.set_xlim(*self.xlim)
+        # axes.set_ylim(*self.ylim)
+        if self.logy:
+            axes.set_ylim(None, math.exp(math.log(max(bottoms))*1.4))
+        else:
+            axes.set_ylim(None, max(bottoms)*1.2)
+        axes.legend(frameon=True, ncol=2)
+        add_decorations(axes, self.luminosity, self.energy)
+
+    def draw(self, axes, save=False, filename=None, **kwargs):
+        self.do_draw(axes, **kwargs)
+        if save:
+            if filename is None:
+                filename = "".join(c for c in self.title if c.isalnum() or c in (' ._+-'))+".png"
+            self.save(filename, **kwargs)
+
+
+class StackHistWithSignificance(StackHist):
+
+    def __init__(self, *args, **kwargs):
+        super().__init__(*args, **kwargs)
+
+    def do_draw(self, axes, bin_significance=True, low_cut_significance=False, high_cut_significance=False):
+        bottom_box, _, top_box = axes.get_position().splity(0.28, 0.30)
+        axes.set_position(top_box)
+        super().do_draw(axes)
+        axes.set_xticks([])
+        rhs_color = '#cc6600'
+
+        bottom = axes.get_figure().add_axes(bottom_box)
+        bottom_rhs = bottom.twinx()
+        bgs = [0]*self.n_bins
+        for (_, _, bins, _) in self.backgrounds:
+            for i, (left, right, value) in enumerate(bins):
+                bgs[i] += value
+
+        sigs = [0]*self.n_bins
+        if bin_significance:
+            xs = []
+            for i, (left, right, value) in enumerate(self.signal[2]):
+                sigs[i] += value
+                xs.append(left)
+            xs, ys = zip(*[(x, sig/(sig+bg)) for x, sig, bg in zip(xs, sigs, bgs) if (sig+bg)>0])
+            bottom.plot(xs, ys, '.k')
+
+        if high_cut_significance:
+            # s/(s+b) for events passing a minimum cut requirement
+            min_bg = [sum(bgs[i:]) for i in range(self.n_bins)]
+            min_sig = [sum(sigs[i:]) for i in range(self.n_bins)]
+            min_xs, min_ys = zip(*[(x, sig/math.sqrt(sig+bg)) for x, sig, bg in zip(xs, min_sig, min_bg)
+                                   if (sig+bg) > 0])
+            bottom_rhs.plot(min_xs, min_ys, '->', color=rhs_color)
+
+        if low_cut_significance:
+            # s/(s+b) for events passing a maximum cut requirement
+            max_bg = [sum(bgs[:i]) for i in range(self.n_bins)]
+            max_sig = [sum(sigs[:i]) for i in range(self.n_bins)]
+            max_xs, max_ys = zip(*[(x, sig/math.sqrt(sig+bg)) for x, sig, bg in zip(xs, max_sig, max_bg)
+                                   if (sig+bg) > 0])
+            bottom_rhs.plot(max_xs, max_ys, '-<', color=rhs_color)
+
+        bottom.set_ylabel(r'$S/(S+B)$')
+        bottom.set_xlim(axes.get_xlim())
+        bottom.set_ylim((0, 1.1))
+        if low_cut_significance or high_cut_significance:
+            bottom_rhs.set_ylabel(r'$S/\sqrt{S+B}$')
+            bottom_rhs.yaxis.label.set_color(rhs_color)
+            bottom_rhs.tick_params(axis='y', colors=rhs_color, size=4, width=1.5)
+        # bottom.grid()
+
+
+if __name__ == '__main__':
+    import matplotlib.pyplot as plt
+    from utils import ResultSet
+
+    rs_TTZ =  ResultSet("TTZ",  "../data/TTZToLLNuNu_treeProducerSusyMultilepton_tree.root")
+    rs_TTW  = ResultSet("TTW",  "../data/TTWToLNu_treeProducerSusyMultilepton_tree.root")
+    rs_TTH  = ResultSet("TTH", "../data/TTHnobb_mWCutfix_ext1_treeProducerSusyMultilepton_tree.root")
+    rs_TTTT = ResultSet("TTTT", "../data/TTTT_ext_treeProducerSusyMultilepton_tree.root")
+
+    sh = StackHist('B-Jet Multiplicity')
+    sh.add_mc_background(rs_TTZ.b_jet_count, 'TTZ', lumi=40)
+    sh.add_mc_background(rs_TTW.b_jet_count, 'TTW', lumi=40)
+    sh.add_mc_background(rs_TTH.b_jet_count, 'TTH', lumi=40)
+    sh.set_mc_signal(rs_TTTT.b_jet_count, 'TTTT', lumi=40, scale=10)
+
+    sh.luminosity = 40
+    sh.energy = 13
+    sh.xlabel = 'B-Jet Count'
+    sh.ylabel = r'\# Events'
+    sh.xlim = (-.5, 9.5)
+    sh.signal_stack = False
+
+    fig = plt.figure()
+    sh.draw(fig.gca())
+    plt.show()
+    # sh.add_data(rs_TTZ.b_jet_count, 'TTZ')

utils.py

@@ -0,0 +1,317 @@
+
+import io
+import sys
+import itertools as it
+from os.path import dirname, join, abspath, normpath
+from math import ceil, floor, sqrt
+from collections import deque
+from IPython.display import Image
+
+import ROOT
+from graph_vals import parse
+
+PRJ_PATH = normpath(join(dirname(abspath(__file__)), "../"))
+EXE_PATH = join(PRJ_PATH, "build/main")
+
+PDG = {1:   'd',   -1:  'd̄',
+       2:   'u',   -2:  'ū',
+       3:   's',   -3:  's̄',
+       4:   'c',   -4:  'c̄',
+       5:   'b',   -5:  'b̄',
+       6:   't',   -6:  't̄',
+
+       11:  'e-',  -11: 'e+',
+       12:  'ν_e', -12: 'ῡ_e',
+
+       13:  'μ-',  -13: 'μ+',
+       14:  'ν_μ', -14: 'ῡ_μ',
+
+       15:  'τ-',  -15: 'τ+',
+       16:  'ν_τ', -16: 'ῡ_τ',
+
+       21:  'g',
+       22:  'γ',
+       23:  'Z0',
+       24:  'W+',  -24: 'W-',
+       25:  'H',
+       }
+
+SINGLE_PLOT_SIZE = (600, 450)
+MAX_WIDTH = 1800
+
+SCALE = .75
+CAN_SIZE_DEF = (int(1600*SCALE), int(1200*SCALE))
+CANVAS = ROOT.TCanvas("c1", "", *CAN_SIZE_DEF)
+ROOT.gStyle.SetPalette(112)  # set the "virdidis" color map
+
+VALUES = {}
+
+
+def clear():
+    CANVAS.Clear()
+    CANVAS.SetCanvasSize(*CAN_SIZE_DEF)
+
+
+def get_color(val, max_val, min_val=0):
+    val = (val-min_val)/(max_val-min_val)
+    val = round(val * (ROOT.gStyle.GetNumberOfColors()-1))
+    col_idx = ROOT.gStyle.GetColorPalette(val)
+    col = ROOT.gROOT.GetColor(col_idx)
+    r = floor(256*col.GetRed())
+    g = floor(256*col.GetGreen())
+    b = floor(256*col.GetBlue())
+    gs = (r + g + b)//3
+    text_color = 'white' if gs < 100 else 'black'
+    return '#{:02x}{:02x}{:02x}'.format(r, g, b), text_color
+
+
+def show_function(dataset, fname):
+    from IPython.display import Markdown
+
+    def md_single(fname_):
+        impl = dataset._function_impl_lookup[fname_]
+        return '*{}*\n-----\n```cpp\n{}\n```\n\n---'.format(fname_, impl)
+    try:
+        return Markdown('\n'.join(md_single(fname_) for fname_ in iter(fname)))
+    except TypeError:
+        return Markdown(md_single(fname))
+
+
+def show_value(dataset, container):
+    if type(container) != str:
+        container = container.GetName().split(':')[1]
+    g, functions = parse(VALUES[container], container)
+    try:
+        return Image(g.create_gif()), show_function(dataset, functions)
+    except Exception as e:
+        print(e)
+        print(g.to_string())
+
+
+class OutputCapture:
+    def __init__(self):
+        self.my_stdout = io.StringIO()
+        self.my_stderr = io.StringIO()
+
+    def get_stdout(self):
+        self.my_stdout.seek(0)
+        return self.my_stdout.read()
+
+    def get_stderr(self):
+        self.my_stderr.seek(0)
+        return self.my_stderr.read()
+
+    def __enter__(self):
+        self.stdout = sys.stdout
+        self.stderr = sys.stderr
+        sys.stdout = self.my_stdout
+        sys.stderr = self.my_stderr
+
+    def __exit__(self, *args):
+        sys.stdout = self.stdout
+        sys.stderr = self.stderr
+        self.stdout = None
+        self.stderr = None
+
+
+def normalize_columns(hist2d):
+    normHist = ROOT.TH2D(hist2d)
+    cols, rows = hist2d.GetNbinsX(), hist2d.GetNbinsY()
+    for col in range(1, cols+1):
+        sum_ = 0
+        for row in range(1, rows+1):
+            sum_ += hist2d.GetBinContent(col, row)
+        if sum_ == 0:
+            continue
+        for row in range(1, rows+1):
+            norm = hist2d.GetBinContent(col, row) / sum_
+            normHist.SetBinContent(col, row, norm)
+    return normHist
+
+
+class ResultSet:
+
+    def __init__(self, sample_name, input_filename):
+        self.sample_name = sample_name
+        self.input_filename = input_filename
+        # self.output_filename = self.input_filename.replace(".root", "_result.root")
+        # self.conditional_recompute()
+        self.load_objects()
+
+        ResultSet.add_collection(self)
+
+    def load_objects(self):
+        file = ROOT.TFile.Open(self.input_filename)
+        l = file.GetListOfKeys()
+        self.map = {}
+        VALUES.update(dict(file.Get("_value_lookup")))
+        for i in range(l.GetSize()):
+            name = l.At(i).GetName()
+            new_name = ":".join((self.sample_name, name))
+            obj = file.Get(name)
+            try:
+                obj.SetName(new_name)
+                obj.SetDirectory(0)  # disconnects Object from file
+            except AttributeError:
+                pass
+            if 'ROOT.vector<int>' in str(type(obj)) and '_count' in name:
+                obj = obj[0]
+            self.map[name] = obj
+            setattr(self, name, obj)
+        file.Close()
+
+        # Now add these histograms into the current ROOT directory (in memory)
+        # and remove old versions if needed
+        for obj in self.map.values():
+            try:
+                old_obj = ROOT.gDirectory.Get(obj.GetName())
+                ROOT.gDirectory.Remove(old_obj)
+                ROOT.gDirectory.Add(obj)
+            except AttributeError:
+                pass
+
+    @classmethod
+    def calc_shape(cls, n_plots):
+        if n_plots*SINGLE_PLOT_SIZE[0] > MAX_WIDTH:
+            shape_x = MAX_WIDTH//SINGLE_PLOT_SIZE[0]
+            shape_y = ceil(n_plots / shape_x)
+            return (shape_x, shape_y)
+        else:
+            return (n_plots, 1)
+
+    def draw(self, shape=None):
+        objs = [obj for obj in self.map.values() if hasattr(obj, "Draw")]
+        if shape is None:
+            n_plots = len(objs)
+            shape = self.calc_shape(n_plots)
+        CANVAS.Clear()
+        CANVAS.SetCanvasSize(shape[0]*SINGLE_PLOT_SIZE[0], shape[1]*SINGLE_PLOT_SIZE[1])
+        CANVAS.Divide(*shape)
+        i = 1
+        for hist in objs:
+            CANVAS.cd(i)
+            try:
+                hist.SetStats(False)
+            except AttributeError:
+                pass
+            if type(hist) in (ROOT.TH1I, ROOT.TH1F, ROOT.TH1D):
+                hist.SetMinimum(0)
+            hist.Draw(self.get_draw_option(hist))
+            i += 1
+        CANVAS.Draw()
+
+    @staticmethod
+    def get_draw_option(obj):
+        obj_type = type(obj)
+        if obj_type in (ROOT.TH1F, ROOT.TH1I, ROOT.TH1D):
+            return ""
+        elif obj_type in (ROOT.TH2F, ROOT.TH2I, ROOT.TH2D):
+            return "COLZ"
+        elif obj_type in (ROOT.TGraph,):
+            return "A*"
+        else:
+            return None
+
+    @classmethod
+    def get_hist_set(cls, attrname):
+        labels, hists = zip(*[(sample_name, getattr(h, attrname))
+                              for sample_name, h in cls.collections.items()])
+        return labels, hists
+
+    @classmethod
+    def add_collection(cls, hc):
+        if not hasattr(cls, "collections"):
+            cls.collections = {}
+        cls.collections[hc.sample_name] = hc
+
+    @classmethod
+    def stack_hist(cls,
+                   hist_name,
+                   title="",
+                   enable_fill=False,
+                   normalize_to=0,
+                   draw=False,
+                   draw_canvas=True,
+                   draw_option="",
+                   make_legend=False,
+                   _stacks={}):
+        labels, hists = cls.get_hist_set(hist_name)
+        if draw_canvas:
+            CANVAS.Clear()
+            CANVAS.SetCanvasSize(SINGLE_PLOT_SIZE[0],
+                                 SINGLE_PLOT_SIZE[1])
+
+        colors = it.cycle([ROOT.kRed, ROOT.kBlue, ROOT.kGreen, ROOT.kYellow])
+        stack = ROOT.THStack(hist_name+"_stack", title)
+        if labels is None:
+            labels = [hist.GetName() for hist in hists]
+        if type(normalize_to) in (int, float):
+            normalize_to = [normalize_to]*len(hists)
+        ens = enumerate(zip(hists, labels, colors, normalize_to))
+        for i, (hist, label, color, norm) in ens:
+            hist_copy = hist
+            hist_copy = hist.Clone(hist.GetName()+"_clone" + draw_option)
+            hist_copy.SetTitle(label)
+            if enable_fill:
+                hist_copy.SetFillColorAlpha(color, 0.75)
+                hist_copy.SetLineColorAlpha(color, 0.75)
+            if norm:
+                integral = hist_copy.Integral()
+                hist_copy.Scale(norm/integral, "nosw2")
+                hist_copy.SetStats(True)
+            stack.Add(hist_copy)
+        if draw:
+            stack.Draw(draw_option)
+            if make_legend:
+                CANVAS.BuildLegend(0.75, 0.75, 0.95, 0.95, "")
+        # prevent stack from getting garbage collected
+        _stacks[stack.GetName()] = stack
+        if draw_canvas:
+            CANVAS.Draw()
+        return stack
+
+    @classmethod
+    def stack_hist_array(cls,
+                         hist_names,
+                         titles,
+                         shape=None, **kwargs):
+        n_hist = len(hist_names)
+        if shape is None:
+            if n_hist <= 4:
+                shape = (1, n_hist)
+            else:
+                shape = (ceil(sqrt(n_hist)),)*2
+        CANVAS.SetCanvasSize(SINGLE_PLOT_SIZE[0]*shape[0],
+                             SINGLE_PLOT_SIZE[1]*shape[1])
+        CANVAS.Divide(*shape)
+        for i, hist_name, title in zip(range(1, n_hist+1), hist_names, titles):
+            CANVAS.cd(i)
+            cls.stack_hist(hist_name, title=title, draw=True,
+                           draw_canvas=False, **kwargs)
+        CANVAS.cd(n_hist).BuildLegend(0.75, 0.75, 0.95, 0.95, "")
+
+    pts = deque([], 50)
+
+    @classmethod
+    def hist_array_single(cls,
+                          hist_name,
+                          title=None,
+                          **kwargs):
+        n_hist = len(cls.collections)
+        shape = cls.calc_shape(n_hist)
+        CANVAS.SetCanvasSize(SINGLE_PLOT_SIZE[0]*shape[0],
+                             SINGLE_PLOT_SIZE[1]*shape[1])
+        CANVAS.Divide(*shape)
+        labels, hists = cls.get_hist_set(hist_name)
+
+        def pave_loc():
+            hist.Get
+        for i, label, hist in zip(range(1, n_hist+1), labels, hists):
+            CANVAS.cd(i)
+            hist.SetStats(False)
+            hist.Draw(cls.get_draw_option(hist))
+
+            pt = ROOT.TPaveText(0.70, 0.87, 0.85, 0.95, "NDC")
+            pt.AddText("Dataset: "+label)
+            pt.Draw()
+            cls.pts.append(pt)