mfx102101026 Fe Kα pipeline steps: event-code mask, spatial background, smalldata override, Zr calibration

XSpect/analysis/spectroscopy.py

@@ -54,6 +54,52 @@ def get_run_shot_properties(run, **kwargs):
     run.get_run_shot_properties()
 
 
+@register_step("make_eventcode_mask")
+def make_eventcode_mask(run, **kwargs):
+    """Build a per-shot boolean mask from an EVR/timing event-code column.
+
+    The smalldata timing group stores a per-shot event-code table (shape
+    n_shots x n_codes, e.g. timing/eventcodes with 288 columns). This step
+    extracts one code's column as a boolean shot mask, so it can be used by
+    union_shots / filter_shots like any other mask.
+
+    Typical use: beam is delivered at 30 Hz inside a 120 Hz DAQ, tagged by
+    event code 198. Gate the analysis on code 198 to drop the 90 Hz of empty
+    frames (the smalldata `xray` flag does NOT distinguish these).
+
+    Parameters from YAML:
+        code:        event-code index (column) to extract, e.g. 198
+        eventcodes:  attribute holding the 2D code table (default "eventcodes")
+        new_key:     output mask attribute name (default: "ec<code>")
+    """
+    code = kwargs.get("code")
+    ec_key = kwargs.get("eventcodes", "eventcodes")
+    new_key = kwargs.get("new_key", None)
+    if code is None:
+        run.update_status("make_eventcode_mask: 'code' is required")
+        return
+    table = getattr(run, ec_key, None)
+    if table is None:
+        run.update_status(
+            f"make_eventcode_mask: '{ec_key}' not loaded (add it to data.keys)"
+        )
+        return
+    table = np.asarray(table)
+    if table.ndim != 2 or code >= table.shape[1]:
+        run.update_status(
+            f"make_eventcode_mask: '{ec_key}' shape {table.shape} cannot index code {code}"
+        )
+        return
+    mask = table[:, code].astype(bool)
+    if new_key is None:
+        new_key = f"ec{code}"
+    setattr(run, new_key, mask)
+    run.update_status(
+        f"make_eventcode_mask: {ec_key}[:, {code}] -> {new_key} "
+        f"({int(mask.sum())}/{mask.size} shots = {100 * mask.mean():.1f}%)"
+    )
+
+
 @register_step("filter_shots")
 def filter_shots(run, **kwargs):
     """Filter a shot mask by thresholding on another key.
@@ -478,6 +524,83 @@ def apply_roi(run, **kwargs):
     run.update_status(f"Applied ROI to {detector_key}")
 
 
+@register_step("subtract_spatial_background")
+def subtract_spatial_background(run, **kwargs):
+    """Subtract a per-line background estimated from flanking regions.
+
+    Designed for a dispersed spectral streak sitting on a smooth, spatially
+    slowly-varying background (e.g. isotropic fluorescence under a Von Hamos
+    emission line). For each line along the dispersion axis, the background
+    level is estimated from one or two "side band" regions flanking the signal
+    (on the cross-dispersion axis) and subtracted from the whole line.
+
+    Works on a 2D image (rows x cols). ``bkg_axis`` is the CROSS-DISPERSION
+    axis, i.e. the axis along which the side bands are taken and which is
+    reduced away by the subsequent reduce_detector_spatial step. The background
+    per line is: mean over the side-band pixels, scaled to the number of pixels
+    in the signal band, then subtracted from the signal band; side-band columns
+    themselves are zeroed so they do not contribute downstream.
+
+    Parameters from YAML:
+        on:            detector key (2D rows x cols)
+        signal:        [start, end] of the signal band on bkg_axis
+        sidebands:     list of [start, end] flanking regions on bkg_axis used to
+                       estimate the background (1 or 2 regions typical)
+        bkg_axis:      cross-dispersion axis (default 1 = columns)
+        new_key:       output key (default: "<on>_bkgsub"); the original is kept
+        estimator:     "mean" or "median" over side-band pixels (default "median")
+    """
+    detector_key = kwargs.get("on")
+    signal = kwargs.get("signal")
+    sidebands = kwargs.get("sidebands")
+    bkg_axis = kwargs.get("bkg_axis", 1)
+    new_key = kwargs.get("new_key", None)
+    estimator = kwargs.get("estimator", "median")
+
+    if detector_key is None or signal is None or not sidebands:
+        run.update_status(
+            "subtract_spatial_background: 'on', 'signal' and 'sidebands' required"
+        )
+        return
+    data = getattr(run, detector_key, None)
+    if data is None:
+        return
+    if data.ndim != 2:
+        run.update_status(
+            f"subtract_spatial_background: expected 2D image, got ndim={data.ndim}"
+        )
+        return
+
+    # Move the cross-dispersion axis to position 1 so we can index columns.
+    work = data if bkg_axis == 1 else data.T  # (lines, cross)
+    est_fn = np.nanmedian if estimator == "median" else np.nanmean
+
+    # Per-line background level (one value per dispersion line) from side bands.
+    side_pixels = []
+    for lo, hi in sidebands:
+        side_pixels.append(work[:, lo:hi])
+    side = np.concatenate(side_pixels, axis=1)
+    bkg_per_pixel = est_fn(side, axis=1, keepdims=True)  # (lines, 1)
+
+    out = work.astype(float).copy()
+    s0, s1 = signal
+    out[:, s0:s1] = out[:, s0:s1] - bkg_per_pixel  # subtract per-pixel bkg
+    # zero everything outside the signal band so downstream reduction only
+    # integrates the background-subtracted signal columns
+    keep = np.zeros(work.shape[1], dtype=bool)
+    keep[s0:s1] = True
+    out[:, ~keep] = 0.0
+
+    result = out if bkg_axis == 1 else out.T
+    if new_key is None:
+        new_key = f"{detector_key}_bkgsub"
+    setattr(run, new_key, result)
+    run.update_status(
+        f"subtract_spatial_background: {detector_key} -> {new_key} "
+        f"signal={signal} sidebands={sidebands} estimator={estimator}"
+    )
+
+
 @register_step("time_binning")
 def time_binning(run, **kwargs):
     """Create time delay bins from laser timing data.

