From 0bb475a4dfd2f849c8d09328f0b8e960b4e7b3e9 Mon Sep 17 00:00:00 2001 From: Johannes Kasimir Date: Fri, 16 Jan 2026 14:24:04 +0100 Subject: [PATCH 1/7] feat: bin mcstas data into artificial pixel grid after loading --- src/ess/beer/clustering.py | 19 +++++++++++++------ src/ess/beer/conversions.py | 17 ++++++++++------- src/ess/beer/io.py | 30 ++++++++++++++++++++---------- 3 files changed, 43 insertions(+), 23 deletions(-) diff --git a/src/ess/beer/clustering.py b/src/ess/beer/clustering.py index 40c848a4..82718afc 100644 --- a/src/ess/beer/clustering.py +++ b/src/ess/beer/clustering.py @@ -12,20 +12,25 @@ def cluster_events_by_streak(da: RawDetector[RunType]) -> StreakClusteredData[Ru da = da.copy(deep=False) t = time_of_arrival( - da.coords['event_time_offset'], - da.coords['tc'].to(unit=da.coords['event_time_offset'].unit), + da.bins.coords['event_time_offset'], + da.coords['tc'].to(unit=da.bins.coords['event_time_offset'].unit), ) approximate_t0 = t0_estimate( da.coords['wavelength_estimate'], da.coords['L0'], da.coords['Ltotal'] ).to(unit=t.unit) - da.coords['coarse_d'] = dspacing_from_tof( + da.bins.coords['coarse_d'] = dspacing_from_tof( tof=t - approximate_t0, Ltotal=da.coords['Ltotal'], two_theta=da.coords['two_theta'], ).to(unit='angstrom') - h = da.hist(coarse_d=1000) + # 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.bins.concat().hist(coarse_d=1000) i_peaks, _ = find_peaks( h.data.values, height=medfilt(h.values, kernel_size=99), distance=3 ) @@ -52,8 +57,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..110e3fb7 100644 --- a/src/ess/beer/conversions.py +++ b/src/ess/beer/conversions.py @@ -92,10 +92,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 = time_of_arrival.copy().to(dtype='float64') + # d.unit = dhkl_list[0].unit + # d[...] = sc.scalar(float('nan'), unit=d.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' @@ -145,9 +149,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 diff --git a/src/ess/beer/io.py b/src/ess/beer/io.py index 629c1d85..a866c2d7 100644 --- a/src/ess/beer/io.py +++ b/src/ess/beer/io.py @@ -183,28 +183,38 @@ 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']) + da = da.bin( + y=sc.linspace('y', -0.5, 0.5, 500, unit='m'), + x=sc.linspace('x', -0.5, 0.5, 500, unit='m'), + ) + da.coords['position'] = ( + da.coords['detector_position'] + + sc.spatial.as_vectors( + sc.scalar(0.0, unit='m'), + sc.midpoints(da.coords['y']), + sc.midpoints(-da.coords['x'] if bank == 1 else da.coords['x']), + ) + ).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 + (da.coords['position'] - da.coords['sample_position']).fields.z / L2 ) - # Save some space - da.coords.pop('x') - da.coords.pop('y') - da.coords.pop('n') - - 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 + ) da.coords["tc"] = ( sc.constants.m_n / sc.constants.h From a694d3e68cdcc679b8a7ce7cca8db8fd90408fc1 Mon Sep 17 00:00:00 2001 From: Johannes Kasimir Date: Fri, 16 Jan 2026 14:26:39 +0100 Subject: [PATCH 2/7] test: fix --- tests/beer/mcstas_reduction_test.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/beer/mcstas_reduction_test.py b/tests/beer/mcstas_reduction_test.py index fd7aa75f..ac62e236 100644 --- a/tests/beer/mcstas_reduction_test.py +++ b/tests/beer/mcstas_reduction_test.py @@ -22,7 +22,7 @@ def test_can_reduce_using_known_peaks_workflow(): 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',), @@ -64,7 +64,7 @@ def test_pulse_shaping_workflow(): 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',), From 021737d480eb1668b92510430d14d7ed331cdd0a Mon Sep 17 00:00:00 2001 From: Johannes Kasimir Date: Fri, 16 Jan 2026 14:29:39 +0100 Subject: [PATCH 3/7] fix --- src/ess/beer/conversions.py | 4 ---- 1 file changed, 4 deletions(-) diff --git a/src/ess/beer/conversions.py b/src/ess/beer/conversions.py index 110e3fb7..69bd5ea9 100644 --- a/src/ess/beer/conversions.py +++ b/src/ess/beer/conversions.py @@ -92,13 +92,9 @@ def _compute_d( sinth = sc.sin(theta) t = time_of_arrival - # d = time_of_arrival.copy().to(dtype='float64') - # d.unit = dhkl_list[0].unit - # d[...] = sc.scalar(float('nan'), unit=d.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( From 5271d6b7798bf9941408fefb69316f552dbcd09d Mon Sep 17 00:00:00 2001 From: Johannes Kasimir Date: Thu, 22 Jan 2026 12:20:53 +0100 Subject: [PATCH 4/7] refactoring and geometry fixes + adding DetectorBank domain type to select detector bank --- .../beer/beer_modulation_mcstas.ipynb | 189 +++++++++++++----- src/ess/beer/clustering.py | 33 ++- src/ess/beer/conversions.py | 24 ++- src/ess/beer/io.py | 133 ++++++++---- src/ess/beer/types.py | 2 + tests/beer/mcstas_reduction_test.py | 14 +- 6 files changed, 272 insertions(+), 123 deletions(-) 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 82718afc..e37395dc 100644 --- a/src/ess/beer/clustering.py +++ b/src/ess/beer/clustering.py @@ -1,29 +1,22 @@ 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.bins.coords['event_time_offset'], - da.coords['tc'].to(unit=da.bins.coords['event_time_offset'].unit), - ) - approximate_t0 = t0_estimate( - da.coords['wavelength_estimate'], da.coords['L0'], da.coords['Ltotal'] - ).to(unit=t.unit) - - da.bins.coords['coarse_d'] = dspacing_from_tof( - tof=t - approximate_t0, - Ltotal=da.coords['Ltotal'], - two_theta=da.coords['two_theta'], - ).to(unit='angstrom') + da = da.transform_coords(['dspacing'], graph=graph) + da.bins.coords['coarse_d'] = da.bins.coords.pop('dspacing').to(unit='angstrom') # We need to keep these coordinates after binning, # adding them to the binned data coords achieves this. diff --git a/src/ess/beer/conversions.py b/src/ess/beer/conversions.py index 69bd5ea9..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, @@ -100,7 +102,6 @@ def _compute_d( 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 @@ -160,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.""" @@ -181,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, @@ -190,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, @@ -222,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, @@ -239,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 a866c2d7..76e220c2 100644 --- a/src/ess/beer/io.py +++ b/src/ess/beer/io.py @@ -3,11 +3,14 @@ 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, @@ -119,6 +122,9 @@ def _load_beer_mcstas(f, bank=1): mca_pos, mcb_pos, mcc_pos, + source_rotation, + bank1_rotation, + bank2_rotation, ) = _load_h5( f, f'NXentry/NXdetector/bank{bank:02}_events_dat_list_p_x_y_n_id_t', @@ -131,6 +137,9 @@ def _load_beer_mcstas(f, bank=1): positions['MCA'], positions['MCB'], positions['MCC'], + '/entry1/instrument/components/0173_sourceMantid/Rotation', + '/entry1/instrument/components/0199_nD_Mantid1/Rotation', + '/entry1/instrument/components/0200_nD_Mantid2/Rotation', ) events = events[()] @@ -188,33 +197,84 @@ def _load_beer_mcstas(f, bank=1): da.coords['y'].unit = 'm' da.coords['t'].unit = 's' + # 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, 500, unit='m'), - x=sc.linspace('x', -0.5, 0.5, 500, unit='m'), + 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. + detector_rotation_angle = np.acos( + bank1_rotation[0, 0] if bank == 1 else bank2_rotation[0, 0] ) da.coords['position'] = ( da.coords['detector_position'] - + sc.spatial.as_vectors( - sc.scalar(0.0, unit='m'), + + ( + sc.spatial.rotation( + value=[ + 0.0, + np.sin(detector_rotation_angle / 2), + 0.0, + np.cos(detector_rotation_angle / 2), + ] + ) + if bank == 1 + else sc.spatial.rotation( + value=[ + 0.0, + np.sin(-detector_rotation_angle / 2), + 0.0, + np.cos(-detector_rotation_angle / 2), + ] + ) + ) + * sc.spatial.as_vectors( + sc.midpoints(da.coords['x']), sc.midpoints(da.coords['y']), - sc.midpoints(-da.coords['x'] if bank == 1 else da.coords['x']), + 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.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['position'] - da.coords['sample_position']).fields.z / L2 + # Define the incident beam by rotating the z-axis by + # the rotation of the "source" in McStas. + beam_rotation_angle = np.acos(source_rotation[0, 0]) + incident_beam = L1 * ( + sc.spatial.rotation( + value=[ + 0.0, + np.sin(beam_rotation_angle / 2), + 0.0, + np.cos(beam_rotation_angle / 2), + ] + ) + * 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 + + # L0 is the total length of the instrument + da.coords['L0'] = L1 + L2 + sc.norm(da.coords['chopper_position']) 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 @@ -222,48 +282,39 @@ 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 x, y, t coordinates are kept because they can be useful. 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( @@ -273,7 +324,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': @@ -286,7 +337,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': @@ -301,7 +352,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 ac62e236..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,11 +17,9 @@ 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.bins.coords # assert dataarray has all coords required to compute dspacing da = da.transform_coords( @@ -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,10 +56,8 @@ 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.bins.coords # assert dataarray has all coords required to compute dspacing da = da.transform_coords( From f19a9460fb2ad7d8fafbaefa7cae0315bd2ce666 Mon Sep 17 00:00:00 2001 From: Johannes Kasimir Date: Fri, 23 Jan 2026 09:40:39 +0100 Subject: [PATCH 5/7] fix: search for component matching name --- src/ess/beer/io.py | 24 +++++++++++++++--------- 1 file changed, 15 insertions(+), 9 deletions(-) diff --git a/src/ess/beer/io.py b/src/ess/beer/io.py index 76e220c2..4fa402fd 100644 --- a/src/ess/beer/io.py +++ b/src/ess/beer/io.py @@ -19,6 +19,14 @@ ) +def _find_h5(group: h5py.Group, matches): + for p in group.keys(): + if matches in 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: @@ -122,9 +130,6 @@ def _load_beer_mcstas(f, bank=1): mca_pos, mcb_pos, mcc_pos, - source_rotation, - bank1_rotation, - bank2_rotation, ) = _load_h5( f, f'NXentry/NXdetector/bank{bank:02}_events_dat_list_p_x_y_n_id_t', @@ -137,10 +142,13 @@ def _load_beer_mcstas(f, bank=1): positions['MCA'], positions['MCB'], positions['MCC'], - '/entry1/instrument/components/0173_sourceMantid/Rotation', - '/entry1/instrument/components/0199_nD_Mantid1/Rotation', - '/entry1/instrument/components/0200_nD_Mantid2/Rotation', ) + source_rotation = _find_h5(f['/entry1/instrument/components'], 'sourceMantid')[ + 'Rotation' + ] + detector_rotation = _find_h5(f['/entry1/instrument/components'], 'nD_Mantid')[ + 'Rotation' + ] events = events[()] da = sc.DataArray( @@ -207,9 +215,7 @@ def _load_beer_mcstas(f, bank=1): # 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. - detector_rotation_angle = np.acos( - bank1_rotation[0, 0] if bank == 1 else bank2_rotation[0, 0] - ) + detector_rotation_angle = np.acos(detector_rotation[0, 0]) da.coords['position'] = ( da.coords['detector_position'] + ( From 6c8e0bc45ceabdf6e3f35c65bae9e11952a6a1a1 Mon Sep 17 00:00:00 2001 From: Johannes Kasimir Date: Fri, 23 Jan 2026 10:20:59 +0100 Subject: [PATCH 6/7] fix --- src/ess/beer/io.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/ess/beer/io.py b/src/ess/beer/io.py index 4fa402fd..c59dc1fd 100644 --- a/src/ess/beer/io.py +++ b/src/ess/beer/io.py @@ -290,7 +290,8 @@ def _load_beer_mcstas(f, bank=1): del da.coords['x'] del da.coords['y'] - # The binned x, y, t coordinates are kept because they can be useful. + # The binned t coordinate is kept because it can be useful + # to understand resolution and to debug tof estimation. return da From 2c31563c483cf5a5a6914ac572b50e0dd070e06d Mon Sep 17 00:00:00 2001 From: Johannes Kasimir Date: Fri, 23 Jan 2026 11:47:39 +0100 Subject: [PATCH 7/7] fix: remove cases --- src/ess/beer/io.py | 57 ++++++++++++++++++---------------------------- 1 file changed, 22 insertions(+), 35 deletions(-) diff --git a/src/ess/beer/io.py b/src/ess/beer/io.py index c59dc1fd..464961ec 100644 --- a/src/ess/beer/io.py +++ b/src/ess/beer/io.py @@ -1,5 +1,6 @@ # SPDX-License-Identifier: BSD-3-Clause # Copyright (c) 2025 Scipp contributors (https://github.com/scipp) +import re from pathlib import Path import h5py @@ -19,9 +20,23 @@ ) +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 matches in p: + if re.match(matches, p): return group[p] else: raise RuntimeError(f'Could not find "{matches}" in {group}.') @@ -143,12 +158,12 @@ def _load_beer_mcstas(f, bank=1): positions['MCB'], positions['MCC'], ) - source_rotation = _find_h5(f['/entry1/instrument/components'], 'sourceMantid')[ - 'Rotation' - ] - detector_rotation = _find_h5(f['/entry1/instrument/components'], 'nD_Mantid')[ + 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( @@ -215,28 +230,9 @@ def _load_beer_mcstas(f, bank=1): # 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. - detector_rotation_angle = np.acos(detector_rotation[0, 0]) da.coords['position'] = ( da.coords['detector_position'] - + ( - sc.spatial.rotation( - value=[ - 0.0, - np.sin(detector_rotation_angle / 2), - 0.0, - np.cos(detector_rotation_angle / 2), - ] - ) - if bank == 1 - else sc.spatial.rotation( - value=[ - 0.0, - np.sin(-detector_rotation_angle / 2), - 0.0, - np.cos(-detector_rotation_angle / 2), - ] - ) - ) + + _rotation_from_y_rotation_matrix(detector_rotation) * sc.spatial.as_vectors( sc.midpoints(da.coords['x']), sc.midpoints(da.coords['y']), @@ -251,17 +247,8 @@ def _load_beer_mcstas(f, bank=1): # Define the incident beam by rotating the z-axis by # the rotation of the "source" in McStas. - beam_rotation_angle = np.acos(source_rotation[0, 0]) incident_beam = L1 * ( - sc.spatial.rotation( - value=[ - 0.0, - np.sin(beam_rotation_angle / 2), - 0.0, - np.cos(beam_rotation_angle / 2), - ] - ) - * sc.vector([0, 0, 1.0]) + _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.