feat: add config for beam steering devices

pxiii_bec/device_configs/beam_steerer.yaml

@@ -0,0 +1,28 @@
+beam_profile:
+  readoutPriority: baseline
+  description: Beam location and size based on AD sample camera analysis
+  deviceClass: pxiii_bec.devices.BeamProfile
+  deviceConfig:
+    prefix: "X06DA-ES-MS:XYP1:"
+  deviceTags:
+    - beam
+  enabled: true
+  readOnly: false
+  softwareTrigger: false
+beam_steering:
+  readoutPriority: on_request
+  description: Beam location and size based on AD sample camera analysis
+  deviceClass: pxiii_bec.devices.BeamSteerer
+  deviceConfig:
+    prefix: ""
+    hfm_step_per_x_px: -0.00013
+    vfm_step_per_y_px: 0.00033
+    sample_loc_x_px: 1000
+    sample_loc_y_px: 1000
+    px_tolerance: 2
+    step_fraction: 0.4
+  deviceTags:
+    - beam
+  enabled: true
+  readOnly: false
+  softwareTrigger: false

pxiii_bec/device_configs/x06da_device_config.yaml

@@ -1,5 +1,6 @@
 base_config:
    - !include ./pxiii-devices.yaml
+   - !include ./beam_steerer.yaml
 #states_config:
 #  - !include ./pxiii-state-devices.yaml
   

pxiii_bec/devices/__init__.py

@@ -4,11 +4,15 @@ Ophyd devices for the PX III beamline, including the MX specific Aerotech A3200
 
 @author: mohacsi_i
 """
+
+from mx_bec.devices.beam_profile import BeamProfile
+from mx_bec.devices.beam_steering import BeamSteerer
+
 from .A3200 import AerotechAbrStage
 from .A3200utils import A3200Axis
+from .NDArrayPreview import NDArrayPreview
+from .PneumaticValve import PneumaticValve
+from .SamCamDetector import SamCamDetector
 from .SmarGonA import SmarGonAxis as SmarGonAxisA
 from .SmarGonB import SmarGonAxis as SmarGonAxisB
 from .StdDaqPreview import StdDaqPreviewDetector
-from .NDArrayPreview import NDArrayPreview
-from .SamCamDetector import SamCamDetector
-from .PneumaticValve import PneumaticValve

pxiii_bec/macros/katscripts.py

@@ -1,14 +1,15 @@
 """Some simple scripts for commissioning"""
 
-import time, os
-from datetime import datetime, timedelta
+import os
+import time
 from dataclasses import dataclass
-import numpy as np
+from datetime import datetime, timedelta
+
 import matplotlib.pyplot as plt
-from scipy.signal import find_peaks as fp
+import numpy as np
 import pandas as pd
 from scipy.optimize import linear_sum_assignment
-
+from scipy.signal import find_peaks as fp
 
 
 def bragg_scan():
@@ -173,7 +174,7 @@ def generate_dcm_lut(start_energy=6000, end_energy=30000, step=300):
     now = datetime.now()
     fnow = now.strftime("%d%m%H%M")
     filename = f"luts/{fnow}_lut.csv"
-    zoom = dev.scam_zoom.read()['scam_zoom']['value']
+    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")
 
@@ -197,8 +198,8 @@ def generate_dcm_lut(start_energy=6000, end_energy=30000, step=300):
         # 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']
+        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"
@@ -252,7 +253,6 @@ def write_beamsize():
         f.write(f"{lens_config},{sx:.4g},{sy:.4g},{i1_reading:.4g}\n")
 
 
-
 def scan_Ge():
     scanbefore = 30  # eV
     scanafter = 30  # eV
@@ -269,9 +269,9 @@ def scan_Ge():
     filedir, fname = fnow()
     filename = f"{filedir}{fname}_ge_edge.csv"
     print(filename)
-    with open(filename, 'w') as f:
+    with open(filename, "w") as f:
         f.write(f"Energy,Bragg,bcu_bpmsum\n")
-    
+
     # scanstart = edge_Ge + scanafter
     # scanend = edge_Ge - scanbefore
     # bragg_start = convert_from_energy(scanstart)["bragg_angle_mrad"]
@@ -286,9 +286,9 @@ def scan_Ge():
     # print(s)
     # umv(dev.ss_f4_x, 2)
 
-    plot_live_data_bec(dev.dcm_bragg,dev.bcu_bpmsum)
+    plot_live_data_bec(dev.dcm_bragg, dev.bcu_bpmsum)
     scanpoints = np.linspace(current_bragg - 0.5, current_bragg + 0.5, 300)
-    start_energy = convert_from_bragg(current_bragg+0.5)["energy_ev"]
+    start_energy = convert_from_bragg(current_bragg + 0.5)["energy_ev"]
     bl_energy(start_energy)
     umv(dev.ss_f4_x, 6)  # Put Ge filter in place
 
@@ -298,7 +298,7 @@ def scan_Ge():
         inten = average_inten()
         xdata.append(energy)
         ydata.append(inten)
-        with open(filename, 'a') as f:
+        with open(filename, "a") as f:
             f.write(f"{energy:.2f},{bragg:.2f},{inten:.2f}\n")
     data = {
         "x_data": np.array(xdata),
@@ -310,23 +310,20 @@ def scan_Ge():
     }
 
     # Define and fit model to scan data
-    fit_params = create_fit_parameters(deriv = True,
-                                       negative = False, 
-                                       model = FitDefaults.MODEL, 
-                                       baseline = FitDefaults.BASELINE)
+    fit_params = create_fit_parameters(
+        deriv=True, negative=False, model=FitDefaults.MODEL, baseline=FitDefaults.BASELINE
+    )
     fit_result = fit(data, fit_params)
 
     # Plot the fitted data if plot = True
- 
+
     plot_fitted_data_bec(data, fit_result)
 
 
-
-
 def cal_samcam_fixed_x():
     # lat = dev.hfm_x.read()["hfm_x"]["value"]
     # pitch = dev.hfm_pitch.read()["hfm_pitch"]["value"]
-    fpitch = dev.dcm_fpitch.read()['dcm_fpitch']['value']
+    fpitch = dev.dcm_fpitch.read()["dcm_fpitch"]["value"]
     scanpoints = np.linspace(3.75, 4.75, 7)
     zoom_high = 800
     zoom_low = 200
@@ -690,9 +687,10 @@ def move_to_x(targetx=0.13):
 #     umv(dev.vfm_y, mot_new)
 
 
-from epics import caget, caput
 import time
 
+from epics import caget, caput
+
 
 def auto_exposure(
     cam="samcam", target=200, tol=5, gain=0.9, min_exp=0.00001, max_exp=0.2, max_iter=10
@@ -802,12 +800,14 @@ 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 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()
@@ -842,7 +842,6 @@ def energy_vs_height(start_energy=11000, end_energy=25000, step=1000):
     # umv(dev.id_gap, 20.0)
 
 
-
 def scan_in_x(xsize, xcen, mon, small, large, width, steps):
     umv(xsize, small)
     go_to_peak(xcen, mon, -(width * 0.5), (width * 0.5), steps=steps, relative=False)
@@ -860,16 +859,7 @@ def trans_fils(inc):
     x2 = dev.ss_f2_x.read()["ss_f2_x"]["value"]
     x3 = dev.ss_f3_x.read()["ss_f3_x"]["value"]
     x4 = dev.ss_f4_x.read()["ss_f4_x"]["value"]
-    umv(
-        dev.ss_f1_x,
-        x1 + inc,
-        dev.ss_f2_x,
-        x2 + inc,
-        dev.ss_f3_x,
-        x3 + inc,
-        dev.ss_f4_x,
-        x4 + inc,
-    )
+    umv(dev.ss_f1_x, x1 + inc, dev.ss_f2_x, x2 + inc, dev.ss_f3_x, x3 + inc, dev.ss_f4_x, x4 + inc)
 
 
 def calc_from_centre(cam="samcam"):
@@ -915,17 +905,19 @@ def get_cam(low_zoom, high_zoom):
 #     yn = ((readings["bpm1"] + readings["bpm3"]) - (readings["bpm2"] + readings["bpm4"])) / total
 #     return xn, yn
 
-def average_inten():
-            vals_inten = []
-            for _ in range(5):
-                vals_inten.append(dev.bcu_bpmsum.read()['bcu_bpmsum']['value'])
-                time.sleep(0.1)
-            inten = np.mean(vals_inten)
-            return inten
 
-def mirror_pitch(mirror = 'vfm'):
-    
-    if mirror == 'vfm':
+def average_inten():
+    vals_inten = []
+    for _ in range(5):
+        vals_inten.append(dev.bcu_bpmsum.read()["bcu_bpmsum"]["value"])
+        time.sleep(0.1)
+    inten = np.mean(vals_inten)
+    return inten
+
+
+def mirror_pitch(mirror="vfm"):
+
+    if mirror == "vfm":
         pitch = dev.vfm_pitch.position
         mot = dev.vfm_pitch
         start = 2.94
@@ -934,9 +926,9 @@ def mirror_pitch(mirror = 'vfm'):
         mot = dev.hfm_pitch
         start = 3.3
     print(f"Scanning {mirror}_pitch from 10-29 keV")
-    pitch_range = [start - 0.25, start - 0.2, start -0.15, start - 0.1, start -0.05]
+    pitch_range = [start - 0.25, start - 0.2, start - 0.15, start - 0.1, start - 0.05]
 
-    umv(mot,start)
+    umv(mot, start)
 
     results = []
     now = datetime.now()
@@ -958,65 +950,58 @@ def mirror_pitch(mirror = 'vfm'):
             vals_diode = []
 
             for _ in range(5):
-                vals_diode.append(dev.xidiode.read()['xidiode']['value'])
+                vals_diode.append(dev.xidiode.read()["xidiode"]["value"])
                 time.sleep(0.1)
 
             diode = np.mean(vals_diode)
 
-            results.append([
-                energy,
-                p,
-                diode,
-            ])
+            results.append([energy, p, diode])
 
             with open(csvfile, "a") as f:
-                f.write(
-                    f"{energy:.2f},{p:.5f},{diode:.5g}\n"
-                )
-    umv(mot,start)
+                f.write(f"{energy:.2f},{p:.5f},{diode:.5g}\n")
+    umv(mot, start)
     bl_energy(12400)
     results = np.array(results)
 
-def find_best_roll(low = 6000, high = 20000):
+
+def find_best_roll(low=6000, high=20000):
     low_scanpoints = np.linspace(2.0, 7.0, 11)
-    high_scanpoints = np.linspace(1.0, 9.0,9)
+    high_scanpoints = np.linspace(1.0, 9.0, 9)
     highs_y, lows_y = [], []
 
     # high energy
     bl_energy(high)
     filename = f"./luts/roll_{high:.0f}.csv"
     with open(filename, "w", encoding="utf-8") as f:
-        f.write("Roll, Y pos\n") 
+        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 = dev.samcam_y.read()["samcam_y"]["value"]
         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") 
+            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") 
+        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']
+        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") 
+            f.write(f"{i:.4g},{low_y:.4g} \n")
 
-    
     y_high = np.array(highs_y)
     y_low = np.array(lows_y)
 
-
     # # Fit linear models: y = m*x + b
     m_high, b_high = np.polyfit(high_scanpoints, y_high, 1)
     m_low, b_low = np.polyfit(low_scanpoints, y_low, 1)
@@ -1030,24 +1015,25 @@ def find_best_roll(low = 6000, high = 20000):
     print(f"  samcam Y = {y_intersect:.2f}")
 
     # # --- Plot ---
-    plt.figure(figsize=(7,5))
-    plt.scatter(high_scanpoints, y_high, color='red', label=f'High energy')
-    plt.scatter(low_scanpoints, y_low, color='blue', label=f'Low energy')
+    plt.figure(figsize=(7, 5))
+    plt.scatter(high_scanpoints, y_high, color="red", label=f"High energy")
+    plt.scatter(low_scanpoints, y_low, color="blue", label=f"Low energy")
 
     x_fit = np.linspace(min(high_scanpoints), max(high_scanpoints), 100)
-    plt.plot(x_fit, m_high*x_fit + b_high, 'r--', label=f'High energy')
-    plt.plot(x_fit, m_low*x_fit + b_low, 'b--', label=f'Low energy')
+    plt.plot(x_fit, m_high * x_fit + b_high, "r--", label=f"High energy")
+    plt.plot(x_fit, m_low * x_fit + b_low, "b--", label=f"Low energy")
 
-    plt.scatter(x_intersect, y_intersect, color='green', s=80, zorder=5, label='Intersection')
+    plt.scatter(x_intersect, y_intersect, color="green", s=80, zorder=5, label="Intersection")
 
-    plt.xlabel('DCM fine roll')
-    plt.ylabel('samcam Y pos')
-    plt.title('Roll Calibration')
+    plt.xlabel("DCM fine roll")
+    plt.ylabel("samcam Y pos")
+    plt.title("Roll Calibration")
     plt.legend()
     plt.grid(True)
     plt.tight_layout()
     plt.show()
 
+
 def vfm_yaw_check():
     slit_size = 0.02
     slit_point1 = -1.2
@@ -1066,18 +1052,18 @@ def vfm_yaw_check():
         # Point 1
         print(f"Moving slits to position 1")
         umv(dev.ss_sl_ycen, slit_point1)
-        scans.line_scan(dev.vfm_x, mirror_start, mirror_end, steps = 20, relative=False)
+        scans.line_scan(dev.vfm_x, mirror_start, mirror_end, steps=20, relative=False)
         fit_result = fit_history(-1, "xidiode", deriv=True)
-        centre = fit_result['centre']
-        print(f'Yaw = {yaw}, position = {slit_point1}, centre = {centre}')
+        centre = fit_result["centre"]
+        print(f"Yaw = {yaw}, position = {slit_point1}, centre = {centre}")
         point1.append(centre)
         # Point 2
         print(f"Moving slits to position 2")
         umv(dev.ss_sl_ycen, slit_point2)
-        scans.line_scan(dev.vfm_x, mirror_start, mirror_end, steps = 20, relative=False)
+        scans.line_scan(dev.vfm_x, mirror_start, mirror_end, steps=20, relative=False)
         fit_result = fit_history(-1, "xidiode", deriv=True)
-        centre = fit_result['centre']
-        print(f'Yaw = {yaw}, position = {slit_point2}, centre = {centre}')
+        centre = fit_result["centre"]
+        print(f"Yaw = {yaw}, position = {slit_point2}, centre = {centre}")
         point2.append(centre)
     print(f"Point 1 values are {point1}")
     print(f"Point 2 values are {point2}")
@@ -1086,7 +1072,6 @@ def vfm_yaw_check():
     y_point1 = np.array(point1)
     y_point2 = np.array(point2)
 
-
     # # Fit linear models: y = m*x + b
     m_point1, b_point1 = np.polyfit(yaw_values, y_point1, 1)
     m_point2, b_point2 = np.polyfit(yaw_values, y_point2, 1)
@@ -1100,34 +1085,35 @@ def vfm_yaw_check():
     print(f"  Y = {y_intersect:.2f}")
 
     # # --- Plot ---
-    plt.figure(figsize=(7,5))
-    plt.scatter(yaw_values, y_point1, color='red', label=f'point1')
-    plt.scatter(yaw_values, y_point2, color='blue', label=f'point2')
+    plt.figure(figsize=(7, 5))
+    plt.scatter(yaw_values, y_point1, color="red", label=f"point1")
+    plt.scatter(yaw_values, y_point2, color="blue", label=f"point2")
 
     x_fit = np.linspace(min(yaw_values), max(yaw_values), 100)
-    plt.plot(x_fit, m_point1*x_fit + b_point1, 'r--', label=f'point1')
-    plt.plot(x_fit, m_point2*x_fit + b_point2, 'b--', label=f'point2')
+    plt.plot(x_fit, m_point1 * x_fit + b_point1, "r--", label=f"point1")
+    plt.plot(x_fit, m_point2 * x_fit + b_point2, "b--", label=f"point2")
 
-    plt.scatter(x_intersect, y_intersect, color='green', s=80, zorder=5, label='Intersection')
+    plt.scatter(x_intersect, y_intersect, color="green", s=80, zorder=5, label="Intersection")
 
-    plt.xlabel('VFM yaw')
-    plt.ylabel('Edge')
-    plt.title('VFM Yaw Calibration')
+    plt.xlabel("VFM yaw")
+    plt.ylabel("Edge")
+    plt.title("VFM Yaw Calibration")
     plt.legend()
     plt.grid(True)
     plt.tight_layout()
     plt.show()
 
-from dataclasses import dataclass
+
 import time
+from dataclasses import dataclass
 
 
 @dataclass
 class SlitAxis:
     name: str
-    centre_motor: object      # e.g. dev.slit1_x
-    size_motor: object        # e.g. dev.slit1_xsize
-    scan_func: callable       # your existing peak scan function
+    centre_motor: object  # e.g. dev.slit1_x
+    size_motor: object  # e.g. dev.slit1_xsize
+    scan_func: callable  # your existing peak scan function
     scan_range: float = 2.0
     open_size: float = 2.0
     min_size: float = -0.4
@@ -1149,7 +1135,6 @@ def read_intensity(detector, n=5, delay=0.1):
     return sum(values) / len(values)
 
 
-
 def reduce_to_transmission(axis, detector, baseline, target=0.15, step=0.05):
     """
     Close slit size until intensity is around target * baseline.
@@ -1161,8 +1146,7 @@ def reduce_to_transmission(axis, detector, baseline, target=0.15, step=0.05):
         transmission = intensity / baseline
 
         print(
-            f"{axis.name}: size={size:.4f}, "
-            f"I={intensity:.3g}, transmission={transmission:.3f}"
+            f"{axis.name}: size={size:.4f}, " f"I={intensity:.3g}, transmission={transmission:.3f}"
         )
 
         if transmission <= target:
@@ -1185,8 +1169,7 @@ def open_to_transmission(axis, detector, baseline, target=0.99, step=0.05):
         transmission = intensity / baseline
 
         print(
-            f"{axis.name}: size={size:.4f}, "
-            f"I={intensity:.3g}, transmission={transmission:.3f}"
+            f"{axis.name}: size={size:.4f}, " f"I={intensity:.3g}, transmission={transmission:.3f}"
         )
 
         if transmission >= target:
@@ -1200,6 +1183,7 @@ def open_to_transmission(axis, detector, baseline, target=0.99, step=0.05):
 
     raise RuntimeError(f"{axis.name}: reached maximum slit size without restoring intensity")
 
+
 # def approx_centre_axis(axis, detector, slitsize=0.5):
 #     umv(axis.size_motor, slitsize)
 #     go_to_peak(axis.centre_motor, detector, -4, 4, 50, confirm=False)
