feat: load 3D detector geometry

docs/user-guide/beer/beer_modulation_mcstas.ipynb

@@ -161,7 +161,7 @@
    "outputs": [],
    "source": [
     "wf = BeerModMcStasWorkflowKnownPeaks()\n",
-    "wf[DetectorBank] = 2\n",
+    "wf[DetectorBank] = DetectorBank.north\n",
     "wf[Filename[SampleRun]] = mcstas_silicon_medium_resolution()\n",
     "da = wf.compute(RawDetector[SampleRun])\n",
     "da.masks.clear()\n",
@@ -187,7 +187,7 @@
     "wf[Filename[SampleRun]] = mcstas_silicon_medium_resolution()\n",
     "\n",
     "results = {}\n",
-    "for bank in (1, 2):\n",
+    "for bank in DetectorBank:\n",
     "    wf[DetectorBank] = bank\n",
     "    da = wf.compute(TofDetector[SampleRun])\n",
     "    results[bank] = (\n",
@@ -218,7 +218,7 @@
     "wf[Filename[SampleRun]] = mcstas_silicon_medium_resolution()\n",
     "\n",
     "results = {}\n",
-    "for bank in (1, 2):\n",
+    "for bank in DetectorBank:\n",
     "    wf[DetectorBank] = bank\n",
     "    da = wf.compute(TofDetector[SampleRun])\n",
     "    results[bank] = (\n",
@@ -277,7 +277,7 @@
    "outputs": [],
    "source": [
     "wf = BeerModMcStasWorkflowKnownPeaks()\n",
-    "wf[DetectorBank] = 1\n",
+    "wf[DetectorBank] = DetectorBank.south\n",
     "wf[Filename[SampleRun]] = mcstas_duplex(8)\n",
     "wf.compute(RawDetector[SampleRun]).hist(two_theta=400, event_time_offset=1000).plot(norm='log', cmin=1.0e-2)"
    ]
@@ -300,7 +300,7 @@
     "wf[DHKLList] = duplex_peaks_array()\n",
     "\n",
     "results = {}\n",
-    "for bank in (1, 2):\n",
+    "for bank in DetectorBank:\n",
     "    wf[DetectorBank] = bank\n",
     "    da = wf.compute(TofDetector[SampleRun])\n",
     "    results[bank] = (\n",
@@ -331,7 +331,7 @@
     "wf[Filename[SampleRun]] = mcstas_duplex(8)\n",
     "\n",
     "results = {}\n",
-    "for bank in (1, 2):\n",
+    "for bank in DetectorBank:\n",
     "    wf[DetectorBank] = bank\n",
     "    da = wf.compute(TofDetector[SampleRun])\n",
     "    results[bank] = (\n",
@@ -390,7 +390,7 @@
    "outputs": [],
    "source": [
     "wf = BeerModMcStasWorkflowKnownPeaks()\n",
-    "wf[DetectorBank] = 1\n",
+    "wf[DetectorBank] = DetectorBank.south\n",
     "wf[Filename[SampleRun]] = mcstas_duplex(9)\n",
     "wf.compute(RawDetector[SampleRun]).hist(two_theta=400, event_time_offset=1000).plot(norm='log', cmin=1.0e-3)"
    ]
@@ -413,7 +413,7 @@
     "wf[DHKLList] = duplex_peaks_array()\n",
     "\n",
     "results = {}\n",
-    "for bank in (1, 2):\n",
+    "for bank in DetectorBank:\n",
     "    wf[DetectorBank] = bank\n",
     "    da = wf.compute(TofDetector[SampleRun])\n",
     "    results[bank] = (\n",
@@ -444,7 +444,7 @@
     "wf[Filename[SampleRun]] = mcstas_duplex(9)\n",
     "\n",
     "results = {}\n",
-    "for bank in (1, 2):\n",
+    "for bank in DetectorBank:\n",
     "    wf[DetectorBank] = bank\n",
     "    da = wf.compute(TofDetector[SampleRun])\n",
     "    results[bank] = (\n",
@@ -503,7 +503,7 @@
    "outputs": [],
    "source": [
     "wf = BeerModMcStasWorkflowKnownPeaks()\n",
-    "wf[DetectorBank] = 1\n",
+    "wf[DetectorBank] = DetectorBank.south\n",
     "wf[Filename[SampleRun]] = mcstas_duplex(10)\n",
     "wf.compute(RawDetector[SampleRun]).hist(two_theta=400, event_time_offset=1000).plot(norm='log', cmin=1.0e-3)"
    ]
@@ -526,7 +526,7 @@
     "wf[DHKLList] = duplex_peaks_array()\n",
     "\n",
     "results = {}\n",
-    "for bank in (1, 2):\n",
+    "for bank in DetectorBank:\n",
     "    wf[DetectorBank] = bank\n",
     "    da = wf.compute(TofDetector[SampleRun])\n",
     "    results[bank] = (\n",
@@ -557,7 +557,7 @@
     "wf[Filename[SampleRun]] = mcstas_duplex(10)\n",
     "\n",
     "results = {}\n",
-    "for bank in (1, 2):\n",
+    "for bank in DetectorBank:\n",
     "    wf[DetectorBank] = bank\n",
     "    da = wf.compute(TofDetector[SampleRun])\n",
     "    results[bank] = (\n",
@@ -616,7 +616,7 @@
    "outputs": [],
    "source": [
     "wf = BeerModMcStasWorkflowKnownPeaks()\n",
-    "wf[DetectorBank] = 1\n",
+    "wf[DetectorBank] = DetectorBank.south\n",
     "wf[Filename[SampleRun]] = mcstas_duplex(16)\n",
     "wf.compute(RawDetector[SampleRun]).hist(two_theta=400, event_time_offset=1000).plot(norm='log', cmin=1.0e-3)"
    ]
@@ -639,7 +639,7 @@
     "wf[DHKLList] = duplex_peaks_array()\n",
     "\n",
     "results = {}\n",
-    "for bank in (1, 2):\n",
+    "for bank in DetectorBank:\n",
     "    wf[DetectorBank] = bank\n",
     "    da = wf.compute(TofDetector[SampleRun])\n",
     "    results[bank] = (\n",
@@ -670,7 +670,7 @@
     "wf[Filename[SampleRun]] = mcstas_duplex(16)\n",
     "\n",
     "results = {}\n",
-    "for bank in (1, 2):\n",
+    "for bank in DetectorBank:\n",
     "    wf[DetectorBank] = bank\n",
     "    da = wf.compute(TofDetector[SampleRun])\n",
     "    results[bank] = (\n",

src/ess/beer/data.py

@@ -41,6 +41,9 @@
         "silicon-mode7-new-model.h5": "md5:d2070d3132722bb551d99b243c62752f",
         "silicon-mode8-new-model.h5": "md5:d6dfdf7e87eccedf4f83c67ec552ca22",
         "silicon-mode9-new-model.h5": "md5:694a17fb616b7f1c20e94d9da113d201",
+        # Simulation with 3D detector model - almost no events
+        # - only used to verify we can load the 3D geometry.
+        "few_neutrons_3d_detector_example.h5": "md5:88cbe29cb539c8acebf9fd7cee9d3c57",
     },
 )
 
