407 add automatic data analysis com gaussian fit 1st2nd derivatives after scans complete

docs/tutorials/3_blueapi_session.rst

@@ -57,6 +57,9 @@ Once the interactive IPython shell is running, you have access to:
 * ``scan_data``: A dictionary containing all collected event data from your scans.
 * ``plot(scan_id=None, cmap="viridis")``: Opens an interactive visualisation window.
 * ``run(script)``: Executes a Python script inside the active session namespace (a shortcut for ``%run -i <script>``).
+* ``fit_scan(scan_id=None, x=None, y=None, guess=None, log_df=None)``: Performs a 1D Gaussian fit, centre-of-mass (COM) calculation, and derivative on scan data.
+* ``set_fit_axes(x=None, y=None, log_df=False)``: Sets preferred X and Y channel names (and optional log derivative flag) for ``fit_scan()`` for post-scan analysis.
+* ``last_fit``: Stores the latest fitted summary parameters (``centre``, ``com``, ``max_slope_x``, ``min_slope_x``, ``zero_crossing``) for interactive use and scripting.
 
 Writing & Running Scan Scripts
 ------------------------------
@@ -81,6 +84,10 @@ Create a ``.py`` file (for example, ``my_scans.py``) containing the plans and Py
     # Stack a step scan afterwards
     pl.step_scan([dev.idd_gap], [[dev.pgm.energy, [200, 300, 10]]])
 
+    # Access the auto-fitted parameters from the last scan
+    print(f"Peak centre: {last_fit.centre}, COM: {last_fit.com}, Zero crossings: {last_fit.zero_crossing}")
+    print(f"Slopes: max={last_fit.max_slope_x}, min={last_fit.min_slope_x}")
+
 Running the Script
 ^^^^^^^^^^^^^^^^^^
 
@@ -124,6 +131,36 @@ Live data is stored in the ``scan_data`` dictionary, which is keyed first by the
     last_scan_id = next(reversed(scan_data))
     print(scan_data[last_scan_id])
 
+Scan Fitting & Auto-Analysis
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+By default (``auto_analyse=True``), ``BlueAPISession`` automatically analyses each scan when it completes:
+
+* Fits a 1D Gaussian plus linear background to the scan data.
+* Calculates the centre of mass (COM) and 1st/2nd numerical derivatives (``d1``, ``d2``).
+* Appends the derived curves (``gaussian``, ``d1``, ``d2``) and scalar metrics (``centre``, ``com``, ``zero_crossing``, ``max_slope_x``, ``min_slope_x``) to ``scan_data[scan_id]``.
+* Updates ``last_fit`` with ``centre``, ``com``, ``max_slope_x``, ``min_slope_x``, and ``zero_crossing``.
+
+You can also configure preferred fit axes or run fits manually on any scan:
+
+.. code-block:: python
+
+    # Set preferred X and Y channels for auto-analysis and fit_scan()
+    set_fit_axes(x="pgm-energy", y="det1")
+
+    # Fit the most recent scan (or pass a specific scan_id, x, y, or initial guess)
+    res = fit_scan()
+    res = fit_scan(scan_id="30", x="pgm-energy", y="det1")
+
+    # Access the latest fitted summary parameters via last_fit
+    print(
+        last_fit.centre,
+        last_fit.com,
+        last_fit.max_slope_x,
+        last_fit.min_slope_x,
+        last_fit.zero_crossing,
+    )
+
 Interactive Plotting
 ^^^^^^^^^^^^^^^^^^^^
 

pyproject.toml

@@ -32,7 +32,7 @@ requires-python = ">=3.11"
 
 [project.optional-dependencies]
 testing = ["jupyter", "pyepics"]
