WorkspaceCreationHelper.cpp

/*********************************************************************************
 *  PLEASE READ THIS!!!!!!!
 *
 *  This collection of functions MAY NOT be used in any test from a package
 *below
 *  DataObjects (e.g. Kernel, Geometry, API).
 *  Conversely, this file MAY NOT be modified to use anything from a package
 *higher
 *  than DataObjects (e.g. any algorithm), even if going via the factory.
 *********************************************************************************/
#include "MantidTestHelpers/WorkspaceCreationHelper.h"
#include "MantidTestHelpers/ComponentCreationHelper.h"
#include "MantidTestHelpers/InstrumentCreationHelper.h"

#include "MantidAPI/Algorithm.h"
#include "MantidAPI/NumericAxis.h"
#include "MantidAPI/Sample.h"
#include "MantidAPI/SpectraAxis.h"
#include "MantidAPI/SpectrumInfo.h"
#include "MantidAPI/WorkspaceGroup.h"
#include "MantidDataObjects/PeaksWorkspace.h"
#include "MantidDataObjects/ScanningWorkspaceBuilder.h"
#include "MantidDataObjects/WorkspaceCreation.h"
#include "MantidGeometry/Crystal/OrientedLattice.h"
#include "MantidGeometry/Instrument/Component.h"
#include "MantidGeometry/Instrument/Detector.h"
#include "MantidGeometry/Instrument/DetectorInfo.h"
#include "MantidGeometry/Instrument/Goniometer.h"
#include "MantidGeometry/Instrument/ParameterMap.h"
#include "MantidGeometry/Objects/ShapeFactory.h"
#include "MantidHistogramData/LinearGenerator.h"
#include "MantidIndexing/IndexInfo.h"
#include "MantidKernel/MersenneTwister.h"
#include "MantidKernel/OptionalBool.h"
#include "MantidKernel/TimeSeriesProperty.h"
#include "MantidKernel/UnitFactory.h"
#include "MantidKernel/VectorHelper.h"
#include "MantidKernel/make_unique.h"
#include "MantidKernel/V3D.h"

#include <cmath>
#include <sstream>
#include <utility>

namespace WorkspaceCreationHelper {
using namespace Mantid;
using namespace Mantid::DataObjects;
using namespace Mantid::Kernel;
using namespace Mantid::API;
using namespace Mantid::Geometry;
using namespace Mantid::HistogramData;
using Mantid::MantidVec;
using Mantid::MantidVecPtr;
using Mantid::Types::Core::DateAndTime;
using Mantid::Types::Event::TofEvent;

MockAlgorithm::MockAlgorithm(size_t nSteps)
    : m_Progress(
          Mantid::Kernel::make_unique<API::Progress>(this, 0.0, 1.0, nSteps)) {}

EPPTableRow::EPPTableRow(const double peakCentre_, const double sigma_,
                         const double height_, const FitStatus fitStatus_)
    : peakCentre(peakCentre_), sigma(sigma_), height(height_),
      fitStatus(fitStatus_) {}

EPPTableRow::EPPTableRow(const int index, const double peakCentre_,
                         const double sigma_, const double height_,
                         const FitStatus fitStatus_)
    : workspaceIndex(index), peakCentre(peakCentre_), sigma(sigma_),
      height(height_), fitStatus(fitStatus_) {}

/**
 * @param name :: The name of the workspace
 * @param ws :: The workspace object
 */
void storeWS(const std::string &name, Mantid::API::Workspace_sptr ws) {
  Mantid::API::AnalysisDataService::Instance().add(name, ws);
}

/**
 * Deletes a workspace
 * @param name :: The name of the workspace
 */
void removeWS(const std::string &name) {
  Mantid::API::AnalysisDataService::Instance().remove(name);
}

/**
  * Creates bin or point based histograms based on the data passed
  * in for Y and E values and the bool specified.
  *
  * @param isHistogram :: Specifies whether the returned histogram
  * should use points or bin edges for the x axis. True gives bin edges.
  * @param yAxis :: Takes an rvalue (move) of the y axis for the new histogram
  * @param eAxis :: Takes an rvalue (move) of the e axis for the new histogram
  *
  * @return :: Returns a histogram with the user specified X axis type
  * and the data the user passed in.
  */
template <typename YType, typename EType>
Histogram createHisto(bool isHistogram, YType &&yAxis, EType &&eAxis) {
  // We don't need to check if y.size() == e.size() as the histogram
  // type does this at construction
  const size_t yValsSize = yAxis.size();
  if (isHistogram) {
    BinEdges xAxis(yValsSize + 1, LinearGenerator(1, 1));
    Histogram histo{std::move(xAxis), std::forward<YType>(yAxis),
                    std::forward<EType>(eAxis)};
    return histo;
  } else {
    Points xAxis(yValsSize, LinearGenerator(1, 1));
    Histogram pointsHisto{std::move(xAxis), std::forward<YType>(yAxis),
                          std::forward<EType>(eAxis)};
    return pointsHisto;
  }
}

Workspace2D_sptr create1DWorkspaceRand(int size, bool isHisto) {

  MersenneTwister randomGen(DateAndTime::getCurrentTime().nanoseconds(), 0,
                            std::numeric_limits<int>::max());

  auto randFunc = [&randomGen] { return randomGen.nextValue(); };
  Counts counts(size, randFunc);
  CountStandardDeviations errorVals(size, randFunc);

  auto retVal = boost::make_shared<Workspace2D>();
  retVal->initialize(1, createHisto(isHisto, counts, errorVals));
  return retVal;
}

Workspace2D_sptr create1DWorkspaceConstant(int size, double value, double error,
                                           bool isHisto) {
  Counts yVals(size, value);
  CountStandardDeviations errVals(size, error);

  auto retVal = boost::make_shared<Workspace2D>();
  retVal->initialize(1, createHisto(isHisto, yVals, errVals));
  return retVal;
}

Workspace2D_sptr create1DWorkspaceConstantWithXerror(int size, double value,
                                                     double error,
                                                     double xError,
                                                     bool isHisto) {
  auto ws = create1DWorkspaceConstant(size, value, error, isHisto);
  auto dx1 = Kernel::make_cow<HistogramData::HistogramDx>(size, xError);
  ws->setSharedDx(0, dx1);
  return ws;
}

Workspace2D_sptr create1DWorkspaceFib(int size, bool isHisto) {
  BinEdges xVals(size + 1, LinearGenerator(1, 1));
  Counts yVals(size, FibSeries<double>());
  CountStandardDeviations errVals(size);

  auto retVal = boost::make_shared<Workspace2D>();
  retVal->initialize(1, createHisto(isHisto, yVals, errVals));
  return retVal;
}

Workspace2D_sptr create2DWorkspace(int nhist, int numBoundaries) {
  return create2DWorkspaceBinned(nhist, numBoundaries);
}

/** Create a Workspace2D where the Y value at each bin is
 * == to the workspace index
 * @param nhist :: # histograms
 * @param numBoundaries :: # of bins
 * @return Workspace2D
 */
Workspace2D_sptr create2DWorkspaceWhereYIsWorkspaceIndex(int nhist,
                                                         int numBoundaries) {
  Workspace2D_sptr out = create2DWorkspaceBinned(nhist, numBoundaries);
  for (int workspaceIndex = 0; workspaceIndex < nhist; workspaceIndex++) {
    std::vector<double> yValues(numBoundaries,
                                static_cast<double>(workspaceIndex));
    out->mutableY(workspaceIndex) = std::move(yValues);
  }

  return out;
}

Workspace2D_sptr create2DWorkspaceThetaVsTOF(int nHist, int nBins) {

  Workspace2D_sptr outputWS = create2DWorkspaceBinned(nHist, nBins);
  auto const newAxis = new NumericAxis(nHist);
  outputWS->replaceAxis(1, newAxis);
  newAxis->unit() = boost::make_shared<Units::Degrees>();
  for (int i = 0; i < nHist; ++i) {
    newAxis->setValue(i, i + 1);
  }

  return outputWS;
}

Workspace2D_sptr
create2DWorkspaceWithValues(int64_t nHist, int64_t nBins, bool isHist,
                            const std::set<int64_t> &maskedWorkspaceIndices,
                            double xVal, double yVal, double eVal) {
  auto x1 = Kernel::make_cow<HistogramData::HistogramX>(
      isHist ? nBins + 1 : nBins, LinearGenerator(xVal, 1.0));
  Counts y1(nBins, yVal);
  CountStandardDeviations e1(nBins, eVal);
  auto retVal = boost::make_shared<Workspace2D>();
  retVal->initialize(nHist, isHist ? nBins + 1 : nBins, nBins);
  for (int i = 0; i < nHist; i++) {
    retVal->setX(i, x1);
    retVal->setCounts(i, y1);
    retVal->setCountStandardDeviations(i, e1);
    retVal->getSpectrum(i).setDetectorID(i);
    retVal->getSpectrum(i).setSpectrumNo(i);
  }
  retVal = maskSpectra(retVal, maskedWorkspaceIndices);
  return retVal;
}

Workspace2D_sptr create2DWorkspaceWithValuesAndXerror(
    int64_t nHist, int64_t nBins, bool isHist, double xVal, double yVal,
    double eVal, double dxVal,
    const std::set<int64_t> &maskedWorkspaceIndices) {
  auto ws = create2DWorkspaceWithValues(
      nHist, nBins, isHist, maskedWorkspaceIndices, xVal, yVal, eVal);
  PointStandardDeviations dx1(nBins, dxVal);
  for (int i = 0; i < nHist; i++) {
    ws->setPointStandardDeviations(i, dx1);
  }
  return ws;
}

Workspace2D_sptr
create2DWorkspace123(int64_t nHist, int64_t nBins, bool isHist,
                     const std::set<int64_t> &maskedWorkspaceIndices) {
  return create2DWorkspaceWithValues(nHist, nBins, isHist,
                                     maskedWorkspaceIndices, 1.0, 2.0, 3.0);
}

Workspace2D_sptr
create2DWorkspace154(int64_t nHist, int64_t nBins, bool isHist,
                     const std::set<int64_t> &maskedWorkspaceIndices) {
  return create2DWorkspaceWithValues(nHist, nBins, isHist,
                                     maskedWorkspaceIndices, 1.0, 5.0, 4.0);
}

Workspace2D_sptr maskSpectra(Workspace2D_sptr workspace,
                             const std::set<int64_t> &maskedWorkspaceIndices) {
  const int nhist = static_cast<int>(workspace->getNumberHistograms());
  if (workspace->getInstrument()->nelements() == 0) {
    // We need detectors to be able to mask them.
    auto instrument = boost::make_shared<Instrument>();
    workspace->setInstrument(instrument);

    std::string xmlShape = "<sphere id=\"shape\"> ";
    xmlShape += "<centre x=\"0.0\"  y=\"0.0\" z=\"0.0\" /> ";
    xmlShape += "<radius val=\"0.05\" /> ";
    xmlShape += "</sphere>";
    xmlShape += "<algebra val=\"shape\" /> ";

    ShapeFactory sFactory;
    boost::shared_ptr<Object> shape = sFactory.createShape(xmlShape);
    for (int i = 0; i < nhist; ++i) {
      Detector *det = new Detector("det", detid_t(i + 1), shape, nullptr);
      det->setPos(i, i + 1, 1);
      instrument->add(det);
      instrument->markAsDetector(det);
    }
    workspace->setInstrument(instrument);
    // Set IndexInfo without explicit spectrum definitions to trigger building
    // default mapping of spectra to detectors in new instrument.
    workspace->setIndexInfo(Indexing::IndexInfo(nhist));
  }

  auto &spectrumInfo = workspace->mutableSpectrumInfo();
  for (const auto index : maskedWorkspaceIndices)
    spectrumInfo.setMasked(index, true);
  return workspace;
}

/**
 * Create a group with nEntries. It is added to the ADS with the given stem
 */
WorkspaceGroup_sptr createWorkspaceGroup(int nEntries, int nHist, int nBins,
                                         const std::string &stem) {
  auto group = boost::make_shared<WorkspaceGroup>();
  AnalysisDataService::Instance().add(stem, group);
  for (int i = 0; i < nEntries; ++i) {
    Workspace2D_sptr ws = create2DWorkspace(nHist, nBins);
    std::ostringstream os;
    os << stem << "_" << i;
    AnalysisDataService::Instance().add(os.str(), ws);
    group->add(os.str());
  }
  return group;
}

/** Create a 2D workspace with this many histograms and bins.
 * Filled with Y = 2.0 and E = M_SQRT2w
 */
Workspace2D_sptr create2DWorkspaceBinned(int nhist, int numVals, double x0,
                                         double deltax) {
  BinEdges x(numVals + 1, LinearGenerator(x0, deltax));
  Counts y(numVals, 2);
  CountStandardDeviations e(numVals, M_SQRT2);
  auto retVal = boost::make_shared<Workspace2D>();
  retVal->initialize(nhist, numVals + 1, numVals);
  for (int i = 0; i < nhist; i++) {
    retVal->setBinEdges(i, x);
    retVal->setCounts(i, y);
    retVal->setCountStandardDeviations(i, e);
  }
  return retVal;
}

/** Create a 2D workspace with this many histograms and bins. The bins are
 * assumed to be non-uniform and given by the input array
 * Filled with Y = 2.0 and E = M_SQRT2w
 */
Workspace2D_sptr create2DWorkspaceBinned(int nhist, const int numBoundaries,
                                         const double xBoundaries[]) {
  BinEdges x(xBoundaries, xBoundaries + numBoundaries);
  const int numBins = numBoundaries - 1;
  Counts y(numBins, 2);
  CountStandardDeviations e(numBins, M_SQRT2);
  auto retVal = boost::make_shared<Workspace2D>();
  retVal->initialize(nhist, numBins + 1, numBins);
  for (int i = 0; i < nhist; i++) {
    retVal->setBinEdges(i, x);
    retVal->setCounts(i, y);
    retVal->setCountStandardDeviations(i, e);
  }
  return retVal;
}

/**
 * Add random noise to the signal
 * @param ws :: The workspace to add the noise to
 * @param noise :: The mean noise level
 * @param lower :: The lower bound of the flucation (default=-0.5)
 * @param upper:: The upper bound of the flucation (default=-0.5)
 */
void addNoise(Mantid::API::MatrixWorkspace_sptr ws, double noise,
              const double lower, const double upper) {
  const size_t seed(12345);
  MersenneTwister randGen(seed, lower, upper);
  for (size_t iSpec = 0; iSpec < ws->getNumberHistograms(); iSpec++) {
    auto &mutableY = ws->mutableY(iSpec);
    auto &mutableE = ws->mutableE(iSpec);
    for (size_t i = 0; i < mutableY.size(); i++) {
      mutableY[i] += noise * randGen.nextValue();
      mutableE[i] += noise;
    }
  }
}

//================================================================================================================
/**
 * Create a test workspace with a fully defined instrument
 * Each spectra will have a cylindrical detector defined 2*cylinder_radius away
 * from the centre of the
 * previous.
 * Data filled with: Y: 2.0, E: M_SQRT2, X: nbins of width 1 starting at 0
 */
Workspace2D_sptr
create2DWorkspaceWithFullInstrument(int nhist, int nbins, bool includeMonitors,
                                    bool startYNegative, bool isHistogram,
                                    const std::string &instrumentName) {
  if (includeMonitors && nhist < 2) {
    throw std::invalid_argument("Attempting to 2 include monitors for a "
                                "workspace with fewer than 2 histograms");
  }

  Workspace2D_sptr space;
  if (isHistogram)
    space = create2DWorkspaceBinned(
        nhist, nbins); // A 1:1 spectra is created by default
  else
    space = create2DWorkspace123(nhist, nbins, false);
  space->setTitle(
      "Test histogram"); // actually adds a property call run_title to the logs
  space->getAxis(0)->setUnit("TOF");
  space->setYUnit("Counts");

  InstrumentCreationHelper::addFullInstrumentToWorkspace(
      *space, includeMonitors, startYNegative, instrumentName);

  return space;
}

//================================================================================================================
/*
 * startTime is in seconds
 */
MatrixWorkspace_sptr create2DDetectorScanWorkspaceWithFullInstrument(
    int nhist, int nbins, size_t nTimeIndexes, size_t startTime,
    size_t firstInterval, bool includeMonitors, bool startYNegative,
    bool isHistogram, const std::string &instrumentName) {

  auto baseWS = create2DWorkspaceWithFullInstrument(
      nhist, nbins, includeMonitors, startYNegative, isHistogram,
      instrumentName);

  auto builder =
      ScanningWorkspaceBuilder(baseWS->getInstrument(), nTimeIndexes, nbins);

  std::vector<double> timeRanges;
  for (size_t i = 0; i < nTimeIndexes; ++i) {
    timeRanges.push_back(double(i + firstInterval));
  }

  builder.setTimeRanges(DateAndTime(int(startTime), 0), timeRanges);

  return builder.buildWorkspace();
}

//================================================================================================================
/** Create an Workspace2D with an instrument that contains
 *RectangularDetector's.
 * Bins will be 0.0, 1.0, to numBins, filled with signal=2.0, M_SQRT2
 *
 * @param numBanks :: number of rectangular banks
 * @param numPixels :: each bank will be numPixels*numPixels
 * @param numBins :: each spectrum will have this # of bins
 * @return The EventWorkspace
 */
Mantid::DataObjects::Workspace2D_sptr
create2DWorkspaceWithRectangularInstrument(int numBanks, int numPixels,
                                           int numBins) {
  Instrument_sptr inst =
      ComponentCreationHelper::createTestInstrumentRectangular(numBanks,
                                                               numPixels);
  Workspace2D_sptr ws =
      create2DWorkspaceBinned(numBanks * numPixels * numPixels, numBins);
  ws->setInstrument(inst);
  ws->getAxis(0)->setUnit("dSpacing");
  for (size_t wi = 0; wi < ws->getNumberHistograms(); wi++) {
    ws->getSpectrum(wi).setDetectorID(detid_t(numPixels * numPixels + wi));
    ws->getSpectrum(wi).setSpectrumNo(specnum_t(wi));
  }

  return ws;
}

//================================================================================================================
/** Create an Eventworkspace with an instrument that contains
 *RectangularDetector's.
 * X axis = 100 histogrammed bins from 0.0 in steps of 1.0.
 * 200 events per pixel.
 *
 * @param numBanks :: number of rectangular banks
 * @param numPixels :: each bank will be numPixels*numPixels
 * @param clearEvents :: if true, erase the events from list
 * @return The EventWorkspace
 */
Mantid::DataObjects::EventWorkspace_sptr
createEventWorkspaceWithFullInstrument(int numBanks, int numPixels,
                                       bool clearEvents) {
  Instrument_sptr inst =
      ComponentCreationHelper::createTestInstrumentRectangular(numBanks,
                                                               numPixels);
  EventWorkspace_sptr ws =
      createEventWorkspace2(numBanks * numPixels * numPixels, 100);
  ws->setInstrument(inst);

  // Set the X axes
  const auto &xVals = ws->x(0);
  const size_t xSize = xVals.size();
  auto ax0 = new NumericAxis(xSize);
  ax0->setUnit("dSpacing");
  for (size_t i = 0; i < xSize; i++) {
    ax0->setValue(i, xVals[i]);
  }
  ws->replaceAxis(0, ax0);

  // re-assign detector IDs to the rectangular detector
  int detID = numPixels * numPixels;
  for (int wi = 0; wi < static_cast<int>(ws->getNumberHistograms()); wi++) {
    ws->getSpectrum(wi).clearDetectorIDs();
    if (clearEvents)
      ws->getSpectrum(wi).clear(true);
    ws->getSpectrum(wi).setDetectorID(detID);
    detID++;
  }
  return ws;
}

Mantid::DataObjects::EventWorkspace_sptr
createEventWorkspaceWithNonUniformInstrument(int numBanks, bool clearEvents) {
  // Number of detectors in a bank as created by createTestInstrumentCylindrical
  const int DETECTORS_PER_BANK(9);

  V3D srcPos(0., 0., -10.), samplePos;
  Instrument_sptr inst =
      ComponentCreationHelper::createTestInstrumentCylindrical(
          numBanks, srcPos, samplePos, 0.0025, 0.005);
  EventWorkspace_sptr ws =
      createEventWorkspace2(numBanks * DETECTORS_PER_BANK, 100);
  ws->setInstrument(inst);

  std::vector<detid_t> detectorIds = inst->getDetectorIDs();

  // Should be equal if DETECTORS_PER_BANK is correct
  assert(detectorIds.size() == ws->getNumberHistograms());

  // Re-assign detector IDs
  for (size_t wi = 0; wi < ws->getNumberHistograms(); wi++) {
    ws->getSpectrum(wi).clearDetectorIDs();
    if (clearEvents)
      ws->getSpectrum(wi).clear(true);
    ws->getSpectrum(wi).setDetectorID(detectorIds[wi]);
  }

  return ws;
}

/** Adds a component to an instrument
 *
 * @param instrument :: instrument to which the component will be added
 * @param position :: position of the component
 * @param name :: name of the component
 * @return a component pointer
 */
ObjComponent *addComponent(Instrument_sptr &instrument, const V3D &position,
                           std::string name = "component") {
  ObjComponent *component = new ObjComponent(name);
  component->setPos(position);
  instrument->add(component);
  return component;
}

/** Adds a source to an instrument
 *
 * @param instrument :: instrument to which the source will be added
 * @param position :: position of the source
 * @param name :: name of the source
 */
void addSource(Instrument_sptr &instrument, const V3D &position,
               std::string name = "source") {
  auto source = addComponent(instrument, position, name);
  instrument->markAsSource(source);
}

/** Adds a sample to an instrument
 *
 * @param instrument :: instrument to which the sample will be added
 * @param position :: position of the sample
 * @param name :: name of the sample
 */
void addSample(Instrument_sptr &instrument, const V3D &position,
               std::string name = "some-surface-holder") {
  auto sample = addComponent(instrument, position, name);
  instrument->markAsSamplePos(sample);
}

/** Adds a monitor to an instrument
 *
 * @param instrument :: instrument to which the monitor will be added
 * @param position :: position of the monitor
 * @param ID :: identification number of the monitor
 * @param name :: name of the monitor
 */
void addMonitor(Instrument_sptr &instrument, const V3D &position, const int ID,
                std::string name = "Monitor") {
  Detector *monitor = new Detector(name, ID, nullptr);
  monitor->setPos(position);
  instrument->add(monitor);
  instrument->markAsMonitor(monitor);
}

/** Adds a detector to an instrument
 *
 * @param instrument :: instrument to which the detector will be added
 * @param position :: position of the detector
 * @param ID :: identification number of the detector
 * @param name :: name of the detector
 */
void addDetector(Instrument_sptr &instrument, const V3D &position, const int ID,
                 std::string name = "point-detector") {
  // Where 0.01 is half detector width etc.
  Detector *detector = new Detector(
      name, ID, ComponentCreationHelper::createCuboid(0.01, 0.02, 0.03),
      nullptr);
  detector->setPos(position);
  instrument->add(detector);
  instrument->markAsDetector(detector);
}

/** Creates a binned 2DWorkspace with title and TOF x-axis and counts y-axis
 *
 * @param startX :: start TOF x-value
 * @param nSpectra :: number of spectra
 * @param nBins :: number of bins
 * @param deltaX :: TOF delta x-value
 * @return a Workspace2D
 */
DataObjects::Workspace2D_sptr reflectometryWorkspace(const double startX,
                                                     const int nSpectra,
                                                     const int nBins,
                                                     const double deltaX) {

  auto workspace = create2DWorkspaceBinned(nSpectra, nBins, startX, deltaX);

  workspace->setTitle(
      "Test histogram"); // actually adds a property call run_title to the logs
  workspace->getAxis(0)->setUnit("TOF");
  workspace->setYUnit("Counts");
  return workspace;
}

/**
 * Create a very small 2D workspace for a virtual reflectometry instrument.
 * @return workspace with instrument attached.
 * @param startX : X Tof start value for the workspace.
 * @param slit1Pos :: slit 1 position
 * @param slit2Pos :: slit 2 position
 * @param vg1 :: vertical slit 1
 * @param vg2 :: vertical slit 2
 * @param sourcePos :: source position
 * @param monitorPos :: monitor position
 * @param samplePos :: sample position
 * @param detectorPos :: detector position
 * @param up :: pointing up axis
 * @param along :: pointing along beam axis
 * @param nSpectra :: number of spectra
 * @param nBins :: number of bins
 * @param deltaX :: TOF delta x-value
 */
MatrixWorkspace_sptr create2DWorkspaceWithReflectometryInstrument(
    double startX, V3D slit1Pos, V3D slit2Pos, double vg1, double vg2,
    V3D sourcePos, V3D monitorPos, V3D samplePos, V3D detectorPos,
    PointingAlong up, PointingAlong along, const int nSpectra, const int nBins,
    const double deltaX) {
  Instrument_sptr instrument = boost::make_shared<Instrument>();
  instrument->setReferenceFrame(
      boost::make_shared<ReferenceFrame>(up, along, Handedness::Left, "0,0,0"));

  addSource(instrument, sourcePos);
  addMonitor(instrument, monitorPos, 1);
  addSample(instrument, samplePos);
  addDetector(instrument, detectorPos, 2);
  auto slit1 = addComponent(instrument, slit1Pos, "slit1");
  auto slit2 = addComponent(instrument, slit2Pos, "slit2");

  auto workspace = reflectometryWorkspace(startX, nSpectra, nBins, deltaX);
  workspace->setInstrument(instrument);

  ParameterMap &pmap = workspace->instrumentParameters();
  pmap.addDouble(slit1, "vertical gap", vg1);
  pmap.addDouble(slit2, "vertical gap", vg2);

  workspace->getSpectrum(0).setDetectorID(2);
  workspace->getSpectrum(1).setDetectorID(1);

  return workspace;
}

/**
* Create a very small 2D workspace for a virtual reflectometry instrument with
* multiple detectors
* @return workspace with instrument attached.
* @param startX : X Tof start value for the workspace.
* @param detSize : optional detector height (default is 0 which puts all
* detectors at the same position)
* @param nSpectra :: number of spectra
* @param nBins :: number of bins
* @param deltaX :: TOF delta x-value
*/
MatrixWorkspace_sptr create2DWorkspaceWithReflectometryInstrumentMultiDetector(
    double startX, const double detSize, const int nSpectra, const int nBins,
    const double deltaX) {
  Instrument_sptr instrument = boost::make_shared<Instrument>();
  instrument->setReferenceFrame(
      boost::make_shared<ReferenceFrame>(Y /*up*/, X /*along*/, Left, "0,0,0"));

  addSource(instrument, V3D(0, 0, 0));
  addSample(instrument, V3D(15, 0, 0));
  addMonitor(instrument, V3D(14, 0, 0), 1);

  // Place the central detector at 45 degrees (i.e. the distance
  // from the sample in Y is the same as the distance in X).
  const double posX = 20;
  const double posY = posX - 15;
  addDetector(instrument, V3D(posX, posY - detSize, 0),
              2);                                 // offset below centre
  addDetector(instrument, V3D(posX, posY, 0), 3); // at centre
  addDetector(instrument, V3D(posX, posY + detSize, 0),
              4); // offset above centre

  auto workspace = reflectometryWorkspace(startX, nSpectra, nBins, deltaX);
  workspace->setInstrument(instrument);
  workspace->getSpectrum(0).setDetectorID(1);
  workspace->getSpectrum(1).setDetectorID(2);
  workspace->getSpectrum(2).setDetectorID(3);
  workspace->getSpectrum(3).setDetectorID(4);
  return workspace;
}

void createInstrumentForWorkspaceWithDistances(
    MatrixWorkspace_sptr workspace, const V3D &samplePosition,
    const V3D &sourcePosition, const std::vector<V3D> &detectorPositions) {
  Instrument_sptr instrument = boost::make_shared<Instrument>();
  instrument->setReferenceFrame(
      boost::make_shared<ReferenceFrame>(Y, X, Left, "0,0,0"));

  addSource(instrument, sourcePosition);
  addSample(instrument, samplePosition, "sample");

  for (int i = 0; i < static_cast<int>(detectorPositions.size()); ++i) {
    std::stringstream buffer;
    buffer << "detector_" << i;
    addDetector(instrument, detectorPositions[i], i, buffer.str());

    // Link it to the workspace
    workspace->getSpectrum(i).addDetectorID(i);
  }
  workspace->setInstrument(instrument);
}

//================================================================================================================
WorkspaceSingleValue_sptr createWorkspaceSingleValue(double value) {
  return boost::make_shared<WorkspaceSingleValue>(value, sqrt(value));
}

WorkspaceSingleValue_sptr createWorkspaceSingleValueWithError(double value,
                                                              double error) {
  return boost::make_shared<WorkspaceSingleValue>(value, error);
}

/** Perform some finalization on event workspace stuff */
void eventWorkspace_Finalize(EventWorkspace_sptr ew) {
  // get a proton charge
  ew->mutableRun().integrateProtonCharge();
}

/** Create event workspace with:
 * 500 pixels
 * 1000 histogrammed bins.
 */
EventWorkspace_sptr createEventWorkspace() {
  return createEventWorkspace(500, 1001, 100, 1000);
}

/** Create event workspace with:
 * numPixels pixels
 * numBins histogrammed bins from 0.0 in steps of 1.0
 * 200 events; two in each bin, at time 0.5, 1.5, etc.
 * PulseTime = 0 second x2, 1 second x2, 2 seconds x2, etc. after 2010-01-01
 */
EventWorkspace_sptr createEventWorkspace2(int numPixels, int numBins) {
  return createEventWorkspace(numPixels, numBins, 100, 0.0, 1.0, 2);
}

/** Create event workspace
 */
EventWorkspace_sptr createEventWorkspace(int numPixels, int numBins,
                                         int numEvents, double x0,
                                         double binDelta, int eventPattern,
                                         int start_at_pixelID) {
  return createEventWorkspaceWithStartTime(
      numPixels, numBins, numEvents, x0, binDelta, eventPattern,
      start_at_pixelID, DateAndTime("2010-01-01T00:00:00"));
}

/**
 * Create event workspace with defined start date time
 */
EventWorkspace_sptr
createEventWorkspaceWithStartTime(int numPixels, int numBins, int numEvents,
                                  double x0, double binDelta, int eventPattern,
                                  int start_at_pixelID, DateAndTime run_start) {

  // add one to the number of bins as this is histogram
  numBins++;

  auto retVal = boost::make_shared<EventWorkspace>();
  retVal->initialize(numPixels, 1, 1);

  // Make fake events
  if (eventPattern) // 0 == no events
  {
    size_t workspaceIndex = 0;
    for (int pix = start_at_pixelID + 0; pix < start_at_pixelID + numPixels;
         pix++) {
      EventList &el = retVal->getSpectrum(workspaceIndex);
      el.setSpectrumNo(pix);
      el.setDetectorID(pix);

      for (int i = 0; i < numEvents; i++) {
        if (eventPattern == 1) // 0, 1 diagonal pattern
          el += TofEvent((pix + i + 0.5) * binDelta, run_start + double(i));
        else if (eventPattern == 2) // solid 2
        {
          el += TofEvent((i + 0.5) * binDelta, run_start + double(i));
          el += TofEvent((i + 0.5) * binDelta, run_start + double(i));
        } else if (eventPattern == 3) // solid 1
        {
          el += TofEvent((i + 0.5) * binDelta, run_start + double(i));
        } else if (eventPattern == 4) // Number of events per bin = pixelId (aka
                                      // workspace index in most cases)
        {
          for (int q = 0; q < pix; q++)
            el += TofEvent((i + 0.5) * binDelta, run_start + double(i));
        }
      }
      workspaceIndex++;
    }
  }

  retVal->setAllX(BinEdges(numBins, LinearGenerator(x0, binDelta)));

  return retVal;
}

// =====================================================================================
/** Create event workspace, with several detector IDs in one event list.
 */
EventWorkspace_sptr
createGroupedEventWorkspace(std::vector<std::vector<int>> groups, int numBins,
                            double binDelta, double xOffset) {

  auto retVal = boost::make_shared<EventWorkspace>();
  retVal->initialize(groups.size(), 2, 1);

  for (size_t g = 0; g < groups.size(); g++) {
    retVal->getSpectrum(g).clearDetectorIDs();
    std::vector<int> dets = groups[g];
    for (auto det : dets) {
      for (int i = 0; i < numBins; i++)
        retVal->getSpectrum(g) += TofEvent((i + 0.5) * binDelta, 1);
      retVal->getSpectrum(g).addDetectorID(det);
    }
  }

  if (xOffset == 0.) {
    retVal->setAllX(BinEdges(numBins, LinearGenerator(0.0, binDelta)));
  } else {
    for (size_t g = 0; g < groups.size(); g++) {
      // Create the x-axis for histogramming.
      const double x0 = xOffset * static_cast<double>(g);
      retVal->setX(
          g, make_cow<HistogramX>(numBins, LinearGenerator(x0, binDelta)));
    }
  }

  return retVal;
}

// =====================================================================================
/** Create an event workspace with randomized TOF and pulsetimes
 *
 * @param numbins :: # of bins to set. This is also = # of events per EventList
 * @param numpixels :: number of pixels
 * @param bin_delta :: a constant offset to shift the bin bounds by
 * @return EventWorkspace
 */
EventWorkspace_sptr createRandomEventWorkspace(size_t numbins, size_t numpixels,
                                               double bin_delta) {
  auto retVal = boost::make_shared<EventWorkspace>();
  retVal->initialize(numpixels, numbins, numbins - 1);

  // and X-axis for references:
  auto pAxis0 = new NumericAxis(numbins);
  // Create the original X axis to histogram on.
  // Create the x-axis for histogramming.
  HistogramData::BinEdges axis(numbins, LinearGenerator(0.0, bin_delta));
  for (int i = 0; i < static_cast<int>(numbins); ++i) {
    pAxis0->setValue(i, axis[i]);
  }
  pAxis0->setUnit("TOF");

  MersenneTwister randomGen(DateAndTime::getCurrentTime().nanoseconds(), 0,
                            std::numeric_limits<int>::max());
  // Make up some data for each pixels
  for (size_t i = 0; i < numpixels; i++) {
    // Create one event for each bin
    EventList &events = retVal->getSpectrum(static_cast<detid_t>(i));
    for (std::size_t ie = 0; ie < numbins; ie++) {
      // Create a list of events, randomize
      events += TofEvent(static_cast<double>(randomGen.nextValue()),
                         static_cast<int64_t>(randomGen.nextValue()));
    }
    events.addDetectorID(detid_t(i));
  }
  retVal->setAllX(axis);
  retVal->replaceAxis(0, pAxis0);

  return retVal;
}

// =====================================================================================
/** Create Workspace2d, with numHist spectra, each with 9 detectors,
 * with IDs 1-9, 10-18, 19-27
 */
MatrixWorkspace_sptr createGroupedWorkspace2D(size_t numHist, int numBins,
                                              double binDelta) {
  Workspace2D_sptr retVal = create2DWorkspaceBinned(static_cast<int>(numHist),
                                                    numBins, 0.0, binDelta);
  retVal->setInstrument(
      ComponentCreationHelper::createTestInstrumentCylindrical(
          static_cast<int>(numHist)));

  for (int g = 0; g < static_cast<int>(numHist); g++) {
    auto &spec = retVal->getSpectrum(g);
    for (int i = 1; i <= 9; i++)
      spec.addDetectorID(g * 9 + i);
    spec.setSpectrumNo(g + 1); // Match detector ID and spec NO
  }
  return boost::dynamic_pointer_cast<MatrixWorkspace>(retVal);
}

// =====================================================================================
// RootOfNumHist == square root of hystohram number;
MatrixWorkspace_sptr
createGroupedWorkspace2DWithRingsAndBoxes(size_t RootOfNumHist, int numBins,
                                          double binDelta) {
  size_t numHist = RootOfNumHist * RootOfNumHist;
  Workspace2D_sptr retVal = create2DWorkspaceBinned(static_cast<int>(numHist),
                                                    numBins, 0.0, binDelta);
  retVal->setInstrument(
      ComponentCreationHelper::createTestInstrumentCylindrical(
          static_cast<int>(numHist)));
  for (int g = 0; g < static_cast<int>(numHist); g++) {
    auto &spec = retVal->getSpectrum(g);
    spec.addDetectorID(
        g + 1); // Legacy comptibilty: Used to be default IDs in Workspace2D.
    for (int i = 1; i <= 9; i++)
      spec.addDetectorID(g * 9 + i);
    spec.setSpectrumNo(g + 1); // Match detector ID and spec NO
  }
  return boost::dynamic_pointer_cast<MatrixWorkspace>(retVal);
}

// not strictly creating a workspace, but really helpful to see what one
// contains
void displayDataY(MatrixWorkspace_const_sptr ws) {
  const size_t numHists = ws->getNumberHistograms();
  for (size_t i = 0; i < numHists; ++i) {
    std::cout << "Histogram " << i << " = ";
    const auto &y = ws->y(i);
    for (size_t j = 0; j < y.size(); ++j) {
      std::cout << y[j] << " ";
    }
    std::cout << '\n';
  }
}
void displayData(MatrixWorkspace_const_sptr ws) { displayDataX(ws); }

// not strictly creating a workspace, but really helpful to see what one
// contains
void displayDataX(MatrixWorkspace_const_sptr ws) {
  const size_t numHists = ws->getNumberHistograms();
  for (size_t i = 0; i < numHists; ++i) {
    std::cout << "Histogram " << i << " = ";
    const auto &x = ws->x(i);
    for (size_t j = 0; j < x.size(); ++j) {
      std::cout << x[j] << " ";
    }
    std::cout << '\n';
  }
}

// not strictly creating a workspace, but really helpful to see what one
// contains
void displayDataE(MatrixWorkspace_const_sptr ws) {
  const size_t numHists = ws->getNumberHistograms();
  for (size_t i = 0; i < numHists; ++i) {
    std::cout << "Histogram " << i << " = ";
    const auto &e = ws->e(i);
    for (size_t j = 0; j < e.size(); ++j) {
      std::cout << e[j] << " ";
    }
    std::cout << '\n';
  }
}

// =====================================================================================
/** Utility function to add a TimeSeriesProperty with a name and value
 *
 * @param runInfo :: Run to add to
 * @param name :: property name
 * @param val :: value
 */
void addTSPEntry(Run &runInfo, std::string name, double val) {
  TimeSeriesProperty<double> *tsp;
  tsp = new TimeSeriesProperty<double>(name);
  tsp->addValue("2011-05-24T00:00:00", val);
  runInfo.addProperty(tsp);
}

// =====================================================================================
/** Sets the OrientedLattice in the crystal as an crystal with given lattice
 *lengths, angles of 90 deg
 *
 * @param ws :: workspace to set
 * @param a :: lattice length
 * @param b :: lattice length
 * @param c :: lattice length
 */
void setOrientedLattice(Mantid::API::MatrixWorkspace_sptr ws, double a,
                        double b, double c) {
  auto latt =