Include wavelength component in Q resolution calculation

.gitattributes

@@ -0,0 +1,2 @@
+# SCM syntax highlighting & preventing 3-way merges
+pixi.lock merge=binary linguist-language=YAML linguist-generated=true

docs/api/lr_reduction.rst

@@ -87,10 +87,10 @@ lr_reduction.time_resolved
    :undoc-members:
    :show-inheritance:
 
-lr_reduction.typing
+lr_reduction.types
 -------------------
 
-.. automodule:: lr_reduction.typing
+.. automodule:: lr_reduction.types
    :members:
    :undoc-members:
    :show-inheritance:

pyproject.toml

@@ -87,6 +87,8 @@ matplotlib = "==3.9.4"
 plot-publisher = ">=1.0,<2.0"
 pyqt = ">=5,<6"
 qtpy = ">=2.4.3"
+sympy = ">=1.14.0,<2"
+jupyter = ">=1.1.1,<2"
 
 # Pypi packages to be installed in the environments with `pip`
 [tool.pixi.pypi-dependencies]
@@ -102,6 +104,7 @@ matplotlib = "==3.9.4"
 plot-publisher = ">=1.0,<2.0"
 pyqt = ">=5,<6"
 qtpy = ">=2.4.3"
+sympy = ">=1.14.0,<2"
 
 [tool.pixi.package]
 name = "lr_reduction"

src/lr_reduction/data_info.py

@@ -37,23 +37,66 @@ def from_workspace(cls, input_workspace: MantidWorkspace):
         sample_logs = SampleLogValues(input_workspace)
         value = cls.REFLECTED_BEAM
         try:
+            coordinate_system = CoordinateSystem.from_workspace(input_workspace)
             # Determine whether this is a direct beam based on the geometry
-            if (("BL4B:CS:Mode:Coordinates" in sample_logs and sample_logs["BL4B:CS:Mode:Coordinates"] == 0) or  # This is a new log for earth-centered
-                    sample_logs["BL4B:CS:ExpPl:OperatingMode"] == "Free Liquid"):  # This is for backward compatibility from before the new log value
-                # Earth-centered coordinate system
+            if coordinate_system == CoordinateSystem.EARTH_CENTERED:
                 thi = sample_logs["thi"]
                 tthd = sample_logs["tthd"]
                 if np.isclose(thi, tthd, atol=0.01):
                     value = cls.DIRECT_BEAM
-            else:
+            elif coordinate_system == CoordinateSystem.BEAM_CENTERED:
                 # Beam-centered coordinate system
                 ths = sample_logs["ths"]
                 tthd = sample_logs["tthd"]
                 if np.fabs(tthd) < 0.001 and np.fabs(ths) < 0.001:
                     value = cls.DIRECT_BEAM
-        except KeyError as e:
+            else:
+                logger.warning("Unknown coordinate system; assuming reflected beam")
+        except RuntimeError as e:
             logger.warning(f"Missing sample log {e}, assuming reflected beam")
         return value
 
     def __str__(self):
         return self.name
+
+
+class CoordinateSystem(IntEnum):
+    """
+    Enum to represent the coordinate system used in the experiment.
+
+    Attributes:
+        UNKNOWN (int): Represents unknown coordinate system.
+        EARTH_CENTERED (int): Represents the Earth-centered coordinate system.
+        BEAM_CENTERED (int): Represents the Beam-centered coordinate system.
+    """
+
+    UNKNOWN = -1
+    EARTH_CENTERED = 0
+    BEAM_CENTERED = 1
+
+    @classmethod
+    def from_workspace(cls, input_workspace: MantidWorkspace):
+        """
+        Determine the coordinate system from the given workspace.
+        """
+        sample_logs = SampleLogValues(input_workspace)
+        value = cls.UNKNOWN
+        try:
+            # This is the new log for coordinate system mode
+            if "BL4B:CS:Mode:Coordinates" in sample_logs:
+                if sample_logs["BL4B:CS:Mode:Coordinates"] == 0:
+                    value = cls.EARTH_CENTERED
+                else:
+                    value = cls.BEAM_CENTERED
+            # Fallback to older method using operating mode
+            elif "BL4B:CS:ExpPl:OperatingMode" in sample_logs:
+                if sample_logs["BL4B:CS:ExpPl:OperatingMode"] == "Free Liquid":
+                    value = cls.EARTH_CENTERED
+                else:
+                    value = cls.BEAM_CENTERED
+        except RuntimeError as e:
+            logger.warning(f"Missing sample log {e}, unable to determine coordinate system")
+        return value
+
+    def __str__(self):
+        return self.name

