fix(par-dsp): handle peaks with no candidate events in evtsel

src/legenddataflowscripts/par/geds/dsp/dplms.py

@@ -21,6 +21,7 @@
     expand_filelist,
     prepare_output_paths,
     require_config_keys,
+    require_peaks_present,
     take_table_rows,
 )
 
@@ -152,6 +153,11 @@ def par_geds_dsp_dplms() -> None:
 
         peaks_rounded = [int(peak) for peak in peaks_kev]
         peaks = lh5.read_as(f"{args.raw_table_name}/peak", args.peak_file, library="np")
+        require_peaks_present(
+            np.unique(peaks),
+            peaks_rounded,
+            f"peak file {args.peak_file} requested by dplms config ({args.config_file})",
+        )
         ids = np.isin(peaks, peaks_rounded)
 
         raw_cal = lh5.read(args.raw_table_name, args.peak_file, idx=ids)

src/legenddataflowscripts/par/geds/dsp/eopt.py

@@ -24,6 +24,7 @@
     check_input_files,
     prepare_output_paths,
     require_config_keys,
+    require_peaks_present,
 )
 
 warnings.filterwarnings(action="ignore", category=RuntimeWarning)
@@ -181,6 +182,11 @@ def par_geds_dsp_eopt() -> None:
 
         peaks_rounded = [int(peak) for peak in peaks_kev]
         peaks = lh5.read_as(f"{args.raw_table_name}/peak", args.peak_file, library="np")
+        require_peaks_present(
+            np.unique(peaks),
+            peaks_rounded,
+            f"peak file {args.peak_file} requested by eopt config ({args.config_file})",
+        )
         ids = np.isin(peaks, peaks_rounded)
         peaks = peaks[ids]
         idx_list = [np.where(peaks == peak)[0] for peak in peaks_rounded]

src/legenddataflowscripts/par/geds/dsp/evtsel.py

