Skip to content
Snippets Groups Projects
Commit fcc6efd0 authored by Dmitry Kalinkin's avatar Dmitry Kalinkin
Browse files

backwards_ecal: drop full calorimeter hit sum plot

parent 34ae6a4e
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
benchmarks/backwards_ecal/backwards_ecal.org
+ 80
− 97
View file @ fcc6efd0
......@@ -91,10 +91,8 @@ energies = [
"10GeV",
"20GeV",
]
filter_name = [
"MCParticles.*",
"*EcalEndcapNRecHits*",
"*EcalEndcapNClusters*",
]
......@@ -128,61 +126,53 @@ sigmas_rel_FWHM_cb = {}
fractions_below = {}
for ix, energy in enumerate(energies):
for use_clusters in [False, True]:
energy_value = float(energy.replace("GeV", "").replace("MeV", "e-3"))
if use_clusters:
clf_label = "leading cluster"
else:
clf_label = "sum all hits"
def clf(events):
if use_clusters:
return ak.drop_none(ak.max(events["EcalEndcapNClusters.energy"], axis=-1)) / energy_value
else:
return ak.sum(events["EcalEndcapNRecHits.energy"], axis=-1) / energy_value
e_pred = clf(e_eval[energy])
plt.sca(axs[ix])
counts, bins, patches = plt.hist(e_pred, weights=np.full_like(e_pred, 1.0 / ak.num(e_pred, axis=0)), bins=np.linspace(0.01, 1.01, 101), label=rf"$e^-$ {clf_label}", hatch=None if use_clusters else r"xxx", alpha=0.8 if use_clusters else 1.)
plt.title(f"{energy}")
energy_value = float(energy.replace("GeV", "").replace("MeV", "e-3"))
clf_label = "leading cluster"
def clf(events):
return ak.drop_none(ak.max(events["EcalEndcapNClusters.energy"], axis=-1)) / energy_value
e_pred = clf(e_eval[energy])
plt.sca(axs[ix])
counts, bins, patches = plt.hist(e_pred, weights=np.full_like(e_pred, 1.0 / ak.num(e_pred, axis=0)), bins=np.linspace(0.01, 1.01, 101), label=rf"$e^-$ {clf_label}")
plt.title(f"{energy}")
e_over_p = (bins[1:] + bins[:-1]) / 2
import scipy.stats
def f(x, n, beta, m, loc, scale):
return n * scipy.stats.crystalball.pdf(x, beta, m, loc, scale)
p0 = (np.sum(counts[10:]), 2., 3., 0.95, 0.05)
e_over_p = (bins[1:] + bins[:-1]) / 2
import scipy.stats
def f(x, n, beta, m, loc, scale):
return n * scipy.stats.crystalball.pdf(x, beta, m, loc, scale)
p0 = (np.sum(counts[10:]), 2., 3., 0.95, 0.05)
try:
import scipy.optimize
par, pcov = scipy.optimize.curve_fit(f, e_over_p[5:], counts[5:], p0=p0, maxfev=10000)
except RuntimeError:
par = None
plt.plot(e_over_p, f(e_over_p, *par), label=rf"Crystal Ball fit", color="tab:green" if use_clusters else "green", lw=0.8)
def summarize_fit(par):
_, _, _, loc_cb, scale_cb = par
# Calculate FWHM
y_max = np.max(f(np.linspace(0., 1., 100), *par))
f_prime = lambda x: f(x, *par) - y_max / 2
x_plus, = scipy.optimize.root(f_prime, loc_cb + scale_cb).x
x_minus, = scipy.optimize.root(f_prime, loc_cb - scale_cb).x
color = "cyan" if use_clusters else "orange"
plt.axvline(x_minus, ls="--", lw=0.75, color=patches[0].get_facecolor(), label=r"$\mu - $FWHM")
plt.axvline(x_plus, ls=":", lw=0.75, color=patches[0].get_facecolor(), label=r"$\mu + $FWHM")
fwhm = (x_plus - x_minus) / loc_cb
sigma_rel_FWHM_cb = fwhm / 2 / np.sqrt(2 * np.log(2))
cutoff_x = loc_cb - fwhm
fraction_below = np.sum(counts[e_over_p < cutoff_x]) / ak.num(e_pred, axis=0)
return sigma_rel_FWHM_cb, fraction_below
sigma_rel_FWHM_cb, fraction_below = summarize_fit(par)
sigmas_rel_FWHM_cb.setdefault(clf_label, {})[energy] = sigma_rel_FWHM_cb
fractions_below.setdefault(clf_label, {})[energy] = fraction_below
try:
import scipy.optimize
par, pcov = scipy.optimize.curve_fit(f, e_over_p[5:], counts[5:], p0=p0, maxfev=10000)
except RuntimeError:
par = None
plt.plot(e_over_p, f(e_over_p, *par), label=rf"Crystal Ball fit", color="tab:green", lw=0.8)
plt.legend()
plt.xlabel("$E/p$", loc="right")
plt.ylabel("Event yield", loc="top")
def summarize_fit(par):
_, _, _, loc_cb, scale_cb = par
# Calculate FWHM
y_max = np.max(f(np.linspace(0., 1., 100), *par))
f_prime = lambda x: f(x, *par) - y_max / 2
x_plus, = scipy.optimize.root(f_prime, loc_cb + scale_cb).x
x_minus, = scipy.optimize.root(f_prime, loc_cb - scale_cb).x
plt.axvline(x_minus, ls="--", lw=0.75, color=patches[0].get_facecolor(), label=r"$\mu - $FWHM")
plt.axvline(x_plus, ls=":", lw=0.75, color=patches[0].get_facecolor(), label=r"$\mu + $FWHM")
fwhm = (x_plus - x_minus) / loc_cb
sigma_rel_FWHM_cb = fwhm / 2 / np.sqrt(2 * np.log(2))
cutoff_x = loc_cb - fwhm
fraction_below = np.sum(counts[e_over_p < cutoff_x]) / ak.num(e_pred, axis=0)
return sigma_rel_FWHM_cb, fraction_below
sigma_rel_FWHM_cb, fraction_below = summarize_fit(par)
sigmas_rel_FWHM_cb.setdefault(clf_label, {})[energy] = sigma_rel_FWHM_cb
fractions_below.setdefault(clf_label, {})[energy] = fraction_below
plt.legend()
plt.xlabel("$E/p$", loc="right")
plt.ylabel("Event yield", loc="top")
fig.savefig(output_dir / f"resolution_plots.pdf", bbox_inches="tight")
fig.savefig(output_dir / f"resolution_plots.png", bbox_inches="tight")
......@@ -248,49 +238,42 @@ axs_roc = np.ravel(np.array(axs_roc))
rocs = {}
for ix, energy in enumerate(energies):
for use_clusters in [False, True]:
energy_value = float(energy.replace("GeV", "").replace("MeV", "e-3"))
if use_clusters:
clf_label = "leading cluster"
else:
clf_label = "sum all hits"
def clf(events):
if use_clusters:
return ak.drop_none(ak.max(events["EcalEndcapNClusters.energy"], axis=-1)) / energy_value
else:
return ak.sum(events["EcalEndcapNRecHits.energy"], axis=-1) / energy_value
e_pred = clf(e_eval[energy])
pi_pred = clf(pi_eval[energy])
for do_log, ax in [(False, axs[ix]), (True, axs_log[ix])]:
plt.sca(ax)
plt.hist(e_pred, weights=np.full_like(e_pred, 1.0 / ak.num(e_pred, axis=0)), bins=np.linspace(0., 1.01, 101), label=rf"$e^-$ {clf_label}", hatch=None if use_clusters else r"xxx", alpha=0.8 if use_clusters else 1.)
plt.hist(pi_pred, weights=np.full_like(pi_pred, 1.0 / ak.num(pi_pred, axis=0)), bins=np.linspace(0., 1.01, 101), label=rf"$\pi^-$ {clf_label}", histtype="step")
plt.title(f"{energy}")
plt.legend()
plt.xlabel("Classifier output")
plt.ylabel("Event yield")
if do_log:
plt.yscale("log")
plt.sca(axs_roc[ix])
fpr, tpr, _ = roc_curve(
np.concatenate([np.ones_like(e_pred), np.zeros_like(pi_pred)]),
np.concatenate([e_pred, pi_pred]),
)
cond = fpr != 0 # avoid infinite rejection (region of large uncertainty)
cond &= tpr != 1 # avoid linear interpolation (region of large uncertainty)
def mk_interp(tpr, fpr):
def interp(eff):
return np.interp(eff, tpr, fpr)
return interp
rocs.setdefault(clf_label, {})[energy] = mk_interp(tpr, fpr)
plt.plot(tpr[cond] * 100, 1 / fpr[cond], label=f"{clf_label}")
plt.yscale("log")
energy_value = float(energy.replace("GeV", "").replace("MeV", "e-3"))
clf_label = "leading cluster"
def clf(events):
return ak.drop_none(ak.max(events["EcalEndcapNClusters.energy"], axis=-1)) / energy_value
e_pred = clf(e_eval[energy])
pi_pred = clf(pi_eval[energy])
for do_log, ax in [(False, axs[ix]), (True, axs_log[ix])]:
plt.sca(ax)
plt.hist(e_pred, weights=np.full_like(e_pred, 1.0 / ak.num(e_pred, axis=0)), bins=np.linspace(0., 1.01, 101), label=rf"$e^-$ {clf_label}")
plt.hist(pi_pred, weights=np.full_like(pi_pred, 1.0 / ak.num(pi_pred, axis=0)), bins=np.linspace(0., 1.01, 101), label=rf"$\pi^-$ {clf_label}", histtype="step")
plt.title(f"{energy}")
plt.legend(loc="lower left")
plt.xlabel("Electron efficiency, %")
plt.ylabel("Pion rejection factor")
plt.legend()
plt.xlabel("Classifier output")
plt.ylabel("Event yield")
if do_log:
plt.yscale("log")
plt.sca(axs_roc[ix])
fpr, tpr, _ = roc_curve(
np.concatenate([np.ones_like(e_pred), np.zeros_like(pi_pred)]),
np.concatenate([e_pred, pi_pred]),
)
cond = fpr != 0 # avoid infinite rejection (region of large uncertainty)
cond &= tpr != 1 # avoid linear interpolation (region of large uncertainty)
def mk_interp(tpr, fpr):
def interp(eff):
return np.interp(eff, tpr, fpr)
return interp
rocs.setdefault(clf_label, {})[energy] = mk_interp(tpr, fpr)
plt.plot(tpr[cond] * 100, 1 / fpr[cond], label=f"{clf_label}")
plt.yscale("log")
plt.title(f"{energy}")
plt.legend(loc="lower left")
plt.xlabel("Electron efficiency, %")
plt.ylabel("Pion rejection factor")
fig.savefig(output_dir / f"pred.pdf", bbox_inches="tight")
fig.savefig(output_dir / f"pred.png", bbox_inches="tight")
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment