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 boost_histogram as bh
0022 import dask_histogram as dh
0023 import numpy as np
0024 import uproot
0025 import vector
0026 #+end_src
0027
0028 #+begin_src jupyter-python :results silent
0029 from dask.distributed import Client
0030
0031 client = Client(os.environ.get("DASK_SCHEDULER", "localhost:8786"))
0032 #+end_src
0033
0034 #+begin_src jupyter-python
0035 vector.register_awkward()
0036
0037 from dask.distributed import WorkerPlugin
0038 class VectorImporter(WorkerPlugin):
0039 idempotent=True
0040
0041 def setup(self, worker):
0042 import vector
0043
0044 vector.register_awkward()
0045
0046 client.register_plugin(VectorImporter())
0047 #+end_src
0048
0049 * Plotting setup
0050
0051 #+begin_src jupyter-python :results silent
0052 import matplotlib as mpl
0053 import matplotlib.pyplot as plt
0054
0055 def setup_presentation_style():
0056 mpl.rcParams.update(mpl.rcParamsDefault)
0057 if int(os.environ.get("USE_PUBLICATION_STYLE", "0")):
0058 import mplhep
0059 plt.style.use([mplhep.style.CMSTex])
0060 plt.rcParams.update({
0061 "text.usetex": True,
0062 "savefig.pad_inches": 0.,
0063 "text.latex.preamble": r"\usepackage{amsmath}",
0064 "pgf.rcfonts": False,
0065 "font.family": "serif",
0066 "font.serif": [],
0067 "font.sans-serif": [],
0068 "mathtext.fontset": "cm",
0069 "figure.figsize": (3.4, 3.4),
0070 "font.size": 10,
0071 "lines.linewidth": 0.25,
0072 "xtick.major.width": 0.4,
0073 "xtick.minor.width": 0.3,
0074 "xtick.major.size": 6.0,
0075 "xtick.minor.size": 3.0,
0076 "ytick.major.width": 0.4,
0077 "ytick.minor.width": 0.3,
0078 "ytick.major.size": 6.0,
0079 "ytick.minor.size": 3.0,
0080 "axes.linewidth": 0.5,
0081 "lines.markersize": 5,
0082 "hatch.linewidth": 0.5
0083 })
0084 else:
0085 plt.style.use("ggplot")
0086 plt.rcParams.update({
0087 "axes.labelsize": 8,
0088 "axes.titlesize": 9,
0089 "figure.titlesize": 9,
0090 "figure.figsize": (4, 3),
0091 "legend.fontsize": 7,
0092 "xtick.labelsize": 8,
0093 "ytick.labelsize": 8,
0094 "xaxis.labellocation": "right",
0095 "yaxis.labellocation": "top",
0096 "pgf.rcfonts": False,
0097 })
0098
0099
0100 setup_presentation_style()
0101 #+end_src
0102
0103 * Parameters
0104
0105 #+begin_src jupyter-python :results silent
0106 DETECTOR_CONFIG=os.environ.get("DETECTOR_CONFIG")
0107 PLOT_TITLE=os.environ.get("PLOT_TITLE")
0108 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.edm4eic.root")
0109
0110 output_dir=Path(os.environ.get("OUTPUT_DIR", "./"))
0111 output_dir.mkdir(parents=True, exist_ok=True)
0112 #+end_src
0113
0114 * Analysis
0115
0116 First, we define a requirement on what phase we will consider for our
0117 analysis. The following function filters single-particle events that
0118 are thrown within $-3.5 < \eta < -2.0$:
0119
0120 #+begin_src jupyter-python
0121 def filter_pointing(events):
0122 part_momentum = ak.zip({
0123 "m": events["MCParticles.mass"],
0124 "px": events["MCParticles.momentum.x"],
0125 "py": events["MCParticles.momentum.y"],
0126 "pz": events["MCParticles.momentum.z"],
0127 }, with_name="Momentum4D")
0128 cond = (part_momentum.eta[:,0] > -3.5) & (part_momentum.eta[:,0] < -2.)
0129 return events[cond]
0130 #+end_src
0131
0132 #+begin_src jupyter-python
0133 energies = [
0134 "100MeV",
0135 "200MeV",
0136 "500MeV",
0137 "1GeV",
0138 "2GeV",
0139 "5GeV",
0140 "10GeV",
0141 "20GeV",
0142 ]
0143 filter_name = [
0144 "MCParticles.*",
0145 "EcalEndcapNClusters.energy",
0146 ]
0147
0148 pi_eval = {}
0149 e_eval = {}
0150
0151 def readlist(path):
0152 with open(path, "rt") as fp:
0153 paths = [line.rstrip() for line in fp.readlines()]
0154 return paths
0155
0156 for energy in energies:
0157 pi_eval[energy] = uproot.dask(
0158 {path: "events" for path in readlist(INPUT_PATH_FORMAT.format(particle="pi-", energy=energy))},
0159 filter_name=filter_name,
0160 steps_per_file=1,
0161 open_files=False,
0162 ).map_partitions(filter_pointing)
0163 e_eval[energy] = uproot.dask(
0164 {path: "events" for path in readlist(INPUT_PATH_FORMAT.format(particle="e-", energy=energy))},
0165 filter_name=filter_name,
0166 steps_per_file=1,
0167 open_files=False,
0168 ).map_partitions(filter_pointing)
0169 #+end_src
0170
0171 ** Energy resolution
0172
0173 The first step in the analysis is plotting distributions of cluster energies.
0174 Only electron distribution is needed, but we also do pi- for the pion rejection
0175 study.
0176
0177 #+begin_src jupyter-python
0178 e_over_p_hist = {}
0179 e_over_p_axis = bh.axis.Regular(110, 0., 1.10)
0180
0181 for ix, energy in enumerate(energies):
0182 for particle_name, dataset in [("e-", e_eval[energy]), ("pi-", pi_eval[energy])]:
0183 energy_value = float(energy.replace("GeV", "").replace("MeV", "e-3"))
0184 def clf(events):
0185 return ak.drop_none(ak.max(events["EcalEndcapNClusters.energy"], axis=-1)) / energy_value
0186 e_over_p_hist[(particle_name, energy)] = dh.factory(
0187 clf(dataset),
0188 axes=(e_over_p_axis,),
0189 )
0190
0191 e_over_p_hist = client.gather(client.compute(e_over_p_hist))
0192 #+end_src
0193
0194 The next step is to plot the histograms together with fits and produce a plot
0195 for resolution vs energy that summarizes them.
0196
0197 #+begin_src jupyter-python
0198 fig, axs = plt.subplots(2, 4, sharex=True, sharey=True, figsize=(15, 6))
0199
0200 axs = np.ravel(np.array(axs))
0201
0202 sigmas_rel_FWHM_cb = {}
0203 fractions_below = {}
0204
0205 def norm_by_sum(arr):
0206 return arr / np.sum(arr)
0207
0208 for ix, energy in enumerate(energies):
0209 plt.sca(axs[ix])
0210 hist = e_over_p_hist[("e-", energy)]
0211 patch = plt.stairs(norm_by_sum(hist.values()), hist.axes[0].edges, label=rf"$e^-$ {PLOT_TITLE}")
0212 plt.title(f"{energy}")
0213
0214 import scipy.stats
0215 def f(x, n, beta, m, loc, scale):
0216 return n * scipy.stats.crystalball.pdf(x, beta, m, loc, scale)
0217 n0 = np.sum(hist.values()[10:])
0218 loc0 = hist.axes[0].centers[5:][np.argmax(hist.values()[5:])]
0219 p0 = (n0, 2., 3., loc0, 0.05)
0220
0221 try:
0222 import scipy.optimize
0223 par, pcov = scipy.optimize.curve_fit(f, hist.axes[0].centers[5:], hist.values()[5:], p0=p0, maxfev=10000)
0224 except RuntimeError:
0225 print("Initial parameters:", p0)
0226 print(hist)
0227 raise
0228 plt.plot(hist.axes[0].centers, f(hist.axes[0].centers, *par), label=rf"Crystal Ball fit", color="tab:green", lw=0.8)
0229
0230 def summarize_fit(par):
0231 _, _, _, loc_cb, scale_cb = par
0232 # Calculate FWHM
0233 y_max = np.max(f(np.linspace(0., 1., 100), *par))
0234 f_prime = lambda x: f(x, *par) - y_max / 2
0235 x_plus, = scipy.optimize.root(f_prime, loc_cb + scale_cb).x
0236 x_minus, = scipy.optimize.root(f_prime, loc_cb - scale_cb).x
0237 plt.axvline(x_minus, ls="--", lw=0.75, color=patch.get_facecolor(), label=r"$\mu - $FWHM")
0238 plt.axvline(x_plus, ls=":", lw=0.75, color=patch.get_facecolor(), label=r"$\mu + $FWHM")
0239 fwhm = (x_plus - x_minus) / loc_cb
0240 sigma_rel_FWHM_cb = fwhm / 2 / np.sqrt(2 * np.log(2))
0241
0242 cutoff_x = loc_cb - fwhm
0243 fraction_below = np.sum(hist.values()[hist.axes[0].centers < cutoff_x]) / np.sum(hist.values())
0244
0245 return sigma_rel_FWHM_cb, fraction_below
0246
0247 sigma_rel_FWHM_cb, fraction_below = summarize_fit(par)
0248 sigmas_rel_FWHM_cb.setdefault(PLOT_TITLE, {})[energy] = sigma_rel_FWHM_cb
0249 fractions_below.setdefault(PLOT_TITLE, {})[energy] = fraction_below
0250
0251 plt.legend()
0252 plt.xlabel("$E/p$", loc="right")
0253 plt.ylabel("Event yield", loc="top")
0254
0255 fig.savefig(output_dir / f"resolution_plots.pdf", bbox_inches="tight")
0256 fig.savefig(output_dir / f"resolution_plots.png", bbox_inches="tight")
0257 plt.show()
0258 plt.close(fig)
0259
0260 plt.figure()
0261 energy_values = np.array([float(energy.replace("GeV", "").replace("MeV", "e-3")) for energy in energies])
0262
0263 for clf_label, sigma_rel_FWHM_cb in sigmas_rel_FWHM_cb.items():
0264 sigma_over_e = np.array([sigma_rel_FWHM_cb[energy] for energy in energies]) * 100 # convert to %
0265
0266 def f(energy, stochastic, constant):
0267 return np.sqrt((stochastic / np.sqrt(energy)) ** 2 + constant ** 2)
0268 cond = energy_values >= 0.5
0269 try:
0270 import scipy.optimize
0271 par, pcov = scipy.optimize.curve_fit(f, energy_values[cond], sigma_over_e[cond], maxfev=10000)
0272 except RuntimeError:
0273 print("energy_values", energy_values[cond])
0274 print("sigma_over_e", sigma_over_e[cond])
0275 raise
0276 stochastic, constant = par
0277
0278 plt.plot(
0279 energy_values,
0280 sigma_over_e,
0281 marker=".",
0282 ls="none",
0283 label=f"{clf_label}"
0284 )
0285 xmin = np.min(energy_values[cond])
0286 xmax = np.max(energy_values[cond])
0287 xs = np.arange(xmin, xmax, 0.1)
0288 plt.plot(
0289 xs,
0290 f(xs, *par),
0291 ls="--",
0292 lw=0.5,
0293 label=fr"Functional fit: ${np.ceil(stochastic * 10) / 10:.1f}\% / \sqrt{{E}} \oplus {np.ceil(constant * 10) / 10:.1f}\%$",
0294 )
0295 xmin = np.min(energy_values)
0296 xmax = np.max(energy_values) * 1.05
0297 xs = np.arange(xmin, xmax + 0.1, 0.1)
0298 plt.fill_between(
0299 xs,
0300 np.sqrt((2 / np.sqrt(xs)) ** 2 + 1 ** 2),
0301 np.sqrt((3 / np.sqrt(xs)) ** 2 + 2 ** 2),
0302 lw=0., alpha=0.2, color="black", label=r"YR requirement $(2-3)\% / \sqrt{E} \oplus (1-2)\%$",
0303 )
0304 plt.xlim(0., xmax)
0305 plt.ylim(top=6.)
0306 plt.legend()
0307 plt.xlabel("Energy, GeV", loc="right")
0308 plt.ylabel(r"$\sigma_{E} / E$ derived from FWHM, \%", loc="top")
0309 plt.savefig(output_dir / f"resolution.pdf", bbox_inches="tight")
0310 plt.savefig(output_dir / f"resolution.png", bbox_inches="tight")
0311 plt.show()
0312 #+end_src
0313
0314 ** Pion rejection
0315
0316 #+begin_src jupyter-python
0317 fig, axs = plt.subplots(2, 4, sharex=True, sharey=True, figsize=(15, 6))
0318 fig_log, axs_log = plt.subplots(2, 4, sharex=True, sharey=True, figsize=(15, 6))
0319 fig_roc, axs_roc = plt.subplots(2, 4, sharex=True, sharey=True, figsize=(15, 6))
0320
0321 axs = np.ravel(np.array(axs))
0322 axs_log = np.ravel(np.array(axs_log))
0323 axs_roc = np.ravel(np.array(axs_roc))
0324
0325 rocs = {}
0326
0327 for ix, energy in enumerate(energies):
0328 energy_value = float(energy.replace("GeV", "").replace("MeV", "e-3"))
0329 def clf(events):
0330 return ak.drop_none(ak.max(events["EcalEndcapNClusters.energy"], axis=-1)) / energy_value
0331 e_pred = clf(e_eval[energy])
0332 pi_pred = clf(pi_eval[energy])
0333
0334 for do_log, ax in [(False, axs[ix]), (True, axs_log[ix])]:
0335 plt.sca(ax)
0336 patch = plt.stairs(norm_by_sum(e_over_p_hist[("e-", energy)].values()), e_over_p_hist[("e-", energy)].axes[0].edges, label=rf"$e^-$ {PLOT_TITLE}")
0337 patch = plt.stairs(norm_by_sum(e_over_p_hist[("pi-", energy)].values()), e_over_p_hist[("pi-", energy)].axes[0].edges, label=rf"$\pi^-$ {PLOT_TITLE}")
0338 plt.title(f"{energy}")
0339 plt.legend()
0340 plt.xlabel("Classifier output")
0341 plt.ylabel("Event yield")
0342 if do_log:
0343 plt.yscale("log")
0344
0345 plt.sca(axs_roc[ix])
0346 fpr = np.cumsum(e_over_p_hist[("pi-", energy)].values()[::-1])
0347 fpr /= fpr[-1]
0348 tpr = np.cumsum(e_over_p_hist[("e-", energy)].values()[::-1])
0349 tpr /= tpr[-1]
0350 cond = fpr != 0 # avoid infinite rejection (region of large uncertainty)
0351 cond &= tpr != 1 # avoid linear interpolation (region of large uncertainty)
0352 def mk_interp(tpr, fpr):
0353 def interp(eff):
0354 return np.interp(eff, tpr, fpr)
0355 return interp
0356 rocs.setdefault(PLOT_TITLE, {})[energy] = mk_interp(tpr, fpr)
0357 plt.plot(tpr[cond] * 100, 1 / fpr[cond], label=f"{PLOT_TITLE}")
0358 plt.yscale("log")
0359 plt.title(f"{energy}")
0360 plt.legend(loc="lower left")
0361 plt.xlabel("Electron efficiency, %")
0362 plt.ylabel("Pion rejection factor")
0363
0364 fig.savefig(output_dir / f"pred.pdf", bbox_inches="tight")
0365 fig.savefig(output_dir / f"pred.png", bbox_inches="tight")
0366 plt.close(fig)
0367 fig_log.savefig(output_dir / f"pred_log.pdf", bbox_inches="tight")
0368 fig_log.savefig(output_dir / f"pred_log.png", bbox_inches="tight")
0369 fig_log.show()
0370 fig_roc.savefig(output_dir / f"roc.pdf", bbox_inches="tight")
0371 fig_roc.savefig(output_dir / f"roc.png", bbox_inches="tight")
0372 fig_roc.show()
0373
0374 plt.figure()
0375 for clf_label, roc in rocs.items():
0376 plt.plot(
0377 [float(energy.replace("GeV", "").replace("MeV", "e-3")) for energy in energies],
0378 [1 / roc[energy](0.95) for energy in energies],
0379 marker=".",
0380 label=f"{clf_label}",
0381 )
0382 xmax = np.max(energy_values) * 1.05
0383 plt.xlim(0., xmax)
0384 plt.yscale("log")
0385 plt.legend()
0386 plt.xlabel("Energy, GeV")
0387 plt.ylabel(r"Pion rejection at 95\%")
0388 plt.savefig(output_dir / f"pion_rej.pdf", bbox_inches="tight")
0389 plt.savefig(output_dir / f"pion_rej.png", bbox_inches="tight")
0390 plt.show()
0391 #+end_src
