fix: keep analyzer crystals and inputs in geometry artifacts

src/ess/livedata/scripts/make_geometry_nexus.py

@@ -39,6 +39,17 @@ def _copy_source_fields(src: h5py.Group, dst: h5py.Group) -> None:
             src.copy(key, dst)
 
 
+def _copy_crystal_fields(src: h5py.Group, dst: h5py.Group) -> None:
+    """Copy an analyzer crystal's placement and lattice spacing.
+
+    Indirect-geometry instruments (BIFROST) derive each detector's final energy
+    from the position and ``d_spacing`` of the crystal its ``inputs`` lead to.
+    """
+    _copy_depends_on(src, dst)
+    if 'd_spacing' in src:
+        src.copy('d_spacing', dst)
+
+
 def _pixel_offsets_from_off(
     src_group: h5py.Group, active_face: str
 ) -> dict[str, tuple[np.ndarray, str]] | None:
@@ -215,29 +226,40 @@ def _visitor(name: str, obj: h5py.Group | h5py.Dataset) -> None:
     return targets
 
 
+def _create_group(fin: h5py.File, fout: h5py.File, path: str) -> h5py.Group:
+    """Create ``fout[path]`` with the attributes and ``inputs`` of ``fin[path]``.
+
+    ``inputs`` names the component(s) upstream in the beam path. BIFROST finds
+    each detector's analyzer crystal by following ``inputs`` from the detector,
+    possibly through intermediate groups. Files following the NeXus standard
+    store it as a dataset, older ESS files as an attribute.
+    """
+    src = fin[path]
+    dst = fout.create_group(path)
+    _copy_attributes(src, dst)
+    if 'inputs' in src:
+        src.copy('inputs', dst)
+    return dst
+
+
 def _ensure_parent_groups(fin: h5py.File, fout: h5py.File, path: str) -> None:
-    """Create parent groups in *fout*, copying attributes from *fin*."""
+    """Create parent groups in *fout* from their counterparts in *fin*."""
     parts = path.split('/')
     current = ''
     for part in parts[:-1]:
         if not part:
             continue
         current = f'{current}/{part}'
         if current not in fout:
-            fout.create_group(current)
-            _copy_attributes(fin[current], fout[current])
+            _create_group(fin, fout, current)
 
 
-def _get_or_create_dst(
-    fin: h5py.File, fout: h5py.File, name: str, src: h5py.Group
-) -> h5py.Group:
-    """Return ``fout[name]``, creating it (and parents) from ``src`` if absent."""
+def _get_or_create_dst(fin: h5py.File, fout: h5py.File, name: str) -> h5py.Group:
+    """Return ``fout[name]``, creating it (and parents) from ``fin`` if absent."""
     _ensure_parent_groups(fin, fout, name)
     if name in fout:
         return fout[name]
-    dst = fout.create_group(name)
-    _copy_attributes(src, dst)
-    return dst
+    return _create_group(fin, fout, name)
 
 
 def _resolve_depends_on_chains(fin: h5py.File, fout: h5py.File) -> set[str]:
@@ -284,6 +306,7 @@ def _resolve_depends_on_chains(fin: h5py.File, fout: h5py.File) -> set[str]:
     'NXsample',
     'NXtransformations',
     'NXdisk_chopper',
+    'NXcrystal',
 )
 
 
@@ -309,7 +332,7 @@ def visit_and_copy(name: str, obj: h5py.Group | h5py.Dataset) -> None:
             nx_class = _nx_class(obj)
             if nx_class not in _HANDLED_NX_CLASSES:
                 return
-            dst = _get_or_create_dst(fin, fout, name, obj)
+            dst = _get_or_create_dst(fin, fout, name)
             if nx_class == 'NXdetector':
                 _copy_detector_fields(
                     obj,
@@ -323,6 +346,8 @@ def visit_and_copy(name: str, obj: h5py.Group | h5py.Dataset) -> None:
                 _copy_source_fields(obj, dst)
             elif nx_class == 'NXsample':
                 _copy_depends_on(obj, dst)
+            elif nx_class == 'NXcrystal':
+                _copy_crystal_fields(obj, dst)
             else:  # NXtransformations or NXdisk_chopper
                 for key in obj:
                     _copy_child(obj, key, dst)

tests/make_geometry_nexus_test.py

@@ -4,6 +4,10 @@
 import h5py
 import numpy as np
 import pytest
+import scipp as sc
+from ess.bifrost.io.nexus import load_analyzer_for_detector
+from ess.reduce.nexus.types import NeXusComponentLocationSpec, SampleRun
+from scippnexus import NXdetector
 
 from ess.livedata.scripts.make_geometry_nexus import write_minimal_geometry
 
@@ -644,3 +648,102 @@ def test_skips_ess_writer_bookkeeping_groups(tmp_path):
     with h5py.File(output, 'r') as f:
         assert 'entry/instrument/monitor_1' in f
         assert '.ESS_monitor' not in f['entry']
+
+
+def _write_analyzer_chain(path, inputs_encoding: str) -> None:
+    """BIFROST-like detector whose ``inputs`` lead to its analyzer crystal.
+
+    The chain passes through an intermediate group that holds nothing but its
+    own placement, as ``channel_*_detector_angle`` does in BIFROST files.
+    """
+    with h5py.File(path, 'w') as f:
+        entry = f.create_group('entry')
+        entry.attrs['NX_class'] = 'NXentry'
+        inst = entry.create_group('instrument')
+        inst.attrs['NX_class'] = 'NXinstrument'
+
+        def add_inputs(group: h5py.Group, upstream: str) -> None:
+            if inputs_encoding == 'dataset':
+                group.create_dataset('inputs', data=f'/entry/instrument/{upstream}')
+            else:
+                group.attrs['inputs'] = upstream
+
+        def add_placement(group: h5py.Group) -> None:
+            group.create_dataset(
+                'depends_on', data=f'{group.name}/transformations/offset'
+            )
+            tr = group.create_group('transformations')
+            tr.attrs['NX_class'] = 'NXtransformations'
+            ds = tr.create_dataset('offset', data=1.0)
+            ds.attrs['transformation_type'] = 'translation'
+            ds.attrs['vector'] = [0.0, 0.0, 1.0]
+            ds.attrs['units'] = 'm'
+            ds.attrs['depends_on'] = '.'
+
+        det = inst.create_group('triplet')
+        det.attrs['NX_class'] = 'NXdetector'
+        det.create_dataset('detector_number', data=np.arange(3))
+        det.create_dataset('x_pixel_offset', data=np.zeros(3))
+        add_placement(det)
+        add_inputs(det, 'detector_angle')
+
+        angle = inst.create_group('detector_angle')
+        angle.attrs['NX_class'] = 'NXcoordinate_system'
+        angle.create_group('transformations').attrs['NX_class'] = 'NXtransformations'
+        add_inputs(angle, 'analyzer')
+
+        analyzer = inst.create_group('analyzer')
+        analyzer.attrs['NX_class'] = 'NXcrystal'
+        analyzer.create_dataset('d_spacing', data=3.355).attrs['units'] = 'angstrom'
+        add_placement(analyzer)
+        # Points at a component the artifact does not keep.
+        add_inputs(analyzer, 'collimator')
+
+
+@pytest.mark.parametrize('inputs_encoding', ['dataset', 'attribute'])
+def test_keeps_inputs_chain_to_analyzer_crystal(tmp_path, inputs_encoding):
+    """BIFROST finds each detector's analyzer by following ``inputs`` from the
+    detector to an ``NXcrystal``, whose position and ``d_spacing`` give the
+    final energy. The artifact must keep both the chain and the crystal.
+    """
+    src = tmp_path / 'input.nxs'
+    _write_analyzer_chain(src, inputs_encoding)
+    output = tmp_path / 'output.nxs'
+    write_minimal_geometry(src, output)
+
+    def inputs(group: h5py.Group) -> str:
+        if inputs_encoding == 'dataset':
+            return group['inputs'][()].decode()
+        return group.attrs['inputs']
+
+    with h5py.File(output, 'r') as f:
+        inst = f['entry/instrument']
+        assert inputs(inst['triplet']).endswith('detector_angle')
+        assert inputs(inst['detector_angle']).endswith('analyzer')
+        analyzer = inst['analyzer']
+        assert analyzer.attrs['NX_class'] == 'NXcrystal'
+        assert inputs(analyzer).endswith('collimator')
+        assert analyzer['d_spacing'][()] == pytest.approx(3.355)
+        assert analyzer['d_spacing'].attrs['units'] == 'angstrom'
+        assert analyzer['depends_on'][()].decode() == (
+            '/entry/instrument/analyzer/transformations/offset'
+        )
+        assert analyzer['transformations/offset'][()] == pytest.approx(1.0)
+
+
+def test_analyzer_lookup_finds_crystal_in_artifact(tmp_path):
+    # The attribute encoding is the one every essspectroscopy version reads.
+    src = tmp_path / 'input.nxs'
+    _write_analyzer_chain(src, 'attribute')
+    output = tmp_path / 'output.nxs'
+    write_minimal_geometry(src, output)
+
+    analyzer = load_analyzer_for_detector(
+        NeXusComponentLocationSpec[NXdetector, SampleRun](
+            filename=output, component_name='triplet'
+        )
+    )
+    assert sc.identical(analyzer['d_spacing'], sc.scalar(3.355, unit='angstrom'))
+    assert analyzer['depends_on'].value == (
+        '/entry/instrument/analyzer/transformations/offset'
+    )