EGamma_ElectronTrackingValidation/plotting/eff_plots.py

#!/usr/bin/env python
from itertools import product
import numpy as np
import matplotlib.pyplot as plt

from uproot import open as root_open
from matplottery.utils import Hist1D, Hist2D, to_html_table
from matplottery.plotter import set_defaults, plot_2d

import matplotboard as mpb

matching_cuts = {
    'new-extra-narrow': [
        dict(
            dPhiMaxHighEt=0.025,
            dPhiMaxHighEtThres=20.0,
            dPhiMaxLowEtGrad=-0.002,
            dRzMaxHighEt=9999.0,
            dRzMaxHighEtThres=0.0,
            dRzMaxLowEtGrad=0.0,
        ),
        dict(
            dPhiMaxHighEt=0.0015,
            dPhiMaxHighEtThres=0.0,
            dPhiMaxLowEtGrad=0.0,
            dRzMaxHighEt=0.025,
            dRzMaxHighEtThres=30.0,
            dRzMaxLowEtGrad=-0.002,
        ),
        dict(
            dPhiMaxHighEt=0.0015,
            dPhiMaxHighEtThres=0.0,
            dPhiMaxLowEtGrad=0.0,
            dRzMaxHighEt=0.025,
            dRzMaxHighEtThres=30.0,
            dRzMaxLowEtGrad=-0.002,
        )
    ],
    'new-default': [
        dict(
            dPhiMaxHighEt=0.05,
            dPhiMaxHighEtThres=20.0,
            dPhiMaxLowEtGrad=-0.002,
            dRzMaxHighEt=9999.0,
            dRzMaxHighEtThres=0.0,
            dRzMaxLowEtGrad=0.0,
        ),
        dict(
            dPhiMaxHighEt=0.003,
            dPhiMaxHighEtThres=0.0,
            dPhiMaxLowEtGrad=0.0,
            dRzMaxHighEt=0.05,
            dRzMaxHighEtThres=30.0,
            dRzMaxLowEtGrad=-0.002,
        ),
        dict(
            dPhiMaxHighEt=0.003,
            dPhiMaxHighEtThres=0.0,
            dPhiMaxLowEtGrad=0.0,
            dRzMaxHighEt=0.05,
            dRzMaxHighEtThres=30.0,
            dRzMaxLowEtGrad=-0.002,
        )
    ],
    'new-wide': [
        dict(
            dPhiMaxHighEt=0.10,
            dPhiMaxHighEtThres=20.0,
            dPhiMaxLowEtGrad=-0.002,
            dRzMaxHighEt=9999.0,
            dRzMaxHighEtThres=0.0,
            dRzMaxLowEtGrad=0.0,
        ),
        dict(
            dPhiMaxHighEt=0.006,
            dPhiMaxHighEtThres=0.0,
            dPhiMaxLowEtGrad=0.0,
            dRzMaxHighEt=0.10,
            dRzMaxHighEtThres=30.0,
            dRzMaxLowEtGrad=-0.002,
        ),
        dict(
            dPhiMaxHighEt=0.006,
            dPhiMaxHighEtThres=0.0,
            dPhiMaxLowEtGrad=0.0,
            dRzMaxHighEt=0.10,
            dRzMaxHighEtThres=30.0,
            dRzMaxLowEtGrad=-0.002,
        )
    ],
    'new-extra-wide': [
        dict(
            dPhiMaxHighEt=0.15,
            dPhiMaxHighEtThres=20.0,
            dPhiMaxLowEtGrad=-0.002,
            dRzMaxHighEt=9999.0,
            dRzMaxHighEtThres=0.0,
            dRzMaxLowEtGrad=0.0,
        ),
        dict(
            dPhiMaxHighEt=0.009,
            dPhiMaxHighEtThres=0.0,
            dPhiMaxLowEtGrad=0.0,
            dRzMaxHighEt=0.15,
            dRzMaxHighEtThres=30.0,
            dRzMaxLowEtGrad=-0.002,
        ),
        dict(
            dPhiMaxHighEt=0.009,
            dPhiMaxHighEtThres=0.0,
            dPhiMaxLowEtGrad=0.0,
            dRzMaxHighEt=0.15,
            dRzMaxHighEtThres=30.0,
            dRzMaxLowEtGrad=-0.002,
        )
    ],
}

