Merge pull request #170 from scipp/drop_variances_in_sans

Drop normalization variances in sans and reflectometry workflows

Author: nvaytet

docs/instruments/external/powgen/powgen_reduction.ipynb

@@ -29,7 +29,6 @@
    },
    "outputs": [],
    "source": [
-    "%matplotlib widget\n",
     "import scipp as sc\n",
     "import scippneutron as scn\n",
     "import plopp as pp\n",
@@ -378,12 +377,53 @@
   },
   {
    "cell_type": "markdown",
-   "id": "a8f39532-e0ee-48cc-8c49-40ac7427ded8",
+   "id": "97ecaa8f-c905-4d02-aad6-dc7ef7c707a8",
    "metadata": {},
    "source": [
     "Currently, `load_and_preprocess_vanadium` uses a crude, preliminary event filtering mechanism which removed all events if the proton charge is too low.\n",
     "This happens with the empty instrument data in this case as can be seen in the log widget at the top of the notebook.\n",
     "\n",
+    "### Removing the variances of the Vanadium data"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "3d41c537-40b1-4ea7-87b1-f0dc3929ef55",
+   "metadata": {},
+   "source": [
+    "<div class=\"alert alert-warning\">\n",
+    "\n",
+    "**Warning**\n",
+    "    \n",
+    "Heybrock et al. (2023) have shown that Scipp's uncertainty propagation is not suited for broadcast operations\n",
+    "such as normalizing event counts by a scalar value, which is the case when normalizing by Vanadium counts.\n",
+    "These operations are forbidden in recent versions of Scipp.\n",
+    "Until an alternative method is found to satisfactorily track the variances in this workflow,\n",
+    "we remove the variances in the Vanadium data.\n",
+    "The issue is tracked [here](https://github.com/scipp/ess/issues/171).\n",
+    "\n",
+    "</div>"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "7ee60cf0-f2f5-426b-89e7-aa137a73b8c5",
+   "metadata": {
+    "tags": []
+   },
+   "outputs": [],
+   "source": [
+    "vana.bins.constituents['data'].variances = None"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "31a1525b-c5a2-4164-862d-5e2152d84c83",
+   "metadata": {},
+   "source": [
+    "### Conversion to d-spacing\n",
+    "\n",
     "Now, we compute d-spacing using the same calibration parameters as before."
    ]
   },
@@ -459,7 +499,9 @@
    "metadata": {},
    "source": [
     "We start simple by combining all spectra using `data.bins.concat('spectrum')`.\n",
-    "Then, we can normalize the same data by vanadium to get a d-spacing distribution."
+    "Then, we can normalize the same data by vanadium to get a d-spacing distribution.\n",
+    "\n",
+    "**Note that because we removed the variances on the Vanadium data, after the following cell, the standard deviations on the result are underestimated.**"
    ]
   },
   {
@@ -624,6 +666,7 @@
    },
    "outputs": [],
    "source": [
+    "%matplotlib widget\n",
     "pp.superplot(normalized)"
    ]
   }

docs/instruments/loki/sans2d_reduction.ipynb

@@ -105,7 +105,7 @@
    },
    "outputs": [],
    "source": [
-    "direct_beam = sc.io.open_hdf5(loki.data.get_path('DIRECT_SANS2D_REAR_34327_4m_8mm_16Feb16.hdf5'))\n",
+    "direct_beam = sc.io.load_hdf5(loki.data.get_path('DIRECT_SANS2D_REAR_34327_4m_8mm_16Feb16.hdf5'))\n",
     "sample = loki.io.load_sans2d(filename=loki.data.get_path('SANS2D00063114.hdf5'))\n",
     "direct = loki.io.load_sans2d(filename=loki.data.get_path('SANS2D00063091.hdf5'))\n",
     "background = loki.io.load_sans2d(filename=loki.data.get_path('SANS2D00063159.hdf5'))\n",

docs/instruments/loki/sans2d_to_I_of_Q.ipynb

@@ -117,7 +117,7 @@
     "# Direct measurement is with the empty sample holder/cuvette\n",
     "direct = loki.io.load_sans2d(filename=loki.data.get_path('SANS2D00063091.hdf5'))\n",
     "# Load direct beam function for main detector\n",
-    "direct_beam = sc.io.open_hdf5(loki.data.get_path('DIRECT_SANS2D_REAR_34327_4m_8mm_16Feb16.hdf5'))\n",
+    "direct_beam = sc.io.load_hdf5(loki.data.get_path('DIRECT_SANS2D_REAR_34327_4m_8mm_16Feb16.hdf5'))\n",
     "# Inspect sample data\n",
     "sample"
    ]
