diff --git a/pxii_bec/device_configs/pxii-devices.yaml b/pxii_bec/device_configs/pxii-devices.yaml
index 2d0ab75..782114e 100644
--- a/pxii_bec/device_configs/pxii-devices.yaml
+++ b/pxii_bec/device_configs/pxii-devices.yaml
@@ -55,8 +55,8 @@ aerotech:
   deviceConfig:
     prefix: http://mx-x10sa-queue-01:5234
   deviceTags:
-  - aerotech
   - motors
+  - aerotech
   enabled: true
   onFailure: buffer
   readoutPriority: baseline
@@ -64,14 +64,13 @@ aerotech:
     u:
       mount: 0
       work: 0
-    work: null
     x:
       in: 0
       out: -10
       safe: -50
     y:
       mount: -0.0
-      work: -0.16
+      work: -0.0
     z:
       mount: 0
       work: 0
@@ -292,7 +291,7 @@ bs_x:
   onFailure: buffer
   readoutPriority: baseline
   userParameter:
-    in: 1.9
+    in: 0.2
 bs_y:
   description: Beamstop Y
   deviceClass: ophyd_devices.EpicsMotor
@@ -304,7 +303,7 @@ bs_y:
   onFailure: buffer
   readoutPriority: baseline
   userParameter:
-    in: -0.6
+    in: -3.154
 bs_z:
   description: Beamstop Z
   deviceClass: ophyd_devices.EpicsMotor
@@ -320,7 +319,7 @@ bs_z:
     max_blout: 70
     min: 13
     safe: 45.1
-    samp: 15.3
+    samp: 16.408
     type: continuous
     work_min: 20
 bsf_f1_y:
@@ -433,6 +432,8 @@ coll_x:
   enabled: true
   onFailure: buffer
   readoutPriority: baseline
+  userParameter:
+    in: -3.080
 coll_y:
   description: Collimator Y
   deviceClass: ophyd_devices.EpicsMotor
@@ -444,7 +445,7 @@ coll_y:
   onFailure: buffer
   readoutPriority: baseline
   userParameter:
-    in: 39.393
+    in: 39.16
     intermediate: 32
     out: 20.002
     park: 1
@@ -602,6 +603,9 @@ det_y:
   enabled: true
   onFailure: buffer
   readoutPriority: baseline
+  userParameter:
+    work: 70
+    xi: 214
 det_z:
   description: Detector Z
   deviceClass: ophyd_devices.EpicsMotor
@@ -614,7 +618,7 @@ det_z:
   readoutPriority: baseline
   userParameter:
     mse: 800
-    vis: 1198
+    vis: 800 
 diag_y:
   description: Scintillator/diode Y
   deviceClass: ophyd_devices.EpicsMotor
@@ -1206,6 +1210,26 @@ mag_ref:
   onFailure: buffer
   readOnly: true
   readoutPriority: baseline
+microscope_stats:
+  description: Stats Signal of MS camera
+  deviceClass: pxii_bec.devices.stats_signal.StatsSignal
+  deviceConfig:
+    prefix: 'X10SA-ES-MS:'
+  enabled: true
+  readoutPriority: async
+  softwareTrigger: true
+samcam_ecc:
+  description: Sample Camera Eccentricity
+  deviceClass: ophyd.EpicsSignalRO
+  deviceConfig:
+    auto_monitor: true
+    read_pv: X10SA-ES-MS:Stats5:Eccentricity_RBV
+  deviceTags:
+  - scam
+  enabled: true
+  onFailure: buffer
+  readOnly: true
+  readoutPriority: monitored
 samcam_exp:
   description: Sample Camera Exposure
   deviceClass: ophyd.EpicsSignal
@@ -1240,6 +1264,30 @@ samcam_max:
   onFailure: buffer
   readOnly: true
   readoutPriority: monitored
+samcam_mean:
+  description: Sample Camera Mean
+  deviceClass: ophyd.EpicsSignalRO
+  deviceConfig:
+    auto_monitor: true
+    read_pv: X10SA-ES-MS:Stats1:MeanValue_RBV
+  deviceTags:
+  - scam
+  enabled: true
+  onFailure: buffer
+  readOnly: true
+  readoutPriority: monitored
+samcam_sigma:
+  description: Sample Camera Sigma
+  deviceClass: ophyd.EpicsSignalRO
+  deviceConfig:
+    auto_monitor: true
+    read_pv: X10SA-ES-MS:Stats1:Sigma_RBV
+  deviceTags:
+  - scam
+  enabled: true
+  onFailure: buffer
+  readOnly: true
+  readoutPriority: monitored
 samcam_x:
   description: Sample Camera X
   deviceClass: ophyd.EpicsSignalRO
@@ -1288,18 +1336,6 @@ samcam_ysig:
   onFailure: buffer
   readOnly: true
   readoutPriority: monitored
-samcam_ecc:
-  description: Sample Camera Eccentricity
-  deviceClass: ophyd.EpicsSignalRO
-  deviceConfig:
-    auto_monitor: true
-    read_pv: X10SA-ES-MS:Stats5:Eccentricity_RBV
-  deviceTags:
-  - scam
-  enabled: true
-  onFailure: buffer
-  readOnly: true
-  readoutPriority: monitored
 scam_zoom:
   description: Sample Camera Zoom
   deviceClass: ophyd_devices.EpicsMotor
@@ -1328,8 +1364,8 @@ smargon:
   deviceConfig:
     prefix: http://x10sa-smargopolo.psi.ch:3000
   deviceTags:
-  - smargon
   - motors
+  - smargon
   enabled: true
   onFailure: buffer
   readoutPriority: baseline
@@ -1424,6 +1460,18 @@ ss_bpmsum:
   onFailure: buffer
   readOnly: true
   readoutPriority: monitored
+transm:
+  description: Beam Transmission
+  deviceClass: ophyd_devices.EpicsSignal
+  deviceConfig:
+    read_pv: 'X10SA-ES-SSFI:TRANSM-SET'
+    write_pv: 'X10SA-ES-SSFI:TRANSM-SET'
+  deviceTags:
+  - transm
+  enabled: true
+  onFailure: retry
+  readOnly: false
+  readoutPriority: baseline
 ss_f1_x:
   description: SS Filter 1 X
   deviceClass: ophyd_devices.EpicsMotor
diff --git a/pxii_bec/device_configs/x10sa_device_config.yaml b/pxii_bec/device_configs/x10sa_device_config.yaml
index 8e54df4..6778a92 100644
--- a/pxii_bec/device_configs/x10sa_device_config.yaml
+++ b/pxii_bec/device_configs/x10sa_device_config.yaml
@@ -1,5 +1,17 @@
 base_config:
   - !include ./pxii-devices.yaml
+
+microscope_stats:
+  readoutPriority: async
+  description: Stats Signal of MS camera
+  deviceClass: pxii_bec.devices.stats_signal.StatsSignal
+  deviceConfig:
+     prefix: 'X10SA-ES-MS:'
+  onFailure: retry
+  enabled: true
+  readOnly: false
+  softwareTrigger: true
+
 # sample_env:
 #   - !include ./pxii-state-devices.yaml
 # id_gap:
diff --git a/pxii_bec/devices/stats_signal.py b/pxii_bec/devices/stats_signal.py
index 92dbc86..f9ac6da 100644
--- a/pxii_bec/devices/stats_signal.py
+++ b/pxii_bec/devices/stats_signal.py
@@ -1,14 +1,18 @@
-from ophyd_devices import AsyncMultiSignal
-
-from ophyd_devices.interfaces.base_classes.psi_device_base import PSIDeviceBase
-from ophyd import Component as Cpt
-from ophyd import EpicsSignalRO, Kind
+import time
 from collections import defaultdict
 from threading import RLock
