Add photon edep export scripts and per-step presentation plots

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>

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")

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")

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")