Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Estia #108

Open
wants to merge 35 commits into
base: main
Choose a base branch
from
Open

Estia #108

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
35 commits
Select commit Hold shift + click to select a range
8c26f15
Create estia-data-reduction.md
jokasimr Nov 27, 2024
b10a2c0
feat: initial setup for estia workflow
jokasimr Dec 10, 2024
f84db80
feat: compute angle without divergence angle per pixel
jokasimr Jan 29, 2025
b643492
fix: dont modify input + dont keep intermediates
jokasimr Feb 5, 2025
baca44b
fix: remove operation modifying input
jokasimr Feb 5, 2025
4e5fb07
feat: estia coordinate transformations + workflow
jokasimr Feb 5, 2025
43e63c9
WIP: estia workflow - changes to the generic part of the workflow
jokasimr Feb 11, 2025
2b32c7e
WIP: estia workflow
jokasimr Feb 11, 2025
aa4b6a0
fix
jokasimr Feb 17, 2025
10ac274
wip
jokasimr Feb 18, 2025
ec0b317
feat: add mcstas h5 file loader
jokasimr Feb 18, 2025
e8b23f6
fix: allow loading ascii and h5 files
jokasimr Feb 18, 2025
c3a6afe
Apply automatic formatting
pre-commit-ci-lite[bot] Feb 17, 2025
dd1b7c0
fix: reshape to expected dimensions
jokasimr Feb 19, 2025
9ac64dc
test: add mcstas test files and tests
jokasimr Feb 20, 2025
95cf51e
minor
jokasimr Feb 27, 2025
4db26a8
docs: removing from this pr because documentation needs some rewriting
jokasimr Feb 27, 2025
df83bfe
refactor: share parts of the workflow between amor and estia
jokasimr Feb 27, 2025
e516537
fix: remove comment and add providers to load
jokasimr Feb 28, 2025
8ab311a
docs: typehints
jokasimr Feb 28, 2025
e03119e
docs: add parameters in docstrings
jokasimr Feb 28, 2025
9481703
fix: mod t to get event_time_offset + remove unnecessary coordinates
jokasimr Feb 28, 2025
ac1063e
docs: update licence headers
jokasimr Feb 28, 2025
8ece26c
docs: add licence header
jokasimr Feb 28, 2025
6539130
fix: add comment explaining missing value and set better default
jokasimr Mar 3, 2025
a427864
docs: add docstring
jokasimr Mar 3, 2025
a27676d
docs: add type hints
jokasimr Mar 3, 2025
1f9f48c
fix: rename test file
jokasimr Mar 3, 2025
40f641f
fix: rename test
jokasimr Mar 3, 2025
2dea495
Merge branch 'main' into estia
jokasimr Mar 4, 2025
d9610b8
fix: remove center of incident angle concept
jokasimr Mar 4, 2025
7b68c2c
Merge branch 'estia' of github.com:scipp/essreflectometry into estia
jokasimr Mar 4, 2025
7efcc91
Update src/ess/estia/conversions.py
jokasimr Mar 4, 2025
fbc77b1
refactor: move masking
jokasimr Mar 4, 2025
9d48bc5
Merge branch 'estia' of github.com:scipp/essreflectometry into estia
jokasimr Mar 4, 2025
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion src/ess/amor/__init__.py
Original file line number Diff line number Diff line change
@@ -1,5 +1,5 @@
# SPDX-License-Identifier: BSD-3-Clause
# Copyright (c) 2023 Scipp contributors (https://github.com/scipp)
# Copyright (c) 2025 Scipp contributors (https://github.com/scipp)
import importlib.metadata

import sciline
Expand Down
38 changes: 12 additions & 26 deletions src/ess/amor/conversions.py
Original file line number Diff line number Diff line change
@@ -1,20 +1,20 @@
# SPDX-License-Identifier: BSD-3-Clause
# Copyright (c) 2023 Scipp contributors (https://github.com/scipp)
# Copyright (c) 2025 Scipp contributors (https://github.com/scipp)
import numpy as np
import scipp as sc

from ..reflectometry.conversions import reflectometry_q
from ..reflectometry.types import (
BeamDivergenceLimits,
CoordTransformationGraph,
WavelengthBins,
YIndexLimits,
ZIndexLimits,
)
from .geometry import Detector
from .types import CoordTransformationGraph