@@ -127,7 +127,7 @@
    "id": "c4569a53-a19a-440e-8aa0-63df29f67d05",
    "metadata": {},
    "source": [
-    "## Inspect monitor data\n",
+    "## Inspecting the monitor data\n",
     "\n",
     "The monitor counts play a major role in normalizing the detector counts.\n",
     "From the two loaded runs, we extract the data from the incident and transmission monitors."
@@ -261,18 +261,18 @@
    },
    "outputs": [],
    "source": [
-    "monitors = sans.i_of_q.preprocess_monitor_data(\n",
+    "processed_monitors = sans.i_of_q.preprocess_monitor_data(\n",
     "    monitors,\n",
     "    non_background_range=non_background_range,\n",
     "    wavelength_bins=wavelength_bins)\n",
     "\n",
     "sample_monitors = {\n",
-    "    'incident': monitors['sample-incident'],\n",
-    "    'transmission': monitors['sample-transmission']\n",
+    "    'incident': processed_monitors['sample-incident'],\n",
+    "    'transmission': processed_monitors['sample-transmission']\n",
     "}\n",
     "direct_monitors = {\n",
-    "    'incident': monitors['direct-incident'],\n",
-    "    'transmission': monitors['direct-transmission']\n",
+    "    'incident': processed_monitors['direct-incident'],\n",
+    "    'transmission': processed_monitors['direct-transmission']\n",
     "}"
    ]
   },
@@ -623,32 +623,6 @@
     "sc.plot(direct_beam)"
    ]
   },
