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 8540a2e119..3b47cd7223 100644 --- a/src/dials/algorithms/integration/boost_python/integration_integrator_ext.cc +++ b/src/dials/algorithms/integration/boost_python/integration_integrator_ext.cc @@ -14,6 +14,7 @@ #include #include #include +#include using namespace boost::python; @@ -349,6 +350,36 @@ namespace dials { namespace algorithms { namespace boost_python { def("max_memory_needed", &max_memory_needed, (arg("reflection table"), arg("frame0"), arg("frame1"), arg("flatten"))); + + def("tof_extract_shoeboxes_to_reflection_table", + &tof_extract_shoeboxes_to_reflection_table, + (arg("reflection_table"), + arg("experiment"), + arg("incident_data"), + arg("empty_data"), + 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"))); + + def("tof_extract_shoeboxes_to_reflection_table2", + &tof_extract_shoeboxes_to_reflection_table2, + (arg("reflection_table"), + arg("experiment"), + arg("incident_data"), + arg("empty_data"), + 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"))); } }}} // namespace dials::algorithms::boost_python diff --git a/src/dials/algorithms/profile_model/gaussian_rs/__init__.py b/src/dials/algorithms/profile_model/gaussian_rs/__init__.py index 1002d80da4..d20c7a7f13 100644 --- a/src/dials/algorithms/profile_model/gaussian_rs/__init__.py +++ b/src/dials/algorithms/profile_model/gaussian_rs/__init__.py @@ -6,6 +6,7 @@ BBoxCalculator2D, BBoxCalculator3D, BBoxCalculatorIface, + BBoxCalculatorTOF, BBoxMultiCalculator, CoordinateSystem, CoordinateSystem2d, @@ -27,6 +28,7 @@ "BBoxCalculator", "BBoxCalculator2D", "BBoxCalculator3D", + "BBoxCalculatorTOF", "BBoxCalculatorIface", "BBoxMultiCalculator", "CoordinateSystem", @@ -52,7 +54,9 @@ def BBoxCalculator(crystal, beam, detector, goniometer, scan, delta_b, delta_m): """Return the relevant bbox calculator.""" - if goniometer is None or scan is None or scan.is_still(): + if scan.has_property("time_of_flight"): + algorithm = BBoxCalculatorTOF(beam, detector, scan, delta_b, delta_m) + elif goniometer is None or scan is None or scan.is_still(): algorithm = BBoxCalculator2D(beam, detector, delta_b, delta_m) else: algorithm = BBoxCalculator3D(beam, detector, goniometer, scan, delta_b, delta_m) diff --git a/src/dials/algorithms/profile_model/gaussian_rs/bbox_calculator.h b/src/dials/algorithms/profile_model/gaussian_rs/bbox_calculator.h index f6d8fe1486..42ed8002f2 100644 --- a/src/dials/algorithms/profile_model/gaussian_rs/bbox_calculator.h +++ b/src/dials/algorithms/profile_model/gaussian_rs/bbox_calculator.h @@ -33,6 +33,7 @@ namespace dials { using dxtbx::model::BeamBase; using dxtbx::model::Detector; using dxtbx::model::Goniometer; + using dxtbx::model::PolychromaticBeam; using dxtbx::model::Scan; using scitbx::vec2; using scitbx::vec3; @@ -349,6 +350,139 @@ namespace dials { double delta_divergence_; }; + class BBoxCalculatorTOF { + public: + /** + * Initialise the bounding box calculation. + * @param beam The beam parameters + * @param detector The detector parameters + * @param delta_divergence The xds delta_divergence parameter + * @param delta_mosaicity The xds delta_mosaicity parameter + */ + BBoxCalculatorTOF(const PolychromaticBeam &beam, + const Detector &detector, + const Scan &scan, + double delta_divergence, + double delta_mosaicity) + : detector_(detector), + scan_(scan), + beam_(beam), + delta_divergence_(delta_divergence), + delta_mosaicity_(delta_mosaicity) { + DIALS_ASSERT(delta_divergence > 0.0); + DIALS_ASSERT(delta_mosaicity >= 0.0); + } + + /** + * Calculate the bbox on the detector image volume for the reflection. + * + * The roi is calculated using the parameters delta_divergence and + * delta_mosaicity. The reflection mask comprises all pixels where: + * |e1| <= delta_d, |e2| <= delta_d, |e3| <= delta_m + * + * We transform the coordinates of the box + * (-delta_d, -delta_d, 0) + * (+delta_d, -delta_d, 0) + * (-delta_d, +delta_d, 0) + * (+delta_d, +delta_d, 0) + * + * to the detector image volume and return the minimum and maximum values + * for the x, y, z image volume coordinates. + * + * @param s1 The diffracted beam vector + * @param frame The predicted frame number + * @returns A 6 element array: (minx, maxx, miny, maxy, minz, maxz) + */ + virtual int6 single(vec3 s0, + vec3 s1, + double frame, + double L1, + std::size_t panel) const { + // Ensure our values are ok + DIALS_ASSERT(s1.length_sq() > 0); + + // Create the coordinate system for the reflection + CoordinateSystemTOF xcs(s0, s1, L1); + + // Get the divergence and mosaicity for this point + double delta_d = delta_divergence_; + double delta_m = delta_mosaicity_; + + // Calculate the beam vectors at the following xds coordinates: + // (-delta_d, -delta_d, 0) + // (+delta_d, -delta_d, 0) + // (-delta_d, +delta_d, 0) + // (+delta_d, +delta_d, 0) + double point = delta_d; + double3 sdash1 = xcs.to_beam_vector(double2(-point, -point)); + double3 sdash2 = xcs.to_beam_vector(double2(+point, -point)); + double3 sdash3 = xcs.to_beam_vector(double2(-point, +point)); + double3 sdash4 = xcs.to_beam_vector(double2(+point, +point)); + + // Get the detector coordinates (px) at the ray intersections + double2 xy1 = detector_[panel].get_ray_intersection_px(sdash1); + double2 xy2 = detector_[panel].get_ray_intersection_px(sdash2); + double2 xy3 = detector_[panel].get_ray_intersection_px(sdash3); + double2 xy4 = detector_[panel].get_ray_intersection_px(sdash4); + + // Return the roi in the following form: + // (minx, maxx, miny, maxy, minz, maxz) + // Min's are rounded down to the nearest integer, Max's are rounded up + double4 x(xy1[0], xy2[0], xy3[0], xy4[0]); + double4 y(xy1[1], xy2[1], xy3[1], xy4[1]); + + int x0 = (int)floor(min(x)); + int x1 = (int)ceil(max(x)); + int y0 = (int)floor(min(y)); + int y1 = (int)ceil(max(y)); + + double z0 = frame - delta_m; + if (z0 < 0) { + z0 = 0; + } + double max_z = scan_.get_array_range()[1]; + double z1 = frame + delta_m; + if (z1 > max_z) { + z1 = max_z; + } + + int6 bbox(x0, x1, y0, y1, z0, z1); + DIALS_ASSERT(bbox[1] > bbox[0]); + DIALS_ASSERT(bbox[3] > bbox[2]); + DIALS_ASSERT(bbox[5] > bbox[4]); + return bbox; + } + + /** + * Calculate the rois for an array of reflections given by the array of + * diffracted beam vectors and rotation angles. + * @param s1 The array of diffracted beam vectors + * @param phi The array of rotation angles. + */ + virtual af::shared array(const af::const_ref > &s0, + const af::const_ref > &s1, + const af::const_ref &frame, + const af::const_ref &L1, + const af::const_ref &panel) const { + DIALS_ASSERT(s1.size() == frame.size()); + DIALS_ASSERT(s1.size() == panel.size()); + DIALS_ASSERT(s0.size() == frame.size()); + DIALS_ASSERT(s0.size() == panel.size()); + af::shared result(s1.size(), af::init_functor_null()); + for (std::size_t i = 0; i < s1.size(); ++i) { + result[i] = single(s0[i], s1[i], frame[i], L1[i], panel[i]); + } + return result; + } + + private: + Detector detector_; + Scan scan_; + PolychromaticBeam beam_; + double delta_divergence_; + double delta_mosaicity_; + }; + /** * Class to help compute bbox for multiple experiments. */ diff --git a/src/dials/algorithms/profile_model/gaussian_rs/boost_python/gaussian_rs_ext.cc b/src/dials/algorithms/profile_model/gaussian_rs/boost_python/gaussian_rs_ext.cc index 2568531f27..71948ed573 100644 --- a/src/dials/algorithms/profile_model/gaussian_rs/boost_python/gaussian_rs_ext.cc +++ b/src/dials/algorithms/profile_model/gaussian_rs/boost_python/gaussian_rs_ext.cc @@ -212,6 +212,21 @@ namespace dials { arg("delta_divergence"), arg("delta_mosaicity")))); + class_("BBoxCalculatorTOF", no_init) + .def( + init( + (arg("beam"), + arg("detector"), + arg("scan"), + arg("delta_divergence"), + arg("delta_mosaicity")))) + .def("__call__", + &BBoxCalculatorTOF::single, + (arg("s0"), arg("s1"), arg("frame"), arg("L1"), arg("panel"))) + .def("__call__", + &BBoxCalculatorTOF::array, + (arg("s0"), arg("s1"), arg("frame"), arg("L1"), arg("panel"))); + class_("BBoxMultiCalculator") .def("append", &BBoxMultiCalculator::push_back) .def("__len__", &BBoxMultiCalculator::size) diff --git a/src/dials/algorithms/profile_model/gaussian_rs/coordinate_system.h b/src/dials/algorithms/profile_model/gaussian_rs/coordinate_system.h index 99f3c1194c..6afdec100b 100644 --- a/src/dials/algorithms/profile_model/gaussian_rs/coordinate_system.h +++ b/src/dials/algorithms/profile_model/gaussian_rs/coordinate_system.h @@ -374,6 +374,132 @@ namespace dials { double zeta_; }; + /** + * Class representing the local reflection coordinate system + */ + class CoordinateSystemTOF { + public: + /** + * Initialise coordinate system. s0 should be the same length as s1. + * These quantities are not checked because this class will be created for + * each reflection and we want to maximize performance. + * @param s0 The incident beam vector + * @param s1 The diffracted beam vector + */ + CoordinateSystemTOF(vec3 s0, vec3 s1, double L1) + : s0_(s0), + s1_(s1), + p_star_(s1 - s0), + L1_(L1), + e1_(s1.cross(s0).normalize()), + e2_(s1.cross(e1_).normalize()), + e3_((s1 + s0).normalize()) {} + + vec3 s0() const { + return s0_; + } + vec3 s1() const { + return s1_; + } + double L1() const { + return L1_; + } + vec3 p_star() const { + return p_star_; + } + vec3 e1_axis() const { + return e1_; + } + vec3 e2_axis() const { + return e2_; + } + vec3 e3_axis() const { + return e3_; + } + + /** + * Transform the beam vector to the reciprocal space coordinate system. + * @param s_dash The beam vector + * @param s0_dash The incident beam vector + * @returns The e1, e2, e3 coordinates + */ + vec2 from_beam_vector(const vec3 &s_dash) const { + double s1_length = s1_.length(); + double s0_length = s0_.length(); + DIALS_ASSERT(s1_length > 0); + DIALS_ASSERT(s0_length > 0); + // vec3 p_star0 = s_dash-s0_dash; + vec3 e1 = e1_ / s1_length; + vec3 e2 = e2_ / s1_length; + // vec3 e3 = (s1_+s0_)/(s1_length + s0_length); + /* + return vec3( + e1 * (s_dash - s1_), + e2 * (s_dash - s1_), + e3 * (p_star0 - p_star_)/p_star_.length()); + */ + return vec2(e1 * (s_dash - s1_), e2 * (s_dash - s1_)); + } + + /** + * Transform the reciprocal space coordinate to get the beam vector. + * @param c12 The e1 and e2 coordinates. + * @returns The beam vector + */ + vec3 to_beam_vector(const vec2 &c12) const { + double radius = s1_.length(); + DIALS_ASSERT(radius > 0); + vec3 scaled_e1 = e1_ * radius; + vec3 scaled_e2 = e2_ * radius; + vec3 normalized_s1 = s1_ / radius; + + vec3 p = c12[0] * scaled_e1 + c12[1] * scaled_e2; + double b = radius * radius - p.length_sq(); + DIALS_ASSERT(b >= 0); + double d = -(normalized_s1 * p) + std::sqrt(b); + return p + d * normalized_s1; + } + + /** + * @param c3 The XDS e3 coordinate + * @param s_dash The beam vector from the e1 and e2 coordinates + * @returns The wavelength of the e3 coordinate (s) + * + * Solved be rearranging c3 = e3(p_star0 - p_star) / s1_length, + * noting that p_star0 = s_dash - s0_dash, + * and s0_dash = unit_s0/wavelength_dash + */ + double to_wavelength(double c3, vec3 s_dash) const { + double p_star_length = p_star_.length(); + DIALS_ASSERT(p_star_length > 0); + vec3 unit_s0 = s0_.normalize(); + return (e3_ * unit_s0) / (e3_ * s_dash - e3_ * p_star_ - c3 * p_star_length); + } + + /** + * Transform the rotation angle to the reciprocal space coordinate system + * @param phi_dash The rotation angle + * @returns The e3 coordinate. + */ + double from_wavelength(double wavelength) const { + double p_star_length = p_star_.length(); + DIALS_ASSERT(p_star_length > 0); + vec3 scaled_e3 = e3_ / p_star_length; + vec3 s0 = s0_.normalize() / wavelength; + vec3 p_star0 = s1_ - s0; + return scaled_e3 * (p_star0 - p_star_); + } + + private: + vec3 s0_; + vec3 s1_; + vec3 p_star_; + vec3 e1_; + vec3 e2_; + vec3 e3_; + double L1_; + }; + }}}} // namespace dials::algorithms::profile_model::gaussian_rs #endif // DIALS_ALGORITHMS_PROFILE_MODEL_GAUSSIAN_RS_COORDINATE_SYSTEM_H diff --git a/src/dials/algorithms/profile_model/gaussian_rs/model.py b/src/dials/algorithms/profile_model/gaussian_rs/model.py index 710e974929..56fa71c75b 100644 --- a/src/dials/algorithms/profile_model/gaussian_rs/model.py +++ b/src/dials/algorithms/profile_model/gaussian_rs/model.py @@ -515,9 +515,20 @@ def compute_bbox( ) # Calculate the bounding boxes of all the reflections - bbox = calculate( - reflections["s1"], reflections["xyzcal.px"].parts()[2], reflections["panel"] - ) + if scan.has_property("time_of_flight"): + bbox = calculate( + reflections["s0_cal"], + reflections["s1"], + reflections["xyzcal.px"].parts()[2], + reflections["L1"], + reflections["panel"], + ) + else: + bbox = calculate( + reflections["s1"], + reflections["xyzcal.px"].parts()[2], + reflections["panel"], + ) # Return the bounding boxes return bbox 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..4a23a86bb3 --- /dev/null +++ b/src/dials/algorithms/scaling/tof_scaling_corrections.h @@ -0,0 +1,571 @@ + +#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 + +namespace dials { namespace algorithms { + + using dxtbx::ImageSequence; + using dxtbx::af::flex_table; + using dxtbx::model::Detector; + using dxtbx::model::Experiment; + using dxtbx::model::PolychromaticBeam; + using dxtbx::model::Scan; + using dxtbx::model::scan_property_types; + using scitbx::deg_as_rad; + 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; + } + + void tof_extract_shoeboxes_to_reflection_table2( + af::reflection_table &reflection_table, + Experiment &experiment, + ImageSequence &incident_data, + ImageSequence &empty_data, + 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) { + double sample_linear_scattering_c = + sample_number_density * sample_scattering_x_section; + double incident_linear_scattering_c = + incident_number_density * incident_scattering_x_section; + double sample_linear_absorption_c = + sample_number_density * sample_absorption_x_section; + double incident_linear_absorption_c = + incident_number_density * 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"); + ImageSequence data = + boost::python::extract(experiment.get_imageset()); + + int n_panels = detector.size(); + int num_images = data.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 > 63 || panel_x < 0) { + continue; + } + for (std::size_t y = 0; y < shoebox.ysize(); ++y) { + int panel_y = bbox[2] + y; + if (panel_y > 63 || panel_y < 0) { + continue; + } + + double incident_pixel_data = i_shoebox.data(x, y, z); + double empty_pixel_data = e_shoebox.data(x, y, z); + double pixel_data = shoebox.data(x, y, z); + + // Subtract empty from incident and sample + pixel_data -= empty_pixel_data; + incident_pixel_data -= empty_pixel_data; + + // Spherical absorption correction + 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); + + double sample_muR = + (sample_linear_scattering_c + (sample_linear_absorption_c / 1.8) * wl) + * sample_radius; + double sample_absorption_correction = + tof_pixel_spherical_absorption_correction( + pixel_data, sample_muR, two_theta, two_theta_idx); + + if (sample_absorption_correction < 1e-5) { + shoebox.data(x, y, z) = 0; + continue; + } + + double incident_muR = + (incident_linear_scattering_c + (incident_linear_absorption_c / 1.8) * wl) + * incident_radius; + double incident_absorption_correction = + tof_pixel_spherical_absorption_correction( + pixel_data, incident_muR, two_theta, two_theta_idx); + + if (incident_absorption_correction < 1e-5) { + shoebox.data(x, y, z) = 0; + continue; + } + + incident_pixel_data /= incident_absorption_correction; + if (incident_pixel_data < 1e-5) { + shoebox.data(x, y, z) = 0; + continue; + } + + pixel_data /= incident_pixel_data; + pixel_data /= sample_absorption_correction; + + // 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); + pixel_data *= lorentz_correction; + + shoebox.data(x, y, z) = (int)pixel_data; + } + } + } + } + } + + void tof_extract_shoeboxes_to_reflection_table(af::reflection_table &reflection_table, + Experiment &experiment, + ImageSequence &incident_data, + ImageSequence &empty_data, + 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) { + auto start = std::chrono::high_resolution_clock::now(); + + double sample_linear_scattering_c = + sample_number_density * sample_scattering_x_section; + double incident_linear_scattering_c = + incident_number_density * incident_scattering_x_section; + double sample_linear_absorption_c = + sample_number_density * sample_absorption_x_section; + double incident_linear_absorption_c = + incident_number_density * 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"); + + ImageSequence data = + boost::python::extract(experiment.get_imageset()); + int n_panels = detector.size(); + int num_images = data.size(); + DIALS_ASSERT(num_images == img_tof.size()); + + ShoeboxProcessor shoebox_processor( + reflection_table, n_panels, 0, num_images, false); + + auto block1_end = std::chrono::high_resolution_clock::now(); + auto duration_block1 = + std::chrono::duration_cast(block1_end - start); + std::cout << "Time for init: " << duration_block1.count() << " microseconds" + << std::endl; + + for (std::size_t img_num = 0; img_num < num_images; ++img_num) { + // Image for each panel + dxtbx::format::Image imgs = data.get_corrected_data(img_num); + dxtbx::format::Image masks = data.get_mask(img_num); + dxtbx::format::Image i_imgs = incident_data.get_corrected_data(img_num); + dxtbx::format::Image e_imgs = empty_data.get_corrected_data(img_num); + + DIALS_ASSERT(imgs.n_tiles() == i_imgs.n_tiles()); + DIALS_ASSERT(imgs.n_tiles() == e_imgs.n_tiles()); + + dxtbx::format::Image corrected_imgs; + double tof = img_tof[img_num] * std::pow(10, -6); // (s) + + for (std::size_t panel_num = 0; panel_num < imgs.n_tiles(); ++panel_num) { + // Get panel data + scitbx::af::versa > panel_data = + imgs.tile(panel_num).data(); + scitbx::af::versa > i_panel_data = + i_imgs.tile(panel_num).data(); + scitbx::af::versa > e_panel_data = + e_imgs.tile(panel_num).data(); + DIALS_ASSERT(panel_data.accessor().all_eq(i_panel_data.accessor())); + DIALS_ASSERT(panel_data.accessor().all_eq(e_panel_data.accessor())); + + scitbx::af::versa > corrected_img_data( + panel_data.accessor()); + + for (std::size_t i = 0; i < panel_data.accessor()[0]; ++i) { + for (std::size_t j = 0; j < panel_data.accessor()[1]; ++j) { + // Get data for pixel + double pixel_data = panel_data(i, j); + double incident_pixel_data = i_panel_data(i, j); + double empty_pixel_data = e_panel_data(i, j); + + // Subtract empty from incident and sample + pixel_data -= empty_pixel_data; + incident_pixel_data -= empty_pixel_data; + + // Spherical absorption correction + double two_theta = detector[panel_num].get_two_theta_at_pixel( + unit_s0, scitbx::vec2(i, j)); + double two_theta_deg = two_theta * (180 / pi); + int two_theta_idx = static_cast(two_theta_deg / 10); + + scitbx::vec3 s1 = + detector[panel_num].get_pixel_lab_coord(scitbx::vec2(i, j)); + 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); + + double sample_muR = + (sample_linear_scattering_c + (sample_linear_absorption_c / 1.8) * wl) + * sample_radius; + double sample_absorption_correction = + tof_pixel_spherical_absorption_correction( + pixel_data, sample_muR, two_theta, two_theta_idx); + if (sample_absorption_correction < 1e-5) { + corrected_img_data(i, j) = 0.0; + continue; + } + double incident_muR = + (incident_linear_scattering_c + (incident_linear_absorption_c / 1.8) * wl) + * incident_radius; + double incident_absorption_correction = + tof_pixel_spherical_absorption_correction( + pixel_data, incident_muR, two_theta, two_theta_idx); + + if (incident_absorption_correction < 1e-5) { + corrected_img_data(i, j) = 0.0; + continue; + } + + incident_pixel_data /= incident_absorption_correction; + if (incident_pixel_data < 1e-5) { + corrected_img_data(i, j) = 0.0; + continue; + } + + pixel_data /= incident_pixel_data; + pixel_data /= sample_absorption_correction; + + // 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); + pixel_data *= lorentz_correction; + corrected_img_data(i, j) = pixel_data; + } + } + corrected_imgs.push_back(dxtbx::format::ImageTile(corrected_img_data)); + } + + af::shared > > output_data( + corrected_imgs.n_tiles()); + af::shared > > output_mask( + corrected_imgs.n_tiles()); + for (std::size_t i = 0; i < output_data.size(); ++i) { + output_data[i] = corrected_imgs.tile(i).data(); + output_mask[i] = masks.tile(i).data(); + } + shoebox_processor.next_data_only( + model::Image(output_data.const_ref(), output_mask.const_ref())); + + auto block2_end = std::chrono::high_resolution_clock::now(); + auto duration_block2 = + std::chrono::duration_cast(block2_end - block1_end); + std::cout << "Time for shoebox assignment: " << duration_block2.count() + << " microseconds" << std::endl; + } + } +}} // namespace dials::algorithms + +#endif /* DIALS_ALGORITHMS_SCALING_TOF_SCALING_CORRECTIONS_H */ \ No newline at end of file diff --git a/src/dials/algorithms/spot_prediction/__init__.py b/src/dials/algorithms/spot_prediction/__init__.py index 2f3a82a242..bdb84e7cb0 100644 --- a/src/dials/algorithms/spot_prediction/__init__.py +++ b/src/dials/algorithms/spot_prediction/__init__.py @@ -34,6 +34,9 @@ "StillsReflectionPredictor", "LaueReflectionPredictor", ] +from dxtbx.model import tof_helpers + +import dials.array_family.flex as flex def ScanStaticReflectionPredictor(experiment, dmin=None, margin=1, padding=0, **kwargs): @@ -169,3 +172,95 @@ def LaueReflectionPredictor(experiment, dmin: float): experiment.crystal.get_space_group().type(), dmin, ) + + +class TOFReflectionPredictor: + def __init__(self, experiment, dmin): + self.experiment = experiment + self.dmin = dmin + self.predictor = dials_algorithms_spot_prediction_ext.LaueReflectionPredictor( + experiment.beam, + experiment.detector, + experiment.goniometer, + experiment.crystal.get_A(), + experiment.crystal.get_unit_cell(), + experiment.crystal.get_space_group().type(), + dmin, + ) + + def post_prediction(self, reflections): + + if "tof_cal" not in reflections: + reflections["tof_cal"] = flex.double(reflections.nrows()) + if "L1" not in reflections: + reflections["L1"] = flex.double(reflections.nrows()) + + tof_cal = flex.double(reflections.nrows()) + L1 = flex.double(reflections.nrows()) + L0 = self.experiment.beam.get_sample_to_source_distance() * 10**-3 # (m) + + panel_numbers = flex.size_t(reflections["panel"]) + expt = self.experiment + + for i_panel in range(len(expt.detector)): + sel = panel_numbers == i_panel + expt_reflections = reflections.select(sel) + x, y, _ = expt_reflections["xyzcal.mm"].parts() + s1 = expt.detector[i_panel].get_lab_coord(flex.vec2_double(x, y)) + expt_L1 = s1.norms() + expt_tof_cal = flex.double(expt_reflections.nrows()) + + for idx in range(len(expt_reflections)): + wavelength = expt_reflections[idx]["wavelength_cal"] + tof = tof_helpers.tof_from_wavelength( + wavelength, L0 + expt_L1[idx] * 10**-3 + ) + expt_tof_cal[idx] = tof + tof_cal.set_selected(sel, expt_tof_cal) + L1.set_selected(sel, expt_L1) + + reflections["tof_cal"] = tof_cal + reflections["L1"] = L1 + + # Filter out predicted reflections outside of experiment range + wavelength_range = expt.beam.get_wavelength_range() + sel = reflections["wavelength_cal"] >= wavelength_range[0] + reflections = reflections.select(sel) + sel = reflections["wavelength_cal"] <= wavelength_range[1] + reflections = reflections.select(sel) + + return reflections + + def for_ub(self, ub): + + reflection_table = self.predictor.for_ub(ub) + reflection_table = self.post_prediction(reflection_table) + + image_range = self.experiment.scan.get_image_range() + frames = list(range(image_range[0], image_range[1])) + x, y, z = reflection_table["xyzcal.px"].parts() + xyz = flex.vec3_double(len(reflection_table)) + + interpolation_tof = self.experiment.scan.get_property("time_of_flight") + interpolation_frames = list(range(len(interpolation_tof))) + tof_to_frame = tof_helpers.tof_to_frame_interpolator( + interpolation_tof, interpolation_frames + ) + L0 = self.experiment.beam.get_sample_to_source_distance() * 10**-3 # (m) + + for i in range(len(reflection_table)): + wavelength = reflection_table["wavelength_cal"][i] + L1 = reflection_table["L1"][i] + tof = ( + tof_helpers.tof_from_wavelength(wavelength, L0 + (L1 * 10**-3)) + * 10**6 + ) # (usec) + frame = tof_to_frame(tof) + xyz[i] = (x[i], y[i], float(frame)) + x, y, z = xyz.parts() + reflection_table["xyzcal.px"] = xyz + sel = z > frames[0] + return reflection_table.select(sel) + + def for_reflection_table(self, reflections, UB): + return self.predictor.for_reflection_table(reflections, UB) diff --git a/src/dials/algorithms/spot_prediction/reflection_predictor.py b/src/dials/algorithms/spot_prediction/reflection_predictor.py index 95058000ec..ed57540258 100644 --- a/src/dials/algorithms/spot_prediction/reflection_predictor.py +++ b/src/dials/algorithms/spot_prediction/reflection_predictor.py @@ -4,6 +4,7 @@ logger = logging.getLogger(__name__) +from dxtbx.model import ExperimentType from libtbx.phil import parse from dials.util import Sorry @@ -62,6 +63,7 @@ def __init__( ScanStaticReflectionPredictor, ScanVaryingReflectionPredictor, StillsReflectionPredictor, + TOFReflectionPredictor, ) from dials.array_family import flex @@ -79,6 +81,15 @@ def __call__(self): result.del_selected(mask) return result + if experiment.get_type() == ExperimentType.TOF: + predictor = TOFReflectionPredictor(experiment=experiment, dmin=dmin) + predict = Predictor( + "ToF prediction", + lambda: predictor.for_ub(experiment.crystal.get_A()), + ) + self._predict = predict + return + # Check prediction to maximum resolution is possible wl = experiment.beam.get_wavelength() if dmin is not None and dmin < 0.5 * wl: