diff --git a/src/ess/dream/__init__.py b/src/ess/dream/__init__.py new file mode 100644 index 00000000..b621a546 --- /dev/null +++ b/src/ess/dream/__init__.py @@ -0,0 +1,23 @@ +# SPDX-License-Identifier: BSD-3-Clause +# Copyright (c) 2023 Scipp contributors (https://github.com/scipp) + +""" +Components for DREAM +""" +import importlib.metadata + +from . import data +from .io import fold_nexus_detectors, load_nexus + +try: + __version__ = importlib.metadata.version(__package__ or __name__) +except importlib.metadata.PackageNotFoundError: + __version__ = "0.0.0" + +del importlib + +__all__ = [ + "data", + "fold_nexus_detectors", + "load_nexus", +] diff --git a/src/ess/dream/data.py b/src/ess/dream/data.py new file mode 100644 index 00000000..8b09dd43 --- /dev/null +++ b/src/ess/dream/data.py @@ -0,0 +1,34 @@ +# SPDX-License-Identifier: BSD-3-Clause +# Copyright (c) 2023 Scipp contributors (https://github.com/scipp) +_version = '1' + +__all__ = ['get_path'] + + +def _make_pooch(): + import pooch + + return pooch.create( + path=pooch.os_cache('ess/powgen'), + env='ESS_DATA_DIR', + base_url='https://public.esss.dk/groups/scipp/ess/dream/{version}/', + version=_version, + registry={ + 'DREAM_nexus_sorted-2023-12-07.nxs': 'md5:22824e14f6eb950d24a720b2a0e2cb66', + }, + ) + + +_pooch = _make_pooch() + + +def get_path(name: str, unzip: bool = False) -> str: + """ + Return the path to a data file bundled with ess.dream. + + This function only works with example data and cannot handle + paths to custom files. + """ + import pooch + + return _pooch.fetch(name, processor=pooch.Unzip() if unzip else None) diff --git a/src/ess/dream/io.py b/src/ess/dream/io.py new file mode 100644 index 00000000..b0ad7e38 --- /dev/null +++ b/src/ess/dream/io.py @@ -0,0 +1,133 @@ +# SPDX-License-Identifier: BSD-3-Clause +# Copyright (c) 2023 Scipp contributors (https://github.com/scipp) + +import os +from typing import Any, Dict, Union + +import scipp as sc +import scippnexus as snx + + +def load_nexus( + path: Union[str, os.PathLike], + *, + load_pixel_shape: bool = False, + entry: str = 'entry', + fold_detectors: bool = True, +) -> sc.DataGroup: + """ + Load an unprocessed DREAM NeXus file. + + See https://confluence.esss.lu.se/pages/viewpage.action?pageId=462000005 + and the ICD DREAM interface specification for details. + + Notes (2023-12-06): + + - Mounting-unit, cassette, and counter roughly correspond to the azimuthal angle + in the mantle detector. However, counter is reversed in the current files. This + may also be the case in the other detectors. + - The endcap detectors have a irregular structure that cannot be fully folded. + This is not a problem but note again that the counter may be reversed. It is + currently not clear if this is a bug in the file. + - The high-resolution detector has a very odd numbering scheme. The SANS detector + is using the same, but is not populated in the current files. The scheme + attempts to follows some sort of physical ordering in space (x,y,z), but it + looks partially messed up. + + Parameters + ---------- + path: + Path to the NeXus file. + load_pixel_shape: + If True, load the pixel shape from the file's NXoff_geometry group. This is + often unused by would slow down file loading. Default is False. + entry: + Name of the entry to load. Default is "entry". + fold_detectors: + If True, fold the detector data to (partially) mimic the logical detector + structure. Default is True. + + Returns + ------- + : + A data group with the loaded file contents. + """ + definitions = snx.base_definitions() + if not load_pixel_shape: + definitions["NXdetector"] = FilteredDetector + + with snx.File(path, definitions=definitions) as f: + dg = f[entry][()] + dg = snx.compute_positions(dg) + return fold_nexus_detectors(dg) if fold_detectors else dg + + +def fold_nexus_detectors(dg: sc.DataGroup) -> sc.DataGroup: + """ + Fold the detector data in a DREAM NeXus file. + + The detector banks in the returned data group will have a multi-dimensional shape, + following the logical structure as far as possible. Note that the full structure + cannot be folded, as some dimensions are irregular. + """ + dg = dg.copy() + dg['instrument'] = dg['instrument'].copy() + instrument = dg['instrument'] + mantle = instrument['mantle_detector'] + mantle['mantle_event_data'] = mantle['mantle_event_data'].fold( + dim='detector_number', + sizes={ + 'wire': 32, + 'mounting_unit': 5, + 'cassette': 6, + 'counter': 2, + 'strip': 256, + }, + ) + for direction in ('backward', 'forward'): + endcap = instrument[f'endcap_{direction}_detector'] + endcap[f'endcap_{direction}_event_data'] = endcap[ + f'endcap_{direction}_event_data' + ].fold( + dim='detector_number', + sizes={ + 'strip': 16, + 'wire': 16, + 'sector': 5 if direction == 'forward' else 11, + 'sumo_cass_ctr': -1, # sumo*cassette*counter, irregular, cannot fold + }, + ) + high_resolution = instrument['high_resolution_detector'] + high_resolution['high_resolution_event_data'] = high_resolution[ + 'high_resolution_event_data' + ].fold( + dim='detector_number', + sizes={ + 'strip': 32, + 'other': -1, # some random order that is hard to follow + }, + ) + sans = instrument['sans_detector'] + sans['sans_event_data'] = sans['sans_event_data'].fold( + dim='detector_number', + sizes={ + 'strip': 32, + 'other': -1, # some random order that is hard to follow + }, + ) + return dg + + +def _skip(name: str, obj: Union[snx.Field, snx.Group]) -> bool: + skip_classes = (snx.NXoff_geometry,) + return isinstance(obj, snx.Group) and (obj.nx_class in skip_classes) + + +class FilteredDetector(snx.NXdetector): + def __init__( + self, attrs: Dict[str, Any], children: Dict[str, Union[snx.Field, snx.Group]] + ): + children = { + name: child for name, child in children.items() if not _skip(name, child) + } + super().__init__(attrs=attrs, children=children) diff --git a/tests/dream/io_test.py b/tests/dream/io_test.py new file mode 100644 index 00000000..71ba2695 --- /dev/null +++ b/tests/dream/io_test.py @@ -0,0 +1,77 @@ +# SPDX-License-Identifier: BSD-3-Clause +# Copyright (c) 2023 Scipp contributors (https://github.com/scipp) +import pytest + +from ess import dream + + +@pytest.fixture +def filename(): + return dream.data.get_path('DREAM_nexus_sorted-2023-12-07.nxs') + + +@pytest.mark.filterwarnings("ignore:Failed to load /entry/instrument/monitor_bunker") +@pytest.mark.filterwarnings("ignore:Failed to load /entry/instrument/monitor_cave") +@pytest.mark.filterwarnings("ignore:Failed to load /entry/instrument/polarizer/rate") +@pytest.mark.filterwarnings("ignore:Failed to load /entry/instrument/sans_detector") +def test_load_nexus_loads_file(filename): + dg = dream.load_nexus(filename) + assert 'instrument' in dg + instr = dg['instrument'] + for name in ( + 'mantle', + 'endcap_backward', + 'endcap_forward', + 'high_resolution', + 'sans', + ): + assert f'{name}_detector' in instr + det = instr[f'{name}_detector'] + assert 'pixel_shape' not in det + + +def test_load_nexus_fails_if_entry_not_found(filename): + with pytest.raises(KeyError): + dream.load_nexus(filename, entry='foo') + + +@pytest.mark.filterwarnings("ignore:Failed to load /entry/instrument/monitor_bunker") +@pytest.mark.filterwarnings("ignore:Failed to load /entry/instrument/monitor_cave") +@pytest.mark.filterwarnings("ignore:Failed to load /entry/instrument/polarizer/rate") +@pytest.mark.filterwarnings("ignore:Failed to load /entry/instrument/sans_detector") +def test_load_nexus_folds_detectors_by_default(filename): + dg = dream.load_nexus(filename) + instr = dg['instrument'] + # sans_detector is not populated in the current files + for name in ('mantle', 'endcap_backward', 'endcap_forward', 'high_resolution'): + det = instr[f'{name}_detector'] + # There may be other dims, but some are irregular and this may be subject to + # change + assert 'strip' in det.dims + + +@pytest.mark.filterwarnings("ignore:Failed to load /entry/instrument/monitor_bunker") +@pytest.mark.filterwarnings("ignore:Failed to load /entry/instrument/monitor_cave") +@pytest.mark.filterwarnings("ignore:Failed to load /entry/instrument/polarizer/rate") +@pytest.mark.filterwarnings("ignore:Failed to load /entry/instrument/sans_detector") +def test_load_nexus_with_disabled_fold(filename): + dg = dream.load_nexus(filename, fold_detectors=False) + instr = dg['instrument'] + for name in ('mantle', 'endcap_backward', 'endcap_forward', 'high_resolution'): + det = instr[f'{name}_detector'] + assert det.dims == ('detector_number',) + + +@pytest.mark.filterwarnings("ignore:Failed to load /entry/instrument/monitor_bunker") +@pytest.mark.filterwarnings("ignore:Failed to load /entry/instrument/monitor_cave") +@pytest.mark.filterwarnings("ignore:Failed to load /entry/instrument/polarizer/rate") +@pytest.mark.filterwarnings("ignore:Failed to load /entry/instrument/sans_detector") +def test_load_nexus_with_pixel_shape(filename): + dg = dream.load_nexus(filename, load_pixel_shape=True) + assert 'instrument' in dg + instr = dg['instrument'] + # sans_detector is not populated in the current files + for name in ('mantle', 'endcap_backward', 'endcap_forward', 'high_resolution'): + assert f'{name}_detector' in instr + det = instr[f'{name}_detector'] + assert 'pixel_shape' in det