EGamma_ElectronTrackingValidation/plotting/eff_plots.py

#!/usr/bin/env python
from itertools import product
from argparse import ArgumentParser
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 = {}
        for proc, wp, config in product(procs, wps, configs):
            samples[(proc, wp, config)] = root_open(f'../hists/{proc}-{wp}-{config}.root')
            # if 'old' in wp:
            #     samples[(proc, wp)] = root_open(f'../hists/with_skipping/{proc}-{wp}.root')
            # elif 'noskip' in wp:
            #     wpkey = wp.replace('-noskip', '')
            #     samples[(proc, wp)] = root_open(f'../hists/no_skipping/{proc}-{wpkey}.root')
            # elif 'skip' in wp:
            #     wpkey = wp.replace('-skip', '')
            #     samples[(proc, wp)] = root_open(f'../hists/with_skipping/{proc}-{wpkey}.root')
    return samples


def calc_window(et, eta, hit, variable, cut_sel):
    idx = min(hit-1, 2)
    cut_sel = cut_sel.replace('-skip', '').replace('-noskip', '')
    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.integral
    den_int = den.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.counts
    edges = h.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:
        if h.errors_up is not None:
            errors = np.vstack((h.errors_down, h.errors_up))
        else:
            errors = h.errors
        plt.errorbar(h.bin_centers, h.counts, yerr=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.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()
    if region == 'BPIX':
        region = 'barrel'
    if region == 'FPIX':
        region = 'forward'
    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()
    show_fr = pfx == "tracking"

    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 = []
    row_labels = []
    for i, proc in enumerate(procs):
        plt.subplot(len(procs), 1, i+1)
        row_labels.append(procs[proc])
        row_labels.extend(['']*(len(wps)*len(configs)-2))
        for wp, config in product(wps, configs):
            sample = samples[(proc, wp, config)]
            sample_name = f'{proc}-{wp}-{config}'
            if wp == 'old-default' and config == 'noskip': continue
            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}')
            if show_fr:
                fr, fr_err = get_num_den(sample, f'fake_rate_no_e_match_v_phi')

                rows.append([wp, config,
                             rf'${eff*100:0.2f}\pm{eff_err*100:0.2f}\%$',
                             rf'${pur*100:0.2f}\pm{pur_err*100:0.2f}\%$',
                             rf'${fr*100:0.2f}\pm{fr_err*100:0.2f}\%$'])
            if config == 'skip-pileup':
                marker = 'o'
            elif config == 'noskip':
                marker = '*'
            else:
                marker = '^'

            # plt.errorbar([pur], [eff], xerr=[pur_err], yerr=[eff_err],
            #              label=sample_name, marker=marker, color=color(proc, wp))
            plt.scatter([pur], [eff], label=sample_name, marker=marker, color=color(proc, wp))

        # center_text(0.3, 0.3, r'$p_T>20$ and $|\eta|<2.5$')
        # plt.axis('equal')
        plt.xlim((0.5, 1.02))
        plt.ylim((0.5, 1.02))
        # plt.xlim((0, 1.02))
        # plt.ylim((0.7, 1.02))
        plt.ylabel('Efficiency')
        plt.grid()
        plt.legend(loc='lower right', ncol=1, fancybox=True, numpoints=1)
    plt.xlabel('Purity')

    col_labels = ['Sample', 'Working Point', 'Config', 'Efficiency', 'Purity']
    if show_fr:
        col_labels.append("Fake Rate")
    return to_html_table(rows, col_labels, row_labels, 'table-condensed')


