From 485c6955e0631c1dbb2c085ac83c5ae0572a19fc Mon Sep 17 00:00:00 2001 From: davidmcdonagh Date: Thu, 18 Jul 2024 11:28:18 +0100 Subject: [PATCH] Added basic tof_integrate.py script. Added scaling methods. --- .../integration_integrator_ext.cc | 24 +- .../integration/fit/tof_line_profile.py | 230 +++++++ .../scaling/tof_scaling_corrections.h | 583 ++++++++++++++++++ src/dials/command_line/tof_integrate.py | 488 +++++++++++++++ 4 files changed, 1316 insertions(+), 9 deletions(-) create mode 100644 src/dials/algorithms/integration/fit/tof_line_profile.py create mode 100644 src/dials/algorithms/scaling/tof_scaling_corrections.h create mode 100644 src/dials/command_line/tof_integrate.py diff --git a/src/dials/algorithms/integration/boost_python/integration_integrator_ext.cc b/src/dials/algorithms/integration/boost_python/integration_integrator_ext.cc index 9e27aae925..f62d06c698 100644 --- a/src/dials/algorithms/integration/boost_python/integration_integrator_ext.cc +++ b/src/dials/algorithms/integration/boost_python/integration_integrator_ext.cc @@ -351,21 +351,27 @@ namespace dials { namespace algorithms { namespace boost_python { &max_memory_needed, (arg("reflection table"), arg("frame0"), arg("frame1"), arg("flatten"))); + class_("TOFCorrectionsData", no_init) + .def(init()); + def("tof_extract_shoeboxes_to_reflection_table", &tof_extract_shoeboxes_to_reflection_table, (arg("reflection_table"), arg("experiment"), + arg("data"), arg("incident_data"), arg("empty_data"), - arg("proton_charges"), - arg("sample_radius"), - arg("sample_scattering_x_section"), - arg("sample_absorption_x_section"), - arg("sample_number_density"), - arg("incident_radius"), - arg("incident_scattering_x_section"), - arg("incident_absorption_x_section"), - arg("incident_number_density"), + arg("corrections_data"), arg("apply_lorentz_correction"), arg("apply_spherical_absorption_correction"))); 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..21daff232b --- /dev/null +++ b/src/dials/algorithms/integration/fit/tof_line_profile.py @@ -0,0 +1,230 @@ +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)), + (100000000000, 1, 100000000, 10000000000, 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_data_for_reflection( + reflection_table, A=200.0, alpha=0.4, beta=0.4, sigma=8.0 +): + + assert len(reflection_table) == 1 + + bg_code = MaskCode.Valid | MaskCode.Background | MaskCode.BackgroundUsed + + shoebox = reflection_table["shoebox"][0] + 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) + summation_intensity = float(np.sum(intensity)) + coords = coords[m] + tof = coords[:, 2] + + summed_values = {} + summed_background_values = {} + + for j in np.unique(tof): + indices = np.where(tof == j) + summed_values[j] = np.sum(intensity[indices]) + summed_background_values[j] = np.sum(background[indices]) + + # Remove background and project onto ToF axis + projected_intensity = np.array(list(summed_values.values())) + projected_background = np.array(list(summed_background_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() + line_profile = l.result() + fit_intensity = integrate.simpson(line_profile, tof) + except ValueError: + return [], [], [], [], -1, -1, -1, -1 + + if n_background > 0: + m_n = n_signal / n_background + else: + m_n = 0.0 + fit_std = np.sqrt(abs(fit_intensity) + abs(background_sum) * (1.0 + m_n)) + summation_std = np.sqrt( + abs(summation_intensity) + abs(background_sum) * (1.0 + m_n) + ) + + return ( + tof, + projected_intensity, + projected_background, + line_profile, + fit_intensity, + fit_std, + summation_intensity, + summation_std, + ) + + +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())) + + fit_intensity = None + 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 + fit_variances[i] = -1 + continue + + if n_background > 0: + m_n = n_signal / n_background + else: + m_n = 0.0 + fit_variance = abs(fit_intensity) + abs(background_sum) * (1.0 + m_n) + fit_variances[i] = fit_variance + + reflections["intensity.prf.value"] = fit_intensities + reflections["intensity.prf.variance"] = fit_variances + reflections.set_flags( + reflections["intensity.prf.value"] < 0, + reflections.flags.failed_during_profile_fitting, + ) + reflections.set_flags( + reflections["intensity.prf.value"] > 0, + reflections.flags.integrated_prf, + ) + return reflections diff --git a/src/dials/algorithms/scaling/tof_scaling_corrections.h b/src/dials/algorithms/scaling/tof_scaling_corrections.h new file mode 100644 index 0000000000..ee4081fa45 --- /dev/null +++ b/src/dials/algorithms/scaling/tof_scaling_corrections.h @@ -0,0 +1,583 @@ + +#ifndef DIALS_ALGORITHMS_SCALING_TOF_SCALING_CORRECTIONS_H +#define DIALS_ALGORITHMS_SCALING_TOF_SCALING_CORRECTIONS_H + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +namespace dials { namespace algorithms { + + using dials::model::Background; + using dials::model::Foreground; + using dxtbx::ImageSequence; + using dxtbx::af::flex_table; + using dxtbx::model::Detector; + using dxtbx::model::Experiment; + using dxtbx::model::Goniometer; + using dxtbx::model::PolychromaticBeam; + using dxtbx::model::Scan; + using dxtbx::model::scan_property_types; + using scitbx::deg_as_rad; + using scitbx::mat3; + using scitbx::vec2; + using scitbx::vec3; + using scitbx::constants::m_n; + using scitbx::constants::pi; + using scitbx::constants::Planck; + + double tof_pixel_spherical_absorption_correction(double pixel_data, + double muR, + double two_theta, + int two_theta_idx) { + // Taken from + // https://github.com/mantidproject/mantid/blob/main/Framework/Crystal/inc/MantidCrystal/AnvredCorrection.h + const double pc[8][19] = {{-6.4910e-07, + -6.8938e-07, + -7.8149e-07, + 8.1682e-08, + 1.8008e-06, + 3.3916e-06, + 4.5095e-06, + 4.7970e-06, + 4.4934e-06, + 3.6700e-06, + 2.5881e-06, + 1.5007e-06, + 3.7669e-07, + -7.9487e-07, + -1.7935e-06, + -2.5563e-06, + -3.1113e-06, + -3.3993e-06, + -3.5091e-06}, + {1.0839e-05, + 1.1582e-05, + 1.1004e-05, + -2.2848e-05, + -8.1974e-05, + -1.3268e-04, + -1.6486e-04, + -1.6839e-04, + -1.5242e-04, + -1.1949e-04, + -7.8682e-05, + -3.7973e-05, + 2.9117e-06, + 4.4823e-05, + 8.0464e-05, + 1.0769e-04, + 1.2753e-04, + 1.3800e-04, + 1.4190e-04}, + {8.7140e-05, + 9.0870e-05, + 1.6706e-04, + 6.9008e-04, + 1.4781e-03, + 2.0818e-03, + 2.3973e-03, + 2.3209e-03, + 1.9935e-03, + 1.4508e-03, + 8.1903e-04, + 1.9608e-04, + -4.1128e-04, + -1.0205e-03, + -1.5374e-03, + -1.9329e-03, + -2.2212e-03, + -2.3760e-03, + -2.4324e-03}, + {-2.9549e-03, + -3.1360e-03, + -4.2431e-03, + -8.1103e-03, + -1.2989e-02, + -1.6012e-02, + -1.6815e-02, + -1.4962e-02, + -1.1563e-02, + -6.8581e-03, + -1.7302e-03, + 3.2400e-03, + 7.9409e-03, + 1.2528e-02, + 1.6414e-02, + 1.9394e-02, + 2.1568e-02, + 2.2758e-02, + 2.3182e-02}, + {1.7934e-02, + 1.9304e-02, + 2.4706e-02, + 3.6759e-02, + 4.8351e-02, + 5.1049e-02, + 4.5368e-02, + 3.0864e-02, + 1.2086e-02, + -1.0254e-02, + -3.2992e-02, + -5.4495e-02, + -7.4205e-02, + -9.2818e-02, + -1.0855e-01, + -1.2068e-01, + -1.2954e-01, + -1.3451e-01, + -1.3623e-01}, + {6.2799e-02, + 6.3892e-02, + 6.4943e-02, + 6.4881e-02, + 7.2169e-02, + 9.5669e-02, + 1.3082e-01, + 1.7694e-01, + 2.2559e-01, + 2.7655e-01, + 3.2483e-01, + 3.6888e-01, + 4.0783e-01, + 4.4330e-01, + 4.7317e-01, + 4.9631e-01, + 5.1334e-01, + 5.2318e-01, + 5.2651e-01}, + {-1.4949e+00, + -1.4952e+00, + -1.4925e+00, + -1.4889e+00, + -1.4867e+00, + -1.4897e+00, + -1.4948e+00, + -1.5025e+00, + -1.5084e+00, + -1.5142e+00, + -1.5176e+00, + -1.5191e+00, + -1.5187e+00, + -1.5180e+00, + -1.5169e+00, + -1.5153e+00, + -1.5138e+00, + -1.5125e+00, + -1.5120e+00}, + {0.0000e+00, + 0.0000e+00, + 0.0000e+00, + 0.0000e+00, + 0.0000e+00, + 0.0000e+00, + 0.0000e+00, + 0.0000e+00, + 0.0000e+00, + 0.0000e+00, + 0.0000e+00, + 0.0000e+00, + 0.0000e+00, + 0.0000e+00, + 0.0000e+00, + 0.0000e+00, + 0.0000e+00, + 0.0000e+00, + 0.0000e+00}}; + double ln_t1 = 0; + double ln_t2 = 0; + for (std::size_t k = 0; k < 8; ++k) { + ln_t1 = ln_t1 * muR + pc[k][two_theta_idx]; + ln_t2 = ln_t2 * muR + pc[k][two_theta_idx + 1]; + } + const double t1 = exp(ln_t1); + const double t2 = exp(ln_t2); + const double sin_theta_1 = pow(sin(deg_as_rad(two_theta_idx * 5.0)), 2); + const double sin_theta_2 = pow(sin(deg_as_rad((two_theta_idx + 1) * 5.0)), 2); + const double l1 = (t1 - t2) / (sin_theta_1 - sin_theta_2); + const double l0 = t1 - l1 * sin_theta_1; + const double correction = 1 / (l0 + l1 * pow(sin(two_theta * .5), 2)); + return correction; + } + + struct TOFCorrectionsData { + double sample_proton_charge; + double incident_proton_charge; + double empty_proton_charge; + double sample_radius; + double sample_scattering_x_section; + double sample_absorption_x_section; + double sample_number_density; + double incident_radius; + double incident_scattering_x_section; + double incident_absorption_x_section; + double incident_number_density; + + TOFCorrectionsData(double sample_proton_charge, + double incident_proton_charge, + double empty_proton_charge, + double sample_radius, + double sample_scattering_x_section, + double sample_absorption_x_section, + double sample_number_density, + double incident_radius, + double incident_scattering_x_section, + double incident_absorption_x_section, + double incident_number_density) + : sample_proton_charge(sample_proton_charge), + incident_proton_charge(incident_proton_charge), + empty_proton_charge(empty_proton_charge), + sample_radius(sample_radius), + sample_scattering_x_section(sample_scattering_x_section), + sample_absorption_x_section(sample_absorption_x_section), + sample_number_density(sample_number_density), + incident_radius(incident_radius), + incident_scattering_x_section(incident_scattering_x_section), + incident_absorption_x_section(incident_absorption_x_section), + incident_number_density(incident_number_density) {} + }; + + void tof_extract_shoeboxes_to_reflection_table( + af::reflection_table &reflection_table, + Experiment &experiment, + ImageSequence &data, + ImageSequence &incident_data, + ImageSequence &empty_data, + TOFCorrectionsData &corrections_data, + bool apply_lorentz_correction, + bool apply_spherical_absorption_correction) { + double sample_linear_scattering_c = corrections_data.sample_number_density + * corrections_data.sample_scattering_x_section; + double incident_linear_scattering_c = + corrections_data.incident_number_density + * corrections_data.incident_scattering_x_section; + double sample_linear_absorption_c = corrections_data.sample_number_density + * corrections_data.sample_absorption_x_section; + double incident_linear_absorption_c = + corrections_data.incident_number_density + * corrections_data.incident_absorption_x_section; + + Detector detector = *experiment.get_detector(); + Scan scan = *experiment.get_scan(); + + std::shared_ptr beam_ptr = experiment.get_beam(); + std::shared_ptr beam = + std::dynamic_pointer_cast(beam_ptr); + DIALS_ASSERT(beam != nullptr); + + vec3 unit_s0 = beam->get_unit_s0(); + double sample_to_source_distance = beam->get_sample_to_source_distance(); + + scitbx::af::shared img_tof = scan.get_property("time_of_flight"); + + int n_panels = detector.size(); + int num_images = data.size(); + vec2 image_size = detector[0].get_image_size(); + DIALS_ASSERT(num_images == img_tof.size()); + + boost::python::dict d; + af::reflection_table i_reflection_table = + dxtbx::af::flex_table_suite::deepcopy(reflection_table, d); + af::reflection_table e_reflection_table = + dxtbx::af::flex_table_suite::deepcopy(reflection_table, d); + + ShoeboxProcessor shoebox_processor( + reflection_table, n_panels, 0, num_images, false); + + ShoeboxProcessor incident_shoebox_processor( + i_reflection_table, n_panels, 0, num_images, false); + + ShoeboxProcessor empty_shoebox_processor( + e_reflection_table, n_panels, 0, num_images, false); + + for (std::size_t img_num = 0; img_num < num_images; ++img_num) { + dxtbx::format::Image img = data.get_corrected_data(img_num); + dxtbx::format::Image mask = data.get_mask(img_num); + + af::shared > > output_data( + n_panels); + af::shared > > output_mask( + n_panels); + + for (std::size_t i = 0; i < output_data.size(); ++i) { + output_data[i] = img.tile(i).data(); + output_mask[i] = mask.tile(i).data(); + } + shoebox_processor.next_data_only( + model::Image(output_data.const_ref(), output_mask.const_ref())); + } + + for (std::size_t img_num = 0; img_num < num_images; ++img_num) { + dxtbx::format::Image img = incident_data.get_corrected_data(img_num); + dxtbx::format::Image mask = incident_data.get_mask(img_num); + + af::shared > > output_data( + n_panels); + af::shared > > output_mask( + n_panels); + + for (std::size_t i = 0; i < output_data.size(); ++i) { + output_data[i] = img.tile(i).data(); + output_mask[i] = mask.tile(i).data(); + } + incident_shoebox_processor.next_data_only( + model::Image(output_data.const_ref(), output_mask.const_ref())); + } + + for (std::size_t img_num = 0; img_num < num_images; ++img_num) { + dxtbx::format::Image img = empty_data.get_corrected_data(img_num); + dxtbx::format::Image mask = empty_data.get_mask(img_num); + + af::shared > > output_data( + n_panels); + af::shared > > output_mask( + n_panels); + + for (std::size_t i = 0; i < output_data.size(); ++i) { + output_data[i] = img.tile(i).data(); + output_mask[i] = mask.tile(i).data(); + } + empty_shoebox_processor.next_data_only( + model::Image(output_data.const_ref(), output_mask.const_ref())); + } + + af::shared > shoeboxes = reflection_table["shoebox"]; + af::shared > e_shoeboxes = i_reflection_table["shoebox"]; + af::shared > i_shoeboxes = e_reflection_table["shoebox"]; + af::const_ref bboxes = reflection_table["bbox"]; + + for (std::size_t i = 0; i < reflection_table.size(); ++i) { + Shoebox<> shoebox = shoeboxes[i]; + Shoebox<> i_shoebox = e_shoeboxes[i]; + Shoebox<> e_shoebox = i_shoeboxes[i]; + int6 bbox = bboxes[i]; + int panel = shoebox.panel; + + for (std::size_t z = 0; z < shoebox.zsize(); ++z) { + int frame_z = bbox[4] + z; + double tof = img_tof[frame_z] * std::pow(10, -6); // (s) + + for (std::size_t x = 0; x < shoebox.xsize(); ++x) { + int panel_x = bbox[0] + x; + if (panel_x > image_size[0] || panel_x < 0) { + continue; + } + for (std::size_t y = 0; y < shoebox.ysize(); ++y) { + int panel_y = bbox[2] + y; + if (panel_y > image_size[1] || panel_y < 0) { + continue; + } + + double incident_pixel_data = i_shoebox.data(z, y, x); + double empty_pixel_data = e_shoebox.data(z, y, x); + double pixel_data = shoebox.data(z, y, x); + /* + if (i==0){ + std::cout<<"TEST pixel data " << pixel_data << " empty data " << + empty_pixel_data << " incident data " << incident_pixel_data << std::endl; + } + */ + + pixel_data /= corrections_data.sample_proton_charge; + incident_pixel_data /= corrections_data.incident_proton_charge; + empty_pixel_data /= corrections_data.empty_proton_charge; + /* + if (i==0){ + std::cout<<"TEST after proton charge pixel data " << pixel_data << + " empty data " << empty_pixel_data << " incident data " << + incident_pixel_data << std::endl; + } + */ + + // Subtract empty from incident and sample + pixel_data -= empty_pixel_data; + incident_pixel_data -= empty_pixel_data; + /* + if (i==0){ + std::cout<<"TEST after empty subtract pixel data " << pixel_data << + " empty data " << empty_pixel_data << " incident data " << + incident_pixel_data << std::endl; + } + */ + + double two_theta = detector[panel].get_two_theta_at_pixel( + unit_s0, scitbx::vec2(panel_x, panel_y)); + double two_theta_deg = two_theta * (180 / pi); + int two_theta_idx = static_cast(two_theta_deg / 10); + + scitbx::vec3 s1 = detector[panel].get_pixel_lab_coord( + scitbx::vec2(panel_x, panel_y)); + double distance = s1.length() + sample_to_source_distance; + distance *= std::pow(10, -3); // (m) + double wl = ((Planck * tof) / (m_n * (distance))) * std::pow(10, 10); + + // Spherical absorption correction + if (apply_spherical_absorption_correction) { + double sample_muR = + (sample_linear_scattering_c + (sample_linear_absorption_c / 1.8) * wl) + * corrections_data.sample_radius; + double sample_absorption_correction = + tof_pixel_spherical_absorption_correction( + pixel_data, sample_muR, two_theta, two_theta_idx); + + /* + if (i==0){ + std::cout<<"TEST sample absorption correction: " << + sample_absorption_correction << std::endl; + } + */ + if (sample_absorption_correction < 1e-5) { + shoebox.data(z, y, x) = 0; + continue; + } + + double incident_muR = (incident_linear_scattering_c + + (incident_linear_absorption_c / 1.8) * wl) + * corrections_data.incident_radius; + double incident_absorption_correction = + tof_pixel_spherical_absorption_correction( + pixel_data, incident_muR, two_theta, two_theta_idx); + + /* + if (i==0){ + std::cout<<"TEST incident absorption correction: " << + incident_absorption_correction << std::endl; + } + */ + if (incident_absorption_correction < 1e-5) { + shoebox.data(z, y, x) = 0; + continue; + } + + incident_pixel_data /= incident_absorption_correction; + /* + if (i==0){ + std::cout<<"TEST incident after correction: " << + incident_pixel_data << std::endl; + } + */ + if (incident_pixel_data < 1e-5) { + shoebox.data(z, y, x) = 0; + continue; + } + + pixel_data /= incident_pixel_data; + /* + if (i==0){ + std::cout<<"TEST pixel data after normalisation: " << pixel_data + << std::endl; + } + */ + pixel_data /= sample_absorption_correction; + /* + if (i==0){ + std::cout<<"TEST pixel data after correction: " << pixel_data << + std::endl; + } + */ + } + + else { + if (incident_pixel_data < 1e-5) { + shoebox.data(z, y, x) = 0; + continue; + } + + pixel_data /= incident_pixel_data; + } + + // Lorentz correction + if (apply_lorentz_correction) { + double sin_two_theta_sq = std::pow(sin(two_theta * .5), 2); + double lorentz_correction = sin_two_theta_sq / std::pow(wl, 4); + // std::cout<<"TEST lorentz two_theta " << two_theta << " sin_two_theta_sq + // " << sin_two_theta_sq << " correction " << lorentz_correction << " + // pixel data " << pixel_data << " empty " << empty_pixel_data << + // std::endl; + pixel_data *= lorentz_correction; + // if (i==0){ + // std::cout<<"TEST lorentz " << lorentz_correction << " " << + // pixel_data << std::endl; + // } + } + + shoebox.data(z, y, x) = double(pixel_data); + } + } + } + } + } + + void tof_calculate_shoebox_foreground(af::reflection_table &reflection_table, + Experiment &experiment, + double foreground_radius) { + af::shared > shoeboxes = reflection_table["shoebox"]; + Scan scan = *experiment.get_scan(); + Detector detector = *experiment.get_detector(); + Goniometer goniometer = *experiment.get_goniometer(); + mat3 setting_rotation = goniometer.get_setting_rotation(); + + std::shared_ptr beam_ptr = experiment.get_beam(); + std::shared_ptr beam = + std::dynamic_pointer_cast(beam_ptr); + DIALS_ASSERT(beam != nullptr); + vec3 unit_s0 = beam->get_unit_s0(); + double sample_to_source_distance = beam->get_sample_to_source_distance(); + + scitbx::af::shared img_tof = scan.get_property("time_of_flight"); + af::const_ref bboxes = reflection_table["bbox"]; + scitbx::af::shared > rlps = reflection_table["rlp"]; + + for (std::size_t i = 0; i < reflection_table.size(); ++i) { + Shoebox<> shoebox = shoeboxes[i]; + af::ref > mask = shoebox.mask.ref(); + int panel = shoebox.panel; + int6 bbox = bboxes[i]; + vec3 rlp = rlps[i]; + + for (std::size_t z = 0; z < shoebox.zsize(); ++z) { + int frame_z = bbox[4] + z; + double tof = img_tof[frame_z] * std::pow(10, -6); // (s) + + for (std::size_t y = 0; y < shoebox.ysize(); ++y) { + int panel_y = bbox[2] + y; + + for (std::size_t x = 0; x < shoebox.xsize(); ++x) { + int panel_x = bbox[0] + x; + vec3 s1 = + detector[panel].get_pixel_lab_coord(vec2(panel_x, panel_y)); + double pixel_distance = s1.length() + sample_to_source_distance; + pixel_distance *= std::pow(10, -3); // (m) + double wl = + ((Planck * tof) / (m_n * (pixel_distance))) * std::pow(10, 10); // (A) + vec3 s0 = unit_s0 / wl; + s1 = s1 / s1.length() * (1 / wl); + vec3 S = s1 - s0; + S = setting_rotation.inverse() * S; + double distance = (S - rlp).length(); + + int mask_value = (distance <= foreground_radius) ? Foreground : Background; + mask(z, y, x) |= mask_value; + } + } + } + } + } + +}} // namespace dials::algorithms + +#endif /* DIALS_ALGORITHMS_SCALING_TOF_SCALING_CORRECTIONS_H */ \ No newline at end of file diff --git a/src/dials/command_line/tof_integrate.py b/src/dials/command_line/tof_integrate.py new file mode 100644 index 0000000000..ee32bfcdc2 --- /dev/null +++ b/src/dials/command_line/tof_integrate.py @@ -0,0 +1,488 @@ +# LIBTBX_SET_DISPATCHER_NAME dials.tof_integrate +from __future__ import annotations + +import logging +import multiprocessing + +import numpy as np + +import cctbx.array_family.flex +import libtbx + +import dials.util.log +from dials.algorithms.integration.fit.tof_line_profile import ( + compute_line_profile_intensity, +) +from dials.algorithms.integration.report import IntegrationReport +from dials.algorithms.profile_model.gaussian_rs import Model as GaussianRSProfileModel +from dials.algorithms.profile_model.gaussian_rs.calculator import ( + ComputeEsdBeamDivergence, +) +from dials.algorithms.shoebox import MaskCode +from dials.array_family import flex +from dials.command_line.integrate import process_reference +from dials.extensions.simple_background_ext import SimpleBackgroundExt +from dials.util.options import ArgumentParser, reflections_and_experiments_from_files +from dials.util.phil import parse +from dials.util.version import dials_version +from dials_algorithms_integration_integrator_ext import ( + TOFCorrectionsData, + tof_calculate_shoebox_foreground, + tof_extract_shoeboxes_to_reflection_table, +) + +logger = logging.getLogger("dials.command_line.simple_integrate") + +phil_scope = parse( + """ +input{ + incident_run = None + .type = str + .help = "Path to incident run to normalize intensities (e.g. Vanadium)." + empty_run = None + .type = str + .help = "Path to empty run to correct empty counts." +} +output { +experiments = 'integrated.expt' + .type = str + .help = "The experiments output filename" +reflections = 'integrated.refl' + .type = str + .help = "The integrated output filename" +output_hkl = True + .type = bool + .help = "Output the integrated intensities as a SHELX hkl file" +hkl = 'integrated.hkl' + .type = str + .help = "The hkl output filename" +phil = 'dials.simple_integrate.phil' + .type = str + .help = "The output phil file" +log = 'tof_integrate.log' + .type = str + .help = "The log filename" +} +method{ +line_profile_fitting = False + .type = bool + .help = "Use integration by profile fitting using a Gaussian" + "convoluted with back-to-back exponential functions" +} +corrections{ + lorentz = True + .type = bool + .help = "Apply the Lorentz correction to target spectrum." + spherical_absorption = True + .type = bool + .help = "Apply a spherical absorption correction." +} +incident_spectrum{ + sample_number_density = 0.0722 + .type = float + .help = "Sample number density for incident run." + "Default is Vanadium used at SXD" + sample_radius = 0.3 + .type = float + .help = "Sample radius incident run." + "Default is Vanadium used at SXD" + scattering_x_section = 5.158 + .type = float + .help = "Sample scattering cross section used for incident run." + "Default is Vanadium used at SXD" + absorption_x_section = 4.4883 + .type = float + .help = "Sample absorption cross section for incident run." + "Default is Vanadium used at SXD" +} +target_spectrum{ + sample_number_density = None + .type = float + .help = "Sample number density for target run." + sample_radius = None + .type = float + .help = "Sample radius target run." + scattering_x_section = None + .type = float + .help = "Sample scattering cross section used for target run." + absorption_x_section = None + .type = float + .help = "Sample absorption cross section for target run." +} +mp{ + nproc = Auto + .type = int(value_min=1) + .help = "Number of processors to use during parallelized steps." + "If set to Auto DIALS will choose automatically." +} +sigma_b = 0.01 + .type = float + .help = "Used to calculate xy bounding box of predicted reflections" +sigma_m = 3 + .type = float + .help = "Used to calculate z bounding box of predicted reflections" +keep_shoeboxes = False + .type = bool + .help = "Retain shoeboxes in output reflection table" +""" +) + +""" +Kabsch 2010 refers to +Kabsch W., Integration, scaling, space-group assignment and +post-refinment, Acta Crystallographica Section D, 2010, D66, 133-144 +Usage: +$ dials.tof_integrate.py refined.expt refined.refl +""" + + +def output_reflections_as_hkl(reflections, filename): + def get_corrected_intensity_and_variance(reflections, idx): + intensity = reflections["intensity.sum.value"][idx] + variance = reflections["intensity.sum.variance"][idx] + return intensity, variance + + def get_line_profile_intensity_and_variance(reflections, idx): + intensity = reflections["intensity.prf.value"][idx] + variance = reflections["intensity.prf.variance"][idx] + return intensity, variance + + def valid_intensity(intensity, variance): + from math import isinf, isnan + + if isnan(intensity) or isinf(intensity): + return False + if isnan(variance) or isinf(variance): + return False + return intensity > 0 and variance > 0 + + with open(filename, "w") as g: + for i in range(len(reflections)): + h, k, l = reflections["miller_index"][i] + batch_number = 1 + intensity, variance = get_corrected_intensity_and_variance(reflections, i) + if not valid_intensity(intensity, variance): + continue + intensity = round(intensity, 2) + sigma = round(np.sqrt(variance), 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(f"{int(h)} {int(k)} {int(l)} {float(intensity)} {float(sigma)} {int(batch_number)} {float(wavelength)}\n") + 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) + ) + ) + # g.write(f"{int(0)} {int(0)} {int(0)} {float(0.00)} {float(0.00)} {int(0)} {float(0.0000)}") + + if "intensity.prf.value" not in reflections: + return + with open("prf_" + filename, "w") as g: + for i in range(len(reflections)): + h, k, l = reflections["miller_index"][i] + batch_number = 1 + intensity, variance = get_line_profile_intensity_and_variance( + reflections, i + ) + if not valid_intensity(intensity, variance): + continue + intensity = round(intensity, 2) + sigma = round(np.sqrt(variance), 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 split_reflections(reflections, n, by_panel=False): + if by_panel: + for i in range(max(reflections["panel"]) + 1): + sel = reflections["panel"] == i + yield reflections.select(sel) + else: + d, r = divmod(len(reflections), n) + for i in range(n): + si = (d + 1) * (i if i < r else r) + d * (0 if i < r else i - r) + yield reflections[si : si + (d + 1 if i < r else d)] + + +def join_reflections(list_of_reflections): + reflections = list_of_reflections[0] + for i in range(1, len(list_of_reflections)): + reflections.extend(list_of_reflections[i]) + return reflections + + +def run(): + + """ + Input setup + """ + + phil = phil_scope.fetch() + + usage = "usage: dials.tof_integrate.py refined.expt refined.refl" + parser = ArgumentParser( + usage=usage, + phil=phil, + epilog=__doc__, + read_experiments=True, + read_reflections=True, + ) + + params, options = parser.parse_args(args=None, show_diff_phil=False) + + dials.util.log.config(verbosity=options.verbose, logfile=params.output.log) + logger.info(dials_version()) + + """ + Load experiment and reflections + """ + + reflections, experiments = reflections_and_experiments_from_files( + params.input.reflections, params.input.experiments + ) + reflections = reflections[0] + + integrated_reflections = run_integrate(params, experiments, reflections) + integrated_reflections.as_msgpack_file(params.output.reflections) + experiments.as_file(params.output.experiments) + if params.output.output_hkl: + output_reflections_as_hkl(integrated_reflections, params.output.hkl) + + +def run_integrate(params, experiments, reflections): + nproc = params.mp.nproc + if nproc is libtbx.Auto: + nproc = multiprocessing.cpu_count() + + reflections, _ = process_reference(reflections) + + """ + Predict reflections using experiment crystal + """ + + min_s0_idx = min( + range(len(reflections["wavelength"])), key=reflections["wavelength"].__getitem__ + ) + min_s0 = reflections["s0"][min_s0_idx] + dmin = None + for experiment in experiments: + expt_dmin = experiment.detector.get_max_resolution(min_s0) + if dmin is None or expt_dmin < dmin: + dmin = expt_dmin + + predicted_reflections = None + for idx, experiment in enumerate(experiments): + + if predicted_reflections is None: + predicted_reflections = flex.reflection_table.from_predictions( + experiment, padding=1.0, dmin=dmin + ) + predicted_reflections["id"] = cctbx.array_family.flex.int( + len(predicted_reflections), idx + ) + predicted_reflections["imageset_id"] = cctbx.array_family.flex.int( + len(predicted_reflections), idx + ) + else: + r = flex.reflection_table.from_predictions( + experiment, padding=1.0, dmin=dmin + ) + r["id"] = cctbx.array_family.flex.int(len(r), idx) + r["imageset_id"] = cctbx.array_family.flex.int(len(r), idx) + predicted_reflections.extend(r) + predicted_reflections.calculate_entering_flags(experiments) + + for i in range(len(experiments)): + predicted_reflections.experiment_identifiers()[i] = experiments[i].identifier + + # Updates flags to set which reflections to use in generating reference profiles + matched, reflections, unmatched = predicted_reflections.match_with_reference( + reflections + ) + # if "idx" in reflections: + # predicted_reflections["idx"] = reflections["idx"] + # sel = predicted_reflections.get_flags(predicted_reflections.flags.reference_spot) + predicted_reflections = predicted_reflections.select(matched) + + """ + Create profile model and add it to experiment. + This is used to predict reflection properties. + """ + + # Filter reflections to use to create the model + used_in_ref = reflections.get_flags(reflections.flags.used_in_refinement) + model_reflections = reflections.select(used_in_ref) + + sigma_b = ComputeEsdBeamDivergence( + experiment.detector, model_reflections, centroid_definition="s1" + ).sigma() + + # sigma_m in 3.1 of Kabsch 2010 + sigma_m = params.sigma_m + # sigma_b = 0.001 + # The Gaussian model given in 2.3 of Kabsch 2010 + for idx, experiment in enumerate(experiments): + experiments[idx].profile = GaussianRSProfileModel( + params=params, n_sigma=2.5, sigma_b=sigma_b, sigma_m=sigma_m + ) + + """ + Compute properties for predicted reflections using profile model, + accessed via experiment.profile_model. These reflection_table + methods are largely just wrappers for profile_model.compute_bbox etc. + Note: I do not think all these properties are needed for integration, + but are all present in the current dials.integrate output. + """ + + predicted_reflections.compute_bbox(experiments) + x1, x2, y1, y2, t1, t2 = predicted_reflections["bbox"].parts() + predicted_reflections = predicted_reflections.select( + (t1 > 0) + & (t2 < max([e.scan.get_image_range()[1] for e in experiments])) + & (x1 > 0) + & (y1 > 0) + ) + + predicted_reflections.compute_d(experiments) + # predicted_reflections.compute_partiality(experiments) + + # Shoeboxes + print("Getting shoebox data") + predicted_reflections["shoebox"] = flex.shoebox( + predicted_reflections["panel"], + predicted_reflections["bbox"], + allocate=False, + flatten=False, + ) + + experiment_cls = experiments[0].imageset.get_format_class() + incident_fmt_class = experiment_cls.get_instance(params.input.incident_run) + empty_fmt_class = experiment_cls.get_instance(params.input.empty_run) + + incident_data = experiment_cls(params.input.incident_run).get_imageset( + params.input.incident_run + ) + empty_data = experiment_cls(params.input.empty_run).get_imageset( + params.input.empty_run + ) + incident_proton_charge = incident_fmt_class.get_proton_charge() + empty_proton_charge = empty_fmt_class.get_proton_charge() + + predicted_reflections.map_centroids_to_reciprocal_space( + experiments, calculated=True + ) + + for expt in experiments: + + expt_proton_charge = experiment_cls.get_instance( + expt.imageset.paths()[0], **expt.imageset.data().get_params() + ).get_proton_charge() + expt_data = expt.imageset + corrections_data = TOFCorrectionsData( + expt_proton_charge, + incident_proton_charge, + empty_proton_charge, + params.target_spectrum.sample_radius, + params.target_spectrum.scattering_x_section, + params.target_spectrum.absorption_x_section, + params.target_spectrum.sample_number_density, + params.incident_spectrum.sample_radius, + params.incident_spectrum.scattering_x_section, + params.incident_spectrum.absorption_x_section, + params.incident_spectrum.sample_number_density, + ) + + tof_extract_shoeboxes_to_reflection_table( + predicted_reflections, + expt, + expt_data, + incident_data, + empty_data, + corrections_data, + params.corrections.lorentz, + params.corrections.spherical_absorption, + ) + + tof_calculate_shoebox_foreground(predicted_reflections, expt, 0.5) + predicted_reflections.is_overloaded(experiments) + predicted_reflections.contains_invalid_pixels() + predicted_reflections["partiality"] = flex.double(len(predicted_reflections), 1.0) + + # Background calculated explicitly to expose underlying algorithm + background_algorithm = SimpleBackgroundExt(params=None, experiments=experiments) + success = background_algorithm.compute_background(predicted_reflections) + predicted_reflections.set_flags( + ~success, predicted_reflections.flags.failed_during_background_modelling + ) + + # Centroids calculated explicitly to expose underlying algorithm + # centroid_algorithm = SimpleCentroidExt(params=None, experiments=experiments) + # centroid_algorithm.compute_centroid(predicted_reflections) + + print("Calculating summed intensities") + predicted_reflections.compute_summed_intensity() + + if params.method.line_profile_fitting: + print("Calculating line profile fitted intensities") + predicted_reflections = compute_line_profile_intensity(predicted_reflections) + + # Filter reflections with a high fraction of masked foreground + valid_foreground_threshold = 1.0 # DIALS default + sboxs = predicted_reflections["shoebox"] + nvalfg = sboxs.count_mask_values(MaskCode.Valid | MaskCode.Foreground) + nforeg = sboxs.count_mask_values(MaskCode.Foreground) + fraction_valid = nvalfg.as_double() / nforeg.as_double() + selection = fraction_valid < valid_foreground_threshold + predicted_reflections.set_flags( + selection, predicted_reflections.flags.dont_integrate + ) + + predicted_reflections["num_pixels.valid"] = sboxs.count_mask_values(MaskCode.Valid) + predicted_reflections["num_pixels.background"] = sboxs.count_mask_values( + MaskCode.Valid | MaskCode.Background + ) + predicted_reflections["num_pixels.background_used"] = sboxs.count_mask_values( + MaskCode.Valid | MaskCode.Background | MaskCode.BackgroundUsed + ) + predicted_reflections["num_pixels.foreground"] = nvalfg + + integration_report = IntegrationReport(experiments, predicted_reflections) + logger.info("") + logger.info(integration_report.as_str(prefix=" ")) + + if not params.keep_shoeboxes: + del predicted_reflections["shoebox"] + return predicted_reflections + + +if __name__ == "__main__": + run()