In [113]:
import numpy as np
import matplotlib as mpl
import mplhep
import mplhep.error_estimation
mpl.style.use(mplhep.style.ROOT)
import matplotlib.pyplot as plt
def poisson_interval_ignore_empty(sumw, sumw2):
interval = mplhep.error_estimation.poisson_interval(sumw, sumw2)
lo, hi = interval[0,...], interval[1,...]
to_ignore = np.isnan(lo)
lo[to_ignore] = 0.0
hi[to_ignore] = 0.0
return np.array([lo,hi])
def histplot(h, scale = None, **kwargs):
for disallowed in ['bins','w2method','w2', 'density']:
if disallowed in kwargs:
raise KeyError(f'Cannot manually pass {disallowed} to our histplot')
values, variances = h.values(), h.variances()
if scale == 'density':
kwargs['density'] = True
elif scale == 'bin sum':
variances[values > 0] = variances[values > 0] / values[values > 0]
values = values / values.sum()
variances = values * variances
elif isinstance(scale, float):
values = scale*values
variances = scale*variances
return mplhep.histplot(
values,
bins = h.axes[0].edges,
w2method = poisson_interval_ignore_empty,
w2 = variances,
**kwargs
)
def show(
beam, *,
stage = 'Preliminary',
ax = None,
filename = None,
display = True,
exp_loc = 0
):
mplhep.label.lumitext(f'{beam}GeV', ax = ax)
mplhep.label.exp_text('LDMX',stage, loc=exp_loc, ax=ax)
if filename is not None:
plt.savefig(
filename,
bbox_inches='tight'
)
if display:
plt.show()
else:
plt.clf()
In [114]:
import pickle
import hist
In [115]:
mass_points_to_plot = [1,1000]
In [128]:
from pathlib import Path
def load_hists(run = None):
run_dir = Path('analysis/coffea/hists/re-ana')
run_opts = [
f
for f in run_dir.iterdir()
if f.suffix == '.pkl'
]
# sort by last modified time
run_opts.sort(key = lambda f: f.stat().st_mtime)
if run is None:
run = run_opts[-1] # use last run by default
elif isinstance(run, int):
run = run_opts[run]
elif isinstance(run, str):
run = (run_dir / (run+'.pkl'))
if not run.is_file():
raise ValueError(f'{run} does not exist.')
elif not isinstance(run, Path):
raise ValueError(f'{run} is not None, an int, str, or Path.')
print(f'Loading {run}')
with open(run, 'rb') as f:
return pickle.load(f)
In [129]:
h = load_hists()
Loading analysis/coffea/hists/re-ana/4-more-rec-bins.pkl
In [158]:
# temporary until coffea run that actually counts the runs is done
h[4]['enriched-nuclear']['nthrown'] = 1_000_000*(10*6000)
h[8]['enriched-nuclear']['nthrown'] = 1_000_000*(10*5000)
h[4]['dimuon']['nthrown'] = 1_000_000*200
h[8]['dimuon']['nthrown'] = 1_000_000*1200
In [159]:
# scaling factors to 1e13 EoT
for beam in [4,8]:
for bkgd in ['enriched-nuclear','dimuon']:
h[beam][bkgd]['scale'] = (1e13/h[beam][bkgd]['nthrown'])
Readout Threshold Comparison¶
In [70]:
for b in [4,8]:
for rot in [None,0.5]:
h[b]['histograms'][rot]['EcalEnergyFrac']['enriched-nuclear','trigger',:].plot(
flow='none',
label=(
'Nominal'
if rot is None else
'Only Hits with Amplitude above half a MIP'
)
)
plt.ylabel('Events')
plt.xlabel(r"$E_\text{ECAL}/E_\text{Beam}$")
plt.ylim(ymax=10_000)
plt.xlim(0,1)
plt.legend(title='Enriched Nuclear')
show(b, filename=f'plots/readout-threshold/energy-after-trigger-{b}.pdf')
Yield¶
In [206]:
# yield_table_rows = []
for b in [4,8]:
for rot in [None,0.5]:
print(b,'GeV', rot)
eef_h = h[b]['histograms'][rot]['EcalEnergyFrac']
bkg_h = (
eef_h['dimuon',...]*h[b]['dimuon']['scale']
+eef_h['enriched-nuclear',...]*h[b]['enriched-nuclear']['scale']
)
for cut in bkg_h.axes[0]:
total = bkg_h[cut,sum]
# row = (b,rot,cut,total.value,total.variance)
print(cut, f'{total.value:.2e}', f'{total.variance:.2e}',end='')
for m in [1,100,1000]: #5,10,50,100,500,1000]:
total_sim = h[b][f'ap{m}']['weightsum']
total_pass_trig = eef_h[f'ap{m}','trigger',sum].value
print(f' {eef_h[f'ap{m}',cut,sum].value/total_sim*100:.1f}',end='')
print()
4 GeV None no-selection 1.37e+08 6.42e+08 100.0 100.0 100.0 trigger 4.60e+07 2.23e+08 57.8 70.6 83.1 analysis_trigger 1.95e+06 1.13e+07 34.9 52.7 71.8 hcal_max_pe 1.15e+03 9.59e+02 34.2 51.8 69.0 ecal_rms 1.22e+02 1.02e+02 27.5 40.5 32.4 4 GeV 0.5 no-selection 1.37e+08 6.42e+08 100.0 100.0 100.0 trigger 4.60e+07 2.23e+08 57.8 70.6 83.1 analysis_trigger 1.98e+06 1.15e+07 35.1 52.9 71.9 hcal_max_pe 1.20e+03 1.00e+03 34.3 51.9 69.1 ecal_rms 1.34e+02 1.12e+02 27.7 40.8 32.8 8 GeV None no-selection 1.48e+08 8.54e+08 100.0 100.0 100.0 trigger 6.10e+07 3.50e+08 66.0 78.9 89.3 analysis_trigger 6.88e+06 3.52e+07 52.0 68.8 84.1 hcal_max_pe 3.18e+01 3.17e+01 50.0 66.5 81.0 8 GeV 0.5 no-selection 1.48e+08 8.54e+08 100.0 100.0 100.0 trigger 6.10e+07 3.50e+08 66.0 78.9 89.3 analysis_trigger 6.95e+06 3.56e+07 52.2 68.9 84.2 hcal_max_pe 3.18e+01 3.17e+01 50.2 66.7 81.1
Recon ECal Energy after Certain Selections¶
In [118]:
def signal_rec_efrac(
beam,
*,
masses = [1,10,100,1000],
selection = 'no-selection',
max_efrac = 1.0,
filename = None,
density = True
):
hh = h[beam]['histograms'][0.5]['EcalEnergyFrac']
for m in masses:
histplot(
hh[f'ap{m}',selection,:hist.loc(max_efrac)],
label=r"$m_{A'} = $"+f'{m}MeV',
scale='density' if density else None
)
plt.yscale('log')
if density:
plt.ylabel('Events (Integral Normalized to 1)')
else:
plt.ylabel('Events')
tt = (3160 / 8000 if beam == 8 else 1500 / 4000)
plt.axvline(tt, color='gray')
plt.text(tt-0.01, 1e-1, 'Below Trigger\nThreshold', color='gray', va='top', ha='right')
plt.xlabel(r'$E_\text{ECAL}/E_\text{Beam}$')
plt.xlim(0.0, max_efrac)
plt.ylim(ymin=None, ymax=2e1)
plt.legend()
show(beam, filename = filename)
signal_rec_efrac(4, filename='plots/selection/energy-signal-no-selection-4gev.pdf')
signal_rec_efrac(8, filename='plots/selection/energy-signal-no-selection-8gev.pdf')
In [162]:
def rec_efrac_density(
beam,
selection,
*,
max_efrac = 0.395,
filename = None
):
hh = h[beam]['histograms'][0.5]['EcalEnergyFrac']
bkgd_h = (hh['dimuon',...]*h[beam]['dimuon']['scale']+hh['enriched-nuclear',...]*h[beam]['enriched-nuclear']['scale'])
histplot(
bkgd_h[selection,:hist.loc(max_efrac)],
label='Background',
scale='density'
)
for m in mass_points_to_plot:
histplot(
hh[f'ap{m}',selection,:hist.loc(max_efrac)],
label=r"$m_{A'} = $"+f'{m}MeV',
scale='density'
)
plt.yscale('log')
plt.ylabel('Events (Integral Normalized to 1)')
plt.xlabel(r'$E_\text{ECAL}/E_\text{Beam}$')
plt.xlim(0.0, max_efrac)
plt.legend()
show(beam, filename = filename)
rec_efrac_density(4, 'trigger', max_efrac = 0.375, filename='plots/selection/energy-after-trigger-4gev.pdf')
rec_efrac_density(8, 'trigger', filename='plots/selection/energy-after-trigger-8gev.pdf')
In [166]:
def load_signal_rates(f = 'prod-rate/eat-nom-rates.csv'):
import pandas as pd
rates = pd.read_csv(f)
rates['rate'] = rates.rate/0.01**2 # remove eps2 dependence
return rates
signal_rates = load_signal_rates().set_index(['target','beam','min_efrac','map'])
In [205]:
max_efrac = 0.375
for beam, selection, ymax in [
(4, 'ecal_rms', None),
(8, 'hcal_max_pe', 5e1)
]:
hh = h[beam]['histograms'][0.5]['EcalEnergyFrac']
bkgd_h = (
hh['dimuon',...]*h[beam]['dimuon']['scale']
+hh['enriched-nuclear',...]*h[beam]['enriched-nuclear']['scale']
)
histplot(
bkgd_h[selection,:hist.loc(max_efrac)],
label='Background',
)
for m, e in [
(1, 5e-7),
(10, 5e-6),
(100, 5e-5),
(1000, 1e-2)
]:
no_selection_weightsum = hh[f'ap{m}','no-selection',...].values()
# to avoid divide by zero, I just replace the zeros with 1 since I know that
# the numerator will also be zero in this case.
no_selection_weightsum[no_selection_weightsum == 0.0] = 1.0
selection_eff = hh[f'ap{m}',selection,...]/no_selection_weightsum
signal_event_yield = (
selection_eff
*signal_rates.loc['eat',beam,0.5,m].iloc[0]
*e**2
*1e13
)
histplot(
signal_event_yield[:hist.loc(max_efrac)],
label=r"$m_{A'} = $"+f'{m}MeV'+r"; $\epsilon = $"+f'{e:.0e}'
)
plt.yscale('log')
plt.ylabel('Events in $10^{13}$ EoT')
plt.xlabel(r'$E_\text{ECAL}/E_\text{Beam}$')
plt.xlim(0.0, max_efrac)
plt.ylim(ymax=ymax)
# plt.ylim(ymin=None, ymax=5e1)
plt.legend()
show(beam, filename = f'plots/selection/energy-after-full-selection-{beam}gev.pdf')
N-1 Plots¶
In [74]:
def nm1(
beam,
selection,
cut_val,
var_name,
*,
filename = None,
rebin = 5
):
hh = h[beam]['histograms'][0.5][f'N-1_{selection}']
(hh['dimuon',:]+hh['enriched-nuclear',:])[hist.rebin(rebin)].plot(
density=True,
label='Background',
flow='sum'
)
for m in mass_points_to_plot:
hh[f'ap{m}',hist.rebin(rebin)].plot(density=True,label=r"$m_{A'} = $"+f'{m}MeV',flow='sum')
plt.xlim(
np.min(hh.axes[1].edges),
np.max(hh.axes[1].edges)
)
plt.xlabel(var_name)
plt.yscale('log')
plt.ylim(ymin=None,ymax=2)
plt.ylabel('Events (Integral Normalized to 1)')
plt.legend(title=r'$N-1$')
plt.axvline(cut_val, color='gray')
show(beam, filename=filename)
nm1(8,'hcal_max_pe',10,r"$\text{max}(\text{PE}_\text{HCal})$", filename='plots/selection/nm1-hcal-max-pe-8gev.pdf')
nm1(4,'hcal_max_pe',10,r"$\text{max}(\text{PE}_\text{HCal})$", filename='plots/selection/nm1-hcal-max-pe-4gev.pdf')
# nm1(4,'ecal_n_hits',28,r'$N_\text{hits-ECal}$', rebin=2, filename='plots/selection/nm1-ecal-nhits-4gev.pdf')
nm1(4,'ecal_rms',20,'ECal Hit RMS [mm]', rebin=1, filename='plots/selection/nm1-ecal-rms-4gev.pdf')
In [75]:
nm1_rms_h = (
h[4]['histograms'][0.5]['N-1_ecal_rms']['dimuon',...]
+h[4]['histograms'][0.5]['N-1_ecal_rms']['enriched-nuclear',...]
)
content = nm1_rms_h.view(flow=True)
cumulants = np.add.accumulate(content['value'])
rel_variance = content['variance']
rel_variance[content['value']>0] /= content['value'][content['value']>0]
cumulant_variance = cumulants*np.add.accumulate(rel_variance)
mplhep.histplot(
cumulants,
bins = np.concatenate(([-1],nm1_rms_h.axes[0].edges,[51])),
w2method = poisson_interval_ignore_empty,
w2 = cumulant_variance,
)
plt.yscale('log')
plt.ylabel('Events with RMS < Cut')
plt.xlabel('ECal Hit RMS Cut / mm')
show(4)
In [ ]: