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/InstrumentCreationHelper.h"

#include "MantidGeometry/Instrument/DetectorInfo.h"
#include "MantidAPI/IAlgorithm.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/Goniometer.h"
#include "MantidGeometry/Instrument/ParameterMap.h"
#include "MantidGeometry/Instrument/ReferenceFrame.h"
#include "MantidGeometry/Objects/ShapeFactory.h"
#include "MantidHistogramData/LinearGenerator.h"
#include "MantidKernel/MersenneTwister.h"
#include "MantidKernel/OptionalBool.h"
#include "MantidKernel/TimeSeriesProperty.h"
#include "MantidKernel/UnitFactory.h"
#include "MantidKernel/VectorHelper.h"
#include "MantidIndexing/IndexInfo.h"

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::move(yAxis), std::move(eAxis)};
    return histo;
  } else {
    Points xAxis(yValsSize, LinearGenerator(1, 1));
    Histogram pointsHisto{std::move(xAxis), std::move(yAxis), std::move(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;
}

void addSource(Instrument_sptr &instrument, const V3D &position,
               std::string name = "source") {
  ObjComponent *source = new ObjComponent(name);
  source->setPos(position);
  instrument->add(source);
  instrument->markAsSource(source);
}

void addSample(Instrument_sptr &instrument, const V3D &position,
               std::string name = "some-surface-holder") {
  ObjComponent *sample = new ObjComponent(name);
  sample->setPos(position);
  instrument->add(sample);
  instrument->markAsSamplePos(sample);
}

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);
}

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);
}

void addSlit(ObjComponent *slit, Instrument_sptr &instrument,
             const V3D &position, std::string name = "slit") {
  slit = new ObjComponent(name);
  slit->setPos(position);
  instrument->add(slit);
}

Workspace2D_sptr reflectometryWorkspace(const double startX, const int nSpectra,
                                        const int nBins, const int 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.
 */
Workspace2D_sptr create2DWorkspaceWithReflectometryInstrument(
    double startX, int nSpectra, int nBins, int deltaX, V3D slit1Pos,
    V3D slit2Pos, double vg1, double vg2) {
  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);
  addDetector(instrument, V3D(20, (20 - 15), 0), 2);

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

  ObjComponent *slit1 = nullptr, *slit2 = nullptr;

  slit1 = new ObjComponent("slit1");
  slit1->setPos(slit1Pos);
  instrument->add(slit1);
  slit2 = new ObjComponent("slit2");
  slit2->setPos(slit2Pos);
  instrument->add(slit2);

  /*
  addSlit(slit1, instrument, slit1Pos, "slit1");
  addSlit(slit2, instrument, slit2Pos, "slit2");
  */

  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)
*/
Workspace2D_sptr create2DWorkspaceWithReflectometryInstrumentMultiDetector(
    double startX, double detSize, int nSpectra, int nBins, int 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);
  // addSlit(instrument, V3D(4, 3, 0), "slit1");
  // addSlit(instrument, V3D(5, 2, 0), "slit2");
  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 =
      Mantid::Kernel::make_unique<OrientedLattice>(a, b, c, 90., 90., 90.);
  ws->mutableSample().setOrientedLattice(latt.release());
}

// =====================================================================================
/** Create a default universal goniometer and set its angles
 *
 * @param ws :: workspace to set
 * @param phi :: +Y rotation angle (deg)
 * @param chi :: +X rotation angle (deg)
 * @param omega :: +Y rotation angle (deg)
 */
void setGoniometer(Mantid::API::MatrixWorkspace_sptr ws, double phi, double chi,
                   double omega) {
  addTSPEntry(ws->mutableRun(), "phi", phi);
  addTSPEntry(ws->mutableRun(), "chi", chi);
  addTSPEntry(ws->mutableRun(), "omega", omega);
  Mantid::Geometry::Goniometer gm;
  gm.makeUniversalGoniometer();
  ws->mutableRun().setGoniometer(gm, true);
}

//
Mantid::API::MatrixWorkspace_sptr
createProcessedWorkspaceWithCylComplexInstrument(size_t numPixels,
                                                 size_t numBins,
                                                 bool has_oriented_lattice) {
  size_t rHist = static_cast<size_t>(std::sqrt(static_cast<double>(numPixels)));
  while (rHist * rHist < numPixels)
    rHist++;

  Mantid::API::MatrixWorkspace_sptr ws =
      createGroupedWorkspace2DWithRingsAndBoxes(rHist, 10, 0.1);
  auto pAxis0 = new NumericAxis(numBins);
  for (size_t i = 0; i < numBins; i++) {
    double dE = -1.0 + static_cast<double>(i) * 0.8;
    pAxis0->setValue(i, dE);
  }
  pAxis0->setUnit("DeltaE");
  ws->replaceAxis(0, pAxis0);
  if (has_oriented_lattice) {
    auto latt =
        Mantid::Kernel::make_unique<OrientedLattice>(1, 1, 1, 90., 90., 90.);
    ws->mutableSample().setOrientedLattice(latt.release());

    addTSPEntry(ws->mutableRun(), "phi", 0);
    addTSPEntry(ws->mutableRun(), "chi", 0);
    addTSPEntry(ws->mutableRun(), "omega", 0);
    Mantid::Geometry::Goniometer gm;
    gm.makeUniversalGoniometer();
    ws->mutableRun().setGoniometer(gm, true);
  }