samples = None


def load_samples():
    global samples
    if samples is None:
        print("loading samples")
        samples = {(proc, wp): root_open(f'../hists/{proc}-{wp}.root') for proc, wp in product(procs, wps)}
    return samples


def calc_window(et, eta, hit, variable, cut_sel):
    idx = min(hit-1, 2)
    cuts = matching_cuts[cut_sel][idx]
    if 'etaBins' in cuts:
        for eta_idx, bin_high in enumerate(cuts['etaBins']):
            if eta < bin_high:
                high_et = cuts[f'{variable}MaxHighEt'][eta_idx]
                high_et_thres = cuts[f'{variable}MaxHighEtThres'][eta_idx]
                low_et_grad = cuts[f'{variable}MaxLowEtGrad'][eta_idx]
                break
        else:  # highest bin
            high_et = cuts[f'{variable}MaxHighEt'][-1]
            high_et_thres = cuts[f'{variable}MaxHighEtThres'][-1]
            low_et_grad = cuts[f'{variable}MaxLowEtGrad'][-1]
    else:
        high_et = cuts[f'{variable}MaxHighEt']
        high_et_thres = cuts[f'{variable}MaxHighEtThres']
        low_et_grad = cuts[f'{variable}MaxLowEtGrad']
    return high_et + min(0, et-high_et_thres)*low_et_grad


def hist_integral_ratio(num, den):
    num_int = num.get_integral()
    den_int = den.get_integral()

    ratio = num_int / den_int
    error = np.sqrt(den_int) / den_int  # TODO: Check this definition of error
    return ratio, error


def center_text(x, y, txt, **kwargs):
    plt.text(x, y, txt,
             horizontalalignment='center', verticalalignment='center',
             transform=plt.gca().transAxes, size=24, **kwargs)


def hist_plot(h: Hist1D, *args, include_errors=False, line_width=1, **kwargs):
    """ Plots a 1D ROOT histogram object using matplotlib """

    counts = h.get_counts()
    edges = h.get_edges()
    left, right = edges[:-1], edges[1:]
    x = np.array([left, right]).T.flatten()
    y = np.array([counts, counts]).T.flatten()

    plt.plot(x, y, *args, linewidth=line_width, **kwargs)
    if include_errors:
        plt.errorbar(h.get_bin_centers(), h.counts, yerr=h.errors,
                     color='k', marker=None, linestyle='None',
                     barsabove=True, elinewidth=.7, capsize=1)


def hist2d_percent_contour(h: Hist1D, percent: float, axis: str):
    values = h.counts
    try:
        axis = axis.lower()
        axis_idx = {'x': 1, 'y': 0}[axis]
    except KeyError:
        raise ValueError('axis must be \'x\' or \'y\'')
    if percent < 0 or percent > 1:
        raise ValueError('percent must be in [0,1]')

    with np.warnings.catch_warnings():
        np.warnings.filterwarnings('ignore', 'invalid value encountered in true_divide')
        values = values / np.sum(values, axis=axis_idx, keepdims=True)
        np.nan_to_num(values, copy=False)
    values = np.cumsum(values, axis=axis_idx)
    idxs = np.argmax(values > percent, axis=axis_idx)

    x_centers, y_centers = h.get_bin_centers()

    if axis == 'x':
        return x_centers[idxs], y_centers
    else:
        return x_centers, y_centers[idxs]


