From 826279bcf99724dd784f0cae4d077f92fe9a1030 Mon Sep 17 00:00:00 2001 From: Johannes Kasimir Date: Mon, 27 Apr 2026 15:40:41 +0200 Subject: [PATCH 1/8] feat: add monitor to estia workflow --- .../src/ess/estia/corrections.py | 51 ++++++++++++++++++- .../essreflectometry/src/ess/estia/mcstas.py | 1 + .../essreflectometry/src/ess/estia/types.py | 8 ++- .../src/ess/estia/workflow.py | 6 ++- .../tests/estia/mcstas_data_test.py | 26 ++++++++++ 5 files changed, 88 insertions(+), 4 deletions(-) diff --git a/packages/essreflectometry/src/ess/estia/corrections.py b/packages/essreflectometry/src/ess/estia/corrections.py index a83dc2f49..9fc222e97 100644 --- a/packages/essreflectometry/src/ess/estia/corrections.py +++ b/packages/essreflectometry/src/ess/estia/corrections.py @@ -1,6 +1,8 @@ # SPDX-License-Identifier: BSD-3-Clause # Copyright (c) 2025 Scipp contributors (https://github.com/scipp) +import ess.reduce import scipp as sc +from ess.reduce.uncertainty import UncertaintyBroadcastMode from ..reflectometry.conversions import ( add_proton_current_coord, @@ -21,6 +23,47 @@ ) from .conversions import add_coords from .maskings import add_masks +from .types import WavelengthMonitor + + +def normalize_by_monitor_histogram( + detector: ReducibleData[RunType], + *, + monitor: WavelengthMonitor[RunType], +) -> ReducibleData[RunType]: + """Normalize detector data by a histogrammed monitor. + + The detector is normalized according to + + .. math:: + + d_i^\\text{Norm} = \\frac{d_i}{m_i} \\Delta \\lambda_i + + Parameters + ---------- + detector: + Input event data in wavelength. + monitor: + A histogrammed monitor in wavelength. + uncertainty_broadcast_mode: + Choose how uncertainties of the monitor are broadcast to the sample data. + + Returns + ------- + : + `detector` normalized by a monitor. + + See also + -------- + ess.reduce.normalization.normalize_by_monitor_histogram: + For details and the actual implementation. + """ + return ess.reduce.normalization.normalize_by_monitor_histogram( + detector=detector, + monitor=monitor, + uncertainty_broadcast_mode=UncertaintyBroadcastMode.drop, + skip_range_check=False, + ) def add_coords_masks_and_apply_corrections( @@ -30,6 +73,7 @@ def add_coords_masks_and_apply_corrections( bdlim: BeamDivergenceLimits, wbins: WavelengthBins, proton_current: ProtonCurrent[RunType], + monitor: WavelengthMonitor[RunType], graph: CoordTransformationGraph[RunType], corrections_to_apply: CorrectionsToApply, ) -> ReducibleData[RunType]: @@ -45,7 +89,10 @@ def add_coords_masks_and_apply_corrections( da = add_proton_current_mask(da) for correction in corrections_to_apply: - da = correction(da) + if correction == 'monitor': + da = normalize_by_monitor_histogram(da, monitor=monitor) + else: + da = correction(da) return ReducibleData[RunType](da) @@ -61,6 +108,6 @@ def assume_time_series_constant_with_zero_default_value_if_empty(da: sc.DataArra return da.mean() if len(da) > 0 else sc.scalar(0.0, unit=da.unit) -default_corrections = {correct_by_proton_current, correct_by_footprint} +default_corrections = {correct_by_footprint, 'proton_current', 'monitor'} providers = (add_coords_masks_and_apply_corrections,) diff --git a/packages/essreflectometry/src/ess/estia/mcstas.py b/packages/essreflectometry/src/ess/estia/mcstas.py index ae21711bb..072c6fdbb 100644 --- a/packages/essreflectometry/src/ess/estia/mcstas.py +++ b/packages/essreflectometry/src/ess/estia/mcstas.py @@ -302,6 +302,7 @@ def use_mcstas_wavelengths_instead_of_estimates_from_time_of_arrival( zlims=zlims, bdlim=bdlim, wbins=wbins, + monitor=None, proton_current=proton_current, graph={ **graph, diff --git a/packages/essreflectometry/src/ess/estia/types.py b/packages/essreflectometry/src/ess/estia/types.py index 8acdf9fcd..efcfb64fd 100644 --- a/packages/essreflectometry/src/ess/estia/types.py +++ b/packages/essreflectometry/src/ess/estia/types.py @@ -1,8 +1,14 @@ # Copyright (c) 2025 Scipp contributors (https://github.com/scipp) -from typing import NewType +from typing import NewType, TypeAlias import scipp as sc +from ess.reduce.nexus.types import CaveMonitor +from ess.reduce.unwrap.types import WavelengthMonitor + +from ess.reflectometry.types import RunType WavelengthResolution = NewType("WavelengthResolution", sc.Variable) AngularResolution = NewType("AngularResolution", sc.Variable) SampleSizeResolution = NewType("SampleSizeResolution", sc.Variable) + +WavelengthMonitor: TypeAlias = WavelengthMonitor[RunType, CaveMonitor] diff --git a/packages/essreflectometry/src/ess/estia/workflow.py b/packages/essreflectometry/src/ess/estia/workflow.py index 847ad883e..1fbbc33bf 100644 --- a/packages/essreflectometry/src/ess/estia/workflow.py +++ b/packages/essreflectometry/src/ess/estia/workflow.py @@ -28,6 +28,7 @@ normalization, orso, ) +from .types import WavelengthMonitor _general_providers = ( *reflectometry_providers, @@ -70,10 +71,12 @@ def mcstas_default_parameters() -> dict: ), SampleRotationOffset[RunType]: sc.scalar(0.0, unit='deg'), CorrectionsToApply: corrections.default_corrections + - {'monitor'} - {correct_by_proton_current}, LookupTableRelativeErrorThreshold: { "multiblade_detector": 0.06, }, + WavelengthMonitor[RunType]: None, } @@ -82,11 +85,12 @@ def default_parameters() -> dict: return { NeXusDetectorName: "multiblade_detector", SampleRotationOffset[RunType]: sc.scalar(0.0, unit='deg'), - CorrectionsToApply: corrections.default_corrections, + CorrectionsToApply: corrections.default_corrections - {'monitor'}, DetectorSpatialResolution: 0.0025 * sc.units.m, LookupTableRelativeErrorThreshold: { "multiblade_detector": float('inf'), }, + WavelengthMonitor[RunType]: None, } diff --git a/packages/essreflectometry/tests/estia/mcstas_data_test.py b/packages/essreflectometry/tests/estia/mcstas_data_test.py index 0edf9e398..24543c3d4 100644 --- a/packages/essreflectometry/tests/estia/mcstas_data_test.py +++ b/packages/essreflectometry/tests/estia/mcstas_data_test.py @@ -8,6 +8,9 @@ import pytest import sciline import scipp as sc +import scipp.testing +from ess.reduce.normalization import normalize_by_monitor_histogram +from ess.reduce.uncertainty import UncertaintyBroadcastMode from orsopy import fileio from ess.estia import EstiaMcStasWorkflow @@ -19,9 +22,11 @@ from ess.estia.mcstas import ( use_mcstas_wavelengths_instead_of_estimates_from_time_of_arrival, ) +from ess.estia.types import WavelengthMonitor from ess.reflectometry import orso from ess.reflectometry.types import ( BeamDivergenceLimits, + CorrectionsToApply, Filename, LookupTableFilename, ProtonCurrent, @@ -97,6 +102,27 @@ def test_compute_reducible_data(estia_mcstas_pipeline: sciline.Pipeline): assert 'Q' in da.bins.coords +def test_compute_reducible_data_with_monitor(estia_mcstas_pipeline: sciline.Pipeline): + wf = estia_mcstas_pipeline + wf[Filename[SampleRun]] = estia_mcstas_sample_run(11) + without_monitor = wf.compute(ReducibleData[SampleRun]) + wf[WavelengthMonitor[SampleRun]] = sc.DataArray( + sc.array(dims=['wavelength'], values=[30.0], variances=[1.0]), + coords={'wavelength': sc.linspace('wavelength', 0, 15, 2, unit='angstrom')}, + ) + corrections = wf.compute(CorrectionsToApply) + wf[CorrectionsToApply] = {*corrections, 'monitor'} + with_monitor = wf.compute(ReducibleData[SampleRun]) + scipp.testing.assert_allclose( + normalize_by_monitor_histogram( + without_monitor, + monitor=wf.compute(WavelengthMonitor[SampleRun]), + uncertainty_broadcast_mode=UncertaintyBroadcastMode.drop, + ), + with_monitor, + ) + + def test_can_compute_reflectivity_curve_exact_wavelengths( estia_mcstas_pipeline: sciline.Pipeline, ): From 98e1efb84cc4d9f43e312a9593e4943efc90e39f Mon Sep 17 00:00:00 2001 From: Johannes Kasimir Date: Mon, 27 Apr 2026 15:46:59 +0200 Subject: [PATCH 2/8] fix --- packages/essreflectometry/src/ess/estia/workflow.py | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/packages/essreflectometry/src/ess/estia/workflow.py b/packages/essreflectometry/src/ess/estia/workflow.py index 1fbbc33bf..515daa390 100644 --- a/packages/essreflectometry/src/ess/estia/workflow.py +++ b/packages/essreflectometry/src/ess/estia/workflow.py @@ -70,9 +70,9 @@ def mcstas_default_parameters() -> dict: sc.scalar(0.75, unit='deg'), ), SampleRotationOffset[RunType]: sc.scalar(0.0, unit='deg'), - CorrectionsToApply: corrections.default_corrections - - {'monitor'} - - {correct_by_proton_current}, + CorrectionsToApply: ( + corrections.default_corrections - {'monitor', correct_by_proton_current} + ), LookupTableRelativeErrorThreshold: { "multiblade_detector": 0.06, }, @@ -90,6 +90,9 @@ def default_parameters() -> dict: LookupTableRelativeErrorThreshold: { "multiblade_detector": float('inf'), }, + # The monitor is missing from the Nexus files, + # so to be able to load Nexus files anyway + # this is set to None. WavelengthMonitor[RunType]: None, } From b60433ac6498c4c87d3f3a2369cca86e20dda9d3 Mon Sep 17 00:00:00 2001 From: Johannes Kasimir Date: Tue, 28 Apr 2026 15:33:56 +0200 Subject: [PATCH 3/8] fix: handle both proton charge with and without variance --- .../essreflectometry/src/ess/amor/workflow.py | 8 +--- .../src/ess/estia/corrections.py | 11 ++---- .../essreflectometry/src/ess/estia/mcstas.py | 7 ++-- .../src/ess/estia/workflow.py | 10 ++--- .../src/ess/reflectometry/conversions.py | 38 ------------------- .../src/ess/reflectometry/corrections.py | 20 +++++++++- .../tests/amor/pipeline_test.py | 2 - .../tests/estia/mcstas_data_test.py | 23 +++++++++++ 8 files changed, 52 insertions(+), 67 deletions(-) diff --git a/packages/essreflectometry/src/ess/amor/workflow.py b/packages/essreflectometry/src/ess/amor/workflow.py index 576d73596..7ef3ccedc 100644 --- a/packages/essreflectometry/src/ess/amor/workflow.py +++ b/packages/essreflectometry/src/ess/amor/workflow.py @@ -1,8 +1,4 @@ # Copyright (c) 2025 Scipp contributors (https://github.com/scipp) -from ..reflectometry.conversions import ( - add_proton_current_coord, - add_proton_current_mask, -) from ..reflectometry.corrections import correct_by_footprint, correct_by_proton_current from ..reflectometry.types import ( BeamDivergenceLimits, @@ -38,9 +34,7 @@ def add_coords_masks_and_apply_corrections( # For some older Amor files there are no entries in the proton current log if len(proton_current) != 0: - da = add_proton_current_coord(da, proton_current) - da = add_proton_current_mask(da) - da = correct_by_proton_current(da) + da = correct_by_proton_current(da, proton_current=proton_current) return ReducibleData[RunType](da) diff --git a/packages/essreflectometry/src/ess/estia/corrections.py b/packages/essreflectometry/src/ess/estia/corrections.py index 9fc222e97..60ea043bf 100644 --- a/packages/essreflectometry/src/ess/estia/corrections.py +++ b/packages/essreflectometry/src/ess/estia/corrections.py @@ -4,10 +4,6 @@ import scipp as sc from ess.reduce.uncertainty import UncertaintyBroadcastMode -from ..reflectometry.conversions import ( - add_proton_current_coord, - add_proton_current_mask, -) from ..reflectometry.corrections import correct_by_proton_current from ..reflectometry.types import ( BeamDivergenceLimits, @@ -84,13 +80,12 @@ def add_coords_masks_and_apply_corrections( da = add_coords(da, graph) da = add_masks(da, ylim, zlims, bdlim, wbins) - if len(proton_current) != 0: - da = add_proton_current_coord(da, proton_current) - da = add_proton_current_mask(da) - for correction in corrections_to_apply: if correction == 'monitor': da = normalize_by_monitor_histogram(da, monitor=monitor) + elif correction == 'proton_current': + if len(proton_current) != 0: + da = correct_by_proton_current(da, proton_current=proton_current) else: da = correction(da) diff --git a/packages/essreflectometry/src/ess/estia/mcstas.py b/packages/essreflectometry/src/ess/estia/mcstas.py index 072c6fdbb..b4bcdbe01 100644 --- a/packages/essreflectometry/src/ess/estia/mcstas.py +++ b/packages/essreflectometry/src/ess/estia/mcstas.py @@ -223,10 +223,11 @@ def load_mcstas( da.bins.coords['event_time_zero'] = ( sc.scalar(0, unit='s') * da.bins.coords['t'] - ).to(unit='ns') + ).to(unit='ns', dtype='int64') + sc.datetime(0, unit='ns') da.bins.coords['event_time_offset'] = ( - sc.scalar(1, unit='s') * da.bins.coords['t'] - ).to(unit='ns') % sc.scalar(1 / 14, unit='s').to(unit='ns') + (sc.scalar(1, unit='s') * da.bins.coords['t']).to(unit='ns') + % sc.scalar(1 / 14, unit='s').to(unit='ns') + ).to(dtype='int32') da.bins.coords.pop('t') if 'L' in da.bins.coords: diff --git a/packages/essreflectometry/src/ess/estia/workflow.py b/packages/essreflectometry/src/ess/estia/workflow.py index 515daa390..796cb1112 100644 --- a/packages/essreflectometry/src/ess/estia/workflow.py +++ b/packages/essreflectometry/src/ess/estia/workflow.py @@ -8,7 +8,6 @@ from ..reflectometry import providers as reflectometry_providers from ..reflectometry import supermirror -from ..reflectometry.corrections import correct_by_proton_current from ..reflectometry.types import ( BeamDivergenceLimits, CorrectionsToApply, @@ -71,7 +70,7 @@ def mcstas_default_parameters() -> dict: ), SampleRotationOffset[RunType]: sc.scalar(0.0, unit='deg'), CorrectionsToApply: ( - corrections.default_corrections - {'monitor', correct_by_proton_current} + corrections.default_corrections - {'monitor', 'proton_current'} ), LookupTableRelativeErrorThreshold: { "multiblade_detector": 0.06, @@ -85,15 +84,12 @@ def default_parameters() -> dict: return { NeXusDetectorName: "multiblade_detector", SampleRotationOffset[RunType]: sc.scalar(0.0, unit='deg'), - CorrectionsToApply: corrections.default_corrections - {'monitor'}, + CorrectionsToApply: corrections.default_corrections + - {'monitor', 'proton_current'}, DetectorSpatialResolution: 0.0025 * sc.units.m, LookupTableRelativeErrorThreshold: { "multiblade_detector": float('inf'), }, - # The monitor is missing from the Nexus files, - # so to be able to load Nexus files anyway - # this is set to None. - WavelengthMonitor[RunType]: None, } diff --git a/packages/essreflectometry/src/ess/reflectometry/conversions.py b/packages/essreflectometry/src/ess/reflectometry/conversions.py index d28663943..4a11b2ff8 100644 --- a/packages/essreflectometry/src/ess/reflectometry/conversions.py +++ b/packages/essreflectometry/src/ess/reflectometry/conversions.py @@ -4,11 +4,6 @@ from scipp.constants import pi from scippneutron._utils import elem_dtype -from .types import ( - ProtonCurrent, - RunType, -) - def reflectometry_q(wavelength: sc.Variable, theta: sc.Variable) -> sc.Variable: """ @@ -31,37 +26,4 @@ def reflectometry_q(wavelength: sc.Variable, theta: sc.Variable) -> sc.Variable: return c * sc.sin(theta.astype(dtype, copy=False)) / wavelength -def add_proton_current_coord( - da: sc.DataArray, - pc: ProtonCurrent[RunType], -) -> sc.DataArray: - """Find the proton current value for each event and - adds it as a coord to the data array.""" - pc_lookup = sc.lookup( - pc.assign_coords(time=pc.coords['time'].to(unit='ns')), - dim='time', - mode='previous', - fill_value=sc.scalar(float('nan'), unit=pc.unit), - ) - # Useful for comparing the proton current to what is typical - da = da.assign_coords(median_proton_current=sc.median(pc).data) - da.coords.set_aligned('median_proton_current', False) - da = da.bins.assign_coords( - proton_current=pc_lookup(da.bins.coords['event_time_zero']) - ) - return da - - -def add_proton_current_mask(da: sc.DataArray) -> sc.DataArray: - """Masks events where the proton current was too low or where - the proton current is nan.""" - # Take inverse and use >= because we want to mask nan values - da = da.bins.assign_masks( - proton_current_too_low=~( - da.bins.coords['proton_current'] >= da.coords['median_proton_current'] / 4 - ) - ) - return da - - providers = () diff --git a/packages/essreflectometry/src/ess/reflectometry/corrections.py b/packages/essreflectometry/src/ess/reflectometry/corrections.py index 13eccf155..58169df19 100644 --- a/packages/essreflectometry/src/ess/reflectometry/corrections.py +++ b/packages/essreflectometry/src/ess/reflectometry/corrections.py @@ -43,9 +43,25 @@ def correct_by_footprint(da: sc.DataArray) -> sc.DataArray: ) -def correct_by_proton_current(da: sc.DataArray) -> sc.DataArray: +def correct_by_proton_current( + da: sc.DataArray, proton_current: sc.DataArray +) -> sc.DataArray: "Corrects the data by the proton current during the time of data collection" - return da / da.bins.coords['proton_current'] + pc_lookup = sc.lookup( + sc.values( + proton_current.assign_coords( + time=proton_current.coords['time'].to(unit='ns') + ) + ), + dim='time', + mode='previous', + fill_value=sc.scalar(float('nan'), unit=proton_current.unit), + ) + median_pc = sc.median(sc.values(proton_current)).data + pc = pc_lookup(da.bins.coords['event_time_zero']) + # The test is "not larger than" because that masks nan values as well. + da = da.bins.assign_masks(proton_current_too_low=~(pc >= median_pc / 4)) + return da / pc def correct_sample_rotation( diff --git a/packages/essreflectometry/tests/amor/pipeline_test.py b/packages/essreflectometry/tests/amor/pipeline_test.py index 4579741a6..a30426a32 100644 --- a/packages/essreflectometry/tests/amor/pipeline_test.py +++ b/packages/essreflectometry/tests/amor/pipeline_test.py @@ -199,12 +199,10 @@ def test_proton_current(amor_pipeline: sciline.Pipeline): ) da_with_proton_current = amor_pipeline.compute(ReducibleData[SampleRun]) - assert "proton_current" in da_with_proton_current.bins.coords assert "proton_current_too_low" in da_with_proton_current.bins.masks assert da_with_proton_current.bins.masks["proton_current_too_low"].any() assert not da_with_proton_current.bins.masks["proton_current_too_low"].all() - assert "proton_current" not in da_without_proton_current.bins.coords assert "proton_current_too_low" not in da_without_proton_current.bins.masks t = ( diff --git a/packages/essreflectometry/tests/estia/mcstas_data_test.py b/packages/essreflectometry/tests/estia/mcstas_data_test.py index 24543c3d4..db6c58b66 100644 --- a/packages/essreflectometry/tests/estia/mcstas_data_test.py +++ b/packages/essreflectometry/tests/estia/mcstas_data_test.py @@ -24,6 +24,7 @@ ) from ess.estia.types import WavelengthMonitor from ess.reflectometry import orso +from ess.reflectometry.corrections import correct_by_proton_current from ess.reflectometry.types import ( BeamDivergenceLimits, CorrectionsToApply, @@ -123,6 +124,28 @@ def test_compute_reducible_data_with_monitor(estia_mcstas_pipeline: sciline.Pipe ) +def test_compute_reducible_data_with_proton_current( + estia_mcstas_pipeline: sciline.Pipeline, +): + wf = estia_mcstas_pipeline + wf[Filename[SampleRun]] = estia_mcstas_sample_run(11) + without_proton_current = wf.compute(ReducibleData[SampleRun]) + wf[ProtonCurrent[SampleRun]] = sc.DataArray( + sc.array(dims=['time'], values=[2.0], variances=[1.0]), + coords={'time': sc.datetimes(dims=['time'], values=[0], unit='s')}, + ) + corrections = wf.compute(CorrectionsToApply) + wf[CorrectionsToApply] = {*corrections, 'proton_current'} + with_proton_current = wf.compute(ReducibleData[SampleRun]) + scipp.testing.assert_allclose( + correct_by_proton_current( + without_proton_current, + proton_current=wf.compute(ProtonCurrent[SampleRun]), + ), + with_proton_current, + ) + + def test_can_compute_reflectivity_curve_exact_wavelengths( estia_mcstas_pipeline: sciline.Pipeline, ): From 8ddb6fedc18410952ec91d5bd1b93f799b2d3acf Mon Sep 17 00:00:00 2001 From: Johannes Kasimir Date: Tue, 28 Apr 2026 15:56:29 +0200 Subject: [PATCH 4/8] test: add test suite for estia nexusfile workflow --- .../tests/estia/nexus_data_test.py | 165 ++++++++++++++++++ 1 file changed, 165 insertions(+) create mode 100644 packages/essreflectometry/tests/estia/nexus_data_test.py diff --git a/packages/essreflectometry/tests/estia/nexus_data_test.py b/packages/essreflectometry/tests/estia/nexus_data_test.py new file mode 100644 index 000000000..220c2f29a --- /dev/null +++ b/packages/essreflectometry/tests/estia/nexus_data_test.py @@ -0,0 +1,165 @@ +# SPDX-License-Identifier: BSD-3-Clause +# Copyright (c) 2024 Scipp contributors (https://github.com/scipp) +from datetime import datetime +from pathlib import Path +from zoneinfo import ZoneInfo + +import numpy as np +import pytest +import sciline +import scipp as sc +import scipp.testing +from ess.reduce.normalization import normalize_by_monitor_histogram +from ess.reduce.uncertainty import UncertaintyBroadcastMode +from orsopy import fileio + +from ess.estia import EstiaWorkflow +from ess.estia.data import ( + estia_mcstas_nexus_reference_example, + estia_mcstas_nexus_sample_example, + estia_wavelength_lookup_table, +) +from ess.estia.types import WavelengthMonitor +from ess.reflectometry import orso, supermirror +from ess.reflectometry.types import ( + BeamDivergenceLimits, + CorrectionsToApply, + Filename, + LookupTableFilename, + ProtonCurrent, + QBins, + ReducibleData, + ReferenceRun, + ReflectivityOverQ, + SampleRun, + WavelengthBins, + YIndexLimits, + ZIndexLimits, +) + + +@pytest.fixture +def estia_pipeline() -> sciline.Pipeline: + wf = EstiaWorkflow() + wf[Filename[ReferenceRun]] = estia_mcstas_nexus_reference_example() + + wf[YIndexLimits] = sc.scalar(35), sc.scalar(64) + wf[ZIndexLimits] = sc.scalar(0), sc.scalar(48 * 32) + wf[BeamDivergenceLimits] = sc.scalar(-1.0, unit='deg'), sc.scalar(1.0, unit='deg') + wf[WavelengthBins] = sc.geomspace('wavelength', 3.5, 12, 2001, unit='angstrom') + wf[QBins] = sc.geomspace('Q', 0.005, 0.1, 200, unit='1/angstrom') + + wf[supermirror.CriticalEdge] = sc.scalar(float('inf'), unit='1/angstrom') + wf[supermirror.Alpha] = sc.scalar(0.25 / 0.088, unit=sc.units.angstrom) + wf[supermirror.MValue] = sc.scalar(5, unit=sc.units.dimensionless) + + wf[LookupTableFilename] = estia_wavelength_lookup_table() + + wf[ProtonCurrent[SampleRun]] = sc.DataArray( + sc.array(dims=('time',), values=[]), + coords={'time': sc.array(dims=('time',), values=[], unit='s')}, + ) + wf[ProtonCurrent[ReferenceRun]] = sc.DataArray( + sc.array(dims=('time',), values=[]), + coords={'time': sc.array(dims=('time',), values=[], unit='s')}, + ) + wf[orso.OrsoCreator] = orso.OrsoCreator( + fileio.base.Person( + name="Max Mustermann", + affiliation="European Spallation Source ERIC", + contact="max.mustermann@ess.eu", + ) + ) + wf[orso.OrsoExperiment] = orso.OrsoExperiment( + fileio.data_source.Experiment( + title='McStas run', + instrument='Estia', + facility='ESS', + start_date=datetime(2025, 3, 20, tzinfo=ZoneInfo("Europe/Stockholm")), + probe='neutron', + ) + ) + wf[orso.OrsoOwner] = orso.OrsoOwner( + fileio.base.Person( + name='John Doe', + contact='john.doe@ess.eu', + affiliation='ESS', + ) + ) + wf[orso.OrsoSample] = orso.OrsoSample(fileio.data_source.Sample.empty()) + wf[WavelengthMonitor[SampleRun]] = None + wf[WavelengthMonitor[ReferenceRun]] = None + return wf + + +def test_compute_reducible_data(estia_pipeline: sciline.Pipeline): + estia_pipeline[Filename[SampleRun]] = estia_mcstas_nexus_sample_example( + 'Ni/Ti-multilayer' + )[0] + da = estia_pipeline.compute(ReducibleData[SampleRun]) + assert da.dims == ('strip', 'blade', 'wire') + assert da.shape == (64, 48, 32) + assert 'position' in da.coords + assert 'sample_rotation' in da.coords + assert 'detector_rotation' in da.coords + assert 'theta' in da.coords + assert 'wavelength' in da.bins.coords + assert 'Q' in da.bins.coords + + +def test_compute_reducible_data_with_monitor(estia_pipeline: sciline.Pipeline): + wf = estia_pipeline + wf[Filename[SampleRun]] = estia_mcstas_nexus_sample_example('Ni/Ti-multilayer')[0] + without_monitor = wf.compute(ReducibleData[SampleRun]) + wf[WavelengthMonitor[SampleRun]] = sc.DataArray( + sc.array(dims=['wavelength'], values=[30.0], variances=[1.0]), + coords={'wavelength': sc.linspace('wavelength', 0, 15, 2, unit='angstrom')}, + ) + corrections = wf.compute(CorrectionsToApply) + wf[CorrectionsToApply] = {*corrections, 'monitor'} + with_monitor = wf.compute(ReducibleData[SampleRun]) + scipp.testing.assert_allclose( + normalize_by_monitor_histogram( + without_monitor, + monitor=wf.compute(WavelengthMonitor[SampleRun]), + uncertainty_broadcast_mode=UncertaintyBroadcastMode.drop, + ), + with_monitor, + ) + + +def test_can_compute_reflectivity_curve(estia_pipeline: sciline.Pipeline): + estia_pipeline[Filename[SampleRun]] = estia_mcstas_nexus_sample_example( + 'Ni/Ti-multilayer' + )[0] + r = estia_pipeline.compute(ReflectivityOverQ) + assert "Q" in r.coords + assert "Q_resolution" in r.coords + + +def test_orso_pipeline(estia_pipeline: sciline.Pipeline): + estia_pipeline[Filename[SampleRun]] = estia_mcstas_nexus_sample_example( + 'Ni/Ti-multilayer' + )[0] + res = estia_pipeline.compute(orso.OrsoIofQDataset) + assert res.info.data_source.experiment.instrument == "Estia" + assert res.info.reduction.software.name == "ess.reflectometry" + assert res.info.reduction.corrections == [ + "chopper ToF correction", + "footprint correction", + "supermirror calibration", + ] + assert res.data.ndim == 2 + assert res.data.shape[1] == 4 + assert np.all(res.data[:, 1] >= 0) + assert np.isfinite(res.data).all() + + +def test_save_reduced_orso_file(estia_pipeline: sciline.Pipeline, output_folder: Path): + estia_pipeline[Filename[SampleRun]] = estia_mcstas_nexus_sample_example( + 'Ni/Ti-multilayer' + )[0] + res = estia_pipeline.compute(orso.OrsoIofQDataset) + fileio.orso.save_orso( + datasets=[res], fname=output_folder / 'estia_reduced_iofq.ort' + ) From 3b4ae3203f9fb5144965127ed5290272c1895391 Mon Sep 17 00:00:00 2001 From: Johannes Kasimir Date: Wed, 29 Apr 2026 14:26:05 +0200 Subject: [PATCH 5/8] fix --- packages/essreflectometry/src/ess/estia/corrections.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/packages/essreflectometry/src/ess/estia/corrections.py b/packages/essreflectometry/src/ess/estia/corrections.py index 60ea043bf..cfa8e0978 100644 --- a/packages/essreflectometry/src/ess/estia/corrections.py +++ b/packages/essreflectometry/src/ess/estia/corrections.py @@ -26,6 +26,7 @@ def normalize_by_monitor_histogram( detector: ReducibleData[RunType], *, monitor: WavelengthMonitor[RunType], + uncertainty_broadcast_mode: UncertaintyBroadcastMode, ) -> ReducibleData[RunType]: """Normalize detector data by a histogrammed monitor. @@ -57,7 +58,7 @@ def normalize_by_monitor_histogram( return ess.reduce.normalization.normalize_by_monitor_histogram( detector=detector, monitor=monitor, - uncertainty_broadcast_mode=UncertaintyBroadcastMode.drop, + uncertainty_broadcast_mode=uncertainty_broadcast_mode, skip_range_check=False, ) @@ -103,6 +104,6 @@ def assume_time_series_constant_with_zero_default_value_if_empty(da: sc.DataArra return da.mean() if len(da) > 0 else sc.scalar(0.0, unit=da.unit) -default_corrections = {correct_by_footprint, 'proton_current', 'monitor'} +default_corrections = {correct_by_footprint, 'proton_current'} providers = (add_coords_masks_and_apply_corrections,) From b0aa006d0e92e0b7f4f2c021477b951776052961 Mon Sep 17 00:00:00 2001 From: Johannes Kasimir Date: Thu, 30 Apr 2026 09:36:22 +0200 Subject: [PATCH 6/8] fix: add uncertainty broadcast mode --- packages/essreflectometry/src/ess/estia/corrections.py | 7 ++++++- packages/essreflectometry/src/ess/estia/mcstas.py | 7 ++++++- packages/essreflectometry/src/ess/estia/workflow.py | 3 +++ 3 files changed, 15 insertions(+), 2 deletions(-) diff --git a/packages/essreflectometry/src/ess/estia/corrections.py b/packages/essreflectometry/src/ess/estia/corrections.py index cfa8e0978..d457ba1f1 100644 --- a/packages/essreflectometry/src/ess/estia/corrections.py +++ b/packages/essreflectometry/src/ess/estia/corrections.py @@ -71,6 +71,7 @@ def add_coords_masks_and_apply_corrections( wbins: WavelengthBins, proton_current: ProtonCurrent[RunType], monitor: WavelengthMonitor[RunType], + uncertainty_broadcast_mode: UncertaintyBroadcastMode, graph: CoordTransformationGraph[RunType], corrections_to_apply: CorrectionsToApply, ) -> ReducibleData[RunType]: @@ -83,7 +84,11 @@ def add_coords_masks_and_apply_corrections( for correction in corrections_to_apply: if correction == 'monitor': - da = normalize_by_monitor_histogram(da, monitor=monitor) + da = normalize_by_monitor_histogram( + da, + monitor=monitor, + uncertainty_broadcast_mode=uncertainty_broadcast_mode, + ) elif correction == 'proton_current': if len(proton_current) != 0: da = correct_by_proton_current(da, proton_current=proton_current) diff --git a/packages/essreflectometry/src/ess/estia/mcstas.py b/packages/essreflectometry/src/ess/estia/mcstas.py index b4bcdbe01..f8465ccc4 100644 --- a/packages/essreflectometry/src/ess/estia/mcstas.py +++ b/packages/essreflectometry/src/ess/estia/mcstas.py @@ -6,6 +6,7 @@ import pandas as pd import scipp as sc from ess.reduce.nexus.types import Position +from ess.reduce.uncertainty import UncertaintyBroadcastMode from scippnexus import NXsample, NXsource from ..reflectometry.load import load_h5 @@ -28,6 +29,7 @@ ) from .beamline import DETECTOR_BANK_SIZES from .corrections import add_coords_masks_and_apply_corrections +from .types import WavelengthMonitor def parse_metadata_ascii(lines) -> dict: @@ -294,6 +296,8 @@ def use_mcstas_wavelengths_instead_of_estimates_from_time_of_arrival( bdlim: BeamDivergenceLimits, wbins: WavelengthBins, proton_current: ProtonCurrent[RunType], + monitor: WavelengthMonitor[RunType], + uncertainty_broadcast_mode: UncertaintyBroadcastMode, graph: CoordTransformationGraph[RunType], corrections_to_apply: CorrectionsToApply, ) -> ReducibleData[RunType]: @@ -303,7 +307,8 @@ def use_mcstas_wavelengths_instead_of_estimates_from_time_of_arrival( zlims=zlims, bdlim=bdlim, wbins=wbins, - monitor=None, + monitor=monitor, + uncertainty_broadcast_mode=uncertainty_broadcast_mode, proton_current=proton_current, graph={ **graph, diff --git a/packages/essreflectometry/src/ess/estia/workflow.py b/packages/essreflectometry/src/ess/estia/workflow.py index 796cb1112..cf144ac3d 100644 --- a/packages/essreflectometry/src/ess/estia/workflow.py +++ b/packages/essreflectometry/src/ess/estia/workflow.py @@ -5,6 +5,7 @@ import scipp as sc import scippnexus as snx from ess.reduce.nexus.types import TransformationTimeFilter +from ess.reduce.uncertainty import UncertaintyBroadcastMode from ..reflectometry import providers as reflectometry_providers from ..reflectometry import supermirror @@ -76,6 +77,7 @@ def mcstas_default_parameters() -> dict: "multiblade_detector": 0.06, }, WavelengthMonitor[RunType]: None, + UncertaintyBroadcastMode: UncertaintyBroadcastMode.drop, } @@ -90,6 +92,7 @@ def default_parameters() -> dict: LookupTableRelativeErrorThreshold: { "multiblade_detector": float('inf'), }, + UncertaintyBroadcastMode: UncertaintyBroadcastMode.drop, } From 0e9fbd50cd1f1ff36bbad5ff44ed0a9132441161 Mon Sep 17 00:00:00 2001 From: Johannes Kasimir Date: Fri, 1 May 2026 11:47:18 +0200 Subject: [PATCH 7/8] fix: now that corrections are configurable, remove guard for empty proton charge --- packages/essreflectometry/src/ess/estia/corrections.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/packages/essreflectometry/src/ess/estia/corrections.py b/packages/essreflectometry/src/ess/estia/corrections.py index d457ba1f1..ebaf4b47c 100644 --- a/packages/essreflectometry/src/ess/estia/corrections.py +++ b/packages/essreflectometry/src/ess/estia/corrections.py @@ -90,8 +90,7 @@ def add_coords_masks_and_apply_corrections( uncertainty_broadcast_mode=uncertainty_broadcast_mode, ) elif correction == 'proton_current': - if len(proton_current) != 0: - da = correct_by_proton_current(da, proton_current=proton_current) + da = correct_by_proton_current(da, proton_current=proton_current) else: da = correction(da) From e172d004df1b02a83fe511a16104198cd87f4886 Mon Sep 17 00:00:00 2001 From: Johannes Kasimir Date: Tue, 5 May 2026 11:09:17 +0200 Subject: [PATCH 8/8] feat: add monitor type for freia (WIP) --- .../src/ess/freia/corrections.py | 65 ++++++++++++++++--- .../essreflectometry/src/ess/freia/types.py | 8 ++- 2 files changed, 63 insertions(+), 10 deletions(-) diff --git a/packages/essreflectometry/src/ess/freia/corrections.py b/packages/essreflectometry/src/ess/freia/corrections.py index a00487372..51ed8222c 100644 --- a/packages/essreflectometry/src/ess/freia/corrections.py +++ b/packages/essreflectometry/src/ess/freia/corrections.py @@ -1,11 +1,9 @@ # SPDX-License-Identifier: BSD-3-Clause # Copyright (c) 2025 Scipp contributors (https://github.com/scipp) +import ess.reduce import scipp as sc +from ess.reduce.uncertainty import UncertaintyBroadcastMode -from ..reflectometry.conversions import ( - add_proton_current_coord, - add_proton_current_mask, -) from ..reflectometry.corrections import correct_by_proton_current from ..reflectometry.types import ( BeamDivergenceLimits, @@ -21,6 +19,48 @@ ) from .conversions import add_coords from .maskings import add_masks +from .types import WavelengthMonitor + + +def normalize_by_monitor_histogram( + detector: ReducibleData[RunType], + *, + monitor: WavelengthMonitor[RunType], + uncertainty_broadcast_mode: UncertaintyBroadcastMode, +) -> ReducibleData[RunType]: + """Normalize detector data by a histogrammed monitor. + + The detector is normalized according to + + .. math:: + + d_i^\\text{Norm} = \\frac{d_i}{m_i} \\Delta \\lambda_i + + Parameters + ---------- + detector: + Input event data in wavelength. + monitor: + A histogrammed monitor in wavelength. + uncertainty_broadcast_mode: + Choose how uncertainties of the monitor are broadcast to the sample data. + + Returns + ------- + : + `detector` normalized by a monitor. + + See also + -------- + ess.reduce.normalization.normalize_by_monitor_histogram: + For details and the actual implementation. + """ + return ess.reduce.normalization.normalize_by_monitor_histogram( + detector=detector, + monitor=monitor, + uncertainty_broadcast_mode=uncertainty_broadcast_mode, + skip_range_check=False, + ) def add_coords_masks_and_apply_corrections( @@ -30,6 +70,8 @@ def add_coords_masks_and_apply_corrections( bdlim: BeamDivergenceLimits, wbins: WavelengthBins, proton_current: ProtonCurrent[RunType], + monitor: WavelengthMonitor[RunType], + uncertainty_broadcast_mode: UncertaintyBroadcastMode, graph: CoordTransformationGraph[RunType], corrections_to_apply: CorrectionsToApply, ) -> ReducibleData[RunType]: @@ -40,12 +82,17 @@ def add_coords_masks_and_apply_corrections( da = add_coords(da, graph) da = add_masks(da, ylim, zlims, bdlim, wbins) - if len(proton_current) != 0: - da = add_proton_current_coord(da, proton_current) - da = add_proton_current_mask(da) - for correction in corrections_to_apply: - da = correction(da) + if correction == 'monitor': + da = normalize_by_monitor_histogram( + da, + monitor=monitor, + uncertainty_broadcast_mode=uncertainty_broadcast_mode, + ) + elif correction == 'proton_current': + da = correct_by_proton_current(da, proton_current=proton_current) + else: + da = correction(da) return ReducibleData[RunType](da) diff --git a/packages/essreflectometry/src/ess/freia/types.py b/packages/essreflectometry/src/ess/freia/types.py index 8acdf9fcd..472dac2d4 100644 --- a/packages/essreflectometry/src/ess/freia/types.py +++ b/packages/essreflectometry/src/ess/freia/types.py @@ -1,8 +1,14 @@ # Copyright (c) 2025 Scipp contributors (https://github.com/scipp) -from typing import NewType +from typing import NewType, TypeAlias import scipp as sc +from ess.reduce.nexus.types import IncidentMonitor +from ess.reduce.unwrap.types import WavelengthMonitor + +from ..reflectometry.types import RunType WavelengthResolution = NewType("WavelengthResolution", sc.Variable) AngularResolution = NewType("AngularResolution", sc.Variable) SampleSizeResolution = NewType("SampleSizeResolution", sc.Variable) + +WavelengthMonitor: TypeAlias = WavelengthMonitor[RunType, IncidentMonitor]