MDHistoWorkspace.cpp

#include "MantidGeometry/MDGeometry/IMDDimension.h"
#include "MantidGeometry/MDGeometry/MDGeometryXMLBuilder.h"
#include "MantidKernel/System.h"
#include "MantidKernel/Utils.h"
#include "MantidKernel/VMD.h"
#include "MantidDataObjects/MDHistoWorkspace.h"
#include "MantidDataObjects/MDHistoWorkspaceIterator.h"
#include "MantidGeometry/MDGeometry/MDHistoDimension.h"
#include "MantidGeometry/MDGeometry/MDDimensionExtents.h"
#include <map>
#include "MantidAPI/IMDWorkspace.h"
#include "MantidAPI/IMDIterator.h"
#include <boost/scoped_array.hpp>
#include <boost/make_shared.hpp>
#include <boost/math/special_functions/fpclassify.hpp>

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

namespace Mantid {
namespace DataObjects {
//----------------------------------------------------------------------------------------------
/** Constructor given the 4 dimensions
 * @param dimX :: X dimension binning parameters
 * @param dimY :: Y dimension binning parameters
 * @param dimZ :: Z dimension binning parameters
 * @param dimT :: T (time) dimension binning parameters
 * @param displayNormalization :: optional display normalization to use as the
 * default.
 */
MDHistoWorkspace::MDHistoWorkspace(
    Mantid::Geometry::MDHistoDimension_sptr dimX,
    Mantid::Geometry::MDHistoDimension_sptr dimY,
    Mantid::Geometry::MDHistoDimension_sptr dimZ,
    Mantid::Geometry::MDHistoDimension_sptr dimT,
    Mantid::API::MDNormalization displayNormalization)
    : IMDHistoWorkspace(), numDimensions(0),
      m_nEventsContributed(std::numeric_limits<uint64_t>::quiet_NaN()),
      m_coordSystem(None), m_displayNormalization(displayNormalization) {
  std::vector<Mantid::Geometry::MDHistoDimension_sptr> dimensions;
  if (dimX)
    dimensions.push_back(dimX);
  if (dimY)
    dimensions.push_back(dimY);
  if (dimZ)
    dimensions.push_back(dimZ);
  if (dimT)
    dimensions.push_back(dimT);
  this->init(dimensions);
}

//----------------------------------------------------------------------------------------------
/** Constructor given a vector of dimensions
 * @param dimensions :: vector of MDHistoDimension; no limit to how many.
 * @param displayNormalization :: optional display normalization to use as the
 * default.
 */
MDHistoWorkspace::MDHistoWorkspace(
    std::vector<Mantid::Geometry::MDHistoDimension_sptr> &dimensions,
    Mantid::API::MDNormalization displayNormalization)
    : IMDHistoWorkspace(), numDimensions(0), m_numEvents(NULL),
      m_nEventsContributed(std::numeric_limits<uint64_t>::quiet_NaN()),
      m_coordSystem(None), m_displayNormalization(displayNormalization) {
  this->init(dimensions);
}

//----------------------------------------------------------------------------------------------
/** Constructor given a vector of dimensions
 * @param dimensions :: vector of MDHistoDimension; no limit to how many.
 * @param displayNormalization :: optional display normalization to use as the
 * default.
 */
MDHistoWorkspace::MDHistoWorkspace(
    std::vector<Mantid::Geometry::IMDDimension_sptr> &dimensions,
    Mantid::API::MDNormalization displayNormalization)
    : IMDHistoWorkspace(), numDimensions(0), m_numEvents(NULL),
      m_nEventsContributed(std::numeric_limits<uint64_t>::quiet_NaN()),
      m_coordSystem(None), m_displayNormalization(displayNormalization) {
  this->init(dimensions);
}

//----------------------------------------------------------------------------------------------
/** Copy constructor
 *
 * @param other :: MDHistoWorkspace to copy from.
 */
MDHistoWorkspace::MDHistoWorkspace(const MDHistoWorkspace &other)
    : IMDHistoWorkspace(other) {
  // Dimensions are copied by the copy constructor of MDGeometry
  this->cacheValues();
  // Allocate the linear arrays
  m_signals = new signal_t[m_length];
  m_errorsSquared = new signal_t[m_length];
  m_numEvents = new signal_t[m_length];
  m_masks = new bool[m_length];
  m_nEventsContributed = other.m_nEventsContributed;
  this->setCoordinateSystem(other.getSpecialCoordinateSystem());
  m_displayNormalization = other.displayNormalization();
  // Now copy all the data
  for (size_t i = 0; i < m_length; ++i) {
    m_signals[i] = other.m_signals[i];
    m_errorsSquared[i] = other.m_errorsSquared[i];
    m_numEvents[i] = other.m_numEvents[i];
    m_masks[i] = other.m_masks[i];
  }
}

//----------------------------------------------------------------------------------------------
/** Destructor
 */
MDHistoWorkspace::~MDHistoWorkspace() {
  delete[] m_signals;
  delete[] m_errorsSquared;
  delete[] m_numEvents;
  delete[] indexMultiplier;
  delete[] m_vertexesArray;
  delete[] m_boxLength;
  delete[] m_indexMaker;
  delete[] m_indexMax;
  delete[] m_origin;
  delete[] m_masks;
}

//----------------------------------------------------------------------------------------------
/** Constructor helper method
 * @param dimensions :: vector of MDHistoDimension; no limit to how many.
 */
void MDHistoWorkspace::init(
    std::vector<Mantid::Geometry::MDHistoDimension_sptr> &dimensions) {
  std::vector<IMDDimension_sptr> dim2;
  for (size_t i = 0; i < dimensions.size(); i++)
    dim2.push_back(boost::dynamic_pointer_cast<IMDDimension>(dimensions[i]));
  this->init(dim2);
  m_nEventsContributed = 0;
}

//----------------------------------------------------------------------------------------------
/** Constructor helper method
 * @param dimensions :: vector of IMDDimension; no limit to how many.
 */
void MDHistoWorkspace::init(
    std::vector<Mantid::Geometry::IMDDimension_sptr> &dimensions) {
  MDGeometry::initGeometry(dimensions);
  this->cacheValues();

  // Allocate the linear arrays
  m_signals = new signal_t[m_length];
  m_errorsSquared = new signal_t[m_length];
  m_numEvents = new signal_t[m_length];
  m_masks = new bool[m_length];
  // Initialize them to NAN (quickly)
  signal_t nan = std::numeric_limits<signal_t>::quiet_NaN();
  this->setTo(nan, nan, nan);
  m_nEventsContributed = 0;
}

//----------------------------------------------------------------------------------------------
/** When all dimensions have been initialized, this caches all the necessary
 * values for later use.
 */
void MDHistoWorkspace::cacheValues() {
  // Copy the dimensions array
  numDimensions = m_dimensions.size();

  // For indexing.
  if (numDimensions < 4)
    indexMultiplier = new size_t[4];
  else
    indexMultiplier = new size_t[numDimensions];

  // For quick indexing, accumulate these values
  // First multiplier
  indexMultiplier[0] = m_dimensions[0]->getNBins();
  for (size_t d = 1; d < numDimensions; d++)
    indexMultiplier[d] = indexMultiplier[d - 1] * m_dimensions[d]->getNBins();

  // This is how many dense data points
  m_length = indexMultiplier[numDimensions - 1];

  // Now fix things for < 4 dimensions. Indices > the number of dimensions will
  // be ignored (*0)
  for (size_t d = numDimensions - 1; d < 4; d++)
    indexMultiplier[d] = 0;

  // Compute the volume of each cell.
  coord_t volume = 1.0;
  for (size_t i = 0; i < numDimensions; ++i)
    volume *= m_dimensions[i]->getBinWidth();
  m_inverseVolume = 1.0f / volume;

  // Continue with the vertexes array
  this->initVertexesArray();
  m_nEventsContributed = 0;
}
//----------------------------------------------------------------------------------------------
/** After initialization, call this to initialize the vertexes array
 * to the vertexes of the 0th box.
 * Will be used by getVertexesArray()
 * */