Commit 71575ab5 authored by Vacaliuc, Bogdan's avatar Vacaliuc, Bogdan
Browse files

the plan that claude developed for creating flying-chi.py

parent dc889cc0
Loading
Loading
Loading
Loading
+494 −0
Original line number Diff line number Diff line
# Refactor CHI_scan.py using flying_z.py's camonitor + timestamp correlation

**Date:** 2026-04-15
**Branch:** bl4b-alignment-integration (tasking) → a new branch on `bl4b-scripts` (submodule, currently `ScriptScan/2026A`)
**Target output:** `bl4b-scripts/flying_chi.py`
**Status:** Plan — awaiting implementation

---

## 1. Goal

Produce a new standalone Spyder-runnable script at `bl4b-scripts/flying_chi.py`
that performs the **same alignment function** as `CHI_scan.py` (find the
`BL4B:Mot:chis` position that maximises σ_X, or minimises σ_Y, or maximises
integrated intensity of the beam on the 2D detector) using the **flying-scan
pattern** established in `bl4b-scripts/flying_z.py`:

- One continuous motor sweep from `startCHI` to `endCHI` instead of an N-point
  step scan.
- `camonitor` callbacks on motor `.RBV`, proton charge, and the detector
  observables — each callback appends a `(timestamp, value)` tuple to a
  per-PV list.
- Post-sweep reduction: `np.interp` every observable onto the `chis.RBV`
  timestamp grid, then run the same CHI_scan.py fit functions on the resulting
  `(chi, observable)` arrays.

Net effect: replace 11 discrete dwell points (each accumulating 1 mC of
proton charge) with one continuous sweep whose sample density is set by the
detector's update cadence — typically yielding 30–100× more (chi, σ) samples
in roughly the same wall-clock time.

---

## 2. Source material

| File | Role |
|---|---|
| `bl4b-scripts/flying_z.py` (159 lines, Apr 14 2026) | Template for the camonitor / interp / fit pattern on `BL4B:Mot:hs`. |
| `bl4b-scripts/CHI_scan.py` | Reference for PV list, slit setup, fit functions, ROI, fit-mode semantics, and the final move-to-fit behaviour. |
| `bl4b-scripts/CHI_measurement.py` | Reference for the FitImage (N-slice line-fit) approach — **deferred** from v1. |
| `tasking/analysis-chis-scan-integration.md` §4, §9, §13 | Confirms ADnED exposes 1D projection waveforms as `BL4B:Det:N1:Det1:XY:Stat:ProfileAverage{X,Y}_RBV`. This is the key insight that lets us `camonitor` shape observables without shipping the full 256×304 image. |

---

## 3. Key design differences from flying_z.py

### 3.1 The observable is a waveform, not a scalar

flying_z.py reduces the detector to a single scalar (`EventRate_RBV`) and a
single normalisation (`PCharge`). CHI alignment needs beam *shape*. To get
that without hauling a 256×304 2D array at every update we use ADnED's
pre-computed 1D projections:

- `BL4B:Det:N1:Det1:XY:Stat:ProfileAverageX_RBV` — 256-element X projection
- `BL4B:Det:N1:Det1:XY:Stat:ProfileAverageY_RBV` — 304-element Y projection

Each profile update is ~1–2 KB of payload. Our callback reduces the waveform
to scalars *in the callback* (weighted mean, σ, integrated area over the ROI)
and stores only the scalars. Memory is O(N_samples) × 4 floats per profile,
not O(N_samples × 300).

### 3.2 `EventUpdatePeriod = 10 ms` is unsafe for profile PVs

flying_z.py pins `EventUpdatePeriod = 10 ms` during the sweep. Q5 diagnostic
(`/SNS/users/6ov/BL4B/2026/04/14/baseline-event-update-period.txt`) shows:

- At the site default `EventUpdatePeriod = 500 ms`, all three observable PVs
  (`EventRate_RBV`, `ProfileAverageX_RBV`, `ProfileAverageY_RBV`) post
  reliably at ~2 updates/s matching the configured period.
