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 "MantidHistogramData/LinearGenerator.h"
#include "MantidAPI/Run.h"
#include "MantidAPI/IAlgorithm.h"
#include "MantidAPI/Algorithm.h"
#include "MantidAPI/Sample.h"
#include "MantidAPI/SpectraAxis.h"
#include "MantidAPI/NumericAxis.h"
#include "MantidAPI/WorkspaceGroup.h"
#include "MantidDataObjects/PeaksWorkspace.h"
#include "MantidGeometry/Instrument/Detector.h"
#include "MantidGeometry/Instrument/ParameterMap.h"
#include "MantidGeometry/Instrument/ReferenceFrame.h"
#include "MantidGeometry/Instrument/Component.h"
#include "MantidGeometry/Objects/ShapeFactory.h"
#include "MantidGeometry/Crystal/OrientedLattice.h"
#include "MantidKernel/MersenneTwister.h"
#include "MantidKernel/TimeSeriesProperty.h"
#include "MantidKernel/UnitFactory.h"
#include "MantidKernel/VectorHelper.h"
#include "MantidKernel/make_unique.h"

#include <cmath>
#include <sstream>

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;

MockAlgorithm::MockAlgorithm(size_t nSteps) {
  m_Progress = Mantid::Kernel::make_unique<API::Progress>(this, 0, 1, nSteps);
}

/**
 * @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);
}

Workspace2D_sptr Create1DWorkspaceRand(int size) {
  MantidVec y1(size);
  MantidVec e1(size);

  MersenneTwister randomGen(DateAndTime::getCurrentTime().nanoseconds(), 0,
                            std::numeric_limits<int>::max());
  auto randFunc = [&randomGen] { return randomGen.nextValue(); };

  std::generate(y1.begin(), y1.end(), randFunc);
  std::generate(e1.begin(), e1.end(), randFunc);
  auto retVal = boost::make_shared<Workspace2D>();
  retVal->initialize(1, size, size);
  retVal->setPoints(0, size, LinearGenerator(1.0, 1.0));
  retVal->dataY(0) = y1;
  retVal->dataE(0) = e1;
  return retVal;
}

Workspace2D_sptr Create1DWorkspaceConstant(int size, double value,
                                           double error) {
  MantidVec y1(size, value);
  MantidVec e1(size, error);
  auto retVal = boost::make_shared<Workspace2D>();
  retVal->initialize(1, size, size);
  retVal->setPoints(0, size, LinearGenerator(1.0, 1.0));
  retVal->dataY(0) = y1;
  retVal->dataE(0) = e1;
  return retVal;
}

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

Workspace2D_sptr Create1DWorkspaceFib(int size) {
  MantidVec y1(size);
  MantidVec e1(size);
  std::generate(y1.begin(), y1.end(), FibSeries<double>());
  auto retVal = boost::make_shared<Workspace2D>();
  retVal->initialize(1, size, size);
  retVal->setPoints(0, size, LinearGenerator(1.0, 1.0));
  retVal->dataY(0) = y1;
  retVal->dataE(0) = e1;
  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 wi = 0; wi < nhist; wi++)
    for (int x = 0; x < numBoundaries; x++)
      out->dataY(wi)[x] = wi * 1.0;
  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), shape, nullptr);
      det->setPos(i, i + 1, 1);
      instrument->add(det);
      instrument->markAsDetector(det);
    }
    workspace->setInstrument(instrument);
  }

  ParameterMap &pmap = workspace->instrumentParameters();
  for (int i = 0; i < nhist; ++i) {
    if (maskedWorkspaceIndices.find(i) != maskedWorkspaceIndices.end()) {
      IDetector_const_sptr det = workspace->getDetector(i);
      pmap.addBool(det.get(), "masked", 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 nbins, double x0,
                                         double deltax) {
  BinEdges x(nbins + 1, LinearGenerator(x0, deltax));
  Counts y(nbins, 2);
  CountStandardDeviations e(nbins, M_SQRT2);
  auto retVal = boost::make_shared<Workspace2D>();
  retVal->initialize(nhist, nbins + 1, nbins);
  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++) {
    Mantid::MantidVec &Y = ws->dataY(iSpec);
    Mantid::MantidVec &E = ws->dataE(iSpec);
    for (size_t i = 0; i < Y.size(); i++) {
      Y[i] += noise * randGen.nextValue();
      E[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;
}

//================================================================================================================
/** 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
  MantidVec x = ws->readX(0);
  auto ax0 = new NumericAxis(x.size());
  ax0->setUnit("dSpacing");
  for (size_t i = 0; i < x.size(); i++) {
    ax0->setValue(i, x[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;
}

/**
 * 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.
 */
