MDNorm.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/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/UnitLabelTypes.h"
#include "MantidKernel/VectorHelper.h"
#include "MantidKernel/VisibleWhenProperty.h"
 
#include <algorithm>
#include <boost/lexical_cast.hpp>
 
namespace Mantid::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>>;
 
// 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 doubles
static bool abs_compare(double a, double b) { return (std::fabs(a) < std::fabs(b)); }
 
bool isMonochromaticSCDInstrument(const IMDEventWorkspace &workspace) {
  if (workspace.getNumExperimentInfo() == 0)
    return false;
  const std::string instrument = workspace.getExperimentInfo(0)->getInstrumentName();
  return instrument == "WAND" || instrument == "HB2C" || instrument == "HB3A" || instrument == "DEMAND";
}
} // namespace
 
// Register the algorithm into the AlgorithmFactory
DECLARE_ALGORITHM(MDNorm)
 
//----------------------------------------------------------------------------------------------
MDNorm::MDNorm() : m_isRLU(false), m_monochromatic(false) {}
 
const std::string MDNorm::name() const { return "MDNorm"; }
 
int MDNorm::version() const { return 1; }
 
const std::string MDNorm::category() const { return "MDAlgorithms\\Normalisation"; }
 
const std::string MDNorm::summary() const {
  return "Bins multidimensional data and calculate the normalization on the "
         "same grid";
}
 
//----------------------------------------------------------------------------------------------
void MDNorm::init() {
  declareProperty(
      std::make_unique<WorkspaceProperty<API::IMDEventWorkspace>>("InputWorkspace", "", Kernel::Direction::Input),
      "An input MDEventWorkspace. Must be in Q_sample frame.");
 
  declareProperty(std::make_unique<WorkspaceProperty<API::IMDEventWorkspace>>(
                      "BackgroundWorkspace", "", Kernel::Direction::Input, PropertyMode::Optional),
                  "An (optional) input MDEventWorkspace for background.  Must be in Q_lab frame.");
 
  // RLU and settings
  declareProperty("RLU", true, "Use reciprocal lattice units. If false, use Q_sample");
  setPropertyGroup("RLU", "Q projections RLU");
 
  auto mustBe3D = std::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(std::make_unique<ArrayProperty<double>>("QDimension0", Q0, mustBe3D),
                  "The first Q projection axis - Default is (1,0,0)");
  setPropertySettings("QDimension0", std::make_unique<Kernel::VisibleWhenProperty>("RLU", IS_EQUAL_TO, "1"));
  setPropertyGroup("QDimension0", "Q projections RLU");
 
  declareProperty(std::make_unique<ArrayProperty<double>>("QDimension1", Q1, mustBe3D),
                  "The second Q projection axis - Default is (0,1,0)");
  setPropertySettings("QDimension1", std::make_unique<Kernel::VisibleWhenProperty>("RLU", IS_EQUAL_TO, "1"));
  setPropertyGroup("QDimension1", "Q projections RLU");
 
  declareProperty(std::make_unique<ArrayProperty<double>>("QDimension2", Q2, mustBe3D),
                  "The thirdtCalculateCover Q projection axis - Default is (0,0,1)");
  setPropertySettings("QDimension2", std::make_unique<Kernel::VisibleWhenProperty>("RLU", IS_EQUAL_TO, "1"));
  setPropertyGroup("QDimension2", "Q projections RLU");
 
  // vanadium
  auto fluxValidator = std::make_shared<CompositeValidator>();
  fluxValidator->add<InstrumentValidator>();
  fluxValidator->add<CommonBinsValidator>();
  auto solidAngleValidator = fluxValidator->clone();
  declareProperty(std::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(std::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");
 
  // monochromatic single crystal diffraction (WAND, DEMAND)
  declareProperty(std::make_unique<WorkspaceProperty<API::IMDEventWorkspace>>("MonoSCDNormalizationWorkspace", "",
                                                                              Direction::Input, PropertyMode::Optional),
                  "An (optional) input MDEventWorkspace containing a pre-computed normalization "
                  "for monochromatic single crystal diffraction (e.g. produced by "
                  "ConvertHFIRSCDtoMDE). Must be in Q_sample frame with the same number of "
                  "dimensions as InputWorkspace. Cannot be used together with "
                  "SolidAngleWorkspace/FluxWorkspace or BackgroundWorkspace.");
  setPropertyGroup("MonoSCDNormalizationWorkspace", "monochromatic-SCD");
 
  // 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";
    std::string defaultName = "";
    if (i < 3) {
      defaultName = "QDimension" + Strings::toString(i);
    }
    declareProperty(std::make_unique<PropertyWithValue<std::string>>(propName, defaultName, Direction::Input),
                    "Name for the " + Strings::toString(i) + "th dimension. Leave blank for NONE.");
    auto atMost3 = std::make_shared<ArrayLengthValidator<double>>(0, 3);
    std::vector<double> temp;
    declareProperty(std::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(std::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(std::make_unique<WorkspaceProperty<IMDHistoWorkspace>>("TemporaryDataWorkspace", "", Direction::Input,
                                                                         PropertyMode::Optional),
                  "An (optional) input MDHistoWorkspace used to accumulate data from "
                  "multiple MDEventWorkspaces. If unspecified a blank "
                  "MDHistoWorkspace will be created.");
  declareProperty(std::make_unique<WorkspaceProperty<IMDHistoWorkspace>>("TemporaryNormalizationWorkspace", "",
                                                                         Direction::Input, PropertyMode::Optional),
                  "An (optional) input MDHistoWorkspace used to accumulate normalization "
                  "from multiple MDEventWorkspaces. If unspecified a blank "
                  "MDHistoWorkspace will be created.");
 
  // temporary background workspace
  declareProperty(std::make_unique<WorkspaceProperty<IMDHistoWorkspace>>("TemporaryBackgroundDataWorkspace", "",
                                                                         Direction::Input, PropertyMode::Optional),
                  "An (optional) input MDHistoWorkspace used to accumulate background from "
                  "multiple background MDEventWorkspaces. If unspecified but "
                  "BackgroundWorkspace is specified, a blank "
                  "MDHistoWorkspace will be created.");
  declareProperty(std::make_unique<WorkspaceProperty<IMDHistoWorkspace>>("TemporaryBackgroundNormalizationWorkspace",
                                                                         "", Direction::Input, PropertyMode::Optional),
                  "An (optional) input MDHistoWorkspace used to accumulate background normalization "
                  "from multiple background MDEventWorkspaces. If unspecified but "
                  "BackgroundWorkspace is specified, a blank "
                  "MDHistoWorkspace will be created.");
 
  setPropertyGroup("TemporaryDataWorkspace", "Temporary workspaces");
  setPropertyGroup("TemporaryNormalizationWorkspace", "Temporary workspaces");
  setPropertyGroup("TemporaryBackgroundDataWorkspace", "Temporary workspaces");
  setPropertyGroup("TemporaryBackgroundNormalizationWorkspace", "Temporary workspaces");
 
  declareProperty(std::make_unique<WorkspaceProperty<API::Workspace>>("OutputWorkspace", "", Kernel::Direction::Output),
                  "A name for the normalized output MDHistoWorkspace.");
  declareProperty(
      std::make_unique<WorkspaceProperty<API::Workspace>>("OutputDataWorkspace", "", Kernel::Direction::Output),
      "A name for the output data MDHistoWorkspace.");
  declareProperty(std::make_unique<WorkspaceProperty<Workspace>>("OutputNormalizationWorkspace", "", Direction::Output),
                  "A name for the output normalization MDHistoWorkspace.");
  declareProperty(std::make_unique<WorkspaceProperty<API::Workspace>>(
                      "OutputBackgroundDataWorkspace", "", Kernel::Direction::Output, PropertyMode::Optional),
                  "A name for the optional output background data MDHistoWorkspace.");
  declareProperty(std::make_unique<WorkspaceProperty<Workspace>>("OutputBackgroundNormalizationWorkspace", "",
                                                                 Direction::Output, PropertyMode::Optional),
                  "A name for the optional output background normalization MDHistoWorkspace.");
 
  m_progress = std::make_unique<API::Progress>(this, 0, 1, 1);
}
 
//----------------------------------------------------------------------------------------------
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");
      }
    }
  }
 
  // Optional background input IMDE
  Mantid::API::IMDEventWorkspace_sptr bkgdWS = this->getProperty("BackgroundWorkspace");
  if (bkgdWS) {
    if (bkgdWS->getNumDims() < 3) {
      // must have at least 3 dimensions
      errorMessage.emplace("BackgroundWorkspace", "The input background workspace must be at least 3D");
    } else {
      // Check first 3 dimension for Q lab,
      for (size_t i = 0; i < 3; i++) {
        if (bkgdWS->getDimension(i)->getMDFrame().name() != Mantid::Geometry::QLab::QLabName) {
          errorMessage.emplace("BackgroundWorkspace", "The input backgound workspace must be in Q_lab");
        }
      }
 
      // Check 4th dimension if input workspace is elastic
      if (inputWS->getNumDims() > 3) {
        if (bkgdWS->getNumDims() <= 3) {
          errorMessage.emplace("BackgroundWorkspace", "The input background workspace must have at 4 dimensions when "
                                                      "input workspace has more than 4 dimensions (inelastic case).");
        } else if (bkgdWS->getDimension(3)->getName() != inputWS->getDimension(3)->getName()) {
          errorMessage.emplace("BackgroundWorkspace", "The input background workspace 4th dimension must be DeltaE "
                                                      "for inelastic case.");
        }
      }
    }
  }
 
  // Check if the vanadium is available for diffraction
  bool diffraction = true;
  if ((inputWS->getNumDims() > 3) && (inputWS->getDimension(3)->getName() == "DeltaE")) {
    diffraction = false;
  }
 
  // Optional pre-computed normalization workspace for monochromatic single crystal diffraction
  // (e.g. WAND, DEMAND). This is an alternative to SolidAngleWorkspace/FluxWorkspace.
  Mantid::API::IMDEventWorkspace_sptr monoNormWS = this->getProperty("MonoSCDNormalizationWorkspace");
  bool monochromatic = bool(monoNormWS);
  API::MatrixWorkspace_const_sptr solidAngleWS = getProperty("SolidAngleWorkspace");
  API::MatrixWorkspace_const_sptr fluxWS = getProperty("FluxWorkspace");
 
  // No background for monochromatic single crystal diffraction: it is subtracted in detector space instead. This
  // message replaces any other problem found with BackgroundWorkspace, since the property is not to be used at all.
  if (bkgdWS && (monochromatic || isMonochromaticSCDInstrument(*inputWS))) {
    errorMessage["BackgroundWorkspace"] = "MDNorm does not support a background for monochromatic single crystal "
                                          "diffraction. Subtract the background in detector space with "
                                          "HFIRSCDSubtractBackground before ConvertHFIRSCDtoMDE.";
  }
 
  if (monochromatic) {
    if (!diffraction) {
      errorMessage.emplace("MonoSCDNormalizationWorkspace", "MonoSCDNormalizationWorkspace can only be used for "
                                                            "diffraction (InputWorkspace must not have a DeltaE "
                                                            "dimension)");
    }
    if (solidAngleWS || fluxWS) {
      errorMessage.emplace("MonoSCDNormalizationWorkspace", "MonoSCDNormalizationWorkspace cannot be used together "
                                                            "with SolidAngleWorkspace/FluxWorkspace");
    }
    if (monoNormWS->getNumDims() < 3) {
      errorMessage.emplace("MonoSCDNormalizationWorkspace", "MonoSCDNormalizationWorkspace must be at least 3D");
    } else {
      if (monoNormWS->getNumDims() != inputWS->getNumDims()) {
        errorMessage.emplace("MonoSCDNormalizationWorkspace", "MonoSCDNormalizationWorkspace must have the same "
                                                              "number of dimensions as InputWorkspace");
      }
      for (size_t i = 0; i < 3; i++) {
        if (monoNormWS->getDimension(i)->getMDFrame().name() != Mantid::Geometry::QSample::QSampleName) {
          errorMessage.emplace("MonoSCDNormalizationWorkspace", "MonoSCDNormalizationWorkspace must be in Q_sample");
        }
      }
    }
  } else if (diffraction) {
    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 (TOF only), or, for monochromatic
  // input, that a wavelength log is present (set by e.g. ConvertHFIRSCDtoMDE)
  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 (monochromatic) {
        if (!currentExptInfo.run().hasProperty("wavelength")) {
          errorMessage.emplace("InputWorkspace", "Missing wavelength log. InputWorkspace does not look like it was "
                                                 "produced by ConvertHFIRSCDtoMDE, as expected when "
                                                 "MonoSCDNormalizationWorkspace is provided");
        }
      } else {
        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.emplace_back(inputWS->getDimension(i)->getName());
  }
  originalDimensionNames.emplace_back("QDimension0");
  originalDimensionNames.emplace_back("QDimension1");
  originalDimensionNames.emplace_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 directional dimension");
      } else {
        // make sure dimension is unique
        auto itSel = std::find(selectedDimensions.begin(), selectedDimensions.end(), dimName);
        if (itSel == selectedDimensions.end()) {
          selectedDimensions.emplace_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");
      }
    }
  }
  // validate accumulation workspaces, if provided
  std::shared_ptr<IMDHistoWorkspace> tempNormWS = this->getProperty("TemporaryNormalizationWorkspace");
  Mantid::API::IMDHistoWorkspace_sptr tempDataWS = this->getProperty("TemporaryDataWorkspace");
 
  // check that either both or neuther accumulation workspaces are provied
  if ((tempNormWS && !tempDataWS) || (!tempNormWS && tempDataWS)) {
    errorMessage.emplace("TemporaryDataWorkspace", "Must provide either no accumulation workspaces or,"
                                                   "both TemporaryNormalizationWorkspaces and TemporaryDataWorkspace");
  }
  // check that both accumulation workspaces are on the same grid
  if (tempNormWS && tempDataWS) {
    size_t numNormDims = tempNormWS->getNumDims();
    size_t numDataDims = tempDataWS->getNumDims();
    if (numNormDims == numDataDims) {
      for (size_t i = 0; i < numNormDims; i++) {
        const auto dim1 = tempNormWS->getDimension(i);
        const auto dim2 = tempDataWS->getDimension(i);
        if ((dim1->getMinimum() != dim2->getMinimum()) || (dim1->getMaximum() != dim2->getMaximum()) ||
            (dim1->getNBins() != dim2->getNBins()) || (dim1->getName() != dim2->getName())) {
          errorMessage.emplace("TemporaryDataWorkspace", "Binning for TemporaryNormalizationWorkspaces "
                                                         "and TemporaryDataWorkspace must be the same.");
          break;
        }
      }
    } else { // accumulation workspaces have different number of dimensions
      errorMessage.emplace("TemporaryDataWorkspace", "TemporaryNormalizationWorkspace and TemporaryDataWorkspace "
                                                     "do not have the same number of dimensions");
    }
  }
 
  // validate accumulated background workspaces
  Mantid::API::IMDHistoWorkspace_sptr tempBkgdDataWS = this->getProperty("TemporaryBackgroundDataWorkspace");
  Mantid::API::IMDHistoWorkspace_sptr tempBkgdNormWS = this->getProperty("TemporaryBackgroundNormalizationWorkspace");
  // check existing criteria: Background, TempBackgroundData and
  // TempBackgroundNormalization must be specified
  if (tempBkgdDataWS && (!bkgdWS || !tempDataWS || !tempBkgdNormWS)) {
    errorMessage.emplace("TemporaryBackgroundDataWorkspace", "TemporaryBackgroundDataWorkspace is specified but at "
                                                             "least one of these is not.");
  } else if (tempBkgdNormWS && (!bkgdWS || !tempNormWS || !tempBkgdDataWS)) {
    errorMessage.emplace("TemporaryBackgroundNormalizationWorkspace", "TemporaryBackgroundNormalizationWorkspace is "
                                                                      "specified but at least one of these is not.");
  } else if (bkgdWS && tempDataWS && !tempBkgdDataWS) {
    errorMessage.emplace("TemporaryDataWorkspace",
                         "With Background is specifed and TemporaryDataWorkspace is specifed, "
                         "TemporaryBackgroundDataWorkspace must be specified.");
  } else if (tempBkgdDataWS && tempNormWS) {
    // check when they both exist
    size_t numBkgdDataDims = tempBkgdDataWS->getNumDims();
    size_t numBkgdNormDims = tempBkgdNormWS->getNumDims();
    size_t numDataDims = tempDataWS->getNumDims();
    if (numBkgdDataDims == numBkgdNormDims && numBkgdDataDims == numDataDims) {
      // On each dimension, compare min, max, NBins and name
      for (size_t idim = 0; idim < numBkgdDataDims; ++idim) {
        const auto dimB = tempBkgdDataWS->getDimension(idim);
        const auto dimN = tempBkgdNormWS->getDimension(idim);
        const auto dimD = tempDataWS->getDimension(idim);
        if ((dimB->getMinimum() != dimN->getMinimum()) || (dimB->getMinimum() != dimD->getMinimum()) ||
            (dimB->getMaximum() != dimN->getMaximum()) || (dimB->getMaximum() != dimD->getMaximum()) ||
            (dimB->getNBins() != dimN->getNBins()) || (dimB->getNBins() != dimD->getNBins()) ||
            (dimB->getName() != dimN->getName()) || (dimB->getName() != dimD->getName())) {
          errorMessage.emplace("TemporaryBackgroundDataWorkspace",
                               "TemporaryBackgroundDataWorkspace, "
                               "TemporaryBackgroundNormalizationWorkspace and "
                               "TemporaryDataWorkspace "
                               "must have same minimum, maximum, number of bins and name.");
          break;
        }
      }
    } else {
      errorMessage.emplace("TemporaryBackgroundDataWorkspace", "TemporaryBackgroundDataWorkspace, "
                                                               "TemporaryBackgroundNormalizationWorkspace and "
                                                               "TemporaryDataWorkspace must have same dimensions");
    }
  }
 
  return errorMessage;
}
 