def theta(wavelength, divergence_angle, L2, sample_rotation, detector_rotation):
def theta(wavelength, pixel_divergence_angle, L2, sample_rotation, detector_rotation):
'''
Angle of reflection.

Expand Down Expand Up @@ -61,14 +61,15 @@ def theta(wavelength, divergence_angle, L2, sample_rotation, detector_rotation):
'''
c = sc.constants.g * sc.constants.m_n**2 / sc.constants.h**2
out = (c * L2 * wavelength**2).to(unit='dimensionless') + sc.sin(
divergence_angle.to(unit='rad', copy=False) + detector_rotation.to(unit='rad')
pixel_divergence_angle.to(unit='rad', copy=False)
+ detector_rotation.to(unit='rad')
)
out = sc.asin(out, out=out)
out -= sample_rotation.to(unit='rad')
return out


def angle_of_divergence(theta, sample_rotation, angle_to_center_of_beam):
def divergence_angle(theta, sample_rotation, angle_to_center_of_beam):
"""
Difference between the incident angle and the center of the incident beam.
Useful for filtering parts of the beam that have too high divergence.
Expand All @@ -84,7 +85,7 @@ def angle_of_divergence(theta, sample_rotation, angle_to_center_of_beam):


def wavelength(
event_time_offset, divergence_angle, L1, L2, chopper_phase, chopper_frequency
event_time_offset, pixel_divergence_angle, L1, L2, chopper_phase, chopper_frequency
):
"Converts event_time_offset to wavelength using the chopper settings."
out = event_time_offset.to(unit="ns", dtype="float64", copy=True)
Expand All @@ -108,37 +109,22 @@ def wavelength(
)
# Correction for path length through guides being different
# depending on incident angle.
out -= (divergence_angle.to(unit="rad") / (np.pi * sc.units.rad)) * tau
out -= (pixel_divergence_angle.to(unit="rad") / (np.pi * sc.units.rad)) * tau
out *= (sc.constants.h / sc.constants.m_n) / (L1 + L2)
return out.to(unit='angstrom', copy=False)


def coordinate_transformation_graph() -> CoordTransformationGraph:
return {
"divergence_angle": "pixel_divergence_angle",
"wavelength": wavelength,
"theta": theta,
"angle_of_divergence": angle_of_divergence,
"divergence_angle": divergence_angle,
"Q": reflectometry_q,
"L1": lambda chopper_distance: sc.abs(chopper_distance),
"L2": lambda distance_in_detector: distance_in_detector + Detector.distance,
}


def add_coords(
da: sc.DataArray,
graph: dict,
) -> sc.DataArray:
"Adds scattering coordinates to the raw detector data."
return da.transform_coords(
("wavelength", "theta", "angle_of_divergence", "Q", "L1", "L2"),
graph,
rename_dims=False,
keep_intermediate=False,
keep_aliases=False,
)


def _not_between(v, a, b):
return (v < a) | (v > b)

Expand All @@ -161,9 +147,9 @@ def add_masks(
)
da = da.bins.assign_masks(
divergence_too_large=_not_between(
da.bins.coords["angle_of_divergence"],
bdlim[0].to(unit=da.bins.coords["angle_of_divergence"].bins.unit),
bdlim[1].to(unit=da.bins.coords["angle_of_divergence"].bins.unit),
da.bins.coords["divergence_angle"],
bdlim[0].to(unit=da.bins.coords["divergence_angle"].bins.unit),
bdlim[1].to(unit=da.bins.coords["divergence_angle"].bins.unit),
),
wavelength=_not_between(
da.bins.coords['wavelength'],
Expand Down
2 changes: 1 addition & 1 deletion src/ess/amor/data.py
Original file line number Diff line number Diff line change
@@ -1,5 +1,5 @@
# SPDX-License-Identifier: BSD-3-Clause
# Copyright (c) 2023 Scipp contributors (https://github.com/scipp)
# Copyright (c) 2025 Scipp contributors (https://github.com/scipp)

from ..reflectometry.types import Filename, ReferenceRun, SampleRun

Expand Down
1 change: 1 addition & 0 deletions src/ess/amor/figures.py
Original file line number Diff line number Diff line change
@@ -1,3 +1,4 @@
# Copyright (c) 2025 Scipp contributors (https://github.com/scipp)
from collections.abc import Sequence

import numpy as np
Expand Down
2 changes: 1 addition & 1 deletion src/ess/amor/load.py
Original file line number Diff line number Diff line change
@@ -1,5 +1,5 @@
# SPDX-License-Identifier: BSD-3-Clause
# Copyright (c) 2023 Scipp contributors (https://github.com/scipp)
# Copyright (c) 2025 Scipp contributors (https://github.com/scipp)
import scipp as sc

from ..reflectometry.load import load_nx
Expand Down
10 changes: 5 additions & 5 deletions src/ess/amor/normalization.py
Original file line number Diff line number Diff line change
@@ -1,3 +1,4 @@
# Copyright (c) 2025 Scipp contributors (https://github.com/scipp)
import scipp as sc

from ..reflectometry.conversions import reflectometry_q
Expand All @@ -8,8 +9,8 @@
supermirror_reflectivity,
)
from ..reflectometry.types import (
CoordTransformationGraph,
DetectorSpatialResolution,
QBins,
ReducedReference,
ReducibleData,
Reference,
Expand All @@ -23,7 +24,6 @@
sample_size_resolution,
wavelength_resolution,
)
from .types import CoordTransformationGraph


def mask_events_where_supermirror_does_not_cover(
Expand Down Expand Up @@ -64,14 +64,12 @@ def mask_events_where_supermirror_does_not_cover(
m=mvalue,
alpha=alpha,
)
sam.bins.masks["supermirror_does_not_cover"] = sc.isnan(R)
return sam
return sam.bins.assign_masks(supermirror_does_not_cover=sc.isnan(R))


def evaluate_reference_at_sample_coords(
reference: ReducedReference,
sample: ReducibleData[SampleRun],
qbins: QBins,
detector_spatial_resolution: DetectorSpatialResolution[SampleRun],
graph: CoordTransformationGraph,
) -> Reference:
Expand Down Expand Up @@ -103,6 +101,8 @@ def evaluate_reference_at_sample_coords(
"Q_resolution": q_resolution,
},
rename_dims=False,
keep_intermediate=False,
keep_aliases=False,
)
return sc.values(ref)

Expand Down
2 changes: 1 addition & 1 deletion src/ess/amor/resolution.py
Original file line number Diff line number Diff line change
@@ -1,5 +1,5 @@
# SPDX-License-Identifier: BSD-3-Clause
# Copyright (c) 2023 Scipp contributors (https://github.com/scipp)
# Copyright (c) 2025 Scipp contributors (https://github.com/scipp)
import scipp as sc

from ..reflectometry.tools import fwhm_to_std
Expand Down
3 changes: 1 addition & 2 deletions src/ess/amor/types.py
Original file line number Diff line number Diff line change
@@ -1,3 +1,4 @@
# Copyright (c) 2025 Scipp contributors (https://github.com/scipp)
from typing import Any, NewType

import sciline
Expand All @@ -9,8 +10,6 @@
AngularResolution = NewType("AngularResolution", sc.Variable)
SampleSizeResolution = NewType("SampleSizeResolution", sc.Variable)

CoordTransformationGraph = NewType("CoordTransformationGraph", dict)


class ChopperFrequency(sciline.Scope[RunType, sc.Variable], sc.Variable):
"""Frequency of the choppers used in the run."""
Expand Down
1 change: 1 addition & 0 deletions src/ess/amor/utils.py
Original file line number Diff line number Diff line change
@@ -1,3 +1,4 @@
# Copyright (c) 2025 Scipp contributors (https://github.com/scipp)
import scipp as sc

from ess.reflectometry.types import (
Expand Down
6 changes: 4 additions & 2 deletions src/ess/amor/workflow.py
Original file line number Diff line number Diff line change
@@ -1,10 +1,13 @@
# Copyright (c) 2025 Scipp contributors (https://github.com/scipp)
from ..reflectometry.conversions import (
add_coords,
add_proton_current_coord,
add_proton_current_mask,
)
from ..reflectometry.corrections import correct_by_footprint, correct_by_proton_current
from ..reflectometry.types import (
BeamDivergenceLimits,
CoordTransformationGraph,
ProtonCurrent,
RawDetectorData,
ReducibleData,
Expand All @@ -13,8 +16,7 @@
YIndexLimits,
ZIndexLimits,
)
from .conversions import add_coords, add_masks
from .types import CoordTransformationGraph
from .conversions import add_masks


def add_coords_masks_and_apply_corrections(
Expand Down
80 changes: 80 additions & 0 deletions src/ess/estia/__init__.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,80 @@
# SPDX-License-Identifier: BSD-3-Clause
# Copyright (c) 2025 Scipp contributors (https://github.com/scipp)
import importlib.metadata

import sciline
import scipp as sc

from ..reflectometry import providers as reflectometry_providers
from ..reflectometry import supermirror
from ..reflectometry.types import (
BeamDivergenceLimits,
BeamSize,
DetectorSpatialResolution,
NeXusDetectorName,
RunType,
SamplePosition,
)
from . import conversions, load, maskings, normalization, resolution, workflow
from .types import (
AngularResolution,
SampleSizeResolution,
WavelengthResolution,
)

try:
__version__ = importlib.metadata.version(__package__ or __name__)
except importlib.metadata.PackageNotFoundError:
__version__ = "0.0.0"


providers = (
*reflectometry_providers,
*load.providers,
*conversions.providers,
*maskings.providers,
*workflow.providers,
*normalization.providers,
)
"""
List of providers for setting up a Sciline pipeline.

