From 64c6bd1cefee8ce836637d47fe4571d7c2381274 Mon Sep 17 00:00:00 2001 From: Lars Bogner Date: Thu, 25 Jun 2026 13:36:54 +0200 Subject: [PATCH] Add photon edep export scripts and per-step presentation plots Co-Authored-By: Claude Opus 4.8 --- analysis/export_photon_edep_by_process.py | 65 +++++++++++++++++++++++ analysis/export_photon_edep_ev.py | 57 ++++++++++++++++++++ analysis/export_presentation_plots.py | 6 +++ 3 files changed, 128 insertions(+) create mode 100644 analysis/export_photon_edep_by_process.py create mode 100644 analysis/export_photon_edep_ev.py diff --git a/analysis/export_photon_edep_by_process.py b/analysis/export_photon_edep_by_process.py new file mode 100644 index 0000000..b83cf90 --- /dev/null +++ b/analysis/export_photon_edep_by_process.py @@ -0,0 +1,65 @@ +"""Truth-only histogram of deposited photon (pdg=22) energy in eV, grouped by +the Geant4 physics process that ended the step. Step-type histograms, one per +process, in both linear and log energy space. + +Data is the raw miniCaloSim steps parquet (truth), loaded lazily with polars — +only the photon rows and the (process, edep) columns are materialized. +""" + +from pathlib import Path + +import matplotlib.pyplot as plt +import numpy as np +import polars as pl + +FILE = "/home/lars/Programming/minicalo-data-exploration/pbwo4_10000events_hits.parquet" +OUT = Path("/home/lars/knowledge-base/meta/attachments") + +MEV_TO_EV = 1e6 + +# Lazy load: keep only photon steps and the two columns we need. +df = ( + pl.scan_parquet(FILE) + .filter(pl.col("pdg") == 22) + .select( + pl.col("process"), + (pl.col("edep") * MEV_TO_EV).alias("edep_ev"), + ) + .collect() +) +print("n photon rows:", df.height) + +# Process order by abundance, so the legend is stable and the busiest on top. +processes = ( + df.group_by("process").len().sort("len", descending=True).get_column("process").to_list() +) +series = {p: df.filter(pl.col("process") == p).get_column("edep_ev").to_numpy() for p in processes} + +all_edep = df.get_column("edep_ev").to_numpy() +lin_bins = np.linspace(0, np.quantile(all_edep, 0.999), 80) +pos = all_edep[all_edep > 0] +log_bins = np.geomspace(max(pos.min(), 1e-3), pos.max(), 80) + + +def plot(bins, xscale: str, fname: str) -> None: + fig, ax = plt.subplots(figsize=(7, 4.5)) + lo, hi = bins[0], bins[-1] + for p in processes: + v = series[p] + # Skip processes with nothing inside the bin range (e.g. all-zero edep + # processes on the log axis), which would add phantom legend entries. + if not np.any((v >= lo) & (v <= hi)): + continue + ax.hist(v, bins=bins, histtype="step", density=True, label=f"{p} (n={v.size})") + ax.set_xscale(xscale) + ax.set_yscale("log") + ax.set_xlabel("deposited energy [eV] (photons, pdg=22)") + ax.set_ylabel("density") + ax.legend(fontsize=8, title="process") + fig.tight_layout() + fig.savefig(OUT / fname, dpi=150, bbox_inches="tight") + print("saved", fname) + + +plot(lin_bins, "linear", "giant-photon-edep-by-process-ev-lin.png") +plot(log_bins, "log", "giant-photon-edep-by-process-ev-log.png") diff --git a/analysis/export_photon_edep_ev.py b/analysis/export_photon_edep_ev.py new file mode 100644 index 0000000..3cbb6d9 --- /dev/null +++ b/analysis/export_photon_edep_ev.py @@ -0,0 +1,57 @@ +"""Recreate of the photon (pdg=22) edep histogram in transformed (local-frame) +coordinates, with the deposited-energy axis in electron volts. + +Same data path as export_photon_edep.py, but edep is converted from the raw +MeV units to eV (x1e6) so the small-deposition photon spike is readable. +""" + +from pathlib import Path + +import matplotlib.pyplot as plt +import numpy as np + +import giant.analysis as a + +FILE = "/home/lars/Programming/giant/pbwo4_10k_9_predicted_local.parquet" +OUT = Path("/home/lars/knowledge-base/meta/attachments") + +MEV_TO_EV = 1e6 + +samples = a.load_predicted_local(FILE, sample_frac=0.15) +mask = samples.pdg == 22 +edep_idx = a.RAW_TARGET_NAMES.index("edep") +real = samples.real_raw[mask, edep_idx] * MEV_TO_EV # MeV -> eV +gen = samples.gen_raw[mask, edep_idx] * MEV_TO_EV +print("n photon rows:", mask.sum()) +print("real [eV]: mean", real.mean(), "max", real.max(), "frac==0", (real == 0).mean()) +print("gen [eV]: mean", gen.mean(), "max", gen.max(), "frac==0", (gen == 0).mean()) + +fig, ax = plt.subplots(figsize=(6, 4)) +bins = np.linspace(0, np.quantile(real, 0.999), 80) +ax.hist(real, bins=bins, alpha=0.6, label="real", density=True) +ax.hist(gen, bins=bins, alpha=0.6, label="gen", density=True) +ax.set_yscale("log") +ax.set_xlabel("deposited energy [eV] (photons, pdg=22)") +ax.set_ylabel("density") +ax.legend() + +fig.tight_layout() +fig.savefig(OUT / "giant-h1024n8d0.1lr3e-4-photon-edep-ev.png", dpi=150, bbox_inches="tight") +print("saved photon-edep-ev") + +# Second version: log-spaced energy axis to expose the low-deposition structure. +fig, ax = plt.subplots(figsize=(6, 4)) +lo = max(min(real[real > 0].min(), gen[gen > 0].min()), 1.0) +hi = max(real.max(), gen.max()) +log_bins = np.geomspace(lo, hi, 80) +ax.hist(real, bins=log_bins, alpha=0.6, label="real", density=True) +ax.hist(gen, bins=log_bins, alpha=0.6, label="gen", density=True) +ax.set_xscale("log") +ax.set_yscale("log") +ax.set_xlabel("deposited energy [eV] (photons, pdg=22)") +ax.set_ylabel("density") +ax.legend() + +fig.tight_layout() +fig.savefig(OUT / "giant-h1024n8d0.1lr3e-4-photon-edep-ev-logx.png", dpi=150, bbox_inches="tight") +print("saved photon-edep-ev-logx") diff --git a/analysis/export_presentation_plots.py b/analysis/export_presentation_plots.py index a8a08dd..c149e0c 100644 --- a/analysis/export_presentation_plots.py +++ b/analysis/export_presentation_plots.py @@ -110,6 +110,12 @@ savefig(fig, "event-total-energy") fig = a.plot_total_length(obs) savefig(fig, "event-total-length") +fig = a.plot_mean_energy_per_step(obs) +savefig(fig, "event-mean-energy-per-step") + +fig = a.plot_mean_length_per_step(obs) +savefig(fig, "event-mean-length-per-step") + fig = a.plot_longitudinal_profile(obs) savefig(fig, "event-longitudinal-profile")