@mpb.decl_fig
def plot_residuals(sample, layer, hit, variable, subdet):
    load_samples()
    proc, wp = sample

    track_matched = samples[sample][f'{variable}_{subdet}_L{layer}_H{hit}_v_Et_TrackMatched']
    # no_match = samples[sample][f'{variable}_{subdet}_L{layer}_H{hit}_v_Et_NoMatch']

    h_real = Hist2D(track_matched)
    # h_fake = Hist2D(no_match)

    def do_plot(h):
        plot_2d(h, cmap='viridis')

        xs, ys = hist2d_percent_contour(h, .90, 'x')
        plt.plot(xs, ys, color='green', label='90\% contour')
        xs, ys = hist2d_percent_contour(h, .995, 'x')
        plt.plot(xs, ys, color='darkgreen', label='99.5\% contour')

        ets = h.edges[1]
        cuts = [calc_window(et, 0, hit, variable, wp) for et in ets]
        plt.plot(cuts, ets, color='red', label='Cut Value')
        plt.xlabel({'dPhi': r'$\delta \phi$ (rads)',
                    'dRz': r'$\delta R/z$ (cm)'}[variable])

    # plt.sca(plt.subplot(1, 2, 1))
    do_plot(h_real)
    plt.title('Truth-Matched Seeds')
    plt.ylabel('$E_T$ (GeV)')

    # plt.sca(plt.subplot(1, 2, 2))
    # do_plot(h_fake)
    # plt.title('Not Truth-Matched Seeds')
    plt.legend(loc='upper right')


@mpb.decl_fig
def plot_residuals_eta(sample, hit, variable):
    load_samples()

    h = Hist2D(samples[sample][f'{variable}_residuals_v_eta_H{hit}'])

    plot_2d(h)

    xs, ys = hist2d_percent_contour(h, .90, 'x')
    plt.plot(xs, ys, color='green', label='90\% contour')
    xs, ys = hist2d_percent_contour(h, .995, 'x')
    plt.plot(xs, ys, color='darkgreen', label='99.5\% contour')

    plt.xlabel({'dPhi': r'$\delta \phi$ (rads)',
                'dRz': r'$\delta R/z$ (cm)'}[variable])
    plt.ylabel(r'$|\eta|$')


@mpb.decl_fig
def plot_hit_vs_layer(sample, region):
    load_samples()
    h = Hist2D(samples[sample][f'hit_vs_layer_{region}'])
    plot_2d(h, colz_fmt='2.0f')
    plt.xlabel('Layer #')
    plt.ylabel('Hit #')


@mpb.decl_fig
def plot_roc_curve(pfx, ext=''):
    load_samples()

    def get_num_den(sample, basename):
        num = Hist1D(sample[f'{basename}_num'])
        den = Hist1D(sample[f'{basename}_den'])
        return hist_integral_ratio(num, den)
    rows = []
    for (proc, wp), sample in samples.items():
        sample_name = f'{proc}-{wp}'
        eff, eff_err = get_num_den(sample, f'{pfx}_eff_v_phi{ext}')
        pur, pur_err = get_num_den(sample, f'{pfx}_pur_v_phi{ext}')
        rows.append([wp,
                     rf'${eff*100:0.2f}\pm{eff_err*100:0.2f}\%$',
                     rf'${pur*100:0.2f}\pm{pur_err*100:0.2f}\%$'])

        plt.errorbar([pur], [eff], xerr=[pur_err], yerr=[eff_err],
                     label=sample_name, marker='o', color=color(proc, wp))
    center_text(0.3, 0.3, r'$p_T>20$ and $|\eta|<2.4$')
    plt.axis('equal')
    plt.xlim((0.5, 1.02))
    plt.ylim((0.5, 1.02))
    plt.xlabel('Purity')
    plt.ylabel('Efficiency')
    plt.grid()
    plt.legend(loc='lower right')

    col_labels = ['Sample', 'Working Point', 'Efficiency', 'Purity']
    row_labels = [r'$Z \rightarrow ee$', '', '', r'$t\bar{t}$', '', '']
    return to_html_table(rows, col_labels, row_labels, 'table-condensed')


