Analysis steps

Registered pipeline steps and reductions. Each is a stateless function reached by name through the registry.

XSpect.analysis.registry

Step and reduction registry for the XSpect pipeline.

Steps are stateless functions with signature: step(run, kwargs) -> None Reductions receive all runs: reduction(runs, kwargs) -> dict

ReductionNotFoundError

Bases: KeyError

Raised when a reduction name is not found in the registry.

Source code in XSpect/analysis/registry.py

class ReductionNotFoundError(KeyError):
    """Raised when a reduction name is not found in the registry."""
    pass

StepNotFoundError

Bases: KeyError

Raised when a step name is not found in the registry.

Source code in XSpect/analysis/registry.py

class StepNotFoundError(KeyError):
    """Raised when a step name is not found in the registry."""
    pass

clear_registry()

Clear all registered steps and reductions. Mainly for testing.

Source code in XSpect/analysis/registry.py

def clear_registry():
    """Clear all registered steps and reductions. Mainly for testing."""
    _STEP_REGISTRY.clear()
    _REDUCTION_REGISTRY.clear()

get_reduction(name)

Look up a registered reduction by name.

Source code in XSpect/analysis/registry.py

def get_reduction(name: str) -> callable:
    """Look up a registered reduction by name."""
    try:
        return _REDUCTION_REGISTRY[name]
    except KeyError:
        raise ReductionNotFoundError(
            f"Reduction '{name}' not found. Available: {list(_REDUCTION_REGISTRY.keys())}"
        )

get_step(name)

Look up a registered step by name.

Source code in XSpect/analysis/registry.py

def get_step(name: str) -> callable:
    """Look up a registered step by name."""
    try:
        return _STEP_REGISTRY[name]
    except KeyError:
        raise StepNotFoundError(f"Step '{name}' not found. Available: {list(_STEP_REGISTRY.keys())}")

list_reductions()

Return all registered reduction names.

Source code in XSpect/analysis/registry.py

def list_reductions() -> list[str]:
    """Return all registered reduction names."""
    return list(_REDUCTION_REGISTRY.keys())

list_steps()

Return all registered step names.

Source code in XSpect/analysis/registry.py

def list_steps() -> list[str]:
    """Return all registered step names."""
    return list(_STEP_REGISTRY.keys())

register_reduction(name)

Decorator that registers a function as a reduction step.

Source code in XSpect/analysis/registry.py

def register_reduction(name: str):
    """Decorator that registers a function as a reduction step."""
    def decorator(func):
        if name in _REDUCTION_REGISTRY:
            raise ValueError(f"Reduction '{name}' is already registered")
        _REDUCTION_REGISTRY[name] = func
        func._reduction_name = name
        return func
    return decorator

register_step(name)

Decorator that registers a function as a pipeline step.

Source code in XSpect/analysis/registry.py

def register_step(name: str):
    """Decorator that registers a function as a pipeline step."""
    def decorator(func):
        if name in _STEP_REGISTRY:
            raise ValueError(f"Step '{name}' is already registered")
        _STEP_REGISTRY[name] = func
        func._step_name = name
        return func
    return decorator

XSpect.analysis.spectroscopy

Core spectroscopy pipeline steps.

Each function is registered via @register_step and follows the contract: step(run, **kwargs) -> None All results written to run.results[key].

apply_roi(run, **kwargs)

Extract ROI without spatial reduction (keep spatial dimension).

Parameters from YAML: on: detector key rois: list of [start, end] pixel ranges combine_rois: bool (default True)

Source code in XSpect/analysis/spectroscopy.py

@register_step("apply_roi")
def apply_roi(run, **kwargs):
    """Extract ROI without spatial reduction (keep spatial dimension).

    Parameters from YAML:
        on: detector key
        rois: list of [start, end] pixel ranges
        combine_rois: bool (default True)
    """
    detector_key = kwargs.get("on")
    rois = kwargs.get("rois", [[0, None]])
    combine = kwargs.get("combine_rois", True)
    if detector_key is None:
        return

    detector = getattr(run, detector_key, None)
    if detector is None:
        return

    if combine:
        mask = np.zeros(detector.shape[-1], dtype=bool)
        for roi in rois:
            end = roi[1] if roi[1] is not None else detector.shape[-1]
            mask[roi[0] : end] = True

        if detector.ndim == 3:
            masked = detector[:, :, mask]
        elif detector.ndim == 2:
            masked = detector[:, mask]
        else:
            masked = detector[mask]
        setattr(run, f"{detector_key}_ROI_1", masked)
    else:
        for idx, roi in enumerate(rois):
            end = roi[1] if roi[1] is not None else detector.shape[-1]
            if detector.ndim == 3:
                chunk = detector[:, :, roi[0] : end]
            else:
                chunk = detector[:, roi[0] : end]
            setattr(run, f"{detector_key}_ROI_{idx + 1}", chunk)

    run.update_status(f"Applied ROI to {detector_key}")

bin_uniques(run, **kwargs)

Bin unique scan variable values for scan-based analysis.

Parameters from YAML: on: key with scan variable values

Source code in XSpect/analysis/spectroscopy.py

@register_step("bin_uniques")
def bin_uniques(run, **kwargs):
    """Bin unique scan variable values for scan-based analysis.

    Parameters from YAML:
        on: key with scan variable values
    """
    key = kwargs.get("on")
    if key is None:
        return

    vals = getattr(run, key, None)
    if vals is None:
        return

    bins = np.unique(vals)
    addon = (bins[-1] - bins[-2]) / 2
    bins2 = np.append(bins, bins[-1] + addon)
    bins_center = np.empty_like(bins2)
    for ii in range(len(bins)):
        if ii == 0:
            bins_center[ii] = bins2[ii] - (bins2[ii + 1] - bins2[ii]) / 2
        else:
            bins_center[ii] = bins2[ii] - (bins2[ii] - bins2[ii - 1]) / 2
    bins_center[-1] = bins2[-1]

    run.scanvar_indices = np.digitize(vals, bins_center)
    run.scanvar_bins = bins_center
    run.update_status(f"Binned uniques on {key}: {len(bins)} unique values")

combine_runs(runs, **kwargs)

Aggregate time-binned data across multiple runs.

Sums the time-binned detector data and bin counts from all runs. Returns normalized difference (laser_on - laser_off) / laser_off.

Parameters from YAML: detector_key: base detector key (e.g., "epix_ROI_1") laser_on_suffix: suffix for laser-on data (default "_simultaneous_laser_time_binned") laser_off_suffix: suffix for laser-off data (default "_xray_not_laser_time_binned")

Source code in XSpect/analysis/spectroscopy.py

@register_reduction("combine_runs")
def combine_runs(runs, **kwargs):
    """Aggregate time-binned data across multiple runs.

    Sums the time-binned detector data and bin counts from all runs.
    Returns normalized difference (laser_on - laser_off) / laser_off.

    Parameters from YAML:
        detector_key: base detector key (e.g., "epix_ROI_1")
        laser_on_suffix: suffix for laser-on data (default "_simultaneous_laser_time_binned")
        laser_off_suffix: suffix for laser-off data (default "_xray_not_laser_time_binned")
    """
    detector_key = kwargs.get("detector_key", "epix_ROI_1")
    laser_on_suffix = kwargs.get("laser_on_suffix", "_simultaneous_laser_time_binned")
    laser_off_suffix = kwargs.get("laser_off_suffix", "_xray_not_laser_time_binned")

    laser_on_key = f"{detector_key}{laser_on_suffix}"
    laser_off_key = f"{detector_key}{laser_off_suffix}"
    laser_on_count_key = f"{detector_key}_simultaneous_laser_bincount"
    laser_off_count_key = f"{detector_key}_xray_not_laser_bincount"

    sum_on = None
    sum_off = None
    count_on = None
    count_off = None

    for r in runs:
        on_data = getattr(r, laser_on_key, None)
        off_data = getattr(r, laser_off_key, None)
        on_count = getattr(r, laser_on_count_key, None)
        off_count = getattr(r, laser_off_count_key, None)

        if on_data is not None:
            sum_on = on_data if sum_on is None else sum_on + on_data
        if off_data is not None:
            sum_off = off_data if sum_off is None else sum_off + off_data
        if on_count is not None:
            count_on = on_count if count_on is None else count_on + on_count
        if off_count is not None:
            count_off = off_count if count_off is None else count_off + off_count

    results = {
        "laser_on_summed": sum_on,
        "laser_off_summed": sum_off,
        "laser_on_count": count_on,
        "laser_off_count": count_off,
    }

    if (
        sum_on is not None
        and sum_off is not None
        and count_on is not None
        and count_off is not None
    ):
        safe_count_on = np.where(count_on > 0, count_on, 1)
        safe_count_off = np.where(count_off > 0, count_off, 1)
        avg_on = (
            sum_on / safe_count_on[:, np.newaxis]
            if sum_on.ndim == 2
            else sum_on / safe_count_on
        )
        avg_off = (
            sum_off / safe_count_off[:, np.newaxis]
            if sum_off.ndim == 2
            else sum_off / safe_count_off
        )

        safe_off = np.where(np.abs(avg_off) > 1e-10, avg_off, 1e-10)
        difference = (avg_on - avg_off) / safe_off

        results["laser_on_average"] = avg_on
        results["laser_off_average"] = avg_off
        results["difference"] = difference

    return results

common_mode_correction(run, **kwargs)

Subtract per-row, per-column, or per-bank common-mode offsets.

Detector electronics add a slowly varying baseline that is shared across all pixels in a readout unit (a row, a column, or an ePix100 bank). The offset is estimated from a reference region that carries no signal (a dark band) and subtracted from the whole unit. Runs on 3D detector data (shots x rows x cols) before shot reduction and preserves the input shape.

Parameters from YAML: on: detector key (3D: shots x rows x cols, or 2D: rows x cols) axis: "row" (default), "column", or "bank". "row" removes a per-row offset (shared across columns), "column" a per-column offset, "bank" a per-bank offset across fixed-width column blocks. method: "median" (default, robust to outliers) or "mean". reference: [start, end] pixel range of the signal-free band used to estimate the offset. Indexed along the axis orthogonal to the correction: for axis="row" it is a column range, for axis="column" (and "bank") it is a row range. Default: full extent. bank_size: int, column width of an ePix100 bank (default 128). Only used when axis="bank".

Source code in XSpect/analysis/spectroscopy.py

