diff --git a/packages/essdiffraction/docs/user-guide/beer/beer_modulation_mcstas.ipynb b/packages/essdiffraction/docs/user-guide/beer/beer_modulation_mcstas.ipynb index e59c8e52e..7d07a471b 100644 --- a/packages/essdiffraction/docs/user-guide/beer/beer_modulation_mcstas.ipynb +++ b/packages/essdiffraction/docs/user-guide/beer/beer_modulation_mcstas.ipynb @@ -713,7 +713,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.11.13" + "version": "3.12.13" } }, "nbformat": 4, diff --git a/packages/essdiffraction/docs/user-guide/beer/beer_powder_mcstas_analytical.ipynb b/packages/essdiffraction/docs/user-guide/beer/beer_powder_mcstas_analytical.ipynb index 1858282d3..c9d64ea9f 100644 --- a/packages/essdiffraction/docs/user-guide/beer/beer_powder_mcstas_analytical.ipynb +++ b/packages/essdiffraction/docs/user-guide/beer/beer_powder_mcstas_analytical.ipynb @@ -80,7 +80,12 @@ "outputs": [], "source": [ "workflow[DetectorBank] = DetectorBank.north\n", - "workflow[DspacingBins] = sc.linspace(\"dspacing\", 0.6, 2.0, 701, unit=\"angstrom\")" + "\n", + "# Select the number of dspacing bins.\n", + "# By default the range of `DspacingBins` (the binning in dspacing) is\n", + "# selected to be the dspacing range covered by the instrument.\n", + "workflow[DspacingNBins] = 1000\n", + "workflow.compute(DspacingBins)" ] }, { @@ -198,8 +203,6 @@ " title=\"Intensity over dspacing\",\n", " xlabel=r\"$d$-spacing [$\\AA$]\",\n", " ylabel=f\"Intensity [{histogram.unit}]\",\n", - " xmin=sc.scalar(0.6, unit=\"angstrom\"),\n", - " xmax=sc.scalar(2.0, unit=\"angstrom\"),\n", ")\n", "\n", "ymin, ymax = fig.ax.get_ylim()\n", @@ -308,8 +311,6 @@ " },\n", " xlabel=r\"$d$-spacing [$\\AA$]\",\n", " ylabel=f\"Intensity [{histogram.unit}]\",\n", - " xmin=sc.scalar(0.6, unit=\"angstrom\"),\n", - " xmax=sc.scalar(2.0, unit=\"angstrom\"),\n", ")\n", "\n", "ymin, ymax = fig.ax.get_ylim()\n", @@ -337,7 +338,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.11.13" + "version": "3.12.13" } }, "nbformat": 4, diff --git a/packages/essdiffraction/src/ess/beer/mcstas/__init__.py b/packages/essdiffraction/src/ess/beer/mcstas/__init__.py index 5c2cf0a14..0ae4dbb8c 100644 --- a/packages/essdiffraction/src/ess/beer/mcstas/__init__.py +++ b/packages/essdiffraction/src/ess/beer/mcstas/__init__.py @@ -11,6 +11,8 @@ from .load import ( chopper_mode_from_mcstas_mode, load_beer_mcstas, + load_beer_mcstas_geometry, + load_beer_mcstas_geometry_provider, load_beer_mcstas_monitor, load_beer_mcstas_monitor_provider, load_beer_mcstas_provider, @@ -26,6 +28,8 @@ 'PulseShapingMode', 'chopper_mode_from_mcstas_mode', 'load_beer_mcstas', + 'load_beer_mcstas_geometry', + 'load_beer_mcstas_geometry_provider', 'load_beer_mcstas_monitor', 'load_beer_mcstas_monitor_provider', 'load_beer_mcstas_provider', diff --git a/packages/essdiffraction/src/ess/beer/mcstas/load.py b/packages/essdiffraction/src/ess/beer/mcstas/load.py index 5ab8d918d..65427976c 100644 --- a/packages/essdiffraction/src/ess/beer/mcstas/load.py +++ b/packages/essdiffraction/src/ess/beer/mcstas/load.py @@ -7,7 +7,7 @@ import mcstastox import scipp as sc import scippnexus as snx -from ess.powder.types import CaveMonitor, RunType, WavelengthMonitor +from ess.powder.types import CaveMonitor, EmptyDetector, RunType, WavelengthMonitor from ess.reduce.nexus.types import DiskChoppers, Position from ess.reduce.unwrap.types import DetectorLtotal @@ -24,6 +24,8 @@ __all__ = [ 'chopper_mode_from_mcstas_mode', 'load_beer_mcstas', + 'load_beer_mcstas_geometry', + 'load_beer_mcstas_geometry_provider', 'load_beer_mcstas_monitor', 'load_beer_mcstas_monitor_provider', 'load_beer_mcstas_provider', @@ -62,6 +64,25 @@ def _detector_components(data: mcstastox.Read, bank: DetectorBank) -> list[str]: ) +def _load_position_table(data: mcstastox.Read, components: list[str]) -> sc.DataArray: + """Load pixel positions for selected detector components.""" + position_tables = [] + for component in components: + positions = data.get_component_global(component) + begin, end = data.pixel_range[component] + position_tables.append( + sc.DataArray( + sc.vectors(dims=['pixel_id'], values=positions, unit='m'), + coords={ + 'pixel_id': sc.arange( + 'pixel_id', int(begin), int(end) + 1, dtype='int32' + ) + }, + ) + ) + return sc.sort(sc.concat(position_tables, dim='pixel_id'), key='pixel_id') + + def _load_events(data: mcstastox.Read, components: list[str]) -> sc.DataArray: """Load weighted events and pixel positions for selected components only. @@ -86,30 +107,22 @@ def _load_events(data: mcstastox.Read, components: list[str]) -> sc.DataArray: ), 't': sc.array(dims=['events'], values=event_data['t'], unit='s'), }, - ).group('pixel_id') - - position_tables = [] - for component in components: - positions = data.get_component_global(component) - begin, end = data.pixel_range[component] - position_tables.append( - sc.DataArray( - sc.vectors(dims=['pixel_id'], values=positions, unit='m'), - coords={ - 'pixel_id': sc.arange( - 'pixel_id', int(begin), int(end) + 1, dtype='int32' - ) - }, - ) - ) + ) - position_table = sc.sort(sc.concat(position_tables, dim='pixel_id'), key='pixel_id') - events.coords['position'] = sc.lookup( - position_table, dim='pixel_id', mode='previous' - )[events.coords['pixel_id']] + position_table = _load_position_table(data, components) + # Group by the complete list of detector pixels, rather than by the IDs that + # occur in the event table. This retains pixels with no events as empty bins and + # keeps detector geometry independent of the measured event population. + events = events.group(position_table.coords['pixel_id']) + events.coords['position'] = position_table.data return events +def _load_mode(data: mcstastox.Read) -> str: + mode = data.file['entry1/simulation/Param/mode'][0] + return mode.decode() if isinstance(mode, bytes) else str(mode) + + def load_beer_mcstas( filename: str | Path, bank: DetectorBank, @@ -132,8 +145,7 @@ def load_beer_mcstas( filename = Path(filename) with mcstastox.Read(filename.parent, filename.name) as data: - mode = data.file['entry1/simulation/Param/mode'][0] - mode = mode.decode() if isinstance(mode, bytes) else str(mode) + mode = _load_mode(data) events = _load_events(data, _detector_components(data, bank)) events.coords['mode'] = sc.scalar(mode) @@ -143,6 +155,40 @@ def load_beer_mcstas( return events +def load_beer_mcstas_geometry( + filename: str | Path, + bank: DetectorBank, +) -> sc.DataArray: + """Load detector geometry from a BEER McStas file without loading events.""" + if not isinstance(bank, DetectorBank): + raise ValueError( + 'bank must be ``DetectorBank.north``, ' + '``DetectorBank.south``, or ``DetectorBank.both``' + ) + + if bank == DetectorBank.both: + return sc.concat( + [ + load_beer_mcstas_geometry(filename, bank) + for bank in (DetectorBank.south, DetectorBank.north) + ], + dim='pixel_id', + ) + + filename = Path(filename) + with mcstastox.Read(filename.parent, filename.name) as data: + positions = _load_position_table(data, _detector_components(data, bank)) + + pixel_id = positions.coords['pixel_id'] + return sc.DataArray( + pixel_id, + coords={ + 'pixel_id': pixel_id, + 'position': positions.data, + }, + ) + + def _to_edges(centers: sc.Variable) -> sc.Variable: interior_edges = sc.midpoints(centers) return sc.concat( @@ -201,6 +247,13 @@ def load_beer_mcstas_provider( return load_beer_mcstas(fname, bank) +def load_beer_mcstas_geometry_provider( + fname: Filename[RunType], bank: DetectorBank +) -> EmptyDetector[RunType]: + """Sciline provider for loading BEER McStas detector geometry.""" + return EmptyDetector[RunType](load_beer_mcstas_geometry(fname, bank)) + + def load_beer_mcstas_monitor_provider( fname: Filename[RunType], ) -> WavelengthMonitor[RunType, CaveMonitor]: @@ -287,22 +340,34 @@ def mcstas_choppers( return simulation_choppers(mode, source_position) +def _mcstas_choppers_from_file( + fname: Filename[RunType], + source_position: Position[snx.NXsource, RunType], +) -> DiskChoppers[RunType]: + """Return BEER choppers without loading detector events.""" + fname = Path(fname) + with mcstastox.Read(fname.parent, fname.name) as data: + mode = chopper_mode_from_mcstas_mode(_load_mode(data)) + return simulation_choppers(mode, source_position) + + def mcstas_detector_ltotal( - da: RawDetector[RunType], + detector: RawDetector[RunType], source_position: Position[snx.NXsource, RunType], sample_position: Position[snx.NXsample, RunType], ) -> DetectorLtotal[RunType]: """Return moderator-to-detector flight path lengths for BEER McStas data.""" source_to_sample = sc.norm(sample_position - source_position) - sample_to_detector = sc.norm(da.coords['position'] - sample_position) + sample_to_detector = sc.norm(detector.coords['position'] - sample_position) return source_to_sample + sample_to_detector mcstas_providers = ( load_beer_mcstas_provider, + load_beer_mcstas_geometry_provider, load_beer_mcstas_monitor_provider, mcstas_source_position, mcstas_sample_position, - mcstas_choppers, + _mcstas_choppers_from_file, mcstas_detector_ltotal, ) diff --git a/packages/essdiffraction/src/ess/powder/__init__.py b/packages/essdiffraction/src/ess/powder/__init__.py index faa956396..62f5d3250 100644 --- a/packages/essdiffraction/src/ess/powder/__init__.py +++ b/packages/essdiffraction/src/ess/powder/__init__.py @@ -7,6 +7,7 @@ import importlib.metadata from . import ( + binning, calibration, conversion, correction, @@ -27,6 +28,7 @@ del importlib providers = ( + *binning.providers, *calibration.providers, *conversion.providers, *correction.providers, @@ -41,6 +43,7 @@ __all__ = [ "RunNormalization", "__version__", + "binning", "calibration", "conversion", "correction", diff --git a/packages/essdiffraction/src/ess/powder/binning.py b/packages/essdiffraction/src/ess/powder/binning.py new file mode 100644 index 000000000..fc877e5e4 --- /dev/null +++ b/packages/essdiffraction/src/ess/powder/binning.py @@ -0,0 +1,61 @@ +# SPDX-License-Identifier: BSD-3-Clause +# Copyright (c) 2026 Scipp contributors (https://github.com/scipp) +"""Automatic bin-edge generation for powder diffraction.""" + +import scipp as sc +import scippneutron as scn + +from ess.reduce.unwrap import ChopperFrameSequence + +from .types import ( + DetectorTwoTheta, + DspacingBins, + DspacingNBins, + RunType, + SampleRun, + WavelengthRange, +) + + +def wavelength_range_from_chopper_frames( + frames: ChopperFrameSequence[RunType], +) -> WavelengthRange[RunType]: + """Return the wavelength envelope transmitted by the chopper cascade.""" + return WavelengthRange[RunType](frames[-1].bounds()['wavelength']) + + +def dspacing_bins_from_wavelength_and_two_theta( + wavelength_range: WavelengthRange[SampleRun], + two_theta: DetectorTwoTheta[SampleRun], + nbins: DspacingNBins, +) -> DspacingBins: + """Make d-spacing bin edges from wavelength and detector geometry. + + Wavelength and two-theta are independent geometry and beamline quantities. + Pairing their opposite extrema and transforming them to d-spacing gives an + envelope that covers every detector pixel. + """ + if nbins < 1: + raise ValueError(f'DspacingNBins must be positive, got {nbins}.') + dspacing_range = scn.conversion.tof.dspacing_from_wavelength( + wavelength=sc.concat( + [wavelength_range.nanmin(), wavelength_range.nanmax()], dim='bound' + ), + two_theta=sc.concat([two_theta.nanmax(), two_theta.nanmin()], dim='bound'), + ) + return DspacingBins( + sc.linspace( + dim='dspacing', + start=dspacing_range[0].value, + stop=dspacing_range[-1].value, + num=nbins + 1, + unit=dspacing_range.unit, + ) + ) + + +providers = ( + wavelength_range_from_chopper_frames, + dspacing_bins_from_wavelength_and_two_theta, +) +"""Sciline providers for automatic powder-diffraction binning.""" diff --git a/packages/essdiffraction/src/ess/powder/types.py b/packages/essdiffraction/src/ess/powder/types.py index 5cc18dea6..f3dd55b4c 100644 --- a/packages/essdiffraction/src/ess/powder/types.py +++ b/packages/essdiffraction/src/ess/powder/types.py @@ -62,6 +62,9 @@ DspacingBins = NewType("DspacingBins", sc.Variable) """Bin edges for d-spacing.""" +DspacingNBins = NewType("DspacingNBins", int) +"""Number of d-spacing bins.""" + OutFilename = NewType("OutFilename", str) """Filename of the output.""" @@ -115,6 +118,10 @@ class DetectorTwoTheta(sciline.Scope[RunType, sc.Variable], sc.Variable): """Scattering angle (two-theta) for each detector pixel.""" +class WavelengthRange(sciline.Scope[RunType, sc.Variable], sc.Variable): + """Minimum and maximum wavelength transmitted by the chopper cascade.""" + + class ElasticCoordTransformGraph(sciline.Scope[RunType, dict], dict): """Graph for transforming coordinates in elastic scattering.""" diff --git a/packages/essdiffraction/tests/beer/mcstas_reduction_test.py b/packages/essdiffraction/tests/beer/mcstas_reduction_test.py index 6006a6bab..15b74cfc0 100644 --- a/packages/essdiffraction/tests/beer/mcstas_reduction_test.py +++ b/packages/essdiffraction/tests/beer/mcstas_reduction_test.py @@ -1,6 +1,7 @@ import importlib import sys +import mcstastox import numpy as np import pytest import scipp as sc @@ -20,11 +21,14 @@ ) from ess.beer.mcstas import ( load_beer_mcstas, + load_beer_mcstas_geometry, load_beer_mcstas_monitor, ) from ess.beer.types import DetectorBank, DHKLList, WavelengthDetector from ess.powder.types import ( + DspacingBins, DspacingDetector, + DspacingNBins, ElasticCoordTransformGraph, SampleRun, ) @@ -135,6 +139,24 @@ def test_powder_mcstas_analytical_workflow_computes_dspacing(): ) +def test_powder_workflow_computes_dspacing_bins_without_loading_events(monkeypatch): + wf = BeerPowderMcStasWorkflow() + wf[Filename[SampleRun]] = mcstas_silicon_new_model(6) + wf[DetectorBank] = DetectorBank.north + wf[DspacingNBins] = 123 + + def fail_if_events_are_loaded(*args, **kwargs): + raise AssertionError('event data must not be used to determine bin edges') + + monkeypatch.setattr(mcstastox.Read, 'get_event_data', fail_if_events_are_loaded) + + bins = wf.compute(DspacingBins) + + assert bins.sizes == {'dspacing': 124} + assert sc.all(sc.isfinite(bins)).value + assert sc.all(bins[1:] > bins[:-1]).value + + @pytest.mark.parametrize( 'fname', [ @@ -165,6 +187,16 @@ def test_load_both_detector_banks(): ) +def test_loaded_detector_includes_pixels_without_events(): + filename = mcstas_silicon_new_model(10) + detector = load_beer_mcstas(filename, DetectorBank.north) + geometry = load_beer_mcstas_geometry(filename, DetectorBank.north) + + assert sc.identical(detector.coords['pixel_id'], geometry.coords['pixel_id']) + assert sc.identical(detector.coords['position'], geometry.coords['position']) + assert sc.any(detector.bins.size() == sc.index(0)).value + + def test_loaded_mcstas_event_variances_are_squared_weights(): da = load_beer_mcstas(mcstas_few_neutrons_3d_detector_example(), DetectorBank.north) weights = da.bins.constituents['data'] diff --git a/packages/essdiffraction/tests/powder/binning_test.py b/packages/essdiffraction/tests/powder/binning_test.py new file mode 100644 index 000000000..c6153cef6 --- /dev/null +++ b/packages/essdiffraction/tests/powder/binning_test.py @@ -0,0 +1,84 @@ +# SPDX-License-Identifier: BSD-3-Clause +# Copyright (c) 2026 Scipp contributors (https://github.com/scipp) + +import pytest +import scipp as sc +from ess.powder.binning import ( + dspacing_bins_from_wavelength_and_two_theta, + wavelength_range_from_chopper_frames, +) +from ess.powder.types import ( + DetectorTwoTheta, + DspacingNBins, + SampleRun, + WavelengthRange, +) +from scippneutron.tof import chopper_cascade + + +def test_wavelength_range_uses_frame_after_last_chopper(): + source = chopper_cascade.FrameSequence.from_source_pulse( + time_min=sc.scalar(0.0, unit='ms'), + time_max=sc.scalar(3.0, unit='ms'), + wavelength_min=sc.scalar(0.5, unit='angstrom'), + wavelength_max=sc.scalar(5.0, unit='angstrom'), + ) + last_subframe = chopper_cascade.Subframe( + time=sc.array(dims=['vertex'], values=[1.0, 2.0], unit='ms'), + wavelength=sc.array(dims=['vertex'], values=[1.2, 3.4], unit='angstrom'), + ) + frames = chopper_cascade.FrameSequence( + [ + *source.frames, + chopper_cascade.Frame( + distance=sc.scalar(10.0, unit='m'), subframes=[last_subframe] + ), + ] + ) + + wavelength_range = wavelength_range_from_chopper_frames(frames) + + assert sc.identical( + wavelength_range, + sc.array(dims=['bound'], values=[1.2, 3.4], unit='angstrom'), + ) + + +def test_dspacing_bins_span_envelope_of_wavelength_and_two_theta(): + wavelength_range = WavelengthRange[SampleRun]( + sc.array(dims=['bound'], values=[1.0, 4.0], unit='angstrom') + ) + two_theta = DetectorTwoTheta[SampleRun]( + sc.array(dims=['pixel'], values=[30.0, 60.0, 90.0], unit='deg') + ) + + bins = dspacing_bins_from_wavelength_and_two_theta( + wavelength_range, two_theta, DspacingNBins(4) + ) + + assert bins.sizes == {'dspacing': 5} + assert sc.allclose( + bins[[0, -1]], + sc.array( + dims=['dspacing'], + values=[ + 1.0 / (2**0.5), + 4.0 / (2 * sc.sin(15.0 * sc.Unit('deg')).value), + ], + unit='angstrom', + ), + ) + + +def test_dspacing_nbins_must_be_positive(): + wavelength_range = WavelengthRange[SampleRun]( + sc.array(dims=['bound'], values=[1.0, 4.0], unit='angstrom') + ) + two_theta = DetectorTwoTheta[SampleRun]( + sc.array(dims=['pixel'], values=[30.0, 90.0], unit='deg') + ) + + with pytest.raises(ValueError, match='DspacingNBins must be positive'): + dspacing_bins_from_wavelength_and_two_theta( + wavelength_range, two_theta, DspacingNBins(0) + )