#!/usr/bin/env python import numpy as np import matplotlib.pyplot as plt from filval.result_set import ResultSet from filval.histogram import hist, hist_integral, hist_rebin, hist_norm, hist2d, hist2d_percent_contour from filval.plotting import (decl_plot, render_plots, hist_plot, hist2d_plot, Plot, generate_dashboard, simple_plot) matching_cuts = { 'extra-narrow-window': [ 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, ) ], 'narrow-window': [ 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, ) ], 'wide-window': [ 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, ) ], 'extra-wide-window': [ 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, ) ], 'nwp-tight-window': [ dict( dPhiMaxHighEt=0.025, dPhiMaxHighEtThres=20.0, dPhiMaxLowEtGrad=-0.002, dRzMaxHighEt=9999.0, dRzMaxHighEtThres=0.0, dRzMaxLowEtGrad=0.0, ), dict( dPhiMaxHighEt=0.005, dPhiMaxHighEtThres=0.0, dPhiMaxLowEtGrad=0.0, dRzMaxHighEt=0.07, dRzMaxHighEtThres=30.0, dRzMaxLowEtGrad=-0.002, ), dict( dPhiMaxHighEt=0.006, dPhiMaxHighEtThres=20.0, dPhiMaxLowEtGrad=-0.0001, dRzMaxHighEt=0.08, dRzMaxHighEtThres=30.0, dRzMaxLowEtGrad=-0.002, ) ], 'nwp-window': [ dict( dPhiMaxHighEt=0.05, dPhiMaxHighEtThres=20.0, dPhiMaxLowEtGrad=-0.002, dRzMaxHighEt=9999.0, dRzMaxHighEtThres=0.0, dRzMaxLowEtGrad=0.0, ), dict( dPhiMaxHighEt=0.005, dPhiMaxHighEtThres=0.0, dPhiMaxLowEtGrad=0.0, dRzMaxHighEt=0.07, dRzMaxHighEtThres=30.0, dRzMaxLowEtGrad=-0.002, ), dict( dPhiMaxHighEt=0.006, dPhiMaxHighEtThres=20.0, dPhiMaxLowEtGrad=-0.0001, dRzMaxHighEt=0.08, dRzMaxHighEtThres=30.0, dRzMaxLowEtGrad=-0.002, ) ], 'nwp-eta-breakdown': [ dict( dPhiMaxHighEt=[0.05, 0.07, 0.06], dPhiMaxHighEtThres=[25.0, 25.0, 25.0], dPhiMaxLowEtGrad=[-0.002, -0.006, -0.002], dRzMaxHighEt=[9999.0, 9999.0, 9999.0], dRzMaxHighEtThres=[0.0, 0.0, 0.0], dRzMaxLowEtGrad=[0.0, 0.0, 0.0], etaBins = [1.1, 1.8] ), dict( dPhiMaxHighEt=[0.0035, 0.006, 0.007], dPhiMaxHighEtThres=[0.0, 0.0, 0.0], dPhiMaxLowEtGrad=[0.0, 0.0, 0.0], dRzMaxHighEt=[0.045, 0.08, 0.045], dRzMaxHighEtThres=[30.0, 30.0, 30.0], dRzMaxLowEtGrad=[-0.002, -0.006, -0.002], etaBins=[1.4, 2.3] ), dict( dPhiMaxHighEt=[0.006, 0.007, 0.007], dPhiMaxHighEtThres=[0.0, 20, 20], dPhiMaxLowEtGrad=[0.0, -0.0002, -0.0002], dRzMaxHighEt=[0.04, 0.10, 0.60], dRzMaxHighEtThres=[25.0, 25.0, 25.0], dRzMaxLowEtGrad=[-0.007, -0.007, -0.007], etaBins=[1.0, 2.0] ) ], } 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 center_text(x, y, txt, **kwargs): plt.text(x, y, txt, horizontalalignment='center', verticalalignment='center', transform=plt.gca().transAxes, **kwargs) def hist_integral_ratio(num, den): num_int = hist_integral(num, times_bin_width=False) den_int = hist_integral(den, times_bin_width=False) ratio = num_int / den_int error = np.sqrt(den_int) / den_int # TODO: Check this definition of error return ratio, error @decl_plot def plot_residuals(rs, layer, hit, variable, subdet, cut_sel=None): h_real = hist2d(getattr(rs, f'{variable}_{subdet}_L{layer}_H{hit}_v_Et_TrackMatched')) h_fake = hist2d(getattr(rs, f'{variable}_{subdet}_L{layer}_H{hit}_v_Et_NoMatch')) def do_plot(h): hist2d_plot(h, colorbar=True) 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') if cut_sel: ets = h[3][:, 0] cuts = [calc_window(et, 0, hit, variable, cut_sel) 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') @decl_plot def plot_residuals_eta(rs, hit, variable): h = hist2d(getattr(rs, f'{variable}_residuals_v_eta_H{hit}')) hist2d_plot(h, colorbar=True) 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$ (GeV)') @decl_plot def plot_seed_eff(rs): r"""## ECAL-Driven Seeding Efficiency The proportion of gen-level electrons originating in the luminous region that have an associated Seed, matched via rechit-simhit associations in the pixel detector. Cuts are on simtrack quantities. """ ax_pt = plt.subplot(221) ax_eta = plt.subplot(222) ax_phi = plt.subplot(223) errors = True plt.sca(ax_pt) hist_plot(hist(rs.seed_eff_v_pt), include_errors=errors) center_text(0.5, 0.3, r'$|\eta|<2.4$') plt.xlabel(r"Sim-Track $p_T$") plt.ylim((0, 1.1)) plt.sca(ax_eta) hist_plot(hist(rs.seed_eff_v_eta), include_errors=errors) center_text(0.5, 0.3, r'$p_T>20$') plt.xlabel(r"Sim-Track $\eta$") plt.ylim((0, 1.1)) plt.sca(ax_phi) hist_plot(hist(rs.seed_eff_v_phi), include_errors=errors) center_text(0.5, 0.3, r'$p_T>20$ and $|\eta|<2.4$') plt.xlabel(r"Sim-Track $\phi$") plt.ylim((0, 1.1)) @decl_plot def plot_tracking_eff(rs): r"""## GSF Tracking Efficiency The proportion of electrons origination in the luminous region from the that have an associated GSF track. Cuts are on simtrack quantities. """ ax_pt = plt.subplot(221) ax_eta = plt.subplot(222) ax_phi = plt.subplot(223) ax_eta_pt = plt.subplot(224) errors = True plt.sca(ax_pt) hist_plot(hist(rs.tracking_eff_v_pt), include_errors=errors) center_text(0.5, 0.3, r'$|\eta|<2.4$') plt.xlabel(r"Sim-Track $p_T$") plt.ylim((0, 1.1)) plt.sca(ax_eta) hist_plot(hist(rs.tracking_eff_v_eta), include_errors=errors) center_text(0.5, 0.3, r'$p_T>20$') plt.xlabel(r"Sim-Track $\eta$") plt.ylim((0, 1.1)) plt.sca(ax_phi) hist_plot(hist(rs.tracking_eff_v_phi), include_errors=errors) center_text(0.5, 0.3, r'$p_T>20$ and $|\eta|<2.4$') plt.xlabel(r"Sim-Track $\phi$") plt.ylim((0, 1.1)) plt.sca(ax_eta_pt) hist2d_plot(hist2d(rs.tracking_eff_v_eta_pt)) plt.xlabel(r"Sim-Track $\eta$") plt.ylabel(r"Sim-Track $p_T$") plt.colorbar() @decl_plot def plot_seed_purity(rs, ext=""): r"""## ECAL-Driven Seed Purity The proportion of ECAL-driven seeds that have a matched gen-level electron originating in the luminous region. Cuts are on seed quantities. """ ax_pt = plt.subplot(221) ax_eta = plt.subplot(222) ax_phi = plt.subplot(223) def get_hist(base_name): return hist(getattr(rs, base_name+ext)) errors = True plt.sca(ax_pt) hist_plot(get_hist("seed_pur_v_pt"), include_errors=errors) center_text(0.5, 0.3, r'$|\eta|<2.4$') plt.xlabel(r"Seed $p_T$") if not ext: plt.ylim((0, 1.1)) plt.sca(ax_eta) hist_plot(get_hist("seed_pur_v_eta"), include_errors=errors) center_text(0.5, 0.3, r'$p_T>20$') plt.xlabel(r"Seed $\eta$") if not ext: plt.ylim((0, 1.1)) plt.sca(ax_phi) hist_plot(get_hist("seed_pur_v_phi"), include_errors=errors) center_text(0.5, 0.3, r'$p_T>20$ and $|\eta|<2.4$') plt.xlabel(r"Seed $\phi$") if not ext: plt.ylim((0, 1.1)) @decl_plot def plot_track_purity(rs, ext=""): r"""## GSF Track Purity The proportion of GSF-tracks w\ ECAL-driven seeds that have a matched gen-level electron originating in the luminous region. Cuts are on GSF track quantities. """ ax_pt = plt.subplot(221) ax_eta = plt.subplot(222) ax_phi = plt.subplot(223) def get_hist( base_name): return hist(getattr(rs, base_name+ext)) errors = True plt.sca(ax_pt) hist_plot(get_hist("tracking_pur_v_pt"), include_errors=errors) center_text(0.5, 0.3, r'$|\eta|<2.4$') plt.xlabel(r"GSF-Track $p_T$") plt.ylim((0, 1.1)) plt.sca(ax_eta) hist_plot(get_hist("tracking_pur_v_eta"), include_errors=errors) center_text(0.5, 0.3, r'$p_T>20$') plt.xlabel(r"GSF-Track $\eta$") plt.ylim((0, 1.1)) plt.sca(ax_phi) hist_plot(get_hist("tracking_pur_v_phi"), include_errors=errors) center_text(0.5, 0.3, r'$p_T>20$ and $|\eta|<2.4$') plt.xlabel(r"GSF-Track $\phi$") plt.ylim((0, 1.1)) @decl_plot def plot_hit_vs_layer(rs, region): h = hist2d(getattr(rs, f'hit_vs_layer_{region}')) hist2d_plot(h, txt_format='{:2.0f}') plt.xlabel('Layer \#') plt.ylabel('Hit \#') def single_cut_plots(cut_sel): rs = ResultSet(f'{cut_sel}', f'../hists/{cut_sel}.root') seed_eff = plot_seed_eff, (rs,) tracking_eff = plot_tracking_eff, (rs,) seed_pur = plot_seed_purity, (rs,) track_pur = plot_track_purity, (rs,) track_pur_seed_match = plot_track_purity, (rs,), dict(ext='2') BPIX_residuals_L1_H1_dPhi = plot_residuals, (rs, 1, 1, 'dPhi', 'BPIX'), dict(cut_sel=cut_sel) BPIX_residuals_L2_H2_dPhi = plot_residuals, (rs, 2, 2, 'dPhi', 'BPIX'), dict(cut_sel=cut_sel) BPIX_residuals_L3_H3_dPhi = plot_residuals, (rs, 3, 3, 'dPhi', 'BPIX'), dict(cut_sel=cut_sel) BPIX_residuals_L1_H1_dRz = plot_residuals, (rs, 1, 1, 'dRz', 'BPIX') BPIX_residuals_L2_H2_dRz = plot_residuals, (rs, 2, 2, 'dRz', 'BPIX'), dict(cut_sel=cut_sel) BPIX_residuals_L3_H3_dRz = plot_residuals, (rs, 3, 3, 'dRz', 'BPIX'), dict(cut_sel=cut_sel) BPIX_residuals_L2_H1_dPhi = plot_residuals, (rs, 2, 1, 'dPhi', 'BPIX'), dict(cut_sel=cut_sel) BPIX_residuals_L3_H2_dPhi = plot_residuals, (rs, 3, 2, 'dPhi', 'BPIX'), dict(cut_sel=cut_sel) BPIX_residuals_L4_H3_dPhi = plot_residuals, (rs, 4, 3, 'dPhi', 'BPIX'), dict(cut_sel=cut_sel) BPIX_residuals_L2_H1_dRz = plot_residuals, (rs, 2, 1, 'dRz', 'BPIX') BPIX_residuals_L3_H2_dRz = plot_residuals, (rs, 3, 2, 'dRz', 'BPIX'), dict(cut_sel=cut_sel) BPIX_residuals_L4_H3_dRz = plot_residuals, (rs, 4, 3, 'dRz', 'BPIX'), dict(cut_sel=cut_sel) FPIX_residuals_L1_H1_dPhi = plot_residuals, (rs, 1, 1, 'dPhi', 'FPIX'), dict(cut_sel=cut_sel) FPIX_residuals_L2_H2_dPhi = plot_residuals, (rs, 2, 2, 'dPhi', 'FPIX'), dict(cut_sel=cut_sel) FPIX_residuals_L3_H3_dPhi = plot_residuals, (rs, 3, 3, 'dPhi', 'FPIX'), dict(cut_sel=cut_sel) FPIX_residuals_L1_H1_dRz = plot_residuals, (rs, 1, 1, 'dRz', 'FPIX') FPIX_residuals_L2_H2_dRz = plot_residuals, (rs, 2, 2, 'dRz', 'FPIX'), dict(cut_sel=cut_sel) FPIX_residuals_L3_H3_dRz = plot_residuals, (rs, 3, 3, 'dRz', 'FPIX'), dict(cut_sel=cut_sel) FPIX_residuals_L1_H2_dPhi = plot_residuals, (rs, 1, 2, 'dPhi', 'FPIX'), dict(cut_sel=cut_sel) FPIX_residuals_L1_H3_dPhi = plot_residuals, (rs, 1, 3, 'dPhi', 'FPIX'), dict(cut_sel=cut_sel) FPIX_residuals_L2_H3_dPhi = plot_residuals, (rs, 2, 3, 'dPhi', 'FPIX'), dict(cut_sel=cut_sel) FPIX_residuals_L1_H2_dRz = plot_residuals, (rs, 1, 2, 'dPhi', 'FPIX'), dict(cut_sel=cut_sel) FPIX_residuals_L1_H3_dRz = plot_residuals, (rs, 1, 3, 'dPhi', 'FPIX'), dict(cut_sel=cut_sel) FPIX_residuals_L2_H3_dRz = plot_residuals, (rs, 2, 3, 'dPhi', 'FPIX'), dict(cut_sel=cut_sel) hit_vs_layer_barrel = plot_hit_vs_layer, (rs, 'barrel') hit_vs_layer_forward = plot_hit_vs_layer, (rs, 'forward') dRz_residuals_v_eta_H1 = plot_residuals_eta, (rs, 1, 'dRz') dRz_residuals_v_eta_H2 = plot_residuals_eta, (rs, 2, 'dRz') dRz_residuals_v_eta_H3 = plot_residuals_eta, (rs, 3, 'dRz') dPhi_residuals_v_eta_H1 = plot_residuals_eta, (rs, 1, 'dPhi') dPhi_residuals_v_eta_H2 = plot_residuals_eta, (rs, 2, 'dPhi') dPhi_residuals_v_eta_H3 = plot_residuals_eta, (rs, 3, 'dPhi') plots = [ Plot(BPIX_residuals_L1_H1_dPhi, 'Phi Residuals Layer 1 Hit 1 - BPIX'), Plot(BPIX_residuals_L2_H2_dPhi, 'Phi Residuals Layer 2 Hit 2 - BPIX'), Plot(BPIX_residuals_L3_H3_dPhi, 'Phi Residuals Layer 3 Hit 3 - BPIX'), Plot(BPIX_residuals_L1_H1_dRz, 'dZ Residuals Layer 1 Hit 1 without cuts - BPIX'), Plot(BPIX_residuals_L2_H2_dRz, 'dZ Residuals Layer 2 Hit 2 - BPIX'), Plot(BPIX_residuals_L3_H3_dRz, 'dZ Residuals Layer 3 Hit 3 - BPIX'), Plot(BPIX_residuals_L2_H1_dPhi, 'Phi Residuals Layer 2 Hit 1 - BPIX'), Plot(BPIX_residuals_L3_H2_dPhi, 'Phi Residuals Layer 3 Hit 2 - BPIX'), Plot(BPIX_residuals_L4_H3_dPhi, 'Phi Residuals Layer 4 Hit 3 - BPIX'), Plot(BPIX_residuals_L2_H1_dRz, 'dZ Residuals Layer 2 Hit 1 without cuts - BPIX'), Plot(BPIX_residuals_L3_H2_dRz, 'dZ Residuals Layer 3 Hit 2 - BPIX'), Plot(BPIX_residuals_L4_H3_dRz, 'dZ Residuals Layer 4 Hit 3 - BPIX'), Plot(FPIX_residuals_L1_H1_dPhi, 'Phi Residuals Layer 1 Hit 1 - FPIX'), Plot(FPIX_residuals_L2_H2_dPhi, 'Phi Residuals Layer 2 Hit 2 - FPIX'), Plot(FPIX_residuals_L3_H3_dPhi, 'Phi Residuals Layer 3 Hit 3 - FPIX'), Plot(FPIX_residuals_L1_H1_dRz, 'dR Residuals Layer 1 Hit 1 without cuts - FPIX'), Plot(FPIX_residuals_L2_H2_dRz, 'dR Residuals Layer 2 Hit 2 - FPIX'), Plot(FPIX_residuals_L3_H3_dRz, 'dR Residuals Layer 3 Hit 3 - FPIX'), Plot(FPIX_residuals_L1_H2_dPhi, 'Phi Residuals Layer 1 Hit 2 - FPIX'), Plot(FPIX_residuals_L1_H3_dPhi, 'Phi Residuals Layer 1 Hit 3 - FPIX'), Plot(FPIX_residuals_L2_H3_dPhi, 'Phi Residuals Layer 2 Hit 3 - FPIX'), Plot(FPIX_residuals_L1_H2_dRz, 'dR Residuals Layer 1 Hit 2 - FPIX'), Plot(FPIX_residuals_L1_H3_dRz, 'dR Residuals Layer 1 Hit 3 - FPIX'), Plot(FPIX_residuals_L2_H3_dRz, 'dR Residuals Layer 2 Hit 3 - FPIX'), Plot(seed_eff, 'ECAL-Driven Seeding Efficiency'), Plot(tracking_eff, 'GSF Tracking Efficiency'), Plot(hit_vs_layer_barrel, 'Hit vs Layer - Barrel'), Plot(hit_vs_layer_forward, 'Hit vs Layer - Forward'), Plot(seed_pur, 'ECAL-Driven Seeding Purity'), Plot(track_pur, 'GSF Track Purity'), Plot(track_pur_seed_match , 'GSF Track Purity (Seed Truth Match)'), simple_plot(rs.gsf_tracks_nmatch_sim_tracks, log='y'), Plot(dRz_residuals_v_eta_H1, 'dRz Hit 1 Residuals v eta'), Plot(dRz_residuals_v_eta_H2, 'dRz Hit 2 Residuals v eta'), Plot(dRz_residuals_v_eta_H3, 'dRz Hit 3 Residuals v eta'), Plot(dPhi_residuals_v_eta_H1, 'dPhi Hit 1 Residuals v eta'), Plot(dPhi_residuals_v_eta_H2, 'dPhi Hit 2 Residuals v eta'), Plot(dPhi_residuals_v_eta_H3, 'dPhi Hit 3 Residuals v eta'), ] render_plots(plots, directory='output/figures/'+rs.sample_name, to_disk=to_disk) if not to_disk: generate_dashboard(plots, 'Seeding Efficiency', output=f'{rs.sample_name}.html', source=__file__, config=rs.config) def eta_region_plots(cut_sel): rs = ResultSet(f'{cut_sel}', f'../hists/{cut_sel}.root') @decl_plot def residual_in_region(var, hit): for region in (1, 2, 3): h = hist2d(getattr(rs, f'{var}_residuals_H{hit}_R{region}')) xs, ys = hist2d_percent_contour(h, .99, 'x') plt.plot(xs, ys, label=f'99\%, Region {region}') plt.legend() plt.xlabel({'dPhi': r'$\delta \phi$ (rads)', 'dRz': r'$\delta R/z$ (cm)'}[var]) plt.ylabel('$E_T$ (GeV)') plots = [] for hit in (1, 2, 3): plt_tup = residual_in_region, ('dPhi', hit) plots.append(Plot(plt_tup, f'dPhi residuals, hit {hit}')) plt_tup = residual_in_region, ('dRz', hit) plots.append(Plot(plt_tup, f'dRz residuals, hit {hit}')) render_plots(plots, directory='output/figures/'+rs.sample_name, to_disk=to_disk) if not to_disk: generate_dashboard(plots, 'Breakdown by Eta Region', output=f'{rs.sample_name}-eta-regions.html', source=__file__, config=rs.config) @decl_plot def plot_seed_roc_curve(rss): def get_num_den(rs, basename): num = hist(getattr(rs, f'{basename}_num')) den = hist(getattr(rs, f'{basename}_den')) return hist_integral_ratio(num, den) for rs in rss: eff, eff_err = get_num_den(rs, 'seed_eff_v_phi') pur, pur_err = get_num_den(rs, 'seed_pur_v_phi') if rs.sample_name == 'old-seeding': plt.errorbar([pur], [eff], xerr=[pur_err], yerr=[eff_err], label=rs.sample_name, color='k', marker='o') else: plt.errorbar([pur], [eff], xerr=[pur_err], yerr=[eff_err], label=rs.sample_name[:-7], marker='o') center_text(0.5, 0.3, r'$p_T>20$ and $|\eta|<2.4$') plt.axis('equal') plt.xlim((0.8, 1.0)) plt.ylim((0.8, 1.0)) plt.xlabel('ECAL-Driven Seeding Purity') plt.ylabel('ECAL-Driven Seeding Efficiency') plt.grid() plt.legend() @decl_plot def plot_seed_eff_all(rss): ax_pt = plt.subplot(221) ax_eta = plt.subplot(222) ax_phi = plt.subplot(223) errors = True for rs in rss: plt.sca(ax_pt) hist_plot(hist(rs.seed_eff_v_pt), include_errors=errors, label=rs.sample_name) plt.sca(ax_eta) hist_plot(hist(rs.seed_eff_v_eta), include_errors=errors, label=rs.sample_name) plt.sca(ax_phi) hist_plot(hist(rs.seed_eff_v_phi), include_errors=errors, label=rs.sample_name) plt.sca(ax_pt) center_text(0.5, 0.3, r'$|\eta|<2.4$') plt.xlabel(r"Sim-Track $p_T$") plt.ylim((0, 1.1)) plt.sca(ax_eta) center_text(0.5, 0.3, r'$p_T>20$') plt.xlabel(r"Sim-Track $\eta$") plt.ylim((0, 1.1)) plt.legend(loc='lower right') plt.sca(ax_phi) center_text(0.5, 0.3, r'$p_T>20$ and $|\eta|<2.4$') plt.xlabel(r"Sim-Track $\phi$") plt.ylim((0, 1.1)) @decl_plot def plot_seed_pur_all(rss): ax_pt = plt.subplot(221) ax_eta = plt.subplot(222) ax_phi = plt.subplot(223) errors = True for rs in rss: plt.sca(ax_pt) hist_plot(hist(rs.seed_pur_v_pt), include_errors=errors, label=rs.sample_name) plt.sca(ax_eta) hist_plot(hist(rs.seed_pur_v_eta), include_errors=errors, label=rs.sample_name) plt.sca(ax_phi) hist_plot(hist(rs.seed_pur_v_phi), include_errors=errors, label=rs.sample_name) plt.sca(ax_pt) center_text(0.5, 0.3, r'$|\eta|<2.4$') plt.xlabel(r"Seed $p_T$") plt.ylim((0, 1.1)) plt.sca(ax_eta) center_text(0.5, 0.3, r'$p_T>20$') plt.xlabel(r"Seed $\eta$") plt.ylim((0, 1.1)) plt.legend(loc='lower right') plt.sca(ax_phi) center_text(0.5, 0.3, r'$p_T>20$ and $|\eta|<2.4$') plt.xlabel(r"Seed $\phi$") plt.ylim((0, 1.1)) @decl_plot def plot_tracking_roc_curve(rss, ext=''): def get_num_den(rs, basename): num = hist(getattr(rs, f'{basename}{ext}_num')) den = hist(getattr(rs, f'{basename}{ext}_den')) return hist_integral_ratio(num, den) for rs in rss: eff, eff_err = get_num_den(rs, 'tracking_eff_v_phi') pur, pur_err = get_num_den(rs, 'tracking_pur_v_phi') if rs.sample_name == 'old-seeding': plt.errorbar([pur], [eff], xerr=[pur_err], yerr=[eff_err], label=rs.sample_name, color='k', marker='o') else: plt.errorbar([pur], [eff], xerr=[pur_err], yerr=[eff_err], label=rs.sample_name[:-7], marker='o') center_text(0.5, 0.3, r'$p_T>20$ and $|\eta|<2.4$') plt.axis('equal') plt.xlim((0.8, 1.0)) plt.ylim((0.8, 1.0)) plt.xlabel('GSF-Track Purity') plt.ylabel('GSF-Track Efficiency') plt.grid() plt.legend() @decl_plot def plot_tracking_eff_all(rss, ext=''): ax_pt = plt.subplot(221) ax_eta = plt.subplot(222) ax_phi = plt.subplot(223) errors = True for rs in rss: plt.sca(ax_pt) hist_plot(hist(getattr(rs, f'tracking_eff_v_pt{ext}')), include_errors=errors, label=rs.sample_name) plt.sca(ax_eta) hist_plot(hist(getattr(rs, f'tracking_eff_v_eta{ext}')), include_errors=errors, label=rs.sample_name) plt.sca(ax_phi) hist_plot(hist(getattr(rs, f'tracking_eff_v_phi{ext}')), include_errors=errors, label=rs.sample_name) plt.sca(ax_pt) center_text(0.5, 0.3, r'$|\eta|<2.4$') plt.xlabel(r"Sim-Track $p_T$") plt.ylim((0, 1.1)) plt.sca(ax_eta) center_text(0.5, 0.3, r'$p_T>20$') plt.xlabel(r"Sim-Track $\eta$") plt.ylim((0, 1.1)) plt.legend(loc='lower right') plt.sca(ax_phi) center_text(0.5, 0.3, r'$p_T>20$ and $|\eta|<2.4$') plt.xlabel(r"Sim-Track $\phi$") plt.ylim((0, 1.1)) @decl_plot def plot_tracking_pur_all(rss, ext=''): ax_pt = plt.subplot(221) ax_eta = plt.subplot(222) ax_phi = plt.subplot(223) errors = True for rs in rss: plt.sca(ax_pt) hist_plot(hist(getattr(rs, f'tracking_pur_v_pt{ext}')), include_errors=errors, label=rs.sample_name) plt.sca(ax_eta) hist_plot(hist(getattr(rs, f'tracking_pur_v_eta{ext}')), include_errors=errors, label=rs.sample_name) plt.sca(ax_phi) hist_plot(hist(getattr(rs, f'tracking_pur_v_phi{ext}')), include_errors=errors, label=rs.sample_name) plt.sca(ax_pt) center_text(0.5, 0.3, r'$|\eta|<2.4$') plt.xlabel(r"GSF-Track $p_T$") plt.ylim((0, 1.1)) plt.sca(ax_eta) center_text(0.5, 0.3, r'$p_T>20$') plt.xlabel(r"GSF-Track $\eta$") plt.ylim((0, 1.1)) plt.legend(loc='lower right') plt.sca(ax_phi) center_text(0.5, 0.3, r'$p_T>20$ and $|\eta|<2.4$') plt.xlabel(r"GSF-Track $\phi$") plt.ylim((0, 1.1)) @decl_plot def plot_ecal_rel_res(rss): for rs in rss: hist_plot(hist(rs.ecal_energy_resolution), label=rs.sample_name) plt.xlabel(r"ECAL $E_T$ relative error") plt.legend() @decl_plot def plot_res_contour(rss, hit_number, var, layers, ext='_TrackMatched'): from itertools import chain _, axs = plt.subplots(2, 3) axs_all = list(chain(*axs)) def do_plot(ax, rs): plt.sca(ax) plt.title(rs.sample_name) h = None for layer in layers: subdet = 'BPIX' if layer[0]=='B' else 'FPIX' h = hist2d(getattr(rs, f'{var}_{subdet}_L{layer[1]}_H{hit_number}_v_Et{ext}')) pass xs, ys = hist2d_percent_contour(h, .99, 'x') plt.plot(xs, ys, label=f'{subdet} - L{layer[1]}') ets = h[3][:, 0] cuts = [calc_window(et, 0, hit_number, var, rs.sample_name) for et in ets] plt.plot(cuts, ets, color='red', label='Cut Value') max_x = 0 for ax, rs in zip(axs_all, rss): try: do_plot(ax, rs) _, x_up = ax.get_xlim() max_x = max((max_x, x_up)) except KeyError: pass plt.sca(axs[0][-1]) plt.legend(loc='best') for ax in axs_all: ax.set_xlim((None, max_x)) @decl_plot def number_of_seeds(rss): from filval.histogram import hist_mean for rs in rss: h = hist_rebin(hist(rs.n_seeds), 50, -0.5, 200.5) h = hist_norm(h) mean = int(hist_mean(h)) hist_plot(h, label=f'{rs.sample_name} ($\\mu={mean:d}$)') plt.xlabel('Number of Seeds Produced') plt.legend() def all_cut_plots(cuts): rss = [ResultSet(f'{cut_sel}', f'../hists/{cut_sel}.root') for cut_sel in cuts] tracking_roc_curve = plot_tracking_roc_curve, (rss,) tracking_eff_all = plot_tracking_eff_all, (rss,) tracking_pur_all = plot_tracking_pur_all, (rss,) tracking_roc_curve2 = plot_tracking_roc_curve, (rss, '2') tracking_eff_all2 = plot_tracking_eff_all, (rss, '2') tracking_pur_all2 = plot_tracking_pur_all, (rss, '2') seed_roc_curve = plot_seed_roc_curve, (rss,) seed_eff_all = plot_seed_eff_all, (rss,) seed_pur_all = plot_seed_pur_all, (rss,) ecal_rel_res = plot_ecal_rel_res, (rss,) res_contour_dphi_H1 = plot_res_contour, (rss, 1, 'dPhi', ['B1', 'B2', 'F1']) res_contour_dphi_H2 = plot_res_contour, (rss, 2, 'dPhi', ['B2', 'B3', 'B4', 'F1', 'F2']) res_contour_dRz_H2 = plot_res_contour, (rss, 2, 'dRz', ['B2', 'B3', 'B4', 'F1', 'F2']) res_contour_dphi_H3 = plot_res_contour, (rss, 3, 'dPhi', ['B3', 'B4''F1', 'F2', 'F3']) res_contour_dRz_H3 = plot_res_contour, (rss, 3, 'dRz', ['B3', 'B4''F1', 'F2', 'F3']) plots = [ Plot(tracking_roc_curve, 'Tracking ROC Curve'), Plot(tracking_eff_all, 'Tracking Efficiency'), Plot(tracking_pur_all, 'Tracking Purity'), Plot(tracking_roc_curve2, 'Tracking ROC Curve (Seed Matched)'), Plot(tracking_eff_all2, 'Tracking Efficiency (Seed Matched)'), Plot(tracking_pur_all2, 'Tracking Purity (Seed Matched)'), Plot(seed_roc_curve, 'Seeding ROC Curve'), Plot(seed_eff_all, 'ECAL-Driven Seeding Efficiency'), Plot(seed_pur_all, 'ECAL-Driven Seeding Purity'), Plot(ecal_rel_res, 'ECAL ET Relative Resolution'), Plot(res_contour_dphi_H1, 'dPhi Residual 99% Contours - Hit 1'), Plot(res_contour_dphi_H2, 'dPhi Residual 99% Contours - Hit 2'), Plot(res_contour_dRz_H2, 'dRz Residual 99% Contours - Hit 2'), Plot(res_contour_dphi_H3, 'dPhi Residual 99% Contours - Hit 3'), Plot(res_contour_dRz_H3, 'dRz Residual 99% Contours - Hit 3'), Plot((number_of_seeds, (rss,)), 'Number of Electron Seeds'), ] render_plots(plots, to_disk=to_disk) if not to_disk: generate_dashboard(plots, 'Comparisons', output='comparisons.html', source=__file__, config=rss[0].config) if __name__ == '__main__': to_disk = False all_cuts = [ 'extra-narrow-window', 'narrow-window', 'wide-window', 'extra-wide-window', # 'nwp-window', # 'nwp-tight-window', # 'nwp-eta-breakdown', ] all_cut_plots(all_cuts+['old-seeding']) # for cut in all_cuts: # single_cut_plots(cut) # eta_region_plots('extra-wide-window')