@register_step("common_mode_correction")
def common_mode_correction(run, **kwargs):
    """Subtract per-row, per-column, or per-bank common-mode offsets.

    Detector electronics add a slowly varying baseline that is shared across
    all pixels in a readout unit (a row, a column, or an ePix100 bank). The
    offset is estimated from a reference region that carries no signal (a dark
    band) and subtracted from the whole unit. Runs on 3D detector data
    (shots x rows x cols) before shot reduction and preserves the input shape.

    Parameters from YAML:
        on: detector key (3D: shots x rows x cols, or 2D: rows x cols)
        axis: "row" (default), "column", or "bank". "row" removes a per-row
            offset (shared across columns), "column" a per-column offset,
            "bank" a per-bank offset across fixed-width column blocks.
        method: "median" (default, robust to outliers) or "mean".
        reference: [start, end] pixel range of the signal-free band used to
            estimate the offset. Indexed along the axis orthogonal to the
            correction: for axis="row" it is a column range, for axis="column"
            (and "bank") it is a row range. Default: full extent.
        bank_size: int, column width of an ePix100 bank (default 128). Only
            used when axis="bank".
    """
    detector_key = kwargs.get("on")
    axis = kwargs.get("axis", "row")
    method = kwargs.get("method", "median")
    reference = kwargs.get("reference", None)
    bank_size = int(kwargs.get("bank_size", 128))
    if detector_key is None:
        return

    data = getattr(run, detector_key, None)
    if data is None:
        run.update_status(f"common_mode_correction: {detector_key} not found")
        return

    data = np.asarray(data, dtype=np.float64)
    if data.ndim == 2:
        # promote to (1, rows, cols) so one code path handles both
        working = data[np.newaxis, ...]
        squeezed = True
    elif data.ndim == 3:
        working = data
        squeezed = False
    else:
        run.update_status(
            f"common_mode_correction: {detector_key} must be 2D or 3D, got {data.ndim}D"
        )
        return

    reducer = np.nanmedian if method == "median" else np.nanmean

    # working is (shots, rows, cols). The reference range slices the axis
    # orthogonal to the correction direction so the offset is estimated only
    # from the dark band.
    if axis == "row":
        # per-row offset shared across columns; reference is a column range
        if reference is not None:
            band = working[:, :, reference[0] : reference[1]]
        else:
            band = working
        offset = reducer(band, axis=2, keepdims=True)  # (shots, rows, 1)
        corrected = working - offset
    elif axis == "column":
        # per-column offset shared across rows; reference is a row range
        if reference is not None:
            band = working[:, reference[0] : reference[1], :]
        else:
            band = working
        offset = reducer(band, axis=1, keepdims=True)  # (shots, 1, cols)
        corrected = working - offset
    elif axis == "bank":
        # per-bank offset: split columns into fixed-width blocks, estimate one
        # offset per block from the reference row band, subtract per block.
        if reference is not None:
            band = working[:, reference[0] : reference[1], :]
        else:
            band = working
        n_cols = working.shape[2]
        corrected = working.copy()
        for start in range(0, n_cols, bank_size):
            end = min(start + bank_size, n_cols)
            offset = reducer(band[:, :, start:end], axis=(1, 2), keepdims=True)
            corrected[:, :, start:end] = working[:, :, start:end] - offset
    else:
        run.update_status(
            f"common_mode_correction: unknown axis '{axis}' (use row|column|bank)"
        )
        return

    if squeezed:
        corrected = corrected[0]
    setattr(run, detector_key, corrected)
    run.update_status(
        f"Common-mode corrected {detector_key} (axis={axis}, method={method})"
    )

filter_detector_adu(run, **kwargs)

Zero out detector pixels below ADU threshold.

Parameters from YAML: on: detector key adu_threshold: float or [min, max] (default 3.0)

Source code in XSpect/analysis/spectroscopy.py

@register_step("filter_detector_adu")
def filter_detector_adu(run, **kwargs):
    """Zero out detector pixels below ADU threshold.

    Parameters from YAML:
        on: detector key
        adu_threshold: float or [min, max] (default 3.0)
    """
    detector_key = kwargs.get("on")
    adu_threshold = kwargs.get("adu_threshold", 3.0)
    if detector_key is None:
        return

    detector_images = getattr(run, detector_key, None)
    if detector_images is None:
        return
    if isinstance(adu_threshold, list):
        filtered = detector_images * (detector_images > adu_threshold[0])
        filtered = filtered * (filtered < adu_threshold[1])
    else:
        filtered = detector_images * (detector_images > adu_threshold)

    setattr(run, detector_key, filtered)
    run.update_status(f"ADU filtered {detector_key} with threshold {adu_threshold}")

filter_detector_variance(run, **kwargs)

Zero out low-variance detector pixels using sklearn VarianceThreshold.

Data-driven alternative to filter_detector_adu: instead of a hand-tuned intensity cutoff, pixels whose value barely changes across shots (dead pixels, static hot pixels, constant background) are dropped. Signal-bearing pixels vary shot to shot and are retained.

Parameters from YAML: on: detector key variance_threshold: float (default 0.0). Pixels with variance <= this across shots are zeroed. 0.0 removes only constant pixels.

Detector data is treated as (shots, features): a 3D array (shots x rows x cols) is flattened per shot, filtered, then reshaped back. Writes the filtered array back to detector_key and stores the retained boolean mask as _variance_mask.

Source code in XSpect/analysis/spectroscopy.py

@register_step("filter_detector_variance")
def filter_detector_variance(run, **kwargs):
    """Zero out low-variance detector pixels using sklearn VarianceThreshold.

    Data-driven alternative to filter_detector_adu: instead of a hand-tuned
    intensity cutoff, pixels whose value barely changes across shots (dead
    pixels, static hot pixels, constant background) are dropped. Signal-bearing
    pixels vary shot to shot and are retained.

    Parameters from YAML:
        on: detector key
        variance_threshold: float (default 0.0). Pixels with variance <= this
            across shots are zeroed. 0.0 removes only constant pixels.

    Detector data is treated as (shots, features): a 3D array
    (shots x rows x cols) is flattened per shot, filtered, then reshaped back.
    Writes the filtered array back to detector_key and stores the retained
    boolean mask as <detector_key>_variance_mask.
    """
    from sklearn.feature_selection import VarianceThreshold

    detector_key = kwargs.get("on")
    variance_threshold = kwargs.get("variance_threshold", 0.0)
    if detector_key is None:
        return

    detector_images = getattr(run, detector_key, None)
    if detector_images is None:
        return

    detector_images = np.asarray(detector_images)
    original_shape = detector_images.shape
    n_shots = original_shape[0]
    flat = detector_images.reshape(n_shots, -1)

    selector = VarianceThreshold(threshold=variance_threshold)
    selector.fit(flat)
    keep_mask = selector.get_support()

    filtered = flat * keep_mask  # broadcast over shots, zero the dropped pixels
    setattr(run, detector_key, filtered.reshape(original_shape))
    setattr(run, f"{detector_key}_variance_mask", keep_mask.reshape(original_shape[1:]))

    n_removed = int(np.sum(~keep_mask))
    run.update_status(
        f"Variance filtered {detector_key} (threshold={variance_threshold}): "
        f"{n_removed}/{keep_mask.size} pixels zeroed"
    )

filter_shots(run, **kwargs)

Filter a shot mask by thresholding on another key.

Parameters from YAML: on: shot mask key (e.g., "xray", "simultaneous") filter_key: key to threshold on (e.g., "ipm") threshold: float or [min, max]

Source code in XSpect/analysis/spectroscopy.py

@register_step("filter_shots")
def filter_shots(run, **kwargs):
    """Filter a shot mask by thresholding on another key.

    Parameters from YAML:
        on: shot mask key (e.g., "xray", "simultaneous")
        filter_key: key to threshold on (e.g., "ipm")
        threshold: float or [min, max]
    """
    shot_mask_key = kwargs.get("on")
    filter_key = kwargs.get("filter_key", "ipm")
    threshold = kwargs.get("threshold", 1.0e4)
    if shot_mask_key is None:
        return

    shot_mask = getattr(run, shot_mask_key)
    count_before = np.sum(shot_mask)
    filter_data = getattr(run, filter_key)
    nan_mask = np.isnan(filter_data)

    if isinstance(threshold, (int, float)):
        filtered = shot_mask * (filter_data > threshold) * (~nan_mask)
    elif len(threshold) == 2:
        filtered = (
            shot_mask
            * (filter_data > threshold[0])
            * (filter_data < threshold[1])
            * (~nan_mask)
        )
    else:
        filtered = shot_mask

    setattr(run, shot_mask_key, filtered)
    count_after = np.sum(filtered)
    run.update_status(
        f"Filtered {shot_mask_key} on {filter_key}: {int(count_before - count_after)} shots removed"
    )

find_rotation_angle(run, **kwargs)

Auto-detect the tilt angle of a dispersed spectral signal.

Delegates to XSpectDetectorProcessor which uses: 1. Canny edge detection on the 8-bit normalised sum image 2. DBSCAN clustering to isolate distinct signal streaks 3. PCA per cluster to find each streak's orientation 4. Mean angle across clusters, clamped to [-90, 90]

The result (degrees) is stored as run.<angle_key> for use by rotate_detector.

Parameters from YAML: on: detector key (3D shots x rows x cols, or 2D rows x cols) low_threshold: Canny lower hysteresis threshold (default 30) high_threshold: Canny upper hysteresis threshold (default 100) angle_key: attribute name to store the result (default: _angle)

Source code in XSpect/analysis/spectroscopy.py

@register_step("find_rotation_angle")
def find_rotation_angle(run, **kwargs):
    """Auto-detect the tilt angle of a dispersed spectral signal.

    Delegates to XSpectDetectorProcessor which uses:
      1. Canny edge detection on the 8-bit normalised sum image
      2. DBSCAN clustering to isolate distinct signal streaks
      3. PCA per cluster to find each streak's orientation
      4. Mean angle across clusters, clamped to [-90, 90]

    The result (degrees) is stored as ``run.<angle_key>`` for use by
    rotate_detector.

    Parameters from YAML:
        on: detector key (3D shots x rows x cols, or 2D rows x cols)
        low_threshold:  Canny lower hysteresis threshold (default 30)
        high_threshold: Canny upper hysteresis threshold (default 100)
        angle_key: attribute name to store the result (default: <on>_angle)
    """
    from XSpect.XSpect_Processor.XSpectDetectorProcessor import XSpectDetectorProcessor

    detector_key = kwargs.get("on")
    low_threshold = kwargs.get("low_threshold", 30)
    high_threshold = kwargs.get("high_threshold", 100)
    angle_key = kwargs.get("angle_key", f"{detector_key}_angle")

    if detector_key is None:
        return
    data = getattr(run, detector_key, None)
    if data is None:
        return

    # Build 2D sum image, clamp negatives
    img = data.sum(axis=0) if data.ndim == 3 else data.copy()
    img = np.clip(img, 0, None).astype(float)

    try:
        proc = XSpectDetectorProcessor(img)
        proc.detect_edges(low_threshold=low_threshold, high_threshold=high_threshold)
        angle = proc.find_optimal_rotation_angle()
    except Exception as e:
        run.update_status(f"find_rotation_angle: failed for {detector_key}: {e}")
        return

    setattr(run, angle_key, float(angle))
    run.update_status(
        f"find_rotation_angle: {detector_key} -> {angle_key} = {angle:.3f} deg"
    )

