diff --git a/pxii_bec/macros/katscripts.py b/pxii_bec/macros/katscripts.py
index ff32868..cbf07a4 100755
--- a/pxii_bec/macros/katscripts.py
+++ b/pxii_bec/macros/katscripts.py
@@ -199,24 +199,134 @@ def scan_fe():
     close_fe_slits(False)
 
 
-def scan_gap(gomax = False, close = True, fine = False):
-    """Scan the gap and monitor on lu_bpmsum or bcu_bpmsum"""
-    close_fe_slits(close)
-    current_gap = dev.id_gap.read()["id_gap"]["value"]
-    if current_gap > 4.6:
-        go_to_peak(
-            dev.id_gap, dev.lu_bpmsum, -0.15, 0.15, 25, relative=True, gap=True, confirm=False, gomax=gomax
-        )
-    else:
-        go_to_peak(
-            dev.id_gap, dev.lu_bpmsum, 4.5, 4.75, 25, relative=False, gap=True, confirm=False, gomax=gomax
-        )
-    if fine:
-        go_to_peak(
-            dev.id_gap, dev.lu_bpmsum, -0.012, 0.012, 6, relative=True, gap=True, confirm=False, gomax=gomax
-        )
-    close_fe_slits(False)
+# def scan_gap(gomax = False, close = True, fine = False):
+#     """Scan the gap and monitor on lu_bpmsum or bcu_bpmsum"""
+#     close_fe_slits(close)
+#     current_gap = dev.id_gap.read()["id_gap"]["value"]
+#     if current_gap > 4.6:
+#         go_to_peak(
+#             dev.id_gap, dev.lu_bpmsum, -0.15, 0.15, 25, relative=True, gap=True, confirm=False, gomax=gomax
+#         )
+#     else:
+#         go_to_peak(
+#             dev.id_gap, dev.lu_bpmsum, 4.5, 4.75, 25, relative=False, gap=True, confirm=False, gomax=gomax
+#         )
+#     if fine:
+#         go_to_peak(
+#             dev.id_gap, dev.lu_bpmsum, -0.012, 0.012, 6, relative=True, gap=True, confirm=False, gomax=gomax
+#         )
+#     close_fe_slits(False)
 
+def scan_gap(gomax=False, close=True, fine=False, min_fwhm_steps=1.5):
+    """
+    Scan the undulator gap and monitor the LU BPM sum.
+
+    If gomax=True, go_to_peak() moves to the highest measured point
+    without applying peak-fit validation.
+
+    Returns
+    -------
+    dict
+        Results from the coarse and, if requested, fine scans.
+        `min_gap_reached` is True when the coarse-scan maximum is at
+        the physical minimum gap of 4.5 mm.
+    """
+
+    MIN_GAP = 4.5
+
+    close_fe_slits(close)
+
+    try:
+        current_gap = get_value(dev.id_gap)
+
+        # --------------------------------------------------------
+        # Coarse scan
+        # --------------------------------------------------------
+
+        if current_gap > 4.65:
+            coarse_result = go_to_peak(
+                dev.id_gap,
+                dev.lu_bpmsum,
+                -0.15,
+                0.15,
+                25,
+                relative=True,
+                gap=True,
+                confirm=False,
+                gomax=gomax,
+                min_fwhm_steps=min_fwhm_steps,
+            )
+
+        else:
+            coarse_result = go_to_peak(
+                dev.id_gap,
+                dev.lu_bpmsum,
+                MIN_GAP,
+                4.75,
+                25,
+                relative=False,
+                gap=True,
+                confirm=False,
+                gomax=gomax,
+                min_fwhm_steps=min_fwhm_steps,
+            )
+
+        # --------------------------------------------------------
+        # Check for physical minimum-gap limit
+        # --------------------------------------------------------
+
+        min_gap_reached = False
+
+        if gomax:
+            min_gap_reached = (
+                coarse_result["at_lower_boundary"]
+                and np.isclose(coarse_result["x_max"], MIN_GAP)
+            )
+
+            if min_gap_reached:
+                print(
+                    "\nMaximum occurs at the minimum gap of "
+                    f"{MIN_GAP:.2f} mm."
+                )
+                print(
+                    "The optimum for this harmonic may lie below "
+                    "the accessible gap."
+                )
+
+        # --------------------------------------------------------
+        # Fine scan
+        # --------------------------------------------------------
+
+        fine_result = None
+
+        if fine and not min_gap_reached:
+            fine_result = go_to_peak(
+                dev.id_gap,
+                dev.lu_bpmsum,
+                -0.012,
+                0.012,
+                6,
+                relative=True,
+                gap=True,
+                confirm=False,
+                gomax=gomax,
+            )
+
+        elif fine and min_gap_reached:
+            print(
+                "Skipping fine gap scan because the coarse maximum "
+                "is at the minimum gap."
+            )
+
+        return {
+            "success": True,
+            "coarse": coarse_result,
+            "fine": fine_result,
+            "min_gap_reached": min_gap_reached,
+        }
+
+    finally:
+        close_fe_slits(False)
 
 def record_pitch():
     """Record the pitch and perp distance to file"""
