diff --git a/docs/user-guide/beer/beer_modulation_mcstas.ipynb b/docs/user-guide/beer/beer_modulation_mcstas.ipynb index d269ced1..39227ea9 100644 --- a/docs/user-guide/beer/beer_modulation_mcstas.ipynb +++ b/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", diff --git a/src/ess/beer/data.py b/src/ess/beer/data.py index 1e144080..21cd6dc5 100644 --- a/src/ess/beer/data.py +++ b/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') diff --git a/src/ess/beer/io.py b/src/ess/beer/io.py index 464961ec..91282dbe 100644 --- a/src/ess/beer/io.py +++ b/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( diff --git a/src/ess/beer/types.py b/src/ess/beer/types.py index 144cfbb7..c4655bf5 100644 --- a/src/ess/beer/types.py +++ b/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]) diff --git a/tests/beer/mcstas_reduction_test.py b/tests/beer/mcstas_reduction_test.py index f19ee808..40caec52 100644 --- a/tests/beer/mcstas_reduction_test.py +++ b/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'