Outline for on the fly scaling corrections.

src/dials/algorithms/integration/boost_python/integration_integrator_ext.cc

@@ -14,6 +14,7 @@
 #include <dials/algorithms/integration/integrator.h>
 #include <dials/algorithms/integration/manager.h>
 #include <dxtbx/array_family/flex_table_suite.h>
+#include <dials/algorithms/scaling/tof_scaling_corrections.h>
 
 using namespace boost::python;
 
@@ -349,6 +350,36 @@ namespace dials { namespace algorithms { namespace boost_python {
  def("max_memory_needed",
  &max_memory_needed,
  (arg("reflection table"), arg("frame0"), arg("frame1"), arg("flatten")));
+
+ def("tof_extract_shoeboxes_to_reflection_table",
+ &tof_extract_shoeboxes_to_reflection_table,
+ (arg("reflection_table"),
+ arg("experiment"),
+ arg("incident_data"),
+ arg("empty_data"),
+ arg("sample_radius"),
+ arg("sample_scattering_x_section"),
+ arg("sample_absorption_x_section"),
+ arg("sample_number_density"),
+ arg("incident_radius"),
+ arg("incident_scattering_x_section"),
+ arg("incident_absorption_x_section"),
+ arg("incident_number_density")));
+
+ def("tof_extract_shoeboxes_to_reflection_table2",
+ &tof_extract_shoeboxes_to_reflection_table2,
+ (arg("reflection_table"),
+ arg("experiment"),
+ arg("incident_data"),
+ arg("empty_data"),
+ arg("sample_radius"),
+ arg("sample_scattering_x_section"),
+ arg("sample_absorption_x_section"),
+ arg("sample_number_density"),
+ arg("incident_radius"),
+ arg("incident_scattering_x_section"),
+ arg("incident_absorption_x_section"),
+ arg("incident_number_density")));
  }
 
 }}} // namespace dials::algorithms::boost_python

src/dials/algorithms/profile_model/gaussian_rs/__init__.py

@@ -6,6 +6,7 @@
  BBoxCalculator2D,
  BBoxCalculator3D,
  BBoxCalculatorIface,
+ BBoxCalculatorTOF,
  BBoxMultiCalculator,
  CoordinateSystem,
  CoordinateSystem2d,
@@ -27,6 +28,7 @@
  "BBoxCalculator",
  "BBoxCalculator2D",
  "BBoxCalculator3D",
+ "BBoxCalculatorTOF",
  "BBoxCalculatorIface",
  "BBoxMultiCalculator",
  "CoordinateSystem",
@@ -52,7 +54,9 @@
 
 def BBoxCalculator(crystal, beam, detector, goniometer, scan, delta_b, delta_m):
  """Return the relevant bbox calculator."""
- if goniometer is None or scan is None or scan.is_still():
+ if scan.has_property("time_of_flight"):
+ algorithm = BBoxCalculatorTOF(beam, detector, scan, delta_b, delta_m)
+ elif goniometer is None or scan is None or scan.is_still():
  algorithm = BBoxCalculator2D(beam, detector, delta_b, delta_m)
  else:
  algorithm = BBoxCalculator3D(beam, detector, goniometer, scan, delta_b, delta_m)

src/dials/algorithms/profile_model/gaussian_rs/bbox_calculator.h