@mpb.decl_fig
def plot_kinematic_eff(pref, proc, config, ext='', ylim=(None, None), norm=None, label_pfx='', incl_sel=True,
                       bins_pt=None, bins_eta=None, bins_phi=None, bins_PU=None,
                       xlim_pt=(None, None), xlim_eta=(None, None), xlim_phi=(None, None),
                       is_ratio=False):
    load_samples()
    # Figure out if this one has v_PU
    has_PU =  f'{pref}_v_PU{ext}' in list(samples.values())[0]

    ax_pt = plt.subplot(221)
    ax_eta = plt.subplot(222)
    ax_phi = plt.subplot(223)
    if has_PU:
        ax_PU = plt.subplot(224)
    errors = True
    for (proc_, wp, config_), sample in samples.items():
        if proc != proc_ or config != config_: continue
        sample_name = f'{proc}-{wp}-{config}'
        l = sample_name
        c = color(proc, wp)

        def do_plot(ax, name, bins):
            plt.sca(ax)
            if is_ratio:
                num = Hist1D(sample[name+"_num"], no_overflow=True)
                den = Hist1D(sample[name+"_den"], no_overflow=True)
                if bins:
                    num.rebin(bins)
                    den.rebin(bins)
                h = num // den
            else:
                h = Hist1D(sample[name], no_overflow=True)
                if norm:
                    h = h / (norm*h.integral)
                if bins is not None:
                    h.rebin(bins)
            hist_plot(h, include_errors=errors, label=l, color=c, linestyle=style(config))

        do_plot(ax_pt, f'{pref}_v_pt{ext}', bins_pt)
        do_plot(ax_eta, f'{pref}_v_eta{ext}', bins_eta)
        do_plot(ax_phi, f'{pref}_v_phi{ext}', bins_phi)
        if has_PU:
            do_plot(ax_PU, f'{pref}_v_PU{ext}', [-0.5, 20, 30, 40, 50, 60, 70, 80, 100.5])

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

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

    plt.sca(ax_phi)
    if not incl_sel: center_text(0.5, 0.15, r'$p_T>20$ and $|\eta|<2.5$')
    plt.xlabel(fr"{label_pfx} $\phi$")
    plt.ylim(ylim)
    plt.xlim(xlim_phi)

    if has_PU:
        plt.sca(ax_PU)
        if not incl_sel: center_text(0.5, 0.15, r'$p_T>20$ and $|\eta|<2.5$')
        plt.xlabel("# Pile-up Interactions")
        plt.ylim(ylim)
        plt.tight_layout()
        plt.sca(ax_eta)
        plt.legend(loc='best', prop={'size': 10})
    else:
        plt.tight_layout()
        plt.legend(loc='upper left', bbox_to_anchor=(0.6, 0.45), bbox_transform=plt.gcf().transFigure,
                   prop={'size': 15})


@mpb.decl_fig
def plot_kinematic_eff2d(pref, proc, config, ext=''):
    from math import sqrt, ceil
    load_samples()

    n_col = ceil(sqrt(len(wps)))
    n_row = ceil(len(wps) / n_col)

    plt_idx = 0
    for (proc_, wp, config_), sample in samples.items():
        if proc != proc_ or config != config_: continue
        sample_name = f'{proc}-{wp}-{config}'
        plt_idx += 1

        plt.subplot(n_row, n_col, plt_idx)
        h = Hist2D(sample[f'{pref}_v_eta_phi{ext}'], no_overflow=True)
        h.rebin(div_bins_x=4)
        plot_2d(h)
        plt.text(0.1, 0.1, sample_name,
                 transform=plt.gca().transAxes, size=15, backgroundcolor='#FFFFFFA0')
        if (plt_idx-1) % n_col == 0:  # left col of plots
            plt.ylabel(r'$\phi$')
        if plt_idx > n_col*(n_row-1):  # bottom row of plots
            plt.xlabel(r'$\eta$')


@mpb.decl_fig
def plot_ecal_rel_res():
    load_samples()
    for sample_name, sample in samples.items():
        h = Hist1D(sample['ecal_energy_resolution'])
        try:
            h = h / h.integral
            hist_plot(h, label=sample_name)
        except ZeroDivisionError:
            pass
    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-noskip', 'new-wide-noskip']):
        do_plot(ax, (proc, wp))

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


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

        mean = np.sum(h.counts * h.bin_centers) / h.integral
        try:
            if norm is not None:
                h = h * (norm / h.integral)
            hist_plot(h, label=f'{sample_name} ($\\mu={mean:.2f}$)',
                      color=color(proc_, wp_), line_width=line_width, linestyle=style(config_))
        except ZeroDivisionError:
            pass
    if xlim:
        plt.xlim(xlim)
    if ylim:
        plt.ylim(ylim)
    plt.xlabel(xlabel)
    plt.ylabel(ylabel)
    plt.legend(fontsize=20)


@mpb.decl_fig
def simple_dist2d(hist_name, proc, wp, config, norm=None, colz_fmt='g', clear_zero=False, **kwargs):
    load_samples()
    sample = samples[(proc, wp, config)]
    h = Hist2D(sample[hist_name])
    if norm is not None:
        h = h * (norm / h.integral)
    if clear_zero:
        h.counts[0, 0] = 0
    plot_2d(h, colz_fmt=colz_fmt, cmap='viridis', **kwargs)