@@ -257,41 +367,310 @@ def beam_centre_from_ss():
     return beampos
 
 
-def generate_dcm_lut(env, planner, start_energy=6000, end_energy=30000, step=300, scint=True):
-    """Generate a lookup table for the dcm from start_energy to end_energy in steps of step"""
-    filename = env.lut_dir / f"dcm_calibration_{filenow()}.csv"
-    zoom = dev.scam_zoom.read()['scam_zoom']['value']
-    with open(filename, "w", encoding="utf-8") as f:
-        f.write(f"energy,harmonic,gap,pitch,samcam_x,samcam_y,zoom={zoom}\n")
+def generate_dcm_lut(
+    start_energy=6000,
+    end_energy=30000,
+    step=300,
+    scan_id_gap=True,
+    scint=True,
+):
+    """Generate a DCM calibration lookup table.
+
+    Parameters
+    ----------
+    start_energy : int
+        Starting energy in eV.
+    end_energy : int
+        Final energy in eV.
+    step : int
+        Energy step in eV.
+    scan_id_gap : bool
+        If True, optimise both undulator gap and DCM pitch.
+        If False, only optimise the DCM pitch.
+    scint : bool
+        If True, record sample-camera beam position.
+
+    Returns
+    -------
+    list
+        List of failed scans with energy, scan type and failure reason.
+    """
+
+    filename = bl.env.lut_dir / f"dcm_calibration_{filenow(date_only=True)}.csv"
+    failure_filename = (
+        bl.env.lut_dir / f"dcm_calibration_failures_{filenow(date_only=True)}.csv"
+    )
+
+    failures = []
+
+    zoom = get_value(dev.scam_zoom)
+
+    # Create calibration file
+    if scan_id_gap:
+        with open(filename, "w", encoding="utf-8") as f:
+            if scint:
+                f.write(
+                    f"energy,harmonic,gap,pitch,"
+                    f"samcam_x,samcam_y,zoom={zoom}\n"
+                )
+            else:
+                f.write(
+                    "energy,harmonic,gap,pitch\n"
+                )
+    else:
+        with open(filename, "w", encoding="utf-8") as f:
+            if scint:
+                f.write(
+                    f"energy,pitch,samcam_x,samcam_y,zoom={zoom}\n"
+                )
+            else:
+                f.write(
+                    "energy,pitch\n"
+                )
+
+    # Create failure log
+    with open(failure_filename, "w", encoding="utf-8") as f:
+        f.write("energy,scan,reason\n")
+    try:
+        for energy in range(start_energy, end_energy+1, step):
+
+            print("\n" + "=" * 60)
+            print(f"Moving to {energy:.0f} eV")
+            print("=" * 60)
+
+            # ------------------------------------------------------------
+            # Move to energy without automatically running mono pitch scan
+            # ------------------------------------------------------------
+
+            try:
+                bl_energy(energy, mono_scan=False)
+
+            except Exception as exc:
+                reason = str(exc)
+
+                print(f"Energy move failed at {energy:.0f} eV: {reason}")
+
+                failure = {
+                    "energy": energy,
+                    "scan": "energy",
+                    "reason": reason,
+                }
+
+                failures.append(failure)
+
+                with open(failure_filename, "a", encoding="utf-8") as f:
+                    f.write(
+                        f"{energy:.0f},energy,{reason}\n"
+                    )
+
+                # Cannot calibrate this energy if the energy move failed
+                continue
+
+            # ------------------------------------------------------------
+            # Undulator gap scan
+            # ------------------------------------------------------------
+
+            if scan_id_gap:
+
+                try:
+                    scan_gap(
+                        close=True,
+                        gomax=True,
+                        fine=True,
+                    )
+
+                except Exception as exc:
+                    reason = str(exc)
+
+                    print(
+                        f"Gap scan failed at {energy:.0f} eV: "
+                        f"{reason}"
+                    )
+
+                    failure = {
+                        "energy": energy,
+                        "scan": "gap",
+                        "reason": reason,
+                    }
+
+                    failures.append(failure)
+
+                    with open(
+                        failure_filename,
+                        "a",
+                        encoding="utf-8",
+                    ) as f:
+                        f.write(
+                            f"{energy:.0f},gap,{reason}\n"
+                        )
+
+                    # Fall back to calculated gap
+                    print(
+                        f"Setting calculated gap for "
+                        f"{energy:.0f} eV"
+                    )
+                    set_gap(energy)
+
+            # ------------------------------------------------------------
+            # DCM pitch scan
+            # ------------------------------------------------------------
+
+            try:
+                result = mono_pitch_scan()
+
+                if not result["success"]:
+
+                    reason = "; ".join(result["problems"])
+
+                    print(
+                        f"Mono pitch scan rejected at "
+                        f"{energy:.0f} eV:"
+                    )
+
+                    for problem in result["problems"]:
+                        print(f"  - {problem}")
+
+                    failure = {
+                        "energy": energy,
+                        "scan": "mono_pitch",
+                        "reason": reason,
+                    }
+
+                    failures.append(failure)
+
+                    with open(
+                        failure_filename,
+                        "a",
+                        encoding="utf-8",
+                    ) as f:
+                        f.write(
+                            f"{energy:.0f},mono_pitch,{reason}\n"
+                        )
+
+            except Exception as exc:
+                reason = str(exc)
+
+                print(
+                    f"Mono pitch scan failed at "
+                    f"{energy:.0f} eV: {reason}"
+                )
+
+                failure = {
+                    "energy": energy,
+                    "scan": "mono_pitch_exception",
+                    "reason": reason,
+                }
+
+                failures.append(failure)
+
+                with open(
+                    failure_filename,
+                    "a",
+                    encoding="utf-8",
+                ) as f:
+                    f.write(
+                        f"{energy:.0f},mono_pitch_exception,"
+                        f"{reason}\n"
+                    )
+
+            # ------------------------------------------------------------
+            # Record calibration values
+            # ------------------------------------------------------------
+
+            pitch = get_value(dev.dcm_pitch)
+
+            if scan_id_gap:
+                gap_value = get_value(dev.id_gap)
+                harm = Gap()
+                h = harm.get_harmonic_by_energy(energy)
+
+            if scint:
+                x, y = read_scintXY()
+
+            with open(filename, "a", encoding="utf-8") as f:
+
+                if scan_id_gap:
+
+                    if scint:
+                        f.write(
+                            f"{energy:.0f},"
+                            f"{h},"
+                            f"{gap_value:.5g},"
+                            f"{pitch:.5g},"
+                            f"{x:.4g},"
+                            f"{y:.4g}\n"
+                        )
+                    else:
+                        f.write(
+                            f"{energy:.0f},"
+                            f"{h},"
+                            f"{gap_value:.5g},"
+                            f"{pitch:.5g}\n"
+                        )
+
+                else:
+
+                    if scint:
+                        f.write(
+                            f"{energy:.0f},"
+                            f"{pitch:.5g},"
+                            f"{x:.4g},"
+                            f"{y:.4g}\n"
+                        )
+                    else:
+                        f.write(
+                            f"{energy:.0f},"
+                            f"{pitch:.5g}\n"
+                        )
+
+    finally:
+        print("\nReturning beamline to safe conditions...")
 
