diff --git a/packages/esssans/docs/user-guide/common/beam-center-finder.ipynb b/packages/esssans/docs/user-guide/common/beam-center-finder.ipynb index 736fc73fe..887497cdb 100644 --- a/packages/esssans/docs/user-guide/common/beam-center-finder.ipynb +++ b/packages/esssans/docs/user-guide/common/beam-center-finder.ipynb @@ -188,6 +188,10 @@ "id": "13", "metadata": {}, "source": [ + "The result is the beam center relative to the sample.\n", + "Its component along the beam is the distance from the sample at which the transverse offset was determined, here the distance of the detector panel.\n", + "That distance is part of the result because the same beam direction corresponds to a proportionally smaller transverse offset on a detector closer to the sample.\n", + "\n", "We can now update our previous figure with the new position of the beam center (pink dot),\n", "which is clearly in the center of the beam/beam stop." ] @@ -443,16 +447,28 @@ "metadata": {}, "outputs": [], "source": [ + "from ess.sans.beam_center_finder import _BeamPlane\n", + "from ess.sans.conversions import ElasticCoordTransformGraph\n", + "\n", "workflow = workflow.copy()\n", "workflow[QBins] = q_bins\n", "workflow[ReturnEvents] = False\n", "workflow[DimsToKeep] = ()\n", "workflow[WavelengthMask] = None\n", "workflow[WavelengthBands] = None\n", + "# Offsets inside the plane normal to the beam only define a beam center together with\n", + "# the distance from the sample at which they are given. We take that distance from the\n", + "# center-of-mass estimate above.\n", + "plane = _BeamPlane.from_beam_center(\n", + " data=workflow.compute(CorrectedDetector[SampleRun, Numerator]),\n", + " graph=workflow.compute(ElasticCoordTransformGraph[SampleRun]),\n", + " beam_center=com,\n", + ")\n", "kwargs = dict( # noqa: C408\n", " workflow=workflow,\n", " detector=detector['data'],\n", " norm=workflow.compute(CleanDirectBeam),\n", + " plane=plane,\n", ")" ] }, diff --git a/packages/esssans/src/ess/isissans/general.py b/packages/esssans/src/ess/isissans/general.py index ee0127052..b8d4968c6 100644 --- a/packages/esssans/src/ess/isissans/general.py +++ b/packages/esssans/src/ess/isissans/general.py @@ -12,7 +12,6 @@ import scippneutron as scn import scippnexus as snx from ess.sans.types import ( - BeamCenter, DetectorIDs, DetectorPixelShape, DetectorPositionOffset, @@ -81,9 +80,9 @@ def default_parameters() -> dict: def to_detector_position_offset( - global_offset: DetectorBankOffset, beam_center: BeamCenter + global_offset: DetectorBankOffset, ) -> DetectorPositionOffset[RunType]: - return DetectorPositionOffset[RunType](global_offset - beam_center) + return DetectorPositionOffset[RunType](global_offset) def to_monitor_position_offset( diff --git a/packages/esssans/src/ess/isissans/sans2d.py b/packages/esssans/src/ess/isissans/sans2d.py index 6f839785a..119ad1626 100644 --- a/packages/esssans/src/ess/isissans/sans2d.py +++ b/packages/esssans/src/ess/isissans/sans2d.py @@ -6,7 +6,7 @@ import scipp as sc from ess.sans import SansWorkflow from ess.sans.parameters import typical_outputs -from ess.sans.types import BeamCenter, DetectorMasks, EmptyDetector, SampleRun +from ess.sans.types import DetectorMasks, EmptyDetector, SampleRun from ess.reduce.workflow import register_workflow @@ -25,33 +25,26 @@ """Sample holder mask""" -def detector_edge_mask( - beam_center: BeamCenter, sample: EmptyDetector[SampleRun] -) -> DetectorEdgeMask: - # These values were determined by hand before the beam center was available. - # We therefore undo the shift introduced by the beam center. - raw_pos = sample.coords['position'] + beam_center - mask_edges = (sc.abs(raw_pos.fields.x) > sc.scalar(0.48, unit='m')) | ( - sc.abs(raw_pos.fields.y) > sc.scalar(0.45, unit='m') +def detector_edge_mask(sample: EmptyDetector[SampleRun]) -> DetectorEdgeMask: + pos = sample.coords['position'] + mask_edges = (sc.abs(pos.fields.x) > sc.scalar(0.48, unit='m')) | ( + sc.abs(pos.fields.y) > sc.scalar(0.45, unit='m') ) return DetectorEdgeMask(mask_edges) def sample_holder_mask( - beam_center: BeamCenter, sample: EmptyDetector[SampleRun], low_counts_threshold: LowCountThreshold, ) -> SampleHolderMask: - # These values were determined by hand before the beam center was available. - # We therefore undo the shift introduced by the beam center. - raw_pos = sample.coords['position'] + beam_center + pos = sample.coords['position'] summed = sample.hist() holder_mask = ( (summed.data < low_counts_threshold) - & (raw_pos.fields.x > sc.scalar(0, unit='m')) - & (raw_pos.fields.x < sc.scalar(0.42, unit='m')) - & (raw_pos.fields.y < sc.scalar(0.05, unit='m')) - & (raw_pos.fields.y > sc.scalar(-0.15, unit='m')) + & (pos.fields.x > sc.scalar(0, unit='m')) + & (pos.fields.x < sc.scalar(0.42, unit='m')) + & (pos.fields.y < sc.scalar(0.05, unit='m')) + & (pos.fields.y > sc.scalar(-0.15, unit='m')) ) return SampleHolderMask(holder_mask) diff --git a/packages/esssans/src/ess/loki/larmor_workflow.py b/packages/esssans/src/ess/loki/larmor_workflow.py index 664feee8f..0ee8d66e4 100644 --- a/packages/esssans/src/ess/loki/larmor_workflow.py +++ b/packages/esssans/src/ess/loki/larmor_workflow.py @@ -108,12 +108,14 @@ def larmor_detector_coord_transform_graph( *, sample_position: Position[snx.NXsample, RunType], source_position: Position[snx.NXsource, RunType], + beam_center: BeamCenter, gravity: GravityVector, ) -> ElasticCoordTransformGraph[RunType]: graph = sans_elastic( correct_for_gravity=correct_for_gravity, sample_position=sample_position, source_position=source_position, + beam_center=beam_center, gravity=gravity, ) return ElasticCoordTransformGraph[RunType]({**graph, **tof.elastic_Q('tof')}) @@ -186,5 +188,9 @@ def LokiAtLarmorTutorialWorkflow() -> sciline.Pipeline: data.loki_tutorial_run_60392() ) workflow[Filename[EmptyBeamRun]] = str(data.loki_tutorial_run_60392()) - workflow[BeamCenter] = sc.vector(value=[-0.02914868, -0.01816138, 0.0], unit='m') + # Beam center determined with `beam_center_from_center_of_mass`, the Z component + # being the distance from the sample at which it was determined. + workflow[BeamCenter] = sc.vector( + value=[-0.02914868, -0.01816138, 4.11610], unit='m' + ) return workflow diff --git a/packages/esssans/src/ess/sans/beam_center_finder.py b/packages/esssans/src/ess/sans/beam_center_finder.py index 30578ba2b..746ea85d2 100644 --- a/packages/esssans/src/ess/sans/beam_center_finder.py +++ b/packages/esssans/src/ess/sans/beam_center_finder.py @@ -2,6 +2,7 @@ # Copyright (c) 2023 Scipp contributors (https://github.com/scipp) import uuid +from dataclasses import dataclass import numpy as np import sciline @@ -120,7 +121,6 @@ def _find_beam_center( > slope * distance_to_center.fields.x - beam_stop_arm_width ) - center.data.fields.z = sc.scalar(0.0, unit=center.data.unit) return center.data @@ -178,13 +178,34 @@ def beam_center_from_center_of_mass_alternative( if beam_stop_arm_width is None: beam_stop_arm_width = sc.scalar(0.02, unit='m') + workflow = _with_default_beam_center(workflow) + data = workflow.compute(CorrectedDetector[SampleRun, Numerator]) + graph = workflow.compute(ElasticCoordTransformGraph[SampleRun]) + com = _find_beam_center(data, beam_stop_radius, beam_stop_arm_width) + return BeamCenter(com - _sample_position(data, graph)) + + +def _sample_position(data: sc.DataArray, graph: dict) -> sc.Variable: + """Sample position as defined by the coordinate transformation graph.""" + return data.transform_coords('sample_position', graph=graph).coords[ + 'sample_position' + ] + + +def _with_default_beam_center(workflow: sciline.Pipeline) -> sciline.Pipeline: + """ + Return a workflow that can compute detector data even without a beam center set. + + Computing the detector data may require a beam center, but the beam center found by + the functions below does not depend on it: it is derived from the pixel positions, + which the beam center does not modify. + """ try: - beam_center = workflow.compute(BeamCenter) + workflow.compute(BeamCenter) except sciline.UnsatisfiedRequirement: - beam_center = sc.vector([0.0, 0.0, 0.0], unit='m') - workflow[BeamCenter] = beam_center - data = workflow.compute(CorrectedDetector[SampleRun, Numerator]) - return _find_beam_center(data, beam_stop_radius, beam_stop_arm_width) + beam_center + workflow = workflow.copy() + workflow[BeamCenter] = sc.vector([0.0, 0.0, 0.0], unit='m') + return workflow def beam_center_from_center_of_mass(workflow: sciline.Pipeline) -> BeamCenter: @@ -194,9 +215,37 @@ def beam_center_from_center_of_mass(workflow: sciline.Pipeline) -> BeamCenter: The center-of-mass is simply the weighted mean of the positions. Areas with low counts are excluded from the center of mass calculation, as they typically fall into asymmetric regions of the detector panel and would thus lead - to a biased result. The beam is assumed to be roughly aligned with the Z axis. - The returned beam center is the component normal to the beam direction, projected - onto the X-Y plane. + to a biased result. + + The result is the center-of-mass relative to the sample, i.e., the transverse offset + of the beam together with the distance from the sample at which that offset was + determined. See :py:class:`ess.sans.types.BeamCenter`. + + On a bank with depth, the reported distance is the intensity-weighted mean depth, + not the geometric center of the bank: layers closer to the sample see more counts + and pull it forward. That is consistent rather than wrong. The transverse position + of the beam is linear in the distance from the sample, so any intensity-weighted + average of points on the beam is again a point on the beam. The weighting moves the + result along the beam, not away from it, and the reported distance says where it + ended up. + + This holds as long as the intensity is symmetric around the beam at every depth. It + is not exactly true: masking, the shadow of the beam stop and, with gravity, the + wavelength-dependent drop all bias the center-of-mass, and they do so differently at + different depths. The resulting uncertainty on the distance is a few centimetres for + a Loki bank, which is of the order of the depth of the bank itself. The distance is + used to extrapolate the transverse offset to other banks, where an error of a few + centimetres out of several metres is negligible, but it means the variation of the + correction *within* one bank is not resolved any better than the noise on it. + + The result is anchored at the sample position, so an error in the sample position at + the time of determination tilts the inferred beam direction by that error divided by + the distance. One millimetre of transverse error at five metres is 0.2 mrad, which + is 0.24 mm of transverse error on a bank one metre from the sample. Errors along the + beam are second order and negligible. Note that this is an error on the beam + direction and does not require the sample to be in the same place for the run the + beam center is applied to: the direction is re-anchored at that run's own sample + position. Parameters ---------- @@ -208,12 +257,7 @@ def beam_center_from_center_of_mass(workflow: sciline.Pipeline) -> BeamCenter: : The beam center position as a vector. """ - try: - beam_center = workflow.compute(BeamCenter) - except sciline.UnsatisfiedRequirement: - beam_center = sc.vector([0.0, 0.0, 0.0], unit='m') - workflow = workflow.copy() - workflow[BeamCenter] = beam_center + workflow = _with_default_beam_center(workflow) data = workflow.compute(CorrectedDetector[SampleRun, Numerator]) graph = workflow.compute(ElasticCoordTransformGraph[SampleRun]) @@ -246,29 +290,50 @@ def beam_center_from_center_of_mass(workflow: sciline.Pipeline) -> BeamCenter: pos = pos[select] com = sc.sum(pos * v) / v.sum() - # We compute the shift between the incident beam direction and the center-of-mass - incident_beam = summed.transform_coords('incident_beam', graph=graph).coords[ - 'incident_beam' - ] - n_beam = incident_beam / sc.norm(incident_beam) - com_shift = com - sc.dot(com, n_beam) * n_beam - xy = [com_shift.fields.x.value, com_shift.fields.y.value] - return beam_center + _offsets_to_vector(data=summed, xy=xy, graph=graph) + # The center-of-mass is a point on the beam, so relative to the sample it is the + # beam center, including the distance at which it was determined. + return BeamCenter(com - _sample_position(summed, graph)) -def _offsets_to_vector(data: sc.DataArray, xy: list[float], graph: dict) -> sc.Variable: +@dataclass(frozen=True) +class _BeamPlane: """ - Convert x,y offsets inside the plane normal to the beam to a vector in absolute - coordinates. + Plane normal to the incident beam in which beam-center offsets are expressed. + + ``axis`` points from the sample to the plane, i.e., it carries the distance at which + the beam center is determined. Offsets alone do not define a beam direction, so + searching for a beam center means searching within such a plane. """ - u = data.coords['position'].unit - # Get two vectors that define the plane normal to the beam - coords = data.transform_coords( - ['cyl_x_unit_vector', 'cyl_y_unit_vector'], graph=graph - ).coords - center = xy[0] * coords['cyl_x_unit_vector'] + xy[1] * coords['cyl_y_unit_vector'] - center.unit = u - return center + + unit_x: sc.Variable + unit_y: sc.Variable + axis: sc.Variable + + @staticmethod + def from_beam_center( + data: sc.DataArray, graph: dict, beam_center: sc.Variable + ) -> '_BeamPlane': + """Plane containing ``beam_center``, normal to the incident beam.""" + coords = data.transform_coords( + ['cyl_x_unit_vector', 'cyl_y_unit_vector', 'incident_beam'], graph=graph + ).coords + incident_beam = coords['incident_beam'] + direction = incident_beam / sc.norm(incident_beam) + return _BeamPlane( + unit_x=coords['cyl_x_unit_vector'], + unit_y=coords['cyl_y_unit_vector'], + axis=sc.dot(beam_center, direction) * direction, + ) + + def offsets(self, beam_center: sc.Variable) -> list[float]: + """Components of ``beam_center`` within the plane.""" + return [sc.dot(beam_center, e).value for e in (self.unit_x, self.unit_y)] + + def beam_center(self, xy: list[float]) -> sc.Variable: + """Beam center given by offsets within the plane.""" + center = xy[0] * self.unit_x + xy[1] * self.unit_y + center.unit = self.axis.unit + return center + self.axis def _iofq_in_quadrants( @@ -276,6 +341,7 @@ def _iofq_in_quadrants( workflow: sciline.Pipeline, detector: sc.DataArray, norm: sc.DataArray, + plane: _BeamPlane, ) -> dict[str, sc.DataArray]: """ Compute the intensity as a function of Q inside 4 quadrants in Phi. @@ -288,6 +354,8 @@ def _iofq_in_quadrants( The raw detector. norm: The denominator data for normalization. + plane: + The plane normal to the beam in which ``xy`` are given. Returns ------- @@ -300,9 +368,9 @@ def _iofq_in_quadrants( phi_bins = sc.linspace('phi', -pi, pi, 5, unit='rad') quadrants = ['south-west', 'south-east', 'north-east', 'north-west'] - graph = workflow.compute(ElasticCoordTransformGraph[SampleRun]) workflow = workflow.copy() - workflow[BeamCenter] = _offsets_to_vector(data=detector, xy=xy, graph=graph) + workflow[BeamCenter] = plane.beam_center(xy) + graph = workflow.compute(ElasticCoordTransformGraph[SampleRun]) calibrated = workflow.compute(CorrectedDetector[SampleRun, Numerator]) with_phi = calibrated.transform_coords( 'phi', graph=graph, keep_intermediate=False, keep_inputs=False @@ -321,8 +389,8 @@ def _iofq_in_quadrants( for i, quad in enumerate(quadrants): # Select pixels based on phi sel = (phi >= phi_bins[i]) & (phi < phi_bins[i + 1]) - # The beam center is applied when computing EmptyDetector, set quadrant - # *before* that step. + # Restrict the raw detector to the quadrant, so the denominator (solid angle) + # covers the same pixels as the numerator. workflow[NeXusComponent[snx.NXdetector, SampleRun]] = sc.DataGroup( data=detector[sel] ) @@ -537,8 +605,13 @@ def beam_center_from_iofq( workflow[QBins] = q_bins # Use center of mass to get initial guess for beam center - com_shift = beam_center_from_center_of_mass(workflow) - logger.info('Initial guess for beam center: %s', com_shift) + com = beam_center_from_center_of_mass(workflow) + logger.info('Initial guess for beam center: %s', com) + + # The refinement below only varies the offsets within the plane normal to the beam. + # The distance at which the beam center is determined is taken from the initial + # guess and kept fixed. + plane = _BeamPlane.from_beam_center(data=data, graph=graph, beam_center=com) coords = data.transform_coords( ['cylindrical_x', 'cylindrical_y'], graph=graph @@ -551,14 +624,14 @@ def beam_center_from_iofq( # Refine using Scipy optimize res = minimize( _cost, - x0=[com_shift.fields.x.value, com_shift.fields.y.value], - args=(workflow, detector, norm), + x0=plane.offsets(com), + args=(workflow, detector, norm, plane), bounds=bounds, method=minimizer, tol=tolerance, ) - center = _offsets_to_vector(data=data, xy=res.x, graph=graph) + center = plane.beam_center(res.x) logger.info('Final beam center value: %s', center) logger.info('Beam center finder minimizer info: %s', res) return center diff --git a/packages/esssans/src/ess/sans/common.py b/packages/esssans/src/ess/sans/common.py index e3de665c1..debfe23e6 100644 --- a/packages/esssans/src/ess/sans/common.py +++ b/packages/esssans/src/ess/sans/common.py @@ -6,8 +6,6 @@ import scipp as sc from scipp.constants import g -from .types import BeamCenter, DetectorPositionOffset, RunType - def gravity_vector() -> sc.Variable: """ @@ -91,10 +89,3 @@ def mask_range( else: out.masks[name] = mask_values return out - - -def beam_center_to_detector_position_offset( - beam_center: BeamCenter, -) -> DetectorPositionOffset[RunType]: - """Convert beam center to detector position offset for all runs.""" - return DetectorPositionOffset[RunType](-beam_center) diff --git a/packages/esssans/src/ess/sans/conversions.py b/packages/esssans/src/ess/sans/conversions.py index 44dadbd2c..8fed5262e 100644 --- a/packages/esssans/src/ess/sans/conversions.py +++ b/packages/esssans/src/ess/sans/conversions.py @@ -1,5 +1,7 @@ # SPDX-License-Identifier: BSD-3-Clause # Copyright (c) 2023 Scipp contributors (https://github.com/scipp) +from collections.abc import Callable + import sciline import scipp as sc import scippnexus as snx @@ -13,6 +15,7 @@ from .common import mask_range from .types import ( + BeamCenter, BinnedQ, BinnedQxQy, CorrectForGravity, @@ -33,6 +36,34 @@ ) +def _incident_beam_from_beam_center( + beam_center: sc.Variable, nominal_incident_beam: sc.Variable +) -> Callable[[], sc.Variable]: + """ + Return a provider for ``incident_beam`` that points along the actual beam. + + Scattering angles must be measured relative to the actual beam. The beam passes + through the sample, but its direction, set by the collimation, in general deviates + from the nominal beam axis defined by source and sample. The beam center is a point + on the actual beam, measured from the sample, so it defines that direction. + + Only the direction is taken from the beam center. The length is that of the nominal + incident beam, since ``L1`` is a property of the beamline and not of the beam + center. Scaling ``beam_center`` therefore leaves the result unchanged. + """ + axis = nominal_incident_beam / sc.norm(nominal_incident_beam) + if sc.dot(beam_center, axis).value <= 0.0: + raise ValueError( + f'Invalid beam center {beam_center}. The beam center is the position of ' + 'the beam relative to the sample, so its component along the beam must be ' + 'the (positive) distance from the sample to the plane in which the beam ' + 'center was determined. A beam center given as a transverse offset alone, ' + 'without that distance, cannot define a beam direction.' + ) + incident_beam = sc.norm(nominal_incident_beam) * beam_center / sc.norm(beam_center) + return lambda: incident_beam + + def cyl_unit_vectors(incident_beam: sc.Variable, gravity: sc.Variable): vectors = beam_aligned_unit_vectors(incident_beam=incident_beam, gravity=gravity) return { @@ -102,6 +133,7 @@ def sans_elastic( *, sample_position: Position[snx.NXsample, RunType], source_position: Position[snx.NXsource, RunType], + beam_center: BeamCenter, gravity: GravityVector, ) -> ElasticCoordTransformGraph[RunType]: """ @@ -142,6 +174,10 @@ def sans_elastic( Position of the sample as a vector. source_position: Position of the source as a vector. + beam_center: + Position of the beam center relative to the sample, i.e., the transverse offset + of the beam together with the distance from the sample at which that offset was + determined. Set to a zero vector to apply no correction. """ # noqa: E501 graph = { **beamline.beamline(scatter=True), @@ -150,6 +186,11 @@ def sans_elastic( 'source_position': lambda: source_position, 'gravity': lambda: gravity, } + # A zero beam center means no correction, leaving the plain beamline graph. + if sc.norm(beam_center).value != 0.0: + graph['incident_beam'] = _incident_beam_from_beam_center( + beam_center, nominal_incident_beam=sample_position - source_position + ) if correct_for_gravity: del graph['two_theta'] graph[('two_theta', 'phi')] = scattering_angles_with_gravity diff --git a/packages/esssans/src/ess/sans/parameters.py b/packages/esssans/src/ess/sans/parameters.py index 561e88d57..cfb5f1efe 100644 --- a/packages/esssans/src/ess/sans/parameters.py +++ b/packages/esssans/src/ess/sans/parameters.py @@ -15,7 +15,7 @@ MultiFilenameParameter, ParamWithOptions, StringParameter, - Vector2dParameter, + Vector3dParameter, parameter_registry, ) @@ -112,7 +112,7 @@ parameter_registry[DirectBeamFilename] = FilenameParameter.from_type( DirectBeamFilename, switchable=True ) -parameter_registry[BeamCenter] = Vector2dParameter.from_type( +parameter_registry[BeamCenter] = Vector3dParameter.from_type( BeamCenter, default=sc.vector([0, 0, 0], unit='m') ) diff --git a/packages/esssans/src/ess/sans/types.py b/packages/esssans/src/ess/sans/types.py index 1c7152070..ee5340097 100644 --- a/packages/esssans/src/ess/sans/types.py +++ b/packages/esssans/src/ess/sans/types.py @@ -137,7 +137,28 @@ class TransmissionRun(Generic[ScatteringRunType]): """Filename of direct beam correction""" BeamCenter = NewType('BeamCenter', sc.Variable) -"""Beam center, may be set directly or computed using beam-center finder""" +"""Position of the beam center relative to the sample. + +The beam passes through the sample, but its direction, set by the collimation, in +general deviates from the nominal beam axis defined by source and sample. Together with +the sample position, the beam center defines that direction: it is a point on the actual +beam, given as an offset from the sample. + +Its component along the nominal beam is therefore the distance from the sample at which +the transverse offset was determined, typically the distance of the detector bank used +for the beam-center determination. Recording that distance is essential, since the same +angular deviation of the beam corresponds to a proportionally smaller transverse offset +on banks closer to the sample. + +Only the direction of this vector enters the correction, its magnitude cancels. A value +determined on one detector bank therefore applies to all banks, and, as long as the beam +direction is set by the collimation alone, also to a detector that has been moved to a +different distance from the sample. The value is tied to the collimation that produced +it, so re-determine the beam center whenever the slits change. + +Set to a zero vector to apply no beam-center correction. May be set directly or computed +using a beam-center finder, see :py:mod:`ess.sans.beam_center_finder`. +""" WavelengthMask = NewType('WavelengthMask', sc.DataArray | None) """Optional wavelength mask""" diff --git a/packages/esssans/src/ess/sans/workflow.py b/packages/esssans/src/ess/sans/workflow.py index 16303599b..fa854f8ad 100644 --- a/packages/esssans/src/ess/sans/workflow.py +++ b/packages/esssans/src/ess/sans/workflow.py @@ -9,7 +9,7 @@ from ess.reduce.parameter import parameter_mappers from ess.reduce.unwrap import GenericUnwrapWorkflow, WavelengthLutMode -from . import common, conversions, i_of_q, masking, normalization +from . import conversions, i_of_q, masking, normalization from .types import ( BackgroundRun, CorrectForGravity, @@ -155,7 +155,6 @@ def with_background_runs( *i_of_q.providers, *masking.providers, *normalization.providers, - common.beam_center_to_detector_position_offset, ) """ List of providers for setting up a Sciline pipeline. diff --git a/packages/esssans/tests/beam_center_finder_test.py b/packages/esssans/tests/beam_center_finder_test.py new file mode 100644 index 000000000..a664f09c6 --- /dev/null +++ b/packages/esssans/tests/beam_center_finder_test.py @@ -0,0 +1,154 @@ +# SPDX-License-Identifier: BSD-3-Clause +# Copyright (c) 2023 Scipp contributors (https://github.com/scipp) + +import sciline +import scipp as sc +import scippnexus as snx +from ess.sans.beam_center_finder import beam_center_from_center_of_mass +from ess.sans.common import gravity_vector +from ess.sans.conversions import ElasticCoordTransformGraph, sans_elastic +from ess.sans.types import ( + BeamCenter, + CorrectedDetector, + CorrectForGravity, + Numerator, + Position, + SampleRun, +) + + +def make_detector( + *, + center: sc.Variable, + distance: sc.Variable, + size: int = 41, + extent: float = 1.0, +) -> sc.DataArray: + """ + Flat square detector with a Gaussian intensity distribution centered on ``center``. + + The detector is normal to the Z axis. The pattern is narrow enough to be contained + well within the detector, so that its center of mass is unbiased. + """ + x = sc.linspace('x', -extent, extent, size, unit='m') + y = sc.linspace('y', -extent, extent, size, unit='m') + position = sc.spatial.as_vectors( + x.broadcast(sizes={'y': size, 'x': size}), + y.broadcast(sizes={'y': size, 'x': size}), + distance.broadcast(sizes={'y': size, 'x': size}), + ) + r2 = (position.fields.x - center.fields.x) ** 2 + ( + position.fields.y - center.fields.y + ) ** 2 + sigma = sc.scalar(0.2, unit='m') + counts = sc.exp(-r2 / (2 * sigma**2)) + # A wavelength dim exercises the summing of non-position dims. + counts = counts.broadcast(sizes={**counts.sizes, 'wavelength': 3}).copy() + return sc.DataArray( + data=counts, + coords={ + 'position': position, + 'wavelength': sc.linspace('wavelength', 1.0, 3.0, 3, unit='angstrom'), + }, + ) + + +def make_workflow( + *, + detector: sc.DataArray, + sample_position: sc.Variable, + source_position: sc.Variable | None = None, +) -> sciline.Pipeline: + if source_position is None: + source_position = sample_position - sc.vector([0, 0, 20.0], unit='m') + workflow = sciline.Pipeline(()) + workflow[CorrectedDetector[SampleRun, Numerator]] = detector + workflow[ElasticCoordTransformGraph[SampleRun]] = sans_elastic( + CorrectForGravity(False), + sample_position=Position[snx.NXsample, SampleRun](sample_position), + source_position=Position[snx.NXsource, SampleRun](source_position), + beam_center=BeamCenter(sc.vector([0, 0, 0], unit='m')), + gravity=gravity_vector(), + ) + return workflow + + +DISTANCE = sc.scalar(5.0, unit='m') + + +def test_center_of_mass_finds_center_of_symmetric_pattern() -> None: + center = sc.vector([0.1, -0.07, 0.0], unit='m') + workflow = make_workflow( + detector=make_detector(center=center, distance=DISTANCE), + sample_position=sc.vector([0, 0, 0], unit='m'), + ) + result = beam_center_from_center_of_mass(workflow) + assert sc.allclose( + result, center + DISTANCE * sc.vector([0, 0, 1]), atol=sc.scalar(1e-3, unit='m') + ) + + +def test_center_of_mass_records_the_distance_at_which_it_was_determined() -> None: + # The beam center only defines a beam direction together with the distance from the + # sample at which it was determined, so the finder must report that distance. + center = sc.vector([0.1, -0.07, 0.0], unit='m') + near = beam_center_from_center_of_mass( + make_workflow( + detector=make_detector(center=center, distance=sc.scalar(2.0, unit='m')), + sample_position=sc.vector([0, 0, 0], unit='m'), + ) + ) + far = beam_center_from_center_of_mass( + make_workflow( + detector=make_detector(center=center, distance=sc.scalar(8.0, unit='m')), + sample_position=sc.vector([0, 0, 0], unit='m'), + ) + ) + tol = sc.scalar(1e-3, unit='m') + assert sc.isclose(near.fields.z, sc.scalar(2.0, unit='m'), atol=tol) + assert sc.isclose(far.fields.z, sc.scalar(8.0, unit='m'), atol=tol) + + +def test_center_of_mass_is_measured_relative_to_the_sample_position() -> None: + # Pattern is centered on the beam axis through the sample, so there is no transverse + # offset of the beam center, despite the sample being off-axis. + sample_position = sc.vector([0.1, -0.07, 0.0], unit='m') + workflow = make_workflow( + detector=make_detector(center=sample_position, distance=DISTANCE), + sample_position=sample_position, + ) + result = beam_center_from_center_of_mass(workflow) + assert sc.allclose( + result, + DISTANCE * sc.vector([0, 0, 1]), + atol=sc.scalar(1e-3, unit='m'), + ) + + +def test_center_of_mass_result_yields_zero_scattering_angle_at_the_pattern_center() -> ( + None +): + center = sc.vector([0.1, -0.07, 0.0], unit='m') + sample_position = sc.vector([0, 0, 0], unit='m') + workflow = make_workflow( + detector=make_detector(center=center, distance=DISTANCE), + sample_position=sample_position, + ) + beam_center = beam_center_from_center_of_mass(workflow) + graph = sans_elastic( + CorrectForGravity(False), + sample_position=Position[snx.NXsample, SampleRun](sample_position), + source_position=Position[snx.NXsource, SampleRun]( + sample_position - sc.vector([0, 0, 20.0], unit='m') + ), + beam_center=BeamCenter(beam_center), + gravity=gravity_vector(), + ) + da = sc.DataArray( + data=sc.scalar(1.0), + coords={'position': center + DISTANCE * sc.vector([0, 0, 1])}, + ) + two_theta = da.transform_coords('two_theta', graph=graph).coords['two_theta'] + assert sc.isclose( + two_theta, sc.scalar(0.0, unit='rad'), atol=sc.scalar(1e-4, unit='rad') + ) diff --git a/packages/esssans/tests/beam_center_test.py b/packages/esssans/tests/beam_center_test.py new file mode 100644 index 000000000..5eef241aa --- /dev/null +++ b/packages/esssans/tests/beam_center_test.py @@ -0,0 +1,131 @@ +# SPDX-License-Identifier: BSD-3-Clause +# Copyright (c) 2023 Scipp contributors (https://github.com/scipp) +"""Tests for how the beam center is applied, as opposed to how it is determined.""" + +import pytest +import scipp as sc +import scippnexus as snx +from ess.sans.common import gravity_vector +from ess.sans.conversions import sans_elastic +from ess.sans.types import BeamCenter, CorrectForGravity, Position, SampleRun + +SAMPLE_POSITION = sc.vector([0.0, 0.0, 0.0], unit='m') +# Offset of the source from the sample, i.e., the nominal beam is along Z. +SOURCE_OFFSET = sc.vector([0.0, 0.0, -20.0], unit='m') +NO_BEAM_CENTER = sc.vector([0.0, 0.0, 0.0], unit='m') +# Beam center determined on a detector 5 m downstream of the sample. +BEAM_CENTER = sc.vector([0.05, -0.03, 5.0], unit='m') + + +def make_graph( + *, + beam_center: sc.Variable, + correct_for_gravity: bool = False, + sample_position: sc.Variable = SAMPLE_POSITION, +) -> dict: + return sans_elastic( + CorrectForGravity(correct_for_gravity), + sample_position=Position[snx.NXsample, SampleRun](sample_position), + source_position=Position[snx.NXsource, SampleRun]( + sample_position + SOURCE_OFFSET + ), + beam_center=BeamCenter(beam_center), + gravity=gravity_vector(), + ) + + +def angle(position: sc.Variable, graph: dict, name: str = 'two_theta') -> sc.Variable: + # A negligible wavelength keeps the gravity drop out of the comparison. + da = sc.DataArray( + data=sc.scalar(1.0), + coords={'position': position, 'wavelength': sc.scalar(1e-6, unit='angstrom')}, + ) + return da.transform_coords(name, graph=graph).coords[name] + + +def assert_zero_angle(value: sc.Variable) -> None: + assert sc.isclose( + value, sc.scalar(0.0, unit='rad'), atol=sc.scalar(1e-9, unit='rad') + ) + + +def test_two_theta_is_zero_on_beam_axis_without_beam_center() -> None: + position = sc.vector([0.0, 0.0, 5.0], unit='m') + assert_zero_angle(angle(position, make_graph(beam_center=NO_BEAM_CENTER))) + + +def test_two_theta_is_zero_at_beam_center() -> None: + position = SAMPLE_POSITION + BEAM_CENTER + assert_zero_angle(angle(position, make_graph(beam_center=BEAM_CENTER))) + + +@pytest.mark.parametrize('fraction', [0.2, 0.5, 0.9, 2.0]) +def test_two_theta_is_zero_along_the_beam_at_any_distance(fraction: float) -> None: + # The beam center measured at one distance defines the beam direction, so pixels + # on the same ray but at a different distance also see zero scattering angle. + position = SAMPLE_POSITION + fraction * BEAM_CENTER + assert_zero_angle(angle(position, make_graph(beam_center=BEAM_CENTER))) + + +def test_beam_center_shifts_two_theta() -> None: + position = SAMPLE_POSITION + BEAM_CENTER + without = angle(position, make_graph(beam_center=NO_BEAM_CENTER)) + with_center = angle(position, make_graph(beam_center=BEAM_CENTER)) + assert without.value > 0.0 + assert with_center.value < without.value + + +def _transverse_offset(direction: sc.Variable) -> sc.Variable: + """Offset of 20 cm along ``direction``.""" + return sc.scalar(0.2, unit='m') * direction / sc.norm(direction) + + +def test_phi_is_zero_horizontally_beside_the_beam_center() -> None: + # Perpendicular to the beam and to gravity, i.e., phi == 0. + right = sc.cross(-gravity_vector(), BEAM_CENTER) + position = SAMPLE_POSITION + BEAM_CENTER + _transverse_offset(right) + assert_zero_angle(angle(position, make_graph(beam_center=BEAM_CENTER), name='phi')) + + +def test_phi_is_ninety_degrees_above_the_beam_center() -> None: + # Perpendicular to the beam, in the plane spanned by the beam and gravity. + right = sc.cross(-gravity_vector(), BEAM_CENTER) + position = ( + SAMPLE_POSITION + BEAM_CENTER + _transverse_offset(sc.cross(BEAM_CENTER, right)) + ) + phi = angle(position, make_graph(beam_center=BEAM_CENTER), name='phi') + assert sc.isclose( + phi, + 0.5 * sc.constants.pi.value * sc.Unit('rad'), + atol=sc.scalar(1e-9, unit='rad'), + ) + + +@pytest.mark.parametrize('correct_for_gravity', [False, True]) +def test_beam_center_is_applied_relative_to_sample_position( + correct_for_gravity: bool, +) -> None: + sample_position = sc.vector([0.1, -0.07, 0.0], unit='m') + position = sample_position + BEAM_CENTER + graph = make_graph( + beam_center=BEAM_CENTER, + correct_for_gravity=correct_for_gravity, + sample_position=sample_position, + ) + assert_zero_angle(angle(position, graph)) + + +@pytest.mark.parametrize( + 'beam_center', + [ + # Offset without a distance, i.e., the pre-existing 2-D convention. + sc.vector([0.05, -0.03, 0.0], unit='m'), + # Upstream of the sample. + sc.vector([0.05, -0.03, -5.0], unit='m'), + ], +) +def test_beam_center_without_a_positive_distance_raises( + beam_center: sc.Variable, +) -> None: + with pytest.raises(ValueError, match='distance'): + make_graph(beam_center=beam_center) diff --git a/packages/esssans/tests/isissans/sans2d_reduction_test.py b/packages/esssans/tests/isissans/sans2d_reduction_test.py index e32af35cc..8229bf1f6 100644 --- a/packages/esssans/tests/isissans/sans2d_reduction_test.py +++ b/packages/esssans/tests/isissans/sans2d_reduction_test.py @@ -197,7 +197,9 @@ def test_pixel_dependent_direct_beam_is_supported(pipeline, uncertainties): assert result.dims == ('Q',) -MANTID_BEAM_CENTER = sc.vector([0.09288, -0.08195, 0], unit='m') +# Transverse offset obtained from Mantid, at the distance of the SANS2D detector, which +# is where it was determined. +MANTID_BEAM_CENTER = sc.vector([0.09288, -0.08195, 4.0], unit='m') def test_beam_center_from_center_of_mass_is_close_to_verified_result(pipeline): @@ -209,7 +211,7 @@ def test_beam_center_from_center_of_mass_is_close_to_verified_result(pipeline): def test_beam_center_from_center_of_mass_independent_of_set_beam_center(pipeline): - pipeline[BeamCenter] = sc.vector([0.1, -0.1, 0], unit='m') + pipeline[BeamCenter] = sc.vector([0.1, -0.1, 3.0], unit='m') center = sans.beam_center_from_center_of_mass(pipeline) assert sc.allclose(center, MANTID_BEAM_CENTER, atol=sc.scalar(3e-3, unit='m')) @@ -262,7 +264,7 @@ def test_beam_center_finder_works_with_direct_beam(pipeline): reason="Beam center finder tests are slow and we are revisiting the approach." ) def test_beam_center_finder_independent_of_set_beam_center(pipeline): - pipeline[BeamCenter] = sc.vector([0.1, -0.1, 0], unit='m') + pipeline[BeamCenter] = sc.vector([0.1, -0.1, 3.0], unit='m') q_bins = sc.linspace('Q', 0.02, 0.3, 71, unit='1/angstrom') center_with_direct_beam = sans.beam_center_finder.beam_center_from_iofq( workflow=pipeline, q_bins=q_bins diff --git a/packages/esssans/tests/loki/iofq_test.py b/packages/esssans/tests/loki/iofq_test.py index 5e1f6d0ad..4d445b867 100644 --- a/packages/esssans/tests/loki/iofq_test.py +++ b/packages/esssans/tests/loki/iofq_test.py @@ -111,9 +111,11 @@ def test_pipeline_can_compute_IofQ_in_event_mode( sc.values(reference.data), sc.values(result.hist().data), # Could be 1e-11, but currently the workflow defaults to float32 data, as - # returned by ScippNexus. + # returned by ScippNexus. The background subtraction cancels most of the + # signal, so a single ulp of the unsubtracted intensity shows up here as a + # much larger deviation. rtol=sc.scalar(1e-7), - atol=sc.scalar(1e-7), + atol=sc.scalar(1e-6), ) if uncertainties == UncertaintyBroadcastMode.drop: # Could both be 1e-14 if using float64 @@ -280,7 +282,7 @@ def test_beam_center_from_center_of_mass_is_close_to_verified_result(larmor_work pipeline, loki.data.loki_tutorial_mask_filenames() ) center = sans.beam_center_from_center_of_mass(pipeline) - reference = sc.vector([-0.0291487, -0.0181614, 0], unit='m') + reference = sc.vector([-0.0291487, -0.0181614, 4.1161], unit='m') assert sc.allclose(center, reference) @@ -294,7 +296,7 @@ def test_beam_center_from_center_of_mass_alternative_is_close_to_verified_result center = sans.beam_center_finder.beam_center_from_center_of_mass_alternative( pipeline ) - reference = sc.vector([-0.0248743, -0.0166965, 0], unit='m') + reference = sc.vector([-0.0248743, -0.0166965, 4.12268], unit='m') assert sc.allclose(center, reference) diff --git a/packages/esssans/tests/loki/reference_IofQ_UncertaintyBroadcastMode.drop.hdf5 b/packages/esssans/tests/loki/reference_IofQ_UncertaintyBroadcastMode.drop.hdf5 index 6a6c950bd..0ca6d0e1c 100644 Binary files a/packages/esssans/tests/loki/reference_IofQ_UncertaintyBroadcastMode.drop.hdf5 and b/packages/esssans/tests/loki/reference_IofQ_UncertaintyBroadcastMode.drop.hdf5 differ diff --git a/packages/esssans/tests/loki/reference_IofQxy_UncertaintyBroadcastMode.drop.hdf5 b/packages/esssans/tests/loki/reference_IofQxy_UncertaintyBroadcastMode.drop.hdf5 index 28855e167..30b3af3d2 100644 Binary files a/packages/esssans/tests/loki/reference_IofQxy_UncertaintyBroadcastMode.drop.hdf5 and b/packages/esssans/tests/loki/reference_IofQxy_UncertaintyBroadcastMode.drop.hdf5 differ diff --git a/packages/esssans/tests/loki/workflow_test.py b/packages/esssans/tests/loki/workflow_test.py index 0d20b94e1..e754620da 100644 --- a/packages/esssans/tests/loki/workflow_test.py +++ b/packages/esssans/tests/loki/workflow_test.py @@ -10,6 +10,7 @@ BackgroundRun, BackgroundSubtractedIofQ, BeamCenter, + EmptyDetector, Filename, IntensityQ, LookupTableFilename, @@ -18,6 +19,7 @@ QBins, ReturnEvents, SampleRun, + SolidAngle, UncertaintyBroadcastMode, WavelengthDetector, ) @@ -57,6 +59,23 @@ def test_loki_larmor_workflow_parameters_with_param_returns_param(): assert parameters.keys() == {ReturnEvents} +def test_beam_center_does_not_move_detector_positions(larmor_workflow): + wf = larmor_workflow() + wf[BeamCenter] = sc.vector([0.0, 0.0, 0.0], unit='m') + reference = wf.compute(EmptyDetector[SampleRun]).coords['position'] + wf[BeamCenter] = sc.vector([0.03, -0.02, 0.0], unit='m') + position = wf.compute(EmptyDetector[SampleRun]).coords['position'] + assert sc.identical(position, reference) + + +def test_solid_angle_is_independent_of_beam_center(larmor_workflow): + wf = larmor_workflow() + wf[BeamCenter] = sc.vector([0.0, 0.0, 0.0], unit='m') + reference = wf.compute(SolidAngle[SampleRun]) + wf[BeamCenter] = sc.vector([0.03, -0.02, 0.0], unit='m') + assert sc.identical(wf.compute(SolidAngle[SampleRun]), reference) + + def test_loki_larmor_workflow_compute_with_single_pixel_mask(larmor_workflow): wf = larmor_workflow(no_masks=False) wf[UncertaintyBroadcastMode] = UncertaintyBroadcastMode.drop