From 273b300c426c13342a9b88d0f355cb76ea83dd0f Mon Sep 17 00:00:00 2001 From: davidmcdonagh Date: Wed, 15 Nov 2023 10:15:07 +0000 Subject: [PATCH] Added method to predict spots for asymmetric unit. --- .../algorithms/indexing/model_evaluation.py | 2 +- .../boost_python/reflection_predictor.cc | 5 +- .../spot_prediction/reflection_predictor.h | 73 +++++++++++++++++-- src/dials/command_line/dials_import.py | 5 -- 4 files changed, 70 insertions(+), 15 deletions(-) diff --git a/src/dials/algorithms/indexing/model_evaluation.py b/src/dials/algorithms/indexing/model_evaluation.py index 9959757e8e..b9dbdddeee 100644 --- a/src/dials/algorithms/indexing/model_evaluation.py +++ b/src/dials/algorithms/indexing/model_evaluation.py @@ -36,7 +36,7 @@ def filter_doubled_cell(solutions): accepted_solutions = [] for i1, s1 in enumerate(solutions): doubled_cell = False - for (m1, m2, m3) in ( + for m1, m2, m3 in ( (2, 1, 1), (1, 2, 1), (1, 1, 2), diff --git a/src/dials/algorithms/spot_prediction/boost_python/reflection_predictor.cc b/src/dials/algorithms/spot_prediction/boost_python/reflection_predictor.cc index 4aec350169..47d7979a7a 100644 --- a/src/dials/algorithms/spot_prediction/boost_python/reflection_predictor.cc +++ b/src/dials/algorithms/spot_prediction/boost_python/reflection_predictor.cc @@ -184,7 +184,10 @@ namespace dials { namespace algorithms { namespace boost_python { .def("__call__", predict_observed_with_panel) .def("__call__", predict_observed_with_panel_list) .def("for_reflection_table", &Predictor::for_reflection_table) - .def("for_reflection_table", &Predictor::for_reflection_table_with_individual_ub); + .def("for_reflection_table", &Predictor::for_reflection_table_with_individual_ub) + .def("all_reflections_for_asu", + &Predictor::all_reflections_for_asu, + (arg("goniometer"), arg("phi"))); } void export_reflection_predictor() { export_scan_static_reflection_predictor(); diff --git a/src/dials/algorithms/spot_prediction/reflection_predictor.h b/src/dials/algorithms/spot_prediction/reflection_predictor.h index b9345e9d84..f50ca0c7c0 100644 --- a/src/dials/algorithms/spot_prediction/reflection_predictor.h +++ b/src/dials/algorithms/spot_prediction/reflection_predictor.h @@ -28,6 +28,7 @@ #include #include #include +#include namespace dials { namespace algorithms { @@ -1356,6 +1357,70 @@ namespace dials { namespace algorithms { return af::reflection_table(); } + af::reflection_table all_reflections_for_asu(Goniometer goniometer, double phi) { + mat3 fixed_rotation = goniometer.get_fixed_rotation(); + mat3 setting_rotation = goniometer.get_setting_rotation(); + vec3 rotation_axis = goniometer.get_rotation_axis(); + mat3 rotation = + scitbx::math::r3_rotation::axis_and_angle_as_matrix(rotation_axis, phi); + vec3 unit_s0 = beam_.get_unit_s0(); + vec2 wavelength_range = beam_.get_wavelength_range(); + + cctbx::miller::index_generator indices = + cctbx::miller::index_generator(unit_cell_, space_group_type_, false, dmin_); + + af::shared indices_arr = indices.to_array(); + + af::reflection_table table; + af::shared wavelength_column; + table["wavelength_cal"] = wavelength_column; + af::shared > s0_column; + table["s0_cal"] = s0_column; + tof_prediction_data predictions(table); + + for (std::size_t i = 0; i < indices_arr.size(); ++i) { + miller_index h = indices_arr[i]; + + vec3 q = setting_rotation * rotation * fixed_rotation * ub_ * h; + + // Calculate the wavelength required to meet the diffraction condition + double wavelength = -2 * ((unit_s0 * q) / (q * q)); + if (wavelength < wavelength_range[0] || wavelength > wavelength_range[1]) { + continue; + } + vec3 s0 = unit_s0 / wavelength; + DIALS_ASSERT(s0.length() > 0); + + // Calculate the Ray (default zero angle and 'entering' as false) + vec3 s1 = s0 + q; + + int panel = detector_.get_panel_intersection(s1); + if (panel == -1) { + continue; + } + + Detector::coord_type coord; + coord.first = panel; + coord.second = detector_[panel].get_ray_intersection(s1); + vec2 mm = coord.second; + vec2 px = detector_[panel].millimeter_to_pixel(mm); + + // Add the reflections to the table + predictions.hkl.push_back(h); + predictions.enter.push_back(false); + predictions.s1.push_back(s1); + predictions.xyz_mm.push_back(vec3(mm[0], mm[1], 0.0)); + predictions.xyz_px.push_back(vec3(px[0], px[1], 0.0)); + predictions.panel.push_back(panel); + predictions.flags.push_back(af::Predicted); + predictions.wavelength_cal.push_back(wavelength); + predictions.s0_cal.push_back(s0); + } + + // Return the reflection table + return table; + } + /** * Predict reflections for UB. Also filters based on ewald sphere proximity. * @param ub The UB matrix @@ -1514,14 +1579,6 @@ namespace dials { namespace algorithms { append_for_ray(p, h, ray, panel, wavelength, s0); } - /** - * Predict for the given Miller index, ray, panel number and delta psi - * @param p The reflection data - * @param h The miller index - * @param ray The ray data - * @param panel The panel number - * @param delpsi The calculated minimum rotation to Ewald sphere - */ void append_for_ray(tof_prediction_data &p, const miller_index &h, const Ray &ray, diff --git a/src/dials/command_line/dials_import.py b/src/dials/command_line/dials_import.py index d3023c1c2b..44ea9e3973 100644 --- a/src/dials/command_line/dials_import.py +++ b/src/dials/command_line/dials_import.py @@ -209,7 +209,6 @@ def _extract_or_read_imagesets(params): # Check we have some filenames if len(experiments) == 0: - # FIXME Should probably make this smarter since it requires editing here # and in dials.import phil scope try: @@ -223,7 +222,6 @@ def _extract_or_read_imagesets(params): # Check if a template has been set and print help if not, otherwise try to # import the images based on the template input if len(params.input.template) > 0: - experiments = ExperimentListFactory.from_templates( params.input.template, image_range=params.geometry.scan.image_range, @@ -530,7 +528,6 @@ def __call__(self, imageset_list): # Loop through imagesets for imageset in imageset_list: - # Set the external lookups imageset = self.update_lookup(imageset, lookup) @@ -794,7 +791,6 @@ def assert_single_sequence(experiments, params): ] if len(sequences) > 1: - # Print some info about multiple sequences diagnose_multiple_sequences(sequences, params) @@ -914,7 +910,6 @@ def do_import( # Print out info for all experiments for experiment in experiments: - # Print some experiment info - override the output of image range # if appropriate image_range = params.geometry.scan.image_range