Support event-mode I(Q)

docs/examples/loki-direct-beam.ipynb

@@ -98,6 +98,7 @@
     "\n",
     "params[CorrectForGravity] = True\n",
     "params[UncertaintyBroadcastMode] = UncertaintyBroadcastMode.upper_bound\n",
+    "params[sans.ReturnEvents] = False\n",
     "\n",
     "params[QBins] = sc.linspace(dim='Q', start=0.01, stop=0.3, num=101, unit='1/angstrom')"
    ]
@@ -413,7 +414,8 @@
    "mimetype": "text/x-python",
    "name": "python",
    "nbconvert_exporter": "python",
-   "pygments_lexer": "ipython3"
+   "pygments_lexer": "ipython3",
+   "version": "3.10.12"
   }
  },
  "nbformat": 4,

docs/examples/sans2d.ipynb

@@ -98,7 +98,8 @@
     "    dims=['wavelength'], values=[0.7, 17.1], unit='angstrom'\n",
     ")\n",
     "pipeline[CorrectForGravity] = True\n",
-    "pipeline[UncertaintyBroadcastMode] = UncertaintyBroadcastMode.upper_bound"
+    "pipeline[UncertaintyBroadcastMode] = UncertaintyBroadcastMode.upper_bound\n",
+    "pipeline[sans.ReturnEvents] = True"
    ]
   },
   {
@@ -158,7 +159,25 @@
    "outputs": [],
    "source": [
     "result = iofq.compute()\n",
-    "result.plot()"
+    "result.hist().plot(scale={'Q': 'log'}, norm='log')"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "28532aa7",
+   "metadata": {},
+   "source": [
+    "As the result was computed in event-mode, we can also use a different $Q$-binning, without re-reducing the data:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "e0748a1a",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "result.hist(Q=60).plot(scale={'Q': 'log'}, norm='log')"
    ]
   },
   {
@@ -182,7 +201,7 @@
     "result_drop = pipeline.compute(BackgroundSubtractedIofQ)\n",
     "# Reset the UnsertaintyBroadcastMode to the old value\n",
     "pipeline[UncertaintyBroadcastMode] = UncertaintyBroadcastMode.upper_bound\n",
-    "sc.DataGroup(upper_bound=result, dropped=result_drop).plot(norm='log')"
+    "sc.DataGroup(upper_bound=result, dropped=result_drop).hist().plot(norm='log')"
    ]
   },
   {
@@ -291,7 +310,7 @@
     ")\n",
     "pipeline.set_param_table(param_table)\n",
     "results = pipeline.compute(sciline.Series[Mode, BackgroundSubtractedIofQ])\n",
-    "sc.DataGroup(results).plot(norm='log')"
+    "sc.DataGroup(results).hist().plot(norm='log')"
    ]
   },
   {
@@ -357,7 +376,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "pp.plot(sc.collapse(result, keep='Q'), norm='log')"
+    "pp.plot(sc.collapse(result.hist(), keep='Q'), norm='log')"
    ]
   }
  ],
@@ -377,7 +396,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.10.13"
+   "version": "3.10.12"
   }
  },
  "nbformat": 4,

docs/examples/zoom.ipynb

@@ -109,7 +109,8 @@
     "    dims=['wavelength'], values=[0.7, 17.1], unit='angstrom'\n",
     ")\n",
     "params[CorrectForGravity] = True\n",
-    "params[UncertaintyBroadcastMode] = UncertaintyBroadcastMode.upper_bound"
+    "params[UncertaintyBroadcastMode] = UncertaintyBroadcastMode.upper_bound\n",
+    "params[sans.ReturnEvents] = False"
    ]
   },
   {

src/esssans/__init__.py

@@ -21,6 +21,7 @@
 )
 from .common import transmission_from_background_run, transmission_from_sample_run
 from .direct_beam import direct_beam
