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,109 @@
+.. _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(\lambda,x,y)}
+
+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": None,
+        "loadOptions": {},
+        "thickness": 1.0,
+        "transmission": {
+          "runNumber": None,
+          "value": None,
+          "errorTolerance": None
+        }
+      },
+      "background": {
+        "runNumber": None,
+        "transmission": {
+          "runNumber": None,
+          "value": None
+        }
+      },
+      "emptyTransmission": {
+        "runNumber": None,
+        "value": None
+      },
+      "mmRadiusForTransmission": None,
+      "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).
+
+.. code-block:: json
+
+    {
+      "sample": {
+        "runNumber": 10010,
+        "loadOptions": {},
+        "thickness": 1.0,
+        "transmission": {
+          "runNumber": None,
+          "value": None,
+          "errorTolerance": 0.05
+        }
+      },
+      "emptyTransmission": {
+        "runNumber": 10005,
+        "value": None
+      },
+      "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 ``None``, then the beam
+       radius is calculated from the sample logs.
+     - ``None``
+   * - ``"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
@@ -1205,6 +1206,9 @@ def reduce_single_configuration(
     absolute_scale = reduction_config["StandardAbsoluteScale"]
     time_slice = reduction_config["useTimeSlice"]
     time_slice_transmission = reduction_config["useTimeSlice"] and reduction_config["useTimeSliceTransmission"]
+    sample_trans_error_tol = None
+    if time_slice_transmission:
+        sample_trans_error_tol = reduction_input["sample"]["transmission"]["errorTolerance"]
     output_dir = reduction_config["outputDir"]
 
     nxbins_main = reduction_config["numMainQxQyBins"]
@@ -1369,6 +1373,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
@@ -1488,7 +1493,7 @@ def _prepare_sample_transmission_ws(_sample_transmission):
                         "background": None,
                     },
                 )
-        except ZeroMonitorCountsError as e:
+        except (ZeroMonitorCountsError, TransmissionErrorToleranceError, TransmissionNanError) as e:
             if time_slice:
                 logger.warning(f"Skipping slice {sample_name}: {e}.")
                 continue

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
+    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
+        Maximum relative error for transmission
     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,11 @@ def calculate_transmission(input_sample, input_reference, radius, radius_unit="m
     -------
     ~mantid.api.MatrixWorkspace
         Workspace containing the raw transmission values
+
+    Raises
+    ------
+    TransmissionToleranceError
+        If there is insufficient statistics to calculate the transmission correction
     """
     if output_workspace is None:
         output_workspace = mtd.unique_hidden_name()
@@ -124,10 +150,26 @@ 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)
+            error_max = transmission_relative_error[i_max]
+            error_max_transmission = zero_angle_transmission_workspace.readY(0)[non_gap_indexes][i_max]
+            raise TransmissionErrorToleranceError(
+                f"Transmission error {error_max} > transmission error tolerance "
+                f"{transmission_error_tolerance} (transmission {error_max_transmission})"
+            )
+
+    # 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(

tests/integration/drtsans/mono/biosans/test_api.py

@@ -320,7 +320,7 @@ def _load_synthetic_run(keyword: str) -> Workspace2D:
 @pytest.mark.datarepo
 def test_reduce_single_configuration_slice_transmission_false(datarepo_dir, temp_directory):
     reduction_input = REDUCTION_INPUT.copy()
-
+    reduction_input = reduction_parameters(parameters_particular=reduction_input, validate=False)
     reduction_input["configuration"]["outputDir"] = temp_directory(prefix="trans_slice_false")
 
     prefix = "sans-backend-test" + str(threading.get_ident()) + "_"
@@ -359,20 +359,25 @@ def test_reduce_single_configuration_slice_transmission_true(datarepo_dir, temp_
     reduction_input = REDUCTION_INPUT.copy()
     reduction_input["configuration"]["useTimeSlice"] = True
     reduction_input["configuration"]["useTimeSliceTransmission"] = True
+    reduction_input["sample"]["transmission"]["errorTolerance"] = 0.05
     reduction_input = reduction_parameters(parameters_particular=reduction_input, validate=False)
     reduction_input["configuration"]["outputDir"] = temp_directory(prefix="trans_slice_true")
     prefix = "sans-backend-test" + str(threading.get_ident()) + "_"
     with amend_config(data_dir=datarepo_dir.biosans):
         loaded = load_all_files(reduction_input, prefix)
     assert len(loaded.sample) == 11  # 600.06/60 = 11 rounded up
     output = reduce_single_configuration(loaded, reduction_input)
-    assert len(output) == 11  # expect 11 out since nothing should be skipped
+    assert len(output) == 10  # the 11th slice is skipped due to too high transmission error
 
     # just need a couple components from reduce
     # but the whole thing needs to be run then a few components pulled
+    # The sample transmission workspace is for the last slice, which is skipped due to too large error
     transmission = calculate_transmission(
         mtd["_sample_trans"],  # pull relevant transmission
         mtd["_empty"],  # pull relevant
+        # The sample transmission workspace in the MDS is for the last (skipped) slice, therefore,
+        # tolerate a high transmission error. The transmission value is OK to compare below.
+        transmission_error_tolerance=3.0,
     )
 
     transmission_val = transmission.extractY()[0][0]