def all_cut_plots(build=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'),
    }
    for proc, config in product(procs, configs):

        figures[f'number_of_seeds_{proc}_{config}'] = simple_dist('n_seeds', rebin=(50, -0.5, 200.5),
                                                                  xlabel='Number of Seeds', proc=proc, config=config)
        figures[f'number_of_good_seeds_{proc}_{config}'] = simple_dist('n_good_seeds', rebin=(50, -0.5, 200.5),
                                                                       xlabel='Number of Seeds', proc=proc, config=config)
        figures[f'number_of_scls_{proc}_{config}'] = simple_dist('n_scl', xlabel='Number of Super-Clusters',
                                                                 xlim=(-0.5, 25.5), proc=proc, config=config)
        figures[f'number_of_good_scls_{proc}_{config}'] = simple_dist('n_good_scl', xlabel='Number of Super-Clusters',
                                                                      xlim=(-0.5, 25.5), proc=proc, config=config)

        figures[f'number_of_sim_els_{proc}_{config}'] = simple_dist('n_good_sim', xlabel='Number of prompt(ish) electrons',
                                                                    xlim=(-0.5, 20.5), proc=proc, config=config)
        figures[f'number_of_gsf_tracks_{proc}_{config}'] = simple_dist('n_gsf_track', xlabel='Number of reco electrons',
                                                                       xlim=(-0.5, 20.5), proc=proc, config=config)

        figures[f'number_of_prompt_{proc}_{config}'] = simple_dist('n_prompt', xlabel='Number of prompt electrons',
                                                                   xlim=(-0.5, 20.5), proc=proc, config=config)
        figures[f'number_of_nonprompt_{proc}_{config}'] = simple_dist('n_nonprompt', xlabel='Number of nonprompt electrons',
                                                                      xlim=(-0.5, 20.5), proc=proc, config=config)

        figures[f'number_of_matched_{proc}_{config}'] = simple_dist('n_matched_dR', xlabel='Number of matched electrons',
                                                                    xlim=(-0.5, 10.5), line_width=4, proc=proc, config=config)
        figures[f'number_of_merged_{proc}_{config}'] = simple_dist('n_merged_dR', xlabel='Number of merged electrons',
                                                                   xlim=(-0.5, 10.5), line_width=4, proc=proc, config=config)
        figures[f'number_of_lost_{proc}_{config}'] = simple_dist('n_lost_dR', xlabel='Number of lost electrons',
                                                                 xlim=(-0.5, 10.5), line_width=4, proc=proc, config=config)
        figures[f'number_of_split_{proc}_{config}'] = simple_dist('n_split_dR', xlabel='Number of split electrons',
                                                                  xlim=(-0.5, 10.5), line_width=4, proc=proc, config=config)
        figures[f'number_of_faked_{proc}_{config}'] = simple_dist('n_faked_dR', xlabel='Number of faked electrons',
                                                                  xlim=(-0.5, 10.5), line_width=4, proc=proc, config=config)
        figures[f'number_of_flipped_{proc}_{config}'] = simple_dist('n_flipped_dR', xlabel='Number of flipped electrons',
                                                                    xlim=(-0.5, 10.5), line_width=4, proc=proc, config=config)
        figures[f'matched_dR_{proc}_{config}'] = simple_dist('matched_dR', xlabel='dR between sim and reco',
                                                             proc=proc, config=config)
        figures[f'matched_dpT_{proc}_{config}'] = simple_dist('matched_dpT', xlabel='dpT between sim and reco',
                                                              proc=proc, config=config)

        for wp in wps:
            figures[f'hit_v_layer_BPIX_{proc}_{wp}_{config}'] = plot_hit_vs_layer((proc, wp, config), 'barrel')
            figures[f'hit_v_layer_FPIX_{proc}_{wp}_{config}'] = plot_hit_vs_layer((proc, wp, config), 'forward')

    for proc, config in product(procs, configs):
        bins_pt = np.linspace(0, 100, 30)
        figures[f'good_sim_kinem_{proc}_{config}'] = plot_kinematic_eff('good_sim', proc, config, norm=1, ylim=(0, None),
                                                                        bins_eta=30, bins_phi=30, bins_pt=bins_pt)
        figures[f'good_sim_kinem2d_{proc}_{config}'] = plot_kinematic_eff2d('good_sim', proc, config)

        figures[f'gsf_track_kinem_{proc}_{config}'] = plot_kinematic_eff('gsf_track', proc, config, norm=1, ylim=(0, None),
                                                                         bins_eta=30, bins_phi=30, bins_pt=bins_pt)
        figures[f'gsf_track_kinem2d_{proc}_{config}'] = plot_kinematic_eff2d('gsf_track', proc, config)

        figures[f'gsf_high_pt1_kinem_{proc}_{config}'] = plot_kinematic_eff('gsf_high_pt1', proc, config,
                                                                            norm=1, ylim=(0, None), bins_eta=30,
                                                                            bins_phi=30, bins_pt=bins_pt)
        figures[f'gsf_high_pt1_kinem2d_{proc}_{config}'] = plot_kinematic_eff2d('gsf_high_pt1', proc, config)

        figures[f'gsf_high_pt2_kinem_{proc}_{config}'] = plot_kinematic_eff('gsf_high_pt2', proc, config,
                                                                            norm=1, ylim=(0, None), bins_eta=30,
                                                                            bins_phi=30, bins_pt=bins_pt)
        figures[f'gsf_high_pt2_kinem2d_{proc}_{config}'] = plot_kinematic_eff2d('gsf_high_pt2', proc, config)


    # for proc, wp, region in product(procs, wps, ('BPIX', 'FPIX')):
    #     figures[f'hit_v_layer_{region}_{wp}_{proc}'] = plot_hit_vs_layer((proc, wp), region)

    for proc, config in product(procs, configs):
        def add_num_den(key, func, args, kwargs):
            key = f'{key}_{proc}_{config}'
            base_ext = kwargs.get('ext', '')
            bins_pt_ = [0, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100, 150, 200, 300]
            kwargs['bins_pt'] = kwargs.get('bins_pt', bins_pt_)
            kwargs['bins_eta'] = kwargs.get('bins_eta', 15)
            kwargs['bins_phi'] = kwargs.get('bins_phi', 15)
            args = *args, proc, config
            figures[key] = func(*args, ylim=(0, 1.1), is_ratio=True, **kwargs)
            kwargs_ = kwargs.copy()
            kwargs_['ext'] = base_ext+'_num'
            figures[key+'_num'] = func(*args, ylim=(0, None), **kwargs_)
            kwargs_ = kwargs.copy()
            kwargs_['ext'] = base_ext+'_den'
            figures[key+'_den'] = func(*args, ylim=(0, None), **kwargs_)

        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('prompt_eff_dR', plot_kinematic_eff, ('prompt_eff',), dict(ext='_dR', incl_sel=False))
        # add_num_den('prompt_pur_dR', plot_kinematic_eff, ('prompt_pur',), dict(ext='_dR', incl_sel=False))
        # add_num_den('nonprompt_eff_dR', plot_kinematic_eff, ('nonprompt_eff',), dict(ext='_dR', incl_sel=False))
        # add_num_den('nonprompt_pur_dR', plot_kinematic_eff, ('nonprompt_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('fake_rate_no_e_match_incl', plot_kinematic_eff, ('fake_rate_no_e_match_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('fake_rate_no_e_match', plot_kinematic_eff, ('fake_rate_no_e_match',), 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(['dy', 'tt'], ['new-default-noskip', 'new-wide-noskip'],
    #                                                    hit_layers, ['dPhi', 'dRz'], ['BPIX', 'FPIX']):
    #     figures[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(['dy', 'tt'], [1, 2, 3], ['dPhi', 'dRz']):
    #     figures[f'resall_H{hit}_{var}_{proc}'] = plot_res_contour(proc, hit, var, rel_layers[hit])

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

    def add_simple_plot(proc, wp, config, plot_name, xlabel, ylabel, is2d=False, xlim=None, ylim=None, clear_zero=False):
        if is2d:
            figures[f'{plot_name}_{proc}_{wp}_{config}'] = \
                simple_dist2d(plot_name, proc, wp, config, colz_fmt='', xlabel=xlabel, ylabel=ylabel, do_profile=True,
                              xlim=xlim, ylim=ylim, clear_zero=clear_zero)
        else:
            figures[f'{plot_name}_{proc}_{wp}_{config}'] = \
                simple_dist(plot_name, proc, wp, config, xlabel=xlabel, ylabel=ylabel, xlim=xlim, ylim=ylim)

    for proc, wp, config, in product(procs, wps, configs):
        # add_simple_plot(proc, wp, 'n_pileup')
        add_simple_plot(proc, wp, config, 'n_initseeds_v_PU', '# Initial Seeds', '# In-Time Pileup', True)
        add_simple_plot(proc, wp, config, 'n_seeds_v_PU', '# Seeds', '# In-Time Pileup', True)
        add_simple_plot(proc, wp, config, 'n_good_seeds_v_PU', '# Seeds', '# In-Time Pileup', True, xlim=(-0.5, 40.5))
        add_simple_plot(proc, wp, config, 'n_tracks_v_PU', '# GSF Tracks', '# In-Time Pileup', True, xlim=(-0.5, 20.5))
        # add_simple_plot(proc, wp, config, 'n_pv_v_PU', '# PV', '# In-Time Pileup', True, xlim=(-0.5, 50.5))  # garbage

        # add_simple_plot(proc, wp, config, 'n_seeds_v_tPU', '# Seeds', '# True Pileup', True)
        # add_simple_plot(proc, wp, config, 'n_good_seeds_v_tPU', '# Seeds', '# True Pileup', True, xlim=(-0.5, 40.5))
        # add_simple_plot(proc, wp, config, 'n_tracks_v_tPU', '# GSF Tracks', '# True Pileup', True, xlim=(-0.5, 20.5))
        # add_simple_plot(proc, wp, config, 'n_pv_v_tPU', '# PV', '# True Pileup', True, xlim=(-0.5, 50.5))

        # add_simple_plot(proc, wp, config, 'n_seeds_v_n_scl', '# Seeds', '# Super Clusters', True,
        #                 xlim=(-0.5, 20.5), ylim=(-0.5, 20.5), clear_zero=False)
        add_simple_plot(proc, wp, config, 'n_good_seeds_v_n_scl', '# Good Seeds', '# Super Clusters', True,
                        xlim=(-0.5, 20.5), ylim=(-0.5, 20.5), clear_zero=False)
    #     add_simple_plot(proc, wp, config, 'n_good_seeds_v_n_good_scl', '# Good Seeds', '# Good Super Clusters', True,
    #                     xlim=(-0.5, 20.5), ylim=(-0.5, 20.5), clear_zero=False)
    #     add_simple_plot(proc, wp, config, 'n_initseeds_v_PU', '# Initial Seeds', '# In-Time Pileup', True,
    #                     xlim=None, ylim=None, clear_zero=False)
    #     add_simple_plot(proc, wp, config, 'n_initseeds_v_tPU', '# Initial Seeds', '# True Pileup', True,
    #                     xlim=None, ylim=None, clear_zero=False)

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

    # mpb.generate_report(figures, 'Fixing Hit Skipping',
    #                     output='report.html',
    #                     body='../docs/reports/report_2018_08_16.md')
    if publish:
        mpb.publish()


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

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

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

def style(config):
    return {
        'noskip': '-',
        'skip': '--',
        'skip-pileup': '-.',
    }[config]


if __name__ == '__main__':
    parser = ArgumentParser('Electron Seeding Efficiency Plot Maker')
    parser.add_argument('--build', action='store_true')
    parser.add_argument('--publish', action='store_true')
    args = parser.parse_args()
    set_defaults()
    mpb.configure(output_dir='seeding_studies',
                  multiprocess=True,
                  publish_remote="cfangmeier@t3.unl.edu",
                  publish_dir="/home/dominguez/cfangmeier/public_html/eg/",
                  publish_url="t3.unl.edu/~cfangmeier/eg/",
                  early_abort=True,
                  )
    procs = {
        'dy': r'Drell-Yan',
        'tt': r'$t\bar{t}$'
    }
    configs = [
        'noskip',
        'skip',
        'skip-pileup',
    ]
    wps = {
        'old-default': 'Old-Method Seeding',
        # 'new-narrow': 'HLT Settings',
        'new-default': 'HLT Settings',
        'new-wide': 'Wide Settings',
        'new-extra-wide': 'Extra Wide Settings',
    }
    all_cut_plots(build=args.build, publish=args.publish)