feat: rectangular approximation for cylinder solid angle

src/esssans/beam_center_finder.py

@@ -17,14 +17,12 @@
 )
 from .i_of_q import merge_spectra
 from .logging import get_logger
-from .normalization import (
-    iofq_denominator,
-    normalize,
-    solid_angle_rectangular_approximation,
-)
+from .normalization import iofq_denominator, normalize, solid_angle
 from .types import (
     BeamCenter,
+    DetectorPixelShape,
     IofQ,
+    LabFrameTransform,
     MaskedData,
     NormWavelengthTerm,
     QBins,
@@ -121,6 +119,8 @@ def _iofq_in_quadrants(
     graph: dict,
     q_bins: Union[int, sc.Variable],
     wavelength_range: sc.Variable,
+    transform: sc.Variable,
+    pixel_shape: sc.DataGroup,
 ) -> Dict[str, sc.DataArray]:
     """
     Compute the intensity as a function of Q inside 4 quadrants in Phi.
@@ -161,13 +161,15 @@ def _iofq_in_quadrants(
         iofq_denominator,
         mask_wavelength,
         detector_to_wavelength,
+        solid_angle,
         calibrate_positions,
-        solid_angle_rectangular_approximation,
     ]
     params = {}
     params[UncertaintyBroadcastMode] = UncertaintyBroadcastMode.upper_bound
     params[WavelengthBands] = wavelength_range
     params[QBins] = q_bins
+    params[DetectorPixelShape[SampleRun]] = pixel_shape
+    params[LabFrameTransform[SampleRun]] = transform
     params[ElasticCoordTransformGraph] = graph
     params[BeamCenter] = _offsets_to_vector(data=data, xy=xy, graph=graph)
 
@@ -266,6 +268,8 @@ def beam_center_from_iofq(
     wavelength_bins: WavelengthBins,
     norm: NormWavelengthTerm[SampleRun],
     q_bins: BeamCenterFinderQBins,
+    transform: LabFrameTransform[SampleRun],
+    pixel_shape: DetectorPixelShape[SampleRun],
     minimizer: Optional[BeamCenterFinderMinimizer],
     tolerance: Optional[BeamCenterFinderTolerance],
 ) -> BeamCenter:
@@ -397,7 +401,7 @@ def beam_center_from_iofq(
     res = minimize(
         _cost,
         x0=[com_shift.fields.x.value, com_shift.fields.y.value],
-        args=(data, norm, graph, q_bins, wavelength_range),
+        args=(data, norm, graph, q_bins, wavelength_range, transform, pixel_shape),
         bounds=bounds,
         method=minimizer,
         tol=tolerance,

src/esssans/data.py

@@ -15,9 +15,11 @@ def _make_pooch():
         version=_version,
         registry={
             'DIRECT_SANS2D_REAR_34327_4m_8mm_16Feb16.hdf5': 'md5:43f4188301d709aa49df0631d03a67cb',  # noqa: E501
+            'SANS2D00063091.nxs': 'md5:05929753ea06eca5fe4be164cb06b4d6',
             'SANS2D00063091.hdf5': 'md5:1fdbe36a496e914336f2f9b5cad9f00e',
             'SANS2D00063114.hdf5': 'md5:536303077b9f55286af5ef6ec5124e1c',
             'SANS2D00063159.hdf5': 'md5:e2f2aea3ed9251e579210718de568203',
+            'SANS2D00063091.SolidAngle_from_mantid.hdf5': 'md5:d57b82db377cb1aea0beac7202713861',  # noqa: E501
         },
     )
 

src/esssans/normalization.py

@@ -4,7 +4,6 @@
 from typing import Optional
 
 import scipp as sc
-import scippneutron as scn
 from scipp.core import concepts
 
 from .types import (
@@ -14,9 +13,11 @@
     CleanMonitor,
     CleanSummedQ,
     Denominator,
+    DetectorPixelShape,
     DirectRun,
     Incident,
     IofQ,
+    LabFrameTransform,
     NormWavelengthTerm,
     Numerator,
     RunType,
@@ -32,36 +33,92 @@
 )
 
 
-def solid_angle_rectangular_approximation(
+def solid_angle(
     data: CalibratedMaskedData[RunType],
+    pixel_shape: DetectorPixelShape[RunType],
+    transform: LabFrameTransform[RunType],
 ) -> SolidAngle[RunType]:
     """
-    Solid angle computed from rectangular pixels with a 'width' and a 'height'.
+    Solid angle for cylindrical pixels.
 
-    Note that this is an approximation which is only valid for small angles
-    between the line of sight and the rectangle normal.
+    Note that the approximation is valid when the distance from sample
+    to pixel is much larger than the length or the radius of the pixels.
 
     Parameters
     ----------
     data:
-        The DataArray that contains the positions for the detector pixels and the
-        sample, as well as `pixel_width` and `pixel_height` as coordinates.
+        The DataArray that contains the positions of the detector pixels in coords.
+        The positions must be in the sample reference frame.
+
+    pixel_shape:
+        Contains the description of the detector pixel shape.
+
+    transform:
+        Transformation from the local coordinate system of the detector
+        to the coordinate system of the sample.
 
     Returns
     -------
     :
         The solid angle of the detector pixels, as viewed from the sample position.
-        Any coords and masks that have a dimension common to the dimensions of the
-        position coordinate are retained in the output.
     """
-    pixel_width = data.coords['pixel_width']
-    pixel_height = data.coords['pixel_height']
-    L2 = scn.L2(data)
-    omega = (pixel_width * pixel_height) / (L2 * L2)
-    dims = set(data.dims) - set(omega.dims)
+    face_1_center, face_1_edge, face_2_center = (
+        pixel_shape['vertices']['vertex', i] for i in range(3)
+    )
+    cylinder_axis = transform * face_2_center - transform * face_1_center
+    radius = sc.norm(face_1_center - face_1_edge)
+    length = sc.norm(cylinder_axis)
+
+    omega = _approximate_solid_angle_for_cylinder_shaped_pixel_of_detector(
+        pixel_position=data.coords['position'],
+        cylinder_axis=cylinder_axis,
+        radius=radius,
+        length=length,
+    )
     return SolidAngle[RunType](
-        concepts.rewrap_reduced_data(prototype=data, data=omega, dim=dims)
+        concepts.rewrap_reduced_data(
+            prototype=data, data=omega, dim=set(data.dims) - set(omega.dims)
+        )
+    )
+
+
+def _approximate_solid_angle_for_cylinder_shaped_pixel_of_detector(
+    pixel_position: sc.Variable,
+    cylinder_axis: sc.Variable,
+    radius: sc.Variable,
+    length: sc.Variable,
+):
+    r"""Computes solid angle of a detector pixel under the assumption that the
+    distance to the detector pixel from the sample is large compared to the radius
+    and to the length of the piece of the detector cylinder that makes up the pixel.
+
+    The cylinder is approximated by a rectangle element contained in the cylinder.
+    The normal of the rectangle is selected to
+    1. be orthogonal to the cylinder axis and
+    2. "maximally parallel" to the pixel position vector.
+    This occurs when the normal :math:`n`, the cylinder axis :math:`c` and the
+    pixel position :math:`r` fall in a shared plane such that
+    .. math::
+        r = (n\cdot r)n + (c \cdot r)c.
+
+    From the above relationship we get the scalar product of the position
+    vector and the normal:
+
+    .. math::
+        |\hat{r}\cdot n| = \sqrt{\frac{r \cdot r - (c \cdot r^2)}{||r||^2}}
+
+    The solid angle contribution from the detector element is approximated as
+    .. math::
+        \Delta\Omega =  A \frac{|\hat{r}\cdot n|}{||r||^2}
+    where :math:`A = 2 R L` is the area of the rectangular element
+    and :math:`\hat{r}` is the normalized pixel position vector.
+    """
+    norm_pp = sc.norm(pixel_position)
+    norm_ca = sc.norm(cylinder_axis)
+    cosalpha = sc.sqrt(
+        1 - (sc.dot(pixel_position, cylinder_axis) / (norm_pp * norm_ca)) ** 2
     )
+    return (2 * radius * length) * cosalpha / norm_pp**2
 
 
 def transmission_fraction(
@@ -272,5 +329,5 @@ def normalize(
     iofq_norm_wavelength_term,
     iofq_denominator,
     normalize,
-    solid_angle_rectangular_approximation,
+    solid_angle,
 )

src/esssans/sans2d.py

@@ -11,9 +11,11 @@
 from .common import gravity_vector
 from .types import (
     DetectorEdgeMask,
+    DetectorPixelShape,
     DirectBeam,
     DirectBeamFilename,
     Filename,
+    LabFrameTransform,
     LoadedFileContents,
     MaskedData,
     MonitorType,
@@ -39,8 +41,6 @@ def pooch_load(filename: Filename[RunType]) -> LoadedFileContents[RunType]:
     data = dg['data']
     if 'gravity' not in data.coords:
         data.coords["gravity"] = gravity_vector()
-    data.coords['pixel_width'] = sc.scalar(0.002033984375, unit='m')
-    data.coords['pixel_height'] = sc.scalar(0.0035, unit='m')
 
     # Some fixes specific for these Sans2d runs
     sample_pos_z_offset = 0.053 * sc.units.m
@@ -133,6 +133,38 @@ def run_title(dg: LoadedFileContents[SampleRun]) -> RunTitle:
     return RunTitle(dg['run_title'].value)
 
 
+def sans2d_tube_detector_pixel_shape() -> DetectorPixelShape[RunType]:
+    # Pixel radius and length
+    # found here:
+    # https://github.com/mantidproject/mantid/blob/main/instrument/SANS2D_Definition_Tubes.xml
+    R = 0.00405
+    L = 0.002033984375
+    pixel_shape = sc.DataGroup(
+        {
+            'vertices': sc.vectors(
+                dims=['vertex'],
+                values=[
+                    # Coordinates in pixel-local coordinate system
+                    # Bottom face center
+                    [0, 0, 0],
+                    # Bottom face edge
+                    [R, 0, 0],
+                    # Top face center
+                    [0, L, 0],
+                ],
+                unit='m',
+            ),
+            'nexus_class': 'NXcylindrical_geometry',
+        }
+    )
+    return pixel_shape
+
+
+def lab_frame_transform() -> LabFrameTransform[RunType]:
+    # Rotate +y to -x
+    return sc.spatial.rotation(value=[0, 0, 1 / 2**0.5, 1 / 2**0.5])
+
+
 providers = (
     pooch_load_direct_beam,
     pooch_load,
@@ -143,6 +175,8 @@ def run_title(dg: LoadedFileContents[SampleRun]) -> RunTitle:
     mask_detectors,
     run_number,
     run_title,
+    lab_frame_transform,
+    sans2d_tube_detector_pixel_shape,
 )
 """
 Providers for loading and masking Sans2d data.

src/esssans/types.py

@@ -114,6 +114,15 @@ class TransmissionFraction(sciline.Scope[RunType, sc.DataArray], sc.DataArray):
 """Direct beam after resampling to required wavelength bins"""
 
 
+class DetectorPixelShape(sciline.Scope[RunType, sc.DataGroup], sc.DataGroup):
+    """Geometry of the detector from description in nexus file."""
+
+
+class LabFrameTransform(sciline.Scope[RunType, sc.Variable], sc.Variable):
+    """Coordinate transformation from detector local coordinates
+    to the sample frame of reference."""
+
+
 class SolidAngle(sciline.Scope[RunType, sc.DataArray], sc.DataArray):
     """Solid angle of detector pixels seen from sample position"""