Warning, /detector_benchmarks/benchmarks/backwards_ecal/backwards_ecal.org is written in an unsupported language. File is not indexed.
0001 #+PROPERTY: header-args:jupyter-python :session /jpy:localhost#8888:backwards_ecal :async yes :results drawer :exports both
0002
0003 #+TITLE: ePIC EEEMCal benchmark
0004 #+AUTHOR: Dmitry Kalinkin
0005 #+OPTIONS: d:t
0006
0007 #+LATEX_CLASS_OPTIONS: [9pt,letter]
0008 #+BIND: org-latex-image-default-width ""
0009 #+BIND: org-latex-image-default-option "scale=0.3"
0010 #+BIND: org-latex-images-centered nil
0011 #+BIND: org-latex-minted-options (("breaklines") ("bgcolor" "black!5") ("frame" "single"))
0012 #+LATEX_HEADER: \usepackage[margin=1in]{geometry}
0013 #+LATEX_HEADER: \setlength{\parindent}{0pt}
0014 #+LATEX: \sloppy
0015
0016 #+begin_src jupyter-python :results silent
0017 import os
0018 from pathlib import Path
0019
0020 import awkward as ak
0021 import numpy as np
0022 import vector
0023 import uproot
0024 from sklearn.metrics import roc_curve
0025
0026 vector.register_awkward()
0027 #+end_src
0028
0029 * Plotting setup
0030
0031 #+begin_src jupyter-python :results silent
0032 import matplotlib as mpl
0033 import matplotlib.pyplot as plt
0034
0035 def setup_presentation_style():
0036 mpl.rcParams.update(mpl.rcParamsDefault)
0037 plt.style.use('ggplot')
0038 plt.rcParams.update({
0039 'axes.labelsize': 8,
0040 'axes.titlesize': 9,
0041 'figure.titlesize': 9,
0042 'figure.figsize': (4, 3),
0043 'legend.fontsize': 7,
0044 'xtick.labelsize': 8,
0045 'ytick.labelsize': 8,
0046 'xaxis.labellocation': 'right',
0047 'yaxis.labellocation': 'top',
0048 'pgf.rcfonts': False,
0049 })
0050
0051 setup_presentation_style()
0052 #+end_src
0053
0054 * Parameters
0055
0056 #+begin_src jupyter-python :results silent
0057 DETECTOR_CONFIG=os.environ.get("DETECTOR_CONFIG")
0058 PLOT_TITLE=os.environ.get("PLOT_TITLE")
0059 INPUT_PATH_FORMAT=os.environ.get("INPUT_PATH_FORMAT", "EPIC/RECO/24.04.0/epic_craterlake/SINGLE/{particle}/{energy}/130to177deg/{particle}_{energy}_130to177deg.{ix:04d}.eicrecon.tree.edm4eic.root")
0060 INPUT_PATH_INDEX_RANGE=list(range(20))
0061
0062 output_dir=Path(os.environ.get("OUTPUT_DIR", "./"))
0063 output_dir.mkdir(parents=True, exist_ok=True)
0064 #+end_src
0065
0066 * Analysis
0067
0068 First, we define a requirement on what phase we will consider for our
0069 analysis. The following function filters single-particle events that
0070 are thrown within $-3.5 < \eta < -2.0$:
0071
0072 #+begin_src jupyter-python
0073 def filter_pointing(events):
0074 part_momentum = ak.zip({
0075 "m": events["MCParticles.mass"],
0076 "px": events["MCParticles.momentum.x"],
0077 "py": events["MCParticles.momentum.y"],
0078 "pz": events["MCParticles.momentum.z"],
0079 }, with_name="Momentum4D")
0080 cond = (part_momentum.eta[:,0] > -3.5) & (part_momentum.eta[:,0] < -2.)
0081 return events[cond]
0082 #+end_src
0083
0084 #+begin_src jupyter-python
0085 energies = [
0086 "100MeV",
0087 "200MeV",
0088 "500MeV",
0089 "1GeV",
0090 "2GeV",
0091 "5GeV",
0092 "10GeV",
0093 "20GeV",
0094 ]
0095
0096 filter_name = [
0097 "MCParticles.*",
0098 "*EcalEndcapNRecHits*",
0099 "*EcalEndcapNClusters*",
0100 ]
0101
0102 pi_eval = {}
0103 e_eval = {}
0104
0105 for energy in energies:
0106 pi_eval[energy] = filter_pointing(uproot.concatenate(
0107 {INPUT_PATH_FORMAT.format(particle="pi-", energy=energy, ix=ix): "events" for ix in INPUT_PATH_INDEX_RANGE},
0108 filter_name=filter_name,
0109 ))
0110 e_eval[energy] = filter_pointing(uproot.concatenate(
0111 {INPUT_PATH_FORMAT.format(particle="e-", energy=energy, ix=ix): "events" for ix in INPUT_PATH_INDEX_RANGE},
0112 filter_name=filter_name,
0113 ))
0114 #+end_src
0115
0116 ** Pion rejection
0117
0118 #+begin_src jupyter-python
0119 fig, axs = plt.subplots(2, 4, sharex=True, sharey=True, figsize=(15, 6))
0120 fig_log, axs_log = plt.subplots(2, 4, sharex=True, sharey=True, figsize=(15, 6))
0121 fig_roc, axs_roc = plt.subplots(2, 4, sharex=True, sharey=True, figsize=(15, 6))
0122
0123 fig.suptitle(PLOT_TITLE)
0124 fig_log.suptitle(PLOT_TITLE)
0125 fig_roc.suptitle(PLOT_TITLE)
0126
0127 axs = np.ravel(np.array(axs))
0128 axs_log = np.ravel(np.array(axs_log))
0129 axs_roc = np.ravel(np.array(axs_roc))
0130
0131 rocs = {}
0132
0133 for ix, energy in enumerate(energies):
0134 for use_clusters in [False, True]:
0135 energy_value = float(energy.replace("GeV", "").replace("MeV", "e-3"))
0136 if use_clusters:
0137 clf_label = "leading cluster"
0138 else:
0139 clf_label = "sum all hits"
0140 def clf(events):
0141 if use_clusters:
0142 return ak.drop_none(ak.max(events["EcalEndcapNClusters.energy"], axis=-1)) / energy_value
0143 else:
0144 return ak.sum(events["EcalEndcapNRecHits.energy"], axis=-1) / energy_value
0145 e_pred = clf(e_eval[energy])
0146 pi_pred = clf(pi_eval[energy])
0147
0148 for do_log, ax in [(False, axs[ix]), (True, axs_log[ix])]:
0149 plt.sca(ax)
0150 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.)
0151 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")
0152 plt.title(f"{energy}")
0153 plt.legend()
0154 plt.xlabel("Classifier output")
0155 plt.ylabel("Event yield")
0156 if do_log:
0157 plt.yscale("log")
0158
0159 plt.sca(axs_roc[ix])
0160 fpr, tpr, _ = roc_curve(
0161 np.concatenate([np.ones_like(e_pred), np.zeros_like(pi_pred)]),
0162 np.concatenate([e_pred, pi_pred]),
0163 )
0164 cond = fpr != 0 # avoid infinite rejection (region of large uncertainty)
0165 cond &= tpr != 1 # avoid linear interpolation (region of large uncertainty)
0166 def mk_interp(tpr, fpr):
0167 def interp(eff):
0168 return np.interp(eff, tpr, fpr)
0169 return interp
0170 rocs.setdefault(clf_label, {})[energy] = mk_interp(tpr, fpr)
0171 plt.plot(tpr[cond] * 100, 1 / fpr[cond], label=f"{clf_label}")
0172 plt.yscale("log")
0173 plt.title(f"{energy}")
0174 plt.legend(loc="lower left")
0175 plt.xlabel("Electron efficiency, %")
0176 plt.ylabel("Pion rejection factor")
0177
0178 fig.savefig(output_dir / f"pred.pdf", bbox_inches="tight")
0179 plt.close(fig)
0180 fig_log.savefig(output_dir / f"pred_log.pdf", bbox_inches="tight")
0181 fig_log.show()
0182 fig_roc.savefig(output_dir / f"roc.pdf", bbox_inches="tight")
0183 fig_roc.show()
0184
0185 plt.figure()
0186 for clf_label, roc in rocs.items():
0187 plt.plot(
0188 [float(energy.replace("GeV", "").replace("MeV", "e-3")) for energy in energies],
0189 [1 / roc[energy](0.95) for energy in energies],
0190 marker=".",
0191 label=f"{clf_label}",
0192 )
0193 plt.yscale("log")
0194 plt.title(INPUT_PATH_FORMAT)
0195 plt.legend()
0196 plt.xlabel("Energy, GeV")
0197 plt.ylabel("Pion rejection at 95%")
0198 plt.savefig(output_dir / f"pion_rej.pdf", bbox_inches="tight")
0199 plt.show()
0200 #+end_src
