diff --git a/.github/workflows/build-pull-request.yaml b/.github/workflows/build-pull-request.yaml index 2a8b63fe7..b04e065b2 100644 --- a/.github/workflows/build-pull-request.yaml +++ b/.github/workflows/build-pull-request.yaml @@ -79,3 +79,4 @@ jobs: docker run --rm --gpus 'device=1' ${{ env.IMAGE_NAME }}:${{ env.IMAGE_TAG }} tests/test_GPURaytrace.sh docker run --rm --gpus 'device=1' ${{ env.IMAGE_NAME }}:${{ env.IMAGE_TAG }} tests/test_GPUPhotonFileSource.sh docker run --rm --gpus 'device=1' ${{ env.IMAGE_NAME }}:${{ env.IMAGE_TAG }} tests/test_GPUPhotonSource_8x8SiPM.sh + docker run --rm --gpus 'device=1' ${{ env.IMAGE_NAME }}:${{ env.IMAGE_TAG }} tests/test_wavelength_shifting.sh diff --git a/config/wls_100k.json b/config/wls_100k.json new file mode 100644 index 000000000..26166e47b --- /dev/null +++ b/config/wls_100k.json @@ -0,0 +1,30 @@ +{ + "torch": { + "gentype": "TORCH", + "trackid": 0, + "matline": 0, + "numphoton": 100000, + + "pos": [0.0, 0.0, 0.0], + "time": 0.0, + + "mom": [0.0, 0.0, 1.0], + "weight": 0.0, + + "pol": [1.0, 0.0, 0.0], + "wavelength": 350.0, + + "zenith": [0.0, 1.0], + "azimuth": [0.0, 1.0], + + "radius": 0.0, + "distance": 0.0, + "mode": 255, + "type": "disc" + }, + + "event": { + "mode": "DebugLite", + "maxslot": 10000000 + } +} diff --git a/config/wls_scatter_viz.json b/config/wls_scatter_viz.json new file mode 100644 index 000000000..76bd9ece2 --- /dev/null +++ b/config/wls_scatter_viz.json @@ -0,0 +1,30 @@ +{ + "torch": { + "gentype": "TORCH", + "trackid": 0, + "matline": 0, + "numphoton": 10000, + + "pos": [0.0, 0.0, -25.0], + "time": 0.0, + + "mom": [0.0, 0.0, 1.0], + "weight": 0.0, + + "pol": [1.0, 0.0, 0.0], + "wavelength": 350.0, + + "zenith": [0.0, 0.3], + "azimuth": [0.0, 1.0], + + "radius": 0.0, + "distance": 0.0, + "mode": 255, + "type": "disc" + }, + + "event": { + "mode": "HitPhoton", + "maxslot": 100000 + } +} diff --git a/config/wls_slab.json b/config/wls_slab.json new file mode 100644 index 000000000..bfd412305 --- /dev/null +++ b/config/wls_slab.json @@ -0,0 +1,30 @@ +{ + "torch": { + "gentype": "TORCH", + "trackid": 0, + "matline": 0, + "numphoton": 100000, + + "pos": [0.0, 0.0, -50.0], + "time": 0.0, + + "mom": [0.0, 0.0, 1.0], + "weight": 0.0, + + "pol": [1.0, 0.0, 0.0], + "wavelength": 400.0, + + "zenith": [0.0, 0.0], + "azimuth": [0.0, 1.0], + + "radius": 0.0, + "distance": 0.0, + "mode": 255, + "type": "disc" + }, + + "event": { + "mode": "DebugLite", + "maxslot": 10000000 + } +} diff --git a/optiphy/ana/wls_diagnostic.py b/optiphy/ana/wls_diagnostic.py new file mode 100644 index 000000000..6975983cb --- /dev/null +++ b/optiphy/ana/wls_diagnostic.py @@ -0,0 +1,290 @@ +#!/usr/bin/env python3 +""" +wls_diagnostic.py : Detailed WLS wavelength distribution comparison GPU vs G4 +============================================================================== + +Compares wavelength and time distributions from GPU (opticks) and G4 hits, +performs KS test, and checks per-pass WLS conversion probability. + +Usage:: + + python ana/wls_diagnostic.py [--input-wavelength 350] +""" +import sys +import os +import argparse +import numpy as np + + +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 + + +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 ks_test_2sample(a, b): + """Two-sample Kolmogorov-Smirnov test (no scipy dependency).""" + na, nb = len(a), len(b) + a_sorted = np.sort(a) + b_sorted = np.sort(b) + all_vals = np.sort(np.concatenate([a_sorted, b_sorted])) + + cdf_a = np.searchsorted(a_sorted, all_vals, side='right') / na + cdf_b = np.searchsorted(b_sorted, all_vals, side='right') / nb + d_stat = np.max(np.abs(cdf_a - cdf_b)) + + # Approximate p-value (asymptotic) + n_eff = np.sqrt(na * nb / (na + nb)) + lam = (n_eff + 0.12 + 0.11 / n_eff) * d_stat + # Kolmogorov distribution approximation + if lam < 0.001: + p_val = 1.0 + else: + p_val = 2.0 * np.exp(-2.0 * lam * lam) + p_val = max(0.0, min(1.0, p_val)) + return d_stat, p_val + + +def print_header(title): + print() + print("=" * 74) + print(f" {title}") + print("=" * 74) + + +def print_hit_summary(gpu_hits, g4_hits, n_photons, input_wl): + print_header("HIT COUNT SUMMARY") + ng, nc = len(gpu_hits), len(g4_hits) + print(f" Input photons: {n_photons:>10d} (wavelength = {input_wl:.0f} nm)") + print(f" GPU hits: {ng:>10d} ({100*ng/n_photons:.2f}%)") + print(f" G4 hits: {nc:>10d} ({100*nc/n_photons:.2f}%)") + if ng > 0 and nc > 0: + ratio = nc / ng + # Significance + p_pool = (ng + nc) / (2 * n_photons) + se = np.sqrt(2 * p_pool * (1 - p_pool) / n_photons) + z = abs(ng/n_photons - nc/n_photons) / se if se > 0 else 0 + print(f" Ratio G4/GPU: {ratio:>10.4f}") + print(f" Z-score: {z:>10.1f} {'** SIGNIFICANT **' if z > 3 else '(within noise)'}") + print() + + +def print_wavelength_comparison(gpu_wl, g4_wl): + print_header("WAVELENGTH DISTRIBUTION COMPARISON") + + print(f"\n {'Statistic':<25s} {'GPU':>12s} {'G4':>12s} {'Diff':>12s}") + print(f" {'-'*25} {'-'*12} {'-'*12} {'-'*12}") + + for name, fn in [("Mean (nm)", np.mean), ("Std (nm)", np.std), + ("Median (nm)", np.median), ("Min (nm)", np.min), + ("Max (nm)", np.max)]: + gv, cv = fn(gpu_wl), fn(g4_wl) + print(f" {name:<25s} {gv:12.2f} {cv:12.2f} {cv-gv:12.2f}") + + # Percentiles + print() + for pct in [5, 25, 75, 95]: + gv = np.percentile(gpu_wl, pct) + cv = np.percentile(g4_wl, pct) + print(f" {'P%d (nm)' % pct:<25s} {gv:12.2f} {cv:12.2f} {cv-gv:12.2f}") + + # KS test + d_stat, p_val = ks_test_2sample(gpu_wl, g4_wl) + print(f"\n KS statistic: {d_stat:.6f}") + print(f" KS p-value: {p_val:.2e}") + if p_val < 0.01: + print(" ** Wavelength distributions are SIGNIFICANTLY DIFFERENT **") + else: + print(" Wavelength distributions are statistically compatible") + print() + + +def print_fine_histogram(gpu_wl, g4_wl, bin_width=10): + print_header(f"WAVELENGTH HISTOGRAM (bin={bin_width}nm)") + + lo = min(gpu_wl.min(), g4_wl.min()) + hi = max(gpu_wl.max(), g4_wl.max()) + bins = np.arange(np.floor(lo / bin_width) * bin_width, + np.ceil(hi / bin_width) * bin_width + bin_width, bin_width) + + gc, _ = np.histogram(gpu_wl, bins=bins) + cc, _ = np.histogram(g4_wl, bins=bins) + + # Normalize to density (per nm per photon) + gpu_dens = gc / (len(gpu_wl) * bin_width) + g4_dens = cc / (len(g4_wl) * bin_width) + + max_dens = max(gpu_dens.max(), g4_dens.max()) + bar_w = 25 + + print(f"\n {'Bin (nm)':<14s} {'GPU':>8s} {'G4':>8s} {'GPU/G4':>7s} GPU|G4") + print(f" {'-'*14} {'-'*8} {'-'*8} {'-'*7} {'-'*51}") + + for i in range(len(bins) - 1): + if gc[i] == 0 and cc[i] == 0: + continue + ratio_str = f"{gc[i]/cc[i]:.2f}" if cc[i] > 0 else " inf" + gb = "#" * int(gpu_dens[i] / max_dens * bar_w) if max_dens > 0 else "" + cb = "=" * int(g4_dens[i] / max_dens * bar_w) if max_dens > 0 else "" + print(f" {bins[i]:5.0f}-{bins[i+1]:5.0f} {gc[i]:8d} {cc[i]:8d} {ratio_str:>7s} {gb:<25s}|{cb:<25s}") + print() + + +def print_time_comparison(gpu_t, g4_t): + print_header("TIME DISTRIBUTION COMPARISON") + + print(f"\n {'Statistic':<25s} {'GPU':>12s} {'G4':>12s} {'Diff':>12s}") + print(f" {'-'*25} {'-'*12} {'-'*12} {'-'*12}") + + for name, fn in [("Mean (ns)", np.mean), ("Std (ns)", np.std), + ("Median (ns)", np.median), ("Min (ns)", np.min), + ("Max (ns)", np.max)]: + gv, cv = fn(gpu_t), fn(g4_t) + print(f" {name:<25s} {gv:12.4f} {cv:12.4f} {cv-gv:12.4f}") + + d_stat, p_val = ks_test_2sample(gpu_t, g4_t) + print(f"\n KS statistic: {d_stat:.6f}") + print(f" KS p-value: {p_val:.2e}") + if p_val < 0.01: + print(" ** Time distributions are SIGNIFICANTLY DIFFERENT **") + else: + print(" Time distributions are statistically compatible") + print() + + +def print_gpu_outcomes(photon): + print_header("GPU PHOTON OUTCOMES (all photons)") + q3 = photon[:, 3, :].view(np.uint32) + flag = q3[:, 0] & 0xFFFF + 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.2f}%") + print() + + +def print_position_comparison(gpu_hits, g4_hits): + print_header("SPATIAL DISTRIBUTION") + + gpu_pos = gpu_hits[:, 0, :3] + g4_pos = g4_hits[:, 0, :3] + gpu_r = np.sqrt(np.sum(gpu_pos**2, axis=1)) + g4_r = np.sqrt(np.sum(g4_pos**2, axis=1)) + + print(f"\n {'Statistic':<25s} {'GPU':>12s} {'G4':>12s} {'Diff':>12s}") + print(f" {'-'*25} {'-'*12} {'-'*12} {'-'*12}") + + for name, gv, cv in [ + ("Mean radius (mm)", gpu_r.mean(), g4_r.mean()), + ("Mean X (mm)", gpu_pos[:, 0].mean(), g4_pos[:, 0].mean()), + ("Mean Y (mm)", gpu_pos[:, 1].mean(), g4_pos[:, 1].mean()), + ("Mean Z (mm)", gpu_pos[:, 2].mean(), g4_pos[:, 2].mean()), + ]: + print(f" {name:<25s} {gv:12.3f} {cv:12.3f} {cv-gv:12.3f}") + print() + + +def print_energy_conservation(gpu_wl, g4_wl, input_wl): + print_header("ENERGY CONSERVATION CHECK") + gpu_viol = np.sum(gpu_wl < input_wl) + g4_viol = np.sum(g4_wl < input_wl) + print(f" Input wavelength: {input_wl:.0f} nm") + print(f" GPU hits with wl < input: {gpu_viol} / {len(gpu_wl)}") + print(f" G4 hits with wl < input: {g4_viol} / {len(g4_wl)}") + if gpu_viol == 0 and g4_viol == 0: + print(" ALL PASS: energy conservation satisfied") + else: + if gpu_viol > 0: + bad = gpu_wl[gpu_wl < input_wl] + print(f" GPU violations: min={bad.min():.1f}nm, mean={bad.mean():.1f}nm") + if g4_viol > 0: + bad = g4_wl[g4_wl < input_wl] + print(f" G4 violations: min={bad.min():.1f}nm, mean={bad.mean():.1f}nm") + print() + + +def main(): + parser = argparse.ArgumentParser(description=__doc__, + formatter_class=argparse.RawDescriptionHelpFormatter) + parser.add_argument("gpu_path", help="GPU opticks event folder") + parser.add_argument("g4_hits", help="G4 hits file (g4_hits.npy)") + parser.add_argument("--input-wavelength", type=float, default=350.0, + help="Input photon wavelength in nm (default: 350)") + parser.add_argument("--bin-width", type=float, default=5.0, + help="Histogram bin width in nm (default: 5)") + args = parser.parse_args() + + gpu_path = resolve_event_path(args.gpu_path) + hit_path = os.path.join(gpu_path, "hit.npy") + photon_path = os.path.join(gpu_path, "photon.npy") + + if not os.path.exists(photon_path): + 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) + + gpu_photon = np.load(photon_path) + gpu_hits = np.load(hit_path) if os.path.exists(hit_path) else np.zeros((0, 4, 4), dtype=np.float32) + g4_hits = np.load(args.g4_hits) + n_photons = len(gpu_photon) + + print(f"\n GPU event path: {gpu_path}") + print(f" G4 hits file: {args.g4_hits}") + + # Hit summary + print_hit_summary(gpu_hits, g4_hits, n_photons, args.input_wavelength) + + if len(gpu_hits) == 0 or len(g4_hits) == 0: + print(" Cannot compare — one side has zero hits.") + return + + gpu_wl = gpu_hits[:, 2, 3] + g4_wl = g4_hits[:, 2, 3] + gpu_t = gpu_hits[:, 0, 3] + g4_t = g4_hits[:, 0, 3] + + # GPU outcomes + print_gpu_outcomes(gpu_photon) + + # Energy conservation + print_energy_conservation(gpu_wl, g4_wl, args.input_wavelength) + + # Wavelength comparison + print_wavelength_comparison(gpu_wl, g4_wl) + print_fine_histogram(gpu_wl, g4_wl, bin_width=args.bin_width) + + # Time comparison + print_time_comparison(gpu_t, g4_t) + + # Spatial + print_position_comparison(gpu_hits, g4_hits) + + +if __name__ == "__main__": + main() diff --git a/tests/geom/wls_scatter_viz.gdml b/tests/geom/wls_scatter_viz.gdml new file mode 100644 index 000000000..db1dc872b --- /dev/null +++ b/tests/geom/wls_scatter_viz.gdml @@ -0,0 +1,111 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/tests/geom/wls_slab.gdml b/tests/geom/wls_slab.gdml new file mode 100644 index 000000000..50eb18af5 --- /dev/null +++ b/tests/geom/wls_slab.gdml @@ -0,0 +1,113 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/tests/test_wavelength_shifting.sh b/tests/test_wavelength_shifting.sh new file mode 100755 index 000000000..74278fec8 --- /dev/null +++ b/tests/test_wavelength_shifting.sh @@ -0,0 +1,257 @@ +#!/bin/bash +# +# test_wavelength_shifting.sh +# ============================ +# End-to-end test: GPU vs G4 wavelength shifting physics +# +# Fires 10000 UV photons (350nm) from outside a WLS sphere into a scattering +# medium. Compares GPU (opticks) and G4 hit wavelength distributions, WLS +# conversion rate, and arrival time distributions using chi-squared test. +# +# Geometry: tests/geom/wls_scatter_viz.gdml +# - WLS sphere r=10mm (absorbs UV, re-emits visible) +# - Scattering medium (Rayleigh, 10mm mean free path) +# - Detector shell r=30mm (100% efficiency) +# +# Usage: +# ./tests/test_wavelength_shifting.sh [seed] +# +# Exit code 0 = PASS, 1 = FAIL +# + +set -e + +SCRIPT_DIR="$(cd "$(dirname "${BASH_SOURCE[0]}")" && pwd)" +REPO_DIR="$(cd "$SCRIPT_DIR/.." && pwd)" +SEED=${1:-42} +NUMPHOTON=10000 +GEOM="$REPO_DIR/tests/geom/wls_scatter_viz.gdml" +CONFIG="wls_scatter_viz" + +source /opt/eic-opticks/eic-opticks-env.sh 2>/dev/null || true +export OPTICKS_MAX_BOUNCE=100 +export OPTICKS_EVENT_MODE=HitPhoton +export OPTICKS_MAX_SLOT=100000 + +echo "==============================================" +echo " WLS Test: GPU vs G4 Wavelength Shifting" +echo "==============================================" +echo " Geometry: $GEOM" +echo " Photons: $NUMPHOTON (350nm UV)" +echo " Seed: $SEED" +echo "" + +# --- GPU run --- +echo "[GPU] Running GPUPhotonSourceMinimal..." +GPU_OUT=$(/opt/eic-opticks/bin/GPUPhotonSourceMinimal \ + -g "$GEOM" -c "$CONFIG" -m "$REPO_DIR/tests/run.mac" -s "$SEED" 2>&1) +GPU_HITS=$(echo "$GPU_OUT" | grep "Opticks: NumHits" | head -1 | awk '{print $NF}') +echo "[GPU] Hits: $GPU_HITS" + +GPU_HIT_FILE="/tmp/MISSING_USER/opticks/GEOM/GEOM/GPUPhotonSourceMinimal/ALL0_no_opticks_event_name/A000/hit.npy" + +# --- G4 run --- +echo "[G4] Running StandAloneGeant4Validation..." +G4_OUT=$(/opt/eic-opticks/bin/StandAloneGeant4Validation \ + -g "$GEOM" -c "$CONFIG" -s "$SEED" 2>&1) +G4_HITS=$(echo "$G4_OUT" | grep "Total accumulated hits" | awk '{print $NF}') +echo "[G4] Hits: $G4_HITS" + +G4_HIT_FILE="g4_hits.npy" + +# --- Compare --- +echo "" +echo "[COMPARE] Analyzing wavelength and time distributions..." +echo "" + +python3 - "$GPU_HIT_FILE" "$G4_HIT_FILE" "$GPU_HITS" "$G4_HITS" << 'PYEOF' +import sys +import numpy as np + +gpu_hit_file = sys.argv[1] +g4_hit_file = sys.argv[2] +gpu_nhits = int(sys.argv[3]) +g4_nhits = int(sys.argv[4]) + +gpu = np.load(gpu_hit_file).reshape(-1, 4, 4) +g4 = np.load(g4_hit_file).reshape(-1, 4, 4) + +gpu_wl = gpu[:, 2, 3] +g4_wl = g4[:, 2, 3] +gpu_time = gpu[:, 0, 3] +g4_time = g4[:, 0, 3] + +PASS = True +ALPHA = 0.001 # significance level (tolerates minor ICDF interpolation difference) + + +def chi2_test(h_obs, h_exp, label): + """Chi-squared test for two histograms. Returns (chi2, ndf, p_value, pass).""" + # Scale expected to match observed total + scale = h_obs.sum() / h_exp.sum() if h_exp.sum() > 0 else 1.0 + h_exp_scaled = h_exp * scale + + # Only use bins with sufficient statistics (>5 expected) + mask = h_exp_scaled > 5 + if mask.sum() < 2: + print(f" {label}: Too few bins with sufficient stats") + return 0, 0, 1.0, True + + obs = h_obs[mask].astype(float) + exp = h_exp_scaled[mask].astype(float) + chi2 = np.sum((obs - exp) ** 2 / exp) + ndf = mask.sum() - 1 + + # p-value from chi2 distribution using Wilson-Hilferty approximation + if ndf > 0: + z = (chi2 / ndf) ** (1.0 / 3) - (1 - 2.0 / (9 * ndf)) + z /= np.sqrt(2.0 / (9 * ndf)) + # Approximate p-value from standard normal + p = 0.5 * (1.0 + math.erf(-z / np.sqrt(2))) + else: + p = 1.0 + + passed = p >= ALPHA + return chi2, ndf, p, passed + + +def ks_test(a, b): + """Two-sample Kolmogorov-Smirnov test.""" + a, b = np.sort(a), np.sort(b) + na, nb = len(a), len(b) + combined = np.concatenate([a, b]) + combined.sort() + cdf_a = np.searchsorted(a, combined, side='right') / na + cdf_b = np.searchsorted(b, combined, side='right') / nb + d = np.max(np.abs(cdf_a - cdf_b)) + en = np.sqrt(na * nb / (na + nb)) + p = min(np.exp(-2.0 * (en * d) ** 2) * 2.0, 1.0) + return d, p + + +# ------------------------------------------------------- +# Test 1: Hit count comparison +# ------------------------------------------------------- +print("=" * 55) +print(" TEST 1: Hit Count") +print("=" * 55) +print(f" GPU: {len(gpu)}") +print(f" G4: {len(g4)}") +import math +sigma = math.sqrt(len(gpu) + len(g4)) +z = abs(len(gpu) - len(g4)) / sigma if sigma > 0 else 0 +print(f" |Z| = {z:.1f}σ") +t1_pass = z < 5 +status = "PASS" if t1_pass else "FAIL" +print(f" Result: {status} (threshold: 5σ)") +PASS = PASS and t1_pass + + +# ------------------------------------------------------- +# Test 2: WLS conversion fraction +# ------------------------------------------------------- +print() +print("=" * 55) +print(" TEST 2: WLS Conversion Fraction") +print("=" * 55) +WLS_THRESHOLD = 380 # nm + +gpu_frac = np.mean(gpu_wl > WLS_THRESHOLD) +g4_frac = np.mean(g4_wl > WLS_THRESHOLD) +frac_diff = abs(gpu_frac - g4_frac) + +print(f" GPU shifted: {100*gpu_frac:.1f}%") +print(f" G4 shifted: {100*g4_frac:.1f}%") +print(f" |Difference|: {100*frac_diff:.2f}%") +t2_pass = frac_diff < 0.03 # 3% tolerance +status = "PASS" if t2_pass else "FAIL" +print(f" Result: {status} (threshold: 3%)") +PASS = PASS and t2_pass + + +# Pre-compute shifted/unshifted arrays +gpu_shifted = gpu_wl[gpu_wl > WLS_THRESHOLD] +g4_shifted = g4_wl[g4_wl > WLS_THRESHOLD] + +# ------------------------------------------------------- +# Test 3: Shifted wavelength spectrum (KS test) +# ------------------------------------------------------- +print() +print("=" * 55) +print(" TEST 3: Shifted Wavelength Spectrum (KS Test)") +print("=" * 55) + +if len(gpu_shifted) > 10 and len(g4_shifted) > 10: + d, p3 = ks_test(gpu_shifted, g4_shifted) + print(f" GPU shifted: N={len(gpu_shifted)}, mean={gpu_shifted.mean():.1f}nm") + print(f" G4 shifted: N={len(g4_shifted)}, mean={g4_shifted.mean():.1f}nm") + print(f" KS D={d:.6f} p={p3:.4f}") + t3_pass = p3 >= ALPHA +else: + print(" Too few shifted photons for KS test") + t3_pass = True + +status = "PASS" if t3_pass else "FAIL" +print(f" Result: {status} (threshold: p > {ALPHA})") +PASS = PASS and t3_pass + + +# ------------------------------------------------------- +# Test 4: Arrival time for shifted photons (KS test) +# ------------------------------------------------------- +print() +print("=" * 55) +print(" TEST 4: Shifted Photon Arrival Time (KS Test)") +print("=" * 55) + +# Compare shifted photon times — these include WLS exponential delay + transport +# With the G4 WLS time profile set to "exponential", distributions should match +gpu_shifted_t = gpu_time[gpu_wl > WLS_THRESHOLD] +g4_shifted_t = g4_time[g4_wl > WLS_THRESHOLD] + +print(f" GPU shifted: N={len(gpu_shifted_t)}, mean={gpu_shifted_t.mean():.3f}ns, std={gpu_shifted_t.std():.3f}ns") +print(f" G4 shifted: N={len(g4_shifted_t)}, mean={g4_shifted_t.mean():.3f}ns, std={g4_shifted_t.std():.3f}ns") +print(f" Std ratio: {gpu_shifted_t.std()/g4_shifted_t.std():.3f} (expect ~1.0)") + +if len(gpu_shifted_t) > 10 and len(g4_shifted_t) > 10: + d_t, p_t = ks_test(gpu_shifted_t, g4_shifted_t) + print(f" KS D={d_t:.6f} p={p_t:.4f}") + t4_pass = p_t >= ALPHA +else: + print(" Too few shifted photons for KS test") + t4_pass = True + +# Also check unshifted time (pure transport, no WLS delay) +gpu_unshifted_t = gpu_time[gpu_wl <= WLS_THRESHOLD] +g4_unshifted_t = g4_time[g4_wl <= WLS_THRESHOLD] +print(f" Unshifted time: GPU mean={gpu_unshifted_t.mean():.3f}ns G4 mean={g4_unshifted_t.mean():.3f}ns") + +status = "PASS" if t4_pass else "FAIL" +print(f" Result: {status} (KS p > {ALPHA})") +PASS = PASS and t4_pass + + +# ------------------------------------------------------- +# Summary +# ------------------------------------------------------- +print() +print("=" * 55) +print(" SUMMARY") +print("=" * 55) +tests = [ + ("Hit count", t1_pass), + ("WLS fraction", t2_pass), + ("Shifted wavelength KS", t3_pass), + ("Shifted time KS", t4_pass), +] +for name, passed in tests: + print(f" {name:>25s}: {'PASS' if passed else 'FAIL'}") + +print() +if PASS: + print(" *** ALL TESTS PASSED ***") + sys.exit(0) +else: + print(" *** SOME TESTS FAILED ***") + sys.exit(1) +PYEOF