src/lr_reduction/event_reduction.py

@@ -12,8 +12,11 @@
 import numpy as np
 from scipy.optimize import brentq
 
+from lr_reduction.data_info import CoordinateSystem
 from lr_reduction.gravity_correction import GravityDirection, gravity_correction
 from lr_reduction.instrument_settings import InstrumentSettings
+from lr_reduction.types import MantidWorkspace
+from lr_reduction.user_defined_function import UserDefinedFunction
 from lr_reduction.utils import mantid_algorithm_exec
 
 from . import background, dead_time_correction
@@ -573,7 +576,6 @@ def to_dict(self):
             norm_run=norm_run,
             time=time.ctime(),
             dq0=dq0,
-            dq_over_q=self.dq_over_q,
             sequence_number=sequence_number,
             sequence_id=sequence_id,
             q_summing=self.q_summing,
@@ -608,6 +610,8 @@ def specular(self, q_summing=False, tof_weighted=False, bck_in_q=False, clean=Fa
             The reflectivity values
         d_refl
             The uncertainties in the reflectivity values
+        dq_over_q_bins
+            The Q resolution associated with the q_bins
         """
         if tof_weighted:
             self.specular_weighted(q_summing=q_summing, bck_in_q=bck_in_q)
@@ -631,12 +635,13 @@ def specular(self, q_summing=False, tof_weighted=False, bck_in_q=False, clean=Fa
             self.d_refl = self.d_refl[idx[:-1]]
             self.q_bins = self.q_bins[idx]
 
-        # Compute Q resolution
-        self.dq_over_q = compute_resolution(self._ws_sc, theta=self.theta, q_summing=q_summing)
-        # TODO: integrate wavelength component to self.dq_over_q
         self.q_summing = q_summing
 
-        return self.q_bins, self.refl, self.d_refl
+        # Compute Q resolution
+        # For now keep the dq_over_q the same length as the q_bins to ensure conversion to mid-points is the same.
+        lambda_bin_list = 4 * np.pi * np.sin(self.theta) / self.q_bins
+        dq_over_q_bins = compute_resolution(self._ws_sc, wl_list=lambda_bin_list, theta=self.theta, q_summing=q_summing)
+        return self.q_bins, self.refl, self.d_refl, dq_over_q_bins
 
     def specular_unweighted(self, q_summing=False, normalize=True):
         """
@@ -1264,71 +1269,98 @@ def emission_time_correction(self, ws, tofs):
         return tofs
 
 
-def compute_resolution(ws, default_dq=0.027, theta=None, q_summing=False):
+def compute_theta_resolution(ws: MantidWorkspace, default_dq: float=0.027, theta: float=None, q_summing: bool=False) -> tuple[float, float]:
     """
-    Compute the Q resolution from the meta data.
+    Compute the theta component of the q resolution.
     Corrected to include both slits.
 
     Parameters
     ----------
     ws : mantid.api.Workspace
         Mantid workspace to extract correction meta-data from.
+    default_dq : float
+        Default dth/th value to use if resolution cannot be computed from logs.
     theta : float
-        Scattering angle in radians
+        Scattering angle in radians. If not provided will read ths or thi from ws.
     q_summing : bool
         If True, the pixel size will be used for the resolution
 
     Returns
     -------
-    float
-        The dQ/Q resolution (FWHM)
+    tuple[float, float]
+        The theta and dtheta values in degrees
     """
     settings = read_settings(ws)
 
     if theta is None:
-        theta = abs(ws.getRun().getProperty("ths").value[0]) * np.pi / 180.0
+        coordinate_system = CoordinateSystem.from_workspace(ws)
+        if coordinate_system == CoordinateSystem.EARTH_CENTERED:
+            Theta_deg = abs(ws.getRun().getProperty("thi").value[0])
+            theta = np.radians(Theta_deg)
+        else:
+            Theta_deg = abs(ws.getRun().getProperty("ths").value[0])
+            theta = np.radians(Theta_deg)
 
     if q_summing:
-        # Q resolution is reported as FWHM, so here we consider this to be
-        # related to the pixel width
+        # Compute pixel contribution to angular resolution
         sdd = 1830
         if settings.sample_detector_distance:
             sdd = settings.sample_detector_distance * 1000
-        pixel_size = 0.7
+        pixel_size = 0.7 # mm, default pixel size for BL4B
         if settings.pixel_width:
             pixel_size = settings.pixel_width
 
-        # All angles here in radians, assuming small angles
-        dq_over_q = np.arcsin(pixel_size / sdd) / theta
-        print("Q summing: %g" % dq_over_q)
-        return dq_over_q
+        det_res = 1.0  # mm, approximate value for detector resolution contribution
+        # Standard deviation of a top-hat function (uniform distribution)
+        sigma_pix = pixel_size / np.sqrt(12)
+        sigma_y = np.sqrt((det_res ** 2 + sigma_pix ** 2))
+
+        dtheta = np.arctan(sigma_y / sdd) # in radians
+
+        # Want to export in degrees
+        dtheta_deg = np.degrees(dtheta)
+        theta_deg = np.degrees(theta)
+
+        #print("Q summing: %g" % dq_over_q) ## Work out some print functions throughout.
+        return theta_deg, dtheta_deg
 
     # We can't compute the resolution if the value of xi is not in the logs.
     # Since it was not always logged, check for it here.
     if not ws.getRun().hasProperty("BL4B:Mot:xi.RBV"):
         # For old data, the resolution can't be computed, so use
         # the standard value for the instrument
         print("Could not find BL4B:Mot:xi.RBV: using supplied dQ/Q")
-        return default_dq
+        theta_deg = np.degrees(theta)
+        return theta_deg, default_dq
 
-    # Xi reference would be the position of xi if the si slit were to be positioned
-    # at the sample. The distance from the sample to si is then xi_reference - xi.
-    xi_reference = 445
-    if ws.getInstrument().hasParameter("xi-reference"):
-        xi_reference = ws.getInstrument().getNumberParameter("xi-reference")[0]
+    # Compute the trapezoidal equivalent sigma of the angular distribution
+    theta_deg = np.degrees(theta)
+    _, _, _, dth_over_th = trapezoidal_distribution_params(
+        ws, Theta_deg=theta_deg
+    )
+    dtheta_deg = dth_over_th * theta_deg
 
-    # Distance between the s1 and the sample
-    s1_sample_distance = 1485
-    if settings.s1_sample_distance is not None:
-        s1_sample_distance = settings.s1_sample_distance * 1000
-
-    # Get slit gap openings.
-    s1h = abs(ws.getRun().getProperty("S1VHeight").value[0])
-    sih = abs(ws.getRun().getProperty("SiVHeight").value[0])
-    xi = abs(ws.getRun().getProperty("BL4B:Mot:xi.RBV").value[0])
-    sample_si_distance = xi_reference - xi
-    slit_distance = s1_sample_distance - sample_si_distance
-    dq_over_q = (s1h + sih) * 0.5 / slit_distance / theta
+    return theta_deg, dtheta_deg
+
+
+def compute_resolution(ws: MantidWorkspace, wl_list: np.ndarray, theta: float=None, q_summing: bool=False) -> np.ndarray:
+    """
+    Compute the q resolution including both theta and lambda terms.
+
+    :param ws: workspace for meta-data. If this is a lambda workspace
+        already it can be used as the lambda imports in place of the wl_list.
+    :param wl_list: Provide a list of wavelengths for the lambda calculation.
+    :param theta: option to overwrite theta input, otherwise uses ths/thi value.
+    :param q_summing: Changes the angular calculation to be based on slits or pixel sizes depending on method of q conversion.
+    """
+    ## Need to check through all the None's etc.
+    ## Need to check separate outputs...
+    theta_deg, delta_th_deg = compute_theta_resolution(ws, theta=theta, q_summing=q_summing)
+
+    wl_resolution_function = read_settings(ws).wavelength_resolution_function
+    lambda_list, delta_lam = compute_wavelength_resolution(wl_list, wl_resolution_function)
+
+    dq_over_q = np.sqrt((delta_lam / lambda_list) ** 2 + (delta_th_deg / theta_deg) ** 2)
     return dq_over_q
 
 
@@ -1347,7 +1379,7 @@ def trapezoidal_distribution_params(ws, Theta_deg=None, FootPrint=None, SlitRati
     ws : mantid.api.Workspace
         Mantid workspace to extract correction meta-data from.
     Theta_deg : float
-        Incident beam angle in degrees.
+        Incident beam angle in degrees. If not provided will read ths or thi from ws.
     FootPrint : float (optional)
         Beam footprint length on sample [mm].
     SlitRatio : float (optional)
@@ -1385,7 +1417,11 @@ def trapezoidal_distribution_params(ws, Theta_deg=None, FootPrint=None, SlitRati
     dS1Si = dS1Samp - dSiSamp
 
     if Theta_deg is None:
-        Theta_deg = abs(ws.getRun().getProperty("ths").value[0])
+        coordinate_system = CoordinateSystem.from_workspace(ws)
+        if coordinate_system == CoordinateSystem.EARTH_CENTERED:
+            Theta_deg = abs(ws.getRun().getProperty("thi").value[0])
+        else:
+            Theta_deg = abs(ws.getRun().getProperty("ths").value[0])
 
     # Slit openings - can be taken from a FootPrint and SlitRatio if provided, or the logs
     if FootPrint is not None and SlitRatio is not None:
@@ -1449,37 +1485,30 @@ def func(x):
 ## Fix q-bin error and auto trimming points
 
 
-def compute_wavelength_resolution(ws):
+def compute_wavelength_resolution(wl_list: np.ndarray, resolution_function_str: str) -> tuple[np.ndarray, np.ndarray]:
     """
-    Compute the wavelength resolution from the meta data.
+    Compute the wavelength resolution using the given resolution function.
 
     Parameters
     ----------
-    ws : mantid.api.Workspace
-        Mantid workspace to extract correction meta-data from.
+    wl_list : np.ndarray
+        Wavelength values to compute the resolution for.
+    resolution_function_str : str
+        User-defined resolution function.
 
     Returns
     -------
     tuple of np.ndarray
         (wavelength, d_lambda):
-            wavelength: the fitted wavelength values
-            d_lambda: the difference between wavelength and the fit
-
-    Raises
-    ------
-    ValueError : if ws does not have exactly one spectrum
+            wavelength: the wavelength values
+            d_lambda: the wavelength resolution values
     """
+    # Parse user-defined resolution function
+    resolution_function = UserDefinedFunction(resolution_function_str)
+    d_lambda = np.array(resolution_function(wl_list))
+    wavelength = np.array(wl_list)
 
-    if ws.spectrumInfo().size() != 1:
-        raise ValueError("Workspace must have only one spectrum")
-
-    settings = read_settings(ws)
-
-    out = api.EvaluateFunction(
-        Function=settings.wavelength_resolution_function, InputWorkspace=ws, OutputWorkspace="out"
-    )
-
-    wavelength = out.readX(1)
-    d_lambda = out.readY(1)
+    # Set any negative values to zero
+    d_lambda[d_lambda < 0] = 0
 
-    return np.array(wavelength), np.array(d_lambda)
+    return wavelength, d_lambda

src/lr_reduction/gravity_correction.py

@@ -2,6 +2,7 @@
 
 import numpy as np
 
+from lr_reduction.data_info import CoordinateSystem
 from lr_reduction.types import MantidWorkspace
 from lr_reduction.utils import workspace_handle
 
@@ -94,15 +95,11 @@ def _theta_in(workspace: MantidWorkspace) -> float:
 
     # Angle calculated from thi and a flag on earth-centered vs beam-centered
     thi = _log_value(run, "BL4B:Mot:thi.RBV")
-    if "BL4B:CS:Mode:Coordinates" in run:
-        if _log_value(run, "BL4B:CS:Mode:Coordinates") == 0:  # Earth-centered=0
-            theta_in = thi
-        else:
-            theta_in = thi - 4.0  # Beamline optics gives -4 deg incline. In future will have PV.
-    elif _log_value(run, "BL4B:CS:ExpPl:OperatingMode", default="") == "Free Liquid":
+    coordinate_system = CoordinateSystem.from_workspace(workspace)
+    if coordinate_system == CoordinateSystem.EARTH_CENTERED:
         theta_in = thi
     else:
-        theta_in = thi - 4.0
+        theta_in = thi - 4.0  # Beamline optics gives -4 deg incline. In future will have PV.
     return abs(theta_in)