Skip to content
8 changes: 1 addition & 7 deletions packages/essreflectometry/src/ess/amor/workflow.py
Original file line number Diff line number Diff line change
@@ -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,
Expand Down Expand Up @@ -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)

Expand Down
67 changes: 57 additions & 10 deletions packages/essreflectometry/src/ess/estia/corrections.py
Original file line number Diff line number Diff line change
@@ -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,
Expand All @@ -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(
Expand All @@ -30,6 +70,8 @@ def add_coords_masks_and_apply_corrections(
bdlim: BeamDivergenceLimits,
wbins: WavelengthBins,
proton_current: ProtonCurrent[RunType],
monitor: WavelengthMonitor[RunType],
Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is a problem. You now have to always load and process both the monitor and proton charge even if you don't need them. Plus, this means that the graph visualization is not enough to understand the data flow.
So I still think normalisations should be selected on the graph level, not as a parameter.

uncertainty_broadcast_mode: UncertaintyBroadcastMode,
graph: CoordTransformationGraph[RunType],
corrections_to_apply: CorrectionsToApply,
) -> ReducibleData[RunType]:
Expand All @@ -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)

Expand All @@ -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'}
Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
default_corrections = {correct_by_footprint, 'proton_current'}
default_corrections = {correct_by_footprint}

because the default is set in the workflows without proton current correction.

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The Estia workflow should have proton current correction. Right now it's unclear how to load that information from the nexus files, but I don't think we should adjust the default parameters because of that.


providers = (add_coords_masks_and_apply_corrections,)
13 changes: 10 additions & 3 deletions packages/essreflectometry/src/ess/estia/mcstas.py
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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:
Expand Down Expand Up @@ -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:
Expand Down Expand Up @@ -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]:
Expand All @@ -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,
Expand Down
8 changes: 7 additions & 1 deletion packages/essreflectometry/src/ess/estia/types.py
Original file line number Diff line number Diff line change
@@ -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]
14 changes: 10 additions & 4 deletions packages/essreflectometry/src/ess/estia/workflow.py
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand All @@ -28,6 +28,7 @@
normalization,
orso,
)
from .types import WavelengthMonitor

_general_providers = (
*reflectometry_providers,
Expand Down Expand Up @@ -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,
}


Expand All @@ -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,
}


Expand Down
65 changes: 56 additions & 9 deletions packages/essreflectometry/src/ess/freia/corrections.py
Original file line number Diff line number Diff line change
@@ -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,
Expand All @@ -21,6 +19,48 @@
)
from .conversions import add_coords
from .maskings import add_masks
from .types import WavelengthMonitor


def normalize_by_monitor_histogram(
Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Do we need this copy of the estia function? Can you move this to reflectometry.corrections?

Copy link
Copy Markdown
Contributor Author

@jokasimr jokasimr May 5, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Freia and Estia don't have the same WavelengthMonitor types. Freia has several monitors, while Estia only has one. I though it was less hassle to just duplicate it for now. This will definitely change when we get real requirements about how Freia wants to do monitor normalization.

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Then the monitor should be parametrized over the MonitorType. So can you replace the custom type with the one define din essreduce?
If you think it's too verbose to always specify the type when using an estia workflow, you can provide an alias for the type there. But we should avoid duplicating code just because typing [Monitor] is too much.

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

But that's what is already done right? I mean, there's a type alias for the essreduce monitor in estia/types.py.

The point is: I have no idea what the monitor normalization for Freia will look like, what monitor it will use, etc. It doesn't make any sense to do anything there right now. It's just a placeholder.

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(
Expand All @@ -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]:
Expand All @@ -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)

Expand Down
8 changes: 7 additions & 1 deletion packages/essreflectometry/src/ess/freia/types.py
Original file line number Diff line number Diff line change
@@ -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]
Loading
Loading