-  {
-   "cell_type": "markdown",
-   "id": "894e8329-8987-4565-b070-5d262af746d8",
-   "metadata": {},
-   "source": [
-    "### Solid Angle\n",
-    "\n",
-    "The `sans.normalization` module also provides a utility to compute the solid angles of rectangular detector pixels:"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "id": "e126e749-3dad-4ad8-a6f9-2d446e0153e4",
-   "metadata": {
-    "tags": []
-   },
-   "outputs": [],
-   "source": [
-    "solid_angle = sans.normalization.solid_angle_of_rectangular_pixels(\n",
-    "    sample,\n",
-    "    pixel_width=sample.coords['pixel_width'],\n",
-    "    pixel_height=sample.coords['pixel_height'])\n",
-    "solid_angle"
-   ]
-  },
   {
    "cell_type": "markdown",
    "id": "92b7cba6-c4aa-458f-8564-a0c265cf4332",
@@ -657,23 +631,39 @@
     "### The denominator term\n",
     "\n",
     "We combine all the terms above to compute the `denominator`.\n",
-    "We then attach to the denominator some coordinates required to perform the conversion to $Q$."
+    "\n",
+    "**Note:**\n",
+    "\n",
+    "[Heybrock et al. (2023)](#heybrock2023) describe how broadcasting operations introduce correlations which are not tracked by Scipp's uncertainty propagation.\n",
+    "In the normalization term above, multiplying $M(\\lambda)T(\\lambda)D(\\lambda)$ (wavelength dependent) by the solid angle $\\Omega(R)$ (pixel dependent) is such a broadcasting operation.\n",
+    "The article however states that for normalization operations, under the limit where the counts of the numerator a much smaller than that of the denominator,\n",
+    "dropping the variances of the denominator is a satisfactory approximation.\n",
+    "\n",
+    "The helper function below internally calculates the solid angles of the detector pixels $\\Omega(R)$,\n",
+    "and performs a verification that our data satisfies that condition.\n",
+    "\n",
+    "We are not able to perform this check earlier in the workflow, because the verification involves comparing integrated counts,\n",
+    "and the integration needs to be performed over the same wavelength range for both the monitors and the detector counts.\n",
+    "To be able to compute the wavelengths of the detector data, we needed to first run the beam center finder."
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "30250d51-0d6f-4265-b372-881900a16e2a",
+   "id": "fc07f7b3-73dd-410d-9c1c-905d579a97cf",
    "metadata": {
     "tags": []
    },
    "outputs": [],
    "source": [
-    "denominator = sans.normalization.compute_denominator(\n",
-    "        direct_beam=direct_beam,\n",
-    "        data_incident_monitor=sample_monitors['incident'],\n",
-    "        transmission_fraction=transmission_fraction,\n",
-    "        solid_angle=solid_angle)\n",
+    "denominator = sans.normalization.iofq_denominator(\n",
+    "        data=sample,\n",
+    "        data_transmission_monitor=sample_monitors['transmission'],\n",
+    "        direct_incident_monitor=direct_monitors['incident'],\n",
+    "        direct_transmission_monitor=direct_monitors['transmission'],\n",
+    "        direct_beam=direct_beam\n",
+    "    )\n",
+    "\n",
     "# Insert a copy of coords needed for conversion to Q.\n",
     "# TODO: can this be avoided by copying the Q coords from the converted numerator?\n",
     "for coord in ['position', 'sample_position', 'source_position']:\n",
@@ -826,7 +816,7 @@
   },
   {
    "cell_type": "markdown",
-   "id": "faf45578-6b6a-4046-82aa-4864a2e8bb8b",
+   "id": "cf7ef2db-6ea0-4ea8-a0c4-04ef96831596",
    "metadata": {},
    "source": [
     "<div id='heenan1997'></div>\n",
@@ -836,6 +826,18 @@
     "[J. Appl. Cryst., 30, 1140-1147](https://doi.org/10.1107/S0021889897002173)"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "id": "faf45578-6b6a-4046-82aa-4864a2e8bb8b",
+   "metadata": {},
+   "source": [
+    "<div id='heybrock2023'></div>\n",
+    "\n",
+    "Heybrock S., Wynen J.-L., Vaytet N., **2023**,\n",
+    "*Systematic underestimation of uncertainties by widespread neutron-scattering data-reduction software*,\n",
+    "Journal of Neutron Research"
+   ]
+  },
   {
    "cell_type": "markdown",
    "id": "4ca6948e-ed72-49de-8529-4faf1474535c",

docs/techniques/sans/i_of_q.rst

@@ -18,6 +18,5 @@ Building blocks
    :toctree: ../../generated
 
    convert_to_q_and_merge_spectra
-   normalization_denominator
    preprocess_monitor_data
    resample_direct_beam

docs/techniques/sans/reference.rst

@@ -21,7 +21,6 @@ Beam center finder
    beam_center_finder.center_of_mass
    beam_center_finder.cost
    beam_center_finder.iofq_in_quadrants
-   beam_center_finder.iofq_normalization
    beam_center_finder.minimize
 
 Conversions
@@ -41,7 +40,7 @@ Normalization
 
    normalization.solid_angle_of_rectangular_pixels
    normalization.transmission_fraction
-   normalization.compute_denominator
+   normalization.iofq_denominator
    normalization.normalize
 
 Common

docs/techniques/sans/sans-beam-center-finder.ipynb

@@ -341,9 +341,12 @@
    },
    "outputs": [],
    "source": [
-    "norm = bcf.iofq_normalization(data=sample,\n",
-    "                              data_monitors=sample_monitors,\n",
-    "                              direct_monitors=direct_monitors)\n",
+    "norm = sans.normalization.iofq_denominator(\n",
+    "        data=sample,\n",
+    "        data_transmission_monitor=sc.values(sample_monitors['transmission']),\n",
+    "        direct_incident_monitor=sc.values(direct_monitors['incident']),\n",
+    "        direct_transmission_monitor=sc.values(direct_monitors['transmission']))\n",
+    "\n",
     "norm"
    ]
   },