@mpb.decl_fig
def plot_kinematic_eff(pref, ext='', xlim=(None, None), ylim=(None, None), norm=None, label_pfx='', incl_sel=True):
    load_samples()
    ax_pt = plt.subplot(221)
    ax_eta = plt.subplot(222)
    ax_phi = plt.subplot(223)
    errors = True
    for (proc, wp), sample in samples.items():
        sample_name = f'{proc}-{wp}'
        l = sample_name
        c = color(proc, wp)

        def do_plot(ax, name):
            plt.sca(ax)
            h = Hist1D(sample[name], no_overflow=True)
            if norm:
                h = h / (norm*h.get_integral())
            hist_plot(h, include_errors=errors, label=l, color=c)

        do_plot(ax_pt, f'{pref}_v_pt{ext}')
        do_plot(ax_eta, f'{pref}_v_eta{ext}')
        do_plot(ax_phi, f'{pref}_v_phi{ext}')

    plt.sca(ax_pt)
    if not incl_sel: center_text(0.5, 0.15, r'$|\eta|<2.4$')
    plt.xlabel(fr"{label_pfx} $p_T$")
    plt.ylim(ylim)
    plt.xlim(xlim)

    plt.sca(ax_eta)
    if not incl_sel: center_text(0.5, 0.15, r'$p_T>20$')
    plt.xlabel(fr"{label_pfx} $\eta$")
    plt.ylim(ylim)
    plt.xlim(xlim)

    plt.sca(ax_phi)
    if not incl_sel: center_text(0.5, 0.15, r'$p_T>20$ and $|\eta|<2.4$')
    plt.xlabel(fr"{label_pfx} $\phi$")
    plt.ylim(ylim)
    plt.xlim(xlim)
    plt.tight_layout()
    plt.legend(loc='upper left', bbox_to_anchor=(0.6, 0.45), bbox_transform=plt.gcf().transFigure)


@mpb.decl_fig
def plot_ecal_rel_res():
    load_samples()
    for sample_name, sample in samples.items():
        h = Hist1D(sample['ecal_energy_resolution'])
        h = h / h.get_integral()
        hist_plot(h, label=sample_name)
    plt.xlabel(r"ECAL $E_T$ relative error")
    plt.legend()


@mpb.decl_fig
def plot_res_contour(proc, hit_number, var, layers, ext='_TrackMatched'):
    load_samples()
    _, axs = plt.subplots(1, 2, sharey=True)

    def do_plot(ax, sample):
        plt.sca(ax)
        wp = sample[1]
        plt.title(wp)
        h = None
        for subdet, layer in layers:
            h = Hist2D(samples[sample][f'{var}_{subdet}_L{layer}_H{hit_number}_v_Et{ext}'])
            xs, ys = hist2d_percent_contour(h, .99, 'x')
            plt.plot(xs, ys, label=f'{subdet} - L{layer}')

        ets = h.edges[1]
        cuts = [calc_window(et, 0, hit_number, var, wp) for et in ets]
        plt.plot(cuts, ets, color='k', label='Cut Value')

    for ax, wp in zip(axs, ['new-default', 'new-wide']):
        do_plot(ax, (proc, wp))

    plt.sca(axs[-1])
    plt.legend(loc='upper right')


@mpb.decl_fig
def simple_dist(hist_name, rebin=(), norm=1, xlabel="", ylabel="", xlim=None, ylim=None, line_width=1):
    load_samples()
    for (proc, wp), sample in samples.items():
        sample_name = f'{proc}-{wp}'
        h = Hist1D(sample[hist_name])
        if rebin:
            h.rebin(*rebin)

        mean = np.sum(h.get_counts() * h.get_bin_centers()) / h.get_integral()
        if norm is not None:
            h = h * (norm / h.get_integral())
        hist_plot(h, label=f'{sample_name} ($\\mu={mean:.2f}$)',
                  color=color(proc, wp), line_width=line_width)
    if xlim:
        plt.xlim(xlim)
    if ylim:
        plt.ylim(ylim)
    plt.xlabel(xlabel)
    plt.ylabel(ylabel)
    plt.legend()


@mpb.decl_fig
def simple_dist2d(hist_name, proc, wp, xlabel="", ylabel="", xlim=None, ylim=None, norm=None):
    load_samples()
    sample = samples[(proc, wp)]
    # sample_name = f'{proc}-{wp}'
    h = Hist2D(sample[hist_name])
    if norm is not None:
        h = h * (norm / h.get_integral())
    plot_2d(h, colz_fmt='g')
    if xlim:
        plt.xlim(xlim)
    if ylim:
        plt.ylim(ylim)
    plt.xlabel(xlabel)
    plt.ylabel(ylabel)


