From 826279bcf99724dd784f0cae4d077f92fe9a1030 Mon Sep 17 00:00:00 2001 From: Johannes Kasimir Date: Mon, 27 Apr 2026 15:40:41 +0200 Subject: [PATCH 01/11] 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 02/11] 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 03/11] 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 04/11] 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 05/11] 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 06/11] 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 07/11] 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 08/11] 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] From 56769616cd698b197e4795ca5ddf92fbadfcbb86 Mon Sep 17 00:00:00 2001 From: Johannes Kasimir Date: Tue, 19 May 2026 16:16:06 +0200 Subject: [PATCH 09/11] fix --- packages/essreflectometry/src/ess/amor/load.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/packages/essreflectometry/src/ess/amor/load.py b/packages/essreflectometry/src/ess/amor/load.py index 75ec37ae2..48f0973dd 100644 --- a/packages/essreflectometry/src/ess/amor/load.py +++ b/packages/essreflectometry/src/ess/amor/load.py @@ -117,7 +117,7 @@ def load_amor_proton_charge( fp: Filename[RunType], ) -> ProtonCharge[RunType]: """Load proton charge log from the NeXus detector group.""" - (pc,) = load_nx(fp, 'NXentry/NXinstrument/NXdetector/proton_charge') + (pc,) = load_nx(fp, 'NXentry/NXinstrument/NXdetector/proton_current') pc = pc['value']['dim_1', 0] pc.data.unit = 'mA/s' return pc From dd0c8304511b0db47fda15651eb64fe335ead37c Mon Sep 17 00:00:00 2001 From: Johannes Kasimir Date: Wed, 20 May 2026 10:20:11 +0200 Subject: [PATCH 10/11] fix --- packages/essreflectometry/src/ess/amor/load.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/packages/essreflectometry/src/ess/amor/load.py b/packages/essreflectometry/src/ess/amor/load.py index 48f0973dd..49647d5de 100644 --- a/packages/essreflectometry/src/ess/amor/load.py +++ b/packages/essreflectometry/src/ess/amor/load.py @@ -116,7 +116,7 @@ def load_amor_detector_rotation(fp: Filename[RunType]) -> DetectorRotation[RunTy def load_amor_proton_charge( fp: Filename[RunType], ) -> ProtonCharge[RunType]: - """Load proton charge log from the NeXus detector group.""" + """Load proton current log from the NeXus detector group.""" (pc,) = load_nx(fp, 'NXentry/NXinstrument/NXdetector/proton_current') pc = pc['value']['dim_1', 0] pc.data.unit = 'mA/s' From d01f604259172bf04d4a387a2505ff95f63562d4 Mon Sep 17 00:00:00 2001 From: Johannes Kasimir Date: Wed, 20 May 2026 14:15:46 +0200 Subject: [PATCH 11/11] feat: constructors for workflows with different monitor normalizations --- .../src/ess/estia/__init__.py | 21 ++++- .../src/ess/estia/corrections.py | 58 +++++++----- .../essreflectometry/src/ess/estia/mcstas.py | 15 +-- .../essreflectometry/src/ess/estia/types.py | 5 +- .../src/ess/estia/workflow.py | 93 +++++++++++++++--- .../src/ess/freia/__init__.py | 21 ++++- .../src/ess/freia/corrections.py | 58 +++++++----- .../src/ess/freia/workflow.py | 91 ++++++++++++++++-- .../src/ess/reflectometry/corrections.py | 94 ++++++++++++++++++- .../src/ess/reflectometry/types.py | 4 + .../tests/estia/mcstas_data_test.py | 56 +++++++---- .../tests/estia/nexus_data_test.py | 47 +++++++--- .../tests/estia/workflow_test.py | 41 ++++++++ .../tests/freia/workflow_test.py | 29 ++++++ 14 files changed, 523 insertions(+), 110 deletions(-) create mode 100644 packages/essreflectometry/tests/estia/workflow_test.py diff --git a/packages/essreflectometry/src/ess/estia/__init__.py b/packages/essreflectometry/src/ess/estia/__init__.py index 742f60f43..ba503d461 100644 --- a/packages/essreflectometry/src/ess/estia/__init__.py +++ b/packages/essreflectometry/src/ess/estia/__init__.py @@ -9,7 +9,18 @@ SampleSizeResolution, WavelengthResolution, ) -from .workflow import EstiaMcStasWorkflow, EstiaWorkflow +from .workflow import ( + EstiaMcStasMonitorHistogramWorkflow, + EstiaMcStasMonitorIntegratedWorkflow, + EstiaMcStasProtonChargeWorkflow, + EstiaMcStasUnnormalizedWorkflow, + EstiaMcStasWorkflow, + EstiaMonitorHistogramWorkflow, + EstiaMonitorIntegratedWorkflow, + EstiaProtonChargeWorkflow, + EstiaUnnormalizedWorkflow, + EstiaWorkflow, +) try: __version__ = importlib.metadata.version(__package__ or __name__) @@ -18,7 +29,15 @@ __all__ = [ "AngularResolution", + "EstiaMcStasMonitorHistogramWorkflow", + "EstiaMcStasMonitorIntegratedWorkflow", + "EstiaMcStasProtonChargeWorkflow", + "EstiaMcStasUnnormalizedWorkflow", "EstiaMcStasWorkflow", + "EstiaMonitorHistogramWorkflow", + "EstiaMonitorIntegratedWorkflow", + "EstiaProtonChargeWorkflow", + "EstiaUnnormalizedWorkflow", "EstiaWorkflow", "SampleSizeResolution", "WavelengthResolution", diff --git a/packages/essreflectometry/src/ess/estia/corrections.py b/packages/essreflectometry/src/ess/estia/corrections.py index fbddf3304..ff2df30ca 100644 --- a/packages/essreflectometry/src/ess/estia/corrections.py +++ b/packages/essreflectometry/src/ess/estia/corrections.py @@ -1,17 +1,18 @@ # SPDX-License-Identifier: BSD-3-Clause # Copyright (c) 2025 Scipp contributors (https://github.com/scipp) -import ess.reduce +import sciline import scipp as sc from ess.reduce.uncertainty import UncertaintyBroadcastMode -from ..reflectometry.corrections import correct_by_proton_charge +from ..reflectometry import corrections as common_corrections +from ..reflectometry.corrections import RunNormalization from ..reflectometry.types import ( BeamDivergenceLimits, CoordTransformationGraph, CorrectionsToApply, - ProtonCharge, ReducibleData, RunType, + RunUnnormalizedData, WavelengthBins, WavelengthDetector, YIndexLimits, @@ -23,7 +24,7 @@ def normalize_by_monitor_histogram( - detector: ReducibleData[RunType], + detector: RunUnnormalizedData[RunType], *, monitor: WavelengthMonitor[RunType], uncertainty_broadcast_mode: UncertaintyBroadcastMode, @@ -55,11 +56,36 @@ def normalize_by_monitor_histogram( ess.reduce.normalization.normalize_by_monitor_histogram: For details and the actual implementation. """ - return ess.reduce.normalization.normalize_by_monitor_histogram( + return common_corrections.normalize_by_monitor_histogram( detector=detector, monitor=monitor, uncertainty_broadcast_mode=uncertainty_broadcast_mode, - skip_range_check=False, + ) + + +def normalize_by_monitor_integrated( + detector: RunUnnormalizedData[RunType], + *, + monitor: WavelengthMonitor[RunType], + uncertainty_broadcast_mode: UncertaintyBroadcastMode, +) -> ReducibleData[RunType]: + """Normalize detector data by an integrated wavelength monitor.""" + return common_corrections.normalize_by_monitor_integrated( + detector=detector, + monitor=monitor, + uncertainty_broadcast_mode=uncertainty_broadcast_mode, + ) + + +def insert_run_normalization( + workflow: sciline.Pipeline, run_norm: RunNormalization +) -> None: + """Insert providers for a specific ESTIA run normalization into a workflow.""" + common_corrections.insert_run_normalization( + workflow, + run_norm, + monitor_histogram_provider=normalize_by_monitor_histogram, + monitor_integrated_provider=normalize_by_monitor_integrated, ) @@ -69,12 +95,9 @@ def add_coords_masks_and_apply_corrections( zlims: ZIndexLimits, bdlim: BeamDivergenceLimits, wbins: WavelengthBins, - proton_charge: ProtonCharge[RunType], - monitor: WavelengthMonitor[RunType], - uncertainty_broadcast_mode: UncertaintyBroadcastMode, graph: CoordTransformationGraph[RunType], corrections_to_apply: CorrectionsToApply, -) -> ReducibleData[RunType]: +) -> RunUnnormalizedData[RunType]: """ Computes coordinates, masks and corrections that are the same for the sample measurement and the reference measurement. @@ -83,18 +106,9 @@ def add_coords_masks_and_apply_corrections( da = add_masks(da, ylim, zlims, bdlim, wbins) for correction in corrections_to_apply: - if correction == 'monitor': - da = normalize_by_monitor_histogram( - da, - monitor=monitor, - uncertainty_broadcast_mode=uncertainty_broadcast_mode, - ) - elif correction == 'proton_charge': - da = correct_by_proton_charge(da, proton_charge=proton_charge) - else: - da = correction(da) + da = correction(da) - return ReducibleData[RunType](da) + return RunUnnormalizedData[RunType](da) def correct_by_footprint(da: sc.DataArray) -> sc.DataArray: @@ -108,6 +122,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_charge'} +default_corrections = {correct_by_footprint} 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 2b128f5cd..061a74a14 100644 --- a/packages/essreflectometry/src/ess/estia/mcstas.py +++ b/packages/essreflectometry/src/ess/estia/mcstas.py @@ -6,7 +6,6 @@ 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 @@ -17,10 +16,9 @@ DetectorLtotal, DetectorRotation, Filename, - ProtonCharge, RawDetector, - ReducibleData, RunType, + RunUnnormalizedData, SampleRotation, SampleRotationOffset, WavelengthBins, @@ -29,7 +27,6 @@ ) 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: @@ -295,28 +292,22 @@ def use_mcstas_wavelengths_instead_of_estimates_from_time_of_arrival( zlims: ZIndexLimits, bdlim: BeamDivergenceLimits, wbins: WavelengthBins, - proton_charge: ProtonCharge[RunType], - monitor: WavelengthMonitor[RunType], - uncertainty_broadcast_mode: UncertaintyBroadcastMode, graph: CoordTransformationGraph[RunType], corrections_to_apply: CorrectionsToApply, -) -> ReducibleData[RunType]: +) -> RunUnnormalizedData[RunType]: out = add_coords_masks_and_apply_corrections( da=da, ylim=ylim, zlims=zlims, bdlim=bdlim, wbins=wbins, - monitor=monitor, - uncertainty_broadcast_mode=uncertainty_broadcast_mode, - proton_charge=proton_charge, graph={ **graph, "wavelength": lambda wavelength_from_mcstas: wavelength_from_mcstas, }, corrections_to_apply=corrections_to_apply, ) - return ReducibleData[RunType](out) + return RunUnnormalizedData[RunType](out) providers = ( diff --git a/packages/essreflectometry/src/ess/estia/types.py b/packages/essreflectometry/src/ess/estia/types.py index efcfb64fd..9c2711d77 100644 --- a/packages/essreflectometry/src/ess/estia/types.py +++ b/packages/essreflectometry/src/ess/estia/types.py @@ -3,7 +3,7 @@ import scipp as sc from ess.reduce.nexus.types import CaveMonitor -from ess.reduce.unwrap.types import WavelengthMonitor +from ess.reduce.unwrap.types import WavelengthMonitor as _WavelengthMonitor from ess.reflectometry.types import RunType @@ -11,4 +11,5 @@ AngularResolution = NewType("AngularResolution", sc.Variable) SampleSizeResolution = NewType("SampleSizeResolution", sc.Variable) -WavelengthMonitor: TypeAlias = WavelengthMonitor[RunType, CaveMonitor] + +WavelengthMonitor: TypeAlias = _WavelengthMonitor[RunType, CaveMonitor] diff --git a/packages/essreflectometry/src/ess/estia/workflow.py b/packages/essreflectometry/src/ess/estia/workflow.py index c1c103cd3..04ad01ee7 100644 --- a/packages/essreflectometry/src/ess/estia/workflow.py +++ b/packages/essreflectometry/src/ess/estia/workflow.py @@ -6,6 +6,7 @@ import scippnexus as snx from ess.reduce.nexus.types import TransformationTimeFilter from ess.reduce.uncertainty import UncertaintyBroadcastMode +from ess.reduce.workflow import register_workflow from ..reflectometry import providers as reflectometry_providers from ..reflectometry import supermirror @@ -28,7 +29,7 @@ normalization, orso, ) -from .types import WavelengthMonitor +from .corrections import RunNormalization, insert_run_normalization _general_providers = ( *reflectometry_providers, @@ -70,13 +71,10 @@ def mcstas_default_parameters() -> dict: sc.scalar(0.75, unit='deg'), ), SampleRotationOffset[RunType]: sc.scalar(0.0, unit='deg'), - CorrectionsToApply: ( - corrections.default_corrections - {'monitor', 'proton_charge'} - ), + CorrectionsToApply: corrections.default_corrections, LookupTableRelativeErrorThreshold: { "multiblade_detector": 0.06, }, - WavelengthMonitor[RunType]: None, UncertaintyBroadcastMode: UncertaintyBroadcastMode.drop, } @@ -86,8 +84,7 @@ def default_parameters() -> dict: return { NeXusDetectorName: "multiblade_detector", SampleRotationOffset[RunType]: sc.scalar(0.0, unit='deg'), - CorrectionsToApply: corrections.default_corrections - - {'monitor', 'proton_charge'}, + CorrectionsToApply: corrections.default_corrections, DetectorSpatialResolution: 0.0025 * sc.units.m, LookupTableRelativeErrorThreshold: { "multiblade_detector": float('inf'), @@ -96,21 +93,31 @@ def default_parameters() -> dict: } -def EstiaMcStasWorkflow() -> sciline.Pipeline: +def EstiaMcStasWorkflow( + *, + run_norm: RunNormalization = RunNormalization.none, + **kwargs, +) -> sciline.Pipeline: """Workflow for reduction of McStas data for the Estia instrument.""" - workflow = beamline.LoadNeXusWorkflow() + workflow = beamline.LoadNeXusWorkflow(**kwargs) for provider in mcstas_providers: workflow.insert(provider) + insert_run_normalization(workflow, run_norm) for name, param in mcstas_default_parameters().items(): workflow[name] = param return workflow -def EstiaWorkflow() -> sciline.Pipeline: +def EstiaWorkflow( + *, + run_norm: RunNormalization = RunNormalization.proton_charge, + **kwargs, +) -> sciline.Pipeline: """Workflow for reduction of data for the Estia instrument.""" - workflow = beamline.LoadNeXusWorkflow() + workflow = beamline.LoadNeXusWorkflow(**kwargs) for provider in providers: workflow.insert(provider) + insert_run_normalization(workflow, run_norm) for name, param in default_parameters().items(): workflow[name] = param @@ -122,3 +129,67 @@ def EstiaWorkflow() -> sciline.Pipeline: corrections.assume_time_series_constant_with_zero_default_value_if_empty ) return workflow + + +@register_workflow +def EstiaMcStasUnnormalizedWorkflow() -> sciline.Pipeline: + """Workflow for Estia McStas data without run normalization.""" + return EstiaMcStasWorkflow(run_norm=RunNormalization.none) + + +@register_workflow +def EstiaMcStasMonitorHistogramWorkflow() -> sciline.Pipeline: + """Workflow for Estia McStas data using histogrammed monitor normalization.""" + return EstiaMcStasWorkflow(run_norm=RunNormalization.monitor_histogram) + + +@register_workflow +def EstiaMcStasMonitorIntegratedWorkflow() -> sciline.Pipeline: + """Workflow for Estia McStas data using integrated monitor normalization.""" + return EstiaMcStasWorkflow(run_norm=RunNormalization.monitor_integrated) + + +@register_workflow +def EstiaMcStasProtonChargeWorkflow() -> sciline.Pipeline: + """Workflow for Estia McStas data using proton charge normalization.""" + return EstiaMcStasWorkflow(run_norm=RunNormalization.proton_charge) + + +@register_workflow +def EstiaUnnormalizedWorkflow() -> sciline.Pipeline: + """Workflow for Estia NeXus data without run normalization.""" + return EstiaWorkflow(run_norm=RunNormalization.none) + + +@register_workflow +def EstiaMonitorHistogramWorkflow() -> sciline.Pipeline: + """Workflow for Estia NeXus data using histogrammed monitor normalization.""" + return EstiaWorkflow(run_norm=RunNormalization.monitor_histogram) + + +@register_workflow +def EstiaMonitorIntegratedWorkflow() -> sciline.Pipeline: + """Workflow for Estia NeXus data using integrated monitor normalization.""" + return EstiaWorkflow(run_norm=RunNormalization.monitor_integrated) + + +@register_workflow +def EstiaProtonChargeWorkflow() -> sciline.Pipeline: + """Workflow for Estia NeXus data using proton charge normalization.""" + return EstiaWorkflow(run_norm=RunNormalization.proton_charge) + + +__all__ = [ + 'EstiaMcStasMonitorHistogramWorkflow', + 'EstiaMcStasMonitorIntegratedWorkflow', + 'EstiaMcStasProtonChargeWorkflow', + 'EstiaMcStasUnnormalizedWorkflow', + 'EstiaMcStasWorkflow', + 'EstiaMonitorHistogramWorkflow', + 'EstiaMonitorIntegratedWorkflow', + 'EstiaProtonChargeWorkflow', + 'EstiaUnnormalizedWorkflow', + 'EstiaWorkflow', + 'default_parameters', + 'mcstas_default_parameters', +] diff --git a/packages/essreflectometry/src/ess/freia/__init__.py b/packages/essreflectometry/src/ess/freia/__init__.py index afac536da..3860491e8 100644 --- a/packages/essreflectometry/src/ess/freia/__init__.py +++ b/packages/essreflectometry/src/ess/freia/__init__.py @@ -9,7 +9,18 @@ SampleSizeResolution, WavelengthResolution, ) -from .workflow import FreiaMcStasWorkflow, FreiaWorkflow +from .workflow import ( + FreiaMcStasMonitorHistogramWorkflow, + FreiaMcStasMonitorIntegratedWorkflow, + FreiaMcStasProtonChargeWorkflow, + FreiaMcStasUnnormalizedWorkflow, + FreiaMcStasWorkflow, + FreiaMonitorHistogramWorkflow, + FreiaMonitorIntegratedWorkflow, + FreiaProtonChargeWorkflow, + FreiaUnnormalizedWorkflow, + FreiaWorkflow, +) try: __version__ = importlib.metadata.version(__package__ or __name__) @@ -18,7 +29,15 @@ __all__ = [ "AngularResolution", + "FreiaMcStasMonitorHistogramWorkflow", + "FreiaMcStasMonitorIntegratedWorkflow", + "FreiaMcStasProtonChargeWorkflow", + "FreiaMcStasUnnormalizedWorkflow", "FreiaMcStasWorkflow", + "FreiaMonitorHistogramWorkflow", + "FreiaMonitorIntegratedWorkflow", + "FreiaProtonChargeWorkflow", + "FreiaUnnormalizedWorkflow", "FreiaWorkflow", "SampleSizeResolution", "WavelengthResolution", diff --git a/packages/essreflectometry/src/ess/freia/corrections.py b/packages/essreflectometry/src/ess/freia/corrections.py index 113bc0532..a6f5c1238 100644 --- a/packages/essreflectometry/src/ess/freia/corrections.py +++ b/packages/essreflectometry/src/ess/freia/corrections.py @@ -1,17 +1,18 @@ # SPDX-License-Identifier: BSD-3-Clause # Copyright (c) 2025 Scipp contributors (https://github.com/scipp) -import ess.reduce +import sciline import scipp as sc from ess.reduce.uncertainty import UncertaintyBroadcastMode -from ..reflectometry.corrections import correct_by_proton_charge +from ..reflectometry import corrections as common_corrections +from ..reflectometry.corrections import RunNormalization from ..reflectometry.types import ( BeamDivergenceLimits, CoordTransformationGraph, CorrectionsToApply, - ProtonCharge, ReducibleData, RunType, + RunUnnormalizedData, WavelengthBins, WavelengthDetector, YIndexLimits, @@ -23,7 +24,7 @@ def normalize_by_monitor_histogram( - detector: ReducibleData[RunType], + detector: RunUnnormalizedData[RunType], *, monitor: WavelengthMonitor[RunType], uncertainty_broadcast_mode: UncertaintyBroadcastMode, @@ -55,11 +56,36 @@ def normalize_by_monitor_histogram( ess.reduce.normalization.normalize_by_monitor_histogram: For details and the actual implementation. """ - return ess.reduce.normalization.normalize_by_monitor_histogram( + return common_corrections.normalize_by_monitor_histogram( detector=detector, monitor=monitor, uncertainty_broadcast_mode=uncertainty_broadcast_mode, - skip_range_check=False, + ) + + +def normalize_by_monitor_integrated( + detector: RunUnnormalizedData[RunType], + *, + monitor: WavelengthMonitor[RunType], + uncertainty_broadcast_mode: UncertaintyBroadcastMode, +) -> ReducibleData[RunType]: + """Normalize detector data by an integrated wavelength monitor.""" + return common_corrections.normalize_by_monitor_integrated( + detector=detector, + monitor=monitor, + uncertainty_broadcast_mode=uncertainty_broadcast_mode, + ) + + +def insert_run_normalization( + workflow: sciline.Pipeline, run_norm: RunNormalization +) -> None: + """Insert providers for a specific FREIA run normalization into a workflow.""" + common_corrections.insert_run_normalization( + workflow, + run_norm, + monitor_histogram_provider=normalize_by_monitor_histogram, + monitor_integrated_provider=normalize_by_monitor_integrated, ) @@ -69,12 +95,9 @@ def add_coords_masks_and_apply_corrections( zlims: ZIndexLimits, bdlim: BeamDivergenceLimits, wbins: WavelengthBins, - proton_charge: ProtonCharge[RunType], - monitor: WavelengthMonitor[RunType], - uncertainty_broadcast_mode: UncertaintyBroadcastMode, graph: CoordTransformationGraph[RunType], corrections_to_apply: CorrectionsToApply, -) -> ReducibleData[RunType]: +) -> RunUnnormalizedData[RunType]: """ Computes coordinates, masks and corrections that are the same for the sample measurement and the reference measurement. @@ -83,18 +106,9 @@ def add_coords_masks_and_apply_corrections( da = add_masks(da, ylim, zlims, bdlim, wbins) for correction in corrections_to_apply: - if correction == 'monitor': - da = normalize_by_monitor_histogram( - da, - monitor=monitor, - uncertainty_broadcast_mode=uncertainty_broadcast_mode, - ) - elif correction == 'proton_charge': - da = correct_by_proton_charge(da, proton_charge=proton_charge) - else: - da = correction(da) + da = correction(da) - return ReducibleData[RunType](da) + return RunUnnormalizedData[RunType](da) def correct_by_footprint(da: sc.DataArray) -> sc.DataArray: @@ -102,6 +116,6 @@ def correct_by_footprint(da: sc.DataArray) -> sc.DataArray: return da / sc.sin(da.coords['theta']) -default_corrections = {correct_by_proton_charge, correct_by_footprint} +default_corrections = {correct_by_footprint} providers = (add_coords_masks_and_apply_corrections,) diff --git a/packages/essreflectometry/src/ess/freia/workflow.py b/packages/essreflectometry/src/ess/freia/workflow.py index a4007ba61..352180e59 100644 --- a/packages/essreflectometry/src/ess/freia/workflow.py +++ b/packages/essreflectometry/src/ess/freia/workflow.py @@ -3,9 +3,10 @@ import sciline import scipp as sc +from ess.reduce.uncertainty import UncertaintyBroadcastMode +from ess.reduce.workflow import register_workflow from ..reflectometry import providers as reflectometry_providers -from ..reflectometry.corrections import correct_by_proton_charge from ..reflectometry.types import ( BeamDivergenceLimits, CorrectionsToApply, @@ -25,6 +26,7 @@ normalization, orso, ) +from .corrections import RunNormalization, insert_run_normalization _general_providers = ( *reflectometry_providers, @@ -62,11 +64,11 @@ def mcstas_default_parameters() -> dict: sc.scalar(0.75, unit='deg'), ), SampleRotationOffset[RunType]: sc.scalar(0.0, unit='deg'), - CorrectionsToApply: corrections.default_corrections - - {correct_by_proton_charge}, + CorrectionsToApply: corrections.default_corrections, LookupTableRelativeErrorThreshold: { "detector": 0.06, }, + UncertaintyBroadcastMode: UncertaintyBroadcastMode.drop, } @@ -80,24 +82,99 @@ def default_parameters() -> dict: LookupTableRelativeErrorThreshold: { "multiblade_detector": float('inf'), }, + UncertaintyBroadcastMode: UncertaintyBroadcastMode.drop, } -def FreiaMcStasWorkflow() -> sciline.Pipeline: +def FreiaMcStasWorkflow( + *, + run_norm: RunNormalization = RunNormalization.none, + **kwargs, +) -> sciline.Pipeline: """Workflow for reduction of McStas data for the Freia instrument.""" - workflow = beamline.LoadNeXusWorkflow() + workflow = beamline.LoadNeXusWorkflow(**kwargs) for provider in mcstas_providers: workflow.insert(provider) + insert_run_normalization(workflow, run_norm) for name, param in mcstas_default_parameters().items(): workflow[name] = param return workflow -def FreiaWorkflow() -> sciline.Pipeline: +def FreiaWorkflow( + *, + run_norm: RunNormalization = RunNormalization.proton_charge, + **kwargs, +) -> sciline.Pipeline: """Workflow for reduction of data for the Freia instrument.""" - workflow = beamline.LoadNeXusWorkflow() + workflow = beamline.LoadNeXusWorkflow(**kwargs) for provider in providers: workflow.insert(provider) + insert_run_normalization(workflow, run_norm) for name, param in default_parameters().items(): workflow[name] = param return workflow + + +@register_workflow +def FreiaMcStasUnnormalizedWorkflow() -> sciline.Pipeline: + """Workflow for Freia McStas data without run normalization.""" + return FreiaMcStasWorkflow(run_norm=RunNormalization.none) + + +@register_workflow +def FreiaMcStasMonitorHistogramWorkflow() -> sciline.Pipeline: + """Workflow for Freia McStas data using histogrammed monitor normalization.""" + return FreiaMcStasWorkflow(run_norm=RunNormalization.monitor_histogram) + + +@register_workflow +def FreiaMcStasMonitorIntegratedWorkflow() -> sciline.Pipeline: + """Workflow for Freia McStas data using integrated monitor normalization.""" + return FreiaMcStasWorkflow(run_norm=RunNormalization.monitor_integrated) + + +@register_workflow +def FreiaMcStasProtonChargeWorkflow() -> sciline.Pipeline: + """Workflow for Freia McStas data using proton charge normalization.""" + return FreiaMcStasWorkflow(run_norm=RunNormalization.proton_charge) + + +@register_workflow +def FreiaUnnormalizedWorkflow() -> sciline.Pipeline: + """Workflow for Freia NeXus data without run normalization.""" + return FreiaWorkflow(run_norm=RunNormalization.none) + + +@register_workflow +def FreiaMonitorHistogramWorkflow() -> sciline.Pipeline: + """Workflow for Freia NeXus data using histogrammed monitor normalization.""" + return FreiaWorkflow(run_norm=RunNormalization.monitor_histogram) + + +@register_workflow +def FreiaMonitorIntegratedWorkflow() -> sciline.Pipeline: + """Workflow for Freia NeXus data using integrated monitor normalization.""" + return FreiaWorkflow(run_norm=RunNormalization.monitor_integrated) + + +@register_workflow +def FreiaProtonChargeWorkflow() -> sciline.Pipeline: + """Workflow for Freia NeXus data using proton charge normalization.""" + return FreiaWorkflow(run_norm=RunNormalization.proton_charge) + + +__all__ = [ + 'FreiaMcStasMonitorHistogramWorkflow', + 'FreiaMcStasMonitorIntegratedWorkflow', + 'FreiaMcStasProtonChargeWorkflow', + 'FreiaMcStasUnnormalizedWorkflow', + 'FreiaMcStasWorkflow', + 'FreiaMonitorHistogramWorkflow', + 'FreiaMonitorIntegratedWorkflow', + 'FreiaProtonChargeWorkflow', + 'FreiaUnnormalizedWorkflow', + 'FreiaWorkflow', + 'default_parameters', + 'mcstas_default_parameters', +] diff --git a/packages/essreflectometry/src/ess/reflectometry/corrections.py b/packages/essreflectometry/src/ess/reflectometry/corrections.py index 628ff1e82..ebd604ff7 100644 --- a/packages/essreflectometry/src/ess/reflectometry/corrections.py +++ b/packages/essreflectometry/src/ess/reflectometry/corrections.py @@ -1,8 +1,30 @@ # Copyright (c) 2025 Scipp contributors (https://github.com/scipp) +import enum + +import ess.reduce +import sciline import scipp as sc +from ess.reduce.uncertainty import UncertaintyBroadcastMode from .tools import fwhm_to_std -from .types import RawSampleRotation, RunType, SampleRotation, SampleRotationOffset +from .types import ( + ProtonCharge, + RawSampleRotation, + ReducibleData, + RunType, + RunUnnormalizedData, + SampleRotation, + SampleRotationOffset, +) + + +class RunNormalization(enum.Enum): + """Type of normalization applied to each run.""" + + none = enum.auto() + monitor_histogram = enum.auto() + monitor_integrated = enum.auto() + proton_charge = enum.auto() def footprint_on_sample( @@ -62,6 +84,76 @@ def correct_by_proton_charge( return da / pc +def no_run_normalization( + detector: RunUnnormalizedData[RunType], +) -> ReducibleData[RunType]: + """Use prepared detector data without applying a run normalization.""" + return ReducibleData[RunType](detector) + + +def normalize_by_monitor_histogram( + detector: RunUnnormalizedData[RunType], + *, + monitor: sc.DataArray, + uncertainty_broadcast_mode: UncertaintyBroadcastMode, +) -> ReducibleData[RunType]: + """Normalize detector data by a histogrammed monitor.""" + return ReducibleData[RunType]( + ess.reduce.normalization.normalize_by_monitor_histogram( + detector=detector, + monitor=monitor, + uncertainty_broadcast_mode=uncertainty_broadcast_mode, + skip_range_check=False, + ) + ) + + +def normalize_by_monitor_integrated( + detector: RunUnnormalizedData[RunType], + *, + monitor: sc.DataArray, + uncertainty_broadcast_mode: UncertaintyBroadcastMode, +) -> ReducibleData[RunType]: + """Normalize detector data by an integrated wavelength monitor.""" + return ReducibleData[RunType]( + ess.reduce.normalization.normalize_by_monitor_integrated( + detector=detector, + monitor=monitor, + uncertainty_broadcast_mode=uncertainty_broadcast_mode, + skip_range_check=False, + ) + ) + + +def normalize_by_proton_charge( + detector: RunUnnormalizedData[RunType], + proton_charge: ProtonCharge[RunType], +) -> ReducibleData[RunType]: + """Normalize detector data by time-dependent proton charge.""" + return ReducibleData[RunType]( + correct_by_proton_charge(detector, proton_charge=proton_charge) + ) + + +def insert_run_normalization( + workflow: sciline.Pipeline, + run_norm: RunNormalization, + *, + monitor_histogram_provider, + monitor_integrated_provider, +) -> None: + """Insert providers for a specific normalization into a workflow.""" + match run_norm: + case RunNormalization.none: + workflow.insert(no_run_normalization) + case RunNormalization.monitor_histogram: + workflow.insert(monitor_histogram_provider) + case RunNormalization.monitor_integrated: + workflow.insert(monitor_integrated_provider) + case RunNormalization.proton_charge: + workflow.insert(normalize_by_proton_charge) + + def correct_sample_rotation( mu: RawSampleRotation[RunType], mu_offset: SampleRotationOffset[RunType] ) -> SampleRotation[RunType]: diff --git a/packages/essreflectometry/src/ess/reflectometry/types.py b/packages/essreflectometry/src/ess/reflectometry/types.py index a96902e5c..0e8613a4e 100644 --- a/packages/essreflectometry/src/ess/reflectometry/types.py +++ b/packages/essreflectometry/src/ess/reflectometry/types.py @@ -43,6 +43,10 @@ class ReducibleData(sciline.Scope[RunType, sc.DataArray], sc.DataArray): """Event data with common coordinates added""" +class RunUnnormalizedData(sciline.Scope[RunType, sc.DataArray], sc.DataArray): + """Detector data prepared for reduction, before run normalization.""" + + ReducedReference = NewType("ReducedReference", sc.DataArray) """Intensity distribution on the detector for a sample with :math`R(Q) = 1`""" diff --git a/packages/essreflectometry/tests/estia/mcstas_data_test.py b/packages/essreflectometry/tests/estia/mcstas_data_test.py index ccdb46797..43d71c202 100644 --- a/packages/essreflectometry/tests/estia/mcstas_data_test.py +++ b/packages/essreflectometry/tests/estia/mcstas_data_test.py @@ -9,11 +9,15 @@ import sciline import scipp as sc import scipp.testing -from ess.reduce.normalization import normalize_by_monitor_histogram +from ess.reduce.normalization import ( + normalize_by_monitor_histogram, + normalize_by_monitor_integrated, +) from ess.reduce.uncertainty import UncertaintyBroadcastMode from orsopy import fileio from ess.estia import EstiaMcStasWorkflow +from ess.estia.corrections import RunNormalization from ess.estia.data import ( estia_mcstas_reference_run, estia_mcstas_sample_run, @@ -27,7 +31,6 @@ from ess.reflectometry.corrections import correct_by_proton_charge from ess.reflectometry.types import ( BeamDivergenceLimits, - CorrectionsToApply, Filename, LookupTableFilename, ProtonCharge, @@ -42,9 +45,10 @@ ) -@pytest.fixture -def estia_mcstas_pipeline() -> sciline.Pipeline: - wf = EstiaMcStasWorkflow() +def _make_estia_mcstas_pipeline( + *, run_norm: RunNormalization = RunNormalization.none +) -> sciline.Pipeline: + wf = EstiaMcStasWorkflow(run_norm=run_norm) wf[Filename[ReferenceRun]] = estia_mcstas_reference_run() wf[YIndexLimits] = sc.scalar(35), sc.scalar(64) @@ -90,6 +94,11 @@ def estia_mcstas_pipeline() -> sciline.Pipeline: return wf +@pytest.fixture +def estia_mcstas_pipeline() -> sciline.Pipeline: + return _make_estia_mcstas_pipeline() + + def test_compute_reducible_data(estia_mcstas_pipeline: sciline.Pipeline): estia_mcstas_pipeline[Filename[SampleRun]] = estia_mcstas_sample_run(11) da = estia_mcstas_pipeline.compute(ReducibleData[SampleRun]) @@ -103,19 +112,31 @@ 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]) +@pytest.mark.parametrize( + ("run_norm", "normalize"), + [ + (RunNormalization.monitor_histogram, normalize_by_monitor_histogram), + (RunNormalization.monitor_integrated, normalize_by_monitor_integrated), + ], +) +def test_compute_reducible_data_with_monitor( + estia_mcstas_pipeline: sciline.Pipeline, + run_norm: RunNormalization, + normalize, +): + filename = estia_mcstas_sample_run(11) + estia_mcstas_pipeline[Filename[SampleRun]] = filename + without_monitor = estia_mcstas_pipeline.compute(ReducibleData[SampleRun]) + + wf = _make_estia_mcstas_pipeline(run_norm=run_norm) + wf[Filename[SampleRun]] = filename 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( + normalize( without_monitor, monitor=wf.compute(WavelengthMonitor[SampleRun]), uncertainty_broadcast_mode=UncertaintyBroadcastMode.drop, @@ -127,15 +148,16 @@ def test_compute_reducible_data_with_monitor(estia_mcstas_pipeline: sciline.Pipe def test_compute_reducible_data_with_proton_charge( estia_mcstas_pipeline: sciline.Pipeline, ): - wf = estia_mcstas_pipeline - wf[Filename[SampleRun]] = estia_mcstas_sample_run(11) - without_proton_charge = wf.compute(ReducibleData[SampleRun]) + filename = estia_mcstas_sample_run(11) + estia_mcstas_pipeline[Filename[SampleRun]] = filename + without_proton_charge = estia_mcstas_pipeline.compute(ReducibleData[SampleRun]) + + wf = _make_estia_mcstas_pipeline(run_norm=RunNormalization.proton_charge) + wf[Filename[SampleRun]] = filename wf[ProtonCharge[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_charge'} with_proton_charge = wf.compute(ReducibleData[SampleRun]) scipp.testing.assert_allclose( correct_by_proton_charge( diff --git a/packages/essreflectometry/tests/estia/nexus_data_test.py b/packages/essreflectometry/tests/estia/nexus_data_test.py index d2de80947..156a0904b 100644 --- a/packages/essreflectometry/tests/estia/nexus_data_test.py +++ b/packages/essreflectometry/tests/estia/nexus_data_test.py @@ -9,11 +9,15 @@ import sciline import scipp as sc import scipp.testing -from ess.reduce.normalization import normalize_by_monitor_histogram +from ess.reduce.normalization import ( + normalize_by_monitor_histogram, + normalize_by_monitor_integrated, +) from ess.reduce.uncertainty import UncertaintyBroadcastMode from orsopy import fileio from ess.estia import EstiaWorkflow +from ess.estia.corrections import RunNormalization from ess.estia.data import ( estia_mcstas_nexus_reference_example, estia_mcstas_nexus_sample_example, @@ -23,7 +27,6 @@ from ess.reflectometry import orso, supermirror from ess.reflectometry.types import ( BeamDivergenceLimits, - CorrectionsToApply, Filename, LookupTableFilename, ProtonCharge, @@ -38,9 +41,10 @@ ) -@pytest.fixture -def estia_pipeline() -> sciline.Pipeline: - wf = EstiaWorkflow() +def _make_estia_pipeline( + *, run_norm: RunNormalization = RunNormalization.none +) -> sciline.Pipeline: + wf = EstiaWorkflow(run_norm=run_norm) wf[Filename[ReferenceRun]] = estia_mcstas_nexus_reference_example() wf[YIndexLimits] = sc.scalar(35), sc.scalar(64) @@ -87,11 +91,14 @@ def estia_pipeline() -> sciline.Pipeline: ) ) wf[orso.OrsoSample] = orso.OrsoSample(fileio.data_source.Sample.empty()) - wf[WavelengthMonitor[SampleRun]] = None - wf[WavelengthMonitor[ReferenceRun]] = None return wf +@pytest.fixture +def estia_pipeline() -> sciline.Pipeline: + return _make_estia_pipeline() + + def test_compute_reducible_data(estia_pipeline: sciline.Pipeline): estia_pipeline[Filename[SampleRun]] = estia_mcstas_nexus_sample_example( 'Ni/Ti-multilayer' @@ -107,19 +114,31 @@ def test_compute_reducible_data(estia_pipeline: sciline.Pipeline): 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]) +@pytest.mark.parametrize( + ("run_norm", "normalize"), + [ + (RunNormalization.monitor_histogram, normalize_by_monitor_histogram), + (RunNormalization.monitor_integrated, normalize_by_monitor_integrated), + ], +) +def test_compute_reducible_data_with_monitor( + estia_pipeline: sciline.Pipeline, + run_norm: RunNormalization, + normalize, +): + filename = estia_mcstas_nexus_sample_example('Ni/Ti-multilayer')[0] + estia_pipeline[Filename[SampleRun]] = filename + without_monitor = estia_pipeline.compute(ReducibleData[SampleRun]) + + wf = _make_estia_pipeline(run_norm=run_norm) + wf[Filename[SampleRun]] = filename 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( + normalize( without_monitor, monitor=wf.compute(WavelengthMonitor[SampleRun]), uncertainty_broadcast_mode=UncertaintyBroadcastMode.drop, diff --git a/packages/essreflectometry/tests/estia/workflow_test.py b/packages/essreflectometry/tests/estia/workflow_test.py new file mode 100644 index 000000000..227716151 --- /dev/null +++ b/packages/essreflectometry/tests/estia/workflow_test.py @@ -0,0 +1,41 @@ +# SPDX-License-Identifier: BSD-3-Clause +# Copyright (c) 2026 Scipp contributors (https://github.com/scipp) +import inspect + +from ess.reduce import workflow as reduce_workflow + +from ess import estia +from ess.estia.corrections import RunNormalization +from ess.reflectometry.types import CorrectionsToApply, NeXusDetectorName + + +def test_estia_workflow_uses_estia_defaults(): + params = estia.workflow.default_parameters() + + assert params[NeXusDetectorName] == "multiblade_detector" + assert params[CorrectionsToApply] == estia.corrections.default_corrections + + +def test_estia_workflows_have_expected_run_normalization_defaults(): + assert ( + inspect.signature(estia.EstiaMcStasWorkflow).parameters["run_norm"].default + is RunNormalization.none + ) + assert ( + inspect.signature(estia.EstiaWorkflow).parameters["run_norm"].default + is RunNormalization.proton_charge + ) + + +def test_estia_workflow_registers_run_normalization_variants(): + for wf in ( + estia.EstiaMcStasUnnormalizedWorkflow, + estia.EstiaMcStasMonitorHistogramWorkflow, + estia.EstiaMcStasMonitorIntegratedWorkflow, + estia.EstiaMcStasProtonChargeWorkflow, + estia.EstiaUnnormalizedWorkflow, + estia.EstiaMonitorHistogramWorkflow, + estia.EstiaMonitorIntegratedWorkflow, + estia.EstiaProtonChargeWorkflow, + ): + assert wf in reduce_workflow.workflow_registry diff --git a/packages/essreflectometry/tests/freia/workflow_test.py b/packages/essreflectometry/tests/freia/workflow_test.py index 9927a6484..41bc7a887 100644 --- a/packages/essreflectometry/tests/freia/workflow_test.py +++ b/packages/essreflectometry/tests/freia/workflow_test.py @@ -1,7 +1,11 @@ # SPDX-License-Identifier: BSD-3-Clause # Copyright (c) 2026 Scipp contributors (https://github.com/scipp) +import inspect + +from ess.reduce import workflow as reduce_workflow from ess import freia +from ess.freia.corrections import RunNormalization from ess.reflectometry.types import CorrectionsToApply, NeXusDetectorName @@ -10,3 +14,28 @@ def test_freia_workflow_uses_freia_defaults(): assert params[NeXusDetectorName] == "multiblade_detector" assert params[CorrectionsToApply] == freia.corrections.default_corrections + + +def test_freia_workflows_have_expected_run_normalization_defaults(): + assert ( + inspect.signature(freia.FreiaMcStasWorkflow).parameters["run_norm"].default + is RunNormalization.none + ) + assert ( + inspect.signature(freia.FreiaWorkflow).parameters["run_norm"].default + is RunNormalization.proton_charge + ) + + +def test_freia_workflow_registers_run_normalization_variants(): + for wf in ( + freia.FreiaMcStasUnnormalizedWorkflow, + freia.FreiaMcStasMonitorHistogramWorkflow, + freia.FreiaMcStasMonitorIntegratedWorkflow, + freia.FreiaMcStasProtonChargeWorkflow, + freia.FreiaUnnormalizedWorkflow, + freia.FreiaMonitorHistogramWorkflow, + freia.FreiaMonitorIntegratedWorkflow, + freia.FreiaProtonChargeWorkflow, + ): + assert wf in reduce_workflow.workflow_registry