- At `EventUpdatePeriod = 10 ms`, only ONE real profile update (beyond the
  initial connection value) landed in a 5 s window. The Stats plugin almost
  certainly has an internal floor that suppresses updates this fast — either
  a 60 Hz pulse-rate coupling or a minimum acquire-time check in ADnED.

**Consequence for flying_chi:** do not blindly override
`EventUpdatePeriod`. Expose it as an optional parameter
(`event_update_period_ms`, default `None` = don't touch). A safe default
known to post reliably is 100 ms or 500 ms — operator picks.

### 3.3 Sweep velocity must be sized to detector cadence, not the motor's VMAX

Sample density across the sweep is
`N = (endCHI − startCHI) / (VELO × EventUpdatePeriod)`.
For a defensible four-parameter Gaussian+bkg fit we want `N ≥ 30`:

```
VELO_max = (endCHI − startCHI) / (30 × EventUpdatePeriod)
```

Examples:

| sweep span | EventUpdatePeriod | max VELO for ≥30 samples |
|---|---|---|
| 3.6° (−1.8 → +1.8) | 500 ms | 0.24 °/s |
| 3.6° | 100 ms | 1.2 °/s |
| 3.6° | 50 ms | 2.4 °/s |

This gets chosen once per beamline by running the optional slew-rate probe
(§5.2) and paired with a confirmed-safe `EventUpdatePeriod`.

### 3.4 VELO / VMAX lifecycle must be save-restore (flying_z doesn't do this)

flying_z.py writes `HS_MaxVel = 2` and `HS_Vel = 2` unconditionally and never
restores them. flying_chi will:

- If `flight_velocity` parameter is `None` → leave `.VELO` / `.VMAX` alone
  (user has pre-configured them).
- If `flight_velocity` is set → snapshot current `.VELO` and `.VMAX`, write
  the flight value, and restore in a `try / finally` so Ctrl-C and fit
  failures don't leak the override.

Same save-restore pattern applies to `.BDST`, `EventUpdatePeriod`, and slit
`Y:Gap` values. Everything the script mutates must be listed in a single
restore block.

---

## 4. File layout

```
flying_chi.py
├── Module docstring (what, when, author, derived-from)
├── User-editable parameters block (match CHI_scan.py names where possible)
├── Fit functions (reused from CHI_scan.py)
├── Helper: sigma_of_profile(profile, mask, coords) → (mean, sigma, total)
├── Helper: measure_slew_rate(...)  — optional, not run by default
├── Main scan body (wrapped in try/finally)
│   ├── PV bindings
│   ├── Snapshot initial values of everything we will mutate
│   ├── Set slits, EventUpdatePeriod, BDST, VELO/VMAX (conditionally)
│   ├── Move chis → startCHI, wait idle
│   ├── Register camonitor callbacks (chi.RBV, PCharge, XPROF, YPROF)
│   ├── StartDiag; CHI.put(endCHI); wait idle; StopDiag
│   ├── Clear callbacks
│   └── finally: restore all snapshotted values
├── Post-process: interpolate observables onto tCHI
├── Fits per mode (FitXdist / FitYdist / FitIdist)
├── Plots (match CHI_scan.py visual style)
└── Optional: if MoveToXFit and fit-center is in range, CHI.put(fit-center)
```

---

## 5. Detailed implementation notes

### 5.1 User-editable parameters (mirror CHI_scan.py naming)

```python
# --- Sweep range ---
startCHI = -1.8
endCHI   = +1.8

# --- Optional flight overrides (None = don't touch) ---
flight_velocity       = None   # °/s; if None, use whatever .VELO is configured
flight_max_velocity   = None   # °/s; if None, leave .VMAX alone
event_update_period_ms = None  # ms; if None, leave EventUpdatePeriod alone
                               # (avoid 10 ms — see §3.2)

# --- Slit setup (same as CHI_scan.py) ---
SlitHeight = 0.25              # mm; Si, S1 Y-gap; S2 gets 2 × SlitHeight

# --- ROI for σ calculation (same as CHI_scan.py) ---
Xrange = [50, 200]             # X pixel ROI (out of 256)
Yrange = [120, 180]            # Y pixel ROI (out of 304)

# --- Fit mode selection (same as CHI_scan.py) ---
FitXdist   = True              # maximise σ_X
FitYdist   = False             # minimise σ_Y
FitIdist   = False             # maximise integrated intensity
MoveToXFit = True              # move chis to fit centre if in range

# --- Sanity-probe thresholds ---
min_profile_updates_per_sec = 1.0   # fail pre-flight if profile cadence is too slow
min_fit_samples             = 30    # warn if sweep produced too few samples to fit
```

FitImage mode is **deferred** to a later revision (would require camonitoring
`BL4B:Det:N1:Det1:XY:Array:ArrayData`, which is 311 KB/update and still needs
per-frame 2D analysis). The code will carry a comment pointing at
`CHI_scan.py` lines 176–224 so the next editor knows where to pick up.

### 5.2 Optional slew-rate probe

```python
def measure_slew_rate(motor_pv='BL4B:Mot:chis',
                      delta=1.0, trial_velocity=None, timeout_s=30):
    """
    One-shot utility to characterise the motor's actual slew rate. NOT called
    from the main flying_chi scan — invoke from the IPython prompt once to
    choose a sensible flight_velocity.

    Returns observed °/s in the steady middle 60% of the move and the maximum
    observed |ΔRBV| between consecutive camonitor samples (useful for sanity
    checking whether EventUpdatePeriod is short enough to resolve the sweep).

    If trial_velocity is None, uses the motor's currently configured .VELO.
    Always restores the starting .VELO on exit.
    """
```

Algorithm: snapshot `.VELO`, optionally set trial velocity, camonitor `.RBV`,
move `delta` relative to current position, wait idle, stop monitor. Fit a
linear slope to the middle 60% of `(t, v)` samples (skip accel/decel ends).

### 5.3 Pre-flight sanity probe (runs inside main scan)

Before registering the scan callbacks, do a short dry run: subscribe to the
two profile PVs for 1.0 s, count updates, raise a clear error if either is
below `min_profile_updates_per_sec`. This catches the "I set
EventUpdatePeriod too low and now the profiles aren't posting" failure mode
from Q5 before the sweep starts.

### 5.4 Callback design

Per-PV callback (mirrors `populate_HS_array` / `populate_RATE_array` /
`populate_PC_array` from flying_z.py):

```python
CHIarray   = []   # (ts, chi_rbv)
PCarray    = []   # (ts, pcharge)
XPROFarray = []   # (ts, mean_x, sigma_x, total_x_in_roi)
YPROFarray = []   # (ts, mean_y, sigma_y, total_y_in_roi)

def cb_chi(*a, **kw):
    CHIarray.append((kw['timestamp'], kw['value']))

def cb_pc(*a, **kw):
    PCarray.append((kw['timestamp'], kw['value']))

def cb_xprof(*a, **kw):
    prof = np.asarray(kw['value'])
    mean, sigma, total = sigma_of_profile(prof, Xmask, Xvec)
    XPROFarray.append((kw['timestamp'], mean, sigma, total))

def cb_yprof(*a, **kw):
    prof = np.asarray(kw['value'])
    mean, sigma, total = sigma_of_profile(prof, Ymask, Yvec)
    YPROFarray.append((kw['timestamp'], mean, sigma, total))
```

`sigma_of_profile` is the vectorised equivalent of CHI_scan.py lines 166–173
(weighted mean, variance, σ over the ROI, with an explicit guard against
zero-sum ROIs → `nan`).

### 5.5 PV registration

```python
# Motor
CHI        = PV('BL4B:Mot:chis')
CHIrbv     = PV('BL4B:Mot:chis.RBV')
CHIstat    = PV('BL4B:Mot:chis_Stat')   # 0 = idle (same convention as HSstat)
CHI_Vel    = PV('BL4B:Mot:chis.VELO')
CHI_MaxVel = PV('BL4B:Mot:chis.VMAX')
CHI_BDST   = PV('BL4B:Mot:chis.BDST')

# Detector / normalisation
StartDiag    = PV('BL4B:Det:N1:Start')
StopDiag     = PV('BL4B:Det:N1:Stop')
UpdatePeriod = PV('BL4B:Det:N1:EventUpdatePeriod')
PC           = PV('BL4B:Det:PCharge')

XPROF_PV = 'BL4B:Det:N1:Det1:XY:Stat:ProfileAverageX_RBV'
YPROF_PV = 'BL4B:Det:N1:Det1:XY:Stat:ProfileAverageY_RBV'

# Slits (same as CHI_scan.py)
SiY, S1Y, S2Y = PV('BL4B:Mot:si:Y:Gap'), PV('BL4B:Mot:s1:Y:Gap'), PV('BL4B:Mot:s2:Y:Gap')
SiYstat, S1Ystat, S2Ystat = (PV(f'BL4B:Mot:{m}:Y:Gap:Status') for m in ('si','s1','s2'))
```

### 5.6 Scan sequence (pseudo-code, wrapped in try/finally)

```python
initial = dict(
    EventUpdatePeriod = UpdatePeriod.get(),
    BDST              = CHI_BDST.get(),
    VELO              = CHI_Vel.get(),
    VMAX              = CHI_MaxVel.get(),
    SiY = SiY.get(),  S1Y = S1Y.get(),  S2Y = S2Y.get(),
)
try:
    StopDiag.put(1)
    if event_update_period_ms is not None: UpdatePeriod.put(event_update_period_ms)
    CHI_BDST.put(0)
    if flight_velocity     is not None: CHI_Vel.put(flight_velocity)
    if flight_max_velocity is not None: CHI_MaxVel.put(flight_max_velocity)

    # Slits
    SiY.put(SlitHeight); S1Y.put(SlitHeight); S2Y.put(2 * SlitHeight)
    wait_all_idle([SiYstat, S1Ystat, S2Ystat])

    # Move to start
    CHI.put(startCHI)
    while CHIstat.get() != 0: time.sleep(0.01)

    # Pre-flight profile-cadence sanity check
    _preflight_profile_cadence_check(XPROF_PV, YPROF_PV, window_s=1.0,
                                     min_hz=min_profile_updates_per_sec)

    # Register scan callbacks
    camonitor('BL4B:Mot:chis.RBV', callback=cb_chi)
    camonitor('BL4B:Det:PCharge',  callback=cb_pc)
    if FitXdist or FitIdist: camonitor(XPROF_PV, callback=cb_xprof)
    if FitYdist or FitIdist: camonitor(YPROF_PV, callback=cb_yprof)

    # Fly
    StartDiag.put(1)
    CHI.put(endCHI)
    while CHIstat.get() != 0: time.sleep(0.01)
    StopDiag.put(1)

    camonitor_clear('BL4B:Mot:chis.RBV')
    camonitor_clear('BL4B:Det:PCharge')
    camonitor_clear(XPROF_PV)
    camonitor_clear(YPROF_PV)
finally:
    # Unconditional restore
    UpdatePeriod.put(initial['EventUpdatePeriod'])
    CHI_BDST.put(initial['BDST'])
    if flight_velocity     is not None: CHI_Vel.put(initial['VELO'])
    if flight_max_velocity is not None: CHI_MaxVel.put(initial['VMAX'])
    SiY.put(initial['SiY']); S1Y.put(initial['S1Y']); S2Y.put(initial['S2Y'])
```

`BDST` is always snapshotted and always restored — we always override it to
zero for the flight (prevents the BDST sign-flip retry behaviour documented
in the parent CLAUDE.md "BDST/MRES backlash trap" section).

### 5.7 Post-process

```python
tCHI, vCHI = map(np.array, zip(*CHIarray))
tPC,  vPC  = map(np.array, zip(*PCarray))

if XPROFarray:
    tX = np.array([r[0] for r in XPROFarray])
    sX = np.array([r[2] for r in XPROFarray])
    iX = np.array([r[3] for r in XPROFarray])
    sX_on_chi = np.interp(tCHI, tX, sX)
    iX_on_chi = np.interp(tCHI, tX, iX)

if YPROFarray:
    tY = np.array([r[0] for r in YPROFarray])
    sY = np.array([r[2] for r in YPROFarray])
    sY_on_chi = np.interp(tCHI, tY, sY)

PC_on_chi = np.interp(tCHI, tPC, vPC)
```

For `FitIdist`, normalise by proton-charge increment so we get intensity per
unit beam current, not cumulative counts:

```python
# Use gradient in case tCHI is not exactly uniform — pyepics timestamps are
# server-side and may skew slightly from the update cadence
dPC_dt = np.gradient(PC_on_chi, tCHI)
intensity_normalised = iX_on_chi / np.maximum(dPC_dt, 1e-9)
```

### 5.8 Fits (reuse CHI_scan.py fit functions verbatim)

```python
def Gauss_ConstantBkg(x, a, b, c, d): return d + a * np.exp(-(x - b)**2 / (2 * c**2))
def GaussDip_ConstantBkg(x, a, b, c, d): return d - a * np.exp(-(x - b)**2 / (2 * c**2))
```

Per-mode fit routine mirrors CHI_scan.py lines 265–326:

| Mode | Input | Fit fn | p0 | Optimal parameter |
|---|---|---|---|---|
| `FitXdist` | `sX_on_chi` vs `vCHI` | `Gauss_ConstantBkg` | `[drop, argmax, 0.5, min]` | `px[1]` = chi at max σ_X |
| `FitYdist` | `sY_on_chi` vs `vCHI` | `GaussDip_ConstantBkg` | `[drop, argmin, 0.5, max]` | `py[1]` = chi at min σ_Y |
| `FitIdist` | `intensity_normalised` vs `vCHI` | `Gauss_ConstantBkg` | `[drop, argmax, 0.5, min]` | `pi[1]` = chi at max I |

Sigma width initial guess `0.5°` matches the sweep span sanity (±1.8° span,
typical rocking curve width << 1°).

### 5.9 Move-to-fit safety

```python
if MoveToXFit and FitXdist:
    fit_centre = px[1]
    if startCHI < fit_centre < endCHI:
        print('Moving chi to', np.round(fit_centre, 3))
        CHI.put(fit_centre)
    else:
        print('Bad fit (centre %.3f outside sweep range)' % fit_centre)
```

Only one `MoveTo*Fit` can be active at a time; flying_chi will enforce this
(unlike CHI_scan.py, which has a single `MoveToXFit` hard-coded to the X
fit). If multiple fit modes are enabled, move only to the one marked as
primary — or leave the move opt-in per mode (`MoveToXFit`, `MoveToYFit`,
`MoveToIFit`) for clarity.

### 5.10 Plotting

Match CHI_scan.py's visual style to avoid surprising the operator:

- One `plt.cla()` before drawing.
- Data as `'ok'` markers (size 2).
- Fit curve as `'-r'`.
- FWHM as a horizontal `hlines` at half-max between `b ± c`.
- Vertical `vlines` at `b` from `min` to `max`.
- Print `f'Chi for {mode}: {b:.3f} +/- {err:.3f}'` where `err = sqrt(cov[1,1])`.

Add one extra plot that flying_z.py pioneered: the σ_X(t) vs t raw trace,
with the chis(t) trace on a twin y-axis. Useful for diagnosing sweep
problems (did the motor stall? did the detector stop updating mid-sweep?).

---

## 6. Testing & validation

flying_chi runs against the live IOC — no unit-testable layer. Validation is
empirical, via the operator and the archiver:

1. **Smoke test with beam off.** Run the sweep; confirm no exceptions,
   confirm all three PVs gave back ≥30 samples each, confirm the restore
   block actually restored `VELO`, `VMAX`, `BDST`, `EventUpdatePeriod`, and
   slit gaps (use `caget` to check post-run state).
2. **Back-to-back comparison with CHI_scan.py** on the same sample/slit
   setup, beam on. Compare the fitted chi centre from each; they should
   agree within the fit error bar (typically 0.05–0.1°).
3. **Archiver cross-check** — after a successful flying_chi run, query the
   archiver for `BL4B:Mot:chis.RBV` and `BL4B:Det:N1:Det1:EventRate_RBV`
   across the sweep window using `setup/archiver-query.sh`. Confirm the
   camonitor samples we captured line up with the archived sample density
   (they should — camonitor sees the same CA monitor updates the archive
   engine does).
4. **Failure-mode test**: set `event_update_period_ms = 10`, run — expect
   the pre-flight cadence check to raise a clear error instead of
   completing an under-sampled sweep.

---

## 7. Open questions (for implementation time)

- **Primary fit when multiple modes enabled.** If the operator enables both
  `FitXdist` and `FitYdist`, and `MoveToXFit` is true, do we move to the X
  fit unconditionally, or to whichever has smaller fit error? Proposal:
  unconditional X fit wins, matching CHI_scan.py. Document this in the
  docstring so it's not a surprise.
- **Dedicated MoveTo flags per mode?** CHI_scan.py has only `MoveToXFit`.
  If operators ever want to land on the σ_Y minimum or the intensity peak
  instead, we'd need `MoveToYFit` / `MoveToIFit`. Defer to user request —
  do not speculatively add.
- **Restoring `.BDST` to zero.** If the original `.BDST` is non-zero (e.g.
  0.05°), restoring it is correct. But if the chis axis should not have a
  non-zero BDST (per the parent CLAUDE.md "air-padded rotary stages have
  negligible mechanical backlash" guidance), we might want to **print a
  warning** when we observe a non-zero initial BDST rather than silently
  restoring it. Propose: print the observed initial BDST value at scan
  start; leave restore as-is; operator can decide.
- **Velocity choice without the slew-rate probe.** First-run operator won't
  have measured the slew rate yet. Can we derive a conservative default
  from `chis.VMAX`? Proposal: if `flight_velocity is None`, simply leave
  `.VELO` alone (operator's last-used value). Don't try to be clever.

---

## 8. Out-of-scope for v1

- `FitImage` mode (N-slice Y-gaussian line-fit). Would need to camonitor
  `XY:Array:ArrayData` (311 KB per update). Defer.
- Integration with the FittingIOC / scantools framework (tracked separately
  in Redmine #4850, documented in
  `tasking/analysis-chis-scan-integration.md`). flying_chi v1 stays a
  standalone Spyder script, callable the way `CHI_scan.py` is today.
- Dual-Gaussian fits (for `Solid_theta_alignment_*_2gauss.py` patterns).
  Separate refactor.
- A unit-test scaffolding around `sigma_of_profile` — trivial pure function,
  add only if repeated elsewhere.

---

## 9. Estimated scope

~250 lines of Python, modelled closely on the two source scripts. Single
commit on a new branch of `bl4b-scripts`. No changes required to any other
repository in this session.

---

## 10. Post-implementation follow-ups to flying_z.py

The review of flying_z.py for this refactor exposed two gaps that should get
their own ticket (not blocking flying_chi):

- **No save/restore of `HS.VELO` and `HS.VMAX`.** Once flying_chi's
  try/finally pattern is validated, port it back to flying_z.py.
- **`EventUpdatePeriod = 10 ms` may be unsafe** depending on how
  `EventRate_RBV` is plumbed in the same Stats plugin. Q5 only probed the
  profile PVs. Worth a separate diagnostic on `EventRate_RBV` alone to
  confirm flying_z's 10 ms setting isn't silently under-sampling.