@@ -33,6 +33,7 @@ namespace dials {
  using dxtbx::model::BeamBase;
  using dxtbx::model::Detector;
  using dxtbx::model::Goniometer;
+ using dxtbx::model::PolychromaticBeam;
  using dxtbx::model::Scan;
  using scitbx::vec2;
  using scitbx::vec3;
@@ -349,6 +350,139 @@ namespace dials {
  double delta_divergence_;
  };
 
+ class BBoxCalculatorTOF {
+ public:
+ /**
+ * Initialise the bounding box calculation.
+ * @param beam The beam parameters
+ * @param detector The detector parameters
+ * @param delta_divergence The xds delta_divergence parameter
+ * @param delta_mosaicity The xds delta_mosaicity parameter
+ */
+ BBoxCalculatorTOF(const PolychromaticBeam &beam,
+ const Detector &detector,
+ const Scan &scan,
+ double delta_divergence,
+ double delta_mosaicity)
+ : detector_(detector),
+ scan_(scan),
+ beam_(beam),
+ delta_divergence_(delta_divergence),
+ delta_mosaicity_(delta_mosaicity) {
+ DIALS_ASSERT(delta_divergence > 0.0);
+ DIALS_ASSERT(delta_mosaicity >= 0.0);
+ }
+
+ /**
+ * Calculate the bbox on the detector image volume for the reflection.
+ *
+ * The roi is calculated using the parameters delta_divergence and
+ * delta_mosaicity. The reflection mask comprises all pixels where:
+ * |e1| <= delta_d, |e2| <= delta_d, |e3| <= delta_m
+ *
+ * We transform the coordinates of the box
+ * (-delta_d, -delta_d, 0)
+ * (+delta_d, -delta_d, 0)
+ * (-delta_d, +delta_d, 0)
+ * (+delta_d, +delta_d, 0)
+ *
+ * to the detector image volume and return the minimum and maximum values
+ * for the x, y, z image volume coordinates.
+ *
+ * @param s1 The diffracted beam vector
+ * @param frame The predicted frame number
+ * @returns A 6 element array: (minx, maxx, miny, maxy, minz, maxz)
+ */
+ virtual int6 single(vec3<double> s0,
+ vec3<double> s1,
+ double frame,
+ double L1,
+ std::size_t panel) const {
+ // Ensure our values are ok
+ DIALS_ASSERT(s1.length_sq() > 0);
+
+ // Create the coordinate system for the reflection
+ CoordinateSystemTOF xcs(s0, s1, L1);
+
+ // Get the divergence and mosaicity for this point
+ double delta_d = delta_divergence_;
+ double delta_m = delta_mosaicity_;
+
+ // Calculate the beam vectors at the following xds coordinates:
+ // (-delta_d, -delta_d, 0)
+ // (+delta_d, -delta_d, 0)
+ // (-delta_d, +delta_d, 0)
+ // (+delta_d, +delta_d, 0)
+ double point = delta_d;
+ double3 sdash1 = xcs.to_beam_vector(double2(-point, -point));
+ double3 sdash2 = xcs.to_beam_vector(double2(+point, -point));
+ double3 sdash3 = xcs.to_beam_vector(double2(-point, +point));
+ double3 sdash4 = xcs.to_beam_vector(double2(+point, +point));
+
+ // Get the detector coordinates (px) at the ray intersections
+ double2 xy1 = detector_[panel].get_ray_intersection_px(sdash1);
+ double2 xy2 = detector_[panel].get_ray_intersection_px(sdash2);
+ double2 xy3 = detector_[panel].get_ray_intersection_px(sdash3);
+ double2 xy4 = detector_[panel].get_ray_intersection_px(sdash4);
+
+ // Return the roi in the following form:
+ // (minx, maxx, miny, maxy, minz, maxz)
+ // Min's are rounded down to the nearest integer, Max's are rounded up
+ double4 x(xy1[0], xy2[0], xy3[0], xy4[0]);
+ double4 y(xy1[1], xy2[1], xy3[1], xy4[1]);
+
+ int x0 = (int)floor(min(x));
+ int x1 = (int)ceil(max(x));
+ int y0 = (int)floor(min(y));
+ int y1 = (int)ceil(max(y));
+
+ double z0 = frame - delta_m;
+ if (z0 < 0) {
+ z0 = 0;
+ }
+ double max_z = scan_.get_array_range()[1];
+ double z1 = frame + delta_m;
+ if (z1 > max_z) {
+ z1 = max_z;
+ }
+
+ int6 bbox(x0, x1, y0, y1, z0, z1);
+ DIALS_ASSERT(bbox[1] > bbox[0]);
+ DIALS_ASSERT(bbox[3] > bbox[2]);
+ DIALS_ASSERT(bbox[5] > bbox[4]);
+ return bbox;
+ }
+
+ /**
+ * Calculate the rois for an array of reflections given by the array of
+ * diffracted beam vectors and rotation angles.
+ * @param s1 The array of diffracted beam vectors
+ * @param phi The array of rotation angles.
+ */
+ virtual af::shared<int6> array(const af::const_ref<vec3<double> > &s0,
+ const af::const_ref<vec3<double> > &s1,
+ const af::const_ref<double> &frame,
+ const af::const_ref<double> &L1,
+ const af::const_ref<std::size_t> &panel) const {
+ DIALS_ASSERT(s1.size() == frame.size());
+ DIALS_ASSERT(s1.size() == panel.size());
+ DIALS_ASSERT(s0.size() == frame.size());
+ DIALS_ASSERT(s0.size() == panel.size());
+ af::shared<int6> result(s1.size(), af::init_functor_null<int6>());
+ for (std::size_t i = 0; i < s1.size(); ++i) {
+ result[i] = single(s0[i], s1[i], frame[i], L1[i], panel[i]);
+ }
+ return result;
+ }
+
+ private:
+ Detector detector_;
+ Scan scan_;
+ PolychromaticBeam beam_;
+ double delta_divergence_;
+ double delta_mosaicity_;
+ };
+
  /**
  * Class to help compute bbox for multiple experiments.
  */

src/dials/algorithms/profile_model/gaussian_rs/boost_python/gaussian_rs_ext.cc

@@ -212,6 +212,21 @@ namespace dials {
  arg("delta_divergence"),
  arg("delta_mosaicity"))));
 
