diff --git a/benchmarks/backwards_ecal/backwards_ecal.org b/benchmarks/backwards_ecal/backwards_ecal.org
index fb0792161abbeb3f3c1ea8837cf9177dc781fb0b..ec359464f7e6563cf3055b530bb595f08ffd72b2 100644
--- a/benchmarks/backwards_ecal/backwards_ecal.org
+++ b/benchmarks/backwards_ecal/backwards_ecal.org
@@ -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")