From e33538b5ea4aa8f4f806f3cd3b417c0002822542 Mon Sep 17 00:00:00 2001 From: x10sa Date: Tue, 25 Aug 2026 17:48:13 +0200 Subject: [PATCH 1/2] beamline staff changes --- .../startup/post_startup.py | 5 +- .../bec_ipython_client/startup/pre_startup.py | 8 + pxii_bec/device_configs/pxii-devices.yaml | 16 +- pxii_bec/macros/APscripts.py | 41 +- pxii_bec/macros/Undu_Helpers.py | 28 +- pxii_bec/macros/bcu_bpm_positions.py | 84 +++ pxii_bec/macros/beamline_planner.py | 1 + pxii_bec/macros/beampos_with_en.py~ | 9 - pxii_bec/macros/build_devices.py | 119 +--- pxii_bec/macros/calibrate_bpm.py | 338 ++++++--- pxii_bec/macros/calibrate_bpm.py.bak | 670 ------------------ pxii_bec/macros/dependencies.py | 4 +- pxii_bec/macros/devices_manager.py | 7 +- pxii_bec/macros/flux_macros.py | 298 ++++++++ ...{init_beamline.py => init_beamline.py.bak} | 0 pxii_bec/macros/katscripts.py | 424 ++++++++--- .../macros/{mx_basics.py => mx_basics.py.bak} | 229 +++++- .../{mx_methods.py => mx_methods.py.bak} | 341 ++++++++- pxii_bec/macros/mxcomm_analysis.py | 315 ++++++++ pxii_bec/macros/mxcomm_data.py | 99 +++ pxii_bec/macros/mxcomm_plot.py | 127 ++++ pxii_bec/macros/mxcomm_scan.py | 439 ++++++++++++ pxii_bec/macros/position_device.py | 174 +---- pxii_bec/macros/pxii_energy.py | 123 +++- pxii_bec/macros/pxii_parameters.py | 270 ++++--- pxii_bec/macros/set_kbox.sh | 7 +- pxii_bec/macros/test.py | 31 - pxii_bec/macros/test_state_transitions.py | 55 ++ pxii_bec/scripts/__init__.py | 0 pxii_bec/scripts/beamline_context.py | 23 + pxii_bec/scripts/init_beamline.py | 106 +++ 31 files changed, 3065 insertions(+), 1326 deletions(-) create mode 100644 pxii_bec/macros/bcu_bpm_positions.py delete mode 100644 pxii_bec/macros/beampos_with_en.py~ delete mode 100644 pxii_bec/macros/calibrate_bpm.py.bak create mode 100644 pxii_bec/macros/flux_macros.py rename pxii_bec/macros/{init_beamline.py => init_beamline.py.bak} (100%) rename pxii_bec/macros/{mx_basics.py => mx_basics.py.bak} (61%) rename pxii_bec/macros/{mx_methods.py => mx_methods.py.bak} (57%) create mode 100644 pxii_bec/macros/mxcomm_analysis.py create mode 100644 pxii_bec/macros/mxcomm_data.py create mode 100644 pxii_bec/macros/mxcomm_plot.py create mode 100755 pxii_bec/macros/mxcomm_scan.py delete mode 100755 pxii_bec/macros/test.py create mode 100755 pxii_bec/macros/test_state_transitions.py create mode 100644 pxii_bec/scripts/__init__.py create mode 100644 pxii_bec/scripts/beamline_context.py create mode 100755 pxii_bec/scripts/init_beamline.py diff --git a/pxii_bec/bec_ipython_client/startup/post_startup.py b/pxii_bec/bec_ipython_client/startup/post_startup.py index eae64ae..e13fb8e 100644 --- a/pxii_bec/bec_ipython_client/startup/post_startup.py +++ b/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() diff --git a/pxii_bec/bec_ipython_client/startup/pre_startup.py b/pxii_bec/bec_ipython_client/startup/pre_startup.py index ea42d71..b835d3a 100644 --- a/pxii_bec/bec_ipython_client/startup/pre_startup.py +++ b/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): """ diff --git a/pxii_bec/device_configs/pxii-devices.yaml b/pxii_bec/device_configs/pxii-devices.yaml index 3358c3a..8a7a824 100644 --- a/pxii_bec/device_configs/pxii-devices.yaml +++ b/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 diff --git a/pxii_bec/macros/APscripts.py b/pxii_bec/macros/APscripts.py index 94a7d59..770b1c1 100644 --- a/pxii_bec/macros/APscripts.py +++ b/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}') diff --git a/pxii_bec/macros/Undu_Helpers.py b/pxii_bec/macros/Undu_Helpers.py index baf31a1..0302550 100755 --- a/pxii_bec/macros/Undu_Helpers.py +++ b/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") diff --git a/pxii_bec/macros/bcu_bpm_positions.py b/pxii_bec/macros/bcu_bpm_positions.py new file mode 100644 index 0000000..165085e --- /dev/null +++ b/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) diff --git a/pxii_bec/macros/beamline_planner.py b/pxii_bec/macros/beamline_planner.py index e96f01d..b62264a 100755 --- a/pxii_bec/macros/beamline_planner.py +++ b/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): diff --git a/pxii_bec/macros/beampos_with_en.py~ b/pxii_bec/macros/beampos_with_en.py~ deleted file mode 100644 index 6d89248..0000000 --- a/pxii_bec/macros/beampos_with_en.py~ +++ /dev/null @@ -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) - - diff --git a/pxii_bec/macros/build_devices.py b/pxii_bec/macros/build_devices.py index f731107..bba4f2a 100755 --- a/pxii_bec/macros/build_devices.py +++ b/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) \ No newline at end of file + return round(data[signal_name]["value"], 8) diff --git a/pxii_bec/macros/calibrate_bpm.py b/pxii_bec/macros/calibrate_bpm.py index ff6c2dc..c0ca7c3 100644 --- a/pxii_bec/macros/calibrate_bpm.py +++ b/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 + + diff --git a/pxii_bec/macros/calibrate_bpm.py.bak b/pxii_bec/macros/calibrate_bpm.py.bak deleted file mode 100644 index 7ad0ad5..0000000 --- a/pxii_bec/macros/calibrate_bpm.py.bak +++ /dev/null @@ -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) - - - diff --git a/pxii_bec/macros/dependencies.py b/pxii_bec/macros/dependencies.py index 6e9e522..2ff2fd2 100755 --- a/pxii_bec/macros/dependencies.py +++ b/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")], diff --git a/pxii_bec/macros/devices_manager.py b/pxii_bec/macros/devices_manager.py index cc201db..720caff 100755 --- a/pxii_bec/macros/devices_manager.py +++ b/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 diff --git a/pxii_bec/macros/flux_macros.py b/pxii_bec/macros/flux_macros.py new file mode 100644 index 0000000..11984e5 --- /dev/null +++ b/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) diff --git a/pxii_bec/macros/init_beamline.py b/pxii_bec/macros/init_beamline.py.bak similarity index 100% rename from pxii_bec/macros/init_beamline.py rename to pxii_bec/macros/init_beamline.py.bak diff --git a/pxii_bec/macros/katscripts.py b/pxii_bec/macros/katscripts.py index 451fb84..513080a 100755 --- a/pxii_bec/macros/katscripts.py +++ b/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) + + + + + + + diff --git a/pxii_bec/macros/mx_basics.py b/pxii_bec/macros/mx_basics.py.bak similarity index 61% rename from pxii_bec/macros/mx_basics.py rename to pxii_bec/macros/mx_basics.py.bak index c089c45..22158ef 100755 --- a/pxii_bec/macros/mx_basics.py +++ b/pxii_bec/macros/mx_basics.py.bak @@ -193,22 +193,131 @@ def fit(data, fit_params): 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) +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) + + # 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": lmfit_result.params["peak_fwhm"].value, - "centre": lmfit_result.best_values["peak_center"], - "height": lmfit_result.params["peak_height"].value, + "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.""" + + 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 plot_fitted_data(data, fit_result): - """Plot the original data and the fitted model.""" + """Plot the original data and display fitting parameters.""" + 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"]) @@ -221,14 +330,21 @@ def plot_fitted_data(data, fit_result): 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 + dock_area = open_docks[dock_area_name] + wf = dock_area.Plot + text_box = dock_area.Results + return wf, text_box @@ -260,12 +376,28 @@ def plot_live_data_bec( wf.y_label = signal_name wf.plot(device_x=motor_name, device_y=signal_name) +def plot_live_data( + wf, + motor_name, + signal_name, + title="Live scan", +): + 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_bec( data, fit_result, + window_name="Fitting", + plot_widget=None ): """ Plot fitted data and display fitting parameters in the specified window. @@ -285,15 +417,22 @@ def plot_fitted_data_bec( 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() + if plot_widget is not None: + wf = plot_widget + else: + wf, text_box = select_bec_window(window_name) + 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}" + f"Chi sq = {fit_result['chi_sq']:.3g}\n" + f"SNR = {fit_result['snr']:.3g}" ) - text_box.set_plain_text(fit_text) + if plot_widget is None: + text_box.set_plain_text(fit_text) + wf.clear_all() wf.title = f"Scan: {data['scan_number']}" wf.x_label = data["motor_name"] @@ -303,3 +442,75 @@ def plot_fitted_data_bec( # wf.Fit.set(symbol_size = 0) wf.get_curve('Fit').set(symbol_size=0) +# def plot_feature_bec( +# data, +# result, +# window_name="Feature", +# ): +# """ +# Plot feature data and indicate the detected edges and centre. +# """ + +# wf, text_box = select_bec_window(window_name) + +# feature_text = ( +# f"Feature centre = {result['centre']:.5g}\n" +# f"Left edge = {result['left_edge']:.5g}\n" +# f"Right edge = {result['right_edge']:.5g}\n" +# f"Width = {result['width']:.5g}\n" +# f"Threshold = {result['threshold']:.4g}" +# ) + +# text_box.set_plain_text(feature_text) + +# wf.clear_all() +# # wf.title = f"Scan: {data['scan_number']}" +# wf.title = "What now??" +# wf.x_label = data["motor_name"] +# wf.y_label = data["signal_name"] + +# x = data["x_data"] +# y = data["y_data"] + +# wf.plot( +# x=x, +# y=y, +# label="Data", +# ) + +# ymin = np.min(y) +# ymax = np.max(y) + +# # Horizontal threshold +# wf.plot( +# x=[np.min(x), np.max(x)], +# y=[result["threshold"], result["threshold"]], +# label="Threshold", +# ) + +# # Left edge +# wf.plot( +# x=[result["left_edge"], result["left_edge"]], +# y=[ymin, ymax], +# label="Left edge", +# ) + +# # Right edge +# wf.plot( +# x=[result["right_edge"], result["right_edge"]], +# y=[ymin, ymax], +# label="Right edge", +# ) + +# # Centre +# wf.plot( +# x=[result["centre"], result["centre"]], +# y=[ymin, ymax], +# label="Centre", +# ) + +# # Remove symbols from guide lines +# wf.get_curve("Threshold").set(symbol_size=0) +# wf.get_curve("Left edge").set(symbol_size=0) +# wf.get_curve("Right edge").set(symbol_size=0) +# wf.get_curve("Centre").set(symbol_size=0) \ No newline at end of file diff --git a/pxii_bec/macros/mx_methods.py b/pxii_bec/macros/mx_methods.py.bak similarity index 57% rename from pxii_bec/macros/mx_methods.py rename to pxii_bec/macros/mx_methods.py.bak index 6a00c64..0775450 100755 --- a/pxii_bec/macros/mx_methods.py +++ b/pxii_bec/macros/mx_methods.py.bak @@ -105,7 +105,10 @@ def go_to_peak( settle: float = FitDefaults.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,7 +134,7 @@ 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) @@ -141,7 +144,15 @@ def go_to_peak( print(f"Signal name is {signal_name}") # wf.plot(x_name=motor_name, y_name=signal_name) if plot: - plot_live_data_bec(motor_name, signal_name) + 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) # Validate and calculate step size step_size = calculate_step_size(start, stop, steps) @@ -192,9 +203,17 @@ def go_to_peak( 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}") + return + # Plot the fitted data if plot = True if plot: - plot_fitted_data_bec(data, fit_result) + plot_fitted_data_bec(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: @@ -202,8 +221,22 @@ def go_to_peak( print(f"Max position is at {value}") move_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_position(data["motor_device"], data["motor_name"], fit_result["centre"], data) + else: + # Safely move the motor to the peak position + move_to_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( @@ -335,3 +368,299 @@ def optimise_kb(mirror): 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) + + +def find_feature( + motor_device, + signal_device, + start: float, + stop: float, + steps: int, + feature: str = "peak", + offset=None, + fraction: float = 0.5, + relative: bool = FitDefaults.RELATIVE_MODE, + plot: bool = True, + settle: float = FitDefaults.SETTLE_TIME, + confirm: bool = True, + window_name: str = "Fitting", + plot_widget=None, +): + """ + Scan a motor, determine the centre of a high- or low-signal region, + and move the motor to that position. + """ + + # 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, + ) + + # Retrieve data + motor_data = ( + scan_result.scan.live_data[mres["live_root"]][mres["name"]].val + ) + + signal_data = ( + scan_result.scan.live_data[sres["root"]][sres["name"]].val + ) + + 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", + } + + # 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"], + feature=feature, + 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}") + + # Move safely to centre + move_to_position( + data["motor_device"], + data["motor_name"], + result["centre"], + data, + ) + + 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 + + plot_widget.plot( + x=[result["centre"], result["centre"]], + y=[ + np.min(data["y_data"]), + np.max(data["y_data"]), + ], + label="Centre", + ) + + + plot_widget.title = f"{motor_name} scan" + plot_widget.x_label = motor_name + plot_widget.y_label = signal_name + + return result + +def find_edge_offset(x, y, fraction=0.5, feature="edge_offset", offset=0.0): + """ + Find a single signal edge and calculate a target position + at a known offset from that edge. + """ + + 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) + diff --git a/pxii_bec/macros/mxcomm_analysis.py b/pxii_bec/macros/mxcomm_analysis.py new file mode 100644 index 0000000..7662d71 --- /dev/null +++ b/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) diff --git a/pxii_bec/macros/mxcomm_data.py b/pxii_bec/macros/mxcomm_data.py new file mode 100644 index 0000000..73c1bb3 --- /dev/null +++ b/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 diff --git a/pxii_bec/macros/mxcomm_plot.py b/pxii_bec/macros/mxcomm_plot.py new file mode 100644 index 0000000..ab9b7b5 --- /dev/null +++ b/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", + ) diff --git a/pxii_bec/macros/mxcomm_scan.py b/pxii_bec/macros/mxcomm_scan.py new file mode 100755 index 0000000..92c49b3 --- /dev/null +++ b/pxii_bec/macros/mxcomm_scan.py @@ -0,0 +1,439 @@ +"""Custom scan routines: +1) a go_to_peak scan, that scans a motor, finds the peak position and moves to peak +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 mxcomm_analysis import ( +# create_fit_parameters, +# get_data_from_history, +# fit, +# 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 +def calculate_step_size(start: float, stop: float, steps: int) -> float: + """ + Provides the function to calculate the step size for dividing a specified range + into a given number of steps. + + Args: + start: The starting value of the range. + stop: The stopping value of the range. + steps: The number of steps to divide the range into. Must be at least 1. + + Raises: + ValueError: If the steps value is less than 1. + + Returns: + The calculated step size as a float, rounded to three decimal places. + """ + if steps < 1: + raise ValueError("Number of steps must be at least 1.") + return round((stop - start) / (steps - 1), 5) + + +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. + + The function verifies if the requested position is within the scan range of the + motor device provided. If the position is outside the range, the motor is + moved to the center of its scan range, an error message is raised, and the + operation is halted. If the position is valid, the motor is moved to the + specified position. + + Parameters: + motor_device: The motor device to be moved. + motor_name: str + The name of the motor as a string + position: float + The desired position to move the motor to. Position should be within + the scan range of the motor determined by the provided data. + data: dict + A dictionary containing "x_data", which is used to determine the + scan range of the motor. + + Raises: + ValueError: Raised if the specified position is outside the valid scan + range determined by "x_data" in the data dictionary. The motor will + return to the center of its scan range in this case. + """ + + motor_min = np.min(data["x_data"]) + motor_max = np.max(data["x_data"]) + motor_centre = (motor_max + motor_min) / 2 + + if not motor_min <= position <= motor_max: + scans.umv(motor_device, motor_centre, relative=False) + msg = ( + f"Position {position: .2f} is outside the scan range of " + f"{motor_min: .2f} to {motor_max: .2f}. " + f"Returning to centre of scan range {motor_centre: .3f}." + ) + raise ValueError(msg) + motor_position = round(position, 4) + scans.umv(motor_device, motor_position, relative=False) + print(f"\n Moving {motor_name} to position {motor_position: .3f}") + +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, + signal_device, + start: float, + stop: float, + steps: int, + relative: bool = ScanDefaults.RELATIVE_MODE, + plot: bool = True, + settle: float = ScanDefaults.SETTLE_TIME, + confirm: bool = True, + gomax: 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 + identifying the optimal position based on signal peak data. + + Parameters: + motor_device: The motor device to be scanned. + signal_device: The signal device to monitor during the scan. + start (float): The starting position of the scan. Ignored if `relative` is True. + stop (float): The ending position of the scan. Ignored if `relative` is True. + steps (int): The number of steps to divide the scan range into. + relative (bool, optional): If True, interpret `start` and `stop` as relative to + the current motor position. Defaults to RELATIVE_MODE constant. + plot (bool, optional): If True, plot the scan data and the fitted results. + Defaults to True. + settle (float, optional): The time in seconds to wait after each step for the + signal to stabilize. Defaults to DEFAULT_SETTLE_TIME constant. + confirm (bool, optional): If True, ask for user confirmation before starting + the scan. Defaults to True. + + Raises: + Exception: Raises exceptions potentially raised by dependent functions or + operations such as plotting, fitting, or motor movement. + + Returns: + fit parameters + """ + # 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}") + + # 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( + 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 + starting_position = mres["value"] + print(f"Starting motor position = {starting_position}") + + if confirm: + if relative: + 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}" + ) + print(f"Relative mode = {relative}") + else: + print( + f"\nScanning from {start: .5g} to {stop: .5g} in {steps} steps of size {step_size}" + ) + 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[motor_name][motor_name].val + # sres = get_value(signal_name, True) + signal_data = scan_result.scan.live_data[sres["root"]][sres["name"]].val + scan_number = "Current" + + 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": scan_number, + } + + # Define and fit model to scan data + 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(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_scan_to_position(data["motor_device"], data["motor_name"], fit_result["x_max"], data) + else: + 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( + history_index: int, + signal_name: str, + deriv: bool = False, + negative: bool = False, + smoothing: bool = False, + model: str = "Voigt", + move_to_peak: bool = False, +): + """ + Retrieve and analyze historical data by fitting a model, optionally moving to + a peak position. + + Parameters: + history_index (int): Index of the historical data set to retrieve. + signal_name (str): Name of the signal to fit. + deriv (bool, optional): Whether to include the derivative in the fitting + procedure. Defaults to False. + model (str, optional): Name of the model to use for fitting. Defaults to + DEFAULT_MODEL. + move_to_peak (bool, optional): Whether to move the motor to the peak position + after fitting. Defaults to False. + + Raises: + KeyError: If required keys are not found in the retrieved data dictionary. + ValueError: If the fitting process fails or produces invalid results. + + Returns: + None + """ + # Retrieve historical data + data = get_data_from_history(history_index, signal_name) + + # Define fitting parameters + fit_params = create_fit_parameters(deriv, negative, model, FitDefaults.BASELINE) + + # Perform fit and plot the data + fit_result = fit(data, fit_params) + plot_fitted_data(data, fit_result) + + # Optionally move the motor to the peak position + if move_to_peak: + move_to_position(data["motor_device"], data["motor_name"], fit_result["centre"], data) + return fit_result + + +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, +): + """ + Scan a motor, determine the centre of a high- or low-signal region, + and move the motor to that position. + """ + + # 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 + ) + + # Retrieve data + motor_data = scan_result.scan.live_data[mres["live_root"]][mres["name"]].val + + signal_data = scan_result.scan.live_data[sres["root"]][sres["name"]].val + + 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", + } + + # 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}") + + # Move safely to centre + move_to_position(data["motor_device"], data["motor_name"], result["centre"], data) + + 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 + + plot_widget.plot( + x=[result["centre"], result["centre"]], + y=[np.min(data["y_data"]), np.max(data["y_data"])], + label="Centre", + ) + + plot_widget.title = f"{motor_name} scan" + plot_widget.x_label = motor_name + plot_widget.y_label = signal_name + + return result diff --git a/pxii_bec/macros/position_device.py b/pxii_bec/macros/position_device.py index f6f7589..a7664de 100755 --- a/pxii_bec/macros/position_device.py +++ b/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 + \ No newline at end of file diff --git a/pxii_bec/macros/pxii_energy.py b/pxii_bec/macros/pxii_energy.py index 541e6e9..1ef07cd 100755 --- a/pxii_bec/macros/pxii_energy.py +++ b/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 + diff --git a/pxii_bec/macros/pxii_parameters.py b/pxii_bec/macros/pxii_parameters.py index c21b95b..44b6b92 100755 --- a/pxii_bec/macros/pxii_parameters.py +++ b/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, -# } diff --git a/pxii_bec/macros/set_kbox.sh b/pxii_bec/macros/set_kbox.sh index c071ac1..83ae5c4 100755 --- a/pxii_bec/macros/set_kbox.sh +++ b/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'..." diff --git a/pxii_bec/macros/test.py b/pxii_bec/macros/test.py deleted file mode 100755 index cc368dc..0000000 --- a/pxii_bec/macros/test.py +++ /dev/null @@ -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)) - diff --git a/pxii_bec/macros/test_state_transitions.py b/pxii_bec/macros/test_state_transitions.py new file mode 100755 index 0000000..64b1944 --- /dev/null +++ b/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") diff --git a/pxii_bec/scripts/__init__.py b/pxii_bec/scripts/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/pxii_bec/scripts/beamline_context.py b/pxii_bec/scripts/beamline_context.py new file mode 100644 index 0000000..3e0de47 --- /dev/null +++ b/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() diff --git a/pxii_bec/scripts/init_beamline.py b/pxii_bec/scripts/init_beamline.py new file mode 100755 index 0000000..a698bae --- /dev/null +++ b/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 -- 2.54.0 From b77ba7b73e21d0f36d8886312c7e357eeafb6845 Mon Sep 17 00:00:00 2001 From: perl_d Date: Tue, 25 Aug 2026 17:55:13 +0200 Subject: [PATCH 2/2] chore: cleanup junk and backup files --- pxii_bec/macros/.nfs000000000001101900000001 | Bin 12288 -> 0 bytes pxii_bec/macros/init_beamline.py.bak | 144 ---- pxii_bec/macros/mx_basics.bak | 313 --------- pxii_bec/macros/mx_basics.py.bak | 516 -------------- pxii_bec/macros/mx_methods.py.bak | 666 ------------------- 5 files changed, 1639 deletions(-) delete mode 100644 pxii_bec/macros/.nfs000000000001101900000001 delete mode 100755 pxii_bec/macros/init_beamline.py.bak delete mode 100755 pxii_bec/macros/mx_basics.bak delete mode 100755 pxii_bec/macros/mx_basics.py.bak delete mode 100755 pxii_bec/macros/mx_methods.py.bak diff --git a/pxii_bec/macros/.nfs000000000001101900000001 b/pxii_bec/macros/.nfs000000000001101900000001 deleted file mode 100644 index 1b30d5e970bdf8e3c8237388300e2274d6b60906..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 12288 zcmeI&F;Buk7zW_O9YvkoDlU*fdo3VH!fYHIF(wYAvC0p)nDYKgNv*2 zm-rRRVn8B^gXX>D&B5{UrT3Jj;kx|nUD;k|DE*NVamy8+ zvQ=@~q3h9MYZ3J)1=Ff>K3iPUwz5{+8eh@YTC1&fSm(vWV7iSZS8T4+`Y5^JRh{`d ho$D(otRD!^Zlc!GY8rDzPn-V+()06$^zIj%>=W7e@wos1 diff --git a/pxii_bec/macros/init_beamline.py.bak b/pxii_bec/macros/init_beamline.py.bak deleted file mode 100755 index edf1a33..0000000 --- a/pxii_bec/macros/init_beamline.py.bak +++ /dev/null @@ -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 diff --git a/pxii_bec/macros/mx_basics.bak b/pxii_bec/macros/mx_basics.bak deleted file mode 100755 index 5bd6f9a..0000000 --- a/pxii_bec/macros/mx_basics.bak +++ /dev/null @@ -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) - diff --git a/pxii_bec/macros/mx_basics.py.bak b/pxii_bec/macros/mx_basics.py.bak deleted file mode 100755 index 22158ef..0000000 --- a/pxii_bec/macros/mx_basics.py.bak +++ /dev/null @@ -1,516 +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) - -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) - - # 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.""" - - 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 plot_fitted_data(data, fit_result): - """Plot the original data and display fitting parameters.""" - - 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: - dock_area = open_docks[dock_area_name] - wf = dock_area.Plot - text_box = dock_area.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_live_data( - wf, - motor_name, - signal_name, - title="Live scan", -): - 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_bec( - data, - fit_result, - window_name="Fitting", - plot_widget=None -): - """ - 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. - """ - if plot_widget is not None: - wf = plot_widget - else: - wf, text_box = select_bec_window(window_name) - - 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}\n" - f"SNR = {fit_result['snr']:.3g}" - ) - if plot_widget is None: - 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) - -# def plot_feature_bec( -# data, -# result, -# window_name="Feature", -# ): -# """ -# Plot feature data and indicate the detected edges and centre. -# """ - -# wf, text_box = select_bec_window(window_name) - -# feature_text = ( -# f"Feature centre = {result['centre']:.5g}\n" -# f"Left edge = {result['left_edge']:.5g}\n" -# f"Right edge = {result['right_edge']:.5g}\n" -# f"Width = {result['width']:.5g}\n" -# f"Threshold = {result['threshold']:.4g}" -# ) - -# text_box.set_plain_text(feature_text) - -# wf.clear_all() -# # wf.title = f"Scan: {data['scan_number']}" -# wf.title = "What now??" -# wf.x_label = data["motor_name"] -# wf.y_label = data["signal_name"] - -# x = data["x_data"] -# y = data["y_data"] - -# wf.plot( -# x=x, -# y=y, -# label="Data", -# ) - -# ymin = np.min(y) -# ymax = np.max(y) - -# # Horizontal threshold -# wf.plot( -# x=[np.min(x), np.max(x)], -# y=[result["threshold"], result["threshold"]], -# label="Threshold", -# ) - -# # Left edge -# wf.plot( -# x=[result["left_edge"], result["left_edge"]], -# y=[ymin, ymax], -# label="Left edge", -# ) - -# # Right edge -# wf.plot( -# x=[result["right_edge"], result["right_edge"]], -# y=[ymin, ymax], -# label="Right edge", -# ) - -# # Centre -# wf.plot( -# x=[result["centre"], result["centre"]], -# y=[ymin, ymax], -# label="Centre", -# ) - -# # Remove symbols from guide lines -# wf.get_curve("Threshold").set(symbol_size=0) -# wf.get_curve("Left edge").set(symbol_size=0) -# wf.get_curve("Right edge").set(symbol_size=0) -# wf.get_curve("Centre").set(symbol_size=0) \ No newline at end of file diff --git a/pxii_bec/macros/mx_methods.py.bak b/pxii_bec/macros/mx_methods.py.bak deleted file mode 100755 index 0775450..0000000 --- a/pxii_bec/macros/mx_methods.py.bak +++ /dev/null @@ -1,666 +0,0 @@ -"""Use the methods in mx_basics to perform: -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 -""" - -from dataclasses import dataclass -import numpy as np - -# from pxiii_parameters import FitDefaults, BPMScans, MirrorConfig - -# from mx_basics import ( -# create_fit_parameters, -# get_data_from_history, -# fit, -# plot_fitted_data_bec, -# plot_live_data_bec, -# ) - - -# Method functions -def calculate_step_size(start: float, stop: float, steps: int) -> float: - """ - Provides the function to calculate the step size for dividing a specified range - into a given number of steps. - - Args: - start: The starting value of the range. - stop: The stopping value of the range. - steps: The number of steps to divide the range into. Must be at least 1. - - Raises: - ValueError: If the steps value is less than 1. - - Returns: - The calculated step size as a float, rounded to three decimal places. - """ - if steps < 1: - raise ValueError("Number of steps must be at least 1.") - return round((stop - start) / (steps - 1), 5) - - -def move_to_position(motor_device, motor_name: str, position: float, data: dict): - """ - Function to move a specified motor device to a given position. - - The function verifies if the requested position is within the scan range of the - motor device provided. If the position is outside the range, the motor is - moved to the center of its scan range, an error message is raised, and the - operation is halted. If the position is valid, the motor is moved to the - specified position. - - Parameters: - motor_device: The motor device to be moved. - motor_name: str - The name of the motor as a string - position: float - The desired position to move the motor to. Position should be within - the scan range of the motor determined by the provided data. - data: dict - A dictionary containing "x_data", which is used to determine the - scan range of the motor. - - Raises: - ValueError: Raised if the specified position is outside the valid scan - range determined by "x_data" in the data dictionary. The motor will - return to the center of its scan range in this case. - """ - - motor_min = np.min(data["x_data"]) - motor_max = np.max(data["x_data"]) - motor_centre = (motor_max + motor_min) / 2 - - if not motor_min <= position <= motor_max: - scans.umv(motor_device, motor_centre, relative=False) - msg = ( - f"Position {position: .2f} is outside the scan range of " - f"{motor_min: .2f} to {motor_max: .2f}. " - f"Returning to centre of scan range {motor_centre: .3f}." - ) - raise ValueError(msg) - motor_position = round(position, 4) - 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 go_to_peak( - motor_device, - signal_device, - start: float, - stop: float, - steps: int, - relative: bool = FitDefaults.RELATIVE_MODE, - plot: bool = True, - settle: float = FitDefaults.SETTLE_TIME, - confirm: bool = True, - gomax: 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 - identifying the optimal position based on signal peak data. - - Parameters: - motor_device: The motor device to be scanned. - signal_device: The signal device to monitor during the scan. - start (float): The starting position of the scan. Ignored if `relative` is True. - stop (float): The ending position of the scan. Ignored if `relative` is True. - steps (int): The number of steps to divide the scan range into. - relative (bool, optional): If True, interpret `start` and `stop` as relative to - the current motor position. Defaults to RELATIVE_MODE constant. - plot (bool, optional): If True, plot the scan data and the fitted results. - Defaults to True. - settle (float, optional): The time in seconds to wait after each step for the - signal to stabilize. Defaults to DEFAULT_SETTLE_TIME constant. - confirm (bool, optional): If True, ask for user confirmation before starting - the scan. Defaults to True. - - Raises: - Exception: Raises exceptions potentially raised by dependent functions or - operations such as plotting, fitting, or motor movement. - - Returns: - fit parameters - """ - # 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}") - # 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) - - # 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}") - - if confirm: - if relative: - scan_start = current_motor_position + start - scan_end = current_motor_position + stop - print( - f"\nScanning from {scan_start: .6g} to {scan_end: .6g} in " - f"{steps} steps of size {step_size}" - ) - print(f"Relative mode = {relative}") - else: - print( - f"\nScanning from {start: .5g} to {stop: .5g} in {steps} steps of size {step_size}" - ) - 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[motor_name][motor_name].val - # sres = get_value(signal_name, True) - signal_data = scan_result.scan.live_data[sres['root']][sres['name']].val - scan_number = "Current" - - 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": scan_number, - } - - # Define and fit model to scan data - 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}") - return - - # Plot the fitted data if plot = True - if plot: - plot_fitted_data_bec(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) - else: - if gap: - if fit_result["centre"] < 4.5: - umv(dev.id_gap, 4.5) - else: - move_to_position(data["motor_device"], data["motor_name"], fit_result["centre"], data) - else: - # Safely move the motor to the peak position - move_to_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( - history_index: int, - signal_name: str, - deriv: bool = False, - negative: bool = False, - smoothing: bool = False, - model: str = FitDefaults.MODEL, - move_to_peak: bool = False, -): - """ - Retrieve and analyze historical data by fitting a model, optionally moving to - a peak position. - - Parameters: - history_index (int): Index of the historical data set to retrieve. - signal_name (str): Name of the signal to fit. - deriv (bool, optional): Whether to include the derivative in the fitting - procedure. Defaults to False. - model (str, optional): Name of the model to use for fitting. Defaults to - DEFAULT_MODEL. - move_to_peak (bool, optional): Whether to move the motor to the peak position - after fitting. Defaults to False. - - Raises: - KeyError: If required keys are not found in the retrieved data dictionary. - ValueError: If the fitting process fails or produces invalid results. - - Returns: - None - """ - # Retrieve historical data - data = get_data_from_history(history_index, signal_name) - - # Define fitting parameters - fit_params = create_fit_parameters(deriv, negative, model, FitDefaults.BASELINE) - - # Perform fit and plot the data - fit_result = fit(data, fit_params) - plot_fitted_data_bec(data, fit_result) - - - # Optionally move the motor to the peak position - if move_to_peak: - move_to_position(data["motor_device"], data["motor_name"], fit_result["centre"], data) - return fit_result - - -def scan_bpm(bpmname): - """ - 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" - - """ - - # 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 - ) - 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) - - 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") - - 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") - - 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") - - 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") - - 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) - - -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. - - Parameters: - mirror: either "hfm" or :vfm" - - """ - - # 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) - - -def find_feature( - motor_device, - signal_device, - start: float, - stop: float, - steps: int, - feature: str = "peak", - offset=None, - fraction: float = 0.5, - relative: bool = FitDefaults.RELATIVE_MODE, - plot: bool = True, - settle: float = FitDefaults.SETTLE_TIME, - confirm: bool = True, - window_name: str = "Fitting", - plot_widget=None, -): - """ - Scan a motor, determine the centre of a high- or low-signal region, - and move the motor to that position. - """ - - # 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, - ) - - # Retrieve data - motor_data = ( - scan_result.scan.live_data[mres["live_root"]][mres["name"]].val - ) - - signal_data = ( - scan_result.scan.live_data[sres["root"]][sres["name"]].val - ) - - 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", - } - - # 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"], - feature=feature, - 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}") - - # Move safely to centre - move_to_position( - data["motor_device"], - data["motor_name"], - result["centre"], - data, - ) - - 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 - - plot_widget.plot( - x=[result["centre"], result["centre"]], - y=[ - np.min(data["y_data"]), - np.max(data["y_data"]), - ], - label="Centre", - ) - - - plot_widget.title = f"{motor_name} scan" - plot_widget.x_label = motor_name - plot_widget.y_label = signal_name - - return result - -def find_edge_offset(x, y, fraction=0.5, feature="edge_offset", offset=0.0): - """ - Find a single signal edge and calculate a target position - at a known offset from that edge. - """ - - 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) - -- 2.54.0