pxii_bec/bec_ipython_client/startup/post_startup.py

@@ -33,5 +33,6 @@ to setup the prompts.
         logger.success("cSAXS session loaded.")
 """
 
-# pylint: disable=invalid-name, unused-import, import-error, undefined-variable, unused-variable, unused-argument, no-name-in-module
-d, planner = init_beamline_environment()
+from pxii_bec.scripts import beamline_context as bl
+
+bl.initialise()

pxii_bec/bec_ipython_client/startup/pre_startup.py

@@ -4,11 +4,19 @@ is started. It can be used to add additional command line arguments.
 """
 
 import os
+import sys
+from pathlib import Path
 
 from bec_lib.service_config import ServiceConfig
 
 import pxii_bec
 
+_package_dir = Path(pxii_bec.__file__).resolve().parent
+_scripts_dir = _package_dir / "scripts"
+
+if str(_scripts_dir) not in sys.path:
+    sys.path.insert(0, str(_scripts_dir))
+
 
 def extend_command_line_args(parser):
     """

pxii_bec/device_configs/pxii-devices.yaml

@@ -65,12 +65,12 @@ aerotech:
       mount: 0
       work: 0
     x:
-      in: 0
+      in: 0.0
       out: -10
       safe: -50
     y:
       mount: -0.0
-      work: -0.031
+      work: -0.027
     z:
       mount: 0
       work: -0.01
@@ -257,6 +257,7 @@ bl_bright:
   onFailure: buffer
   readoutPriority: baseline
   userParameter:
+    move_method: signal
     'off': 0.0
     'on': 1.3
     tol: 0.01
@@ -302,7 +303,7 @@ bs_x:
   onFailure: buffer
   readoutPriority: baseline
   userParameter:
-    in: 0.3
+    in: 0.361
 bs_y:
   description: Beamstop Y
   deviceClass: ophyd_devices.EpicsMotor
@@ -314,7 +315,7 @@ bs_y:
   onFailure: buffer
   readoutPriority: baseline
   userParameter:
-    in: -3.35
+    in: 0.258
 bs_z:
   description: Beamstop Z
   deviceClass: ophyd_devices.EpicsMotor
@@ -444,7 +445,7 @@ coll_x:
   onFailure: buffer
   readoutPriority: baseline
   userParameter:
-    in: -1.275
+    in: -1.297
 coll_y:
   description: Collimator Y
   deviceClass: ophyd_devices.EpicsMotor
@@ -456,7 +457,7 @@ coll_y:
   onFailure: buffer
   readoutPriority: baseline
   userParameter:
-    in: 39.06
+    in: 39.002
     intermediate: 32
     out: 20.002
     park: 1
@@ -829,6 +830,7 @@ fl_bright:
   onFailure: buffer
   readoutPriority: baseline
   userParameter:
+    move_method: signal
     'off': 0.0
     'on': 3.0
     tol: 0.01
@@ -1385,8 +1387,8 @@ smargon:
   deviceConfig:
     prefix: http://x10sa-smargopolo.psi.ch:3000
   deviceTags:
-  - smargon
   - motors
+  - smargon
   enabled: true
   onFailure: buffer
   readoutPriority: baseline

pxii_bec/macros/APscripts.py

@@ -95,6 +95,12 @@ def angle(e, *hkl):
 
     return a
 
+##############
+def resol(dstart, dend, nsteps):
+    res = (Dstart - Dend) / nsteps
+    print(f"resolution is: {res}")
+    return res 
+
 
 #####################
 ### setenergy
@@ -215,11 +221,38 @@ def rock(**kwargs):
     peak = GaussianModel()
     background = LinearModel()
     model = peak + background
-    maxy = max(data_y)
-    indmax = np.argmax(data_y)
-    xm = data_x[indmax]
 
-    p = model.make_params(amplitude=maxy, center=xm)
+    # Estimate background from the edges of the data
+    n = len(data_x)
+    edge = max(5, n // 10)
+
+    xb = np.r_[data_x[:edge], data_x[-edge:]]
+    yb = np.r_[data_y[:edge], data_y[-edge:]]
+
+    slope, intercept = np.polyfit(xb, yb, 1)
+
+    # Background-subtracted signal
+    bg = slope * data_x + intercept
+    ys = data_y - bg
+
+    # Gaussian initial estimates
+    center = data_x[np.argmax(ys)]
+    amplitude = np.trapezoid(np.maximum(ys, 0), data_x)
+
+    p = model.make_params(
+        amplitude=amplitude,
+        center=center,
+        sigma=(data_x[-1] - data_x[0]) / 20,
+        slope=slope,
+        intercept=intercept
+    )
+
+    
+   # maxy = max(data_y)
+   # indmax = np.argmax(data_y)
+   # xm = data_x[indmax]
+   # p = model.make_params(amplitude=maxy, center=xm)
+   
     g = model.fit(data_y, p, x=data_x)
     print(f'Center: {g.params["center"].value:.5f} FWHM: {g.params["fwhm"].value:.5f}')
 

pxii_bec/macros/Undu_Helpers.py

@@ -214,7 +214,7 @@ def gap_harm(e=12.4):
 
 # test with estart = 6, end_en = 7.5 # should be 3 files , estimated time: 90 min 
 
-def long_gscan(estart=6, end_en=30.5, g_low=4.5, g_high=9.0, nsteps=1500):
+def long_gscan(estart=6, end_en=30.5, g_low=4.5, g_high=9.0, nsteps=1000, sett = 1): # resol = 4.5 um 
     import time
     import numpy as np
 
@@ -227,7 +227,7 @@ def long_gscan(estart=6, end_en=30.5, g_low=4.5, g_high=9.0, nsteps=1500):
     )
     fe_h_size0 = dev.fe_sl_xsize.user_readback.get()
     fe_v_size0 = dev.fe_sl_ysize.user_readback.get()
-    umv(dev.fe_sl_xsize, 0.5)
+    umv(dev.fe_sl_xsize, 0.3)   # 25 urad 
     umv(dev.fe_sl_ysize, 0.3)
     
     dock_area = bec.gui.new("LongGapScan", geometry = [4000, 900, 1200, 700])
@@ -252,7 +252,7 @@ def long_gscan(estart=6, end_en=30.5, g_low=4.5, g_high=9.0, nsteps=1500):
 
     while en < end_en:
         if en >= 17:
-            nsteps = 1000
+            nsteps = 778  # same resol as with 100 steps for 4.5 
             g_high = 7.5 
         #sete(en)
         #time.sleep(0.2)
@@ -265,7 +265,7 @@ def long_gscan(estart=6, end_en=30.5, g_low=4.5, g_high=9.0, nsteps=1500):
         
         print(f"setting energy to {en}")
         time.sleep(0.2)
-        ds = scans.line_scan(dev.id_gap, g_low, g_high, steps=nsteps, exp_time=0.1, relative=False)
+        ds = scans.line_scan(dev.id_gap, g_low, g_high, steps=nsteps, settling_time= sett, exp_time=0.1, relative=False)
         gap_data = ds.scan.live_data.id_gap.id_gap.val
         bpm_data = ds.scan.live_data.lu_bpmsum.lu_bpmsum.val
         wr.plot(x=gap_data, y=bpm_data)
@@ -282,7 +282,7 @@ def long_gscan(estart=6, end_en=30.5, g_low=4.5, g_high=9.0, nsteps=1500):
         
         print(f"setting energy to {en}")
         time.sleep(0.2)
-        ds = scans.line_scan(dev.id_gap, g_high, g_low, steps=nsteps, exp_time=0.1, relative=False)
+        ds = scans.line_scan(dev.id_gap, g_high, g_low, steps=nsteps, settling_time= sett, exp_time=0.1, relative=False)
         gap_data = ds.scan.live_data.id_gap.id_gap.val
         bpm_data = ds.scan.live_data.lu_bpmsum.lu_bpmsum.val
         wr.plot(x=gap_data, y=bpm_data)
@@ -307,7 +307,7 @@ def long_gscan(estart=6, end_en=30.5, g_low=4.5, g_high=9.0, nsteps=1500):
 
 
 ##################################
-def gscan(centre=0, gomax=0, detune=0):
+def gscan(centre=0, gomax=0, detune=0,sett=1):
     """
     Scan the ID GAP and go
     to the max of lu_bpm intensity
@@ -330,13 +330,21 @@ def gscan(centre=0, gomax=0, detune=0):
 
     g0 = dev.id_gap.readback.get()
     deltag = 0.05
-    ds = scans.line_scan(dev.id_gap, -deltag, deltag, steps=30, exp_time=0.5, relative=True)
+    steps = 30
+    resol = 2*deltag/steps
+    print(f"Gap resol is {resol} mm ")
+    
+    ds = scans.line_scan(dev.id_gap, -deltag, deltag, steps=steps,  settling_time=sett,exp_time=0.5, relative=True)
     gap_data = ds.scan.live_data.id_gap.id_gap.val
     bpm_data = ds.scan.live_data.lu_bpmsum.lu_bpmsum.val
+    
+    gap_data=np.array(gap_data)
+    bpm_data = np.array(bpm_data)
+    
 
-    # maxy = max(bpm_data)
-    # indmax = np.argmax(bpm_data)
-    # gm = gap_data[indmax]
+    maxy = max(bpm_data)
+    indmax = np.argmax(bpm_data)
+    gm = gap_data[indmax]
 
     gcen, xm = fit_plot(gap_data, bpm_data, model="voigt")
 

pxii_bec/macros/bcu_bpm_positions.py

@@ -0,0 +1,84 @@
+"""
+Output the beam position from the BCU XBPM.
+"""
+
+import time
+import numpy as np
+
+
+def compute_norm(n_readings=5):
+    """Average normalized BPM X/Y from repeated sets of four channels."""
+
+    xn_values = []
+    yn_values = []
+
+    for _ in range(n_readings):
+        bpm1 = dev.bcu_bpm1.read()["bcu_bpm1"]["value"]
+        bpm2 = dev.bcu_bpm2.read()["bcu_bpm2"]["value"]
+        bpm3 = dev.bcu_bpm3.read()["bcu_bpm3"]["value"]
+        bpm4 = dev.bcu_bpm4.read()["bcu_bpm4"]["value"]
+
+        total = bpm1 + bpm2 + bpm3 + bpm4
+
+        if total == 0:
+            continue
+
+        xn = ((bpm1 + bpm2) - (bpm3 + bpm4)) / total
+        yn = ((bpm1 + bpm3) - (bpm2 + bpm4)) / total
+
+        xn_values.append(xn)
+        yn_values.append(yn)
+
+        time.sleep(0.1)
+
+    return float(np.mean(xn_values)), float(np.mean(yn_values))
+
+
+def get_bcu_beampos():
+    """Determine beam position"""
+    bpm_xn, bpm_yn = compute_norm()
+    bpm_x = dev.bcu_bpm_x.position
+    bpm_y = dev.bcu_bpm_y.position
+
+    beam_x = beam_offset_um(bpm_xn, bpm_x, BPM_X_Calibration)
+
+    beam_y = beam_offset_um(bpm_yn, bpm_y, BPM_Y_Calibration)
+
+    return beam_x, beam_y
+
+
+def beam_offset_um(bpm_signal, bpm_motor_position, calibration):
+    """
+    Calculate beam position relative to the sample-camera crosshair.
+
+    Returns position in microns.
+    """
+
+    if calibration.crosshair_signal is None:
+        raise ValueError("Camera crosshair reference has not been calibrated.")
+
+    if calibration.crosshair_bpm_position is None:
+        raise ValueError("Crosshair BPM motor position is not defined.")
+
+    if abs(bpm_signal) > calibration.signal_limit:
+        raise ValueError(
+            f"BPM signal {bpm_signal:.3f} is outside the "
+            f"calibration range ±{calibration.signal_limit:.3f}"
+        )
+
+    # Has the BPM itself moved since the crosshair was calibrated?
+    motor_shift = bpm_motor_position - calibration.crosshair_bpm_position
+
+    # Convert current BPM signal to equivalent BPM position
+    current_local = np.polyval(calibration.coeffs, bpm_signal)
+
+    # Convert the BPM signal measured at the camera crosshair
+    # to equivalent BPM position
+    crosshair_local = np.polyval(calibration.coeffs, calibration.crosshair_signal)
+
+    signal_shift = current_local - crosshair_local
+
+    # BPM motion and apparent beam motion have opposite signs
+    beam_shift_mm = motor_shift - signal_shift
+
+    return float(beam_shift_mm * 1000)

pxii_bec/macros/beamline_planner.py

@@ -20,6 +20,7 @@ class BeamlineState(str, Enum):
     BEAMSTOP_ALIGNMENT = "beamstop_alignment"
     MAINTENANCE = "maintenance"
     XTAL_SNAPSHOT = "xtal_snapshot"
+    DET_DIODE = "det_diode"
 
 
 class TemperatureMode(str, Enum):

pxii_bec/macros/beampos_with_en.py~

@@ -1,9 +0,0 @@
-def beampos_with_en(planner):
-    deltat = 5 
-    record_beampos_until("23:00",planner,interval_minutes=deltat)
-    bl_energy(16000)
-    record_beampos_until("05:00",planner,interval_minutes=deltat)
-    bl_energy(12400)    
-    record_beampos_until("11:00",planner,interval_minutes=deltat)
-
-    

pxii_bec/macros/build_devices.py

@@ -1,8 +1,13 @@
 """Build the sample environment devices"""
 
 import os
+import sys
+from pathlib import Path
+
 import yaml
 
+from pxii_bec.scripts import beamline_context as bl
+
 # from position_device import PositionDevice
 
 
@@ -26,7 +31,7 @@ def motor_resolver(bec_name):
     raise ValueError(f"Cannot resolve motor for '{bec_name}'")
 
 
-def build_devices(yaml_file, mock_devices):
+def build_devices(yaml_file):
     """Build devices from the beamline states yaml"""
 
     state_devices = {}
@@ -35,7 +40,6 @@ def build_devices(yaml_file, mock_devices):
         data = yaml.safe_load(f)
 
     for bec_name, cfg in data.items():
-
         user = cfg.get("userParameter")
 
         # Skip devices without user parameters
@@ -43,11 +47,9 @@ def build_devices(yaml_file, mock_devices):
             continue
 
         if bec_name in ("aerotech", "smargon"):
-
             tol = user.get("tol", 0.01)
 
             for axis in ["x", "y", "z", "u", "chi", "phi"]:
-
                 positions = user.get(axis)
 
                 if positions is None:
@@ -61,17 +63,16 @@ def build_devices(yaml_file, mock_devices):
                     positions=positions,
                     tol=tol,
                     allow_arbitrary=True,
-                    use_mock=component_name in mock_devices,
                 )
 
                 state_devices[component_name] = pos_dev
 
         else:
-
             tol = user.get("tol", 0.1)
             allow_arbitrary = user.get("type") == "continuous"
+            move_method = user.get("move_method", "motor")
 
-            positions = {k: v for k, v in user.items() if k not in ("type", "tol")}
+            positions = {k: v for k, v in user.items() if k not in ("type", "tol", "move_method")}
 
             pos_dev = PositionDevice(
                 bec_name=bec_name,
@@ -79,7 +80,7 @@ def build_devices(yaml_file, mock_devices):
                 positions=positions,
                 tol=tol,
                 allow_arbitrary=allow_arbitrary,
-                use_mock=bec_name in mock_devices,
+                move_method=move_method,
             )
 
             state_devices[bec_name] = pos_dev
@@ -87,7 +88,6 @@ def build_devices(yaml_file, mock_devices):
     return state_devices
 
 
