diff --git a/docs/user-guide/common/beam-center-finder.ipynb b/docs/user-guide/common/beam-center-finder.ipynb index 67693e42..6027a246 100644 --- a/docs/user-guide/common/beam-center-finder.ipynb +++ b/docs/user-guide/common/beam-center-finder.ipynb @@ -451,7 +451,7 @@ "kwargs = dict( # noqa: C408\n", " workflow=workflow,\n", " detector=detector['data'],\n", - " norm=workflow.compute(NormWavelengthTerm[SampleRun]),\n", + " norm=workflow.compute(CleanDirectBeam),\n", ")" ] }, diff --git a/docs/user-guide/isis/sans2d.ipynb b/docs/user-guide/isis/sans2d.ipynb index e98bbecd..808d3320 100644 --- a/docs/user-guide/isis/sans2d.ipynb +++ b/docs/user-guide/isis/sans2d.ipynb @@ -319,7 +319,10 @@ " WavelengthMonitor[BackgroundRun, Incident],\n", " WavelengthMonitor[BackgroundRun, Transmission],\n", ")\n", - "parts = (CleanSummedQ[SampleRun, Numerator], CleanSummedQ[SampleRun, Denominator])\n", + "parts = (\n", + " WavelengthScaledQ[SampleRun, Numerator],\n", + " WavelengthScaledQ[SampleRun, Denominator],\n", + ")\n", "iofqs = (IofQ[SampleRun], IofQ[BackgroundRun], BackgroundSubtractedIofQ)\n", "keys = (*monitors, MaskedData[SampleRun], *parts, *iofqs)\n", "\n", @@ -335,12 +338,12 @@ "\n", "wavelength = workflow.compute(WavelengthBins)\n", "display(\n", - " results[CleanSummedQ[SampleRun, Numerator]]\n", + " results[WavelengthScaledQ[SampleRun, Numerator]]\n", " .hist(wavelength=wavelength)\n", " .transpose()\n", " .plot(norm='log')\n", ")\n", - "display(results[CleanSummedQ[SampleRun, Denominator]].plot(norm='log'))\n", + "display(results[WavelengthScaledQ[SampleRun, Denominator]].plot(norm='log'))\n", "parts = {str(key): results[key].sum('wavelength') for key in parts}\n", "display(sc.plot(parts, norm='log'))\n", "\n", @@ -429,7 +432,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.10.14" + "version": "3.10.12" } }, "nbformat": 4, diff --git a/docs/user-guide/isis/zoom.ipynb b/docs/user-guide/isis/zoom.ipynb index 99449935..f572ff3e 100644 --- a/docs/user-guide/isis/zoom.ipynb +++ b/docs/user-guide/isis/zoom.ipynb @@ -233,7 +233,10 @@ " WavelengthMonitor[SampleRun, Incident],\n", " WavelengthMonitor[SampleRun, Transmission],\n", ")\n", - "parts = (CleanSummedQ[SampleRun, Numerator], CleanSummedQ[SampleRun, Denominator])\n", + "parts = (\n", + " WavelengthScaledQ[SampleRun, Numerator],\n", + " WavelengthScaledQ[SampleRun, Denominator],\n", + ")\n", "iofqs = (IofQ[SampleRun],)\n", "keys = (*monitors, MaskedData[SampleRun], *parts, *iofqs)\n", "\n", @@ -249,12 +252,12 @@ "\n", "wavelength = workflow.compute(WavelengthBins)\n", "display(\n", - " results[CleanSummedQ[SampleRun, Numerator]]\n", + " results[WavelengthScaledQ[SampleRun, Numerator]]\n", " .hist(wavelength=wavelength)\n", " .transpose()\n", " .plot(norm='log')\n", ")\n", - "display(results[CleanSummedQ[SampleRun, Denominator]].plot(norm='log'))\n", + "display(results[WavelengthScaledQ[SampleRun, Denominator]].plot(norm='log'))\n", "parts = {str(key): results[key] for key in parts}\n", "parts = {\n", " key: val.sum('wavelength') if val.bins is None else val.hist()\n", @@ -309,7 +312,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.10.12" + "version": "3.10.14" } }, "nbformat": 4, diff --git a/src/ess/sans/beam_center_finder.py b/src/ess/sans/beam_center_finder.py index 18372ade..f1cb3ff9 100644 --- a/src/ess/sans/beam_center_finder.py +++ b/src/ess/sans/beam_center_finder.py @@ -14,12 +14,12 @@ from .logging import get_logger from .types import ( BeamCenter, + CleanDirectBeam, DetectorBankSizes, DimsToKeep, IofQ, MaskedData, NeXusDetector, - NormWavelengthTerm, QBins, ReturnEvents, SampleRun, @@ -175,9 +175,7 @@ def _iofq_in_quadrants( workflow[NeXusDetector[SampleRun]] = sc.DataGroup(data=detector[sel]) # MaskedData would be computed automatically, but we did it above already workflow[MaskedData[SampleRun]] = calibrated[sel] - workflow[NormWavelengthTerm[SampleRun]] = ( - norm if norm.dims == ('wavelength',) else norm[sel] - ) + workflow[CleanDirectBeam] = norm if norm.dims == ('wavelength',) else norm[sel] out[quad] = workflow.compute(IofQ[SampleRun]) return out @@ -364,7 +362,7 @@ def beam_center_from_iofq( keys = ( NeXusDetector[SampleRun], MaskedData[SampleRun], - NormWavelengthTerm[SampleRun], + CleanDirectBeam, ElasticCoordTransformGraph, ) workflow = workflow.copy() @@ -373,7 +371,7 @@ def beam_center_from_iofq( results = workflow.compute(keys) detector = results[NeXusDetector[SampleRun]]['data'] data = results[MaskedData[SampleRun]] - norm = results[NormWavelengthTerm[SampleRun]] + norm = results[CleanDirectBeam] graph = results[ElasticCoordTransformGraph] # Avoid reloading the detector diff --git a/src/ess/sans/conversions.py b/src/ess/sans/conversions.py index 1e4bb823..ce157dd9 100644 --- a/src/ess/sans/conversions.py +++ b/src/ess/sans/conversions.py @@ -9,22 +9,30 @@ ) from scippneutron.conversion.graph import beamline, tof +from ess.reduce.uncertainty import broadcast_uncertainties + from .common import mask_range from .types import ( CleanQ, CleanQxy, + CleanSummedQ, + CleanSummedQxy, CleanWavelength, - CleanWavelengthMasked, CorrectForGravity, + Denominator, IofQPart, MaskedData, + MonitorTerm, MonitorType, Numerator, RunType, ScatteringRunType, TofMonitor, + UncertaintyBroadcastMode, WavelengthMask, WavelengthMonitor, + WavelengthScaledQ, + WavelengthScaledQxy, ) @@ -163,12 +171,44 @@ def detector_to_wavelength( ) -def mask_wavelength( - da: CleanWavelength[ScatteringRunType, IofQPart], mask: WavelengthMask -) -> CleanWavelengthMasked[ScatteringRunType, IofQPart]: +def mask_wavelength_q( + da: CleanSummedQ[ScatteringRunType, Numerator], mask: WavelengthMask +) -> WavelengthScaledQ[ScatteringRunType, Numerator]: + if mask is not None: + da = mask_range(da, mask=mask) + return WavelengthScaledQ[ScatteringRunType, Numerator](da) + + +def mask_wavelength_qxy( + da: CleanSummedQxy[ScatteringRunType, Numerator], mask: WavelengthMask +) -> WavelengthScaledQxy[ScatteringRunType, Numerator]: + if mask is not None: + da = mask_range(da, mask=mask) + return WavelengthScaledQxy[ScatteringRunType, Numerator](da) + + +def mask_and_scale_wavelength_q( + da: CleanSummedQ[ScatteringRunType, Denominator], + mask: WavelengthMask, + wavelength_term: MonitorTerm[ScatteringRunType], + uncertainties: UncertaintyBroadcastMode, +) -> WavelengthScaledQ[ScatteringRunType, Denominator]: + da = da * broadcast_uncertainties(wavelength_term, prototype=da, mode=uncertainties) + if mask is not None: + da = mask_range(da, mask=mask) + return WavelengthScaledQ[ScatteringRunType, Denominator](da) + + +def mask_and_scale_wavelength_qxy( + da: CleanSummedQxy[ScatteringRunType, Denominator], + mask: WavelengthMask, + wavelength_term: MonitorTerm[ScatteringRunType], + uncertainties: UncertaintyBroadcastMode, +) -> WavelengthScaledQxy[ScatteringRunType, Denominator]: + da = da * broadcast_uncertainties(wavelength_term, prototype=da, mode=uncertainties) if mask is not None: da = mask_range(da, mask=mask) - return CleanWavelengthMasked[ScatteringRunType, IofQPart](da) + return WavelengthScaledQxy[ScatteringRunType, Denominator](da) def _compute_Q( @@ -186,7 +226,7 @@ def _compute_Q( def compute_Q( - data: CleanWavelengthMasked[ScatteringRunType, IofQPart], + data: CleanWavelength[ScatteringRunType, IofQPart], graph: ElasticCoordTransformGraph, ) -> CleanQ[ScatteringRunType, IofQPart]: """ @@ -198,7 +238,7 @@ def compute_Q( def compute_Qxy( - data: CleanWavelengthMasked[ScatteringRunType, IofQPart], + data: CleanWavelength[ScatteringRunType, IofQPart], graph: ElasticCoordTransformGraph, ) -> CleanQxy[ScatteringRunType, IofQPart]: """ @@ -214,7 +254,10 @@ def compute_Qxy( sans_monitor, monitor_to_wavelength, detector_to_wavelength, - mask_wavelength, + mask_wavelength_q, + mask_wavelength_qxy, + mask_and_scale_wavelength_q, + mask_and_scale_wavelength_qxy, compute_Q, compute_Qxy, ) diff --git a/src/ess/sans/direct_beam.py b/src/ess/sans/direct_beam.py index 3fd2fe8c..db2bd86a 100644 --- a/src/ess/sans/direct_beam.py +++ b/src/ess/sans/direct_beam.py @@ -8,7 +8,6 @@ from .types import ( BackgroundRun, BackgroundSubtractedIofQ, - CleanSummedQ, Denominator, DirectBeam, Numerator, @@ -16,6 +15,7 @@ SampleRun, WavelengthBands, WavelengthBins, + WavelengthScaledQ, ) @@ -103,16 +103,16 @@ def direct_beam(*, workflow: Pipeline, I0: sc.Variable, niter: int = 5) -> list[ wavelength_bins = workflow.compute(WavelengthBins) parts = ( - CleanSummedQ[SampleRun, Numerator], - CleanSummedQ[SampleRun, Denominator], - CleanSummedQ[BackgroundRun, Numerator], - CleanSummedQ[BackgroundRun, Denominator], + WavelengthScaledQ[SampleRun, Numerator], + WavelengthScaledQ[SampleRun, Denominator], + WavelengthScaledQ[BackgroundRun, Numerator], + WavelengthScaledQ[BackgroundRun, Denominator], ) parts = workflow.compute(parts) # Convert events to histograms to make normalization (in every iteration) cheap for key in [ - CleanSummedQ[SampleRun, Numerator], - CleanSummedQ[BackgroundRun, Numerator], + WavelengthScaledQ[SampleRun, Numerator], + WavelengthScaledQ[BackgroundRun, Numerator], ]: parts[key] = parts[key].hist(wavelength=wavelength_bins) @@ -121,8 +121,8 @@ def direct_beam(*, workflow: Pipeline, I0: sc.Variable, niter: int = 5) -> list[ parts = {key: sc.values(result) for key, result in parts.items()} for key, part in parts.items(): workflow[key] = part - sample0 = parts[CleanSummedQ[SampleRun, Denominator]] - background0 = parts[CleanSummedQ[BackgroundRun, Denominator]] + sample0 = parts[WavelengthScaledQ[SampleRun, Denominator]] + background0 = parts[WavelengthScaledQ[BackgroundRun, Denominator]] results = [] @@ -158,8 +158,8 @@ def direct_beam(*, workflow: Pipeline, I0: sc.Variable, niter: int = 5) -> list[ db.coords['wavelength'] = sc.midpoints( db.coords['wavelength'], dim='wavelength' ) - workflow[CleanSummedQ[SampleRun, Denominator]] = sample0 * db - workflow[CleanSummedQ[BackgroundRun, Denominator]] = background0 * db + workflow[WavelengthScaledQ[SampleRun, Denominator]] = sample0 * db + workflow[WavelengthScaledQ[BackgroundRun, Denominator]] = background0 * db results.append( { diff --git a/src/ess/sans/i_of_q.py b/src/ess/sans/i_of_q.py index bb4d2995..8348b6fc 100644 --- a/src/ess/sans/i_of_q.py +++ b/src/ess/sans/i_of_q.py @@ -112,7 +112,12 @@ def resample_direct_beam( The direct beam function resampled to the requested resolution. """ if direct_beam is None: - return CleanDirectBeam(None) + return CleanDirectBeam( + sc.DataArray( + sc.ones(dims=wavelength_bins.dims, shape=[len(wavelength_bins) - 1]), + coords={'wavelength': wavelength_bins}, + ) + ) if sc.identical(direct_beam.coords['wavelength'], wavelength_bins): return direct_beam if direct_beam.variances is not None: diff --git a/src/ess/sans/normalization.py b/src/ess/sans/normalization.py index d38eaadf..311619c9 100644 --- a/src/ess/sans/normalization.py +++ b/src/ess/sans/normalization.py @@ -9,8 +9,6 @@ CalibratedDetector, CleanDirectBeam, CleanMonitor, - CleanSummedQ, - CleanSummedQxy, CleanWavelength, Denominator, DetectorMasks, @@ -18,12 +16,15 @@ EmptyBeamRun, Incident, IofQ, + IofQPart, IofQxy, LabFrameTransform, MaskedSolidAngle, - NormWavelengthTerm, + MonitorTerm, Numerator, ProcessedWavelengthBands, + ReducedQ, + ReducedQxy, ReturnEvents, ScatteringRunType, SolidAngle, @@ -32,6 +33,8 @@ TransmissionRun, WavelengthBands, WavelengthBins, + WavelengthScaledQ, + WavelengthScaledQxy, ) @@ -169,28 +172,18 @@ def transmission_fraction( return TransmissionFraction[ScatteringRunType](frac) -def iofq_norm_wavelength_term( +def norm_monitor_term( incident_monitor: CleanMonitor[ScatteringRunType, Incident], transmission_fraction: TransmissionFraction[ScatteringRunType], - direct_beam: CleanDirectBeam, - uncertainties: UncertaintyBroadcastMode, -) -> NormWavelengthTerm[ScatteringRunType]: +) -> MonitorTerm[ScatteringRunType]: """ - Compute the wavelength-dependent contribution to the denominator term for the I(Q) - normalization. - Keeping this as a separate function allows us to compute it once during the - iterations for finding the beam center, while the solid angle is recomputed - for each iteration. - - This is basically: - ``incident_monitor * transmission_fraction * direct_beam`` - If the direct beam is not supplied, it is assumed to be 1. + Compute the monitor-dependent contribution to the denominator term of I(Q). + + This is basically ``incident_monitor * transmission_fraction``. - Because the multiplication between the ``incident_monitor * transmission_fraction`` - (pixel-independent) and the direct beam (potentially pixel-dependent) consists of a - broadcast operation which would introduce correlations, variances of the direct - beam are dropped or replaced by an upper-bound estimation, depending on the - configured mode. + Keeping the monitor term separate from the detector term allows us to compute + the latter only once when repeatedly processing chunks of events in streamed data + processing, e.g., for live data reduction. Parameters ---------- @@ -198,94 +191,44 @@ def iofq_norm_wavelength_term( The incident monitor data (depends on wavelength). transmission_fraction: The transmission fraction (depends on wavelength). - direct_beam: - The direct beam function (depends on wavelength). - uncertainties: - The mode for broadcasting uncertainties. See - :py:class:`ess.reduce.uncertainty.UncertaintyBroadcastMode` for details. Returns ------- : - Wavelength-dependent term - (incident_monitor * transmission_fraction * direct_beam) to be used for - the denominator of the SANS I(Q) normalization. - Used by :py:func:`iofq_denominator`. + Monitor-dependent factor of the normalization term for the SANS I(Q). """ out = incident_monitor * transmission_fraction - if direct_beam is not None: - # Make wavelength the inner dim - dims = list(direct_beam.dims) - dims.remove('wavelength') - dims.append('wavelength') - direct_beam = direct_beam.transpose(dims) - out = direct_beam * broadcast_uncertainties( - out, prototype=direct_beam, mode=uncertainties - ) # Convert wavelength coordinate to midpoints for future histogramming out.coords['wavelength'] = sc.midpoints(out.coords['wavelength']) - return NormWavelengthTerm[ScatteringRunType](out) + return MonitorTerm[ScatteringRunType](out) -def iofq_denominator( - wavelength_term: NormWavelengthTerm[ScatteringRunType], +def norm_detector_term( solid_angle: MaskedSolidAngle[ScatteringRunType], + direct_beam: CleanDirectBeam, uncertainties: UncertaintyBroadcastMode, ) -> CleanWavelength[ScatteringRunType, Denominator]: """ - Compute the denominator term for the I(Q) normalization. - - In a SANS experiment, the scattering cross section :math:`I(Q)` is defined as - (`Heenan et al. 1997 `_): - - .. math:: + Compute the detector-dependent contribution to the denominator term of I(Q). - I(Q) = \\frac{\\partial\\Sigma{Q}}{\\partial\\Omega} = \\frac{A_{H} \\Sigma_{R,\\lambda\\subset Q} C(R, \\lambda)}{A_{M} t \\Sigma_{R,\\lambda\\subset Q}M(\\lambda)T(\\lambda)D(\\lambda)\\Omega(R)} - - where :math:`A_{H}` is the area of a mask (which avoids saturating the detector) - placed between the monitor of area :math:`A_{M}` and the main detector. - :math:`\\Omega` is the detector solid angle, and :math:`C` is the count rate on the - main detector, which depends on the position :math:`R` and the wavelength. - :math:`t` is the sample thickness, :math:`M` represents the incident monitor count - rate for the sample run, and :math:`T` is known as the transmission fraction. - - Note that the incident monitor used to compute the transmission fraction is not - necessarily the same as :math:`M`, as the transmission fraction is usually computed - from a separate 'transmission' run (in the 'sample' run, the transmission monitor is - commonly moved out of the beam path, to avoid polluting the sample detector signal). - - Finally, :math:`D` is the 'direct beam function', and is defined as - - .. math:: - - D(\\lambda) = \\frac{\\eta(\\lambda)}{\\eta_{M}(\\lambda)} \\frac{A_{H}}{A_{M}} - - where :math:`\\eta` and :math:`\\eta_{M}` are the detector and monitor - efficiencies, respectively. - - Hence, in order to normalize the main detector counts :math:`C`, we need compute the - transmission fraction :math:`T(\\lambda)`, the direct beam function - :math:`D(\\lambda)` and the solid angle :math:`\\Omega(R)`. + This is basically ``solid_angle * direct_beam``. + If the direct beam is not supplied, it is assumed to be 1. - The denominator is then simply: - :math:`M_{\\lambda} T_{\\lambda} D_{\\lambda} \\Omega_{R}`, - which is equivalent to ``wavelength_term * solid_angle``. - The ``wavelength_term`` includes all but the ``solid_angle`` and is computed by - :py:func:`iofq_norm_wavelength_term_sample` or - :py:func:`iofq_norm_wavelength_term_background`. + Keeping the monitor term separate from the detector term allows us to compute the + the latter only once when repeatedly processing chunks of events in streamed data + processing, e.g., for live data reduction. - Because the multiplication between the wavelength dependent terms - and the pixel dependent term (solid angle) consists of a broadcast operation which - would introduce correlations, variances are dropped or replaced by an upper-bound - estimation, depending on the configured mode. + Because the multiplication between the pixel-dependent solid angle and the + (potentially) pixel-independent direct beam constitutes a broadcast operation which + would introduce correlations, variances of the direct beam are dropped or replaced + by an upper-bound estimation, depending on the configured mode. Parameters ---------- - wavelength_term: - The term that depends on wavelength, computed by - :py:func:`iofq_norm_wavelength_term`. solid_angle: The solid angle of the detector pixels, as viewed from the sample position. + direct_beam: + The direct beam function (depends on wavelength). uncertainties: The mode for broadcasting uncertainties. See :py:class:`ess.reduce.uncertainty.UncertaintyBroadcastMode` for details. @@ -293,12 +236,19 @@ def iofq_denominator( Returns ------- : - The denominator for the SANS I(Q) normalization. - """ # noqa: E501 - denominator = solid_angle * broadcast_uncertainties( - wavelength_term, prototype=solid_angle, mode=uncertainties + Detector-dependent factor of the normalization term for the SANS I(Q). + """ + # Make wavelength the inner dim + dims = list(direct_beam.dims) + dims.remove('wavelength') + dims.append('wavelength') + direct_beam = direct_beam.transpose(dims) + out = solid_angle * broadcast_uncertainties( + direct_beam, prototype=solid_angle, mode=uncertainties ) - return CleanWavelength[ScatteringRunType, Denominator](denominator) + # Convert wavelength coordinate to midpoints for future histogramming + out.coords['wavelength'] = sc.midpoints(out.coords['wavelength']) + return CleanWavelength[ScatteringRunType, Denominator](out) def process_wavelength_bands( @@ -342,13 +292,46 @@ def _normalize( denominator: sc.DataArray, return_events: ReturnEvents, uncertainties: UncertaintyBroadcastMode, - wavelength_bands: ProcessedWavelengthBands, ) -> sc.DataArray: """ - Perform normalization of counts as a function of Q. - If the numerator contains events, we use the sc.lookup function to perform the + Perform normalization of counts as a function of Q to obtain I(Q) or I(Qx, Qy). division. + In a SANS experiment, the scattering cross section :math:`I(Q)` is defined as + (`Heenan et al. 1997 `_): + + .. math:: + + I(Q) = \\frac{\\partial\\Sigma{Q}}{\\partial\\Omega} = \\frac{A_{H} \\Sigma_{R,\\lambda\\subset Q} C(R, \\lambda)}{A_{M} t \\Sigma_{R,\\lambda\\subset Q}M(\\lambda)T(\\lambda)D(\\lambda)\\Omega(R)} + + where :math:`A_{H}` is the area of a mask (which avoids saturating the detector) + placed between the monitor of area :math:`A_{M}` and the main detector. + :math:`\\Omega` is the detector solid angle, and :math:`C` is the count rate on the + main detector, which depends on the position :math:`R` and the wavelength. + :math:`t` is the sample thickness, :math:`M` represents the incident monitor count + rate for the sample run, and :math:`T` is known as the transmission fraction. + + Note that the incident monitor used to compute the transmission fraction is not + necessarily the same as :math:`M`, as the transmission fraction is usually computed + from a separate 'transmission' run (in the 'sample' run, the transmission monitor is + commonly moved out of the beam path, to avoid polluting the sample detector signal). + + Finally, :math:`D` is the 'direct beam function', and is defined as + + .. math:: + + D(\\lambda) = \\frac{\\eta(\\lambda)}{\\eta_{M}(\\lambda)} \\frac{A_{H}}{A_{M}} + + where :math:`\\eta` and :math:`\\eta_{M}` are the detector and monitor + efficiencies, respectively. + + Hence, in order to normalize the main detector counts :math:`C`, we need compute the + transmission fraction :math:`T(\\lambda)`, the direct beam function + :math:`D(\\lambda)` and the solid angle :math:`\\Omega(R)`. + + The denominator is then simply: + :math:`M_{\\lambda} T_{\\lambda} D_{\\lambda} \\Omega_{R}`. + Parameters ---------- numerator: @@ -359,11 +342,6 @@ def _normalize( contain histogrammed data. return_events: Whether to return the result as event data or histogrammed data. - wavelength_bands: - Defines bands in wavelength that can be used to separate different wavelength - ranges that contribute to different regions in Q space. Note that this needs to - be defined, so if all wavelengths should be used, this should simply be a start - and end edges that encompass the entire wavelength range. uncertainties: The mode for broadcasting uncertainties. See :py:class:`ess.reduce.uncertainty.UncertaintyBroadcastMode` for details. @@ -371,35 +349,8 @@ def _normalize( Returns ------- : - The input data normalized by the supplied denominator. - """ - wav = 'wavelength' - wavelength_bounds = sc.sort(wavelength_bands.flatten(to=wav), wav) - if numerator.bins is not None: - # 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. - numerator = numerator.bin(wavelength=wavelength_bounds) - - def _reduce(da: sc.DataArray) -> sc.DataArray: - if da.sizes[wav] == 1: # Can avoid costly event-data da.bins.concat - return da.squeeze(wav) - return da.sum(wav) if da.bins is None else da.bins.concat(wav) - - num_parts = [] - denom_parts = [] - for wav_range in sc.collapse(wavelength_bands, keep=wav).values(): - num_parts.append(_reduce(numerator[wav, wav_range[0] : wav_range[1]])) - denom_parts.append(_reduce(denominator[wav, wav_range[0] : wav_range[1]])) - band_dim = (set(wavelength_bands.dims) - {'wavelength'}).pop() - if len(num_parts) == 1: - numerator = num_parts[0] - denominator = denom_parts[0] - else: - numerator = sc.concat(num_parts, band_dim) - denominator = sc.concat(denom_parts, band_dim) - numerator.coords[wav] = wavelength_bands.squeeze() - denominator.coords[wav] = wavelength_bands.squeeze() - + Normalized I(Q) or I(Qx, Qy). + """ # noqa: E501 if return_events and numerator.bins is not None: # Naive event-mode normalization is not correct if norm-term has variances. # See https://doi.org/10.3233/JNR-220049 for context. @@ -414,12 +365,65 @@ def _reduce(da: sc.DataArray) -> sc.DataArray: return numerator +def _do_reduce(da: sc.DataArray) -> sc.DataArray: + wav = 'wavelength' + if da.sizes[wav] == 1: # Can avoid costly event-data da.bins.concat + return da.squeeze(wav) + return da.sum(wav) if da.bins is None else da.bins.concat(wav) + + +def _reduce(part: sc.DataArray, /, *, bands: ProcessedWavelengthBands) -> sc.DataArray: + """ + Reduce data by summing or concatenating along the wavelength dimension. + + Parameters + ---------- + data: + Numerator or denominator data to be reduced. + wavelength_bands: + Defines bands in wavelength that can be used to separate different wavelength + ranges that contribute to different regions in Q space. Note that this needs to + be defined, so if all wavelengths should be used, this should simply be a start + and end edges that encompass the entire wavelength range. + + Returns + ------- + : + Q-dependent data, ready for normalization. + """ + wav = 'wavelength' + if part.bins is not None: + # 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)) + parts = [ + _do_reduce(part[wav, wav_range[0] : wav_range[1]]) + for wav_range in sc.collapse(bands, keep=wav).values() + ] + band_dim = (set(bands.dims) - {'wavelength'}).pop() + reduced = parts[0] if len(parts) == 1 else sc.concat(parts, band_dim) + return reduced.assign_coords(wavelength=bands.squeeze()) + + +def reduce_q( + data: WavelengthScaledQ[ScatteringRunType, IofQPart], + bands: ProcessedWavelengthBands, +) -> ReducedQ[ScatteringRunType, IofQPart]: + return ReducedQ[ScatteringRunType, IofQPart](_reduce(data, bands=bands)) + + +def reduce_qxy( + data: WavelengthScaledQxy[ScatteringRunType, IofQPart], + bands: ProcessedWavelengthBands, +) -> ReducedQxy[ScatteringRunType, IofQPart]: + return ReducedQxy[ScatteringRunType, IofQPart](_reduce(data, bands=bands)) + + def normalize_q( - numerator: CleanSummedQ[ScatteringRunType, Numerator], - denominator: CleanSummedQ[ScatteringRunType, Denominator], + numerator: ReducedQ[ScatteringRunType, Numerator], + denominator: ReducedQ[ScatteringRunType, Denominator], return_events: ReturnEvents, uncertainties: UncertaintyBroadcastMode, - wavelength_bands: ProcessedWavelengthBands, ) -> IofQ[ScatteringRunType]: return IofQ[ScatteringRunType]( _normalize( @@ -427,17 +431,15 @@ def normalize_q( denominator=denominator, return_events=return_events, uncertainties=uncertainties, - wavelength_bands=wavelength_bands, ) ) def normalize_qxy( - numerator: CleanSummedQxy[ScatteringRunType, Numerator], - denominator: CleanSummedQxy[ScatteringRunType, Denominator], + numerator: ReducedQxy[ScatteringRunType, Numerator], + denominator: ReducedQxy[ScatteringRunType, Denominator], return_events: ReturnEvents, uncertainties: UncertaintyBroadcastMode, - wavelength_bands: ProcessedWavelengthBands, ) -> IofQxy[ScatteringRunType]: return IofQxy[ScatteringRunType]( _normalize( @@ -445,19 +447,22 @@ def normalize_qxy( denominator=denominator, return_events=return_events, uncertainties=uncertainties, - wavelength_bands=wavelength_bands, ) ) +reduce_q.__doc__ = _reduce.__doc__ +reduce_qxy.__doc__ = _reduce.__doc__ normalize_q.__doc__ = _normalize.__doc__ normalize_qxy.__doc__ = _normalize.__doc__ providers = ( transmission_fraction, - iofq_norm_wavelength_term, - iofq_denominator, + norm_monitor_term, + norm_detector_term, + reduce_q, + reduce_qxy, normalize_q, normalize_qxy, process_wavelength_bands, diff --git a/src/ess/sans/types.py b/src/ess/sans/types.py index 86fb8019..392d13a3 100644 --- a/src/ess/sans/types.py +++ b/src/ess/sans/types.py @@ -154,8 +154,8 @@ class TransmissionFraction( """Transmission fraction""" -CleanDirectBeam = NewType('CleanDirectBeam', sc.DataArray | None) -"""Direct beam after resampling to required wavelength bins""" +CleanDirectBeam = NewType('CleanDirectBeam', sc.DataArray) +"""Direct beam after resampling to required wavelength bins, else and array of ones.""" class DetectorPixelShape(sciline.Scope[ScatteringRunType, sc.DataGroup], sc.DataGroup): @@ -192,8 +192,8 @@ class MaskedData(sciline.Scope[ScatteringRunType, sc.DataArray], sc.DataArray): """Raw data with pixel-specific masks applied""" -class NormWavelengthTerm(sciline.Scope[ScatteringRunType, sc.DataArray], sc.DataArray): - """Normalization term (numerator) for IofQ before scaling with solid-angle.""" +class MonitorTerm(sciline.Scope[ScatteringRunType, sc.DataArray], sc.DataArray): + """Monitor-dependent factor of the Normalization term (numerator) for IofQ.""" class CleanWavelength( @@ -208,10 +208,16 @@ class CleanWavelength( """ -class CleanWavelengthMasked( +class WavelengthScaledQ( sciline.ScopeTwoParams[ScatteringRunType, IofQPart, sc.DataArray], sc.DataArray ): - """Result of applying wavelength masking to :py:class:`CleanWavelength`""" + """Result of applying wavelength scaling/masking to :py:class:`CleanSummedQ`""" + + +class WavelengthScaledQxy( + sciline.ScopeTwoParams[ScatteringRunType, IofQPart, sc.DataArray], sc.DataArray +): + """Result of applying wavelength scaling/masking to :py:class:`CleanSummedQxy`""" class CleanQ( @@ -239,6 +245,18 @@ class CleanSummedQxy( Qy bins""" +class ReducedQ( + sciline.ScopeTwoParams[ScatteringRunType, IofQPart, sc.DataArray], sc.DataArray +): + """Result of reducing :py:class:`CleanSummedQ` over the wavelength dimensions""" + + +class ReducedQxy( + sciline.ScopeTwoParams[ScatteringRunType, IofQPart, sc.DataArray], sc.DataArray +): + """Result of reducing :py:class:`CleanSummedQxy` over the wavelength dimensions""" + + class IofQ(sciline.Scope[ScatteringRunType, sc.DataArray], sc.DataArray): """I(Q)""" diff --git a/tests/loki/common.py b/tests/loki/common.py index 4d4002cb..5ecd2b91 100644 --- a/tests/loki/common.py +++ b/tests/loki/common.py @@ -3,6 +3,7 @@ import sciline import scipp as sc +import ess.loki.data # noqa: F401 from ess import loki from ess.sans.types import ( BackgroundRun, diff --git a/tests/loki/iofq_test.py b/tests/loki/iofq_test.py index f401a923..e979ce07 100644 --- a/tests/loki/iofq_test.py +++ b/tests/loki/iofq_test.py @@ -17,7 +17,7 @@ BackgroundSubtractedIofQxy, BeamCenter, CleanSummedQ, - CleanWavelengthMasked, + CleanWavelength, CorrectForGravity, Denominator, DimsToKeep, @@ -280,13 +280,13 @@ def test_phi_with_gravity(): pipeline = make_workflow() pipeline[BeamCenter] = _compute_beam_center() pipeline[CorrectForGravity] = False - data_no_grav = pipeline.compute( - CleanWavelengthMasked[SampleRun, Numerator] - ).flatten(to='pixel') + data_no_grav = pipeline.compute(CleanWavelength[SampleRun, Numerator]).flatten( + to='pixel' + ) graph_no_grav = pipeline.compute(ElasticCoordTransformGraph) pipeline[CorrectForGravity] = True data_with_grav = ( - pipeline.compute(CleanWavelengthMasked[SampleRun, Numerator]) + pipeline.compute(CleanWavelength[SampleRun, Numerator]) .flatten(to='pixel') .hist(wavelength=sc.linspace('wavelength', 1.0, 12.0, 101, unit='angstrom')) )