MDNorm.cpp

// 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
// SPDX - License - Identifier: GPL - 3.0 +

#include "MantidMDAlgorithms/MDNorm.h"
#include "MantidAPI/CommonBinsValidator.h"
#include "MantidAPI/IMDEventWorkspace.h"
#include "MantidAPI/InstrumentValidator.h"
#include "MantidAPI/Run.h"
#include "MantidAPI/Sample.h"
#include "MantidAPI/SpectrumInfo.h"
#include "MantidDataObjects/MDHistoWorkspace.h"
#include "MantidGeometry/Crystal/OrientedLattice.h"
#include "MantidGeometry/Crystal/PointGroupFactory.h"
#include "MantidGeometry/Crystal/SpaceGroupFactory.h"
#include "MantidGeometry/Crystal/SymmetryOperationFactory.h"
#include "MantidGeometry/Instrument.h"
#include "MantidGeometry/MDGeometry/HKL.h"
#include "MantidGeometry/MDGeometry/MDFrameFactory.h"
#include "MantidGeometry/MDGeometry/QSample.h"
#include "MantidKernel/ArrayLengthValidator.h"
#include "MantidKernel/ArrayProperty.h"
#include "MantidKernel/CompositeValidator.h"
#include "MantidKernel/ConfigService.h"
#include "MantidKernel/Exception.h"
#include "MantidKernel/Strings.h"
#include "MantidKernel/VisibleWhenProperty.h"
#include <boost/lexical_cast.hpp>

