From 566f7c3fa9514e8d8cc64287242a857d7f634f48 Mon Sep 17 00:00:00 2001 From: davidmcdonagh Date: Fri, 15 Dec 2023 00:02:24 +0000 Subject: [PATCH] Fixed frame not being in indexed reflections. Fixed dmin being calculated for all experiments. Fixed predicted experiments not having ids set correctly. Fixed line profile values having prf columns. --- .../integration/fit/tof_line_profile.py | 10 +- .../algorithms/spot_prediction/__init__.py | 6 +- .../boost_python/reflection_predictor.cc | 1 + .../spot_prediction/reflection_predictor.h | 9 +- src/dials/array_family/flex_ext.py | 1 + .../command_line/simple_tof_integrate.py | 99 +++++++++++-------- 6 files changed, 79 insertions(+), 47 deletions(-) diff --git a/src/dials/algorithms/integration/fit/tof_line_profile.py b/src/dials/algorithms/integration/fit/tof_line_profile.py index 9a6d026ffd..8c3c340475 100644 --- a/src/dials/algorithms/integration/fit/tof_line_profile.py +++ b/src/dials/algorithms/integration/fit/tof_line_profile.py @@ -219,10 +219,14 @@ def compute_line_profile_intensity(reflections): fit_variance = abs(fit_intensity) + abs(background_sum) * (1.0 + m_n) fit_variances[i] = fit_variance - reflections["line_profile_intensity"] = fit_intensities - reflections["line_profile_variance"] = fit_variances + reflections["intensity.prf.value"] = fit_intensities + reflections["intensity.prf.variance"] = fit_variances reflections.set_flags( - reflections["line_profile_intensity"] < 0, + 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/spot_prediction/__init__.py b/src/dials/algorithms/spot_prediction/__init__.py index 4546e4dfc1..2a511b2fd4 100644 --- a/src/dials/algorithms/spot_prediction/__init__.py +++ b/src/dials/algorithms/spot_prediction/__init__.py @@ -116,6 +116,7 @@ def __init__(self, experiment, dmin): experiment.beam, experiment.detector, experiment.goniometer, + experiment.sequence, experiment.crystal.get_A(), experiment.crystal.get_unit_cell(), experiment.crystal.get_space_group().type(), @@ -147,9 +148,8 @@ def post_prediction(self, reflections): for idx in range(len(expt_reflections)): wavelength = expt_reflections[idx]["wavelength_cal"] - expt_tof_cal[idx] = expt.beam.get_tof_from_wavelength( - wavelength, expt_L1[idx] - ) + tof = expt.beam.get_tof_from_wavelength(wavelength, expt_L1[idx]) + expt_tof_cal[idx] = tof tof_cal.set_selected(sel, expt_tof_cal) L1.set_selected(sel, expt_L1) 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 aef90c7727..3a6880f63b 100644 --- a/src/dials/algorithms/spot_prediction/boost_python/reflection_predictor.cc +++ b/src/dials/algorithms/spot_prediction/boost_python/reflection_predictor.cc @@ -175,6 +175,7 @@ namespace dials { namespace algorithms { namespace boost_python { .def(init, const cctbx::uctbx::unit_cell&, const cctbx::sgtbx::space_group_type&, diff --git a/src/dials/algorithms/spot_prediction/reflection_predictor.h b/src/dials/algorithms/spot_prediction/reflection_predictor.h index b0b3beb810..309bdd8f77 100644 --- a/src/dials/algorithms/spot_prediction/reflection_predictor.h +++ b/src/dials/algorithms/spot_prediction/reflection_predictor.h @@ -42,6 +42,7 @@ namespace dials { namespace algorithms { using dxtbx::model::plane_ray_intersection; using dxtbx::model::PolyBeam; using dxtbx::model::Scan; + using dxtbx::model::TOFSequence; using scitbx::constants::pi; using scitbx::constants::pi_180; using scitbx::constants::two_pi; @@ -1337,6 +1338,7 @@ namespace dials { namespace algorithms { TOFReflectionPredictor(const PolyBeam &beam, const Detector &detector, const Goniometer &goniometer, + const TOFSequence &sequence, mat3 ub, const cctbx::uctbx::unit_cell &unit_cell, const cctbx::sgtbx::space_group_type &space_group_type, @@ -1344,6 +1346,7 @@ namespace dials { namespace algorithms { : beam_(beam), detector_(detector), goniometer_(goniometer), + sequence_(sequence), ub_(ub), unit_cell_(unit_cell), space_group_type_(space_group_type), @@ -1595,13 +1598,16 @@ namespace dials { namespace algorithms { std::size_t panel = impact.first; vec2 mm = impact.second; vec2 px = detector_[panel].millimeter_to_pixel(mm); + double L1 = ray.s1.length() * std::pow(10, -3); + double tof = beam_.get_tof_from_wavelength(wavelength, L1); + double frame = sequence_.get_frame_from_tof(tof); // Add the reflections to the table p.hkl.push_back(h); p.enter.push_back(ray.entering); p.s1.push_back(ray.s1); p.xyz_mm.push_back(vec3(mm[0], mm[1], 0.0)); - p.xyz_px.push_back(vec3(px[0], px[1], 0.0)); + p.xyz_px.push_back(vec3(px[0], px[1], frame)); p.panel.push_back(panel); p.flags.push_back(af::Predicted); p.wavelength_cal.push_back(wavelength); @@ -1631,6 +1637,7 @@ namespace dials { namespace algorithms { PolyBeam beam_; Detector detector_; Goniometer goniometer_; + TOFSequence sequence_; mat3 ub_; cctbx::uctbx::unit_cell unit_cell_; cctbx::sgtbx::space_group_type space_group_type_; diff --git a/src/dials/array_family/flex_ext.py b/src/dials/array_family/flex_ext.py index 286682f47b..6d098264ed 100644 --- a/src/dials/array_family/flex_ext.py +++ b/src/dials/array_family/flex_ext.py @@ -689,6 +689,7 @@ def match( x *= scale[0] y *= scale[1] z *= scale[2] + xyz = cctbx.array_family.flex.vec3_double(x, y, z) x, y, z = ref.parts() diff --git a/src/dials/command_line/simple_tof_integrate.py b/src/dials/command_line/simple_tof_integrate.py index eb80eb0d37..c24bea0b78 100644 --- a/src/dials/command_line/simple_tof_integrate.py +++ b/src/dials/command_line/simple_tof_integrate.py @@ -296,32 +296,34 @@ def print_data(reflections, panel): def output_reflections_as_hkl(reflections, filename): - def get_corrected_intensity_and_sigma(reflections, idx): + def get_corrected_intensity_and_variance(reflections, idx): intensity = reflections["intensity.sum.value"][idx] variance = reflections["intensity.sum.variance"][idx] - return intensity, np.sqrt(variance) + return intensity, variance - def get_line_profile_intensity_and_sigma(reflections, idx): - intensity = reflections["line_profile_intensity"][idx] - variance = reflections["line_profile_variance"][idx] - return intensity, np.sqrt(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): + def valid_intensity(intensity, variance): from math import isinf, isnan if isnan(intensity) or isinf(intensity): return False - return intensity > 0 + 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, sigma = get_corrected_intensity_and_sigma(reflections, i) - if not valid_intensity(intensity): + intensity, variance = get_corrected_intensity_and_variance(reflections, i) + if not valid_intensity(intensity, variance): continue intensity = round(intensity, 2) - sigma = round(sigma, 2) + sigma = round(np.sqrt(variance), 2) wavelength = round(reflections["wavelength_cal"][i], 4) g.write( "" @@ -342,15 +344,17 @@ def valid_intensity(intensity): ) ) - with open("lp_" + filename, "w") as g: + with open("prf_" + filename, "w") as g: for i in range(len(reflections)): h, k, l = reflections["miller_index"][i] batch_number = 1 - intensity, sigma = get_line_profile_intensity_and_sigma(reflections, i) - if not valid_intensity(intensity): + intensity, variance = get_line_profile_intensity_and_variance( + reflections, i + ) + if not valid_intensity(intensity, variance): continue intensity = round(intensity, 2) - sigma = round(sigma, 2) + sigma = round(np.sqrt(variance), 2) wavelength = round(reflections["wavelength_cal"][i], 4) g.write( "" @@ -526,9 +530,6 @@ def run(): ) reflections = reflections[0] - reflections["id"] = cctbx.array_family.flex.int(len(reflections), 0) - reflections["imageset_id"] = cctbx.array_family.flex.int(len(reflections), 0) - integrated_reflections = run_simple_integrate(params, experiments, reflections) integrated_reflections.as_msgpack_file(params.output.reflections) experiments.as_file(params.output.experiments) @@ -540,8 +541,6 @@ def run_simple_integrate(params, experiments, reflections): nproc = 11 pool = multiprocessing.Pool(nproc) - experiment = experiments[0] - reflections, _ = process_reference(reflections) """ @@ -552,16 +551,34 @@ def run_simple_integrate(params, experiments, reflections): range(len(reflections["wavelength"])), key=reflections["wavelength"].__getitem__ ) min_s0 = reflections["s0"][min_s0_idx] - dmin = experiment.detector.get_max_resolution(min_s0) - predicted_reflections = flex.reflection_table.from_predictions( - experiment, padding=1.0, dmin=dmin - ) - predicted_reflections["id"] = cctbx.array_family.flex.int( - len(predicted_reflections), 0 - ) - predicted_reflections["imageset_id"] = cctbx.array_family.flex.int( - len(predicted_reflections), 0 - ) + 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) + # Updates flags to set which reflections to use in generating reference profiles matched, reflections, unmatched = predicted_reflections.match_with_reference( reflections @@ -584,9 +601,10 @@ def run_simple_integrate(params, experiments, reflections): sigma_m = 3 sigma_b = 0.01 # The Gaussian model given in 2.3 of Kabsch 2010 - experiment.profile = GaussianRSProfileModel( - params=params, n_sigma=3, sigma_b=sigma_b, sigma_m=sigma_m - ) + for idx, experiment in enumerate(experiments): + experiments[idx].profile = GaussianRSProfileModel( + params=params, n_sigma=3, sigma_b=sigma_b, sigma_m=sigma_m + ) """ Compute properties for predicted reflections using profile model, @@ -599,7 +617,7 @@ def run_simple_integrate(params, experiments, reflections): predicted_reflections.compute_bbox(experiments) x1, x2, y1, y2, t1, t2 = predicted_reflections["bbox"].parts() predicted_reflections = predicted_reflections.select( - t2 < experiment.sequence.get_image_range()[1] + t2 < experiments[0].sequence.get_image_range()[1] ) predicted_reflections.compute_d(experiments) predicted_reflections.compute_partiality(experiments) @@ -615,16 +633,17 @@ def run_simple_integrate(params, experiments, reflections): # Get actual shoebox values and the reflections for each image shoebox_processor = ShoeboxProcessor( predicted_reflections, - len(experiment.detector), + len(experiments[0].detector), 0, - len(experiment.imageset), + sum([len(experiment.imageset) for experiment in experiments]), False, ) - for i in range(len(experiment.imageset)): - image = experiment.imageset.get_corrected_data(i) - mask = experiment.imageset.get_mask(i) - shoebox_processor.next_data_only(make_image(image, mask)) + for experiment in experiments: + for i in range(len(experiment.imageset)): + image = experiment.imageset.get_corrected_data(i) + mask = experiment.imageset.get_mask(i) + shoebox_processor.next_data_only(make_image(image, mask)) predicted_reflections.is_overloaded(experiments) predicted_reflections.compute_mask(experiments)