//----------------------------------------------------------------------------------------------
void MDNorm::exec() {
  m_convention = Kernel::ConfigService::Instance().getString("Q.convention");
  m_hIntegrated = m_kIntegrated = m_lIntegrated = false;
  // 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 (const 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");
  Mantid::API::IMDEventWorkspace_sptr monoNormInputWS = this->getProperty("MonoSCDNormalizationWorkspace");
  m_monochromatic = bool(monoNormInputWS);
  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 = normalize(m_samplePos - source->getPos());
  if ((m_inputWS->getNumDims() > 3) && (m_inputWS->getDimension(3)->getName() == "DeltaE")) {
    // DeltaE in input MDE: it cannot be diffraction!
    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.");
    }
  }
 
  // Calculate (BinMD) input sample MDE to MDH and create noramlization MDH from
  // it
  auto outputDataWS = binInputWS(symmetryOps);
  if (m_monochromatic) {
    // Monochromatic single crystal diffraction (WAND, DEMAND): the normalization is a
    // pre-computed MDEventWorkspace, binned identically to the data instead of being
    // calculated from solid angle/flux trajectories.
    m_normWS = binMonoSCDNormalizationWS(symmetryOps);
  } else {
    createNormalizationWS(*outputDataWS);
  }
  this->setProperty("OutputNormalizationWorkspace", m_normWS);
  this->setProperty("OutputDataWorkspace", outputDataWS);
 
  // Background
  m_backgroundWS = this->getProperty("BackgroundWorkspace");
  DataObjects::MDHistoWorkspace_sptr outputBackgroundDataWS(nullptr);
  // Outputs for background related
  if (m_backgroundWS) {
    outputBackgroundDataWS = binBackgroundWS(symmetryOps);
    createBackgroundNormalizationWS(*outputBackgroundDataWS);
    this->setProperty("OutputBackgroundNormalizationWorkspace", m_bkgdNormWS);
    this->setProperty("OutputBackgroundDataWorkspace", outputBackgroundDataWS);
  }
 
  m_numExptInfos = outputDataWS->getNumExperimentInfo();
  if (!m_monochromatic) {
    m_signalArray = std::vector<std::atomic<signal_t>>(m_normWS->getNPoints());
    if (m_backgroundWS)
      m_bkgdSignalArray = std::vector<std::atomic<signal_t>>(m_bkgdNormWS->getNPoints());
    // loop over all experiment infos, computing the normalization from solid angle/flux
    // trajectories (TOF only; for monochromatic input, m_normWS was already binned above)
    m_progress->resetNumSteps(m_numExptInfos * m_numSymmOps, 0.3, 0.9);
    for (uint16_t expInfoIndex = 0; expInfoIndex < m_numExptInfos; expInfoIndex++) {
      const auto &currentExptInfo = *(m_inputWS->getExperimentInfo(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) {
        for (const auto &so : symmetryOps) {
          if (currentExptInfo.run().hasProperty("useLogTimes")) {
            calculateNormContinuous(otherValues, expInfoIndex, &so);
          } else {
            calculateNormalization(otherValues, so, expInfoIndex);
          }
          m_progress->report();
        }
      } else {
        g_log.warning("Binning limits are outside the limits of the MDWorkspace. "
                      "Not applying normalization.");
      }
    }
    if (m_accumulate) {
      std::transform(m_signalArray.cbegin(), m_signalArray.cend(), m_normWS->getSignalArray(),
                     m_normWS->mutableSignalArray(),
                     [](const std::atomic<signal_t> &a, const signal_t &b) { return a + b; });
      if (m_backgroundWS)
        std::transform(m_bkgdSignalArray.cbegin(), m_bkgdSignalArray.cend(), m_bkgdNormWS->getSignalArray(),
                       m_bkgdNormWS->mutableSignalArray(),
                       [](const std::atomic<signal_t> &a, const signal_t &b) { return a + b; });
    } else {
      std::copy(m_signalArray.cbegin(), m_signalArray.cend(), m_normWS->mutableSignalArray());
      if (m_backgroundWS)
        std::copy(m_bkgdSignalArray.cbegin(), m_bkgdSignalArray.cend(), m_bkgdNormWS->mutableSignalArray());
    }
  }
 
  API::IMDWorkspace_sptr out(nullptr);
 
  if (m_backgroundWS) {
    // Normalize binned (BinMD) sample workspace with background
    out = divideMD(outputDataWS, m_normWS, getPropertyValue("OutputWorkspace"), 0.97, 0.98);
 
    // Normalize background
    const std::string normedBkgdWSName("_normedBkgd");
    API::IMDWorkspace_sptr outbkgd = divideMD(outputBackgroundDataWS, m_bkgdNormWS, normedBkgdWSName, 0.98, 0.99);
 
    // Clean workspace
    auto minusMD = createChildAlgorithm("MinusMD", 0.99, 1.00);
    // set up
    minusMD->setProperty("LHSWorkspace", out);
    minusMD->setProperty("RHSWorkspace", outbkgd);
    minusMD->setPropertyValue("OutputWorkspace", getPropertyValue("OutputWorkspace"));
    // run and return
    minusMD->executeAsChildAlg();
    out = minusMD->getProperty("OutputWorkspace");
 
  } else {
    // Normalize binned (BinMD) sample workspace without background
    out = divideMD(outputDataWS, m_normWS, getPropertyValue("OutputWorkspace"), 0.97, 1.);
  }
 
  // Set output workspace
  this->setProperty("OutputWorkspace", out);
}
 