@@ -1218,25 +1202,19 @@ def centre_axis(axis, detector, close_fraction=0.10, open_fraction=0.95):
     # if baseline <= 0:
     #     raise RuntimeError(f"{axis.name}: invalid baseline intensity")
 
-    reduce_to_transmission(
-        axis,
-        detector,
-        baseline,
-        target=close_fraction,
-    )
+    reduce_to_transmission(axis, detector, baseline, target=close_fraction)
 
-    go_to_peak(axis.centre_motor, detector, axis.centre_motor.position - axis.scan_range,
-               axis.centre_motor.position + axis.scan_range,
-               50,
-               confirm=False)
+    go_to_peak(
+        axis.centre_motor,
+        detector,
+        axis.centre_motor.position - axis.scan_range,
+        axis.centre_motor.position + axis.scan_range,
+        50,
+        confirm=False,
+    )
     print(f"{axis.name}: moving centre to peak position")
 
-    open_to_transmission(
-        axis,
-        detector,
-        baseline,
-        target=open_fraction,
-    )
+    open_to_transmission(axis, detector, baseline, target=open_fraction)
 
     final_intensity = read_intensity(detector)
     print(
@@ -1244,6 +1222,7 @@ def centre_axis(axis, detector, close_fraction=0.10, open_fraction=0.95):
         f"transmission = {final_intensity / baseline:.3f}"
     )
 
