diff --git a/src/ess/dream/io/geant4.py b/src/ess/dream/io/geant4.py index d5060e4a..d9aec0ed 100644 --- a/src/ess/dream/io/geant4.py +++ b/src/ess/dream/io/geant4.py @@ -10,10 +10,10 @@ CalibrationData, CalibrationFilename, Filename, + NeXusDetector, NeXusDetectorDimensions, NeXusDetectorName, RawDetector, - RawDetectorData, RawSample, RawSource, ReducibleDetectorData, @@ -64,16 +64,16 @@ def load_geant4_csv(file_path: Filename[RunType]) -> AllRawDetectors[RunType]: def extract_geant4_detector( detectors: AllRawDetectors[RunType], detector_name: NeXusDetectorName -) -> RawDetector[RunType]: +) -> NeXusDetector[RunType]: """Extract a single detector from a loaded GEANT4 simulation.""" - return RawDetector[RunType](detectors["instrument"][detector_name]) + return NeXusDetector[RunType](detectors["instrument"][detector_name]) def extract_geant4_detector_data( - detector: RawDetector[RunType], -) -> RawDetectorData[RunType]: + detector: NeXusDetector[RunType], +) -> RawDetector[RunType]: """Extract the histogram or event data from a loaded GEANT4 detector.""" - return RawDetectorData[RunType](extract_detector_data(detector)) + return RawDetector[RunType](extract_detector_data(detector)) def _load_raw_events(file_path: str) -> sc.DataArray: @@ -176,7 +176,7 @@ def get_sample_position(raw_sample: RawSample[RunType]) -> SamplePosition[RunTyp def patch_detector_data( - detector_data: RawDetectorData[RunType], + detector_data: RawDetector[RunType], source_position: SourcePosition[RunType], sample_position: SamplePosition[RunType], ) -> ReducibleDetectorData[RunType]: @@ -188,7 +188,7 @@ def patch_detector_data( def geant4_detector_dimensions( - data: RawDetectorData[SampleRun], + data: RawDetector[SampleRun], ) -> NeXusDetectorDimensions: # For geant4 data, we group by detector identifier, so the data already has # logical dimensions, so we simply return the dimensions of the detector. diff --git a/src/ess/dream/io/nexus.py b/src/ess/dream/io/nexus.py index 98fb618f..3e9755a3 100644 --- a/src/ess/dream/io/nexus.py +++ b/src/ess/dream/io/nexus.py @@ -13,21 +13,7 @@ but it is not possible to reshape the data into all the logical dimensions. """ -import scipp as sc -from ess.reduce import nexus - -from ess.powder.types import ( - Filename, - LoadedNeXusDetector, - NeXusDetectorName, - RawDetectorData, - RawSample, - RawSource, - ReducibleDetectorData, - RunType, - SamplePosition, - SourcePosition, -) +from ess import powder DETECTOR_BANK_SIZES = { "endcap_backward_detector": { @@ -51,95 +37,19 @@ "strip": 256, "counter": 2, }, - "high_resolution_detector": { - "strip": 32, - "other": -1, - }, + "high_resolution_detector": {"strip": 32, "other": -1}, "sans_detector": lambda x: x.fold( dim="detector_number", - sizes={ - "strip": 32, - "other": -1, - }, + sizes={"strip": 32, "other": -1}, ), } -def load_nexus_sample(file_path: Filename[RunType]) -> RawSample[RunType]: - return RawSample[RunType](nexus.load_sample(file_path)) - - -def dummy_load_sample(file_path: Filename[RunType]) -> RawSample[RunType]: - """ - In test files there is not always a sample, so we need a dummy. - """ - return RawSample[RunType]( - sc.DataGroup({'position': sc.vector(value=[0, 0, 0], unit='m')}) - ) - - -def load_nexus_source(file_path: Filename[RunType]) -> RawSource[RunType]: - return RawSource[RunType](nexus.load_source(file_path)) - - -def load_nexus_detector( - file_path: Filename[RunType], detector_name: NeXusDetectorName -) -> LoadedNeXusDetector[RunType]: - out = nexus.load_detector(file_path=file_path, detector_name=detector_name) - out.pop("pixel_shape", None) - return LoadedNeXusDetector[RunType](out) - - -def get_source_position( - raw_source: RawSource[RunType], -) -> SourcePosition[RunType]: - return SourcePosition[RunType](raw_source["position"]) - - -def get_sample_position( - raw_sample: RawSample[RunType], -) -> SamplePosition[RunType]: - return SamplePosition[RunType](raw_sample["position"]) - - -def get_detector_data( - detector: LoadedNeXusDetector[RunType], - detector_name: NeXusDetectorName, -) -> RawDetectorData[RunType]: - da = nexus.extract_detector_data(detector) - if detector_name in DETECTOR_BANK_SIZES: - da = da.fold(dim="detector_number", sizes=DETECTOR_BANK_SIZES[detector_name]) - return RawDetectorData[RunType](da) - - -def patch_detector_data( - detector_data: RawDetectorData[RunType], - source_position: SourcePosition[RunType], - sample_position: SamplePosition[RunType], -) -> ReducibleDetectorData[RunType]: - """ - Patch a detector data object with source and sample positions. - Also adds variances to the event data if they are missing. - """ - out = detector_data.copy(deep=False) - if out.bins is not None: - content = out.bins.constituents["data"] - if content.variances is None: - content.variances = content.values - out.coords["sample_position"] = sample_position - out.coords["source_position"] = source_position - return ReducibleDetectorData[RunType](out) +def dream_detector_bank_sizes() -> powder.types.DetectorBankSizes | None: + return powder.types.DetectorBankSizes(DETECTOR_BANK_SIZES) -providers = ( - load_nexus_sample, - load_nexus_source, - load_nexus_detector, - get_source_position, - get_sample_position, - get_detector_data, - patch_detector_data, -) +providers = (*powder.nexus.providers, dream_detector_bank_sizes) """ -Providers for loading and processing DREAM NeXus data. +Providers for loading and processing NeXus data. """ diff --git a/src/ess/powder/__init__.py b/src/ess/powder/__init__.py index 50ecb914..38ba4624 100644 --- a/src/ess/powder/__init__.py +++ b/src/ess/powder/__init__.py @@ -15,6 +15,7 @@ smoothing, ) from .masking import with_pixel_mask_filenames +from . import nexus try: __version__ = importlib.metadata.version(__package__ or __name__) @@ -39,6 +40,7 @@ "filtering", "grouping", "masking", + "nexus", "providers", "smoothing", "with_pixel_mask_filenames", diff --git a/src/ess/powder/nexus.py b/src/ess/powder/nexus.py new file mode 100644 index 00000000..46728900 --- /dev/null +++ b/src/ess/powder/nexus.py @@ -0,0 +1,289 @@ +# SPDX-License-Identifier: BSD-3-Clause +# Copyright (c) 2024 Scipp contributors (https://github.com/scipp) + +"""NeXus input/output for ESS powder reduction.""" + +from typing import Any + +import scipp as sc +import scippnexus as snx +from ess.reduce import nexus + +from ess.powder.types import ( + DetectorBankSizes, + DetectorEventData, + Filename, + MonitorEventData, + MonitorType, + NeXusDetector, + NeXusDetectorName, + NeXusMonitor, + NeXusMonitorName, + RawDetector, + RawMonitor, + RawMonitorData, + RawSample, + RawSource, + ReducibleDetectorData, + RunType, + SamplePosition, + SourcePosition, +) + + +def load_nexus_sample(file_path: Filename[RunType]) -> RawSample[RunType]: + return RawSample[RunType](nexus.load_sample(file_path)) + + +def dummy_load_sample(file_path: Filename[RunType]) -> RawSample[RunType]: + """ + In test files there is not always a sample, so we need a dummy. + """ + return RawSample[RunType]( + sc.DataGroup({'position': sc.vector(value=[0, 0, 0], unit='m')}) + ) + + +def load_nexus_source(file_path: Filename[RunType]) -> RawSource[RunType]: + return RawSource[RunType](nexus.load_source(file_path)) + + +def load_nexus_detector( + file_path: Filename[RunType], detector_name: NeXusDetectorName +) -> NeXusDetector[RunType]: + """ + Load detector from NeXus, but with event data replaced by placeholders. + + Currently the placeholder is the detector number, but this may change in the future. + + The returned object is a scipp.DataGroup, as it may contain additional information + about the detector that cannot be represented as a single scipp.DataArray. Most + downstream code will only be interested in the contained scipp.DataArray so this + needs to be extracted. However, other processing steps may require the additional + information, so it is kept in the DataGroup. + + Loading thus proceeds in three steps: + + 1. This function loads the detector, but replaces the event data with placeholders. + 2. :py:func:`get_detector_array` drops the additional information, returning only + the contained scipp.DataArray, reshaped to the logical detector shape. + This will generally contain coordinates as well as pixel masks. + 3. :py:func:`assemble_detector_data` replaces placeholder data values with the + event data, and adds source and sample positions. + """ + definitions = snx.base_definitions() + definitions["NXdetector"] = _StrippedDetector + dg = nexus.load_detector( + file_path=file_path, + detector_name=detector_name, + definitions=definitions, + ) + # The name is required later, e.g., for determining logical detector shape + dg['detector_name'] = detector_name + return NeXusDetector[RunType](dg) + + +def load_nexus_monitor( + file_path: Filename[RunType], monitor_name: NeXusMonitorName[MonitorType] +) -> NeXusMonitor[RunType, MonitorType]: + """ + Load monitor from NeXus, but with event data replaced by placeholders. + + Currently the placeholder is a size-0 array, but this may change in the future. + + The returned object is a scipp.DataGroup, as it may contain additional information + about the monitor that cannot be represented as a single scipp.DataArray. Most + downstream code will only be interested in the contained scipp.DataArray so this + needs to be extracted. However, other processing steps may require the additional + information, so it is kept in the DataGroup. + + Loading thus proceeds in three steps: + + 1. This function loads the monitor, but replaces the event data with placeholders. + 2. :py:func:`get_monitor_array` drops the additional information, returning only + the contained scipp.DataArray. + This will generally contain coordinates as well as pixel masks. + 3. :py:func:`assemble_monitor_data` replaces placeholder data values with the + event data, and adds source and sample positions. + """ + definitions = snx.base_definitions() + definitions["NXmonitor"] = _StrippedMonitor + monitor = nexus.load_monitor( + file_path=file_path, monitor_name=monitor_name, definitions=definitions + ) + return NeXusMonitor[RunType, MonitorType](monitor) + + +def get_source_position( + raw_source: RawSource[RunType], +) -> SourcePosition[RunType]: + return SourcePosition[RunType](raw_source["position"]) + + +def get_sample_position( + raw_sample: RawSample[RunType], +) -> SamplePosition[RunType]: + return SamplePosition[RunType](raw_sample["position"]) + + +def get_detector_signal_array( + detector: NeXusDetector[RunType], + bank_sizes: DetectorBankSizes | None = None, +) -> RawDetector[RunType]: + """ + Extract the data array corresponding to a detector's signal field. + + The returned data array includes coords and masks pertaining directly to the + signal values array, but not additional information about the detector. The + data array is reshaped to the logical detector shape, which by folding the data + array along the detector_number dimension. + """ + da = nexus.extract_detector_data(detector) + if (sizes := (bank_sizes or {}).get(detector['detector_name'])) is not None: + da = da.fold(dim="detector_number", sizes=sizes) + return RawDetector[RunType](da) + + +def get_monitor_signal_array( + monitor: NeXusMonitor[RunType, MonitorType], + source_position: SourcePosition[RunType], +) -> RawMonitor[RunType, MonitorType]: + """ + Extract the data array corresponding to a monitor's signal field. + + The returned data array includes coords pertaining directly to the + signal values array, but not additional information about the monitor. + """ + return RawMonitor[RunType, MonitorType]( + nexus.extract_monitor_data(monitor).assign_coords( + position=monitor['position'], source_position=source_position + ) + ) + + +def assemble_detector_data( + detector: RawDetector[RunType], + event_data: DetectorEventData[RunType], + source_position: SourcePosition[RunType], + sample_position: SamplePosition[RunType], +) -> ReducibleDetectorData[RunType]: + """ + Assemble a detector data array with event data and source- and sample-position. + + Also adds variances to the event data if they are missing. + """ + grouped = nexus.group_event_data( + event_data=event_data, detector_number=detector.coords['detector_number'] + ) + return ReducibleDetectorData[RunType]( + _add_variances(grouped) + .assign_coords(source_position=source_position, sample_position=sample_position) + .assign_coords(detector.coords) + .assign_masks(detector.masks) + ) + + +def assemble_monitor_data( + monitor: RawMonitor[RunType, MonitorType], + event_data: MonitorEventData[RunType, MonitorType], +) -> RawMonitorData[RunType, MonitorType]: + """ + Assemble a monitor data array with event data. + + Also adds variances to the event data if they are missing. + """ + da = event_data.assign_coords(monitor.coords).assign_masks(monitor.masks) + return RawMonitorData[RunType, MonitorType](_add_variances(da=da)) + + +def _drop( + children: dict[str, snx.Field | snx.Group], classes: tuple[snx.NXobject, ...] +) -> dict[str, snx.Field | snx.Group]: + return { + name: child + for name, child in children.items() + if not (isinstance(child, snx.Group) and (child.nx_class in classes)) + } + + +class _StrippedDetector(snx.NXdetector): + """Detector definition without large geometry or event data for ScippNexus. + + Drops NXoff_geometry and NXevent_data groups, data is replaced by detector_number. + """ + + def __init__( + self, attrs: dict[str, Any], children: dict[str, snx.Field | snx.Group] + ): + children = _drop(children, (snx.NXoff_geometry, snx.NXevent_data)) + children['data'] = children['detector_number'] + super().__init__(attrs=attrs, children=children) + + +class _DummyField: + """Dummy field that can replace snx.Field in NXmonitor.""" + + def __init__(self): + self.attrs = {} + self.sizes = {'event_time_zero': 0} + self.dims = ('event_time_zero',) + self.shape = (0,) + + def __getitem__(self, key: Any) -> sc.Variable: + return sc.empty(dims=self.dims, shape=self.shape, unit=None) + + +class _StrippedMonitor(snx.NXmonitor): + """Monitor definition without event data for ScippNexus. + + Drops NXevent_data group, data is replaced by a dummy field. + """ + + def __init__( + self, attrs: dict[str, Any], children: dict[str, snx.Field | snx.Group] + ): + children = _drop(children, (snx.NXevent_data,)) + children['data'] = _DummyField() + super().__init__(attrs=attrs, children=children) + + +def load_detector_event_data( + file_path: Filename[RunType], detector_name: NeXusDetectorName +) -> DetectorEventData[RunType]: + da = nexus.load_event_data(file_path=file_path, component_name=detector_name) + return DetectorEventData[RunType](da) + + +def load_monitor_event_data( + file_path: Filename[RunType], monitor_name: NeXusMonitorName[MonitorType] +) -> MonitorEventData[RunType, MonitorType]: + da = nexus.load_event_data(file_path=file_path, component_name=monitor_name) + return MonitorEventData[RunType, MonitorType](da) + + +def _add_variances(da: sc.DataArray) -> sc.DataArray: + out = da.copy(deep=False) + if out.bins is not None: + content = out.bins.constituents['data'] + if content.variances is None: + content.variances = content.values + return out + + +providers = ( + assemble_detector_data, + assemble_monitor_data, + get_detector_signal_array, + get_monitor_signal_array, + get_sample_position, + get_source_position, + load_detector_event_data, + load_monitor_event_data, + load_nexus_detector, + load_nexus_monitor, + load_nexus_sample, + load_nexus_source, +) +""" +Providers for loading and processing NeXus data. +""" diff --git a/src/ess/powder/types.py b/src/ess/powder/types.py index 658bcc1a..5697475d 100644 --- a/src/ess/powder/types.py +++ b/src/ess/powder/types.py @@ -28,9 +28,18 @@ RunType = TypeVar("RunType", EmptyInstrumentRun, SampleRun, VanadiumRun) """TypeVar used for specifying the run.""" +# 1.2 Monitor types +Monitor1 = NewType('Monitor1', int) +"""Placeholder for monitor 1.""" +Monitor2 = NewType('Monitor2', int) +"""Placeholder for monitor 2.""" +MonitorType = TypeVar('MonitorType', Monitor1, Monitor2) +"""TypeVar used for identifying a monitor""" # 2 Workflow parameters +DetectorBankSizes = NewType("DetectorBankSizes", dict[str, dict[str, int | Any]]) + CalibrationFilename = NewType("CalibrationFilename", str | None) """Filename of the instrument calibration file.""" @@ -38,6 +47,11 @@ NeXusDetectorName = NewType("NeXusDetectorName", str) """Name of detector entry in NeXus file""" + +class NeXusMonitorName(sciline.Scope[MonitorType, str], str): + """Name of Incident|Transmission monitor in NeXus file""" + + DspacingBins = NewType("DSpacingBins", sc.Variable) """Bin edges for d-spacing.""" @@ -126,9 +140,46 @@ class FocussedDataDspacingTwoTheta(sciline.Scope[RunType, sc.DataArray], sc.Data """Data that has been normalized by a vanadium run, and grouped into 2theta bins.""" -class LoadedNeXusDetector(sciline.Scope[RunType, sc.DataGroup], sc.DataGroup): - """Detector data, loaded from a NeXus file, containing not only neutron events - but also pixel shape information, transformations, ...""" +class NeXusDetector(sciline.Scope[RunType, sc.DataGroup], sc.DataGroup): + """ + Detector loaded from a NeXus file, without event data. + + Contains detector numbers, pixel shape information, transformations, ... + """ + + +class NeXusMonitor( + sciline.ScopeTwoParams[RunType, MonitorType, sc.DataGroup], sc.DataGroup +): + """ + Monitor loaded from a NeXus file, without event data. + + Contains detector numbers, pixel shape information, transformations, ... + """ + + +class DetectorEventData(sciline.Scope[RunType, sc.DataArray], sc.DataArray): + """Event data loaded from a detector in a NeXus file""" + + +class MonitorEventData( + sciline.ScopeTwoParams[RunType, MonitorType, sc.DataArray], sc.DataArray +): + """Event data loaded from a monitor in a NeXus file""" + + +class RawMonitor( + sciline.ScopeTwoParams[RunType, MonitorType, sc.DataArray], sc.DataArray +): + """Raw monitor data""" + + +class RawMonitorData( + sciline.ScopeTwoParams[RunType, MonitorType, sc.DataArray], sc.DataArray +): + """Raw monitor data where variances and necessary coordinates + (e.g. source position) have been added, and where optionally some + user configuration was applied to some of the coordinates.""" class MaskedData(sciline.Scope[RunType, sc.DataArray], sc.DataArray): @@ -156,11 +207,7 @@ class RawDataAndMetadata(sciline.Scope[RunType, sc.DataGroup], sc.DataGroup): """Raw data and associated metadata.""" -class RawDetector(sciline.Scope[RunType, sc.DataGroup], sc.DataGroup): - """Full raw data for a detector.""" - - -class RawDetectorData(sciline.Scope[RunType, sc.DataArray], sc.DataArray): +class RawDetector(sciline.Scope[RunType, sc.DataArray], sc.DataArray): """Data (events / histogram) extracted from a RawDetector.""" diff --git a/tests/dream/io/geant4_test.py b/tests/dream/io/geant4_test.py index 74522663..70adecb0 100644 --- a/tests/dream/io/geant4_test.py +++ b/tests/dream/io/geant4_test.py @@ -11,7 +11,7 @@ import scipp.testing from ess.dream import data, load_geant4_csv -from ess.powder.types import Filename, NeXusDetectorName, RawDetectorData, SampleRun +from ess.powder.types import Filename, NeXusDetectorName, RawDetector, SampleRun @pytest.fixture(scope="module") @@ -180,6 +180,6 @@ def test_geant4_in_pipeline(file_path, file): NeXusDetectorName: NeXusDetectorName("mantle"), }, ) - detector = pipeline.compute(RawDetectorData[SampleRun]) + detector = pipeline.compute(RawDetector[SampleRun]) expected = load_geant4_csv(file)["instrument"]["mantle"]["events"] sc.testing.assert_identical(detector, expected) diff --git a/tests/dream/io/nexus_test.py b/tests/dream/io/nexus_test.py index 2e684279..7b28e4f1 100644 --- a/tests/dream/io/nexus_test.py +++ b/tests/dream/io/nexus_test.py @@ -2,14 +2,18 @@ # Copyright (c) 2023 Scipp contributors (https://github.com/scipp) import pytest import sciline -from ess import dream +import scipp as sc +from ess import dream, powder import ess.dream.data # noqa: F401 from ess.dream import nexus from ess.powder.types import ( Filename, + Monitor1, NeXusDetectorName, - RawDetectorData, + NeXusMonitorName, + RawDetector, + RawMonitor, ReducibleDetectorData, SampleRun, ) @@ -20,15 +24,7 @@ @pytest.fixture() def providers(): - return ( - nexus.dummy_load_sample, - nexus.load_nexus_source, - nexus.load_nexus_detector, - nexus.get_sample_position, - nexus.get_source_position, - nexus.get_detector_data, - nexus.patch_detector_data, - ) + return (*nexus.providers, powder.nexus.dummy_load_sample) @pytest.fixture( @@ -50,7 +46,7 @@ def params(request): def test_can_load_nexus_detector_data(providers, params): pipeline = sciline.Pipeline(params=params, providers=providers) - result = pipeline.compute(RawDetectorData[SampleRun]) + result = pipeline.compute(RawDetector[SampleRun]) assert ( set(result.dims) == hr_sans_dims if params[NeXusDetectorName] @@ -61,6 +57,18 @@ def test_can_load_nexus_detector_data(providers, params): else bank_dims ) + assert sc.identical(result.data, result.coords['detector_number']) + + +def test_can_load_nexus_monitor_data(providers): + pipeline = sciline.Pipeline(providers=providers) + pipeline[Filename[SampleRun]] = dream.data.get_path( + 'DREAM_nexus_sorted-2023-12-07.nxs' + ) + pipeline[NeXusMonitorName[Monitor1]] = 'monitor_cave' + result = pipeline.compute(RawMonitor[SampleRun, Monitor1]) + assert result.sizes == {'event_time_zero': 0} + def test_load_fails_with_bad_detector_name(providers): params = { @@ -69,10 +77,10 @@ def test_load_fails_with_bad_detector_name(providers): } pipeline = sciline.Pipeline(params=params, providers=providers) with pytest.raises(KeyError, match='bad_detector'): - pipeline.compute(RawDetectorData[SampleRun]) + pipeline.compute(RawDetector[SampleRun]) -def test_patch_nexus_detector_data(providers, params): +def test_assemble_nexus_detector_data(providers, params): pipeline = sciline.Pipeline(params=params, providers=providers) result = pipeline.compute(ReducibleDetectorData[SampleRun]) assert ( diff --git a/tests/powder/filtering_test.py b/tests/powder/filtering_test.py index 95661cae..a32968a0 100644 --- a/tests/powder/filtering_test.py +++ b/tests/powder/filtering_test.py @@ -39,12 +39,8 @@ def make_data_with_pulse_time(rng, n_event) -> sc.DataArray: ), }, ) - return sc.binning.make_binned( - events, - edges=[ - sc.array(dims=['tof'], values=[10, 500, 1000], unit='us', dtype='int64') - ], - groups=[sc.arange('spectrum', 0, 10, unit=None, dtype='int64')], + return events.group(sc.arange('spectrum', 0, 10, unit=None, dtype='int64')).bin( + tof=sc.array(dims=['tof'], values=[10, 500, 1000], unit='us', dtype='int64') )