@@ -133,6 +133,52 @@ def get_out_data(
     return out_tbl, int(np.count_nonzero(final_mask))
 
 
+def build_peak_dicts(peaks_kev, kev_widths, masks, log=None):
+    """Build the per-peak accumulator state for the file x peak read loop.
+
+    Peaks for which *masks* holds no candidate event are dropped: with an
+    empty index array nothing is ever read into the peak's buffer, and the
+    end-of-file flush would then hand ``None`` to
+    :func:`dspeed.build_dsp`.  This happens for channels whose DAQ trigger
+    threshold sits above a configured gamma line.
+
+    Parameters
+    ----------
+    peaks_kev : sequence of float
+        Nominal gamma-line energies in keV.
+    kev_widths : sequence of (float, float)
+        Lower/upper selection widths in keV, one pair per peak.
+    masks : dict
+        Mapping of peak energy to the array of candidate event indices.
+    log : logging.Logger, optional
+        Logger used to warn about skipped peaks.
+
+    Returns
+    -------
+    dict
+        Mapping of peak energy to its accumulator state, with peaks that have
+        no candidate events omitted.
+    """
+    pk_dicts = {}
+    for peak, kev_width in zip(peaks_kev, kev_widths, strict=False):
+        if len(masks[peak]) == 0:
+            if log is not None:
+                msg = (
+                    f"no events in the rough energy range for {peak}, "
+                    "skipping this peak"
+                )
+                log.warning(msg)
+            continue
+        pk_dicts[peak] = {
+            "idxs": (masks[peak],),
+            "n_rows_read": 0,
+            "obj_buf_start": 0,
+            "obj_buf": None,
+            "kev_width": kev_width,
+        }
+    return pk_dicts
+
+
 def par_geds_dsp_evtsel() -> None:
     """Select calibration peak events for DSP parameter optimisation.
 
@@ -363,15 +409,7 @@ def par_geds_dsp_evtsel() -> None:
         # the whole file x peak loop below
         del input_data, tb_data
 
-        pk_dicts = {}
-        for peak, kev_width in zip(peaks_kev, kev_widths, strict=False):
-            pk_dicts[peak] = {
-                "idxs": (masks[peak],),
-                "n_rows_read": 0,
-                "obj_buf_start": 0,
-                "obj_buf": None,
-                "kev_width": kev_width,
-            }
+        pk_dicts = build_peak_dicts(peaks_kev, kev_widths, masks, log=log)
 
         for file in raw_files:
             log.debug(Path(file).name)
@@ -403,7 +441,13 @@ def par_geds_dsp_evtsel() -> None:
                         log.debug(msg)
                         peak_dict["obj_buf_start"] += n_rows_read_i
                     if peak_dict["n_rows_read"] >= 10000 or file == raw_files[-1]:
-                        if "e_lower_lim" not in peak_dict:
+                        if (
+                            peak_dict["obj_buf"] is None
+                            or len(peak_dict["obj_buf"]) == 0
+                        ):
+                            # nothing accumulated for this peak yet, nothing to flush
+                            pass
+                        elif "e_lower_lim" not in peak_dict:
                             tb_out = build_dsp(
                                 raw_in=peak_dict["obj_buf"],
                                 dsp_config=dsp_config,
@@ -501,10 +545,7 @@ def par_geds_dsp_evtsel() -> None:
                             peak_dict["n_events"] = n_wfs
                             msg = f"found {peak_dict['n_events']} events for {peak}"
                             log.debug(msg)
-                        elif (
-                            peak_dict["obj_buf"] is not None
-                            and len(peak_dict["obj_buf"]) > 0
-                        ):
+                        else:
                             tb_out = build_dsp(
                                 raw_in=peak_dict["obj_buf"],
                                 dsp_config=dsp_config,

src/legenddataflowscripts/utils/__init__.py

@@ -1,7 +1,12 @@
 from __future__ import annotations
 
 from .alias_table import alias_table
-from .cfgtools import get_channel_config, get_rule_config, require_config_keys
+from .cfgtools import (
+    get_channel_config,
+    get_rule_config,
+    require_config_keys,
+    require_peaks_present,
+)
 from .convert_np import convert_dict_np_to_float
 from .discharge import get_is_recovering_mask
 from .files import (
@@ -30,5 +35,6 @@
     "parse_json_arg",
     "prepare_output_paths",
     "require_config_keys",
+    "require_peaks_present",
     "take_table_rows",
 ]

src/legenddataflowscripts/utils/cfgtools.py

@@ -75,6 +75,33 @@ def require_config_keys(config: Mapping, keys: Iterable[str], context: str) -> N
         raise ValueError(msg)
 
 
+def require_peaks_present(
+    available: Iterable, required: Iterable, context: str
+) -> None:
+    """Raise ValueError listing all ``required`` peaks absent from ``available``.
+
+    The peak files written by ``par-geds-dsp-evtsel`` tag every row with the
+    nominal gamma-line energy in a ``peak`` column, so a peak for which no
+    events were selected is simply absent rather than empty.  Consumers filter
+    on that column and would otherwise operate silently on zero rows.
+
+    Parameters
+    ----------
+    available : iterable
+        Peak labels actually present in the peak file.
+    required : iterable
+        Peak labels the consumer's configuration asks for.
+    context : str
+        Free-form string naming the peak file and requesting config in the
+        error message.
+    """
+    present = set(available)
+    missing = sorted({peak for peak in required if peak not in present})
+    if missing:
+        msg = f"{context} is missing required peak(s) {missing}"
+        raise ValueError(msg)
+
+
 def get_rule_config(configs_path, rule_name, timestamp, datatype):
     """Resolve the dataflow config for one Snakemake rule.
 

tests/test_peak_selection.py

@@ -0,0 +1,83 @@
+"""Tests for the peak-selection helpers used by the dsp par scripts.
+
+Covers ``build_peak_dicts`` (which drops peaks with no candidate events, so a
+channel whose DAQ threshold sits above a configured gamma line does not hand
+``None`` to ``build_dsp``) and ``require_peaks_present`` (which makes the peak
+file consumers fail loudly instead of filtering down to zero rows).
+"""
+
+from __future__ import annotations
+
+import numpy as np
+import pytest
+
+from legenddataflowscripts.par.geds.dsp.evtsel import build_peak_dicts
+from legenddataflowscripts.utils import require_peaks_present
+
+PEAKS = [583.191, 727.33, 860.564, 1620.5, 2614.553]
+KEV_WIDTHS = [(25, 30), (25, 35), (25, 40), (15, 55), (70, 70)]
+
+
+def _masks(counts):
+    return {
+        peak: np.arange(count, dtype=int)
+        for peak, count in zip(PEAKS, counts, strict=True)
+    }
+
+
+def test_all_peaks_populated():
+    counts = [33123, 10172, 6538, 2519, 12999]
+    pk_dicts = build_peak_dicts(PEAKS, KEV_WIDTHS, _masks(counts))
+
+    assert list(pk_dicts) == PEAKS
+    for peak, kev_width, count in zip(PEAKS, KEV_WIDTHS, counts, strict=True):
+        entry = pk_dicts[peak]
+        assert entry["kev_width"] == kev_width
+        assert len(entry["idxs"][0]) == count
+        assert entry["n_rows_read"] == 0
+        assert entry["obj_buf_start"] == 0
+        assert entry["obj_buf"] is None
+
+
+def test_empty_peaks_are_skipped():
+    # V01404A in p18: the DAQ threshold sits above the three lowest lines
+    counts = [0, 0, 0, 7315, 12112]
+    pk_dicts = build_peak_dicts(PEAKS, KEV_WIDTHS, _masks(counts))
+
+    assert list(pk_dicts) == [1620.5, 2614.553]
+    assert len(pk_dicts[1620.5]["idxs"][0]) == 7315
+    assert pk_dicts[1620.5]["kev_width"] == (15, 55)
+
+
+def test_no_peaks_at_all():
+    assert build_peak_dicts(PEAKS, KEV_WIDTHS, _masks([0] * 5)) == {}
+
+
+def test_skipped_peaks_are_logged():
+    class _Recorder:
+        def __init__(self):
+            self.messages = []
+
+        def warning(self, msg):
+            self.messages.append(msg)
+
+    log = _Recorder()
+    build_peak_dicts(PEAKS, KEV_WIDTHS, _masks([0, 0, 0, 7315, 12112]), log=log)
+
+    assert len(log.messages) == 3
+    assert all("skipping this peak" in msg for msg in log.messages)
+    assert "583.191" in log.messages[0]
+
+
+def test_require_peaks_present_passes():
+    require_peaks_present(np.array([1620, 2614]), [2614], "peak file")
+
+
+def test_require_peaks_present_lists_all_missing():
+    with pytest.raises(ValueError, match=r"missing required peak\(s\) \[583, 727\]"):
+        require_peaks_present(np.array([1620, 2614]), [583, 727, 2614], "peak file")
+
+
+def test_require_peaks_present_names_the_context():
+    with pytest.raises(ValueError, match=r"my-peaks\.lh5"):
+        require_peaks_present(np.array([2614]), [583], "peak file my-peaks.lh5")