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 a83dc2f49..ebaf4b47c 100644 --- a/packages/essreflectometry/src/ess/estia/corrections.py +++ b/packages/essreflectometry/src/ess/estia/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) @@ -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'} 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..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: @@ -223,10 +225,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: @@ -293,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]: @@ -302,6 +307,8 @@ def use_mcstas_wavelengths_instead_of_estimates_from_time_of_arrival( zlims=zlims, bdlim=bdlim, wbins=wbins, + monitor=monitor, + uncertainty_broadcast_mode=uncertainty_broadcast_mode, 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..cf144ac3d 100644 --- a/packages/essreflectometry/src/ess/estia/workflow.py +++ b/packages/essreflectometry/src/ess/estia/workflow.py @@ -5,10 +5,10 @@ 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 -from ..reflectometry.corrections import correct_by_proton_current from ..reflectometry.types import ( BeamDivergenceLimits, CorrectionsToApply, @@ -28,6 +28,7 @@ normalization, orso, ) +from .types import WavelengthMonitor _general_providers = ( *reflectometry_providers, @@ -69,11 +70,14 @@ 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_current}, + CorrectionsToApply: ( + corrections.default_corrections - {'monitor', 'proton_current'} + ), LookupTableRelativeErrorThreshold: { "multiblade_detector": 0.06, }, + WavelengthMonitor[RunType]: None, + UncertaintyBroadcastMode: UncertaintyBroadcastMode.drop, } @@ -82,11 +86,13 @@ 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', 'proton_current'}, DetectorSpatialResolution: 0.0025 * sc.units.m, LookupTableRelativeErrorThreshold: { "multiblade_detector": float('inf'), }, + UncertaintyBroadcastMode: UncertaintyBroadcastMode.drop, } 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] 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 0edf9e398..db6c58b66 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,12 @@ 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.corrections import correct_by_proton_current from ess.reflectometry.types import ( BeamDivergenceLimits, + CorrectionsToApply, Filename, LookupTableFilename, ProtonCurrent, @@ -97,6 +103,49 @@ 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_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, ): 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' + )