From b3bf0c505769e0a729a7bb17b0ec08fa1f8939f6 Mon Sep 17 00:00:00 2001 From: Simon Heybrock Date: Thu, 7 Nov 2024 15:18:42 +0100 Subject: [PATCH 01/11] Draft of Q-resolution computation --- src/ess/sans/conversions.py | 29 ++++++++++++++++++++ src/ess/sans/qresolution.py | 54 +++++++++++++++++++++++++++++++++++++ src/ess/sans/types.py | 8 +++--- src/ess/sans/workflow.py | 3 ++- 4 files changed, 90 insertions(+), 4 deletions(-) create mode 100644 src/ess/sans/qresolution.py diff --git a/src/ess/sans/conversions.py b/src/ess/sans/conversions.py index ce157dd9..5cd5278c 100644 --- a/src/ess/sans/conversions.py +++ b/src/ess/sans/conversions.py @@ -3,6 +3,7 @@ from typing import NewType import scipp as sc +from scipp.constants import pi from scippneutron.conversion.beamline import ( beam_aligned_unit_vectors, scattering_angles_with_gravity, @@ -25,6 +26,7 @@ MonitorTerm, MonitorType, Numerator, + Resolution, RunType, ScatteringRunType, TofMonitor, @@ -171,6 +173,10 @@ def detector_to_wavelength( ) +ModeratorTimeSpread = NewType("ModeratorTimeSpread", sc.DataArray) +"""Moderator time-spread as a function of wavelength.""" + + def mask_wavelength_q( da: CleanSummedQ[ScatteringRunType, Numerator], mask: WavelengthMask ) -> WavelengthScaledQ[ScatteringRunType, Numerator]: @@ -179,6 +185,28 @@ def mask_wavelength_q( return WavelengthScaledQ[ScatteringRunType, Numerator](da) +def mask_and_compute_resolution_q( + pixel_term: CleanSummedQ[ScatteringRunType, Resolution], + mask: WavelengthMask, + moderator_time_spread: ModeratorTimeSpread, +) -> WavelengthScaledQ[ScatteringRunType, Resolution]: + """ + Compute the masked Q-resolution in (Q, lambda) space. + + CleanSummedQ has been summed over pixels but not over wavelengths. This is exactly + what is required for performing the remaining scaling and addition of the moderator + term to obtain the Q-resolution. The result is still in (Q, lambda) space. + """ + if mask is not None: + pixel_term = mask_range(pixel_term, mask=mask) + lambda2 = pixel_term.coords['wavelength'] ** 2 + Q2 = pixel_term.coords['Q'] ** 2 + resolution = ( + pi**2 / (3 * lambda2**2) * pixel_term + Q2 * moderator_time_spread**2 / lambda2 + ) + return WavelengthScaledQ[ScatteringRunType, Resolution](resolution) + + def mask_wavelength_qxy( da: CleanSummedQxy[ScatteringRunType, Numerator], mask: WavelengthMask ) -> WavelengthScaledQxy[ScatteringRunType, Numerator]: @@ -256,6 +284,7 @@ def compute_Qxy( detector_to_wavelength, mask_wavelength_q, mask_wavelength_qxy, + mask_and_compute_resolution_q, mask_and_scale_wavelength_q, mask_and_scale_wavelength_qxy, compute_Q, diff --git a/src/ess/sans/qresolution.py b/src/ess/sans/qresolution.py new file mode 100644 index 00000000..ab71e16f --- /dev/null +++ b/src/ess/sans/qresolution.py @@ -0,0 +1,54 @@ +# SPDX-License-Identifier: BSD-3-Clause +# Copyright (c) 2024 Scipp contributors (https://github.com/scipp) +"""Q-resolution calculation for SANS data.""" + +from typing import NewType + +import scipp as sc + +from .types import CleanQ, Denominator, Resolution, SampleRun + +DeltaR = NewType("DeltaR", sc.Variable) +"""Virtual ring width on the detector.""" + +SampleApertureRadius = NewType("SampleApertureRadius", sc.Variable) +"""Sample aperture radius, R2.""" +SourceApertureRadius = NewType("SourceApertureRadius", sc.Variable) +"""Source aperture radius, R1.""" +SigmaModerator = NewType("SigmaModerator", sc.DataArray) +"""Moderator time spread as a function of wavelength.""" +CollimationLength = NewType("CollimationLength", sc.Variable) +"""Collimation length.""" + +ModeratorTimeSpread = NewType("ModeratorTimeSpread", sc.DataArray) +"""Moderator time-spread as a function of wavelength.""" + +QResolutionDetectorTerm = NewType("QResolutionDetectorTerm", sc.DataArray) + + +def pixel_term( + detector: CleanQ[SampleRun, Denominator], + delta_r: DeltaR, + sample_aperture: SampleApertureRadius, + source_aperture: SourceApertureRadius, + collimation_length: CollimationLength, +) -> CleanQ[SampleRun, Resolution]: + """ + Calculate the pixel term for Q-resolution. + + We compute this based on CleanQ[SampleRun, Denominator]. This ensures that + + 1. We get the correct Q-value, based on wavelength binning used elsewhere. + 2. Masks are included. + 3. We do not depend on neutron data, by using the denominator instead of the + numerator. + """ + L2 = detector.coords["L2"] + L3 = sc.reciprocal(sc.reciprocal(collimation_length) + sc.reciprocal(L2)) + # + result = detector.copy(deep=False) + result.data = ( + 3 * ((source_aperture / collimation_length) ** 2 + (sample_aperture / L3) ** 2) + + (delta_r / L2) ** 2 + ) + return CleanQ[SampleRun, Resolution](result) diff --git a/src/ess/sans/types.py b/src/ess/sans/types.py index a5e21bca..4c9a9488 100644 --- a/src/ess/sans/types.py +++ b/src/ess/sans/types.py @@ -40,11 +40,13 @@ UncertaintyBroadcastMode = _UncertaintyBroadcastMode # 1.3 Numerator and denominator of IofQ -Numerator = NewType('Numerator', sc.DataArray) +Numerator = NewType('Numerator', int) """Numerator of IofQ""" -Denominator = NewType('Denominator', sc.DataArray) +Denominator = NewType('Denominator', int) """Denominator of IofQ""" -IofQPart = TypeVar('IofQPart', Numerator, Denominator) +Resolution = NewType('Resolution', int) +"""Q-Resolution""" +IofQPart = TypeVar('IofQPart', Numerator, Denominator, Resolution) """TypeVar used for specifying Numerator or Denominator of IofQ""" # 1.4 Entry paths in NeXus files diff --git a/src/ess/sans/workflow.py b/src/ess/sans/workflow.py index 15a1a567..f03b9240 100644 --- a/src/ess/sans/workflow.py +++ b/src/ess/sans/workflow.py @@ -9,7 +9,7 @@ from ess.reduce.nexus.workflow import GenericNeXusWorkflow from ess.reduce.parameter import parameter_mappers -from . import common, conversions, i_of_q, masking, normalization +from . import common, conversions, i_of_q, masking, normalization, qresolution from .types import ( BackgroundRun, CleanSummedQ, @@ -151,6 +151,7 @@ def with_background_runs( *masking.providers, *normalization.providers, common.beam_center_to_detector_position_offset, + qresolution.pixel_term, ) """ List of providers for setting up a Sciline pipeline. From 3865e0a2bd43398025886c66ee996f78cadc6a2f Mon Sep 17 00:00:00 2001 From: Simon Heybrock Date: Thu, 7 Nov 2024 15:52:50 +0100 Subject: [PATCH 02/11] WIP --- src/ess/sans/conversions.py | 29 -------------- src/ess/sans/qresolution.py | 78 ++++++++++++++++++++++++++++++++++--- src/ess/sans/types.py | 4 +- src/ess/sans/workflow.py | 2 +- 4 files changed, 74 insertions(+), 39 deletions(-) diff --git a/src/ess/sans/conversions.py b/src/ess/sans/conversions.py index 5cd5278c..ce157dd9 100644 --- a/src/ess/sans/conversions.py +++ b/src/ess/sans/conversions.py @@ -3,7 +3,6 @@ from typing import NewType import scipp as sc -from scipp.constants import pi from scippneutron.conversion.beamline import ( beam_aligned_unit_vectors, scattering_angles_with_gravity, @@ -26,7 +25,6 @@ MonitorTerm, MonitorType, Numerator, - Resolution, RunType, ScatteringRunType, TofMonitor, @@ -173,10 +171,6 @@ def detector_to_wavelength( ) -ModeratorTimeSpread = NewType("ModeratorTimeSpread", sc.DataArray) -"""Moderator time-spread as a function of wavelength.""" - - def mask_wavelength_q( da: CleanSummedQ[ScatteringRunType, Numerator], mask: WavelengthMask ) -> WavelengthScaledQ[ScatteringRunType, Numerator]: @@ -185,28 +179,6 @@ def mask_wavelength_q( return WavelengthScaledQ[ScatteringRunType, Numerator](da) -def mask_and_compute_resolution_q( - pixel_term: CleanSummedQ[ScatteringRunType, Resolution], - mask: WavelengthMask, - moderator_time_spread: ModeratorTimeSpread, -) -> WavelengthScaledQ[ScatteringRunType, Resolution]: - """ - Compute the masked Q-resolution in (Q, lambda) space. - - CleanSummedQ has been summed over pixels but not over wavelengths. This is exactly - what is required for performing the remaining scaling and addition of the moderator - term to obtain the Q-resolution. The result is still in (Q, lambda) space. - """ - if mask is not None: - pixel_term = mask_range(pixel_term, mask=mask) - lambda2 = pixel_term.coords['wavelength'] ** 2 - Q2 = pixel_term.coords['Q'] ** 2 - resolution = ( - pi**2 / (3 * lambda2**2) * pixel_term + Q2 * moderator_time_spread**2 / lambda2 - ) - return WavelengthScaledQ[ScatteringRunType, Resolution](resolution) - - def mask_wavelength_qxy( da: CleanSummedQxy[ScatteringRunType, Numerator], mask: WavelengthMask ) -> WavelengthScaledQxy[ScatteringRunType, Numerator]: @@ -284,7 +256,6 @@ def compute_Qxy( detector_to_wavelength, mask_wavelength_q, mask_wavelength_qxy, - mask_and_compute_resolution_q, mask_and_scale_wavelength_q, mask_and_scale_wavelength_qxy, compute_Q, diff --git a/src/ess/sans/qresolution.py b/src/ess/sans/qresolution.py index ab71e16f..93866d58 100644 --- a/src/ess/sans/qresolution.py +++ b/src/ess/sans/qresolution.py @@ -5,8 +5,18 @@ from typing import NewType import scipp as sc +from scipp.constants import pi -from .types import CleanQ, Denominator, Resolution, SampleRun +from .common import mask_range +from .types import ( + CleanQ, + Denominator, + DimsToKeep, + QBins, + SampleRun, + WavelengthMask, + WavelengthScaledQ, +) DeltaR = NewType("DeltaR", sc.Variable) """Virtual ring width on the detector.""" @@ -19,11 +29,14 @@ """Moderator time spread as a function of wavelength.""" CollimationLength = NewType("CollimationLength", sc.Variable) """Collimation length.""" - ModeratorTimeSpread = NewType("ModeratorTimeSpread", sc.DataArray) """Moderator time-spread as a function of wavelength.""" -QResolutionDetectorTerm = NewType("QResolutionDetectorTerm", sc.DataArray) + +QResolutionPixelTerm = NewType("QResolutionPixelTerm", sc.DataArray) +QResolutionPixelTermGroupedQ = NewType("QResolutionPixelTermGroupedQ", sc.DataArray) +QResolutionByWavelength = NewType("QResolutionByWavelength", sc.DataArray) +QResolution = NewType("QResolution", sc.DataArray) def pixel_term( @@ -32,7 +45,7 @@ def pixel_term( sample_aperture: SampleApertureRadius, source_aperture: SourceApertureRadius, collimation_length: CollimationLength, -) -> CleanQ[SampleRun, Resolution]: +) -> QResolutionPixelTerm: """ Calculate the pixel term for Q-resolution. @@ -45,10 +58,63 @@ def pixel_term( """ L2 = detector.coords["L2"] L3 = sc.reciprocal(sc.reciprocal(collimation_length) + sc.reciprocal(L2)) - # result = detector.copy(deep=False) result.data = ( 3 * ((source_aperture / collimation_length) ** 2 + (sample_aperture / L3) ** 2) + (delta_r / L2) ** 2 ) - return CleanQ[SampleRun, Resolution](result) + # TODO + # Give different name, as we do not actually feed into bin_in_q + return QResolutionPixelTerm(result) + + +# TODO What is the point of naming the input CleanQ and output +# CleanSummedQ? It is not sharing functions, so use a different name +def groupby_q_max( + data: QResolutionPixelTerm, + q_bins: QBins, + dims_to_keep: DimsToKeep, +) -> QResolutionPixelTermGroupedQ: + out = data.groupby('Q', bins=q_bins).max('detector_number') + return QResolutionPixelTermGroupedQ(out) + + +def mask_and_compute_resolution_q( + pixel_term: QResolutionPixelTermGroupedQ, + mask: WavelengthMask, + moderator_time_spread: ModeratorTimeSpread, +) -> QResolutionByWavelength: + """ + Compute the masked Q-resolution in (Q, lambda) space. + + CleanSummedQ has been summed over pixels but not over wavelengths. This is exactly + what is required for performing the remaining scaling and addition of the moderator + term to obtain the Q-resolution. The result is still in (Q, lambda) space. + """ + if mask is not None: + pixel_term = mask_range(pixel_term, mask=mask) + lambda2 = pixel_term.coords['wavelength'] ** 2 + Q2 = pixel_term.coords['Q'] ** 2 + resolution = ( + pi**2 / (3 * lambda2**2) * pixel_term + Q2 * moderator_time_spread**2 / lambda2 + ) + return QResolutionByWavelength(resolution) + + +def reduce_resolution_q( + data: QResolutionByWavelength, + bands: ProcessedWavelengthBands, +) -> QResolution: + # TODO + # We do not want to use reduce_q, but use a max again! + # but reduce_q is quite complex, can we reuse it but use binned data + # so it does not sum? or call the underlying helper function! + pass + + +providers = ( + pixel_term, + groupby_q_max, + mask_and_compute_resolution_q, + reduce_resolution_q, +) diff --git a/src/ess/sans/types.py b/src/ess/sans/types.py index 4c9a9488..f5ff061c 100644 --- a/src/ess/sans/types.py +++ b/src/ess/sans/types.py @@ -44,9 +44,7 @@ """Numerator of IofQ""" Denominator = NewType('Denominator', int) """Denominator of IofQ""" -Resolution = NewType('Resolution', int) -"""Q-Resolution""" -IofQPart = TypeVar('IofQPart', Numerator, Denominator, Resolution) +IofQPart = TypeVar('IofQPart', Numerator, Denominator) """TypeVar used for specifying Numerator or Denominator of IofQ""" # 1.4 Entry paths in NeXus files diff --git a/src/ess/sans/workflow.py b/src/ess/sans/workflow.py index f03b9240..287d097a 100644 --- a/src/ess/sans/workflow.py +++ b/src/ess/sans/workflow.py @@ -151,7 +151,7 @@ def with_background_runs( *masking.providers, *normalization.providers, common.beam_center_to_detector_position_offset, - qresolution.pixel_term, + *qresolution.providers, ) """ List of providers for setting up a Sciline pipeline. From 09e130f7791c08fcd16976d3217f6d6aeebeba34 Mon Sep 17 00:00:00 2001 From: Simon Heybrock Date: Thu, 7 Nov 2024 15:56:31 +0100 Subject: [PATCH 03/11] Outline remaining work --- src/ess/sans/qresolution.py | 15 +++++++-------- 1 file changed, 7 insertions(+), 8 deletions(-) diff --git a/src/ess/sans/qresolution.py b/src/ess/sans/qresolution.py index 93866d58..2622fd3c 100644 --- a/src/ess/sans/qresolution.py +++ b/src/ess/sans/qresolution.py @@ -8,14 +8,15 @@ from scipp.constants import pi from .common import mask_range +from .normalization import _reduce from .types import ( CleanQ, Denominator, DimsToKeep, + ProcessedWavelengthBands, QBins, SampleRun, WavelengthMask, - WavelengthScaledQ, ) DeltaR = NewType("DeltaR", sc.Variable) @@ -75,6 +76,8 @@ def groupby_q_max( q_bins: QBins, dims_to_keep: DimsToKeep, ) -> QResolutionPixelTermGroupedQ: + # TODO Handle multi dim and dims_to_keep! + # Can we use common helper function from bin_in_q? out = data.groupby('Q', bins=q_bins).max('detector_number') return QResolutionPixelTermGroupedQ(out) @@ -102,14 +105,10 @@ def mask_and_compute_resolution_q( def reduce_resolution_q( - data: QResolutionByWavelength, - bands: ProcessedWavelengthBands, + data: QResolutionByWavelength, bands: ProcessedWavelengthBands ) -> QResolution: - # TODO - # We do not want to use reduce_q, but use a max again! - # but reduce_q is quite complex, can we reuse it but use binned data - # so it does not sum? or call the underlying helper function! - pass + # TODO Add op argument to allow injection of different reduction functions + return QResolution(_reduce(data, bands, op=sc.max)) providers = ( From 0c4056ba047642a9fcdc5feca1f002ef78e9e8b3 Mon Sep 17 00:00:00 2001 From: Simon Heybrock Date: Thu, 7 Nov 2024 16:05:17 +0100 Subject: [PATCH 04/11] Progress --- src/ess/sans/qresolution.py | 17 ++++++++++------- 1 file changed, 10 insertions(+), 7 deletions(-) diff --git a/src/ess/sans/qresolution.py b/src/ess/sans/qresolution.py index 2622fd3c..c469d672 100644 --- a/src/ess/sans/qresolution.py +++ b/src/ess/sans/qresolution.py @@ -2,6 +2,7 @@ # Copyright (c) 2024 Scipp contributors (https://github.com/scipp) """Q-resolution calculation for SANS data.""" +import uuid from typing import NewType import scipp as sc @@ -72,14 +73,16 @@ def pixel_term( # TODO What is the point of naming the input CleanQ and output # CleanSummedQ? It is not sharing functions, so use a different name def groupby_q_max( - data: QResolutionPixelTerm, - q_bins: QBins, - dims_to_keep: DimsToKeep, + data: QResolutionPixelTerm, q_bins: QBins, dims_to_keep: DimsToKeep ) -> QResolutionPixelTermGroupedQ: - # TODO Handle multi dim and dims_to_keep! - # Can we use common helper function from bin_in_q? - out = data.groupby('Q', bins=q_bins).max('detector_number') - return QResolutionPixelTermGroupedQ(out) + dims_to_flatten = [dim for dim in data.dims if dim not in dims_to_keep] + dim = uuid.uuid4().hex() + return QResolutionPixelTermGroupedQ( + data.transpose((*dims_to_flatten, *dims_to_keep)) + .flatten(dims_to_flatten, to=dim) + .groupby('Q', bins=q_bins) + .max(dim) + ) def mask_and_compute_resolution_q( From 65738611608e06dcfa4203a86cb35166bc856bf7 Mon Sep 17 00:00:00 2001 From: Simon Heybrock Date: Mon, 11 Nov 2024 06:43:58 +0100 Subject: [PATCH 05/11] Fix methodology --- src/ess/sans/common.py | 2 +- src/ess/sans/normalization.py | 2 +- src/ess/sans/qresolution.py | 68 +++++++++++++++-------------------- 3 files changed, 31 insertions(+), 41 deletions(-) diff --git a/src/ess/sans/common.py b/src/ess/sans/common.py index e3de665c..837ea97c 100644 --- a/src/ess/sans/common.py +++ b/src/ess/sans/common.py @@ -61,7 +61,7 @@ def mask_range( coord = ( da.bins.constituents['data'].coords[dim] - if da.bins is not None + if da.bins is not None and dim in da.bins.coords else da.coords[dim] ) edges = edges.to(unit=coord.unit) diff --git a/src/ess/sans/normalization.py b/src/ess/sans/normalization.py index 4b859ac1..40b17991 100644 --- a/src/ess/sans/normalization.py +++ b/src/ess/sans/normalization.py @@ -392,7 +392,7 @@ def _reduce(part: sc.DataArray, /, *, bands: ProcessedWavelengthBands) -> sc.Dat Q-dependent data, ready for normalization. """ wav = 'wavelength' - if part.bins is not None: + if part.bins is not None and wav in part.bins.coords: # If in event mode the desired wavelength binning has not been applied, we need # it for splitting by bands, or restricting the range in case of a single band. part = part.bin(wavelength=sc.sort(bands.flatten(to=wav), wav)) diff --git a/src/ess/sans/qresolution.py b/src/ess/sans/qresolution.py index c469d672..2c0902c1 100644 --- a/src/ess/sans/qresolution.py +++ b/src/ess/sans/qresolution.py @@ -2,13 +2,13 @@ # Copyright (c) 2024 Scipp contributors (https://github.com/scipp) """Q-resolution calculation for SANS data.""" -import uuid from typing import NewType import scipp as sc from scipp.constants import pi from .common import mask_range +from .conversions import ElasticCoordTransformGraph from .normalization import _reduce from .types import ( CleanQ, @@ -35,21 +35,23 @@ """Moderator time-spread as a function of wavelength.""" -QResolutionPixelTerm = NewType("QResolutionPixelTerm", sc.DataArray) -QResolutionPixelTermGroupedQ = NewType("QResolutionPixelTermGroupedQ", sc.DataArray) +QResolutionByPixel = NewType("QResolutionByPixel", sc.DataArray) +QResolutionByQ = NewType("QResolutionPixelTermGroupedQ", sc.DataArray) QResolutionByWavelength = NewType("QResolutionByWavelength", sc.DataArray) QResolution = NewType("QResolution", sc.DataArray) -def pixel_term( +def q_resolution_by_pixel( detector: CleanQ[SampleRun, Denominator], delta_r: DeltaR, sample_aperture: SampleApertureRadius, source_aperture: SourceApertureRadius, collimation_length: CollimationLength, -) -> QResolutionPixelTerm: + moderator_time_spread: ModeratorTimeSpread, + graph: ElasticCoordTransformGraph, +) -> QResolutionByPixel: """ - Calculate the pixel term for Q-resolution. + Calculate the Q-resolution per pixel. We compute this based on CleanQ[SampleRun, Denominator]. This ensures that @@ -58,37 +60,31 @@ def pixel_term( 3. We do not depend on neutron data, by using the denominator instead of the numerator. """ + detector = detector.transform_coords("L2", graph=graph, keep_inputs=False) L2 = detector.coords["L2"] L3 = sc.reciprocal(sc.reciprocal(collimation_length) + sc.reciprocal(L2)) result = detector.copy(deep=False) - result.data = ( + pixel_term = ( 3 * ((source_aperture / collimation_length) ** 2 + (sample_aperture / L3) ** 2) + (delta_r / L2) ** 2 ) - # TODO - # Give different name, as we do not actually feed into bin_in_q - return QResolutionPixelTerm(result) - - -# TODO What is the point of naming the input CleanQ and output -# CleanSummedQ? It is not sharing functions, so use a different name -def groupby_q_max( - data: QResolutionPixelTerm, q_bins: QBins, dims_to_keep: DimsToKeep -) -> QResolutionPixelTermGroupedQ: - dims_to_flatten = [dim for dim in data.dims if dim not in dims_to_keep] - dim = uuid.uuid4().hex() - return QResolutionPixelTermGroupedQ( - data.transpose((*dims_to_flatten, *dims_to_keep)) - .flatten(dims_to_flatten, to=dim) - .groupby('Q', bins=q_bins) - .max(dim) + inv_lambda2 = sc.reciprocal(detector.coords['wavelength'] ** 2) + Q2 = detector.coords['Q'] ** 2 + result.data = (pi**2 / 3) * inv_lambda2 * pixel_term + Q2 * ( + moderator_time_spread**2 * inv_lambda2 ) + return QResolutionByPixel(result) -def mask_and_compute_resolution_q( - pixel_term: QResolutionPixelTermGroupedQ, - mask: WavelengthMask, - moderator_time_spread: ModeratorTimeSpread, +def q_resolution_by_q( + data: QResolutionByPixel, q_bins: QBins, dims_to_keep: DimsToKeep +) -> QResolutionByQ: + dims = [dim for dim in data.dims if dim not in (*dims_to_keep, 'wavelength')] + return QResolutionByQ(data.bin(Q=q_bins, dim=dims)) + + +def mask_qresolution_in_wavelength( + resolution: QResolutionByQ, mask: WavelengthMask ) -> QResolutionByWavelength: """ Compute the masked Q-resolution in (Q, lambda) space. @@ -98,25 +94,19 @@ def mask_and_compute_resolution_q( term to obtain the Q-resolution. The result is still in (Q, lambda) space. """ if mask is not None: - pixel_term = mask_range(pixel_term, mask=mask) - lambda2 = pixel_term.coords['wavelength'] ** 2 - Q2 = pixel_term.coords['Q'] ** 2 - resolution = ( - pi**2 / (3 * lambda2**2) * pixel_term + Q2 * moderator_time_spread**2 / lambda2 - ) + resolution = mask_range(resolution, mask=mask) return QResolutionByWavelength(resolution) def reduce_resolution_q( data: QResolutionByWavelength, bands: ProcessedWavelengthBands ) -> QResolution: - # TODO Add op argument to allow injection of different reduction functions - return QResolution(_reduce(data, bands, op=sc.max)) + return QResolution(_reduce(data, bands=bands)) providers = ( - pixel_term, - groupby_q_max, - mask_and_compute_resolution_q, + q_resolution_by_pixel, + q_resolution_by_q, + mask_qresolution_in_wavelength, reduce_resolution_q, ) From 1ffeefb5adb236d26121eaf6e99914c758791657 Mon Sep 17 00:00:00 2001 From: Simon Heybrock Date: Mon, 11 Nov 2024 10:49:37 +0100 Subject: [PATCH 06/11] Fix moderator term --- src/ess/sans/qresolution.py | 108 +++++++++++++++++++++++++++++------- src/ess/sans/types.py | 1 + 2 files changed, 88 insertions(+), 21 deletions(-) diff --git a/src/ess/sans/qresolution.py b/src/ess/sans/qresolution.py index 2c0902c1..91d11db5 100644 --- a/src/ess/sans/qresolution.py +++ b/src/ess/sans/qresolution.py @@ -4,19 +4,23 @@ from typing import NewType +import numpy as np import scipp as sc from scipp.constants import pi from .common import mask_range from .conversions import ElasticCoordTransformGraph +from .i_of_q import resample_direct_beam from .normalization import _reduce from .types import ( + CalibratedBeamline, CleanQ, Denominator, DimsToKeep, ProcessedWavelengthBands, QBins, SampleRun, + WavelengthBins, WavelengthMask, ) @@ -28,27 +32,89 @@ SourceApertureRadius = NewType("SourceApertureRadius", sc.Variable) """Source aperture radius, R1.""" SigmaModerator = NewType("SigmaModerator", sc.DataArray) -"""Moderator time spread as a function of wavelength.""" +""" +Moderator wavelength spread as a function of wavelength. + +This is derived from ModeratorTimeSpread. +""" CollimationLength = NewType("CollimationLength", sc.Variable) """Collimation length.""" +ModeratorTimeSpreadFilename = NewType("ModeratorTimeSpreadFilename", str) ModeratorTimeSpread = NewType("ModeratorTimeSpread", sc.DataArray) """Moderator time-spread as a function of wavelength.""" QResolutionByPixel = NewType("QResolutionByPixel", sc.DataArray) -QResolutionByQ = NewType("QResolutionPixelTermGroupedQ", sc.DataArray) QResolutionByWavelength = NewType("QResolutionByWavelength", sc.DataArray) QResolution = NewType("QResolution", sc.DataArray) +def load_isis_moderator_time_spread( + filename: ModeratorTimeSpreadFilename, +) -> ModeratorTimeSpread: + """ + Load moderator time spread from an ISIS moderator file. + + Files looks as follows: + + .. code-block:: text + + Fri 08-Aug-2015, LET exptl data (FWHM/2.35) [...] + + 61 0 0 0 1 61 0 + 0 0 0 0 + 3 (F12.5,2E14.6) + 0.00000 2.257600E+01 0.000000E+00 + 0.50000 2.677152E+01 0.000000E+00 + 1.00000 3.093920E+01 0.000000E+00 + 1.50000 3.507903E+01 0.000000E+00 + 2.00000 3.919100E+01 0.000000E+00 + + The first column is the wavelength in Angstrom, the second is the time spread in + microseconds. The third column is the error on the time spread, which we ignore. + """ + wavelength, time_spread = np.loadtxt( + filename, skiprows=5, usecols=(0, 1), unpack=True + ) + wav = 'wavelength' + return ModeratorTimeSpread( + sc.DataArray( + data=sc.array(dims=[wav], values=time_spread, unit='us'), + coords={wav: sc.array(dims=[wav], values=wavelength, unit='angstrom')}, + ) + ) + + +def moderator_time_spread_to_wavelength_spread( + moderator_time_spread: ModeratorTimeSpread, + beamline: CalibratedBeamline[SampleRun], + graph: ElasticCoordTransformGraph, + wavelength_bins: WavelengthBins, +) -> SigmaModerator: + """ + Convert the moderator time spread to a wavelength spread. + """ + dtof = resample_direct_beam(moderator_time_spread, wavelength_bins=wavelength_bins) + # We would like to "transform" the *data*, but we only have transform_coords, so + # there is some back and forth between data and coords here. + dummy = beamline.broadcast(sizes={**beamline.sizes, **dtof.sizes}) + dummy.data = sc.empty(sizes=dummy.sizes) + da = dummy.assign_coords(tof=dtof.data).transform_coords( + 'wavelength', graph=graph, keep_inputs=False + ) + da.data = da.coords.pop('wavelength') + return SigmaModerator(da.assign_coords(wavelength=wavelength_bins)) + + def q_resolution_by_pixel( detector: CleanQ[SampleRun, Denominator], delta_r: DeltaR, sample_aperture: SampleApertureRadius, source_aperture: SourceApertureRadius, collimation_length: CollimationLength, - moderator_time_spread: ModeratorTimeSpread, + sigma_moderator: SigmaModerator, graph: ElasticCoordTransformGraph, + wavelength_bins: WavelengthBins, ) -> QResolutionByPixel: """ Calculate the Q-resolution per pixel. @@ -60,7 +126,7 @@ def q_resolution_by_pixel( 3. We do not depend on neutron data, by using the denominator instead of the numerator. """ - detector = detector.transform_coords("L2", graph=graph, keep_inputs=False) + detector = detector.transform_coords(('L1', 'L2'), graph=graph, keep_inputs=False) L2 = detector.coords["L2"] L3 = sc.reciprocal(sc.reciprocal(collimation_length) + sc.reciprocal(L2)) result = detector.copy(deep=False) @@ -70,21 +136,18 @@ def q_resolution_by_pixel( ) inv_lambda2 = sc.reciprocal(detector.coords['wavelength'] ** 2) Q2 = detector.coords['Q'] ** 2 - result.data = (pi**2 / 3) * inv_lambda2 * pixel_term + Q2 * ( - moderator_time_spread**2 * inv_lambda2 - ) - return QResolutionByPixel(result) - - -def q_resolution_by_q( - data: QResolutionByPixel, q_bins: QBins, dims_to_keep: DimsToKeep -) -> QResolutionByQ: - dims = [dim for dim in data.dims if dim not in (*dims_to_keep, 'wavelength')] - return QResolutionByQ(data.bin(Q=q_bins, dim=dims)) - - -def mask_qresolution_in_wavelength( - resolution: QResolutionByQ, mask: WavelengthMask + result.data = (pi**2 / 3) * inv_lambda2 * pixel_term + delta_lambda = wavelength_bins[1:] - wavelength_bins[:-1] + sigma_lambda2 = delta_lambda**2 / 12 + sigma_moderator**2 + result += Q2 * (sigma_lambda2 * inv_lambda2) + return QResolutionByPixel(sc.sqrt(result)) + + +def q_resolution_by_wavelength( + data: QResolutionByPixel, + q_bins: QBins, + dims_to_keep: DimsToKeep, + mask: WavelengthMask, ) -> QResolutionByWavelength: """ Compute the masked Q-resolution in (Q, lambda) space. @@ -93,6 +156,8 @@ def mask_qresolution_in_wavelength( what is required for performing the remaining scaling and addition of the moderator term to obtain the Q-resolution. The result is still in (Q, lambda) space. """ + dims = [dim for dim in data.dims if dim not in (*dims_to_keep, 'wavelength')] + resolution = data.bin(Q=q_bins, dim=dims) if mask is not None: resolution = mask_range(resolution, mask=mask) return QResolutionByWavelength(resolution) @@ -105,8 +170,9 @@ def reduce_resolution_q( providers = ( + load_isis_moderator_time_spread, + moderator_time_spread_to_wavelength_spread, q_resolution_by_pixel, - q_resolution_by_q, - mask_qresolution_in_wavelength, + q_resolution_by_wavelength, reduce_resolution_q, ) diff --git a/src/ess/sans/types.py b/src/ess/sans/types.py index f5ff061c..e17c719a 100644 --- a/src/ess/sans/types.py +++ b/src/ess/sans/types.py @@ -16,6 +16,7 @@ from ess.reduce.uncertainty import UncertaintyBroadcastMode as _UncertaintyBroadcastMode BackgroundRun = reduce_t.BackgroundRun +CalibratedBeamline = reduce_t.CalibratedBeamline CalibratedDetector = reduce_t.CalibratedDetector CalibratedMonitor = reduce_t.CalibratedMonitor DetectorData = reduce_t.DetectorData From c17bd6a2a8bc06b9262619be4070610d80e103fa Mon Sep 17 00:00:00 2001 From: Simon Heybrock Date: Mon, 11 Nov 2024 11:31:55 +0100 Subject: [PATCH 07/11] Cleanup and docs --- src/ess/sans/qresolution.py | 132 +++++++++++++++++++++++++++--------- 1 file changed, 101 insertions(+), 31 deletions(-) diff --git a/src/ess/sans/qresolution.py b/src/ess/sans/qresolution.py index 91d11db5..9d833d09 100644 --- a/src/ess/sans/qresolution.py +++ b/src/ess/sans/qresolution.py @@ -14,8 +14,7 @@ from .normalization import _reduce from .types import ( CalibratedBeamline, - CleanQ, - Denominator, + DetectorMasks, DimsToKeep, ProcessedWavelengthBands, QBins, @@ -90,57 +89,95 @@ def moderator_time_spread_to_wavelength_spread( beamline: CalibratedBeamline[SampleRun], graph: ElasticCoordTransformGraph, wavelength_bins: WavelengthBins, + masks: DetectorMasks, ) -> SigmaModerator: """ - Convert the moderator time spread to a wavelength spread. + Convert the moderator time spread to a wavelength spread and mask pixels. + + This resamples the raw moderator time spread to the wavelength bins used in the data + reduction. The result is a DataArray with the time spread as a function of pixel and + wavelength. + + Parameters + ---------- + moderator_time_spread: + Moderator time spread. + beamline: + Beamline geometry information required for converting time spread to wavelength + spread. The latter depends on the pixel position via sample-detector distance + L2. + graph: + Coordinate transformation graph for elastic scattering for computing wavelength. + wavelength_bins: + Binning in wavelength. + masks: + Detector masks. + + Returns + ------- + : + Wavelength spread of the moderator. """ dtof = resample_direct_beam(moderator_time_spread, wavelength_bins=wavelength_bins) # We would like to "transform" the *data*, but we only have transform_coords, so # there is some back and forth between data and coords here. dummy = beamline.broadcast(sizes={**beamline.sizes, **dtof.sizes}) dummy.data = sc.empty(sizes=dummy.sizes) - da = dummy.assign_coords(tof=dtof.data).transform_coords( - 'wavelength', graph=graph, keep_inputs=False + da = ( + dummy.assign_coords(tof=dtof.data) + .assign_masks(masks) + .transform_coords('wavelength', graph=graph, keep_inputs=False) ) da.data = da.coords.pop('wavelength') return SigmaModerator(da.assign_coords(wavelength=wavelength_bins)) def q_resolution_by_pixel( - detector: CleanQ[SampleRun, Denominator], - delta_r: DeltaR, - sample_aperture: SampleApertureRadius, + sigma_moderator: SigmaModerator, source_aperture: SourceApertureRadius, + sample_aperture: SampleApertureRadius, collimation_length: CollimationLength, - sigma_moderator: SigmaModerator, + delta_r: DeltaR, graph: ElasticCoordTransformGraph, - wavelength_bins: WavelengthBins, ) -> QResolutionByPixel: """ - Calculate the Q-resolution per pixel. - - We compute this based on CleanQ[SampleRun, Denominator]. This ensures that - - 1. We get the correct Q-value, based on wavelength binning used elsewhere. - 2. Masks are included. - 3. We do not depend on neutron data, by using the denominator instead of the - numerator. + Calculate the Q-resolution per pixel and wavelength. + + Parameters + ---------- + sigma_moderator: + Moderator wavelength spread as a function of pixel and wavelength. + source_aperture: + Source aperture radius, R1. + sample_aperture: + Sample aperture radius, R2. + collimation_length: + Collimation length. + delta_r: + Virtual ring width on the detector. + graph: + Coordinate transformation graph for elastic scattering. + + Returns + ------- + : + Q-resolution per pixel and wavelength. """ - detector = detector.transform_coords(('L1', 'L2'), graph=graph, keep_inputs=False) - L2 = detector.coords["L2"] + wavelength_bins = sigma_moderator.coords['wavelength'] + delta_lambda = wavelength_bins[1:] - wavelength_bins[:-1] + result = sigma_moderator.assign_coords( + wavelength=sc.midpoints(wavelength_bins) + ).transform_coords('Q', graph=graph, keep_inputs=True) + L2 = result.coords["L2"] L3 = sc.reciprocal(sc.reciprocal(collimation_length) + sc.reciprocal(L2)) - result = detector.copy(deep=False) pixel_term = ( 3 * ((source_aperture / collimation_length) ** 2 + (sample_aperture / L3) ** 2) + (delta_r / L2) ** 2 ) - inv_lambda2 = sc.reciprocal(detector.coords['wavelength'] ** 2) - Q2 = detector.coords['Q'] ** 2 - result.data = (pi**2 / 3) * inv_lambda2 * pixel_term - delta_lambda = wavelength_bins[1:] - wavelength_bins[:-1] - sigma_lambda2 = delta_lambda**2 / 12 + sigma_moderator**2 - result += Q2 * (sigma_lambda2 * inv_lambda2) - return QResolutionByPixel(sc.sqrt(result)) + sigma_lambda2 = result.data**2 + delta_lambda**2 / 12 + result.data = (pi**2 / 3) * pixel_term + result.coords['Q'] ** 2 * sigma_lambda2 + inv_lambda = sc.reciprocal(result.coords['wavelength']) + return QResolutionByPixel(sc.sqrt(result) * inv_lambda) def q_resolution_by_wavelength( @@ -152,9 +189,21 @@ def q_resolution_by_wavelength( """ Compute the masked Q-resolution in (Q, lambda) space. - CleanSummedQ has been summed over pixels but not over wavelengths. This is exactly - what is required for performing the remaining scaling and addition of the moderator - term to obtain the Q-resolution. The result is still in (Q, lambda) space. + Parameters + ---------- + data: + Q-resolution per pixel and wavelength. + q_bins: + Binning in Q. + dims_to_keep: + Detector dimensions that should not be reduced. + mask: + Wavelength mask. + + Returns + ------- + : + Q-resolution in (Q, lambda) space. """ dims = [dim for dim in data.dims if dim not in (*dims_to_keep, 'wavelength')] resolution = data.bin(Q=q_bins, dim=dims) @@ -166,6 +215,27 @@ def q_resolution_by_wavelength( def reduce_resolution_q( data: QResolutionByWavelength, bands: ProcessedWavelengthBands ) -> QResolution: + """ + Q-resolution reduced over wavelength dimension. + + The result depends only Q, or (Q, wavelength_band) if wavelength bands are used. + + Note that the result is *binned data*, with each bin entry representing a specific + input pixel and wavelength. The final solution can be obtained, e.g., by computing + `resolution.bins.max()`, assuming `resolution is the output of this function. + + Parameters + ---------- + data: + Q-resolution in (Q, lambda) space. + bands: + Wavelength bands. + + Returns + ------- + : + Reduced Q-resolution. + """ return QResolution(_reduce(data, bands=bands)) From fd26aa6b73f53cf5638786bfe16e19c95a162135 Mon Sep 17 00:00:00 2001 From: Simon Heybrock Date: Mon, 11 Nov 2024 14:43:53 +0100 Subject: [PATCH 08/11] Do not use expensive event-data mask handling when not necessary --- src/ess/sans/common.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/ess/sans/common.py b/src/ess/sans/common.py index 837ea97c..4073f22e 100644 --- a/src/ess/sans/common.py +++ b/src/ess/sans/common.py @@ -66,7 +66,7 @@ def mask_range( ) edges = edges.to(unit=coord.unit) lu = sc.DataArray(data=mask.data, coords={dim: edges}) - if da.bins is not None: + if da.bins is not None and dim in da.bins.coords: if dim not in da.coords: underlying = da.bins.coords[dim] new_bins = np.union1d( From c1d8037897d1006085135935df8297550ee010d9 Mon Sep 17 00:00:00 2001 From: Simon Heybrock Date: Mon, 11 Nov 2024 14:54:45 +0100 Subject: [PATCH 09/11] Small speedup. --- src/ess/sans/qresolution.py | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/src/ess/sans/qresolution.py b/src/ess/sans/qresolution.py index 9d833d09..6d61f716 100644 --- a/src/ess/sans/qresolution.py +++ b/src/ess/sans/qresolution.py @@ -174,8 +174,12 @@ def q_resolution_by_pixel( 3 * ((source_aperture / collimation_length) ** 2 + (sample_aperture / L3) ** 2) + (delta_r / L2) ** 2 ) - sigma_lambda2 = result.data**2 + delta_lambda**2 / 12 - result.data = (pi**2 / 3) * pixel_term + result.coords['Q'] ** 2 * sigma_lambda2 + # Written in multiple steps to avoid allocation of intermediate arrays. This is + # makes this step about twice as fast (but computing Q above dominates anyway). + result.data *= result.data + result.data += delta_lambda**2 / 12 + result.data *= result.coords['Q'] ** 2 + result.data += (pi**2 / 3) * pixel_term inv_lambda = sc.reciprocal(result.coords['wavelength']) return QResolutionByPixel(sc.sqrt(result) * inv_lambda) From ff331281b84f30ae55bb555e5a805a71a82b2061 Mon Sep 17 00:00:00 2001 From: Simon Heybrock Date: Tue, 12 Nov 2024 11:27:11 +0100 Subject: [PATCH 10/11] Do not correct for gravity in moderator spread conversion to wavelength --- src/ess/sans/qresolution.py | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/src/ess/sans/qresolution.py b/src/ess/sans/qresolution.py index 6d61f716..d67d5e1d 100644 --- a/src/ess/sans/qresolution.py +++ b/src/ess/sans/qresolution.py @@ -9,7 +9,7 @@ from scipp.constants import pi from .common import mask_range -from .conversions import ElasticCoordTransformGraph +from .conversions import ElasticCoordTransformGraph, sans_elastic from .i_of_q import resample_direct_beam from .normalization import _reduce from .types import ( @@ -87,7 +87,6 @@ def load_isis_moderator_time_spread( def moderator_time_spread_to_wavelength_spread( moderator_time_spread: ModeratorTimeSpread, beamline: CalibratedBeamline[SampleRun], - graph: ElasticCoordTransformGraph, wavelength_bins: WavelengthBins, masks: DetectorMasks, ) -> SigmaModerator: @@ -106,8 +105,6 @@ def moderator_time_spread_to_wavelength_spread( Beamline geometry information required for converting time spread to wavelength spread. The latter depends on the pixel position via sample-detector distance L2. - graph: - Coordinate transformation graph for elastic scattering for computing wavelength. wavelength_bins: Binning in wavelength. masks: @@ -118,6 +115,10 @@ def moderator_time_spread_to_wavelength_spread( : Wavelength spread of the moderator. """ + # We want to convert a small time-of-flight spread to a wavelength spread. As this + # time spread is assumed to be symmetric around the base time-of-flight we do not + # want to correct for gravity here as it would skew the result. + graph = sans_elastic(correct_for_gravity=False) dtof = resample_direct_beam(moderator_time_spread, wavelength_bins=wavelength_bins) # We would like to "transform" the *data*, but we only have transform_coords, so # there is some back and forth between data and coords here. From a0cb4cc5f46df0529369fc37bc17ee2c53f0654d Mon Sep 17 00:00:00 2001 From: Simon Heybrock Date: Tue, 12 Nov 2024 12:00:01 +0100 Subject: [PATCH 11/11] Add information on limitations of the Q-resolution implementation. --- src/ess/sans/qresolution.py | 11 +++++++++-- 1 file changed, 9 insertions(+), 2 deletions(-) diff --git a/src/ess/sans/qresolution.py b/src/ess/sans/qresolution.py index d67d5e1d..d6df2d0e 100644 --- a/src/ess/sans/qresolution.py +++ b/src/ess/sans/qresolution.py @@ -144,6 +144,12 @@ def q_resolution_by_pixel( """ Calculate the Q-resolution per pixel and wavelength. + This is a Gaussian approximation to the Q-resolution (Mildner and Carpenter). It is + likely not sufficient for the long ESS pulse. Based on communication with Andrew + Jackson this "approximation works OK when you have all your Q points from a single + planar detector (as on SANS2D) and was implemented by Richard Heenan as a 'hack' to + get a resolution value of some kind". + Parameters ---------- sigma_moderator: @@ -181,8 +187,9 @@ def q_resolution_by_pixel( result.data += delta_lambda**2 / 12 result.data *= result.coords['Q'] ** 2 result.data += (pi**2 / 3) * pixel_term - inv_lambda = sc.reciprocal(result.coords['wavelength']) - return QResolutionByPixel(sc.sqrt(result) * inv_lambda) + result = sc.sqrt(result) + result /= result.coords['wavelength'] + return QResolutionByPixel(result) def q_resolution_by_wavelength(