diff --git a/docs/user-guide/beer/beer_modulation_mcstas.ipynb b/docs/user-guide/beer/beer_modulation_mcstas.ipynb index 83fe0155..d269ced1 100644 --- a/docs/user-guide/beer/beer_modulation_mcstas.ipynb +++ b/docs/user-guide/beer/beer_modulation_mcstas.ipynb @@ -15,6 +15,8 @@ "metadata": {}, "outputs": [], "source": [ + "%matplotlib ipympl\n", + "import plopp as pp\n", "import scippneutron as scn\n", "\n", "from ess.beer import BeerModMcStasWorkflow, BeerModMcStasWorkflowKnownPeaks\n", @@ -159,8 +161,11 @@ "outputs": [], "source": [ "wf = BeerModMcStasWorkflowKnownPeaks()\n", + "wf[DetectorBank] = 2\n", "wf[Filename[SampleRun]] = mcstas_silicon_medium_resolution()\n", - "wf.compute(RawDetector[SampleRun])['bank1'].hist(two_theta=1000, event_time_offset=1000).plot(norm='log', cmin=1.0e-3)" + "da = wf.compute(RawDetector[SampleRun])\n", + "da.masks.clear()\n", + "da.hist(two_theta=400, event_time_offset=1000).plot(norm='log', cmin=1.0e-3)" ] }, { @@ -179,9 +184,19 @@ "outputs": [], "source": [ "wf[DHKLList] = silicon_peaks_array()\n", - "da = wf.compute(TofDetector[SampleRun])\n", - "da = da.transform_coords(('dspacing',), graph=scn.conversion.graph.tof.elastic('tof'),)\n", - "ground_truth_peak_positions(da.hist(dspacing=dspacing, dim=da.dims).plot(), silicon_peaks_array())" + "wf[Filename[SampleRun]] = mcstas_silicon_medium_resolution()\n", + "\n", + "results = {}\n", + "for bank in (1, 2):\n", + " wf[DetectorBank] = bank\n", + " da = wf.compute(TofDetector[SampleRun])\n", + " results[bank] = (\n", + " da\n", + " .transform_coords(('dspacing',), graph=scn.conversion.graph.tof.elastic('tof'),)\n", + " .hist(dspacing=dspacing, dim=da.dims)\n", + " )\n", + "\n", + "ground_truth_peak_positions(pp.plot(results), silicon_peaks_array())" ] }, { @@ -201,9 +216,19 @@ "source": [ "wf = BeerModMcStasWorkflow()\n", "wf[Filename[SampleRun]] = mcstas_silicon_medium_resolution()\n", - "da = wf.compute(TofDetector[SampleRun])\n", - "da = da.transform_coords(('dspacing',), graph=scn.conversion.graph.tof.elastic('tof'),)\n", - "ground_truth_peak_positions(da.hist(dspacing=dspacing, dim=da.dims).plot(), silicon_peaks_array())" + "\n", + "results = {}\n", + "for bank in (1, 2):\n", + " wf[DetectorBank] = bank\n", + " da = wf.compute(TofDetector[SampleRun])\n", + " results[bank] = (\n", + " da\n", + " .transform_coords(('dspacing',), graph=scn.conversion.graph.tof.elastic('tof'),)\n", + " .hist(dspacing=dspacing, dim=da.dims)\n", + " )\n", + "\n", + "\n", + "ground_truth_peak_positions(pp.plot(results), silicon_peaks_array())" ] }, { @@ -223,11 +248,9 @@ "metadata": {}, "outputs": [], "source": [ - "da = da['bank2'].copy()\n", "da.masks.clear()\n", "da.bins.masks['too_far_from_center'] = ~da.bins.masks.pop('too_far_from_center')\n", - "\n", - "da.bins.concat().hist(two_theta=1000, event_time_offset=1000).plot(norm='log', cmin=1.0e-3)" + "da.bins.concat().hist(two_theta=400, event_time_offset=400).plot(norm='log', cmin=1.0e-3)" ] }, { @@ -254,8 +277,9 @@ "outputs": [], "source": [ "wf = BeerModMcStasWorkflowKnownPeaks()\n", + "wf[DetectorBank] = 1\n", "wf[Filename[SampleRun]] = mcstas_duplex(8)\n", - "wf.compute(RawDetector[SampleRun])['bank1'].hist(two_theta=1000, event_time_offset=1000).plot(norm='log', cmin=1.0e-2)" + "wf.compute(RawDetector[SampleRun]).hist(two_theta=400, event_time_offset=1000).plot(norm='log', cmin=1.0e-2)" ] }, { @@ -274,9 +298,18 @@ "outputs": [], "source": [ "wf[DHKLList] = duplex_peaks_array()\n", - "da = wf.compute(TofDetector[SampleRun])\n", - "da = da.transform_coords(('dspacing',), graph=scn.conversion.graph.tof.elastic('tof'),)\n", - "ground_truth_peak_positions(da.hist(dspacing=dspacing, dim=da.dims).plot(), duplex_peaks_array())" + "\n", + "results = {}\n", + "for bank in (1, 2):\n", + " wf[DetectorBank] = bank\n", + " da = wf.compute(TofDetector[SampleRun])\n", + " results[bank] = (\n", + " da\n", + " .transform_coords(('dspacing',), graph=scn.conversion.graph.tof.elastic('tof'),)\n", + " .hist(dspacing=dspacing, dim=da.dims)\n", + " )\n", + "\n", + "ground_truth_peak_positions(pp.plot(results), duplex_peaks_array())" ] }, { @@ -296,9 +329,18 @@ "source": [ "wf = BeerModMcStasWorkflow()\n", "wf[Filename[SampleRun]] = mcstas_duplex(8)\n", - "da = wf.compute(TofDetector[SampleRun])\n", - "da = da.transform_coords(('dspacing',), graph=scn.conversion.graph.tof.elastic('tof'),)\n", - "ground_truth_peak_positions(da.hist(dspacing=dspacing, dim=da.dims).plot(), duplex_peaks_array())" + "\n", + "results = {}\n", + "for bank in (1, 2):\n", + " wf[DetectorBank] = bank\n", + " da = wf.compute(TofDetector[SampleRun])\n", + " results[bank] = (\n", + " da\n", + " .transform_coords(('dspacing',), graph=scn.conversion.graph.tof.elastic('tof'),)\n", + " .hist(dspacing=dspacing, dim=da.dims)\n", + " )\n", + "\n", + "ground_truth_peak_positions(pp.plot(results), duplex_peaks_array())" ] }, { @@ -318,11 +360,10 @@ "metadata": {}, "outputs": [], "source": [ - "da = da['bank2'].copy()\n", "da.masks.clear()\n", "da.bins.masks['too_far_from_center'] = ~da.bins.masks.pop('too_far_from_center')\n", "\n", - "da.bins.concat().hist(two_theta=1000, event_time_offset=1000).plot(norm='log', cmin=1.0e-3)" + "da.bins.concat().hist(two_theta=400, event_time_offset=1000).plot(norm='log', cmin=1.0e-3)" ] }, { @@ -349,8 +390,9 @@ "outputs": [], "source": [ "wf = BeerModMcStasWorkflowKnownPeaks()\n", + "wf[DetectorBank] = 1\n", "wf[Filename[SampleRun]] = mcstas_duplex(9)\n", - "wf.compute(RawDetector[SampleRun])['bank1'].hist(two_theta=1000, event_time_offset=1000).plot(norm='log', cmin=1.0e-3)" + "wf.compute(RawDetector[SampleRun]).hist(two_theta=400, event_time_offset=1000).plot(norm='log', cmin=1.0e-3)" ] }, { @@ -369,9 +411,18 @@ "outputs": [], "source": [ "wf[DHKLList] = duplex_peaks_array()\n", - "da = wf.compute(TofDetector[SampleRun])\n", - "da = da.transform_coords(('dspacing',), graph=scn.conversion.graph.tof.elastic('tof'),)\n", - "ground_truth_peak_positions(da.hist(dspacing=dspacing, dim=da.dims).plot(), duplex_peaks_array())" + "\n", + "results = {}\n", + "for bank in (1, 2):\n", + " wf[DetectorBank] = bank\n", + " da = wf.compute(TofDetector[SampleRun])\n", + " results[bank] = (\n", + " da\n", + " .transform_coords(('dspacing',), graph=scn.conversion.graph.tof.elastic('tof'),)\n", + " .hist(dspacing=dspacing, dim=da.dims)\n", + " )\n", + "\n", + "ground_truth_peak_positions(pp.plot(results), duplex_peaks_array())" ] }, { @@ -391,9 +442,18 @@ "source": [ "wf = BeerModMcStasWorkflow()\n", "wf[Filename[SampleRun]] = mcstas_duplex(9)\n", - "da = wf.compute(TofDetector[SampleRun])\n", - "da = da.transform_coords(('dspacing',), graph=scn.conversion.graph.tof.elastic('tof'),)\n", - "ground_truth_peak_positions(da.hist(dspacing=dspacing, dim=da.dims).plot(), duplex_peaks_array())" + "\n", + "results = {}\n", + "for bank in (1, 2):\n", + " wf[DetectorBank] = bank\n", + " da = wf.compute(TofDetector[SampleRun])\n", + " results[bank] = (\n", + " da\n", + " .transform_coords(('dspacing',), graph=scn.conversion.graph.tof.elastic('tof'),)\n", + " .hist(dspacing=dspacing, dim=da.dims)\n", + " )\n", + "\n", + "ground_truth_peak_positions(pp.plot(results), duplex_peaks_array())" ] }, { @@ -413,11 +473,10 @@ "metadata": {}, "outputs": [], "source": [ - "da = da['bank2'].copy()\n", "da.masks.clear()\n", "da.bins.masks['too_far_from_center'] = ~da.bins.masks.pop('too_far_from_center')\n", "\n", - "da.bins.concat().hist(two_theta=1000, event_time_offset=1000).plot(norm='log', cmin=1.0e-3)" + "da.bins.concat().hist(two_theta=400, event_time_offset=1000).plot(norm='log', cmin=1.0e-3)" ] }, { @@ -444,8 +503,9 @@ "outputs": [], "source": [ "wf = BeerModMcStasWorkflowKnownPeaks()\n", + "wf[DetectorBank] = 1\n", "wf[Filename[SampleRun]] = mcstas_duplex(10)\n", - "wf.compute(RawDetector[SampleRun])['bank1'].hist(two_theta=1000, event_time_offset=1000).plot(norm='log', cmin=1.0e-3)" + "wf.compute(RawDetector[SampleRun]).hist(two_theta=400, event_time_offset=1000).plot(norm='log', cmin=1.0e-3)" ] }, { @@ -464,9 +524,18 @@ "outputs": [], "source": [ "wf[DHKLList] = duplex_peaks_array()\n", - "da = wf.compute(TofDetector[SampleRun])\n", - "da = da.transform_coords(('dspacing',), graph=scn.conversion.graph.tof.elastic('tof'),)\n", - "ground_truth_peak_positions(da.hist(dspacing=dspacing, dim=da.dims).plot(), duplex_peaks_array())" + "\n", + "results = {}\n", + "for bank in (1, 2):\n", + " wf[DetectorBank] = bank\n", + " da = wf.compute(TofDetector[SampleRun])\n", + " results[bank] = (\n", + " da\n", + " .transform_coords(('dspacing',), graph=scn.conversion.graph.tof.elastic('tof'),)\n", + " .hist(dspacing=dspacing, dim=da.dims)\n", + " )\n", + "\n", + "ground_truth_peak_positions(pp.plot(results), duplex_peaks_array())" ] }, { @@ -486,9 +555,18 @@ "source": [ "wf = BeerModMcStasWorkflow()\n", "wf[Filename[SampleRun]] = mcstas_duplex(10)\n", - "da = wf.compute(TofDetector[SampleRun])\n", - "da = da.transform_coords(('dspacing',), graph=scn.conversion.graph.tof.elastic('tof'),)\n", - "ground_truth_peak_positions(da.hist(dspacing=dspacing, dim=da.dims).plot(), duplex_peaks_array())" + "\n", + "results = {}\n", + "for bank in (1, 2):\n", + " wf[DetectorBank] = bank\n", + " da = wf.compute(TofDetector[SampleRun])\n", + " results[bank] = (\n", + " da\n", + " .transform_coords(('dspacing',), graph=scn.conversion.graph.tof.elastic('tof'),)\n", + " .hist(dspacing=dspacing, dim=da.dims)\n", + " )\n", + "\n", + "ground_truth_peak_positions(pp.plot(results), duplex_peaks_array())" ] }, { @@ -508,11 +586,10 @@ "metadata": {}, "outputs": [], "source": [ - "da = da['bank2'].copy()\n", "da.masks.clear()\n", "da.bins.masks['too_far_from_center'] = ~da.bins.masks.pop('too_far_from_center')\n", "\n", - "da.bins.concat().hist(two_theta=1000, event_time_offset=1000).plot(norm='log', cmin=1.0e-3)" + "da.bins.concat().hist(two_theta=400, event_time_offset=1000).plot(norm='log', cmin=1.0e-3)" ] }, { @@ -539,8 +616,9 @@ "outputs": [], "source": [ "wf = BeerModMcStasWorkflowKnownPeaks()\n", + "wf[DetectorBank] = 1\n", "wf[Filename[SampleRun]] = mcstas_duplex(16)\n", - "wf.compute(RawDetector[SampleRun])['bank1'].hist(two_theta=1000, event_time_offset=1000).plot(norm='log', cmin=1.0e-3)" + "wf.compute(RawDetector[SampleRun]).hist(two_theta=400, event_time_offset=1000).plot(norm='log', cmin=1.0e-3)" ] }, { @@ -559,9 +637,18 @@ "outputs": [], "source": [ "wf[DHKLList] = duplex_peaks_array()\n", - "da = wf.compute(TofDetector[SampleRun])\n", - "da = da.transform_coords(('dspacing',), graph=scn.conversion.graph.tof.elastic('tof'),)\n", - "ground_truth_peak_positions(da.hist(dspacing=dspacing, dim=da.dims).plot(), duplex_peaks_array())" + "\n", + "results = {}\n", + "for bank in (1, 2):\n", + " wf[DetectorBank] = bank\n", + " da = wf.compute(TofDetector[SampleRun])\n", + " results[bank] = (\n", + " da\n", + " .transform_coords(('dspacing',), graph=scn.conversion.graph.tof.elastic('tof'),)\n", + " .hist(dspacing=dspacing, dim=da.dims)\n", + " )\n", + "\n", + "ground_truth_peak_positions(pp.plot(results), duplex_peaks_array())" ] }, { @@ -581,9 +668,18 @@ "source": [ "wf = BeerModMcStasWorkflow()\n", "wf[Filename[SampleRun]] = mcstas_duplex(16)\n", - "da = wf.compute(TofDetector[SampleRun])\n", - "da = da.transform_coords(('dspacing',), graph=scn.conversion.graph.tof.elastic('tof'),)\n", - "ground_truth_peak_positions(da.hist(dspacing=dspacing, dim=da.dims).plot(), duplex_peaks_array())" + "\n", + "results = {}\n", + "for bank in (1, 2):\n", + " wf[DetectorBank] = bank\n", + " da = wf.compute(TofDetector[SampleRun])\n", + " results[bank] = (\n", + " da\n", + " .transform_coords(('dspacing',), graph=scn.conversion.graph.tof.elastic('tof'),)\n", + " .hist(dspacing=dspacing, dim=da.dims)\n", + " )\n", + "\n", + "ground_truth_peak_positions(pp.plot(results), duplex_peaks_array())" ] }, { @@ -603,11 +699,10 @@ "metadata": {}, "outputs": [], "source": [ - "da = da['bank2'].copy()\n", "da.masks.clear()\n", "da.bins.masks['too_far_from_center'] = ~da.bins.masks.pop('too_far_from_center')\n", "\n", - "da.bins.concat().hist(two_theta=1000, event_time_offset=1000).plot(norm='log', cmin=1.0e-3)" + "da.bins.concat().hist(two_theta=400, event_time_offset=1000).plot(norm='log', cmin=1.0e-3)" ] } ], @@ -627,7 +722,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.12.7" + "version": "3.11.13" } }, "nbformat": 4, diff --git a/src/ess/beer/clustering.py b/src/ess/beer/clustering.py index 40c848a4..e37395dc 100644 --- a/src/ess/beer/clustering.py +++ b/src/ess/beer/clustering.py @@ -1,31 +1,29 @@ import scipp as sc -from scippneutron.conversion.tof import dspacing_from_tof from scipy.signal import find_peaks, medfilt -from .conversions import t0_estimate, time_of_arrival -from .types import RawDetector, RunType, StreakClusteredData +from .conversions import tof_from_t0_estimate_graph +from .types import ( + GeometryCoordTransformGraph, + RawDetector, + RunType, + StreakClusteredData, +) -def cluster_events_by_streak(da: RawDetector[RunType]) -> StreakClusteredData[RunType]: - if isinstance(da, sc.DataGroup): - return sc.DataGroup({k: cluster_events_by_streak(v) for k, v in da.items()}) - da = da.copy(deep=False) +def cluster_events_by_streak( + da: RawDetector[RunType], gg: GeometryCoordTransformGraph +) -> StreakClusteredData[RunType]: + graph = tof_from_t0_estimate_graph(gg) - t = time_of_arrival( - da.coords['event_time_offset'], - da.coords['tc'].to(unit=da.coords['event_time_offset'].unit), - ) - approximate_t0 = t0_estimate( - da.coords['wavelength_estimate'], da.coords['L0'], da.coords['Ltotal'] - ).to(unit=t.unit) + da = da.transform_coords(['dspacing'], graph=graph) + da.bins.coords['coarse_d'] = da.bins.coords.pop('dspacing').to(unit='angstrom') - da.coords['coarse_d'] = dspacing_from_tof( - tof=t - approximate_t0, - Ltotal=da.coords['Ltotal'], - two_theta=da.coords['two_theta'], - ).to(unit='angstrom') + # We need to keep these coordinates after binning, + # adding them to the binned data coords achieves this. + for coord in ('two_theta', 'Ltotal'): + da.bins.coords[coord] = sc.bins_like(da, da.coords[coord]) - h = da.hist(coarse_d=1000) + h = da.bins.concat().hist(coarse_d=1000) i_peaks, _ = find_peaks( h.data.values, height=medfilt(h.values, kernel_size=99), distance=3 ) @@ -52,8 +50,10 @@ def cluster_events_by_streak(da: RawDetector[RunType]) -> StreakClusteredData[Ru ) ] has_peak = peaks.bin(coarse_d=filtered_valleys).bins.size().data - b = da.bin(coarse_d=filtered_valleys).assign_masks( - no_peak=has_peak != sc.scalar(1, unit=None) + b = ( + da.bins.concat() + .bin(coarse_d=filtered_valleys) + .assign_masks(no_peak=has_peak != sc.scalar(1, unit=None)) ) b = b.drop_coords(('coarse_d',)) b = b.bins.drop_coords(('coarse_d',)) diff --git a/src/ess/beer/conversions.py b/src/ess/beer/conversions.py index 0475b101..78fd6080 100644 --- a/src/ess/beer/conversions.py +++ b/src/ess/beer/conversions.py @@ -1,8 +1,10 @@ import scipp as sc import scipp.constants +from scippneutron.conversion import graph from .types import ( DHKLList, + GeometryCoordTransformGraph, ModulationPeriod, PulseLength, RawDetector, @@ -79,7 +81,7 @@ def _linear_regression_by_bin( return b1, b0 -def _compute_d( +def _compute_d_given_list_of_peaks( time_of_arrival: sc.Variable, theta: sc.Variable, dhkl_list: sc.Variable, @@ -92,15 +94,14 @@ def _compute_d( sinth = sc.sin(theta) t = time_of_arrival - d = sc.empty(dims=sinth.dims, shape=sinth.shape, unit=dhkl_list[0].unit) - d[:] = sc.scalar(float('nan'), unit=dhkl_list[0].unit) - dtfound = sc.empty(dims=sinth.dims, shape=sinth.shape, dtype='float64', unit=t.unit) - dtfound[:] = sc.scalar(float('nan'), unit=t.unit) + d = sc.full_like( + time_of_arrival, value=float('nan'), unit=dhkl_list[0].unit, dtype='float64' + ) + dtfound = sc.full_like(time_of_arrival, value=float('nan'), dtype='float64') const = (2 * sinth * L0 / (scipp.constants.h / scipp.constants.m_n)).to( unit=f'{time_of_arrival.unit}/angstrom' ) - for dhkl in dhkl_list: dt = sc.abs(t - dhkl * const) dt_in_range = dt < pulse_length / 2 @@ -145,9 +146,8 @@ def _tof_from_dhkl( c = (-2 * 1.0 / (scipp.constants.h / scipp.constants.m_n)).to( unit=f'{time_of_arrival.unit}/m/angstrom' ) - out = sc.sin(theta) - out *= c - out *= coarse_dhkl + out = c * coarse_dhkl + out *= sc.sin(theta) out *= Ltotal out += time_of_arrival out -= time0 @@ -161,11 +161,16 @@ def _tof_from_dhkl( return out +def geometry_graph() -> GeometryCoordTransformGraph: + return graph.beamline.beamline(scatter=True) + + def tof_from_known_dhkl_graph( mod_period: ModulationPeriod, pulse_length: PulseLength, time0: WavelengthDefinitionChopperDelay, dhkl_list: DHKLList, + gg: GeometryCoordTransformGraph, ) -> TofCoordTransformGraph: """Graph computing ``tof`` in modulation chopper modes using list of peak positions.""" @@ -182,7 +187,7 @@ def _compute_coarse_dspacing( The error happens because data arrays do not allow coordinates with dimensions not present on the data. """ - return _compute_d( + return _compute_d_given_list_of_peaks( time_of_arrival=time_of_arrival, theta=theta, pulse_length=pulse_length, @@ -191,6 +196,8 @@ def _compute_coarse_dspacing( ) return { + **gg, + **graph.tof.elastic("tof"), 'pulse_length': lambda: pulse_length, 'mod_period': lambda: mod_period, 'time0': lambda: time0, @@ -223,9 +230,13 @@ def _tof_from_t0( return time_of_arrival - t0 -def tof_from_t0_estimate_graph() -> TofCoordTransformGraph: +def tof_from_t0_estimate_graph( + gg: GeometryCoordTransformGraph, +) -> TofCoordTransformGraph: """Graph for computing ``tof`` in pulse shaping chopper modes.""" return { + **gg, + **graph.tof.elastic("tof"), 't0': t0_estimate, 'tof': _tof_from_t0, 'time_of_arrival': time_of_arrival, @@ -240,11 +251,13 @@ def compute_tof( convert_from_known_peaks_providers = ( + geometry_graph, tof_from_known_dhkl_graph, compute_tof, ) convert_pulse_shaping = ( + geometry_graph, tof_from_t0_estimate_graph, compute_tof, ) -providers = (compute_tof_in_each_cluster,) +providers = (compute_tof_in_each_cluster, geometry_graph) diff --git a/src/ess/beer/io.py b/src/ess/beer/io.py index 629c1d85..464961ec 100644 --- a/src/ess/beer/io.py +++ b/src/ess/beer/io.py @@ -1,13 +1,17 @@ # SPDX-License-Identifier: BSD-3-Clause # Copyright (c) 2025 Scipp contributors (https://github.com/scipp) +import re from pathlib import Path import h5py +import numpy as np import scipp as sc import scipp.constants from .types import ( + DetectorBank, Filename, + GeometryCoordTransformGraph, ModulationPeriod, RawDetector, SampleRun, @@ -16,6 +20,28 @@ ) +def _rotation_from_y_rotation_matrix(rot): + '''Assuming the rotation is around the y-axis + this function creates a rotation operator from the rotation matrix.''' + angle = np.atan2(rot[2, 0], rot[0, 0]) + return sc.spatial.rotation( + value=[ + 0.0, + np.sin(angle / 2), + 0.0, + np.cos(angle / 2), + ] + ) + + +def _find_h5(group: h5py.Group, matches): + for p in group.keys(): + if re.match(matches, p): + return group[p] + else: + raise RuntimeError(f'Could not find "{matches}" in {group}.') + + def _load_h5(group: h5py.Group | str, *paths: str): if isinstance(group, str): with h5py.File(group) as group: @@ -132,6 +158,12 @@ def _load_beer_mcstas(f, bank=1): positions['MCB'], positions['MCC'], ) + beam_rotation = _find_h5(f['/entry1/instrument/components'], '.*sourceMantid.*')[ + 'Rotation' + ] + detector_rotation = _find_h5( + f['/entry1/instrument/components'], f'.*nD_Mantid_?{bank}.*' + )['Rotation'] events = events[()] da = sc.DataArray( @@ -183,28 +215,59 @@ def _load_beer_mcstas(f, bank=1): list(map(float, da.coords.pop('position').value.split(' '))), unit='m' ) + da.coords.pop('n') da.coords['x'].unit = 'm' da.coords['y'].unit = 'm' da.coords['t'].unit = 's' - z = sc.norm(da.coords['detector_position'] - da.coords['sample_position']) + # Bin detector panel into rectangular "pixels" + # similar in size to the physical detector pixels. + da = da.bin( + y=sc.linspace('y', -0.5, 0.5, 501, unit='m'), + x=sc.linspace('x', -0.5, 0.5, 201, unit='m'), + ) + + # Compute the position of each pixel in the global coordinate system. + # The detector local coordinate system is rotatated by the detector rotation, + # and translated to the location of the detector in the global coordinate system. + da.coords['position'] = ( + da.coords['detector_position'] + + _rotation_from_y_rotation_matrix(detector_rotation) + * sc.spatial.as_vectors( + sc.midpoints(da.coords['x']), + sc.midpoints(da.coords['y']), + sc.scalar(0.0, unit='m'), + ) + # We need the dimension order of the positions to be the same + # as the dimension order of the binned data array. + ).transpose(da.dims) + L1 = sc.norm(da.coords['sample_position'] - da.coords['chopper_position']) - L2 = sc.sqrt(da.coords['x'] ** 2 + da.coords['y'] ** 2 + z**2) + L2 = sc.norm(da.coords['position'] - da.coords['sample_position']) - # Source is assumed to be at the origin - da.coords['L0'] = L1 + L2 + sc.norm(da.coords['chopper_position']) - da.coords['Ltotal'] = L1 + L2 - da.coords['two_theta'] = sc.acos( - (-da.coords['x'] if bank == 1 else da.coords['x']) / L2 + # Define the incident beam by rotating the z-axis by + # the rotation of the "source" in McStas. + incident_beam = L1 * ( + _rotation_from_y_rotation_matrix(beam_rotation) * sc.vector([0, 0, 1.0]) ) + # Create a source position that gives us the incident beam + # direction and length that we want. + # In practice this should be hardcoded or determined from + # some entry in the Nexus file. + da.coords['source_position'] = da.coords['sample_position'] - incident_beam - # Save some space - da.coords.pop('x') - da.coords.pop('y') - da.coords.pop('n') + # L0 is the total length of the instrument + da.coords['L0'] = L1 + L2 + sc.norm(da.coords['chopper_position']) - t = da.coords.pop('t') - da.coords['event_time_offset'] = t % sc.scalar(1 / 14, unit='s').to(unit=t.unit) + t = da.bins.coords['t'] + da.bins.coords['event_time_offset'] = t % sc.scalar(1 / 14, unit='s').to( + unit=t.unit + ) + # Estimate of the time the neutron passed the virtual source chopper. + # Used in pulse shaping mode to determine the wavelength. + # Used in modulation mode automatic-peak-finding reduction to estimate d. + # In practice this will probably be replaced by the regular tof workflow. + # But I'm not 100% sure. da.coords["tc"] = ( sc.constants.m_n / sc.constants.h @@ -212,48 +275,40 @@ def _load_beer_mcstas(f, bank=1): * da.coords['L0'].min().to(unit='angstrom') ).to(unit='s') - sc.scalar(1 / 14, unit='s') / 2 + del da.coords['x'] + del da.coords['y'] + # The binned t coordinate is kept because it can be useful + # to understand resolution and to debug tof estimation. return da -def load_beer_mcstas(f: str | Path | h5py.File) -> sc.DataGroup: +def _not_between(x, a, b): + return (x < a) | (b < x) + + +def load_beer_mcstas(f: str | Path | h5py.File, bank: int) -> sc.DataArray: '''Load beer McStas data from a file to a data group with one data array for each bank. ''' if isinstance(f, str | Path): with h5py.File(f) as ff: - return load_beer_mcstas(ff) - - return sc.DataGroup( - { - 'bank1': _load_beer_mcstas(f, bank=1), - 'bank2': _load_beer_mcstas(f, bank=2), - } - ) + return load_beer_mcstas(ff, bank=bank) - -def _not_between(x, a, b): - return (x < a) | (b < x) + return _load_beer_mcstas(f, bank=bank) def load_beer_mcstas_provider( - fname: Filename[SampleRun], two_theta_limits: TwoThetaLimits + fname: Filename[SampleRun], + bank: DetectorBank, + two_theta_limits: TwoThetaLimits, + graph: GeometryCoordTransformGraph, ) -> RawDetector[SampleRun]: - da = load_beer_mcstas(fname) - da = ( - sc.DataGroup( - { - k: v.assign_masks( - two_theta=_not_between(v.coords['two_theta'], *two_theta_limits) - ) - for k, v in da.items() - } - ) - if isinstance(da, sc.DataGroup) - else da.assign_masks( - two_theta=_not_between(da.coords['two_theta'], *two_theta_limits) - ) + da = load_beer_mcstas(fname, bank) + da = da.transform_coords(['two_theta'], graph=graph) + da = da.assign_masks( + two_theta=_not_between(da.coords['two_theta'], *two_theta_limits) ) - return RawDetector[SampleRun](da) + return da def mcstas_chopper_delay_from_mode( @@ -263,7 +318,7 @@ def mcstas_chopper_delay_from_mode( use in the docs currently. Eventually we will want to determine this from the chopper information in the files, but that information is not in the simulation output.''' - mode = next(iter(d.coords['mode'] for d in da.values())).value + mode = da.coords['mode'].value if mode in ('7', '8', '9', '10'): return sc.scalar(0.0024730158730158727, unit='s') if mode == '16': @@ -276,7 +331,7 @@ def mcstas_chopper_delay_from_mode_new_simulations( ) -> WavelengthDefinitionChopperDelay: '''Celine has a new simulation with some changes to the chopper placement(?). For those simulations we need to adapt the chopper delay values.''' - mode = next(iter(d.coords['mode'] for d in da.values())).value + mode = da.coords['mode'].value if mode == '7': return sc.scalar(0.001370158730158727, unit='s') if mode == '8': @@ -291,7 +346,7 @@ def mcstas_chopper_delay_from_mode_new_simulations( def mcstas_modulation_period_from_mode(da: RawDetector[SampleRun]) -> ModulationPeriod: - mode = next(iter(d.coords['mode'] for d in da.values())).value + mode = da.coords['mode'].value if mode in ('7', '8'): return sc.scalar(1.0 / (8 * 70), unit='s') if mode == '9': diff --git a/src/ess/beer/types.py b/src/ess/beer/types.py index 0c773643..144cfbb7 100644 --- a/src/ess/beer/types.py +++ b/src/ess/beer/types.py @@ -25,10 +25,12 @@ class StreakClusteredData(sciline.Scope[RunType, sc.DataArray], sc.DataArray): SampleRun = SampleRun TofDetector = TofDetector +DetectorBank = NewType('DetectorBank', int) TwoThetaLimits = NewType("TwoThetaLimits", tuple[sc.Variable, sc.Variable]) TofCoordTransformGraph = NewType("TofCoordTransformGraph", dict) +GeometryCoordTransformGraph = NewType("GeometryCoordTransformGraph", dict) PulseLength = NewType("PulseLength", sc.Variable) """Length of the neutron source pulse in time.""" diff --git a/tests/beer/mcstas_reduction_test.py b/tests/beer/mcstas_reduction_test.py index fd7aa75f..f19ee808 100644 --- a/tests/beer/mcstas_reduction_test.py +++ b/tests/beer/mcstas_reduction_test.py @@ -9,7 +9,7 @@ BeerModMcStasWorkflowKnownPeaks, ) from ess.beer.data import duplex_peaks_array, mcstas_duplex, mcstas_silicon_new_model -from ess.beer.types import DHKLList +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,12 +17,10 @@ def test_can_reduce_using_known_peaks_workflow(): wf = BeerModMcStasWorkflowKnownPeaks() wf[DHKLList] = duplex_peaks_array() + wf[DetectorBank] = 1 wf[Filename[SampleRun]] = mcstas_duplex(7) da = wf.compute(TofDetector[SampleRun]) - assert 'bank1' in da - assert 'bank2' in da - da = da['bank1'] - assert 'tof' in da.coords + assert 'tof' in da.bins.coords # assert dataarray has all coords required to compute dspacing da = da.transform_coords( ('dspacing',), @@ -40,10 +38,8 @@ 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 da = wf.compute(TofDetector[SampleRun]) - assert 'bank1' in da - assert 'bank2' in da - da = da['bank1'] da = da.transform_coords( ('dspacing',), graph=scn.conversion.graph.tof.elastic('tof'), @@ -60,11 +56,9 @@ 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 da = wf.compute(TofDetector[SampleRun]) - assert 'bank1' in da - assert 'bank2' in da - da = da['bank1'] - assert 'tof' in da.coords + assert 'tof' in da.bins.coords # assert dataarray has all coords required to compute dspacing da = da.transform_coords( ('dspacing',),