eic-opticks/tests/g4trap_compare.py

0001 #!/usr/bin/env python3
0002 """Compare G4 vs Opticks hit files: hit count + per-axis chi2."""
0003 import re, sys, numpy as np
0004
0005 PAT = re.compile(r'^\s*([\-\d.eE+]+)\s+([\-\d.eE+]+)\s+\(([^)]+)\)\s+\(([^)]+)\)')
0006
0007 def load(p):
0008     rows = []
0009     for line in open(p):
0010         m = PAT.match(line)
0011         if not m: continue
0012         rows.append((float(m.group(1)), float(m.group(2)),
0013                      *[float(x) for x in m.group(3).split(',')],
0014                      *[float(x) for x in m.group(4).split(',')]))
0015     return np.array(rows, float) if rows else np.empty((0,8))
0016
0017 def chi2_1d(a, b, bins):
0018     ha, _ = np.histogram(a, bins=bins)
0019     hb, _ = np.histogram(b, bins=bins)
0020     m = (ha + hb) > 0
0021     return float(np.sum((ha[m] - hb[m])**2 / (ha[m] + hb[m]))), int(m.sum())
0022
0023 g_path, o_path, label = sys.argv[1], sys.argv[2], sys.argv[3]
0024 tolerance_count = float(sys.argv[4]) if len(sys.argv) > 4 else 5.0
0025 tolerance_chi2  = float(sys.argv[5]) if len(sys.argv) > 5 else 5.0
0026
0027 g = load(g_path); o = load(o_path)
0028 n_g, n_o = len(g), len(o)
0029 rel = abs(n_o - n_g) / ((n_o + n_g) / 2) * 100 if n_o + n_g > 0 else 0
0030
0031 print(f"=== {label} ===")
0032 print(f"  G4 hits:      {n_g}")
0033 print(f"  Opticks hits: {n_o}")
0034 print(f"  rel diff:     {rel:.3f}%   (tol={tolerance_count}%)")
0035
0036 fail = []
0037 if rel > tolerance_count:
0038     fail.append(f"count rel-diff {rel:.2f}% > {tolerance_count}%")
0039 if n_g == 0 or n_o == 0:
0040     print("  no hits, skip distributions")
0041 else:
0042     for col, name, bins in [
0043         (2, 'x',  np.linspace(min(g[:,2].min(), o[:,2].min()), max(g[:,2].max(), o[:,2].max()), 31)),
0044         (3, 'y',  np.linspace(min(g[:,3].min(), o[:,3].min()), max(g[:,3].max(), o[:,3].max()), 31)),
0045         (5, 'dx', np.linspace(-1, 1, 21)),
0046         (6, 'dy', np.linspace(-1, 1, 21)),
0047         (7, 'dz', np.linspace(-1, 1, 21)),
0048     ]:
0049         chi2, ndf = chi2_1d(g[:, col], o[:, col], bins)
0050         ratio = chi2 / max(ndf, 1)
0051         marker = "FAIL" if ratio > tolerance_chi2 else "ok"
0052         print(f"  {name:>2}: chi2/ndf = {chi2:7.2f}/{ndf} = {ratio:5.2f}  {marker}")
0053         if ratio > tolerance_chi2:
0054             fail.append(f"{name} chi2/ndf {ratio:.2f}")
0055
0056 if fail:
0057     print(f"  FAILED: {', '.join(fail)}")
0058     sys.exit(1)
0059 print("  PASS")