get_run_shot_properties(run, **kwargs)

Load xray/laser/simultaneous boolean masks from lightStatus.

Source code in XSpect/analysis/spectroscopy.py

@register_step("get_run_shot_properties")
def get_run_shot_properties(run, **kwargs):
    """Load xray/laser/simultaneous boolean masks from lightStatus."""
    run.get_run_shot_properties()

hitfinding(run, **kwargs)

Filter shots by total signal intensity (hit detection).

Keeps shots whose per-shot detector sum exceeds a threshold. Two modes (applied after any ADU filtering upstream):

min_sum (preferred for sparse/XES data): Keep shots where sum(detector) > min_sum. Use min_sum=1.0 to reject shots that are completely zero after the ADU threshold — i.e. true dark shots with no photons.

cutoff_multiplier (legacy, relative threshold): threshold = median(sums) - cutoff_multiplier * std(sums) Works well when signal shots are the majority; breaks down when most shots are dark (median ≈ 0 → threshold is negative).

If both are supplied, min_sum takes precedence.

Parameters from YAML: on: detector key (3D: shots × rows × cols) min_sum: absolute ADU floor per shot (default None) cutoff_multiplier: std multiplier for relative threshold (default 1.0)

Source code in XSpect/analysis/spectroscopy.py

@register_step("hitfinding")
def hitfinding(run, **kwargs):
    """Filter shots by total signal intensity (hit detection).

    Keeps shots whose per-shot detector sum exceeds a threshold.
    Two modes (applied after any ADU filtering upstream):

      min_sum (preferred for sparse/XES data):
        Keep shots where sum(detector) > min_sum.
        Use min_sum=1.0 to reject shots that are completely zero after
        the ADU threshold — i.e. true dark shots with no photons.

      cutoff_multiplier (legacy, relative threshold):
        threshold = median(sums) - cutoff_multiplier * std(sums)
        Works well when signal shots are the majority; breaks down when
        most shots are dark (median ≈ 0 → threshold is negative).

    If both are supplied, min_sum takes precedence.

    Parameters from YAML:
        on:                 detector key (3D: shots × rows × cols)
        min_sum:            absolute ADU floor per shot (default None)
        cutoff_multiplier:  std multiplier for relative threshold (default 1.0)
    """
    detector_key = kwargs.get("on")
    min_sum = kwargs.get("min_sum", None)
    cutoff_multiplier = kwargs.get("cutoff_multiplier", 1.0)
    if detector_key is None:
        return

    detector = getattr(run, detector_key, None)
    if detector is None or detector.ndim < 3:
        return

    sum_images = np.sum(detector, axis=(1, 2))

    if min_sum is not None:
        threshold = float(min_sum)
    else:
        threshold = np.median(sum_images) - cutoff_multiplier * np.std(sum_images)

    hits = np.where(sum_images > threshold)[0]

    if len(hits) == 0:
        setattr(
            run, detector_key, np.zeros((0,) + detector.shape[1:], dtype=detector.dtype)
        )
    elif len(hits) < detector.shape[0]:
        setattr(run, detector_key, detector[hits])
    run.update_status(
        f"Hitfinding on {detector_key}: {len(hits)}/{detector.shape[0]} shots kept "
        f"(threshold={threshold:.1f}, mode={'min_sum' if min_sum is not None else 'relative'})"
    )

load_detector(run, **kwargs)

Load 3D detector data (shots x rows x cols) with optional ROI and transpose.

Parameters from YAML: keys: list of HDF5 paths friendly_names: list of attribute names transpose: bool (default False) rois: list of [start, end] pixel ranges (default None) combine_rois: bool (default True)

Source code in XSpect/analysis/spectroscopy.py

@register_step("load_detector")
def load_detector(run, **kwargs):
    """Load 3D detector data (shots x rows x cols) with optional ROI and transpose.

    Parameters from YAML:
        keys: list of HDF5 paths
        friendly_names: list of attribute names
        transpose: bool (default False)
        rois: list of [start, end] pixel ranges (default None)
        combine_rois: bool (default True)
    """
    keys = kwargs.get("keys", [])
    friendly_names = kwargs.get("friendly_names", [])
    transpose = kwargs.get("transpose", False)
    rois = kwargs.get("rois", None)
    combine = kwargs.get("combine_rois", True)
    if not keys:
        return
    run.load_run_key_delayed(
        keys, friendly_names, transpose=transpose, rois=rois, combine=combine
    )

load_run_keys(run, **kwargs)

Load scalar/1D keys from HDF5 into run attributes.

Reads the keys defined in the pipeline's data section. Expects run to have run_file, start_index, end_index set.

Source code in XSpect/analysis/spectroscopy.py

@register_step("load_run_keys")
def load_run_keys(run, **kwargs):
    """Load scalar/1D keys from HDF5 into run attributes.

    Reads the keys defined in the pipeline's data section.
    Expects run to have run_file, start_index, end_index set.
    """
    keys = kwargs.get("keys", [])
    friendly_names = kwargs.get("friendly_names", [])
    if not keys:
        return
    run.load_run_keys(keys, friendly_names)

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")

Source code in XSpect/analysis/spectroscopy.py

@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}%)"
    )

purge_keys(run, **kwargs)

Delete specified attributes from run to free memory.

Parameters from YAML: keys: list of attribute names to purge

Source code in XSpect/analysis/spectroscopy.py

@register_step("purge_keys")
def purge_keys(run, **kwargs):
    """Delete specified attributes from run to free memory.

    Parameters from YAML:
        keys: list of attribute names to purge
    """
    keys = kwargs.get("keys", [])
    for key in keys:
        if hasattr(run, key):
            setattr(run, key, None)
    run.update_status(f"Purged keys: {keys}")

reduce_detector_shots(run, **kwargs)

Collapse the shot dimension using sum/mean.

Parameters from YAML: on: detector key reduction: "sum" or "mean" (default "sum") purge: bool (default True)

Source code in XSpect/analysis/spectroscopy.py

@register_step("reduce_detector_shots")
def reduce_detector_shots(run, **kwargs):
    """Collapse the shot dimension using sum/mean.

    Parameters from YAML:
        on: detector key
        reduction: "sum" or "mean" (default "sum")
        purge: bool (default True)
    """
    detector_key = kwargs.get("on")
    reduction_name = kwargs.get("reduction", "sum")
    purge = kwargs.get("purge", True)
    if detector_key is None:
        return

    reduction_fn = np.nansum if reduction_name == "sum" else np.nanmean
    detector = getattr(run, detector_key, None)
    if detector is None:
        return

    reduced = reduction_fn(detector, axis=0)
    setattr(run, f"{detector_key}_reduced", reduced)
    if purge:
        setattr(run, detector_key, None)
    run.update_status(f"Shot reduction: {detector_key} -> {detector_key}_reduced")

reduce_detector_spatial(run, **kwargs)

Reduce spatial dimension of detector using ROIs.

For 3D data (shots x rows x cols), the ROI selects along one spatial axis (default: axis 1 = cross-dispersion) and reduces by summing/averaging, producing (shots x remaining_spatial).

Parameters from YAML: on: detector key rois: list of [start, end] pixel ranges combine_rois: bool (default True) reduction: "sum" or "mean" (default "sum") axis: which spatial axis to apply ROI on (default 1 for 3D, -1 for 2D) purge: bool (default True) - delete original after reduction

Source code in XSpect/analysis/spectroscopy.py

@register_step("reduce_detector_spatial")
def reduce_detector_spatial(run, **kwargs):
    """Reduce spatial dimension of detector using ROIs.

    For 3D data (shots x rows x cols), the ROI selects along one spatial axis
    (default: axis 1 = cross-dispersion) and reduces by summing/averaging,
    producing (shots x remaining_spatial).

    Parameters from YAML:
        on: detector key
        rois: list of [start, end] pixel ranges
        combine_rois: bool (default True)
        reduction: "sum" or "mean" (default "sum")
        axis: which spatial axis to apply ROI on (default 1 for 3D, -1 for 2D)
        purge: bool (default True) - delete original after reduction
    """
    detector_key = kwargs.get("on")
    rois = kwargs.get("rois", [[0, None]])
    combine = kwargs.get("combine_rois", True)
    reduction_name = kwargs.get("reduction", "sum")
    roi_axis = kwargs.get("axis", None)
    purge = kwargs.get("purge", True)
    if detector_key is None:
        return

    reduction_fn = np.nansum if reduction_name == "sum" else np.nanmean
    detector = getattr(run, detector_key, None)
    if detector is None:
        return

    if roi_axis is None:
        roi_axis = 1 if detector.ndim == 3 else -1

    axis_size = detector.shape[roi_axis]

    # If the detector was cropped at import (row_range), the ROIs in the YAML
    # are in ABSOLUTE (full-frame) coordinates but the array now starts at the
    # crop origin. Translate ROIs into cropped-array coordinates automatically
    # so the YAML can always use absolute row numbers.
    # The offset is keyed by the loaded detector name (e.g. "epix"); derived
    # keys like "epix_reduced" inherit it by stripping known suffixes.
    row_offset = 0
    offsets = getattr(run, "_row_offset", None)
    if offsets:
        candidate = detector_key
        for suffix in ("_reduced", "_ROI_1", "_ROI_2"):
            candidate = candidate.replace(suffix, "")
        row_offset = offsets.get(detector_key, offsets.get(candidate, 0))
    if row_offset:
        adjusted = []
        for roi in rois:
            start = roi[0] - row_offset
            end = (roi[1] - row_offset) if roi[1] is not None else None
            adjusted.append([max(0, start), end])
        run.update_status(
            f"reduce_detector_spatial: applied row_range offset {row_offset} "
            f"to ROIs {rois} -> {adjusted}"
        )
        rois = adjusted

    if combine:
        mask = np.zeros(axis_size, dtype=bool)
        for roi in rois:
            end = roi[1] if roi[1] is not None else axis_size
            mask[roi[0] : end] = True

        idx = [slice(None)] * detector.ndim
        idx[roi_axis] = mask
        masked_data = detector[tuple(idx)]

        reduced = reduction_fn(masked_data, axis=roi_axis)
        setattr(run, f"{detector_key}_ROI_1", reduced)
        run.update_status(f"Spatial reduction: {detector_key} -> {detector_key}_ROI_1")
    else:
        for roi_idx, roi in enumerate(rois):
            end = roi[1] if roi[1] is not None else axis_size
            idx = [slice(None)] * detector.ndim
            idx[roi_axis] = slice(roi[0], end)
            chunk = detector[tuple(idx)]
            reduced = reduction_fn(chunk, axis=roi_axis)
            setattr(run, f"{detector_key}_ROI_{roi_idx + 1}", reduced)
            run.update_status(
                f"Spatial reduction: {detector_key} -> {detector_key}_ROI_{roi_idx + 1}"
            )

    if purge:
        setattr(run, detector_key, None)