def all_cut_plots(refresh=True, publish=False):

    figures = {
        'tracking_roc_curve': (plot_roc_curve, ('tracking',)),
        'tracking_roc_curve_dR': (plot_roc_curve, ('tracking',), {'ext': '_dR'}),
        'seeding_roc_curve': (plot_roc_curve, ('seed',)),

        'number_of_seeds': (simple_dist, ('n_seeds',),
                            dict(xlabel='Number of Seeds', rebin=(50, -0.5, 200.5))),
        'number_of_good_seeds': (simple_dist, ('n_good_seeds',),
                                 dict(xlabel='Number of Good Seeds', rebin=(50, -0.5, 200.5))),
        'number_of_scls': (simple_dist, ('n_scl',),
                           dict(xlabel='Number of Super-Clusters', xlim=(-0.5, 25.5))),
        'number_of_good_scls': (simple_dist, ('n_good_scl',),
                                dict(xlabel='Number of Super-Clusters', xlim=(-0.5, 25.5))),

        'number_of_sim_els': (simple_dist, ('n_good_sim',),
                              dict(xlabel='Number of prompt(ish) electrons', xlim=(-0.5, 20.5))),
        'number_of_gsf_tracks': (simple_dist, ('n_gsf_track',),
                                 dict(xlabel='Number of reco electrons', xlim=(-0.5, 20.5))),

        'number_of_matched': (simple_dist, ('n_matched',),
                              dict(xlabel='Number of matched electrons', xlim=(-0.5, 10.5), line_width=4)),
        'number_of_merged': (simple_dist, ('n_merged',),
                             dict(xlabel='Number of merged electrons', xlim=(-0.5, 10.5), line_width=4)),
        'number_of_lost': (simple_dist, ('n_lost',),
                           dict(xlabel='Number of lost electrons', xlim=(-0.5, 10.5), line_width=4)),
        'number_of_split': (simple_dist, ('n_split',),
                            dict(xlabel='Number of split electrons', xlim=(-0.5, 10.5), line_width=4)),
        'number_of_faked': (simple_dist, ('n_faked',),
                            dict(xlabel='Number of faked electrons', xlim=(-0.5, 10.5), line_width=4)),
        'number_of_flipped': (simple_dist, ('n_flipped',),
                              dict(xlabel='Number of flipped electrons', xlim=(-0.5, 10.5), line_width=4)),
        'matched_dR': (simple_dist, ('matched_dR',),
                       dict(xlabel='dR between sim and reco')),
        'matched_dpT': (simple_dist, ('matched_dpT',),
                        dict(xlabel='dpT between sim and reco')),

        'number_of_matched_dR': (simple_dist, ('n_matched_dR',),
                                 dict(xlabel='Number of matched electrons - dR Matched', xlim=(-0.5, 10.5),
                                      line_width=4)),
        'number_of_merged_dR': (simple_dist, ('n_merged_dR',),
                                dict(xlabel='Number of merged electrons - dR Matched', xlim=(-0.5, 10.5),
                                     line_width=4)),
        'number_of_lost_dR': (simple_dist, ('n_lost_dR',),
                              dict(xlabel='Number of lost electrons - dR Matched', xlim=(-0.5, 10.5),
                                   line_width=4)),
        'number_of_split_dR': (simple_dist, ('n_split_dR',),
                               dict(xlabel='Number of split electrons - dR Matched', xlim=(-0.5, 10.5),
                                    line_width=4)),
        'number_of_faked_dR': (simple_dist, ('n_faked_dR',),
                               dict(xlabel='Number of faked electrons - dR Matched', xlim=(-0.5, 10.5),
                                    line_width=4)),
        'number_of_flipped_dR': (simple_dist, ('n_flipped_dR',),
                                 dict(xlabel='Number of flipped electrons - dR Matched', xlim=(-0.5, 10.5),
                                      line_width=4)),
        'matched_dR_dR': (simple_dist, ('matched_dR_dR',),
                          dict(xlabel='dR between sim and reco - dR Matched')),
        'matched_dpT_dR': (simple_dist, ('matched_dpT_dR',),
                           dict(xlabel='dpT between sim and reco - dR Matched')),

        'sim_pt': (simple_dist, ('sim_pt',),
                   dict(xlabel='Sim Track $p_T$', xlim=(0, None))),
        'sim_eta': (simple_dist, ('sim_eta',),
                    dict(xlabel='Sim Track $eta$', rebin=(20, -3, 3))),
        'sim_phi': (simple_dist, ('sim_phi',),
                    dict(xlabel='Sim Track $phi$', rebin=(20, -3.14, 3.14), ylim=(0, None))),
        'reco_pt': (simple_dist, ('reco_pt',),
                    dict(xlabel='Reco Track $p_T$', xlim=(0, None))),
        'reco_eta': (simple_dist, ('reco_eta',),
                     dict(xlabel='Reco Track $eta$', rebin=(20, -3, 3))),
        'reco_phi': (simple_dist, ('reco_phi',),
                     dict(xlabel='Reco Track $phi$', rebin=(20, -3.14, 3.14), ylim=(0, None))),

        'tm_corr': (simple_dist2d, ('tm_corr', 'zee', 'old-default'),
                    dict(xlabel='Seed Matched', ylabel='Track Matched', norm=1)),

        'ecal_rel_res': plot_ecal_rel_res,
        'hit_v_layer_BPIX_new-default_zee': (plot_hit_vs_layer, (('zee', 'new-default'), 'barrel')),
        'hit_v_layer_FPIX_new-default_zee': (plot_hit_vs_layer, (('zee', 'new-default'), 'forward')),
        'hit_v_layer_BPIX_new-default_tt': (plot_hit_vs_layer, (('tt', 'new-default'), 'barrel')),
        'hit_v_layer_FPIX_new-default_tt': (plot_hit_vs_layer, (('tt', 'new-default'), 'forward')),
        'hit_v_layer_BPIX_new-wide_zee': (plot_hit_vs_layer, (('zee', 'new-wide'), 'barrel')),
        'hit_v_layer_FPIX_new-wide_zee': (plot_hit_vs_layer, (('zee', 'new-wide'), 'forward')),
        'hit_v_layer_BPIX_new-wide_tt': (plot_hit_vs_layer, (('tt', 'new-wide'), 'barrel')),
        'hit_v_layer_FPIX_new-wide_tt': (plot_hit_vs_layer, (('tt', 'new-wide'), 'forward')),


        'good_sim_kinem': (plot_kinematic_eff, ('good_sim',),
                           dict(norm=1, ylim=(0, None))),
        'gsf_track_kinem': (plot_kinematic_eff, ('gsf_track',),
                            dict(norm=1, ylim=(0, None))),
        'seed_kinem': (plot_kinematic_eff, ('seed',),
                       dict(norm=1, ylim=(0, None))),
        'scl_kinem': (plot_kinematic_eff, ('scl',),
                      dict(norm=1, ylim=(0, None))),
    }

    def add_num_den(key, func, args, kwargs):
        figures[key] = (func, args, dict(**kwargs, ylim=(0, 1.1)))
        base_ext = kwargs.get('ext', '')
        kwargs_ = kwargs.copy()
        kwargs_['ext'] = base_ext+'_num'
        figures[key+'_num'] = (func, args, kwargs_)
        kwargs_ = kwargs.copy()
        kwargs_['ext'] = base_ext+'_den'
        figures[key+'_den'] = (func, args, kwargs_)

    add_num_den('tracking_eff', plot_kinematic_eff, ('tracking_eff',), dict(incl_sel=False))
    add_num_den('tracking_pur', plot_kinematic_eff, ('tracking_pur',), dict(incl_sel=False))
    add_num_den('tracking_eff_dR', plot_kinematic_eff, ('tracking_eff',), dict(ext='_dR', incl_sel=False))
    add_num_den('tracking_pur_dR', plot_kinematic_eff, ('tracking_pur',), dict(ext='_dR', incl_sel=False))
    add_num_den('seeding_eff', plot_kinematic_eff, ('seed_eff',), dict(incl_sel=False))
    add_num_den('seeding_pur', plot_kinematic_eff, ('seed_pur',), dict(incl_sel=False))
    add_num_den('fake_rate_incl', plot_kinematic_eff, ('fake_rate_incl',), {})
    add_num_den('partial_fake_rate_incl', plot_kinematic_eff, ('partial_fake_rate_incl',), {})
    add_num_den('full_fake_rate_incl', plot_kinematic_eff, ('full_fake_rate_incl',), {})
    add_num_den('clean_fake_rate_incl', plot_kinematic_eff, ('clean_fake_rate_incl',), {})
    add_num_den('fake_rate', plot_kinematic_eff, ('fake_rate',), dict(incl_sel=False))
    add_num_den('partial_fake_rate', plot_kinematic_eff, ('partial_fake_rate',), dict(incl_sel=False))
    add_num_den('full_fake_rate', plot_kinematic_eff, ('full_fake_rate',), dict(incl_sel=False))
    add_num_den('clean_fake_rate', plot_kinematic_eff, ('clean_fake_rate',), dict(incl_sel=False))

    hit_layers = [(1, 1), (1, 2), (2, 2), (2, 3), (3, 3), (3, 4)]
    for proc, wp, (hit, layer), var, subdet in product(['zee', 'tt'], ['new-default', 'new-wide'],
                                                       hit_layers,
                                                       ['dPhi', 'dRz'], ['BPIX', 'FPIX']):
        figures.update({
            f'res_{subdet}_L{layer}_H{hit}_{var}_{proc}_{wp}':
                (plot_residuals, ((proc, wp), layer, hit, var, subdet))})

    rel_layers = {1: [('BPIX', 1), ('BPIX', 2), ('FPIX', 1), ('FPIX', 2)],
                  2: [('BPIX', 2), ('BPIX', 3), ('FPIX', 2), ('FPIX', 3)],
                  3: [('BPIX', 3), ('BPIX', 4), ('FPIX', 3)], }
    for proc, hit, var in product(['zee', 'tt'], [1, 2, 3], ['dPhi', 'dRz']):
        figures.update({
            f'resall_H{hit}_{var}_{proc}': (plot_res_contour, (proc, hit, var, rel_layers[hit]))})

    for proc, wp, hit, var in product(['zee', 'tt'], ['new-default', 'new-wide'], [1, 2, 3], ['dPhi', 'dRz']):
        figures.update({
            f'res_v_eta_H{hit}_{var}_{proc}_{wp}': (plot_residuals_eta, ((proc, wp), hit, var))})

    mpb.render(figures, refresh=refresh)
    mpb.generate_report(figures, 'Electron Seeding Studies',
                        output=f'hists.html',
                        source=__file__)

    mpb.generate_report(figures, 'Update',
                        output='report.html',
                        body='../docs/reports/report_2018_05_30.md')
    if publish:
        mpb.publish()


