Experiment updates

## experiment.py
 - replaced _scan_numbers with all_scan_numbers
 - allow scans to be loaded from NexusScan objects
 - add get_scan_number()
 - add get_nearby_scan_numbers()
 - improved find_scans
 - add get_value_changes()

Add examples and tests.

All tests pass.

## xas/container_functions.py
 - add find_similar_measurements

Author: DanPorter

docs/examples/examples.md

@@ -55,7 +55,7 @@ for scan in exp[-5:]:
 
 ### Plotting
 ```python
-from mmg_toolbox import data_file_reader
+from mmg_toolbox import data_file_reader, Experiment
 from mmg_toolbox.plotting.matplotlib import set_plot_defaults
 set_plot_defaults()  # changes the matplotlib defaults to larger fonts, thicker lines etc.
 
@@ -66,4 +66,7 @@ scan.plot.plot(xaxis='axes', yaxis='signal')
 # multi-line plot
 scan.plot.plot(xaxis='axes', yaxis=['signal', 'signal2'])
 
+# Plotting multiple files
+exp = Experiment('/dls/i06/2025/mm1234-1', instrument='i06') 
+exp.plot.multi_lines(*range(12345, 12355), value='Ta')
 ```

examples/example_2d_surface_scans.py

@@ -0,0 +1,82 @@
+"""
+mmg_toolbox example
+Example script to load multiple files from a 2D grid scan
+"""
+
+import numpy as np
+import matplotlib.pyplot as plt
+from mmg_toolbox import Experiment
+
+
+exp = Experiment('/dls/science/groups/das/ExampleData/i16/2dscans_2026/', instrument='i16')
+
+# Add custom ROI
+exp.add_roi('new_roi', cen_i=95, cen_j=200, wid_i=40, wid_j=45, image_name='pil3_100k')
+
+scan_numbers = range(1144247, 1144272)
+scans = exp.scans(*scan_numbers)
+ttl = exp.generate_scans_title(*scan_numbers)
+
+xx, yy, image = exp.generate_mesh(*scan_numbers, axes='sx', values='sy', signal='signal')
+
+plt.figure()
+plt.pcolormesh(xx, yy, image)
+plt.axis('image')
+cb = plt.colorbar(label='Intensity')
+plt.xlabel('sx',fontsize=22)
+plt.ylabel('sy',fontsize=22)
+plt.title(ttl,fontsize=18)
+
+"""
+# perform fits and extract data
+sx_values = np.zeros(len(scans))
+sy_values = np.zeros(len(scans))
+area_values = np.zeros(len(scans))
+for n, scan in enumerate(scans):
+    s = scan.map.get_image_shape()
+    scan.map.add_roi('pilsum', s[0]//2, s[1]//2, s[0], s[1])  # add ROI for whole detector
+    result = scan.fit.multi_peak_fit(
+        xaxis='axes',  # default scan axes
+        yaxis='pilsum_total / Transmission / (rc / 300)', #'signal',  # default scan values
+        npeaks=None, # determine the number of peaks automatically
+        min_peak_power=None,
+        peak_distance_idx=6,
+        model='Gaussian',
+        background='Slope',
+        plot_result=True
+    )
+    sx_values[n] = scan('mean(sx)')
+    sy_values[n] = scan('mean(sy)')
+    area_values[n] = result.amplitude  # total area under all peaks
+    #area_values[n] = result.fwhm  # width
+    #area_values[n] = result.center  # centre
+
+
+# Determine the repeat length of the scans
+delta_x = np.abs(np.diff(sx_values))
+ch_idx_x = np.where(delta_x > delta_x.max()*0.9) # find biggest changes
+ch_delta_x = np.diff(ch_idx_x)
+rep_len_x = np.round(np.mean(ch_delta_x))
+delta_y = np.abs(np.diff(sy_values))
+ch_idx_y = np.where(delta_y > delta_y.max()*0.9) # find biggest changes
+ch_delta_y = np.diff(ch_idx_y)
+rep_len_y = np.round(np.mean(ch_delta_y))
+print('Scans in sx are repeating every {} iterations'.format(rep_len_x))
+print('Scans in sy are repeating every {} iterations'.format(rep_len_y))
+rep_len = int(max(rep_len_x, rep_len_y))
+
+# Reshape into square arrays
+# If this is problematic, look at scipy.interpolate.griddata
+sx_squareA = sx_values[:rep_len*(len(sx_values)//rep_len)].reshape(-1,rep_len)
+sy_squareA = sy_values[:rep_len*(len(sy_values)//rep_len)].reshape(-1,rep_len)
+area_squareA = area_values[:rep_len*(len(area_values)//rep_len)].reshape(-1,rep_len)
+
+plt.figure(figsize=(12, 10), dpi=60)
+plt.pcolormesh(sx_squareA, sy_squareA, area_squareA, shading='auto')
+plt.axis('image')
+cb = plt.colorbar(label='Amplitude')
+plt.xlabel('sx',fontsize=22)
+plt.ylabel('sy',fontsize=22)
+plt.title(ttl,fontsize=18)
+cb.set_label('Integrated Area',fontsize=22)
+"""