reduce_detector_temporal(run, **kwargs)

Bin detector data along time dimension using timing indices.

Parameters from YAML: on: detector key (2D: shots x pixels) timing_bin_key: attribute name for bin indices (default "timing_bin_indices") average: bool (default False) - if True, divide by bin count

Source code in XSpect/analysis/spectroscopy.py

@register_step("reduce_detector_temporal")
def reduce_detector_temporal(run, **kwargs):
    """Bin detector data along time dimension using timing indices.

    Parameters from YAML:
        on: detector key (2D: shots x pixels)
        timing_bin_key: attribute name for bin indices (default "timing_bin_indices")
        average: bool (default False) - if True, divide by bin count
    """
    detector_key = kwargs.get("on")
    timing_bin_key = kwargs.get("timing_bin_key", "timing_bin_indices")
    average = kwargs.get("average", False)
    if detector_key is None:
        return

    detector = getattr(run, detector_key, None)
    timing_indices = getattr(run, timing_bin_key, None)
    time_bins = getattr(run, "time_bins", None)

    if detector is None or timing_indices is None or time_bins is None:
        run.update_status(f"reduce_detector_temporal: missing data for {detector_key}")
        return

    n_bins = len(time_bins)

    if detector.ndim == 1:
        binned = np.zeros(n_bins)
        bincount = np.zeros(n_bins)
        for i in range(len(detector)):
            idx = timing_indices[i] - 1
            if 0 <= idx < n_bins and not np.isnan(detector[i]):
                binned[idx] += detector[i]
                bincount[idx] += 1
    elif detector.ndim == 2:
        n_pixels = detector.shape[1]
        binned = np.zeros((n_bins, n_pixels))
        bincount = np.zeros(n_bins)
        for i in range(detector.shape[0]):
            idx = timing_indices[i] - 1
            if 0 <= idx < n_bins:
                row = detector[i]
                if not np.all(np.isnan(row)):
                    binned[idx] += np.where(np.isnan(row), 0.0, row)
                    bincount[idx] += 1
    else:
        run.update_status(f"reduce_detector_temporal: unsupported ndim={detector.ndim}")
        return

    if average:
        safe_count = np.where(bincount > 0, bincount, 1)
        if binned.ndim == 2:
            binned = binned / safe_count[:, np.newaxis]
        else:
            binned = binned / safe_count

    setattr(run, f"{detector_key}_time_binned", binned)
    setattr(run, f"{detector_key}_bincount", bincount)
    run.update_status(
        f"Temporal reduction: {detector_key} -> {detector_key}_time_binned ({n_bins} bins)"
    )

rotate_detector(run, **kwargs)

Rotate detector images by a fixed or auto-detected angle (scipy.ndimage.rotate).

Applied per-shot. For 3D data (shots x rows x cols), rotates in the (rows, cols) plane. For 2D data (rows x cols), rotates directly. Skipped when the resolved angle is 0.

Parameters from YAML: on: detector key angle: rotation in degrees (positive = CCW). Mutually exclusive with angle_key. angle_key: name of a run attribute holding the angle (set by find_rotation_angle). Takes precedence over angle if both are provided. axes: rotation plane as [ax1, ax2] (default [1, 2] for 3D, [0, 1] for 2D) reshape: bool (default False) — keep original array shape after rotation

Source code in XSpect/analysis/spectroscopy.py

@register_step("rotate_detector")
def rotate_detector(run, **kwargs):
    """Rotate detector images by a fixed or auto-detected angle (scipy.ndimage.rotate).

    Applied per-shot. For 3D data (shots x rows x cols), rotates in the
    (rows, cols) plane. For 2D data (rows x cols), rotates directly.
    Skipped when the resolved angle is 0.

    Parameters from YAML:
        on: detector key
        angle: rotation in degrees (positive = CCW). Mutually exclusive with angle_key.
        angle_key: name of a run attribute holding the angle (set by find_rotation_angle).
                   Takes precedence over angle if both are provided.
        axes: rotation plane as [ax1, ax2] (default [1, 2] for 3D, [0, 1] for 2D)
        reshape: bool (default False) — keep original array shape after rotation
    """
    detector_key = kwargs.get("on")
    angle_key = kwargs.get("angle_key", None)
    axes = kwargs.get("axes", None)
    if detector_key is None:
        return

    # Resolve angle: dynamic key takes precedence over static value
    if angle_key is not None:
        angle = getattr(run, angle_key, None)
        if angle is None:
            run.update_status(
                f"rotate_detector: angle_key '{angle_key}' not found, skipping"
            )
            return
    else:
        angle = kwargs.get("angle", 0)

    if angle == 0:
        return

    data = getattr(run, detector_key, None)
    if data is None:
        return

    if axes is None:
        axes = [0, 1] if data.ndim == 2 else [1, 2]

    # Static analysis path (2D) uses axes=[0,1]; per-shot 3D uses [1,2]
    # The old XESBatchAnalysisRotation uses [0,1] for the static case
    if data.ndim == 3 and axes == [0, 1]:
        axes = [1, 2]

    reshape = kwargs.get("reshape", False)
    rotated = rotate(data, angle=angle, axes=axes, reshape=reshape)
    setattr(run, detector_key, rotated)
    run.update_status(
        f"Rotated {detector_key} by {angle:.2f} degrees (axes={axes}, reshape={reshape})"
    )

separate_shots(run, **kwargs)

Extract shots matching first mask but NOT second (A and not B).

Parameters from YAML: on: detector/data key filter_keys: [include_mask, exclude_mask] new_key: optional output key name

Source code in XSpect/analysis/spectroscopy.py

@register_step("separate_shots")
def separate_shots(run, **kwargs):
    """Extract shots matching first mask but NOT second (A and not B).

    Parameters from YAML:
        on: detector/data key
        filter_keys: [include_mask, exclude_mask]
        new_key: optional output key name
    """
    detector_key = kwargs.get("on")
    filter_keys = kwargs.get("filter_keys", [])
    new_key = kwargs.get("new_key", None)
    if detector_key is None or len(filter_keys) < 2:
        return

    detector = getattr(run, detector_key)
    include_mask = getattr(run, filter_keys[0]).astype(bool)
    exclude_mask = getattr(run, filter_keys[1]).astype(bool)
    separation_mask = include_mask & (~exclude_mask)

    filtered_data = detector[separation_mask]

    if new_key is None:
        new_key = f"{detector_key}_{filter_keys[0]}_not_{filter_keys[1]}"
    setattr(run, new_key, filtered_data)
    run.update_status(
        f"Separated shots: {detector_key} ({filter_keys[0]} not {filter_keys[1]}) -> {new_key} ({int(np.sum(separation_mask))} shots)"
    )

subtract_polynomial_background(run, **kwargs)

Subtract a polynomial baseline fit along the spatial axis.

Fits a low-order polynomial to signal-free regions along one axis and subtracts it, removing smooth scattering/fluorescence background while preserving peak area. The peak region is excluded from the fit either by naming the background regions explicitly (background) or by masking the peak (peak_mask). Non-destructive: writes <on>_bkgsub.

Works on a 1D spectrum or a 2D array (bins x pixels). The fit reuses the vectorized weighted-polynomial projection from patch_pixels: the offset vector depends only on the sample positions and weights, so it is built once and applied to every row with a single matrix multiply.

Parameters from YAML: on: spectrum key (1D pixels, or 2D bins x pixels). axis: spatial axis to fit along (default: last axis). order: polynomial degree (default 2). background: list of [start, end] pixel ranges to fit (signal-free). If given, only these ranges anchor the fit. peak_mask: pixel range(s) to EXCLUDE from the fit (the emission peaks). Accepts a single [start, end] or a list of [start, end] ranges, so several dispersed lines (e.g. Kalpha and Kbeta on one detector) can all be masked at once. Used when naming the peaks is easier than the background. Ignored if background is given.

Source code in XSpect/analysis/spectroscopy.py

@register_step("subtract_polynomial_background")
def subtract_polynomial_background(run, **kwargs):
    """Subtract a polynomial baseline fit along the spatial axis.

    Fits a low-order polynomial to signal-free regions along one axis and
    subtracts it, removing smooth scattering/fluorescence background while
    preserving peak area. The peak region is excluded from the fit either by
    naming the background regions explicitly (``background``) or by masking
    the peak (``peak_mask``). Non-destructive: writes ``<on>_bkgsub``.

    Works on a 1D spectrum or a 2D array (bins x pixels). The fit reuses the
    vectorized weighted-polynomial projection from patch_pixels: the offset
    vector depends only on the sample positions and weights, so it is built
    once and applied to every row with a single matrix multiply.

    Parameters from YAML:
        on: spectrum key (1D pixels, or 2D bins x pixels).
        axis: spatial axis to fit along (default: last axis).
        order: polynomial degree (default 2).
        background: list of [start, end] pixel ranges to fit (signal-free).
            If given, only these ranges anchor the fit.
        peak_mask: pixel range(s) to EXCLUDE from the fit (the emission peaks).
            Accepts a single [start, end] or a list of [start, end] ranges, so
            several dispersed lines (e.g. Kalpha and Kbeta on one detector) can
            all be masked at once. Used when naming the peaks is easier than the
            background. Ignored if ``background`` is given.
    """
    detector_key = kwargs.get("on")
    axis = kwargs.get("axis", None)
    order = int(kwargs.get("order", 2))
    background = kwargs.get("background", None)
    peak_mask = kwargs.get("peak_mask", None)
    if detector_key is None:
        return

    data = getattr(run, detector_key, None)
    if data is None:
        run.update_status(f"subtract_polynomial_background: {detector_key} not found")
        return

    data = np.asarray(data, dtype=np.float64)
    if axis is None:
        axis = data.ndim - 1

    n_pixels = data.shape[axis]
    x = np.arange(n_pixels, dtype=np.float64)

    # weights select which pixels anchor the fit: 1 for background, 0 for peak.
    weights = np.ones(n_pixels, dtype=np.float64)
    if background is not None:
        weights[:] = 0.0
        for rng in background:
            weights[rng[0] : rng[1]] = 1.0
    elif peak_mask is not None:
        # accept a single [start, end] or a list of ranges (multiple lines)
        mask_ranges = peak_mask
        if len(peak_mask) == 2 and np.isscalar(peak_mask[0]):
            mask_ranges = [peak_mask]
        for rng in mask_ranges:
            weights[rng[0] : rng[1]] = 0.0

    # Bring the fit axis to front so every other axis is a batch dimension.
    moved = np.moveaxis(data, axis, 0)  # (n_pixels, ...)
    flat = moved.reshape(n_pixels, -1).copy()  # (n_pixels, N_rows)
    nan_mask = np.isnan(flat)
    flat[nan_mask] = 0.0

    if np.sum(weights > 0.5) < order + 1:
        run.update_status(
            f"subtract_polynomial_background: too few background pixels "
            f"({int(np.sum(weights > 0.5))}) for order {order}; skipped"
        )
        return

    # Per-row weights: the base background/peak mask, zeroed wherever a row has
    # a NaN so those points never enter that row's fit. NaN positions differ
    # per row, so the normal equations are solved per row (batched). The base
    # weight is squared to match numpy.polyfit's 1/sigma convention
    # (minimise sum(w**2 * r**2)).
    V = np.vander(x, order + 1)  # (n_pixels, order+1)
    w_full = (weights**2)[:, np.newaxis] * (~nan_mask)  # (n_pixels, N_rows)

    # A[r] = V^T diag(w_r) V ; b[r] = V^T diag(w_r) flat_r, batched over rows r.
    A = np.einsum("pk,pr,pl->rkl", V, w_full, V)  # (N_rows, order+1, order+1)
    b = np.einsum("pk,pr->rk", V, w_full * flat)  # (N_rows, order+1)
    coeffs = np.linalg.solve(A, b)  # (N_rows, order+1)
    baseline = (V @ coeffs.T)  # (n_pixels, N_rows)

    subtracted = flat - baseline
    subtracted[nan_mask] = np.nan  # keep original NaN positions
    result = np.moveaxis(subtracted.reshape(moved.shape), 0, axis)

    setattr(run, f"{detector_key}_bkgsub", result)
    run.update_status(
        f"Polynomial background subtracted {detector_key} -> "
        f"{detector_key}_bkgsub (order={order}, axis={axis})"
    )

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: "_bkgsub"); the original is kept estimator: "mean" or "median" over side-band pixels (default "median")