def color(proc, wp):
    from matplotlib.colors import XKCD_COLORS

    def f(name):
        return XKCD_COLORS['xkcd:'+name]

    return {
        ('zee', 'new-extra-narrow'): f('kelly green'),
        ('tt',  'new-extra-narrow'): f('grey green'),
        ('zee', 'new-default'):      f('red'),
        ('tt',  'new-default'):      f('pink'),
        ('zee', 'new-wide'):         f('strong blue'),
        ('tt',  'new-wide'):         f('bright blue'),
        ('zee', 'new-extra-wide'):   f('indigo'),
        ('tt',  'new-extra-wide'):   f('bright purple'),
        ('zee', 'old-default'):      f('black'),
        ('tt',  'old-default'):      f('blue grey'),
    }[(proc, wp)]


if __name__ == '__main__':
    set_defaults()
    mpb.configure(output_dir='seeding_studies',
                  multiprocess=True,
                  publish_remote="caleb@fangmeier.tech",
                  publish_dir="/var/www/eg",
                  publish_url="eg.fangmeier.tech/",
                  )
    procs = {
        'zee': r'$Z\rightarrow e^+e_-}$',
        'tt': r'$t\bar{t}$'
    }
    wps = {
        'new-default': 'HLT Settings',
        'new-wide': 'Wide Settings',
        'old-default': 'Old Seeding',
    }
    all_cut_plots(refresh=True, publish=True)