-    for energy in range(start_energy, end_energy, step):
-        print(f"Moving to {energy: .0f}eV")
-        bl_energy(energy, mono_scan=False)
         try:
-            scan_gap(close=True, gomax=True, fine=True)
+            bl_energy(12400)
         except Exception as exc:
-            print(f"Gap scan failed at {energy:.0f} eV: {exc}")
-            set_gap(energy)  # set fallback gap here
+            print(f"WARNING: Could not return to 12400 eV: {exc}")
 
-        mono_pitch_scan()
-        gap = dev.id_gap.read()["id_gap"]["value"]
-        harm = Gap()
-        h = harm.get_harmonic_by_energy(energy)
-        pitch = dev.dcm_pitch.read()["dcm_pitch"]["value"]
-        if scint:
-            x, y = read_samcam_scint(planner)
-            with open(filename, "a", encoding="utf-8") as f:
-                f.write(
-                    f"{energy:.0f},{h},{gap:.5g},{pitch:.5g},{x:.4g},{y:.4g}\n"
-                )
-        else:
-            with open(filename, "a", encoding="utf-8") as f:
-                f.write(
-                    f"{energy:.0f},{h},{gap:.5g},{pitch:.5g}\n"
-                )
+        try:
+            umv(dev.id_gap, 20)
+        except Exception as exc:
+            print(f"WARNING: Could not open undulator gap: {exc}")
+
+
+    # ------------------------------------------------------------
+    # Summary
+    # ------------------------------------------------------------
+
+    print("\n" + "=" * 60)
+    print("DCM calibration complete")
+    print("=" * 60)
+
+    if failures:
+
+        print(f"\n{len(failures)} scan failure(s):\n")
+
+        for failure in failures:
+            print(
+                f"{failure['energy']:5.0f} eV  "
+                f"{failure['scan']:22s}  "
+                f"{failure['reason']}"
+            )
+
+        print(
+            f"\nFailure log written to:\n"
+            f"{failure_filename}"
+        )
+
+    else:
+        print("\nNo scan failures.")
+
+    print(
+        f"\nCalibration data written to:\n"
+        f"{filename}"
+    )
+
+    return failures
 