inline API::IMDWorkspace_sptr MDNorm::divideMD(const API::IMDHistoWorkspace_sptr &lhs,
                                               const API::IMDHistoWorkspace_sptr &rhs, const std::string &outputwsname,
                                               const double &startProgress, const double &endProgress) {
  auto divideMD = createChildAlgorithm("DivideMD", startProgress, endProgress);
  divideMD->setProperty("LHSWorkspace", lhs);
  divideMD->setProperty("RHSWorkspace", rhs);
  divideMD->setPropertyValue("OutputWorkspace", outputwsname);
  divideMD->executeAsChildAlg();
  // API::IMDWorkspace_sptr
  API::IMDWorkspace_sptr out = divideMD->getProperty("OutputWorkspace");
 
  return out;
}
 
std::string MDNorm::QDimensionNameQSample(int i) {
  if (i == 0)
    return std::string("Q_sample_x");
  else if (i == 1)
    return std::string("Q_sample_y");
  else if (i == 2)
    return std::string("Q_sample_z");
  else
    throw std::invalid_argument("Index must be 0, 1, or 2 for QDimensionNameQSample");
}
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();
}
 
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.emplace_back("QDimension0");
  originalDimensionNames.emplace_back("QDimension1");
  originalDimensionNames.emplace_back("QDimension2");
  for (size_t i = 3; i < m_inputWS->getNumDims(); i++) {
    originalDimensionNames.emplace_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, an isotropic bound used below to set the default HKL bin extents
  // (only consumed when a Q dimension's binning is left automatic/step-only, see the
  // m_isRLU branch below). For TOF input this comes from the physical wavelength/TOF
  // window each detector reaches (MDNorm_low/MDNorm_high logs, set by
  // CropWorkspaceForMDNorm). Monochromatic input has no such trajectory -- each event
  // is a single measured Q_sample point -- so maxQ is instead estimated from the
  // workspace's own data-occupied extents: the box-tree's tight bounding box (not the
  // nominal/arbitrary MinValues/MaxValues used at MD-conversion time), converted to an
  // isotropic bound by taking the modulus of its farthest corner from the origin.
  double maxQ = 0.;
  if (!m_monochromatic) {
    auto &exptInfo0 = *(m_inputWS->getExperimentInfo(static_cast<uint16_t>(0)));
    auto upperLimitsVector =
        (*(dynamic_cast<Kernel::PropertyWithValue<std::vector<double>> *>(exptInfo0.getLog("MDNorm_high"))))();
    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.");
      }
      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);
    }
  } else {
    // getMinimumExtents() walks the box tree to the given depth and unions the extents of
    // populated boxes found there; its default depth of 2 can be much looser than the box
    // tree's actual depth for sparser inputs (measured up to ~3.7x tighter per-axis bounds by
    // depth 4 on real WAND/DEMAND data), while costing at most tens of milliseconds even for
    // O(1e7)-event workspaces, so a fixed depth of 4 is used here instead of the default.
    auto dataExtents = m_inputWS->getMinimumExtents(4);
    double qx = std::max(std::fabs(dataExtents[0].getMin()), std::fabs(dataExtents[0].getMax()));
    double qy = std::max(std::fabs(dataExtents[1].getMin()), std::fabs(dataExtents[1].getMax()));
    double qz = std::max(std::fabs(dataExtents[2].getMin()), std::fabs(dataExtents[2].getMax()));
    maxQ = std::sqrt(qx * qx + qy * qy + qz * qz);
  }
  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.emplace_back(1.f);
        } else {
          propertyValue << ",0";
          transformation.emplace_back(0.f);
        }
      }
      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;
}
 