+
 def open_slits(size=3):
     slits = define_slits()
     for slit in slits:
@@ -1268,12 +1247,14 @@ def centre_all_slits(slits, detector):
             print(f"WARNING: failed to align {slit.name}: {exc}")
             print("Continuing with next slit.")
 
+
 def centre_one(name, detector):
     slits = define_slits()
     for slit in slits:
         if slit.name == name:
             centre_slit(slit, detector)
 
+
 def define_slits():
     slits = [
         Slit(
@@ -1321,9 +1302,10 @@ def define_slits():
                 scan_func=go_to_peak,
             ),
         ),
-    ] 
+    ]
     return slits
 
+
 def get_mirror_data(history_index: int):
     """Read data from the BEC history and return the X and Y data as arrays."""
     motor1 = "vfm_bu"
@@ -1333,14 +1315,13 @@ def get_mirror_data(history_index: int):
     md = scan.metadata["bec"]
     x1_data = scan.devices[motor1][motor1].read()["value"]
     x2_data = scan.devices[motor2][motor2].read()["value"]
-    y1_data = scan.devices['microscope_stats']['microscope_stats_data_sigma_y'].read()['value']
-    y2_data = scan.devices['microscope_stats']['microscope_stats_data_eccentricity'].read()['value']
+    y1_data = scan.devices["microscope_stats"]["microscope_stats_data_sigma_y"].read()["value"]
+    y2_data = scan.devices["microscope_stats"]["microscope_stats_data_eccentricity"].read()["value"]
     np.savetxt("luts/x1_data", x1_data)
     np.savetxt("luts/x2_data", x2_data)
     np.savetxt("luts/y1_data", y1_data)
     np.savetxt("luts/y2_data", y2_data)
 
