MDNormSCD.cpp

#include "MantidMDAlgorithms/MDNormSCD.h"

#include "MantidAPI/CommonBinsValidator.h"
#include "MantidAPI/InstrumentValidator.h"
#include "MantidAPI/WorkspaceUnitValidator.h"
#include "MantidDataObjects/EventWorkspace.h"
#include "MantidDataObjects/MDEventWorkspace.h"
#include "MantidDataObjects/MDHistoWorkspace.h"
#include "MantidGeometry/Instrument.h"
#include "MantidKernel/CompositeValidator.h"
#include "MantidKernel/TimeSeriesProperty.h"
#include "MantidKernel/VectorHelper.h"

namespace Mantid {
namespace MDAlgorithms {

using Mantid::Kernel::Direction;
using Mantid::API::WorkspaceProperty;
using namespace Mantid::DataObjects;
using namespace Mantid::API;
using namespace Mantid::Kernel;

namespace {
// function to  compare two intersections (h,k,l,Momentum) by Momentum
bool compareMomentum(const Mantid::Kernel::VMD &v1,
                     const Mantid::Kernel::VMD &v2) {
  return (v1[3] < v2[3]);
}
}

// Register the algorithm into the AlgorithmFactory
DECLARE_ALGORITHM(MDNormSCD)

//----------------------------------------------------------------------------------------------
/**
 * Constructor
 */
MDNormSCD::MDNormSCD()
    : m_normWS(), m_inputWS(), m_hmin(0.0f), m_hmax(0.0f), m_kmin(0.0f),
      m_kmax(0.0f), m_lmin(0.0f), m_lmax(0.0f), m_hIntegrated(true),
      m_kIntegrated(true), m_lIntegrated(true), m_rubw(3, 3), m_kiMin(0.0),
      m_kiMax(EMPTY_DBL()), m_hIdx(-1), m_kIdx(-1), m_lIdx(-1), m_hX(), m_kX(),
      m_lX(), m_samplePos(), m_beamDir() {}

/// Algorithm's version for identification. @see Algorithm::version
int MDNormSCD::version() const { return 1; }

/// Algorithm's category for identification. @see Algorithm::category
const std::string MDNormSCD::category() const {
  return "MDAlgorithms\\Normalisation";
}

/// Algorithm's summary for use in the GUI and help. @see Algorithm::summary
const std::string MDNormSCD::summary() const {
  return "Calculate normalization for an MDEvent workspace for single crystal "
         "diffraction.";
}

/// Algorithm's name for use in the GUI and help. @see Algorithm::name
const std::string MDNormSCD::name() const { return "MDNormSCD"; }

//----------------------------------------------------------------------------------------------
/**
  * Initialize the algorithm's properties.
  */
void MDNormSCD::init() {
  declareProperty(new WorkspaceProperty<IMDEventWorkspace>("InputWorkspace", "",
                                                           Direction::Input),
                  "An input MDWorkspace.");

  std::string dimChars = getDimensionChars();
  // --------------- Axis-aligned properties
  // ---------------------------------------
  for (size_t i = 0; i < dimChars.size(); i++) {
    std::string dim(" ");
    dim[0] = dimChars[i];
    std::string propName = "AlignedDim" + dim;
    declareProperty(
        new PropertyWithValue<std::string>(propName, "", Direction::Input),
        "Binning parameters for the " + Strings::toString(i) +
            "th dimension.\n"
            "Enter it as a comma-separated list of values with the format: "
            "'name,minimum,maximum,number_of_bins'. Leave blank for NONE.");
  }

  auto fluxValidator = boost::make_shared<CompositeValidator>();
  fluxValidator->add<WorkspaceUnitValidator>("Momentum");
  fluxValidator->add<InstrumentValidator>();
  fluxValidator->add<CommonBinsValidator>();
  auto solidAngleValidator = fluxValidator->clone();

  declareProperty(new WorkspaceProperty<>("FluxWorkspace", "", Direction::Input,
                                          fluxValidator),
                  "An input workspace containing momentum dependent flux.");
  declareProperty(new WorkspaceProperty<>("SolidAngleWorkspace", "",
                                          Direction::Input,
                                          solidAngleValidator),
                  "An input workspace containing momentum integrated vanadium "
                  "(a measure of the solid angle).");

  declareProperty(new WorkspaceProperty<Workspace>("OutputWorkspace", "",
                                                   Direction::Output),
                  "A name for the output data MDHistoWorkspace.");
  declareProperty(new WorkspaceProperty<Workspace>(
                      "OutputNormalizationWorkspace", "", Direction::Output),
                  "A name for the output normalization MDHistoWorkspace.");
}

//----------------------------------------------------------------------------------------------
/**
 * Execute the algorithm.
 */
void MDNormSCD::exec() {
  cacheInputs();
  auto outputWS = binInputWS();
  convention = Kernel::ConfigService::Instance().getString("Q.convention");
  setProperty<Workspace_sptr>("OutputWorkspace", outputWS);
  createNormalizationWS(*outputWS);
  setProperty("OutputNormalizationWorkspace", m_normWS);

  // Check for other dimensions if we could measure anything in the original
  // data
  bool skipNormalization = false;
  const std::vector<coord_t> otherValues =
      getValuesFromOtherDimensions(skipNormalization);
  const auto affineTrans =
      findIntergratedDimensions(otherValues, skipNormalization);
  cacheDimensionXValues();

  if (!skipNormalization) {
    calculateNormalization(otherValues, affineTrans);
  } else {
    g_log.warning("Binning limits are outside the limits of the MDWorkspace. "
                  "Not applying normalization.");
  }
}

/**
 * Set up starting values for cached variables
 */
void MDNormSCD::cacheInputs() {
  m_inputWS = getProperty("InputWorkspace");
  if (inputEnergyMode() != "Elastic") {
    throw std::invalid_argument("Invalid energy transfer mode. Algorithm "
                                "currently only supports elastic data.");
  }
  // Min/max dimension values
  const auto hdim(m_inputWS->getDimension(0)), kdim(m_inputWS->getDimension(1)),
      ldim(m_inputWS->getDimension(2));
  m_hmin = hdim->getMinimum();
  m_kmin = kdim->getMinimum();
  m_lmin = ldim->getMinimum();
  m_hmax = hdim->getMaximum();
  m_kmax = kdim->getMaximum();
  m_lmax = ldim->getMaximum();

  const auto &exptInfoZero = *(m_inputWS->getExperimentInfo(0));
  auto source = exptInfoZero.getInstrument()->getSource();
  auto sample = exptInfoZero.getInstrument()->getSample();
  if (source == nullptr || sample == nullptr) {
    throw Kernel::Exception::InstrumentDefinitionError(
        "Instrument not sufficiently defined: failed to get source and/or "
        "sample");
  }
  m_samplePos = sample->getPos();
  m_beamDir = m_samplePos - source->getPos();
  m_beamDir.normalize();
}

/**
 * Currently looks for the ConvertToMD algorithm in the history
 * @return A string donating the energy transfer mode of the input workspace
 */
std::string MDNormSCD::inputEnergyMode() const {
  const auto &hist = m_inputWS->getHistory();
  const size_t nalgs = hist.size();
  const auto &lastAlgorithm = hist.lastAlgorithm();

  std::string emode("");
  if (lastAlgorithm->name() == "ConvertToMD") {
    emode = lastAlgorithm->getPropertyValue("dEAnalysisMode");
  } else if ((lastAlgorithm->name() == "Load" ||
              hist.lastAlgorithm()->name() == "LoadMD") &&
             hist.getAlgorithmHistory(nalgs - 2)->name() == "ConvertToMD") {
    // get dEAnalysisMode
    PropertyHistories histvec =
        hist.getAlgorithmHistory(nalgs - 2)->getProperties();
    for (auto &hist : histvec) {
      if (hist->name() == "dEAnalysisMode") {
        emode = hist->value();
        break;
      }
    }
  } else {
    throw std::invalid_argument("The last algorithm in the history of the "
                                "input workspace is not ConvertToMD");
  }
  return emode;
}

/**
 * Runs the BinMD algorithm on the input to provide the output workspace
 * All slicing algorithm properties are passed along
 * @return MDHistoWorkspace as a result of the binning
 */
MDHistoWorkspace_sptr MDNormSCD::binInputWS() {
  const auto &props = getProperties();
  IAlgorithm_sptr binMD = createChildAlgorithm("BinMD", 0.0, 0.3);
  binMD->setPropertyValue("AxisAligned", "1");
  for (auto prop : props) {
    const auto &propName = prop->name();
    if (propName != "FluxWorkspace" && propName != "SolidAngleWorkspace" &&
        propName != "OutputNormalizationWorkspace") {
      binMD->setPropertyValue(propName, prop->value());
    }
  }
  binMD->executeAsChildAlg();
  Workspace_sptr outputWS = binMD->getProperty("OutputWorkspace");
  return boost::dynamic_pointer_cast<MDHistoWorkspace>(outputWS);
}

/**
 * Create & cached the normalization workspace
 * @param dataWS The binned workspace that will be used for the data
 */
void MDNormSCD::createNormalizationWS(const MDHistoWorkspace &dataWS) {
  // Copy the MDHisto workspace, and change signals and errors to 0.
  m_normWS.reset(dataWS.clone().release());
  m_normWS->setTo(0., 0., 0.);
}

/**
 * Retrieve logged values from non-HKL dimensions
 * @param skipNormalization [InOut] Updated to false if any values are outside
 * range measured by input workspace
 * @return A vector of values from other dimensions to be include in normalized
 * MD position calculation
 */
std::vector<coord_t>
MDNormSCD::getValuesFromOtherDimensions(bool &skipNormalization) const {
  const auto &runZero = m_inputWS->getExperimentInfo(0)->run();

  std::vector<coord_t> otherDimValues;
  for (size_t i = 3; i < m_inputWS->getNumDims(); i++) {
    const auto dimension = m_inputWS->getDimension(i);
    float dimMin = static_cast<float>(dimension->getMinimum());
    float dimMax = static_cast<float>(dimension->getMaximum());
    auto *dimProp = dynamic_cast<Kernel::TimeSeriesProperty<double> *>(
        runZero.getProperty(dimension->getName()));
    if (dimProp) {
      coord_t value = static_cast<coord_t>(dimProp->firstValue());
      otherDimValues.push_back(value);
      // in the original MD data no time was spent measuring between dimMin and
      // dimMax
      if (value < dimMin || value > dimMax) {
        skipNormalization = true;
      }
    }
  }
  return otherDimValues;
}

/**
 * Checks the normalization workspace against the indices of the original
 * dimensions.
 * If not found, the corresponding dimension is integrated
 * @param otherDimValues Values from non-HKL dimensions
 * @param skipNormalization [InOut] Sets the flag true if normalization values
 * are outside of original inputs
 * @return Affine trasform matrix
 */
Kernel::Matrix<coord_t>
MDNormSCD::findIntergratedDimensions(const std::vector<coord_t> &otherDimValues,
                                     bool &skipNormalization) {
  // Get indices of the original dimensions in the output workspace,
  // and if not found, the corresponding dimension is integrated
  Kernel::Matrix<coord_t> affineMat =
      m_normWS->getTransformFromOriginal(0)->makeAffineMatrix();

  const size_t nrm1 = affineMat.numRows() - 1;
  const size_t ncm1 = affineMat.numCols() - 1;
  for (size_t row = 0; row < nrm1; row++) // affine matrix, ignore last row
  {
    const auto dimen = m_normWS->getDimension(row);
    const auto dimMin(dimen->getMinimum()), dimMax(dimen->getMaximum());
    if (affineMat[row][0] == 1.) {
      m_hIntegrated = false;
      m_hIdx = row;
      if (m_hmin < dimMin)
        m_hmin = dimMin;
      if (m_hmax > dimMax)
        m_hmax = dimMax;
      if (m_hmin > dimMax || m_hmax < dimMin) {
        skipNormalization = true;
      }
    }
    if (affineMat[row][1] == 1.) {
      m_kIntegrated = false;
      m_kIdx = row;
      if (m_kmin < dimMin)
        m_kmin = dimMin;
      if (m_kmax > dimMax)
        m_kmax = dimMax;
      if (m_kmin > dimMax || m_kmax < dimMin) {
        skipNormalization = true;
      }
    }
    if (affineMat[row][2] == 1.) {
      m_lIntegrated = false;
      m_lIdx = row;
      if (m_lmin < dimMin)
        m_lmin = dimMin;
      if (m_lmax > dimMax)
        m_lmax = dimMax;
      if (m_lmin > dimMax || m_lmax < dimMin) {
        skipNormalization = true;
      }
    }

    for (size_t col = 3; col < ncm1; col++) // affine matrix, ignore last column
    {
      if (affineMat[row][col] == 1.) {
        double val = otherDimValues.at(col - 3);
        if (val > dimMax || val < dimMin) {
          skipNormalization = true;
        }
      }
    }
  }

  return affineMat;
}

/**
 * Stores the X values from each H,K,L dimension as member variables
 */
void MDNormSCD::cacheDimensionXValues() {
  if (!m_hIntegrated) {
    auto &hDim = *m_normWS->getDimension(m_hIdx);
    m_hX.resize(hDim.getNBins());
    for (size_t i = 0; i < m_hX.size(); ++i) {
      m_hX[i] = hDim.getX(i);
    }
  }
  if (!m_kIntegrated) {
    auto &kDim = *m_normWS->getDimension(m_kIdx);
    m_kX.resize(kDim.getNBins());
    for (size_t i = 0; i < m_kX.size(); ++i) {
      m_kX[i] = kDim.getX(i);
    }
  }
  if (!m_lIntegrated) {
    auto &lDim = *m_normWS->getDimension(m_lIdx);
    m_lX.resize(lDim.getNBins());
    for (size_t i = 0; i < m_lX.size(); ++i) {
      m_lX[i] = lDim.getX(i);
    }
  }
}

/**
 * Computed the normalization for the input workspace. Results are stored in
 * m_normWS
 * @param otherValues
 * @param affineTrans
 */
void MDNormSCD::calculateNormalization(
    const std::vector<coord_t> &otherValues,
    const Kernel::Matrix<coord_t> &affineTrans) {
  API::MatrixWorkspace_const_sptr integrFlux = getProperty("FluxWorkspace");
  integrFlux->getXMinMax(m_kiMin, m_kiMax);
  API::MatrixWorkspace_const_sptr solidAngleWS =
      getProperty("SolidAngleWorkspace");

  const auto &exptInfoZero = *(m_inputWS->getExperimentInfo(0));
  typedef Kernel::PropertyWithValue<std::vector<double>> VectorDoubleProperty;
  auto *rubwLog =
      dynamic_cast<VectorDoubleProperty *>(exptInfoZero.getLog("RUBW_MATRIX"));
  if (!rubwLog) {
    throw std::runtime_error(
        "Wokspace does not contain a log entry for the RUBW matrix."
        "Cannot continue.");
  } else {
    Kernel::DblMatrix rubwValue(
        (*rubwLog)()); // includes the 2*pi factor but not goniometer for now :)
    m_rubw = exptInfoZero.run().getGoniometerMatrix() * rubwValue;
    m_rubw.Invert();
  }
  const double protonCharge = exptInfoZero.run().getProtonCharge();

  auto instrument = exptInfoZero.getInstrument();
  std::vector<detid_t> detIDs = instrument->getDetectorIDs(true);
  // Prune out those that are part of a group and simply leave the head of the
  // group
  detIDs = removeGroupedIDs(exptInfoZero, detIDs);

  // Mappings
  const int64_t ndets = static_cast<int64_t>(detIDs.size());
  const detid2index_map fluxDetToIdx =
      integrFlux->getDetectorIDToWorkspaceIndexMap();
  const detid2index_map solidAngDetToIdx =
      solidAngleWS->getDetectorIDToWorkspaceIndexMap();

  std::unique_ptr<API::Progress> prog(new API::Progress(this, 0.3, 1.0, ndets));
  PARALLEL_FOR1(integrFlux)
  for (int64_t i = 0; i < ndets; i++) {
    PARALLEL_START_INTERUPT_REGION

    const auto detID = detIDs[i];
    double theta(0.0), phi(0.0);
    bool skip(false);
    try {
      auto spectrum = getThetaPhi(detID, exptInfoZero, theta, phi);
      if (spectrum->isMonitor() || spectrum->isMasked())
        continue;
    } catch (
        std::exception &) // detector might not exist or has no been included
                          // in grouping
    {
      skip = true; // Intel compiler has a problem with continue inside a catch
                   // inside openmp...
    }
    if (skip)
      continue;

    // Intersections
    auto intersections = calculateIntersections(theta, phi);
    if (intersections.empty())
      continue;

    // get the flux spetrum number
    size_t wsIdx = fluxDetToIdx.find(detID)->second;
    // Get solid angle for this contribution
    double solid =
        solidAngleWS->readY(solidAngDetToIdx.find(detID)->second)[0] *
        protonCharge;

    // -- calculate integrals for the intersection --
    // momentum values at intersections
    auto intersectionsBegin = intersections.begin();
    std::vector<double> xValues(intersections.size()),
        yValues(intersections.size());
    {
      // copy momenta to xValues
      auto x = xValues.begin();
      for (auto it = intersectionsBegin; it != intersections.end(); ++it, ++x) {
        *x = (*it)[3];
      }
    }
    // calculate integrals at momenta from xValues by interpolating between
    // points in spectrum sp
    // of workspace integrFlux. The result is stored in yValues
    calcIntegralsForIntersections(xValues, *integrFlux, wsIdx, yValues);

    // Compute final position in HKL
    const size_t vmdDims = intersections.front().size();
    // pre-allocate for efficiency and copy non-hkl dim values into place
    std::vector<coord_t> pos(vmdDims + otherValues.size());
    std::copy(otherValues.begin(), otherValues.end(),
              pos.begin() + vmdDims - 1);
    pos.push_back(1.);

    for (auto it = intersectionsBegin + 1; it != intersections.end(); ++it) {
      const auto &curIntSec = *it;
      const auto &prevIntSec = *(it - 1);
      // the full vector isn't used so compute only what is necessary
      double delta = curIntSec[3] - prevIntSec[3];
      if (delta < 1e-07)
        continue; // Assume zero contribution if difference is small

      // Average between two intersections for final position
      std::transform(curIntSec.getBareArray(),
                     curIntSec.getBareArray() + vmdDims - 1,
                     prevIntSec.getBareArray(), pos.begin(),
                     VectorHelper::SimpleAverage<coord_t>());
      std::vector<coord_t> posNew = affineTrans * pos;
      size_t linIndex = m_normWS->getLinearIndexAtCoord(posNew.data());
      if (linIndex == size_t(-1))
        continue;

      // index of the current intersection
      size_t k = static_cast<size_t>(std::distance(intersectionsBegin, it));
      // signal = integral between two consecutive intersections
      double signal = (yValues[k] - yValues[k - 1]) * solid;

      PARALLEL_CRITICAL(updateMD) {
        signal += m_normWS->getSignalAt(linIndex);
        m_normWS->setSignalAt(linIndex, signal);
      }
    }
    prog->report();

    PARALLEL_END_INTERUPT_REGION
  }
  PARALLEL_CHECK_INTERUPT_REGION
}

/**
 * Linearly interpolate between the points in integrFlux at xValues and save the
 * results in yValues.
 * @param xValues :: X-values at which to interpolate
 * @param integrFlux :: A workspace with the spectra to interpolate
 * @param sp :: A workspace index for a spectrum in integrFlux to interpolate.
 * @param yValues :: A vector to save the results.
 */
void MDNormSCD::calcIntegralsForIntersections(
    const std::vector<double> &xValues, const API::MatrixWorkspace &integrFlux,
    size_t sp, std::vector<double> &yValues) const {
  assert(xValues.size() == yValues.size());

  // the x-data from the workspace
  const auto &xData = integrFlux.readX(sp);
  const double xStart = xData.front();
  const double xEnd = xData.back();

  // the values in integrFlux are expected to be integrals of a non-negative
  // function
  // ie they must make a non-decreasing function
  const auto &yData = integrFlux.readY(sp);
  size_t spSize = yData.size();

  const double yMin = 0.0;
  const double yMax = yData.back();

  size_t nData = xValues.size();
  // all integrals below xStart must be 0
  if (xValues[nData - 1] < xStart) {
    std::fill(yValues.begin(), yValues.end(), yMin);
    return;
  }

  // all integrals above xEnd must be equal tp yMax
  if (xValues[0] > xEnd) {
    std::fill(yValues.begin(), yValues.end(), yMax);
    return;
  }

  size_t i = 0;
  // integrals below xStart must be 0
  while (i < nData - 1 && xValues[i] < xStart) {
    yValues[i] = yMin;
    i++;
  }
  size_t j = 0;
  for (; i < nData; i++) {
    // integrals above xEnd must be equal tp yMax
    if (j >= spSize - 1) {
      yValues[i] = yMax;
    } else {
      double xi = xValues[i];
      while (j < spSize - 1 && xi > xData[j])
        j++;
      // if x falls onto an interpolation point return the corresponding y
      if (xi == xData[j]) {
        yValues[i] = yData[j];
      } else if (j == spSize - 1) {
        // if we get above xEnd it's yMax
        yValues[i] = yMax;
      } else if (j > 0) {
        // interpolate between the consecutive points
        double x0 = xData[j - 1];
        double x1 = xData[j];
        double y0 = yData[j - 1];
        double y1 = yData[j];
        yValues[i] = y0 + (y1 - y0) * (xi - x0) / (x1 - x0);
      } else // j == 0
      {
        yValues[i] = yMin;
      }
    }
  }
}

/**
 * Checks for IDs that are actually part of the same group and just keeps one
 * from the group.
 * For a 1:1 map, none will be removed.
 * @param exptInfo An ExperimentInfo object that defines the grouping
 * @param detIDs A list of existing IDs
 * @return A new list of IDs
 */
std::vector<detid_t>
MDNormSCD::removeGroupedIDs(const ExperimentInfo &exptInfo,
                            const std::vector<detid_t> &detIDs) {
  const size_t ntotal = detIDs.size();
  std::vector<detid_t> singleIDs;
  singleIDs.reserve(ntotal / 2); // reserve half. In the case of 1:1 it will
                                 // double to the correct size once
  std::set<detid_t> groupedIDs;

  for (auto curID : detIDs) {
    if (groupedIDs.count(curID) == 1)
      continue; // Already been processed

    try {
      const auto &members = exptInfo.getGroupMembers(curID);
      singleIDs.push_back(members.front());
      std::copy(members.begin() + 1, members.end(),
                std::inserter(groupedIDs, groupedIDs.begin()));
    } catch (std::runtime_error &) {
      singleIDs.push_back(curID);
    }
  }

  g_log.debug() << "Found " << singleIDs.size() << " spectra from  "
                << detIDs.size() << " IDs\n";
  return singleIDs;
}

/**
 * Get the theta and phi angles for the given ID. If the detector was part of a
 * group,
 * as defined in the ExperimentInfo object, then the theta/phi are for the whole
 * set.
 * @param detID A reference to a single ID
 * @param exptInfo A reference to the ExperimentInfo that defines that
 * spectrum->detector mapping
 * @param theta [Output] Set to the theta angle for the detector (set)
 * @param phi [Output] Set to the phi angle for the detector (set)
 * @return A poiner to the Detector object for this spectrum as a whole
 *         (may be a single pixel or group)
 */
Geometry::IDetector_const_sptr
MDNormSCD::getThetaPhi(const detid_t detID, const ExperimentInfo &exptInfo,
                       double &theta, double &phi) {
  const auto spectrum = exptInfo.getDetectorByID(detID);
  theta = spectrum->getTwoTheta(m_samplePos, m_beamDir);
  phi = spectrum->getPhi();
  return spectrum;
}

/**
 * Calculate the points of intersection for the given detector with cuboid
 * surrounding the
 * detector position in HKL
 * @param theta Polar angle withd detector
 * @param phi Azimuthal angle with detector
 * @return A list of intersections in HKL space
 */
std::vector<Kernel::VMD> MDNormSCD::calculateIntersections(const double theta,
                                                           const double phi) {
  V3D q(-sin(theta) * cos(phi), -sin(theta) * sin(phi), 1. - cos(theta));
  q = m_rubw * q;
  if (convention == "Crystallography") {
    q *= -1;
  }

  double hStart = q.X() * m_kiMin, hEnd = q.X() * m_kiMax;
  double kStart = q.Y() * m_kiMin, kEnd = q.Y() * m_kiMax;
  double lStart = q.Z() * m_kiMin, lEnd = q.Z() * m_kiMax;

  double eps = 1e-7;

  auto hNBins = m_hX.size();
  auto kNBins = m_kX.size();
  auto lNBins = m_lX.size();
  std::vector<Kernel::VMD> intersections;
  intersections.reserve(hNBins + kNBins + lNBins + 8);

  // calculate intersections with planes perpendicular to h
  if (fabs(hStart - hEnd) > eps) {
    double fmom = (m_kiMax - m_kiMin) / (hEnd - hStart);
    double fk = (kEnd - kStart) / (hEnd - hStart);
    double fl = (lEnd - lStart) / (hEnd - hStart);
    if (!m_hIntegrated) {
      for (size_t i = 0; i < hNBins; i++) {
        double hi = m_hX[i];
        if ((hi >= m_hmin) && (hi <= m_hmax) &&
            ((hStart - hi) * (hEnd - hi) < 0)) {
          // if hi is between hStart and hEnd, then ki and li will be between
          // kStart, kEnd and lStart, lEnd and momi will be between m_kiMin and
          // KnincidemtmMax
          double ki = fk * (hi - hStart) + kStart;
          double li = fl * (hi - hStart) + lStart;
          if ((ki >= m_kmin) && (ki <= m_kmax) && (li >= m_lmin) &&
              (li <= m_lmax)) {
            double momi = fmom * (hi - hStart) + m_kiMin;
            Mantid::Kernel::VMD v(hi, ki, li, momi);
            intersections.push_back(v);
          }
        }
      }
    }

    double momhMin = fmom * (m_hmin - hStart) + m_kiMin;
    if ((momhMin > m_kiMin) && (momhMin < m_kiMax)) {
      // khmin and lhmin
      double khmin = fk * (m_hmin - hStart) + kStart;
      double lhmin = fl * (m_hmin - hStart) + lStart;
      if ((khmin >= m_kmin) && (khmin <= m_kmax) && (lhmin >= m_lmin) &&
          (lhmin <= m_lmax)) {
        Mantid::Kernel::VMD v(m_hmin, khmin, lhmin, momhMin);
        intersections.push_back(v);
      }
    }
    double momhMax = fmom * (m_hmax - hStart) + m_kiMin;
    if ((momhMax > m_kiMin) && (momhMax < m_kiMax)) {
      // khmax and lhmax
      double khmax = fk * (m_hmax - hStart) + kStart;
      double lhmax = fl * (m_hmax - hStart) + lStart;
      if ((khmax >= m_kmin) && (khmax <= m_kmax) && (lhmax >= m_lmin) &&
          (lhmax <= m_lmax)) {
        Mantid::Kernel::VMD v(m_hmax, khmax, lhmax, momhMax);
        intersections.push_back(v);
      }
    }
  }

  // calculate intersections with planes perpendicular to k
  if (fabs(kStart - kEnd) > eps) {
    double fmom = (m_kiMax - m_kiMin) / (kEnd - kStart);
    double fh = (hEnd - hStart) / (kEnd - kStart);
    double fl = (lEnd - lStart) / (kEnd - kStart);
    if (!m_kIntegrated) {
      for (size_t i = 0; i < kNBins; i++) {
        double ki = m_kX[i];
        if ((ki >= m_kmin) && (ki <= m_kmax) &&
            ((kStart - ki) * (kEnd - ki) < 0)) {
          // if ki is between kStart and kEnd, then hi and li will be between
          // hStart, hEnd and lStart, lEnd
          double hi = fh * (ki - kStart) + hStart;
          double li = fl * (ki - kStart) + lStart;
          if ((hi >= m_hmin) && (hi <= m_hmax) && (li >= m_lmin) &&
              (li <= m_lmax)) {
            double momi = fmom * (ki - kStart) + m_kiMin;
            Mantid::Kernel::VMD v(hi, ki, li, momi);
            intersections.push_back(v);
          }
        }
      }
    }

    double momkMin = fmom * (m_kmin - kStart) + m_kiMin;
    if ((momkMin > m_kiMin) && (momkMin < m_kiMax)) {
      // hkmin and lkmin
      double hkmin = fh * (m_kmin - kStart) + hStart;
      double lkmin = fl * (m_kmin - kStart) + lStart;
      if ((hkmin >= m_hmin) && (hkmin <= m_hmax) && (lkmin >= m_lmin) &&
          (lkmin <= m_lmax)) {
        Mantid::Kernel::VMD v(hkmin, m_kmin, lkmin, momkMin);
        intersections.push_back(v);
      }
    }
    double momkMax = fmom * (m_kmax - kStart) + m_kiMin;
    if ((momkMax > m_kiMin) && (momkMax < m_kiMax)) {
      // hkmax and lkmax
      double hkmax = fh * (m_kmax - kStart) + hStart;
      double lkmax = fl * (m_kmax - kStart) + lStart;
      if ((hkmax >= m_hmin) && (hkmax <= m_hmax) && (lkmax >= m_lmin) &&
          (lkmax <= m_lmax)) {
        Mantid::Kernel::VMD v(hkmax, m_kmax, lkmax, momkMax);
        intersections.push_back(v);
      }
    }
  }

  // calculate intersections with planes perpendicular to l
  if (fabs(lStart - lEnd) > eps) {
    double fmom = (m_kiMax - m_kiMin) / (lEnd - lStart);
    double fh = (hEnd - hStart) / (lEnd - lStart);
    double fk = (kEnd - kStart) / (lEnd - lStart);
    if (!m_lIntegrated) {
      for (size_t i = 0; i < lNBins; i++) {
        double li = m_lX[i];
        if ((li >= m_lmin) && (li <= m_lmax) &&
            ((lStart - li) * (lEnd - li) < 0)) {
          // if li is between lStart and lEnd, then hi and ki will be between
          // hStart, hEnd and kStart, kEnd
          double hi = fh * (li - lStart) + hStart;
          double ki = fk * (li - lStart) + kStart;
          if ((hi >= m_hmin) && (hi <= m_hmax) && (ki >= m_kmin) &&
              (ki <= m_kmax)) {
            double momi = fmom * (li - lStart) + m_kiMin;
            Mantid::Kernel::VMD v(hi, ki, li, momi);
            intersections.push_back(v);
          }
        }
      }
    }

    double momlMin = fmom * (m_lmin - lStart) + m_kiMin;
    if ((momlMin > m_kiMin) && (momlMin < m_kiMax)) {
      // hlmin and klmin
      double hlmin = fh * (m_lmin - lStart) + hStart;
      double klmin = fk * (m_lmin - lStart) + kStart;
      if ((hlmin >= m_hmin) && (hlmin <= m_hmax) && (klmin >= m_kmin) &&
          (klmin <= m_kmax)) {
        Mantid::Kernel::VMD v(hlmin, klmin, m_lmin, momlMin);
        intersections.push_back(v);
      }
    }
    double momlMax = fmom * (m_lmax - lStart) + m_kiMin;
    if ((momlMax > m_kiMin) && (momlMax < m_kiMax)) {
      // khmax and lhmax
      double hlmax = fh * (m_lmax - lStart) + hStart;
      double klmax = fk * (m_lmax - lStart) + kStart;
      if ((hlmax >= m_hmin) && (hlmax <= m_hmax) && (klmax >= m_kmin) &&
          (klmax <= m_kmax)) {
        Mantid::Kernel::VMD v(hlmax, klmax, m_lmax, momlMax);
        intersections.push_back(v);
      }
    }
  }

  // add endpoints
  if ((hStart >= m_hmin) && (hStart <= m_hmax) && (kStart >= m_kmin) &&
      (kStart <= m_kmax) && (lStart >= m_lmin) && (lStart <= m_lmax)) {
    Mantid::Kernel::VMD v(hStart, kStart, lStart, m_kiMin);
    intersections.push_back(v);
  }
  if ((hEnd >= m_hmin) && (hEnd <= m_hmax) && (kEnd >= m_kmin) &&
      (kEnd <= m_kmax) && (lEnd >= m_lmin) && (lEnd <= m_lmax)) {
    Mantid::Kernel::VMD v(hEnd, kEnd, lEnd, m_kiMax);
    intersections.push_back(v);
  }

  // sort intersections by momentum
  typedef std::vector<Mantid::Kernel::VMD>::iterator IterType;
  std::stable_sort<IterType, bool (*)(const Mantid::Kernel::VMD &,
                                      const Mantid::Kernel::VMD &)>(
      intersections.begin(), intersections.end(), compareMomentum);

  return intersections;
}

} // namespace MDAlgorithms
} // namespace Mantid