This provides a default Estia workflow including providers for loadings files.
"""


def default_parameters() -> dict:
return {
supermirror.MValue: sc.scalar(5, unit=sc.units.dimensionless),
supermirror.CriticalEdge: 0.022 * sc.Unit("1/angstrom"),
supermirror.Alpha: sc.scalar(0.25 / 0.088, unit=sc.units.angstrom),
BeamSize[RunType]: 2.0 * sc.units.mm,
DetectorSpatialResolution[RunType]: 0.0025 * sc.units.m,
SamplePosition[RunType]: sc.vector([0, 0, 0], unit="m"),
NeXusDetectorName[RunType]: "detector",
BeamDivergenceLimits: (
sc.scalar(-0.75, unit='deg'),
sc.scalar(0.75, unit='deg'),
),
}


def EstiaWorkflow() -> sciline.Pipeline:
"""
Workflow with default parameters for the Estia instrument.
"""
return sciline.Pipeline(providers=providers, params=default_parameters())


__all__ = [
"AngularResolution",
"EstiaWorkflow",
"SampleSizeResolution",
"WavelengthResolution",
"conversions",
"default_parameters",
"load",
"providers",
"resolution",
"supermirror",
]
83 changes: 83 additions & 0 deletions src/ess/estia/conversions.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,83 @@
# SPDX-License-Identifier: BSD-3-Clause
# Copyright (c) 2025 Scipp contributors (https://github.com/scipp)
import scipp as sc

from ..reflectometry.conversions import reflectometry_q
from ..reflectometry.types import (
CoordTransformationGraph,
)


def theta(
divergence_angle: sc.Variable,
sample_rotation: sc.Variable,
):
'''
Angle of reflection.

Computes the angle between the scattering direction of
the neutron and the sample surface.

Parameters
------------
divergence_angle:
Divergence angle of the scattered beam.
sample_rotation:
Rotation of the sample from to its zero position.

Returns
-----------
The reflection angle of the neutron.
'''
Copy link
Member

Choose a reason for hiding this comment

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

Can you also add descriptions of the parameters in the docstring?

return divergence_angle + sample_rotation.to(unit=divergence_angle.unit)


def divergence_angle(
position: sc.Variable,
sample_position: sc.Variable,
detector_rotation: sc.Variable,
):
"""
Angle between the scattering ray and
the ray that travels parallel to the sample surface
when the sample rotation is zero.