+from .types import BackgroundSubtractedIofQ, IofQ, ReturnEvents, SampleRun
 
 providers = (
     *conversions.providers,
@@ -41,6 +42,7 @@
 del importlib
 
 __all__ = [
+    'ReturnEvents',
     'beam_center_finder',
     'common',
     'conversions',

src/esssans/beam_center_finder.py

@@ -27,6 +27,7 @@
     MaskedData,
     NormWavelengthTerm,
     QBins,
+    ReturnEvents,
     SampleRun,
     UncertaintyBroadcastMode,
     WavelengthBins,
@@ -176,6 +177,7 @@ def _iofq_in_quadrants(
     ]
     params = {}
     params[UncertaintyBroadcastMode] = UncertaintyBroadcastMode.upper_bound
+    params[ReturnEvents] = False
     params[WavelengthBins] = wavelength_bins
     params[QBins] = q_bins
     params[DetectorPixelShape[SampleRun]] = pixel_shape

src/esssans/direct_beam.py

@@ -132,7 +132,7 @@ def direct_beam(pipeline: Pipeline, I0: sc.Variable, niter: int = 5) -> List[dic
         # parameters, nor given by any providers, so it will be considered flat.
         # TODO: Should we have a check that DirectBeam cannot be computed from the
         # pipeline?
-        iofq = pipeline.compute(BackgroundSubtractedIofQ)
+        iofq = sc.values(pipeline.compute(BackgroundSubtractedIofQ))
         iofq_full = iofq['band', -1]
         iofq_bands = iofq['band', :-1]
 

src/esssans/i_of_q.py

@@ -24,6 +24,7 @@
     MonitorType,
     NonBackgroundWavelengthRange,
     QBins,
+    ReturnEvents,
     RunType,
     SampleRun,
     UncertaintyBroadcastMode,
@@ -314,8 +315,12 @@ def _dense_merge_spectra(
 
 
 def subtract_background(
-    sample: IofQ[SampleRun], background: IofQ[BackgroundRun]
+    sample: IofQ[SampleRun],
+    background: IofQ[BackgroundRun],
+    return_events: ReturnEvents,
 ) -> BackgroundSubtractedIofQ:
+    if return_events and sample.bins is not None and background.bins is not None:
+        return sample.bins.concatenate(-background)
     if sample.bins is not None:
         sample = sample.bins.sum()
     if background.bins is not None:

src/esssans/io.py

@@ -15,8 +15,11 @@ def save_background_subtracted_iofq(
     run_number: RunNumber,
     run_title: RunTitle,
 ) -> None:
-    """Save background subtracted IofQ as an NXcanSAS file."""
-    da = iofq.copy(deep=False)
+    """Save background-subtracted I(Q) histogram as an NXcanSAS file."""
+    if iofq.bins is None:
+        da = iofq.copy(deep=False)
+    else:
+        da = iofq.hist()
     if da.coords.is_edges('Q'):
         da.coords['Q'] = sc.midpoints(da.coords['Q'])
     with snx.File(out_filename, 'w') as f:

src/esssans/loki/io.py

@@ -37,6 +37,10 @@ def _patch_data(
     da: sc.DataArray, sample_position: sc.Variable, source_position: sc.Variable
 ) -> sc.DataArray:
     out = da.copy(deep=False)
+    if out.bins is not None:
+        content = out.bins.constituents['data']
+        if content.variances is None:
+            content.variances = content.values
     out.coords['sample_position'] = sample_position
     out.coords['source_position'] = source_position
     out.coords['gravity'] = gravity_vector()

src/esssans/normalization.py

@@ -20,6 +20,7 @@
     LabFrameTransform,
     NormWavelengthTerm,
     Numerator,
+    ReturnEvents,
     RunType,
     SolidAngle,
     Transmission,
@@ -28,6 +29,7 @@
     UncertaintyBroadcastMode,
 )
 from .uncertainty import (
+    broadcast_to_events_with_upper_bound_variances,
     broadcast_with_upper_bound_variances,
     drop_variances_if_broadcast,
 )
@@ -298,6 +300,8 @@ def iofq_denominator(
 def normalize(
     numerator: CleanSummedQ[RunType, Numerator],
     denominator: CleanSummedQ[RunType, Denominator],
+    return_events: ReturnEvents,
+    uncertainties: UncertaintyBroadcastMode,
 ) -> IofQ[RunType]:
     """
     Perform normalization of counts as a function of Q.
@@ -312,17 +316,27 @@ def normalize(
     denominator:
         The divisor for the normalization operation. This cannot be event data, it must
         contain histogrammed data.
+    return_events:
+        Whether to return the result as event data or histogrammed data.
 
     Returns
     -------
     :
         The input data normalized by the supplied denominator.
     """
-    if denominator.variances is not None and numerator.bins is not None:
-        # Event-mode normalization is not correct of norm-term has variances.
+    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.
+        if denominator.variances is not None:
+            if uncertainties == UncertaintyBroadcastMode.drop:
+                denominator = sc.values(denominator)
+            else:
+                denominator = broadcast_to_events_with_upper_bound_variances(
+                    denominator, events=numerator
+                )
+    else:
         numerator = numerator.hist()
-    if numerator.bins is not None:
+    if numerator.bins is not None and denominator.bins is None:
         da = numerator.bins / sc.lookup(func=denominator, dim='Q')
     else:
         da = numerator / denominator

src/esssans/types.py

@@ -82,6 +82,9 @@ class TransmissionRun(sciline.Scope[RunType, int], int):
 See https://doi.org/10.3233/JNR-220049 for context.
 """
 
+ReturnEvents = NewType('ReturnEvents', bool)
+"""Whether to return events in the output I(Q)"""
+
 WavelengthBins = NewType('WavelengthBins', sc.Variable)
 """Wavelength binning"""
 

src/esssans/uncertainty.py

@@ -51,3 +51,32 @@ def _no_variance_broadcast(
     return (data.variances is None) or all(
         data.sizes.get(dim) == size for dim, size in sizes.items()
     )
+
+
+def broadcast_to_events_with_upper_bound_variances(
+    da: sc.DataArray, *, events: sc.DataArray
+) -> sc.DataArray:
+    """
+    Upper-bound estimate for errors from normalization in event-mode.
+
+    Count the number of events in each bin of the input data array. Then scale the
+    variances by the number of events in each bin. An explicit broadcast is performed
+    to bypass Scipp's safety check on broadcasting variances.
+
+    Details will be published in an upcoming publication by Simon Heybrock et al.
+    """
+    if da.variances is None:
+        return da
+    constituents = events.bins.constituents
+    Q = constituents['data'].coords['Q']
+    constituents['data'] = sc.DataArray(
+        sc.ones(sizes=Q.sizes, dtype='float32'), coords={'Q': Q}
+    )
+    # Combine all bins of the events that correspond to the same bin in the input data
+    to_concat = set(events.dims) - set(da.dims)
+    binned = sc.DataArray(sc.bins(**constituents).bins.concat(to_concat))
+    counts = binned.hist(Q=da.coords['Q'])
+    da = da.copy()
+    da.variances *= counts.values
+    da.data = sc.bins_like(events, da.data)
+    return da

tests/loki/common.py

@@ -13,6 +13,7 @@
     EmptyBeamRun,
     FileList,
     QBins,
+    ReturnEvents,
     SampleRun,
     TransmissionRun,
     UncertaintyBroadcastMode,
@@ -65,6 +66,7 @@ def make_params(
     params[BeamStopRadius] = sc.scalar(0.042, unit='m')
     params[CorrectForGravity] = True
     params[UncertaintyBroadcastMode] = UncertaintyBroadcastMode.upper_bound
+    params[ReturnEvents] = False
 
     params[QBins] = sc.linspace(
         dim='Q', start=0.01, stop=0.3, num=101, unit='1/angstrom'

tests/loki/iofq_test.py

@@ -13,6 +13,7 @@
     BackgroundSubtractedIofQ,
     BeamCenter,
     DimsToKeep,
+    ReturnEvents,
     UncertaintyBroadcastMode,
     WavelengthBands,
 )
@@ -39,6 +40,43 @@ def test_pipeline_can_compute_IofQ(uncertainties):
     assert result.dims == ('Q',)
 
 
+@pytest.mark.parametrize(
+    'uncertainties',
+    [UncertaintyBroadcastMode.drop, UncertaintyBroadcastMode.upper_bound],
+)
+@pytest.mark.parametrize(
+    'target', [sans.IofQ[sans.SampleRun], sans.BackgroundSubtractedIofQ]
+)
+def test_pipeline_can_compute_IofQ_in_event_mode(uncertainties, target):
+    params = make_params()
+    params[UncertaintyBroadcastMode] = uncertainties
+    pipeline = sciline.Pipeline(loki_providers(), params=params)
+    reference = pipeline.compute(target)
+    pipeline[ReturnEvents] = True
+    result = pipeline.compute(target)
+    assert result.bins is not None
+    assert sc.allclose(
+        sc.values(reference.data),
+        sc.values(result.hist().data),
+        rtol=sc.scalar(1e-11),
+        atol=sc.scalar(1e-11),
+    )
+    if uncertainties == UncertaintyBroadcastMode.drop:
+        assert sc.allclose(
+            sc.variances(reference).data,
+            sc.variances(result.hist()).data,
+            rtol=sc.scalar(1e-14),
+            atol=sc.scalar(1e-14),
+        )
+    else:
+        assert sc.allclose(
+            sc.variances(reference).data,
+            sc.variances(result.hist()).data,
+            rtol=sc.scalar(1e-9),
+            atol=sc.scalar(1e-9),
+        )
+
+
 def test_pipeline_can_compute_IofQ_in_wavelength_slices():
     params = make_params()
     band = np.linspace(1.0, 13.0, num=11)

tests/sans2d/reduction_test.py

@@ -22,6 +22,7 @@
     NonBackgroundWavelengthRange,
     QBins,
     RawData,
+    ReturnEvents,
     SampleRun,
     SolidAngle,
     TransmissionRun,
@@ -61,6 +62,7 @@ def make_params() -> dict:
     )
     params[CorrectForGravity] = True
     params[UncertaintyBroadcastMode] = UncertaintyBroadcastMode.upper_bound
+    params[ReturnEvents] = False
     return params
 
 
@@ -140,7 +142,7 @@ def test_uncertainty_broadcast_mode_drop_yields_smaller_variances():
     )
     params[UncertaintyBroadcastMode] = UncertaintyBroadcastMode.drop
     pipeline = sciline.Pipeline(sans2d_providers(), params=params)
-    drop = pipeline.compute(IofQ[SampleRun]).hist().data
+    drop = pipeline.compute(IofQ[SampleRun]).data
     params[UncertaintyBroadcastMode] = UncertaintyBroadcastMode.upper_bound
     pipeline = sciline.Pipeline(sans2d_providers(), params=params)
     upper_bound = pipeline.compute(IofQ[SampleRun]).data