diff --git a/sbnd/sbnd_abhat/SUMMARY.txt b/sbnd/sbnd_abhat/SUMMARY.txt index 10f9349..8b9aa68 100644 --- a/sbnd/sbnd_abhat/SUMMARY.txt +++ b/sbnd/sbnd_abhat/SUMMARY.txt @@ -1,5 +1,9 @@ SBND WCT 0.36 SCECorrection validation -- 50 crossing-muon events -points: 228404 (East 70503, West 157901) -per-point residual reco-map: rms E=1.83e-04 W=2.23e-04 cm, max 8.42e-04 cm -pooled mean|Dx|: East=0.4021 West=0.5338 W/E=1.327 +points: 226721 paired by q (East 69913, West 156808); 226721/228762 matched, 2041 dropped (ambiguous q) +=== Per-point reco-vs-map closure === +Dx rms (cm): E = 1.83e-04 W = 2.23e-04 max = 5.05e-04 +Dy rms (cm): E = 2.63e-04 W = 1.65e-04 max = 9.91e-04 [Stage 3a] +Dz rms (cm): E = 3.35e-04 W = 2.99e-04 max = 9.50e-04 [Stage 3a] +=> reproduces the TH3 backward map to interpolation precision in all three components +pooled mean|Dx|: East=0.4026 West=0.5340 W/E=1.326 references: map vol-avg 1.276, 0.33-era reco 1.271 diff --git a/sbnd/sbnd_abhat/clus.jsonnet b/sbnd/sbnd_abhat/clus.jsonnet index 73dc82e..db21c98 100644 --- a/sbnd/sbnd_abhat/clus.jsonnet +++ b/sbnd/sbnd_abhat/clus.jsonnet @@ -21,9 +21,17 @@ local LeventNo = std.parseInt(initial_eventNo); local common_coords = ["x", "y", "z"]; local common_corr_coords = ["x_t0cor", "y", "z"]; -local common_sce_coords = ["x_sce", "y", "z"]; +local common_sce_coords = ["x_sce", "y_sce", "z_sce"]; +local sce_field = { + type: "SCEFieldTH3", + name: "sbnd_dualmap", + data: { + sce_map_file: "/cvmfs/sbnd.opensciencegrid.org/products/sbnd/sbnd_data/v01_42_00/SCEoffsets/SCEoffsets_SBND_E500_dualmap_CV_voxelTH3.root", + }, +}; + local dvm = { overall: { FV_xmin: -202.5 * wc.cm, @@ -51,7 +59,7 @@ local dvm = { FV_xmax: -0.45 * wc.cm, FV_xmin_margin: 2 * wc.cm, FV_xmax_margin: 2 * wc.cm, - sce_map_file: "/cvmfs/sbnd.opensciencegrid.org/products/sbnd/sbnd_data/v01_42_00/SCEoffsets/SCEoffsets_SBND_E500_dualmap_CV_voxelTH3.root", + sce_field: wc.tn(sce_field), }, a1f0pA: $.a0f0pA + { FV_xmin: 0.45 * wc.mm, @@ -81,7 +89,7 @@ local pctransforms(dv) = { type: "PCTransformSet", name: dv.name, data: { detector_volumes: wc.tn(dv) }, - uses: [dv] + uses: [dv, sce_field] }; @@ -343,7 +351,15 @@ local clus_all_apa ( detector: "sbnd", algorithm: "sce", pcname: "3d", - coords: common_sce_coords, // ["x_sce","y","z"] + coords: ["x_sce", "y", "z"], // original y,z kept for (y,z,q) pairing + individual: false + }, + { + name: "sce3d", // full 3D SCE-corrected coords (transverse closure) + detector: "sbnd", + algorithm: "sce3d", + pcname: "3d", + coords: ["x_sce", "y_sce", "z_sce"], individual: false } ], diff --git a/sbnd/sbnd_abhat/make_transverse_profiles.py b/sbnd/sbnd_abhat/make_transverse_profiles.py new file mode 100644 index 0000000..9c66314 --- /dev/null +++ b/sbnd/sbnd_abhat/make_transverse_profiles.py @@ -0,0 +1,266 @@ +#!/usr/bin/env python3 +""" +make_transverse_profiles.py -- Stage 3 (map-level), SBND SCE migration. + +Map-level SCE spatial-offset PROFILES vs drift distance, in Lane Kashur's +technote convention: + + * FORWARD displacement (Delta = reco - true) -> read TrueFwd_Displacement_* + * horizontal axis = DRIFT DISTANCE (0 = anode, ~200 cm = cathode), not signed x + * West = orange, East = blue (matches technote Fig 4) + +Produces one figure (08_transverse_profiles.png) with three panels: + + (a) Delta x vs drift -- CROSS-CHECK. Should reproduce technote Fig 4 + (peak ~0.9 cm near mid-drift, West curve above East). If it does, the + geometry / drift-distance handling is validated and panels (b),(c) are + trustworthy. + (b) Delta y vs drift -- conditioned near top (y>0) and bottom (y<0) walls. + Expect top negative / bottom positive, |Delta y| growing toward the + cathode (large drift). This is the SBND analog of MicroBooNE's + Delta_y_top(x) / Delta_y_bottom(x) (NOTE-1018, Sec 6.1). + (c) Delta z vs drift -- conditioned near the upstream (z->0) and + downstream (z->~500) ends, where Delta z is non-negligible. + +This reads the map ONLY (no BEE / reco output), so it gives the map +EXPECTATION. The reco-vs-map closure (Stage 3a + plot 09) comes later. + +Run on the gpvm with the SCE env loaded (ROOT + matplotlib available): + + source $SCE_TOP/restore_sce_env_036_sce.sh + python make_transverse_profiles.py + # optional: SCE_MAP_FILE=/some/other/map.root python make_transverse_profiles.py +""" + +import os +import numpy as np +import ROOT +import matplotlib +matplotlib.use("Agg") +import matplotlib.pyplot as plt + +# ---------------------------------------------------------------------------- +# Configuration (everything you might want to tune lives here) +# ---------------------------------------------------------------------------- +MAP_FILE = os.environ.get( + "SCE_MAP_FILE", + "/cvmfs/sbnd.opensciencegrid.org/products/sbnd/sbnd_data/v01_42_00/" + "SCEoffsets/SCEoffsets_SBND_E500_dualmap_CV_voxelTH3.root", +) + +# Forward-displacement TH3 names. Lane's Fig 7/8 are the FORWARD map, so we +# match that. (TrueBkwd_* exist too; do NOT mix the two within SBND.) +FWD = { + ("x", "E"): "TrueFwd_Displacement_X_E", ("x", "W"): "TrueFwd_Displacement_X_W", + ("y", "E"): "TrueFwd_Displacement_Y_E", ("y", "W"): "TrueFwd_Displacement_Y_W", + ("z", "E"): "TrueFwd_Displacement_Z_E", ("z", "W"): "TrueFwd_Displacement_Z_W", +} + +# Geometry (cm). Cathode at x = 0, anodes at x = +/- ANODE_ABS_X. +# East TPC = x < 0 (APA 0), West TPC = x > 0 (APA 1). +ANODE_ABS_X = 200.0 +TPC_XRANGE = {"E": (-ANODE_ABS_X, 0.0), "W": (0.0, ANODE_ABS_X)} + +# Conditioning regions (cm). +TOP_Y_MIN = 150.0 # "near top wall" (Delta y) +BOT_Y_MAX = -150.0 # "near bottom wall" (Delta y) +UPSTREAM_Z_MAX = 60.0 # near z = 0 end (Delta z) +DOWNSTREAM_Z_MIN = 440.0 # near z ~ 500 end (Delta z) + +# Bulk region used for the Delta x cross-check (mimics Lane's edge-avoiding +# calibration selection so the cross-check tracks Fig 4 rather than the corners). +BULK_ABS_Y_MAX = 150.0 +BULK_Z_MIN, BULK_Z_MAX = 50.0, 450.0 + +# Drift-distance binning. 20 cm matches existing plot 06; drop to 10 to mirror +# Lane Fig 4 more finely once you trust the script. +DRIFT_BINS = np.arange(0.0, ANODE_ABS_X + 1e-6, 20.0) + +# Lane colors (colorblind-safe). +COL = {"W": "#E69F00", "E": "#0072B2"} # West orange, East blue +TPC_LABEL = {"W": "West", "E": "East"} + +OUTDIR = os.environ.get("SCE_PLOT_OUT", "validation_plots") +OUTNAME = "08_transverse_profiles.png" + + +# ---------------------------------------------------------------------------- +# Helpers +# ---------------------------------------------------------------------------- +def load_voxels(h): + """Flatten a TH3 to (x, y, z, value) arrays of bin centers / contents [cm].""" + ax, ay, az = h.GetXaxis(), h.GetYaxis(), h.GetZaxis() + nx, ny, nz = h.GetNbinsX(), h.GetNbinsY(), h.GetNbinsZ() + xc = np.fromiter((ax.GetBinCenter(i) for i in range(1, nx + 1)), float, nx) + yc = np.fromiter((ay.GetBinCenter(j) for j in range(1, ny + 1)), float, ny) + zc = np.fromiter((az.GetBinCenter(k) for k in range(1, nz + 1)), float, nz) + val = np.empty((nx, ny, nz), dtype=float) + for i in range(nx): + for j in range(ny): + for k in range(nz): + val[i, j, k] = h.GetBinContent(i + 1, j + 1, k + 1) + X, Y, Z = np.meshgrid(xc, yc, zc, indexing="ij") + return X.ravel(), Y.ravel(), Z.ravel(), val.ravel() + + +def drift_profile(x, val, mask): + """Median (+ MAD) of val vs drift distance = ANODE_ABS_X - |x|, within mask.""" + d = ANODE_ABS_X - np.abs(x) + d, v = d[mask], val[mask] + cen, med, mad = [], [], [] + for lo, hi in zip(DRIFT_BINS[:-1], DRIFT_BINS[1:]): + b = (d >= lo) & (d < hi) + if b.sum() < 5: # need a few voxels for a stable median + continue + m = float(np.median(v[b])) + cen.append(0.5 * (lo + hi)) + med.append(m) + mad.append(float(np.median(np.abs(v[b] - m)))) + return np.array(cen), np.array(med), np.array(mad) + + +def watermark(ax): + ax.text(0.03, 0.95, "SBND Simulation\nPreliminary", transform=ax.transAxes, + ha="left", va="top", color="0.45", fontsize=10) + + +# ---------------------------------------------------------------------------- +# Load +# ---------------------------------------------------------------------------- +if not os.path.exists(MAP_FILE): + raise SystemExit(f"[fatal] SCE map not found:\n {MAP_FILE}\n" + f"Set SCE_MAP_FILE to the correct path.") + +f = ROOT.TFile.Open(MAP_FILE) +if not f or f.IsZombie(): + raise SystemExit(f"[fatal] could not open {MAP_FILE}") + +data = {} # (comp, tpc) -> (x,y,z,val) +for (comp, tpc), name in FWD.items(): + h = f.Get(name) + if not h: + raise SystemExit(f"[fatal] histogram '{name}' not in file. " + f"Check names with: rootls -l {MAP_FILE}") + ax = h.GetXaxis() + print(f" {name:32s} nbins=({h.GetNbinsX()},{h.GetNbinsY()},{h.GetNbinsZ()})" + f" x:[{ax.GetXmin():.1f},{ax.GetXmax():.1f}] cm") + data[(comp, tpc)] = load_voxels(h) + +os.makedirs(OUTDIR, exist_ok=True) + +# ---------------------------------------------------------------------------- +# Plot +# ---------------------------------------------------------------------------- +fig, axes = plt.subplots(1, 3, figsize=(16.5, 5.2)) +fig.suptitle("SBND SCE dualmap -- forward spatial offsets vs drift distance " + "(map expectation)", fontsize=13) + +# (a) Delta x cross-check ---------------------------------------------------- +axx = axes[0] +for tpc in ("W", "E"): + x, y, z, v = data[("x", tpc)] + xlo, xhi = TPC_XRANGE[tpc] + bulk = ((x >= xlo) & (x <= xhi) & (np.abs(y) <= BULK_ABS_Y_MAX) + & (z >= BULK_Z_MIN) & (z <= BULK_Z_MAX)) + # plot |Delta x| to match Lane Fig 4 (which shows the offset as positive) + c, m, e = drift_profile(x, np.abs(v), bulk) + axx.errorbar(c, m, yerr=e, marker="s", ms=5, lw=1.6, capsize=2, + color=COL[tpc], label=f"{TPC_LABEL[tpc]} TPC") +axx.set_title(r"(a) $|\Delta x|$ vs drift — cross-check vs Fig 4") +axx.set_xlabel("Drift Distance [cm]") +axx.set_ylabel(r"Spatial Offset $|\Delta x|$ [cm]") +axx.set_xlim(0, ANODE_ABS_X) +axx.grid(alpha=0.3) +axx.legend() +watermark(axx) + +# (b) Delta y top / bottom --------------------------------------------------- +axy = axes[1] +for tpc in ("W", "E"): + x, y, z, v = data[("y", tpc)] + xlo, xhi = TPC_XRANGE[tpc] + inx = (x >= xlo) & (x <= xhi) + top = inx & (y >= TOP_Y_MIN) + bot = inx & (y <= BOT_Y_MAX) + c, m, e = drift_profile(x, v, top) + axy.errorbar(c, m, yerr=e, marker="^", ms=5, lw=1.6, capsize=2, ls="-", + color=COL[tpc], label=f"{TPC_LABEL[tpc]} top") + c, m, e = drift_profile(x, v, bot) + axy.errorbar(c, m, yerr=e, marker="v", ms=5, lw=1.6, capsize=2, ls="--", + color=COL[tpc], label=f"{TPC_LABEL[tpc]} bottom") +axy.axhline(0, color="0.6", lw=0.8) +axy.set_title(r"(b) $\Delta y$ vs drift — near top / bottom walls") +axy.set_xlabel("Drift Distance [cm]") +axy.set_ylabel(r"Spatial Offset $\Delta y$ [cm]") +axy.set_xlim(0, ANODE_ABS_X) +axy.grid(alpha=0.3) +axy.legend(fontsize=8) +watermark(axy) + +# (c) Delta z upstream / downstream ----------------------------------------- +axz = axes[2] +for tpc in ("W", "E"): + x, y, z, v = data[("z", tpc)] + xlo, xhi = TPC_XRANGE[tpc] + inx = (x >= xlo) & (x <= xhi) + up = inx & (z <= UPSTREAM_Z_MAX) + dn = inx & (z >= DOWNSTREAM_Z_MIN) + c, m, e = drift_profile(x, v, up) + axz.errorbar(c, m, yerr=e, marker="^", ms=5, lw=1.6, capsize=2, ls="-", + color=COL[tpc], label=f"{TPC_LABEL[tpc]} upstream") + c, m, e = drift_profile(x, v, dn) + axz.errorbar(c, m, yerr=e, marker="v", ms=5, lw=1.6, capsize=2, ls="--", + color=COL[tpc], label=f"{TPC_LABEL[tpc]} downstream") +axz.axhline(0, color="0.6", lw=0.8) +axz.set_title(r"(c) $\Delta z$ vs drift — near upstream / downstream ends") +axz.set_xlabel("Drift Distance [cm]") +axz.set_ylabel(r"Spatial Offset $\Delta z$ [cm]") +axz.set_xlim(0, ANODE_ABS_X) +axz.grid(alpha=0.3) +axz.legend(fontsize=8) +watermark(axz) + +fig.tight_layout(rect=(0, 0, 1, 0.96)) +outpath = os.path.join(OUTDIR, OUTNAME) +fig.savefig(outpath, dpi=140) +print(f"\n[ok] wrote {outpath}") + +# ---------------------------------------------------------------------------- +# Console summary (sanity numbers) +# ---------------------------------------------------------------------------- +''' +print("\nPeak |median| offset per component (over plotted drift bins):") +for comp, label in (("x", "|dx| bulk"), ("y", "dy top/bot"), ("z", "dz up/dn")): + for tpc in ("E", "W"): + x, y, z, v = data[(comp, tpc)] + xlo, xhi = TPC_XRANGE[tpc] + inx = (x >= xlo) & (x <= xhi) + if comp == "x": + mask = inx & (np.abs(y) <= BULK_ABS_Y_MAX) & (z >= BULK_Z_MIN) & (z <= BULK_Z_MAX) + vals = np.abs(v) + elif comp == "y": + mask = inx & ((y >= TOP_Y_MIN) | (y <= BOT_Y_MAX)) + vals = v + else: + mask = inx & ((z <= UPSTREAM_Z_MAX) | (z >= DOWNSTREAM_Z_MIN)) + vals = v + _, m, _ = drift_profile(x, vals, mask) + peak = np.max(np.abs(m)) if m.size else float("nan") + print(f" {label:12s} {TPC_LABEL[tpc]:5s}: {peak:6.3f} cm") +''' +print("\nPeak |median| offset per component / region (over plotted drift bins):") +def _peak(x, v, mask): + _, m, _ = drift_profile(x, v, mask) + return np.max(np.abs(m)) if m.size else float("nan") + +for tpc in ("E", "W"): + xlo, xhi = TPC_XRANGE[tpc] + x, y, z, v = data[("x", tpc)]; inx = (x >= xlo) & (x <= xhi) + bulk = inx & (np.abs(y) <= BULK_ABS_Y_MAX) & (z >= BULK_Z_MIN) & (z <= BULK_Z_MAX) + print(f" |dx| bulk {TPC_LABEL[tpc]:5s}: {_peak(x, np.abs(v), bulk):6.3f} cm") + x, y, z, v = data[("y", tpc)]; inx = (x >= xlo) & (x <= xhi) + print(f" dy top {TPC_LABEL[tpc]:5s}: {_peak(x, v, inx & (y >= TOP_Y_MIN)):6.3f} cm") + print(f" dy bottom {TPC_LABEL[tpc]:5s}: {_peak(x, v, inx & (y <= BOT_Y_MAX)):6.3f} cm") + x, y, z, v = data[("z", tpc)]; inx = (x >= xlo) & (x <= xhi) + print(f" dz upstream {TPC_LABEL[tpc]:5s}: {_peak(x, v, inx & (z <= UPSTREAM_Z_MAX)):6.3f} cm") + print(f" dz downstream {TPC_LABEL[tpc]:5s}: {_peak(x, v, inx & (z >= DOWNSTREAM_Z_MIN)):6.3f} cm") diff --git a/sbnd/sbnd_abhat/make_validation_3d.py b/sbnd/sbnd_abhat/make_validation_3d.py new file mode 100644 index 0000000..5241d96 --- /dev/null +++ b/sbnd/sbnd_abhat/make_validation_3d.py @@ -0,0 +1,165 @@ +#!/usr/bin/env python3 +""" +make_validation_3d.py -- Stage 3a closure, SBND SCE migration. + +Pairs the img and sce3d bee sets on q (SCE-invariant) to recover, per point: + pre-correction (x_t0, y, z) [img] and (x_sce, y_sce, z_sce) [sce3d]. +Then for x/y/z: + reco Delta = corrected - pre (what the WCT reco applied) + map Delta = SIGN * TrueBkwd(x_t0,y,z) (independent map lookup) + +Outputs: + 09_pairing_check.png reco-vs-map residual in X (q-pairing validity; ~um) + 10_transverse_residual.png reco-vs-map residual in Y and Z (~um => plumbing OK) + 11_transverse_closure.png reco Delta y/z profiles vs drift, map overlay + (Lane convention: drift distance, West orange/East blue) +""" +import zipfile, json, os +import numpy as np +import matplotlib; matplotlib.use("Agg") +import matplotlib.pyplot as plt +import ROOT + +SCE_TOP = os.environ["SCE_TOP"] +ZIP = os.path.join(SCE_TOP, "sce_test", "mabc-all-apa.zip") +OUT = os.path.join(SCE_TOP, "sce_test", "validation_plots"); os.makedirs(OUT, exist_ok=True) +MAP = ("/cvmfs/sbnd.opensciencegrid.org/products/sbnd/sbnd_data/v01_42_00/" + "SCEoffsets/SCEoffsets_SBND_E500_dualmap_CV_voxelTH3.root") +SIGN = -1.0 # matches SCEFieldTH3 sign convention +ANODE = 200.0 +COL = {"W": "#E69F00", "E": "#0072B2"} + +f = ROOT.TFile.Open(MAP) +H = {(c, t): f.Get(f"TrueBkwd_Displacement_{c.upper()}_{t}") + for c in ("x", "y", "z") for t in ("E", "W")} + +def clampax(a, v): + lo = a.GetBinCenter(1); hi = a.GetBinCenter(a.GetNbins()); pad = 1e-3 * a.GetBinWidth(1) + return lo + pad if v < lo + pad else (hi - pad if v > hi - pad else v) +def mapval(h, x, y, z): + return h.Interpolate(clampax(h.GetXaxis(), x), clampax(h.GetYaxis(), y), clampax(h.GetZaxis(), z)) + +z = zipfile.ZipFile(ZIP) +evs = sorted({n.split("/")[1] for n in z.namelist() if n.startswith("data/")}, key=int) + +X0=[]; Y0=[]; Z0=[]; TP=[] +DXr=[]; DYr=[]; DZr=[]; DXm=[]; DYm=[]; DZm=[] +n_img = n_sce = n_match = n_amb = 0 +for e in evs: + try: + img = json.load(z.open("data/%s/%s-img-global.json" % (e, e))) + sce = json.load(z.open("data/%s/%s-sce3d-global.json" % (e, e))) + except KeyError: + continue + n_img += len(img["q"]); n_sce += len(sce["q"]) + qmap = {} + for i, q in enumerate(img["q"]): + qmap.setdefault(q, []).append(i) + for j, q in enumerate(sce["q"]): + idx = qmap.get(q) + if not idx or len(idx) != 1: # ambiguous / unmatched -> drop + n_amb += 1; continue + i = idx[0] + x0, y0, z0 = img["x"][i], img["y"][i], img["z"][i] + xs, ys, zs = sce["x"][j], sce["y"][j], sce["z"][j] + t = "E" if x0 < 0 else "W" + X0.append(x0); Y0.append(y0); Z0.append(z0); TP.append(0 if x0 < 0 else 1) + DXr.append(xs - x0); DYr.append(ys - y0); DZr.append(zs - z0) + DXm.append(SIGN * mapval(H[("x", t)], x0, y0, z0)) + DYm.append(SIGN * mapval(H[("y", t)], x0, y0, z0)) + DZm.append(SIGN * mapval(H[("z", t)], x0, y0, z0)) + n_match += 1 + +X0=np.array(X0); Y0=np.array(Y0); Z0=np.array(Z0); TP=np.array(TP) +DXr=np.array(DXr); DYr=np.array(DYr); DZr=np.array(DZr) +DXm=np.array(DXm); DYm=np.array(DYm); DZm=np.array(DZm) +E = TP == 0; W = TP == 1 +rx, ry, rz = DXr - DXm, DYr - DYm, DZr - DZm # residuals [cm] + +print(f"img pts {n_img} sce3d pts {n_sce} matched {n_match} dropped(ambiguous) {n_amb}") +print(f"x-residual rms: E={rx[E].std()*1e4:.2f}um W={rx[W].std()*1e4:.2f}um " + f"max={np.abs(rx).max()*1e4:.2f}um <-- if ~um, q-pairing is VALID") +print(f"y-residual rms: E={ry[E].std()*1e4:.2f}um W={ry[W].std()*1e4:.2f}um max={np.abs(ry).max()*1e4:.2f}um") +print(f"z-residual rms: E={rz[E].std()*1e4:.2f}um W={rz[W].std()*1e4:.2f}um max={np.abs(rz).max()*1e4:.2f}um") + +# unified validation record (supersedes the Stage-1 SUMMARY) +nev = len(evs) +wx = np.abs(DXr); poolE, poolW = wx[E].mean(), wx[W].mean() +with open(os.path.join(OUT, "SUMMARY.txt"), "w") as fh: + fh.write("SBND WCT 0.36 SCECorrection validation -- %d crossing-muon events\n" % nev) + fh.write("points: %d paired by q (East %d, West %d); %d/%d matched, %d dropped (ambiguous q)\n" + % (n_match, E.sum(), W.sum(), n_match, n_sce, n_amb)) + fh.write("=== Per-point reco-vs-map closure ===\n") + fh.write("Dx rms (cm): E = %.2e W = %.2e max = %.2e\n" % (rx[E].std(), rx[W].std(), np.abs(rx).max())) + fh.write("Dy rms (cm): E = %.2e W = %.2e max = %.2e [Stage 3a]\n" % (ry[E].std(), ry[W].std(), np.abs(ry).max())) + fh.write("Dz rms (cm): E = %.2e W = %.2e max = %.2e [Stage 3a]\n" % (rz[E].std(), rz[W].std(), np.abs(rz).max())) + fh.write("=> reproduces the TH3 backward map to interpolation precision in all three components\n") + fh.write("pooled mean|Dx|: East=%.4f West=%.4f W/E=%.3f\n" % (poolE, poolW, poolW/poolE)) + fh.write("references: map vol-avg 1.276, 0.33-era reco 1.271\n") +print("wrote", os.path.join(OUT, "SUMMARY.txt")) + +def watermark(ax): + ax.text(0.03, 0.95, "SBND Simulation\nPreliminary", transform=ax.transAxes, + ha="left", va="top", color="0.45", fontsize=9) + +# --- 09: X residual (q-pairing validity) --- +fig, ax = plt.subplots(figsize=(7, 4.5)) +for sel, lab in [(E, "East (APA 0)"), (W, "West (APA 1)")]: + ax.hist(rx[sel]*1e4, bins=48, range=(-6, 6), histtype="step", lw=1.6, label=lab) +ax.set_yscale("log"); ax.set_ylim(bottom=0.5) +ax.set_xlabel(r"reco $\Delta x-$map $\Delta x$ [$\mu$m]"); ax.set_ylabel("points (log)") +ax.set_title("q-pairing check (img <-> sce3d): X closure") +ax.legend(loc="upper right") +ax.text(0.02, 0.96, "rms E=%.2f, W=%.2f um\nmax=%.2f um\nmatched=%d dropped=%d" + % (rx[E].std()*1e4, rx[W].std()*1e4, np.abs(rx).max()*1e4, n_match, n_amb), + transform=ax.transAxes, va="top", fontsize=9, + bbox=dict(boxstyle="round,pad=0.35", fc="white", ec="0.7")) +fig.tight_layout(); fig.savefig(OUT+"/09_pairing_check.png", dpi=130); plt.close(fig) + +# --- 10: transverse residual (Y, Z) --- +fig, axs = plt.subplots(1, 2, figsize=(12, 4.5)) +for ax, r, c in [(axs[0], ry, "y"), (axs[1], rz, "z")]: + for sel, lab in [(E, "East"), (W, "West")]: + ax.hist(r[sel]*1e4, bins=48, range=(-6, 6), histtype="step", lw=1.6, label=lab) + ax.set_yscale("log"); ax.set_ylim(bottom=0.5) + ax.set_xlabel(rf"reco $\Delta {c}-$map $\Delta {c}$ [$\mu$m]") + ax.set_title(rf"$\Delta {c}$ closure (reco vs map)") + ax.legend(loc="upper right") + ax.text(0.02, 0.96, "rms E=%.2f, W=%.2f um" % + (r[E].std()*1e4, r[W].std()*1e4), transform=ax.transAxes, + va="top", fontsize=9, bbox=dict(boxstyle="round,pad=0.3", fc="white", ec="0.7")) +axs[0].set_ylabel("points (log)") +fig.suptitle("SCE transverse reco-vs-map per-point residual", fontsize=12) +fig.tight_layout(); fig.savefig(OUT+"/10_transverse_residual.png", dpi=130); plt.close(fig) + +# --- 11: transverse closure profiles vs drift (Lane convention), map overlay --- +d = ANODE - np.abs(X0) +DB = np.arange(0.0, 200.001, 20.0) +def prof(arr, sel): + cen, med = [], [] + for lo, hi in zip(DB[:-1], DB[1:]): + s = sel & (d >= lo) & (d < hi) + if s.sum() < 20: continue + cen.append(0.5*(lo+hi)); med.append(np.median(arr[s])) + return np.array(cen), np.array(med) + +fig, (ay, az) = plt.subplots(1, 2, figsize=(13, 5)) +for ax, dr, dm, comp, regs in [ + (ay, DYr, DYm, "y", [("top", Y0 > 150, "^", "-"), ("bottom", Y0 < -150, "v", "--")]), + (az, DZr, DZm, "z", [("upstream", Z0 < 60, "^", "-"), ("downstream", Z0 > 440, "v", "--")])]: + for tpc, sel0 in [("W", W), ("E", E)]: + for reg, rmask, mk, ls in regs: + sel = sel0 & rmask + c, mr = prof(dr, sel); _, mm = prof(dm, sel) + if c.size: + ax.plot(c, mr, ls=ls, marker=mk, ms=5, color=COL[tpc], + label=f"{tpc} {reg} (reco)") + ax.plot(c, mm, ls=":", color=COL[tpc], alpha=0.5, lw=2.5) # map overlay + ax.axhline(0, color="0.6", lw=0.8) + ax.set_xlabel("Drift Distance [cm]"); ax.set_xlim(0, 200) + ax.set_ylabel(rf"Spatial Offset $\Delta {comp}$ [cm]") + ax.set_title(rf"$\Delta {comp}$ closure: markers=reco, dotted=map"); ax.grid(alpha=0.3) + ax.legend(fontsize=8); watermark(ax) +fig.suptitle("SBND SCE transverse closure -- reco (markers) vs map (dotted) vs drift", fontsize=12) +fig.tight_layout(rect=(0, 0, 1, 0.96)); fig.savefig(OUT+"/11_transverse_closure.png", dpi=130); plt.close(fig) +print("wrote", sorted(n for n in os.listdir(OUT) if n.startswith(("09","10","11")))) diff --git a/sbnd/sbnd_abhat/make_validation_plots.py b/sbnd/sbnd_abhat/make_validation_plots_stage1.py similarity index 100% rename from sbnd/sbnd_abhat/make_validation_plots.py rename to sbnd/sbnd_abhat/make_validation_plots_stage1.py diff --git a/sbnd/sbnd_abhat/pics/08_transverse_profiles.png b/sbnd/sbnd_abhat/pics/08_transverse_profiles.png new file mode 100644 index 0000000..676b0ec Binary files /dev/null and b/sbnd/sbnd_abhat/pics/08_transverse_profiles.png differ diff --git a/sbnd/sbnd_abhat/pics/09_pairing_check.png b/sbnd/sbnd_abhat/pics/09_pairing_check.png new file mode 100644 index 0000000..ed15ceb Binary files /dev/null and b/sbnd/sbnd_abhat/pics/09_pairing_check.png differ diff --git a/sbnd/sbnd_abhat/pics/10_transverse_residual.png b/sbnd/sbnd_abhat/pics/10_transverse_residual.png new file mode 100644 index 0000000..486334e Binary files /dev/null and b/sbnd/sbnd_abhat/pics/10_transverse_residual.png differ diff --git a/sbnd/sbnd_abhat/pics/11_transverse_closure.png b/sbnd/sbnd_abhat/pics/11_transverse_closure.png new file mode 100644 index 0000000..ddb87ef Binary files /dev/null and b/sbnd/sbnd_abhat/pics/11_transverse_closure.png differ