From f98b1cfe3c92ad9a2ebafa3a2eaa0175ba4f5482 Mon Sep 17 00:00:00 2001 From: Johannes Kasimir Date: Fri, 24 Oct 2025 11:01:22 +0200 Subject: [PATCH 1/5] fix: convert unlimited TOA to expected range --- src/ess/dream/io/geant4.py | 35 +++++++++++++++++++++++++++++------ 1 file changed, 29 insertions(+), 6 deletions(-) diff --git a/src/ess/dream/io/geant4.py b/src/ess/dream/io/geant4.py index ae1776d6..db850b2d 100644 --- a/src/ess/dream/io/geant4.py +++ b/src/ess/dream/io/geant4.py @@ -196,6 +196,27 @@ def _to_edges(centers: sc.Variable) -> sc.Variable: ) +def _rebin_to_tofrange(da): + '''Rebins a monitor TOA histogram to the range of values + expected from a real instrument at ESS. + That is, to the range [0, 1/14] s. + + Strategy: + 1. Create a grid that alingns with the pulse times. + 2. Rebin TOA on that grid. + 3. Fold pulses into new dimension and sum over pulses. + ''' + dim = da.dim + unit = da.coords[dim].unit + period = (1.0 / sc.scalar(14.0, unit='Hz')).to(unit=unit) + N = sc.sum(da.coords[dim] < period).value + K = (da.coords[dim].max() // period + 1).to(dtype='int') + grid = sc.linspace(dim, sc.scalar(0.0, unit=unit), K * period, K * (N - 1) + 1) + out = da.rebin(tof=grid).fold(dim, sizes={'_': K, dim: N - 1}).sum('_') + out.coords[dim] = sc.linspace(dim, sc.scalar(0.0, unit=unit), period, N) + return out + + def load_mcstas_monitor( file_path: MonitorFilename[RunType], position: CaveMonitorPosition, @@ -221,14 +242,16 @@ def load_mcstas_monitor( tof, counts, err = np.loadtxt(file_path, usecols=(0, 1, 2), unpack=True) tof = _to_edges(sc.array(dims=["tof"], values=tof, unit="us")) + data = sc.DataArray( + sc.array(dims=["tof"], values=counts, variances=err**2, unit="counts"), + coords={ + "tof": tof, + }, + ) + data = _rebin_to_tofrange(data) return NeXusComponent[CaveMonitor, RunType]( sc.DataGroup( - data=sc.DataArray( - sc.array(dims=["tof"], values=counts, variances=err**2, unit="counts"), - coords={ - "tof": tof, - }, - ), + data=data, position=position, ) ) From 8994b88e039fe926cc82a98f07b269a824397b76 Mon Sep 17 00:00:00 2001 From: Johannes Kasimir Date: Mon, 27 Oct 2025 12:06:08 +0100 Subject: [PATCH 2/5] fix: add masking in monitor correction --- src/ess/powder/correction.py | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/src/ess/powder/correction.py b/src/ess/powder/correction.py index bfbc6d94..c760c4f8 100644 --- a/src/ess/powder/correction.py +++ b/src/ess/powder/correction.py @@ -59,11 +59,15 @@ def normalize_by_monitor_histogram( ) lut = sc.lookup(norm, dim="wavelength") if detector.bins is None: - result = ( - detector / lut[sc.midpoints(detector.coords['wavelength'], dim='dspacing')] - ) + c = lut[sc.midpoints(detector.coords['wavelength'], dim='dspacing')] + result = detector / c + result.masks['monitor_intensity_is_zero'] = c == sc.scalar(0.0, unit=c.unit) else: - result = detector.bins / lut + c = lut(detector.bins.coords['wavelength']) + result = detector / c + result.bins.masks['monitor_intensity_is_zero'] = c == sc.scalar( + 0.0, unit=c.unit + ) return ScaledCountsDspacing[RunType](result) From efdb5ef5b7b8bc9ad762680c433d6baa60e3e67c Mon Sep 17 00:00:00 2001 From: Johannes Kasimir Date: Mon, 27 Oct 2025 17:11:31 +0100 Subject: [PATCH 3/5] fix: use dim coord not 'tof' --- src/ess/dream/io/geant4.py | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/src/ess/dream/io/geant4.py b/src/ess/dream/io/geant4.py index db850b2d..bd74e6e1 100644 --- a/src/ess/dream/io/geant4.py +++ b/src/ess/dream/io/geant4.py @@ -212,7 +212,7 @@ def _rebin_to_tofrange(da): N = sc.sum(da.coords[dim] < period).value K = (da.coords[dim].max() // period + 1).to(dtype='int') grid = sc.linspace(dim, sc.scalar(0.0, unit=unit), K * period, K * (N - 1) + 1) - out = da.rebin(tof=grid).fold(dim, sizes={'_': K, dim: N - 1}).sum('_') + out = da.rebin({dim: grid}).fold(dim, sizes={'_': K, dim: N - 1}).sum('_') out.coords[dim] = sc.linspace(dim, sc.scalar(0.0, unit=unit), period, N) return out @@ -244,9 +244,7 @@ def load_mcstas_monitor( tof = _to_edges(sc.array(dims=["tof"], values=tof, unit="us")) data = sc.DataArray( sc.array(dims=["tof"], values=counts, variances=err**2, unit="counts"), - coords={ - "tof": tof, - }, + coords={"tof": tof}, ) data = _rebin_to_tofrange(data) return NeXusComponent[CaveMonitor, RunType]( From 72a5df3c048742afe26db25a2f6c1efece927a7c Mon Sep 17 00:00:00 2001 From: Johannes Kasimir Date: Mon, 27 Oct 2025 17:11:53 +0100 Subject: [PATCH 4/5] test: add masks to expected --- tests/powder/correction_test.py | 13 +++++++++++-- 1 file changed, 11 insertions(+), 2 deletions(-) diff --git a/tests/powder/correction_test.py b/tests/powder/correction_test.py index f1d77492..c51ebbeb 100644 --- a/tests/powder/correction_test.py +++ b/tests/powder/correction_test.py @@ -347,7 +347,12 @@ def test_normalize_by_monitor_histogram_expected_results(): monitor=WavelengthMonitor[SampleRun, CaveMonitor](monitor), uncertainty_broadcast_mode=UncertaintyBroadcastMode.fail, ) - expected = ScaledCountsDspacing[SampleRun](detector / monitor.data) + # Simple way to get all false array same shape as detector data. + # Is there a better way? + all_false = sc.isinf(detector.bins.data) + expected = ScaledCountsDspacing[SampleRun]( + detector / monitor.data + ).bins.assign_masks({'monitor_intensity_is_zero': all_false}) sc.testing.assert_identical(normalized, expected) @@ -362,12 +367,16 @@ def test_normalize_by_monitor_histogram_ignores_monitor_values_out_of_range(): 'wavelength': sc.array(dims=['wavelength'], values=[0.0, 3, 4], unit='Å') }, ) + # Simple way to get all false array same shape as detector data. + all_false = sc.isinf(detector.bins.data) normalized = normalize_by_monitor_histogram( CountsWavelength[SampleRun](detector), monitor=WavelengthMonitor[SampleRun, CaveMonitor](monitor), uncertainty_broadcast_mode=UncertaintyBroadcastMode.fail, ) - expected = ScaledCountsDspacing[SampleRun](detector / sc.scalar(4.0, unit='counts')) + expected = ScaledCountsDspacing[SampleRun]( + detector / sc.scalar(4.0, unit='counts') + ).bins.assign_masks({'monitor_intensity_is_zero': all_false}) sc.testing.assert_identical(normalized, expected) From fa6731b00c4c1dd94bed071055d51c6e67867d6e Mon Sep 17 00:00:00 2001 From: Johannes Kasimir Date: Tue, 28 Oct 2025 13:53:46 +0100 Subject: [PATCH 5/5] test: loaded monitor has correct 'tof' range --- tests/dream/io/geant4_test.py | 27 +++++++++++++++++++++++++-- 1 file changed, 25 insertions(+), 2 deletions(-) diff --git a/tests/dream/io/geant4_test.py b/tests/dream/io/geant4_test.py index 4d17c973..428fb475 100644 --- a/tests/dream/io/geant4_test.py +++ b/tests/dream/io/geant4_test.py @@ -10,8 +10,19 @@ import scipp.testing import scippnexus as snx -from ess.dream import DreamGeant4ProtonChargeWorkflow, data, load_geant4_csv -from ess.powder.types import Filename, NeXusComponent, NeXusDetectorName, SampleRun +from ess.dream import DreamGeant4ProtonChargeWorkflow, data, io, load_geant4_csv +from ess.powder.types import ( + CaveMonitorPosition, + Filename, + NeXusComponent, + NeXusDetectorName, + SampleRun, +) + + +@pytest.fixture(scope="module") +def monitor_file_path(): + return data.simulated_monitor_diamond_sample() @pytest.fixture(scope="module") @@ -184,3 +195,15 @@ def test_geant4_in_workflow(file_path, file): expected = load_geant4_csv(file)["instrument"]["mantle"]["events"] expected.coords["detector"] = sc.scalar("mantle") sc.testing.assert_identical(detector, expected) + + +def test_geant4_monitor_in_expected_range(monitor_file_path): + dummy_position = CaveMonitorPosition(sc.vector([0, 0, 0.0], unit='m')) + mon = io.geant4.load_mcstas_monitor(monitor_file_path, dummy_position)['data'] + assert ( + mon.coords['tof'] >= sc.scalar(0.0, unit='s').to(unit=mon.coords['tof'].unit) + ).all() + assert ( + mon.coords['tof'] + <= sc.scalar(1 / 14.0, unit='s').to(unit=mon.coords['tof'].unit) + ).all()