Integration improvements

src/dials/algorithms/integration/fit/tof_line_profile.py

@@ -50,7 +50,10 @@ def fit(self):
  xdata=self.tof,
  ydata=self.intensities,
  p0=self.params,
- bounds=((1, 0, 0, 1, min(self.tof)), (1000, 1, 10, 1000, max(self.tof))),
+ bounds=(
+ (1, 0, 0, 1, min(self.tof)),
+ (100000000, 1, 100000, 10000000, max(self.tof)),
+ ),
  max_nfev=10000000,
  )
  self.params = params
@@ -64,6 +67,88 @@ def calc_intensity(self):
  return integrate.simpson(predicted, self.tof)
 
 
+def compute_line_profile_data_for_reflection(reflection_table):
+
+ assert len(reflection_table) == 1
+ A = 200.0
+ alpha = 0.4
+ beta = 0.4
+ sigma = 8.0
+
+ bg_code = MaskCode.Valid | MaskCode.Background | MaskCode.BackgroundUsed
+
+ shoebox = reflection_table["shoebox"][0]
+ data = flumpy.to_numpy(shoebox.data).ravel()
+ background = flumpy.to_numpy(shoebox.background).ravel()
+ mask = flumpy.to_numpy(shoebox.mask).ravel()
+ coords = flumpy.to_numpy(shoebox.coords())
+ m = mask & MaskCode.Foreground == MaskCode.Foreground
+ bg_m = mask & bg_code == bg_code
+ n_background = np.sum(np.bitwise_and(~m, bg_m))
+
+ m = np.bitwise_and(m, mask & MaskCode.Valid == MaskCode.Valid)
+ m = np.bitwise_and(m, mask & MaskCode.Overlapped == 0)
+
+ n_signal = np.sum(m)
+
+ background = background[m]
+ intensity = data[m] - background
+ background_sum = np.sum(background)
+ summation_intensity = float(np.sum(intensity))
+ coords = coords[m]
+ tof = coords[:, 2]
+
+ summed_values = {}
+ summed_background_values = {}
+
+ for j in np.unique(tof):
+ indices = np.where(tof == j)
+ summed_values[j] = np.sum(intensity[indices])
+ summed_background_values[j] = np.sum(background[indices])
+
+ # Remove background and project onto ToF axis
+ projected_intensity = np.array(list(summed_values.values()))
+ projected_background = np.array(list(summed_background_values.values()))
+ tof = np.array(list(summed_values.keys()))
+
+ try:
+ T = tof[np.argmax(projected_intensity)]
+ l = BackToBackExponential(
+ tof=tof,
+ intensities=projected_intensity,
+ A=A,
+ alpha=alpha,
+ beta=beta,
+ sigma=sigma,
+ T=T,
+ )
+ l.fit()
+ line_profile = l.result()
+ fit_intensity = integrate.simpson(line_profile, tof)
+ except ValueError:
+ return [], [], [], [], -1, -1, -1, -1
+
+ if n_background > 0:
+ m_n = n_signal / n_background
+ else:
+ m_n = 0.0
+ fit_std = np.sqrt(abs(fit_intensity) + abs(background_sum) * (1.0 + m_n))
+ summation_std = np.sqrt(
+ abs(summation_intensity) + abs(background_sum) * (1.0 + m_n)
+ )
+
+ return (
+ tof,
+ projected_intensity,
+ projected_background,
+ line_profile,
+ fit_intensity,
+ fit_std,
+ summation_intensity,
+ summation_std,
+ )
+
+
 def compute_line_profile_intensity(reflections):
 
  A = 200.0
@@ -107,6 +192,7 @@ def compute_line_profile_intensity(reflections):
  projected_intensity = np.array(list(summed_values.values()))
  tof = np.array(list(summed_values.keys()))
 
+ fit_intensity = None
  try:
  T = tof[np.argmax(projected_intensity)]
  l = BackToBackExponential(
@@ -123,13 +209,20 @@ def compute_line_profile_intensity(reflections):
  fit_intensities[i] = fit_intensity
  except ValueError:
  fit_intensities[i] = -1
+ fit_variances[i] = -1
+ continue
 
  if n_background > 0:
  m_n = n_signal / n_background
  else:
  m_n = 0.0
- fit_variances[i] = abs(fit_intensity) + abs(background_sum) * (1.0 + m_n)
+ 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.set_flags(
+ reflections["line_profile_intensity"] < 0,
+ reflections.flags.failed_during_profile_fitting,
+ )
  return reflections

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

@@ -452,7 +452,7 @@ namespace dials {
  double z2 = sequence_.get_frame_from_tof(tof2);
  double2 z(z1, z2);
 
- int6 bbox(x0, x1, y0, y1, frame - 5, frame + 7);
+ int6 bbox(x0, x1, y0, y1, frame - 30, frame + 30);
  DIALS_ASSERT(bbox[1] > bbox[0]);
  DIALS_ASSERT(bbox[3] > bbox[2]);
  DIALS_ASSERT(bbox[5] > bbox[4]);

src/dials/algorithms/scaling/boost_python/scaling_ext.cc

@@ -19,6 +19,8 @@ namespace dials_scaling { namespace boost_python {
  void export_create_sph_harm_lookup_table();
  void export_gaussian_smoother_first_fixed();
  void export_limit_outlier_weights();
+ void export_tof_spherical_absorption_correction();
+ void export_tof_lorentz_correction();
 
  BOOST_PYTHON_MODULE(dials_scaling_ext) {
  export_elementwise_square();
@@ -35,6 +37,8 @@ namespace dials_scaling { namespace boost_python {
  export_create_sph_harm_lookup_table();
  export_gaussian_smoother_first_fixed();
  export_limit_outlier_weights();
+ export_tof_spherical_absorption_correction();
+ export_tof_lorentz_correction();
  }
 
 }} // namespace dials_scaling::boost_python

src/dials/algorithms/scaling/boost_python/scaling_helper.cc

@@ -1,6 +1,7 @@
 #include <boost/python.hpp>
 #include <boost/python/def.hpp>
 #include <dials/algorithms/scaling/scaling_helper.h>
+#include <dials/algorithms/scaling/tof_absorption_correction.h>
 
 namespace dials_scaling { namespace boost_python {
 
@@ -101,4 +102,16 @@ namespace dials_scaling { namespace boost_python {
  &GaussianSmootherFirstFixed::multi_value_weight_first_fixed);
  }
 
+ void export_tof_spherical_absorption_correction() {
+ def("tof_spherical_absorption_correction",
+ &dials::algorithms::tof_spherical_absorption_correction,
+ (arg("spectra"), arg("muR"), arg("two_thetas"), arg("two_theta_idxs")));
+ }
+ void export_tof_lorentz_correction() {
+ def(
+ "tof_lorentz_correction",
+ &dials::algorithms::tof_lorentz_correction,
+ (arg("spectra"), arg("L0"), arg("L1"), arg("tof"), arg("two_theta_spectra_sq")));
+ }
+
 }} // namespace dials_scaling::boost_python