Fixed frame not being in indexed reflections. Fixed dmin being calcul…

…ated for all experiments. Fixed predicted experiments not having ids set correctly. Fixed line profile values having prf columns.

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

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)
 

src/dials/algorithms/spot_prediction/boost_python/reflection_predictor.cc

@@ -175,6 +175,7 @@ namespace dials { namespace algorithms { namespace boost_python {
  .def(init<const PolyBeam&,
  const Detector&,
  const Goniometer&,
+ const TOFSequence&,
  mat3<double>,
  const cctbx::uctbx::unit_cell&,
  const cctbx::sgtbx::space_group_type&,

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,13 +1338,15 @@ namespace dials { namespace algorithms {
  TOFReflectionPredictor(const PolyBeam &beam,
  const Detector &detector,
  const Goniometer &goniometer,
+ const TOFSequence &sequence,
  mat3<double> 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),
+ 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<double> mm = impact.second;
  vec2<double> 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<double>(mm[0], mm[1], 0.0));
- p.xyz_px.push_back(vec3<double>(px[0], px[1], 0.0));
+ p.xyz_px.push_back(vec3<double>(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<double> ub_;
  cctbx::uctbx::unit_cell unit_cell_;
  cctbx::sgtbx::space_group_type space_group_type_;

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()

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)