-    bl_energy(12400)
-    umv(dev.id_gap, 20.0)
 
 
 def write_lut():
@@ -631,15 +1010,21 @@ def correlate_fpitch():
     plt.show()
 
 
-def correlate_hfm_x():
+def correlate_hfm_xd(meas="bcu"):
     x_correlation = []
-    lat = dev.hfm_x.read()["hfm_x"]["value"]
-    scanpoints = np.linspace(lat - 0.1, lat + 0.1, 10)
+    startx = get_value(dev.hfm_xd)
+    scanpoints = np.linspace(startx - 0.02, startx + 0.02, 10)
     for pos in scanpoints:
-        umv(dev.hfm_x, pos)
-        xpos, ypos = compute_pos()
+        umv(dev.hfm_xd, pos)
+        if meas == "bcu":
+            xpos, ypos = get_bcu_beampos()
+        elif meas == "scint":
+            xpos, ypos = read_scintXY(1000)
+        else:
+            print("define measurement device, bcu or scint")
+            return
         x_correlation.append(xpos)
-    umv(dev.hfm_x, lat)
+    umv(dev.hfm_xd, startx)
 
     x = np.array(scanpoints)
     y1 = np.array(x_correlation)
@@ -656,45 +1041,61 @@ def correlate_hfm_x():
     plt.plot(
         x_fit, xslope * x_fit + xint, label=f"y = {xslope:.5g}*x + ({xint:.5g})", color="lightcoral"
     )
-    plt.xlabel("HFM lateral (mm)")
-    plt.ylabel("beam X (from BCU BPM)")
+    plt.xlabel("HFM Downstream X (mm)")
+    if meas == "bcu":
+        plt.ylabel("beam X (from BCU BPM)")
+    elif meas == "scint":
+      plt.ylabel("beam X (from  Scintillator)")  
     plt.legend()
     plt.grid(True)
     plt.tight_layout()
     plt.show()
 
