From e3c191d1fdae53b1ddb655bb97562ef85f70fdb8 Mon Sep 17 00:00:00 2001 From: Gabor Galgoczi Date: Tue, 31 Mar 2026 17:33:59 +0000 Subject: [PATCH 01/17] add photon path visualization script with wavelength coloring Plots GPU photon paths from record.npy colored by wavelength. Supports custom photon selection, sphere overlays, and wavelength colorbar. Useful for visualizing WLS conversion and Rayleigh scattering in optical simulations. --- optiphy/ana/plot_photon_paths.py | 161 +++++++++++++++++++++++++++++++ 1 file changed, 161 insertions(+) create mode 100644 optiphy/ana/plot_photon_paths.py diff --git a/optiphy/ana/plot_photon_paths.py b/optiphy/ana/plot_photon_paths.py new file mode 100644 index 000000000..bf6a09245 --- /dev/null +++ b/optiphy/ana/plot_photon_paths.py @@ -0,0 +1,161 @@ +#!/usr/bin/env python3 +"""Plot GPU photon paths colored by wavelength from record.npy. + +Usage: + python optiphy/ana/plot_photon_paths.py [photon_indices] [--output path.png] + +Examples: + # Plot first 10 hit photons + python optiphy/ana/plot_photon_paths.py /tmp/$USER/opticks/GEOM/GEOM/GPUPhotonSourceMinimal/ALL0_no_opticks_event_name/A000 + + # Plot specific photons by index + python optiphy/ana/plot_photon_paths.py /tmp/$USER/opticks/.../A000 2,19,6 + + # Custom output path + python optiphy/ana/plot_photon_paths.py /tmp/$USER/opticks/.../A000 2,19,6 --output my_plot.png +""" +import argparse +import numpy as np +import matplotlib +matplotlib.use('Agg') +import matplotlib.pyplot as plt +from mpl_toolkits.mplot3d import Axes3D +import matplotlib.colors as mcolors +from matplotlib.cm import ScalarMappable + + +def wl_to_rgb(wl): + """Convert wavelength (nm) to RGB tuple. Covers 300-780nm.""" + r = g = b = 0.0 + if 300 <= wl < 380: + t = (wl - 300) / (380 - 300) + r = 0.4 * (1 - t) + 0.5 * t + g = 0 + b = 0.4 * (1 - t) + 1.0 * t + elif 380 <= wl < 440: + r = -(wl - 440) / (440 - 380); g = 0; b = 1 + elif 440 <= wl < 490: + r = 0; g = (wl - 440) / (490 - 440); b = 1 + elif 490 <= wl < 510: + r = 0; g = 1; b = -(wl - 510) / (510 - 490) + elif 510 <= wl < 580: + r = (wl - 510) / (580 - 510); g = 1; b = 0 + elif 580 <= wl < 645: + r = 1; g = -(wl - 645) / (645 - 580); b = 0 + elif 645 <= wl <= 780: + r = 1; g = 0; b = 0 + else: + r = g = b = 0.3 + return (max(0, min(1, r)), max(0, min(1, g)), max(0, min(1, b))) + + +def get_steps(record, pidx): + """Return number of valid steps for photon pidx.""" + rec_p = record[pidx] + rf = rec_p.reshape(rec_p.shape[0], -1) + return int(np.sum(np.any(rf != 0, axis=1))) + + +def plot_photon_paths(event_dir, photon_indices=None, output="photon_paths.png", + sphere_radii=None, title=None, lim=None): + record = np.load(f"{event_dir}/record.npy") + photon = np.load(f"{event_dir}/photon.npy") + + q3 = photon[:, 3, :].copy().view(np.uint32) + flags = q3[:, 0] & 0xFFFF + hit_idx = np.where(flags == 0x40)[0] + + if photon_indices is None: + photon_indices = hit_idx[:10] + + fig = plt.figure(figsize=(12, 10)) + ax = fig.add_subplot(111, projection='3d') + + wl_min, wl_max = 800, 300 + for pidx in photon_indices: + ns = get_steps(record, pidx) + if ns < 2: + continue + rec_p = record[pidx] + x = rec_p[:ns, 0, 0] + y = rec_p[:ns, 0, 1] + z = rec_p[:ns, 0, 2] + wl = rec_p[:ns, 2, 3] + + wl_min = min(wl_min, wl.min()) + wl_max = max(wl_max, wl.max()) + + for s in range(ns - 1): + color = wl_to_rgb(float(wl[s])) + ax.plot([x[s], x[s + 1]], [y[s], y[s + 1]], [z[s], z[s + 1]], + color=color, alpha=0.9, linewidth=2.5) + + ax.scatter(x[0], y[0], z[0], c=[wl_to_rgb(float(wl[0]))], s=60, + marker='o', edgecolors='black', linewidths=0.8, zorder=5) + ax.scatter(x[-1], y[-1], z[-1], c='red', s=100, marker='*', zorder=5) + + # Draw spheres if requested + if sphere_radii: + u = np.linspace(0, 2 * np.pi, 60) + v = np.linspace(0, np.pi, 30) + sphere_colors = ['mediumpurple', 'lightgreen', 'lightyellow', 'lightcoral'] + sphere_alphas = [0.1, 0.05, 0.05, 0.05] + for i, r in enumerate(sphere_radii): + xs = r * np.outer(np.cos(u), np.sin(v)) + ys = r * np.outer(np.sin(u), np.sin(v)) + zs = r * np.outer(np.ones_like(u), np.cos(v)) + ci = min(i, len(sphere_colors) - 1) + ax.plot_surface(xs, ys, zs, alpha=sphere_alphas[ci], color=sphere_colors[ci]) + + # Wavelength colorbar + wl_range = np.linspace(wl_min, wl_max, 256) + colors = [wl_to_rgb(w) for w in wl_range] + cmap = mcolors.ListedColormap(colors) + norm = mcolors.Normalize(vmin=wl_min, vmax=wl_max) + sm = ScalarMappable(cmap=cmap, norm=norm) + sm.set_array([]) + plt.colorbar(sm, ax=ax, shrink=0.5, pad=0.08, label='Wavelength (nm)') + + ax.set_xlabel('X (mm)') + ax.set_ylabel('Y (mm)') + ax.set_zlabel('Z (mm)') + if title: + ax.set_title(title) + if lim: + ax.set_xlim(-lim, lim) + ax.set_ylim(-lim, lim) + ax.set_zlim(-lim, lim) + ax.view_init(elev=20, azim=135) + plt.tight_layout() + plt.savefig(output, dpi=180) + print(f"Saved {output}") + + +def main(): + parser = argparse.ArgumentParser(description=__doc__, + formatter_class=argparse.RawDescriptionHelpFormatter) + parser.add_argument("event_dir", help="Path to opticks event folder containing record.npy") + parser.add_argument("indices", nargs='?', default=None, + help="Comma-separated photon indices (default: first 10 hits)") + parser.add_argument("--output", "-o", default="photon_paths.png", help="Output image path") + parser.add_argument("--spheres", default=None, + help="Comma-separated sphere radii to draw (e.g. 10,30)") + parser.add_argument("--title", "-t", default=None, help="Plot title") + parser.add_argument("--lim", type=float, default=None, + help="Axis limit in mm (symmetric)") + args = parser.parse_args() + + indices = None + if args.indices: + indices = [int(x) for x in args.indices.split(',')] + + spheres = None + if args.spheres: + spheres = [float(x) for x in args.spheres.split(',')] + + plot_photon_paths(args.event_dir, indices, args.output, + sphere_radii=spheres, title=args.title, lim=args.lim) + + +if __name__ == "__main__": + main() From 0a52f1f38fae83ce258ea3d43ea50cab4f4ada2a Mon Sep 17 00:00:00 2001 From: Gabor Galgoczi Date: Wed, 1 Apr 2026 12:47:36 +0000 Subject: [PATCH 02/17] add GPU vs G4 comparison script with simulation runner and plots Runs GPURaytrace, collects GPU and G4 hits, and generates 6 comparison plots: hit count, shifted wavelength, full wavelength, bulk arrival time, full arrival time, and 3D hit positions. All with sqrt(N) error bars. Can also run on pre-existing .npy files without re-simulating. --- optiphy/ana/run_and_compare.py | 231 +++++++++++++++++++++++++++++++++ 1 file changed, 231 insertions(+) create mode 100755 optiphy/ana/run_and_compare.py diff --git a/optiphy/ana/run_and_compare.py b/optiphy/ana/run_and_compare.py new file mode 100755 index 000000000..461d01548 --- /dev/null +++ b/optiphy/ana/run_and_compare.py @@ -0,0 +1,231 @@ +#!/usr/bin/env python3 +"""Run GPU and G4 simulations and compare hit distributions. + +Runs GPURaytrace with a given GDML and config, then plots: + 1. Hit count with sqrt(N) error bars + 2. WLS-shifted wavelength distribution + 3. Full wavelength distribution + 4. Arrival time (bulk, truncated) + 5. Arrival time (full range, no overflow) + 6. 3D hit position scatter for GPU and G4 + +Usage: + python optiphy/ana/run_and_compare.py -g det.gdml -s 42 [--outdir plots] + + # Skip simulation, use existing .npy files: + python optiphy/ana/run_and_compare.py --gpu-hits gpu_hits.npy --g4-hits g4_hits.npy +""" +import argparse +import os +import subprocess +import sys +import math + +import numpy as np +import matplotlib +matplotlib.use('Agg') +import matplotlib.pyplot as plt +from mpl_toolkits.mplot3d import Axes3D + + +def run_simulation(gdml, config, macro, seed): + """Run GPURaytrace and return (gpu_hits_path, g4_hits_path, gpu_nhits, g4_nhits).""" + cmd = ["/opt/eic-opticks/bin/GPURaytrace", + "-g", gdml, "-m", macro, "-s", str(seed)] + if config: + cmd += ["-c", config] + + print(f"Running: {' '.join(cmd)}") + result = subprocess.run(cmd, capture_output=True, text=True, timeout=600) + output = result.stdout + result.stderr + + gpu_nhits = g4_nhits = 0 + for line in output.split('\n'): + if 'Opticks: NumHits:' in line: + gpu_nhits = int(line.strip().split()[-1]) + if 'Geant4: NumHits:' in line: + g4_nhits = int(line.strip().split()[-1]) + + print(f" GPU: {gpu_nhits} hits, G4: {g4_nhits} hits") + return "gpu_hits.npy", "g4_hits.npy", gpu_nhits, g4_nhits + + +def load_hits(path, expected_cols=None): + """Load hit array and reshape to (N, ?, 4).""" + a = np.load(path) + if a.ndim == 2: + ncols = a.shape[1] // 4 + a = a.reshape(-1, ncols, 4) + return a + + +def plot_with_errors(ax, data1, data2, bins, label1, label2, xlabel): + """Plot two histograms as points with sqrt(N) error bars.""" + h1, edges = np.histogram(data1, bins=bins) + h2, _ = np.histogram(data2, bins=bins) + centers = (edges[:-1] + edges[1:]) / 2 + width = (edges[1] - edges[0]) * 0.35 + ax.errorbar(centers - width / 2, h1, yerr=np.sqrt(np.maximum(h1, 1)), + fmt='o', color='dodgerblue', markersize=4, capsize=2, + linewidth=1, label=label1) + ax.errorbar(centers + width / 2, h2, yerr=np.sqrt(np.maximum(h2, 1)), + fmt='s', color='orangered', markersize=4, capsize=2, + linewidth=1, label=label2) + ax.set_xlabel(xlabel) + ax.set_ylabel('Counts') + ax.legend() + + +def make_plots(gpu, g4, outdir, title_extra=""): + os.makedirs(outdir, exist_ok=True) + + gpu_wl = gpu[:, 2, 3] + g4_wl = g4[:, 2, 3] + gpu_t = gpu[:, 0, 3] + g4_t = g4[:, 0, 3] + gpu_pos = gpu[:, 0, :3] + g4_pos = g4[:, 0, :3] + + diff = 100 * (len(gpu) / len(g4) - 1) if len(g4) > 0 else 0 + z_score = (len(gpu) - len(g4)) / math.sqrt(len(gpu) + len(g4)) if (len(gpu) + len(g4)) > 0 else 0 + header = f"GPU={len(gpu)} G4={len(g4)} ({diff:+.1f}%, {z_score:+.1f}σ)" + if title_extra: + header = f"{title_extra}\n{header}" + + # 1. Hit count + fig, ax = plt.subplots(figsize=(6, 5)) + vals = [len(gpu), len(g4)] + errs = [math.sqrt(v) for v in vals] + ax.errorbar([0], [vals[0]], yerr=[errs[0]], fmt='o', markersize=12, + capsize=8, linewidth=2, color='dodgerblue', label='GPU') + ax.errorbar([1], [vals[1]], yerr=[errs[1]], fmt='s', markersize=12, + capsize=8, linewidth=2, color='orangered', label='G4') + ax.set_xticks([0, 1]) + ax.set_xticklabels(['GPU', 'G4']) + ax.set_ylabel('Hits') + ax.set_title(f'Hit Count\n{header}') + ax.set_xlim(-0.5, 1.5) + ax.legend() + for i, (v, e) in enumerate(zip(vals, errs)): + ax.text(i + 0.15, v, f'{v}±{e:.0f}', va='center', fontsize=10) + plt.tight_layout() + plt.savefig(f'{outdir}/hits.png', dpi=150) + plt.close() + + # 2. WLS-shifted wavelength + fig, ax = plt.subplots(figsize=(8, 5)) + gpu_s = gpu_wl[gpu_wl > 380] + g4_s = g4_wl[g4_wl > 380] + plot_with_errors(ax, gpu_s, g4_s, np.arange(380, 550, 15), + f'GPU ({len(gpu_s)})', f'G4 ({len(g4_s)})', + 'Wavelength (nm)') + ax.set_title(f'WLS-shifted Wavelength (>380nm)\n{header}') + plt.tight_layout() + plt.savefig(f'{outdir}/wavelength_shifted.png', dpi=150) + plt.close() + + # 3. Full wavelength + fig, ax = plt.subplots(figsize=(8, 5)) + plot_with_errors(ax, gpu_wl, g4_wl, np.arange(330, 550, 15), + f'GPU ({len(gpu)})', f'G4 ({len(g4)})', + 'Wavelength (nm)') + ax.set_title(f'Full Wavelength\n{header}') + plt.tight_layout() + plt.savefig(f'{outdir}/wavelength_full.png', dpi=150) + plt.close() + + # 4. Arrival time (bulk, truncated at 99th percentile) + fig, ax = plt.subplots(figsize=(8, 5)) + t_cut = max(np.percentile(gpu_t, 99), np.percentile(g4_t, 99)) + t_bins = np.linspace(0, t_cut, 30) + gpu_over = (gpu_t > t_cut).sum() + g4_over = (g4_t > t_cut).sum() + plot_with_errors(ax, gpu_t[gpu_t <= t_cut], g4_t[g4_t <= t_cut], t_bins, + f'GPU (overflow={gpu_over})', f'G4 (overflow={g4_over})', + 'Time (ns)') + ax.set_title(f'Arrival Time (t < {t_cut:.0f}ns)\n{header}') + plt.tight_layout() + plt.savefig(f'{outdir}/time_bulk.png', dpi=150) + plt.close() + + # 5. Arrival time (full range, no overflow) + fig, ax = plt.subplots(figsize=(8, 5)) + t_max = max(gpu_t.max(), g4_t.max()) * 1.05 + t_bins_full = np.linspace(0, t_max, 50) + plot_with_errors(ax, gpu_t, g4_t, t_bins_full, + f'GPU ({len(gpu)})', f'G4 ({len(g4)})', + 'Time (ns)') + ax.set_title(f'Arrival Time (full range)\n{header}') + plt.tight_layout() + plt.savefig(f'{outdir}/time_full.png', dpi=150) + plt.close() + + # 6. 3D hit positions + fig = plt.figure(figsize=(14, 6)) + + ax1 = fig.add_subplot(121, projection='3d') + ax1.scatter(gpu_pos[:, 0], gpu_pos[:, 1], gpu_pos[:, 2], + c='dodgerblue', s=3, alpha=0.5) + ax1.set_xlabel('X (mm)') + ax1.set_ylabel('Y (mm)') + ax1.set_zlabel('Z (mm)') + ax1.set_title(f'GPU hit positions ({len(gpu)})') + + ax2 = fig.add_subplot(122, projection='3d') + ax2.scatter(g4_pos[:, 0], g4_pos[:, 1], g4_pos[:, 2], + c='orangered', s=3, alpha=0.5) + ax2.set_xlabel('X (mm)') + ax2.set_ylabel('Y (mm)') + ax2.set_zlabel('Z (mm)') + ax2.set_title(f'G4 hit positions ({len(g4)})') + + plt.suptitle(f'3D Hit Positions\n{header}') + plt.tight_layout() + plt.savefig(f'{outdir}/positions_3d.png', dpi=150) + plt.close() + + # Print summary + print(f"\nSummary: {header}") + print(f" Wavelength: GPU mean={gpu_wl.mean():.1f}nm G4 mean={g4_wl.mean():.1f}nm") + print(f" Time: GPU mean={gpu_t.mean():.2f}ns G4 mean={g4_t.mean():.2f}ns") + print(f" WLS shifted: GPU {100*(gpu_wl>380).mean():.1f}% G4 {100*(g4_wl>380).mean():.1f}%") + print(f"\nPlots saved to {outdir}/:") + for f in ['hits.png', 'wavelength_shifted.png', 'wavelength_full.png', + 'time_bulk.png', 'time_full.png', 'positions_3d.png']: + print(f" {f}") + + +def main(): + parser = argparse.ArgumentParser(description=__doc__, + formatter_class=argparse.RawDescriptionHelpFormatter) + parser.add_argument("-g", "--gdml", default="det.gdml", help="GDML geometry file") + parser.add_argument("-c", "--config", default=None, help="Config name (e.g. det_debug)") + parser.add_argument("-m", "--macro", default="tests/run_genstep.mac", help="G4 macro file") + parser.add_argument("-s", "--seed", type=int, default=42, help="Random seed") + parser.add_argument("--outdir", default="plots", help="Output directory for plots") + parser.add_argument("--title", default="", help="Extra title text") + parser.add_argument("--gpu-hits", default=None, help="Skip sim, use existing GPU hits .npy") + parser.add_argument("--g4-hits", default=None, help="Skip sim, use existing G4 hits .npy") + args = parser.parse_args() + + if args.gpu_hits and args.g4_hits: + print(f"Using existing files: {args.gpu_hits}, {args.g4_hits}") + else: + run_simulation(args.gdml, args.config, args.macro, args.seed) + args.gpu_hits = "gpu_hits.npy" + args.g4_hits = "g4_hits.npy" + + gpu = load_hits(args.gpu_hits) + g4 = load_hits(args.g4_hits) + + # Normalize to (N, 4, 4) — take first 4 rows if more + if gpu.shape[1] > 4: + gpu = gpu[:, :4, :] + if g4.shape[1] > 4: + g4 = g4[:, :4, :] + + make_plots(gpu, g4, args.outdir, title_extra=args.title) + + +if __name__ == "__main__": + main() From 38f2e263735d7c18dde21cd0f5ca1d3a10c23cc3 Mon Sep 17 00:00:00 2001 From: Gabor Galgoczi Date: Wed, 1 Apr 2026 15:07:54 +0000 Subject: [PATCH 03/17] add step count distribution plot to GPU vs G4 comparison Shows number of optical boundary/surface steps each detected photon takes before reaching the SiPM. GPU steps from record.npy, G4 steps from extended hit array (when available). Uses sqrt(N) error bars. --- optiphy/ana/run_and_compare.py | 91 +++++++++++++++++++++++++++++++--- 1 file changed, 83 insertions(+), 8 deletions(-) diff --git a/optiphy/ana/run_and_compare.py b/optiphy/ana/run_and_compare.py index 461d01548..8ceaf220e 100755 --- a/optiphy/ana/run_and_compare.py +++ b/optiphy/ana/run_and_compare.py @@ -76,7 +76,7 @@ def plot_with_errors(ax, data1, data2, bins, label1, label2, xlabel): ax.legend() -def make_plots(gpu, g4, outdir, title_extra=""): +def make_plots(gpu, g4, outdir, title_extra="", g4_full=None, g4_raw_shape=None): os.makedirs(outdir, exist_ok=True) gpu_wl = gpu[:, 2, 3] @@ -184,14 +184,86 @@ def make_plots(gpu, g4, outdir, title_extra=""): plt.savefig(f'{outdir}/positions_3d.png', dpi=150) plt.close() + # 7. Step count distribution (GPU from record.npy, G4 from extended hit array) + gpu_steps = None + g4_steps = None + + # GPU: try to load record.npy and count steps per hit + record_path = os.path.join(os.path.dirname(os.environ.get('OPTICKS_EVTDIR', '')), + 'record.npy') + # Try common paths + for rpath in [record_path, + '/tmp/MISSING_USER/opticks/GEOM/GEOM/GPURaytrace/ALL0_no_opticks_event_name/A000/record.npy', + '/tmp/MISSING_USER/opticks/GEOM/GEOM/GPUPhotonSourceMinimal/ALL0_no_opticks_event_name/A000/record.npy']: + if os.path.exists(rpath): + try: + record = np.load(rpath) + photon_all = np.load(rpath.replace('record.npy', 'photon.npy')) + q3_all = photon_all[:, 3, :].copy().view(np.uint32) + flags_all = q3_all[:, 0] & 0xFFFF + hit_mask = (flags_all == 0x40) + hit_indices = np.where(hit_mask)[0] + rec_flat = record[hit_indices].reshape(len(hit_indices), record.shape[1], -1) + gpu_steps = np.sum(np.any(rec_flat != 0, axis=2), axis=1) + print(f" GPU step counts loaded from {rpath}") + except Exception: + pass + break + + # G4: check if extended hit array has step count in row 3, col 3 + if g4_raw_shape is not None and g4_raw_shape[1] >= 5: + g4_steps = g4_full[:, 3, 3].astype(int) + print(f" G4 step counts loaded from extended hit array") + + if gpu_steps is not None or g4_steps is not None: + fig, ax = plt.subplots(figsize=(8, 5)) + s_max = 0 + if gpu_steps is not None: + s_max = max(s_max, np.percentile(gpu_steps, 99)) + if g4_steps is not None: + s_max = max(s_max, np.percentile(g4_steps, 99)) + s_bins = np.linspace(0, s_max * 1.1, 30) + + if gpu_steps is not None and g4_steps is not None: + plot_with_errors(ax, gpu_steps, g4_steps, s_bins, + f'GPU ({len(gpu_steps)})', f'G4 ({len(g4_steps)})', + 'Steps to detection') + elif gpu_steps is not None: + h, edges = np.histogram(gpu_steps, bins=s_bins) + centers = (edges[:-1] + edges[1:]) / 2 + ax.errorbar(centers, h, yerr=np.sqrt(np.maximum(h, 1)), + fmt='o', color='dodgerblue', markersize=4, capsize=2, + label=f'GPU ({len(gpu_steps)})') + ax.set_xlabel('Steps to detection') + ax.set_ylabel('Counts') + ax.legend() + elif g4_steps is not None: + h, edges = np.histogram(g4_steps, bins=s_bins) + centers = (edges[:-1] + edges[1:]) / 2 + ax.errorbar(centers, h, yerr=np.sqrt(np.maximum(h, 1)), + fmt='s', color='orangered', markersize=4, capsize=2, + label=f'G4 ({len(g4_steps)})') + ax.set_xlabel('Steps to detection') + ax.set_ylabel('Counts') + ax.legend() + + ax.set_title(f'Steps to Detection\n{header}') + plt.tight_layout() + plt.savefig(f'{outdir}/step_count.png', dpi=150) + plt.close() + # Print summary print(f"\nSummary: {header}") print(f" Wavelength: GPU mean={gpu_wl.mean():.1f}nm G4 mean={g4_wl.mean():.1f}nm") print(f" Time: GPU mean={gpu_t.mean():.2f}ns G4 mean={g4_t.mean():.2f}ns") print(f" WLS shifted: GPU {100*(gpu_wl>380).mean():.1f}% G4 {100*(g4_wl>380).mean():.1f}%") + if gpu_steps is not None: + print(f" GPU steps: mean={gpu_steps.mean():.0f} median={np.median(gpu_steps):.0f} max={gpu_steps.max()}") + if g4_steps is not None: + print(f" G4 steps: mean={g4_steps.mean():.0f} median={np.median(g4_steps):.0f} max={g4_steps.max()}") print(f"\nPlots saved to {outdir}/:") for f in ['hits.png', 'wavelength_shifted.png', 'wavelength_full.png', - 'time_bulk.png', 'time_full.png', 'positions_3d.png']: + 'time_bulk.png', 'time_full.png', 'positions_3d.png', 'step_count.png']: print(f" {f}") @@ -216,15 +288,18 @@ def main(): args.g4_hits = "g4_hits.npy" gpu = load_hits(args.gpu_hits) - g4 = load_hits(args.g4_hits) + g4_raw = load_hits(args.g4_hits) + g4_raw_shape = g4_raw.shape + + # Keep full G4 array for step count extraction + g4_full = g4_raw.copy() if g4_raw.shape[1] >= 5 else None # Normalize to (N, 4, 4) — take first 4 rows if more - if gpu.shape[1] > 4: - gpu = gpu[:, :4, :] - if g4.shape[1] > 4: - g4 = g4[:, :4, :] + gpu = gpu[:, :4, :] if gpu.shape[1] > 4 else gpu + g4 = g4_raw[:, :4, :] if g4_raw.shape[1] > 4 else g4_raw - make_plots(gpu, g4, args.outdir, title_extra=args.title) + make_plots(gpu, g4, args.outdir, title_extra=args.title, + g4_full=g4_full, g4_raw_shape=g4_raw_shape) if __name__ == "__main__": From 55989eef4c26d864f5b433864c630afa26bc684d Mon Sep 17 00:00:00 2001 From: Gabor Galgoczi Date: Thu, 2 Apr 2026 14:57:02 +0000 Subject: [PATCH 04/17] add DebugLite mode hint to run_and_compare output Print a note when running simulations that DebugLite mode is needed in the config JSON for step count plots, while HitPhoton is sufficient for hit-only analysis. --- optiphy/ana/run_and_compare.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/optiphy/ana/run_and_compare.py b/optiphy/ana/run_and_compare.py index 8ceaf220e..5575e3382 100755 --- a/optiphy/ana/run_and_compare.py +++ b/optiphy/ana/run_and_compare.py @@ -30,6 +30,9 @@ def run_simulation(gdml, config, macro, seed): """Run GPURaytrace and return (gpu_hits_path, g4_hits_path, gpu_nhits, g4_nhits).""" + print("NOTE: For step count plots, set \"mode\": \"DebugLite\" in the config JSON file.") + print(" For hit-only analysis, \"HitPhoton\" is sufficient.") + print() cmd = ["/opt/eic-opticks/bin/GPURaytrace", "-g", gdml, "-m", macro, "-s", str(seed)] if config: From 6c0d3ef4b7daf3147085f29c5452ef7a871d66cf Mon Sep 17 00:00:00 2001 From: Gabor Galgoczi Date: Sat, 4 Apr 2026 02:34:45 +0000 Subject: [PATCH 05/17] add benchmark script for apex.gdml speedup measurement --- examples/benchmark_apex.sh | 78 ++++++++++++++++++++++++++++++++++++++ 1 file changed, 78 insertions(+) create mode 100755 examples/benchmark_apex.sh diff --git a/examples/benchmark_apex.sh b/examples/benchmark_apex.sh new file mode 100755 index 000000000..2a4de6f83 --- /dev/null +++ b/examples/benchmark_apex.sh @@ -0,0 +1,78 @@ +#!/bin/bash +# benchmark_apex.sh — Measure GPU vs G4 speedup on apex.gdml +# +# Usage: +# ./examples/benchmark_apex.sh + +GDML="apex.gdml" +MACRO="tests/run.mac" +EPS="0.00001" +EPS0="0.0006" +OUTDIR="plots" +CONFIG="det_debug" + +if [ ! -f "$GDML" ]; then + echo "ERROR: $GDML not found. Run from the eic-opticks root directory." + exit 1 +fi + +echo "=== apex.gdml Benchmark ===" +echo "eps=$EPS, eps0=$EPS0" +echo "Running..." + +LOGFILE=$(mktemp /tmp/bench_XXXXXX.txt) +OPTICKS_MAX_BOUNCE=1000 \ +OPTICKS_PROPAGATE_EPSILON=$EPS \ +OPTICKS_PROPAGATE_EPSILON0=$EPS0 \ +GPURaytrace -g "$GDML" -m "$MACRO" -c "$CONFIG" &> "$LOGFILE" || true + +GPU_TIME=$(grep "Simulation time:" "$LOGFILE" | awk '{print $3}') +G4_LINE=$(grep "^ User=" "$LOGFILE" | tail -1) +G4_CPU=$(echo "$G4_LINE" | grep -oP 'User=\K[0-9.]+') +G4_WALL=$(echo "$G4_LINE" | grep -oP 'Real=\K[0-9.]+') +NPHOTONS=$(grep "NumCollected:" "$LOGFILE" | tail -1 | awk '{print $NF}') +GPU_HITS=$(grep "Opticks: NumHits:" "$LOGFILE" | awk '{print $NF}') +G4_HITS=$(grep "Geant4: NumHits:" "$LOGFILE" | awk '{print $NF}') + +if [ -z "$GPU_TIME" ] || [ -z "$G4_CPU" ]; then + echo "ERROR: Could not parse timing from output" + tail -30 "$LOGFILE" + rm -f "$LOGFILE" + exit 1 +fi + +python3 -c " +gpu = float('$GPU_TIME') +g4_cpu = float('$G4_CPU') +g4_wall = float('$G4_WALL') +nphotons = int('$NPHOTONS') +gpu_hits = int('$GPU_HITS') +g4_hits = int('$G4_HITS') +hit_diff = (gpu_hits - g4_hits) / g4_hits * 100 if g4_hits > 0 else 0 + +print() +print(f'Photons: {nphotons:>10,}') +print(f'GPU sim time: {gpu:>10.4f} s') +print(f'G4 CPU time: {g4_cpu:>10.2f} s') +print(f'G4 wall time: {g4_wall:>10.2f} s') +print() +print(f'Speedup (CPU): {g4_cpu/gpu:>10.0f}x') +print(f'Speedup (wall): {g4_wall/gpu:>10.0f}x') +print() +print(f'GPU rate: {nphotons/gpu/1e6:>10.1f} M photons/s') +print(f'G4 rate: {nphotons/g4_cpu/1e3:>10.1f} k photons/s') +print() +print(f'GPU hits: {gpu_hits:>10}') +print(f'G4 hits: {g4_hits:>10}') +print(f'Hit diff: {hit_diff:>+9.1f}%') +" + +rm -f "$LOGFILE" + +# Generate comparison plots if hit files exist +if [ -f "gpu_hits.npy" ] && [ -f "g4_hits.npy" ]; then + echo "" + echo "=== Generating comparison plots ===" + python3 optiphy/ana/run_and_compare.py --gpu-hits gpu_hits.npy --g4-hits g4_hits.npy --outdir "$OUTDIR" 2>&1 | tail -15 + echo "Plots saved to $OUTDIR/" +fi From 4f03e49cafed7dbaa262a8d5b4e417b60ab93311 Mon Sep 17 00:00:00 2001 From: Gabor Galgoczi Date: Thu, 26 Mar 2026 21:50:52 +0000 Subject: [PATCH 06/17] feat: add GPU vs G4 hit comparison script Side-by-side comparison of opticks GPU and standalone G4 simulation hits. Prints stats table (wavelength, time, position), statistical significance test, and optional wavelength/time distribution histograms (--histograms). --- ana/compare_gpu_g4.py | 286 ++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 286 insertions(+) create mode 100755 ana/compare_gpu_g4.py diff --git a/ana/compare_gpu_g4.py b/ana/compare_gpu_g4.py new file mode 100755 index 000000000..87fe210e6 --- /dev/null +++ b/ana/compare_gpu_g4.py @@ -0,0 +1,286 @@ +#!/usr/bin/env python +""" +compare_gpu_g4.py : Compare GPU (opticks) vs G4 (standalone) simulation hits +============================================================================= + +Reads GPU hit/photon arrays from an opticks event folder and G4 hits from +g4_hits.npy, then prints a side-by-side comparison table. + +Usage:: + + python ana/compare_gpu_g4.py + + # Auto-resolves A000 subfolder: + python ana/compare_gpu_g4.py /tmp/$USER/opticks/GEOM/GEOM/GPUPhotonSourceMinimal/ALL0_no_opticks_event_name g4_hits.npy +""" +import sys +import os +import argparse +import numpy as np + +FLAG_ENUM = { + 0x0004: "TORCH", 0x0008: "BULK_ABSORB", 0x0010: "BULK_REEMIT", + 0x0020: "BULK_SCATTER", 0x0040: "SURFACE_DETECT", 0x0080: "SURFACE_ABSORB", + 0x0100: "SURFACE_DREFLECT", 0x0200: "SURFACE_SREFLECT", + 0x0400: "BOUNDARY_REFLECT", 0x0800: "BOUNDARY_TRANSMIT", + 0x1000: "NAN_ABORT", 0x2000: "EFFICIENCY_COLLECT", 0x8000: "MISS", +} + + +def resolve_event_path(path): + if os.path.exists(os.path.join(path, "photon.npy")): + return path + a000 = os.path.join(path, "A000") + if os.path.exists(os.path.join(a000, "photon.npy")): + return a000 + if os.path.isdir(path): + for d in sorted(os.listdir(path)): + dp = os.path.join(path, d) + if os.path.isdir(dp) and os.path.exists(os.path.join(dp, "photon.npy")): + return dp + return path + + +def hit_stats(hits, label): + """Compute statistics dict from a (N, 4, 4) hit array.""" + n = len(hits) + if n == 0: + return dict(label=label, n=0) + wl = hits[:, 2, 3] + t = hits[:, 0, 3] + pos = hits[:, 0, :3] + r = np.sqrt(np.sum(pos ** 2, axis=1)) + return dict( + label=label, n=n, + wl_min=wl.min(), wl_max=wl.max(), wl_mean=wl.mean(), wl_std=wl.std(), + t_min=t.min(), t_max=t.max(), t_mean=t.mean(), t_std=t.std(), + r_min=r.min(), r_max=r.max(), r_mean=r.mean(), + x_mean=pos[:, 0].mean(), y_mean=pos[:, 1].mean(), z_mean=pos[:, 2].mean(), + ) + + +def print_comparison_table(gpu, g4, n_photons): + """Print side-by-side comparison.""" + w = 14 # column width + + print("=" * 70) + print("GPU vs G4 COMPARISON") + print("=" * 70) + + print(f"\n {'':30s} {'GPU':>{w}s} {'G4':>{w}s} {'Diff':>{w}s}") + print(f" {'-'*30} {'-'*w} {'-'*w} {'-'*w}") + + def row(name, gv, cv, fmt=".1f", diff_fmt=None): + if diff_fmt is None: + diff_fmt = fmt + gs = f"{gv:{fmt}}" if gv is not None else "—" + cs = f"{cv:{fmt}}" if cv is not None else "—" + if gv is not None and cv is not None: + d = cv - gv + ds = f"{d:{diff_fmt}}" + else: + ds = "—" + print(f" {name:30s} {gs:>{w}s} {cs:>{w}s} {ds:>{w}s}") + + row("Hits", gpu["n"], g4["n"], "d") + if n_photons and n_photons > 0: + row("Hit rate (%)", 100.0 * gpu["n"] / n_photons, 100.0 * g4["n"] / n_photons, ".2f") + + if gpu["n"] > 0 and g4["n"] > 0: + ratio = g4["n"] / gpu["n"] + print(f" {'Ratio G4/GPU':30s} {'':>{w}s} {'':>{w}s} {ratio:>{w}.3f}") + + if gpu["n"] == 0 or g4["n"] == 0: + print("\n Cannot compare distributions — one side has zero hits.") + return + + print() + row("Wavelength min (nm)", gpu["wl_min"], g4["wl_min"]) + row("Wavelength max (nm)", gpu["wl_max"], g4["wl_max"]) + row("Wavelength mean (nm)", gpu["wl_mean"], g4["wl_mean"]) + row("Wavelength std (nm)", gpu["wl_std"], g4["wl_std"]) + + print() + row("Time min (ns)", gpu["t_min"], g4["t_min"], ".3f") + row("Time max (ns)", gpu["t_max"], g4["t_max"], ".3f") + row("Time mean (ns)", gpu["t_mean"], g4["t_mean"], ".3f") + row("Time std (ns)", gpu["t_std"], g4["t_std"], ".3f") + + print() + row("Radius min (mm)", gpu["r_min"], g4["r_min"], ".2f") + row("Radius max (mm)", gpu["r_max"], g4["r_max"], ".2f") + row("Radius mean (mm)", gpu["r_mean"], g4["r_mean"], ".2f") + + print() + row("Mean X (mm)", gpu["x_mean"], g4["x_mean"], ".2f") + row("Mean Y (mm)", gpu["y_mean"], g4["y_mean"], ".2f") + row("Mean Z (mm)", gpu["z_mean"], g4["z_mean"], ".2f") + + # Statistical significance + print() + if n_photons and n_photons > 0: + p_pool = (gpu["n"] + g4["n"]) / (2 * n_photons) + std = np.sqrt(p_pool * (1 - p_pool) / n_photons) + if std > 0: + z = abs(gpu["n"] / n_photons - g4["n"] / n_photons) / (std * np.sqrt(2)) + expected_fluct = std * np.sqrt(2) * n_photons + print(f" {'Z-score (hit count)':30s} {z:>{w}.1f}") + print(f" {'Expected 1σ fluctuation':30s} {expected_fluct:>{w}.0f} hits") + if z > 3: + print(f" ** Statistically significant difference (>{3}σ) **") + else: + print(f" Within statistical expectations (<3σ)") + print() + + +def print_gpu_outcomes(photon): + """Print GPU photon outcome summary.""" + q3 = photon[:, 3, :].view(np.uint32) + flag = q3[:, 0] & 0xFFFF + + print("=" * 70) + print("GPU PHOTON OUTCOMES") + print("=" * 70) + + n = len(flag) + vals, counts = np.unique(flag, return_counts=True) + order = np.argsort(-counts) + + print(f"\n {'Flag':<22s} {'Count':>8s} {'%':>7s}") + print(f" {'-'*22} {'-'*8} {'-'*7}") + for idx in order: + f = vals[idx] + c = counts[idx] + name = FLAG_ENUM.get(f, f"0x{f:04x}") + print(f" {name:<22s} {c:8d} {100*c/n:6.1f}%") + print() + + +def print_wavelength_histograms(gpu_hits, g4_hits): + """Print overlaid wavelength histograms.""" + if len(gpu_hits) == 0 or len(g4_hits) == 0: + return + + gpu_wl = gpu_hits[:, 2, 3] + g4_wl = g4_hits[:, 2, 3] + + wl_min = min(gpu_wl.min(), g4_wl.min()) + wl_max = max(gpu_wl.max(), g4_wl.max()) + bins = np.arange(max(100, np.floor(wl_min / 25) * 25), + min(800, np.ceil(wl_max / 25) * 25 + 25), 25) + + gpu_counts, _ = np.histogram(gpu_wl, bins=bins) + g4_counts, _ = np.histogram(g4_wl, bins=bins) + + # Normalize to same total for shape comparison + gpu_norm = gpu_counts / len(gpu_hits) * 1000 + g4_norm = g4_counts / len(g4_hits) * 1000 + + print("=" * 70) + print("WAVELENGTH DISTRIBUTION (per 1000 hits)") + print("=" * 70) + print(f"\n {'Bin (nm)':<14s} {'GPU':>8s} {'G4':>8s} {'GPU':^20s} {'G4':^20s}") + print(f" {'-'*14} {'-'*8} {'-'*8} {'-'*20} {'-'*20}") + + max_bar = 20 + scale = max(gpu_norm.max(), g4_norm.max()) + if scale == 0: + scale = 1 + + for i in range(len(bins) - 1): + if gpu_counts[i] == 0 and g4_counts[i] == 0: + continue + gpu_bar = "#" * int(gpu_norm[i] / scale * max_bar) + g4_bar = "#" * int(g4_norm[i] / scale * max_bar) + print(f" {bins[i]:5.0f}-{bins[i+1]:5.0f} {gpu_norm[i]:8.1f} {g4_norm[i]:8.1f}" + f" {gpu_bar:<20s} {g4_bar:<20s}") + print() + + +def print_time_histograms(gpu_hits, g4_hits): + """Print overlaid time histograms.""" + if len(gpu_hits) == 0 or len(g4_hits) == 0: + return + + gpu_t = gpu_hits[:, 0, 3] + g4_t = g4_hits[:, 0, 3] + + t_max = max(gpu_t.max(), g4_t.max()) + bin_size = max(1.0, np.ceil(t_max / 15)) + bins = np.arange(0, t_max + bin_size, bin_size) + + gpu_counts, _ = np.histogram(gpu_t, bins=bins) + g4_counts, _ = np.histogram(g4_t, bins=bins) + + gpu_norm = gpu_counts / len(gpu_hits) * 1000 + g4_norm = g4_counts / len(g4_hits) * 1000 + + print("=" * 70) + print("TIME DISTRIBUTION (per 1000 hits)") + print("=" * 70) + print(f"\n {'Bin (ns)':<14s} {'GPU':>8s} {'G4':>8s} {'GPU':^20s} {'G4':^20s}") + print(f" {'-'*14} {'-'*8} {'-'*8} {'-'*20} {'-'*20}") + + max_bar = 20 + scale = max(gpu_norm.max(), g4_norm.max()) + if scale == 0: + scale = 1 + + for i in range(len(bins) - 1): + if gpu_counts[i] == 0 and g4_counts[i] == 0: + continue + gpu_bar = "#" * int(gpu_norm[i] / scale * max_bar) + g4_bar = "#" * int(g4_norm[i] / scale * max_bar) + print(f" {bins[i]:5.1f}-{bins[i+1]:5.1f} {gpu_norm[i]:8.1f} {g4_norm[i]:8.1f}" + f" {gpu_bar:<20s} {g4_bar:<20s}") + print() + + +def main(): + parser = argparse.ArgumentParser( + description="Compare GPU (opticks) vs G4 (standalone) simulation hits", + formatter_class=argparse.RawDescriptionHelpFormatter, + epilog=__doc__, + ) + parser.add_argument("gpu_path", help="Path to GPU opticks event folder") + parser.add_argument("g4_hits", help="Path to G4 hits file (g4_hits.npy)") + parser.add_argument("--histograms", action="store_true", + help="Show wavelength and time distribution histograms") + + args = parser.parse_args() + + gpu_path = resolve_event_path(args.gpu_path) + if not os.path.exists(os.path.join(gpu_path, "photon.npy")): + print(f"Error: photon.npy not found in {gpu_path}") + sys.exit(1) + if not os.path.exists(args.g4_hits): + print(f"Error: {args.g4_hits} not found") + sys.exit(1) + + # Load GPU arrays + gpu_hits = np.load(os.path.join(gpu_path, "hit.npy")) if os.path.exists(os.path.join(gpu_path, "hit.npy")) else np.zeros((0, 4, 4), dtype=np.float32) + gpu_photon = np.load(os.path.join(gpu_path, "photon.npy")) + n_photons = len(gpu_photon) + + # Load G4 hits + g4_hits = np.load(args.g4_hits) + + print(f"\nGPU event: {gpu_path}") + print(f"G4 hits: {args.g4_hits}") + print(f"Total photons: {n_photons}\n") + + # Compute stats + gpu_stats = hit_stats(gpu_hits, "GPU") + g4_stats = hit_stats(g4_hits, "G4") + + # Print tables + print_comparison_table(gpu_stats, g4_stats, n_photons) + print_gpu_outcomes(gpu_photon) + + if args.histograms: + print_wavelength_histograms(gpu_hits, g4_hits) + print_time_histograms(gpu_hits, g4_hits) + + +if __name__ == "__main__": + main() From 9ae937fb6156aa12d326022347ede9a406dfbc56 Mon Sep 17 00:00:00 2001 From: Gabor Galgoczi Date: Thu, 26 Mar 2026 21:51:45 +0000 Subject: [PATCH 07/17] move compare_gpu_g4.py to optiphy/ana --- {ana => optiphy/ana}/compare_gpu_g4.py | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename {ana => optiphy/ana}/compare_gpu_g4.py (100%) diff --git a/ana/compare_gpu_g4.py b/optiphy/ana/compare_gpu_g4.py similarity index 100% rename from ana/compare_gpu_g4.py rename to optiphy/ana/compare_gpu_g4.py From b2da948b04776edccd27e078a132585c20f2c33e Mon Sep 17 00:00:00 2001 From: Gabor Galgoczi Date: Thu, 26 Mar 2026 21:50:41 +0000 Subject: [PATCH 08/17] feat: add DUNE example detector geometry for validation testing LAr scintillator tile with two-stage WLS readout (pTP + blue WLS acrylic), 30 SiPMs along edge, Vikuiti reflective foil wrapping. Includes EFFICIENCY and SensDet tags for SiPM hit collection. --- tests/geom/DUNE_example_detector.gdml | 373 ++++++++++++++++++++++++++ 1 file changed, 373 insertions(+) create mode 100644 tests/geom/DUNE_example_detector.gdml diff --git a/tests/geom/DUNE_example_detector.gdml b/tests/geom/DUNE_example_detector.gdml new file mode 100644 index 000000000..11ae55814 --- /dev/null +++ b/tests/geom/DUNE_example_detector.gdml @@ -0,0 +1,373 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + From 50251dbc5cff74c77b4f3ecb35d20e3d9499f661 Mon Sep 17 00:00:00 2001 From: Gabor Galgoczi Date: Sun, 29 Mar 2026 16:23:03 +0000 Subject: [PATCH 09/17] add savephotonhistory config flag and GPU/G4 hit .npy saving When savephotonhistory=true in config JSON, saves full SEvt arrays (photon, record, seq, hit) plus gpu_hits.npy and g4_hits.npy for distribution comparison. G4 hits collected via thread-safe accumulator in EndOfEventAction. GPURaytrace now accepts -c config flag. --- src/GPURaytrace.cpp | 14 ++++++++ src/GPURaytrace.h | 84 +++++++++++++++++++++++++++------------------ src/config.cpp | 3 ++ src/config.h | 2 ++ 4 files changed, 70 insertions(+), 33 deletions(-) diff --git a/src/GPURaytrace.cpp b/src/GPURaytrace.cpp index c97c7ec11..ddc3f494a 100644 --- a/src/GPURaytrace.cpp +++ b/src/GPURaytrace.cpp @@ -12,6 +12,7 @@ #include "sysrap/OPTICKS_LOG.hh" +#include "config.h" #include "GPURaytrace.h" #include "G4RunManager.hh" @@ -68,6 +69,11 @@ int main(int argc, char **argv) .nargs(1) .store_into(macro_name); + program.add_argument("-c", "--config") + .help("config file name (without .json extension)") + .default_value(string("")) + .nargs(1); + program.add_argument("-i", "--interactive") .help("whether to open an interactive window with a viewer") .flag() @@ -108,6 +114,14 @@ int main(int argc, char **argv) G4App *g4app = new G4App(gdml_file); + // Load config and apply savephotonhistory flag if provided + string config_name = program.get("--config"); + if (!config_name.empty()) + { + gphox::Config cfg(config_name); + g4app->run_act_->fSavePhotonHistory = cfg.savephotonhistory; + } + ActionInitialization *actionInit = new ActionInitialization(g4app); run_mgr->SetUserInitialization(actionInit); run_mgr->SetUserInitialization(g4app->det_cons_); diff --git a/src/GPURaytrace.h b/src/GPURaytrace.h index f8443a167..958777d83 100644 --- a/src/GPURaytrace.h +++ b/src/GPURaytrace.h @@ -1,3 +1,4 @@ +#include #include #include #include @@ -47,6 +48,8 @@ namespace { G4Mutex genstep_mutex = G4MUTEX_INITIALIZER; +G4Mutex g4hits_mutex = G4MUTEX_INITIALIZER; +std::vector> g4_accumulated_hits; } bool IsSubtractionSolid(G4VSolid *solid) @@ -303,7 +306,7 @@ struct PrimaryGenerator : G4VUserPrimaryGeneratorAction G4PrimaryVertex *vertex = new G4PrimaryVertex(position_mm, time_ns); G4PrimaryParticle *particle = new G4PrimaryParticle(G4Electron::Definition()); - particle->SetKineticEnergy(5 * GeV); + particle->SetKineticEnergy(10 * MeV); particle->SetMomentumDirection(direction); vertex->SetPrimary(particle); event->AddPrimaryVertex(vertex); @@ -331,7 +334,26 @@ struct EventAction : G4UserEventAction for (G4int i = 0; i < hce->GetNumberOfCollections(); i++) { G4VHitsCollection *hc = hce->GetHC(i); - if (hc) + if (!hc) continue; + + PhotonHitsCollection *phc = dynamic_cast(hc); + if (phc) + { + G4AutoLock lock(&g4hits_mutex); + for (size_t j = 0; j < phc->entries(); j++) + { + PhotonHit* hit = (*phc)[j]; + float wl = 1239.84198f / static_cast(hit->fenergy); + g4_accumulated_hits.push_back({ + float(hit->fposition.x()), float(hit->fposition.y()), float(hit->fposition.z()), float(hit->ftime), + float(hit->fdirection.x()), float(hit->fdirection.y()), float(hit->fdirection.z()), 0.f, + float(hit->fpolarization.x()), float(hit->fpolarization.y()), float(hit->fpolarization.z()), wl, + 0.f, 0.f, 0.f, 0.f + }); + } + fTotalG4Hits += phc->entries(); + } + else { fTotalG4Hits += hc->GetSize(); } @@ -348,6 +370,7 @@ struct EventAction : G4UserEventAction struct RunAction : G4UserRunAction { EventAction *fEventAction; + bool fSavePhotonHistory{false}; RunAction(EventAction *eventAction) : fEventAction(eventAction) { @@ -379,44 +402,39 @@ struct RunAction : G4UserRunAction std::cout << "Opticks: NumCollected: " << sev->GetNumPhotonCollected(0) << std::endl; std::cout << "Opticks: NumHits: " << num_hits << std::endl; std::cout << "Geant4: NumHits: " << fEventAction->GetTotalG4Hits() << std::endl; - std::ofstream outFile("opticks_hits_output.txt"); - if (!outFile.is_open()) - { - std::cerr << "Error opening output file!" << std::endl; - return; - } - for (int idx = 0; idx < int(num_hits); idx++) + if (fSavePhotonHistory) { - sphoton hit; - sev->getHit(hit, idx); - G4ThreeVector position = G4ThreeVector(hit.pos.x, hit.pos.y, hit.pos.z); - G4ThreeVector direction = G4ThreeVector(hit.mom.x, hit.mom.y, hit.mom.z); - G4ThreeVector polarization = G4ThreeVector(hit.pol.x, hit.pol.y, hit.pol.z); - int theCreationProcessid; - if (OpticksPhoton::HasCerenkovFlag(hit.flagmask)) - { - theCreationProcessid = 0; - } - else if (OpticksPhoton::HasScintillationFlag(hit.flagmask)) + // Save full SEvt (photon, record, seq, hit) when DebugLite/DebugHeavy + sev->save(); + std::cout << "SEvt::save() complete" << std::endl; + + // Save GPU hits as .npy (sphoton layout: N x 4 x 4 float32) { - theCreationProcessid = 1; + NP* gpu_h = NP::Make(num_hits, 4, 4); + for (unsigned idx = 0; idx < num_hits; idx++) + { + sphoton hit; + sev->getHit(hit, idx); + memcpy(gpu_h->bytes() + idx * sizeof(sphoton), &hit, sizeof(sphoton)); + } + gpu_h->save("gpu_hits.npy"); + std::cout << "Saved GPU hits: " << num_hits << " to gpu_hits.npy" << std::endl; } - else + + // Save G4 hits as .npy (same layout: N x 4 x 4 float32) { - theCreationProcessid = -1; + G4AutoLock lock(&g4hits_mutex); + size_t ng4 = g4_accumulated_hits.size(); + if (ng4 > 0) + { + NP* g4h = NP::Make(ng4, 4, 4); + memcpy(g4h->bytes(), g4_accumulated_hits.data(), ng4 * 16 * sizeof(float)); + g4h->save("g4_hits.npy"); + std::cout << "Saved G4 hits: " << ng4 << " to g4_hits.npy" << std::endl; + } } - // std::cout << "Adding hit from Opticks:" << hit.wavelength << " " << position << " " << direction - // << " - // " - // << polarization << std::endl; - outFile << hit.time << " " << hit.wavelength << " " << "(" << position.x() << ", " << position.y() - << ", " << position.z() << ") " << "(" << direction.x() << ", " << direction.y() << ", " - << direction.z() << ") " << "(" << polarization.x() << ", " << polarization.y() << ", " - << polarization.z() << ") " << "CreationProcessID=" << theCreationProcessid << std::endl; } - - outFile.close(); } } }; diff --git a/src/config.cpp b/src/config.cpp index 844244017..26cb221e0 100644 --- a/src/config.cpp +++ b/src/config.cpp @@ -136,6 +136,9 @@ void Config::ReadConfig(std::string filepath) SEventConfig::SetEventMode( string(event_["mode"]).c_str() ); SEventConfig::SetMaxSlot( event_["maxslot"] ); + + if (event_.contains("savephotonhistory")) + savephotonhistory = event_["savephotonhistory"].get(); } catch (nlohmann::json::exception& e) { std::string errmsg{"Failed reading config parameters from " + filepath + "\n" + e.what()}; diff --git a/src/config.h b/src/config.h index 1fc5c838d..3dc2c1ebf 100644 --- a/src/config.h +++ b/src/config.h @@ -26,6 +26,8 @@ class Config storch torch; + bool savephotonhistory{false}; + private: std::string Locate(std::string filename) const; From cfe53b9188e06a3c9a0d1eb5e0200bdc3981233c Mon Sep 17 00:00:00 2001 From: Gabor Galgoczi Date: Tue, 7 Apr 2026 14:38:20 +0000 Subject: [PATCH 10/17] remove X-offset dodge between GPU and G4 data points in validation plots MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Plot both datasets at true bin centers instead of shifting ±17.5% of bin width for readability. --- optiphy/ana/run_and_compare.py | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/optiphy/ana/run_and_compare.py b/optiphy/ana/run_and_compare.py index 5575e3382..1cd599cd2 100755 --- a/optiphy/ana/run_and_compare.py +++ b/optiphy/ana/run_and_compare.py @@ -67,11 +67,10 @@ def plot_with_errors(ax, data1, data2, bins, label1, label2, xlabel): h1, edges = np.histogram(data1, bins=bins) h2, _ = np.histogram(data2, bins=bins) centers = (edges[:-1] + edges[1:]) / 2 - width = (edges[1] - edges[0]) * 0.35 - ax.errorbar(centers - width / 2, h1, yerr=np.sqrt(np.maximum(h1, 1)), + ax.errorbar(centers, h1, yerr=np.sqrt(np.maximum(h1, 1)), fmt='o', color='dodgerblue', markersize=4, capsize=2, linewidth=1, label=label1) - ax.errorbar(centers + width / 2, h2, yerr=np.sqrt(np.maximum(h2, 1)), + ax.errorbar(centers, h2, yerr=np.sqrt(np.maximum(h2, 1)), fmt='s', color='orangered', markersize=4, capsize=2, linewidth=1, label=label2) ax.set_xlabel(xlabel) From fa816eef5c4ffb4592582a5afe268b17e08f1e2d Mon Sep 17 00:00:00 2001 From: Gabor Galgoczi Date: Tue, 7 Apr 2026 14:55:29 +0000 Subject: [PATCH 11/17] revert primary energy to 5 GeV for CI raindrop test The 10 MeV change was for APEX validation and belongs on the physics branch, not analysis-tools. --- src/GPURaytrace.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/GPURaytrace.h b/src/GPURaytrace.h index 958777d83..c01f8a263 100644 --- a/src/GPURaytrace.h +++ b/src/GPURaytrace.h @@ -306,7 +306,7 @@ struct PrimaryGenerator : G4VUserPrimaryGeneratorAction G4PrimaryVertex *vertex = new G4PrimaryVertex(position_mm, time_ns); G4PrimaryParticle *particle = new G4PrimaryParticle(G4Electron::Definition()); - particle->SetKineticEnergy(10 * MeV); + particle->SetKineticEnergy(5 * GeV); particle->SetMomentumDirection(direction); vertex->SetPrimary(particle); event->AddPrimaryVertex(vertex); From 027719569970e562c46a351db172abd773aaf24e Mon Sep 17 00:00:00 2001 From: Gabor Galgoczi Date: Thu, 9 Apr 2026 22:55:17 +0000 Subject: [PATCH 12/17] add apex.gdml detector geometry for validation APEX liquid argon scintillator tile with WLS light guide, Vikuiti reflective foils, and 31 SiPMs. Used for GPU vs G4 benchmarking. --- apex.gdml | 379 ++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 379 insertions(+) create mode 100644 apex.gdml diff --git a/apex.gdml b/apex.gdml new file mode 100644 index 000000000..49d36656b --- /dev/null +++ b/apex.gdml @@ -0,0 +1,379 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + From 1bfcb760cb9881459f252d69a0535343fa29cb3c Mon Sep 17 00:00:00 2001 From: Gabor Galgoczi Date: Sun, 12 Apr 2026 23:07:36 +0000 Subject: [PATCH 13/17] Apply clang-format Microsoft style to flagged files --- src/GPURaytrace.cpp | 2 +- src/GPURaytrace.h | 24 +++---- src/config.cpp | 158 ++++++++++++++++++++++---------------------- src/config.h | 28 ++++---- 4 files changed, 105 insertions(+), 107 deletions(-) diff --git a/src/GPURaytrace.cpp b/src/GPURaytrace.cpp index ddc3f494a..f0590db8d 100644 --- a/src/GPURaytrace.cpp +++ b/src/GPURaytrace.cpp @@ -12,8 +12,8 @@ #include "sysrap/OPTICKS_LOG.hh" -#include "config.h" #include "GPURaytrace.h" +#include "config.h" #include "G4RunManager.hh" #include "G4RunManagerFactory.hh" diff --git a/src/GPURaytrace.h b/src/GPURaytrace.h index c01f8a263..3ea2e4a5d 100644 --- a/src/GPURaytrace.h +++ b/src/GPURaytrace.h @@ -50,7 +50,7 @@ namespace G4Mutex genstep_mutex = G4MUTEX_INITIALIZER; G4Mutex g4hits_mutex = G4MUTEX_INITIALIZER; std::vector> g4_accumulated_hits; -} +} // namespace bool IsSubtractionSolid(G4VSolid *solid) { @@ -334,22 +334,22 @@ struct EventAction : G4UserEventAction for (G4int i = 0; i < hce->GetNumberOfCollections(); i++) { G4VHitsCollection *hc = hce->GetHC(i); - if (!hc) continue; + if (!hc) + continue; - PhotonHitsCollection *phc = dynamic_cast(hc); + PhotonHitsCollection *phc = dynamic_cast(hc); if (phc) { G4AutoLock lock(&g4hits_mutex); for (size_t j = 0; j < phc->entries(); j++) { - PhotonHit* hit = (*phc)[j]; + PhotonHit *hit = (*phc)[j]; float wl = 1239.84198f / static_cast(hit->fenergy); - g4_accumulated_hits.push_back({ - float(hit->fposition.x()), float(hit->fposition.y()), float(hit->fposition.z()), float(hit->ftime), - float(hit->fdirection.x()), float(hit->fdirection.y()), float(hit->fdirection.z()), 0.f, - float(hit->fpolarization.x()), float(hit->fpolarization.y()), float(hit->fpolarization.z()), wl, - 0.f, 0.f, 0.f, 0.f - }); + g4_accumulated_hits.push_back( + {float(hit->fposition.x()), float(hit->fposition.y()), float(hit->fposition.z()), + float(hit->ftime), float(hit->fdirection.x()), float(hit->fdirection.y()), + float(hit->fdirection.z()), 0.f, float(hit->fpolarization.x()), + float(hit->fpolarization.y()), float(hit->fpolarization.z()), wl, 0.f, 0.f, 0.f, 0.f}); } fTotalG4Hits += phc->entries(); } @@ -411,7 +411,7 @@ struct RunAction : G4UserRunAction // Save GPU hits as .npy (sphoton layout: N x 4 x 4 float32) { - NP* gpu_h = NP::Make(num_hits, 4, 4); + NP *gpu_h = NP::Make(num_hits, 4, 4); for (unsigned idx = 0; idx < num_hits; idx++) { sphoton hit; @@ -428,7 +428,7 @@ struct RunAction : G4UserRunAction size_t ng4 = g4_accumulated_hits.size(); if (ng4 > 0) { - NP* g4h = NP::Make(ng4, 4, 4); + NP *g4h = NP::Make(ng4, 4, 4); memcpy(g4h->bytes(), g4_accumulated_hits.data(), ng4 * 16 * sizeof(float)); g4h->save("g4_hits.npy"); std::cout << "Saved G4 hits: " << ng4 << " to g4_hits.npy" << std::endl; diff --git a/src/config.cpp b/src/config.cpp index 26cb221e0..8c898cdce 100644 --- a/src/config.cpp +++ b/src/config.cpp @@ -6,15 +6,16 @@ #include #include -#include #include +#include #include "sysrap/SEventConfig.hh" #include "config.h" #include "config_path.h" -namespace gphox { +namespace gphox +{ using namespace std; @@ -28,10 +29,9 @@ bool FileExists(const std::string &path) return std::filesystem::exists(path, ec) && !ec; } -Config::Config(std::string config_name) : - name{std::getenv("GPHOX_CONFIG") ? std::getenv("GPHOX_CONFIG") : config_name} +Config::Config(std::string config_name) : name{std::getenv("GPHOX_CONFIG") ? std::getenv("GPHOX_CONFIG") : config_name} { - ReadConfig(Locate(name + ".json")); + ReadConfig(Locate(name + ".json")); } std::string Config::PtxPath(const std::string &ptx_name) @@ -54,96 +54,96 @@ std::string Config::PtxPath(const std::string &ptx_name) std::string Config::Locate(std::string filename) const { - std::vector search_paths; + std::vector search_paths; + + const std::string user_dir{std::getenv("GPHOX_CONFIG_DIR") ? std::getenv("GPHOX_CONFIG_DIR") : ""}; - const std::string user_dir{std::getenv("GPHOX_CONFIG_DIR") ? std::getenv("GPHOX_CONFIG_DIR") : ""}; + if (user_dir.empty()) + { + std::string paths(GPHOX_CONFIG_SEARCH_PATHS); - if (user_dir.empty()) - { - std::string paths(GPHOX_CONFIG_SEARCH_PATHS); + size_t last = 0; + size_t next = 0; + while ((next = paths.find(':', last)) != std::string::npos) + { + search_paths.push_back(paths.substr(last, next - last)); + last = next + 1; + } - size_t last = 0; - size_t next = 0; - while ((next = paths.find(':', last)) != std::string::npos) + search_paths.push_back(paths.substr(last)); + } + else { - search_paths.push_back(paths.substr(last, next-last)); - last = next + 1; + search_paths.push_back(user_dir); } - search_paths.push_back(paths.substr(last)); - } - else - { - search_paths.push_back(user_dir); - } - - struct stat buffer; - std::string filepath{""}; - for (std::string path : search_paths) - { - std::string fpath{path + "/" + filename}; - if (stat(fpath.c_str(), &buffer) == 0) + struct stat buffer; + std::string filepath{""}; + for (std::string path : search_paths) { - filepath = fpath; - break; + std::string fpath{path + "/" + filename}; + if (stat(fpath.c_str(), &buffer) == 0) + { + filepath = fpath; + break; + } } - } - if (filepath.empty()) - { - std::string errmsg{"Could not find config file \"" + filename + "\" in "}; - for (std::string path : search_paths) errmsg += (path + ":"); - throw std::runtime_error(errmsg); - } + if (filepath.empty()) + { + std::string errmsg{"Could not find config file \"" + filename + "\" in "}; + for (std::string path : search_paths) + errmsg += (path + ":"); + throw std::runtime_error(errmsg); + } - return filepath; + return filepath; } - /** * Expects a valid filepath. */ void Config::ReadConfig(std::string filepath) { - nlohmann::json json; - - try { - std::ifstream ifs(filepath); - ifs >> json; - - nlohmann::json torch_ = json["torch"]; - - torch = { - .gentype = OpticksGenstep_::Type(torch_["gentype"]), - .trackid = torch_["trackid"], - .matline = torch_["matline"], - .numphoton = torch_["numphoton"], - .pos = make_float3(torch_["pos"][0], torch_["pos"][1], torch_["pos"][2]), - .time = torch_["time"], - .mom = normalize(make_float3(torch_["mom"][0], torch_["mom"][1], torch_["mom"][2])), - .weight = torch_["weight"], - .pol = make_float3(torch_["pol"][0], torch_["pol"][1], torch_["pol"][2]), - .wavelength = torch_["wavelength"], - .zenith = make_float2(torch_["zenith"][0], torch_["zenith"][1]), - .azimuth = make_float2(torch_["azimuth"][0], torch_["azimuth"][1]), - .radius = torch_["radius"], - .distance = torch_["distance"], - .mode = torch_["mode"], - .type = storchtype::Type(torch_["type"]) - }; - - nlohmann::json event_ = json["event"]; - - SEventConfig::SetEventMode( string(event_["mode"]).c_str() ); - SEventConfig::SetMaxSlot( event_["maxslot"] ); - - if (event_.contains("savephotonhistory")) - savephotonhistory = event_["savephotonhistory"].get(); - } - catch (nlohmann::json::exception& e) { - std::string errmsg{"Failed reading config parameters from " + filepath + "\n" + e.what()}; - throw std::runtime_error{errmsg}; - } -} + nlohmann::json json; + try + { + std::ifstream ifs(filepath); + ifs >> json; + + nlohmann::json torch_ = json["torch"]; + + torch = {.gentype = OpticksGenstep_::Type(torch_["gentype"]), + .trackid = torch_["trackid"], + .matline = torch_["matline"], + .numphoton = torch_["numphoton"], + .pos = make_float3(torch_["pos"][0], torch_["pos"][1], torch_["pos"][2]), + .time = torch_["time"], + .mom = normalize(make_float3(torch_["mom"][0], torch_["mom"][1], torch_["mom"][2])), + .weight = torch_["weight"], + .pol = make_float3(torch_["pol"][0], torch_["pol"][1], torch_["pol"][2]), + .wavelength = torch_["wavelength"], + .zenith = make_float2(torch_["zenith"][0], torch_["zenith"][1]), + .azimuth = make_float2(torch_["azimuth"][0], torch_["azimuth"][1]), + .radius = torch_["radius"], + .distance = torch_["distance"], + .mode = torch_["mode"], + .type = storchtype::Type(torch_["type"])}; + + nlohmann::json event_ = json["event"]; + + SEventConfig::SetEventMode(string(event_["mode"]).c_str()); + SEventConfig::SetMaxSlot(event_["maxslot"]); + + if (event_.contains("savephotonhistory")) + savephotonhistory = event_["savephotonhistory"].get(); + } + catch (nlohmann::json::exception &e) + { + std::string errmsg{"Failed reading config parameters from " + filepath + "\n" + e.what()}; + throw std::runtime_error{errmsg}; + } } + +} // namespace gphox diff --git a/src/config.h b/src/config.h index 3dc2c1ebf..28e26416c 100644 --- a/src/config.h +++ b/src/config.h @@ -7,31 +7,29 @@ #include "sysrap/srng.h" #include "sysrap/storch.h" -namespace gphox { - +namespace gphox +{ /** * Provides access to all configuration types and data. */ class Config { - public: - - Config(std::string config_name = "dev"); - - static std::string PtxPath(const std::string &ptx_name = "CSGOptiX7.ptx"); + public: + Config(std::string config_name = "dev"); - /// A unique name associated with this Config - std::string name; + static std::string PtxPath(const std::string &ptx_name = "CSGOptiX7.ptx"); - storch torch; + /// A unique name associated with this Config + std::string name; - bool savephotonhistory{false}; + storch torch; - private: + bool savephotonhistory{false}; - std::string Locate(std::string filename) const; - void ReadConfig(std::string filepath); + private: + std::string Locate(std::string filename) const; + void ReadConfig(std::string filepath); }; -} +} // namespace gphox From bccc7187345dcc1ffb7b191799603afb21debefc Mon Sep 17 00:00:00 2001 From: Gabor Galgoczi Date: Sun, 12 Apr 2026 23:24:45 +0000 Subject: [PATCH 14/17] Move benchmark script from examples/ to benchmarks/ --- {examples => benchmarks}/benchmark_apex.sh | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename {examples => benchmarks}/benchmark_apex.sh (100%) diff --git a/examples/benchmark_apex.sh b/benchmarks/benchmark_apex.sh similarity index 100% rename from examples/benchmark_apex.sh rename to benchmarks/benchmark_apex.sh From 9b0e6e0e099de8beaa54d11ad1d7c6772105e1ca Mon Sep 17 00:00:00 2001 From: Gabor Galgoczi Date: Sun, 12 Apr 2026 23:30:13 +0000 Subject: [PATCH 15/17] Move apex.gdml to tests/geom/ --- apex.gdml => tests/geom/apex.gdml | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename apex.gdml => tests/geom/apex.gdml (100%) diff --git a/apex.gdml b/tests/geom/apex.gdml similarity index 100% rename from apex.gdml rename to tests/geom/apex.gdml From 3458524a3fb2f79a34f2473b47b3dc193ba3482c Mon Sep 17 00:00:00 2001 From: Gabor Galgoczi Date: Mon, 13 Apr 2026 00:02:00 +0000 Subject: [PATCH 16/17] Update benchmark_apex.sh path after apex.gdml move to tests/geom/ --- benchmarks/benchmark_apex.sh | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/benchmarks/benchmark_apex.sh b/benchmarks/benchmark_apex.sh index 2a4de6f83..4ff805b43 100755 --- a/benchmarks/benchmark_apex.sh +++ b/benchmarks/benchmark_apex.sh @@ -2,9 +2,9 @@ # benchmark_apex.sh — Measure GPU vs G4 speedup on apex.gdml # # Usage: -# ./examples/benchmark_apex.sh +# ./benchmarks/benchmark_apex.sh -GDML="apex.gdml" +GDML="tests/geom/apex.gdml" MACRO="tests/run.mac" EPS="0.00001" EPS0="0.0006" From fde7d90bbc3462106044ce2185cbbabfa90a8ff8 Mon Sep 17 00:00:00 2001 From: Gabor Galgoczi Date: Tue, 14 Apr 2026 12:10:57 +0000 Subject: [PATCH 17/17] refactor(benchmarks): extract speedup formatter to optiphy/ana/print_speedup.py Move the inline python3 -c "..." block out of benchmarks/benchmark_apex.sh into optiphy/ana/print_speedup.py. Shell now only runs the benchmark and parses log values; Python module owns the formatting and is testable standalone. Behavior verified identical against same input values. --- benchmarks/benchmark_apex.sh | 27 ++---------------------- optiphy/ana/print_speedup.py | 41 ++++++++++++++++++++++++++++++++++++ 2 files changed, 43 insertions(+), 25 deletions(-) create mode 100644 optiphy/ana/print_speedup.py diff --git a/benchmarks/benchmark_apex.sh b/benchmarks/benchmark_apex.sh index 4ff805b43..cab100fde 100755 --- a/benchmarks/benchmark_apex.sh +++ b/benchmarks/benchmark_apex.sh @@ -41,31 +41,8 @@ if [ -z "$GPU_TIME" ] || [ -z "$G4_CPU" ]; then exit 1 fi -python3 -c " -gpu = float('$GPU_TIME') -g4_cpu = float('$G4_CPU') -g4_wall = float('$G4_WALL') -nphotons = int('$NPHOTONS') -gpu_hits = int('$GPU_HITS') -g4_hits = int('$G4_HITS') -hit_diff = (gpu_hits - g4_hits) / g4_hits * 100 if g4_hits > 0 else 0 - -print() -print(f'Photons: {nphotons:>10,}') -print(f'GPU sim time: {gpu:>10.4f} s') -print(f'G4 CPU time: {g4_cpu:>10.2f} s') -print(f'G4 wall time: {g4_wall:>10.2f} s') -print() -print(f'Speedup (CPU): {g4_cpu/gpu:>10.0f}x') -print(f'Speedup (wall): {g4_wall/gpu:>10.0f}x') -print() -print(f'GPU rate: {nphotons/gpu/1e6:>10.1f} M photons/s') -print(f'G4 rate: {nphotons/g4_cpu/1e3:>10.1f} k photons/s') -print() -print(f'GPU hits: {gpu_hits:>10}') -print(f'G4 hits: {g4_hits:>10}') -print(f'Hit diff: {hit_diff:>+9.1f}%') -" +python3 optiphy/ana/print_speedup.py \ + "$GPU_TIME" "$G4_CPU" "$G4_WALL" "$NPHOTONS" "$GPU_HITS" "$G4_HITS" rm -f "$LOGFILE" diff --git a/optiphy/ana/print_speedup.py b/optiphy/ana/print_speedup.py new file mode 100644 index 000000000..fd039db3c --- /dev/null +++ b/optiphy/ana/print_speedup.py @@ -0,0 +1,41 @@ +""" +print_speedup.py +================ +Format and print the GPU vs G4 speedup summary used by +benchmarks/benchmark_apex.sh after a single benchmark run. + +Usage: + python3 print_speedup.py +""" + +import sys + + +def main(argv): + gpu = float(argv[1]) + g4_cpu = float(argv[2]) + g4_wall = float(argv[3]) + nphotons = int(argv[4]) + gpu_hits = int(argv[5]) + g4_hits = int(argv[6]) + hit_diff = (gpu_hits - g4_hits) / g4_hits * 100 if g4_hits > 0 else 0 + + print() + print(f'Photons: {nphotons:>10,}') + print(f'GPU sim time: {gpu:>10.4f} s') + print(f'G4 CPU time: {g4_cpu:>10.2f} s') + print(f'G4 wall time: {g4_wall:>10.2f} s') + print() + print(f'Speedup (CPU): {g4_cpu/gpu:>10.0f}x') + print(f'Speedup (wall): {g4_wall/gpu:>10.0f}x') + print() + print(f'GPU rate: {nphotons/gpu/1e6:>10.1f} M photons/s') + print(f'G4 rate: {nphotons/g4_cpu/1e3:>10.1f} k photons/s') + print() + print(f'GPU hits: {gpu_hits:>10}') + print(f'G4 hits: {g4_hits:>10}') + print(f'Hit diff: {hit_diff:>+9.1f}%') + + +if __name__ == "__main__": + main(sys.argv)