Parameters
------------
position:
Detector position where the neutron was detected.
sample_position:
Position of the sample.
detector_rotation:
Rotation of the detector from its zero position.
Returns
----------
The divergence angle of the scattered beam.
"""
p = position - sample_position.to(unit=position.unit)
return sc.atan2(y=p.fields.x, x=p.fields.z) - detector_rotation.to(unit='rad')


def wavelength(
event_time_offset,
# Other inputs
):
"Converts event_time_offset to wavelength"
# Use frame unwrapping from scippneutron
raise NotImplementedError()


def coordinate_transformation_graph() -> CoordTransformationGraph:
return {
"wavelength": wavelength,
"theta": theta,
"divergence_angle": divergence_angle,
"Q": reflectometry_q,
"L1": lambda source_position, sample_position: sc.norm(
sample_position - source_position
), # + extra correction for guides?
"L2": lambda position, sample_position: sc.norm(position - sample_position),
Comment on lines +76 to +79
Copy link
Member

Choose a reason for hiding this comment

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

Can we use/import this from scippneutron, thus removing the need to re-implement here?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

I know this is something that's come up before, but I really don't think it's useful to import these two tiny definitions from a separate package.

  1. If they change they are not going to change "globally" so that we want to apply those changes to every instrument simultaneously.
  2. It's also something the user or instrument scientists want's to see and understand, and finding where this is done in scippneutron is hard.

Besides, it's plausible that the L1 definition will change to take into account the non-straight (focused) beam.
That's done for Amor.

}


providers = (coordinate_transformation_graph,)
Loading