+def correlate_vfm_ydr(meas="bcu"):
+    y_correlation = []
+    starty = get_value(dev.vfm_ydr)
+    if meas == "bcu":
+        scanpoints = np.linspace(starty - 0.005, starty + 0.005, 10)
+    elif meas == "scint":
+        scanpoints = np.linspace(starty - 0.02, starty + 0.02, 10)
+    
+    for pos in scanpoints:
+        umv(dev.vfm_ydr, pos)
+        if meas == "bcu":
+            xpos, ypos = get_bcu_beampos()
+        elif meas == "scint":
+            xpos, ypos = read_scintXY(1000)
+        else:
+            print("define measurement device, bcu or scint")
+            return
+        y_correlation.append(ypos)
+    umv(dev.vfm_ydr, starty)
 
-# def correlate_vfm_y():
-#     y_correlation = []
-#     vert = dev.vfm_y.read()["vfm_y"]["value"]
-#     scanpoints = np.linspace(vert - 0.06, vert + 0.06, 10)
-#     for pos in scanpoints:
-#         umv(dev.vfm_y, pos)
-#         xpos, ypos = compute_pos()
-#         y_correlation.append(ypos)
-#     umv(dev.vfm_y, vert)
+    x = np.array(scanpoints)
+    y1 = np.array(y_correlation)
 
-#     x = np.array(scanpoints)
-#     y1 = np.array(y_correlation)
+    xslope, xint = np.polyfit(x, y1, 1)
 