-    
     # plt.imshow(y1_data, cmap="hot")
     # plt.clf
     # return {
@@ -1364,21 +1345,19 @@ def record_beampos_until(stop_time_str, planner, interval_minutes=10):
     fdir = "/sls/x06da/config/bec/production/pxiii_bec/pxiii_bec/macros/"
     fname = f"luts/beampos/beampos_{fnow}.csv"
     filename = fdir + fname
-    scam_zoom = dev.scam_zoom.read()['scam_zoom']['value']
+    scam_zoom = dev.scam_zoom.read()["scam_zoom"]["value"]
     umv(dev.scam_zoom, 1000)
-    #dev.transm.put(0.1)
+    # dev.transm.put(0.1)
     beam_size = "25 x 20"
 
     energy = get_current_energy()
 
-
     # Write the header to the file
     with open(filename, "w", encoding="utf-8") as f:
         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")
 
-
     # Parse the stop time (HH:MM)
     stop_hour, stop_minute = map(int, stop_time_str.split(":"))
     end_time = now.replace(hour=stop_hour, minute=stop_minute, second=0, microsecond=0)
@@ -1405,8 +1384,6 @@ def record_beampos_until(stop_time_str, planner, interval_minutes=10):
 
     print(f"Finished — reached {end_time:%H:%M}.")
     planner.move_to(BeamlineState.MANUAL_SAMPLE_EXCHANGE)
