optiphy/ana/wls_diagnostic.py

0001 #!/usr/bin/env python3
0002 """
0003 wls_diagnostic.py : Detailed WLS wavelength distribution comparison GPU vs G4
0004 ==============================================================================
0005
0006 Compares wavelength and time distributions from GPU (opticks) and G4 hits,
0007 performs KS test, and checks per-pass WLS conversion probability.
0008
0009 Usage::
0010
0011     python ana/wls_diagnostic.py <gpu_event_folder> <g4_hits.npy> [--input-wavelength 350]
0012 """
0013 import sys
0014 import os
0015 import argparse
0016 import numpy as np
0017
0018
0019 def resolve_event_path(path):
0020     if os.path.exists(os.path.join(path, "photon.npy")):
0021         return path
0022     a000 = os.path.join(path, "A000")
0023     if os.path.exists(os.path.join(a000, "photon.npy")):
0024         return a000
0025     if os.path.isdir(path):
0026         for d in sorted(os.listdir(path)):
0027             dp = os.path.join(path, d)
0028             if os.path.isdir(dp) and os.path.exists(os.path.join(dp, "photon.npy")):
0029                 return dp
0030     return path
0031
0032
0033 FLAG_ENUM = {
0034     0x0004: "TORCH", 0x0008: "BULK_ABSORB", 0x0010: "BULK_REEMIT",
0035     0x0020: "BULK_SCATTER", 0x0040: "SURFACE_DETECT", 0x0080: "SURFACE_ABSORB",
0036     0x0100: "SURFACE_DREFLECT", 0x0200: "SURFACE_SREFLECT",
0037     0x0400: "BOUNDARY_REFLECT", 0x0800: "BOUNDARY_TRANSMIT",
0038     0x1000: "NAN_ABORT", 0x2000: "EFFICIENCY_COLLECT", 0x8000: "MISS",
0039 }
0040
0041
0042 def ks_test_2sample(a, b):
0043     """Two-sample Kolmogorov-Smirnov test (no scipy dependency)."""
0044     na, nb = len(a), len(b)
0045     a_sorted = np.sort(a)
0046     b_sorted = np.sort(b)
0047     all_vals = np.sort(np.concatenate([a_sorted, b_sorted]))
0048
0049     cdf_a = np.searchsorted(a_sorted, all_vals, side='right') / na
0050     cdf_b = np.searchsorted(b_sorted, all_vals, side='right') / nb
0051     d_stat = np.max(np.abs(cdf_a - cdf_b))
0052
0053     # Approximate p-value (asymptotic)
0054     n_eff = np.sqrt(na * nb / (na + nb))
0055     lam = (n_eff + 0.12 + 0.11 / n_eff) * d_stat
0056     # Kolmogorov distribution approximation
0057     if lam < 0.001:
0058         p_val = 1.0
0059     else:
0060         p_val = 2.0 * np.exp(-2.0 * lam * lam)
0061         p_val = max(0.0, min(1.0, p_val))
0062     return d_stat, p_val
0063
0064
0065 def print_header(title):
0066     print()
0067     print("=" * 74)
0068     print(f"  {title}")
0069     print("=" * 74)
0070
0071
0072 def print_hit_summary(gpu_hits, g4_hits, n_photons, input_wl):
0073     print_header("HIT COUNT SUMMARY")
0074     ng, nc = len(gpu_hits), len(g4_hits)
0075     print(f"  Input photons:     {n_photons:>10d}   (wavelength = {input_wl:.0f} nm)")
0076     print(f"  GPU hits:          {ng:>10d}   ({100*ng/n_photons:.2f}%)")
0077     print(f"  G4  hits:          {nc:>10d}   ({100*nc/n_photons:.2f}%)")
0078     if ng > 0 and nc > 0:
0079         ratio = nc / ng
0080         # Significance
0081         p_pool = (ng + nc) / (2 * n_photons)
0082         se = np.sqrt(2 * p_pool * (1 - p_pool) / n_photons)
0083         z = abs(ng/n_photons - nc/n_photons) / se if se > 0 else 0
0084         print(f"  Ratio G4/GPU:      {ratio:>10.4f}")
0085         print(f"  Z-score:           {z:>10.1f}   {'** SIGNIFICANT **' if z > 3 else '(within noise)'}")
0086     print()
0087
0088
0089 def print_wavelength_comparison(gpu_wl, g4_wl):
0090     print_header("WAVELENGTH DISTRIBUTION COMPARISON")
0091
0092     print(f"\n  {'Statistic':<25s} {'GPU':>12s} {'G4':>12s} {'Diff':>12s}")
0093     print(f"  {'-'*25} {'-'*12} {'-'*12} {'-'*12}")
0094
0095     for name, fn in [("Mean (nm)", np.mean), ("Std (nm)", np.std),
0096                      ("Median (nm)", np.median), ("Min (nm)", np.min),
0097                      ("Max (nm)", np.max)]:
0098         gv, cv = fn(gpu_wl), fn(g4_wl)
0099         print(f"  {name:<25s} {gv:12.2f} {cv:12.2f} {cv-gv:12.2f}")
0100
0101     # Percentiles
0102     print()
0103     for pct in [5, 25, 75, 95]:
0104         gv = np.percentile(gpu_wl, pct)
0105         cv = np.percentile(g4_wl, pct)
0106         print(f"  {'P%d (nm)' % pct:<25s} {gv:12.2f} {cv:12.2f} {cv-gv:12.2f}")
0107
0108     # KS test
0109     d_stat, p_val = ks_test_2sample(gpu_wl, g4_wl)
0110     print(f"\n  KS statistic:      {d_stat:.6f}")
0111     print(f"  KS p-value:        {p_val:.2e}")
0112     if p_val < 0.01:
0113         print("  ** Wavelength distributions are SIGNIFICANTLY DIFFERENT **")
0114     else:
0115         print("  Wavelength distributions are statistically compatible")
0116     print()
0117
0118
0119 def print_fine_histogram(gpu_wl, g4_wl, bin_width=10):
0120     print_header(f"WAVELENGTH HISTOGRAM (bin={bin_width}nm)")
0121
0122     lo = min(gpu_wl.min(), g4_wl.min())
0123     hi = max(gpu_wl.max(), g4_wl.max())
0124     bins = np.arange(np.floor(lo / bin_width) * bin_width,
0125                      np.ceil(hi / bin_width) * bin_width + bin_width, bin_width)
0126
0127     gc, _ = np.histogram(gpu_wl, bins=bins)
0128     cc, _ = np.histogram(g4_wl, bins=bins)
0129
0130     # Normalize to density (per nm per photon)
0131     gpu_dens = gc / (len(gpu_wl) * bin_width)
0132     g4_dens = cc / (len(g4_wl) * bin_width)
0133
0134     max_dens = max(gpu_dens.max(), g4_dens.max())
0135     bar_w = 25
0136
0137     print(f"\n  {'Bin (nm)':<14s} {'GPU':>8s} {'G4':>8s} {'GPU/G4':>7s}  GPU|G4")
0138     print(f"  {'-'*14} {'-'*8} {'-'*8} {'-'*7}  {'-'*51}")
0139
0140     for i in range(len(bins) - 1):
0141         if gc[i] == 0 and cc[i] == 0:
0142             continue
0143         ratio_str = f"{gc[i]/cc[i]:.2f}" if cc[i] > 0 else "  inf"
0144         gb = "#" * int(gpu_dens[i] / max_dens * bar_w) if max_dens > 0 else ""
0145         cb = "=" * int(g4_dens[i] / max_dens * bar_w) if max_dens > 0 else ""
0146         print(f"  {bins[i]:5.0f}-{bins[i+1]:5.0f}   {gc[i]:8d} {cc[i]:8d} {ratio_str:>7s}  {gb:<25s}|{cb:<25s}")
0147     print()
0148
0149
0150 def print_time_comparison(gpu_t, g4_t):
0151     print_header("TIME DISTRIBUTION COMPARISON")
0152
0153     print(f"\n  {'Statistic':<25s} {'GPU':>12s} {'G4':>12s} {'Diff':>12s}")
0154     print(f"  {'-'*25} {'-'*12} {'-'*12} {'-'*12}")
0155
0156     for name, fn in [("Mean (ns)", np.mean), ("Std (ns)", np.std),
0157                      ("Median (ns)", np.median), ("Min (ns)", np.min),
0158                      ("Max (ns)", np.max)]:
0159         gv, cv = fn(gpu_t), fn(g4_t)
0160         print(f"  {name:<25s} {gv:12.4f} {cv:12.4f} {cv-gv:12.4f}")
0161
0162     d_stat, p_val = ks_test_2sample(gpu_t, g4_t)
0163     print(f"\n  KS statistic:      {d_stat:.6f}")
0164     print(f"  KS p-value:        {p_val:.2e}")
0165     if p_val < 0.01:
0166         print("  ** Time distributions are SIGNIFICANTLY DIFFERENT **")
0167     else:
0168         print("  Time distributions are statistically compatible")
0169     print()
0170
0171
0172 def print_gpu_outcomes(photon):
0173     print_header("GPU PHOTON OUTCOMES (all photons)")
0174     q3 = photon[:, 3, :].view(np.uint32)
0175     flag = q3[:, 0] & 0xFFFF
0176     n = len(flag)
0177     vals, counts = np.unique(flag, return_counts=True)
0178     order = np.argsort(-counts)
0179
0180     print(f"\n  {'Flag':<22s} {'Count':>8s} {'%':>7s}")
0181     print(f"  {'-'*22} {'-'*8} {'-'*7}")
0182     for idx in order:
0183         f = vals[idx]
0184         c = counts[idx]
0185         name = FLAG_ENUM.get(f, f"0x{f:04x}")
0186         print(f"  {name:<22s} {c:8d} {100*c/n:6.2f}%")
0187     print()
0188
0189
0190 def print_position_comparison(gpu_hits, g4_hits):
0191     print_header("SPATIAL DISTRIBUTION")
0192
0193     gpu_pos = gpu_hits[:, 0, :3]
0194     g4_pos = g4_hits[:, 0, :3]
0195     gpu_r = np.sqrt(np.sum(gpu_pos**2, axis=1))
0196     g4_r = np.sqrt(np.sum(g4_pos**2, axis=1))
0197
0198     print(f"\n  {'Statistic':<25s} {'GPU':>12s} {'G4':>12s} {'Diff':>12s}")
0199     print(f"  {'-'*25} {'-'*12} {'-'*12} {'-'*12}")
0200
0201     for name, gv, cv in [
0202         ("Mean radius (mm)", gpu_r.mean(), g4_r.mean()),
0203         ("Mean X (mm)", gpu_pos[:, 0].mean(), g4_pos[:, 0].mean()),
0204         ("Mean Y (mm)", gpu_pos[:, 1].mean(), g4_pos[:, 1].mean()),
0205         ("Mean Z (mm)", gpu_pos[:, 2].mean(), g4_pos[:, 2].mean()),
0206     ]:
0207         print(f"  {name:<25s} {gv:12.3f} {cv:12.3f} {cv-gv:12.3f}")
0208     print()
0209
0210
0211 def print_energy_conservation(gpu_wl, g4_wl, input_wl):
0212     print_header("ENERGY CONSERVATION CHECK")
0213     gpu_viol = np.sum(gpu_wl < input_wl)
0214     g4_viol = np.sum(g4_wl < input_wl)
0215     print(f"  Input wavelength:          {input_wl:.0f} nm")
0216     print(f"  GPU hits with wl < input:  {gpu_viol} / {len(gpu_wl)}")
0217     print(f"  G4  hits with wl < input:  {g4_viol} / {len(g4_wl)}")
0218     if gpu_viol == 0 and g4_viol == 0:
0219         print("  ALL PASS: energy conservation satisfied")
0220     else:
0221         if gpu_viol > 0:
0222             bad = gpu_wl[gpu_wl < input_wl]
0223             print(f"  GPU violations: min={bad.min():.1f}nm, mean={bad.mean():.1f}nm")
0224         if g4_viol > 0:
0225             bad = g4_wl[g4_wl < input_wl]
0226             print(f"  G4  violations: min={bad.min():.1f}nm, mean={bad.mean():.1f}nm")
0227     print()
0228
0229
0230 def main():
0231     parser = argparse.ArgumentParser(description=__doc__,
0232                                      formatter_class=argparse.RawDescriptionHelpFormatter)
0233     parser.add_argument("gpu_path", help="GPU opticks event folder")
0234     parser.add_argument("g4_hits", help="G4 hits file (g4_hits.npy)")
0235     parser.add_argument("--input-wavelength", type=float, default=350.0,
0236                         help="Input photon wavelength in nm (default: 350)")
0237     parser.add_argument("--bin-width", type=float, default=5.0,
0238                         help="Histogram bin width in nm (default: 5)")
0239     args = parser.parse_args()
0240
0241     gpu_path = resolve_event_path(args.gpu_path)
0242     hit_path = os.path.join(gpu_path, "hit.npy")
0243     photon_path = os.path.join(gpu_path, "photon.npy")
0244
0245     if not os.path.exists(photon_path):
0246         print(f"Error: photon.npy not found in {gpu_path}")
0247         sys.exit(1)
0248     if not os.path.exists(args.g4_hits):
0249         print(f"Error: {args.g4_hits} not found")
0250         sys.exit(1)
0251
0252     gpu_photon = np.load(photon_path)
0253     gpu_hits = np.load(hit_path) if os.path.exists(hit_path) else np.zeros((0, 4, 4), dtype=np.float32)
0254     g4_hits = np.load(args.g4_hits)
0255     n_photons = len(gpu_photon)
0256
0257     print(f"\n  GPU event path: {gpu_path}")
0258     print(f"  G4 hits file:   {args.g4_hits}")
0259
0260     # Hit summary
0261     print_hit_summary(gpu_hits, g4_hits, n_photons, args.input_wavelength)
0262
0263     if len(gpu_hits) == 0 or len(g4_hits) == 0:
0264         print("  Cannot compare — one side has zero hits.")
0265         return
0266
0267     gpu_wl = gpu_hits[:, 2, 3]
0268     g4_wl = g4_hits[:, 2, 3]
0269     gpu_t = gpu_hits[:, 0, 3]
0270     g4_t = g4_hits[:, 0, 3]
0271
0272     # GPU outcomes
0273     print_gpu_outcomes(gpu_photon)
0274
0275     # Energy conservation
0276     print_energy_conservation(gpu_wl, g4_wl, args.input_wavelength)
0277
0278     # Wavelength comparison
0279     print_wavelength_comparison(gpu_wl, g4_wl)
0280     print_fine_histogram(gpu_wl, g4_wl, bin_width=args.bin_width)
0281
0282     # Time comparison
0283     print_time_comparison(gpu_t, g4_t)
0284
0285     # Spatial
0286     print_position_comparison(gpu_hits, g4_hits)
0287
0288
0289 if __name__ == "__main__":
0290     main()