-#     yslope, yint = np.polyfit(x, y1, 1)
+    print(f"Fit for ypos is: y = {xslope:.5g}.x + {xint:.5g}")
+    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, xslope * x_fit + xint, label=f"y = {xslope:.5g}*x + ({xint:.5g})", color="lightcoral"
+    )
+    plt.xlabel("VFM Downstream Ring Y (mm)")
+    if meas == "bcu":
+        plt.ylabel("beam Y (from BCU BPM)")
+    if meas == "scint":
+        plt.ylabel("beam Y (from Scintillator)")
+    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():
@@ -1080,15 +1481,20 @@ def read_scintXY(zoom=800):
     a sepcified zoom. Default zoom is 800"""
     # Ensure shutter and detector cover are closed before changing states
     dev.bcu_shutter.put(0)
-    bl.d['det_cov'].move('close')
+    if bl.d['det_cov'].pos != 'close':
+        bl.d['det_cov'].move('close')
     try:
         # set transmission to 10%
-        dev.transm.put(0.1)
+        current_transm = get_value(dev.transm)
+        if abs(current_transm - 0.1) > 0.01:
+            dev.transm.put(0.1)
         # set camera zoom
-        umv(dev.scam_zoom, zoom)
+        if get_value(dev.scam_zoom) != zoom:
+           umv(dev.scam_zoom, zoom)
         # move to beam visualisation
 
-        bl.planner.move_to(BeamlineState.BEAM_VISUALISATION)
+        if not bl.planner.is_state(BeamlineState.BEAM_VISUALISATION):
+            bl.planner.move_to(BeamlineState.BEAM_VISUALISATION)
         # open shutter
         dev.bcu_shutter.put(1)
         # set autoexposure
@@ -1100,78 +1506,7 @@ def read_scintXY(zoom=800):
     finally:
         dev.bcu_shutter.put(0)
 
-def find_best_roll(low = 6000, high = 20000):
-    # low_scanpoints = np.linspace(2.0, 7.0, 11)
-    low_scanpoints = np.linspace(2.5, 6.5, 9)
-    high_scanpoints = np.linspace(1.0, 9.0,9)
-    highs_y, lows_y = [], []
-    filedir, _ = filenow(dir="test")
 
-    # low energy
-    bl_energy(low)
-    # filename = f"{filedir}/roll_{low:.0f}.csv"
-    filename = env.lut_dir / f"roll_{low:.0f}.csv"
-    with open(filename, "w", encoding="utf-8") as f:
-        f.write("Roll, Y pos\n") 
-    for i in low_scanpoints:
-        umv(dev.dcm_froll, i)
-        time.sleep(0.2)
-        _, low_y = read_samcam_scint()
-        print(f"Ypos at sample cam is {low_y}")
-        lows_y.append(low_y)
-        with open(filename, "a", encoding="utf-8") as f:
-            f.write(f"{i:.4g},{low_y:.4g} \n") 
-
-    # high energy
-    bl_energy(high)
-    # filename = f"{filedir}/roll_{high:.0f}.csv"
-    filename = env_lut_dir / f"roll_{high:.0f}.csv"
-    with open(filename, "w", encoding="utf-8") as f:
-        f.write("Roll, Y pos\n") 
-    auto_exposure("samcam", max_iter=25)
-    for i in high_scanpoints:
-        umv(dev.dcm_froll, i)
-        time.sleep(0.2)
-        _, high_y = read_samcam_scint()
-        print(f"Ypos at sample cam is {high_y}")
-        highs_y.append(high_y)
-        with open(filename, "a", encoding="utf-8") as f:
-            f.write(f"{i:.4g},{high_y:.4g} \n") 
-    
-    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)
-
-    # # 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"  DCM Roll = {x_intersect:.5f}")
-    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')
-
-    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.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.legend()
-    plt.grid(True)
-    plt.tight_layout()
-    plt.show()
 
 def vfm_yaw_check():
     slit_size = 0.02
@@ -1739,3 +2074,417 @@ def read_i2():
 
 def quick_test():
     bl.planner.closest_states()
+
+import time
+import numpy as np
+
+
+def _read_signal(signal):
+    """Return the value from a BEC/ophyd signal."""
+    data = signal.read()
+    key = next(iter(data))
+    return data[key]["value"]
+
+
+def _smooth_profile(profile, window=7):
+    """Simple moving-average smoothing."""
+    profile = np.asarray(profile, dtype=float)
+
+    if window <= 1:
+        return profile
+
+    kernel = np.ones(window) / window
+    return np.convolve(profile, kernel, mode="same")
+
+
+def find_dark_feature(
+    profile,
+    search_min,
+    search_max,
+    smooth_window=7,
+    threshold_fraction=0.35,
+):
+    """
+    Find the centre and edges of a dark feature in a 1-D intensity profile.
+
+    Parameters
+    ----------
+    profile : array-like
+        Intensity profile.
+
+    search_min, search_max : int
+        Pixel range in which to search.
+
+    smooth_window : int
+        Width of moving-average smoothing.
+
+    threshold_fraction : float
+        Fraction of the maximum dip depth used to define the feature edges.
+
+    Returns
+    -------
+    dict
+        centre, low_edge, high_edge, width, depth
+    """
+
+    profile = np.asarray(profile, dtype=float)
+    smooth = _smooth_profile(profile, smooth_window)
+
+    x = np.arange(len(smooth))
+
+    mask = (x >= search_min) & (x <= search_max)
+    xx = x[mask]
+    yy = smooth[mask]
+
+    if len(yy) == 0:
+        raise ValueError("Search range contains no profile points")
+
+    # Estimate local background from the brighter part of the profile.
+    background = np.percentile(yy, 80)
+
+    # How far each pixel lies below the background.
+    darkness = np.clip(background - yy, 0, None)
+
+    peak_index = np.argmax(darkness)
+    peak_depth = darkness[peak_index]
+
+    if peak_depth <= 0:
+        raise RuntimeError("No dark feature found")
+
+    threshold = threshold_fraction * peak_depth
+
+    # Start at deepest point and walk outwards until the dip disappears.
+    left = peak_index
+    while left > 0 and darkness[left] > threshold:
+        left -= 1
+
+    right = peak_index
+    while right < len(darkness) - 1 and darkness[right] > threshold:
+        right += 1
+
+    low_edge = xx[left]
+    high_edge = xx[right]
+
+    # Darkness-weighted centre is more stable than a single minimum pixel.
+    feature_mask = np.arange(len(xx))
+    feature_mask = (feature_mask >= left) & (feature_mask <= right)
+
+    weights = darkness[feature_mask]
+    coords = xx[feature_mask]
+
+    centre = np.average(coords, weights=weights)
+
+    return {
+        "centre": float(centre),
+        "low_edge": float(low_edge),
+        "high_edge": float(high_edge),
+        "width": float(high_edge - low_edge),
+        "depth": float(peak_depth),
+    }
+
+
+def measure_omega_runout(
+    omega_motor,
+    profile_x_signal,
+    profile_y_signal,
+    cursor_x_signal,
+    *,
+    angles=None,
+    x_search=(900, 1200),
+    y_search=(850, 1150),
+    settle=1.0,
+    cursor_settle=0.2,
+    smooth_window=7,
+    threshold_fraction=0.35,
+):
+    """
+    Measure apparent needle position while rotating omega.
+
+    CursorY should already be positioned so that the X profile crosses
+    the needle.
+
+    At each omega:
+        1. Read ProfileCursorX
+        2. Find needle X
+        3. Move CursorX to that position
+        4. Read ProfileCursorY
+        5. Find needle Y
+
+    Returns
+    -------
+    list[dict]
+        One result dictionary per omega position.
+    """
+
+    if angles is None:
+        angles = np.arange(0, 361, 45)
+
+    results = []
+
+    for angle in angles:
+        print(f"\nOmega = {angle:.1f}°")
+
+        umv(omega_motor,(float(angle)))
+        time.sleep(settle)
+
+        # ------------------------------------------------------
+        # X position
+        # ------------------------------------------------------
+        profile_x = _read_signal(profile_x_signal)
+
+        x_result = find_dark_feature(
+            profile_x,
+            search_min=x_search[0],
+            search_max=x_search[1],
+            smooth_window=smooth_window,
+            threshold_fraction=threshold_fraction,
+        )
+
+        needle_x = x_result["centre"]
+
+        # Put vertical cursor through the detected needle position.
+        cursor_x_signal.put(int(round(needle_x)))
+
+        time.sleep(cursor_settle)
+
+        # ------------------------------------------------------
+        # Y position
+        # ------------------------------------------------------
+        profile_y = _read_signal(profile_y_signal)
+
+        y_result = find_dark_feature(
+            profile_y,
+            search_min=y_search[0],
+            search_max=y_search[1],
+            smooth_window=smooth_window,
+            threshold_fraction=threshold_fraction,
+        )
+
+        needle_y = y_result["centre"]
+
+        result = {
+            "omega": float(angle),
+
+            "x": needle_x,
+            "x_width": x_result["width"],
+            "x_depth": x_result["depth"],
+
+            "y": needle_y,
+            "y_width": y_result["width"],
+            "y_depth": y_result["depth"],
+        }
+
+        results.append(result)
+
+        print(
+            f"  X = {needle_x:8.2f} px"
+            f"    width = {x_result['width']:6.1f}"
+        )
+        print(
+            f"  Y = {needle_y:8.2f} px"
+            f"    width = {y_result['width']:6.1f}"
+        )
+
+    # ----------------------------------------------------------
+    # Summary
+    # ----------------------------------------------------------
+
+    xs = np.asarray([r["x"] for r in results])
+    ys = np.asarray([r["y"] for r in results])
+
+    print("\nOmega runout")
+    print("------------")
+    print(f"X peak-to-peak : {np.ptp(xs):.2f} px")
+    print(f"Y peak-to-peak : {np.ptp(ys):.2f} px")
+
+    if len(results) > 1:
+        return_error = np.hypot(
+            xs[-1] - xs[0],
+            ys[-1] - ys[0],
+        )
+
+        print(
+            f"0° -> 360° error: {return_error:.2f} px"
+        )
+
+    return results
+
+
+def dcm_pitch_calibration(start_energy=6000, end_energy=31000, step=3000):
+    filename = bl.env.lut_dir / f"dcm_pitch_{filenow(date_only=True)}.csv"
+    with open(filename, "w", encoding="utf-8") as f:
+        f.write("energy,pitch\n")
+    for energy in range(start_energy, end_energy, step):
+    # energy_list = [29000,26000,23000,20000,17000,14000,11000,8000,6000]
+    # for energy in energy_list:
+        bl_energy(energy, mono_scan=False)
+        result = go_to_peak(dev.dcm_pitch, 
+                            dev.lu_bpmsum, 
+                            -0.1,
+                            0.4, 
+                            40, 
+                            min_fwhm_steps=1, 
+                            confirm=False
+                            )
+        pitch_new = get_value(dev.dcm_pitch)
+        with open(filename, "a", encoding="utf-8") as f:
+            f.write(f"{energy:.0f},{pitch_new:.5g}\n")
+        
+def scan_roll(
+    energy,
+    start,
+    stop,
+    steps=9,
+    settle=0.2,
+    overwrite=False,
+):
+    """Measure scintillator Y position as a function of DCM fine roll."""
+
+    scanpoints = np.linspace(start, stop, steps)
+    filename = bl.env.lut_dir / f"roll_{energy:.0f}.csv"
+
+    if filename.exists() and not overwrite:
+        print(f"{filename.name} already exists.")
+        print("Use overwrite=True if you want to repeat this energy.")
+        return filename
+
+    bl_energy(energy)
+
+    with open(filename, "w", encoding="utf-8") as f:
+        f.write("Roll,Y pos\n")
+
+        for roll in scanpoints:
+            umv(dev.dcm_froll, roll)
+            time.sleep(settle)
+
+            _, ypos = read_scintXY()
+
+            print(
+                f"{energy / 1000:.1f} keV | "
+                f"roll = {roll:.3f} | "
+                f"Y = {ypos:.2f}"
+            )
+
+            f.write(f"{roll:.6f},{ypos:.4f}\n")
+
+    return filename
+
+def analyse_roll(low=6000, high=25000):
+    """Find the DCM fine-roll value where low/high-energy beam positions intersect."""
+
+    low_file = bl.env.lut_dir / f"roll_{low:.0f}.csv"
+    high_file = bl.env.lut_dir / f"roll_{high:.0f}.csv"
+
+    low_data = np.loadtxt(
+        low_file,
+        delimiter=",",
+        skiprows=1,
+    )
+
+    high_data = np.loadtxt(
+        high_file,
+        delimiter=",",
+        skiprows=1,
+    )
+
+    low_roll = low_data[:, 0]
+    low_y = low_data[:, 1]
+
+    high_roll = high_data[:, 0]
+    high_y = high_data[:, 1]
+
+    # Linear fits: y = mx + b
+    m_low, b_low = np.polyfit(low_roll, low_y, 1)
+    m_high, b_high = np.polyfit(high_roll, high_y, 1)
+
+    # Intersection
+    x_intersect = (b_low - b_high) / (m_high - m_low)
+    y_intersect = m_low * x_intersect + b_low
+
+    print("\nIntersection:")
+    print(f"  DCM Roll = {x_intersect:.5f}")
+    print(f"  samcam Y = {y_intersect:.2f}")
+
+    # Plot
+    plt.figure(figsize=(7, 5))
+
+    plt.scatter(
+        low_roll,
+        low_y,
+        label=f"{low / 1000:.1f} keV",
+    )
+
+    plt.scatter(
+        high_roll,
+        high_y,
+        label=f"{high / 1000:.1f} keV",
+    )
+
+    x_min = min(low_roll.min(), high_roll.min())
+    x_max = max(low_roll.max(), high_roll.max())
+
+    x_fit = np.linspace(x_min, x_max, 200)
+
+    plt.plot(
+        x_fit,
+        m_low * x_fit + b_low,
+        "--",
+        label=f"{low / 1000:.1f} keV fit",
+    )
+
+    plt.plot(
+        x_fit,
+        m_high * x_fit + b_high,
+        "--",
+        label=f"{high / 1000:.1f} keV fit",
+    )
+
+    plt.scatter(
+        x_intersect,
+        y_intersect,
+        s=80,
+        zorder=5,
+        label="Intersection",
+    )
+
+    plt.xlabel("DCM fine roll")
+    plt.ylabel("samcam Y pos")
+    plt.title("DCM Roll Calibration")
+    plt.legend()
+    plt.grid()
+    plt.tight_layout()
+    plt.show()
+
+    return x_intersect
+
+def find_best_roll(
+    low=6000,
+    high=25000,
+    low_range=(4.0, 7.0),
+    high_range=(1.0, 8.0),
+    steps=9,
+):
+    """Collect missing roll scans and calculate optimum roll."""
+
+    low_file = bl.env.lut_dir / f"roll_{low:.0f}.csv"
+    high_file = bl.env.lut_dir / f"roll_{high:.0f}.csv"
+
+    if not low_file.exists():
+        scan_roll(
+            low,
+            *low_range,
+            steps=steps,
+        )
+    else:
+        print(f"Using existing {low_file.name}")
+
+    if not high_file.exists():
+        scan_roll(
+            high,
+            *high_range,
+            steps=steps,
+        )
+    else:
+        print(f"Using existing {high_file.name}")
+
+    return analyse_roll(low, high)
\ No newline at end of file