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