XSpect/analysis/xes.py

@@ -109,6 +109,138 @@ def make_energy_axis(run, **kwargs):
     )
 
 
+# Elemental emission line energies (eV). Source: standard X-ray data tables
+# (e.g. Bearden / NIST). Used by calibrate_energy to fit the geometry against a
+# known foil line. Add elements/lines as needed.
+EMISSION_LINES = {
+    "Ti": {"Ka1": 4510.84, "Ka2": 4504.86, "Kb1": 4931.81},
+    "V": {"Ka1": 4952.20, "Ka2": 4944.64, "Kb1": 5427.29},
+    "Cr": {"Ka1": 5414.72, "Ka2": 5405.509, "Kb1": 5946.71},
+    "Mn": {"Ka1": 5898.75, "Ka2": 5887.65, "Kb1": 6490.45},
+    "Fe": {"Ka1": 6403.84, "Ka2": 6390.84, "Kb1": 7057.98},
+    "Co": {"Ka1": 6930.32, "Ka2": 6915.30, "Kb1": 7649.43},
+    "Ni": {"Ka1": 7478.15, "Ka2": 7460.89, "Kb1": 8264.66},
+    "Cu": {"Ka1": 8047.78, "Ka2": 8027.83, "Kb1": 8905.29},
+    "Zn": {"Ka1": 8638.86, "Ka2": 8615.78, "Kb1": 9572.0},
+    "Ga": {"Ka1": 9251.74, "Ka2": 9224.82, "Kb1": 10264.2},
+    "Ge": {"Ka1": 9886.42, "Ka2": 9855.32, "Kb1": 10982.1},
+    "As": {"Ka1": 10543.72, "Ka2": 10507.99, "Kb1": 11726.2},
+    "Se": {"Ka1": 11222.4, "Ka2": 11181.4, "Kb1": 12495.9},
+    "Br": {"Ka1": 11924.2, "Ka2": 11877.6, "Kb1": 13291.4},
+    "Kr": {"Ka1": 12649.0, "Ka2": 12598.0, "Kb1": 14112.0},
+    "Rb": {"Ka1": 13395.3, "Ka2": 13335.8, "Kb1": 14961.3},
+    "Sr": {"Ka1": 14165.0, "Ka2": 14097.9, "Kb1": 15835.7},
+    "Y": {"Ka1": 14958.4, "Ka2": 14882.9, "Kb1": 16737.8},
+    "Zr": {"Ka1": 15775.1, "Ka2": 15690.9, "Kb1": 17667.8},
+}
+
+
+@register_step("calibrate_energy")
+def calibrate_energy(run, **kwargs):
+    """Calibrate the vonHamos crystal_detector_distance (A) against a known line.
+
+    The energy axis (make_energy_axis) is
+        ll   = pixel*mm_per_pixel/2 - (max(gl)-min(gl))/4
+        E(p) = hc / (2 d sin(arctan(R / (ll(p) + A))))
+    Everything except A is fixed by the spectrometer. Given a foil emission line
+    of known energy E0 whose peak falls at pixel p0 in a measured spectrum, A is
+    solved in closed form:
+        s     = hc / (2 d E0)
+        A     = R * sqrt(1 - s^2) / s - ll(p0)
+    This step finds p0 as the argmax of the spectrum (optionally within a pixel
+    window), solves A, and rebuilds `{name}_energy` so the line sits exactly at
+    E0. The fitted A is stored as `{name}_calibrated_A` (and returned in the
+    result dict via run.results) so it can be reused for subsequent, non-foil
+    measurements on the same spectrometer.
+
+    Parameters from YAML
+    --------------------
+    on            : spectrum key to locate the peak in (e.g. epix_reduced_ROI_1)
+    element       : element symbol for the foil, e.g. "Mn" (uses EMISSION_LINES)
+    line          : which line, "Ka1" (default), "Ka2", or "Kb1"
+    energy        : explicit line energy in eV (overrides element/line lookup)
+    crystal_radius, d_spacing, mm_per_pixel : geometry (same as make_energy_axis)
+    n_pixels      : pixel count (defaults to len of the `on` spectrum)
+    peak_window   : optional [lo, hi] pixel range to search for the peak
+    name          : output prefix (default "xes"); writes {name}_energy and
+                    {name}_calibrated_A
+    """
+    spec_key = kwargs.get("on")
+    element = kwargs.get("element")
+    line = kwargs.get("line", "Ka1")
+    energy = kwargs.get("energy", None)
+    R = kwargs.get("crystal_radius")
+    d = kwargs.get("d_spacing")
+    mm_per_pixel = kwargs.get("mm_per_pixel", 0.05)
+    peak_window = kwargs.get("peak_window", None)
+    name = kwargs.get("name", "xes")
+
+    if spec_key is None or R is None or d is None:
+        run.update_status("calibrate_energy: 'on', crystal_radius, d_spacing required")
+        return
+
+    # Resolve the reference energy
+    if energy is None:
+        if element is None or element not in EMISSION_LINES:
+            run.update_status(
+                f"calibrate_energy: provide 'energy' or a known 'element' "
+                f"(got element={element!r})"
+            )
+            return
+        if line not in EMISSION_LINES[element]:
+            run.update_status(
+                f"calibrate_energy: line {line!r} not tabulated for {element}"
+            )
+            return
+        energy = EMISSION_LINES[element][line]
+
+    spec = getattr(run, spec_key, None)
+    if spec is None:
+        run.update_status(f"calibrate_energy: spectrum '{spec_key}' not found")
+        return
+    spec = np.asarray(spec, dtype=np.float64)
+    if spec.ndim > 1:
+        # collapse any leading axes; energy runs along the last (dispersion) axis
+        spec = spec.reshape(-1, spec.shape[-1]).sum(axis=0)
+
+    n_pixels = kwargs.get("n_pixels", spec.shape[-1])
+
+    # Locate the peak pixel (optionally restricted to a window)
+    if peak_window is not None:
+        lo, hi = int(peak_window[0]), int(peak_window[1])
+        p0 = lo + int(np.argmax(spec[lo:hi]))
+    else:
+        p0 = int(np.argmax(spec))
+
+    hc = 12398.42  # eV * Angstrom
+    gl = np.arange(n_pixels, dtype=np.float64) * mm_per_pixel
+    ll = gl / 2.0 - (np.amax(gl) - np.amin(gl)) / 4.0
+
+    # Closed-form solve for A so that E(p0) == energy
+    s = hc / (2.0 * d * energy)
+    if not (0.0 < s < 1.0):
+        run.update_status(
+            f"calibrate_energy: unphysical sin(theta)={s:.4f} for E={energy} "
+            f"(check d_spacing)"
+        )
+        return
+    A = R * np.sqrt(1.0 - s * s) / s - ll[p0]
+
+    energy_axis = hc / (2.0 * d * np.sin(np.arctan(R / (ll + A))))
+    setattr(run, f"{name}_energy", energy_axis)
+    setattr(run, f"{name}_calibrated_A", float(A))
+    # expose in results for reuse as the calibration on later measurements
+    if hasattr(run, "results") and isinstance(run.results, dict):
+        run.results[f"{name}_calibrated_A"] = float(A)
+        run.results[f"{name}_energy"] = energy_axis
+
+    run.update_status(
+        f"calibrate_energy: {name} line {element or ''} {line}={energy:.2f} eV at "
+        f"pixel {p0} -> A={A:.4f} mm; axis {energy_axis.min():.1f}-"
+        f"{energy_axis.max():.1f} eV"
+    )
+
+
 @register_step("patch_pixels")
 def patch_pixels(run, **kwargs):
     """Repair bad pixels using polynomial fitting from neighbors.

XSpect/controller/config_parser.py

@@ -28,6 +28,9 @@ class ExperimentConfig:
     hutch: str
     experiment_id: str
     lcls_run: int
+    smalldata_dir: "str | None" = None  # Optional override for the smalldata
+    # directory. If set, XSpect looks for {exp}_Run{run:04d}.h5 here instead of
+    # the default /sdf/data/lcls/ds/{hutch}/{exp}/hdf5/smalldata search list.
 
 
 @dataclass(frozen=True)
@@ -107,6 +110,7 @@ def _parse_experiment(raw: dict) -> ExperimentConfig:
         hutch=str(raw["hutch"]),
         experiment_id=str(raw["experiment_id"]),
         lcls_run=int(raw["lcls_run"]),
+        smalldata_dir=(str(raw["smalldata_dir"]) if raw.get("smalldata_dir") else None),
     )
 
 

XSpect/controller/pipeline.py

@@ -281,7 +281,12 @@ def _create_experiment(self):
         """Create experiment object from config. Wraps the directory lookup failure gracefully."""
         cfg = self.config.experiment
         try:
-            exp = spectroscopy_experiment(cfg.lcls_run, cfg.hutch, cfg.experiment_id)
+            exp = spectroscopy_experiment(
+                cfg.lcls_run,
+                cfg.hutch,
+                cfg.experiment_id,
+                smalldata_dir=cfg.smalldata_dir,
+            )
         except Exception:
             exp = _MockExperiment(cfg.lcls_run, cfg.hutch, cfg.experiment_id)
         return exp

XSpect/model/experiment.py

@@ -2,7 +2,7 @@
 
 
 class experiment:
-    def __init__(self, lcls_run, hutch, experiment_id):
+    def __init__(self, lcls_run, hutch, experiment_id, smalldata_dir=None):
         """
         Initializes an experiment instance.
 