-    
-
 
 
 def record_beam(filename):
@@ -1414,9 +1391,9 @@ def record_beam(filename):
     Record the beam position to file.
     """
     # Set transmission and bpm position
-    #dev.transm.put(0.1)
-    #umv(dev.bcu_bpm_x, 0.43)
-    #umv(dev.bcu_bpm_y, -1.543 )
+    # dev.transm.put(0.1)
+    # umv(dev.bcu_bpm_x, 0.43)
+    # umv(dev.bcu_bpm_y, -1.543 )
 
     # Open shutter and record positions
     dev.bcu_shutter.put(1)
@@ -1424,6 +1401,13 @@ def record_beam(filename):
     time.sleep(2)
     now = datetime.now()
 
+    # TODO: introduce set threshold based on histogram
+    # PV: Stats5:CentroidThreshold
+    # Get histogram, find index of 97th %ile, use it as threshold
+    # histo -> cumsum -> final value * 0.97 -> index where cumsum exceeds this -> this value should be the threshold
+
+    # Alternatively, use only green channel, fit gaussian
+
     # Append data to the file
     with open(filename, "a", encoding="utf-8") as f:
         fnow = now.strftime("%H:%M:%S")
@@ -1431,11 +1415,11 @@ def record_beam(filename):
         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()
-        #f.write(
+        # x_bpm, y_bpm = compute_norm()
+        # 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.write(f"{fnow},{xpos_cam: .5g},{ypos_cam: .4g}\n")
     # Close shutter
     dev.bcu_shutter.put(0)
@@ -1464,7 +1448,6 @@ def analyse_beampos(file):
     data["Scin_x"] = (data["Scin_x"] - ref_scinx) * calc_scam_microns(1)
     data["Scin_y"] = (data["Scin_y"] - ref_sciny) * calc_scam_microns(1)
 
-
     # Extract Y-position columns (all columns except the first column)
     y_positions = data.iloc[:, 1:]  # Exclude the first column
 
@@ -1474,7 +1457,7 @@ def analyse_beampos(file):
     # Flatten axs for easy iteration (if axs is a 2D array of axes)
     axs = axs.flatten()
     # Specify which columns to color differently
-    #red_columns = ["bpm_x", "scin_x"]
+    # red_columns = ["bpm_x", "scin_x"]
     red_columns = ["scin_x"]
 
     for idx, column in enumerate(y_positions.columns):
@@ -1485,9 +1468,7 @@ def analyse_beampos(file):
 
         # Plot the line and scatter with the chosen color
         ax.plot(meas_time, y_positions[column], color=color, linewidth=0.5)  # Line plot
-        ax.scatter(
-            meas_time, y_positions[column], color=color, label=column
-        )  # Scatter plot
+        ax.scatter(meas_time, y_positions[column], color=color, label=column)  # Scatter plot
 
         # Set x-axis tick positions and labels
         tick_positions = range(0, len(meas_time), 7)  # Adjust the step size as needed
@@ -1504,14 +1485,15 @@ def analyse_beampos(file):
     # Show the resulting figure
     plt.show()
 
+
 import time
-import numpy as np
+
 import matplotlib.pyplot as plt
+import numpy as np
 from epics import caget
 
 
-
-def scan_scin_focus(planner, type = 'fine'):
+def scan_scin_focus(planner, type="fine"):
     dev.bcu_shutter.put(0)
     planner.move_to(BeamlineState.BEAM_VISUALISATION)
     dev.bcu_shutter.put(1)
@@ -1520,10 +1502,9 @@ def scan_scin_focus(planner, type = 'fine'):
     min = 0.1
     max = 9.9
 
-
-    if type == 'fine':
+    if type == "fine":
         steps = 100
-    elif type == 'coarse':
+    elif type == "coarse":
         steps = 30
     else:
         steps = 60