+ class_<BBoxCalculatorTOF>("BBoxCalculatorTOF", no_init)
+ .def(
+ init<const PolychromaticBeam&, const Detector&, const Scan&, double, double>(
+ (arg("beam"),
+ arg("detector"),
+ arg("scan"),
+ arg("delta_divergence"),
+ arg("delta_mosaicity"))))
+ .def("__call__",
+ &BBoxCalculatorTOF::single,
+ (arg("s0"), arg("s1"), arg("frame"), arg("L1"), arg("panel")))
+ .def("__call__",
+ &BBoxCalculatorTOF::array,
+ (arg("s0"), arg("s1"), arg("frame"), arg("L1"), arg("panel")));
+
  class_<BBoxMultiCalculator>("BBoxMultiCalculator")
  .def("append", &BBoxMultiCalculator::push_back)
  .def("__len__", &BBoxMultiCalculator::size)

src/dials/algorithms/profile_model/gaussian_rs/coordinate_system.h

@@ -374,6 +374,132 @@ namespace dials {
  double zeta_;
  };
 
+ /**
+ * Class representing the local reflection coordinate system
+ */
+ class CoordinateSystemTOF {
+ public:
+ /**
+ * Initialise coordinate system. s0 should be the same length as s1.
+ * These quantities are not checked because this class will be created for
+ * each reflection and we want to maximize performance.
+ * @param s0 The incident beam vector
+ * @param s1 The diffracted beam vector
+ */
+ CoordinateSystemTOF(vec3<double> s0, vec3<double> s1, double L1)
+ : s0_(s0),
+ s1_(s1),
+ p_star_(s1 - s0),
+ L1_(L1),
+ e1_(s1.cross(s0).normalize()),
+ e2_(s1.cross(e1_).normalize()),
+ e3_((s1 + s0).normalize()) {}
+
+ vec3<double> s0() const {
+ return s0_;
+ }
+ vec3<double> s1() const {
+ return s1_;
+ }
+ double L1() const {
+ return L1_;
+ }
+ vec3<double> p_star() const {
+ return p_star_;
+ }
+ vec3<double> e1_axis() const {
+ return e1_;
+ }
+ vec3<double> e2_axis() const {
+ return e2_;
+ }
+ vec3<double> e3_axis() const {
+ return e3_;
+ }
+
+ /**
+ * Transform the beam vector to the reciprocal space coordinate system.
+ * @param s_dash The beam vector
+ * @param s0_dash The incident beam vector
+ * @returns The e1, e2, e3 coordinates
+ */
+ vec2<double> from_beam_vector(const vec3<double> &s_dash) const {
+ double s1_length = s1_.length();
+ double s0_length = s0_.length();
+ DIALS_ASSERT(s1_length > 0);
+ DIALS_ASSERT(s0_length > 0);
+ // vec3<double> p_star0 = s_dash-s0_dash;
+ vec3<double> e1 = e1_ / s1_length;
+ vec3<double> e2 = e2_ / s1_length;
+ // vec3<double> e3 = (s1_+s0_)/(s1_length + s0_length);
+ /*
+ return vec3<double>(
+ e1 * (s_dash - s1_),
+ e2 * (s_dash - s1_),
+ e3 * (p_star0 - p_star_)/p_star_.length());
+ */
+ return vec2<double>(e1 * (s_dash - s1_), e2 * (s_dash - s1_));
+ }
+
+ /**
+ * Transform the reciprocal space coordinate to get the beam vector.
+ * @param c12 The e1 and e2 coordinates.
+ * @returns The beam vector
+ */
+ vec3<double> to_beam_vector(const vec2<double> &c12) const {
+ double radius = s1_.length();
+ DIALS_ASSERT(radius > 0);
+ vec3<double> scaled_e1 = e1_ * radius;
+ vec3<double> scaled_e2 = e2_ * radius;
+ vec3<double> normalized_s1 = s1_ / radius;
+
+ vec3<double> p = c12[0] * scaled_e1 + c12[1] * scaled_e2;
+ double b = radius * radius - p.length_sq();
+ DIALS_ASSERT(b >= 0);
+ double d = -(normalized_s1 * p) + std::sqrt(b);
+ return p + d * normalized_s1;
+ }
+
+ /**
+ * @param c3 The XDS e3 coordinate
+ * @param s_dash The beam vector from the e1 and e2 coordinates
+ * @returns The wavelength of the e3 coordinate (s)
+ *
+ * Solved be rearranging c3 = e3(p_star0 - p_star) / s1_length,
+ * noting that p_star0 = s_dash - s0_dash,
+ * and s0_dash = unit_s0/wavelength_dash
+ */
+ double to_wavelength(double c3, vec3<double> s_dash) const {
+ double p_star_length = p_star_.length();
+ DIALS_ASSERT(p_star_length > 0);
+ vec3<double> unit_s0 = s0_.normalize();
+ return (e3_ * unit_s0) / (e3_ * s_dash - e3_ * p_star_ - c3 * p_star_length);
+ }
+
+ /**
+ * Transform the rotation angle to the reciprocal space coordinate system
+ * @param phi_dash The rotation angle
+ * @returns The e3 coordinate.
+ */
+ double from_wavelength(double wavelength) const {
+ double p_star_length = p_star_.length();
+ DIALS_ASSERT(p_star_length > 0);
+ vec3<double> scaled_e3 = e3_ / p_star_length;
+ vec3<double> s0 = s0_.normalize() / wavelength;
+ vec3<double> p_star0 = s1_ - s0;
+ return scaled_e3 * (p_star0 - p_star_);
+ }
+
+ private:
+ vec3<double> s0_;
+ vec3<double> s1_;
+ vec3<double> p_star_;
+ vec3<double> e1_;
+ vec3<double> e2_;
+ vec3<double> e3_;
+ double L1_;
+ };
+
 }}}} // namespace dials::algorithms::profile_model::gaussian_rs
 
 #endif // DIALS_ALGORITHMS_PROFILE_MODEL_GAUSSIAN_RS_COORDINATE_SYSTEM_H

src/dials/algorithms/profile_model/gaussian_rs/model.py

@@ -515,9 +515,20 @@ def compute_bbox(
  )
 
  # Calculate the bounding boxes of all the reflections
- bbox = calculate(
- reflections["s1"], reflections["xyzcal.px"].parts()[2], reflections["panel"]
- )
+ if scan.has_property("time_of_flight"):
+ bbox = calculate(
+ reflections["s0_cal"],
+ reflections["s1"],
+ reflections["xyzcal.px"].parts()[2],
+ reflections["L1"],
+ reflections["panel"],
+ )
+ else:
+ bbox = calculate(
+ reflections["s1"],
+ reflections["xyzcal.px"].parts()[2],
+ reflections["panel"],
+ )
 
  # Return the bounding boxes
  return bbox