Warning, /detector_benchmarks/benchmarks/backgrounds/ecal_backwards.org is written in an unsupported language. File is not indexed.
0001 #+PROPERTY: header-args:jupyter-python :session /jpy:localhost#8888:eeemcal :async yes :results drawer :exports both
0002
0003 #+TITLE: ePIC EEEMCal background rates 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 from pyHepMC3 import HepMC3
0019 #+end_src
0020
0021 #+begin_src jupyter-python :results silent
0022 from dask.distributed import Client
0023
0024 client = Client(os.environ.get("DASK_SCHEDULER", "localhost:8786"))
0025 #+end_src
0026
0027 * Plotting setup
0028
0029 #+begin_src jupyter-python :results silent
0030 import matplotlib as mpl
0031 import matplotlib.pyplot as plt
0032
0033 def setup_presentation_style():
0034 mpl.rcParams.update(mpl.rcParamsDefault)
0035 plt.style.use('ggplot')
0036 plt.rcParams.update({
0037 'axes.labelsize': 8,
0038 'axes.titlesize': 9,
0039 'figure.titlesize': 9,
0040 'figure.figsize': (4, 3),
0041 'legend.fontsize': 7,
0042 'legend.loc': 'upper right',
0043 'xtick.labelsize': 8,
0044 'ytick.labelsize': 8,
0045 'xaxis.labellocation': 'right',
0046 'yaxis.labellocation': 'top',
0047 'pgf.rcfonts': False,
0048 })
0049
0050 setup_presentation_style()
0051 #+end_src
0052
0053 * Analysis
0054
0055 ** Input files
0056
0057 #+begin_src jupyter-python :results silent
0058 ELECTRON_BEAM_GAS_GEN=os.environ.get("ELECTRON_BEAM_GAS_GEN", "../beam_gas_ep_10GeV_foam_emin10keV_10Mevt_vtx.hepmc")
0059 ELECTRON_BEAM_GAS_SIM=os.environ.get("ELECTRON_BEAM_GAS_SIM", "electron_beam_gas.edm4hep.root")
0060 PHYSICS_PROCESS_SIM=os.environ.get("PHYSICS_PROCESS_SIM", "pythia8NCDIS_10x100_minQ2=1_beamEffects_xAngle=-0.025_hiDiv_1.edm4hep.root")
0061 PROTON_BEAM_GAS_GEN=os.environ.get("PROTON_BEAM_GAS_GEN", "100GeV.hepmc")
0062 PROTON_BEAM_GAS_SIM=os.environ.get("PROTON_BEAM_GAS_SIM", "proton+beam_gas_ep.edm4hep.root")
0063
0064 output_dir=Path(os.environ.get("OUTPUT_DIR", "./"))
0065 output_dir.mkdir(parents=True, exist_ok=True)
0066 #+end_src
0067
0068 #+begin_src jupyter-python :results silent
0069 builder = ak.ArrayBuilder()
0070 n = 0
0071 f = HepMC3.ReaderPlugin(PROTON_BEAM_GAS_GEN, "libHepMC3rootIO.so", "newReaderRootTreefile")
0072 event = HepMC3.GenEvent()
0073 while f.read_event(event):
0074 builder.begin_list()
0075 assert event.length_unit().name == "MM"
0076 for vertex in event.vertices():
0077 with builder.record("Vector4D"):
0078 builder.field("x")
0079 builder.real(vertex.position().x())
0080 builder.field("y")
0081 builder.real(vertex.position().y())
0082 builder.field("z")
0083 builder.real(vertex.position().z())
0084 builder.field("t")
0085 builder.real(vertex.position().t())
0086 builder.end_list()
0087 n += 1
0088 if n > 10000: break
0089
0090 vertices_proton_beam_gas = builder.snapshot()
0091
0092 builder = ak.ArrayBuilder()
0093 n = 0
0094 f = HepMC3.ReaderPlugin(ELECTRON_BEAM_GAS_GEN, "libHepMC3rootIO.so", "newReaderRootTreefile")
0095 event = HepMC3.GenEvent()
0096 while f.read_event(event):
0097 builder.begin_list()
0098 assert event.length_unit().name == "MM"
0099 for vertex in event.vertices():
0100 with builder.record("Vector4D"):
0101 builder.field("x")
0102 builder.real(vertex.position().x())
0103 builder.field("y")
0104 builder.real(vertex.position().y())
0105 builder.field("z")
0106 builder.real(vertex.position().z())
0107 builder.field("t")
0108 builder.real(vertex.position().t())
0109 builder.end_list()
0110 n += 1
0111 if n > 10000: break
0112
0113 vertices_electron_beam_gas = builder.snapshot()
0114 #+end_src
0115
0116 #+begin_src jupyter-python :results silent
0117 def filter_name(name):
0118 return "Hits." in name or "MCParticles." in name
0119
0120 datasets = {
0121 # https://wiki.bnl.gov/EPIC/index.php?title=Electron_Beam_Gas
0122 "electron beam gas 10 GeV": {
0123 "vertices": vertices_electron_beam_gas,
0124 "events": uproot.dask({ELECTRON_BEAM_GAS_SIM: "events"}, filter_name=filter_name, open_files=False, steps_per_file=32),
0125 "cmap": "cool",
0126 "rate": 3.2e6, # Hz
0127 },
0128 "DIS 10x100, $Q^2 > 1$ GeV${}^2$": {
0129 "events": uproot.dask({PHYSICS_PROCESS_SIM: "events"}, filter_name=filter_name, open_files=False, steps_per_file=32),
0130 "rate": 184e3, # Hz
0131 },
0132 # https://wiki.bnl.gov/EPIC/index.php?title=Hadron_Beam_Gas
0133 "proton beam gas 100 GeV": {
0134 "vertices": vertices_proton_beam_gas,
0135 "events": uproot.dask({PROTON_BEAM_GAS_SIM: "events"}, filter_name=filter_name, open_files=False, steps_per_file=32),
0136 "cmap": "summer",
0137 "rate": 31.9e3, # Hz
0138 },
0139 }
0140
0141 for ds in datasets.values():
0142 ds["events"].eager_compute_divisions()
0143 #+end_src
0144
0145 ** Vertex distributions
0146
0147 #+begin_src jupyter-python
0148 for label, ds in datasets.items():
0149 if "vertices" not in ds: continue
0150 vs = ds["vertices"]
0151 weight = ds["rate"] / ak.num(ds["vertices"], axis=0)
0152 plt.hist(ak.ravel(vs.t[:,0]), bins=100, histtype="step", label=label)
0153 plt.minorticks_on()
0154 plt.xlabel("vertex[0].t, mm")
0155 plt.legend()
0156 plt.savefig(output_dir / "vertex_time_distribution.png", bbox_inches="tight")
0157 plt.show()
0158 plt.clf()
0159
0160 for label, ds in datasets.items():
0161 if "vertices" not in ds: continue
0162 vs = ds["vertices"]
0163 weight = ds["rate"] / ak.num(ds["vertices"], axis=0)
0164 plt.hist(ak.ravel(vs.z[:,0]), bins=100, histtype="step", label=label)
0165 plt.minorticks_on()
0166 plt.xlabel("vertex[0].z, mm")
0167 plt.legend()
0168 plt.savefig(output_dir / "vertex_z_distribution.png", bbox_inches="tight")
0169 plt.show()
0170 plt.clf()
0171
0172 for label, ds in datasets.items():
0173 if "vertices" not in ds: continue
0174 vs = ds["vertices"]
0175 cmap = ds["cmap"]
0176 weight = ds["rate"] / ak.num(ds["vertices"], axis=0)
0177 plt.hist2d(vs.z[:,0].to_numpy(), vs.x[:,0].to_numpy(), bins=(100, np.linspace(-130, 130, 160)), cmin=1e-30, label=label, cmap=cmap)
0178 plt.plot([], color=mpl.colormaps[cmap](0.5), label=label)
0179 plt.minorticks_on()
0180 plt.xlabel("vertex[0].z, mm")
0181 plt.ylabel("vertex[0].x, mm")
0182 plt.legend()
0183 plt.savefig(output_dir / "vertex_xz_distribution.png", bbox_inches="tight")
0184 plt.show()
0185 plt.clf()
0186
0187 for ix, (label, ds) in enumerate(datasets.items()):
0188 if "vertices" not in ds: continue
0189 vs = ds["vertices"]
0190 cmap = ds["cmap"]
0191 weight = ds["rate"] / ak.num(ds["vertices"], axis=0)
0192 plt.hist2d(vs.z[:,0].to_numpy(), vs.y[:,0].to_numpy(), bins=(100, 100), cmin=1e-30, cmap=cmap)
0193 plt.colorbar()
0194 plt.minorticks_on()
0195 plt.xlabel("vertex[0].z, mm")
0196 plt.ylabel("vertex[0].y, mm")
0197 plt.title(label)
0198 plt.savefig(output_dir / f"vertex_yz_distribution_{ix}.png", bbox_inches="tight")
0199 plt.show()
0200 plt.clf()
0201 #+end_src
0202
0203 ** Simulation results
0204
0205 #+begin_src jupyter-python
0206 for collection_name in ["EcalEndcapNHits", "EcalEndcapPHits"]:
0207 for dataset_ix, (label, ds) in enumerate(datasets.items()):
0208 events = ds["events"]
0209
0210 energy_sums = ak.sum(events[f"{collection_name}.energy"].head(10000), axis=1)
0211 event_id = ak.argmax(energy_sums)
0212 xs = events[f"{collection_name}.position.x"].head(event_id + 1)[event_id].to_numpy()
0213 ys = events[f"{collection_name}.position.y"].head(event_id + 1)[event_id].to_numpy()
0214
0215 bin_widths = [None, None]
0216 for ix, vals in enumerate([xs, ys]):
0217 centers = np.unique(vals)
0218 diffs = centers[1:] - centers[:-1]
0219 EPSILON = 1e-5
0220 bin_widths[ix] = np.min(diffs[diffs > EPSILON]) if np.sum(diffs > EPSILON) > 0 else 1.
0221 print(f"bin_widths[{ix}]", bin_widths[ix])
0222
0223 bins = {
0224 "EcalEndcapNHits": [np.arange(-750., 750., bin_width) for bin_width in bin_widths],
0225 "EcalEndcapPHits": [np.arange(-1800., 1800., bin_width) for bin_width in bin_widths],
0226 }[collection_name]
0227
0228 plt.hist2d(
0229 xs,
0230 ys,
0231 weights=events[f"{collection_name}.energy"].head(event_id + 1)[event_id].to_numpy(),
0232 bins=bins,
0233 cmin=1e-10,
0234 )
0235 plt.colorbar().set_label("energy, GeV", loc="top")
0236 plt.title(f"{label}, event_id={event_id}\n{collection_name}")
0237 plt.xlabel("hit x, mm")
0238 plt.ylabel("hit y, mm")
0239 plt.savefig(output_dir / f"{collection_name}_event_display_{dataset_ix}.png", bbox_inches="tight")
0240 plt.show()
0241 plt.clf()
0242 #+end_src
0243
0244 ** Discovering number of cells
0245
0246 Using HyperLogLog algorithm would be faster here, or actually load
0247 DD4hep geometry and count sensitive volumes.
0248
0249 #+begin_src jupyter-python
0250 def unique(array):
0251 if ak.backend(array) == "typetracer":
0252 ak.typetracer.touch_data(array)
0253 return array
0254 return ak.from_numpy(np.unique(ak.to_numpy(ak.ravel(array))))
0255 unique_delayed = dask.delayed(unique)
0256 len_delayed = dask.delayed(len)
0257
0258 cellID_for_r = dict()
0259
0260 for collection_name in ["EcalEndcapNHits", "EcalEndcapPHits"]:
0261 r_axis = {
0262 "EcalEndcapNHits": bh.axis.Regular(75, 0., 750.),
0263 "EcalEndcapPHits": bh.axis.Regular(90, 0., 1800.),
0264 }[collection_name]
0265 ds = datasets["DIS 10x100, $Q^2 > 1$ GeV${}^2$"]
0266 events = ds["events"]
0267
0268 r = np.hypot(
0269 ak.ravel(events[f"{collection_name}.position.x"]),
0270 ak.ravel(events[f"{collection_name}.position.y"]),
0271 )
0272 cellID = ak.ravel(events[f"{collection_name}.cellID"])
0273
0274 cellID_for_r[collection_name] = np.array(client.gather(client.compute([
0275 len_delayed(unique_delayed(
0276 cellID[(r >= r_min) & (r < r_max)].map_partitions(unique)
0277 ))
0278 for r_min, r_max in zip(r_axis.edges[:-1], r_axis.edges[1:])
0279 ])))
0280
0281 print(cellID_for_r[collection_name])
0282 print(sum(cellID_for_r[collection_name]))
0283
0284 plt.stairs(
0285 cellID_for_r[collection_name],
0286 r_axis.edges,
0287 )
0288
0289 plt.title(f"{collection_name}")
0290 plt.legend()
0291 plt.xlabel("r, mm")
0292 dr = (r_axis.edges[1] - r_axis.edges[0])
0293 plt.ylabel(f"Number of towers per {dr} mm slice in $r$")
0294 plt.savefig(output_dir / f"{collection_name}_num_towers.png", bbox_inches="tight")
0295 plt.show()
0296 plt.clf()
0297 #+end_src
0298
0299 ** Plotting the rates
0300
0301 #+begin_src jupyter-python
0302 for collection_name in ["EcalEndcapNHits", "EcalEndcapPHits"]:
0303 r_axis = {
0304 "EcalEndcapNHits": bh.axis.Regular(75, 0., 750.),
0305 "EcalEndcapPHits": bh.axis.Regular(90, 0., 1800.),
0306 }[collection_name]
0307 for edep_min in [0.005, 0.015, 0.050]: # GeV
0308 for label, ds in datasets.items():
0309 events = ds["events"]
0310 weight = ds["rate"] / len(events)
0311
0312 r = np.hypot(
0313 ak.ravel(events[f"{collection_name}.position.x"]),
0314 ak.ravel(events[f"{collection_name}.position.y"]),
0315 )
0316 edep = ak.ravel(events[f"{collection_name}.energy"])
0317 r = r[edep > edep_min]
0318
0319 hist = dh.factory(
0320 r,
0321 axes=(r_axis,),
0322 ).compute()
0323 plt.stairs(
0324 hist.values() * weight / cellID_for_r[collection_name],
0325 hist.axes[0].edges,
0326 label=f"{label}",
0327 )
0328
0329 plt.title(f"for $E_{{dep.}} >$ {edep_min * 1000} MeV\n{collection_name}")
0330 plt.legend()
0331 plt.xlabel("r, mm")
0332 plt.ylabel("rate per tower, Hz")
0333 plt.yscale("log")
0334 plt.savefig(output_dir / f"{collection_name}_hit_rate_vs_r_edep_min_{edep_min:.3f}.png", bbox_inches="tight")
0335 plt.show()
0336 plt.clf()
0337 #+end_src
0338
0339 #+begin_src jupyter-python
0340 for collection_name in ["EcalEndcapNHits", "EcalEndcapPHits"]:
0341 for totedep_min in [-1, 0, 0.1, 0.5, 1.0, 5.0, 10.]: # GeV
0342 for label, ds in datasets.items():
0343 events = ds["events"]
0344 weight = ds["rate"] / len(events)
0345
0346 z = ds["events"]["MCParticles.vertex.z"][:,1]
0347 totedep = ak.sum(events[f"{collection_name}.energy"], axis=1)
0348 z = z[totedep > totedep_min]
0349
0350 hist = dh.factory(
0351 z,
0352 axes=(bh.axis.Regular(250, -7500., 17500.),),
0353 ).compute()
0354 plt.stairs(
0355 hist.values() * weight,
0356 hist.axes[0].edges,
0357 label=f"{label}",
0358 )
0359
0360 plt.title(rf"for events with $E_{{\mathrm{{dep. tot.}}}}$ $>$ {totedep_min} GeV" + f"\n{collection_name}")
0361 plt.legend()
0362 plt.xlabel("$z$ of the first interaction vertex, mm")
0363 plt.ylabel("rate, Hz")
0364 plt.yscale("log")
0365 plt.savefig(output_dir / f"{collection_name}_hit_rate_vs_z_totedep_min_{totedep_min:.1f}.png", bbox_inches="tight")
0366 plt.show()
0367 plt.clf()
0368 #+end_src
0369
0370 #+begin_src jupyter-python
0371 num_towers_cache = {
0372 "LumiSpecCAL": 200,
0373 "LumiDirectPCAL": 1,
0374 "ZDCHcal": 1470,
0375 "LFHCAL": 578338,
0376 "ZDC_WSi_": 187043,
0377 "EcalBarrelScFi": 124205483,
0378 "EcalEndcapP": 15037,
0379 "ZDCEcal": 400,
0380 "EcalEndcapPInsert": 536,
0381 "HcalEndcapPInsert": 20251,
0382 "B0ECal": 131,
0383 "HcalEndcapN": 13800,
0384 "HcalBarrel": 7680,
0385 "EcalBarrelImaging": 5765469,
0386 "EcalEndcapN": 2988,
0387 "ZDC_PbSi_": 44344,
0388 }
0389
0390 fig_cmb = plt.figure()
0391 ax_cmb = fig_cmb.gca()
0392
0393 for edep_min in [0]: # GeV
0394 for dataset_ix, (x_offset, (ds_label, ds)) in enumerate(zip(np.linspace(-0.3, 0.3, len(datasets)), datasets.items())):
0395 events = ds["events"]
0396 weight = ds["rate"] / len(events)
0397
0398 labels = []
0399 values = []
0400 norms = []
0401
0402 for branch_name in events.fields:
0403 if ".energy" not in branch_name: continue
0404 if "ZDC_SiliconPix_Hits" in branch_name: continue
0405
0406 edep = ak.ravel(events[branch_name])
0407
0408 #cellID = ak.ravel(events[branch_name.replace(".energy", ".cellID")])
0409 #num_towers = len(unique_delayed(
0410 # cellID.map_partitions(unique)
0411 #).compute())
0412
0413 key = branch_name.replace("Hits.energy", "")
0414 if key not in num_towers_cache:
0415 print(f"The \"{key}\" not in num_towers_cache. Skipping.")
0416 continue
0417 num_towers = num_towers_cache[key]
0418
0419 labels.append(branch_name.replace("Hits.energy", ""))
0420 values.append(ak.count(edep[edep > edep_min]))
0421 norms.append(num_towers if num_towers != 0 else np.nan)
0422
0423 fig_cur = plt.figure()
0424 ax_cur = fig_cur.gca()
0425
0426 values, = dask.compute(values)
0427 for ax, per_tower, offset, width in [
0428 (ax_cmb, True, x_offset, 2 * 0.3 / (len(datasets) - 1)),
0429 (ax_cur, False, 0, 2 * 0.3),
0430 ]:
0431 ax.bar(
0432 np.arange(len(labels)) + offset,
0433 weight * np.array(values) / (np.array(norms) if per_tower else np.ones_like(norms)),
0434 width=width,
0435 tick_label=labels,
0436 label=ds_label,
0437 color=f"C{dataset_ix}",
0438 )
0439
0440 plt.sca(ax_cur)
0441 plt.legend()
0442 plt.title(f"for $E_{{dep.}} >$ {edep_min * 1000} MeV")
0443 plt.ylabel("rate, Hz")
0444 plt.yscale("log")
0445 plt.xticks(rotation=90, ha='right')
0446 fig_cur.savefig(f"rates_edep_min_{edep_min}_{dataset_ix}.png", bbox_inches="tight")
0447 plt.show()
0448 plt.close(fig_cur)
0449
0450 plt.sca(ax_cmb)
0451 plt.legend()
0452 plt.title(f"for $E_{{dep.}} >$ {edep_min * 1000} MeV")
0453 plt.ylabel("rate per tower, Hz")
0454 plt.yscale("log")
0455 plt.xticks(rotation=90, ha='right')
0456 fig_cmb.savefig(f"rates_edep_min_{edep_min}.png", bbox_inches="tight")
0457 plt.show()
0458 plt.clf()
0459 #+end_src
