Skip to content

Commit

Permalink
Added method to predict spots for asymmetric unit.
Browse files Browse the repository at this point in the history
  • Loading branch information
toastisme committed Nov 15, 2023
1 parent 7b973b9 commit 273b300
Show file tree
Hide file tree
Showing 4 changed files with 70 additions and 15 deletions.
2 changes: 1 addition & 1 deletion src/dials/algorithms/indexing/model_evaluation.py
Original file line number Diff line number Diff line change
Expand Up @@ -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),
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -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();
Expand Down
73 changes: 65 additions & 8 deletions src/dials/algorithms/spot_prediction/reflection_predictor.h
Original file line number Diff line number Diff line change
Expand Up @@ -28,6 +28,7 @@
#include <dials/algorithms/spot_prediction/stills_ray_predictor.h>
#include <dials/algorithms/spot_prediction/ray_intersection.h>
#include <iostream>
#include <cctbx/miller/index_generator.h>

namespace dials { namespace algorithms {

Expand Down Expand Up @@ -1356,6 +1357,70 @@ namespace dials { namespace algorithms {
return af::reflection_table();
}

af::reflection_table all_reflections_for_asu(Goniometer goniometer, double phi) {
mat3<double> fixed_rotation = goniometer.get_fixed_rotation();
mat3<double> setting_rotation = goniometer.get_setting_rotation();
vec3<double> rotation_axis = goniometer.get_rotation_axis();
mat3<double> rotation =
scitbx::math::r3_rotation::axis_and_angle_as_matrix(rotation_axis, phi);
vec3<double> unit_s0 = beam_.get_unit_s0();
vec2<double> wavelength_range = beam_.get_wavelength_range();

cctbx::miller::index_generator indices =
cctbx::miller::index_generator(unit_cell_, space_group_type_, false, dmin_);

af::shared<miller_index> indices_arr = indices.to_array();

af::reflection_table table;
af::shared<double> wavelength_column;
table["wavelength_cal"] = wavelength_column;
af::shared<vec3<double> > 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<double> 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<double> s0 = unit_s0 / wavelength;
DIALS_ASSERT(s0.length() > 0);

// Calculate the Ray (default zero angle and 'entering' as false)
vec3<double> 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<double> mm = coord.second;
vec2<double> 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<double>(mm[0], mm[1], 0.0));
predictions.xyz_px.push_back(vec3<double>(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
Expand Down Expand Up @@ -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,
Expand Down
5 changes: 0 additions & 5 deletions src/dials/command_line/dials_import.py
Original file line number Diff line number Diff line change
Expand Up @@ -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:
Expand All @@ -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,
Expand Down Expand Up @@ -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)

Expand Down Expand Up @@ -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)

Expand Down Expand Up @@ -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
Expand Down

0 comments on commit 273b300

Please sign in to comment.