@@ -74,6 +77,10 @@ def mcstas_silicon_new_model(mode: int) -> Path:
     return _registry.get_path(f'silicon-mode{mode}-new-model.h5')
 
 
+def mcstas_few_neutrons_3d_detector_example():
+    return _registry.get_path('few_neutrons_3d_detector_example.h5')
+
+
 def duplex_peaks() -> Path:
     return _registry.get_path('duplex-dhkl.tab')
 

src/ess/beer/io.py

@@ -39,7 +39,7 @@ def _find_h5(group: h5py.Group, matches):
         if re.match(matches, p):
             return group[p]
     else:
-        raise RuntimeError(f'Could not find "{matches}" in {group}.')
+        raise ValueError(f'Could not find "{matches}" in {group}.')
 
 
 def _load_h5(group: h5py.Group | str, *paths: str):
@@ -127,7 +127,7 @@ def _effective_chopper_position_from_mode(
         raise ValueError(f'Unkonwn chopper mode {mode}.')
 
 
-def _load_beer_mcstas(f, bank=1):
+def _load_beer_mcstas(f, north_or_south=None, *, number):
     positions = {
         name: f'/entry1/instrument/components/{key}/Position'
         for key in f['/entry1/instrument/components']
@@ -147,8 +147,16 @@ def _load_beer_mcstas(f, bank=1):
         mcc_pos,
     ) = _load_h5(
         f,
-        f'NXentry/NXdetector/bank{bank:02}_events_dat_list_p_x_y_n_id_t',
-        f'NXentry/NXdetector/bank{bank:02}_events_dat_list_p_x_y_n_id_t/events',
+        (
+            f'NXentry/NXdetector/bank_{north_or_south}{number}_events_dat_list_p_x_y_n_id_t'
+            if north_or_south is not None
+            else f'NXentry/NXdetector/bank{number:02}_events_dat_list_p_x_y_n_id_t'
+        ),
+        (
+            f'NXentry/NXdetector/bank_{north_or_south}{number}_events_dat_list_p_x_y_n_id_t/events'
+            if north_or_south is not None
+            else f'NXentry/NXdetector/bank{number:02}_events_dat_list_p_x_y_n_id_t/events'  # noqa: E501
+        ),
         'NXentry/simulation/Param',
         positions['sampleMantid'],
         positions['PSC1'],
@@ -162,7 +170,10 @@ def _load_beer_mcstas(f, bank=1):
         'Rotation'
     ]
     detector_rotation = _find_h5(
-        f['/entry1/instrument/components'], f'.*nD_Mantid_?{bank}.*'
+        f['/entry1/instrument/components'],
+        f'.*nD_Mantid_?{north_or_south}_{number}.*'
+        if north_or_south is not None
+        else f'.*nD_Mantid_?{number}.*',
     )['Rotation']
 
     events = events[()]
@@ -286,15 +297,74 @@ def _not_between(x, a, b):
     return (x < a) | (b < x)
 
 
-def load_beer_mcstas(f: str | Path | h5py.File, bank: int) -> sc.DataArray:
+def load_beer_mcstas(f: str | Path | h5py.File, bank: DetectorBank) -> sc.DataArray:
     '''Load beer McStas data from a file to a
     data group with one data array for each bank.
     '''
+    if not isinstance(bank, DetectorBank):
+        raise ValueError(
+            '"bank" must be either ``DetectorBank.north`` or ``DetectorBank.south``'
+        )
+
     if isinstance(f, str | Path):
         with h5py.File(f) as ff:
             return load_beer_mcstas(ff, bank=bank)
 
-    return _load_beer_mcstas(f, bank=bank)
+    try:
+        _find_h5(f['/entry1/instrument/components'], '.*nD_Mantid_?south_1.*')
+    except ValueError:
+        # The file did not have a detector named 'south'-something.
+        # Load old 2D structure where banks were not named 'north' and 'south'.
+        return _load_beer_mcstas(
+            f, north_or_south=None, number=1 if bank == DetectorBank.south else 2
+        )
+
+    return sc.concat(
+        [
+            _load_beer_mcstas(f, north_or_south=bank.name, number=number)
+            for number in range(1, 13)
+        ],
+        dim='panel',
+    )
+
+
+def load_beer_mcstas_monitor(f: str | Path | h5py.File):
+    if isinstance(f, str | Path):
+        with h5py.File(f) as ff:
+            return load_beer_mcstas_monitor(ff)
+    (
+        monitor,
+        wavelengths,
+        data,
+        errors,
+        ncount,
+    ) = _load_h5(
+        f,
+        'NXentry/NXdetector/Lmon_hereon_dat',
+        'NXentry/NXdetector/Lmon_hereon_dat/Wavelength__AA_',
+        'NXentry/NXdetector/Lmon_hereon_dat/data',
+        'NXentry/NXdetector/Lmon_hereon_dat/errors',
+        'NXentry/NXdetector/Lmon_hereon_dat/ncount',
+    )
+    da = sc.DataArray(
+        sc.array(
+            dims=['wavelength'], values=data[:], variances=errors[:], unit='counts'
+        ),
+        coords={
+            'wavelength': sc.array(
+                dims=['wavelength'], values=wavelengths[:], unit='angstrom'
+            ),
+            'ncount': sc.array(dims=['wavelength'], values=ncount[:], unit='counts'),
+        },
+    )
+    for name, value in monitor.attrs.items():
+        if name in ('position',):
+            da.coords[name] = sc.scalar(value.decode())
+
+    da.coords['position'] = sc.vector(
+        list(map(float, da.coords.pop('position').value.split(' '))), unit='m'
+    )
+    return da
 
 
 def load_beer_mcstas_provider(

src/ess/beer/types.py

@@ -7,6 +7,7 @@
 pipeline.
 """
 
+from enum import Enum
 from typing import NewType
 
 import sciline
@@ -25,7 +26,11 @@ class StreakClusteredData(sciline.Scope[RunType, sc.DataArray], sc.DataArray):
 SampleRun = SampleRun
 TofDetector = TofDetector
 
-DetectorBank = NewType('DetectorBank', int)
+
+class DetectorBank(Enum):
+    north = 'north'
+    south = 'south'
+
 
 TwoThetaLimits = NewType("TwoThetaLimits", tuple[sc.Variable, sc.Variable])
 

tests/beer/mcstas_reduction_test.py

@@ -8,7 +8,13 @@
     BeerModMcStasWorkflow,
     BeerModMcStasWorkflowKnownPeaks,
 )
-from ess.beer.data import duplex_peaks_array, mcstas_duplex, mcstas_silicon_new_model
+from ess.beer.data import (
+    duplex_peaks_array,
+    mcstas_duplex,
+    mcstas_few_neutrons_3d_detector_example,
+    mcstas_silicon_new_model,
+)
+from ess.beer.io import load_beer_mcstas, load_beer_mcstas_monitor
 from ess.beer.types import DetectorBank, DHKLList
 from ess.reduce.nexus.types import Filename, SampleRun
 from ess.reduce.time_of_flight.types import TofDetector
@@ -17,7 +23,7 @@
 def test_can_reduce_using_known_peaks_workflow():
     wf = BeerModMcStasWorkflowKnownPeaks()
     wf[DHKLList] = duplex_peaks_array()
-    wf[DetectorBank] = 1
+    wf[DetectorBank] = DetectorBank.north
     wf[Filename[SampleRun]] = mcstas_duplex(7)
     da = wf.compute(TofDetector[SampleRun])
     assert 'tof' in da.bins.coords
@@ -38,7 +44,7 @@ def test_can_reduce_using_known_peaks_workflow():
 def test_can_reduce_using_unknown_peaks_workflow():
     wf = BeerModMcStasWorkflow()
     wf[Filename[SampleRun]] = mcstas_duplex(7)
-    wf[DetectorBank] = 1
+    wf[DetectorBank] = DetectorBank.north
     da = wf.compute(TofDetector[SampleRun])
     da = da.transform_coords(
         ('dspacing',),
@@ -56,7 +62,7 @@ def test_can_reduce_using_unknown_peaks_workflow():
 def test_pulse_shaping_workflow():
     wf = BeerMcStasWorkflowPulseShaping()
     wf[Filename[SampleRun]] = mcstas_silicon_new_model(6)
-    wf[DetectorBank] = 1
+    wf[DetectorBank] = DetectorBank.north
     da = wf.compute(TofDetector[SampleRun])
     assert 'tof' in da.bins.coords
     # assert dataarray has all coords required to compute dspacing
@@ -71,3 +77,17 @@ def test_pulse_shaping_workflow():
         sc.scalar(1.6374, unit='angstrom'),
         atol=sc.scalar(1e-2, unit='angstrom'),
     )
+
+
+def test_can_load_3d_detector():
+    load_beer_mcstas(mcstas_few_neutrons_3d_detector_example(), DetectorBank.north)
+    da = load_beer_mcstas(mcstas_few_neutrons_3d_detector_example(), DetectorBank.south)
+    assert 'panel' in da.dims
+
+
+def test_can_load_monitor():
+    da = load_beer_mcstas_monitor(mcstas_few_neutrons_3d_detector_example())
+    assert 'wavelength' in da.coords
+    assert 'position' in da.coords
+    assert da.coords['position'].dtype == sc.DType.vector3
+    assert da.coords['position'].unit == 'm'