MatrixWorkspace_sptr
create2DWorkspaceWithReflectometryInstrument(double startX) {
  Instrument_sptr instrument = boost::make_shared<Instrument>();
  instrument->setReferenceFrame(
      boost::make_shared<ReferenceFrame>(Y /*up*/, X /*along*/, Left, "0,0,0"));

  ObjComponent *source = new ObjComponent("source");
  source->setPos(V3D(0, 0, 0));
  instrument->add(source);
  instrument->markAsSource(source);

  Detector *monitor = new Detector("Monitor", 1, nullptr);
  monitor->setPos(14, 0, 0);
  instrument->add(monitor);
  instrument->markAsMonitor(monitor);

  ObjComponent *sample = new ObjComponent("some-surface-holder");
  sample->setPos(V3D(15, 0, 0));
  instrument->add(sample);
  instrument->markAsSamplePos(sample);

  // Where 0.01 is half detector width etc.
  Detector *det = new Detector(
      "point-detector", 2,
      ComponentCreationHelper::createCuboid(0.01, 0.02, 0.03), nullptr);
  det->setPos(20, (20 - sample->getPos().X()), 0);
  instrument->add(det);
  instrument->markAsDetector(det);

  const int nSpectra = 2;
  const int nBins = 100;
  const double deltaX = 2000; // TOF
  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");

  workspace->setInstrument(instrument);
  workspace->getSpectrum(0).setDetectorID(det->getID());
  workspace->getSpectrum(1).setDetectorID(monitor->getID());
  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"));

  ObjComponent *source = new ObjComponent("source");
  source->setPos(sourcePosition);
  instrument->add(source);
  instrument->markAsSource(source);

  ObjComponent *sample = new ObjComponent("sample");
  source->setPos(samplePosition);
  instrument->add(sample);
  instrument->markAsSamplePos(sample);

  workspace->setInstrument(instrument);

  for (int i = 0; i < static_cast<int>(detectorPositions.size()); ++i) {
    std::stringstream buffer;
    buffer << "detector_" << i;
    Detector *det = new Detector(buffer.str(), i, nullptr);
    det->setPos(detectorPositions[i]);
    instrument->add(det);
    instrument->markAsDetector(det);

    // Link it to the workspace
    workspace->getSpectrum(i).clearDetectorIDs();
    workspace->getSpectrum(i).addDetectorID(det->getID());
  }
}

//================================================================================================================
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);
    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 helpfull to see what one
// contains
void DisplayDataY(const MatrixWorkspace_sptr ws) {
  const size_t numHists = ws->getNumberHistograms();
  for (size_t i = 0; i < numHists; ++i) {
    std::cout << "Histogram " << i << " = ";
    for (size_t j = 0; j < ws->blocksize(); ++j) {
      std::cout << ws->readY(i)[j] << " ";
    }
    std::cout << '\n';
  }
}
void DisplayData(const MatrixWorkspace_sptr ws) { DisplayDataX(ws); }

// not strictly creating a workspace, but really helpfull to see what one
// contains
void DisplayDataX(const MatrixWorkspace_sptr ws) {
  const size_t numHists = ws->getNumberHistograms();
  for (size_t i = 0; i < numHists; ++i) {
    std::cout << "Histogram " << i << " = ";
    for (size_t j = 0; j < ws->blocksize(); ++j) {
      std::cout << ws->readX(i)[j] << " ";
    }
    std::cout << '\n';
  }
}

// not strictly creating a workspace, but really helpfull to see what one
// contains
void DisplayDataE(const MatrixWorkspace_sptr ws) {
  const size_t numHists = ws->getNumberHistograms();
  for (size_t i = 0; i < numHists; ++i) {
    std::cout << "Histogram " << i << " = ";
    for (size_t j = 0; j < ws->blocksize(); ++j) {
      std::cout << ws->readE(i)[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);
  }

  return ws;
}

/// Create a workspace with all components needed for inelastic analysis and 3
/// detectors in specific places
/// @param L2        -- the sample to detector flight path
/// @param polar     -- the detector polar angle
/// @param azimutal  -- the detector azimuthal
/// @param numBins   -- the number of histogram bins for the workspace
/// @param Emin      -- minimal energy transfer
/// @param Emax      -- maxinal energy transfer
/// @param Ei        -- input beam energy
Mantid::API::MatrixWorkspace_sptr
createProcessedInelasticWS(const std::vector<double> &L2,
                           const std::vector<double> &polar,
                           const std::vector<double> &azimutal, size_t numBins,
                           double Emin, double Emax, double Ei) {
  // not used but interface needs it
  std::set<int64_t> maskedWorkspaceIndices;
  size_t numPixels = L2.size();

  Mantid::API::MatrixWorkspace_sptr ws =
      Create2DWorkspaceWithValues(uint64_t(numPixels), uint64_t(numBins), true,
                                  maskedWorkspaceIndices, 0, 1, 0.1);

  // detectors at L2, sample at 0 and source at -L2_min
  ws->setInstrument(
      ComponentCreationHelper::createCylInstrumentWithDetInGivenPositions(
          L2, polar, azimutal));

  for (int g = 0; g < static_cast<int>(numPixels); g++) {
    auto &spec = ws->getSpectrum(g);
    // we just made (in createCylInstrumentWithDetInGivenPosisions) det ID-s to
    // start from 1
    spec.setDetectorID(g + 1);
    // and this is absolutely different nummer, corresponding to det ID just by
    // chance ? -- some uncertainties remain
    spec.setSpectrumNo(g + 1);
    // spec->setSpectrumNo(g+1);
    //   spec->addDetectorID(g*9);
    //   spec->setSpectrumNo(g+1); // Match detector ID and spec NO
  }

  double dE = (Emax - Emin) / double(numBins);
  for (size_t j = 0; j < numPixels; j++) {

    MantidVec &E_transfer = ws->dataX(j);
    for (size_t i = 0; i <= numBins; i++) {
      double E = Emin + static_cast<double>(i) * dE;
      E_transfer[i] = E;
    }
  }
  // set axis, correspondent to the X-values
  auto pAxis0 = new NumericAxis(numBins);
  MantidVec &E_transfer = ws->dataX(0);
  for (size_t i = 0; i < numBins; i++) {
    double E = 0.5 * (E_transfer[i] + E_transfer[i + 1]);
    pAxis0->setValue(i, E);
  }

  pAxis0->setUnit("DeltaE");

  ws->replaceAxis(0, pAxis0);

  // define oriented lattice which requested for processed ws
  auto latt =
      Mantid::Kernel::make_unique<OrientedLattice>(1, 1, 1, 90., 90., 90.);
  ws->mutableSample().setOrientedLattice(latt.release());

  ws->mutableRun().addProperty(
      new PropertyWithValue<std::string>("deltaE-mode", "Direct"), true);
  ws->mutableRun().addProperty(new PropertyWithValue<double>("Ei", Ei), true);
  // these properties have to be different -> specific for processed ws, as time
  // now should be reconciled
  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);

  return ws;
}

/*
 * Create an EventWorkspace from a source EventWorkspace.
 * The new workspace should be exactly the same as the source workspace but
 * without any events
 */
Mantid::DataObjects::EventWorkspace_sptr
createEventWorkspace3(Mantid::DataObjects::EventWorkspace_const_sptr sourceWS,
                      std::string wsname, API::Algorithm *alg) {
  UNUSED_ARG(wsname);
  // 1. Initialize:use dummy numbers for arguments, for event workspace it
  // doesn't matter
  Mantid::DataObjects::EventWorkspace_sptr outputWS =
      Mantid::DataObjects::EventWorkspace_sptr(
          new DataObjects::EventWorkspace());
  outputWS->initialize(sourceWS->getInstrument()->getDetectorIDs(true).size(),
                       1, 1);

  // 2. Set the units
  outputWS->getAxis(0)->unit() = UnitFactory::Instance().create("TOF");
  outputWS->setYUnit("Counts");
  outputWS->setTitle("Empty_Title");

  // 3. Add the run_start property:
  int runnumber = sourceWS->getRunNumber();
  outputWS->mutableRun().addProperty("run_number", runnumber);

  std::string runstartstr = sourceWS->run().getProperty("run_start")->value();
  outputWS->mutableRun().addProperty("run_start", runstartstr);

  // 4. Instrument
  Mantid::API::Algorithm_sptr loadInst =
      alg->createChildAlgorithm("LoadInstrument");
  // Now execute the Child Algorithm. Catch and log any error, but don't stop.
  loadInst->setPropertyValue("InstrumentName",
                             sourceWS->getInstrument()->getName());
  loadInst->setProperty<MatrixWorkspace_sptr>("Workspace", outputWS);
  loadInst->setProperty("RewriteSpectraMap",
                        Mantid::Kernel::OptionalBool(true));
  loadInst->executeAsChildAlg();
  // Populate the instrument parameters in this workspace - this works around a
  // bug
  outputWS->populateInstrumentParameters();

  // 6. Build spectrum and event list
  // a) We want to pad out empty pixels.
  detid2det_map detector_map;
  outputWS->getInstrument()->getDetectors(detector_map);

  // b) determine maximum pixel id
  detid2det_map::iterator it;
  detid_t detid_max = 0; // seems like a safe lower bound
  for (it = detector_map.begin(); it != detector_map.end(); ++it)
    if (it->first > detid_max)
      detid_max = it->first;

  // c) Pad all the pixels and Set to zero
  size_t workspaceIndex = 0;
  for (it = detector_map.begin(); it != detector_map.end(); ++it) {
    if (!it->second->isMonitor()) {