@@ -14,17 +14,26 @@ def __init__(self, lcls_run, hutch, experiment_id):
             Hutch name. Example: xcs
         experiment_id : str
             Experiment identifier. Example: xcsl1004021
+        smalldata_dir : str, optional
+            Explicit path to the directory holding the {exp}_Run{run:04d}.h5
+            smalldata files. If provided, this overrides the default location
+            search (useful for locally reprocessed data). Example:
+            /sdf/data/lcls/ds/mfx/mfx101592326/results/lbgee/hdf5/smalldata
         """
         self.lcls_run = lcls_run
         self.hutch = hutch
         self.experiment_id = experiment_id
+        self.smalldata_dir = smalldata_dir
         self.get_experiment_directory()
 
     def get_experiment_directory(self):
         """
         Determines and returns the directory of the experiment based on the hutch and experiment ID.
         It attempts the various paths LCLS has had over the years with recent S3DF paths being the first attempt.
 
+        If `smalldata_dir` was provided at construction, it takes priority over
+        the default location search.
+
         Returns
         -------
         str
@@ -35,15 +44,26 @@ def get_experiment_directory(self):
         Exception
             If the directory cannot be found.
         """
+        # Explicit override (e.g. locally reprocessed smalldata) takes priority.
+        if self.smalldata_dir:
+            if os.path.exists(self.smalldata_dir) and os.listdir(self.smalldata_dir):
+                self.experiment_directory = self.smalldata_dir
+                return self.smalldata_dir
+            raise Exception(
+                f"smalldata_dir '{self.smalldata_dir}' does not exist or is empty."
+            )
+
         experiment_directories = [
-            '/sdf/data/lcls/ds/%s/%s/hdf5/smalldata',
-            '/reg/data/drpsrcf/%s/%s/scratch/hdf5/smalldata',
-            '/cds/data/drpsrcf/%s/%s/scratch/hdf5/smalldata',
-            '/reg/d/psdm/%s/%s/hdf5/smalldata'
+            "/sdf/data/lcls/ds/%s/%s/hdf5/smalldata",
+            "/reg/data/drpsrcf/%s/%s/scratch/hdf5/smalldata",
+            "/cds/data/drpsrcf/%s/%s/scratch/hdf5/smalldata",
+            "/reg/d/psdm/%s/%s/hdf5/smalldata",
         ]
         for directory in experiment_directories:
             experiment_directory = directory % (self.hutch, self.experiment_id)
-            if os.path.exists(experiment_directory) and os.listdir(experiment_directory):
+            if os.path.exists(experiment_directory) and os.listdir(
+                experiment_directory
+            ):
                 self.experiment_directory = experiment_directory
                 return experiment_directory
         raise Exception("Unable to find experiment directory.")
@@ -53,6 +73,7 @@ class spectroscopy_experiment(experiment):
     """
     A class to represent a spectroscopy experiment.
     """
+
     def __init__(self, *args, **kwargs):
         super().__init__(*args, **kwargs)