Source code in XSpect/analysis/spectroscopy.py

@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}"
    )

time_binning(run, **kwargs)

Create time delay bins from laser timing data.

Parameters from YAML: bins: "auto", list, or [min, max, num_points] lxt_key: str or null (default "lxt_ttc"; null means use encoder directly) fast_delay_key: str (default "encoder") tt_correction_key: str (default "time_tool_correction") resolution: bin width for auto mode (default 50e-15, i.e. 50 fs)

Source code in XSpect/analysis/spectroscopy.py

@register_step("time_binning")
def time_binning(run, **kwargs):
    """Create time delay bins from laser timing data.

    Parameters from YAML:
        bins: "auto", list, or [min, max, num_points]
        lxt_key: str or null (default "lxt_ttc"; null means use encoder directly)
        fast_delay_key: str (default "encoder")
        tt_correction_key: str (default "time_tool_correction")
        resolution: bin width for auto mode (default 50e-15, i.e. 50 fs)
    """
    bins_spec = kwargs.get("bins")
    lxt_key = kwargs.get("lxt_key", "lxt_ttc")
    fast_delay_key = kwargs.get("fast_delay_key", "encoder")
    tt_correction_key = kwargs.get("tt_correction_key", "time_tool_correction")
    resolution = kwargs.get("resolution", 50e-15)

    if bins_spec is None:
        return

    resolution = float(resolution)

    lxt_raw = getattr(run, lxt_key, None) if lxt_key else None
    fast_delay_raw = getattr(run, fast_delay_key, None) if fast_delay_key else None
    tt_correction_raw = (
        getattr(run, tt_correction_key, None) if tt_correction_key else None
    )

    lxt = np.asarray(lxt_raw, dtype=np.float64) if lxt_raw is not None else None
    fast_delay = (
        np.asarray(fast_delay_raw, dtype=np.float64)
        if fast_delay_raw is not None
        else None
    )
    tt_correction = (
        np.asarray(tt_correction_raw, dtype=np.float64)
        if tt_correction_raw is not None
        else None
    )

    if lxt is not None:
        delays = lxt.copy()
        if fast_delay is not None:
            if np.mean(np.abs(fast_delay)) > 1e-3:
                delays = delays + fast_delay
        if tt_correction is not None:
            delays = delays + tt_correction
    elif fast_delay is not None:
        delays = fast_delay.copy()
        if tt_correction is not None:
            delays = delays + tt_correction
    elif tt_correction is not None:
        delays = tt_correction.copy()
    else:
        run.update_status("time_binning: no timing keys found")
        return

    if bins_spec == "auto":
        rounded = np.round(delays / resolution) * resolution
        bins = np.unique(rounded)
    elif isinstance(bins_spec, list) and len(bins_spec) == 3:
        bins = np.linspace(bins_spec[0], bins_spec[1], int(bins_spec[2]))
    else:
        bins = np.array(bins_spec)

    bin_edges, _ = _center_binning(delays, bins)
    indices = np.digitize(delays, bin_edges)

    run.delays = delays
    run.time_bins = bins
    run.time_bins_centered = (bins[:-1] + bins[1:]) / 2 if len(bins) > 1 else bins
    run.timing_bin_indices = indices
    run.update_status(
        f"Time binning complete: {len(bins)} bins from {bins[0]:.4g} to {bins[-1]:.4g}"
    )

union_shots(run, **kwargs)

Combine shots matching multiple boolean masks (logical AND).

Parameters from YAML: on: detector/data key to filter filter_keys: list of mask attribute names to AND together new_key: optional output key name (default: auto-generated)

Source code in XSpect/analysis/spectroscopy.py

@register_step("union_shots")
def union_shots(run, **kwargs):
    """Combine shots matching multiple boolean masks (logical AND).

    Parameters from YAML:
        on: detector/data key to filter
        filter_keys: list of mask attribute names to AND together
        new_key: optional output key name (default: auto-generated)
    """
    detector_key = kwargs.get("on")
    filter_keys = kwargs.get("filter_keys", [])
    new_key = kwargs.get("new_key", None)
    if detector_key is None or not filter_keys:
        return

    detector = getattr(run, detector_key, None)
    if detector is None:
        return
    combined_mask = np.ones(detector.shape[0], dtype=bool)
    for fk in filter_keys:
        mask = getattr(run, fk)
        combined_mask = combined_mask * mask.astype(bool)

    if detector.ndim == 1:
        filtered_data = detector[combined_mask]
    elif detector.ndim == 2:
        filtered_data = detector[combined_mask]
    elif detector.ndim == 3:
        filtered_data = detector[combined_mask]
    else:
        filtered_data = detector[combined_mask]

    if new_key is None:
        new_key = f"{detector_key}_{'_'.join(filter_keys)}"
    setattr(run, new_key, filtered_data)
    run.update_status(
        f"Union shots: {detector_key} with masks {filter_keys} -> {new_key} ({int(np.sum(combined_mask))} shots)"
    )

XSpect.analysis.xes

XES (X-ray Emission Spectroscopy) specific pipeline steps.

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

Source code in XSpect/analysis/xes.py

@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"
    )

make_energy_axis(run, **kwargs)

Generate energy axis from vonHamos spectrometer geometry.

Uses the same formula as XSpect_Visualization.make_energy_axis: gl = pixel_indices * mm_per_pixel ll = gl/2 - (max(gl) - min(gl))/4 energy = 12398.42 / (2 * d * sin(arctan(R / (ll + A)))) energy = energy[::-1]

Parameters from YAML: detector_key: key to get pixel count from (or use n_pixels directly) n_pixels: number of pixels (overrides detector_key shape) crystal_detector_distance: A (mm) crystal_radius: R (mm) d_spacing: d (Angstrom) mm_per_pixel: pixel pitch (default 0.05 mm for ePix) name: output attribute name prefix (default: "xes")

Source code in XSpect/analysis/xes.py

@register_step("make_energy_axis")
def make_energy_axis(run, **kwargs):
    """Generate energy axis from vonHamos spectrometer geometry.

    Uses the same formula as XSpect_Visualization.make_energy_axis:
        gl = pixel_indices * mm_per_pixel
        ll = gl/2 - (max(gl) - min(gl))/4
        energy = 12398.42 / (2 * d * sin(arctan(R / (ll + A))))
        energy = energy[::-1]

    Parameters from YAML:
        detector_key: key to get pixel count from (or use n_pixels directly)
        n_pixels: number of pixels (overrides detector_key shape)
        crystal_detector_distance: A (mm)
        crystal_radius: R (mm)
        d_spacing: d (Angstrom)
        mm_per_pixel: pixel pitch (default 0.05 mm for ePix)
        name: output attribute name prefix (default: "xes")
    """
    detector_key = kwargs.get("detector_key", None)
    n_pixels = kwargs.get("n_pixels", None)
    A = kwargs.get("crystal_detector_distance")
    R = kwargs.get("crystal_radius")
    d = kwargs.get("d_spacing")
    mm_per_pixel = kwargs.get("mm_per_pixel", 0.05)
    name = kwargs.get("name", "xes")

    if A is None or R is None or d is None:
        run.update_status("make_energy_axis: missing geometry parameters")
        return

    if n_pixels is None and detector_key is not None:
        det = getattr(run, detector_key, None)
        if det is not None:
            n_pixels = det.shape[-1] if det.ndim >= 1 else 100
        else:
            n_pixels = 100

    if n_pixels is None:
        n_pixels = 100

    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
    energy_axis = hc / (2.0 * d * np.sin(np.arctan(R / (ll + A))))

    setattr(run, f"{name}_energy", energy_axis)
    run.update_status(
        f"Energy axis: {name}_energy ({n_pixels} pixels, center {energy_axis[n_pixels // 2]:.1f} eV)"
    )

normalize_xes(run, **kwargs)

Normalize XES spectra by row sum over a pixel range.

Divides each time bin's spectrum by its total intensity over the specified pixel range, producing area-normalized spectra.

Parameters from YAML: on: detector key (time_binned data, shape: bins x pixels) pixel_range: [start, end] pixels for normalization (default: full range)

Source code in XSpect/analysis/xes.py

@register_step("normalize_xes")
def normalize_xes(run, **kwargs):
    """Normalize XES spectra by row sum over a pixel range.

    Divides each time bin's spectrum by its total intensity over the
    specified pixel range, producing area-normalized spectra.

    Parameters from YAML:
        on: detector key (time_binned data, shape: bins x pixels)
        pixel_range: [start, end] pixels for normalization (default: full range)
    """
    detector_key = kwargs.get("on")
    pixel_range = kwargs.get("pixel_range", None)
    if detector_key is None:
        return

    data = getattr(run, detector_key, None)
    if data is None:
        run.update_status(f"normalize_xes: {detector_key} not found")
        return

    if data.ndim == 2:
        if pixel_range:
            norm_slice = data[:, pixel_range[0] : pixel_range[1]]
        else:
            norm_slice = data

        row_sums = np.sum(norm_slice, axis=1)
        safe_sums = np.where(row_sums > 0, row_sums, 1.0)
        normalized = data / safe_sums[:, np.newaxis]

        setattr(run, f"{detector_key}_normalized", normalized)

        # Propagate std if it exists
        std_key = detector_key.replace("_time_binned", "_std")
        std_data = getattr(run, std_key, None)
        if std_data is not None:
            std_normalized = std_data / safe_sums[:, np.newaxis]
            setattr(run, f"{detector_key}_normalized_std", std_normalized)

    elif data.ndim == 1:
        total = np.sum(data)
        if total > 0:
            normalized = data / total
        else:
            normalized = data
        setattr(run, f"{detector_key}_normalized", normalized)

    run.update_status(f"XES normalization: {detector_key} -> {detector_key}_normalized")

