diff --git a/benchmarks/benchmark_apex.sh b/benchmarks/benchmark_apex.sh new file mode 100755 index 000000000..cab100fde --- /dev/null +++ b/benchmarks/benchmark_apex.sh @@ -0,0 +1,55 @@ +#!/bin/bash +# benchmark_apex.sh — Measure GPU vs G4 speedup on apex.gdml +# +# Usage: +# ./benchmarks/benchmark_apex.sh + +GDML="tests/geom/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 optiphy/ana/print_speedup.py \ + "$GPU_TIME" "$G4_CPU" "$G4_WALL" "$NPHOTONS" "$GPU_HITS" "$G4_HITS" + +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 diff --git a/optiphy/ana/compare_gpu_g4.py b/optiphy/ana/compare_gpu_g4.py new file mode 100755 index 000000000..87fe210e6 --- /dev/null +++ b/optiphy/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() 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() 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) diff --git a/optiphy/ana/run_and_compare.py b/optiphy/ana/run_and_compare.py new file mode 100755 index 000000000..1cd599cd2 --- /dev/null +++ b/optiphy/ana/run_and_compare.py @@ -0,0 +1,308 @@ +#!/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).""" + 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: + 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 + 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, 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="", g4_full=None, g4_raw_shape=None): + 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() + + # 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', 'step_count.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_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 + 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, + g4_full=g4_full, g4_raw_shape=g4_raw_shape) + + +if __name__ == "__main__": + main() diff --git a/src/GPURaytrace.cpp b/src/GPURaytrace.cpp index c97c7ec11..f0590db8d 100644 --- a/src/GPURaytrace.cpp +++ b/src/GPURaytrace.cpp @@ -13,6 +13,7 @@ #include "sysrap/OPTICKS_LOG.hh" #include "GPURaytrace.h" +#include "config.h" #include "G4RunManager.hh" #include "G4RunManagerFactory.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..3ea2e4a5d 100644 --- a/src/GPURaytrace.h +++ b/src/GPURaytrace.h @@ -1,3 +1,4 @@ +#include #include #include #include @@ -47,7 +48,9 @@ namespace { G4Mutex genstep_mutex = G4MUTEX_INITIALIZER; -} +G4Mutex g4hits_mutex = G4MUTEX_INITIALIZER; +std::vector> g4_accumulated_hits; +} // namespace bool IsSubtractionSolid(G4VSolid *solid) { @@ -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..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,93 +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"] ); - } - 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 1fc5c838d..28e26416c 100644 --- a/src/config.h +++ b/src/config.h @@ -7,29 +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"); + public: + Config(std::string config_name = "dev"); - static std::string PtxPath(const std::string &ptx_name = "CSGOptiX7.ptx"); + static std::string PtxPath(const std::string &ptx_name = "CSGOptiX7.ptx"); - /// A unique name associated with this Config - std::string name; + /// A unique name associated with this Config + std::string name; - storch torch; + 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 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 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/tests/geom/apex.gdml b/tests/geom/apex.gdml new file mode 100644 index 000000000..49d36656b --- /dev/null +++ b/tests/geom/apex.gdml @@ -0,0 +1,379 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +