diff --git a/bin/debug/README.md b/bin/debug/README.md new file mode 100644 index 00000000..4525c319 --- /dev/null +++ b/bin/debug/README.md @@ -0,0 +1,35 @@ +# Debug scripts + +Throwaway diagnostics for investigating beamline problems. **Not** part of the plugin: +nothing here is imported by `csaxs_bec`, and none of it needs BEC or a venv beyond the +libraries each script names. + +Kept on the `debug/scan-artefacts` branch rather than `main` so that a half-finished +diagnostic never reaches a production checkout. Promote anything that turns out to be +generally useful; delete the rest once the question it answered is settled. + +## scan_artefact_check.py + +Desk check for the July 2026 commissioning scan artefacts reported by Andreas Apseros +and Kazu Hirosawa. Read-only — opens BEC master files and writes nothing. + +```sh +python3 scan_artefact_check.py --root /sls/x12sa/data/p23080/raw/data \ + --scan 473 --scan 450 --scan 411 --scan 254 +``` + +Answers two questions without beam: + +- **end-of-line intensity** — is the last point of every line systematically weaker, and + is the *first* point affected too (which would mean a different mechanism)? +- **dropped line motion** — did the fast axis actually travel on every line? A line whose + span collapsed is motion that did not happen while the scan still reported success. + +Include a known-clean scan (254) so a null result distinguishes "this scan was fine" from +"the analysis is wrong". + +Run `--list` first if anything looks off; it dumps the file structure. + +> The motion check reads the fast-axis *readback*. If that axis is an open-loop stepper +> with no encoder, the readback is the commanded position and lost steps are invisible: +> a collapsed span still proves failure, but a healthy span proves nothing. diff --git a/bin/debug/scan_artefact_check.py b/bin/debug/scan_artefact_check.py new file mode 100755 index 00000000..8d6082ea --- /dev/null +++ b/bin/debug/scan_artefact_check.py @@ -0,0 +1,323 @@ +#!/usr/bin/env python3 +"""Desk check for the July 2026 commissioning scan artefacts. + +Read-only. Opens BEC master files and answers two questions without needing beam: + + Bug 1 -- is the last point of every line systematically weaker? + Quantifies it per monitored signal, and checks whether the *first* + point is affected too (which would point at a different mechanism). + + Bug 2 -- did the fast axis actually move on every line? + Compares the per-line travel of the fast-axis readback. A line where + motion silently failed shows up as a collapsed span while the scan + still reported success. + +Usage +----- + # what is in the file (run this first if anything looks odd) + python3 scan_artefact_check.py --list S00473_master.h5 + + # analyse one or more scans; a clean scan alongside a bad one is ideal + python3 scan_artefact_check.py S00473_master.h5 S00254_master.h5 + + # or by scan number, resolving the standard cSAXS layout + python3 scan_artefact_check.py --root /sls/x12sa/data/p23080/raw/data \ + --scan 473 --scan 450 --scan 254 + +Needs only h5py and numpy. + +Caveat worth knowing before trusting Bug 2's answer: if the fast axis is an +open-loop stepper with no encoder, its "readback" is really the commanded +position, and lost steps are invisible here. In that case a collapsed span +proves motion failed, but a healthy span does *not* prove it succeeded. +""" + +from __future__ import annotations + +import argparse +import glob +import os +import sys + +import h5py +import numpy as np + +DEVICES_PATH = "/entry/collection/devices" +READOUT_PATH = "/entry/collection/readout_groups" + + +# --------------------------------------------------------------------------- # +# file handling +# --------------------------------------------------------------------------- # + + +def resolve_scan(root: str, number: int) -> str | None: + """Find S_master.h5 under the standard S00000-00999/S/ layout.""" + pattern = os.path.join(root, "S*", f"S{number:05d}", f"S{number:05d}_master.h5") + hits = glob.glob(pattern) + if not hits: # some deployments omit the bucket directory + hits = glob.glob(os.path.join(root, f"S{number:05d}", f"S{number:05d}_master.h5")) + return hits[0] if hits else None + + +def walk_datasets(h5: h5py.File, top: str = "/") -> list[tuple[str, tuple]]: + """Every dataset in the file, as (path, shape). Follows soft links only.""" + found: list[tuple[str, tuple]] = [] + + def visit(name, obj): + if isinstance(obj, h5py.Dataset): + found.append((name, obj.shape)) + + h5[top].visititems(visit) + return found + + +def numeric_signals(h5: h5py.File) -> dict[str, np.ndarray]: + """Collect 1-D numeric per-point signals, keyed by 'device/signal'.""" + out: dict[str, np.ndarray] = {} + if DEVICES_PATH not in h5: + return out + for dev_name, dev in h5[DEVICES_PATH].items(): + if not isinstance(dev, h5py.Group): + continue + for sig_name, sig in dev.items(): + if not isinstance(sig, h5py.Group) or "value" not in sig: + continue + data = sig["value"] + if not isinstance(data, h5py.Dataset) or data.ndim != 1: + continue + if not np.issubdtype(data.dtype, np.number): + continue + arr = np.asarray(data[()], dtype=float) + if arr.size > 1: + out[f"{dev_name}/{sig_name}"] = arr + return out + + +def read_metadata(h5: h5py.File) -> dict: + """Pull the scan bookkeeping we need, tolerating layout differences.""" + meta: dict = {} + for path in ("/entry/collection/metadata", "/entry/collection"): + if path not in h5: + continue + node = h5[path] + for key in ( + "scan_name", + "scan_number", + "num_points", + "num_lines", + "frames_per_trigger", + ): + if key in node: + try: + val = node[key][()] + meta[key] = val.item() if hasattr(val, "item") else val + except Exception: # pylint: disable=broad-except + pass + for key, val in getattr(node, "attrs", {}).items(): + meta.setdefault(key, val) + return {k: (v.decode() if isinstance(v, bytes) else v) for k, v in meta.items()} + + +def monitored_names(h5: h5py.File) -> set[str]: + grp = f"{READOUT_PATH}/monitored" + return set(h5[grp].keys()) if grp in h5 else set() + + +# --------------------------------------------------------------------------- # +# geometry +# --------------------------------------------------------------------------- # + + +def infer_grid(n: int, meta: dict) -> tuple[int, int] | None: + """(num_lines, points_per_line) for n samples, from metadata or factoring.""" + lines = meta.get("num_lines") + frames = meta.get("frames_per_trigger") + if lines and frames and int(lines) * int(frames) == n: + return int(lines), int(frames) + if frames and n % int(frames) == 0: + return n // int(frames), int(frames) + if lines and n % int(lines) == 0: + return int(lines), n // int(lines) + return None + + +def monotonicity(block: np.ndarray) -> float: + """How steadily a signal ramps within a line: 1.0 = perfect ramp, 0 = noise. + + This is what separates a scanned axis from an intensity. Judging by + within-line *range* alone is not enough -- a noisy detector channel can have + a larger spread than the axis it is plotted against, and then gets mistaken + for the fast axis (which is exactly what the first version of this script + did). Correlation against the point index does not care about scale. + """ + idx = np.arange(block.shape[1], dtype=float) + scores = [] + for row in block: + if np.ptp(row) <= 0 or not np.all(np.isfinite(row)): + continue # a stalled line tells us nothing about what kind of signal this is + scores.append(abs(np.corrcoef(idx, row)[0, 1])) + return float(np.mean(scores)) if scores else 0.0 + + +def pick_fast_axis(signals: dict[str, np.ndarray], grid: tuple[int, int]) -> str | None: + """The fast axis is the signal that ramps most cleanly *within* each line.""" + lines, per_line = grid + best, best_score = None, 0.0 + for name, arr in signals.items(): + if arr.size != lines * per_line: + continue + block = arr.reshape(lines, per_line) + if float(np.median(np.ptp(block, axis=1))) <= 0: + continue # constant within lines -> slow axis or a setting + score = monotonicity(block) + if score > best_score: + best, best_score = name, score + return best if best_score > 0.9 else None + + +# --------------------------------------------------------------------------- # +# the two checks +# --------------------------------------------------------------------------- # + + +def check_end_of_line(signals, grid, monitored, fast_axis) -> None: + lines, per_line = grid + print(f"\n BUG 1 -- end-of-line intensity ({lines} lines x {per_line} points)") + if per_line < 4: + print(" too few points per line to judge; skipping") + return + print(f" {'signal':<34} {'last/bulk':>10} {'first/bulk':>11} verdict") + print(" " + "-" * 74) + for name, arr in sorted(signals.items()): + if arr.size != lines * per_line: + continue + dev = name.split("/")[0] + if monitored and dev not in monitored: + continue + if fast_axis and name == fast_axis: + continue + block = arr.reshape(lines, per_line) + if monotonicity(block) > 0.9: + continue # a position ramp, not an intensity + bulk = np.mean(block[:, 1:-1]) + if not np.isfinite(bulk) or abs(bulk) < 1e-30: + continue + last = np.mean(block[:, -1]) / bulk + first = np.mean(block[:, 0]) / bulk + # only remark on signals that plausibly carry intensity + if abs(last - 1) < 0.02 and abs(first - 1) < 0.02: + verdict = "flat" + elif last < 0.98 and first >= 0.98: + verdict = "LAST POINT LOW <-- Bug 1" + elif last < 0.98 and first < 0.98: + verdict = "both ends low (different mechanism)" + else: + verdict = "" + print(f" {name:<34} {last:>10.3f} {first:>11.3f} {verdict}") + print(" last/bulk = mean of final point / mean of interior points, over all lines.") + print(" A value near 1.00 is healthy; Bug 1 predicts a systematic shortfall.") + + +def check_line_motion(signals, grid, fast_axis) -> None: + lines, per_line = grid + print("\n BUG 2 -- did the fast axis move on every line?") + if fast_axis is None: + print(" could not identify a fast axis; re-run with --list and tell me the names") + return + block = signals[fast_axis].reshape(lines, per_line) + spans = np.ptp(block, axis=1) + median = float(np.median(spans)) + print(f" fast axis : {fast_axis}") + print(f" median line span : {median:.6g}") + if median <= 0: + print(" the fast axis does not move at all in this file -- check the signal choice") + return + ratio = spans / median + bad = np.where(ratio < 0.9)[0] + dead = np.where(ratio < 0.1)[0] + print(f" lines below 90% of median span : {len(bad)}") + print(f" lines below 10% (motion absent) : {len(dead)}") + if len(bad): + print(f" {'line':>6} {'span':>14} {'vs median':>10}") + for i in bad[:25]: + print(f" {i:>6} {spans[i]:>14.6g} {ratio[i]:>9.2f}x") + if len(bad) > 25: + print(f" ... and {len(bad) - 25} more") + print("\n Lines with a collapsed span are motion that did not happen while the") + print(" scan still reported success -- the signature we are looking for.") + else: + print(" all lines travelled a consistent distance; no dropped motion visible here.") + starts, ends = block[:, 0], block[:, -1] + print(f" line start scatter : {np.std(starts):.3g} line end scatter : {np.std(ends):.3g}") + print(" (large scatter at one end suggests the return move races the next line)") + + +# --------------------------------------------------------------------------- # + + +def analyse(path: str, do_list: bool) -> None: + print("=" * 80) + print(os.path.basename(path)) + print("=" * 80) + with h5py.File(path, "r") as h5: + meta = read_metadata(h5) + interesting = {k: meta[k] for k in + ("scan_name", "scan_number", "num_points", "num_lines", + "frames_per_trigger") if k in meta} + print(" metadata:", interesting if interesting else "(none of the expected keys found)") + + if do_list: + print("\n datasets:") + for name, shape in walk_datasets(h5): + print(f" {shape!s:>16} {name}") + return + + signals = numeric_signals(h5) + if not signals: + print(" no 1-D numeric signals found -- re-run with --list") + return + n = max(len(v) for v in signals.values()) + grid = infer_grid(n, meta) + if grid is None: + print(f" could not work out the line geometry for {n} points.") + print(" Pass --lines / --points, or send me the --list output.") + return + mon = monitored_names(h5) + fast = pick_fast_axis(signals, grid) + check_line_motion(signals, grid, fast) + check_end_of_line(signals, grid, mon, fast) + + +def main(argv=None) -> int: + ap = argparse.ArgumentParser(description=__doc__.splitlines()[0]) + ap.add_argument("files", nargs="*", help="master .h5 files") + ap.add_argument("--scan", type=int, action="append", default=[], help="scan number") + ap.add_argument("--root", default="/sls/x12sa/data/p23080/raw/data", + help="data root used with --scan") + ap.add_argument("--list", action="store_true", help="dump the file structure and stop") + ap.add_argument("--lines", type=int, help="override number of lines") + ap.add_argument("--points", type=int, help="override points per line") + args = ap.parse_args(argv) + + paths = list(args.files) + for num in args.scan: + found = resolve_scan(args.root, num) + if found: + paths.append(found) + else: + print(f"scan {num}: no master file under {args.root}", file=sys.stderr) + if not paths: + ap.error("give at least one file or --scan") + + for path in paths: + try: + analyse(path, args.list) + except Exception as exc: # pylint: disable=broad-except + print(f" FAILED on {path}: {type(exc).__name__}: {exc}") + print() + return 0 + + +if __name__ == "__main__": + sys.exit(main())