ConvertSpectrumAxis2.cpp

//----------------------------------------------------------------------
// Includes
//----------------------------------------------------------------------
#include "MantidAlgorithms/ConvertSpectrumAxis2.h"
#include "MantidAPI/NumericAxis.h"
#include "MantidKernel/UnitConversion.h"
#include "MantidKernel/UnitFactory.h"
#include "MantidAPI/WorkspaceValidators.h"
#include "MantidAPI/Run.h"
#include <boost/function.hpp>
#include <boost/bind.hpp>
#include <cfloat>
#include "MantidKernel/BoundedValidator.h"
#include "MantidKernel/ListValidator.h"

namespace Mantid {
namespace Algorithms {
// Register the algorithm into the AlgorithmFactory
DECLARE_ALGORITHM(ConvertSpectrumAxis2)
using namespace Kernel;
using namespace API;
using namespace Geometry;

ConvertSpectrumAxis2::ConvertSpectrumAxis2()
    : API::Algorithm(), m_inputWS(), m_indexMap(), m_nBins(0), m_nxBins(0),
      m_nHist(0) {}

void ConvertSpectrumAxis2::init() {
  // Validator for Input Workspace
  auto wsVal = boost::make_shared<CompositeValidator>();
  wsVal->add<HistogramValidator>();
  wsVal->add<SpectraAxisValidator>();
  wsVal->add<InstrumentValidator>();

  declareProperty(
      new WorkspaceProperty<>("InputWorkspace", "", Direction::Input, wsVal),
      "The name of the input workspace.");
  declareProperty(
      new WorkspaceProperty<>("OutputWorkspace", "", Direction::Output),
      "The name to use for the output workspace.");
  std::vector<std::string> targetOptions(6);
  targetOptions[0] = "Theta";
  targetOptions[1] = "SignedTheta";
  targetOptions[2] = "ElasticQ";
  targetOptions[3] = "ElasticQSquared";
  targetOptions[4] = "theta";
  targetOptions[5] = "signed_theta";

  declareProperty(
      "Target", "", boost::make_shared<StringListValidator>(targetOptions),
      "The unit to which spectrum axis is converted to - \"theta\" (for the "
      "angle in degrees), Q or Q^2, where elastic Q is evaluated at EFixed. "
      "Note that 'theta' and 'signed_theta' are there for compatibility "
      "purposes; they are the same as 'Theta' and 'SignedTheta' respectively");
  std::vector<std::string> eModeOptions;
  eModeOptions.push_back("Direct");
  eModeOptions.push_back("Indirect");
  declareProperty("EMode", "Direct",
                  boost::make_shared<StringListValidator>(eModeOptions),
                  "Some unit conversions require this value to be set "
                  "(\"Direct\" or \"Indirect\")");
  auto mustBePositive = boost::make_shared<BoundedValidator<double>>();
  mustBePositive->setLower(0.0);
  declareProperty("EFixed", EMPTY_DBL(), mustBePositive,
                  "Value of fixed energy in meV : EI (EMode=Direct) or EF "
                  "(EMode=Indirect))");
}

void ConvertSpectrumAxis2::exec() {
  // Get the input workspace.
  m_inputWS = getProperty("InputWorkspace");
  // Assign value to the member variable storing the number of histograms.
  m_nHist = m_inputWS->getNumberHistograms();

  // Assign values to the member variables to store number of bins.
  m_nBins = m_inputWS->blocksize();

  const bool isHist = m_inputWS->isHistogramData();

  if (isHist) {
    m_nxBins = m_nBins + 1;
  } else {
    m_nxBins = m_nBins;
  }

  // The unit to convert to.
  const std::string unitTarget = getProperty("Target");

  // Call the functions to convert to the different forms of theta or Q.
  if (unitTarget == "theta" || unitTarget == "Theta" ||
      unitTarget == "signed_theta" || unitTarget == "SignedTheta") {
    createThetaMap(unitTarget);
  } else if (unitTarget == "ElasticQ" || unitTarget == "ElasticQSquared") {
    createElasticQMap(unitTarget);
  }

  // Create an output workspace and set the property for it.
  MatrixWorkspace_sptr outputWS = createOutputWorkspace(unitTarget);
  setProperty("OutputWorkspace", outputWS);
}

void ConvertSpectrumAxis2::createThetaMap(const std::string &targetUnit) {
  // Set up binding to member funtion. Avoids condition as part of loop over
  // nHistograms.
  boost::function<double(IDetector_const_sptr)> thetaFunction;
  if (targetUnit.compare("signed_theta") == 0 ||
      targetUnit.compare("SignedTheta") == 0) {
    thetaFunction =
        boost::bind(&MatrixWorkspace::detectorSignedTwoTheta, m_inputWS, _1);
  } else if (targetUnit == "theta" || targetUnit == "Theta") {
    thetaFunction =
        boost::bind(&MatrixWorkspace::detectorTwoTheta, m_inputWS, _1);
  }

  bool warningGiven = false;
  for (size_t i = 0; i < m_nHist; ++i) {
    try {
      IDetector_const_sptr det = m_inputWS->getDetector(i);
      // Invoke relevant member function.
      m_indexMap.insert(std::make_pair(thetaFunction(det) * 180.0 / M_PI, i));
    } catch (Exception::NotFoundError &) {
      if (!warningGiven)
        g_log.warning("The instrument definition is incomplete - spectra "
                      "dropped from output");
      warningGiven = true;
    }
  }
}

void ConvertSpectrumAxis2::createElasticQMap(const std::string &targetUnit) {
  IComponent_const_sptr source = m_inputWS->getInstrument()->getSource();
  IComponent_const_sptr sample = m_inputWS->getInstrument()->getSample();

  const std::string emodeStr = getProperty("EMode");
  int emode = 0;
  if (emodeStr == "Direct")
    emode = 1;
  else if (emodeStr == "Indirect")
    emode = 2;

  for (size_t i = 0; i < m_nHist; i++) {
    IDetector_const_sptr detector = m_inputWS->getDetector(i);
    double twoTheta(0.0), efixed(0.0);
    if (!detector->isMonitor()) {
      twoTheta = m_inputWS->detectorTwoTheta(detector) / 2.0;
      efixed = getEfixed(detector, m_inputWS, emode); // get efixed
    } else {
      twoTheta = 0.0;
      efixed = DBL_MIN;
    }

    // Convert to MomentumTransfer
    double elasticQInAngstroms = Kernel::UnitConversion::run(twoTheta, efixed);

    if (targetUnit == "ElasticQ") {
      m_indexMap.insert(std::make_pair(elasticQInAngstroms, i));
    } else if (targetUnit == "ElasticQSquared") {
      // The QSquared value.
      double elasticQSquaredInAngstroms =
          elasticQInAngstroms * elasticQInAngstroms;

      m_indexMap.insert(std::make_pair(elasticQSquaredInAngstroms, i));
    }
  }
}

MatrixWorkspace_sptr
ConvertSpectrumAxis2::createOutputWorkspace(const std::string &targetUnit) {
  // Create the output workspace. Can not re-use the input one because the
  // spectra are re-ordered.
  MatrixWorkspace_sptr outputWorkspace = WorkspaceFactory::Instance().create(
      m_inputWS, m_indexMap.size(), m_nxBins, m_nBins);

  // Now set up a new numeric axis holding the theta values corresponding to
  // each spectrum.
  NumericAxis *const newAxis = new NumericAxis(m_indexMap.size());
  outputWorkspace->replaceAxis(1, newAxis);

  // Set the units of the axis.
  if (targetUnit == "theta" || targetUnit == "Theta" ||
      targetUnit == "signed_theta" || targetUnit == "SignedTheta") {
    newAxis->unit() = boost::shared_ptr<Unit>(new Units::Degrees);
  } else if (targetUnit == "ElasticQ") {
    newAxis->unit() = UnitFactory::Instance().create("MomentumTransfer");
  } else if (targetUnit == "ElasticQSquared") {
    newAxis->unit() = UnitFactory::Instance().create("QSquared");
  }

  std::multimap<double, size_t>::const_iterator it;
  size_t currentIndex = 0;
  for (it = m_indexMap.begin(); it != m_indexMap.end(); ++it) {
    // Set the axis value.
    newAxis->setValue(currentIndex, it->first);
    // Copy over the data.
    outputWorkspace->dataX(currentIndex) = m_inputWS->dataX(it->second);
    outputWorkspace->dataY(currentIndex) = m_inputWS->dataY(it->second);
    outputWorkspace->dataE(currentIndex) = m_inputWS->dataE(it->second);
    // We can keep the spectrum numbers etc.
    outputWorkspace->getSpectrum(currentIndex)
        ->copyInfoFrom(*m_inputWS->getSpectrum(it->second));
    ++currentIndex;
  }
  return outputWorkspace;
}

double ConvertSpectrumAxis2::getEfixed(IDetector_const_sptr detector,
                                       MatrixWorkspace_const_sptr inputWS,
                                       int emode) const {
  double efixed(0);
  double efixedProp = getProperty("Efixed");
  if (efixedProp != EMPTY_DBL()) {
    efixed = efixedProp;
    g_log.debug() << "Detector: " << detector->getID() << " Efixed: " << efixed
                  << "\n";
  } else {
    if (emode == 1) {
      if (inputWS->run().hasProperty("Ei")) {
        efixed = inputWS->run().getLogAsSingleValue("Ei");
      } else {
        throw std::invalid_argument("Could not retrieve Efixed from the "
                                    "workspace. Please provide a value.");
      }
    } else if (emode == 2) {
      std::vector<double> efixedVec = detector->getNumberParameter("Efixed");
      if (efixedVec.empty()) {
        int detid = detector->getID();
        IDetector_const_sptr detectorSingle =
            inputWS->getInstrument()->getDetector(detid);
        efixedVec = detectorSingle->getNumberParameter("Efixed");
      }
      if (!efixedVec.empty()) {
        efixed = efixedVec.at(0);
        g_log.debug() << "Detector: " << detector->getID()
                      << " EFixed: " << efixed << "\n";
      } else {
        g_log.warning() << "Efixed could not be found for detector "
                        << detector->getID() << ", please provide a value\n";
        throw std::invalid_argument("Could not retrieve Efixed from the "
                                    "detector. Please provide a value.");
      }
    }
  }
  return efixed;
}

} // namespace Algorithms
} // namespace Mantid