MatrixWorkspace.cpp

#include "MantidAPI/MatrixWorkspace.h"
#include "MantidAPI/BinEdgeAxis.h"
#include "MantidAPI/NumericAxis.h"
#include "MantidAPI/SpectraAxis.h"
#include "MantidAPI/MatrixWorkspaceMDIterator.h"
#include "MantidAPI/SpectrumDetectorMapping.h"
#include "MantidGeometry/Instrument/Detector.h"
#include "MantidGeometry/Instrument/DetectorGroup.h"
#include "MantidGeometry/Instrument/NearestNeighboursFactory.h"
#include "MantidGeometry/Instrument/ReferenceFrame.h"
#include "MantidGeometry/MDGeometry/MDFrame.h"
#include "MantidGeometry/MDGeometry/GeneralFrame.h"
#include "MantidKernel/TimeSeriesProperty.h"
#include "MantidKernel/MDUnit.h"

#include <numeric>
#include <boost/math/special_functions/fpclassify.hpp>

using Mantid::Kernel::DateAndTime;
using Mantid::Kernel::TimeSeriesProperty;
using NeXus::NXcompression;
using Mantid::Kernel::Strings::toString;

namespace Mantid {
namespace API {
using std::size_t;
using namespace Geometry;
using Kernel::V3D;

namespace {
/// static logger
Kernel::Logger g_log("MatrixWorkspace");
}

const std::string MatrixWorkspace::xDimensionId = "xDimension";
const std::string MatrixWorkspace::yDimensionId = "yDimension";

/// Default constructor
MatrixWorkspace::MatrixWorkspace(
    Mantid::Geometry::INearestNeighboursFactory *nnFactory)
    : IMDWorkspace(), ExperimentInfo(), m_axes(), m_isInitialized(false),
      m_YUnit(), m_YUnitLabel(), m_isDistribution(false),
      m_isCommonBinsFlagSet(false), m_isCommonBinsFlag(false), m_masks(),
      m_indexCalculator(),
      m_nearestNeighboursFactory(
          (nnFactory == NULL) ? new NearestNeighboursFactory : nnFactory),
      m_nearestNeighbours() {}

MatrixWorkspace::MatrixWorkspace(const MatrixWorkspace &other)
    : IMDWorkspace(other), ExperimentInfo(other) {
  m_axes.resize(other.m_axes.size());
  for (size_t i = 0; i < m_axes.size(); ++i)
    m_axes[i] = other.m_axes[i]->clone(this);

  m_isInitialized = other.m_isInitialized;
  m_YUnit = other.m_YUnit;
  m_YUnitLabel = other.m_YUnitLabel;
  m_isDistribution = other.m_isDistribution;
  m_isCommonBinsFlagSet = other.m_isCommonBinsFlagSet;
  m_isCommonBinsFlag = other.m_isCommonBinsFlag;
  m_masks = other.m_masks;
  m_indexCalculator = other.m_indexCalculator;
  // I think it is necessary to create our own copy of the factory, since we do
  // not know who owns the factory in other and how its lifetime is controlled.
  m_nearestNeighboursFactory.reset(new NearestNeighboursFactory);
  // m_nearestNeighbours seem to be built automatically when needed, so we do
  // not copy here.

  // TODO: Do we need to init m_monitorWorkspace?

  // This call causes copying of m_parmap (ParameterMap). The constructor of
  // ExperimentInfo just kept a shared_ptr to the same map as in other, which
  // is not enough as soon as the maps in one of the workspaces it edited.
  instrumentParameters();
}

/// Destructor
// RJT, 3/10/07: The Analysis Data Service needs to be able to delete
// workspaces, so I moved this from protected to public.
MatrixWorkspace::~MatrixWorkspace() {
  for (unsigned int i = 0; i < m_axes.size(); ++i) {
    delete m_axes[i];
  }
}

/// @returns A human-readable string of the current state
const std::string MatrixWorkspace::toString() const {
  std::ostringstream os;
  os << id() << "\n"
     << "Title: " << getTitle() << "\n"
     << "Histograms: " << getNumberHistograms() << "\n"
     << "Bins: " << blocksize() << "\n";

  if (isHistogramData())
    os << "Histogram\n";
  else
    os << "Data points\n";

  os << "X axis: ";
  if (axes() > 0) {
    Axis *ax = getAxis(0);
    if (ax && ax->unit())
      os << ax->unit()->caption() << " / " << ax->unit()->label().ascii();
    else
      os << "Not set";
  } else {
    os << "N/A";
  }
  os << "\n"
     << "Y axis: " << YUnitLabel() << "\n";

  os << "Distribution: " << (isDistribution() ? "True" : "False") << "\n";

  os << ExperimentInfo::toString();
  return os.str();
}

/** Initialize the workspace. Calls the protected init() method, which is
* implemented in each type of
*  workspace. Returns immediately if the workspace is already initialized.
*  @param NVectors :: The number of spectra in the workspace (only relevant for
* a 2D workspace
*  @param XLength :: The number of X data points/bin boundaries in each vector
* (must all be the same)
*  @param YLength :: The number of data/error points in each vector (must all be
* the same)
*/
void MatrixWorkspace::initialize(const std::size_t &NVectors,
                                 const std::size_t &XLength,
                                 const std::size_t &YLength) {
  // Check validity of arguments
  if (NVectors == 0 || XLength == 0 || YLength == 0) {
    throw std::out_of_range(
        "All arguments to init must be positive and non-zero");
  }

  // Bypass the initialization if the workspace has already been initialized.
  if (m_isInitialized)
    return;

  // Invoke init() method of the derived class inside a try/catch clause
  try {
    this->init(NVectors, XLength, YLength);
  } catch (std::runtime_error &) {
    throw;
  }

  m_indexCalculator = MatrixWSIndexCalculator(this->blocksize());
  // Indicate that this workspace has been initialized to prevent duplicate
  // attempts.
  m_isInitialized = true;
}

//---------------------------------------------------------------------------------------
/** Set the title of the workspace
*
*  @param t :: The title
*/
void MatrixWorkspace::setTitle(const std::string &t) {
  Workspace::setTitle(t);

  // A MatrixWorkspace contains uniquely one Run object, hence for this
  // workspace
  // keep the Run object run_title property the same as the workspace title
  Run &run = mutableRun();
  run.addProperty("run_title", t, true);
}

//---------------------------------------------------------------------------------------
/** Get the workspace title
*
*  @return The title
*/
const std::string MatrixWorkspace::getTitle() const {
  if (run().hasProperty("run_title")) {
    std::string title = run().getProperty("run_title")->value();
    return title;
  } else
    return Workspace::getTitle();
}

void MatrixWorkspace::updateSpectraUsing(const SpectrumDetectorMapping &map) {
  for (size_t j = 0; j < getNumberHistograms(); ++j) {
    auto spec = getSpectrum(j);
    try {
      if (map.indexIsSpecNumber())
        spec->setDetectorIDs(
            map.getDetectorIDsForSpectrumNo(spec->getSpectrumNo()));
      else
        spec->setDetectorIDs(map.getDetectorIDsForSpectrumIndex(j));
    } catch (std::out_of_range &e) {
      // Get here if the spectrum number is not in the map.
      spec->clearDetectorIDs();
      g_log.debug(e.what());
      g_log.debug() << "Spectrum number " << spec->getSpectrumNo()
                    << " not in map.\n";
    }
  }
}

//---------------------------------------------------------------------------------------
/**
* Rebuild the default spectra mapping for a workspace. If a non-empty
* instrument is set then the default maps each detector to a spectra with
* the same ID. If an empty instrument is set then a 1:1 map from 1->NHistograms
* is created.
* @param includeMonitors :: If false the monitors are not included
*/
void MatrixWorkspace::rebuildSpectraMapping(const bool includeMonitors) {
  if (sptr_instrument->nelements() == 0) {
    return;
  }

  std::vector<detid_t> pixelIDs =
      this->getInstrument()->getDetectorIDs(!includeMonitors);

  try {
    size_t index = 0;
    std::vector<detid_t>::const_iterator iend = pixelIDs.end();
    for (std::vector<detid_t>::const_iterator it = pixelIDs.begin(); it != iend;
         ++it) {
      // The detector ID
      const detid_t detId = *it;
      // By default: Spectrum number = index +  1
      const specid_t specNo = specid_t(index + 1);

      if (index < this->getNumberHistograms()) {
        ISpectrum *spec = getSpectrum(index);
        spec->setSpectrumNo(specNo);
        spec->setDetectorID(detId);
      }

      index++;
    }

    m_nearestNeighbours.reset();

  } catch (std::runtime_error &) {
    throw;
  }
}

//---------------------------------------------------------------------------------------
/**
* Handles the building of the NearestNeighbours object, if it has not already
* been
* populated for this parameter map.
* @param ignoreMaskedDetectors :: flag indicating that masked detectors should
* be ignored. True to ignore detectors.
*/
void MatrixWorkspace::buildNearestNeighbours(
    const bool ignoreMaskedDetectors) const {
  if (!m_nearestNeighbours) {
    boost::shared_ptr<const Instrument> inst = this->getInstrument();
    if (inst) {
      SpectrumDetectorMapping spectraMap(this);
      m_nearestNeighbours.reset(m_nearestNeighboursFactory->create(
          inst, spectraMap.getMapping(), ignoreMaskedDetectors));
    } else {
      throw Mantid::Kernel::Exception::NullPointerException(
          "ParameterMap: buildNearestNeighbours. Can't obtain instrument.",
          "instrument");
    }
  }
}

/*
Allow the NearestNeighbours list to be cleaned and rebuilt. Certain algorithms
require this in order to exclude/include
detectors from previously being considered.
*/
void MatrixWorkspace::rebuildNearestNeighbours() {
  /*m_nearestNeighbours should now be NULL. This will trigger rebuilding on
  subsequent first call to getNeighbours
  ,which peforms a lazy evaluation on the nearest neighbours map */
  m_nearestNeighbours.reset();
}

//---------------------------------------------------------------------------------------
/** Queries the NearestNeighbours object for the selected detector.
* NOTE! getNeighbours(spectrumNumber, radius) is MUCH faster.
*
* @param comp :: pointer to the querying detector
* @param radius :: distance from detector on which to filter results
* @param ignoreMaskedDetectors :: flag indicating that masked detectors should
*be ignored. True to ignore detectors.
* @return map of DetectorID to distance for the nearest neighbours
*/
std::map<specid_t, V3D>
MatrixWorkspace::getNeighbours(const Geometry::IDetector *comp,
                               const double radius,
                               const bool ignoreMaskedDetectors) const {
  if (!m_nearestNeighbours) {
    buildNearestNeighbours(ignoreMaskedDetectors);
  }
  // Find the spectrum number
  std::vector<specid_t> spectra;
  this->getSpectraFromDetectorIDs(std::vector<detid_t>(1, comp->getID()),
                                  spectra);
  if (spectra.empty()) {
    throw Kernel::Exception::NotFoundError("MatrixWorkspace::getNeighbours - "
                                           "Cannot find spectrum number for "
                                           "detector",
                                           comp->getID());
  }
  std::map<specid_t, V3D> neighbours =
      m_nearestNeighbours->neighboursInRadius(spectra[0], radius);
  return neighbours;
}

//---------------------------------------------------------------------------------------
/** Queries the NearestNeighbours object for the selected spectrum number.
*
* @param spec :: spectrum number of the detector you are looking at
* @param radius :: distance from detector on which to filter results
* @param ignoreMaskedDetectors :: flag indicating that masked detectors should
*be ignored. True to ignore detectors.
* @return map of DetectorID to distance for the nearest neighbours
*/
std::map<specid_t, V3D>
MatrixWorkspace::getNeighbours(specid_t spec, const double radius,
                               bool ignoreMaskedDetectors) const {
  if (!m_nearestNeighbours) {
    buildNearestNeighbours(ignoreMaskedDetectors);
  }
  std::map<specid_t, V3D> neighbours =
      m_nearestNeighbours->neighboursInRadius(spec, radius);
  return neighbours;
}

//---------------------------------------------------------------------------------------
/** Queries the NearestNeighbours object for the selected spectrum number.
*
* @param spec :: spectrum number of the detector you are looking at
* @param nNeighbours :: unsigned int, number of neighbours to include.
* @param ignoreMaskedDetectors :: flag indicating that masked detectors should
*be ignored. True to ignore detectors.
* @return map of DetectorID to distance for the nearest neighbours
*/
std::map<specid_t, V3D>
MatrixWorkspace::getNeighboursExact(specid_t spec, const int nNeighbours,
                                    bool ignoreMaskedDetectors) const {
  if (!m_nearestNeighbours) {
    SpectrumDetectorMapping spectraMap(this);
    m_nearestNeighbours.reset(m_nearestNeighboursFactory->create(
        nNeighbours, this->getInstrument(), spectraMap.getMapping(),
        ignoreMaskedDetectors));
  }
  std::map<specid_t, V3D> neighbours = m_nearestNeighbours->neighbours(spec);
  return neighbours;
}

//---------------------------------------------------------------------------------------
/** Return a map where:
*    KEY is the Spectrum #
*    VALUE is the Workspace Index
*/
spec2index_map MatrixWorkspace::getSpectrumToWorkspaceIndexMap() const {
  SpectraAxis *ax = dynamic_cast<SpectraAxis *>(this->m_axes[1]);
  if (!ax)
    throw std::runtime_error("MatrixWorkspace::getSpectrumToWorkspaceIndexMap: "
                             "axis[1] is not a SpectraAxis, so I cannot "
                             "generate a map.");
  spec2index_map map;
  try {
    ax->getSpectraIndexMap(map);
  } catch (std::runtime_error &) {
    throw std::runtime_error(
        "MatrixWorkspace::getSpectrumToWorkspaceIndexMap: no elements!");
  }
  return map;
}

//---------------------------------------------------------------------------------------
/** Return a vector where:
*    The index into the vector = spectrum number + offset
*    The value at that index = the corresponding Workspace Index
*
*  @param out :: vector set to above definition
*  @param offset :: add this to the detector ID to get the index into the
*vector.
*/
void MatrixWorkspace::getSpectrumToWorkspaceIndexVector(
    std::vector<size_t> &out, specid_t &offset) const {
  SpectraAxis *ax = dynamic_cast<SpectraAxis *>(this->m_axes[1]);
  if (!ax)
    throw std::runtime_error("MatrixWorkspace::getSpectrumToWorkspaceIndexMap: "
                             "axis[1] is not a SpectraAxis, so I cannot "
                             "generate a map.");

  // Find the min/max spectra IDs
  specid_t min = std::numeric_limits<specid_t>::max(); // So that any number
                                                       // will be less than this
  specid_t max = -std::numeric_limits<specid_t>::max(); // So that any number
                                                        // will be greater than
                                                        // this
  size_t length = ax->length();
  for (size_t i = 0; i < length; i++) {
    specid_t spec = ax->spectraNo(i);
    if (spec < min)
      min = spec;
    if (spec > max)
      max = spec;
  }

  // Offset so that the "min" value goes to index 0
  offset = -min;

  // Resize correctly
  out.resize(max - min + 1, 0);

  // Make the vector
  for (size_t i = 0; i < length; i++) {
    specid_t spec = ax->spectraNo(i);
    out[spec + offset] = i;
  }
}

//---------------------------------------------------------------------------------------
/** Does the workspace has any grouped detectors?
*  @return true if the workspace has any grouped detectors, otherwise false
*/
bool MatrixWorkspace::hasGroupedDetectors() const {
  bool retVal = false;

  // Loop through the workspace index
  for (size_t workspaceIndex = 0; workspaceIndex < this->getNumberHistograms();
       workspaceIndex++) {
    auto detList = getSpectrum(workspaceIndex)->getDetectorIDs();
    if (detList.size() > 1) {
      retVal = true;
      break;
    }
  }
  return retVal;
}

//---------------------------------------------------------------------------------------
/** Return a map where:
*    KEY is the DetectorID (pixel ID)
*    VALUE is the Workspace Index
*  @param throwIfMultipleDets :: set to true to make the algorithm throw an
* error
*         if there is more than one detector for a specific workspace index.
*  @throw runtime_error if there is more than one detector per spectrum (if
* throwIfMultipleDets is true)
*  @return Index to Index Map object. THE CALLER TAKES OWNERSHIP OF THE MAP AND
* IS RESPONSIBLE FOR ITS DELETION.
*/
detid2index_map MatrixWorkspace::getDetectorIDToWorkspaceIndexMap(
    bool throwIfMultipleDets) const {
  detid2index_map map;

  // Loop through the workspace index
  for (size_t workspaceIndex = 0; workspaceIndex < this->getNumberHistograms();
       ++workspaceIndex) {
    auto detList = getSpectrum(workspaceIndex)->getDetectorIDs();

    if (throwIfMultipleDets) {
      if (detList.size() > 1) {
        throw std::runtime_error(
            "MatrixWorkspace::getDetectorIDToWorkspaceIndexMap(): more than 1 "
            "detector for one histogram! I cannot generate a map of detector "
            "ID to workspace index.");
      }

      // Set the KEY to the detector ID and the VALUE to the workspace index.
      if (detList.size() == 1)
        map[*detList.begin()] = workspaceIndex;
    } else {
      // Allow multiple detectors per workspace index
      for (auto it = detList.begin(); it != detList.end(); ++it)
        map[*it] = workspaceIndex;
    }

    // Ignore if the detector list is empty.
  }

  return map;
}

//---------------------------------------------------------------------------------------
/** Return a vector where:
*    The index into the vector = DetectorID (pixel ID) + offset
*    The value at that index = the corresponding Workspace Index
*
*  @param out :: vector set to above definition
*  @param offset :: add this to the detector ID to get the index into the
*vector.
*  @param throwIfMultipleDets :: set to true to make the algorithm throw an
*error
*         if there is more than one detector for a specific workspace index.
*  @throw runtime_error if there is more than one detector per spectrum (if
*throwIfMultipleDets is true)
*/
void MatrixWorkspace::getDetectorIDToWorkspaceIndexVector(
    std::vector<size_t> &out, detid_t &offset, bool throwIfMultipleDets) const {
  // Make a correct initial size
  out.clear();
  detid_t minId = 0;
  detid_t maxId = 0;
  this->getInstrument()->getMinMaxDetectorIDs(minId, maxId);
  offset = -minId;
  const int outSize = maxId - minId + 1;
  // Allocate at once
  out.resize(outSize, std::numeric_limits<size_t>::max());

  for (size_t workspaceIndex = 0; workspaceIndex < getNumberHistograms();
       ++workspaceIndex) {
    // Get the list of detectors from the WS index
    const std::set<detid_t> &detList =
        this->getSpectrum(workspaceIndex)->getDetectorIDs();

    if (throwIfMultipleDets && (detList.size() > 1))
      throw std::runtime_error(
          "MatrixWorkspace::getDetectorIDToWorkspaceIndexVector(): more than 1 "
          "detector for one histogram! I cannot generate a map of detector ID "
          "to workspace index.");

    // Allow multiple detectors per workspace index, or,
    // If only one is allowed, then this has thrown already
    for (std::set<detid_t>::const_iterator it = detList.begin();
         it != detList.end(); ++it) {
      int index = *it + offset;
      if (index < 0 || index >= outSize) {
        g_log.debug() << "MatrixWorkspace::getDetectorIDToWorkspaceIndexVector("
                         "): detector ID found (" << *it
                      << " at workspace index " << workspaceIndex
                      << ") is invalid." << std::endl;
      } else
        // Save it at that point.
        out[index] = workspaceIndex;
    }

  } // (for each workspace index)
}

//---------------------------------------------------------------------------------------
/** Converts a list of spectrum numbers to the corresponding workspace indices.
*  Not a very efficient operation, but unfortunately it's sometimes required.
*
*  @param spectraList :: The list of spectrum numbers required
*  @param indexList ::   Returns a reference to the vector of indices (empty if
*not a Workspace2D)
*/
void MatrixWorkspace::getIndicesFromSpectra(
    const std::vector<specid_t> &spectraList,
    std::vector<size_t> &indexList) const {
  // Clear the output index list
  indexList.clear();
  indexList.reserve(this->getNumberHistograms());

  std::vector<specid_t>::const_iterator iter = spectraList.begin();
  while (iter != spectraList.end()) {
    for (size_t i = 0; i < this->getNumberHistograms(); ++i) {
      if (this->getSpectrum(i)->getSpectrumNo() == *iter) {
        indexList.push_back(i);
        break;
      }
    }
    ++iter;
  }
}

//---------------------------------------------------------------------------------------
/** Given a spectrum number, find the corresponding workspace index
*
* @param specNo :: spectrum number wanted
* @return the workspace index
* @throw runtime_error if not found.
*/
size_t
MatrixWorkspace::getIndexFromSpectrumNumber(const specid_t specNo) const {
  for (size_t i = 0; i < this->getNumberHistograms(); ++i) {
    if (this->getSpectrum(i)->getSpectrumNo() == specNo)
      return i;
  }
  throw std::runtime_error("Could not find spectrum number in any spectrum.");
}

//---------------------------------------------------------------------------------------
/** Converts a list of detector IDs to the corresponding workspace indices.
*
     *  Note that only known detector IDs are converted (so an empty vector will
*be returned
     *  if none of the IDs are recognised), and that the returned workspace
*indices are
     *  effectively a set (i.e. there are no duplicates).
     *
*  @param detIdList :: The list of detector IDs required
*  @param indexList :: Returns a reference to the vector of indices
*/
void MatrixWorkspace::getIndicesFromDetectorIDs(
    const std::vector<detid_t> &detIdList,
    std::vector<size_t> &indexList) const {
  std::map<detid_t, std::set<size_t>> detectorIDtoWSIndices;
  for (size_t i = 0; i < getNumberHistograms(); ++i) {
    auto detIDs = getSpectrum(i)->getDetectorIDs();
    for (auto it = detIDs.begin(); it != detIDs.end(); ++it) {
      detectorIDtoWSIndices[*it].insert(i);
    }
  }

  indexList.clear();
  indexList.reserve(detIdList.size());
  for (size_t j = 0; j < detIdList.size(); ++j) {
    auto wsIndices = detectorIDtoWSIndices.find(detIdList[j]);
    if (wsIndices != detectorIDtoWSIndices.end()) {
      for (auto it = wsIndices->second.begin(); it != wsIndices->second.end();
           ++it) {
        indexList.push_back(*it);
      }
    }
  }
}

//---------------------------------------------------------------------------------------
/** Converts a list of detector IDs to the corresponding spectrum numbers. Might
*be slow!
*
* @param detIdList :: The list of detector IDs required
* @param spectraList :: Returns a reference to the vector of spectrum numbers.
*                       0 for not-found detectors
*/
void MatrixWorkspace::getSpectraFromDetectorIDs(
    const std::vector<detid_t> &detIdList,
    std::vector<specid_t> &spectraList) const {
  std::vector<detid_t>::const_iterator it_start = detIdList.begin();
  std::vector<detid_t>::const_iterator it_end = detIdList.end();

  spectraList.clear();

  // Try every detector in the list
  std::vector<detid_t>::const_iterator it;
  for (it = it_start; it != it_end; ++it) {
    bool foundDet = false;
    specid_t foundSpecNum = 0;

    // Go through every histogram
    for (size_t i = 0; i < this->getNumberHistograms(); i++) {
      if (this->getSpectrum(i)->hasDetectorID(*it)) {
        foundDet = true;
        foundSpecNum = this->getSpectrum(i)->getSpectrumNo();
        break;
      }
    }

    if (foundDet)
      spectraList.push_back(foundSpecNum);
  } // for each detector ID in the list
}

double MatrixWorkspace::getXMin() const {
  double xmin;
  double xmax;
  this->getXMinMax(xmin, xmax); // delegate to the proper code
  return xmin;
}

double MatrixWorkspace::getXMax() const {
  double xmin;
  double xmax;
  this->getXMinMax(xmin, xmax); // delegate to the proper code
  return xmax;
}

void MatrixWorkspace::getXMinMax(double &xmin, double &xmax) const {
  // set to crazy values to start
  xmin = std::numeric_limits<double>::max();
  xmax = -1.0 * xmin;
  size_t numberOfSpectra = this->getNumberHistograms();

  // determine the data range
  for (size_t workspaceIndex = 0; workspaceIndex < numberOfSpectra;
       workspaceIndex++) {
    const MantidVec &dataX = this->readX(workspaceIndex);
    const double xfront = dataX.front();
    const double xback = dataX.back();
    if (boost::math::isfinite(xfront) && boost::math::isfinite(xback)) {
      if (xfront < xmin)
        xmin = xfront;
      if (xback > xmax)
        xmax = xback;
    }
  }
}

//---------------------------------------------------------------------------------------
/** Integrate all the spectra in the matrix workspace within the range given.
* Default implementation, can be overridden by base classes if they know
*something smarter!
*
* @param out :: returns the vector where there is one entry per spectrum in the
*workspace. Same
*            order as the workspace indices.
* @param minX :: minimum X bin to use in integrating.
* @param maxX :: maximum X bin to use in integrating.
* @param entireRange :: set to true to use the entire range. minX and maxX are
*then ignored!
*/
void MatrixWorkspace::getIntegratedSpectra(std::vector<double> &out,
                                           const double minX, const double maxX,
                                           const bool entireRange) const {
  out.resize(this->getNumberHistograms(), 0.0);

  // Run in parallel if the implementation is threadsafe
  PARALLEL_FOR_IF(this->threadSafe())
  for (int wksp_index = 0;
       wksp_index < static_cast<int>(this->getNumberHistograms());
       wksp_index++) {
    // Get Handle to data
    const Mantid::MantidVec &x = this->readX(wksp_index);
    const Mantid::MantidVec &y = this->readY(wksp_index);
    // If it is a 1D workspace, no need to integrate
    if ((x.size() <= 2) && (y.size() >= 1)) {
      out[wksp_index] = y[0];
    } else {
      // Iterators for limits - whole range by default
      Mantid::MantidVec::const_iterator lowit, highit;
      lowit = x.begin();
      highit = x.end() - 1;

      // But maybe we don't want the entire range?
      if (!entireRange) {
        // If the first element is lower that the xmin then search for new lowit
        if ((*lowit) < minX)
          lowit = std::lower_bound(x.begin(), x.end(), minX);
        // If the last element is higher that the xmax then search for new lowit
        if ((*highit) > maxX)
          highit = std::upper_bound(lowit, x.end(), maxX);
      }

      // Get the range for the y vector
      Mantid::MantidVec::difference_type distmin =
          std::distance(x.begin(), lowit);
      Mantid::MantidVec::difference_type distmax =
          std::distance(x.begin(), highit);
      double sum(0.0);
      if (distmin <= distmax) {
        // Integrate
        sum = std::accumulate(y.begin() + distmin, y.begin() + distmax, 0.0);
      }
      // Save it in the vector
      out[wksp_index] = sum;
    }
  }
}

/** Get the effective detector for the given spectrum
*  @param  workspaceIndex The workspace index for which the detector is required
*  @return A single detector object representing the detector(s) contributing
*          to the given spectrum number. If more than one detector contributes
then
*          the returned object's concrete type will be DetectorGroup.
*  @throw  Kernel::Exception::NotFoundError If the Instrument is missing or the
requested workspace index does not have any associated detectors
*/
Geometry::IDetector_const_sptr
MatrixWorkspace::getDetector(const size_t workspaceIndex) const {
  const ISpectrum *spec = this->getSpectrum(workspaceIndex);
  if (!spec)
    throw Kernel::Exception::NotFoundError("MatrixWorkspace::getDetector(): "
                                           "NULL spectrum found at the given "
                                           "workspace index.",
                                           "");

  const std::set<detid_t> &dets = spec->getDetectorIDs();
  Instrument_const_sptr localInstrument = getInstrument();
  if (!localInstrument) {
    g_log.debug() << "No instrument defined.\n";
    throw Kernel::Exception::NotFoundError("Instrument not found", "");
  }

  const size_t ndets = dets.size();
  if (ndets == 1) {
    // If only 1 detector for the spectrum number, just return it
    return localInstrument->getDetector(*dets.begin());
  } else if (ndets == 0) {
    throw Kernel::Exception::NotFoundError("MatrixWorkspace::getDetector(): No "
                                           "detectors for this workspace "
                                           "index.",
                                           "");
  }
  // Else need to construct a DetectorGroup and return that
  std::vector<Geometry::IDetector_const_sptr> dets_ptr =
      localInstrument->getDetectors(dets);
  return Geometry::IDetector_const_sptr(
      new Geometry::DetectorGroup(dets_ptr, false));
}

/** Returns the signed 2Theta scattering angle for a detector
*  @param det :: A pointer to the detector object (N.B. might be a
* DetectorGroup)
*  @return The scattering angle (0 < theta < pi)
*  @throws InstrumentDefinitionError if source or sample is missing, or they are
* in the same place
*/
double MatrixWorkspace::detectorSignedTwoTheta(
    Geometry::IDetector_const_sptr det) const {
  Instrument_const_sptr instrument = getInstrument();

  Geometry::IComponent_const_sptr source = instrument->getSource();
  Geometry::IComponent_const_sptr sample = instrument->getSample();
  if (source == NULL || sample == NULL) {
    throw Kernel::Exception::InstrumentDefinitionError(
        "Instrument not sufficiently defined: failed to get source and/or "
        "sample");
  }

  const Kernel::V3D samplePos = sample->getPos();
  const Kernel::V3D beamLine = samplePos - source->getPos();

  if (beamLine.nullVector()) {
    throw Kernel::Exception::InstrumentDefinitionError(
        "Source and sample are at same position!");
  }
  // Get the instrument up axis.
  const V3D &instrumentUpAxis =
      instrument->getReferenceFrame()->vecPointingUp();
  return det->getSignedTwoTheta(samplePos, beamLine, instrumentUpAxis);
}

/** Returns the 2Theta scattering angle for a detector
*  @param det :: A pointer to the detector object (N.B. might be a
* DetectorGroup)
*  @return The scattering angle (0 < theta < pi)
*  @throws InstrumentDefinitionError if source or sample is missing, or they are
* in the same place
*/
double
MatrixWorkspace::detectorTwoTheta(Geometry::IDetector_const_sptr det) const {
  Geometry::IComponent_const_sptr source = getInstrument()->getSource();
  Geometry::IComponent_const_sptr sample = getInstrument()->getSample();
  if (source == NULL || sample == NULL) {
    throw Kernel::Exception::InstrumentDefinitionError(
        "Instrument not sufficiently defined: failed to get source and/or "
        "sample");
  }

  const Kernel::V3D samplePos = sample->getPos();
  const Kernel::V3D beamLine = samplePos - source->getPos();

  if (beamLine.nullVector()) {
    throw Kernel::Exception::InstrumentDefinitionError(
        "Source and sample are at same position!");
  }

  return det->getTwoTheta(samplePos, beamLine);
}

//---------------------------------------------------------------------------------------
/** Add parameters to the instrument parameter map that are defined in
* instrument
*   definition file and for which logfile data are available. Logs must be
* loaded
*   before running this method.
*/
void MatrixWorkspace::populateInstrumentParameters() {
  ExperimentInfo::populateInstrumentParameters();
  // Clear out the nearestNeighbors so that it gets recalculated
  this->m_nearestNeighbours.reset();
}

//----------------------------------------------------------------------------------------------------
/// @return The number of axes which this workspace has
int MatrixWorkspace::axes() const { return static_cast<int>(m_axes.size()); }

//----------------------------------------------------------------------------------------------------
/** Get a pointer to a workspace axis
*  @param axisIndex :: The index of the axis required
*  @throw IndexError If the argument given is outside the range of axes held by
* this workspace
*  @return Pointer to Axis object
*/
Axis *MatrixWorkspace::getAxis(const std::size_t &axisIndex) const {
  if (axisIndex >= m_axes.size()) {
    throw Kernel::Exception::IndexError(
        axisIndex, m_axes.size(),
        "Argument to getAxis is invalid for this workspace");
  }

  return m_axes[axisIndex];
}

/** Replaces one of the workspace's axes with the new one provided.
*  @param axisIndex :: The index of the axis to replace
*  @param newAxis :: A pointer to the new axis. The class will take ownership.
*  @throw IndexError If the axisIndex given is outside the range of axes held by
* this workspace
*  @throw std::runtime_error If the new axis is not of the correct length
* (within one of the old one)
*/
void MatrixWorkspace::replaceAxis(const std::size_t &axisIndex,
                                  Axis *const newAxis) {
  // First check that axisIndex is in range
  if (axisIndex >= m_axes.size()) {
    throw Kernel::Exception::IndexError(
        axisIndex, m_axes.size(),
        "Value of axisIndex is invalid for this workspace");
  }
  // If we're OK, then delete the old axis and set the pointer to the new one
  delete m_axes[axisIndex];
  m_axes[axisIndex] = newAxis;
}

//----------------------------------------------------------------------------------------------------
/// Returns the units of the data in the workspace
std::string MatrixWorkspace::YUnit() const { return m_YUnit; }

/// Sets a new unit for the data (Y axis) in the workspace
void MatrixWorkspace::setYUnit(const std::string &newUnit) {
  m_YUnit = newUnit;
}

/// Returns a caption for the units of the data in the workspace
std::string MatrixWorkspace::YUnitLabel() const {
  std::string retVal;
  if (!m_YUnitLabel.empty())
    retVal = m_YUnitLabel;
  else {
    retVal = m_YUnit;
    // If this workspace a distribution & has at least one axis & this axis has
    // its unit set
    // then append that unit to the string to be returned
    if (!retVal.empty() && this->isDistribution() && this->axes() &&
        this->getAxis(0)->unit()) {
      retVal = retVal + " per " + this->getAxis(0)->unit()->label().ascii();
    }
  }

  return retVal;
}

/// Sets a new caption for the data (Y axis) in the workspace
void MatrixWorkspace::setYUnitLabel(const std::string &newLabel) {
  m_YUnitLabel = newLabel;
}

//----------------------------------------------------------------------------------------------------
/** Are the Y-values in this workspace dimensioned?
* TODO: For example: ????
* @return whether workspace is a distribution or not
*/
const bool &MatrixWorkspace::isDistribution() const { return m_isDistribution; }

/** Set the flag for whether the Y-values are dimensioned
*  @return whether workspace is now a distribution
*/
bool &MatrixWorkspace::isDistribution(bool newValue) {
  m_isDistribution = newValue;
  return m_isDistribution;
}

/**
*  Whether the workspace contains histogram data
*  @return whether the workspace contains histogram data
*/
bool MatrixWorkspace::isHistogramData() const {
  return (readX(0).size() == blocksize() ? false : true);
}

/**
*  Whether the workspace contains common X bins
*  @return whether the workspace contains common X bins
*/
bool MatrixWorkspace::isCommonBins() const {
  if (!m_isCommonBinsFlagSet) {
    m_isCommonBinsFlag = true;

    // there being only one or zero histograms is accepted as not being an error
    if (blocksize() || getNumberHistograms() > 1) {
      // otherwise will compare some of the data, to save time just check two
      // the first and the last
      const size_t lastSpec = getNumberHistograms() - 1;
      // Quickest check is to see if they are actually the same vector
      if (&(readX(0)[0]) != &(readX(lastSpec)[0])) {
        // Now check numerically
        const double first =
            std::accumulate(readX(0).begin(), readX(0).end(), 0.);
        const double last =
            std::accumulate(readX(lastSpec).begin(), readX(lastSpec).end(), 0.);
        if (std::abs(first - last) / std::abs(first + last) > 1.0E-9) {
          m_isCommonBinsFlag = false;
        }

        // handle Nan's and inf's
        if ((boost::math::isinf(first) != boost::math::isinf(last)) ||
            (boost::math::isnan(first) != boost::math::isnan(last))) {
          m_isCommonBinsFlag = false;
        }
      }
    }
    m_isCommonBinsFlagSet = true;
  }

  return m_isCommonBinsFlag;
}
//----------------------------------------------------------------------------------------------------
/**
* Mask a given workspace index, setting the data and error values to zero