-client = ["ipython", "matplotlib", "blueapi"]
+client = ["ipython", "matplotlib", "blueapi", "lmfit"]
 server = ["pyserial", "zhinst-core", "blueapi"]
 [dependency-groups]
 dev = [

src/sm_bluesky/common/clients/blueapi_session.py

@@ -3,21 +3,47 @@
 import time
 from collections.abc import Callable
 from concurrent.futures import ThreadPoolExecutor
+from dataclasses import dataclass
 from datetime import datetime
 from functools import wraps
 from pathlib import Path
 from typing import Any
 
 import click
 import matplotlib.pyplot as plt
+import numpy as np
 from blueapi.client import BlueapiClient
 from blueapi.client.event_bus import AnyEvent
 from blueapi.config import ApplicationConfig
 from blueapi.core import DataEvent
+from lmfit.model import ModelResult
+
+from sm_bluesky.common.math import (
+    DerivativeResult,
+    calculate_com,
+    calculate_derivatives,
+    fit_gaussian,
+)
 
 from .interactive_plot import DEFAULT_CMAP, InteractivePlotWindow
 
 
+@dataclass
+class FitResult:
+    gaussian: ModelResult
+    com: float
+    df: DerivativeResult
+
+
+@dataclass
+class LastFit:
+    centre: float | None = None
+    com: float | None = None
+    max_slope_x: float | None = None
+    min_slope_x: float | None = None
+    zero_crossing: list[float] | None = None
+
+
 class BlueAPISession:
     def __init__(
         self, config: ApplicationConfig, instrument_session: str | None = None
@@ -26,6 +52,11 @@ def __init__(
         self.data: dict[str, dict[str, list[Any]]] = {}
         self.current_scan_id: Any | None = None
         self._active_plot_windows: list[InteractivePlotWindow] = []
+        self.auto_analyse: bool = True
+        self.preferred_x: str | None = None
+        self.preferred_y: str | None = None
+        self.log_df: bool = False
+        self.last_fit: LastFit = LastFit()
         click.echo("Connecting to BlueAPI service...")
         self.bc = BlueapiClient.from_config(self.config)
         self.bc.login()
@@ -114,6 +145,9 @@ def start_shell(self) -> None:
                 "scan_data": self.data,
                 "plot": self.plot,
                 "run": self.run,
+                "fit_scan": self.fit_scan,
+                "set_fit_axes": self.set_fit_axes,
+                "last_fit": self.last_fit,
             },
         )
 
@@ -160,6 +194,56 @@ def plot(self, scan_id: Any | None = None, cmap: str = DEFAULT_CMAP) -> None:
         window = InteractivePlotWindow(self.data, scan_id, cmap=cmap)
         self._active_plot_windows.append(window)
 
+    def set_fit_axes(
+        self,
+        x: str | None = None,
+        y: str | None = None,
+        log_df: bool = False,
+    ) -> None:
+        """Set preferred X/Y channels and log-derivative flag for fit_scan."""
+        self.preferred_x = x
+        self.preferred_y = y
+        self.log_df = log_df
+
+    def fit_scan(
+        self,
+        scan_id: Any | None = None,
+        x: str | None = None,
+        y: str | None = None,
+        guess: Any | None = None,
+        log_df: bool | None = None,
+    ) -> FitResult:
+        """Perform Gaussian fit and Centre of Mass calculation on scan data."""
+        if scan_id is None:
+            scan_id = next(reversed(self.data))
+        scan = self.data[str(scan_id)]
+        x_col = x or (self.preferred_x if self.preferred_x in scan else None)
+        y_col = y or (self.preferred_y if self.preferred_y in scan else None)
+        vals = list(scan.values())
+        if len(vals) >= 2:
+            x_data = scan[x_col] if x_col else vals[0]
+            y_data = scan[y_col] if y_col else vals[1]
+        elif len(vals) == 1:
+            y_data = scan[y_col] if y_col else vals[0]
+            x_data = scan[x_col] if x_col else list(range(1, len(y_data) + 1))
+        else:
+            raise ValueError("Scan has no data")
+        use_log = self.log_df if log_df is None else log_df
+        g_res = fit_gaussian(x_data, y_data, guess)
+        com_res = calculate_com(x_data, y_data)
+        df_res = calculate_derivatives(x_data, y_data, log=use_log)
+        click.echo(
+            f"Fit (scan {scan_id}): centre={g_res.best_values['center']:.4f}, "
+            f"FWHM={g_res.params['fwhm'].value:.4f}, COM={com_res:.4f}, "
+            f"zero_crossing={df_res.zero_crossings}"
+        )
+        self.last_fit.centre = float(g_res.best_values["center"])
+        self.last_fit.com = com_res
+        self.last_fit.max_slope_x = df_res.max_slope_x
+        self.last_fit.min_slope_x = df_res.min_slope_x
+        self.last_fit.zero_crossing = df_res.zero_crossings
+        return FitResult(gaussian=g_res, com=com_res, df=df_res)
+
     def print_plans(self) -> None:
         click.echo("\nPlans available:")
         for plan in self.bc.plans:
@@ -239,6 +323,32 @@ def feedback(event: AnyEvent) -> None:
                         f"uid: {uid}): {status}"
                     )
                     if self.current_scan_id is not None:
+                        if self.auto_analyse:
+                            try:
+                                fit_res = self.fit_scan(self.current_scan_id)
+                                scan_entry = self.data[self.current_scan_id]
+                                if fit_res.gaussian.success:
+                                    scan_entry["gaussian"] = (
+                                        fit_res.gaussian.best_fit.tolist()
+                                    )
+                                    scan_entry["centre"] = [
+                                        fit_res.gaussian.best_values["center"]
+                                    ]
+                                if not np.isnan(fit_res.com):
+                                    scan_entry["com"] = [fit_res.com]
+                                if fit_res.df.d1.size > 0:
+                                    scan_entry.update(
+                                        {
+                                            "d1": fit_res.df.d1.tolist(),
+                                            "d2": fit_res.df.d2.tolist(),
+                                            "zero_crossing": fit_res.df.zero_crossings,
+                                            "max_slope_x": [fit_res.df.max_slope_x],
+                                            "min_slope_x": [fit_res.df.min_slope_x],
+                                        }
+                                    )
+
+                            except Exception as e:
+                                click.echo(f"Auto-fit failed: {e}")
                         self._save_scan_to_cache(str(self.current_scan_id))
                     self.current_scan_id = None
 

src/sm_bluesky/common/clients/interactive_plot.py

@@ -26,6 +26,14 @@ class PlotMode(StrEnum):
 
 DEFAULT_CMAP = "viridis"
 NO_FIELDS = "(no data)"
+DERIVED_Y_FIELDS = ("gaussian", "d1", "d2")
+SCALAR_META_FIELDS = {
+    "centre",
+    "com",
+    "zero_crossing",
+    "max_slope_x",
+    "min_slope_x",
+}
 
 
 class InteractivePlotWindow:
@@ -143,7 +151,10 @@ def get_active_fields(self) -> list[str]:
         fields = set()
         for s_id in self.active_scans:
             fields.update(self.data_source[s_id].keys())
-        return sorted(fields)
+        fields -= SCALAR_META_FIELDS
+        device_fields = sorted(fields - set(DERIVED_Y_FIELDS))
+        derived_fields = [f for f in DERIVED_Y_FIELDS if f in fields]
+        return device_fields + derived_fields
 
     def get_visible_scans(self) -> list[Any]:
         """Return up to 10 active and most recent scan IDs in scan order."""
@@ -172,7 +183,8 @@ def get_2d_fields(self) -> list[str]:
         s_id = self.get_2d_scan()
         if s_id is None:
             return []
-        return list(self.data_source[s_id].keys())
+        excluded = SCALAR_META_FIELDS | set(DERIVED_Y_FIELDS)
+        return [k for k in self.data_source[s_id].keys() if k not in excluded]
 
     def _clear_widgets(self, *names: str) -> None:
         """Remove and disconnect the specified Matplotlib widgets."""
@@ -252,7 +264,8 @@ def _build_controls_1d(self, old_x: str | None = None) -> None:
         """Build X-axis radio buttons and Y-channel checkboxes for 1D mode."""
         self.ax_x.set_title("X Axis", fontsize=10, fontweight="bold")
         fields = self.get_active_fields()
-        x_options = ["Point number"] + fields
+        x_fields = [f for f in fields if f not in DERIVED_Y_FIELDS]
+        x_options = ["Point number"] + x_fields
 
         active_x_idx = 0
         if old_x in x_options:

src/sm_bluesky/common/math/__init__.py

@@ -1,3 +1,9 @@
+from .analysis import (
+    DerivativeResult,
+    calculate_com,
+    calculate_derivatives,
+    fit_gaussian,
+)
 from .conversion import cal_range_num, step_size_to_step_num, to_float
 from .rebinning import (
     binned_axis,
@@ -8,6 +14,10 @@
 )
 
 __all__ = [
+    "DerivativeResult",
+    "calculate_com",
+    "calculate_derivatives",
+    "fit_gaussian",
     "cal_range_num",
     "edges_from_centres",
     "step_size_to_step_num",