diff --git a/src/dials/algorithms/integration/fit/tof_line_profile.py b/src/dials/algorithms/integration/fit/tof_line_profile.py new file mode 100644 index 0000000000..b8080d4777 --- /dev/null +++ b/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 diff --git a/src/dials/command_line/simple_tof_integrate.py b/src/dials/command_line/simple_tof_integrate.py index 8af38c39fa..676d80191e 100644 --- a/src/dials/command_line/simple_tof_integrate.py +++ b/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