Added basic line profile to tof integration.

src/dials/algorithms/integration/fit/tof_line_profile.py

@@ -0,0 +1,135 @@
+from __future__ import annotations
+
+import numpy as np
+from scipy import integrate
+from scipy.optimize import curve_fit
+from scipy.special import erfc
+
+import cctbx.array_family.flex
+from dxtbx import flumpy
+
+from dials.algorithms.shoebox import MaskCode
+
+
+class BackToBackExponential:
+
+    """
+    https://www.nature.com/articles/srep36628.pdf
+    """
+
+    def __init__(self, tof, intensities, A, alpha, beta, sigma, T):
+        self.intensities = intensities
+        self.tof = tof
+        self.params = (A, alpha, beta, sigma, T)
+        self.cov = None
+
+    def func(self, tof, A, alpha, beta, sigma, T):
+        dT = tof - T
+        sigma2 = np.square(sigma)
+        sigma_sqrt = np.sqrt(2 * sigma2)
+
+        u = alpha * 0.5 * (alpha * sigma2 + 2 * dT)
+        v = beta * 0.5 * (beta * sigma2 - 2 * dT)
+        y = (alpha * sigma2 + dT) / sigma_sqrt
+        z = (beta * sigma2 - dT) / sigma_sqrt
+
+        N = (alpha * beta) / (2 * (alpha + beta))
+        exp_u = np.exp(u)
+        exp_v = np.exp(v)
+        erfc_y = erfc(y)
+        erfc_z = erfc(z)
+
+        result = A
+        result *= N
+        result *= exp_u * erfc_y + exp_v * erfc_z
+        return result
+
+    def fit(self):
+        params, cov = curve_fit(
+            f=self.func,
+            xdata=self.tof,
+            ydata=self.intensities,
+            p0=self.params,
+            bounds=((1, 0, 0, 1, min(self.tof)), (1000, 1, 10, 1000, max(self.tof))),
+            max_nfev=10000000,
+        )
+        self.params = params
+        self.cov = cov
+
+    def result(self):
+        return self.func(self.tof, *(self.params))
+
+    def calc_intensity(self):
+        predicted = self.result()
+        return integrate.simpson(predicted, self.tof)
+
+
+def compute_line_profile_intensity(reflections):
+
+    A = 200.0
+    alpha = 0.4
+    beta = 0.4
+    sigma = 8.0
+
+    bg_code = MaskCode.Valid | MaskCode.Background | MaskCode.BackgroundUsed
+
+    fit_intensities = cctbx.array_family.flex.double(len(reflections))
+    fit_variances = cctbx.array_family.flex.double(len(reflections))
+
+    for i in range(len(reflections)):
+        shoebox = reflections[i]["shoebox"]
+        data = flumpy.to_numpy(shoebox.data).ravel()
+        background = flumpy.to_numpy(shoebox.background).ravel()
+        mask = flumpy.to_numpy(shoebox.mask).ravel()
+        coords = flumpy.to_numpy(shoebox.coords())
+        m = mask & MaskCode.Foreground == MaskCode.Foreground
+        bg_m = mask & bg_code == bg_code
+        n_background = np.sum(np.bitwise_and(~m, bg_m))
+
+        m = np.bitwise_and(m, mask & MaskCode.Valid == MaskCode.Valid)
+        m = np.bitwise_and(m, mask & MaskCode.Overlapped == 0)
+
+        n_signal = np.sum(m)
+
+        background = background[m]
+        intensity = data[m] - background
+        background_sum = np.sum(background)
+        coords = coords[m]
+        tof = coords[:, 2]
+
+        summed_values = {}
+
+        for j in np.unique(tof):
+            indices = np.where(tof == j)
+            summed_values[j] = np.sum(intensity[indices])
+
+        # Remove background and project onto ToF axis
+        projected_intensity = np.array(list(summed_values.values()))
+        tof = np.array(list(summed_values.keys()))
+
+        try:
+            T = tof[np.argmax(projected_intensity)]
+            l = BackToBackExponential(
+                tof=tof,
+                intensities=projected_intensity,
+                A=A,
+                alpha=alpha,
+                beta=beta,
+                sigma=sigma,
+                T=T,
+            )
+            l.fit()
+            fit_intensity = l.calc_intensity()
+            fit_intensities[i] = fit_intensity
+        except ValueError:
+            fit_intensities[i] = -1
+
+        if n_background > 0:
+            m_n = n_signal / n_background
+        else:
+            m_n = 0.0
+        fit_variances[i] = abs(fit_intensity) + abs(background_sum) * (1.0 + m_n)
+
+    reflections["line_profile_intensity"] = fit_intensities
+    reflections["line_profile_variance"] = fit_variances
+    return reflections

src/dials/command_line/simple_tof_integrate.py

@@ -17,6 +17,9 @@
 
 import dials.util.log
 import dials_array_family_flex_ext
+from dials.algorithms.integration.fit.tof_line_profile import (
+    compute_line_profile_intensity,
+)
 from dials.algorithms.integration.report import IntegrationReport, ProfileModelReport
 from dials.algorithms.profile_model.gaussian_rs import GaussianRSProfileModeller
 from dials.algorithms.profile_model.gaussian_rs import Model as GaussianRSProfileModel
@@ -298,6 +301,11 @@ def get_corrected_intensity_and_sigma(reflections, idx):
         variance = reflections["intensity.sum.variance"][idx]
         return intensity, np.sqrt(variance)
 
+    def get_line_profile_intensity_and_sigma(reflections, idx):
+        intensity = reflections["line_profile_intensity"][idx]
+        variance = reflections["line_profile_variance"][idx]
+        return intensity, np.sqrt(variance)
+
     def valid_intensity(intensity):
         from math import isinf, isnan
 
@@ -334,6 +342,35 @@ def valid_intensity(intensity):
             )
         )
 
+    with open("lp_" + filename, "w") as g:
+        for i in range(len(reflections)):
+            h, k, l = reflections["miller_index"][i]
+            batch_number = 1
+            intensity, sigma = get_line_profile_intensity_and_sigma(reflections, i)
+            if not valid_intensity(intensity):
+                continue
+            intensity = round(intensity, 2)
+            sigma = round(sigma, 2)
+            wavelength = round(reflections["wavelength_cal"][i], 4)
+            g.write(
+                ""
+                + "{:4d}{:4d}{:4d}{:8.1f}{:8.2f}{:4d}{:8.4f}\n".format(
+                    int(h),
+                    int(k),
+                    int(l),
+                    float(intensity),
+                    float(sigma),
+                    int(batch_number),
+                    float(wavelength),
+                )
+            )
+        g.write(
+            ""
+            + "{:4d}{:4d}{:4d}{:8.1f}{:9.2f}{:4d}{:8.4f}\n".format(
+                int(0), int(0), int(0), float(0.00), float(0.00), int(0), float(0.0000)
+            )
+        )
+
 
 def output_expt_as_ins(expt, filename):
     def LATT_SYMM(s, space_group, decimal=False):
@@ -605,6 +642,7 @@ def run_simple_integrate(params, experiments, reflections):
     centroid_algorithm.compute_centroid(predicted_reflections)
 
     predicted_reflections.compute_summed_intensity()
+    predicted_reflections = compute_line_profile_intensity(predicted_reflections)
 
     # Filter reflections with a high fraction of masked foreground
     valid_foreground_threshold = 1.0  # DIALS default