-import time
+
+import numpy as np
+from bec_lib.logger import bec_logger
+from ophyd import Component as Cpt
+from ophyd import EpicsSignalRO, Kind, Signal
+from ophyd_devices import AsyncMultiSignal, StatusBase, TransitionStatus
+from ophyd_devices.interfaces.base_classes.psi_device_base import PSIDeviceBase
+
+logger = bec_logger.logger
 
 
-class StatsSignal(PSIDeviceBase, AsyncMultiSignal):
+class StatsSignal(PSIDeviceBase):
 
     data = Cpt(
         AsyncMultiSignal,
@@ -27,7 +31,7 @@ class StatsSignal(PSIDeviceBase, AsyncMultiSignal):
 
     def on_init(self):
         self._rlock = RLock()
-        self._cached_average_data = defaultdict(float)
+        self._cached_average_data = defaultdict(list)
 
     def on_connected(self):
         self._cached_average_data.clear()
@@ -37,26 +41,44 @@ class StatsSignal(PSIDeviceBase, AsyncMultiSignal):
     def _update_average_data(self, value, obj, **kwargs):
         if value is None:
             return
+        if self.scan_info.msg is None:
+            return
+        if self.scan_info.msg.status != "open":
+            return
+
         sig_name = obj.name
+        logger.info(f"Received update from signal {sig_name} with value{ value}")
 
         with self._rlock:
-            self._cached_average_data[sig_name] += value
+            self._cached_average_data[sig_name].append(value)
 
     def on_trigger(self):
         with self._rlock:
             self._cached_average_data.clear()
-        time.sleep(
-            self.scan_info.msg.scan_parameters.get("exp_time")
-        )  # Wait for the exposure time to finish
+        status = StatusBase(obj=self)  # TransitionStatus(self.acquiring, transitions=[1, 0])
+        status.add_callback(self._publish_data)
+        time.sleep(self.scan_info.msg.scan_parameters.get("exp_time", 0))
+        status.set_finished()
+        return status
 
-        with self._rlock:
-            sigma_x = self._cached_average_data.get("sigma_x_raw", 0.0)
-            sigma_y = self._cached_average_data.get("sigma_y_raw", 0.0)
-            eccentricity = self._cached_average_data.get("eccentricity_raw", 0.0)
+    def _publish_data(self, status, **kwargs):
+        if status.done:
+            with self._rlock:
 
-        ret = {
-            "sigma_x": {"value": sigma_x},
-            "sigma_y": {"value": sigma_y},
-            "eccentricity": {"value": eccentricity},
-        }
-        self.data.put(ret)
+                sigma_x = float(
+                    np.average(self._cached_average_data.get(self.sigma_x_raw.name, [0]))
+                )
+                sigma_y = float(
+                    np.average(self._cached_average_data.get(self.sigma_y_raw.name, [0]))
+                )
+                eccentricity = float(
+                    np.average(self._cached_average_data.get(self.eccentricity_raw.name, [0]))
+                )
+
+            ret = {
+                "sigma_x": {"value": sigma_x},
+                "sigma_y": {"value": sigma_y},
+                "eccentricity": {"value": eccentricity},
+            }
+            logger.info(f"Setting data for {self.name} data signal to {ret}.")
+            self.data.put(ret)
diff --git a/pxii_bec/macros/APscripts.py b/pxii_bec/macros/APscripts.py
index 871c3c0..1dd036e 100644
--- a/pxii_bec/macros/APscripts.py
+++ b/pxii_bec/macros/APscripts.py
@@ -313,9 +313,10 @@ def justfit(data_x, data_y, model="gauss", ibg=0):
     # diagnostics
     # print(f'Gfit: {g.params}')
 
-    print(f'Center of {model} fit: {g.pars["center"].value:.5f}')
-    print(f'Sigma: {g.pars["sigma"].value:.5f}')
-    print(f'FWHM: {g.pars["fwhm"].value:.5f}')
+    #print(f'Center of {model} fit: {g.pars["center"].value:.5f}')
+    #print(f'Sigma: {g.pars["sigma"].value:.5f}')
+    #print(f'FWHM: {g.pars["fwhm"].value:.5f}')
+    print(g.fit_report())
     print(f"Position of maximum: {xm}")
 
     return g, xm
@@ -338,13 +339,21 @@ def fit_plothist(hindex: int, signal_name: str, model="gauss", ibg=0):
     from lmfit.models import LinearModel, GaussianModel, VoigtModel, QuadraticModel, LorentzianModel
     import matplotlib.pyplot as plt
 
-    h = bec.history[hindex]
+    if hindex < 0:
+        h = bec.history[hindex]
+        print('by scan number')
+    else:
+        h = bec.history.get_by_scan_number(hindex)
+        print('by scan number')
     md = h.metadata["bec"]
     scanvar = list(md["args"].keys())[0]  # string, returns the variable of the last performed scan
     # data = h.devices[device_name][signal].read()["value"]
     # data_x = h.devices.dcm_pitch.dcm_pitch.read()["value"]
     # data_y = h.devices.lu_bpmsum.lu_bpmsum.read()["value"]
-    # data_x = h.devices[device_name][device_name].read()["value"]
+
+    ## that should work in any case:
+    device_name = scanvar
+    #data_x = h.devices[device_name][device_name].read()["value"]
 
     data_xd = md["positions"]  # last scan knows which device ... however, double array
     data_x = data_xd.flatten()
@@ -352,6 +361,10 @@ def fit_plothist(hindex: int, signal_name: str, model="gauss", ibg=0):
     data_y = h.devices[signal_name][signal_name].read()["value"]
     print("dy = ", data_y)
 
+    prefix='v_'
+    if model == 'gauss':
+         prefix='g_'
+    
     g, xm = justfit(data_x, data_y, model=model, ibg=ibg)
 
     plt.ion()
@@ -362,7 +375,9 @@ def fit_plothist(hindex: int, signal_name: str, model="gauss", ibg=0):
     plt.ylabel(signal_name)
     plt.show()
 
-    gcen = g.pars["center"].value
+    # g.params
+    returncenvar = prefix+"center"
+    gcen = g.params[returncenvar].value
 
     return gcen, xm, data_x, data_y
 
@@ -499,7 +514,7 @@ def fit_plot(data_x, data_y, model="gauss", ibg=1, fitrange=0, fitclick=0):
 
     fit_id = pre+'_center'
     
-    center = g.pars[fit_id].value
+    center = g.params[fit_id].value
     return center, xm 
 
 
@@ -550,7 +565,7 @@ def save_data(hindex: int, device_name: str, signal_name: str, isave = 1):
 #########################################
 ### just retrieve the saved data from csv
 ########################################
-def read_data(filename: str):
+def read_data_csv(filename: str):
     """
     Get data stored in a CSV file
 
@@ -872,7 +887,7 @@ def colliscan(direction: str, range=0.3, nsteps=30, stime=0.5, centre=1):
         settling time: default = 0.5
         move to centre or not [0/1]
 
-    Example: colliscan("h", centre=1)
+    Example: colliscan("h", range=0.2, nsteps=30, stime=0.5,centre=1)
 
     """
     import sys
@@ -972,13 +987,14 @@ def slitscan(device_location: str, direction: str, range: 1, nsteps=50, centre=0
         default_v = 1.0  ## 1.2
         det = dev.lu_bpmsum
         if direction in ["x","h"]:
-            mot = dev.fe_sxcen
-            size = dev.fe_sxsize
+            mot = dev.fe_sl_xcen
+            size = dev.fe_sl_xsize
             s_closed = 0.1
             s_open = default_h
         else:
-            mot = dev.fe_sycen
-            size = dev.fe_sysize
+            mot = dev.fe_sl_ycen
+            size = dev.fe_sl_ysize
+            
             s_closed = 0.1
             s_open = default_v
 
@@ -1072,9 +1088,9 @@ def slitscan(device_location: str, direction: str, range: 1, nsteps=50, centre=0
     # FE slits ================================
     data_y = s.scan.live_data.lu_bpmsum.lu_bpmsum.val
 
-    if mot.name == "fe_sxcen":
+    if mot.name == "fe_sl_xcen":
         data_x = s.scan.live_data.fe_sxcen.fe_sxcen.val
-    if mot.name == "fe_sycen":
+    if mot.name == "fe_sl_ycen":
         data_x = s.scan.live_data.fe_sycen.fe_sycen.val
     # BSF slits ================================
     if mot.name == "bsf_sl_xcen":
@@ -1486,7 +1502,7 @@ def scan_eg(erange, nsteps=50, fit=True):
         data_e = energy_data
         data_y = bpm_data
 
-        dirname = "/home/gac-x10sa/Data/"
+        dirname = "/sls/x10sa/config/commissioning/Data"
         # writing output to simple data file for later analysis:
         combined = np.column_stack((data_x, data_e, data_y))
         tname = f"Scan_{timestamp}.txt"
@@ -1533,3 +1549,33 @@ def scan_window(wname="Scan", fit=True):
     
     ### if done, remove 
     #w2.remove()
+################################################
+### beam centre on detector with distance 
+################################################
+
+#/sls/x10sa/config/commissioning/Data
+
+def readdist(filename):
+    with open(filename, "r") as f:
+        combined_data = np.loadtxt(f, delimiter=",")
+        
+    ind1 = 0
+    ind2 = 1
+    ind3 = 2
+    detdist = combined_data[:, ind1]
+    bposx = combined_data[:, ind2]
+    bposy = combined_data[:, ind3]
+    
+    plt.ion()
+    plt.figure(1)
+    plt.plot(detdist,bposx, ".")
+    plt.title("Bpos X vs Detectordist")
+    plt.xlabel("dist / mm")
+    plt.ylabel("bpos X in det pix ")
+    
+    plt.figure(2)
+    plt.plot(detdist,bposy, ".")
+    plt.title("Bpos Y vs Detectordist")
+    plt.xlabel("dist / mm")
+    plt.ylabel("bpos Y in det pix ")
+    plt.show()
diff --git a/pxii_bec/macros/KB_Class.py b/pxii_bec/macros/KB_Class.py
index 7c5dccb..51e2c7f 100644
--- a/pxii_bec/macros/KB_Class.py
+++ b/pxii_bec/macros/KB_Class.py
@@ -13,10 +13,10 @@ class KBMirror:
             self.P2 = 0.583490
             self.P3 = -0.372950
             self.L1 = -0.955798
-            self.curv_approx_fact = 1.0121
-            self.s1 = 24.15   #24.35
+            self.curv_approx_fact =  0.9550 # 0.9550 # origin: 1.0121 
+            self.s1 = 24.15   
             self.s2 = 1.45
-            self.asym0 = 0.08
+            self.asym0 = 0.114 ## origin: 0.0.08  + 0.034  # user asymm, eg 0.4750, will be added ( - sign) 
 
         elif self.mirror_type == "PXII_HFM":
             self.C0 = 0.060453
@@ -24,10 +24,10 @@ class KBMirror:
             self.P2 = 0.586057
             self.P3 = -0.088093
             self.L1 = -0.999445
-            self.curv_approx_fact = 0.995
-            self.s1 = 24.90   #25.10
+            self.curv_approx_fact = 0.8810   # origin: # 0.995 
+            self.s1 = 24.90   
             self.s2 = 0.7
-            self.asym0 = -0.0107
+            self.asym0 = 0.06455    ## origin: -0.0107 + 0.07525          # added user asymm 2.75 (-sign)  
 
         else:
             raise ValueError(f"Unknown mirror type: {self.mirror_type}")
diff --git a/pxii_bec/macros/calculator.py b/pxii_bec/macros/calculator.py
index b10a22a..a6278e8 100755
--- a/pxii_bec/macros/calculator.py
+++ b/pxii_bec/macros/calculator.py
@@ -17,6 +17,10 @@ class Constants:
     # d-spacings
     d_spacing = {120: 3.13481, 298: 3.13562}
 
+    # offset
+    # bragg_offset = -0.1457
+    bragg_offset = 0
+
 
 def speed_of_light_ang():
     """
@@ -116,7 +120,9 @@ def calculate_wavelength_from_angle(bragg_angle_mrad: float, temp=120) -> float:
         The calculated wavelength as a float value.
     """
     d = Constants.d_spacing[temp]
-    return 2 * d * np.sin(bragg_angle_mrad / 1000)
+    fudge = bragg_angle_mrad + Constants.bragg_offset
+    # return 2 * d * np.sin((bragg_angle_mrad / 1000)
+    return 2 * d * np.sin(fudge / 1000)
 
 
 def calculate_energy_from_wavelength(wavelength: float) -> float:
diff --git a/pxii_bec/macros/calibrate_bpm.py b/pxii_bec/macros/calibrate_bpm.py
new file mode 100644
index 0000000..f041273
--- /dev/null
+++ b/pxii_bec/macros/calibrate_bpm.py
@@ -0,0 +1,484 @@
+"""
+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
+import numpy as np
+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
+
+
+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)
+   
+    # 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"]
+    # print(f"bpm1 = {bpm1:.3f}, bpm2 = {bpm2:.3f}, bpm3 = {bpm3:.3f}, bpm4 = {bpm4:.3f}")
+    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 xn, 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)
+
+
+# Read data from BPM calibration files
+# def fit_bpm_data(x_y):
+#     """
+#     Read data from bpm calibration files.
+#     """
+#     filedir = "luts/bpm_calib/"
+#     rundate= "02072130_"
+#     name = f"bcu_bpm_{x_y}.csv"
+#     filename = filedir + rundate + name
+#     data = np.loadtxt(filename, delimiter=",", skiprows=1)
+#     xdata = data[:, 0]
+#     ydata = data[:, 1]
+#     return xdata, ydata
+
+
+
+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.65
+    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("HFM Lat Position")
+    plt.ylabel("SamCam X")
+    plt.title("Zoom Calibration")
+    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"The flux is {flux: .3g} ph/s")
+
+    return float(flux)
+
+def calibrate_bpm_flux():
+    # d['det_z'].move('mse')
+    # d['det_y'].move('xi')
+    # d['det_z'].move(200)
+
+    flux = flux(dev.xidiode.read()['xidiode']['value'])
+
+
+def zero_xidiode():
+    vals = []
+    for i in range(1,101):
+        reading = dev.xidiode.read()['xidiode']['value']
+        vals.append(reading)
+
+    zero = np.mean(vals)
+    print(zero)
+
+
+    
diff --git a/pxii_bec/macros/chmirrh.py b/pxii_bec/macros/chmirrh.py
index 7131fb5..94b58ce 100644
--- a/pxii_bec/macros/chmirrh.py
+++ b/pxii_bec/macros/chmirrh.py
@@ -2,60 +2,9 @@ import numpy as np
 import math, time 
 import matplotlib.pyplot as plt
 import scipy
-#########################################
-### just retrieve the saved data from csv
-########################################
-def read_data(filename: str):
-    """
-    Get data stored in a CSV file
-
-    Args:
-        filename (str): the csv file, eg, of a Scan
-        e.g., dat = read_data("/home/e18747/SLS2/Data/gaps/gaps10.txt")
-    """
-
-    # dirname = '/home/gac-x10sa/Data/'
-    # read output to simple data file for later analysis:
-
-    with open(filename, "r") as f:
-        combined_data = np.loadtxt(f, delimiter=",")  # no header and fmt for load
-
-    rows, cols = combined_data.shape
-
-    print("No of Rows =", rows)
-    print("No of Colums =", cols)
-    ind1 = 0
-    ind2 = 1
-    if cols > 2:
-        print(
-            "only first 2 colums (data[:, 0] and  data[:, 1]) are plotted, please consider the other columns as well!"
-        )
-        index1 = 0
-        index2 = 1
-        index1, index2 = input(
-            "please enter the columns to be plotted vs first one [eg: 1,2] :  "
-        ).split(",")
-        ind1 = int(index1)
-        ind2 = int(index2)
-
-        if ind2 > cols - 1:
-            print(" no valid index")
-
-    data_x = combined_data[:, ind1]
-    data_y = combined_data[:, ind2]
-
-    plt.ion()
-    plt.figure()
-    plt.plot(data_x, data_y, ".")
-
-    plt.title(f"Scan from {filename}")
-    plt.xlabel(f"data column {ind1}")
-    plt.ylabel(f"data column {ind2}")
-
-    plt.show()
-
-    return combined_data
-
+########################################################
+### scan FE slits and check BPMs before and after mirror 
+########################################################
 
 def gaussfit(x,y,p):
 # define a gaussian fitting function where
@@ -95,22 +44,42 @@ def gaussfit(x,y,p):
 
     return bestfit, p1
 
-def scan_fe_sl():
-    fe_v_size0 = dev.fe_sl_ysize.user_readback.get()
-    umv(dev.fe_sl_xsize, 0.1)
-    s = scans.line_scan(dev.fe_sl_ycen,-1,1, exp_time=02., steps=21, relative=True)
-    time.sleep(0.1)
-    umv(dev.fe_sl_xsize, fe_v_size0)
+def scan_fe_sl(dir = "y"):
+    if dir in ["y","v"]:
+        sdir ="vertical"
+        fe_size0 = dev.fe_sl_ysize.user_readback.get()
+        print('Closing vertical FE slits to 0.1')  
+        umv(dev.fe_sl_ysize, 0.1)
+        s = scans.line_scan(dev.fe_sl_ycen,-1,1, exp_time=02., steps=21, relative=True)
+        time.sleep(0.1)
+        umv(dev.fe_sl_ysize, fe_size0)
+        print(f'Opening {sdir} FE slits to {fe_size0}')
+    if dir in ["x","h"]:
+        sdir ="horizontal"
+        fe_size0 = dev.fe_sl_xsize.user_readback.get()
+        print('Closing horizontal FE slits to 0.1')  
+        umv(dev.fe_sl_xsize, 0.1)
+        s = scans.line_scan(dev.fe_sl_xcen,-1,1, exp_time=02., steps=21, relative=True)
+        time.sleep(0.1)
+        umv(dev.fe_sl_xsize, fe_size0)
+        print(f'Opening {sdir} FE slits to {fe_size0}')
+    else:
+        print('no valid direction, doing nothing') 
     
-def mirrh():
+def mirrh(dir='y'):
     
-    #dat_ss = read_data("/home/anuschka/SS_pitch2p94.csv")
-    #dat_bcu = read_data("/home/anuschka/BCU_pitch2p94.csv")
+    #dat_ss = read_data_csv("/home/anuschka/SS_pitch2p94.csv")
+    #dat_bcu = read_data_csv("/home/anuschka/BCU_pitch2p94.csv")
 
-    scan_fe_sl()
+    scan_fe_sl(dir=dir)
     ind = -1
-    dat_ss = save_data(ind, 'fe_sl_ycen', 'ss_bpmsum', isave = 0)
-    dat_bcu = save_data(ind, 'fe_sl_ycen', 'bcu_bpmsum', isave = 0)
+    if dir in ["x","h"]:
+        slitmot = 'fe_sl_xcen'
+    if dir in ["y","v"]:
+        slitmot = 'fe_sl_ycen'        
+    
+    dat_ss = save_data(ind, slitmot, 'ss_bpmsum', isave = 0)
+    dat_bcu = save_data(ind, slitmot, 'bcu_bpmsum', isave = 0)
     
 
     data_xs = dat_ss[0]      # csv would be: dat_ss[:, 0]
@@ -124,7 +93,7 @@ def mirrh():
     #plt.plot(data_xs, data_ys, "*")
     #plt.plot(data_xb, data_yb, "+")
 
-    data_ys_norm = data_ys/max(data_ys)
+    data_ys_norm = data_ys/max(data_ys) # looks better than data_ys/sum(data_ys)
     data_yb_norm = data_yb/max(data_yb)
 
     plt.figure()
@@ -138,7 +107,7 @@ def mirrh():
     xmi = data_xs[xm]
     #print('xm, xmi', xm, xmi) 
 
-    pf1 = [1, xmi , np.std(data_ys_norm),0]
+    pf1 = [1, xmi , np.std(data_ys_norm), np.mean(y[0:5])]     # bg = 0 
     fit_s, ps = gaussfit(data_xs,data_ys_norm, pf1) 
 
     y= data_yb_norm
@@ -146,7 +115,7 @@ def mirrh():
     xmi = data_xb[xm]
     #print('xm, xmi', xm, xmi) 
 
-    pf2 = [1, xmi, np.std(data_yb_norm),-11] 
+    pf2 = [1, xmi, np.std(data_yb_norm), np.mean(y[0:5])]  # bg = -11
     fit_b, pb = gaussfit(data_xb,data_yb_norm, pf2) 
 
 
diff --git a/pxii_bec/macros/chmirrh.py~ b/pxii_bec/macros/chmirrh.py~
deleted file mode 100644
index 664ab60..0000000
--- a/pxii_bec/macros/chmirrh.py~
+++ /dev/null
@@ -1,142 +0,0 @@
-import numpy as np 
-import math
-import matplotlib.pyplot as plt
-import scipy
-#########################################
-### just retrieve the saved data from csv
-########################################
-def read_data(filename: str):
-    """
-    Get data stored in a CSV file
-
-    Args:
-        filename (str): the csv file, eg, of a Scan
-        e.g., dat = read_data("/home/e18747/SLS2/Data/gaps/gaps10.txt")
-    """
-
-    # dirname = '/home/gac-x10sa/Data/'
-    # read output to simple data file for later analysis:
-
-    with open(filename, "r") as f:
-        combined_data = np.loadtxt(f, delimiter=",")  # no header and fmt for load
-
-    rows, cols = combined_data.shape
-
-    print("No of Rows =", rows)
-    print("No of Colums =", cols)
-    ind1 = 0
-    ind2 = 1
-    if cols > 2:
-        print(
-            "only first 2 colums (data[:, 0] and  data[:, 1]) are plotted, please consider the other columns as well!"
-        )
-        index1 = 0
-        index2 = 1
-        index1, index2 = input(
-            "please enter the columns to be plotted vs first one [eg: 1,2] :  "
-        ).split(",")
-        ind1 = int(index1)
-        ind2 = int(index2)
-
-        if ind2 > cols - 1:
-            print(" no valid index")
-
-    data_x = combined_data[:, ind1]
-    data_y = combined_data[:, ind2]
-
-    plt.ion()
-    plt.figure()
-    plt.plot(data_x, data_y, ".")
-
-    plt.title(f"Scan from {filename}")
-    plt.xlabel(f"data column {ind1}")
-    plt.ylabel(f"data column {ind2}")
-
-    plt.show()
-
-    return combined_data
-
-
-def gaussfit(x,y,p):
-# define a gaussian fitting function where
-# p[0] = amplitude
-# p[1] = mean
-# p[2] = sigma
-# p[3] = const offset 
-    
-    fitfct = lambda p, x0: p[0]*np.exp(-(x0-p[1])**2/(2.0*p[2]**2))+p[3]
-    #fitfct = lambda p, x0: p[0]*np.exp(-(x0-p[1])**2/(2.0*p[2]**2))
-    errfct = lambda p, x0,y0: fitfct(p,x0)-y0  
-
-    # guess fit parameters
-    ind0 = np.where(y==max(y))[0] 
-    p0=scipy.c_[max(y), x[ind0], np.std(y), np.mean(y[0:5])]
-    #p0=scipy.c_[max(y), np.where(y==max(y))[0], np.std(y)]
-                  # array generation [[]], and taking only one dimens
-    # fit a gaussian to the correlation function
-
-    print('p0 = ' , p0)
-    p1, success = scipy.optimize.leastsq(errfct, p0.copy()[0],
-                                     args=(x,y))
-
-    plt.figure()
-    bestfit=fitfct(p1,x)
-    #print('************************', x,y)
-    plt.plot(x,y,'r*')    
-    plt.plot(x,bestfit,'r-')
-
-    # different fit fct
-    #popt, pcov = curve_fit(fitfct, x, y)
-    #newfit=fitfct(popt,x)
-    #plt.plot(x,newfit,'g+')
-
-    
-    plt.show()
-    #time.sleep(1)
-    #plt.close()
-    return bestfit, p1
-
-def mirh():
-    
-    dat_ss = read_data("/home/anuschka/SS_pitch2p94.csv")
-    dat_bcu = read_data("/home/anuschka/BCU_pitch2p94.csv")
-
-    data_xs = dat_ss[:, 0]
-    data_ys = dat_ss[:, 1]
-    
-    data_xb = dat_bcu[:, 0]
-    data_yb = dat_bcu[:, 1]
-
-    plt.figure()
-
-    plt.plot(data_xs, data_ys, "*")
-    plt.plot(data_xb, data_yb, "+")
-
-    data_ys_norm = data_ys/max(data_ys)
-    data_yb_norm = data_yb/max(data_yb)
-
-    plt.figure()
-
-    plt.plot(data_xs, data_ys_norm, "*")
-    plt.plot(data_xb, data_yb_norm, "+")
-
-
-    pf1 = [1, 0.8 , np.std(data_ys_norm),0]
-
-    y= data_ys_norm
-    xm = np.where(y==max(y))[0]
-    xmi = data_xs[xm]
-    print('xm, xmi', xm, xmi) 
-    
-    fit_s, ps = gaussfit(data_xs,data_ys_norm, pf1) 
-
-    pf2 = [1, 0.8, np.std(data_yb_norm),-11] 
-    fit_b, pb = gaussfit(data_xb,data_yb_norm, pf2) 
-
-
-    plt.figure()
-
-    plt.plot(data_xs, fit_s-ps[3], "*")
-    plt.plot(data_xb, fit_b-pb[3], "+")
-  
-    
diff --git a/pxii_bec/macros/katscripts.py b/pxii_bec/macros/katscripts.py
index 9dff871..9b76ed9 100755
--- a/pxii_bec/macros/katscripts.py
+++ b/pxii_bec/macros/katscripts.py
@@ -1,7 +1,7 @@
 """Some simple scripts for commissioning"""
 
 import time, os
-from datetime import datetime
+from datetime import datetime, timedelta
 from dataclasses import dataclass
 import numpy as np
 import matplotlib.pyplot as plt
@@ -10,15 +10,6 @@ import pandas as pd
 from scipy.optimize import linear_sum_assignment
 
 
-def goto1a():
-    """Set beamline to 1A position"""
-    bragg_angle = 160.117
-    perp_dist = 0.0389
-    gap = 5.319
-    umv(dev.dcm_bragg, bragg_angle)
-    umv(dev.dcm_perp, perp_dist)
-    umv(dev.id_gap, gap)
-
 
 def bragg_scan():
     """Energy scans from 6 to 18 keV at various gap positions"""
@@ -59,26 +50,10 @@ def energy_checks(start_energy=8000, end_energy=20000):
     """Check gap, pitch and intensities at different energies"""
     f = open("luts/test.csv", "a", encoding="utf-8")
 
-    gaps = {
-        8000: 8.04,
-        9000: 9.55,
-        10000: 7.45,
-        11000: 8.46,
-        12000: 7.09,
-        13000: 7.77,
-        14000: 8.61,
-        15000: 9.57,
-        16000: 8.09,
-        17000: 7.18,
-        18000: 7.67,
-        19000: 8.24,
-        20000: 7.52,
-    }
     for energy in range(start_energy, end_energy, 1000):
         f = open("luts/test.csv", "a", encoding="utf-8")
         print(f"Moving to {energy} eV")
-        umv(dev.id_gap, gaps[energy])
-        move_and_set_dcm(energy)
+        bl_energy
         if energy < 16000:
             close_fe_slits(True)
             go_to_peak(
@@ -277,86 +252,75 @@ def write_beamsize():
         f.write(f"{lens_config},{sx:.4g},{sy:.4g},{i1_reading:.4g}\n")
 
 
-def beamsize(size):
-    """WIP to dial up different beam sizes"""
-    # Set mirror trans and pitch to default settings for 12.4 keV
-    # vlaues from a 2025
-    # umv(dev.vfm_vert,2.136)
-    # umv(dev.vfm_pitch, 2.695)
-    # umv(dev.hfm_lat, 1.64)
-    # umv(dev.hfm_pitch, 2.9844)
-    umv(dev.vfm_vert, 2.5301)
-    umv(dev.vfm_pitch, 2.60)
-    umv(dev.hfm_lat, 1.6586)
-    umv(dev.hfm_pitch, 2.9964)
-    print(f"changing beam size to {size}")
-
-    # values = {"xs": {"vbu": -0.9095, "vbd": -1.246, "hbu": -2.9375, "hbd": -5.698},
-    #           "s": {"vbu": -0.8365, "vbd": -1.22, "hbu": -2.804, "hbd": -5.473},
-    #           "m": {"vbu": -0.7865, "vbd": -1.17, "hbu": -2.604, "hbd": -5.273},
-    #           "l": {"vbu": -0.7459, "vbd": -1.13, "hbu": -2.504, "hbd": -5.173}}
-    # values = {"xs": {"vbu": -0.9095, "vbd": -1.246, "hbu": -2.8075, "hbd": -5.668},
-    #           "s": {"vbu": -0.8606, "vbd": -1.195, "hbu": -2.6975, "hbd": -5.5585},
-    #           "m": {"vbu": -0.8410, "vbd": -1.465, "hbu": -2.5575, "hbd": -5.418},
-    #           "l": {"vbu": -0.762, "vbd": -1.094, "hbu": -2.4582, "hbd": -5.3187}}
-    values = {
-        "xs": {"vbu": -0.741, "vbd": -1.206, "hbu": -2.869, "hbd": -5.819},
-        "s": {"vbu": -0.7101, "vbd": -1.176, "hbu": -2.7133, "hbd": -5.5535},
-        "m": {"vbu": -0.6601, "vbd": -1.126, "hbu": -2.5133, "hbd": -5.3535},
-        "l": {"vbu": -0.6196, "vbd": -1.086, "hbu": -2.4133, "hbd": -5.2535},
-    }
-    umv(
-        dev.vfm_bu,
-        values[size]["vbu"],
-        dev.vfm_bd,
-        values[size]["vbd"],
-        dev.hfm_bu,
-        values[size]["hbu"],
-        dev.hfm_bd,
-        values[size]["hbd"],
-    )
-    print("Setting camera exposure time")
-    auto_exposure()
-    xsig = dev.samcam_xsig.read()["samcam_xsig"]["value"]
-    ysig = dev.samcam_ysig.read()["samcam_ysig"]["value"]
-    print(f"Beam has Xsigma = {xsig:.3g}, Ysigma = {ysig:.3g}")
-
-    # if size == 'm':
-    #     move_to_x()
-    #     print(f"Beam has Xsigma = {xsig:.3g}, Ysigma = {ysig:.3g}")
-    #     return
-    # if size == 'l':
-    #     move_to_x()
-    #     move_to_y()
-    #     umv(dev.hfm_pitch, 2.9684)
-    #     print(f"Beam has Xsigma = {xsig:.3g}, Ysigma = {ysig:.3g}")
-    #     return
-    # else:
-    #     move_to_x()
-    #     move_to_y()
-    #     print(f"Beam has Xsigma = {xsig:.3g}, Ysigma = {ysig:.3g}")
-
 
 def scan_Ge():
     scanbefore = 30  # eV
     scanafter = 30  # eV
-    stepsize = 0.5  # eV
-    edge_Ge = 11103.1  # eV
-
-    scanstart = edge_Ge + scanafter
-    scanend = edge_Ge - scanbefore
-    bragg_start = convert_from_energy(scanstart)["bragg_angle_mrad"]
-    bragg_end = convert_from_energy(scanend)["bragg_angle_mrad"]
-    stepnumber = abs(int((scanend - scanstart) / stepsize))
-
-    print(f"Changing energy to {scanstart}")
+    stepsize = 0.2  # eV
+    # Edge sgould be at 11103.1
+    actual_edge_Ge = 11103.1  # eV
+    actual_bragg = convert_from_energy(actual_edge_Ge)["bragg_angle_mrad"]
+    current_edge = 11095
+    current_bragg = convert_from_energy(current_edge)["bragg_angle_mrad"]
+    print(f"Ge edge should be at a bragg angle of {actual_bragg:.5g} mrad\n")
+    print(f"Current estimate of edge is {current_bragg:.5g} mrad\n")
+    xdata = []
+    ydata = []
+    filedir, fname = fnow()
+    filename = f"{filedir}{fname}_ge_edge.csv"
+    print(filename)
+    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"]
+    # bragg_end = convert_from_energy(scanend)["bragg_angle_mrad"]
+    # stepnumber = abs(int((scanend - scanstart) / stepsize))
+    # print(f"Changing energy to {scanstart}")
     # bl_energy(scanstart)
-    # scan_gap()
-    umv(dev.bsc_f4_x, 6)  # Put Ge filter in place
-    print(f"scanning from {bragg_start} to {bragg_end} with {stepnumber} steps")
-    s = scans.line_scan(dev.dcm_bragg, bragg_start, bragg_end, steps=stepnumber, relative=False)
-    print(s)
-    umv(dev.bsc_f4_x, 2)
+    # # scan_gap()
+    # umv(dev.ss_f4_x, 6)  # Put Ge filter in place
+    # print(f"scanning from {bragg_start} to {bragg_end} with {stepnumber} steps")
+    # s = scans.line_scan(dev.dcm_bragg, bragg_start, bragg_end, steps=stepnumber, relative=False)
+    # print(s)
+    # umv(dev.ss_f4_x, 2)
+
+    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"]
+    bl_energy(start_energy)
+    umv(dev.ss_f4_x, 6)  # Put Ge filter in place
+
+    for bragg in scanpoints:
+        scans.umv(dev.dcm_bragg, bragg, relative=False)
+        energy = convert_from_bragg(bragg)["energy_ev"]
+        inten = average_inten()
+        xdata.append(energy)
+        ydata.append(inten)
+        with open(filename, 'a') as f:
+            f.write(f"{energy:.2f},{bragg:.2f},{inten:.2f}\n")
+    data = {
+        "x_data": np.array(xdata),
+        "y_data": np.array(ydata),
+        "motor_name": "dcm_bragg",
+        "signal_name": "bcu_bpmsum",
+        "motor_device": dev.dcm_bragg,
+        "scan_number": "Current",
+    }
+
+    # Define and fit model to scan data
+    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():
@@ -434,79 +398,79 @@ def cal_samcam_fixed_x():
     plt.show()
 
 
-def cal_samcam_fixed_y():
-    vert = dev.vfm_vert.read()["vfm_vert"]["value"]
-    pitch = dev.vfm_pitch.read()["vfm_pitch"]["value"]
-    scanpoints = np.linspace(vert - 0.1, vert + 0.1, 7)
-    zoom_high = 800
-    zoom_low = 200
-    highs, lows = [], []
-    high_vert, low_vert = [], []
+# def cal_samcam_fixed_y():
+#     vert = dev.vfm_vert.read()["vfm_vert"]["value"]
+#     pitch = dev.vfm_pitch.read()["vfm_pitch"]["value"]
+#     scanpoints = np.linspace(vert - 0.1, vert + 0.1, 7)
+#     zoom_high = 800
+#     zoom_low = 200
+#     highs, lows = [], []
+#     high_vert, low_vert = [], []
 
-    # high zoom
-    umv(dev.scam_zoom, zoom_high)
-    auto_exposure("samcam")
-    for i in scanpoints:
-        umv(dev.vfm_vert, i)
-        umv(dev.vfm_pitch, pitch)
-        time.sleep(0.2)
-        high_y = dev.samcam_y.read()["samcam_y"]["value"]
-        high_x = dev.vfm_vert.read()["vfm_vert"]["value"]
-        print(f"Ypos at zoom {zoom_high} is {high_y}")
-        highs.append(high_y)
-        high_vert.append(high_x)
+#     # high zoom
+#     umv(dev.scam_zoom, zoom_high)
+#     auto_exposure("samcam")
+#     for i in scanpoints:
+#         umv(dev.vfm_vert, i)
+#         umv(dev.vfm_pitch, pitch)
+#         time.sleep(0.2)
+#         high_y = dev.samcam_y.read()["samcam_y"]["value"]
+#         high_x = dev.vfm_vert.read()["vfm_vert"]["value"]
+#         print(f"Ypos at zoom {zoom_high} is {high_y}")
+#         highs.append(high_y)
+#         high_vert.append(high_x)
 
-    umv(dev.scam_zoom, zoom_low)
-    auto_exposure("samcam")
-    for i in scanpoints:
-        umv(dev.vfm_vert, i)
-        umv(dev.vfm_pitch, pitch)
-        time.sleep(0.2)
-        low_y = dev.samcam_y.read()["samcam_y"]["value"]
-        low_x = dev.vfm_vert.read()["vfm_vert"]["value"]
-        print(f"Ypos at zoom {zoom_low} is {low_y}")
-        lows.append(low_y)
-        low_vert.append(low_x)
+#     umv(dev.scam_zoom, zoom_low)
+#     auto_exposure("samcam")
+#     for i in scanpoints:
+#         umv(dev.vfm_vert, i)
+#         umv(dev.vfm_pitch, pitch)
+#         time.sleep(0.2)
+#         low_y = dev.samcam_y.read()["samcam_y"]["value"]
+#         low_x = dev.vfm_vert.read()["vfm_vert"]["value"]
+#         print(f"Ypos at zoom {zoom_low} is {low_y}")
+#         lows.append(low_y)
+#         low_vert.append(low_x)
 
-    umv(dev.vfm_vert, vert)
-    umv(dev.vfm_pitch, pitch)
+#     umv(dev.vfm_vert, vert)
+#     umv(dev.vfm_pitch, pitch)
 
-    # Convert to numpy arrays
-    x_high = np.array(high_vert)
-    x_low = np.array(low_vert)
-    y_high = np.array(highs)
-    y_low = np.array(lows)
+#     # Convert to numpy arrays
+#     x_high = np.array(high_vert)
+#     x_low = np.array(low_vert)
+#     y_high = np.array(highs)
+#     y_low = np.array(lows)
 
-    # Fit linear models: y = m*x + b
-    m_high, b_high = np.polyfit(x_high, y_high, 1)
-    m_low, b_low = np.polyfit(x_low, y_low, 1)
+#     # Fit linear models: y = m*x + b
+#     m_high, b_high = np.polyfit(x_high, y_high, 1)
+#     m_low, b_low = np.polyfit(x_low, y_low, 1)
 
-    # Intersection point
-    x_intersect = (b_low - b_high) / (m_high - m_low)
-    y_intersect = m_high * x_intersect + b_high
+#     # Intersection point
+#     x_intersect = (b_low - b_high) / (m_high - m_low)
+#     y_intersect = m_high * x_intersect + b_high
 
-    print(f"\nIntersection at:")
-    print(f"  VFM_vert = {x_intersect:.5f}")
-    print(f"  SamCam_y = {y_intersect:.2f}")
+#     print(f"\nIntersection at:")
+#     print(f"  VFM_vert = {x_intersect:.5f}")
+#     print(f"  SamCam_y = {y_intersect:.2f}")
 
-    # --- Plot ---
-    plt.figure(figsize=(7, 5))
-    plt.scatter(x_high, y_high, color="red", label=f"Zoom = {zoom_high} data")
-    plt.scatter(x_low, y_low, color="blue", label=f"Zoom = {zoom_low} data")
+#     # --- Plot ---
+#     plt.figure(figsize=(7, 5))
+#     plt.scatter(x_high, y_high, color="red", label=f"Zoom = {zoom_high} data")
+#     plt.scatter(x_low, y_low, color="blue", label=f"Zoom = {zoom_low} data")
 
-    x_fit = np.linspace(min(x_high), max(x_high), 100)
-    plt.plot(x_fit, m_high * x_fit + b_high, "r--", label=f"Fit: Zoom {zoom_high}")
-    plt.plot(x_fit, m_low * x_fit + b_low, "b--", label=f"Fit: Zoom {zoom_low}")
+#     x_fit = np.linspace(min(x_high), max(x_high), 100)
+#     plt.plot(x_fit, m_high * x_fit + b_high, "r--", label=f"Fit: Zoom {zoom_high}")
+#     plt.plot(x_fit, m_low * x_fit + b_low, "b--", label=f"Fit: Zoom {zoom_low}")
 
-    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 Vert Position")
-    plt.ylabel("SamCam Y")
-    plt.title("Zoom Calibration")
-    plt.legend()
-    plt.grid(True)
-    plt.tight_layout()
-    plt.show()
+#     plt.xlabel("VFM Vert Position")
+#     plt.ylabel("SamCam Y")
+#     plt.title("Zoom Calibration")
+#     plt.legend()
+#     plt.grid(True)
+#     plt.tight_layout()
+#     plt.show()
 
 
 def correlate_froll():
@@ -620,37 +584,37 @@ def correlate_hfm_lat():
     plt.show()
 
 
-def correlate_vfm_vert():
-    y_correlation = []
-    vert = dev.vfm_vert.read()["vfm_vert"]["value"]
-    scanpoints = np.linspace(vert - 0.06, vert + 0.06, 10)
-    for pos in scanpoints:
-        umv(dev.vfm_vert, pos)
-        xpos, ypos = compute_pos()
-        y_correlation.append(ypos)
-    umv(dev.vfm_vert, vert)
+# def correlate_vfm_vert():
+#     y_correlation = []
+#     vert = dev.vfm_vert.read()["vfm_vert"]["value"]
+#     scanpoints = np.linspace(vert - 0.06, vert + 0.06, 10)
+#     for pos in scanpoints:
+#         umv(dev.vfm_vert, pos)
+#         xpos, ypos = compute_pos()
+#         y_correlation.append(ypos)
+#     umv(dev.vfm_vert, vert)
 
-    x = np.array(scanpoints)
-    y1 = np.array(y_correlation)
+#     x = np.array(scanpoints)
+#     y1 = np.array(y_correlation)
 
-    yslope, yint = np.polyfit(x, y1, 1)
+#     yslope, yint = np.polyfit(x, y1, 1)
 
-    print(f"Fit for xpos is: y = {yslope}.x + {yint}")
-    plt.figure(figsize=(7, 5))
-    # plt.scatter(x, y1, label='xpos')
-    plt.scatter(x, y1, label="beam y position", color="rebeccapurple")
+#     print(f"Fit for xpos is: y = {yslope}.x + {yint}")
+#     plt.figure(figsize=(7, 5))
+#     # plt.scatter(x, y1, label='xpos')
+#     plt.scatter(x, y1, label="beam y position", color="rebeccapurple")
 
-    x_fit = np.linspace(min(x), max(x), 100)
-    # plt.plot(x_fit, xslope*x_fit + xint, 'r--', label=f'y = {xslope}*x + {xint}')
-    plt.plot(
-        x_fit, yslope * x_fit + yint, label=f"y = {yslope:.5g}*x + ({yint:.5g})", color="lightcoral"
-    )
-    plt.xlabel("VFM vertical (mm)")
-    plt.ylabel("beam Y(from BCU BPM)")
-    plt.legend()
-    plt.grid(True)
-    plt.tight_layout()
-    plt.show()
+#     x_fit = np.linspace(min(x), max(x), 100)
+#     # plt.plot(x_fit, xslope*x_fit + xint, 'r--', label=f'y = {xslope}*x + {xint}')
+#     plt.plot(
+#         x_fit, yslope * x_fit + yint, label=f"y = {yslope:.5g}*x + ({yint:.5g})", color="lightcoral"
+#     )
+#     plt.xlabel("VFM vertical (mm)")
+#     plt.ylabel("beam Y(from BCU BPM)")
+#     plt.legend()
+#     plt.grid(True)
+#     plt.tight_layout()
+#     plt.show()
 
 
 def correlate_hfm_pitch():
@@ -703,21 +667,21 @@ def move_to_x(targetx=0.13):
     umv(dev.hfm_lat, mot_new)
 
 
-def move_to_y(targety=-0.033):
-    # Calibration and drift compensation
-    # target from previous run = -0.283
-    slope = -1.2779
-    intercept = 2.4374
+# def move_to_y(targety=-0.033):
+#     # Calibration and drift compensation
+#     # target from previous run = -0.283
+#     slope = -1.2779
+#     intercept = 2.4374
 
-    xpos, ypos = compute_pos()
+#     xpos, ypos = compute_pos()
 
-    beam_drift = ypos - targety  # beam moved +0.1 mm (100 µm)
-    mot_current = dev.vfm_vert.read()["vfm_vert"]["value"]  # current motor position (example)
+#     beam_drift = ypos - targety  # beam moved +0.1 mm (100 µm)
+#     mot_current = dev.vfm_vert.read()["vfm_vert"]["value"]  # current motor position (example)
 
-    # Compute correction
-    delta = -beam_drift / slope
-    mot_new = mot_current + delta
-    umv(dev.vfm_vert, mot_new)
+#     # Compute correction
+#     delta = -beam_drift / slope
+#     mot_new = mot_current + delta
+#     umv(dev.vfm_vert, mot_new)
 
 
 from epics import caget, caput
@@ -792,34 +756,34 @@ def auto_exposure(
     return exposure
 
 
-def get_new_beamsizes(mirror):
-    filename = f"luts/{mirror}_size.csv"
-    if mirror == "hfm":
-        bu = dev.hfm_bu
-        bd = dev.hfm_bd
-        step = 0.2
-    else:
-        bu = dev.vfm_bu
-        bd = dev.vfm_bd
-        step = 0.05
-    beamsize("xs")
-    move_targetx_with_hfm_lat()
-    move_targety_with_vfm_vert()
-    with open(filename, "w") as f:
-        for _ in range(4):
-            curr_bu = bu.user_readback.get()
-            curr_bd = bd.user_readback.get()
-            auto_exposure()
-            umv(bu, curr_bu + step)
-            umv(bd, curr_bd + step)
-            move_targetx_with_hfm_lat()
-            move_targety_with_vfm_vert()
-            xpos = dev.samcam_x.read()["samcam_x"]["value"]
-            ypos = dev.samcam_y.read()["samcam_y"]["value"]
-            xsig = dev.samcam_xsig.read()["samcam_xsig"]["value"]
-            ysig = dev.samcam_ysig.read()["samcam_ysig"]["value"]
-            f.write(f"{curr_bu}, {curr_bd}, {xpos}, {ypos}, {xsig}, {ysig}\n")
-            print(f"{curr_bu}, {curr_bd}, {xpos}, {ypos}, {xsig}, {ysig}")
+# def get_new_beamsizes(mirror):
+#     filename = f"luts/{mirror}_size.csv"
+#     if mirror == "hfm":
+#         bu = dev.hfm_bu
+#         bd = dev.hfm_bd
+#         step = 0.2
+#     else:
+#         bu = dev.vfm_bu
+#         bd = dev.vfm_bd
+#         step = 0.05
+#     beamsize("xs")
+#     move_targetx_with_hfm_lat()
+#     move_targety_with_vfm_vert()
+#     with open(filename, "w") as f:
+#         for _ in range(4):
+#             curr_bu = bu.user_readback.get()
+#             curr_bd = bd.user_readback.get()
+#             auto_exposure()
+#             umv(bu, curr_bu + step)
+#             umv(bd, curr_bd + step)
+#             move_targetx_with_hfm_lat()
+#             move_targety_with_vfm_vert()
+#             xpos = dev.samcam_x.read()["samcam_x"]["value"]
+#             ypos = dev.samcam_y.read()["samcam_y"]["value"]
+#             xsig = dev.samcam_xsig.read()["samcam_xsig"]["value"]
+#             ysig = dev.samcam_ysig.read()["samcam_ysig"]["value"]
+#             f.write(f"{curr_bu}, {curr_bd}, {xpos}, {ypos}, {xsig}, {ysig}\n")
+#             print(f"{curr_bu}, {curr_bd}, {xpos}, {ypos}, {xsig}, {ysig}")
 
 
 def filters(pos="out"):
@@ -832,6 +796,11 @@ 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"""
@@ -867,111 +836,6 @@ def energy_vs_height(start_energy=11000, end_energy=25000, step=1000):
     # umv(dev.id_gap, 20.0)
 
 
-@dataclass
-class SlitParameters:
-    small: float
-    large: float
-    scanwidth: float
-    steps: any
-    x_size: any
-    x_cen: any
-    y_size: any
-    y_cen: any
-    mon: any
-
-
-# Generic scan function
-def scan_slits(name="s1", x=True, y=True):
-    # Define parameters for all slits
-    s1_params = SlitParameters(
-        small=0.075,
-        large=2.5,
-        scanwidth=2.5,
-        steps=40,
-        x_size=dev.s1_xsize,
-        y_size=dev.s1_ysize,
-        x_cen=dev.s1_xcen,
-        y_cen=dev.s1_ycen,
-        mon=dev.lu_bpmsum,
-    )
-    s2_params = SlitParameters(
-        small=0.075,
-        large=3.0,
-        scanwidth=2.5,
-        steps=50,
-        x_size=dev.s2_xsize,
-        y_size=dev.s2_ysize,
-        x_cen=dev.s2_xcen,
-        y_cen=dev.s2_ycen,
-        mon=dev.bcu_bpmsum,
-    )
-    s3_params = SlitParameters(
-        small=0.05,
-        large=0.5,
-        scanwidth=0.7,
-        steps=50,
-        x_size=dev.s3_xsize,
-        y_size=dev.s3_ysize,
-        x_cen=dev.s3_xcen,
-        y_cen=dev.s3_ycen,
-        mon=dev.i1,
-    )
-
-    # Map slit names to their parameters
-    slit_parameters = {"s1": s1_params, "s2": s2_params, "s3": s3_params}
-
-    if name not in slit_parameters:
-        raise ValueError(f"Unknown slit name: {name}")
-
-    params = slit_parameters[name]
-    print("Moving filters out")
-    filters("out")
-
-    print("Moving i1 diode into position")
-    SE.i1.mvin()
-
-    if x and y:
-        scan_in_x(
-            xsize=params.x_size,
-            xcen=params.x_cen,
-            mon=params.mon,
-            small=params.small,
-            large=params.large,
-            width=params.scanwidth,
-            steps=params.steps,
-        )
-        scan_in_y(
-            ysize=params.y_size,
-            ycen=params.y_cen,
-            mon=params.mon,
-            small=params.small,
-            large=params.large,
-            width=params.scanwidth,
-            steps=params.steps,
-        )
-    if x and not y:
-        scan_in_x(
-            xsize=params.x_size,
-            xcen=params.x_cen,
-            mon=params.mon,
-            small=params.small,
-            large=params.large,
-            width=params.scanwidth,
-            steps=params.steps,
-        )
-    else:
-        scan_in_y(
-            ysize=params.y_size,
-            ycen=params.y_cen,
-            mon=params.mon,
-            small=params.small,
-            large=params.large,
-            width=params.scanwidth,
-            steps=params.steps,
-        )
-
-    filters("in")
-
 
 def scan_in_x(xsize, xcen, mon, small, large, width, steps):
     umv(xsize, small)
@@ -1015,35 +879,43 @@ def calc_from_centre(cam="samcam"):
         return read, diff_px
 
 
-def get_cam():
-    umv(dev.scam_zoom, 100)
+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}")
-    umv(dev.scam_zoom, 500)
+    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}")
 
 
-def read_currents():
-    """Read the individual bpm channels"""
-    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"]
-    # print(f"bpm1 = {bpm1:.3f}, bpm2 = {bpm2:.3f}, bpm3 = {bpm3:.3f}, bpm4 = {bpm4:.3f}")
-    return {"bpm1": bpm1, "bpm2": bpm2, "bpm3": bpm3, "bpm4": bpm4}
+# def read_currents():
+#     """Read the individual bpm channels"""
+#     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"]
+#     # print(f"bpm1 = {bpm1:.3f}, bpm2 = {bpm2:.3f}, bpm3 = {bpm3:.3f}, bpm4 = {bpm4:.3f}")
+#     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
+# 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
+
+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'):
     
@@ -1471,4 +1343,184 @@ def get_mirror_data(history_index: int):
     #     "y_data": y_data,
     #     "motor_name": motor_name,
     #     "scan_number": scan_number,
-    # }
\ No newline at end of file
+    # }
+
+
+def record_beampos_until(stop_time_str, 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")
+    filename = f"luts/beampos/beam_pos_record_{fnow}.csv"
+    scam_zoom = dev.scam_zoom.read()['scam_zoom']['value']
+    umv(dev.scam_zoom, 1000)
+    beam_size = "25 x 25"
+
+    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)
+
+    # If that time already passed today, set it for tomorrow
+    if end_time <= now:
+        end_time += timedelta(days=1)
+
+    print(f"Starting now ({now:%Y-%m-%d %H:%M}), running until {end_time:%Y-%m-%d %H:%M}")
+    print(f"Interval: {interval_minutes} minutes\n")
+
+    while datetime.now() < end_time:
+        record_beam(filename)
+        next_run = datetime.now() + timedelta(minutes=interval_minutes)
+        if next_run > end_time:
+            break
+        print(f"Next run at {next_run:%H:%M:%S}")
+        time.sleep(interval_minutes * 60)
+
+    print(f"Finished — reached {end_time:%H:%M}.")
+
+
+
+def record_beam(filename):
+    """
+    Record the beam position to file.
+    """
+
+    # Perform the usual recording steps
+    umv(dev.ss_f4_x, 10)
+    auto_exposure(cam="samcam")
+    time.sleep(2)
+    now = datetime.now()
+
+    # Append data to the file
+    with open(filename, "a", encoding="utf-8") as f:
+        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}")
+        # xpos, ypos = compute_pos()
+        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"
+        )
+
+    umv(dev.ss_f4_x, 16)
+
+
+def analyse_beampos(file):
+    """Plot the recorded beam position data"""
+    # Load the CSV file into a Pandas DataFrame
+    filedir = "luts/beampos/"
+    filename = filedir + file
+    data = pd.read_csv(filename, skiprows=1)
+    print(data.head())
+
+    # ref_bpmx = 0.110
+    # ref_bmpy = -0.175
+    ref_scinx = 1024
+    ref_sciny = 980
+    scam_zoom = 116
+
+    # Extract the first column (time) as the x-axis
+    meas_time = data["Time"]
+
+    # subtract the reference positions from the measures position
+    # data["bpm_x"] = data["bpm_x"] - ref_bpmx
+    # data["bpm_y"] = data["bpm_y"] - ref_bmpy
+    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
+
+    fig, axs = plt.subplots(
+        2, 2, figsize=(12, 8), layout="constrained"
+    )  # Create a 2x2 grid for subplots
+    # 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"]
+
+    for idx, column in enumerate(y_positions.columns):
+        ax = axs[idx]  # Select the current subplot
+
+        # Determine the color based on the column name
+        color = "red" if column in red_columns else "blue"
+
+        # 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
+
+        # Set x-axis tick positions and labels
+        tick_positions = range(0, len(meas_time), 7)  # Adjust the step size as needed
+        ax.set_xticks(tick_positions)
+        ax.set_xticklabels([meas_time[i] for i in tick_positions], rotation=45)
+
+        # Add labels, title, legend, and grid to each subplot
+        ax.set_xlabel("Time")
+        ax.set_ylabel("beam displacement (microns)")
+        ax.set_title(f"Beam position at {column}")
+        ax.legend()
+        ax.grid(True)
+
+    # Show the resulting figure
+    plt.show()
+
+import time
+import numpy as np
+import matplotlib.pyplot as plt
+from epics import caget
+
+
+def scan_scin_focus():
+    z0 = dev.diag_z.position
+
+    z_positions = np.linspace(-4,4,30)
+    contrasts = []
+
+
+
+    for z in z_positions:
+        umv(dev.diag_z, z)
+        time.sleep(0.1)
+
+        sigma = dev.samcam_sigma.read()['samcam_sigma']['value']
+        mean = dev.samcam_mean.read()['samcam_mean']['value']
+
+        contrast = sigma / mean if mean > 0 else 0
+        contrasts.append(contrast)
+
+        print(f"z = {z:.4f}, sigma = {sigma:.3f}, mean = {mean:.3f}, contrast = {contrast:.5f}")
+
+    contrasts = np.asarray(contrasts)
+
+    best_idx = np.argmax(contrasts)
+    best_z = z_positions[best_idx]
+
+    print(f"\nBest scintillator focus z = {best_z:.4f}")
+
+    umv(dev.diag_z,best_z)
+
+    plt.plot(z_positions, contrasts, "o-")
+    plt.axvline(best_z, linestyle="--")
+    plt.xlabel("Scintillator z")
+    plt.ylabel("Sigma / Mean")
+    plt.show()
+
+    return best_z, z_positions, contrasts
\ No newline at end of file
diff --git a/pxii_bec/macros/mx_methods.py b/pxii_bec/macros/mx_methods.py
index 38cc6ed..9e012d9 100755
--- a/pxii_bec/macros/mx_methods.py
+++ b/pxii_bec/macros/mx_methods.py
@@ -105,6 +105,7 @@ def go_to_peak(
     settle: float = FitDefaults.SETTLE_TIME,
     confirm: bool = True,
     gomax: bool = False,
+    negative: bool = False
 ):
     """
     Go to the peak of a signal by scanning a motor within a specified range and
@@ -177,7 +178,10 @@ def go_to_peak(
     }
 
     # Define and fit model to scan data
-    fit_params = create_fit_parameters(False, FitDefaults.MODEL, FitDefaults.BASELINE)
+    fit_params = create_fit_parameters(deriv = False,
+                                       negative = False, 
+                                       model = FitDefaults.MODEL, 
+                                       baseline = FitDefaults.BASELINE)
     fit_result = fit(data, fit_params)
 
     # Plot the fitted data if plot = True
@@ -199,6 +203,7 @@ def fit_history(
     signal_name: str,
     deriv: bool = False,
     negative: bool = False,
+    smoothing: bool = False,
     model: str = FitDefaults.MODEL,
     move_to_peak: bool = False,
 ):