examples/example_experiment.py

@@ -3,16 +3,35 @@
 Example script using the Experiment class to get files from data folders
 """
 
+import numpy as np
 import matplotlib.pyplot as plt
 from mmg_toolbox import Experiment
 
 datadir1 = '/dls/science/groups/das/ExampleData/i16/azimuths'
 datadir2 = '/dls/science/groups/das/ExampleData/hdfmap_tests/i16/cm37262-1'
 exp = Experiment(datadir1, datadir2)
-all_scans = exp.all_scan_numbers()
+all_scan_nos = exp.all_scan_numbers()
 
-print('\n'.join(exp.scans_str(*all_scans)))
+# Print information - all scans
+print('\n'.join(exp.scans_str(*all_scan_nos)))
+# single scan
+print(exp.scan_str(-1))  # -1 is a shortcut for the latest scan, 0 is the first scan
 
+# Load scans
+scan = exp.scan(1032510)  # NexusDataHolder object (loads data)
+print(repr(scan))
+scan1, scan2, scan3 = exp.scans(*range(-3, 0))  # NexusScan object (lazy loader)
+
+# Search for scans
+scan_nos = exp.get_nearby_scan_numbers(1032510)
+scans = exp.find_scans(*scan_nos, **{'_cmd': 'scan energy'})
+
+# Find values that change between scans
+values = exp.get_value_changes(*scan_nos)
+print('\n'.join(f"{name} : {np.std(array):.2f}" for name, array in values.items()))
+
+
+# Plotting
 rng = range(1032510, 1032521)
 # exp.plot(*rng)
 exp.plot.surface_3d(*rng)
@@ -21,9 +40,9 @@
 
 exp.plot.detail(-1)
 
-exp.plot.surface_2d(*all_scans[:5], xaxis='eta_fly', signal='signal', values='psi')
-exp.plot.surface_3d(*all_scans[:5], xaxis='eta_fly', signal='signal', values='psi')
+exp.plot.surface_2d(*all_scan_nos[:5], xaxis='eta_fly', signal='signal', values='psi')
+exp.plot.surface_3d(*all_scan_nos[:5], xaxis='eta_fly', signal='signal', values='psi')
 
-exp.plot.lines_3d(*all_scans[:5], xaxis='eta_fly', signal='signal', values='psi')
+exp.plot.lines_3d(*all_scan_nos[:5], xaxis='eta_fly', signal='signal', values='psi')
 
 plt.show()

examples/example_i16_xmcd.py

@@ -7,28 +7,53 @@
 import numpy as np
 import matplotlib.pyplot as plt
 from mmg_toolbox import Experiment
-from mmg_toolbox.xas import average_polarised_scans, average_scans
+from mmg_toolbox.xas import average_polarised_scans
 
 
-exp = Experiment('/dls/i16/data/2026/cm44164-9', instrument='i16')
+exp = Experiment('/dls/science/groups/das/ExampleData/i16/vortex_2026', instrument='i16')
 
-spectra = exp.load_xas(1145818)
+### Single Energy Vortex spectrum ###
+spectrum, = exp.load_xas(1146115, element_edge='Ir L32')  # scan x 1 3 1 xsp3 1
+m = spectrum.metadata
+print(f"Spectrum: {spectrum.name} T={m.temp:.1f} K, B={m.mag_field:+3.1f} T, pol='{m.pol}'")
+
+# Processing
+processed_spectrum = spectrum.trim(1000, 20000).remove_background('flat')
+fig = processed_spectrum.create_figure()
+
+### Incident Energy Scan with Vortex ###
+scan, = exp.scans(1146116)
+# Add ROI window - xsp3 has 2 channels, 4096 bins. create window around bin 915 == 9150 eV
+scan.map.add_roi('new_window', 1, 915, 2, 100, 'xsp3')
+spectrum = scan.xas_spectra()
+
+fig2 = spectrum.create_figure()  # show the different Windows
 
-for spectrum in spectra:
-    m = spectrum.metadata
-    print(f"Spectrum: {spectrum.name} T={m.temp:.1f} K, B={m.mag_field:+3.1f} T, pol='{m.pol}'")
+
+### Merge and compare scans
+exp.add_roi('my_window', 1, 915, 2, 100, 'xsp3')  # Add ROI to exp so it is added to every scan
+spectra = exp.load_xas(*range(1146116, 1146124))
 
 # Processing
 spectra = [
-    spectrum.trim(1000, 30000).remove_background('flat')
+    spectrum.trim(5, 0).remove_background('flat')
     for spectrum in spectra
 ]
 
-av = average_scans(*spectra)
+# Find opposite polarisations, sum similar polarisations
+pol1, pol2 = average_polarised_scans(*spectra)
+print(f"Found {len(pol1.parents)} spectra with '{pol1.metadata.pol}' polarisation," +
+      " and {len(pol2.parents)} spectra with '{pol2.metadata.pol}' polarisation")
+
+# Align spectra using cross-correlation
+pol1 = pol1.align_spectra('my_window')
+pol2 = pol2.align_spectra('my_window')
+
+# Calculate subtraction
+xmld = pol1 - pol2
 
-print(repr(av))
+fig3 = xmld.create_sum_rules_figure(figsize=[12, 10], dpi=100)
 
-av.create_figure()
 plt.show()
 
 

mmg_toolbox/beamline_metadata/config.py

@@ -32,7 +32,7 @@ class C:
 }
 
 REPLACE_NAMES = {
-    # NEW_NAME: EXPRESSION
+    # NEW_NAME: EXPRESSION (cannot include NEW_NAME)  #TODO: fix this in HdfMap
     '_t': '(count_time|counttime|t?(1.0))',
     '_cmd': DEFAULT_SCAN_DESCRIPTION.strip('{}'),
 }

mmg_toolbox/plotting/exp_plot_manager.py

@@ -5,7 +5,7 @@
 import numpy as np
 import matplotlib.pyplot as plt
 import hdfmap
-from ..utils.experiment import Experiment
+from ..utils.experiment import Experiment, ScanFile
 from .matplotlib import (
     set_plot_defaults, generate_subplots, plot_lines, plot_2d_surface, plot_3d_surface, plot_3d_lines,
     FIG_SIZE, FIG_DPI, DEFAULT_CMAP, Axes3D
@@ -32,7 +32,7 @@ def __init__(self, experiment: Experiment):
     def __call__(self, *args, **kwargs):
         return self.plot(*args, **kwargs)
 
-    def plot(self, *scan_files: int | str, hdf_map: hdfmap.NexusMap | None = None,
+    def plot(self, *scan_files: ScanFile, hdf_map: hdfmap.NexusMap | None = None,
              xaxis: str = 'axes', yaxis: str | list[str] = 'signal', axes: plt.Axes | None = None, **kwargs) -> plt.Axes:
         """
         Create matplotlib figure with a line plot from a or several scans
@@ -57,7 +57,7 @@ def plot(self, *scan_files: int | str, hdf_map: hdfmap.NexusMap | None = None,
         axes.set_title(self.exp.generate_scans_title(*scan_files))
         return axes
 
-    def image(self,  scan_file: int | str = -1, hdf_map: hdfmap.NexusMap | None = None,
+    def image(self,  scan_file: ScanFile = -1, hdf_map: hdfmap.NexusMap | None = None,
               index: int | tuple | slice | None = None, xaxis: str = 'axes',
               axes: plt.Axes | None = None, clim: tuple[float, float] | None = None,
               cmap: str = DEFAULT_CMAP, colorbar: bool = False, **kwargs) -> plt.Axes | None:
@@ -79,7 +79,7 @@ def image(self,  scan_file: int | str = -1, hdf_map: hdfmap.NexusMap | None = No
             return scan.plot.image(index, xaxis, axes, clim, cmap, colorbar, **kwargs)
         return None
 
-    def detail(self, scan_file: int | str = -1, hdf_map: hdfmap.NexusMap | None = None,
+    def detail(self, scan_file: ScanFile = -1, hdf_map: hdfmap.NexusMap | None = None,
                xaxis: str = 'axes', yaxis: str | list[str] = 'signal',
                index: int | tuple | slice | None = None, clim: tuple[float, float] | None = None,
                cmap: str = DEFAULT_CMAP, **kwargs) -> plt.Figure:
@@ -107,8 +107,9 @@ def detail(self, scan_file: int | str = -1, hdf_map: hdfmap.NexusMap | None = No
         lt.set_title(None)
 
         # Top right - image plot
-        scan.plot.image(index, xaxis, cmap=cmap, clim=clim, axes=rt)
-        rt.set_title(None)
+        if scan.map.image_data:
+            scan.plot.image(index, xaxis, cmap=cmap, clim=clim, axes=rt)
+            rt.set_title(None)
 
         # Bottom-Left - details
         details = self.exp.scan_str(scan_file, hdf_map=hdf_map)
@@ -122,7 +123,7 @@ def detail(self, scan_file: int | str = -1, hdf_map: hdfmap.NexusMap | None = No
             rb.text(-0.2, 1, fit_report, ha='left', va='top', multialignment="left", fontsize=12, wrap=True)
         return fig
 
-    def multi_lines(self, *scan_files: int | str, hdf_map: hdfmap.NexusMap | None = None,
+    def multi_lines(self, *scan_files: ScanFile, hdf_map: hdfmap.NexusMap | None = None,
                     xaxis: str = 'axes', yaxis: str = 'signal', value: str | None = None,
                     axes: plt.Axes | None = None, **kwargs) -> plt.Axes:
         """
@@ -162,7 +163,7 @@ def multi_lines(self, *scan_files: int | str, hdf_map: hdfmap.NexusMap | None =
         plt.colorbar(sm, ax=axes, label=value_label)
         return axes
 
-    def multi_plot(self, *scan_files: int | str, hdf_map: hdfmap.NexusMap | None = None,
+    def multi_plot(self, *scan_files: ScanFile, hdf_map: hdfmap.NexusMap | None = None,
                    xaxis: str = 'axes', yaxis: str | list[str] = 'signal', value: str | None = None,
                    subplots: tuple[int, int] = (4, 4), **kwargs) ->  list[tuple[plt.Figure, plt.Axes]]:
         """
@@ -190,7 +191,29 @@ def multi_plot(self, *scan_files: int | str, hdf_map: hdfmap.NexusMap | None = N
             ax.set_title(ttl_exp)
         return fig_ax
 
-    def surface_2d(self, *scan_files: int | str, hdf_map: hdfmap.NexusMap | None = None,
+    def multi_image(self,  *scan_files: ScanFile, hdf_map: hdfmap.NexusMap | None = None,
+                    index: int | tuple | slice | None = None, xaxis: str = 'axes',
+                    clim: tuple[float, float] | None = None, cmap: str = DEFAULT_CMAP,
+                    colorbar: bool = False, **kwargs) -> list[tuple[plt.Figure, plt.Axes]]:
+        """
+        Plot scan images in matplotlib figure (if available)
+        :param scan_files: scan number or filename (multiple allowed)
+        :param hdf_map: hdfmap object or None
+        :param index: int, detector image index, 0-length of scan, if None, use centre index
+        :param xaxis: name or address of xaxis dataset
+        :param clim: [min, max] colormap cut-offs (None for auto)
+        :param cmap: str colormap name (None for auto)
+        :param colorbar: False/ True add colorbar to plot
+        :param kwargs: additional arguments for plot_detector_image
+        :return: [(Figure, Axes)] for each scan
+        """
+        fig_axes = generate_subplots(len(scan_files))
+        for scan, (fig, axes) in zip(scan_files, fig_axes):
+            self.image(scan, axes=axes, hdf_map=hdf_map, index=index, xaxis=xaxis,
+                       clim=clim, cmap=cmap, colorbar=colorbar, **kwargs)
+        return fig_axes
+
+    def surface_2d(self, *scan_files: ScanFile, hdf_map: hdfmap.NexusMap | None = None,
                    xaxis: str = 'axes', signal: str = 'signal', values: str | None = None,
                    axes: plt.Axes | None = None, clim: tuple[float, float] | None = None,
                    axlim: str = 'image', **kwargs) -> plt.Axes:
@@ -226,7 +249,7 @@ def surface_2d(self, *scan_files: int | str, hdf_map: hdfmap.NexusMap | None = N
         axes.figure.colorbar(surf, ax=axes, label=signal)
         return axes
 
-    def lines_3d(self, *scan_files: int | str, hdf_map: hdfmap.NexusMap | None = None,
+    def lines_3d(self, *scan_files: ScanFile, hdf_map: hdfmap.NexusMap | None = None,
                  xaxis: str = 'axes', signal: str = 'signal', values: str | None = None,
                  axes: Axes3D | None = None, legend: bool = False, **kwargs) -> Axes3D:
         """
@@ -265,7 +288,7 @@ def lines_3d(self, *scan_files: int | str, hdf_map: hdfmap.NexusMap | None = Non
             axes.legend()
         return axes
 
-    def surface_3d(self, *scan_files: int | str, hdf_map: hdfmap.NexusMap | None = None,
+    def surface_3d(self, *scan_files: ScanFile, hdf_map: hdfmap.NexusMap | None = None,
                    xaxis: str = 'axes', signal: str = 'signal', values: str | None = None,
                    axes: Axes3D | None = None, clim: tuple[float, float] | None = None,
                    axlim: str = 'image', **kwargs) -> Axes3D:
@@ -302,7 +325,7 @@ def surface_3d(self, *scan_files: int | str, hdf_map: hdfmap.NexusMap | None = N
         axes.set_zlabel(signal_label)
         return axes
 
-    def metadata(self, *scan_files: int | str, values: str | list[str], hdf_map: hdfmap.NexusMap | None = None,
+    def metadata(self, *scan_files: ScanFile, values: str | list[str], hdf_map: hdfmap.NexusMap | None = None,
                  axes: plt.Axes | None = None) -> plt.Axes:
         """
         Create matplotlib figure with plot of metadata vs scan number