void MDNorm::createBackgroundNormalizationWS(const DataObjects::MDHistoWorkspace &bkgdDataWS) {
 
  // requiring background workspace is specified
  if (!m_backgroundWS) {
    return;
  }
 
  // Copy the MDHisto workspace, and change signals and errors to 0.
  std::shared_ptr<IMDHistoWorkspace> tmp = this->getProperty("TemporaryBackgroundNormalizationWorkspace");
  m_bkgdNormWS = std::dynamic_pointer_cast<MDHistoWorkspace>(tmp);
  if (!m_bkgdNormWS) {
    m_bkgdNormWS = bkgdDataWS.clone();
    m_bkgdNormWS->setTo(0., 0., 0.);
  }
}
 
void MDNorm::validateBinningForTemporaryDataWorkspace(const std::map<std::string, std::string> &parameters,
                                                      const Mantid::API::IMDHistoWorkspace_sptr &tempDataWS) {
 
  // parse the paramters map and get extents from tempDataWS
  const std::string numBinsStr = parameters.at("OutputBins");
  const std::string extentsStr = parameters.at("OutputExtents");
  const std::vector<size_t> numBins = VectorHelper::splitStringIntoVector<size_t>(numBinsStr);
  const std::vector<double> extents = VectorHelper::splitStringIntoVector<double>(extentsStr);
 
  // make sure the number of dimensions is the same for both workspaces
  size_t numDimsTemp = tempDataWS->getNumDims();
  if ((numBins.size() != numDimsTemp) || (extents.size() != numDimsTemp * 2)) {
    std::stringstream errorMessage;
    errorMessage << "The number of dimensions in the output and ";
    errorMessage << "TemporaryDataWorkspace are not the same.";
    throw(std::invalid_argument(errorMessage.str()));
  }
 
  // compare the extents and number of bins
  for (size_t i = 0; i < numDimsTemp; i++) {
    auto ax = tempDataWS->getDimension(i);
    if (numBins[i] != ax->getNBins()) {
      std::stringstream errorMessage;
      errorMessage << "The number of bins output and number of bins in ";
      errorMessage << "TemporaryDataWorkspace are not the same along ";
      errorMessage << "dimension " << i;
      throw(std::invalid_argument(errorMessage.str()));
    }
    if (std::abs(extents[2 * i] - ax->getMinimum()) > 1.e-5) {
      std::stringstream errorMessage;
      errorMessage << "The minimum binning value for the output and ";
      errorMessage << "TemporaryDataWorkspace are not the same along ";
      errorMessage << "dimension " << i;
      throw(std::invalid_argument(errorMessage.str()));
    }
    if (std::abs(extents[2 * i + 1] - ax->getMaximum()) > 1.e-5) {
      std::stringstream errorMessage;
      errorMessage << "The maximum binning value for the output and ";
      errorMessage << "TemporaryDataWorkspace are not the same along ";
      errorMessage << "dimension " << i;
      throw(std::invalid_argument(errorMessage.str()));
    }
  }
 
  // sort out which axes are dimensional and check names
  size_t parametersIndex = 0;
  std::vector<size_t> dimensionIndex(numDimsTemp + 1, 3); // stores h, k, l or Qx, Qy, Qz dimensions
  for (const auto &p : parameters) {
    auto key = p.first;
    auto value = p.second;
    // value starts with QDimension0, then other stuff
    // do not use ==
    if (value.find("QDimension0") != std::string::npos) {
      dimensionIndex[0] = parametersIndex;
      const std::string dimXName = tempDataWS->getDimension(parametersIndex)->getName();
      if (m_isRLU) { // hkl
        if (dimXName != QDimensionName(m_Q0Basis)) {
          std::stringstream errorMessage;
          std::stringstream debugMessage;
          errorMessage << "TemporaryDataWorkspace does not have the  ";
          errorMessage << "correct name for dimension " << parametersIndex;
          debugMessage << "QDimension0 Names: Output will be: " << QDimensionName(m_Q0Basis);
          debugMessage << " TemporaryDataWorkspace: " << dimXName;
          g_log.warning(debugMessage.str());
          throw(std::invalid_argument(errorMessage.str()));
        }
      } else {
        if (dimXName != QDimensionNameQSample(0)) {
          std::stringstream errorMessage;
          std::stringstream debugMessage;
          errorMessage << "TemporaryDataWorkspace does not have the  ";
          errorMessage << "correct name for dimension " << parametersIndex;
          debugMessage << "QDimension0 Names: Output will be: " << QDimensionNameQSample(0);
          debugMessage << " TemporaryDataWorkspace: " << dimXName;
          g_log.warning(debugMessage.str());
          throw(std::invalid_argument(errorMessage.str()));
        }
      }
    } else if (value.find("QDimension1") != std::string::npos) {
      dimensionIndex[1] = parametersIndex;
      const std::string dimYName = tempDataWS->getDimension(parametersIndex)->getName();
      if (m_isRLU) { // hkl
        if (dimYName != QDimensionName(m_Q1Basis)) {
          std::stringstream errorMessage;
          std::stringstream debugMessage;
          errorMessage << "TemporaryDataWorkspace does not have the  ";
          errorMessage << "correct name for dimension " << parametersIndex;
          debugMessage << "QDimension1 Names: Output will be: " << QDimensionName(m_Q1Basis);
          debugMessage << " TemporaryDataWorkspace: " << dimYName;
          g_log.warning(debugMessage.str());
          throw(std::invalid_argument(errorMessage.str()));
        }
      } else {
        if (dimYName != QDimensionNameQSample(1)) {
          std::stringstream errorMessage;
          std::stringstream debugMessage;
          errorMessage << "TemporaryDataWorkspace does not have the  ";
          errorMessage << "correct name for dimension " << parametersIndex;
          debugMessage << "QDimension1 Names: Output will be: " << QDimensionNameQSample(1);
          debugMessage << " TemporaryDataWorkspace: " << dimYName;
          g_log.warning(debugMessage.str());
          throw(std::invalid_argument(errorMessage.str()));
        }
      }
    } else if (value.find("QDimension2") != std::string::npos) {
      dimensionIndex[2] = parametersIndex;
      const std::string dimZName = tempDataWS->getDimension(parametersIndex)->getName();
      if (m_isRLU) { // hkl
        if (dimZName != QDimensionName(m_Q2Basis)) {
          std::stringstream errorMessage;
          std::stringstream debugMessage;
          errorMessage << "TemporaryDataWorkspace does not have the  ";
          errorMessage << "correct name for dimension " << parametersIndex;
          debugMessage << "QDimension2 Names: Output will be: " << QDimensionName(m_Q2Basis);
          debugMessage << " TemporaryDataWorkspace: " << dimZName;
          g_log.warning(debugMessage.str());
          throw(std::invalid_argument(errorMessage.str()));
        }
      } else {
        if (dimZName != QDimensionNameQSample(2)) {
          std::stringstream errorMessage;
          std::stringstream debugMessage;
          errorMessage << "TemporaryDataWorkspace does not have the  ";
          errorMessage << "correct name for dimension " << parametersIndex;
          debugMessage << "QDimension2 Names: Output will be: " << QDimensionNameQSample(2);
          debugMessage << " TemporaryDataWorkspace: " << dimZName;
          g_log.warning(debugMessage.str());
          throw(std::invalid_argument(errorMessage.str()));
        }
      }
 
    } else if ((key != "OutputBins") && (key != "OutputExtents")) {
      // make sure the names of non-directional dimensions are the same
      const std::string nameData = tempDataWS->getDimension(parametersIndex)->getName();
      if (value.find(nameData) != 0) {
        g_log.error() << "Dimension " << nameData
                      << " from the temporary workspace"
                         " is not one of the binning dimensions, "
                         " or dimensions are in the wrong order."
                      << std::endl;
        throw(std::invalid_argument("Beside the Q dimensions, "
                                    "TemporaryDataWorkspace does not have the "
                                    "same dimension names as OutputWorkspace."));
      }
    }
    parametersIndex++;
  }
  const auto it = std::find_if(dimensionIndex.cbegin(), dimensionIndex.cend(),
                               [numDimsTemp](const auto &idx) { return idx > numDimsTemp; });
  if (it != dimensionIndex.cend())
    throw(std::invalid_argument("Cannot find at least one of QDimension0, "
                                "QDimension1, or QDimension2"));
}
 