patch_pixels(run, **kwargs)

Repair bad pixels using polynomial fitting from neighbors.

Replicates XSpect_Analysis.patch_pixel with mode='polynomial': fits a weighted polynomial to surrounding pixels (excluding the bad pixel region) and evaluates at the bad pixel location.

Supports automatic detection of ASIC panel-gap spikes via auto_detect: true. Uses a two-method approach:

Method 1 — Global ratio: the column profile (sum over shots and rows) is compared to a median-filtered baseline. Columns with ratio > threshold (bright spike) or ratio < 1/threshold (dead gap within signal region) are flagged. Catches extreme outliers.

Method 2 — Z-score + width filter: computes a robust z-score (residual / MAD-based local sigma) for each column. Bright and dark outliers (|z| > nsigma) are separately filtered to keep only narrow clusters (≤ max_gap_width columns). This catches subtle ASIC gap columns (partially reduced intensity) that the ratio method misses — the ASIC gap pattern is a dark stripe (the actual gap) flanked by bright charge-sharing neighbors. Filtering bright and dark separately avoids merging them into one wide cluster.

The final set is the union of both methods, merged with any manually specified pixels.

Parameters from YAML: on: detector key pixels: list of pixel indices to patch (merged with auto if both given) auto_detect: bool (default False) — auto-find spike/dead-gap columns threshold: float (default 5.0) — ratio threshold for Method 1 (spike: ratio > threshold; dead: ratio < 1/threshold) nsigma: float (default 5.0) — z-score threshold for Method 2 max_gap_width: int (default 4) — maximum cluster width to keep in Method 2; wider clusters are rejected as signal gradients smooth_window: int (default 31) — median-filter window for baseline; must be >= the width of the widest expected spike mode: "polynomial", "interpolate", or "zero" (default "polynomial") axis: which axis the pixel indices refer to (default: 0 for 2D, last axis for 3D) patch_range: pixels on each side of bad pixel to exclude (default 4) poly_range: additional pixels beyond patch_range for fitting (default 6) deg: polynomial degree (default 1)

Source code in XSpect/analysis/xes.py

@register_step("patch_pixels")
def patch_pixels(run, **kwargs):
    """Repair bad pixels using polynomial fitting from neighbors.

    Replicates XSpect_Analysis.patch_pixel with mode='polynomial':
    fits a weighted polynomial to surrounding pixels (excluding the bad
    pixel region) and evaluates at the bad pixel location.

    Supports automatic detection of ASIC panel-gap spikes via
    ``auto_detect: true``.  Uses a two-method approach:

    **Method 1 — Global ratio:** the column profile (sum over shots and
    rows) is compared to a median-filtered baseline.  Columns with
    ratio > ``threshold`` (bright spike) or ratio < 1/threshold (dead
    gap within signal region) are flagged.  Catches extreme outliers.

    **Method 2 — Z-score + width filter:** computes a robust z-score
    (residual / MAD-based local sigma) for each column.  Bright and dark
    outliers (|z| > ``nsigma``) are separately filtered to keep only
    narrow clusters (≤ ``max_gap_width`` columns).  This catches subtle
    ASIC gap columns (partially reduced intensity) that the ratio method
    misses — the ASIC gap pattern is a dark stripe (the actual gap)
    flanked by bright charge-sharing neighbors.  Filtering bright and
    dark separately avoids merging them into one wide cluster.

    The final set is the union of both methods, merged with any manually
    specified ``pixels``.

    Parameters from YAML:
        on:            detector key
        pixels:        list of pixel indices to patch (merged with auto if both given)
        auto_detect:   bool (default False) — auto-find spike/dead-gap columns
        threshold:     float (default 5.0) — ratio threshold for Method 1
                       (spike: ratio > threshold; dead: ratio < 1/threshold)
        nsigma:        float (default 5.0) — z-score threshold for Method 2
        max_gap_width: int (default 4) — maximum cluster width to keep in
                       Method 2; wider clusters are rejected as signal gradients
        smooth_window: int (default 31) — median-filter window for baseline;
                       must be >= the width of the widest expected spike
        mode:          "polynomial", "interpolate", or "zero" (default "polynomial")
        axis:          which axis the pixel indices refer to
                       (default: 0 for 2D, last axis for 3D)
        patch_range:   pixels on each side of bad pixel to exclude (default 4)
        poly_range:    additional pixels beyond patch_range for fitting (default 6)
        deg:           polynomial degree (default 1)
    """
    from scipy.ndimage import median_filter as _mf

    detector_key = kwargs.get("on")
    pixels = list(kwargs.get("pixels", []))
    mode = kwargs.get("mode", "polynomial")
    axis = kwargs.get("axis", None)
    patch_range = kwargs.get("patch_range", 4)
    poly_range = kwargs.get("poly_range", 6)
    deg = kwargs.get("deg", 1)
    auto_detect = kwargs.get("auto_detect", False)

    if detector_key is None:
        return

    data = getattr(run, detector_key, None)
    if data is None:
        return

    if axis is None:
        # Default: pixel indices refer to columns (last axis)
        # 1D: axis=0; 2D (rows×cols): axis=1; 3D (shots×rows×cols): axis=2
        axis = data.ndim - 1

    # ------------------------------------------------------------------ #
    # Auto-detection: find ASIC panel-gap spikes / dead-gap columns       #
    # ------------------------------------------------------------------ #
    manual_pixels = list(pixels)  # preserve user-specified pixels
    if auto_detect:
        from scipy.ndimage import label as _label

        smooth_win = int(kwargs.get("smooth_window", 31))
        threshold = float(kwargs.get("threshold", 5.0))
        nsigma = float(kwargs.get("nsigma", 5.0))
        max_gap_width = int(kwargs.get("max_gap_width", 4))

        # Collapse shot dimension if 3D → 2D (rows × cols)
        if data.ndim == 3:
            img = np.clip(data, 0, None).sum(axis=0).astype(float)
        elif data.ndim == 2:
            img = np.clip(data, 0, None).astype(float)
        else:
            img = None

        if img is not None and img.ndim == 2:
            col_profile = img.sum(axis=0)
        else:
            col_profile = np.clip(data, 0, None).astype(float)

        n_cols = len(col_profile)
        baseline = _mf(col_profile, size=smooth_win, mode="nearest")

        with np.errstate(divide="ignore", invalid="ignore"):
            ratio = np.where(baseline > 0, col_profile / baseline, 0.0)
        in_signal = baseline > baseline.max() * 0.05

        # ── Method 1: Global ratio (extreme outliers) ────────────────
        ratio_spikes = set(np.where(ratio > threshold)[0])
        ratio_dead = set(np.where((ratio < 1.0 / threshold) & in_signal)[0])

        # ── Method 2: Z-score + width filter (subtle ASIC gaps) ──────
        # The column profile z-score catches columns whose intensity
        # deviates significantly from the local baseline.  Width filtering
        # then keeps only narrow defects (≤ max_gap_width cols), rejecting
        # broad signal gradients.  Bright and dark outliers are filtered
        # separately because ASIC gaps produce a dark stripe flanked by
        # bright charge-sharing neighbors; together they form a wide
        # cluster, but independently each is narrow.
        residual = col_profile - baseline
        abs_residual = np.abs(residual)
        local_mad = _mf(abs_residual, size=smooth_win * 2 + 1, mode="nearest")
        local_sigma = np.maximum(local_mad * 1.4826, baseline * 0.01)

        with np.errstate(divide="ignore", invalid="ignore"):
            z_score = np.where(local_sigma > 0, residual / local_sigma, 0.0)

        def _narrow_clusters(mask, max_w):
            """Keep only contiguous flagged regions ≤ max_w wide."""
            labeled_arr, n = _label(mask)
            out = np.zeros_like(mask)
            for i in range(1, n + 1):
                cluster = np.where(labeled_arr == i)[0]
                if (cluster.max() - cluster.min() + 1) <= max_w:
                    out[cluster] = True
            return out

        bright_mask = (z_score > nsigma) & in_signal
        dark_mask = (z_score < -nsigma) & in_signal

        zscore_spikes = set(np.where(_narrow_clusters(bright_mask, max_gap_width))[0])
        zscore_dead = set(np.where(_narrow_clusters(dark_mask, max_gap_width))[0])

        # ── Combine (union of both methods) ──────────────────────────
        spike_cols = sorted(ratio_spikes | zscore_spikes)
        dead_cols = sorted(ratio_dead | zscore_dead)
        auto_pixels = sorted(set(spike_cols + dead_cols))

        # Merge manual pixels with auto-detected ones
        pixels = sorted(set(manual_pixels + auto_pixels))

        # Store diagnostics on the run object for inspection
        setattr(run, f"{detector_key}_auto_patched_pixels", auto_pixels)
        setattr(run, f"{detector_key}_manual_patched_pixels", manual_pixels)
        setattr(run, f"{detector_key}_auto_spike_cols", spike_cols)
        setattr(run, f"{detector_key}_auto_dead_cols", dead_cols)
        setattr(run, f"{detector_key}_col_profile", col_profile)
        setattr(run, f"{detector_key}_col_baseline", baseline)
        setattr(run, f"{detector_key}_col_zscore", z_score)
        run.update_status(
            f"patch_pixels auto_detect: "
            f"{len(spike_cols)} spikes + {len(dead_cols)} dead = "
            f"{len(auto_pixels)} auto + {len(manual_pixels)} manual = "
            f"{len(pixels)} total on {detector_key} "
            f"(threshold={threshold}, nsigma={nsigma}, "
            f"max_gap_width={max_gap_width})"
        )

    if not pixels:
        return

    n_pixels = data.shape[axis]
    bad_set = set(pixels)  # all detected bad columns (for exclusion)

    for pixel in pixels:
        if pixel < 0 or pixel >= n_pixels:
            continue

        if mode == "zero":
            slc = [slice(None)] * data.ndim
            slc[axis] = pixel
            data[tuple(slc)] = 0
        elif mode in ("polynomial", "interpolate"):
            start = pixel - patch_range - poly_range
            end = pixel + patch_range + poly_range + 1
            actual_start = max(start, 0)
            actual_end = min(end, n_pixels)
            slc_region = [slice(None)] * data.ndim
            slc_region[axis] = slice(actual_start, actual_end)
            region = np.moveaxis(data[tuple(slc_region)], axis, 0)

            patch_x = np.arange(actual_start, actual_end)
            weights = np.ones(len(patch_x))

            # Build exclusion mask: zero-weight the target pixel's
            # neighborhood AND the ±patch_range neighborhood of every
            # other bad pixel in the region.  This ensures the fit
            # anchors only on data well outside the entire ASIC feature
            # (gap + charge-sharing halo), preventing overshoot from
            # elevated neighbors.
            for bad_col in bad_set:
                idx_in_region = bad_col - actual_start
                lo = max(idx_in_region - patch_range, 0)
                hi = min(idx_in_region + patch_range + 1, len(weights))
                if lo < len(weights) and hi > 0:
                    weights[lo:hi] = 0.0

            # If too few non-zero weights remain, expand the fit region
            good_count = np.sum(weights > 0.5)
            if good_count < 4:
                # Expand region symmetrically until we get enough clean pixels
                expand = poly_range
                while good_count < 4 and expand < n_pixels // 2:
                    expand += poly_range
                    exp_start = max(pixel - patch_range - expand, 0)
                    exp_end = min(pixel + patch_range + expand + 1, n_pixels)
                    slc_region = [slice(None)] * data.ndim
                    slc_region[axis] = slice(exp_start, exp_end)
                    region = np.moveaxis(data[tuple(slc_region)], axis, 0)
                    patch_x = np.arange(exp_start, exp_end)
                    weights = np.ones(len(patch_x))
                    # Zero-weight ±patch_range around ALL bad pixels
                    for bad_col in bad_set:
                        idx_in_region = bad_col - exp_start
                        lo = max(idx_in_region - patch_range, 0)
                        hi = min(idx_in_region + patch_range + 1, len(weights))
                        if lo < len(weights) and hi > 0:
                            weights[lo:hi] = 0.0
                    good_count = np.sum(weights > 0.5)
                    actual_start, actual_end = exp_start, exp_end

            if data.ndim == 1:
                coeffs = np.polyfit(patch_x, region, deg, w=weights)
                new_val = np.polyval(coeffs, pixel)
            else:
                # Vectorized weighted polynomial fit across all (shot, row)
                # pairs simultaneously.  The projection vector `proj` maps
                # the region slice to the interpolated value at `pixel` —
                # it only depends on patch_x, weights, and deg (same for
                # all columns of `flat`), so we compute it once and apply
                # via a single matrix-vector multiply.
                other_shape = region.shape[1:]
                flat = region.reshape(region.shape[0], -1).astype(
                    np.float64
                )  # (n_region, N_total)

                # Build Vandermonde matrix and solve for projection vector.
                # Square the weights to match numpy.polyfit convention
                # (polyfit treats w as 1/sigma, minimising sum(w**2 * r**2)).
                X = np.vander(patch_x.astype(np.float64), deg + 1)
                W = (weights**2).astype(np.float64)
                XtW = X.T * W[np.newaxis, :]  # (deg+1, n_region)
                XtWX = XtW @ X  # (deg+1, deg+1)
                XtWX_inv_XtW = np.linalg.solve(XtWX, XtW)  # (deg+1, n_region)
                x_eval = np.vander(np.array([float(pixel)]), deg + 1)[0]  # (deg+1,)
                proj = x_eval @ XtWX_inv_XtW  # (n_region,)

                # Single dot product replaces the per-column polyfit loop
                new_vals = proj @ flat  # (N_total,)
                new_val = new_vals.reshape(other_shape)

            slc = [slice(None)] * data.ndim
            slc[axis] = pixel
            data[tuple(slc)] = new_val

    setattr(run, detector_key, data)
    run.update_status(
        f"Patched {len(pixels)} pixels on {detector_key} (axis={axis}, mode={mode})"
    )

