Warning, /detector_benchmarks/benchmarks/ecal_gaps/ecal_gaps.org is written in an unsupported language. File is not indexed.
0001 #+PROPERTY: header-args:jupyter-python :session /jpy:localhost#8888:gap :async yes :results drawer :exports both
0002
0003 #+TITLE: ePIC ECal gap study
0004 #+AUTHOR: Dmitry Kalinkin
0005 #+OPTIONS: d:t
0006
0007 #+begin_src jupyter-python :results silent
0008 import os
0009 from pathlib import Path
0010
0011 import awkward as ak
0012 import boost_histogram as bh
0013 import dask
0014 import dask_awkward as dak
0015 import dask_histogram as dh
0016 import numpy as np
0017 import uproot
0018 #+end_src
0019
0020 #+begin_src jupyter-python :results slient
0021 from dask.distributed import Client
0022 client = Client(os.environ.get("DASK_SCHEDULER", "localhost:8786"))
0023 #+end_src
0024
0025 * Plotting setup
0026
0027 #+begin_src jupyter-python :results silent
0028 import matplotlib as mpl
0029 import matplotlib.pyplot as plt
0030
0031 def setup_presentation_style():
0032 mpl.rcParams.update(mpl.rcParamsDefault)
0033 plt.style.use('ggplot')
0034 plt.rcParams.update({
0035 'axes.labelsize': 8,
0036 'axes.titlesize': 9,
0037 'figure.titlesize': 9,
0038 'figure.figsize': (4, 3),
0039 'legend.fontsize': 7,
0040 'xtick.labelsize': 8,
0041 'ytick.labelsize': 8,
0042 'xaxis.labellocation': 'right',
0043 'yaxis.labellocation': 'top',
0044 'pgf.rcfonts': False,
0045 })
0046
0047 setup_presentation_style()
0048 #+end_src
0049
0050 ** Settings
0051
0052 #+begin_src jupyter-python :results silent
0053 DETECTOR_CONFIG=os.environ.get("DETECTOR_CONFIG")
0054
0055 output_dir=Path(os.environ.get("OUTPUT_DIR", "./"))
0056 output_dir.mkdir(parents=True, exist_ok=True)
0057 #+end_src
0058
0059 * Analysis
0060
0061 #+begin_src jupyter-python :results silent
0062 def filter_name(name):
0063 return (
0064 "PARAMETERS" not in name
0065 and all(coll not in name for coll in ["_intMap", "_floatMap", "_stringMap", "_doubleMap"])
0066 )
0067
0068 import functools
0069
0070 axis_eta = bh.axis.Regular(200, -4, 4)
0071 axis_eta_coarse = bh.axis.Regular(100, -4, 4)
0072
0073 @functools.cache
0074 def get_events(particle="e-", energy="20GeV", num_files=1):
0075 events = uproot.dask(
0076 {}
0077 | {f"sim_output/ecal_gaps/{DETECTOR_CONFIG}/{particle}/{energy}/3to50deg/{particle}_{energy}_3to50deg.{INDEX:04d}.eicrecon.tree.edm4eic.root": "events" for INDEX in range(num_files)}
0078 | {f"sim_output/ecal_gaps/{DETECTOR_CONFIG}/{particle}/{energy}/45to135deg/{particle}_{energy}_45to135deg.{INDEX:04d}.eicrecon.tree.edm4eic.root": "events" for INDEX in range(num_files)}
0079 | {f"sim_output/ecal_gaps/{DETECTOR_CONFIG}/{particle}/{energy}/130to177deg/{particle}_{energy}_130to177deg.{INDEX:04d}.eicrecon.tree.edm4eic.root": "events" for INDEX in range(num_files)}
0080 ,
0081 filter_name=filter_name, open_files=False, steps_per_file=1,
0082 )
0083
0084 pt = np.hypot(events["MCParticles.momentum.x"][:,0], events["MCParticles.momentum.y"][:,0])
0085 theta = np.arctan2(pt, events["MCParticles.momentum.z"][:,0])
0086 eta = -np.log(np.tan(theta / 2))
0087 p = np.hypot(pt, events["MCParticles.momentum.z"][:,0])
0088
0089 hist_norm = dh.factory(
0090 eta,
0091 axes=(
0092 axis_eta,
0093 ),
0094 ).compute()
0095 weight_lut = 1 / hist_norm.values(flow=True)
0096
0097 def get_weight(array):
0098 if ak.backend(array) == "typetracer":
0099 ak.typetracer.touch_data(array)
0100 return array
0101 return ak.from_numpy(weight_lut[axis_eta.index(ak.to_numpy(array))])
0102
0103 weights = eta.map_partitions(get_weight)
0104
0105 events["p_thrown"] = p
0106 events["eta_thrown"] = eta
0107 events["weights"] = weights
0108
0109 return events
0110 #+end_src
0111
0112 #+begin_src jupyter-python
0113 particle = "e-"
0114
0115 for energy in ["500MeV", "5GeV", "20GeV"]:
0116 events = get_events(particle=particle, energy=energy)
0117
0118 for calo_name in ["EcalEndcapN", "EcalBarrelScFi", "EcalBarrelImaging", "EcalEndcapP"]:
0119 cond = ak.sum(events[f"{calo_name}RecHits.energy"], axis=1) > 10e-3 # GeV
0120
0121 hist, profile = client.gather(client.compute([
0122 dh.factory(
0123 events["eta_thrown"][cond],
0124 (ak.sum(events[f"{calo_name}RecHits.energy"], axis=1) / events["p_thrown"])[cond],
0125 axes=(
0126 axis_eta,
0127 bh.axis.Regular(100, 0.1, 1.1),
0128 ),
0129 weights=events["weights"][cond],
0130 ),
0131 dh.factory(
0132 events["eta_thrown"][cond],
0133 sample=(ak.sum(events[f"{calo_name}RecHits.energy"], axis=1) / events["p_thrown"])[cond],
0134 storage=bh.storage.WeightedMean(),
0135 axes=(
0136 axis_eta_coarse,
0137 ),
0138 weights=events["weights"][cond],
0139 ),
0140 ]))
0141
0142 cmap = plt.get_cmap(name="rainbow", lut=None).copy()
0143 cmap.set_under("none")
0144
0145 plt.pcolormesh(
0146 hist.axes[0].edges,
0147 hist.axes[1].edges,
0148 hist.values().T,
0149 cmap=cmap,
0150 norm=mpl.colors.LogNorm(
0151 vmin=np.min(hist.values()[hist.values() > 0]),
0152 ),
0153 )
0154 plt.colorbar()
0155 std = np.sqrt(profile.variances())
0156 cond = profile.values() > std
0157 plt.errorbar(profile.axes[0].centers[cond], profile.values()[cond], yerr=std[cond], marker=".", markersize=2, color="black", ls="none", lw=0.6, capsize=1.)
0158 plt.xlabel(r"$\eta_{thrown}$")
0159 plt.ylabel(r"$\sum E_{\mathrm{dep.}} / p_{\mathrm{thrown}}$")
0160 plt.title(f"{energy} {particle} in {calo_name}")
0161 plt.minorticks_on()
0162 plt.savefig(output_dir / f"ecal_gap_{particle}_{energy}_{calo_name}.png", bbox_inches="tight")
0163 plt.show()
0164 plt.clf()
0165 #+end_src
0166
0167 #+begin_src jupyter-python
0168 particle = "e-"
0169
0170 for energy in ["500MeV", "5GeV", "20GeV"]:
0171 events = get_events(particle=particle, energy=energy)
0172
0173 calos = ["EcalEndcapN", "EcalBarrelScFi", "EcalEndcapP"]
0174 total_energy = sum([
0175 ak.sum(events[f"{calo_name}RecHits.energy"], axis=1)
0176 for calo_name in calos
0177 ])
0178
0179 hist, profile = client.gather(client.compute([
0180 dh.factory(
0181 events["eta_thrown"],
0182 total_energy / events["p_thrown"],
0183 axes=(
0184 axis_eta,
0185 bh.axis.Regular(100, 0.1, 1.1),
0186 ),
0187 weights=events["weights"],
0188 ),
0189 dh.factory(
0190 events["eta_thrown"],
0191 sample=total_energy / events["p_thrown"],
0192 storage=bh.storage.WeightedMean(),
0193 axes=(
0194 axis_eta_coarse,
0195 ),
0196 weights=events["weights"],
0197 ),
0198 ]))
0199
0200 cmap = plt.get_cmap(name="rainbow", lut=None).copy()
0201 cmap.set_under("none")
0202
0203 plt.pcolormesh(
0204 hist.axes[0].edges,
0205 hist.axes[1].edges,
0206 hist.values().T,
0207 cmap=cmap,
0208 norm=mpl.colors.LogNorm(
0209 vmin=np.min(hist.values()[hist.values() > 0]),
0210 ),
0211 )
0212 plt.colorbar()
0213 std = np.sqrt(profile.variances())
0214 cond = profile.values() > std
0215 plt.errorbar(profile.axes[0].centers[cond], profile.values()[cond], yerr=std[cond], marker=".", markersize=2, color="black", ls="none", lw=0.6, capsize=1.)
0216 plt.xlabel(r"$\eta_{thrown}$")
0217 plt.ylabel(r"$\sum E_{\mathrm{dep.}} / p_{\mathrm{thrown}}$")
0218 plt.title(f"{energy} {particle}\n" + "+".join(calos))
0219 plt.minorticks_on()
0220 plt.savefig(output_dir / f"ecal_gap_{particle}_{energy}_sum_all.png", bbox_inches="tight")
0221 plt.show()
0222 plt.clf()
0223 #+end_src
