Skip time slices for which the transmission calculation fails

docs/index.rst

@@ -13,9 +13,10 @@ User Guide
    :maxdepth: 2
 
    /user/calibration
-   /user/reduction_output
+   /user/corrections/index
    /user/slicing
    /user/binning
+   /user/reduction_output
    /drtsans/reduction_parameters
    /drtsans/reduction_scripts
    release_notes

docs/release_notes.rst

@@ -23,6 +23,7 @@ Release Notes
 
 **Of interest to the User**:
 
+- PR #993: Skip slices with too high transmission error when using time sliced sample transmission run
 - PR #325: Migrates repository from GitLab to GitHub
 - MR #1185: Rename output directory for coherent and incoherent-inelastic corrections as "info"
 - MR #1184: Document wedge binning
@@ -37,7 +38,7 @@ Release Notes
 - MR #1167: Allow separate configurations for inelastic incoherence correction per frame in frame skipping mode
 - MR #1166: Added option to _fitSpectrum to auto-find one wedge and mirror it
 - MR #1162: When reducing `gpsans` data with `direct_beam` scaling, the `direct_beam_scaling` parameter is now logged during
-the reduction process and stored in the output Nexus file at `reduction_information/special_parameters/direct_beam_scaling/value`.
+  the reduction process and stored in the output Nexus file at `reduction_information/special_parameters/direct_beam_scaling/value`.
 - MR #1161: Add a parameters removeAlgorithmHistory to write less data and speed up I/O during reduction
 - MR #1160: Expose pixel detector rescalings to the instrument API's
 - MR #1159: Separate configuration for elastic normalization and inelastic incoherence correction

docs/user/corrections/index.rst

@@ -0,0 +1,9 @@
+.. _user.corrections.index:
+
+Corrections
+===========
+
+.. toctree::
+   :maxdepth: 2
+
+   /user/corrections/transmission

docs/user/corrections/transmission.rst

@@ -0,0 +1,112 @@
+.. _user.corrections.transmission:
+
+Transmission
+============
+
+The measured intensity must be corrected to account for sample transmission:
+
+.. math::
+
+    I'_{sample}(x,y,\lambda) = \frac{I_{sample}(x,y,\lambda)}{T(x,y,\lambda)}
+
+The transmission correction is calculated using a transmission run and an empty beam transmission
+run (reference).
+It is possible to specify the transmission value directly in the transmission parameter ``"value"``,
+however, this is mostly used for diagnostic purposes.
+
+.. code-block:: json
+
+    {
+      "sample": {
+        "runNumber": 10010,
+        "loadOptions": {},
+        "thickness": 1.0,
+        "transmission": {
+          "runNumber": null,
+          "value": null,
+          "errorTolerance": null
+        }
+      },
+      "background": {
+        "runNumber": null,
+        "transmission": {
+          "runNumber": null,
+          "value": null
+        }
+      },
+      "emptyTransmission": {
+        "runNumber": null,
+        "value": null
+      },
+      "mmRadiusForTransmission": null,
+      "useTimeSliceTransmission": false,
+      "useThetaDepTransCorrection": true
+    }
+
+Time Slice Transmission (Bio-SANS only)
+---------------------------------------
+
+When using time slicing (``"useTimeSlice": true``), users can optionally calculate the transmission
+correction using the time sliced sample run by setting the parameter
+``"useTimeSliceTransmission": true``. The sample transmission run number is ignored when
+``"useTimeSliceTransmission": true``. The time slice transmission option can be used when the sample
+transmission is expected to change over time.
+
+Time slices for which the transmission calculation fails will be skipped. The transmission
+calculation can fail due to all transmission values being NaN or if the transmission error is
+higher than the allowed relative transmission error, which is configurable in the sample
+transmission parameter ``"errorTolerance"`` (default: 0.01). For example, the last time slice may
+be shorter and, therefore, include fewer counts, resulting in large statistical errors in the
+transmission calculation.
+
+.. code-block:: json
+
+    {
+      "sample": {
+        "runNumber": 10010,
+        "loadOptions": {},
+        "thickness": 1.0,
+        "transmission": {
+          "runNumber": null,
+          "value": null,
+          "errorTolerance": 0.05
+        }
+      },
+      "emptyTransmission": {
+        "runNumber": 10005,
+        "value": null
+      },
+      "useTimeSlice": true,
+      "timeSliceInterval": 100.0,
+      "useTimeSliceTransmission": true
+    }
+
+Parameters
+----------
+
+.. list-table::
+   :widths: 25 65 10
+   :header-rows: 1
+
+   * - Parameter
+     - Description
+     - Default
+   * - ``"mmRadiusForTransmission"``
+     - Beam radius within which the transmission will be calculated. If ``null``, then the beam
+       radius is calculated from the sample logs.
+     - ``null``
+   * - ``"useThetaDepTransCorrection"``
+     - If ``true``, a theta dependent transmission correction will be applied, which takes into
+       account the effect of the scattering angle on the transmission.
+     - ``true``
+   * - ``"useTimeSliceTransmission"``
+     - (`Only for Bio-SANS and when` ``"useTimeSlice": true``.) If ``true``, the transmission
+       correction will be calculated using the time sliced sample run itself instead of a separate
+       sample transmission run. This is useful when the sample transmission is expected to change
+       over time. Slices with relative transmission error larger than
+       ``"transmissionErrorTolerance"`` will be skipped.
+     - ``false``
+   * - ``"errorTolerance"``
+     - (`Only for Bio-SANS and when` ``"useTimeSlice": true`` `and` ``"useTimeSliceTransmission": true``.)
+       Maximum relative transmission error.
+     - 0.01

docs/user/reduction_output.rst

@@ -5,8 +5,8 @@ Reduction Output
 
 
 
-CANSAS
-------
+canSAS format
+-------------
 
 Collective Action for Nomadic Small-Angle Scatterers (canSAS) provides standards and tools for the small-angle scattering user community. See more at https://www.cansas.org
 

src/drtsans/configuration/schema/common.json

@@ -197,9 +197,13 @@
           },
           "value": {
             "$ref": "common.json#/definitions/transmissionValueTypes"
+          },
+          "errorTolerance": {
+            "$ref": "common.json#/definitions/transmissionValueTypes",
+            "description": "Maximum allowed relative error in the calculated transmission"
           }
         },
-        "maxProperties": 2,
+        "maxProperties": 3,
         "required": ["runNumber", "value"],
         "description": "The transmission for the sample, run number or value (0 < value <=1). Can be empty",
         "examples": ["0.9", 1.0]

src/drtsans/mono/biosans/api.py

@@ -31,6 +31,7 @@
 from drtsans.stitch import stitch_binned_profiles
 from drtsans.reductionlog import savereductionlog
 from drtsans.thickness_normalization import normalize_by_thickness
+from drtsans.transmission import TransmissionErrorToleranceError, TransmissionNanError
 
 # third party imports
 from mantid.dataobjects import Workspace2D
@@ -67,6 +68,9 @@
 SI_WINDOW_NOMINAL_DISTANCE_METER = 0.071
 SAMPLE_SI_META_NAME = "CG3:CS:SampleToSi"
 
+# default transmission error tolerance
+DEFAULT_TRANSMISSION_ERROR_TOLERANCE = 0.01
+
 # setup logger
 # NOTE: If logging information is not showing up, please check the mantid log level.
 #       If problem persists, please visit:
@@ -1205,6 +1209,14 @@ def reduce_single_configuration(
     absolute_scale = reduction_config["StandardAbsoluteScale"]
     time_slice = reduction_config["useTimeSlice"]
     time_slice_transmission = reduction_config["useTimeSlice"] and reduction_config["useTimeSliceTransmission"]
+    # Transmission tolerance errors are only checked when slicing the transmission run,
+    # because this is the only observed use-case of poor statistics for transmission calculations
+    if time_slice_transmission:
+        sample_trans_error_tol = reduction_input["sample"]["transmission"]["errorTolerance"]
+        if sample_trans_error_tol is None:
+            sample_trans_error_tol = DEFAULT_TRANSMISSION_ERROR_TOLERANCE
+    else:
+        sample_trans_error_tol = None
     output_dir = reduction_config["outputDir"]
 
     nxbins_main = reduction_config["numMainQxQyBins"]
@@ -1369,6 +1381,7 @@ def _prepare_sample_transmission_ws(_sample_transmission):
             empty_trans_ws,
             radius=transmission_radius,
             radius_unit="mm",
+            transmission_error_tolerance=sample_trans_error_tol,
         )
 
     sample_trans_ws = None
@@ -1391,10 +1404,14 @@ def _prepare_sample_transmission_ws(_sample_transmission):
         if len(loaded_ws.sample) > 1:
             output_suffix = f"_{i}"
 
-        try:
-            if time_slice_transmission:
+        if time_slice_transmission:
+            try:
                 _, sample_trans_ws = _prepare_sample_transmission_ws(raw_sample_ws)
+            except (ZeroMonitorCountsError, TransmissionErrorToleranceError, TransmissionNanError) as e:
+                logger.warning(f"Skipping slice {sample_name}: {e}.")
+                continue
 
+        try:
             processed_data_main, trans_main = process_single_configuration(
                 raw_sample_ws,
                 sample_trans_ws=sample_trans_ws,

src/drtsans/mono/normalization.py

@@ -103,8 +103,10 @@ def normalize_by_monitor(input_workspace, output_workspace=None):
 
     Raises
     ------
-    RuntimeError
+    NoMonitorMetadataError
         No monitor metadata found in the sample logs of the input workspace
+    ZeroMonitorCountsError
+        No monitor counts in the input workspace
     """
     metadata_entry_names = [
         "monitor",  # created by load_events

src/drtsans/mono/transmission.py

@@ -7,7 +7,47 @@
 __all__ = ["calculate_transmission", "apply_transmission_correction"]
 
 
-def calculate_transmission(input_sample, input_reference, radius=None, radius_unit="mm", output_workspace=None):
+def calculate_transmission(
+    input_sample,
+    input_reference,
+    radius=None,
+    radius_unit="mm",
+    transmission_error_tolerance=None,
+    output_workspace=None,
+):
+    """
+    Calculate the raw transmission coefficients at zero scattering angle
+    from already prepared sample and reference data.
+
+    Parameters
+    ----------
+    input_sample: str, ~mantid.api.MatrixWorkspace, ~mantid.api.IEventWorkspace
+        Prepared sample workspace (possibly obtained with an attenuated beam)
+    input_reference: str, ~mantid.api.MatrixWorkspace, ~mantid.api.IEventWorkspace
+        Prepared direct beam workspace (possibly obtained with an attenuated beam)
+    radius: float
+        Radius around the beam center for pixel integration, in millimeters.
+        If None, radius will be obtained or calculated using `input_reference` workspace.
+    radius_unit: str
+        Either 'mm' or 'm', and only used in conjunction with option `radius`.
+    transmission_error_tolerance: float | None
+        Maximum relative error for transmission. If None, the error will not be checked.
+    output_workspace: str
+        Name of the output workspace containing the raw transmission values.
+        If None, a hidden random name will be provided.
+
+    Returns
+    -------
+    ~mantid.api.MatrixWorkspace
+        Workspace containing the raw transmission values
+
+    Raises
+    ------
+    TransmissionNanError
+        If all transmission values are NaN
+    TransmissionToleranceError
+        If there is insufficient statistics to calculate the transmission correction
+    """
     if radius is None:
         logger.information("Calculating beam radius from sample logs")
         radius = beam_radius(input_reference, unit="mm")
@@ -17,6 +57,7 @@ def calculate_transmission(input_sample, input_reference, radius=None, radius_un
         input_reference,
         radius,
         radius_unit=radius_unit,
+        transmission_error_tolerance=transmission_error_tolerance,
         output_workspace=output_workspace,
     )
 

src/drtsans/transmission.py

@@ -25,10 +25,29 @@
 from drtsans.mask_utils import circular_mask_from_beam_center, masked_detectors
 
 # Symbols to be exported
-__all__ = ["apply_transmission_correction", "calculate_transmission"]
+__all__ = [
+    "apply_transmission_correction",
+    "calculate_transmission",
+    "TransmissionErrorToleranceError",
+    "TransmissionNanError",
+]
 
 
-def calculate_transmission(input_sample, input_reference, radius, radius_unit="mm", output_workspace=None):
+class TransmissionErrorToleranceError(Exception):
+    """Exception raised when the transmission error is larger than the transmission error tolerance"""
+
+    pass
+
+
+class TransmissionNanError(Exception):
+    """Exception raised when all transmission values are NaN"""
+
+    pass
+
+
+def calculate_transmission(
+    input_sample, input_reference, radius, radius_unit="mm", transmission_error_tolerance=None, output_workspace=None
+):
     """
     Calculate the raw transmission coefficients at zero scattering angle
     from already prepared sample and reference data.
@@ -49,10 +68,12 @@ def calculate_transmission(input_sample, input_reference, radius, radius_unit="m
     input_reference: str, ~mantid.api.MatrixWorkspace, ~mantid.api.IEventWorkspace
         Prepared direct beam workspace (possibly obtained with an attenuated beam)
     radius: float
-        Radius around the bean center for pixel integration, in millimeters.
+        Radius around the beam center for pixel integration, in millimeters.
         If None, radius will be obtained or calculated using `input_reference` workspace.
     radius_unit: str
         Either 'mm' or 'm', and only used in conjunction with option `radius`.
+    transmission_error_tolerance: float | None
+        Maximum relative error for transmission. If None, the error will not be checked.
     output_workspace: str
         Name of the output workspace containing the raw transmission values.
         If None, a hidden random name will be provided.
@@ -61,6 +82,13 @@ def calculate_transmission(input_sample, input_reference, radius, radius_unit="m
     -------
     ~mantid.api.MatrixWorkspace
         Workspace containing the raw transmission values
+
+    Raises
+    ------
+    TransmissionNanError
+        If all transmission values are NaN.
+    TransmissionToleranceError
+        If there are insufficient statistics to calculate the transmission correction.
     """
     if output_workspace is None:
         output_workspace = mtd.unique_hidden_name()
@@ -124,10 +152,28 @@ def calculate_transmission(input_sample, input_reference, radius, radius_unit="m
     # Notify of incorrect calculation of zero angle transmission
     # Will happen if the beam centers have been totally masked
     if bool(np.all(np.isnan(zero_angle_transmission_workspace.readY(0)))) is True:
-        raise RuntimeError("Transmission at zero-angle is NaN")
+        raise TransmissionNanError("Transmission at zero-angle is NaN")
 
-    # Notify of average transmission value
     non_gap_indexes = np.isfinite(zero_angle_transmission_workspace.readY(0))
+
+    if transmission_error_tolerance:
+        # Verify that errors are below the transmission error tolerance
+        transmission_relative_error = (
+            zero_angle_transmission_workspace.readE(0)[non_gap_indexes]
+            / zero_angle_transmission_workspace.readY(0)[non_gap_indexes]
+        )
+        if not np.all(transmission_relative_error < transmission_error_tolerance):
+            i_max = np.argmax(transmission_relative_error)
+            rel_error = transmission_relative_error[i_max]
+            transmission = zero_angle_transmission_workspace.readY(0)[non_gap_indexes][i_max]
+            abs_error = zero_angle_transmission_workspace.readE(0)[non_gap_indexes][i_max]
+            raise TransmissionErrorToleranceError(
+                f"transmission_error / transmission_value ({abs_error:.4f} / {transmission:.4f} = "
+                f"{rel_error:.4f}) > transmission_relative_error_tolerance "
+                f"({transmission_error_tolerance:.4f})"
+            )
+
+    # Notify of average transmission value
     average_zero_angle_transmission = np.mean(zero_angle_transmission_workspace.readY(0)[non_gap_indexes])
     average_zero_angle_transmission_error = np.linalg.norm(zero_angle_transmission_workspace.readE(0)[non_gap_indexes])
     message = "Average zero angle transmission = {0} +/- {1}".format(