XSpect.analysis.xas

XAS (X-ray Absorption Spectroscopy) specific pipeline steps.

ccm_binning(run, **kwargs)

Digitize CCM data into energy bin indices.

Parameters from YAML: ccm_key: attribute with CCM energy values (default "ccm") ccm_bins_key: attribute with bin edges (default "ccm_bins")

Source code in XSpect/analysis/xas.py

@register_step("ccm_binning")
def ccm_binning(run, **kwargs):
    """Digitize CCM data into energy bin indices.

    Parameters from YAML:
        ccm_key: attribute with CCM energy values (default "ccm")
        ccm_bins_key: attribute with bin edges (default "ccm_bins")
    """
    ccm_key = kwargs.get("ccm_key", "ccm")
    ccm_bins_key = kwargs.get("ccm_bins_key", "ccm_bins")

    ccm_data = getattr(run, ccm_key, None)
    ccm_bins = getattr(run, ccm_bins_key, None)

    if ccm_data is None or ccm_bins is None:
        run.update_status(f"ccm_binning: missing {ccm_key} or {ccm_bins_key}")
        return

    indices = np.digitize(ccm_data, ccm_bins)
    run.ccm_bin_indices = indices
    run.update_status(f"CCM binning complete: {len(ccm_bins)} bins")

make_ccm_axis(run, **kwargs)

Generate CCM (Channel Cut Monochromator) energy bins.

Parameters from YAML: energies: "auto", explicit list, or [min, max, num_points] ccm_key: attribute with CCM values (for auto mode, default "ccm") resolution: rounding resolution in same units as ccm data (for auto mode, default 0.001)

Source code in XSpect/analysis/xas.py

@register_step("make_ccm_axis")
def make_ccm_axis(run, **kwargs):
    """Generate CCM (Channel Cut Monochromator) energy bins.

    Parameters from YAML:
        energies: "auto", explicit list, or [min, max, num_points]
        ccm_key: attribute with CCM values (for auto mode, default "ccm")
        resolution: rounding resolution in same units as ccm data (for auto mode, default 0.001)
    """
    energies_spec = kwargs.get("energies")
    ccm_key = kwargs.get("ccm_key", "ccm")
    resolution = kwargs.get("resolution", 0.001)
    if energies_spec is None:
        return

    # If ccm_bins was pre-computed on the parent run and injected into this
    # batch, skip re-derivation so all batches share the same global axis.
    if getattr(run, "ccm_bins", None) is not None:
        return

    if energies_spec == "auto":
        ccm_data = getattr(run, ccm_key, None)
        if ccm_data is None:
            run.update_status(f"make_ccm_axis: missing {ccm_key} for auto binning")
            return
        rounded = np.round(ccm_data / resolution) * resolution
        energies = np.unique(rounded)
    elif isinstance(energies_spec, list) and len(energies_spec) == 3:
        energies = np.linspace(
            energies_spec[0], energies_spec[1], int(energies_spec[2])
        )
    else:
        energies = np.array(energies_spec)

    if len(energies) < 2:
        addon = resolution / 2
    else:
        addon = (energies[-1] - energies[-2]) / 2
    bins2 = np.append(energies, energies[-1] + addon)
    bins_center = np.empty_like(bins2)
    for ii in range(len(energies)):
        if ii == 0:
            bins_center[ii] = bins2[ii] - (bins2[ii + 1] - bins2[ii]) / 2
        else:
            bins_center[ii] = bins2[ii] - (bins2[ii] - bins2[ii - 1]) / 2
    bins_center[-1] = bins2[-1]

    run.ccm_bins = bins_center
    run.ccm_energies = energies
    run.update_status(
        f"CCM axis: {len(energies)} energy points from {energies[0]:.1f} to {energies[-1]:.1f}"
    )

reduce_detector_ccm(run, **kwargs)

Bin detector data along energy (CCM) dimension.

Parameters from YAML: on: detector key (1D or 2D) ccm_bin_key: bin indices attribute (default "ccm_bin_indices") average: bool (default False)

Source code in XSpect/analysis/xas.py

@register_step("reduce_detector_ccm")
def reduce_detector_ccm(run, **kwargs):
    """Bin detector data along energy (CCM) dimension.

    Parameters from YAML:
        on: detector key (1D or 2D)
        ccm_bin_key: bin indices attribute (default "ccm_bin_indices")
        average: bool (default False)
    """
    detector_key = kwargs.get("on")
    ccm_bin_key = kwargs.get("ccm_bin_key", "ccm_bin_indices")
    average = kwargs.get("average", False)
    if detector_key is None:
        return

    detector = getattr(run, detector_key, None)
    ccm_indices = getattr(run, ccm_bin_key, None)
    ccm_bins = getattr(run, "ccm_bins", None)
    ccm_energies = getattr(run, "ccm_energies", None)

    if detector is None or ccm_indices is None or ccm_bins is None:
        run.update_status(f"reduce_detector_ccm: missing data for {detector_key}")
        return

    # ccm_bins has n_energies+1 elements (shifted edge array from make_ccm_axis).
    # Use n_energies so the output axis aligns 1-to-1 with ccm_energies.
    n_bins = len(ccm_energies) if ccm_energies is not None else len(ccm_bins)

    if detector.ndim == 1:
        binned = np.zeros(n_bins)
        bincount = np.zeros(n_bins)
        for i in range(len(detector)):
            idx = ccm_indices[i] - 1
            if 0 <= idx < n_bins and not np.isnan(detector[i]):
                binned[idx] += detector[i]
                bincount[idx] += 1
    elif detector.ndim == 2:
        n_pixels = detector.shape[1]
        binned = np.zeros((n_bins, n_pixels))
        bincount = np.zeros(n_bins)
        for i in range(detector.shape[0]):
            idx = ccm_indices[i] - 1
            if 0 <= idx < n_bins:
                row = detector[i]
                if not np.all(np.isnan(row)):
                    binned[idx] += np.where(np.isnan(row), 0.0, row)
                    bincount[idx] += 1
    elif detector.ndim == 3:
        n_rows = detector.shape[1]
        n_cols = detector.shape[2]
        binned = np.zeros((n_bins, n_rows, n_cols))
        bincount = np.zeros(n_bins)
        for i in range(detector.shape[0]):
            idx = ccm_indices[i] - 1
            if 0 <= idx < n_bins:
                frame = detector[i]
                if not np.all(np.isnan(frame)):
                    binned[idx] += np.where(np.isnan(frame), 0.0, frame)
                    bincount[idx] += 1
    else:
        run.update_status(f"reduce_detector_ccm: unsupported ndim={detector.ndim}")
        return

    if average:
        safe_count = np.where(bincount > 0, bincount, 1)
        if binned.ndim == 2:
            binned = binned / safe_count[:, np.newaxis]
        else:
            binned = binned / safe_count

    setattr(run, f"{detector_key}_energy_binned", binned)
    setattr(run, f"{detector_key}_energy_bincount", bincount)
    run.update_status(
        f"CCM reduction: {detector_key} -> {detector_key}_energy_binned ({n_bins} bins)"
    )

