Integrate3DEvents.cpp

Go to the documentation of this file.

// Mantid Repository : https://github.com/mantidproject/mantid
//
// Copyright © 2018 ISIS Rutherford Appleton Laboratory UKRI,
//   NScD Oak Ridge National Laboratory, European Spallation Source,
//   Institut Laue - Langevin & CSNS, Institute of High Energy Physics, CAS
// SPDX - License - Identifier: GPL - 3.0 +
#include "MantidMDAlgorithms/Integrate3DEvents.h"
#include "MantidDataObjects/NoShape.h"
#include "MantidDataObjects/PeakShapeEllipsoid.h"
#include "MantidGeometry/Crystal/IndexingUtils.h"
 
#include <boost/math/special_functions/round.hpp>
#include <cmath>
#include <fstream>
#include <memory>
#include <numeric>
#include <tuple>
 
extern "C" {
#include <cstdio>
#include <utility>
 
#include <gsl/gsl_eigen.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
}
 
using namespace Mantid::DataObjects;
namespace Mantid::MDAlgorithms {
 
using namespace std;
using Mantid::Kernel::DblMatrix;
using Mantid::Kernel::V3D;
 
Integrate3DEvents::Integrate3DEvents(
    const std::vector<std::pair<std::pair<double, double>, Mantid::Kernel::V3D>> &peak_q_list, Kernel::DblMatrix UBinv,
    double radius, const bool useOnePercentBackgroundCorrection)
    : m_UBinv(std::move(UBinv)), m_radius(radius), maxOrder(0), crossterm(false),
      m_useOnePercentBackgroundCorrection(useOnePercentBackgroundCorrection) {
  for (size_t it = 0; it != peak_q_list.size(); ++it) {
    int64_t hkl_key = getHklKey(peak_q_list[it].second);
    if (hkl_key != 0) // only save if hkl != (0,0,0)
      m_peak_qs[hkl_key] = peak_q_list[it].second;
  }
}
 
Integrate3DEvents::Integrate3DEvents(
    const std::vector<std::pair<std::pair<double, double>, Mantid::Kernel::V3D>> &peak_q_list,
    std::vector<V3D> const &hkl_list, std::vector<V3D> const &mnp_list, Kernel::DblMatrix UBinv,
    Kernel::DblMatrix ModHKL, double radius_m, double radius_s, int MaxO, const bool CrossT,
    const bool useOnePercentBackgroundCorrection)
    : m_UBinv(std::move(UBinv)), m_ModHKL(std::move(ModHKL)), m_radius(radius_m), s_radius(radius_s), maxOrder(MaxO),
      crossterm(CrossT), m_useOnePercentBackgroundCorrection(useOnePercentBackgroundCorrection) {
  for (size_t it = 0; it != peak_q_list.size(); ++it) {
    int64_t hklmnp_key =
        getHklMnpKey(boost::math::iround<double>(hkl_list[it][0]), boost::math::iround<double>(hkl_list[it][1]),
                     boost::math::iround<double>(hkl_list[it][2]), boost::math::iround<double>(mnp_list[it][0]),
                     boost::math::iround<double>(mnp_list[it][1]), boost::math::iround<double>(mnp_list[it][2]));
    if (hklmnp_key != 0) // only save if hkl != (0,0,0)
      m_peak_qs[hklmnp_key] = peak_q_list[it].second;
  }
}
 
void Integrate3DEvents::addEvents(std::vector<std::pair<std::pair<double, double>, V3D>> const &event_qs,
                                  bool hkl_integ) {
  if (!maxOrder)
    for (const auto &event_q : event_qs)
      addEvent(event_q, hkl_integ);
  else
    for (const auto &event_q : event_qs)
      addModEvent(event_q, hkl_integ);
}
 
std::pair<std::shared_ptr<const Geometry::PeakShape>, std::tuple<double, double, double>>
Integrate3DEvents::integrateStrongPeak(const IntegrationParameters &params, const V3D &peak_q, double &inti,
                                       double &sigi) {
 
  inti = 0.0; // default values, in case something
  sigi = 0.0; // is wrong with the peak.
  auto result = getEvents(peak_q);
  if (!result)
    return std::make_pair(std::make_shared<NoShape>(), make_tuple(0., 0., 0.));
 
  const auto &events = *result;
  if (events.empty())
    return std::make_pair(std::make_shared<NoShape>(), make_tuple(0., 0., 0.));
 
  DblMatrix cov_matrix(3, 3);
  makeCovarianceMatrix(events, cov_matrix, params.regionRadius);
 
  std::vector<V3D> eigen_vectors;
  std::vector<double> eigen_values;
  getEigenVectors(cov_matrix, eigen_vectors, eigen_values);
 
  std::vector<double> sigmas(3);
  for (int i = 0; i < 3; i++) {
    sigmas[i] = sqrt(eigen_values[i]);
  }
 
  bool invalid_peak =
      std::any_of(sigmas.cbegin(), sigmas.cend(), [](const double sigma) { return std::isnan(sigma) || sigma <= 0; });
 
  if (invalid_peak)
    return std::make_pair(std::make_shared<NoShape>(), make_tuple(0., 0., 0.));
 
  const auto max_sigma = *std::max_element(sigmas.begin(), sigmas.end());
  if (max_sigma == 0)
    return std::make_pair(std::make_shared<NoShape>(), make_tuple(0., 0., 0.));
 
  auto rValues = calculateRadiusFactors(params, max_sigma);
  auto &r1 = std::get<0>(rValues), r2 = std::get<1>(rValues), r3 = std::get<2>(rValues);
 
  std::vector<double> abcBackgroundOuterRadii, abcBackgroundInnerRadii;
  std::vector<double> peakRadii;
  for (int i = 0; i < 3; i++) {
    abcBackgroundOuterRadii.emplace_back(r3 * sigmas[i]);
    abcBackgroundInnerRadii.emplace_back(r2 * sigmas[i]);
    peakRadii.emplace_back(r1 * sigmas[i]);
  }
 
  const auto isPeakOnDetector = correctForDetectorEdges(rValues, params.E1Vectors, peak_q, peakRadii,
                                                        abcBackgroundInnerRadii, abcBackgroundOuterRadii);
 
  if (!isPeakOnDetector)
    return std::make_pair(std::make_shared<NoShape>(), make_tuple(0.0, 0.0, 0.));
 
  const auto backgrd = numInEllipsoidBkg(events, eigen_vectors, abcBackgroundOuterRadii, abcBackgroundInnerRadii,
                                         m_useOnePercentBackgroundCorrection);
  const auto core = numInEllipsoid(events, eigen_vectors, sigmas);
  const auto peak = numInEllipsoid(events, eigen_vectors, peakRadii);
  const auto ratio = pow(r1, 3) / (pow(r3, 3) - pow(r2, 3));
 
  inti = peak.first - ratio * backgrd.first;
  sigi = sqrt(peak.second + ratio * ratio * backgrd.second);
 
  // compute the fraction of peak within the standard core
  const auto total = (core.first + peak.first) - ratio * backgrd.first;
  const auto frac = std::min(1.0, std::abs(inti / total));
  // compute the uncertainty in the fraction
  const auto df_ds_core = (1 - frac) / peak.first;
  const auto df_ds_peak = frac / peak.first;
  const auto fracError = sqrt(peak.first * pow(df_ds_core, 2) + core.first * pow(df_ds_peak, 2));
 
  // create the peaks shape for the strong peak
  const auto shape = std::make_shared<const PeakShapeEllipsoid>(eigen_vectors, peakRadii, abcBackgroundInnerRadii,
                                                                abcBackgroundOuterRadii, Mantid::Kernel::QLab,
                                                                "IntegrateEllipsoidsTwoStep");
 
  return std::make_pair(shape, std::make_tuple(frac, fracError, max_sigma));
}
 
std::shared_ptr<const Geometry::PeakShape>
Integrate3DEvents::integrateWeakPeak(const IntegrationParameters &params, PeakShapeEllipsoid_const_sptr shape,
                                     const std::tuple<double, double, double> &libPeak, const V3D &center, double &inti,
                                     double &sigi) {
 
  inti = 0.0; // default values, in case something
  sigi = 0.0; // is wrong with the peak.
 
  auto result = getEvents(center);
  if (!result)
    return std::make_shared<NoShape>();
 
  const auto &events = *result;
 
  const auto &directions = shape->directions();
  auto abcBackgroundInnerRadii = shape->abcRadiiBackgroundInner();
  auto abcBackgroundOuterRadii = shape->abcRadiiBackgroundOuter();
  auto abcRadii = shape->abcRadii();
 
  const auto max_sigma = std::get<2>(libPeak);
  auto rValues = calculateRadiusFactors(params, max_sigma);
 
  const auto isPeakOnDetector = correctForDetectorEdges(rValues, params.E1Vectors, center, abcRadii,
                                                        abcBackgroundInnerRadii, abcBackgroundOuterRadii);
 
  if (!isPeakOnDetector)
    return shape;
 
  const double r1 = std::get<0>(rValues), r2 = std::get<1>(rValues), r3 = std::get<2>(rValues);
 
  // integrate
  std::pair<double, double> backgrd = numInEllipsoidBkg(events, directions, abcBackgroundOuterRadii,
                                                        abcBackgroundInnerRadii, m_useOnePercentBackgroundCorrection);
  std::pair<double, double> peak_w_back = numInEllipsoid(events, directions, abcRadii);
  double ratio = pow(r1, 3) / (pow(r3, 3) - pow(r2, 3));
 
  const auto frac = std::get<0>(libPeak);
  const auto fracError = std::get<1>(libPeak);
 
  inti = peak_w_back.first - ratio * backgrd.first;
 
  // correct for fractional intensity
  sigi = sigi / pow(inti, 2);
  sigi += pow((fracError / frac), 2);
 
  inti = inti * frac;
  sigi = sqrt(sigi) * inti;
 
  // scale integration shape by fractional amount
  for (size_t i = 0; i < abcRadii.size(); ++i) {
    abcRadii[i] *= frac;
    abcBackgroundInnerRadii[i] *= frac;
    abcBackgroundOuterRadii[i] *= frac;
  }
 
  return std::make_shared<const PeakShapeEllipsoid>(shape->directions(), abcRadii, abcBackgroundInnerRadii,
                                                    abcBackgroundOuterRadii, Mantid::Kernel::QLab,
                                                    "IntegrateEllipsoidsTwoStep");
}
 
double Integrate3DEvents::estimateSignalToNoiseRatio(const IntegrationParameters &params, const V3D &center,
                                                     bool forceSpherical, double sphericityTol) {
 
  auto result = getEvents(center);
  if (!result)
    return .0;
 
  const auto &events = *result;
  if (events.empty())
    return .0;
 
  DblMatrix cov_matrix(3, 3);
  makeCovarianceMatrix(events, cov_matrix, params.regionRadius);
 
  std::vector<V3D> eigen_vectors;
  std::vector<double> eigen_values;
  getEigenVectors(cov_matrix, eigen_vectors, eigen_values);
 
  std::vector<double> sigmas(3);
  for (int i = 0; i < 3; i++) {
    sigmas[i] = sqrt(eigen_values[i]);
  }
 
  const auto max_sigma = *std::max_element(sigmas.begin(), sigmas.end());
  const auto min_sigma = *std::min_element(sigmas.begin(), sigmas.end());
  if (max_sigma == 0)
    return .0;
 
  auto rValues = calculateRadiusFactors(params, max_sigma);
  auto &r1 = std::get<0>(rValues), r2 = std::get<1>(rValues), r3 = std::get<2>(rValues);
  std::vector<double> abcBackgroundOuterRadii, abcBackgroundInnerRadii;
  std::vector<double> peakRadii;
  if (forceSpherical) {
    // test for spherically symmeteric peak (within tolerance)
    if ((max_sigma - min_sigma) / max_sigma > sphericityTol)
      return .0;
    for (int i = 0; i < 3; i++) {
      abcBackgroundOuterRadii.emplace_back(r3 * max_sigma);
      abcBackgroundInnerRadii.emplace_back(r2 * max_sigma);
      peakRadii.emplace_back(r1 * max_sigma);
    }
  } else {
    for (int i = 0; i < 3; i++) {
      abcBackgroundOuterRadii.emplace_back(r3 * sigmas[i]);
      abcBackgroundInnerRadii.emplace_back(r2 * sigmas[i]);
      peakRadii.emplace_back(r1 * sigmas[i]);
    }
  }
 
  // Background / Peak / Background
  std::pair<double, double> backgrd = numInEllipsoidBkg(events, eigen_vectors, abcBackgroundOuterRadii,
                                                        abcBackgroundInnerRadii, m_useOnePercentBackgroundCorrection);
 
  std::pair<double, double> peak_w_back = numInEllipsoid(events, eigen_vectors, peakRadii);
 
  double ratio = pow(r1, 3) / (pow(r3, 3) - pow(r2, 3));
  auto inti = peak_w_back.first - ratio * backgrd.first;
  auto sigi = sqrt(peak_w_back.second + ratio * ratio * backgrd.second);
 
  return inti / sigi;
}
 
const std::vector<std::pair<std::pair<double, double>, V3D>> *Integrate3DEvents::getEvents(const V3D &peak_q) {
  auto hkl_key = getHklKey(peak_q);
  if (maxOrder)
    hkl_key = getHklMnpKey(peak_q);
 
  if (hkl_key == 0)
    return nullptr;
 
  const auto pos = m_event_lists.find(hkl_key);
 
  if (m_event_lists.end() == pos)
    return nullptr;
 
  if (pos->second.size() < 3) // if there are not enough events
    return nullptr;
 
  return &(pos->second);
}
 
bool Integrate3DEvents::correctForDetectorEdges(std::tuple<double, double, double> &radii,
                                                const std::vector<V3D> &E1Vecs, const V3D &peak_q,
                                                const std::vector<double> &axesRadii,
                                                const std::vector<double> &bkgInnerRadii,
                                                const std::vector<double> &bkgOuterRadii) {
 
  if (E1Vecs.empty())
    return true;
 
  const auto &r1 = std::get<0>(radii);
  auto &r2 = std::get<1>(radii);
  auto &r3 = std::get<2>(radii);
  auto h3 = 1.0 - detectorQ(E1Vecs, peak_q, bkgOuterRadii);
  // scaled from area of circle minus segment when r normalized to 1
  auto m3 = std::sqrt(1.0 - (std::acos(1.0 - h3) - (1.0 - h3) * std::sqrt(2.0 * h3 - h3 * h3)) / M_PI);
  auto h1 = 1.0 - detectorQ(E1Vecs, peak_q, axesRadii);
  // Do not use peak if edge of detector is inside integration radius
  if (h1 > 0.0)
    return false;
 
  r3 *= m3;
  if (r2 != r1) {
    auto h2 = 1.0 - detectorQ(E1Vecs, peak_q, bkgInnerRadii);
    // scaled from area of circle minus segment when r normalized to 1
    auto m2 = std::sqrt(1.0 - (std::acos(1.0 - h2) - (1.0 - h2) * std::sqrt(2.0 * h2 - h2 * h2)) / M_PI);
    r2 *= m2;
  }
 
  return true;
}
 
Mantid::Geometry::PeakShape_const_sptr
Integrate3DEvents::ellipseIntegrateEvents(const std::vector<V3D> &E1Vec, V3D const &peak_q, bool specify_size,
                                          double peak_radius, double back_inner_radius, double back_outer_radius,
                                          std::vector<double> &axes_radii, double &inti, double &sigi) {
  inti = 0.0; // default values, in case something
  sigi = 0.0; // is wrong with the peak.
 
  int64_t hkl_key = getHklKey(peak_q);
 
  if (hkl_key == 0) {
    return std::make_shared<NoShape>();
  }
 
  auto pos = m_event_lists.find(hkl_key);
  if (m_event_lists.end() == pos)
    return std::make_shared<NoShape>();
  ;
 
  const std::vector<std::pair<std::pair<double, double>, V3D>> &some_events = pos->second;
 
  if (some_events.size() < 3) // if there are not enough events to
  {                           // find covariance matrix, return
    return std::make_shared<NoShape>();
  }
 
  DblMatrix cov_matrix(3, 3);
  makeCovarianceMatrix(some_events, cov_matrix, m_radius);
 
  std::vector<V3D> eigen_vectors;
  std::vector<double> eigen_values;
  getEigenVectors(cov_matrix, eigen_vectors, eigen_values);
 
  std::vector<double> sigmas(3);
  for (int i = 0; i < 3; i++) {
    sigmas[i] = sqrt(eigen_values[i]);
  }
 
  bool invalid_peak =
      std::any_of(sigmas.cbegin(), sigmas.cend(), [](const double sigma) { return std::isnan(sigma) || sigma <= 0; });
 
  if (invalid_peak)                     // if data collapses to a line or
  {                                     // to a plane, the volume of the
    return std::make_shared<NoShape>(); // ellipsoids will be zero.
  }
 
  return ellipseIntegrateEvents(E1Vec, peak_q, some_events, eigen_vectors, sigmas, specify_size, peak_radius,
                                back_inner_radius, back_outer_radius, axes_radii, inti, sigi);
}
 
Mantid::Geometry::PeakShape_const_sptr
Integrate3DEvents::ellipseIntegrateModEvents(const std::vector<V3D> &E1Vec, V3D const &peak_q, V3D const &hkl,
                                             V3D const &mnp, bool specify_size, double peak_radius,
                                             double back_inner_radius, double back_outer_radius,
                                             std::vector<double> &axes_radii, double &inti, double &sigi) {
  inti = 0.0; // default values, in case something
  sigi = 0.0; // is wrong with the peak.
 
  int64_t hkl_key = getHklMnpKey(boost::math::iround<double>(hkl[0]), boost::math::iround<double>(hkl[1]),
                                 boost::math::iround<double>(hkl[2]), boost::math::iround<double>(mnp[0]),
                                 boost::math::iround<double>(mnp[1]), boost::math::iround<double>(mnp[2]));
 
  if (hkl_key == 0) {
    return std::make_shared<NoShape>();
  }
 
  auto pos = m_event_lists.find(hkl_key);
  if (m_event_lists.end() == pos)
    return std::make_shared<NoShape>();
  ;
 
  const std::vector<std::pair<std::pair<double, double>, V3D>> &some_events = pos->second;
 
  if (some_events.size() < 3) // if there are not enough events to
  {                           // find covariance matrix, return
    return std::make_shared<NoShape>();
  }
 
  DblMatrix cov_matrix(3, 3);
  if (hkl_key % 1000 == 0)
    makeCovarianceMatrix(some_events, cov_matrix, m_radius);
  else
    makeCovarianceMatrix(some_events, cov_matrix, s_radius);
 
  std::vector<V3D> eigen_vectors;
  std::vector<double> eigen_values;
  getEigenVectors(cov_matrix, eigen_vectors, eigen_values);
 
  std::vector<double> sigmas(3);
  for (int i = 0; i < 3; i++)
    sigmas[i] = sqrt(eigen_values[i]);
 
  bool invalid_peak =
      std::any_of(sigmas.cbegin(), sigmas.cend(), [](const double sigma) { return std::isnan(sigma) || sigma <= 0; });
 
  if (invalid_peak)                     // if data collapses to a line or
  {                                     // to a plane, the volume of the
    return std::make_shared<NoShape>(); // ellipsoids will be zero.
  }
 
  return ellipseIntegrateEvents(E1Vec, peak_q, some_events, eigen_vectors, sigmas, specify_size, peak_radius,
                                back_inner_radius, back_outer_radius, axes_radii, inti, sigi);
}
std::pair<double, double>
Integrate3DEvents::numInEllipsoid(std::vector<std::pair<std::pair<double, double>, V3D>> const &events,
                                  std::vector<V3D> const &directions, std::vector<double> const &sizes) {
 
  std::pair<double, double> count(0, 0);
  for (const auto &event : events) {
    double sum = 0;
    for (size_t k = 0; k < 3; k++) {
      double comp = event.second.scalar_prod(directions[k]) / sizes[k];
      sum += comp * comp;
    }
    if (sum <= 1) {
      count.first += event.first.first;   // count
      count.second += event.first.second; // error squared (add in quadrature)
    }
  }
 
  return count;
}
std::pair<double, double>
Integrate3DEvents::numInEllipsoidBkg(std::vector<std::pair<std::pair<double, double>, V3D>> const &events,
                                     std::vector<V3D> const &directions, std::vector<double> const &sizes,
                                     std::vector<double> const &sizesIn, const bool useOnePercentBackgroundCorrection) {
  std::pair<double, double> count(0, 0);
  std::vector<std::pair<double, double>> eventVec;
  for (const auto &event : events) {
    double sum = 0;
    double sumIn = 0;
    for (size_t k = 0; k < 3; k++) {
      double comp = event.second.scalar_prod(directions[k]) / sizes[k];
      sum += comp * comp;
      comp = event.second.scalar_prod(directions[k]) / sizesIn[k];
      sumIn += comp * comp;
    }
    if (sum <= 1 && sumIn >= 1)
      eventVec.emplace_back(event.first);
  }
 
  auto endIndex = eventVec.size();
  if (useOnePercentBackgroundCorrection) {
    // Remove top 1% of background
    std::sort(eventVec.begin(), eventVec.end(),
              [](const std::pair<double, double> &a, const std::pair<double, double> &b) { return a.first < b.first; });
    endIndex = static_cast<size_t>(0.99 * static_cast<double>(endIndex));
  }
 
  for (size_t k = 0; k < endIndex; ++k) {
    count.first += eventVec[k].first;
    count.second += eventVec[k].second;
  }
 
  return count;
}
void Integrate3DEvents::makeCovarianceMatrix(std::vector<std::pair<std::pair<double, double>, V3D>> const &events,
                                             DblMatrix &matrix, double radius) {
  double totalCounts;
  for (int row = 0; row < 3; row++) {
    for (int col = 0; col < 3; col++) {
      totalCounts = 0;
      double sum = 0;
      for (const auto &event : events) {
        if (event.second.norm() <= radius) {
          totalCounts += event.first.first;
          sum += event.first.first * event.second[row] * event.second[col];
        }
      }
      if (totalCounts > 1)
        matrix[row][col] = sum / (totalCounts - 1);
      else
        matrix[row][col] = sum;
    }
  }
}
 
void Integrate3DEvents::getEigenVectors(DblMatrix const &cov_matrix, std::vector<V3D> &eigen_vectors,
                                        std::vector<double> &eigen_values) {
  unsigned int size = 3;
 
  gsl_matrix *matrix = gsl_matrix_alloc(size, size);
  gsl_vector *eigen_val = gsl_vector_alloc(size);
  gsl_matrix *eigen_vec = gsl_matrix_alloc(size, size);
  gsl_eigen_symmv_workspace *wkspace = gsl_eigen_symmv_alloc(size);
 
  // copy the matrix data into the gsl matrix
  for (size_t row = 0; row < size; row++)
    for (size_t col = 0; col < size; col++) {
      gsl_matrix_set(matrix, row, col, cov_matrix[row][col]);
    }
 
  gsl_eigen_symmv(matrix, eigen_val, eigen_vec, wkspace);
 
  // copy the resulting eigen vectors to output vector
  for (size_t col = 0; col < size; col++) {
    eigen_vectors.emplace_back(gsl_matrix_get(eigen_vec, 0, col), gsl_matrix_get(eigen_vec, 1, col),
                               gsl_matrix_get(eigen_vec, 2, col));
    eigen_values.emplace_back(gsl_vector_get(eigen_val, col));
  }
 
  gsl_matrix_free(matrix);
  gsl_vector_free(eigen_val);
  gsl_matrix_free(eigen_vec);
  gsl_eigen_symmv_free(wkspace);
}
 
int64_t Integrate3DEvents::getHklKey(int h, int k, int l) {
  int64_t key(0);
 
  if (h != 0 || k != 0 || l != 0)
    key = 1000000000000 * h + 100000000 * k + 10000 * l;
 
  return key;
}
int64_t Integrate3DEvents::getHklMnpKey(int h, int k, int l, int m, int n, int p) {
  int64_t key(0);
 
  if (h != 0 || k != 0 || l != 0 || m != 0 || n != 0 || p != 0)
    key = 1000000000000 * h + 100000000 * k + 10000 * l + 100 * m + 10 * n + p;
 
  return key;
}
int64_t Integrate3DEvents::getHklKey2(V3D const &hkl) {
  int h = boost::math::iround<double>(hkl[0]);
  int k = boost::math::iround<double>(hkl[1]);
  int l = boost::math::iround<double>(hkl[2]);
  return getHklKey(h, k, l);
}
int64_t Integrate3DEvents::getHklMnpKey2(V3D const &hkl) {
  V3D modvec1 = V3D(m_ModHKL[0][0], m_ModHKL[1][0], m_ModHKL[2][0]);
  V3D modvec2 = V3D(m_ModHKL[0][1], m_ModHKL[1][1], m_ModHKL[2][1]);
  V3D modvec3 = V3D(m_ModHKL[0][2], m_ModHKL[1][2], m_ModHKL[2][2]);
  if (Geometry::IndexingUtils::ValidIndex(hkl, m_radius)) {
    int h = boost::math::iround<double>(hkl[0]);
    int k = boost::math::iround<double>(hkl[1]);
    int l = boost::math::iround<double>(hkl[2]);
 
    return getHklMnpKey(h, k, l, 0, 0, 0);
  } else if (!crossterm) {
    if (modvec1 != V3D(0, 0, 0))
      for (int order = -maxOrder; order <= maxOrder; order++) {
        if (order == 0)
          continue; // exclude order 0
        V3D hkl1(hkl);
        hkl1[0] -= order * modvec1[0];
        hkl1[1] -= order * modvec1[1];
        hkl1[2] -= order * modvec1[2];
        if (Geometry::IndexingUtils::ValidIndex(hkl1, s_radius)) {
          int h = boost::math::iround<double>(hkl1[0]);
          int k = boost::math::iround<double>(hkl1[1]);
          int l = boost::math::iround<double>(hkl1[2]);
          return getHklMnpKey(h, k, l, order, 0, 0);
        }
      }
    if (modvec2 != V3D(0, 0, 0))
      for (int order = -maxOrder; order <= maxOrder; order++) {
        if (order == 0)
          continue; // exclude order 0
        V3D hkl1(hkl);
        hkl1[0] -= order * modvec2[0];
        hkl1[1] -= order * modvec2[1];
        hkl1[2] -= order * modvec2[2];
        if (Geometry::IndexingUtils::ValidIndex(hkl1, s_radius)) {
          int h = boost::math::iround<double>(hkl1[0]);
          int k = boost::math::iround<double>(hkl1[1]);
          int l = boost::math::iround<double>(hkl1[2]);
          return getHklMnpKey(h, k, l, 0, order, 0);
        }
      }
    if (modvec3 != V3D(0, 0, 0))
      for (int order = -maxOrder; order <= maxOrder; order++) {
        if (order == 0)
          continue; // exclude order 0
        V3D hkl1(hkl);
        hkl1[0] -= order * modvec3[0];
        hkl1[1] -= order * modvec3[1];
        hkl1[2] -= order * modvec3[2];
        if (Geometry::IndexingUtils::ValidIndex(hkl1, s_radius)) {
          int h = boost::math::iround<double>(hkl1[0]);
          int k = boost::math::iround<double>(hkl1[1]);
          int l = boost::math::iround<double>(hkl1[2]);
          return getHklMnpKey(h, k, l, 0, 0, order);
        }
      }
  } else {
    int maxOrder1 = maxOrder;
    if (modvec1 == V3D(0, 0, 0))
      maxOrder1 = 0;
    int maxOrder2 = maxOrder;
    if (modvec2 == V3D(0, 0, 0))
      maxOrder2 = 0;
    int maxOrder3 = maxOrder;
    if (modvec3 == V3D(0, 0, 0))
      maxOrder3 = 0;
    for (int m = -maxOrder1; m <= maxOrder1; m++)
      for (int n = -maxOrder2; n <= maxOrder2; n++)
        for (int p = -maxOrder3; p <= maxOrder3; p++) {
          if (m == 0 && n == 0 && p == 0)
            continue; // exclude 0,0,0
          V3D hkl1(hkl);
          V3D mnp = V3D(m, n, p);
          hkl1 -= m_ModHKL * mnp;
          if (Geometry::IndexingUtils::ValidIndex(hkl1, s_radius)) {
            int h = boost::math::iround<double>(hkl1[0]);
            int k = boost::math::iround<double>(hkl1[1]);
            int l = boost::math::iround<double>(hkl1[2]);
            return getHklMnpKey(h, k, l, m, n, p);
          }
        }
  }
  return 0;
}
int64_t Integrate3DEvents::getHklKey(V3D const &q_vector) {
  V3D hkl = m_UBinv * q_vector;
  int h = boost::math::iround<double>(hkl[0]);
  int k = boost::math::iround<double>(hkl[1]);
  int l = boost::math::iround<double>(hkl[2]);
  return getHklKey(h, k, l);
}
 
int64_t Integrate3DEvents::getHklMnpKey(V3D const &q_vector) {
  V3D hkl = m_UBinv * q_vector;
 
  V3D modvec1 = V3D(m_ModHKL[0][0], m_ModHKL[1][0], m_ModHKL[2][0]);
  V3D modvec2 = V3D(m_ModHKL[0][1], m_ModHKL[1][1], m_ModHKL[2][1]);
  V3D modvec3 = V3D(m_ModHKL[0][2], m_ModHKL[1][2], m_ModHKL[2][2]);
  if (Geometry::IndexingUtils::ValidIndex(hkl, m_radius)) {
    int h = boost::math::iround<double>(hkl[0]);
    int k = boost::math::iround<double>(hkl[1]);
    int l = boost::math::iround<double>(hkl[2]);
 
    return getHklMnpKey(h, k, l, 0, 0, 0);
  } else if (!crossterm) {
    if (modvec1 != V3D(0, 0, 0))
      for (int order = -maxOrder; order <= maxOrder; order++) {
        if (order == 0)
          continue; // exclude order 0
        V3D hkl1(hkl);
        hkl1[0] -= order * modvec1[0];
        hkl1[1] -= order * modvec1[1];
        hkl1[2] -= order * modvec1[2];
        if (Geometry::IndexingUtils::ValidIndex(hkl1, s_radius)) {
          int h = boost::math::iround<double>(hkl1[0]);
          int k = boost::math::iround<double>(hkl1[1]);
          int l = boost::math::iround<double>(hkl1[2]);
          return getHklMnpKey(h, k, l, order, 0, 0);
        }
      }
    if (modvec2 != V3D(0, 0, 0))
      for (int order = -maxOrder; order <= maxOrder; order++) {
        if (order == 0)
          continue; // exclude order 0
        V3D hkl1(hkl);
        hkl1[0] -= order * modvec2[0];
        hkl1[1] -= order * modvec2[1];
        hkl1[2] -= order * modvec2[2];
        if (Geometry::IndexingUtils::ValidIndex(hkl1, s_radius)) {
          int h = boost::math::iround<double>(hkl1[0]);
          int k = boost::math::iround<double>(hkl1[1]);
          int l = boost::math::iround<double>(hkl1[2]);
          return getHklMnpKey(h, k, l, 0, order, 0);
        }
      }
    if (modvec3 != V3D(0, 0, 0))
      for (int order = -maxOrder; order <= maxOrder; order++) {
        if (order == 0)
          continue; // exclude order 0
        V3D hkl1(hkl);
        hkl1[0] -= order * modvec3[0];
        hkl1[1] -= order * modvec3[1];
        hkl1[2] -= order * modvec3[2];
        if (Geometry::IndexingUtils::ValidIndex(hkl1, s_radius)) {
          int h = boost::math::iround<double>(hkl1[0]);
          int k = boost::math::iround<double>(hkl1[1]);
          int l = boost::math::iround<double>(hkl1[2]);
          return getHklMnpKey(h, k, l, 0, 0, order);
        }
      }
  } else {
    int maxOrder1 = maxOrder;
    if (modvec1 == V3D(0, 0, 0))
      maxOrder1 = 0;
    int maxOrder2 = maxOrder;
    if (modvec2 == V3D(0, 0, 0))
      maxOrder2 = 0;
    int maxOrder3 = maxOrder;
    if (modvec3 == V3D(0, 0, 0))
      maxOrder3 = 0;
    for (int m = -maxOrder1; m <= maxOrder1; m++)
      for (int n = -maxOrder2; n <= maxOrder2; n++)
        for (int p = -maxOrder3; p <= maxOrder3; p++) {
          if (m == 0 && n == 0 && p == 0)
            continue; // exclude 0,0,0
          V3D hkl1(hkl);
          V3D mnp = V3D(m, n, p);
          hkl1 -= m_ModHKL * mnp;
          if (Geometry::IndexingUtils::ValidIndex(hkl1, s_radius)) {
            int h = boost::math::iround<double>(hkl1[0]);
            int k = boost::math::iround<double>(hkl1[1]);
            int l = boost::math::iround<double>(hkl1[2]);
            return getHklMnpKey(h, k, l, m, n, p);
          }
        }
  }
  return 0;
}
 
void Integrate3DEvents::addEvent(std::pair<std::pair<double, double>, V3D> event_Q, bool hkl_integ) {
  int64_t hkl_key;
  if (hkl_integ)
    hkl_key = getHklKey2(event_Q.second);
  else
    hkl_key = getHklKey(event_Q.second);
 
  if (hkl_key == 0) // don't keep events associated with 0,0,0
    return;
 
  auto peak_it = m_peak_qs.find(hkl_key);
  if (peak_it != m_peak_qs.end()) {
    if (!peak_it->second.nullVector()) {
      if (hkl_integ)
        event_Q.second = event_Q.second - m_UBinv * peak_it->second;
      else
        event_Q.second = event_Q.second - peak_it->second;
      if (event_Q.second.norm() < m_radius) {
        m_event_lists[hkl_key].emplace_back(event_Q);
      }
    }
  }
}
 
void Integrate3DEvents::addModEvent(std::pair<std::pair<double, double>, V3D> event_Q, bool hkl_integ) {
  int64_t hklmnp_key;
 
  if (hkl_integ)
    hklmnp_key = getHklMnpKey2(event_Q.second);
  else
    hklmnp_key = getHklMnpKey(event_Q.second);
 
  if (hklmnp_key == 0) // don't keep events associated with 0,0,0
    return;
 
  auto peak_it = m_peak_qs.find(hklmnp_key);
  if (peak_it != m_peak_qs.end()) {
    if (!peak_it->second.nullVector()) {
      if (hkl_integ)
        event_Q.second = event_Q.second - m_UBinv * peak_it->second;
      else
        event_Q.second = event_Q.second - peak_it->second;
 
      if (hklmnp_key % 10000 == 0) {
        if (event_Q.second.norm() < m_radius)
          m_event_lists[hklmnp_key].emplace_back(event_Q);
      } else if (event_Q.second.norm() < s_radius) {
        m_event_lists[hklmnp_key].emplace_back(event_Q);
      }
    }
  }
}
 
PeakShapeEllipsoid_const_sptr Integrate3DEvents::ellipseIntegrateEvents(
    const std::vector<V3D> &E1Vec, V3D const &peak_q,
    std::vector<std::pair<std::pair<double, double>, Mantid::Kernel::V3D>> const &ev_list,
    std::vector<V3D> const &directions, std::vector<double> const &sigmas, bool specify_size, double peak_radius,
    double back_inner_radius, double back_outer_radius, std::vector<double> &axes_radii, double &inti, double &sigi) {
  // r1, r2 and r3 will give the sizes of the major axis of
  // the peak ellipsoid, and of the inner and outer surface
  // of the background ellipsoidal shell, respectively.
  // They will specify the size as the number of standard
  // deviations in the direction of each of the pricipal
  // axes that the ellipsoid will extend from the center.
  double r1, r2, r3;
 
  double max_sigma = sigmas[0];
  for (int i = 1; i < 3; i++) {
    if (sigmas[i] > max_sigma) {
      max_sigma = sigmas[i];
    }
  }
 
  if (specify_size) {
    r1 = peak_radius / max_sigma;       // scale specified sizes by 1/max_sigma
    r2 = back_inner_radius / max_sigma; // so when multiplied by the individual
    r3 = back_outer_radius / max_sigma; // sigmas in different directions, the
  }                                     // major axis has the specified size
  else {
    r1 = 3;
    r2 = 3;
    r3 = r2 * 1.25992105; // A factor of 2 ^ (1/3) will make the background
    // shell volume equal to the peak region volume.
 
    // if necessary restrict the background ellipsoid
    // to lie within the specified sphere, and adjust
    // the other sizes, proportionally
    if (r3 * max_sigma > m_radius) {
      r3 = m_radius / max_sigma;
      r1 = r3 * 0.79370053f; // This value for r1 and r2 makes the background
      r2 = r1;               // shell volume equal to the peak region volume.
    }
  }
 
  axes_radii.clear();
  std::vector<double> abcBackgroundOuterRadii;
  std::vector<double> abcBackgroundInnerRadii;
  std::vector<double> abcRadii;
  for (int i = 0; i < 3; i++) {
    abcBackgroundOuterRadii.emplace_back(r3 * sigmas[i]);
    abcBackgroundInnerRadii.emplace_back(r2 * sigmas[i]);
    abcRadii.emplace_back(r1 * sigmas[i]);
    axes_radii.emplace_back(r1 * sigmas[i]);
  }
 
  if (!E1Vec.empty()) {
    double h3 = 1.0 - detectorQ(E1Vec, peak_q, abcBackgroundOuterRadii);
    // scaled from area of circle minus segment when r normalized to 1
    double m3 = std::sqrt(1.0 - (std::acos(1.0 - h3) - (1.0 - h3) * std::sqrt(2.0 * h3 - h3 * h3)) / M_PI);
    double h1 = 1.0 - detectorQ(E1Vec, peak_q, axes_radii);
    // Do not use peak if edge of detector is inside integration radius
    if (h1 > 0.0)
      return std::make_shared<const PeakShapeEllipsoid>(directions, abcRadii, abcBackgroundInnerRadii,
                                                        abcBackgroundOuterRadii, Mantid::Kernel::QLab,
                                                        "IntegrateEllipsoids");
    r3 *= m3;
    if (r2 != r1) {
      double h2 = 1.0 - detectorQ(E1Vec, peak_q, abcBackgroundInnerRadii);
      // scaled from area of circle minus segment when r normalized to 1
      double m2 = std::sqrt(1.0 - (std::acos(1.0 - h2) - (1.0 - h2) * std::sqrt(2.0 * h2 - h2 * h2)) / M_PI);
      r2 *= m2;
    }
  }
 
  std::pair<double, double> backgrd = numInEllipsoidBkg(ev_list, directions, abcBackgroundOuterRadii,
                                                        abcBackgroundInnerRadii, m_useOnePercentBackgroundCorrection);
 
  std::pair<double, double> peak_w_back = numInEllipsoid(ev_list, directions, axes_radii);
 
  double ratio = pow(r1, 3) / (pow(r3, 3) - pow(r2, 3));
 
  inti = peak_w_back.first - ratio * backgrd.first;
  sigi = sqrt(peak_w_back.second + ratio * ratio * backgrd.second);
 
  // Make the shape and return it.
  return std::make_shared<const PeakShapeEllipsoid>(directions, abcRadii, abcBackgroundInnerRadii,
                                                    abcBackgroundOuterRadii, Mantid::Kernel::QLab,
                                                    "IntegrateEllipsoids");
}
double Integrate3DEvents::detectorQ(const std::vector<V3D> &E1Vec, const V3D &QLabFrame, const std::vector<double> &r) {
  double quot = 1.0;
  for (auto &E1 : E1Vec) {
    V3D distv = QLabFrame - E1 * (QLabFrame.scalar_prod(E1)); // distance to the trajectory as a vector
    double quot0 = distv.norm() / *(std::min_element(r.begin(), r.end()));
    if (quot0 < quot) {
      quot = quot0;
    }
  }
  return quot;
}
 
std::tuple<double, double, double> Integrate3DEvents::calculateRadiusFactors(const IntegrationParameters &params,
                                                                             double max_sigma) const {
  double r1 = 0, r2 = 0, r3 = 0;
 
  if (!params.specifySize) {
    r1 = 3;
    r2 = 3;
    r3 = r2 * 1.25992105; // A factor of 2 ^ (1/3) will make the background
    // shell volume equal to the peak region volume.
 
    // if necessary restrict the background ellipsoid
    // to lie within the specified sphere, and adjust
    // the other sizes, proportionally
    if (r3 * max_sigma > params.regionRadius) {
      r3 = params.regionRadius / max_sigma;
      r1 = r3 * 0.79370053f; // This value for r1 and r2 makes the background
      r2 = r1;               // shell volume equal to the peak region volume.
    }
  } else {
    // scale specified sizes by 1/max_sigma
    // so when multiplied by the individual
    // sigmas in different directions, the
    r1 = params.peakRadius / max_sigma;
    r2 = params.backgroundInnerRadius / max_sigma;
    r3 = params.backgroundOuterRadius / max_sigma;
  }
 
  return std::make_tuple(r1, r2, r3);
}
 
} // namespace Mantid::MDAlgorithms