inline DblMatrix MDNorm::buildSymmetryMatrix(const Geometry::SymmetryOperation &so) {
  // 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);
  soMatrix.Invert();
 
  return soMatrix;
}
 
// projection: input/output
// requiring: m_hIdx, m_kIndex, m_lIdx, meidx, m_dEintegrated, m_Q0Basis,
// mQ1Basis,
inline void MDNorm::determineBasisVector(const size_t &qindex, const std::string &value,
                                         const Mantid::Kernel::DblMatrix &Qtransform, std::vector<double> &projection,
                                         std::stringstream &basisVector, std::vector<size_t> &qDimensionIndices) {
  if (value.find("QDimension0") != std::string::npos) {
    m_hIdx = qindex;
    if (!m_isRLU) {
      projection[0] = 1.;
      basisVector << QDimensionNameQSample(0) << ",A^{-1}";
    } else {
      qDimensionIndices.emplace_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 << QDimensionNameQSample(1) << ",A^{-1}";
    } else {
      qDimensionIndices.emplace_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 << QDimensionNameQSample(2) << ",A^{-1}";
    } else {
      qDimensionIndices.emplace_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;
  }
}
 
inline void MDNorm::setQUnit(const std::vector<size_t> &qDimensionIndices,
                             const Mantid::DataObjects::MDHistoWorkspace_sptr &outputMDHWS) {
  Mantid::Geometry::MDFrameArgument argument(Mantid::Geometry::HKL::HKLName, Mantid::Kernel::Units::Symbol::RLU);
  auto mdFrameFactory = Mantid::Geometry::makeMDFrameFactoryChain();
  Mantid::Geometry::MDFrame_uptr hklFrame = mdFrameFactory->create(argument);
  for (size_t i : qDimensionIndices) {
    auto mdHistoDimension = std::const_pointer_cast<Mantid::Geometry::MDHistoDimension>(
        std::dynamic_pointer_cast<const Mantid::Geometry::MDHistoDimension>(outputMDHWS->getDimension(i)));
    mdHistoDimension->setMDFrame(*hklFrame);
  }
  // add W_matrix
  auto ei = outputMDHWS->getExperimentInfo(0);
  ei->mutableRun().addProperty("W_MATRIX", m_W.getVector(), true);
}
 
DataObjects::MDHistoWorkspace_sptr
MDNorm::binBackgroundWS(const std::vector<Geometry::SymmetryOperation> &symmetryOps) {
  // Create output background data histogram MD workspace
  // Either from TemporaryBackgroundDataWorkspace
  // Or from scratch
  Mantid::API::IMDHistoWorkspace_sptr tempBkgdDataWS = this->getProperty("TemporaryBackgroundDataWorkspace");
  Mantid::API::Workspace_sptr outputWS;
 
  // check that our input matches the temporary workspaces
  std::map<std::string, std::string> parameters = getBinParameters();
  if (tempBkgdDataWS) {
    validateBinningForTemporaryDataWorkspace(parameters, tempBkgdDataWS);
  }
 
  // For each symmetry operation, do binning MD once
  std::vector<size_t> qDimensionIndices;
  uint16_t numexpinfo = static_cast<uint16_t>(m_inputWS->getNumExperimentInfo());
  if (m_numSymmOps != symmetryOps.size())
    throw std::runtime_error("Symmetry operation number m_umSymops is wrong!");
 
  for (uint16_t i_expinfo = 0; i_expinfo < numexpinfo; ++i_expinfo) {
 
    auto rotMatrix = m_inputWS->getExperimentInfo(i_expinfo)->run().getGoniometerMatrix();
 
    // Reset symmetry operation index
    double soIndex = 0;
 
    for (const auto &so : symmetryOps) {
      // Q transformation matrix: From Q_lab to HKL or Q_sample
      // Building symmetric operation matrix
      DblMatrix soMatrix = buildSymmetryMatrix(so);
      // Calculate Q transform matrix
      DblMatrix Qtransform;
      if (m_isRLU) {
        Qtransform = rotMatrix * m_UB * soMatrix * m_W;
      } else {
        Qtransform = rotMatrix * soMatrix * m_W;
      }
 
      // Set up BinMD for this symmetry opeation
      double progress_fraction = 1. / static_cast<double>(symmetryOps.size() * numexpinfo);
      auto binMD =
          createChildAlgorithm("BinMD", soIndex * 0.3 * progress_fraction, (soIndex + 1) * 0.3 * progress_fraction);
 
      binMD->setPropertyValue("AxisAligned", "0");
      binMD->setProperty("InputWorkspace", m_backgroundWS);
      binMD->setProperty("TemporaryDataWorkspace", tempBkgdDataWS);
      binMD->setPropertyValue("NormalizeBasisVectors", "0");
      // Set the output Workspace directly to Algorithm's
      // OutputBackgroundDataWorkspace
      binMD->setPropertyValue("OutputWorkspace", getPropertyValue("OutputBackgroundDataWorkspace"));
      // set binning properties
      size_t qindex = 0;
      for (const auto &p : parameters) {
        auto key = p.first;
        auto value = p.second;
        std::stringstream basisVector;
        std::vector<double> projection(m_inputWS->getNumDims(), 0.);
        // value is a string that can start with QDimension0, etc, but contain
        // other stuff. Do not use ==
        determineBasisVector(qindex, value, Qtransform, projection, basisVector, qDimensionIndices);
 
        if (!basisVector.str().empty()) {
          // reconstruct from calculated basis vector
          for (auto proji : projection) {
            proji = std::abs(proji) > 1e-10 ? proji : 0.0;
            basisVector << "," << proji;
          }
          value = basisVector.str();
        }
 
        binMD->setPropertyValue(key, value);
        qindex++;
      }
      // execute algorithm
      binMD->executeAsChildAlg();
 
      // set the temporary workspace to be the output workspace, so it keeps
      // adding different symmetries AND
      // FIXME in future another ExpInfo
      outputWS = binMD->getProperty("OutputWorkspace");
      tempBkgdDataWS = std::dynamic_pointer_cast<MDHistoWorkspace>(outputWS);
      tempBkgdDataWS->clearOriginalWorkspaces();
      tempBkgdDataWS->clearTransforms();
    }
  }
  auto outputMDHWS = std::dynamic_pointer_cast<MDHistoWorkspace>(outputWS);
  // set MDUnits for Q dimensions
  if (m_isRLU) {
    setQUnit(qDimensionIndices, outputMDHWS);
  }
 
  outputMDHWS->setDisplayNormalization(Mantid::API::NoNormalization);
  return outputMDHWS;
}
 
DataObjects::MDHistoWorkspace_sptr MDNorm::binInputWS(const std::vector<Geometry::SymmetryOperation> &symmetryOps) {
  std::map<std::string, std::string> parameters = getBinParameters();
  return binMDEventWorkspace(m_inputWS, "TemporaryDataWorkspace", "OutputDataWorkspace", symmetryOps, parameters);
}
 
DataObjects::MDHistoWorkspace_sptr
MDNorm::binMonoSCDNormalizationWS(const std::vector<Geometry::SymmetryOperation> &symmetryOps) {
  std::map<std::string, std::string> parameters = getBinParameters();
  Mantid::API::IMDEventWorkspace_sptr monoNormInputWS = this->getProperty("MonoSCDNormalizationWorkspace");
  return binMDEventWorkspace(monoNormInputWS, "TemporaryNormalizationWorkspace", "OutputNormalizationWorkspace",
                             symmetryOps, parameters);
}
 
DataObjects::MDHistoWorkspace_sptr
MDNorm::binMDEventWorkspace(const API::IMDEventWorkspace_sptr &ws, const std::string &temporaryWSPropertyName,
                            const std::string &outputWSPropertyName,
                            const std::vector<Geometry::SymmetryOperation> &symmetryOps,
                            const std::map<std::string, std::string> &parameters) {
  Mantid::API::IMDHistoWorkspace_sptr tempWS = this->getProperty(temporaryWSPropertyName);
  Mantid::API::Workspace_sptr outputWS;
 
  // check that our input matches the temporary workspace
  if (tempWS)
    validateBinningForTemporaryDataWorkspace(parameters, tempWS);
 
  double soIndex = 0;
  std::vector<size_t> qDimensionIndices;
  for (const auto &so : symmetryOps) {
    // calculate dimensions for binning
    DblMatrix soMatrix = buildSymmetryMatrix(so);
    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());
    auto binMD = createChildAlgorithm("BinMD", soIndex * 0.3 * fraction, (soIndex + 1) * 0.3 * fraction);
    binMD->setPropertyValue("AxisAligned", "0");
    binMD->setProperty("InputWorkspace", ws);
    binMD->setProperty("TemporaryDataWorkspace", tempWS);
    binMD->setPropertyValue("NormalizeBasisVectors", "0");
    binMD->setPropertyValue("OutputWorkspace", getPropertyValue(outputWSPropertyName));
    // set binning properties
    size_t qindex = 0;
    for (const auto &p : parameters) {
      auto value = p.second;
      std::stringstream basisVector;
      std::vector<double> projection(ws->getNumDims(), 0.);
      // value is a string that can start with QDimension0, etc, but contain
      // other stuff. Do not use ==
      determineBasisVector(qindex, value, Qtransform, projection, basisVector, qDimensionIndices);
 
      if (!basisVector.str().empty()) {
        // reconstruct from calculated basis vector
        for (auto proji : projection) {
          proji = std::abs(proji) > 1e-10 ? proji : 0.0;
          basisVector << "," << proji;
        }
        value = basisVector.str();
      }
 
      binMD->setPropertyValue(p.first, 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
    tempWS = std::dynamic_pointer_cast<MDHistoWorkspace>(outputWS);
    tempWS->clearOriginalWorkspaces();
    tempWS->clearTransforms();
    soIndex += 1;
  }
 
  auto outputMDHWS = std::dynamic_pointer_cast<MDHistoWorkspace>(outputWS);
  // set MDUnits for Q dimensions
  if (m_isRLU) {
    setQUnit(qDimensionIndices, outputMDHWS);
  }
 
  outputMDHWS->setDisplayNormalization(Mantid::API::NoNormalization);
  return outputMDHWS;
}
 
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);
    auto inputDimMin = static_cast<float>(dimension->getMinimum());
    auto 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.emplace_back(value);
      if (value < inputDimMin || value > inputDimMax) {
        skipNormalization = true;
      }
      if ((!isIntegrated) && (value < outputDimMin || value > outputDimMax)) {
        skipNormalization = true;
      }
    }
  }
  return otherDimValues;
}
 
} // namespace Mantid::MDAlgorithms