reduce_detector_ccm_temporal(run, **kwargs)

2D binning: both time AND energy (CCM) dimensions.

Produces a 3D array: (time_bins x energy_bins x pixels) or 2D: (time_bins x energy_bins) for scalar detectors.

Parameters from YAML: on: detector key timing_bin_key: time bin indices (default "timing_bin_indices") ccm_bin_key: CCM bin indices (default "ccm_bin_indices") average: bool (default False)

Source code in XSpect/analysis/xas.py

@register_step("reduce_detector_ccm_temporal")
def reduce_detector_ccm_temporal(run, **kwargs):
    """2D binning: both time AND energy (CCM) dimensions.

    Produces a 3D array: (time_bins x energy_bins x pixels) or
    2D: (time_bins x energy_bins) for scalar detectors.

    Parameters from YAML:
        on: detector key
        timing_bin_key: time bin indices (default "timing_bin_indices")
        ccm_bin_key: CCM bin indices (default "ccm_bin_indices")
        average: bool (default False)
    """
    detector_key = kwargs.get("on")
    timing_bin_key = kwargs.get("timing_bin_key", "timing_bin_indices")
    ccm_bin_key = kwargs.get("ccm_bin_key", "ccm_bin_indices")
    average = kwargs.get("average", False)
    if detector_key is None:
        return

    detector = getattr(run, detector_key, None)
    timing_indices = getattr(run, timing_bin_key, None)
    ccm_indices = getattr(run, ccm_bin_key, None)
    time_bins = getattr(run, "time_bins", None)
    ccm_bins = getattr(run, "ccm_bins", None)

    if any(
        x is None for x in [detector, timing_indices, ccm_indices, time_bins, ccm_bins]
    ):
        run.update_status(
            f"reduce_detector_ccm_temporal: missing data for {detector_key}"
        )
        return

    n_time = len(time_bins)
    # Use ccm_energies length (n bins) not ccm_bins length (n+1 edges)
    ccm_energies = getattr(run, "ccm_energies", None)
    n_energy = len(ccm_energies) if ccm_energies is not None else len(ccm_bins) - 1

    if detector.ndim == 1:
        binned = np.zeros((n_time, n_energy))
        bincount = np.zeros((n_time, n_energy))
        for i in range(len(detector)):
            t_idx = timing_indices[i] - 1
            e_idx = ccm_indices[i] - 1
            if (
                0 <= t_idx < n_time
                and 0 <= e_idx < n_energy
                and not np.isnan(detector[i])
            ):
                binned[t_idx, e_idx] += detector[i]
                bincount[t_idx, e_idx] += 1
    elif detector.ndim == 2:
        n_pixels = detector.shape[1]
        binned = np.zeros((n_time, n_energy, n_pixels))
        bincount = np.zeros((n_time, n_energy))
        for i in range(detector.shape[0]):
            t_idx = timing_indices[i] - 1
            e_idx = ccm_indices[i] - 1
            if 0 <= t_idx < n_time and 0 <= e_idx < n_energy:
                row = detector[i]
                if not np.all(np.isnan(row)):
                    binned[t_idx, e_idx] += np.where(np.isnan(row), 0.0, row)
                    bincount[t_idx, e_idx] += 1
    else:
        run.update_status(
            f"reduce_detector_ccm_temporal: unsupported ndim={detector.ndim}"
        )
        return

    if average:
        safe_count = np.where(bincount > 0, bincount, 1)
        if binned.ndim == 3:
            binned = binned / safe_count[:, :, np.newaxis]
        else:
            binned = binned / safe_count

    setattr(run, f"{detector_key}_time_energy_binned", binned)
    setattr(run, f"{detector_key}_time_energy_bincount", bincount)
    run.update_status(
        f"CCM+temporal reduction: {detector_key} -> {detector_key}_time_energy_binned ({n_time}x{n_energy})"
    )

XSpect.analysis.droplet

Droplet-to-photon reconstruction pipeline step.

Reads per-shot sparse photon coordinates from the fixed-length or variable-length droplet2photon arrays stored in smalldata HDF5 files and scatters them back onto dense 2D images, producing a (N_shots, rows, cols) stack for downstream steps.

Supported HDF5 layouts

Fixed-length (e.g. epix100_0): /droplet_droplet2phot_sparse_{row,col,data,tile} shape (Nshots, nData) — zero-padded; valid entries have data != 0

Variable-length (e.g. epix100_1): /var_droplet_droplet2phot_sparse/{row,col,data} flat arrays + /var_droplet_droplet2phot_sparse_len (per-shot counts)

Batch awareness

When the pipeline batches a run, the batch manager injects abs_start_index and abs_end_index onto each batch run so the step knows which HDF5 rows to read. In the non-batched path these attributes are absent and the step falls back to run.start_index / run.end_index.

ROI-direct scattering

When a roi is given, photons are filtered to the ROI window and scattered straight into the small cropped image (_scatter_roi) — the full 704x768 panel is never allocated. This is bit-identical to scattering onto the full panel and cropping, but uses ~30x less memory and runs ~2x faster for the typical (60, 300) XES window.

droplet_reconstruction(run, **kwargs)

Reconstruct per-shot 2D images from droplet2photon sparse data.

Reads sparse photon-position arrays directly from the source smalldata HDF5 file and scatters them onto dense images. The resulting (N_shots, rows, cols) array is stored on the run object under new_key and is compatible with all downstream steps (filter_detector_adu, patch_pixels, rotate_detector, reduce_detector_spatial, etc.).

The step reads from run.run_file. In batch mode the batch manager injects abs_start_index and abs_end_index onto the batch run so the correct HDF5 rows are read; in the non-batched path these attributes are absent and the step derives the range from run.start_index / run.end_index.

Parameters from YAML

det : str — detector group in the HDF5, e.g. "epix100_0" new_key : str — attribute name to store the reconstructed stack roi : list — [row0, row1, col0, col1] crop region (omit for full panel) panel_shape : list — [rows, cols] of the full panel (default [704, 768])

Example YAML step

• step: droplet_reconstruction det: epix100_0 new_key: epix_spec roi: [270, 330, 400, 700] # -> (60, 300) output matching ROI_area

Source code in XSpect/analysis/droplet.py

@register_step("droplet_reconstruction")
def droplet_reconstruction(run, **kwargs):
    """Reconstruct per-shot 2D images from droplet2photon sparse data.

    Reads sparse photon-position arrays directly from the source smalldata HDF5
    file and scatters them onto dense images.  The resulting (N_shots, rows, cols)
    array is stored on the run object under ``new_key`` and is compatible with
    all downstream steps (``filter_detector_adu``, ``patch_pixels``,
    ``rotate_detector``, ``reduce_detector_spatial``, etc.).

    The step reads from ``run.run_file``.  In batch mode the batch manager
    injects ``abs_start_index`` and ``abs_end_index`` onto the batch run so
    the correct HDF5 rows are read; in the non-batched path these attributes
    are absent and the step derives the range from ``run.start_index`` /
    ``run.end_index``.

    Parameters from YAML
    --------------------
    det         : str   — detector group in the HDF5, e.g. "epix100_0"
    new_key     : str   — attribute name to store the reconstructed stack
    roi         : list  — [row0, row1, col0, col1] crop region (omit for full panel)
    panel_shape : list  — [rows, cols] of the full panel (default [704, 768])

    Example YAML step
    -----------------
    - step: droplet_reconstruction
      det: epix100_0
      new_key: epix_spec
      roi: [270, 330, 400, 700]   # -> (60, 300) output matching ROI_area
    """
    det = kwargs.get("det")
    new_key = kwargs.get("new_key")
    roi_kwarg = kwargs.get("roi", None)
    panel_shape_kwarg = kwargs.get("panel_shape", None)

    if det is None or new_key is None:
        run.update_status("droplet_reconstruction: 'det' and 'new_key' are required")
        return

    run_file = getattr(run, "run_file", None)
    if run_file is None:
        run.update_status("droplet_reconstruction: run_file not set, skipping")
        return

    # ------------------------------------------------------------------
    # Resolve absolute HDF5 shot indices
    # In batch mode the batch manager sets abs_start_index / abs_end_index.
    # In the non-batched path those attrs are absent; fall back to the
    # run's own start_index / end_index.
    # ------------------------------------------------------------------
    abs_start = getattr(run, "abs_start_index", None)
    if abs_start is None:
        abs_start = getattr(run, "start_index", 0)

    abs_end = getattr(run, "abs_end_index", None)
    if abs_end is None:
        end_idx = getattr(run, "end_index", -1)
        if end_idx == -1:
            # Need to look up the total shot count from the HDF5
            try:
                with h5py.File(run_file, "r") as _fh:
                    fixed_key = f"{det}/droplet_droplet2phot_sparse_data"
                    var_key = f"{det}/var_droplet_droplet2phot_sparse_len"
                    if fixed_key in _fh:
                        abs_end = int(_fh[fixed_key].shape[0])
                    elif var_key in _fh:
                        abs_end = int(_fh[var_key].shape[0])
                    else:
                        # Fall back to total_shots if known
                        abs_end = abs_start + int(getattr(run, "total_shots", 0))
            except Exception as exc:
                run.update_status(
                    f"droplet_reconstruction: could not determine shot count: {exc}"
                )
                return
        else:
            abs_end = int(end_idx)

    if abs_end <= abs_start:
        run.update_status(
            f"droplet_reconstruction: empty range [{abs_start}, {abs_end}), skipping"
        )
        return

    panel_shape = (
        tuple(panel_shape_kwarg) if panel_shape_kwarg is not None else PANEL_SHAPE
    )
    roi = tuple(roi_kwarg) if roi_kwarg is not None else None

    n_shots = abs_end - abs_start

    t0 = time.time()
    try:
        with h5py.File(run_file, "r") as fh:
            fixed_key = f"{det}/droplet_droplet2phot_sparse_data"
            var_key = f"{det}/var_droplet_droplet2phot_sparse/data"

            if fixed_key in fh:
                images = _reconstruct_fixed(
                    fh, det, abs_start, abs_end, roi, panel_shape
                )
            elif var_key in fh:
                images = _reconstruct_var(fh, det, abs_start, abs_end, roi, panel_shape)
            else:
                available = list(fh[det].keys()) if det in fh else []
                run.update_status(
                    f"droplet_reconstruction: no sparse data found for '{det}'. "
                    f"Available keys under '{det}': {available}"
                )
                return
    except Exception as exc:
        run.update_status(f"droplet_reconstruction: HDF5 read failed: {exc}")
        return

    elapsed = time.time() - t0
    setattr(run, new_key, images)
    run.update_status(
        f"droplet_reconstruction: '{det}' -> '{new_key}' "
        f"shape={images.shape} roi={roi} "
        f"shots=[{abs_start},{abs_end}) "
        f"elapsed={elapsed:.1f}s"
    )