namespace Mantid {
namespace MDAlgorithms {

using namespace Mantid::Kernel;
using namespace Mantid::API;
using namespace Mantid::Geometry;
using namespace Mantid::DataObjects;

namespace {
using VectorDoubleProperty = Kernel::PropertyWithValue<std::vector<double>>;
// function to  compare two intersections (h,k,l,Momentum) by Momentum
bool compareMomentum(const std::array<double, 4> &v1,
                     const std::array<double, 4> &v2) {
  return (v1[3] < v2[3]);
}

// k=sqrt(energyToK * E)
constexpr double energyToK = 8.0 * M_PI * M_PI *
                             PhysicalConstants::NeutronMass *
                             PhysicalConstants::meV * 1e-20 /
                             (PhysicalConstants::h * PhysicalConstants::h);

// compare absolute values of integers
static bool abs_compare(int a, int b) { return (std::abs(a) < std::abs(b)); }
} // namespace

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

//----------------------------------------------------------------------------------------------
/**
 * Constructor
 */
MDNorm::MDNorm()
    : m_normWS(), m_inputWS(), m_isRLU(false), m_UB(3, 3, true),
      m_W(3, 3, true), m_transformation(), m_hX(), m_kX(), m_lX(), m_eX(),
      m_hIdx(-1), m_kIdx(-1), m_lIdx(-1), m_eIdx(-1), m_numExptInfos(0),
      m_Ei(0.0), m_diffraction(true), m_accumulate(false),
      m_dEIntegrated(false), m_samplePos(), m_beamDir(), convention("") {}

/// Algorithms name for identification. @see Algorithm::name
const std::string MDNorm::name() const { return "MDNorm"; }

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

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

/// Algorithm's summary for use in the GUI and help. @see Algorithm::summary
const std::string MDNorm::summary() const {
  return "Bins multidimensional data and calculate the normalization on the "
         "same grid";
}

//----------------------------------------------------------------------------------------------
/** Initialize the algorithm's properties.
 */
void MDNorm::init() {
  declareProperty(make_unique<WorkspaceProperty<API::IMDEventWorkspace>>(
                      "InputWorkspace", "", Kernel::Direction::Input),
                  "An input MDEventWorkspace. Must be in Q_sample frame.");

  // RLU and settings
  declareProperty("RLU", true,
                  "Use reciprocal lattice units. If false, use Q_sample");
  setPropertyGroup("RLU", "Q projections RLU");

  auto mustBe3D = boost::make_shared<Kernel::ArrayLengthValidator<double>>(3);
  std::vector<double> Q0(3, 0.), Q1(3, 0), Q2(3, 0);
  Q0[0] = 1.;
  Q1[1] = 1.;
  Q2[2] = 1.;

  declareProperty(
      make_unique<ArrayProperty<double>>("QDimension0", Q0, mustBe3D),
      "The first Q projection axis - Default is (1,0,0)");
  setPropertySettings("QDimension0", make_unique<Kernel::VisibleWhenProperty>(
                                         "RLU", IS_EQUAL_TO, "1"));
  setPropertyGroup("QDimension0", "Q projections RLU");

  declareProperty(
      make_unique<ArrayProperty<double>>("QDimension1", Q1, mustBe3D),
      "The second Q projection axis - Default is (0,1,0)");
  setPropertySettings("QDimension1", make_unique<Kernel::VisibleWhenProperty>(
                                         "RLU", IS_EQUAL_TO, "1"));
  setPropertyGroup("QDimension1", "Q projections RLU");

  declareProperty(
      make_unique<ArrayProperty<double>>("QDimension2", Q2, mustBe3D),
      "The thirdtCalculateCover Q projection axis - Default is (0,0,1)");
  setPropertySettings("QDimension2", make_unique<Kernel::VisibleWhenProperty>(
                                         "RLU", IS_EQUAL_TO, "1"));
  setPropertyGroup("QDimension2", "Q projections RLU");

  // vanadium
  auto fluxValidator = boost::make_shared<CompositeValidator>();
  fluxValidator->add<InstrumentValidator>();
  fluxValidator->add<CommonBinsValidator>();
  auto solidAngleValidator = fluxValidator->clone();
  declareProperty(
      make_unique<WorkspaceProperty<>>(
          "SolidAngleWorkspace", "", Direction::Input,
          API::PropertyMode::Optional, solidAngleValidator),
      "An input workspace containing integrated vanadium "
      "(a measure of the solid angle).\n"
      "Mandatory for diffraction, optional for direct geometry inelastic");
  declareProperty(
      make_unique<WorkspaceProperty<>>("FluxWorkspace", "", Direction::Input,
                                       API::PropertyMode::Optional,
                                       fluxValidator),
      "An input workspace containing momentum dependent flux.\n"
      "Mandatory for diffraction. No effect on direct geometry inelastic");
  setPropertyGroup("SolidAngleWorkspace", "Vanadium normalization");
  setPropertyGroup("FluxWorkspace", "Vanadium normalization");

  // Define slicing
  for (std::size_t i = 0; i < 6; i++) {
    std::string propName = "Dimension" + Strings::toString(i) + "Name";
    std::string propBinning = "Dimension" + Strings::toString(i) + "Binning";
    declareProperty(Kernel::make_unique<PropertyWithValue<std::string>>(
                        propName, "", Direction::Input),
                    "Name for the " + Strings::toString(i) +
                        "th dimension. Leave blank for NONE.");
    auto atMost3 = boost::make_shared<ArrayLengthValidator<double>>(0, 3);
    std::vector<double> temp;
    declareProperty(
        Kernel::make_unique<ArrayProperty<double>>(propBinning, temp, atMost3),
        "Binning for the " + Strings::toString(i) + "th dimension.\n" +
            "- Leave blank for complete integration\n" +
            "- One value is interpreted as step\n"
            "- Two values are interpreted integration interval\n" +
            "- Three values are interpreted as min, step, max");
    setPropertyGroup(propName, "Binning");
    setPropertyGroup(propBinning, "Binning");
  }

  // symmetry operations
  declareProperty(Kernel::make_unique<PropertyWithValue<std::string>>(
                      "SymmetryOperations", "", Direction::Input),
                  "If specified the symmetry will be applied, "
                  "can be space group name, point group name, or list "
                  "individual symmetries.");

  // temporary workspaces
  declareProperty(make_unique<WorkspaceProperty<IMDHistoWorkspace>>(
                      "TemporaryDataWorkspace", "", Direction::Input,
                      PropertyMode::Optional),
                  "An input MDHistoWorkspace used to accumulate data from "
                  "multiple MDEventWorkspaces. If unspecified a blank "
                  "MDHistoWorkspace will be created.");
  declareProperty(make_unique<WorkspaceProperty<IMDHistoWorkspace>>(
                      "TemporaryNormalizationWorkspace", "", Direction::Input,
                      PropertyMode::Optional),
                  "An input MDHistoWorkspace used to accumulate normalization "
                  "from multiple MDEventWorkspaces. If unspecified a blank "
                  "MDHistoWorkspace will be created.");
  setPropertyGroup("TemporaryDataWorkspace", "Temporary workspaces");
  setPropertyGroup("TemporaryNormalizationWorkspace", "Temporary workspaces");

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

//----------------------------------------------------------------------------------------------
/// Validate the input workspace @see Algorithm::validateInputs
std::map<std::string, std::string> MDNorm::validateInputs() {
  std::map<std::string, std::string> errorMessage;

  // Check for input workspace frame
  Mantid::API::IMDEventWorkspace_sptr inputWS =
      this->getProperty("InputWorkspace");
  if (inputWS->getNumDims() < 3) {
    errorMessage.emplace("InputWorkspace",
                         "The input workspace must be at least 3D");
  } else {
    for (size_t i = 0; i < 3; i++) {
      if (inputWS->getDimension(i)->getMDFrame().name() !=
          Mantid::Geometry::QSample::QSampleName) {
        errorMessage.emplace("InputWorkspace",
                             "The input workspace must be in Q_sample");
      }
    }
  }
  // Check if the vanadium is available for diffraction
  bool diffraction = true;
  if ((inputWS->getNumDims() > 3) &&
      (inputWS->getDimension(3)->getName() == "DeltaE")) {
    diffraction = false;
  }
  if (diffraction) {
    API::MatrixWorkspace_const_sptr solidAngleWS =
        getProperty("SolidAngleWorkspace");
    API::MatrixWorkspace_const_sptr fluxWS = getProperty("FluxWorkspace");
    if (solidAngleWS == nullptr) {
      errorMessage.emplace("SolidAngleWorkspace",
                           "SolidAngleWorkspace is required for diffraction");
    }
    if (fluxWS == nullptr) {
      errorMessage.emplace("FluxWorkspace",
                           "FluxWorkspace is required for diffraction");
    }
  }
  // Check for property MDNorm_low and MDNorm_high
  size_t nExperimentInfos = inputWS->getNumExperimentInfo();
  if (nExperimentInfos == 0) {
    errorMessage.emplace("InputWorkspace",
                         "There must be at least one experiment info");
  } else {
    for (size_t iExpInfo = 0; iExpInfo < nExperimentInfos; iExpInfo++) {
      auto &currentExptInfo =
          *(inputWS->getExperimentInfo(static_cast<uint16_t>(iExpInfo)));
      if (!currentExptInfo.run().hasProperty("MDNorm_low")) {
        errorMessage.emplace("InputWorkspace", "Missing MDNorm_low log. Please "
                                               "use CropWorkspaceForMDNorm "
                                               "before converting to MD");
      }
      if (!currentExptInfo.run().hasProperty("MDNorm_high")) {
        errorMessage.emplace("InputWorkspace",
                             "Missing MDNorm_high log. Please use "
                             "CropWorkspaceForMDNorm before converting to MD");
      }
    }
  }
  // check projections and UB
  if (getProperty("RLU")) {
    DblMatrix W = DblMatrix(3, 3);
    std::vector<double> Q0Basis = getProperty("QDimension0");
    std::vector<double> Q1Basis = getProperty("QDimension1");
    std::vector<double> Q2Basis = getProperty("QDimension2");
    W.setColumn(0, Q0Basis);
    W.setColumn(1, Q1Basis);
    W.setColumn(2, Q2Basis);
    if (fabs(W.determinant()) < 1e-5) {
      errorMessage.emplace("QDimension0",
                           "The projection dimensions are coplanar or zero");
      errorMessage.emplace("QDimension1",
                           "The projection dimensions are coplanar or zero");
      errorMessage.emplace("QDimension2",
                           "The projection dimensions are coplanar or zero");
    }
    if (!inputWS->getExperimentInfo(0)->sample().hasOrientedLattice()) {
      errorMessage.emplace("InputWorkspace",
                           "There is no oriented lattice "
                           "associated with the input workspace. "
                           "Use SetUB algorithm");
    }
  }
  // check dimension names
  std::vector<std::string> originalDimensionNames;
  for (size_t i = 3; i < inputWS->getNumDims(); i++) {
    originalDimensionNames.push_back(inputWS->getDimension(i)->getName());
  }
  originalDimensionNames.push_back("QDimension0");
  originalDimensionNames.push_back("QDimension1");
  originalDimensionNames.push_back("QDimension2");
  std::vector<std::string> selectedDimensions;
  for (std::size_t i = 0; i < 6; i++) {
    std::string propName = "Dimension" + Strings::toString(i) + "Name";
    std::string dimName = getProperty(propName);
    std::string binningName = "Dimension" + Strings::toString(i) + "Binning";
    std::vector<double> binning = getProperty(binningName);
    if (!dimName.empty()) {
      auto it = std::find(originalDimensionNames.begin(),
                          originalDimensionNames.end(), dimName);
      if (it == originalDimensionNames.end()) {
        errorMessage.emplace(propName,
                             "Name '" + dimName +
                                 "' is not one of the "
                                 "original workspace names or a Q dimension");
      } else {
        // make sure dimension is unique
        auto itSel = std::find(selectedDimensions.begin(),
                               selectedDimensions.end(), dimName);
        if (itSel == selectedDimensions.end()) {
          selectedDimensions.push_back(dimName);
        } else {
          errorMessage.emplace(propName,
                               "Name '" + dimName + "' was already selected");
        }
      }
    } else {
      if (!binning.empty()) {
        errorMessage.emplace(
            binningName,
            "There should be no binning if the dimension name is empty");
      }
    }
  }
  // since Q dimensions can be non - orthogonal, all must be present
  if ((std::find(selectedDimensions.begin(), selectedDimensions.end(),
                 "QDimension0") == selectedDimensions.end()) ||
      (std::find(selectedDimensions.begin(), selectedDimensions.end(),
                 "QDimension1") == selectedDimensions.end()) ||
      (std::find(selectedDimensions.begin(), selectedDimensions.end(),
                 "QDimension2") == selectedDimensions.end())) {
    for (std::size_t i = 0; i < 6; i++) {
      std::string propName = "Dimension" + Strings::toString(i) + "Name";
      errorMessage.emplace(
          propName,
          "All of QDimension0, QDimension1, QDimension2 must be present");
    }
  }
  // symmetry operations
  std::string symOps = this->getProperty("SymmetryOperations");
  if (!symOps.empty()) {
    bool isSpaceGroup =
        Geometry::SpaceGroupFactory::Instance().isSubscribed(symOps);
    bool isPointGroup =
        Geometry::PointGroupFactory::Instance().isSubscribed(symOps);
    if (!isSpaceGroup && !isPointGroup) {
      try {
        Geometry::SymmetryOperationFactory::Instance().createSymOps(symOps);
      } catch (const Mantid::Kernel::Exception::ParseError &) {
        errorMessage.emplace("SymmetryOperations",
                             "The input is not a space group, a point group, "
                             "or a list of symmetry operations");
      }
    }
  }
  return errorMessage;
}

//----------------------------------------------------------------------------------------------
/** Execute the algorithm.
 */
void MDNorm::exec() {
  convention = Kernel::ConfigService::Instance().getString("Q.convention");
  // symmetry operations
  std::string symOps = this->getProperty("SymmetryOperations");
  std::vector<Geometry::SymmetryOperation> symmetryOps;
  if (symOps.empty()) {
    symOps = "x,y,z";
  }
  if (Geometry::SpaceGroupFactory::Instance().isSubscribed(symOps)) {
    auto spaceGroup =
        Geometry::SpaceGroupFactory::Instance().createSpaceGroup(symOps);
    auto pointGroup = spaceGroup->getPointGroup();
    symmetryOps = pointGroup->getSymmetryOperations();
  } else if (Geometry::PointGroupFactory::Instance().isSubscribed(symOps)) {
    auto pointGroup =
        Geometry::PointGroupFactory::Instance().createPointGroup(symOps);
    symmetryOps = pointGroup->getSymmetryOperations();
  } else {
    symmetryOps =
        Geometry::SymmetryOperationFactory::Instance().createSymOps(symOps);
  }
  g_log.debug() << "Symmetry operations\n";
  for (auto so : symmetryOps) {
    g_log.debug() << so.identifier() << "\n";
  }
  m_numSymmOps = symmetryOps.size();

  m_isRLU = getProperty("RLU");
  // get the workspaces
  m_inputWS = this->getProperty("InputWorkspace");
  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();
  if ((m_inputWS->getNumDims() > 3) &&
      (m_inputWS->getDimension(3)->getName() == "DeltaE")) {
    m_diffraction = false;
    if (exptInfoZero.run().hasProperty("Ei")) {
      Kernel::Property *eiprop = exptInfoZero.run().getProperty("Ei");
      m_Ei = boost::lexical_cast<double>(eiprop->value());
      if (m_Ei <= 0) {
        throw std::invalid_argument(
            "Ei stored in the workspace is not positive");
      }
    } else {
      throw std::invalid_argument("Could not find Ei value in the workspace.");
    }
  }
  auto outputDataWS = binInputWS(symmetryOps);

  createNormalizationWS(*outputDataWS);
  this->setProperty("OutputNormalizationWorkspace", m_normWS);
  this->setProperty("OutputDataWorkspace", outputDataWS);

  m_numExptInfos = outputDataWS->getNumExperimentInfo();
  // loop over all experiment infos
  for (uint16_t expInfoIndex = 0; expInfoIndex < m_numExptInfos;
       expInfoIndex++) {
    // Check for other dimensions if we could measure anything in the original
    // data
    bool skipNormalization = false;
    const std::vector<coord_t> otherValues =
        getValuesFromOtherDimensions(skipNormalization, expInfoIndex);

    cacheDimensionXValues();

    if (!skipNormalization) {
      size_t symmOpsIndex = 0;
      for (const auto &so : symmetryOps) {
        calculateNormalization(otherValues, so, expInfoIndex, symmOpsIndex);
        symmOpsIndex++;
      }

    } else {
      g_log.warning("Binning limits are outside the limits of the MDWorkspace. "
                    "Not applying normalization.");
    }
    // if more than one experiment info, keep accumulating
    m_accumulate = true;
  }
  IAlgorithm_sptr divideMD = createChildAlgorithm("DivideMD", 0.99, 1.);
  divideMD->setProperty("LHSWorkspace", outputDataWS);
  divideMD->setProperty("RHSWorkspace", m_normWS);
  divideMD->setPropertyValue("OutputWorkspace",
                             getPropertyValue("OutputWorkspace"));
  divideMD->executeAsChildAlg();
  API::IMDWorkspace_sptr out = divideMD->getProperty("OutputWorkspace");
  this->setProperty("OutputWorkspace", out);
}

/**
 *
 * @param projection - a vector with 3 elements, containing a
 *   description of the projection ("1,-1,0" for "[H,-H,0]")
 * @return string containing the name
 */
std::string MDNorm::QDimensionName(std::vector<double> projection) {
  std::vector<double>::iterator result;
  result = std::max_element(projection.begin(), projection.end(), abs_compare);
  std::vector<char> symbol{'H', 'K', 'L'};
  char character = symbol[std::distance(projection.begin(), result)];
  std::stringstream name;
  name << "[";
  for (size_t i = 0; i < 3; i++) {
    if (projection[i] == 0) {
      name << "0";
    } else if (projection[i] == 1) {
      name << character;
    } else if (projection[i] == -1) {
      name << "-" << character;
    } else {
      name << std::defaultfloat << std::setprecision(3) << projection[i]
           << character;
    }
    if (i != 2) {
      name << ",";
    }
  }
  name << "]";
  return name.str();
}

/**
 * Calculate binning parameters
 * @return map of parameters to be passed to BinMD (non axis-aligned)
 */
std::map<std::string, std::string> MDNorm::getBinParameters() {
  std::map<std::string, std::string> parameters;
  std::stringstream extents;
  std::stringstream bins;
  std::vector<std::string> originalDimensionNames;
  originalDimensionNames.push_back("QDimension0");
  originalDimensionNames.push_back("QDimension1");
  originalDimensionNames.push_back("QDimension2");
  for (size_t i = 3; i < m_inputWS->getNumDims(); i++) {
    originalDimensionNames.push_back(m_inputWS->getDimension(i)->getName());
  }

  if (m_isRLU) {
    m_Q0Basis = getProperty("QDimension0");
    m_Q1Basis = getProperty("QDimension1");
    m_Q2Basis = getProperty("QDimension2");
    m_UB =
        m_inputWS->getExperimentInfo(0)->sample().getOrientedLattice().getUB() *
        2 * M_PI;
  }

  std::vector<double> W(m_Q0Basis);
  W.insert(W.end(), m_Q1Basis.begin(), m_Q1Basis.end());
  W.insert(W.end(), m_Q2Basis.begin(), m_Q2Basis.end());
  m_W = DblMatrix(W);
  m_W.Transpose();

  // Find maximum Q
  auto &exptInfo0 = *(m_inputWS->getExperimentInfo(static_cast<uint16_t>(0)));
  auto upperLimitsVector =
      (*(dynamic_cast<Kernel::PropertyWithValue<std::vector<double>> *>(
          exptInfo0.getLog("MDNorm_high"))))();
  double maxQ;
  if (m_diffraction) {
    maxQ = 2. * (*std::max_element(upperLimitsVector.begin(),
                                   upperLimitsVector.end()));
  } else {
    double Ei;
    double maxDE =
        *std::max_element(upperLimitsVector.begin(), upperLimitsVector.end());
    auto loweLimitsVector =
        (*(dynamic_cast<Kernel::PropertyWithValue<std::vector<double>> *>(
            exptInfo0.getLog("MDNorm_low"))))();
    double minDE =
        *std::min_element(loweLimitsVector.begin(), loweLimitsVector.end());
    if (exptInfo0.run().hasProperty("Ei")) {
      Kernel::Property *eiprop = exptInfo0.run().getProperty("Ei");
      Ei = boost::lexical_cast<double>(eiprop->value());
      if (Ei <= 0) {
        throw std::invalid_argument(
            "Ei stored in the workspace is not positive");
      }
    } else {
      throw std::invalid_argument("Could not find Ei value in the workspace.");
    }
    const double energyToK = 8.0 * M_PI * M_PI *
                             PhysicalConstants::NeutronMass *
                             PhysicalConstants::meV * 1e-20 /
                             (PhysicalConstants::h * PhysicalConstants::h);
    double ki = std::sqrt(energyToK * Ei);
    double kfmin = std::sqrt(energyToK * (Ei - minDE));
    double kfmax = std::sqrt(energyToK * (Ei - maxDE));

    maxQ = ki + std::max(kfmin, kfmax);
  }
  size_t basisVectorIndex = 0;
  std::vector<coord_t> transformation;
  for (std::size_t i = 0; i < 6; i++) {
    std::string propName = "Dimension" + Strings::toString(i) + "Name";
    std::string binningName = "Dimension" + Strings::toString(i) + "Binning";
    std::string dimName = getProperty(propName);
    std::vector<double> binning = getProperty(binningName);
    std::string bv = "BasisVector";
    if (!dimName.empty()) {
      std::string property = bv + Strings::toString(basisVectorIndex);
      std::stringstream propertyValue;
      propertyValue << dimName;
      // get the index in the original workspace
      auto dimIndex =
          std::distance(originalDimensionNames.begin(),
                        std::find(originalDimensionNames.begin(),
                                  originalDimensionNames.end(), dimName));
      auto dimension = m_inputWS->getDimension(dimIndex);
      propertyValue << "," << dimension->getMDUnits().getUnitLabel().ascii();
      for (size_t j = 0; j < originalDimensionNames.size(); j++) {
        if (j == static_cast<size_t>(dimIndex)) {
          propertyValue << ",1";
          transformation.push_back(1.);
        } else {
          propertyValue << ",0";
          transformation.push_back(0.);
        }
      }
      parameters.emplace(property, propertyValue.str());
      // get the extents an number of bins
      coord_t dimMax = dimension->getMaximum();
      coord_t dimMin = dimension->getMinimum();
      if (m_isRLU) {
        Mantid::Geometry::OrientedLattice ol;
        ol.setUB(m_UB * m_W); // note that this is already multiplied by 2Pi
        if (dimIndex == 0) {
          dimMax = static_cast<coord_t>(ol.a() * maxQ);
          dimMin = -dimMax;
        } else if (dimIndex == 1) {
          dimMax = static_cast<coord_t>(ol.b() * maxQ);
          dimMin = -dimMax;
        } else if (dimIndex == 2) {
          dimMax = static_cast<coord_t>(ol.c() * maxQ);
          dimMin = -dimMax;
        }
      }
      if (binning.size() == 0) {
        // only one bin, integrating from min to max
        extents << dimMin << "," << dimMax << ",";
        bins << 1 << ",";
      } else if (binning.size() == 2) {
        // only one bin, integrating from min to max
        extents << binning[0] << "," << binning[1] << ",";
        bins << 1 << ",";
      } else if (binning.size() == 1) {
        auto step = binning[0];
        double nsteps = (dimMax - dimMin) / step;
        if (nsteps + 1 - std::ceil(nsteps) >= 1e-4) {
          nsteps = std::ceil(nsteps);
        } else {
          nsteps = std::floor(nsteps);
        }
        bins << static_cast<int>(nsteps) << ",";
        extents << dimMin << "," << dimMin + nsteps * step << ",";
      } else if (binning.size() == 3) {
        dimMin = static_cast<coord_t>(binning[0]);
        auto step = binning[1];
        dimMax = static_cast<coord_t>(binning[2]);
        double nsteps = (dimMax - dimMin) / step;
        if (nsteps + 1 - std::ceil(nsteps) >= 1e-4) {
          nsteps = std::ceil(nsteps);
        } else {
          nsteps = std::floor(nsteps);
        }
        bins << static_cast<int>(nsteps) << ",";
        extents << dimMin << "," << dimMin + nsteps * step << ",";
      }
      basisVectorIndex++;
    }
  }
  parameters.emplace("OutputExtents", extents.str());
  parameters.emplace("OutputBins", bins.str());
  m_transformation = Mantid::Kernel::Matrix<coord_t>(
      transformation,
      static_cast<size_t>((transformation.size()) / m_inputWS->getNumDims()),
      m_inputWS->getNumDims());
  return parameters;
}

/**
 * Create & cached the normalization workspace
 * @param dataWS The binned workspace that will be used for the data
 */
void MDNorm::createNormalizationWS(
    const DataObjects::MDHistoWorkspace &dataWS) {
  // Copy the MDHisto workspace, and change signals and errors to 0.
  boost::shared_ptr<IMDHistoWorkspace> tmp =
      this->getProperty("TemporaryNormalizationWorkspace");
  m_normWS = boost::dynamic_pointer_cast<MDHistoWorkspace>(tmp);
  if (!m_normWS) {
    m_normWS = dataWS.clone();
    m_normWS->setTo(0., 0., 0.);
  } else {
    m_accumulate = true;
  }
}

/**
 * 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
 */
DataObjects::MDHistoWorkspace_sptr
MDNorm::binInputWS(std::vector<Geometry::SymmetryOperation> symmetryOps) {
  Mantid::API::IMDHistoWorkspace_sptr tempDataWS =
      this->getProperty("TemporaryDataWorkspace");
  Mantid::API::Workspace_sptr outputWS;
  std::map<std::string, std::string> parameters = getBinParameters();
  double soIndex = 0;
  std::vector<size_t> qDimensionIndices;
  for (auto so : symmetryOps) {
    // calculate dimensions for binning
    DblMatrix soMatrix(3, 3);
    auto v = so.transformHKL(V3D(1, 0, 0));
    soMatrix.setColumn(0, v);
    v = so.transformHKL(V3D(0, 1, 0));
    soMatrix.setColumn(1, v);
    v = so.transformHKL(V3D(0, 0, 1));
    soMatrix.setColumn(2, v);

    DblMatrix Qtransform;
    if (m_isRLU) {
      Qtransform = m_UB * soMatrix * m_W;
    } else {
      Qtransform = soMatrix * m_W;
    }

    // bin the data
    double fraction = 1. / static_cast<double>(symmetryOps.size());
    IAlgorithm_sptr binMD = createChildAlgorithm(
        "BinMD", soIndex * 0.3 * fraction, (soIndex + 1) * 0.3 * fraction);
    binMD->setPropertyValue("AxisAligned", "0");
    binMD->setProperty("InputWorkspace", m_inputWS);
    binMD->setProperty("TemporaryDataWorkspace", tempDataWS);
    binMD->setPropertyValue("NormalizeBasisVectors", "0");
    binMD->setPropertyValue("OutputWorkspace",
                            getPropertyValue("OutputDataWorkspace"));
    // set binning properties
    size_t qindex = 0;
    for (auto const &p : parameters) {
      auto key = p.first;
      auto value = p.second;

      std::stringstream basisVector;
      std::vector<double> projection(m_inputWS->getNumDims(), 0.);
      if (value.find("QDimension0") != std::string::npos) {
        m_hIdx = qindex;
        if (!m_isRLU) {
          projection[0] = 1.;
          basisVector << "Q_sample_x,A^{-1}";
        } else {
          qDimensionIndices.push_back(qindex);
          projection[0] = Qtransform[0][0];
          projection[1] = Qtransform[1][0];
          projection[2] = Qtransform[2][0];
          basisVector << QDimensionName(m_Q0Basis) << ", r.l.u.";
        }
      } else if (value.find("QDimension1") != std::string::npos) {
        m_kIdx = qindex;
        if (!m_isRLU) {
          projection[1] = 1.;
          basisVector << "Q_sample_y,A^{-1}";
        } else {
          qDimensionIndices.push_back(qindex);
          projection[0] = Qtransform[0][1];
          projection[1] = Qtransform[1][1];
          projection[2] = Qtransform[2][1];
          basisVector << QDimensionName(m_Q1Basis) << ", r.l.u.";
        }
      } else if (value.find("QDimension2") != std::string::npos) {
        m_lIdx = qindex;
        if (!m_isRLU) {
          projection[2] = 1.;
          basisVector << "Q_sample_z,A^{-1}";
        } else {
          qDimensionIndices.push_back(qindex);
          projection[0] = Qtransform[0][2];
          projection[1] = Qtransform[1][2];
          projection[2] = Qtransform[2][2];
          basisVector << QDimensionName(m_Q2Basis) << ", r.l.u.";
        }
      } else if (value.find("DeltaE") != std::string::npos) {
        m_eIdx = qindex;
        m_dEIntegrated = false;
      }
      if (!basisVector.str().empty()) {
        for (auto const &proji : projection) {
          basisVector << "," << proji;
        }
        value = basisVector.str();
      }
      if (value.find("DeltaE") != std::string::npos) {
        m_eIdx = qindex;
      }
      g_log.debug() << "Binning parameter " << key << " value: " << value
                    << "\n";
      binMD->setPropertyValue(key, value);
      qindex++;
    }
    // execute algorithm
    binMD->executeAsChildAlg();
    outputWS = binMD->getProperty("OutputWorkspace");

    // set the temporary workspace to be the output workspace, so it keeps
    // adding different symmetries
    tempDataWS = boost::dynamic_pointer_cast<MDHistoWorkspace>(outputWS);
    soIndex += 1;
  }

  auto outputMDHWS = boost::dynamic_pointer_cast<MDHistoWorkspace>(outputWS);
  // set MDUnits for Q dimensions
  if (m_isRLU) {
    Mantid::Geometry::MDFrameArgument argument(Mantid::Geometry::HKL::HKLName,
                                               "r.l.u.");
    auto mdFrameFactory = Mantid::Geometry::makeMDFrameFactoryChain();
    Mantid::Geometry::MDFrame_uptr hklFrame = mdFrameFactory->create(argument);
    for (size_t i : qDimensionIndices) {
      auto mdHistoDimension = boost::const_pointer_cast<
          Mantid::Geometry::MDHistoDimension>(
          boost::dynamic_pointer_cast<const Mantid::Geometry::MDHistoDimension>(
              outputMDHWS->getDimension(i)));
      mdHistoDimension->setMDFrame(*hklFrame);
    }
  }

  outputMDHWS->setDisplayNormalization(Mantid::API::NoNormalization);
  return outputMDHWS;
}

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

  std::vector<coord_t> otherDimValues;
  for (size_t i = 3; i < m_inputWS->getNumDims(); i++) {
    const auto dimension = m_inputWS->getDimension(i);
    coord_t inputDimMin = static_cast<float>(dimension->getMinimum());
    coord_t inputDimMax = static_cast<float>(dimension->getMaximum());
    coord_t outputDimMin(0), outputDimMax(0);
    bool isIntegrated = true;

    for (size_t j = 0; j < m_transformation.numRows(); j++) {
      if (m_transformation[j][i] == 1) {
        isIntegrated = false;
        outputDimMin = m_normWS->getDimension(j)->getMinimum();
        outputDimMax = m_normWS->getDimension(j)->getMaximum();
      }
    }
    if (dimension->getName() == "DeltaE") {
      if ((inputDimMax < outputDimMin) || (inputDimMin > outputDimMax)) {
        skipNormalization = true;
      }
    } else {
      coord_t value = static_cast<coord_t>(currentRun.getLogAsSingleValue(
          dimension->getName(), Mantid::Kernel::Math::TimeAveragedMean));
      otherDimValues.push_back(value);
      if (value < inputDimMin || value > inputDimMax) {
        skipNormalization = true;
      }
      if ((!isIntegrated) && (value < outputDimMin || value > outputDimMax)) {
        skipNormalization = true;
      }
    }
  }
  return otherDimValues;
}

/**
 * Stores the X values from each H,K,L, and optionally DeltaE dimension as
 * member variables
 */
void MDNorm::cacheDimensionXValues() {
  auto &hDim = *m_normWS->getDimension(m_hIdx);
  m_hX.resize(hDim.getNBoundaries());
  for (size_t i = 0; i < m_hX.size(); ++i) {
    m_hX[i] = hDim.getX(i);
  }
  auto &kDim = *m_normWS->getDimension(m_kIdx);
  m_kX.resize(kDim.getNBoundaries());
  for (size_t i = 0; i < m_kX.size(); ++i) {
    m_kX[i] = kDim.getX(i);
  }

  auto &lDim = *m_normWS->getDimension(m_lIdx);
  m_lX.resize(lDim.getNBoundaries());
  for (size_t i = 0; i < m_lX.size(); ++i) {
    m_lX[i] = lDim.getX(i);
  }

  if ((!m_diffraction) && (!m_dEIntegrated)) {
    // NOTE: store k final instead
    auto &eDim = *m_normWS->getDimension(m_eIdx);
    m_eX.resize(eDim.getNBoundaries());
    for (size_t i = 0; i < m_eX.size(); ++i) {
      double temp = m_Ei - eDim.getX(i);
      temp = std::max(temp, 0.);
      m_eX[i] = std::sqrt(energyToK * temp);
    }
  }
}

/**
 * Computed the normalization for the input workspace. Results are stored in
 * m_normWS
 * @param otherValues - values for dimensions other than Q or DeltaE
 * @param so - symmetry operation
 * @param expInfoIndex - current experiment info index
 * @param soIndex - the index of symmetry operation (for progress purposes)
 */
void MDNorm::calculateNormalization(const std::vector<coord_t> &otherValues,
                                    Geometry::SymmetryOperation so,
                                    uint16_t expInfoIndex, size_t soIndex) {
  const auto &currentExptInfo = *(m_inputWS->getExperimentInfo(expInfoIndex));
  std::vector<double> lowValues, highValues;
  auto *lowValuesLog = dynamic_cast<VectorDoubleProperty *>(
      currentExptInfo.getLog("MDNorm_low"));
  lowValues = (*lowValuesLog)();
  auto *highValuesLog = dynamic_cast<VectorDoubleProperty *>(
      currentExptInfo.getLog("MDNorm_high"));
  highValues = (*highValuesLog)();

  DblMatrix R = currentExptInfo.run().getGoniometerMatrix();
  DblMatrix soMatrix(3, 3);
  auto v = so.transformHKL(V3D(1, 0, 0));
  soMatrix.setColumn(0, v);
  v = so.transformHKL(V3D(0, 1, 0));
  soMatrix.setColumn(1, v);
  v = so.transformHKL(V3D(0, 0, 1));
  soMatrix.setColumn(2, v);
  soMatrix.Invert();
  DblMatrix Qtransform = R * m_UB * soMatrix * m_W;
  Qtransform.Invert();
  const double protonCharge = currentExptInfo.run().getProtonCharge();
  const auto &spectrumInfo = currentExptInfo.spectrumInfo();

  // Mappings
  const int64_t ndets = static_cast<int64_t>(spectrumInfo.size());
  detid2index_map fluxDetToIdx;
  detid2index_map solidAngDetToIdx;
  bool haveSA = false;
  API::MatrixWorkspace_const_sptr solidAngleWS =
      getProperty("SolidAngleWorkspace");
  API::MatrixWorkspace_const_sptr integrFlux = getProperty("FluxWorkspace");
  if (solidAngleWS != nullptr) {
    haveSA = true;
    solidAngDetToIdx = solidAngleWS->getDetectorIDToWorkspaceIndexMap();
  }
  if (m_diffraction) {
    fluxDetToIdx = integrFlux->getDetectorIDToWorkspaceIndexMap();
  }

  const size_t vmdDims = (m_diffraction) ? 3 : 4;
  std::vector<std::atomic<signal_t>> signalArray(m_normWS->getNPoints());
  std::vector<std::array<double, 4>> intersections;
  std::vector<double> xValues, yValues;
  std::vector<coord_t> pos, posNew;

  double progStep = 0.7 / static_cast<double>(m_numExptInfos * m_numSymmOps);
  double progIndex = static_cast<double>(soIndex + expInfoIndex * m_numSymmOps);
  auto prog =
      make_unique<API::Progress>(this, 0.3 + progStep * progIndex,
                                 0.3 + progStep * (1. + progIndex), ndets);

  bool safe = true;
  if (m_diffraction) {
    safe = Kernel::threadSafe(*integrFlux);
  }
  // cppcheck-suppress syntaxError
PRAGMA_OMP(parallel for private(intersections, xValues, yValues, pos, posNew) if (safe))
for (int64_t i = 0; i < ndets; i++) {
  PARALLEL_START_INTERUPT_REGION

  if (!spectrumInfo.hasDetectors(i) || spectrumInfo.isMonitor(i) ||
      spectrumInfo.isMasked(i)) {
    continue;
  }

  const auto &detector = spectrumInfo.detector(i);
  double theta = detector.getTwoTheta(m_samplePos, m_beamDir);
  double phi = detector.getPhi();
  // If the dtefctor is a group, this should be the ID of the first detector
  const auto detID = detector.getID();

  // Intersections
  this->calculateIntersections(intersections, theta, phi, Qtransform,
                               lowValues[i], highValues[i]);
  if (intersections.empty())
    continue;
  // Get solid angle for this contribution
  double solid = protonCharge;
  if (haveSA) {
    solid =
        solidAngleWS->y(solidAngDetToIdx.find(detID)->second)[0] * protonCharge;
  }

  if (m_diffraction) {
    // -- calculate integrals for the intersection --
    // momentum values at intersections
    auto intersectionsBegin = intersections.begin();
    // copy momenta to xValues
    xValues.resize(intersections.size());
    yValues.resize(intersections.size());
    auto x = xValues.begin();
    for (auto it = intersectionsBegin; it != intersections.end(); ++it, ++x) {
      *x = (*it)[3];
    }
    // get the flux spetrum number
    size_t wsIdx = fluxDetToIdx.find(detID)->second;
    // 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
  // pre-allocate for efficiency and copy non-hkl dim values into place