-
 def update_position(
     device,
     position,
@@ -109,19 +109,14 @@ def update_position(
     # Nested device (e.g. aerotech.x.in)
     # ------------------------------------------------------------------
     if axis is not None:
-
         if axis not in params:
-            raise KeyError(
-                f"Unknown axis '{axis}'. "
-                f"Valid axes: {list(params.keys())}"
-            )
+            raise KeyError(f"Unknown axis '{axis}'. Valid axes: {list(params.keys())}")
 
         axis_params = params[axis].copy()
 
         if position not in axis_params:
             raise KeyError(
-                f"Unknown position '{position}'. "
-                f"Valid positions: {list(axis_params.keys())}"
+                f"Unknown position '{position}'. Valid positions: {list(axis_params.keys())}"
             )
 
         old_value = axis_params[position]
@@ -139,36 +134,22 @@ def update_position(
 
         axis_params[position] = value
 
-        device.update_user_parameter(
-            {axis: axis_params}
-        )
+        device.update_user_parameter({axis: axis_params})
 
-        print(
-            f"Updated {device.name}.{axis}.{position}: "
-            f"{old_value} -> {value}"
-        )
+        print(f"Updated {device.name}.{axis}.{position}: {old_value} -> {value}")
 
         return
 
     # ------------------------------------------------------------------
     # Flat device (e.g. bs_z)
     # ------------------------------------------------------------------
-    nested_keys = [
-        k for k, v in params.items()
-        if isinstance(v, dict)
-    ]
+    nested_keys = [k for k, v in params.items() if isinstance(v, dict)]
 
     if nested_keys:
-        raise ValueError(
-            f"{device.name} has nested parameters "
-            f"({nested_keys}); specify axis=..."
-        )
+        raise ValueError(f"{device.name} has nested parameters ({nested_keys}); specify axis=...")
 
     if position not in params:
-        raise KeyError(
-            f"Unknown position '{position}'. "
-            f"Valid positions: {list(params.keys())}"
-        )
+        raise KeyError(f"Unknown position '{position}'. Valid positions: {list(params.keys())}")
 
     old_value = params[position]
 
@@ -183,31 +164,22 @@ def update_position(
             f"Use force=True if intentional."
         )
 
-    device.update_user_parameter(
-        {position: value}
-    )
+    device.update_user_parameter({position: value})
 
-    print(
-        f"Updated {device.name}.{position}: "
-        f"{old_value} -> {value}"
-    )
+    print(f"Updated {device.name}.{position}: {old_value} -> {value}")
 
 
 def save_and_reload():
     """
     Save current configuration to YAML and rebuild
-    the beamline environment.
+    the beamline environment using the same deployment.
     """
 
-    beamline = os.environ["BEAMLINE"]
+    print(f"(saving changes to {bl.env.devices_file})")
 
-    env = Environment(beamline=beamline)
+    bec.config.save_current_session(bl.env.devices_file)
 
-    print(f"(saving changes to {env.devices_file})")
-
-    bec.config.save_current_session(env.devices_file)
-
-    return init_beamline_environment()
+    bl.reload()
 
 
 def update_positions(updates):
@@ -230,7 +202,6 @@ def update_positions(updates):
     """
 
     for update in updates:
-
         update_position(
             device=update["device"],
             position=update["position"],
@@ -258,25 +229,20 @@ def save_current_position(
     params = device.user_parameter
 
     if axis is not None:
-
         if axis not in params:
-            raise KeyError(
-                f"Unknown axis '{axis}'. "
-                f"Valid axes: {list(params.keys())}"
-            )
+            raise KeyError(f"Unknown axis '{axis}'. Valid axes: {list(params.keys())}")
 
         axis_params = params[axis].copy()
 
         if position not in axis_params:
             raise KeyError(
-                f"Unknown position '{position}'. "
-                f"Valid positions: {list(axis_params.keys())}"
+                f"Unknown position '{position}'. Valid positions: {list(axis_params.keys())}"
             )
 
         old_value = axis_params[position]
         print(f"Old value is {old_value}")
 
-        current_value = get_current_value(device, axis=None)
+        current_value = get_current_value(device, axis)
 
         delta = abs(current_value - old_value)
 
@@ -291,35 +257,21 @@ def save_current_position(
 
         axis_params[position] = current_value
 
-        device.update_user_parameter(
-            {axis: axis_params}
-        )
+        device.update_user_parameter({axis: axis_params})
 
-        print(
-            f"Updated {device.name}.{axis}.{position}: "
-            f"{old_value} -> {current_value}"
-        )
+        print(f"Updated {device.name}.{axis}.{position}: {old_value} -> {current_value}")
 
         return
     # ------------------------------------------------------------------
     # Flat device (e.g. bs_z)
     # ------------------------------------------------------------------
-    nested_keys = [
-        k for k, v in params.items()
-        if isinstance(v, dict)
-    ]
+    nested_keys = [k for k, v in params.items() if isinstance(v, dict)]
 
     if nested_keys:
-        raise ValueError(
-            f"{device.name} has nested parameters "
-            f"({nested_keys}); specify axis=..."
-        )
+        raise ValueError(f"{device.name} has nested parameters ({nested_keys}); specify axis=...")
 
     if position not in params:
-        raise KeyError(
-            f"Unknown position '{position}'. "
-            f"Valid positions: {list(params.keys())}"
-        )
+        raise KeyError(f"Unknown position '{position}'. Valid positions: {list(params.keys())}")
 
     old_value = params[position]
     current_value = get_current_value(device, axis=None)
@@ -335,14 +287,9 @@ def save_current_position(
             f"Use force=True if intentional."
         )
 
-    device.update_user_parameter(
-        {position: current_value}
-    )
+    device.update_user_parameter({position: current_value})
 
-    print(
-        f"Updated {device.name}.{position}: "
-        f"{old_value} -> {current_value}"
-    )
+    print(f"Updated {device.name}.{position}: {old_value} -> {current_value}")
 
 
 def get_current_value(device, axis=None):
@@ -359,4 +306,4 @@ def get_current_value(device, axis=None):
 
     signal_name = next(iter(data))
 
-    return round(data[signal_name]["value"],3)
+    return round(data[signal_name]["value"], 8)

pxii_bec/macros/calibrate_bpm.py

@@ -1,142 +1,296 @@
 """
-BPM calibration using inverse PCHIP spline.
-
-Calibration concept:
-    BPM signal -> beam position relative to BPM
-
-Final beam position:
-    beam_absolute = bpm_motor_position + local_beam_position
-    beam_relative_to_camera_crosshair = beam_absolute - camera_crosshair_position
+BPM calibration.
 """
 
 import time
-import os
-from datetime import datetime, timedelta
-from pathlib import Path
+from datetime import datetime
 import numpy as np
 import pandas as pd
 import matplotlib.pyplot as plt
-from pyparsing import results
-from scipy.optimize import curve_fit
-import scipy.special as special
-from scipy.interpolate import PchipInterpolator, CubicSpline
 
 
-def read_currents():
-    """Read the individual bpm channels"""
-    vals_bpm1 = []
-    vals_bpm2 = []
-    vals_bpm3 = []
-    vals_bpm4 = []
+# def compute_norm(n_readings=5):
+#     """Average normalized BPM X/Y from repeated sets of four channels."""
 
-    for _ in range(5):
-        vals_bpm1.append(dev.bcu_bpm1.read()["bcu_bpm1"]["value"])
-        time.sleep(0.1)
-    for _ in range(5):
-        vals_bpm2.append(dev.bcu_bpm2.read()["bcu_bpm2"]["value"])
-        time.sleep(0.1)
-    for _ in range(5):
-        vals_bpm3.append(dev.bcu_bpm3.read()["bcu_bpm3"]["value"])
-        time.sleep(0.1)
-    for _ in range(5):
-        vals_bpm4.append(dev.bcu_bpm4.read()["bcu_bpm4"]["value"])
-        time.sleep(0.1)
-    bpm1 = np.mean(vals_bpm1)
-    bpm2 = np.mean(vals_bpm2)
-    bpm3 = np.mean(vals_bpm3)
-    bpm4 = np.mean(vals_bpm4)
-   
-    return {"bpm1": bpm1, "bpm2": bpm2, "bpm3": bpm3, "bpm4": bpm4}
+#     xn_values = []
+#     yn_values = []
+
+#     for _ in range(n_readings):
+#         bpm1 = dev.bcu_bpm1.read()["bcu_bpm1"]["value"]
+#         bpm2 = dev.bcu_bpm2.read()["bcu_bpm2"]["value"]
+#         bpm3 = dev.bcu_bpm3.read()["bcu_bpm3"]["value"]
+#         bpm4 = dev.bcu_bpm4.read()["bcu_bpm4"]["value"]
+
+#         total = bpm1 + bpm2 + bpm3 + bpm4
+
+#         if total == 0:
+#             continue
+
+#         xn = ((bpm1 + bpm2) - (bpm3 + bpm4)) / total
+#         yn = ((bpm1 + bpm3) - (bpm2 + bpm4)) / total
+
+#         xn_values.append(xn)
+#         yn_values.append(yn)
+
+#         time.sleep(0.1)
+
+#     return float(np.mean(xn_values)), float(np.mean(yn_values))
 
 
-def compute_norm():
-    """Compute normalised x and y positions"""
-    readings = read_currents()
-    total = readings["bpm1"] + readings["bpm2"] + readings["bpm3"] + readings["bpm4"]
-    # print(f"Sum is {total:.2f}")
-    xn = (
-        (readings["bpm1"] + readings["bpm2"]) - (readings["bpm3"] + readings["bpm4"])
-    ) / total
-    yn = (
-        (readings["bpm1"] + readings["bpm3"]) - (readings["bpm2"] + readings["bpm4"])
-    ) / total
-    # return xn, yn, readings["bpm1"], readings["bpm2"], readings["bpm3"], readings["bpm4"]
-    return float(xn), float(yn)
-
-
-# def compute_pos():
-#     """Compute the beam x and y positions"""
-#     cal = {
-#         "bcu_bpm": {
-#             "x_slope": 7.778,
-#             "x_off": -1.019,
-#             "y_slope": -8.625,
-#             "y_off": -1.935,
-#         }
-#     }
-#     xn, yn, bpm1, bpm2, bpm3, bpm4 = compute_norm()
-#     # print(f"Normalised positions are {xn}, {yn}")
-#     x_pos = (xn - cal["bcu_bpm"]["x_off"]) / cal["bcu_bpm"]["x_slope"]
-#     y_pos = (yn - cal["bcu_bpm"]["y_off"]) / cal["bcu_bpm"]["y_slope"]
-#     return x_pos, y_pos
-
-
-def run_calibration():
+def run_calibration(env):
     """Scan bpm in x and y and record normalised x and y positions"""
+    centred_x = 0.444
+    centred_y = -1.5404
+    now = datetime.now()
+    fnow = filenow()
     print(f"moving bpm to {centred_x}, {centred_y}")
     umv(dev.bcu_bpm_y, centred_y)
     umv(dev.bcu_bpm_x, centred_x)
-    
+
     # # Calibrate in X
-    umv(dev.ss_sl_xsize,0.02)
+    umv(dev.ss_sl_xsize, 0.02)
     umv(dev.ss_sl_xsize, 3.0)
-    filename = f"luts/bpm_calib/{fnow}_bcu_bpm_x.csv"
+    filename = env.lut_dir / "bpm_calib" / f"bcu_bpm_x_{fnow}.csv"
     xn_data = []
-    x_positions = np.linspace(centred_x - 0.6, centred_x + 0.6, 101)
+    x_positions = np.linspace(centred_x - 0.4, centred_x + 0.4, 101)
     with open(filename, "w", encoding="utf-8") as f:
         f.write(f"Xpos,Xnorm,BPM1,BPM2,BPM3,BPM4\n")
     for x in x_positions:
         umv(dev.bcu_bpm_x, x)
-        print(f"Moving BPM to {x} mm")
+        print(f"Moving BPM X to {x} mm")
         time.sleep(0.2)
-        xn, yn, bpm1, bpm2, bpm3, bpm4 = compute_norm()
+        xn, yn = compute_norm()
         xn_data.append(xn)
         with open(filename, "a", encoding="utf-8") as f:
-            f.write(f"{x},{xn},{bpm1},{bpm2},{bpm3},{bpm4}\n")
+            f.write(f"{x},{xn}\n")
 
-    # Calibrate in Y
+    # # Calibrate in Y
     umv(dev.bcu_bpm_x, centred_x)
     umv(dev.ss_sl_xsize, 3.0)
     umv(dev.ss_sl_ysize, 0.03)
-    filename = f"luts/bpm_calib/{fnow}_bcu_bpm_y.csv"
-    y_positions = np.linspace(centred_y - 0.1, centred_y + 0.1, 251)
+    filename = env.lut_dir / "bpm_calib" / f"bcu_bpm_y_{fnow}.csv"
+    y_positions = np.linspace(centred_y - 0.05, centred_y + 0.05, 101)
     yn_data = []
     with open(filename, "w", encoding="utf-8") as f:
         f.write(f"Ypos,Ynorm,BPM1,BPM2,BPM3,BPM4\n")
     for y in y_positions:
         umv(dev.bcu_bpm_y, y)
-        print(f"Moving BPM to {y} mm")
+        print(f"Moving BPM Y to {y} mm")
         time.sleep(0.2)
-        xn, yn, bpm1, bpm2, bpm3, bpm4 = compute_norm()
+        xn, yn = compute_norm()
         yn_data.append(yn)
         with open(filename, "a", encoding="utf-8") as f:
-            f.write(f"{y},{yn},{bpm1},{bpm2},{bpm3},{bpm4}\n")
+            f.write(f"{y},{yn}\n")
 
     umv(dev.bcu_bpm_x, centred_x)
     umv(dev.bcu_bpm_y, centred_y)
     umv(dev.ss_sl_ysize, 3.0)
 
+def read_bpm_calibration_data(
+    env,
+    x_y,
+    rundate,
+    bpm_name="bcu_bpm",
+):
+    """Read BPM calibration data from CSV."""
 
+    filename = env.lut_dir / f"bpm_calib/{bpm_name}_{x_y}_{rundate}.csv"
 
-def fit_bpm_data(x_y, rundate="02072130_", bpm_name="bcu_bpm"):
-    CALIB_DIR = "luts/bpm_calib/"
-    filename = os.path.join(CALIB_DIR, f"{rundate}{bpm_name}_{x_y}.csv")
-    data = np.loadtxt(filename, delimiter=",", skiprows=1)
+    data = np.loadtxt(
+        filename,
+        delimiter=",",
+        skiprows=1,
+    )
 
-    xdata = data[:, 0]   # scan / motor position
-    ydata = data[:, 1]   # normalised BPM signal
+    position = data[:, 0]
+    bpm_norm = data[:, 1]
 
-    return xdata, ydata
+    return position, bpm_norm
+
+# def fit_bpm_linear(env, filename, x_y, plot=True):
+#     """
+#     Fit the linear region of a BPM calibration scan.
+
+#     The fit is:
+
+#         normalised_signal = slope * bpm_motor_position + intercept
+#     """
+
+#     file = env.lut_dir / "bpm_calib" / filename
+
+#     bpm_centre = {
+#         "x": 0.441,
+#         "y": -1.5405,
+#     }
+
+#     linear_limit = {
+#         "x": 0.59,
+#         "y": 0.5,
+#     }
+#     data = pd.read_csv(file)
+
+#     if x_y == "x":
+#         pos = data["Xpos"].to_numpy()
+#         signal = data["Xnorm"].to_numpy()
+#     elif x_y == "y":
+#         pos = data["Ypos"].to_numpy()
+#         signal = data["Ynorm"].to_numpy()
+#     else:
+#         raise ValueError("x_y must be 'x' or 'y'")
+
+#     limit = linear_limit[x_y]
+
+#     mask = np.abs(signal) < limit
+
+#     pos_fit = pos[mask]
+#     signal_fit = signal[mask]
+
+#     slope, intercept = np.polyfit(pos_fit, signal_fit, 1)
+
+#     fitted = slope * pos_fit + intercept
+
+#     ss_res = np.sum((signal_fit - fitted) ** 2)
+#     ss_tot = np.sum((signal_fit - np.mean(signal_fit)) ** 2)
+
+#     r_squared = 1 - ss_res / ss_tot
+
+#     zero_position = -intercept / slope
+
+#     result = {
+#         "axis": x_y,
+#         "slope": float(slope),
+#         "intercept": float(intercept),
+#         "linear_limit": limit,
+#         "r_squared": float(r_squared),
+#         "zero_position": float(zero_position),
+#         "centred_bpm_position": bpm_centre[x_y],
+#     }
+#     print(result)
+#     print(f"class BPM_{x_y.upper()}_Calibration:")
+#     print(f"    slope = {slope:.6f}")
+#     print(f"    intercept = {intercept:.6f}")
+#     print(f"    linear_limit = {limit:.2f}")
+#     print(f"    centred_position = {bpm_centre[x_y]:.3f}")
+#     print(f"    crosshair_signal = None")
+#     print(f"    crosshair_bpm_position = None")
+
+#     if plot:
+#         xline = np.linspace(pos_fit.min(), pos_fit.max(), 200)
+
+#         plt.figure(figsize=(7, 5))
+#         plt.scatter(pos, signal, alpha=0.4, label="all data")
+#         plt.scatter(pos_fit, signal_fit, label="linear region")
+#         plt.plot(
+#             xline,
+#             slope * xline + intercept,
+#             label="linear fit",
+#         )
+
+#         plt.axvline(
+#             bpm_centre[x_y],
+#             linestyle="--",
+#             label="centred BPM position",
+#         )
+
+#         plt.xlabel(f"BPM {x_y} position (mm)")
+#         plt.ylabel(f"BPM {x_y} normalised signal")
+#         plt.title(f"BPM {x_y.upper()} calibration")
+#         plt.grid(True)
+#         plt.legend()
+#         plt.tight_layout()
+#         plt.show()
+
+#     return result   
 
 
 
+
+def polynomial_calibration(
+    env,
+    x_y,
+    rundate,
+    bpm_name="bcu_bpm",
+    degree=5,
+    signal_limit=0.9,
+    plot=True,
+):
+
+    # BPM centre during calibration run
+    bpm_centre = {"x": 0.441,"y": -1.5405,}
+
+
+    xdata, ydata = read_bpm_calibration_data(
+        env,
+        x_y,
+        rundate=rundate,
+        bpm_name=bpm_name,
+    )
+
+    # Keep useful BPM range
+    mask = np.abs(ydata) < signal_limit
+
+    position = xdata[mask]
+    signal = ydata[mask]
+
+    # Fit position as a function of BPM signal
+    coeffs = np.polyfit(signal, position, degree)
+
+    fitted_position = np.polyval(coeffs, signal)
+
+    residual_um = (position - fitted_position) * 1000
+
+    rms_error_um = np.sqrt(np.mean(residual_um**2))
+    max_error_um = np.max(np.abs(residual_um))
+
+    ss_res = np.sum((position - fitted_position) ** 2)
+    ss_tot = np.sum((position - np.mean(position)) ** 2)
+    r_squared = 1 - ss_res / ss_tot
+
+    print(f"{x_y.upper()} calibration")
+    print(f"Polynomial degree: {degree}")
+    print(f"Signal range: ±{signal_limit}")
+    print(f"R²: {r_squared:.6f}")
+    print(f"RMS error: {rms_error_um:.2f} µm")
+    print(f"Max error: {max_error_um:.2f} µm")
+    print(f"Coefficients: {coeffs}")
+
+    print(f"class BPM_{x_y.upper()}_Calibration:")
+    print(f"    coeffs = {coeffs}")
+    print(f"    signal_limit = {signal_limit:.2f}")
+    print(f"    centred_position = {bpm_centre[x_y]:.3f}")
+    print(f"    crosshair_signal = None")
+    print(f"    crosshair_bpm_position = None")
+
+    if plot:
+        signal_fit = np.linspace(
+            signal.min(),
+            signal.max(),
+            500,
+        )
+
+        position_fit = np.polyval(coeffs, signal_fit)
+
+        plt.figure(figsize=(7, 5))
+        plt.scatter(
+            ydata,
+            xdata,
+            alpha=0.4,
+            label="Calibration data",
+        )
+        plt.plot(
+            signal_fit,
+            position_fit,
+            label=f"{degree}th-order polynomial",
+        )
+
+        plt.xlabel(f"BPM {x_y} normalized signal")
+        plt.ylabel(f"BPM {x_y} position (mm)")
+        plt.title(f"BPM {x_y.upper()} calibration")
+        plt.grid(True)
+        plt.legend()
+        plt.tight_layout()
+        plt.show()
+
+    return coeffs
+
+

pxii_bec/macros/calibrate_bpm.py.bak

@@ -1,670 +0,0 @@
-"""
-BPM calibration using inverse PCHIP spline.
-
-Calibration concept:
-    BPM signal -> beam position relative to BPM
-
-Final beam position:
-    beam_absolute = bpm_motor_position + local_beam_position
-    beam_relative_to_camera_crosshair = beam_absolute - camera_crosshair_position
-"""
-
-import time
-import os
-from datetime import datetime, timedelta
-from pathlib import Path
-import numpy as np
-import pandas as pd
-import matplotlib.pyplot as plt
-from pyparsing import results
-from scipy.optimize import curve_fit
-import scipy.special as special
-from scipy.interpolate import PchipInterpolator, CubicSpline
-
-
-def read_currents():
-    """Read the individual bpm channels"""
-    vals_bpm1 = []
-    vals_bpm2 = []
-    vals_bpm3 = []
-    vals_bpm4 = []
-
-    for _ in range(5):
-        vals_bpm1.append(dev.bcu_bpm1.read()["bcu_bpm1"]["value"])
-        time.sleep(0.1)
-    for _ in range(5):
-        vals_bpm2.append(dev.bcu_bpm2.read()["bcu_bpm2"]["value"])
-        time.sleep(0.1)
-    for _ in range(5):
-        vals_bpm3.append(dev.bcu_bpm3.read()["bcu_bpm3"]["value"])
-        time.sleep(0.1)
-    for _ in range(5):
-        vals_bpm4.append(dev.bcu_bpm4.read()["bcu_bpm4"]["value"])
-        time.sleep(0.1)
-    bpm1 = np.mean(vals_bpm1)
-    bpm2 = np.mean(vals_bpm2)
-    bpm3 = np.mean(vals_bpm3)
-    bpm4 = np.mean(vals_bpm4)
-   
-    return {"bpm1": bpm1, "bpm2": bpm2, "bpm3": bpm3, "bpm4": bpm4}
-
-
-def compute_norm():
-    """Compute normalised x and y positions"""
-    readings = read_currents()
-    total = readings["bpm1"] + readings["bpm2"] + readings["bpm3"] + readings["bpm4"]
-    # print(f"Sum is {total:.2f}")
-    xn = (
-        (readings["bpm1"] + readings["bpm2"]) - (readings["bpm3"] + readings["bpm4"])
-    ) / total
-    yn = (
-        (readings["bpm1"] + readings["bpm3"]) - (readings["bpm2"] + readings["bpm4"])
-    ) / total
-    # return xn, yn, readings["bpm1"], readings["bpm2"], readings["bpm3"], readings["bpm4"]
-    return float(xn), float(yn)
-
-
-# def compute_pos():
-#     """Compute the beam x and y positions"""
-#     cal = {
-#         "bcu_bpm": {
-#             "x_slope": 7.778,
-#             "x_off": -1.019,
-#             "y_slope": -8.625,
-#             "y_off": -1.935,
-#         }
-#     }
-#     xn, yn, bpm1, bpm2, bpm3, bpm4 = compute_norm()
-#     # print(f"Normalised positions are {xn}, {yn}")
-#     x_pos = (xn - cal["bcu_bpm"]["x_off"]) / cal["bcu_bpm"]["x_slope"]
-#     y_pos = (yn - cal["bcu_bpm"]["y_off"]) / cal["bcu_bpm"]["y_slope"]
-#     return x_pos, y_pos
-
-
-def run_calibration():
-    """Scan bpm in x and y and record normalised x and y positions"""
-    # centred_x = 0.749
-    # centred_y = -0.660
-    now = datetime.now()
-    fnow = now.strftime("%d%m%H%M")
-    centred_x = 0.46
-    centred_y = -1.51
-    print(f"moving bpm to {centred_x}, {centred_y}")
-    umv(dev.bcu_bpm_y, centred_y)
-    umv(dev.bcu_bpm_x, centred_x)
-    
-    # # Calibrate in X
-    umv(dev.ss_sl_xsize,0.02)
-    umv(dev.ss_sl_xsize, 3.0)
-    filename = f"luts/bpm_calib/{fnow}_bcu_bpm_x.csv"
-    xn_data = []
-    x_positions = np.linspace(centred_x - 0.6, centred_x + 0.6, 101)
-    with open(filename, "w", encoding="utf-8") as f:
-        f.write(f"Xpos,Xnorm,BPM1,BPM2,BPM3,BPM4\n")
-    for x in x_positions:
-        umv(dev.bcu_bpm_x, x)
-        print(f"Moving BPM to {x} mm")
-        time.sleep(0.2)
-        xn, yn, bpm1, bpm2, bpm3, bpm4 = compute_norm()
-        xn_data.append(xn)
-        with open(filename, "a", encoding="utf-8") as f:
-            f.write(f"{x},{xn},{bpm1},{bpm2},{bpm3},{bpm4}\n")
-
-    # Calibrate in Y
-    umv(dev.bcu_bpm_x, centred_x)
-    umv(dev.ss_sl_xsize, 3.0)
-    umv(dev.ss_sl_ysize, 0.03)
-    filename = f"luts/bpm_calib/{fnow}_bcu_bpm_y.csv"
-    y_positions = np.linspace(centred_y - 0.1, centred_y + 0.1, 251)
-    yn_data = []
-    with open(filename, "w", encoding="utf-8") as f:
-        f.write(f"Ypos,Ynorm,BPM1,BPM2,BPM3,BPM4\n")
-    for y in y_positions:
-        umv(dev.bcu_bpm_y, y)
-        print(f"Moving BPM to {y} mm")
-        time.sleep(0.2)
-        xn, yn, bpm1, bpm2, bpm3, bpm4 = compute_norm()
-        yn_data.append(yn)
-        with open(filename, "a", encoding="utf-8") as f:
-            f.write(f"{y},{yn},{bpm1},{bpm2},{bpm3},{bpm4}\n")
-
-    umv(dev.bcu_bpm_x, centred_x)
-    umv(dev.bcu_bpm_y, centred_y)
-    umv(dev.ss_sl_ysize, 3.0)
-
-
-
-def fit_bpm_data(x_y, rundate="02072130_", bpm_name="bcu_bpm"):
-    CALIB_DIR = "luts/bpm_calib/"
-    filename = os.path.join(CALIB_DIR, f"{rundate}{bpm_name}_{x_y}.csv")
-    data = np.loadtxt(filename, delimiter=",", skiprows=1)
-
-    xdata = data[:, 0]   # scan / motor position
-    ydata = data[:, 1]   # normalised BPM signal
-
-    return xdata, ydata
-
-def linear_calibration(
-        x_y,
-        bpm_motor_position,
-        crosshair_bpm_signal,
-        rundate = "05071923_",
-        bpm_name = "bcu_bpm",
-        sanity_limit = 0.2,
-        plot = True
-):
-    calib_dir = "luts/bpm_calib"
-    xdata, ydata = fit_bpm_data(x_y, rundate = rundate, bpm_name = bpm_name)
-    # Remove non-linear region
-    if x_y == 'x':
-        mask_limit = 0.59
-    else:
-        mask_limit = 0.9
-    mask = np.abs(ydata) < mask_limit
-    xdata = xdata[mask]
-    ydata = ydata[mask]
-
-    m, c = np.polyfit(xdata, ydata, 1)
-
-    # --- Plot ---
-    plt.figure(figsize=(7, 5))
-    plt.scatter(xdata, ydata, color="red", label=f"BPM {x_y} data")
-
-    x_fit = np.linspace(min(xdata), max(xdata), 100)
-    plt.plot(x_fit, m * x_fit + c, "r--", label=f"Fit: Linear")
-
-    plt.xlabel("label")
-    plt.ylabel("label")
-    plt.title("title")
-    plt.legend()
-    plt.grid(True)
-    plt.tight_layout()
-    plt.show()
-
-def make_inverse_spline_calibration(
-    x_y,
-    bpm_motor_position,
-    crosshair_bpm_signal,
-    rundate="05071923_",
-    bpm_name="bcu_bpm",
-    sanity_limit=0.2,
-    plot=True,
-):
-    """
-    Create and save inverse BPM calibration.
-
-    Parameters
-    ----------
-    x_y : str
-        "x" or "y"
-    bpm_motor_position : float
-        Current physical BPM motor position during calibration.
-    crosshair_bpm_signal : float
-        Normalised BPM reading when the beam is on the sample-camera crosshair.
-    mask_limit : float
-        Ignore saturated BPM readings outside +/- mask_limit.
-    sanity_limit : float
-        Warn later if BPM motor has moved by more than this.
-    """
-    CALIB_DIR = "luts/bpm_calib/"
-    xdata, ydata = fit_bpm_data(x_y, rundate=rundate, bpm_name=bpm_name)
-
-    # Remove saturated BPM response
-    if x_y == 'x':
-        mask_limit = 0.98
-    else:
-        mask_limit = 0.91
-    mask = np.abs(ydata) < mask_limit
-    xdata = xdata[mask]
-    ydata = ydata[mask]
-
-    # For inverse calibration: BPM signal -> beam position
-    order = np.argsort(ydata)
-    bpm_signal = ydata[order]
-    local_position = xdata[order]
-
-    # Remove duplicate BPM signal values
-    bpm_signal, unique_idx = np.unique(bpm_signal, return_index=True)
-    local_position = local_position[unique_idx]
-
-    spline = PchipInterpolator(bpm_signal, local_position)
-
-    # Local beam position when beam is visually on camera crosshair
-    crosshair_local_position = float(spline(crosshair_bpm_signal))
-
-    # Absolute beamline coordinate of camera crosshair
-    crosshair_absolute_position = bpm_motor_position + crosshair_local_position
-
-    save_file = os.path.join(CALIB_DIR, f"{bpm_name}_{x_y}_inverse_spline.npz")
-
-    np.savez(
-        save_file,
-        bpm_signal=bpm_signal,
-        local_position=local_position,
-        bpm_motor_position=bpm_motor_position,
-        crosshair_bpm_signal=crosshair_bpm_signal,
-        crosshair_local_position=crosshair_local_position,
-        crosshair_absolute_position=crosshair_absolute_position,
-        mask_limit=mask_limit,
-        sanity_limit=sanity_limit,
-    )
-
-    print(f"Saved {x_y} BPM calibration to:")
-    print(save_file)
-    print()
-    print(f"Calibration BPM motor position: {bpm_motor_position:.6f}")
-    print(f"Crosshair BPM signal:          {crosshair_bpm_signal:.6f}")
-    print(f"Crosshair local position:      {crosshair_local_position:.6f}")
-    print(f"Crosshair absolute position:   {crosshair_absolute_position:.6f}")
-
-    if plot:
-        signal_fit = np.linspace(np.min(bpm_signal), np.max(bpm_signal), 1000)
-        pos_fit = spline(signal_fit)
-
-        plt.figure()
-        plt.scatter(ydata, xdata, label="calibration data")
-        plt.plot(signal_fit, pos_fit, label="inverse PCHIP spline", color='red')
-        plt.axvline(crosshair_bpm_signal, linestyle="--", label="camera crosshair")
-        plt.xlabel(f"{bpm_name}_{x_y} normalised signal")
-        plt.ylabel(f"{x_y} local beam position")
-        plt.title(f"{bpm_name}_{x_y} inverse calibration")
-        plt.grid(True)
-        plt.legend()
-        plt.show()
-
-    return save_file
-
-class BPMCalibration:
-    def __init__(self, filename):
-        data = np.load(filename)
-
-        self.bpm_signal = data["bpm_signal"]
-        self.local_position = data["local_position"]
-
-        self.bpm_motor_position_calib = float(data["bpm_motor_position"])
-        self.crosshair_absolute_position = float(data["crosshair_absolute_position"])
-        self.sanity_limit = float(data["sanity_limit"])
-
-        self.spline = PchipInterpolator(
-            self.bpm_signal,
-            self.local_position,
-        )
-
-    def beam_position(self, bpm_signal, bpm_motor_position, warn=True):
-        """
-        Return beam position relative to the sample-camera crosshair.
-        """
-
-        motor_shift = bpm_motor_position - self.bpm_motor_position_calib
-
-        if warn and abs(motor_shift) > self.sanity_limit:
-            print(
-                "Warning: BPM motor has moved since calibration: "
-                f"delta = {motor_shift:.3f} mm"
-            )
-
-        local_beam_position = float(self.spline(bpm_signal))
-        absolute_beam_position = bpm_motor_position + local_beam_position
-
-        relative_to_crosshair = (
-            absolute_beam_position - self.crosshair_absolute_position
-        )
-
-        return relative_to_crosshair
-
-def get_beampos():
-    x_cal = BPMCalibration("luts/bpm_calib/bcu_bpm_x_inverse_spline.npz")
-    y_cal = BPMCalibration("luts/bpm_calib/bcu_bpm_y_inverse_spline.npz")
-
-    bpm_x_signal, bpm_y_signal = compute_norm()
-
-    beam_x = x_cal.beam_position(
-        bpm_signal=bpm_x_signal,
-        bpm_motor_position=dev.bcu_bpm_x.position,
-    )
-
-    beam_y = y_cal.beam_position(
-        bpm_signal=bpm_y_signal,
-        bpm_motor_position=dev.bcu_bpm_y.position,
-    )
-
-    beam_x_mic = 1000 * beam_x
-    beam_y_mic = 1000 * beam_y
-
-    print(f"Beam relative to camera crosshair:")
-    print(f"X = {beam_x_mic:.4f} um")
-    print(f"Y = {beam_y_mic:.4f} um")
-
-    return beam_x_mic, beam_y_mic
-
-#!/usr/bin/env python
-"""Calculate photon flux from a Si diode current."""
-
-import math
-from dataclasses import dataclass
-
-import numpy as np
-
-
-@dataclass
-class PXIIDiodeConfig:
-    """Diode/beamline-specific constants."""
-
-    si_thickness_um: float = 12.0
-    al_thickness_um: float = 20.0
-    diode_offset_mm: float = 15.0
-
-    si_density: float = 2.33        # g/cm^3
-    al_density: float = 2.699       # g/cm^3
-    air_density: float = 1.205e-3   # g/cm^3
-
-    eps_si: float = 3.62            # eV per electron-hole pair
-
-
-def flux(
-    current_a: float,
-    energy_kev: float | None = None,
-    det_z: float | None = None,
-    config = PXIIDiodeConfig,
-) -> float:
-    """
-    Calculate photon flux from Si diode current.
-
-    Parameters
-    ----------
-    current_a:
-        Diode current in A.
-    energy_kev:
-        Beam energy in keV. Will use the current energy if 
-        not specified.
-    det_z:
-        Detector distance (dev.det_z) (mm). Will use the current
-        detector distance if not specified.
-    config:
-        DiodeConfig containing Si/Al thicknesses and material constants.
-
-    Returns
-    -------
-    float
-        Flux in photons/s.
-    """
-    if energy_kev is None:
-        energy_kev = get_current_energy() / 1000
-
-    if det_z is None:
-        det_z = dev.det_z.position
-
-    air_path_mm = det_z + config.diode_offset_mm
-
-    si_thickness_cm = config.si_thickness_um / 10_000.0
-    al_thickness_cm = config.al_thickness_um / 10_000.0
-
-    log_energy = math.log10(energy_kev)
-
-    # Photoelectric mass attenuation coefficient fits
-    si_poly = np.poly1d([0.0789, -0.477, -2.238, 4.158])
-    al_poly = np.poly1d([0.0638, -0.413, -2.349, 4.106])
-    air_poly = np.poly1d([0.928, -2.348, -1.026, 3.153])
-
-    a_si = 10 ** si_poly(log_energy)
-    a_al = 10 ** al_poly(log_energy)
-    a_air = 10 ** air_poly(log_energy)
-
-    # Fraction of beam absorbed in active Si layer
-    si_absorbed_fraction = -np.expm1(
-        -a_si * config.si_density * si_thickness_cm
-    )
-
-    if si_absorbed_fraction <= 0:
-        raise ValueError("Calculated Si absorption fraction is zero or negative.")
-
-    # Flux before Al/air correction
-    flux = (
-        1000.0
-        * current_a
-        * config.eps_si
-        / 1.602
-        / energy_kev
-        / si_absorbed_fraction
-        * 1.0e13
-    )
-
-    # Transmission through Al and air
-    al_transmission = np.exp(
-        -a_al * config.al_density * al_thickness_cm
-    )
-
-    air_transmission = np.exp(
-        -a_air * config.air_density * air_path_mm / 10.0
-    )
-
-    flux /= al_transmission
-    flux /= air_transmission
-
-
-    print(f"Flux = {flux: .3g} ph/s")
-
-    return float(flux)
-
-def calibrate_bpm_flux(ave = 5):
-    # d['det_z'].move('mse')
-    # d['det_y'].move('xi')
-    # d['det_z'].move(200)
-
-    dir = "luts/bpm_calib/"
-    filename = dir + "bpm_flux_20260708b.csv"
-    with open(filename, 'w', encoding="utf-8") as f:
-        f.write("Energy,Transmission,BPMSum,Flux\n")
-    print("Energy,Transmission,BPMSum,Flux\n")
-    energies = [12400,7000,8000,9000,10000,11000,13000,14000,15000,16000,18000,20000,25000]
-    transmissions = [1, 0.9, 0.8, 0.7, 0.6, 0.5,0.4,0.3,0.2,0.1]
-
-    for energy in energies:
-        print(f"Moving to energy {energy}")
-        bl_energy(energy)
-
-        for set_trans in transmissions:
-            print(f"Setting transmission to {set_trans}")
-            dev.transm.put(set_trans)
-            time.sleep(5)
-
-            flux_values = []
-            bpm_values = []
-
-            print(f"Measuring flux and bpmsum")
-            for _ in range(ave):
-                diode_current = dev.xidiode.read()['xidiode']['value']
-                fl = flux(diode_current)
-                bpmsum = dev.bcu_bpmsum.read()['bcu_bpmsum']['value']
-                flux_values.append(fl)
-                bpm_values. append(bpmsum)
-
-            actual_trans = dev.transm.read()['transm']['value']
-            flux_ave = np.mean(flux_values)
-            bpmsum_ave = np.mean(bpm_values)
-            print(f"{energy:.5g},{actual_trans:.3g},{bpmsum_ave:.4g},{flux_ave:.4f}\n")
-
-            with open(filename, 'a') as f:
-                f.write(f"{energy:.5g},{actual_trans:.3g},{bpmsum_ave:.4g},{flux_ave:.4f}\n")
-    bl_energy(12400)
-    umv(dev.det_y, 70)
-    umv(dev.id_gap, 20)
-
-def analyse_bpm_flux_calibration(
-    filename="luts/bpm_calib/bpm_flux_20260708b.csv",
-    force_zero_intercept=False,
-):
-    data = pd.read_csv(filename)
-
-    results = []
-    with open("luts/bpm_calib/bpm_flux_fits.csv", "w") as f:
-        f.write("Energy,Slope,Intercept\n")
-    for energy, group in data.groupby("Energy"):
-        x = group["BPMSum"].to_numpy(dtype=float)
-        print(f"x is {x}")
-        y = group["Flux"].to_numpy(dtype=float)
-
-        if force_zero_intercept:
-            # y = m*x
-            slope = np.sum(x * y) / np.sum(x * x)
-            intercept = 0.0
-        else:
-            # y = m*x + c
-            slope, intercept = np.polyfit(x, y, deg=1)
-
-        y_fit = slope * x + intercept
-
-        # --- Plot ---
-        plt.figure(figsize=(7, 5))
-        plt.scatter(x, y, color="red")
-
-        x_fit = np.linspace(min(x), max(x), 100)
-        plt.plot(x_fit, slope * x_fit + intercept, 
-        "r--", label=f"Fit: Slope = {slope:.5g}, Intercept = {intercept:.5g}")
-
-        plt.xlabel("BPMSUm")
-        plt.ylabel("Flux")
-        plt.title(f"{energy} eV")
-        plt.legend()
-        plt.grid(True)
-        plt.tight_layout()
-        plt.savefig(f"luts/bpm_calib/{energy}_bpm_flux")
-        # plt.show()
-
-        residuals = y - y_fit
-        ss_res = np.sum(residuals**2)
-        ss_tot = np.sum((y - np.mean(y))**2)
-
-        r2 = 1 - ss_res / ss_tot if ss_tot > 0 else np.nan
-
-        results.append(
-            {
-                "Energy": energy,
-                "Slope": f"{slope:.5g}",
-                "Intercept": f"{intercept:.5g}",
-                "R2": f"{r2:.3g}",
-                # "N": len(group),
-            }
-        )
-
-    results = pd.DataFrame(results)
-
-    results.to_csv("luts/bpm_calib/bpm_flux_fits.csv", index=False)
-    dict = results.to_dict()
-    print(dict)
-    return results
-
-def linear_bpm_flux_calib():
-    filename = 'luts/bpm_calib/fits-25.csv'
-    df = pd.read_csv(filename)
-    energies = df["Energy"].to_numpy()
-    intercepts = df["Intercept"].to_numpy()
-    slopes = df["Slope"].to_numpy()
-    # --- Plot ---
-    plt.figure(figsize=(7, 5))
-    plt.scatter(energies, slopes, color="red")
-
-    x_fit = np.linspace(min(energies), max(energies), 100)
-    coeffs = np.polyfit(energies, slopes,1)
-    y_fit = np.poly1d(coeffs)
-    plt.plot(x_fit, y_fit(x_fit)),    
-    plt.xlabel("Energy eV")
-    plt.ylabel("Slope")
-    plt.title(f"Energies vs slopes, coeffs = {coeffs}")
-    # plt.legend()
-    plt.grid(True)
-    plt.tight_layout()
-    plt.savefig(f"luts/bpm_calib/slopes_vs_energy")
-    plt.show()
-    return coeffs
-
-
-# def get_flux_calibration(energy):
-#     # filename = "luts/bpm_calib/bpm_flux_fits.csv"
-#     filename = 'luts/bpm_calib/fits-25.csv'
-#     default_slope = 4.42e8
-#     default_intercept = -1.66e11
-
-#     try:
-#         df = pd.read_csv(filename)
-#         energies = df["Energy"].to_numpy()
-#         intercepts = df["Intercept"].to_numpy()
-#         slopes = df["Slope"].to_numpy()
-#         slope_interp = CubicSpline(energies, slopes)
-#         intercept_interp = CubicSpline(energies, intercepts)
-#         slope = slope_interp(energy)
-#         intercept = intercept_interp(energy)
-#         print(f"Using data from calibration file")
-#          # --- Plot ---
-#         plt.figure(figsize=(7, 5))
-#         plt.scatter(energies, slopes, color="red")
-
-#         x_fit = np.linspace(min(energies), max(energies), 100)
-#         coeffs = np.polyfit(energies, slopes,1)
-#         print("linear co-efficients: ", coeffs)
-#         y_fit = np.poly1d(coeffs)
-#         plt.plot(x_fit, y_fit(x_fit)), 
-#         # "r--", label=f"Fit: Slope = {slope:.5g}, Intercept = {intercept:.5g}")
-      
-#         plt.xlabel("Energy")
-#         plt.ylabel("Slope")
-#         plt.title(f"{energy} eV")
-#         # plt.legend()
-#         plt.grid(True)
-#         plt.tight_layout()
-#         plt.show()
-#     except FileNotFoundError:
-#         print("Cannot access calibration file, using default parameters")
-#         slope = default_slope
-#         intercept = default_intercept
-    
-#     return coeffs
-
-
-def bpm2flux(energy=None, bpmsum=None, transm=None):
-    """
-    Calculate the flux from BCU BPM reading.
-    Option to provide energy, bpmsum reading and transmission. 
-    If these are not given, the current values are used.
-
-    Returns the flux at the sample (after attenuators) and a calculated
-    value of the full flux if attenuation was 0.
-    
-    """
-    calibration_min_energy = 6500
-    calibration_max_energy = 21000
-
-    if energy is None:
-        energy = get_current_energy()
-
-    if bpmsum is None:
-        bpmsum = dev.bcu_bpmsum.read()["bcu_bpmsum"]["value"]
-
-    if transm is None:
-        transm = dev.transm.read()['transm']['value']
-
-    if not calibration_min_energy <= energy <= calibration_max_energy:
-        raise ValueError(
-            f"Energy {energy:.0f} eV is outside the calibrated range "
-            f"{calibration_min_energy:.0f}–"
-            f"{calibration_max_energy:.0f} eV"
-        )
-
-    coeffs = [
-        1.06436277e5,
-        -6.38943746e8,
-    ]
-
-    calibration_factor = np.polyval(coeffs, energy)
-    flux = calibration_factor * bpmsum
-    full_flux = flux/transm
-    if transm == 1:
-        print(f"Transmission = {transm*100:.5g} %, Flux = {flux: .3g} ph/s,\n")
-    else:
-        print(f"Transmission = {transm*100: .5g} %, Flux = {flux: .3g} ph/s,\n"
-            f"Flux with 100 % transmission would be {full_flux: .3g} ph/s")
-
-    return float(flux), float(full_flux)
-
-
-    

pxii_bec/macros/dependencies.py

@@ -9,8 +9,8 @@ def planner_deps():
         ("aerotech_x", "out"): [("diag_y", "out"), ("bs_z", "safe")],
         ("diag_y", "scint"): [("aerotech_x", "out"), ("bs_z", "safe"), ("cryo_pos", "out")],
         ("diag_y", "i1"): [("aerotech_x", "out"), ("bs_z", "safe"), ("cryo_pos", "out")],
-        ("bs_pos", "out"): [("bs_z", "safe")],
-        ("bs_pos", "in"): [("bs_z", "safe")],
+        ("bs_pos", "out"): [("bs_z", "safe"),("bl_pos", "in")],
+        ("bs_pos", "in"): [("bs_z", "safe"),("bl_pos", "in")],
         ("diag_y", "out"): [("bs_z", "safe")],
         ("diag_y", "park"): [("bs_z", "safe")],
         ("coll_y", "out"): [("bs_z", "safe")],

pxii_bec/macros/devices_manager.py

@@ -10,14 +10,11 @@ class DeviceManager:
     """Class for building devices and attaching safety guards and policies."""
 
     @staticmethod
-    # def initialize_devices(state_devices_file, rest_devices_file, mock_devices):
-    def initialize_devices(devices_file, mock_devices):
+    def initialize_devices(devices_file):
         """
         Initializes sample environment devices from the specified file.
         """
-        devices = build_devices(devices_file, mock_devices)
-        # rest_devices = build_devices(rest_devices_file, mock_devices)
-        # devices.update(rest_devices)
+        devices = build_devices(devices_file)
         attach_guards(devices)
         attach_policies(devices)
         return devices

pxii_bec/macros/flux_macros.py

@@ -0,0 +1,298 @@
+"""
+Macros for measurimng flux from Si diode. Also for estimating the flux from
+the BCU BPM reading
+"""
+
+import numpy as np
+from dataclasses import dataclass
+
+
+@dataclass
+class PXIIDiodeConfig:
+    """Diode/beamline-specific constants."""
+
+    si_thickness_um: float = 12.0
+    al_thickness_um: float = 20.0
+    diode_offset_mm: float = 15.0
+
+    si_density: float = 2.33  # g/cm^3
+    al_density: float = 2.699  # g/cm^3
+    air_density: float = 1.205e-3  # g/cm^3
+
+    eps_si: float = 3.62  # eV per electron-hole pair
+
+
+def flux(
+    current_a: float,
+    energy_kev: float | None = None,
+    det_z: float | None = None,
+    config=PXIIDiodeConfig,
+) -> float:
+    """
+    Calculate photon flux from Si diode current.
+
+    Parameters
+    ----------
+    current_a:
+        Diode current in A.
+    energy_kev:
+        Beam energy in keV. Will use the current energy if
+        not specified.
+    det_z:
+        Detector distance (dev.det_z) (mm). Will use the current
+        detector distance if not specified.
+    config:
+        DiodeConfig containing Si/Al thicknesses and material constants.
+
+    Returns
+    -------
+    float
+        Flux in photons/s.
+    """
+    if energy_kev is None:
+        energy_kev = get_current_energy() / 1000
+
+    if det_z is None:
+        det_z = dev.det_z.position
+
+    air_path_mm = det_z + config.diode_offset_mm
+
+    si_thickness_cm = config.si_thickness_um / 10_000.0
+    al_thickness_cm = config.al_thickness_um / 10_000.0
+
+    log_energy = math.log10(energy_kev)
+
+    # Photoelectric mass attenuation coefficient fits
+    si_poly = np.poly1d([0.0789, -0.477, -2.238, 4.158])
+    al_poly = np.poly1d([0.0638, -0.413, -2.349, 4.106])
+    air_poly = np.poly1d([0.928, -2.348, -1.026, 3.153])
+
+    a_si = 10 ** si_poly(log_energy)
+    a_al = 10 ** al_poly(log_energy)
+    a_air = 10 ** air_poly(log_energy)
+
+    # Fraction of beam absorbed in active Si layer
+    si_absorbed_fraction = -np.expm1(-a_si * config.si_density * si_thickness_cm)
+
+    if si_absorbed_fraction <= 0:
+        raise ValueError("Calculated Si absorption fraction is zero or negative.")
+
+    # Flux before Al/air correction
+    flux = 1000.0 * current_a * config.eps_si / 1.602 / energy_kev / si_absorbed_fraction * 1.0e13
+
+    # Transmission through Al and air
+    al_transmission = np.exp(-a_al * config.al_density * al_thickness_cm)
+
+    air_transmission = np.exp(-a_air * config.air_density * air_path_mm / 10.0)
+
+    flux /= al_transmission
+    flux /= air_transmission
+
+    print(f"Flux = {flux: .3g} ph/s")
+
+    return float(flux)
+
+
+def calibrate_bpm_flux(ave=5):
+    # d['det_z'].move('mse')
+    # d['det_y'].move('xi')
+    # d['det_z'].move(200)
+
+    dir = "luts/bpm_calib/"
+    filename = dir + "bpm_flux_20260708b.csv"
+    with open(filename, "w", encoding="utf-8") as f:
+        f.write("Energy,Transmission,BPMSum,Flux\n")
+    print("Energy,Transmission,BPMSum,Flux\n")
+    energies = [
+        12400,
+        7000,
+        8000,
+        9000,
+        10000,
+        11000,
+        13000,
+        14000,
+        15000,
+        16000,
+        18000,
+        20000,
+        25000,
+    ]
+    transmissions = [1, 0.9, 0.8, 0.7, 0.6, 0.5, 0.4, 0.3, 0.2, 0.1]
+
+    for energy in energies:
+        print(f"Moving to energy {energy}")
+        bl_energy(energy)
+
+        for set_trans in transmissions:
+            print(f"Setting transmission to {set_trans}")
+            dev.transm.put(set_trans)
+            time.sleep(5)
+
+            flux_values = []
+            bpm_values = []
+
+            print(f"Measuring flux and bpmsum")
+            for _ in range(ave):
+                diode_current = dev.xidiode.read()["xidiode"]["value"]
+                fl = flux(diode_current)
+                bpmsum = dev.bcu_bpmsum.read()["bcu_bpmsum"]["value"]
+                flux_values.append(fl)
+                bpm_values.append(bpmsum)
+
+            actual_trans = dev.transm.read()["transm"]["value"]
+            flux_ave = np.mean(flux_values)
+            bpmsum_ave = np.mean(bpm_values)
+            print(f"{energy:.5g},{actual_trans:.3g},{bpmsum_ave:.4g},{flux_ave:.4f}\n")
+
+            with open(filename, "a") as f:
+                f.write(f"{energy:.5g},{actual_trans:.3g},{bpmsum_ave:.4g},{flux_ave:.4f}\n")
+    bl_energy(12400)
+    umv(dev.det_y, 70)
+    umv(dev.id_gap, 20)
+
+
+def analyse_bpm_flux_calibration(
+    filename="luts/bpm_calib/bpm_flux_20260708b.csv", force_zero_intercept=False
+):
+    data = pd.read_csv(filename)
+
+    results = []
+    with open("luts/bpm_calib/bpm_flux_fits.csv", "w") as f:
+        f.write("Energy,Slope,Intercept\n")
+    for energy, group in data.groupby("Energy"):
+        x = group["BPMSum"].to_numpy(dtype=float)
+        print(f"x is {x}")
+        y = group["Flux"].to_numpy(dtype=float)
+
+        if force_zero_intercept:
+            # y = m*x
+            slope = np.sum(x * y) / np.sum(x * x)
+            intercept = 0.0
+        else:
+            # y = m*x + c
+            slope, intercept = np.polyfit(x, y, deg=1)
+
+        y_fit = slope * x + intercept
+
+        # --- Plot ---
+        plt.figure(figsize=(7, 5))
+        plt.scatter(x, y, color="red")
+
+        x_fit = np.linspace(min(x), max(x), 100)
+        plt.plot(
+            x_fit,
+            slope * x_fit + intercept,
+            "r--",
+            label=f"Fit: Slope = {slope:.5g}, Intercept = {intercept:.5g}",
+        )
+
+        plt.xlabel("BPMSUm")
+        plt.ylabel("Flux")
+        plt.title(f"{energy} eV")
+        plt.legend()
+        plt.grid(True)
+        plt.tight_layout()
+        plt.savefig(f"luts/bpm_calib/{energy}_bpm_flux")
+        # plt.show()
+
+        residuals = y - y_fit
+        ss_res = np.sum(residuals**2)
+        ss_tot = np.sum((y - np.mean(y)) ** 2)
+
+        r2 = 1 - ss_res / ss_tot if ss_tot > 0 else np.nan
+
+        results.append(
+            {
+                "Energy": energy,
+                "Slope": f"{slope:.5g}",
+                "Intercept": f"{intercept:.5g}",
+                "R2": f"{r2:.3g}",
+                # "N": len(group),
+            }
+        )
+
+    results = pd.DataFrame(results)
+
+    results.to_csv("luts/bpm_calib/bpm_flux_fits.csv", index=False)
+    dict = results.to_dict()
+    print(dict)
+    return results
+
+
+def linear_bpm_flux_calib():
+    filename = "luts/bpm_calib/fits-25.csv"
+    df = pd.read_csv(filename)
+    energies = df["Energy"].to_numpy()
+    intercepts = df["Intercept"].to_numpy()
+    slopes = df["Slope"].to_numpy()
+    # --- Plot ---
+    plt.figure(figsize=(7, 5))
+    plt.scatter(energies, slopes, color="red")
+
+    x_fit = np.linspace(min(energies), max(energies), 100)
+    coeffs = np.polyfit(energies, slopes, 1)
+    y_fit = np.poly1d(coeffs)
+    plt.plot(x_fit, y_fit(x_fit)),
+    plt.xlabel("Energy eV")
+    plt.ylabel("Slope")
+    plt.title(f"Energies vs slopes, coeffs = {coeffs}")
+    # plt.legend()
+    plt.grid(True)
+    plt.tight_layout()
+    plt.savefig(f"luts/bpm_calib/slopes_vs_energy")
+    plt.show()
+    return coeffs
+
+
+def bpm2flux(energy=None, bpmsum=None, transm=None):
+    """
+    Calculate the flux from BCU BPM reading.
+    Option to provide energy, bpmsum reading and transmission.
+    If these are not given, the current values are used.
+
+    Returns the flux at the sample (after attenuators) and a calculated
+    value of the full flux if attenuation was 0.
+
+    """
+    calibration_min_energy = 6500
+    calibration_max_energy = 21000
+
+    if energy is None:
+        energy = get_current_energy()
+
+    if bpmsum is None:
+        bpmsum = dev.bcu_bpmsum.read(cached=True)["bcu_bpmsum"]["value"]
+        # bpmsum =  read_av(dev.bcu_bpmsum,10) # do not add to traffic jam
+
+    if transm is None:
+        transm = dev.transm.read()["transm"]["value"]
+
+    if not calibration_min_energy <= energy <= calibration_max_energy:
+        raise ValueError(
+            f"Energy {energy:.0f} eV is outside the calibrated range "
+            f"{calibration_min_energy:.0f}–"
+            f"{calibration_max_energy:.0f} eV"
+        )
+
+    # coeffs = [
+    #    1.06436277e5,
+    #    -6.38943746e8,
+    # ]
+    coeff = [-0.00026213, 0.01171065, -0.05379571, 0.04478851]
+    if energy > 24000:
+        coeff = [4.83704287e-05, -3.17174523e-03, 7.33131161e-02, -5.91926617e-01, 1.66572494e00]
+
+    # calibration_factor = np.polyval(coeffs, energy)
+    calibration_factor = np.polyval(coeff, energy / 1000) * 1e9
+    flux = calibration_factor * bpmsum
+    full_flux = flux / transm
+    if transm == 1:
+        print(f"Transmission = {transm*100:.5g} %, Flux = {flux: .3g} ph/s,\n")
+    else:
+        print(
+            f"Transmission = {transm*100: .5g} %, Flux = {flux: .3g} ph/s,\n"
+            f"Flux with 100 % transmission would be {full_flux: .3g} ph/s"
+        )
+
+    return float(flux), float(full_flux)

pxii_bec/macros/init_beamline.py

@@ -1,144 +0,0 @@
-"""Initialise sample environment devices and beamline states"""
-
-import os
-from dataclasses import dataclass
-
-# from devices_manager import DeviceManager
-# from beamline_state_manager import DefineStatesManager
-# from dependencies import planner_deps
-# from beamline_planner import StateChangePlanner
-
-
-@dataclass
-class Environment:
-    """Define the device and state yaml files for each beamline"""
-
-    device_mocks: dict[str, bool] = None
-    beamline: str = None
-
-    beamline_configs = {
-        # "x06sa" : {
-        #     "devdir": "/sls/x10sa/config/bec/production/pxi_bec/pxi_bec/device_configs/",
-        #     "state": "pxi-state-devices.yaml",
-        #     "rest": "pxi-rest-devices.yaml",
-        #     "states": "pxi-beamline-states.yaml",
-        # },
-        "x10sa": {
-            "devdir": "/sls/x10sa/config/bec/production/pxii_bec/pxii_bec/device_configs/",
-            "devices": "pxii-devices.yaml",
-            "states": "pxii-beamline-states.yaml",
-        },
-        "x06da": {
-            "devdir": "/sls/x06da/config/bec/production/pxiii_bec/pxiii_bec/device_configs/",
-            "state": "pxiii-devices.yaml",
-            "states": "pxiii-beamline-states.yaml",
-        },
-    }
-
-    def __post_init__(self):
-
-        if self.device_mocks is None:
-            self.device_mocks = {
-                "aerotech_x": "real",
-                "aerotech_y": "real",
-                "aerotech_z": "real",
-                "aerotech_u": "real",
-                "bl_bright": "real",
-                "bl_pos": "real",
-                "bs_pos": "real",
-                "bs_z": "real",
-                "coll_y": "real",
-                "cryo_pos": "real",
-                "det_cov": "real",
-                "diag_y": "real",
-                "fl_bright": "real",
-                "smargon_x": "real",
-                "smargon_y": "real",
-                "smargon_z": "real",
-                "smargon_chi": "real",
-                "smargon_phi": "real",
-                "xrf_pos": "mock",
-            }
-
-        mocks = sorted(name for name, backend in self.device_mocks.items() if backend == "mock")
-
-        reals = sorted(name for name, backend in self.device_mocks.items() if backend == "real")
-
-        # print(f"Mock devices ({len(mocks)}): {mocks}")
-        # print(f"Real devices ({len(reals)}): {reals}")
-
-        try:
-            cfg = self.beamline_configs[self.beamline]
-        except KeyError:
-            raise ValueError(f"Unknown beamline '{self.beamline}'")
-
-        devdir = cfg["devdir"]
-
-        self.devices_file = devdir + cfg["devices"]
-        self.states_file = devdir + cfg["states"]
-
-    mock_start_positions = {
-        "aerotech_x": "in",
-        "aerotech_y": "mount",
-        "aerotech_z": "mount",
-        "aerotech_u": "mount",
-        "bl_bright": "off",
-        "bl_pos": "out",
-        "bs_pos": "in",
-        "bs_z": "safe",
-        "coll_y": "park",
-        "cryo_pos": "in",
-        "det_cov": "close",
-        "diag_y": "park",
-        "fl_bright": "off",
-        "smargon_x": "mount",
-        "smargon_y": "mount",
-        "smargon_z": "mount",
-        "smargon_chi": "mount",
-        "smargon_phi": "mount",
-        "xrf_pos": "out",
-    }
-
-    @property
-    def mock_devices(self):
-        mock_names = set()
-        for name, device in self.device_mocks.items():
-            if device == "mock":
-                mock_names.add(name)
-        return mock_names
-
-
-def init_beamline_environment():
-    """
-    Initializes the beamline with real or mock devices.
-    """
-
-    beamline = os.environ["BEAMLINE"]
-
-    env = Environment(beamline=beamline)
-
-    # Initialize devices
-    device_manager = DeviceManager()
-    devices = device_manager.initialize_devices(env.devices_file, env.mock_devices)
-
-    # Initialize states
-    state_manager = DefineStatesManager()
-    states, allow_modifiers = state_manager.initialize_states(env.states_file)
-
-    # Setup dependencies
-    deps = planner_deps()
-
-    # Setup planner
-    planner = StateChangePlanner(devices, states, allow_modifiers, deps)
-    print("Initializing beamline state planner and devices")
-
-    # Set mock motor starting positions
-
-    for mock_dev, pos in env.mock_start_positions.items():
-        if mock_dev in env.mock_devices:
-            try:
-                devices[mock_dev].set_position(pos)
-            except Exception as e:
-                print(f"{e} error when setting {mock_dev}")
-
-    return devices, planner

pxii_bec/macros/katscripts.py

@@ -39,11 +39,11 @@ def move_xeye(direction):
 
 
 def close_fe_slits(close=True):
-    """Close FE slits to (0.3,0.3) or open the FE slits to (2.0,1.0)"""
+    """Close FE slits to (0.3,0.3) or open the FE slits to (0.9,0.9)"""
     if close:
         umv(dev.fe_sl_xsize, 0.3, dev.fe_sl_ysize, 0.3)
     else:
-        umv(dev.fe_sl_xsize, 1.5, dev.fe_sl_ysize, 1.0)
+        umv(dev.fe_sl_xsize, 0.9, dev.fe_sl_ysize, 0.9)
 
 
 def energy_checks(start_energy=8000, end_energy=20000):
@@ -106,27 +106,31 @@ def energy_checks(start_energy=8000, end_energy=20000):
 def scan_fe():
     """Scan the FE slits in X and Y"""
     print("Scanning FE slits in X")
-    umv(dev.fe_sxsize, 0.05, dev.fe_sysize, 2)
-    go_to_peak(dev.fe_sxcen, dev.lu_bpmsum, -1.5, 1.5, 20, confirm=False, relative=True)
+    umv(dev.fe_sl_xsize, 0.05, dev.fe_sl_ysize, 1)
+    go_to_peak(dev.fe_sl_xcen, dev.lu_bpmsum, -1.5, 1.5, 20, confirm=False, relative=True)
     print("Scanning FE slits in Y")
-    umv(dev.fe_sxsize, 0.3, dev.fe_sysize, 0.05)
-    go_to_peak(dev.fe_sycen, dev.lu_bpmsum, -1.2, 1.2, 15, confirm=False, relative=True)
+    umv(dev.fe_sl_xsize, 1, dev.fe_sl_ysize, 0.05)
+    go_to_peak(dev.fe_sl_ycen, dev.lu_bpmsum, -1.2, 1.2, 15, confirm=False, relative=True)
     close_fe_slits(False)
 
 
-def scan_gap():
+def scan_gap(gomax = False, close = True, fine = False):
     """Scan the gap and monitor on lu_bpmsum or bcu_bpmsum"""
-    # close_fe_slits(True)
+    close_fe_slits(close)
     current_gap = dev.id_gap.read()["id_gap"]["value"]
     if current_gap > 4.6:
         go_to_peak(
-            dev.id_gap, dev.lu_bpmsum, -0.1, 0.1, 25, relative=True, confirm=False, gomax=False
+            dev.id_gap, dev.lu_bpmsum, -0.15, 0.15, 25, relative=True, gap=True, confirm=False, gomax=gomax
         )
     else:
         go_to_peak(
-            dev.id_gap, dev.lu_bpmsum, 4.5, 4.75, 25, relative=False, confirm=False, gomax=False
+            dev.id_gap, dev.lu_bpmsum, 4.5, 4.75, 25, relative=False, gap=True, confirm=False, gomax=gomax
         )
-    # close_fe_slits(False)
+    if fine:
+        go_to_peak(
+            dev.id_gap, dev.lu_bpmsum, -0.012, 0.012, 6, relative=True, gap=True, confirm=False, gomax=gomax
+        )
+    close_fe_slits(False)
 
 
 def record_pitch():
@@ -168,20 +172,18 @@ def beam_centre_from_ss():
     return beampos
 
 
-def generate_dcm_lut(start_energy=6000, end_energy=30000, step=300):
+def generate_dcm_lut(env, planner, start_energy=6000, end_energy=30000, step=300, scint=True):
     """Generate a lookup table for the dcm from start_energy to end_energy in steps of step"""
-    now = datetime.now()
-    fnow = now.strftime("%d%m%H%M")
-    filename = f"luts/{fnow}_lut.csv"
+    filename = env.lut_dir / f"dcm_calibration_{filenow()}.csv"
     zoom = dev.scam_zoom.read()['scam_zoom']['value']
     with open(filename, "w", encoding="utf-8") as f:
-        f.write(f"energy,harmonic,gap,pitch,roll,perp,samcam_x,samcam_y,zoom={zoom}\n")
+        f.write(f"energy,harmonic,gap,pitch,samcam_x,samcam_y,zoom={zoom}\n")
 
     for energy in range(start_energy, end_energy, step):
         print(f"Moving to {energy: .0f}eV")
         bl_energy(energy, mono_scan=False)
         try:
-            scan_gap()
+            scan_gap(close=True, gomax=True, fine=True)
         except Exception as exc:
             print(f"Gap scan failed at {energy:.0f} eV: {exc}")
             set_gap(energy)  # set fallback gap here
@@ -191,18 +193,18 @@ def generate_dcm_lut(start_energy=6000, end_energy=30000, step=300):
         harm = Gap()
         h = harm.get_harmonic_by_energy(energy)
         pitch = dev.dcm_pitch.read()["dcm_pitch"]["value"]
-        roll = dev.dcm_froll.read()["dcm_froll"]["value"]
-        perp = dev.dcm_perp.read()["dcm_perp"]["value"]
-        # auto_exposure(cam="ss_xicam", max_iter=25)
-        # x = dev.ss_xicam_x.read()["ss_xicam_x"]["value"]
-        # y = dev.ss_xicam_y.read()["ss_xicam_y"]["value"]
-        auto_exposure(cam="samcam", max_iter=20)
-        x = dev.samcam_x.read()['samcam_x']['value']
-        y = dev.samcam_y.read()['samcam_y']['value']
-        with open(filename, "a", encoding="utf-8") as f:
-            f.write(
-                f"{energy:.0f},{h},{gap:.5g},{pitch:.5g},{roll:.4g},{perp:.5g},{x:.4g},{y:.4g}\n"
-            )
+        if scint:
+            x, y = read_samcam_scint(planner)
+            with open(filename, "a", encoding="utf-8") as f:
+                f.write(
+                    f"{energy:.0f},{h},{gap:.5g},{pitch:.5g},{x:.4g},{y:.4g}\n"
+                )
+        else:
+            with open(filename, "a", encoding="utf-8") as f:
+                f.write(
+                    f"{energy:.0f},{h},{gap:.5g},{pitch:.5g}\n"
+                )
+
     bl_energy(12400)
     umv(dev.id_gap, 20.0)
 
@@ -751,10 +753,10 @@ def auto_exposure(
         x = dev.ss_xicam_y.read()["ss_xicam_y"]["value"]
         if cam == "samcam":
             dev.samcam_exp.put(new_exp)
-            print(f" → Adjusting exposure to {new_exp:.4f} s")
+            # print(f" → Adjusting exposure to {new_exp:.4f} s")
         else:
             dev.ss_xicam_exp.put(new_exp)
-            print(f" → Adjusting exposure to {new_exp:.4f} s")
+            # print(f" → Adjusting exposure to {new_exp:.4f} s")
 
         exposure = new_exp
 
@@ -802,18 +804,18 @@ def filters(pos="out"):
         if SE.scin.checkpos():
             auto_exposure(cam="samcam", max_iter=20)
 
-def fnow():
-    now = datetime.now()
-    fnow = now.strftime("%Y%m%d-%H%M")
-    filedir = "luts/"
-    return filedir, fnow
+def filenow(date_only = False):
+    """Return current time formatted for filenames."""
+    if date_only:
+        return datetime.now().strftime("%Y%m%d")
+    else:
+        return datetime.now().strftime("%Y%m%d-%H%M")  
+
 
 def energy_vs_height(start_energy=11000, end_energy=25000, step=1000):
     """Correlate beam height with energy from start_energy to end_energy in steps of step"""
-    now = datetime.now()
-    # fnow = now.strftime("%d%m%H%M_")
-    fnow = now.strftime("%Y-%m-%d-%H%M")
-    filename = f"luts/{fnow}_energy-height-ss.csv"
+    filedir, fnow = fnow(dir = "test")
+    filename = f"{filedir}/{fnow}_energy-height-ss.csv"
     with open(filename, "w", encoding="utf-8") as f:
         # f.write("energy,bcu_norm_y,bcu_y,samcam_ypix,samcam_y\n")
         f.write("energy,y_px,diffy_mic\n")
@@ -887,13 +889,37 @@ def calc_from_centre(cam="samcam"):
 
 def get_cam(low_zoom, high_zoom):
     umv(dev.scam_zoom, low_zoom)
-    x = dev.samcam_x.read()["samcam_x"]["value"]
-    y = dev.samcam_y.read()["samcam_y"]["value"]
-    print(f"Low zoom: x = {x}, y = {y}")
+    auto_exposure(cam="samcam")
+    time.sleep(1)
+    low_x = dev.samcam_x.read()["samcam_x"]["value"]
+    low_y = dev.samcam_y.read()["samcam_y"]["value"]
+    print(f"Low zoom: x = {low_x}, y = {low_y}")
     umv(dev.scam_zoom, high_zoom)
-    x = dev.samcam_x.read()["samcam_x"]["value"]
-    y = dev.samcam_y.read()["samcam_y"]["value"]
-    print(f"High zoom: x = {x}, y = {y}")
+    auto_exposure(cam="samcam")
+    time.sleep(1)
+    high_x = dev.samcam_x.read()["samcam_x"]["value"]
+    high_y = dev.samcam_y.read()["samcam_y"]["value"]
+    xdiff = high_x - low_x
+    ydiff = high_y - low_y
+    print(f"High zoom: x = {high_x}, y = {high_y}")
+    print(f"X Diff = {xdiff:.3f}, Y Diff = {ydiff:.3f}")
+    if -15 <= xdiff <= 15:
+        print("X is OK")
+    elif xdiff > 15:
+        print("Adjust X tilt in clockwise direction")
+    else:
+        print("Adjust X tilt in anticlockwise direction")
+    # if xdiff > 0:
+    #     print("Adjust X tilt in clockwise direction")
+    # else:
+    #     print("Adjust X tilt in anticlockwise direction")
+    if -15 <= ydiff <= 15:
+        print("Y is OK")   
+    elif ydiff > 15:
+            print("Adjust Y tilt in anticlockwise direction")
+    else:
+        print("Adjust Y tilt in clockwise direction")
+
 
 
 # def read_currents():
@@ -977,41 +1003,61 @@ def mirror_pitch(mirror = 'vfm'):
     bl_energy(12400)
     results = np.array(results)
 
+def read_samcam_scint(planner,zoom=800):
+    # set transmission tp 10%
+    dev.transm.put(0.1)
+    # set camera zoom
+    umv(dev.scam_zoom, zoom)
+    # move to beam visualisation, ensure shutter is closed first
+    dev.bcu_shutter.put(0)
+    planner.move_to(BeamlineState.BEAM_VISUALISATION)
+    # open shutter
+    dev.bcu_shutter.put(1)
+    # set autoexposure
+    auto_exposure(cam="samcam", max_iter=25)
+    time.sleep(0.2)
+    x = dev.samcam_x.read()['samcam_x']['value']
+    y = dev.samcam_y.read()['samcam_y']['value']
+    dev.bcu_shutter.put(0)
+    return x, y
+
 def find_best_roll(low = 6000, high = 20000):
-    low_scanpoints = np.linspace(2.0, 7.0, 11)
+    # low_scanpoints = np.linspace(2.0, 7.0, 11)
+    low_scanpoints = np.linspace(2.5, 6.5, 9)
     high_scanpoints = np.linspace(1.0, 9.0,9)
     highs_y, lows_y = [], []
+    filedir, _ = filenow(dir="test")
+
+    # low energy
+    bl_energy(low)
+    # filename = f"{filedir}/roll_{low:.0f}.csv"
+    filename = env.lut_dir / f"roll_{low:.0f}.csv"
+    with open(filename, "w", encoding="utf-8") as f:
+        f.write("Roll, Y pos\n") 
+    for i in low_scanpoints:
+        umv(dev.dcm_froll, i)
+        time.sleep(0.2)
+        _, low_y = read_samcam_scint()
+        print(f"Ypos at sample cam is {low_y}")
+        lows_y.append(low_y)
+        with open(filename, "a", encoding="utf-8") as f:
+            f.write(f"{i:.4g},{low_y:.4g} \n") 
 
     # high energy
     bl_energy(high)
-    filename = f"./luts/roll_{high:.0f}.csv"
+    # filename = f"{filedir}/roll_{high:.0f}.csv"
+    filename = env_lut_dir / f"roll_{high:.0f}.csv"
     with open(filename, "w", encoding="utf-8") as f:
         f.write("Roll, Y pos\n") 
     auto_exposure("samcam", max_iter=25)
     for i in high_scanpoints:
         umv(dev.dcm_froll, i)
         time.sleep(0.2)
-        high_y = dev.samcam_y.read()['samcam_y']['value']
+        _, high_y = read_samcam_scint()
         print(f"Ypos at sample cam is {high_y}")
         highs_y.append(high_y)
         with open(filename, "a", encoding="utf-8") as f:
             f.write(f"{i:.4g},{high_y:.4g} \n") 
-
-    # low energy
-    bl_energy(low)
-    filename = f"./luts/roll_{low:.0f}.csv"
-    with open(filename, "w", encoding="utf-8") as f:
-        f.write("Roll, Y pos\n") 
-    auto_exposure("samcam", max_iter=25)
-    for i in low_scanpoints:
-        umv(dev.dcm_froll, i)
-        time.sleep(0.2)
-        low_y = dev.samcam_y.read()['samcam_y']['value']
-        print(f"Ypos at sample cam is {low_y}")
-        lows_y.append(low_y)
-        with open(filename, "a", encoding="utf-8") as f:
-            f.write(f"{i:.4g},{low_y:.4g} \n") 
-
     
     y_high = np.array(highs_y)
     y_low = np.array(lows_y)
@@ -1352,21 +1398,16 @@ def get_mirror_data(history_index: int):
     # }
 
 
-def record_beampos_until(stop_time_str, planner, interval_minutes=10):
+def record_beampos_until(stop_time_str, planner, env, interval_minutes=10):
     """
     Runs `record_beam()` every `interval_minutes` minutes
     until the next occurrence of the given stop time (HH:MM).
     """
 
     # Prepare the file
-    now = datetime.now()
-    fnow = now.strftime("%d%m%H%M")
-    fdir = "/sls/x10sa/config/bec/production/pxii_bec/pxii_bec/macros/"
-    fname = f"luts/beampos/beampos_{fnow}.csv"
-    filename = fdir + fname
+    filename = env.lut_dir / "beampos" / f"beampos_{filenow()}.csv"
     scam_zoom = dev.scam_zoom.read()['scam_zoom']['value']
-    umv(dev.scam_zoom, 1000)
-    dev.transm.put(0.1)
+    umv(dev.scam_zoom, 800)
     beam_size = "25 x 25"
 
     energy = get_current_energy()
@@ -1377,7 +1418,7 @@ def record_beampos_until(stop_time_str, planner, interval_minutes=10):
         f.write(f"Beam size: {beam_size}, Energy: {energy:.0f} eV, Zoom: {scam_zoom}\n")
         # f.write("Time,BPM_x,BPM_y,Scin_x,Scin_y,BPM_x,BPM_y\n")
         f.write("Time,Scin_x,Scin_y,BPM_x,BPM_y\n")
-
+    now = datetime.now()
 
     # Parse the stop time (HH:MM)
     stop_hour, stop_minute = map(int, stop_time_str.split(":"))
@@ -1415,8 +1456,8 @@ def record_beam(filename):
     """
     # Set transmission and bpm position
     dev.transm.put(0.1)
-    umv(dev.bcu_bpm_x, 0.43)
-    umv(dev.bcu_bpm_y, -1.543 )
+    umv(dev.bcu_bpm_x, BPM_X_Calibration.centred_position)
+    umv(dev.bcu_bpm_y, BPM_Y_Calibration.centred_position)
 
     # Open shutter and record positions
     dev.bcu_shutter.put(1)
@@ -1429,12 +1470,12 @@ def record_beam(filename):
         fnow = now.strftime("%H:%M:%S")
         xpos_cam = dev.samcam_x.read()["samcam_x"]["value"]
         ypos_cam = dev.samcam_y.read()["samcam_y"]["value"]
-        # x_bpm, y_bpm = get_beampos()
-        print(f"Time is {fnow}, X position is {xpos_cam:.5g}, Y position is {ypos_cam:.5g}")
-        x_bpm, y_bpm = compute_norm()
+        xpos_bpm, ypos_bpm = get_bcu_beampos()
+        print(f"Time is {fnow}, ScinX position is {xpos_cam:.5g}, ScinY position is {ypos_cam:.5g}")
+        print(f"BPM X is: {xpos_bpm}, BPM Y is {ypos_bpm}")
         f.write(
             # f"{fnow}, {xpos: .5g}, {ypos: .4g}, {xpos_cam: .5g}, {ypos_cam: .4g}\n"
-            f"{fnow},{xpos_cam: .5g},{ypos_cam: .4g},{x_bpm: .5g},{y_bpm: .5g}\n"
+            f"{fnow},{xpos_cam: .5g},{ypos_cam: .4g},{xpos_bpm: .5g},{ypos_bpm: .5g}\n"
         )
     # Close shutter
     dev.bcu_shutter.put(0)
@@ -1443,7 +1484,7 @@ def record_beam(filename):
 def analyse_beampos(file):
     """Plot the recorded beam position data"""
     # Load the CSV file into a Pandas DataFrame
-    filedir = "/sls/x10sa/config/bec/production/pxii_bec/pxii_bec/macros/luts/beampos/"
+    filedir = "/sls/x10sa/config/bec/test/pxii_bec/pxii_bec/macros/luts/beampos/"
     filename = filedir + file
     data = pd.read_csv(filename, skiprows=1)
     print(data.head())
@@ -1534,3 +1575,224 @@ def scan_scin_focus(planner, type = 'fine'):
 
     go_to_peak(dev.diag_z, dev.scam, min, max, steps, relative=True, gomax=True, confirm=False)
     dev.bcu_shutter.put(0)
+
+def daily_setup(energy=12400):
+
+    wf, results, log, progress = select_daily_setup_gui()
+
+    messages = []
+
+    progress_state = {
+    "Set energy": "todo",
+    "Measure flux": "todo",
+    "Measure beam position": "todo",
+    "Align collimator": "todo",
+    "Check beamstop": "todo",
+    "Align rotation axis": "todo",
+    }
+
+    update_log(log, messages, "Starting daily beamline setup")
+
+
+    # Step 1: set beamline 
+    set_progress(progress, progress_state, "Set energy", "running")
+    current_energy = get_current_energy()
+    if abs(energy - current_energy) > 1:
+        bl_energy(energy, window_name="DailySetup")
+    else:
+        gap_actual = dev.id_gap.position
+        gap_req = new_calc_gap(energy)
+        if abs(gap_actual-gap_req) > 0.05:
+            umv(dev.id_gap, gap_req)
+        mono_pitch_scan(window_name="DailySetup")
+        
+    set_progress(progress, progress_state, "Set energy", "done")
+    update_log(log, messages, f"Beamline energy is {energy} eV")
+
+    # Step 2: measure flux
+    set_progress(progress, progress_state, "Measure flux", "running")
+    flux1,flux2 = bpm2flux()
+    update_log(log, messages, f"Flux at 100% transmission is {flux2:.3g} ph/s")
+    set_progress(progress, progress_state, "Measure flux", "done")
+
+    # Step 3: measure beam position
+
+
+def update_log(log_box, messages, message):
+    messages.append(message)
+    log_box.set_plain_text("\n".join(messages))
+
+def update_progress(progress_box, progress_state):
+    status = {
+        "todo": "\u25CB",
+        "running": "\u25B6",
+        "done": "\u2713",
+        "warning": "\u26A0",
+        "failed": "\u2717",
+    }
+
+
+    text = "\n".join(
+        f"{status[state]} {label}"
+        for label, state in progress_state.items()
+    )
+    progress_box.set_plain_text(text)
+
+def set_progress(progress_box, progress_state, step, status):
+    """Change the status of one step and update the display."""
+    progress_state[step] = status
+    update_progress(progress_box, progress_state)
+
+
+def select_daily_setup_gui():
+    """Return widgets belonging to the Daily Setup GUI."""
+
+    window_name = "DailySetup"
+    open_windows = bec.gui.windows
+
+    if open_windows.get(window_name) is None:
+        gui = bec.gui.new(window_name)
+
+        plot = gui.new(
+            widget = "Waveform",
+            object_name = "Plot"
+        )
+
+        results = gui.new(
+            widget = "TextBox",
+            object_name = "Results",
+            where = "bottom",
+        )
+
+        log = gui.new(
+            widget = "TextBox",
+            object_name = "Info",
+            where = "bottom",
+        )
+
+        progress = gui.new(
+            widget = "TextBox",
+            object_name = "Progress",
+            where = "left",  
+        )
+
+    else:
+        gui = open_windows[window_name]
+        plot = gui.Plot
+        log = gui.Info
+        results = gui.Results
+        progress = gui.Progress
+
+    return plot, results, log, progress
+
+def bpm_errors(ch1, ch2, ch3, ch4):
+    total = ch1 + ch2 + ch3 + ch4
+
+    if total <= 0:
+        raise ValueError("BPM total intensity is zero or negative")
+
+    x_error = ((ch1 - ch3) + (ch2 - ch4)) / total
+    y_error = ((ch1 - ch2) + (ch3 - ch4)) / total
+
+    return x_error, y_error
+
+
+def find_bpm_centre(positions, errors):
+    """Fit the BPM difference signal and return its zero crossing."""
+    
+    slope, intercept = np.polyfit(positions, errors, 1)
+
+    if abs(slope) < 1e-6:
+        raise ValueError("BPM signal is too flat to determine centre")
+
+    centre = -intercept / slope
+
+    return centre, slope, intercept
+
+def centre_bcu_bpm_x():
+    positions = []
+    errors = []
+    startx = dev.bcu_bpm_x.position
+    starty = dev.bcu_bpm_y.position
+    scanpoints_x = np.linspace(startx-0.15, startx+0.15, 20)
+    for x in scanpoints_x:
+        umv(dev.bcu_bpm_x, x)
+        result = read_currents()
+        ch1 = float(result["bpm1"])
+        ch2 = float(result["bpm2"])
+        ch3 = float(result["bpm3"])
+        ch4 = float(result["bpm4"])
+        x_error, _ = bpm_errors(ch1, ch2, ch3, ch4)
+
+        print(f"BPM_X = {float(x)}, Xerror = {x_error}")
+
+        positions.append(x)
+        errors.append(x_error)
+
+    centre, slope, intercept = find_bpm_centre(positions,errors)
+    print(f"BPM centre in X is at {centre} mm")
+    umv(dev.bcu_bpm_x, centre)
+    return centre
+
+def centre_bcu_bpm_y():
+    positions = []
+    errors = []
+    starty = dev.bcu_bpm_y.position
+    scanpoints_y = np.linspace(starty-0.15, starty+0.15, 20)
+    for y in scanpoints_y:
+        umv(dev.bcu_bpm_y, y)
+        result = read_currents()
+        ch1 = float(result["bpm1"])
+        ch2 = float(result["bpm2"])
+        ch3 = float(result["bpm3"])
+        ch4 = float(result["bpm4"])
+        _, y_error = bpm_errors(ch1, ch2, ch3, ch4)
+
+        print(f"BPM_Y = {float(y)}, Yerror = {y_error}")
+
+        positions.append(y)
+        errors.append(y_error)
+
+    centre, slope, intercept = find_bpm_centre(positions,errors)
+    print(f"BPM centre in Y is at {centre} mm")
+    umv(dev.bcu_bpm_y, centre)
+    return centre
+
+def gap_check():
+    gaps = []
+    for repeat in range(6000, 30000, 2000):
+        gap = float(new_calc_gap(repeat))
+        gaps.append(gap)
+    print(gaps)
+
+def read_i1():
+    return dev.i1.read()['i1']['value']
+
+def cen_coll(d, planner):
+    """ Centre collimator by scanning against i1 diode """
+    dev.bcu_shutter.put(0)
+    planner.move_to(BeamlineState.FLUX_MEASUREMENT)
+    # Centre in Y
+    dev.bcu_shutter.put(1)
+    start_inten = read_i1()
+    print(f"Starting intensity is {start_inten:.4g}")
+    d['coll_y'].move('in')
+    coll_start_inten = read_i1()
+    print(f"Intensity with collimator in before Y scan: {coll_start_inten}")
+    go_to_peak(dev.coll_y, dev.i1, -0.4, 0.4, 20, relative=True, confirm=False)
+    coll_end_inten = read_i1()
+    print(f"Intensity with collimator in after Y scan: {coll_end_inten}")
+    # Centre in X
+    coll_start_inten = read_i1()
+    print(f"Intensity with collimator in before X scan: {coll_start_inten}")
+    go_to_peak(dev.coll_x, dev.i1, -0.5, 0.5, 20, relative=True, confirm=False)
+    coll_end_inten = read_i1()
+    print(f"Intensity with collimator in after X scan: {coll_end_inten}")
+    dev.bcu_shutter.put(0)
+
+
+
+
+    
+
+

pxii_bec/macros/mx_basics.bak

@@ -1,313 +0,0 @@
-"""Get data from an h5 file or BEC history and perform fitting."""
-
-import numpy as np
-from lmfit.models import (
-    GaussianModel,
-    LorentzianModel,
-    VoigtModel,
-    ConstantModel,
-    LinearModel,
-)
-from scipy.ndimage import gaussian_filter1d
-import h5py
-import matplotlib.pyplot as plt
-
-
-def create_fit_parameters(
-    deriv: bool = False,
-    negative: bool = False,
-    model: str = "Voigt",
-    baseline: str = "Linear",
-    smoothing: None = None,
-):
-    """Store the fit parameters in a dictionary."""
-    # map input model to lmfit model name
-    model_mappings = {
-        "Gaussian": GaussianModel,
-        "Lorentzian": LorentzianModel,
-        "Voigt": VoigtModel,
-        "Constant": ConstantModel,
-        "Linear": LinearModel,
-    }
-    return {
-        "deriv": deriv,
-        "negative": negative,
-        "model": model_mappings[model],
-        "baseline": model_mappings[baseline],
-        "smoothing": smoothing,
-    }
-
-
-def get_data_from_h5(signal_name: str = "lu_bpmsum"):
-    """Get data from an h5 file."""
-    root, name = get_signal_values(signal_name)
-    with h5py.File("scan_676.h5", "r") as f:
-        entry = f["entry"]["collection"]
-        y_data = entry["devices"][root][name]["value"][:]
-        motor_data = entry["metadata"]["bec"]
-        motor_name = motor_data["scan_motors"][0].decode()
-        scan_number = motor_data["scan_number"][()]
-        x_data = entry["devices"][motor_name][motor_name]["value"][:]
-    return {
-        "x_data": x_data,
-        "y_data": y_data,
-        "signal_name": signal_name,
-        "motor_name": motor_name,
-        "scan_number": str(scan_number),
-    }
-
-def get_history(scan):
-    if isinstance(scan, int):
-        if scan < 0:
-            return bec.history[scan]
-        return bec.history.get_by_scan_number(scan)
-        # return bec.history.get_by_dataset_number(scan)
-    return scan
-
-
-def get_signal_values(signal_name):
-
-    try:
-        mot = getattr(dev,signal_name)
-        root_name = signal_name
-    except:
-        root_name = signal_name.split("_", maxsplit=1)[0]
-    
-    return root_name, signal_name
-
-
-def get_data_from_history(
-    scan: int,
-    signal_name: str = "lu_bpmsum",
-):
-    """Read data from the BEC history and return the X and Y data as arrays."""
-
-    data = get_history(scan)
-    md = data.metadata["bec"]
-    motor_name = md["scan_report_devices"][0].decode()
-    mot_name = motor_name.replace(".","_")
-    print(data.devices)
-    if "." in motor_name:
-        root_name = motor_name.split(".")[0]
-        x_data = data.devices[root_name][mot_name].read()["value"]
-    else:
-        x_data = data.devices[mot_name][mot_name].read()["value"]
-
-    scan_number = md["scan_number"]
-    
-    # y_data = data.devices[signal_name][signal_name].read()["value"]
-    root, name = get_signal_values(signal_name)
-    y_data = data.devices[root][name].read()["value"]
-    return {
-        "signal_name": signal_name,
-        "x_data": x_data,
-        "y_data": y_data,
-        "motor_name": mot_name,
-        "scan_number": float(scan_number),
-    }
-
-
-def process_data(data, fit_params):
-    """
-    Process the signal data for fitting based on derivative or smoothing.
-    """
-    smoothing = fit_params["smoothing"]
-    deriv = fit_params["deriv"]
-    negative = fit_params["negative"]
-
-    signal_name = data["signal_name"]
-    y_data = data["y_data"]
-
-    if smoothing and smoothing > 0.01:
-        fitting_data = gaussian_filter1d(y_data, smoothing)
-        signal_name = f"Smoothed {signal_name}"
-    else:
-        fitting_data = y_data.copy()
-
-    if deriv:
-        fitting_data = np.gradient(fitting_data)
-        signal_name = f"Derivative of {signal_name}"
-
-    if negative:
-        fitting_data = -fitting_data
-        signal_name = f"Negative {signal_name}"
-
-    data.update(
-        {
-            "y_to_fit": fitting_data,
-            "signal_name": signal_name,
-        }
-    )
-    return data
-
-    # if deriv:
-    #     if smoothing:
-    #         y_smooth = gaussian_filter1d(y_data, smoothing)
-    #         fitting_data = np.gradient(y_smooth)
-    #         signal_name = f"Derivative of smoothed {signal_name}"
-    #     else:
-    #         fitting_data = np.gradient(y_data)
-    #         signal_name = f"Derivative of {signal_name}"
-    # elif smoothing and smoothing > 0.01:
-    #     fitting_data = gaussian_filter1d(y_data, smoothing)
-    #     signal_name = f"Smoothed {signal_name}"
-    # else:
-    #     fitting_data = y_data
-
-    # updated_data = {
-    #     "y_to_fit": fitting_data,
-    #     "signal_name": signal_name,
-    # }
-    # data.update(updated_data)
-    # return data
-
-
-def fit(data, fit_params):
-    """Fit a signal to a model and return the fitting results."""
-    # Create the model
-    peak_model = fit_params["model"](prefix="peak_")
-    baseline_model = fit_params["baseline"](prefix="base_")
-    full_model = peak_model + baseline_model
-
-    # Prepare data
-    processed_data = process_data(data, fit_params)
-    params = full_model.make_params()
-    y_min = np.min(processed_data["y_to_fit"])
-
-    # Configure baseline parameters
-    if fit_params["baseline"] == ConstantModel:
-        params["base_c"].set(value=y_min)
-    elif fit_params["baseline"] == LinearModel:
-        params["base_intercept"].set(value=y_min)
-        params["base_slope"].set(value=0)
-
-    # Add peak-specific parameters
-    params.update(
-        peak_model.guess(processed_data["y_to_fit"], x=processed_data["x_data"])
-    )
-
-    # Perform the fitting
-    lmfit_result = full_model.fit(
-        processed_data["y_to_fit"], params, x=processed_data["x_data"]
-    )
-
-    # Find the X that gives the max Y
-    max_index = np.argmax(processed_data["y_to_fit"])
-    x_max = processed_data["x_data"][max_index]
-
-    # Generate data for a smoothed fit curve
-    fit_xdata = np.linspace(np.min(data["x_data"]), np.max(data["x_data"]), 500)
-    fit_ydata = lmfit_result.eval(x=fit_xdata, params=lmfit_result.params)
-
-    # Collect results
-    return {
-        "model": fit_params["model"].__name__,
-        "fwhm": lmfit_result.params["peak_fwhm"].value,
-        "centre": lmfit_result.best_values["peak_center"],
-        "height": lmfit_result.params["peak_height"].value,
-        "chi_sq": lmfit_result.chisqr,
-        "lmfit_result": lmfit_result,
-        "x_max": x_max,
-        "fit_xdata": fit_xdata,
-        "fit_ydata": fit_ydata,
-    }
-
-
-def plot_fitted_data(data, fit_result):
-    """Plot the original data and the fitted model."""
-    plt.plot(data["x_data"], data["y_to_fit"], label="Data")
-    plt.plot(fit_result['fit_xdata'], fit_result['fit_ydata'], label="Fit")
-    plt.xlabel(data["motor_name"])
-    plt.ylabel(data["signal_name"])
-    plt.title(f"Scan {data['scan_number']}, fitted with {fit_result['model']}")
-    plt.grid(True)
-    plt.legend()
-    plt.show()
-
-
-def select_bec_window(dock_area_name="Fitting"):
-    """Check to see if the fitting results dock is already open and re-create it if not"""
-    open_docks = bec.gui.windows
-    if open_docks.get(dock_area_name) is None:
-        dock_area = bec.gui.new(dock_area_name)
-        # wf = dock_area.new("Plot").new(bec.gui.available_widgets.Waveform)
-        wf = dock_area.new(widget='Waveform', object_name='Plot')
-        text_box = dock_area.new(widget='TextBox', object_name="Results", where="bottom")
-    else:
-        wf = bec.gui.Fitting.Plot
-        text_box = bec.gui.Fitting.Results
-    return wf, text_box
-
-
-def plot_live_data_bec(
-    motor_name,
-    signal_name,
-    window_name="Fitting"
-):
-    """
-    Plotting live data for motor and signal using BEC.
-
-    This function plots live data from a specified motor and signal.
-    It clears the current plot window, sets its title, labels the axes
-    with the provided motor and signal names, and initializes live plotting
-    on the given signal against the motor.
-
-    Args:
-        motor_name (str): The name of the motor to be used as the x-axis.
-        signal_name (str): The name of the signal to be used as the y-axis.
-
-    Returns:
-        None
-    """
-    root, name = get_signal_values(signal_name)
-    wf, text_box = select_bec_window(window_name)
-    text_box.set_plain_text("Plotting live data")
-    wf.clear_all()
-    wf.title = "Scan: Live scan"
-    wf.x_label = motor_name
-    wf.y_label = signal_name
-    wf.plot(device_x=motor_name, device_y=root)
-
-
-
-
-def plot_fitted_data_bec(
-    data,
-    fit_result,
-):
-    """
-    Plot fitted data and display fitting parameters in the specified window.
-
-    This function selects a BEC window and plots the original data along with the
-    fitted function. Additionally, it displays the fitting results in a text
-    box within the same window for better visualization of the fit results.
-
-    Parameters:
-    data : dict
-        Dictionary containing the original dataset, where 'x_data' and 'y_to_fit'
-        hold the independent variable and the dependent variable, respectively,
-        'scan_number' represents the scan number, 'motor_name' and 'signal_name'
-        provide axis labels.
-    fit_result : dict
-        Dictionary containing the results of the fit, including parameters such
-        as 'centre', 'fwhm', 'height', and the fitted model stored under
-        'lmfit_result', with its 'best_fit' attribute representing the fitted data.
-    """
-    wf, text_box = select_bec_window()
-    fit_text = (
-        f"Fit parameters: Centre = {fit_result['centre']:.5g}, "
-        f"FWHM = {fit_result['fwhm']:.5f}, "
-        f"Height = {fit_result['height']:.4g}\n"
-        f"Model = {fit_result['model']}\n"
-        f"Chi sq = {fit_result['chi_sq']:.3g}"
-    )
-    text_box.set_plain_text(fit_text)
-    wf.clear_all()
-    wf.title = f"Scan: {data['scan_number']}"
-    wf.x_label = data["motor_name"]
-    wf.y_label = data["signal_name"]
-    wf.plot(x=data["x_data"], y=data["y_to_fit"], label="Data")
-    wf.plot(x=fit_result["fit_xdata"], y=fit_result["fit_ydata"], label="Fit")
-    # wf.Fit.set(symbol_size = 0)
-    wf.get_curve('Fit').set(symbol_size=0)
-

pxii_bec/macros/mx_basics.py

@@ -1,305 +0,0 @@
-"""Get data from an h5 file or BEC history and perform fitting."""
-import re
-import numpy as np
-from lmfit.models import (
-    GaussianModel,
-    LorentzianModel,
-    VoigtModel,
-    ConstantModel,
-    LinearModel,
-)
-from scipy.ndimage import gaussian_filter1d
-import h5py
-import matplotlib.pyplot as plt
-
-
-def create_fit_parameters(
-    deriv: bool = False,
-    negative: bool = False,
-    model: str = "Voigt",
-    baseline: str = "Linear",
-    smoothing: None = None,
-):
-    """Store the fit parameters in a dictionary."""
-    # map input model to lmfit model name
-    model_mappings = {
-        "Gaussian": GaussianModel,
-        "Lorentzian": LorentzianModel,
-        "Voigt": VoigtModel,
-        "Constant": ConstantModel,
-        "Linear": LinearModel,
-    }
-    return {
-        "deriv": deriv,
-        "negative": negative,
-        "model": model_mappings[model],
-        "baseline": model_mappings[baseline],
-        "smoothing": smoothing,
-    }
-
-
-# def get_data_from_h5(signal_name: str = "lu_bpmsum"):
-#     """Get data from an h5 file."""
-#     root, name = get_signal_values(signal_name)
-#     with h5py.File("scan_676.h5", "r") as f:
-#         entry = f["entry"]["collection"]
-#         y_data = entry["devices"][root][name]["value"][:]
-#         motor_data = entry["metadata"]["bec"]
-#         motor_name = motor_data["scan_motors"][0].decode()
-#         scan_number = motor_data["scan_number"][()]
-#         x_data = entry["devices"][motor_name][motor_name]["value"][:]
-#     return {
-#         "x_data": x_data,
-#         "y_data": y_data,
-#         "signal_name": signal_name,
-#         "motor_name": motor_name,
-#         "scan_number": str(scan_number),
-#     }
-
-def get_history(scan):
-    if isinstance(scan, int):
-        if scan < 0:
-            return bec.history[scan]
-        return bec.history.get_by_scan_number(scan)
-        # return bec.history.get_by_dataset_number(scan)
-    return scan
-
-
-def get_value(dev_name, all=False):
-    if dev_name == dev.scam:
-        root = 'scam'
-        live_root = 'scam'
-        name = 'scam_contrast'
-        value = dev_name.read()[name]['value']
-    else:
-        try:
-            name = live_root = root = dev_name.name
-            value = dev_name.read()[name]['value']
-        except:
-            print("non-standard device")
-            root = dev_name.dotted_name
-            name = dev_name.full_name
-            live_root = dev_name.parent.name
-            value = dev_name.read()[name]['value']
- 
-
-    if all:
-        return {"value": value, "root": root, "name": name, "live_root": live_root}
-
-    return value
-
-
-def get_data_from_history(
-    scan: int,
-    signal_dev,
-):
-    """Read data from the BEC history and return the X and Y data as arrays."""
-
-    data = get_history(scan) 
-    md = data.metadata["bec"]
-    motor_name = md["scan_report_devices"][0].decode()
-    if '.' in motor_name:
-        root = motor_name.split(".")[0]
-        name = motor_name.replace(".","_")
-    else:
-        root = motor_name
-        name = motor_name
-
-
-    x_data = data.devices[root][name].read()["value"]
-    y_res = get_value(signal_dev, True)
-    y_data = data.devices[y_res['root']][y_res['name']].read()["value"]
-    signal_name = y_res['name']
-
-    scan_number = md["scan_number"]
-
-    return {
-        "signal_name": signal_name,
-        "x_data": x_data,
-        "y_data": y_data,
-        "motor_name": motor_name,
-        "scan_number": float(scan_number),
-    }
-
-
-def process_data(data, fit_params):
-    """
-    Process the signal data for fitting based on derivative or smoothing.
-    """
-    smoothing = fit_params["smoothing"]
-    deriv = fit_params["deriv"]
-    negative = fit_params["negative"]
-
-    signal_name = data["signal_name"]
-    y_data = data["y_data"]
-
-    if smoothing and smoothing > 0.01:
-        fitting_data = gaussian_filter1d(y_data, smoothing)
-        signal_name = f"Smoothed {signal_name}"
-    else:
-        fitting_data = y_data.copy()
-
-    if deriv:
-        fitting_data = np.gradient(fitting_data)
-        signal_name = f"Derivative of {signal_name}"
-
-    if negative:
-        fitting_data = -fitting_data
-        signal_name = f"Negative {signal_name}"
-
-    data.update(
-        {
-            "y_to_fit": fitting_data,
-            "signal_name": signal_name,
-        }
-    )
-    return data
-
-
-def fit(data, fit_params):
-    """Fit a signal to a model and return the fitting results."""
-    # Create the model
-    peak_model = fit_params["model"](prefix="peak_")
-    baseline_model = fit_params["baseline"](prefix="base_")
-    full_model = peak_model + baseline_model
-
-    # Prepare data
-    processed_data = process_data(data, fit_params)
-    params = full_model.make_params()
-    y_min = np.min(processed_data["y_to_fit"])
-
-    # Configure baseline parameters
-    if fit_params["baseline"] == ConstantModel:
-        params["base_c"].set(value=y_min)
-    elif fit_params["baseline"] == LinearModel:
-        params["base_intercept"].set(value=y_min)
-        params["base_slope"].set(value=0)
-
-    # Add peak-specific parameters
-    params.update(
-        peak_model.guess(processed_data["y_to_fit"], x=processed_data["x_data"])
-    )
-
-    # Perform the fitting
-    lmfit_result = full_model.fit(
-        processed_data["y_to_fit"], params, x=processed_data["x_data"]
-    )
-
-    # Find the X that gives the max Y
-    max_index = np.argmax(processed_data["y_to_fit"])
-    x_max = processed_data["x_data"][max_index]
-
-    # Generate data for a smoothed fit curve
-    fit_xdata = np.linspace(np.min(data["x_data"]), np.max(data["x_data"]), 500)
-    fit_ydata = lmfit_result.eval(x=fit_xdata, params=lmfit_result.params)
-
-    # Collect results
-    return {
-        "model": fit_params["model"].__name__,
-        "fwhm": lmfit_result.params["peak_fwhm"].value,
-        "centre": lmfit_result.best_values["peak_center"],
-        "height": lmfit_result.params["peak_height"].value,
-        "chi_sq": lmfit_result.chisqr,
-        "lmfit_result": lmfit_result,
-        "x_max": x_max,
-        "fit_xdata": fit_xdata,
-        "fit_ydata": fit_ydata,
-    }
-
-
-def plot_fitted_data(data, fit_result):
-    """Plot the original data and the fitted model."""
-    plt.plot(data["x_data"], data["y_to_fit"], label="Data")
-    plt.plot(fit_result['fit_xdata'], fit_result['fit_ydata'], label="Fit")
-    plt.xlabel(data["motor_name"])
-    plt.ylabel(data["signal_name"])
-    plt.title(f"Scan {data['scan_number']}, fitted with {fit_result['model']}")
-    plt.grid(True)
-    plt.legend()
-    plt.show()
-
-
-def select_bec_window(dock_area_name="Fitting"):
-    """Check to see if the fitting results dock is already open and re-create it if not"""
-    open_docks = bec.gui.windows
-    if open_docks.get(dock_area_name) is None:
-        dock_area = bec.gui.new(dock_area_name)
-        wf = dock_area.new(widget='Waveform', object_name='Plot')
-        text_box = dock_area.new(widget='TextBox', object_name="Results", where="bottom")
-    else:
-        wf = bec.gui.Fitting.Plot
-        text_box = bec.gui.Fitting.Results
-    return wf, text_box
-
-
-def plot_live_data_bec(
-    motor_name,
-    signal_name,
-    window_name="Fitting"
-):
-    """
-    Plotting live data for motor and signal using BEC.
-
-    This function plots live data from a specified motor and signal.
-    It clears the current plot window, sets its title, labels the axes
-    with the provided motor and signal names, and initializes live plotting
-    on the given signal against the motor.
-
-    Args:
-        motor_name (str): The name of the motor to be used as the x-axis.
-        signal_name (str): The name of the signal to be used as the y-axis.
-
-    Returns:
-        None
-    """
-    wf, text_box = select_bec_window(window_name)
-    text_box.set_plain_text("Plotting live data")
-    wf.clear_all()
-    wf.title = "Scan: Live scan"
-    wf.x_label = motor_name
-    wf.y_label = signal_name
-    wf.plot(device_x=motor_name, device_y=signal_name)
-
-
-
-
-def plot_fitted_data_bec(
-    data,
-    fit_result,
-):
-    """
-    Plot fitted data and display fitting parameters in the specified window.
-
-    This function selects a BEC window and plots the original data along with the
-    fitted function. Additionally, it displays the fitting results in a text
-    box within the same window for better visualization of the fit results.
-
-    Parameters:
-    data : dict
-        Dictionary containing the original dataset, where 'x_data' and 'y_to_fit'
-        hold the independent variable and the dependent variable, respectively,
-        'scan_number' represents the scan number, 'motor_name' and 'signal_name'
-        provide axis labels.
-    fit_result : dict
-        Dictionary containing the results of the fit, including parameters such
-        as 'centre', 'fwhm', 'height', and the fitted model stored under
-        'lmfit_result', with its 'best_fit' attribute representing the fitted data.
-    """
-    wf, text_box = select_bec_window()
-    fit_text = (
-        f"Fit parameters: Centre = {fit_result['centre']:.5g}, "
-        f"FWHM = {fit_result['fwhm']:.5f}, "
-        f"Height = {fit_result['height']:.4g}\n"
-        f"Model = {fit_result['model']}\n"
-        f"Chi sq = {fit_result['chi_sq']:.3g}"
-    )
-    text_box.set_plain_text(fit_text)
-    wf.clear_all()
-    wf.title = f"Scan: {data['scan_number']}"
-    wf.x_label = data["motor_name"]
-    wf.y_label = data["signal_name"]
-    wf.plot(x=data["x_data"], y=data["y_to_fit"], label="Data")
-    wf.plot(x=fit_result["fit_xdata"], y=fit_result["fit_ydata"], label="Fit")
-    # wf.Fit.set(symbol_size = 0)
-    wf.get_curve('Fit').set(symbol_size=0)
-

pxii_bec/macros/mxcomm_analysis.py

@@ -0,0 +1,315 @@
+"""Functions for fitting scan data"""
+
+from dataclasses import dataclass
+import numpy as np
+from lmfit.models import GaussianModel, LorentzianModel, VoigtModel, ConstantModel, LinearModel
+from scipy.ndimage import gaussian_filter1d
+
+
+@dataclass(frozen=True)
+class FitDefaults:
+    """Default values for fitting routines"""
+
+    # Constants for default models, baselines, and parameters
+    MODEL = "Voigt"
+    BASELINE = "Linear"
+
+
+def create_fit_parameters(
+    deriv: bool = False,
+    negative: bool = False,
+    model: str = "Voigt",
+    baseline: str = "Linear",
+    smoothing: None = None,
+):
+    """Store the fit parameters in a dictionary
+
+    Parameters: optional derivative, negative peak, fit model,
+    baseline model, optional smoothing
+    Returns: dictionary of fit parameters
+    """
+    # map input model to lmfit model name
+    model_mappings = {
+        "Gaussian": GaussianModel,
+        "Lorentzian": LorentzianModel,
+        "Voigt": VoigtModel,
+        "Constant": ConstantModel,
+        "Linear": LinearModel,
+    }
+    return {
+        "deriv": deriv,
+        "negative": negative,
+        "model": model_mappings[model],
+        "baseline": model_mappings[baseline],
+        "smoothing": smoothing,
+    }
+
+
+def process_data(data, fit_params):
+    """
+    Process the signal data for fitting based on derivative or smoothing
+
+    Parameters: data and fitting parameter dictionary
+    Returns: processed data
+    """
+    smoothing = fit_params["smoothing"]
+    deriv = fit_params["deriv"]
+    negative = fit_params["negative"]
+
+    signal_name = data["signal_name"]
+    y_data = data["y_data"]
+
+    if smoothing and smoothing > 0.01:
+        fitting_data = gaussian_filter1d(y_data, smoothing)
+        signal_name = f"Smoothed {signal_name}"
+    else:
+        fitting_data = y_data.copy()
+
+    if deriv:
+        fitting_data = np.gradient(fitting_data)
+        signal_name = f"Derivative of {signal_name}"
+
+    if negative:
+        fitting_data = -fitting_data
+        signal_name = f"Negative {signal_name}"
+
+    data.update({"y_to_fit": fitting_data, "signal_name": signal_name})
+    return data
+
+
+def fit(data, fit_params):
+    """Fit a signal to a model and return the fitting results
+
+    Parameters: data and fitting parameters dictionary
+    Returns: Fit results and quality indicators
+    """
+
+    # Create the model
+    peak_model = fit_params["model"](prefix="peak_")
+    baseline_model = fit_params["baseline"](prefix="base_")
+    full_model = peak_model + baseline_model
+
+    # Prepare data
+    processed_data = process_data(data, fit_params)
+    params = full_model.make_params()
+    y_min = np.min(processed_data["y_to_fit"])
+
+    # Configure baseline parameters
+    if fit_params["baseline"] == ConstantModel:
+        params["base_c"].set(value=y_min)
+    elif fit_params["baseline"] == LinearModel:
+        params["base_intercept"].set(value=y_min)
+        params["base_slope"].set(value=0)
+
+    # Add peak-specific parameters
+    params.update(peak_model.guess(processed_data["y_to_fit"], x=processed_data["x_data"]))
+
+    # Perform the fitting
+    lmfit_result = full_model.fit(processed_data["y_to_fit"], params, x=processed_data["x_data"])
+
+    # Find the X that gives the max Y
+    max_index = np.argmax(processed_data["y_to_fit"])
+    x_max = processed_data["x_data"][max_index]
+
+    # Generate data for a smoothed fit curve
+    fit_xdata = np.linspace(np.min(data["x_data"]), np.max(data["x_data"]), 500)
+    fit_ydata = lmfit_result.eval(x=fit_xdata, params=lmfit_result.params)
+
+    # Fit diagnostics
+    residuals = lmfit_result.residual
+    noise = np.std(residuals)
+
+    height = lmfit_result.params["peak_height"].value
+    centre = lmfit_result.best_values["peak_center"]
+    fwhm = lmfit_result.params["peak_fwhm"].value
+
+    centre_stderr = lmfit_result.params["peak_center"].stderr
+    height_stderr = lmfit_result.params["peak_height"].stderr
+
+    snr = height / noise if noise > 0 else np.inf
+
+    # Collect results
+    return {
+        "model": fit_params["model"].__name__,
+        "fwhm": fwhm,
+        "centre": centre,
+        "height": height,
+        # Fit quality
+        "chi_sq": lmfit_result.chisqr,
+        "redchi": lmfit_result.redchi,
+        "aic": lmfit_result.aic,
+        "bic": lmfit_result.bic,
+        "noise": noise,
+        "snr": snr,
+        "centre_stderr": centre_stderr,
+        "height_stderr": height_stderr,
+        # Existing information
+        "lmfit_result": lmfit_result,
+        "x_max": x_max,
+        "fit_xdata": fit_xdata,
+        "fit_ydata": fit_ydata,
+    }
+
+
+def check_peak_fit(x_data, result):
+    """Check whether a fitted peak is reliable
+    Parameters: X data, results from the fit step
+    Returns: success or problems"""
+
+    x_data = np.asarray(x_data)
+
+    x_min = np.min(x_data)
+    x_max = np.max(x_data)
+    scan_range = x_max - x_min
+    step = np.median(np.diff(np.sort(x_data)))
+
+    centre = result["centre"]
+    fwhm = result["fwhm"]
+
+    problems = []
+
+    if not x_min + 0.1 * scan_range < centre < x_max - 0.1 * scan_range:
+        problems.append("peak centre is too close to scan edge")
+
+    if fwhm < 3 * step:
+        problems.append("peak is narrower than 3 scan steps")
+
+    if fwhm > 0.7 * scan_range:
+        problems.append("peak is too broad for scan range")
+
+    if result["snr"] < 5:
+        problems.append(f"low SNR ({result['snr']:.1f})")
+
+    centre_stderr = result["centre_stderr"]
+
+    if centre_stderr is None:
+        problems.append("centre uncertainty could not be estimated")
+    elif centre_stderr > 0.2 * fwhm:
+        problems.append(f"centre uncertainty is large " f"({centre_stderr:.3g})")
+
+    return {"success": not problems, "problems": problems}
+
+
+def find_edge_offset(x, y, fraction=0.5, offset=0.0):
+    """
+    Find a single signal edge and calculate a target position
+    at a known offset from that edge.
+    Parameters: xdata, ydata, fraction of edge, description, offset
+    Returns: edge, centre, threshold
+    """
+
+    x = np.asarray(x)
+    y = np.asarray(y)
+
+    y_min = np.min(y)
+    y_max = np.max(y)
+
+    threshold = y_min + fraction * (y_max - y_min)
+
+    # Find threshold crossings
+    crossings = np.where(np.diff(np.sign(y - threshold)) != 0)[0]
+
+    if len(crossings) == 0:
+        raise RuntimeError("No edge found in scan")
+
+    # If we expect only one useful edge, use the strongest/first one
+    i = crossings[0]
+
+    # Linear interpolation for more accurate edge position
+    x1, x2 = x[i], x[i + 1]
+    y1, y2 = y[i], y[i + 1]
+
+    edge = x1 + (threshold - y1) * (x2 - x1) / (y2 - y1)
+
+    centre = edge + offset
+
+    return {"edge": edge, "centre": centre, "threshold": threshold}
+
+
+def find_feature_centre(x, y, feature="peak", fraction=0.5):
+    """
+    Find the centre of a peak-like or dip-like region from its two edges.
+
+    Parameters
+    ----------
+    x, y : array-like
+        Scan positions and measured signal.
+
+    mode : {"peak", "dip"}
+        "peak":
+            Feature is the high-signal region, e.g. collimator scan.
+
+        "dip":
+            Feature is the low-signal region, e.g. beamstop scan.
+
+    fraction : float
+        Fraction of the signal range used for the edge threshold.
+        0.5 means the 50% crossing points.
+
+    Returns
+    -------
+    dict
+        centre
+        left_edge
+        right_edge
+        width
+        threshold
+    """
+
+    x = np.asarray(x, dtype=float)
+    y = np.asarray(y, dtype=float)
+
+    y_low = np.min(y)
+    y_high = np.max(y)
+
+    threshold = y_low + fraction * (y_high - y_low)
+
+    if feature == "peak":
+        inside = y >= threshold
+
+    elif feature == "dip":
+        inside = y <= threshold
+
+    else:
+        raise ValueError("mode must be 'peak' or 'dip'")
+
+    # Find transitions into and out of the feature
+    transitions = np.diff(inside.astype(int))
+
+    entering = np.where(transitions == 1)[0]
+    leaving = np.where(transitions == -1)[0]
+
+    if len(entering) == 0 or len(leaving) == 0:
+        raise ValueError("Could not find both feature edges")
+
+    # Use first entrance and last exit
+    left_idx = entering[0]
+    right_idx = leaving[-1]
+
+    left_edge = interpolate_crossing(
+        x[left_idx], y[left_idx], x[left_idx + 1], y[left_idx + 1], threshold
+    )
+
+    right_edge = interpolate_crossing(
+        x[right_idx], y[right_idx], x[right_idx + 1], y[right_idx + 1], threshold
+    )
+
+    centre = 0.5 * (left_edge + right_edge)
+    width = right_edge - left_edge
+
+    return {
+        "centre": centre,
+        "left_edge": left_edge,
+        "right_edge": right_edge,
+        "width": width,
+        "threshold": threshold,
+    }
+
+
+def interpolate_crossing(x1, y1, x2, y2, threshold):
+    """Linearly interpolate the x position where y crosses threshold."""
+
+    if y2 == y1:
+        return 0.5 * (x1 + x2)
+
+    return x1 + (threshold - y1) * (x2 - x1) / (y2 - y1)

pxii_bec/macros/mxcomm_data.py

@@ -0,0 +1,99 @@
+"""Get data from scan number, BEC history or H5 files"""
+
+import h5py
+
+# def get_data_from_h5(signal_name: str = "lu_bpmsum"):
+#     """Get data from an h5 file."""
+#     root, name = get_signal_values(signal_name)
+#     with h5py.File("scan_676.h5", "r") as f:
+#         entry = f["entry"]["collection"]
+#         y_data = entry["devices"][root][name]["value"][:]
+#         motor_data = entry["metadata"]["bec"]
+#         motor_name = motor_data["scan_motors"][0].decode()
+#         scan_number = motor_data["scan_number"][()]
+#         x_data = entry["devices"][motor_name][motor_name]["value"][:]
+#     return {
+#         "x_data": x_data,
+#         "y_data": y_data,
+#         "signal_name": signal_name,
+#         "motor_name": motor_name,
+#         "scan_number": str(scan_number),
+#     }
+
+
+def get_history(scan):
+    """Get scan container from the bec history
+    Parameters: scan number or relative number
+    Returns: scan container
+
+    """
+
+    if isinstance(scan, int):
+        if scan < 0:
+            return bec.history[scan]
+        return bec.history.get_by_scan_number(scan)
+        # return bec.history.get_by_dataset_number(scan)
+    return scan
+
+
+def get_data_from_history(scan: int, signal_dev):
+    """Get scan data from the BEC history
+    Parameters: scan number or relative number
+    Returns:    the X and Y data as arrays
+                signal name, motor name and scan number
+
+    """
+
+    data = get_history(scan)
+    md = data.metadata["bec"]
+    motor_name = md["scan_report_devices"][0].decode()
+    if "." in motor_name:
+        root = motor_name.split(".")[0]
+        name = motor_name.replace(".", "_")
+    else:
+        root = motor_name
+        name = motor_name
+
+    x_data = data.devices[root][name].read()["value"]
+    y_res = get_value(signal_dev, True)
+    y_data = data.devices[y_res["root"]][y_res["name"]].read()["value"]
+    signal_name = y_res["name"]
+
+    scan_number = md["scan_number"]
+
+    return {
+        "signal_name": signal_name,
+        "x_data": x_data,
+        "y_data": y_data,
+        "motor_name": motor_name,
+        "scan_number": float(scan_number),
+    }
+
+
+def get_value(dev_name, details=False):
+    """Get the value of a device
+    Parameters: device, details: bool
+    Returns: value or full details
+
+    """
+
+    if dev_name == dev.scam:
+        root = "scam"
+        live_root = "scam"
+        name = "scam_contrast"
+        value = dev_name.read()[name]["value"]
+    else:
+        try:
+            name = live_root = root = dev_name.name
+            value = dev_name.read()[name]["value"]
+        except:
+            print("non-standard device")
+            root = dev_name.dotted_name
+            name = dev_name.full_name
+            live_root = dev_name.parent.name
+            value = dev_name.read()[name]["value"]
+
+    if details:
+        return {"value": value, "root": root, "name": name, "live_root": live_root}
+
+    return value

pxii_bec/macros/mxcomm_plot.py

@@ -0,0 +1,127 @@
+"""Plotting routines for MX scans and fits in BEC."""
+
+
+def select_bec_window(window_name="Fitting"):
+    """
+    Return the plot and results widgets for a BEC window.
+
+    If the requested window does not already exist, it is created with
+    a Waveform plot and a TextBox for results.
+    """
+
+    open_docks = bec.gui.windows
+
+    if open_docks.get(window_name) is None:
+        dock_area = bec.gui.new(window_name)
+
+        plot_widget = dock_area.new(widget="Waveform", object_name="Plot")
+
+        text_widget = dock_area.new(widget="TextBox", object_name="Results", where="bottom")
+
+    else:
+        dock_area = open_docks[window_name]
+        plot_widget = dock_area.Plot
+        text_widget = dock_area.Results
+
+    return plot_widget, text_widget
+
+
+def plot_live_data(
+    motor_name, signal_name, window_name="Fitting", plot_widget=None, title="Live scan"
+):
+    """
+    Start live plotting of a signal against a motor in BEC.
+
+    If plot_widget is supplied, the data are plotted in that widget.
+    Otherwise, the standard BEC fitting window is used.
+    """
+    if isinstance(plot_widget, str):
+        dock_area = bec.gui.windows[window_name]
+        wf = getattr(dock_area, plot_widget)
+
+    elif plot_widget is not None:
+        wf = plot_widget
+    else:
+        wf, text_widget = select_bec_window(window_name)
+        text_widget.set_plain_text("Plotting live data")
+
+    wf.clear_all()
+    wf.title = title
+    wf.x_label = motor_name
+    wf.y_label = signal_name
+
+    wf.plot(device_x=motor_name, device_y=signal_name)
+
+
+def plot_fitted_data(data, fit_result, window_name="Fitting", plot_widget=None, title=None):
+    """
+    Plot completed scan data and its fitted curve in BEC.
+
+    If plot_widget is supplied, the fit is plotted in that widget.
+    Otherwise, the standard BEC fitting window is used and the fitting
+    results are displayed in its results text box.
+    """
+    if isinstance(plot_widget, str):
+        dock_area = bec.gui.windows[window_name]
+        wf = getattr(dock_area, plot_widget)
+
+    elif plot_widget is not None:
+        wf = plot_widget
+
+    else:
+        wf, text_widget = select_bec_window(window_name)
+
+        fit_text = (
+            f"Centre = {fit_result['centre']:.5g}\n"
+            f"FWHM = {fit_result['fwhm']:.5f}\n"
+            f"Height = {fit_result['height']:.4g}\n"
+            f"Model = {fit_result['model']}\n"
+            f"Chi sq = {fit_result['chi_sq']:.3g}\n"
+            f"SNR = {fit_result['snr']:.3g}"
+        )
+
+        text_widget.set_plain_text(fit_text)
+
+    wf.clear_all()
+    wf.title = f"Scan: {data['scan_number']}"
+    wf.x_label = data["motor_name"]
+    wf.y_label = data["signal_name"]
+
+    wf.plot(x=data["x_data"], y=data["y_to_fit"], label="Data")
+
+    wf.plot(x=fit_result["fit_xdata"], y=fit_result["fit_ydata"], label="Fit")
+
+    # Display the fitted curve as a line without symbols.
+    wf.get_curve("Fit").set(symbol_size=0)
+
+def plot_scan_data(
+    data,
+    window_name="Fitting",
+    plot_widget=None,
+    title=None,
+):
+    """Plot completed scan data without a fitted curve."""
+
+    if isinstance(plot_widget, str):
+        dock_area = bec.gui.windows[window_name]
+        wf = getattr(dock_area, plot_widget)
+
+    elif plot_widget is not None:
+        wf = plot_widget
+
+    else:
+        wf, text_widget = select_bec_window(window_name)
+        text_widget.set_plain_text("Fit rejected")
+
+    # Clear the live plot/subscription
+    wf.clear_all()
+
+    wf.title = title or f"Scan: {data['scan_number']}"
+    wf.x_label = data["motor_name"]
+    wf.y_label = data["signal_name"]
+
+    wf.plot(
+        x=data["x_data"],
+        y=data["y_data"],
+        label="Data",
+    )

pxii_bec/macros/mxcomm_scan.py

@@ -1,20 +1,38 @@
-"""Use the methods in mx_basics to perform:
+"""Custom scan routines:
 1) a go_to_peak scan, that scans a motor, finds the peak position and moves to peak
-2) fits data from a bec history file
+2) a find feature scan, that finds a non-Gaussian peak or dip or an edge
 """
 
 from dataclasses import dataclass
 import numpy as np
 
-# from pxiii_parameters import FitDefaults, BPMScans, MirrorConfig
 
-# from mx_basics import (
+# from mxcomm_analysis import (
 #     create_fit_parameters,
 #     get_data_from_history,
 #     fit,
-#     plot_fitted_data_bec,
+#     FitDefaults,
+#     check_peak_fit,
+#     get_value,
+#     find_feature_centre,
+#     find_edge_offset)
+# from mxcomm_plot import (
 #     plot_live_data_bec,
+#     plot_live_data,
 # )
+# from pxii_bec import (
+#     scans,
+#     umv,
+#     dev
+# )
+
+@dataclass(frozen=True)
+class ScanDefaults:
+    """Default values for fitting routines"""
+
+    # Preferred scan parameters
+    SETTLE_TIME = 0.1
+    RELATIVE_MODE = True
 
 
 # Method functions
@@ -39,7 +57,7 @@ def calculate_step_size(start: float, stop: float, steps: int) -> float:
     return round((stop - start) / (steps - 1), 5)
 
 
-def move_to_position(motor_device, motor_name: str, position: float, data: dict):
+def move_to_scan_position(motor_device, motor_name: str, position: float, data: dict):
     """
     Function to move a specified motor device to a given position.
 
@@ -82,17 +100,9 @@ def move_to_position(motor_device, motor_name: str, position: float, data: dict)
     scans.umv(motor_device, motor_position, relative=False)
     print(f"\n Moving {motor_name} to position {motor_position: .3f}")
 
-
-@dataclass(frozen=True)
-class FitDefaults:
-    """Default values for fitting routines"""
-
-    # Constants for default models, baselines, and parameters
-    MODEL = "Voigt"
-    BASELINE = "Linear"
-    SETTLE_TIME = 0.1
-    RELATIVE_MODE = True
-
+def restore_position(motor_device, position):
+    """Return a motor to its position before the scan."""
+    scans.umv(motor_device, position, relative=False)    
 
 def go_to_peak(
     motor_device,
@@ -100,12 +110,15 @@ def go_to_peak(
     start: float,
     stop: float,
     steps: int,
-    relative: bool = FitDefaults.RELATIVE_MODE,
+    relative: bool = ScanDefaults.RELATIVE_MODE,
     plot: bool = True,
-    settle: float = FitDefaults.SETTLE_TIME,
+    settle: float = ScanDefaults.SETTLE_TIME,
     confirm: bool = True,
     gomax: bool = False,
-    negative: bool = False
+    negative: bool = False,
+    gap: bool = False,
+    window_name: str = "Fitting",
+    plot_widget=None,
 ):
     """
     Go to the peak of a signal by scanning a motor within a specified range and
@@ -131,29 +144,40 @@ def go_to_peak(
             operations such as plotting, fitting, or motor movement.
 
     Returns:
-        None
+        fit parameters
     """
     # Get motor and signal names
     mres = get_value(motor_device, True)
-    motor_name = mres['root']
+    motor_name = mres["root"]
     sres = get_value(signal_device, True)
-    signal_name = sres['root']
+    signal_name = sres["root"]
     print(f"Signal name is {signal_name}")
-    # wf.plot(x_name=motor_name, y_name=signal_name)
+
+    # if plot:
+    #     if plot_widget is not None:
+    #         plot_live_data(plot_widget, motor_name, signal_name, title=f"{motor_name} scan")
+    #     else:
+    #         plot_live_data_bec(motor_name, signal_name, window_name=window_name)
+
     if plot:
-        plot_live_data_bec(motor_name, signal_name)
+        plot_live_data(
+            motor_name, 
+            signal_name,
+            window_name=window_name,
+            plot_widget=plot_widget, 
+            title=f"{motor_name} scan")
 
     # Validate and calculate step size
     step_size = calculate_step_size(start, stop, steps)
 
     # Confirm the scan range
-    current_motor_position = mres['value']
-    print(f"Starting motor position = {current_motor_position}")
+    starting_position = mres["value"]
+    print(f"Starting motor position = {starting_position}")
 
     if confirm:
         if relative:
-            scan_start = current_motor_position + start
-            scan_end = current_motor_position + stop
+            scan_start = starting_position + start
+            scan_end = starting_position + stop
             print(
                 f"\nScanning from {scan_start: .6g} to {scan_end: .6g} in "
                 f"{steps} steps of size {step_size}"
@@ -165,15 +189,15 @@ def go_to_peak(
             )
             print(f"Relative mode = {relative}")
         input("Press Enter to continue...")
-    
+
     # Perform the scan
     scan_result = scans.line_scan(
         motor_device, start, stop, steps=steps, relative=relative, settling_time=settle
     )
-    motor_data = scan_result.scan.live_data[mres['live_root']][mres['name']].val
+    motor_data = scan_result.scan.live_data[mres["live_root"]][mres["name"]].val
     # # motor_data = scan_result.scan.live_data[motor_name][motor_name].val
     # sres = get_value(signal_name, True)
-    signal_data = scan_result.scan.live_data[sres['root']][sres['name']].val
+    signal_data = scan_result.scan.live_data[sres["root"]][sres["name"]].val
     scan_number = "Current"
 
     data = {
@@ -186,24 +210,64 @@ def go_to_peak(
     }
 
     # Define and fit model to scan data
-    fit_params = create_fit_parameters(deriv = False,
-                                       negative = False, 
-                                       model = FitDefaults.MODEL, 
-                                       baseline = FitDefaults.BASELINE)
+    fit_params = create_fit_parameters(
+        deriv=False, negative=False, model=FitDefaults.MODEL, baseline=FitDefaults.BASELINE
+    )
     fit_result = fit(data, fit_params)
 
+    check = check_peak_fit(data["x_data"], fit_result)
+
+    if not check["success"]:
+        print("Peak fit rejected:")
+        for problem in check["problems"]:
+            print(f"  - {problem}")
+
+        if plot:
+            plot_scan_data(
+                data,
+                window_name=window_name,
+                plot_widget=plot_widget,
+                title=f"{motor_name} scan - FIT REJECTED",
+            )
+        print(f"Returning motor to start position, {starting_position:.5g}")
+        restore_position(
+            motor_device,
+            starting_position,
+        )
+
+        return None
+
+
+
+
     # Plot the fitted data if plot = True
     if plot:
-        plot_fitted_data_bec(data, fit_result)
+        plot_fitted_data(data, fit_result, window_name=window_name, plot_widget=plot_widget)
 
     # If gomax is set then move to the maximum value, rather than the fit centre
     if gomax:
         value = fit_result["x_max"]
         print(f"Max position is at {value}")
-        move_to_position(data["motor_device"], data["motor_name"], fit_result["x_max"], data)
+        move_scan_to_position(data["motor_device"], data["motor_name"], fit_result["x_max"], data)
     else:
-        # Safely move the motor to the peak position
-        move_to_position(data["motor_device"], data["motor_name"], fit_result["centre"], data)
+        if gap:
+            if fit_result["centre"] < 4.5:
+                umv(dev.id_gap, 4.5)
+            else:
+                move_to_scan_position(
+                    data["motor_device"], data["motor_name"], fit_result["centre"], data
+                )
+        else:
+            # Safely move the motor to the peak position
+            move_to_scan_position(data["motor_device"], data["motor_name"], fit_result["centre"], data)
+
+    return {
+        "centre": fit_result["centre"],
+        "fwhm": float(fit_result["fwhm"]),
+        "height": float(fit_result["height"]),
+        "chi_sq": float(fit_result["chi_sq"]),
+        "snr": float(fit_result["snr"]),
+    }
 
 
 def fit_history(
@@ -212,7 +276,7 @@ def fit_history(
     deriv: bool = False,
     negative: bool = False,
     smoothing: bool = False,
-    model: str = FitDefaults.MODEL,
+    model: str = "Voigt",
     move_to_peak: bool = False,
 ):
     """
@@ -244,8 +308,7 @@ def fit_history(
 
     # Perform fit and plot the data
     fit_result = fit(data, fit_params)
-    plot_fitted_data_bec(data, fit_result)
-    
+    plot_fitted_data(data, fit_result)
 
     # Optionally move the motor to the peak position
     if move_to_peak:
@@ -253,85 +316,124 @@ def fit_history(
     return fit_result
 
 
-def scan_bpm(bpmname):
+def find_feature(
+    motor_device,
+    signal_device,
+    start: float,
+    stop: float,
+    steps: int,
+    feature: str = "peak",
+    offset=None,
+    fraction: float = 0.5,
+    relative: bool = ScanDefaults.RELATIVE_MODE,
+    plot: bool = True,
+    settle: float = ScanDefaults.SETTLE_TIME,
+    confirm: bool = True,
+    window_name: str = "Fitting",
+    plot_widget=None,
+):
     """
-    Runs a grid scan of a BPM in x and y,  and plots each channel
-    as a heatmap.
-
-    Parameters:
-        bpmname: the name of the bpm to be scanned e.g. "fe"
-
+    Scan a motor, determine the centre of a high- or low-signal region,
+    and move the motor to that position.
     """
 
-    # Open a dock area and set up the heatmaps
-    dock_area = bec.gui.new("XBPM_Scan")
-    wf5 = dock_area.new("Sum").new(bec.gui.available_widgets.Heatmap)
-    wf1 = dock_area.new("Ch1", relative_to="Sum", position="bottom").new(
-        bec.gui.available_widgets.Heatmap
+    # Get motor and signal names
+    mres = get_value(motor_device, True)
+    motor_name = mres["root"]
+
+    sres = get_value(signal_device, True)
+    signal_name = sres["root"]
+
+    print(f"Signal name is {signal_name}")
+
+    # Calculate step size
+    step_size = calculate_step_size(start, stop, steps)
+
+    current_motor_position = mres["value"]
+    print(f"Starting motor position = {current_motor_position}")
+
+    # Confirm scan range
+    if confirm:
+        if relative:
+            scan_start = current_motor_position + start
+            scan_end = current_motor_position + stop
+        else:
+            scan_start = start
+            scan_end = stop
+
+        print(
+            f"\nScanning from {scan_start:.6g} to {scan_end:.6g} "
+            f"in {steps} steps of size {step_size}"
+        )
+        print(f"Relative mode = {relative}")
+        print(f"Looking for a {feature} feature")
+
+        input("Press Enter to continue...")
+
+    # -------------------------------------------------
+    # Set up live plot BEFORE starting the scan
+    # -------------------------------------------------
+    if plot:
+        if plot_widget is not None:
+            plot_live_data(plot_widget, motor_name, signal_name, title=f"{motor_name} scan")
+    else:
+        plot_live_data_bec(motor_name, signal_name, window_name=window_name)
+
+    # Perform scan
+    scan_result = scans.line_scan(
+        motor_device, start, stop, steps=steps, relative=relative, settling_time=settle
     )
-    wf3 = dock_area.new("Ch3", relative_to="Ch1", position="right").new(
-        bec.gui.available_widgets.Heatmap
-    )
-    wf4 = dock_area.new("Ch4", relative_to="Ch3", position="bottom").new(
-        bec.gui.available_widgets.Heatmap
-    )
-    wf2 = dock_area.new("Ch2", relative_to="Ch1", position="bottom").new(
-        bec.gui.available_widgets.Heatmap
-    )
-    wfscan = dock_area.new("ScanControl").new(bec.gui.available_widgets.ScanControl)
 
-    cfg = getattr(BPMScans, bpmname)
+    # Retrieve data
+    motor_data = scan_result.scan.live_data[mres["live_root"]][mres["name"]].val
 
-    wf1.x_label = cfg["x_name"]
-    wf1.y_label = cfg["y_name"]
-    wf1.plot(x_name=cfg["x_name"], y_name=cfg["y_name"], z_name=cfg["z1_name"], color_map="plasma")
+    signal_data = scan_result.scan.live_data[sres["root"]][sres["name"]].val
 
-    wf2.x_label = cfg["x_name"]
-    wf2.y_label = cfg["y_name"]
-    wf2.plot(x_name=cfg["x_name"], y_name=cfg["y_name"], z_name=cfg["z2_name"], color_map="plasma")
+    data = {
+        "x_data": np.array(motor_data),
+        "y_data": np.array(signal_data),
+        "motor_name": motor_name,
+        "signal_name": signal_name,
+        "motor_device": motor_device,
+        "scan_number": "Current",
+    }
 
-    wf3.x_label = cfg["x_name"]
-    wf3.y_label = cfg["y_name"]
-    wf3.plot(x_name=cfg["x_name"], y_name=cfg["y_name"], z_name=cfg["z3_name"], color_map="plasma")
+    # Find feature centre
+    if feature in ("peak", "dip"):
+        result = find_feature_centre(
+            data["x_data"], data["y_data"], feature=feature, fraction=fraction
+        )
+    elif feature == "edge_offset":
+        result = find_edge_offset(data["x_data"], data["y_data"], offset=offset)
+    else:
+        raise ValueError(f"Unknown feature type {feature}")
+    if feature in ("peak", "dip"):
+        print(f"\nLeft edge  = {result['left_edge']:.4f}")
+        print(f"Right edge = {result['right_edge']:.4f}")
+        print(f"Width      = {result['width']:.4f}")
+        print(f"Centre     = {result['centre']:.4f}")
+    elif feature == "edge_offset":
+        print(f"Edge = {result['edge']:.4f}")
+        print(f"Centre     = {result['centre']:.4f}")
 
-    wf4.x_label = cfg["x_name"]
-    wf4.y_label = cfg["y_name"]
-    wf4.plot(x_name=cfg["x_name"], y_name=cfg["y_name"], z_name=cfg["z4_name"], color_map="plasma")
+    # Move safely to centre
+    move_to_position(data["motor_device"], data["motor_name"], result["centre"], data)
 
-    wf5.x_label = cfg["x_name"]
-    wf5.y_label = cfg["y_name"]
-    wf5.plot(x_name=cfg["x_name"], y_name=cfg["y_name"], z_name=cfg["z5_name"], color_map="plasma")
-    # Run the scan
-    x_mot = cfg["x_device"]
-    y_mot = cfg["y_device"]
-    # scans.grid_scan(x_mot, -0.5, 0.5, 20, y_mot, -0.5, 0.5, 20,
-    #                     exp_time=0.5, relative=False, snaked=True)
+    if plot and plot_widget is not None:
+        # Replace the live subscription with static completed scan data
+        plot_widget.clear_all()
 
+        plot_widget.plot(x=data["x_data"], y=data["y_data"], label="Scan data")
+        # Mark centre
 
-def optimise_kb(mirror):
-    """
-    Runs a grid scan of a the upstream and downstream benders,
-    and plots a heatmap of the sample camera x or y sigma.
+        plot_widget.plot(
+            x=[result["centre"], result["centre"]],
+            y=[np.min(data["y_data"]), np.max(data["y_data"])],
+            label="Centre",
+        )
 
-    Parameters:
-        mirror: either "hfm" or :vfm"
+        plot_widget.title = f"{motor_name} scan"
+        plot_widget.x_label = motor_name
+        plot_widget.y_label = signal_name
 
-    """
-
-    # Open a dock area and set up the heatmaps
-    dock_area = bec.gui.new(mirror)
-    wf1 = dock_area.new("Heatmap").new(bec.gui.available_widgets.Heatmap)
-
-    wfscan = dock_area.new("ScanControl").new(bec.gui.available_widgets.ScanControl)
-
-    cfg = getattr(MirrorConfig, mirror)
-
-    wf1.x_label = cfg["bu_name"]
-    wf1.y_label = cfg["bd_name"]
-    wf1.plot(x_name=cfg["bu_name"], y_name=cfg["bd_name"], z_name=cfg["z_name"], color_map="plasma")
-
-    # Run the scan
-    x_mot = cfg["x_device"]
-    y_mot = cfg["y_device"]
-    # scans.grid_scan(x_mot, -0.02, 0.02, 11, y_mot, -0.02, 0.02, 11,
-    #                     exp_time=0.5, relative=True, snaked=True)
+    return result

pxii_bec/macros/position_device.py

@@ -6,129 +6,6 @@ from typing import Callable, List, Dict, Optional, Union
 import time
 
 
-class SimpleStatus:
-    """Makes a mock motor return a status"""
-
-    def __init__(self, motor, target, delay=0.0, success=True, name=""):
-        self.motor = motor
-        self.target = target
-        self.delay = delay
-        self._success = success
-        self.name = name
-        self._done = False
-
-    def wait(self, timeout=None):
-        start = time.time()
-
-        while True:
-            # simulate motion completion
-            if not self._done:
-                if time.time() - start >= self.delay:
-                    if self._success:
-                        self.motor.position = self.target
-                    self._done = True
-
-            if self._done:
-                if not self._success:
-                    raise RuntimeError(f"Motor {self.name} failed")
-                return True
-
-            if timeout is not None and (time.time() - start) > timeout:
-                raise TimeoutError(f"Timeout waiting for {self.name}")
-
-            time.sleep(0.01)
-
-
-class MotorAdapter:
-    """Motor adapter for setting up mock/real motors"""
-
-    def move(self, pos: float):
-        """Move the motor to the given position"""
-        raise NotImplementedError
-
-    def move_with_status(self, pos: float):
-        """Move the motor to the given position with a status"""
-        raise NotImplementedError
-
-    def set_fail(self, value: bool):
-        """Put the motor into a failure state"""
-        raise NotImplementedError
-
-    @property
-    def actual(self) -> float:
-        """The actual position of the motor"""
-        raise NotImplementedError
-
-
-class MockMotorAdapter(MotorAdapter):
-    """Motor adapter for mock motors"""
-
-    def __init__(self, name: str):
-        self._motor = PositionDevice.MockMotor(name)
-
-    def move(self, pos: float):
-        """Move the motor to the given position"""
-        self._motor.move(pos)
-
-    def move_with_status(self, pos: float):
-        """Move the motor to the given position with a status"""
-        if self._motor.fail:
-            return SimpleStatus(
-                motor=self._motor,
-                target=pos,
-                delay=self._motor.delay,
-                success=False,
-                name=self._motor.name,
-            )
-
-        # don't update position immediately
-        return SimpleStatus(
-            motor=self._motor,
-            target=pos,
-            delay=self._motor.delay,
-            success=True,
-            name=self._motor.name,
-        )
-
-    def set_fail(self, value: bool):
-        """Put the motor into a failure state"""
-        self._motor.fail = value
-
-    def set_delay(self, value: float):
-        self._motor.delay = value
-
-    @property
-    def actual(self) -> float:
-        """The actual position of the motor"""
-        return self._motor.position
-
-
-class RealMotorAdapter(MotorAdapter):
-    """Motor adapter for real motors"""
-
-    def __init__(self, mot, name):
-        self._motor = mot
-        self._name = name
-
-    def move(self, pos: float):
-        """Move the motor to the given position"""
-        scans.umv(self._motor, pos, relative=False)
-
-    def move_with_status(self, pos: float):
-        """Move the motor to the given position with a status"""
-        return scans.mv(self._motor, pos, relative=False)
-        # scans.mv(self._motor, pos, relative=False)
-
-    def set_fail(self, value: bool):
-        """Put the motor into a failure state"""
-        raise NotImplementedError
-
-    @property
-    def actual(self) -> float:
-        """The actual position of the motor"""
-        return self._motor.read()[self._name]["value"]
-
-
 @dataclass
 class PositionDevice:
     """Generic device that moves between named or numeric positions"""
@@ -141,15 +18,9 @@ class PositionDevice:
     guards: List[Callable[[], None]] = field(default_factory=list)
     policy: Optional[Callable[[float], None]] = None
     allow_arbitrary: bool = False
-    use_mock: bool = True
+    move_method: str = "motor"
 
     def __post_init__(self):
-        if self.use_mock:
-            self.mot = MockMotorAdapter(self.bec_name)
-        else:
-            self.mot = RealMotorAdapter(self.mot_device, self.bec_name)
-            # self.mot = getattr(dev, self.bec_name)
-
         # Normalize position names
         self.positions = {k.lower(): v for k, v in self.positions.items()}
 
@@ -185,7 +56,10 @@ class PositionDevice:
         if self.policy:
             self.policy(pos)
 
-        self.mot.move(pos)
+        if self.move_method == "signal":
+            self.mot_device.put(pos)
+        else:
+            scans.umv(self.mot_device, pos, relative=False)
 
     def mv(self, target: Union[str, float]):
         """move devices with a timeout"""
@@ -197,20 +71,19 @@ class PositionDevice:
         if self.policy:
             self.policy(pos)
 
-        # status.wait(self.device_timeout)
-        return self.mot.move_with_status(pos)
-        # self.mot.move_with_status(pos)
+        return scans.mv(self.mot_device, pos, relative=False)
 
-    def set_position(self, target: Union[str, float]):
-        """Only to be used for testing purposes, bypasses guards"""
-        pos = self._resolve_target(target)
-        self.mot.move(pos)
+
+    # def set_position(self, target: Union[str, float]):
+    #     """Only to be used for testing purposes, bypasses guards"""
+    #     pos = self._resolve_target(target)
+    #     self.mot.move(pos)
 
     @property
     def actual(self) -> float:
         """Return the actual position of the device."""
         # return self.mot.read()[self.bec_name]["value"]
-        return self.mot.actual
+        return self.mot_device.read()[self.bec_name]["value"]
 
     @property
     def pos(self) -> str:
@@ -230,25 +103,4 @@ class PositionDevice:
         pos = self._resolve_target(target)
         return abs(self.actual - pos) <= self.tol
 
-    # -------------------------
-    # Mock Motor Implementation
-    # -------------------------
-
-    class MockMotor:
-        """Mock motor implementation"""
-
-        def __init__(self, name: str):
-            self.name = name
-            self.position = 0.0
-            self._target = self.position
-            self.fail = False
-            self.delay = 0
-
-        #
-        def move(self, pos: float):
-            """Move the motor to the given position."""
-            if self.fail:
-                return
-                # raise RuntimeError(f"Motor {self.name} failed")
-            time.sleep(self.delay)
-            self.position = pos
+ 

pxii_bec/macros/pxii_energy.py

@@ -78,7 +78,7 @@ def set_mirror_stripe(energy_ev):
     print(f"Selected mirror stripe: {selected_stripe}")
 
 
-def mono_pitch_scan(plot=True):
+def mono_pitch_scan(plot=True, window_name="Fitting", plot_widget=None):
     """Scan the monochromator pitch and move to the peak."""
     # Move to the calculated pitch value for the current energy
     print("Starting Mono Pitch Scan.")
@@ -94,7 +94,7 @@ def mono_pitch_scan(plot=True):
     # Go to peak using default parameters from EnergyDefaults
     if plot:
         print("Scanning monochromator pitch and moving to peak, with plotting.")
-        go_to_peak(
+        result = go_to_peak(
             EnergyDefaults.mono_pitch,
             EnergyDefaults.signals["sig1"],
             -EnergyDefaults.pitch_scan["halfwidth"],
@@ -104,10 +104,12 @@ def mono_pitch_scan(plot=True):
             settle=0.01,
             plot=True,
             confirm=False,
+            window_name=window_name,
+            plot_widget=plot_widget
         )
     else:
         print("Scanning monochromator pitch and moving to peak, without plotting.")
-        go_to_peak(
+        result = go_to_peak(
             EnergyDefaults.mono_pitch,
             EnergyDefaults.signals["sig1"],
             -EnergyDefaults.pitch_scan["halfwidth"],
@@ -119,6 +121,8 @@ def mono_pitch_scan(plot=True):
             confirm=False,
         )
 
+    return result
+
 
 # Specific functions - need to be edited for each beamline
 
@@ -212,7 +216,12 @@ def move_dcm_motors(energy_ev):
     )
 
 
-def bl_energy(energy_ev, move_gap=True, mono_scan=True, plot=True):
+def bl_energy(energy_ev, 
+              move_gap=True, 
+              mono_scan=True, 
+              plot=True, 
+              window_name="Fitting",
+              plot_widget=None):
     """
     Adjusts the beamline's energy to the specified energy in electron volts (eV).
     The function validates the target energy, checks the current energy, and makes
@@ -251,37 +260,87 @@ def bl_energy(energy_ev, move_gap=True, mono_scan=True, plot=True):
     # Step 4: Perform DCM pitch scan and move to peak.
     if mono_scan:
         if plot:
-            mono_pitch_scan(plot=True)
+            result = mono_pitch_scan(plot=True,
+                                     window_name=window_name,
+                                     plot_widget=plot_widget)
+            return result
         else:
             mono_pitch_scan(plot=False)
 
-def new_calc_gap(energy_ev):
-    """
-    Calculates the gap based on the given energy
-    """
-    gapval = Gap()
-    harmonic_values = gapval.get_harmonic_values(energy_ev)
-    if harmonic_values is not None:
-        fit = np.poly1d(harmonic_values)(energy_ev)
-        gap_value = fit * 0.9994
-        return gap_value
-    return None
+    
 
-def set_gap(energy_ev, print_result=True):
-    """
-    set the gap
-    """
-    gap = new_calc_gap(energy_ev)
-    # else:
-    #     gap = calc_gap(energy_kev)
-    if gap is not None and Harmonics.min_gap_value <= gap <= 20:
-        umv(dev.id_gap, gap)
-        if print_result:
-            print(f"Setting gap to {gap:.3f} mm")
-    else:
-        msg = (
-            f"Position {gap: .3f} mm is outside the allowed scan range of "
-            f"4.5 to 20 mm. Not changing the gap."
+# def new_calc_gap(energy_ev):
+#     """
+#     Calculates the gap based on the given energy
+#     """
+#     gapval = Gap()
+#     harmonic_values = gapval.get_harmonic_values(energy_ev)
+#     if harmonic_values is not None:
+#         fit = np.poly1d(harmonic_values)(energy_ev)
+#         gap_value = fit * 0.9994
+#         return gap_value
+#     return None
+
+# def set_gap(energy_ev, print_result=True):
+#     """
+#     set the gap
+#     """
+#     gap = new_calc_gap(energy_ev)
+#     # else:
+#     #     gap = calc_gap(energy_kev)
+#     if gap is not None and Harmonics.min_gap_value <= gap <= 20:
+#         umv(dev.id_gap, gap)
+#         if print_result:
+#             print(f"Setting gap to {gap:.3f} mm")
+#     else:
+#         umv(dev.id_gap, Harmonics.min_gap_value)
+#         # msg = (
+#         #     f"Position {gap: .3f} mm is outside the allowed scan range of "
+#         #     f"4.5 to 20 mm. Not changing the gap."
+#         # )
+#         # raise ValueError(msg)
+#     return gap
+
+def new_calc_gap(energy_ev, print_result=True):
+    """Calculate the undulator gap for the given energy."""
+
+    harmonic_values = Gap.get_harmonic_values(energy_ev)
+
+    if harmonic_values is None:
+        return None
+
+    gap = np.poly1d(harmonic_values)(energy_ev)
+    gap *= Gap.correction_factor
+
+    if print_result:
+        harmonic = Gap.get_harmonic_by_energy(energy_ev)
+        print(
+            f"Calculated gap is {gap:.3f} mm "
+            f"using {harmonic}"
         )
-        raise ValueError(msg)
+
     return gap
+
+def set_gap(energy_ev, print_result=False):
+    """Set the undulator gap."""
+
+    gap = new_calc_gap(energy_ev)
+
+    if gap is not None and Gap.minimum_gap <= gap <= 20:
+        umv(dev.id_gap, gap)
+
+        if print_result:
+            harmonic = Gap.get_harmonic_by_energy(energy_ev)
+            print(
+                f"Setting gap to {gap:.3f} mm "
+                f"using {harmonic}"
+            )
+
+    else:
+      raise ValueError(
+        f"Cannot calculate a valid gap for {energy_ev:.0f} eV. "
+        f"Calculated gap: {gap}"
+    )
+
+    return gap
+

pxii_bec/macros/pxii_parameters.py

@@ -2,7 +2,7 @@
 
 from dataclasses import dataclass
 import numpy as np
-
+from enum import Enum
 
 @dataclass(frozen=True)
 class EnergyDefaults:
@@ -18,7 +18,7 @@ class EnergyDefaults:
     mono_fpitch = dev.dcm_fpitch
     mono_perp = dev.dcm_perp
     mono_roll = dev.dcm_froll
-    mono_roll_value = 4.56
+    mono_roll_value = 4.75
     LUT_table = "luts/energy_lut.csv"
     stripe_thresholds = {"silicon": 9000, "rhodium": 20000, "platinum": 40000}
     pitch_scan = {"halfwidth": 0.075, "steps": 20}
@@ -29,72 +29,120 @@ class Calibration:
     """Calibration parameters for PXII optics"""
 
     # pitch = np.array([4.61823701e-14, -1.97330772e-09, 2.89694543e-05, -5.34468669e00])
-    pitch_calib = [-5.89030403e-19, 
-                   4.85665619e-14, 
-                   -1.49364459e-09, 
-                   2.08652765e-05,
-                   -5.59843072e+00]
+    # pitch_calib = [-5.89030403e-19, 
+    #                4.85665619e-14, 
+    #                -1.49364459e-09, 
+    #                2.08652765e-05,
+    #                -5.59843072e+00]
+    pitch_calib = [-4.09746242e-18,
+                   2.32320520e-13,
+                   -4.95650128e-09,
+                   4.85168715e-05,
+                   -5.67418203e+00]
+    
     # roll = np.array([2.28291039e-03, -2.41928101e01])
 
 
+from dataclasses import dataclass
+import numpy as np
+
+
 @dataclass(frozen=True)
 class Gap:
-    """Fit parameters to calculate gap from harmonics"""
+    """Fit parameters to calculate undulator gap from energy."""
+
+    minimum_gap = 4.5
+    maximum_gap = 20.0
+    correction_factor = 0.9994
 
     harmonics = {
-        "H3": np.array([9.409e-04, 2.797e-01]),
-        "H5": np.array([5.097e-04, 8.065e-01]),
-        "H7": np.array([3.623e-04, 8.252e-01]),
-        "H9": np.array([2.795e-04, 8.595e-01]),
-        "H11": np.array([2.309e-04, 8.167e-01]),
-        "H13": np.array([1.928e-04, 8.748e-01]),
-        "H15": np.array([1.67e-04, 8.743e-01]),
-        "H17": np.array([1.484e-04, 8.546e-01]),
-        "H19": np.array([1.275e-04, 1.0085]),
-    }
-    # Define harmonic ranges as a constant
-    harmonic_ranges = {
-        "H3": (4900, 7000),
-        "H5": (7000, 10000),
-        "H7": (10000, 13000),
-        "H9": (13000, 16000),
-        "H11": (16000, 19000),
-        "H13": (19000, 22000),
-        "H15": (22000, 25000),
-        "H17": (25000, 27500),
-        "H19": (27500, float("inf")),
+        "H3": np.array([9.57400000e-04, 1.62100000e-01]),
+        "H5": np.array([5.03757576e-04, 8.53992727e-01]),
+        "H7": np.array([3.57624242e-04, 8.74460000e-01]),
+        "H9": np.array([2.76218182e-04, 9.00565455e-01]),
+        "H11": np.array([2.25662626e-04, 9.04630909e-01]),
+        "H13": np.array([1.89440404e-04, 9.35049697e-01]),
+        "H15": np.array([1.64179798e-04, 9.35255758e-01]),
+        "H17": np.array([1.44460317e-04, 9.47441667e-01]),
+        "H19": np.array([1.26071429e-04, 1.04085357e00]),
     }
 
-    @staticmethod
-    def get_harmonic_by_energy(energy_ev: float):
+    @classmethod
+    def minimum_energy(cls, harmonic):
         """
-        Determines the harmonic key based on the provided energy.
-
-        Args:
-            energy_ev (float): The energy value (eV).
-
-        Returns:
-            Optional[str]: The harmonic name (e.g., 'H3', 'H7')
-            if the range matches, None otherwise.
+        Return the lowest energy that can be reached with this harmonic
+        without the corrected gap going below the minimum gap.
         """
-        for harmonic, (low, high) in Gap.harmonic_ranges.items():
-            if low < energy_ev <= high:
-                return harmonic
-        return None
+        slope, intercept = cls.harmonics[harmonic]
 
-    def get_harmonic_values(self, energy_ev: float):
+        # new_calc_gap() applies the correction factor afterwards,
+        # so account for it here as well.
+        raw_minimum_gap = cls.minimum_gap / cls.correction_factor
+
+        return (raw_minimum_gap - intercept) / slope
+
+    @classmethod
+    def get_harmonic_by_energy(cls, energy_ev):
         """
-        Retrieves the harmonic array based on the energy value.
+        Return the appropriate harmonic for the requested energy.
 
-        Args:
-            energy_ev (float): The energy value (eV).
-
-        Returns:
-            Optional[np.array]: The corresponding array of harmonic values
-            if the range matches, None otherwise.
+        Switch to the next harmonic as soon as it can produce the
+        requested energy without requiring a gap below minimum_gap.
         """
-        harmonic = self.get_harmonic_by_energy(energy_ev)
-        return self.harmonics.get(harmonic) if harmonic else None
+        selected = None
+
+        for harmonic in cls.harmonics:
+            if energy_ev >= cls.minimum_energy(harmonic):
+                selected = harmonic
+            else:
+                break
+
+        return selected
+
+    @classmethod
+    def get_harmonic_values(cls, energy_ev):
+        """
+        Return the fit coefficients for the appropriate harmonic.
+        """
+        harmonic = cls.get_harmonic_by_energy(energy_ev)
+
+        if harmonic is None:
+            return None
+
+        return cls.harmonics[harmonic]
+
+    @classmethod
+    def calculate_gap(cls, energy_ev):
+        """
+        Calculate the corrected undulator gap for the requested energy.
+
+        Returns
+        -------
+        tuple
+            (harmonic, gap)
+        """
+        harmonic = cls.get_harmonic_by_energy(energy_ev)
+
+        if harmonic is None:
+            raise ValueError(
+                f"{energy_ev:.0f} eV is below the minimum calibrated energy."
+            )
+
+        slope, intercept = cls.harmonics[harmonic]
+
+        gap = slope * energy_ev + intercept
+        gap *= cls.correction_factor
+
+        if not cls.minimum_gap <= gap <= cls.maximum_gap:
+            raise ValueError(
+                f"Calculated gap {gap:.3f} mm for {energy_ev:.0f} eV "
+                f"using {harmonic} is outside the allowed range "
+                f"{cls.minimum_gap:.1f}-{cls.maximum_gap:.1f} mm."
+            )
+
+        return harmonic, gap
+    
+
 
 
 @dataclass(frozen=True)
@@ -126,72 +174,60 @@ class CamConversion:
     a = 0.5208
     b = 0.002586
 
+# @dataclass(frozen=True)
+class BPM_X_Calibration:
+    coeffs = (0.11079256, -0.02493623, -0.00799871, -0.00482772,  0.16569134,  0.44192518)
+    signal_limit = 0.90
+    centred_position = 0.441
+    crosshair_signal = 0.000731
+    crosshair_bpm_position = 0.444
 
-@dataclass(frozen=True)
-class BPMScans:
-    """Define the names of the motors and bpm channels"""
-
-    # fe = {
-    #     "x_name": dev.fe_bpm_x.name,
-    #     "y_name": dev.fe_bpm_y.name,
-    #     "z1_name": dev.fe_bpm1.name,
-    #     "z2_name": dev.fe_bpm2.name,
-    #     "z3_name": dev.fe_bpm3.name,
-    #     "z4_name": dev.fe_bpm3.name,
-    #     "z5_name": dev.fe_bpmsum.name,
-    #     "x_device": dev.fe_bpm_x,
-    #     "y_device": dev.fe_bpm_y,
-    # }
-    lu = {
-        "x_name": dev.lu_bpm_x.name,
-        "y_name": dev.lu_bpm_y.name,
-        "z1_name": dev.lu_bpm1.name,
-        "z2_name": dev.lu_bpm2.name,
-        "z3_name": dev.lu_bpm3.name,
-        "z4_name": dev.lu_bpm4.name,
-        "z5_name": dev.lu_bpmsum.name,
-        "x_device": dev.lu_bpm_x,
-        "y_device": dev.lu_bpm_y,
-    }
-    ss = {
-        "x_name": dev.ss_bpm_x.name,
-        "y_name": dev.ss_bpm_y.name,
-        "z1_name": dev.ss_bpm1.name,
-        "z2_name": dev.ss_bpm2.name,
-        "z3_name": dev.ss_bpm3.name,
-        "z4_name": dev.ss_bpm4.name,
-        "z5_name": dev.ss_bpmsum.name,
-        "x_device": dev.ss_bpm_x,
-        "y_device": dev.ss_bpm_y,
-    }
-    bcu = {
-        "x_name": dev.bcu_bpm_x.name,
-        "y_name": dev.bcu_bpm_y.name,
-        "z1_name": dev.bcu_bpm1.name,
-        "z2_name": dev.bcu_bpm2.name,
-        "z3_name": dev.bcu_bpm3.name,
-        "z4_name": dev.bcu_bpm4.name,
-        "z5_name": dev.bcu_bpmsum.name,
-        "x_device": dev.bcu_bpm_x,
-        "y_device": dev.bcu_bpm_y,
-    }
 
 
 # @dataclass(frozen=True)
-# class MirrorConfig:
-#     """Define the names of the mirror channels"""
+class BPM_Y_Calibration:
+    coeffs = (-0.14645507, 0.010016, 0.08258468, -0.00547233, -0.02715706, -1.54419084)
+    signal_limit = 0.90
+    centred_position = -1.540
+    crosshair_signal = -0.01383
+    crosshair_bpm_position = -1.540
+
+# @dataclass(frozen=True)
+class Beamstop:
+    x_offset = -0.5
+
+
+class CameraZoom(Enum):
+    X1 = ("1x", 1)
+    X2 = ("2x", 280)
+    X3 = ("3x", 500)
+    X5_8 = ("5.8x", 700)
+    X7_5 = ("7.5x", 800)
+    X12_5 = ("12.5x", 1000)
+
+    def __init__(self, label, position):
+        self.label= label
+        self.position = position
+
+    @classmethod
+    def from_position(cls, position, tolerance=2):
+        """Return zoom level according to motor position"""
+        for zoom in cls:
+            if abs(position - zoom.position) <= tolerance:
+                return zoom
+        raise ValueError(
+            f"No camera zoom levle found near motor position {position}"
+        )
+
+    @classmethod
+    def from_label(cls, label):
+        """Return zoom level corresponding to zoom label"""
+        for zoom in cls:
+            if zoom.label.lower() == label.lower():
+                return zoom
+        raise ValueError(f"Unknown camera zoom: {label}")
+
+
+
+
 
-#     hfm = {
-#         "bu_name": dev.hfm_bu.name,
-#         "bd_name": dev.hfm_bd.name,
-#         "z_name": dev.samcam_xsig.name,
-#         "x_device": dev.hfm_bu,
-#         "y_device": dev.hfm_bd,
-#     }
-#     vfm = {
-#         "bu_name": dev.vfm_bu.name,
-#         "bd_name": dev.vfm_bd.name,
-#         "z_name": dev.samcam_ysig.name,
-#         "x_device": dev.vfm_bu,
-#         "y_device": dev.vfm_bd,
-#     }

pxii_bec/macros/set_kbox.sh

@@ -46,24 +46,27 @@ main() {
     echo "Device: $DEVICE_NAME"
     echo "Value : $SET_VALUE"
 
-    # --- Your logic here ---
-    # Example placeholder:
 
     if [[ $DEVICE_NAME == "colli_in" ]]; then
 	echo "caput X10SA-ES-COL:POS-SET-SEQ.DO2 $SET_VALUE"
+	caput X10SA-ES-COL:POS-SET-SEQ.DO2 $SET_VALUE
     fi
     if [[ $DEVICE_NAME == "colli_out" ]]; then
 	echo "caput X10SA-ES-COL:POS-SET-SEQ.DO1 $SET_VALUE"
+	caput X10SA-ES-COL:POS-SET-SEQ.DO1 $SET_VALUE
     fi
     # 
     if [[ $DEVICE_NAME == "scinti_in" ]]; then
 	echo "caput X10SA-ES-SCL:POS-SET-SEQ.DO2 $SET_VALUE"
+	caput X10SA-ES-SCL:POS-SET-SEQ.DO2 $SET_VALUE
     fi
     if [[ $DEVICE_NAME == "diode_in" ]]; then
 	echo "caput X10SA-ES-SCL:POS-SET-SEQ.DO3 $SET_VALUE"
+	caput X10SA-ES-SCL:POS-SET-SEQ.DO3 $SET_VALUE
     fi
     if [[ $DEVICE_NAME == "scinti_out" || $DEVICE_NAME == "diode_out" ]]; then
 	echo "caput X10SA-ES-SCL:POS-SET-SEQ.DO1 $SET_VALUE"
+	caput X10SA-ES-SCL:POS-SET-SEQ.DO1 $SET_VALUE
     fi
     # 
     echo "Setting device '$DEVICE_NAME' to '$SET_VALUE'..."

pxii_bec/macros/test.py

@@ -1,31 +0,0 @@
-"""Test state changes"""
-
-import random
-import time
-
-def state_test():
-    d, planner = init_beamline_environment()
-    tested_states = []
-    states_to_test = [BeamlineState.BEAMSTOP_ALIGNMENT,
-                      BeamlineState.DATA_COLLECTION,
-                      BeamlineState.MAINTENANCE,
-                      BeamlineState.ROBOT_SAMPLE_EXCHANGE,
-                      BeamlineState.SAMPLE_ALIGNMENT,
-                    #   BeamlineState.DC_XRF,
-                      BeamlineState.FLUX_MEASUREMENT,
-                      BeamlineState.XTAL_SNAPSHOT,
-                      BeamlineState.MANUAL_SAMPLE_EXCHANGE,
-                      BeamlineState.BEAM_VISUALISATION]
-    previous = None
-    for i in range(20):
-        choices = [s for s in states_to_test if s != previous]
-        current = random.choice(choices)
-        print(f"Testing {current}")
-        planner.move_to(current)
-        time.sleep(2)
-        tested_states.append(current.name)
-        previous = current
-    planner.move_to(BeamlineState.MANUAL_SAMPLE_EXCHANGE)
-
-    print("\n".join(tested_states))
-

pxii_bec/macros/test_state_transitions.py

@@ -0,0 +1,55 @@
+"""Test state changes"""
+
+import random
+import sys
+import time
+from pathlib import Path
+
+sys.path.insert(0, str(Path(__file__).parent.parent / "scripts"))
+from pxii_bec.scripts import beamline_context as bl
+
+
+def state_test(number_tries=20):
+    tested_states = []
+    states_to_test = [
+        BeamlineState.BEAMSTOP_ALIGNMENT,
+        # BeamlineState.DATA_COLLECTION,
+        BeamlineState.MAINTENANCE,
+        BeamlineState.ROBOT_SAMPLE_EXCHANGE,
+        BeamlineState.SAMPLE_ALIGNMENT,
+        # BeamlineState.DC_XRF,
+        # BeamlineState.FLUX_MEASUREMENT,
+        # BeamlineState.XTAL_SNAPSHOT,
+        BeamlineState.MANUAL_SAMPLE_EXCHANGE,
+        BeamlineState.BEAM_VISUALISATION,
+    ]
+    # BeamlineState.DET_DIODE]
+    previous = None
+    for i in range(number_tries):
+        choices = [s for s in states_to_test if s != previous]
+        current = random.choice(choices)
+        print(f"Testing {current}")
+        bl.planner.move_to(current)
+        time.sleep(2)
+        tested_states.append(current.name)
+        previous = current
+    bl.planner.move_to(BeamlineState.MANUAL_SAMPLE_EXCHANGE)
+
+    print("\n".join(tested_states))
+
+
+def beamstop_test():
+    # time.sleep(10)
+    # bl.d['bs_pos'].move('out')
+    # time.sleep(5)
+    # bl.d['bs_pos'].move('in')
+    # for i in range(6):
+    #     bl.planner.move_to(BeamlineState.ROBOT_SAMPLE_EXCHANGE)
+    #     bl.planner.move_to(BeamlineState.BEAMSTOP_ALIGNMENT)
+    #     time.sleep(1)
+    # bl.planner.move_to(BeamlineState.MANUAL_SAMPLE_EXCHANGE)
+    # bl.planner.move_to(BeamlineState.BEAM_VISUALISATION)
+    time.sleep(10)
+    bl.d["bl_pos"].move("out")
+    time.sleep(10)
+    bl.d["bl_pos"].move("in")

pxii_bec/scripts/beamline_context.py

@@ -0,0 +1,23 @@
+"""Shared beamline context.
+Provides access to the current beamline environment,
+positioned devices, and beamline state planner."""
+
+from .init_beamline import init_beamline_environment
+
+env = None
+d = None
+planner = None
+
+
+def initialise():
+    """Initialise the beamline environment."""
+    global env, d, planner
+
+    env, d, planner = init_beamline_environment()
+
+    print("Beamline environment, devices and planner have been initialised")
+
+
+def reload():
+    """Re-initialise the beamline environment"""
+    initialise()

pxii_bec/scripts/init_beamline.py

@@ -0,0 +1,106 @@
+"""Initialise sample environment devices and beamline states."""
+
+import os
+import importlib
+from dataclasses import dataclass
+from pathlib import Path
+
+
+@dataclass
+class Environment:
+    """Define configuration files for a beamline and BEC deployment."""
+
+    beamline: str
+
+    beamline_configs = {
+        "x06sa": {
+            "package": "pxi_bec",
+            "devices": "pxi-devices.yaml",
+            "states": "pxi-beamline-states.yaml",
+        },
+        "x10sa": {
+            "package": "pxii_bec",
+            "devices": "pxii-devices.yaml",
+            "states": "pxii-beamline-states.yaml",
+        },
+        "x06da": {
+            "package": "pxiii_bec",
+            "devices": "pxiii-devices.yaml",
+            "states": "pxiii-beamline-states.yaml",
+        },
+    }
+
+    valid_deployments = ("test", "production")
+
+    def __post_init__(self):
+
+        if self.beamline not in self.beamline_configs:
+            raise ValueError(
+                f"Unknown beamline '{self.beamline}'. " f"Choose from {list(self.beamline_configs)}"
+            )
+
+        cfg = self.beamline_configs[self.beamline]
+
+        # Import whichever deplyment BEV has already made available
+        package = importlib.import_module(cfg["package"])
+
+        self.bec_dir = Path(package.__file__).resolve().parent
+
+        self.device_config_dir = self.bec_dir / "device_configs"
+        self.lut_dir = self.bec_dir / "macros" / "luts"
+
+        self.devices_file = self.device_config_dir / cfg["devices"]
+        self.states_file = self.device_config_dir / cfg["states"]
+
+        # Determine which deployment is being used
+        self.deployment = self._get_deployment()
+
+        if not self.devices_file.exists():
+            raise FileNotFoundError(f"Device configuration not found: {self.devices_file}")
+
+        if not self.states_file.exists():
+            raise FileNotFoundError(f"State configuration not found: {self.states_file}")
+
+    def _get_deployment(self):
+        """Determine deployment from the loaded package path."""
+
+        for deployment in ("test", "production"):
+            if deployment in self.bec_dir.parts:
+                return deployment
+
+        raise RuntimeError(f"Could not determine BEC deployment from {self.bec_dir}")
+
+
+def init_beamline_environment():
+    """
+    Initialise beamline positioned devices and state planner.
+    """
+
+    beamline = os.environ["BEAMLINE"]
+
+    env = Environment(
+        beamline=beamline,
+    )
+
+    print("=" * 60)
+    print(f"Beamline       : {beamline}")
+    print(f"BEC deployment : {env.deployment.upper()}")
+    print(f"Devices file   : {env.devices_file}")
+    print(f"States file    : {env.states_file}")
+    print("=" * 60)
+
+    # Initialise positioned devices
+    device_manager = DeviceManager()
+    devices = device_manager.initialize_devices(env.devices_file)
+
+    # Initialise beamline states
+    state_manager = DefineStatesManager()
+    states, allow_modifiers = state_manager.initialize_states(env.states_file)
+
+    # Set up dependencies
+    deps = planner_deps()
+
+    # Set up planner
+    planner = StateChangePlanner(devices, states, allow_modifiers, deps)
+
+    return env, devices, planner