Warning, /detector_benchmarks/benchmarks/backwards_ecal/backwards_ecal_dis.org is written in an unsupported language. File is not indexed.
0001 #+TITLE: ePIC EEEMCal benchmark in DIS events
0002 #+AUTHOR: detector_benchmarks contributors
0003 #+OPTIONS: d:t
0004
0005 This analysis measures the response of the backward electromagnetic
0006 calorimeter to the scattered electron in full DIS events. Unlike the
0007 single-particle benchmark, the highest-energy cluster in the event is not
0008 necessarily produced by the electron. We therefore select the scattered
0009 electron from =MCScatteredElectrons= and match the reconstructed cluster
0010 nearest to its endpoint on the EEEMCal surface.
0011
0012 The analysis is deliberately independent of how the input sample was
0013 produced. The benchmark workflow runs it on five reconstructed 10x100 GeV
0014 neutral-current DIS files with minimum $Q^2=1$ GeV$^2$. It reuses the DIS
0015 simulation defined by =tracking_performances_dis= and performs an EEEMCal-
0016 specific reconstruction that writes the collections needed below.
0017
0018 The same program can later be run on signal-only and beam-background-overlay
0019 samples. Comparison plots can then be made from the two saved =summary.json=
0020 and =distributions.npz= files without repeating the event processing.
0021
0022 Run the wired benchmark with:
0023
0024 #+begin_example
0025 snakemake --cores 1 backwards_ecal_dis_run_locally
0026 #+end_example
0027
0028 To analyze another already-reconstructed sample manually, generate the Python
0029 program and provide a file list:
0030
0031 #+begin_example
0032 snakemake benchmarks/backwards_ecal/backwards_ecal_dis.org2py.py
0033 env \
0034 INPUT_FILE_LIST=/path/to/dis_files.list \
0035 OUTPUT_DIR=results/backwards_ecal_dis/my_sample \
0036 SAMPLE_LABEL="DIS, 10x100 GeV" \
0037 python benchmarks/backwards_ecal/backwards_ecal_dis.org2py.py
0038 #+end_example
0039
0040 The file list must contain one EDM4eic ROOT file per line. =MATCH_RADIUS_MM=
0041 and =ENERGY_BIN_EDGES= may be overridden in the environment. The default
0042 50 mm matching radius is intentionally recorded in the output and should be
0043 reviewed using the matching-distance plot before it becomes a fixed benchmark
0044 requirement.
0045
0046 * Setup
0047
0048 #+begin_src jupyter-python
0049 import json
0050 import os
0051 from pathlib import Path
0052
0053 import awkward as ak
0054 import matplotlib
0055 matplotlib.use("Agg")
0056 import matplotlib.pyplot as plt
0057 import numpy as np
0058 import scipy.optimize
0059 import scipy.stats
0060 import uproot
0061
0062
0063 INPUT_FILE_LIST = Path(os.environ["INPUT_FILE_LIST"])
0064 OUTPUT_DIR = Path(os.environ.get("OUTPUT_DIR", "results/backwards_ecal_dis/manual"))
0065 SAMPLE_LABEL = os.environ.get("SAMPLE_LABEL", "DIS")
0066 MATCH_RADIUS_MM = float(os.environ.get("MATCH_RADIUS_MM", "50"))
0067 ENERGY_BIN_EDGES = np.asarray([
0068 float(value)
0069 for value in os.environ.get(
0070 "ENERGY_BIN_EDGES", "0.5,1,2,3,4,5,6,8,10,15,20"
0071 ).split(",")
0072 ])
0073
0074 ETA_MIN = float(os.environ.get("ETA_MIN", "-3.0"))
0075 ETA_MAX = float(os.environ.get("ETA_MAX", "-1.8"))
0076 ENDPOINT_Z_MIN_MM = float(os.environ.get("ENDPOINT_Z_MIN_MM", "-1900"))
0077 ENDPOINT_Z_MAX_MM = float(os.environ.get("ENDPOINT_Z_MAX_MM", "-1700"))
0078
0079 OUTPUT_DIR.mkdir(parents=True, exist_ok=True)
0080
0081
0082 def read_file_list(path):
0083 with path.open() as stream:
0084 files = [
0085 line.strip()
0086 for line in stream
0087 if line.strip() and not line.lstrip().startswith("#")
0088 ]
0089 if not files:
0090 raise RuntimeError(f"No input files found in {path}")
0091 return files
0092
0093
0094 files = read_file_list(INPUT_FILE_LIST)
0095 branches = [
0096 "MCScatteredElectrons_objIdx.index",
0097 "MCParticles.momentum.*",
0098 "MCParticles.endpoint.*",
0099 "EcalEndcapNClusters.energy",
0100 "EcalEndcapNClusters.position.*",
0101 ]
0102 #+end_src
0103
0104 * Object selection
0105
0106 The =MCScatteredElectrons= collection is a subset collection, represented in
0107 the ROOT file by indices into =MCParticles=. We retain candidates in the
0108 backward acceptance and, if more than one remains, use the one with the most
0109 negative longitudinal momentum. This reproduces the selection used in the
0110 original EEEMCal DIS study.
0111
0112 The endpoint requirement ensures that the selected electron actually reaches
0113 the EEEMCal. The closest reconstructed cluster is accepted only inside the
0114 configured matching radius. Events without such a cluster remain in the
0115 efficiency denominator but not in the response or residual distributions.
0116
0117 #+begin_src jupyter-python
0118 def reduce_length_one(array):
0119 """Reduce a zero-or-one element jagged array without using ak.firsts."""
0120 return ak.max(array, axis=-1, mask_identity=True)
0121
0122
0123 def select_scattered_electron(events):
0124 indices = events["MCScatteredElectrons_objIdx.index"]
0125 px = events["MCParticles.momentum.x"][indices]
0126 py = events["MCParticles.momentum.y"][indices]
0127 pz = events["MCParticles.momentum.z"][indices]
0128
0129 momentum = np.sqrt(px * px + py * py + pz * pz)
0130 eta = np.arcsinh(pz / np.sqrt(px * px + py * py))
0131 accepted = (eta > ETA_MIN) & (eta < ETA_MAX)
0132 accepted_pz = ak.mask(pz, accepted)
0133 best = ak.argmin(
0134 accepted_pz, axis=-1, keepdims=True, mask_identity=True
0135 )
0136 selected_indices = indices[best]
0137
0138 def mc_value(name):
0139 return reduce_length_one(events[name][selected_indices])
0140
0141 return {
0142 "energy": reduce_length_one(momentum[best]),
0143 "eta": reduce_length_one(eta[best]),
0144 "endpoint_x": mc_value("MCParticles.endpoint.x"),
0145 "endpoint_y": mc_value("MCParticles.endpoint.y"),
0146 "endpoint_z": mc_value("MCParticles.endpoint.z"),
0147 }
0148
0149
0150 def match_cluster(events, electron):
0151 cluster_x = events["EcalEndcapNClusters.position.x"]
0152 cluster_y = events["EcalEndcapNClusters.position.y"]
0153 distance = np.sqrt(
0154 (cluster_x - electron["endpoint_x"]) ** 2
0155 + (cluster_y - electron["endpoint_y"]) ** 2
0156 )
0157 closest = ak.argmin(distance, axis=-1, keepdims=True, mask_identity=True)
0158
0159 return {
0160 "energy": reduce_length_one(events["EcalEndcapNClusters.energy"][closest]),
0161 "x": reduce_length_one(cluster_x[closest]),
0162 "y": reduce_length_one(cluster_y[closest]),
0163 "distance": reduce_length_one(distance[closest]),
0164 }
0165 #+end_src
0166
0167 * Event processing
0168
0169 Only one scalar record per selected electron is retained, so processing is
0170 performed in chunks and does not require all EDM data to fit in memory.
0171
0172 #+begin_src jupyter-python
0173 selected_chunks = []
0174 matched_chunks = []
0175 total_events = 0
0176
0177 sources = {path: "events" for path in files}
0178 for events in uproot.iterate(
0179 sources,
0180 filter_name=branches,
0181 step_size="100 MB",
0182 library="ak",
0183 ):
0184 total_events += len(events)
0185 electron = select_scattered_electron(events)
0186 cluster = match_cluster(events, electron)
0187
0188 endpoint_accepted = (
0189 (electron["endpoint_z"] > ENDPOINT_Z_MIN_MM)
0190 & (electron["endpoint_z"] < ENDPOINT_Z_MAX_MM)
0191 )
0192 matched = endpoint_accepted & (cluster["distance"] < MATCH_RADIUS_MM)
0193
0194 selected_chunks.append({
0195 key: ak.to_numpy(ak.drop_none(ak.mask(value, endpoint_accepted)))
0196 for key, value in electron.items()
0197 })
0198 matched_chunks.append({
0199 "truth_energy": ak.to_numpy(
0200 ak.drop_none(ak.mask(electron["energy"], matched))
0201 ),
0202 "truth_eta": ak.to_numpy(
0203 ak.drop_none(ak.mask(electron["eta"], matched))
0204 ),
0205 "response": ak.to_numpy(
0206 ak.drop_none(ak.mask(cluster["energy"] / electron["energy"], matched))
0207 ),
0208 "dx": ak.to_numpy(
0209 ak.drop_none(
0210 ak.mask(cluster["x"] - electron["endpoint_x"], matched)
0211 )
0212 ),
0213 "dy": ak.to_numpy(
0214 ak.drop_none(
0215 ak.mask(cluster["y"] - electron["endpoint_y"], matched)
0216 )
0217 ),
0218 "distance": ak.to_numpy(
0219 ak.drop_none(ak.mask(cluster["distance"], endpoint_accepted))
0220 ),
0221 })
0222
0223
0224 def concatenate(chunks, key):
0225 arrays = [chunk[key] for chunk in chunks if len(chunk[key])]
0226 return np.concatenate(arrays) if arrays else np.asarray([], dtype=float)
0227
0228
0229 truth_energy_selected = concatenate(selected_chunks, "energy")
0230 truth_eta_selected = concatenate(selected_chunks, "eta")
0231 truth_energy_matched = concatenate(matched_chunks, "truth_energy")
0232 truth_eta_matched = concatenate(matched_chunks, "truth_eta")
0233 response = concatenate(matched_chunks, "response")
0234 dx = concatenate(matched_chunks, "dx")
0235 dy = concatenate(matched_chunks, "dy")
0236 closest_distance = concatenate(matched_chunks, "distance")
0237
0238 np.savez_compressed(
0239 OUTPUT_DIR / "distributions.npz",
0240 truth_energy_selected=truth_energy_selected,
0241 truth_eta_selected=truth_eta_selected,
0242 truth_energy_matched=truth_energy_matched,
0243 truth_eta_matched=truth_eta_matched,
0244 response=response,
0245 dx=dx,
0246 dy=dy,
0247 closest_distance=closest_distance,
0248 energy_bin_edges=ENERGY_BIN_EDGES,
0249 )
0250 #+end_src
0251
0252 * Resolution summaries
0253
0254 Energy resolution follows the single-particle benchmark convention: fit the
0255 =E/p= distribution with a Crystal Ball function and convert its FWHM to an
0256 equivalent Gaussian sigma. Position resolution is the half-width of the
0257 central 68% interval, as in the position-resolution addition to that
0258 benchmark. A bin with insufficient entries or a failed fit is reported as
0259 =null= in JSON rather than aborting the full analysis.
0260
0261 #+begin_src jupyter-python
0262 def crystal_ball_resolution(values):
0263 values = values[np.isfinite(values)]
0264 if len(values) < 30:
0265 return np.nan
0266
0267 counts, edges = np.histogram(values, bins=110, range=(0.0, 1.10))
0268 centers = 0.5 * (edges[:-1] + edges[1:])
0269
0270 def model(x, normalization, beta, m, location, scale):
0271 return normalization * scipy.stats.crystalball.pdf(
0272 x, beta, m, loc=location, scale=scale
0273 )
0274
0275 fit_slice = slice(5, None)
0276 location0 = centers[fit_slice][np.argmax(counts[fit_slice])]
0277 parameters0 = (max(np.sum(counts[10:]) * 0.02, 1.0), 2.0, 3.0, location0, 0.02)
0278 try:
0279 parameters, _ = scipy.optimize.curve_fit(
0280 model,
0281 centers[fit_slice],
0282 counts[fit_slice],
0283 p0=parameters0,
0284 bounds=([0.0, 0.05, 1.01, 0.1, 0.001], [np.inf, 20.0, 100.0, 1.1, 0.5]),
0285 maxfev=20000,
0286 )
0287 location = parameters[3]
0288 grid = np.linspace(0.0, 1.1, 10000)
0289 curve = model(grid, *parameters)
0290 above_half_maximum = grid[curve >= 0.5 * np.max(curve)]
0291 if location <= 0.0 or len(above_half_maximum) < 2:
0292 return np.nan
0293 fwhm = above_half_maximum[-1] - above_half_maximum[0]
0294 return fwhm / (2.0 * np.sqrt(2.0 * np.log(2.0))) / location
0295 except (RuntimeError, ValueError):
0296 return np.nan
0297
0298
0299 def central_68_half_width(values):
0300 values = values[np.isfinite(values)]
0301 if len(values) < 20:
0302 return np.nan
0303 low, high = np.quantile(values, [0.16, 0.84])
0304 return 0.5 * (high - low)
0305
0306
0307 def optional_number(value):
0308 return float(value) if np.isfinite(value) else None
0309
0310
0311 bin_centers = 0.5 * (ENERGY_BIN_EDGES[:-1] + ENERGY_BIN_EDGES[1:])
0312 energy_resolution = []
0313 position_resolution_x = []
0314 position_resolution_y = []
0315 matching_efficiency = []
0316 selected_counts = []
0317 matched_counts = []
0318
0319 for low, high in zip(ENERGY_BIN_EDGES[:-1], ENERGY_BIN_EDGES[1:]):
0320 selected_in_bin = (truth_energy_selected >= low) & (truth_energy_selected < high)
0321 matched_in_bin = (truth_energy_matched >= low) & (truth_energy_matched < high)
0322 denominator = int(np.count_nonzero(selected_in_bin))
0323 numerator = int(np.count_nonzero(matched_in_bin))
0324
0325 selected_counts.append(denominator)
0326 matched_counts.append(numerator)
0327 matching_efficiency.append(numerator / denominator if denominator else np.nan)
0328 energy_resolution.append(crystal_ball_resolution(response[matched_in_bin]))
0329 position_resolution_x.append(central_68_half_width(dx[matched_in_bin]))
0330 position_resolution_y.append(central_68_half_width(dy[matched_in_bin]))
0331
0332 summary = {
0333 "sample_label": SAMPLE_LABEL,
0334 "input_file_list": str(INPUT_FILE_LIST),
0335 "number_of_files": len(files),
0336 "number_of_events": total_events,
0337 "number_of_endpoint_selected_electrons": len(truth_energy_selected),
0338 "number_of_matched_electrons": len(truth_energy_matched),
0339 "selection": {
0340 "eta": [ETA_MIN, ETA_MAX],
0341 "endpoint_z_mm": [ENDPOINT_Z_MIN_MM, ENDPOINT_Z_MAX_MM],
0342 "match_radius_mm": MATCH_RADIUS_MM,
0343 },
0344 "energy_bin_edges_GeV": ENERGY_BIN_EDGES.tolist(),
0345 "selected_counts": selected_counts,
0346 "matched_counts": matched_counts,
0347 "matching_efficiency": [optional_number(x) for x in matching_efficiency],
0348 "energy_resolution": [optional_number(x) for x in energy_resolution],
0349 "position_resolution_x_mm": [optional_number(x) for x in position_resolution_x],
0350 "position_resolution_y_mm": [optional_number(x) for x in position_resolution_y],
0351 }
0352
0353 with (OUTPUT_DIR / "summary.json").open("w") as stream:
0354 json.dump(summary, stream, indent=2)
0355 #+end_src
0356
0357 * Plots
0358
0359 #+begin_src jupyter-python
0360 fig, axes = plt.subplots(1, 3, figsize=(12, 3.5))
0361
0362 axes[0].hist(response, bins=110, range=(0.0, 1.1), histtype="step")
0363 axes[0].set_xlabel(r"$E_{cluster}/p_{truth}$")
0364 axes[0].set_ylabel("Events")
0365
0366 axes[1].hist(dx, bins=200, range=(-50.0, 50.0), histtype="step", label="x")
0367 axes[1].hist(dy, bins=200, range=(-50.0, 50.0), histtype="step", label="y")
0368 axes[1].set_xlabel("Cluster - endpoint [mm]")
0369 axes[1].set_ylabel("Events")
0370 axes[1].legend()
0371
0372 axes[2].hist(closest_distance, bins=200, range=(0.0, 200.0), histtype="step")
0373 axes[2].axvline(MATCH_RADIUS_MM, color="black", linestyle="--", label="match cut")
0374 axes[2].set_xlabel("Closest cluster distance [mm]")
0375 axes[2].set_ylabel("Endpoint-selected events")
0376 axes[2].legend()
0377
0378 fig.suptitle(SAMPLE_LABEL)
0379 fig.tight_layout()
0380 fig.savefig(OUTPUT_DIR / "distributions.png", dpi=150)
0381 fig.savefig(OUTPUT_DIR / "distributions.pdf")
0382 plt.close(fig)
0383
0384 fig, axes = plt.subplots(1, 3, figsize=(12, 3.5))
0385
0386 axes[0].plot(bin_centers, 100.0 * np.asarray(energy_resolution), marker="o")
0387 axes[0].set_xlabel("Truth electron energy [GeV]")
0388 axes[0].set_ylabel(r"$\sigma_E/E$ from FWHM [%]")
0389
0390 axes[1].plot(bin_centers, position_resolution_x, marker="o", label="x")
0391 axes[1].plot(bin_centers, position_resolution_y, marker="o", label="y")
0392 axes[1].set_xlabel("Truth electron energy [GeV]")
0393 axes[1].set_ylabel("Central 68% half-width [mm]")
0394 axes[1].legend()
0395
0396 axes[2].plot(bin_centers, matching_efficiency, marker="o")
0397 axes[2].set_xlabel("Truth electron energy [GeV]")
0398 axes[2].set_ylabel("Cluster matching efficiency")
0399 axes[2].set_ylim(0.0, 1.05)
0400
0401 fig.suptitle(SAMPLE_LABEL)
0402 fig.tight_layout()
0403 fig.savefig(OUTPUT_DIR / "performance.png", dpi=150)
0404 fig.savefig(OUTPUT_DIR / "performance.pdf")
0405 plt.close(fig)
0406
0407 print(json.dumps(summary, indent=2))
0408 #+end_src
