diff --git a/bin/debug/scan_artefact_check.py b/bin/debug/scan_artefact_check.py index efc5ba51..3921ff3b 100755 --- a/bin/debug/scan_artefact_check.py +++ b/bin/debug/scan_artefact_check.py @@ -285,82 +285,103 @@ def pick_fast_axis(signals: dict[str, np.ndarray], grid: tuple[int, int]) -> str # --------------------------------------------------------------------------- # -def check_end_of_line(signals, grid, monitored, fast_axis, req=None) -> None: +def check_end_of_line(signals, grid, monitored, fast_axis, req=None, + motors=frozenset()) -> None: lines, per_line = grid print(f"\n BUG 1 -- end-of-line intensity ({lines} lines x {per_line} points)") - if per_line < 4: + if per_line < 6: print(" too few points per line to judge; skipping") return exp = (req or {}).get("exp_time") exp = float(exp) if isinstance(exp, (int, float)) and exp > 0 else None - head = f" {'signal':<34} {'last/bulk':>10} {'first/bulk':>11}" + k = max(1, min(3, (per_line - 2) // 2)) + + head = f" {'signal':<30} {'first/near':>11} {'last/near':>10} {'trend':>7}" if exp: head += f" {'lost @start':>12}" print(head + " verdict") - print(" " + "-" * (86 if exp else 74)) + print(" " + "-" * (len(head) + 22)) deficits = [] for name, arr in sorted(signals.items()): if arr.size != lines * per_line: continue - dev = name.split("/")[0] 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: + if name.split("/")[0] in motors: + continue # a positioner, identified by having a setpoint + block = np.asarray(arr, dtype=float).reshape(lines, per_line) + if not motors and monotonicity(block) > 0.9: + continue # no positioner list to go on, so fall back to the shape + + # Each end is compared against its own neighbours, not the whole-line + # mean. A line with an intensity trend along it makes the whole-line + # mean a biased baseline and manufactures a deficit at the low end that + # has nothing to do with the first point -- scan 473 rises 20% along + # each line and reported three times the true figure because of it. + with np.errstate(divide="ignore", invalid="ignore"): + r_first = block[:, 0] / block[:, 1:1 + k].mean(axis=1) + r_last = block[:, -1] / block[:, -1 - k:-1].mean(axis=1) + r_first = r_first[np.isfinite(r_first)] + r_last = r_last[np.isfinite(r_last)] + if r_first.size < 3 or r_last.size < 3: 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: + first, last = float(np.mean(r_first)), float(np.mean(r_last)) + # Scatter over lines, so a channel with too few counts to say anything + # cannot contribute a number. mca8 sits at 2 counts and was polluting + # the summary with a deficit indistinguishable from Poisson noise. + sem = float(np.std(r_first)) / np.sqrt(r_first.size) + interior = block[:, 1:-1] + base = float(np.mean(interior[:, :k])) + trend = float(np.mean(interior[:, -k:])) / base if base else float("nan") + + real = (1.0 - first) > max(3.0 * sem, 0.02) + if abs(last - 1) < 0.02 and not real: verdict = "flat" - elif last < 0.98 and first >= 0.98: - verdict = "LAST POINT LOW <-- Bug 1" - elif last < 0.98 and first < 0.98: + elif real and last < 0.98: verdict = "both ends low (different mechanism)" - elif first < 0.98 and last >= 0.98: + elif real: verdict = "FIRST POINT LOW" + elif last < 0.98: + verdict = "LAST POINT LOW" else: verdict = "" - row = f" {name:<34} {last:>10.3f} {first:>11.3f}" + row = f" {name:<30} {first:>11.3f} {last:>10.3f} {trend:>7.3f}" if exp: - # A shortfall at the first point of every line, expressed as the - # exposure it is missing. If the cause is the shutter opening late - # against the first gate, this is a fixed number of milliseconds and - # is therefore the SAME across scans with different exposure times -- - # while first/bulk is not. That is the test. - lost = (1.0 - first) * exp - row += f" {lost * 1e3:>10.2f}ms" if first < 0.98 else f" {'-':>12}" - if first < 0.98: - deficits.append((name, lost)) + if real: + lost = (1.0 - first) * exp + deficits.append(lost) + row += f" {lost * 1e3:>10.2f}ms" + else: + row += f" {'-':>12}" print(row + f" {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.") + print(f" first/near = first point / mean of the next {k}; last/near likewise.") + print(" trend = interior level at the end of a line / at its start; far from") + print(" 1.00 means the line itself drifts, which is a separate effect.") if deficits: - vals = np.array([d for _, d in deficits]) + vals = np.array(deficits) print(f"\n lost exposure at the first point of every line:" f" {np.median(vals) * 1e3:.2f} ms" - f" (over {len(deficits)} channel(s))") - print(" If this is the shutter opening late against the first gate, the SAME") - print(" millisecond figure must come back from scans taken at other exposure") - print(f" times -- this one ran at {exp * 1e3:.0f} ms/point. A constant ms means") - print(" timing; a constant fraction means something proportional to the dose.") + f" (over {len(deficits)} channel(s), this scan ran at {exp * 1e3:.0f} ms/point)") + print(" A shutter opening late against the first gate costs a fixed number of") + print(" milliseconds, so this figure must repeat across scans taken at other") + print(" exposure times. A constant fraction instead would mean it scales with") + print(" dose and the shutter is innocent.") # End-of-*scan* behaviour is a separate question from end-of-*line*. print(f"\n Trend across the scan (last tenth of lines vs first tenth)") for name, arr in sorted(signals.items()): if arr.size != lines * per_line: continue - block = arr.reshape(lines, per_line) - if monotonicity(block) > 0.9: + if name.split("/")[0] in motors: continue - k = max(1, lines // 10) - head, tail = np.mean(block[:k]), np.mean(block[-k:]) - if not np.isfinite(head) or abs(head) < 1e-30: + block = np.asarray(arr, dtype=float).reshape(lines, per_line) + if not motors and monotonicity(block) > 0.9: continue - rel = tail / head + j = max(1, lines // 10) + first_tenth, last_tenth = np.mean(block[:j]), np.mean(block[-j:]) + if not np.isfinite(first_tenth) or abs(first_tenth) < 1e-30: + continue + rel = last_tenth / first_tenth flag = " <-- degrades towards the end" if rel < 0.95 else "" if abs(rel - 1) > 0.02: print(f" {name:<44} {rel:>7.3f}{flag}") @@ -1147,7 +1168,7 @@ def analyse(path: str, do_list: bool, force_lines=None, force_points=None, # line boundaries. check_line_motion_from_motor(signals, motors, meta, read_request_inputs(h5)) req = read_request_inputs(h5) - check_end_of_line(signals, grid, mon, None, req) + check_end_of_line(signals, grid, mon, None, req, motors) check_stationary_from_data(signals, grid, motors, req) if trace_lines: times = numeric_timestamps(h5)