From 4418bb7112c95bc7e955364f82867302277ef902 Mon Sep 17 00:00:00 2001 From: davidmcdonagh Date: Thu, 14 Dec 2023 10:19:17 +0000 Subject: [PATCH] Added goniometer rotations to ToF workflow. --- src/dials/algorithms/integration/processor.py | 2 +- .../refinement/parameterisation/configure.py | 4 +- .../prediction/managed_predictors.py | 2 +- src/dials/algorithms/refinement/refiner.py | 2 +- .../refinement/reflection_manager.py | 6 +- src/dials/algorithms/refinement/target.py | 2 +- .../algorithms/spot_prediction/__init__.py | 1 + .../boost_python/ray_predictor.cc | 6 +- .../boost_python/reflection_predictor.cc | 1 + .../spot_prediction/ray_predictor.h | 91 ++++++++++--------- .../spot_prediction/reflection_predictor.h | 7 +- src/dials/command_line/integrate.py | 2 +- src/dials/command_line/symmetry.py | 2 +- 13 files changed, 71 insertions(+), 57 deletions(-) diff --git a/src/dials/algorithms/integration/processor.py b/src/dials/algorithms/integration/processor.py index cf091e6500..57a0fb1fce 100644 --- a/src/dials/algorithms/integration/processor.py +++ b/src/dials/algorithms/integration/processor.py @@ -930,7 +930,7 @@ def summary(self): Get a summary of the processing """ # Compute the task table - if self.experiments.all_stills() or self.experiments.is_single_tof_experiment(): + if self.experiments.all_stills() or self.experiments.all_tof_experiments(): rows = [["#", "Group", "Frame From", "Frame To", "# Reflections"]] for i in range(len(self)): job = self.manager.job(i) diff --git a/src/dials/algorithms/refinement/parameterisation/configure.py b/src/dials/algorithms/refinement/parameterisation/configure.py index 0f84a34f43..1acc6453f6 100644 --- a/src/dials/algorithms/refinement/parameterisation/configure.py +++ b/src/dials/algorithms/refinement/parameterisation/configure.py @@ -395,7 +395,7 @@ def _set_n_intervals(smoother_params, analysis, scan, exp_ids): def _parameterise_beams(options, experiments, analysis): beam_params = [] - if experiments.is_single_tof_experiment(): + if experiments.all_tof_experiments(): return beam_params sv_beam = options.scan_varying and not options.beam.force_static @@ -810,7 +810,7 @@ def build_prediction_parameterisation( # Build the prediction equation parameterisation if do_stills: # doing stills - if experiments.is_single_tof_experiment(): + if experiments.all_tof_experiments(): PredParam = TOFPredictionParameterisation elif options.sparse: if options.spherical_relp_model: diff --git a/src/dials/algorithms/refinement/prediction/managed_predictors.py b/src/dials/algorithms/refinement/prediction/managed_predictors.py index 18a6bb6306..ed36ff3cb1 100644 --- a/src/dials/algorithms/refinement/prediction/managed_predictors.py +++ b/src/dials/algorithms/refinement/prediction/managed_predictors.py @@ -235,7 +235,7 @@ class ExperimentsPredictorFactory: @staticmethod def from_experiments(experiments, force_stills=False, spherical_relp=False): - if experiments.is_single_tof_experiment(): + if experiments.all_tof_experiments(): return TOFExperimentsPredictor(experiments) # Determine whether or not to use a stills predictor diff --git a/src/dials/algorithms/refinement/refiner.py b/src/dials/algorithms/refinement/refiner.py index afe4e633d8..f13878996d 100644 --- a/src/dials/algorithms/refinement/refiner.py +++ b/src/dials/algorithms/refinement/refiner.py @@ -373,7 +373,7 @@ def _build_components(cls, params, reflections, experiments): obs["x_resid"] = x_calc - x_obs obs["y_resid"] = y_calc - y_obs - if experiments.is_single_tof_experiment(): + if experiments.all_tof_experiments(): obs["wavelength_resid"] = obs["wavelength_cal"] - obs["wavelength"] obs["wavelength_resid2"] = (obs["wavelength_cal"] - obs["wavelength"]) ** 2 diff --git a/src/dials/algorithms/refinement/reflection_manager.py b/src/dials/algorithms/refinement/reflection_manager.py index 9db98edd24..285ad2f7d4 100644 --- a/src/dials/algorithms/refinement/reflection_manager.py +++ b/src/dials/algorithms/refinement/reflection_manager.py @@ -229,7 +229,7 @@ def from_parameters_reflections_experiments( flex.set_random_seed(params.random_seed) logger.debug("Random seed set to %d", params.random_seed) - if experiments.is_single_tof_experiment(): + if experiments.all_tof_experiments(): refman = TOFReflectionManager # check whether we deal with stills or scans @@ -269,7 +269,7 @@ def from_parameters_reflections_experiments( if params.outlier.algorithm in ("null", None): outlier_detector = None else: - if experiments.is_single_tof_experiment(): + if experiments.all_tof_experiments(): colnames = ["x_resid", "y_resid", "wavelength_resid"] elif do_stills: colnames = ["x_resid", "y_resid"] @@ -362,7 +362,7 @@ def __init__( self._axes = [ matrix.col(g.get_rotation_axis()) if g else None for g in goniometers ] - if experiments.is_single_tof_experiment(): + if experiments.all_tof_experiments(): self._s0vecs = [None] else: self._s0vecs = [matrix.col(e.beam.get_s0()) for e in self._experiments] diff --git a/src/dials/algorithms/refinement/target.py b/src/dials/algorithms/refinement/target.py index 64af1c8ff0..01304592d9 100644 --- a/src/dials/algorithms/refinement/target.py +++ b/src/dials/algorithms/refinement/target.py @@ -80,7 +80,7 @@ def from_parameters_and_experiments( # Determine whether the target is in X, Y, Phi space or just X, Y to choose # the right Target to instantiate - if experiments.is_single_tof_experiment(): + if experiments.all_tof_experiments(): targ = TOFLeastSquaresPositionalResidualWithRmsdCutoff elif do_stills: if do_sparse: diff --git a/src/dials/algorithms/spot_prediction/__init__.py b/src/dials/algorithms/spot_prediction/__init__.py index 693aabeea3..4546e4dfc1 100644 --- a/src/dials/algorithms/spot_prediction/__init__.py +++ b/src/dials/algorithms/spot_prediction/__init__.py @@ -115,6 +115,7 @@ def __init__(self, experiment, dmin): self.predictor = TOFReflectionPredictor( experiment.beam, experiment.detector, + experiment.goniometer, experiment.crystal.get_A(), experiment.crystal.get_unit_cell(), experiment.crystal.get_space_group().type(), diff --git a/src/dials/algorithms/spot_prediction/boost_python/ray_predictor.cc b/src/dials/algorithms/spot_prediction/boost_python/ray_predictor.cc index e3a1d64b1d..fd4535c2f6 100644 --- a/src/dials/algorithms/spot_prediction/boost_python/ray_predictor.cc +++ b/src/dials/algorithms/spot_prediction/boost_python/ray_predictor.cc @@ -62,9 +62,9 @@ namespace dials { namespace algorithms { namespace boost_python { void export_tof_ray_predictor() { // Create and return the wrapper for the spot predictor object class_("TOFRayPredictor", no_init) - .def(init >((arg("unit_s0")))) - .def( - "__call__", &TOFRayPredictor::operator(), (arg("miller_index"), arg("UB"))); + .def(init, mat3, mat3 >( + (arg("unit_s0"), arg("fixed_rotation"), arg("setting_rotation")))) + .def("__call__", &TOFRayPredictor::operator(), (arg("miller_index"), arg("UB"))); } }}} // namespace dials::algorithms::boost_python 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 47d7979a7a..aef90c7727 100644 --- a/src/dials/algorithms/spot_prediction/boost_python/reflection_predictor.cc +++ b/src/dials/algorithms/spot_prediction/boost_python/reflection_predictor.cc @@ -174,6 +174,7 @@ namespace dials { namespace algorithms { namespace boost_python { class_("TOFReflectionPredictor", no_init) .def(init, const cctbx::uctbx::unit_cell&, const cctbx::sgtbx::space_group_type&, diff --git a/src/dials/algorithms/spot_prediction/ray_predictor.h b/src/dials/algorithms/spot_prediction/ray_predictor.h index dd0f4c4144..1f7de9bd26 100644 --- a/src/dials/algorithms/spot_prediction/ray_predictor.h +++ b/src/dials/algorithms/spot_prediction/ray_predictor.h @@ -88,7 +88,7 @@ namespace dials { namespace algorithms { vec2 phi; try { phi = calculate_rotation_angles_(pstar0); - } catch (error const&) { + } catch (error const &) { return rays; } @@ -124,7 +124,7 @@ namespace dials { namespace algorithms { vec2 phi; try { phi = calculate_rotation_angles_(pstar0); - } catch (error const&) { + } catch (error const &) { return rays; } @@ -251,52 +251,59 @@ namespace dials { namespace algorithms { }; /** - * Class to predict s1 rays for ToF data + * Class to predict s1 rays for ToF data */ - class TOFRayPredictor{ - public: - typedef cctbx::miller::index<> miller_index; + class TOFRayPredictor { + public: + typedef cctbx::miller::index<> miller_index; - TOFRayPredictor(const vec3 unit_s0) : unit_s0_(unit_s0){ - DIALS_ASSERT(unit_s0_.length() > 0.0); - } + TOFRayPredictor(const vec3 unit_s0, + mat3 fixed_rotation, + mat3 setting_rotation) + : unit_s0_(unit_s0), + fixed_rotation_(fixed_rotation), + setting_rotation_(setting_rotation) - /** - * For a given miller index and UB matrix, calculates the predicted s1 ray. - * The TOFRayPredictor wavelength and s0 variables are updated during the calculation, - * so that they can be monitored for convergence. - * @param h The miller index - * @param ub The UB matrix - * @returns Ray - */ - Ray operator()(const miller_index &h, - const mat3 &ub){ - - // Calculate the reciprocal lattice vector - vec3 q = ub * h; - - // Calculate the wavelength required to meet the diffraction condition - wavelength_ = -2*((unit_s0_ * q)/(q *q)); - s0_ = unit_s0_ / wavelength_; - DIALS_ASSERT(s0_.length() > 0); - - // Calculate the Ray (default zero angle and 'entering' as false) - vec3 s1 = s0_ + q; - return Ray(s1, 0.0, false); - } + { + DIALS_ASSERT(unit_s0_.length() > 0.0); + } - double get_wavelength() const{ - return wavelength_; - } + /** + * For a given miller index and UB matrix, calculates the predicted s1 ray. + * The TOFRayPredictor wavelength and s0 variables are updated during the + * calculation, so that they can be monitored for convergence. + * @param h The miller index + * @param ub The UB matrix + * @returns Ray + */ + Ray operator()(const miller_index &h, const mat3 &ub) { + // Calculate the reciprocal lattice vector + vec3 q = setting_rotation_ * fixed_rotation_ * ub * h; + + // Calculate the wavelength required to meet the diffraction condition + wavelength_ = -2 * ((unit_s0_ * q) / (q * q)); + s0_ = unit_s0_ / wavelength_; + DIALS_ASSERT(s0_.length() > 0); + + // Calculate the Ray (default zero angle and 'entering' as false) + vec3 s1 = s0_ + q; + return Ray(s1, 0.0, false); + } - vec3 get_s0() const{ - return s0_; - } + double get_wavelength() const { + return wavelength_; + } + + vec3 get_s0() const { + return s0_; + } - private: - const vec3 unit_s0_; - double wavelength_; - vec3 s0_; + private: + const vec3 unit_s0_; + double wavelength_; + vec3 s0_; + mat3 fixed_rotation_; + mat3 setting_rotation_; }; }} // namespace dials::algorithms diff --git a/src/dials/algorithms/spot_prediction/reflection_predictor.h b/src/dials/algorithms/spot_prediction/reflection_predictor.h index f50ca0c7c0..b0b3beb810 100644 --- a/src/dials/algorithms/spot_prediction/reflection_predictor.h +++ b/src/dials/algorithms/spot_prediction/reflection_predictor.h @@ -1336,17 +1336,21 @@ namespace dials { namespace algorithms { */ TOFReflectionPredictor(const PolyBeam &beam, const Detector &detector, + const Goniometer &goniometer, mat3 ub, const cctbx::uctbx::unit_cell &unit_cell, const cctbx::sgtbx::space_group_type &space_group_type, const double &dmin) : beam_(beam), detector_(detector), + goniometer_(goniometer), ub_(ub), unit_cell_(unit_cell), space_group_type_(space_group_type), dmin_(dmin), - predict_ray_(beam.get_unit_s0()) {} + predict_ray_(beam.get_unit_s0(), + goniometer.get_fixed_rotation(), + goniometer.get_setting_rotation()) {} /** * Predict all reflection. @@ -1626,6 +1630,7 @@ namespace dials { namespace algorithms { protected: PolyBeam beam_; Detector detector_; + Goniometer goniometer_; mat3 ub_; cctbx::uctbx::unit_cell unit_cell_; cctbx::sgtbx::space_group_type space_group_type_; diff --git a/src/dials/command_line/integrate.py b/src/dials/command_line/integrate.py index c82a930449..8e875df38c 100644 --- a/src/dials/command_line/integrate.py +++ b/src/dials/command_line/integrate.py @@ -486,7 +486,7 @@ def run_integration(params, experiments, reference=None): logger.info("\n".join(("", "=" * 80, ""))) logger.info(heading("Predicting reflections")) - if experiments.is_single_tof_experiment() and not params.prediction.d_min: + if experiments.all_tof_experiments() and not params.prediction.d_min: min_s0_idx = min( range(len(reference["wavelength"])), key=reference["wavelength"].__getitem__ diff --git a/src/dials/command_line/symmetry.py b/src/dials/command_line/symmetry.py index 6c9647f2f5..7edbf13ba2 100644 --- a/src/dials/command_line/symmetry.py +++ b/src/dials/command_line/symmetry.py @@ -269,7 +269,7 @@ def eliminate_sys_absent(experiments, reflections): def get_subset_for_symmetry(experiments, reflection_tables, exclude_images=None): """Select an image range for symmetry analysis, or just select the first 360 degrees of data.""" - if experiments.is_single_tof_experiment(): + if experiments.all_tof_experiments(): return reflection_tables refls_for_sym = [] if exclude_images: