Newer
Older
/*********************************************************************************
* 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/MatrixWorkspace.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"
Gigg, Martyn Anthony
committed
#include "MantidGeometry/Instrument/Detector.h"
#include "MantidGeometry/Instrument/DetectorInfo.h"
#include "MantidGeometry/Instrument/Goniometer.h"
Gigg, Martyn Anthony
committed
#include "MantidGeometry/Instrument/ParameterMap.h"
#include "MantidGeometry/Instrument/ReferenceFrame.h"
Gigg, Martyn Anthony
committed
#include "MantidGeometry/Objects/ShapeFactory.h"
#include "MantidHistogramData/LinearGenerator.h"
#include "MantidIndexing/IndexInfo.h"
#include "MantidKernel/MersenneTwister.h"
Janik Zikovsky
committed
#include "MantidKernel/TimeSeriesProperty.h"
#include "MantidKernel/UnitFactory.h"
#include "MantidKernel/VectorHelper.h"
#include "MantidKernel/make_unique.h"
#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;
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);
}
Gigg, Martyn Anthony
committed
/**
* Deletes a workspace
* @param name :: The name of the workspace
*/
void removeWS(const std::string &name) {
Mantid::API::AnalysisDataService::Instance().remove(name);
}
Gigg, Martyn Anthony
committed
/**
* 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) {
Federico Montesino Pouzols
committed
MersenneTwister randomGen(DateAndTime::getCurrentTime().nanoseconds(), 0,
std::numeric_limits<int>::max());
Federico Montesino Pouzols
committed
auto randFunc = [&randomGen] { return randomGen.nextValue(); };
Counts counts(size, randFunc);
CountStandardDeviations errorVals(size, randFunc);
Federico Montesino Pouzols
committed
auto retVal = boost::make_shared<Workspace2D>();
retVal->initialize(1, createHisto(isHisto, counts, errorVals));
Gigg, Martyn Anthony
committed
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));
Janik Zikovsky
committed
Workspace2D_sptr create1DWorkspaceConstantWithXerror(int size, double value,
double error,
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));
Gigg, Martyn Anthony
committed
Workspace2D_sptr create2DWorkspace(int nhist, int numBoundaries) {
return create2DWorkspaceBinned(nhist, numBoundaries);
Janik Zikovsky
committed
/** 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,
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);
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);
Gigg, Martyn Anthony
committed
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);
Janik Zikovsky
committed
}
retVal = maskSpectra(retVal, maskedWorkspaceIndices);
return retVal;
}
Gigg, Martyn Anthony
committed
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;
}
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);
}
Gigg, Martyn Anthony
committed
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);
}
Janik Zikovsky
committed
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);
Gigg, Martyn Anthony
committed
std::string xmlShape = "<sphere id=\"shape\"> ";
xmlShape += R"(<centre x="0.0" y="0.0" z="0.0" /> )";
xmlShape += "<radius val=\"0.05\" /> ";
xmlShape += "</sphere>";
xmlShape += "<algebra val=\"shape\" /> ";
ShapeFactory sFactory;
auto 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);
Gigg, Martyn Anthony
committed
}
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));
Gigg, Martyn Anthony
committed
}
auto &spectrumInfo = workspace->mutableSpectrumInfo();
for (const auto index : maskedWorkspaceIndices)
spectrumInfo.setMasked(index, true);
Gigg, Martyn Anthony
committed
/**
* Create a group with nEntries. It is added to the ADS with the given stem
*/
WorkspaceGroup_sptr createWorkspaceGroup(int nEntries, int nHist, int nBins,
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());
Gigg, Martyn Anthony
committed
}
Gigg, Martyn Anthony
committed
/** Create a 2D workspace with this many histograms and bins.
Workspace2D_sptr create2DWorkspaceBinned(int nhist, int numVals, double x0,
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);
Gigg, Martyn Anthony
committed
}
Gigg, Martyn Anthony
committed
/** Create a 2D workspace with this many histograms and bins. The bins are
* assumed to be non-uniform and given by the input array
Workspace2D_sptr create2DWorkspaceBinned(int nhist, const int numBoundaries,
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);
/**
* 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;
Janik Zikovsky
committed
}
}
Janik Zikovsky
committed
//================================================================================================================
/**
* 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)
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);
//================================================================================================================
/*
* 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));
//================================================================================================================
/** 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();
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;
ComponentCreationHelper::createTestInstrumentCylindrical(
numBanks, srcPos, samplePos, 0.0025, 0.005);
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();
ws->getSpectrum(wi).clear(true);
ws->getSpectrum(wi).setDetectorID(detectorIds[wi]);
/** 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
ObjComponent *addComponent(Instrument_sptr &instrument, const V3D &position,
ObjComponent *component = new ObjComponent(name);
component->setPos(position);
instrument->add(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,
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,
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,
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,
// 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");
/**
* 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 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,
const int nSpectra, const int nBins, const double deltaX) {
Instrument_sptr instrument = boost::make_shared<Instrument>();
instrument->setReferenceFrame(boost::make_shared<ReferenceFrame>(
PointingAlong::Y, PointingAlong::X, Handedness::Left, "0,0,0"));
addSource(instrument, sourcePos, "source");
addMonitor(instrument, monitorPos, 1, "Monitor");
addSample(instrument, samplePos, "some-surface-holder");
addDetector(instrument, detectorPos, 2, "point-detector");
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);
* 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,
Instrument_sptr instrument = boost::make_shared<Instrument>();
instrument->setReferenceFrame(boost::make_shared<ReferenceFrame>(
PointingAlong::Y /*up*/, PointingAlong::X /*along*/, Handedness::Left,
"0,0,0"));
addSource(instrument, V3D(0, 0, 0), "source");
addSample(instrument, V3D(15, 0, 0), "some-surface-holder");
addMonitor(instrument, V3D(14, 0, 0), 1, "Monitor");
// 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,
"point-detector"); // offset below centre
addDetector(instrument, V3D(posX, posY, 0), 3, "point-detector"); // at centre
addDetector(instrument, V3D(posX, posY + detSize, 0), 4,
"point-detector"); // 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);
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, "source");
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());
workspace->getSpectrum(i).addDetectorID(i);
Gigg, Martyn Anthony
committed
}
workspace->setInstrument(instrument);
Gigg, Martyn Anthony
committed
//================================================================================================================
WorkspaceSingleValue_sptr createWorkspaceSingleValue(double value) {
return boost::make_shared<WorkspaceSingleValue>(value, sqrt(value));
Gigg, Martyn Anthony
committed
WorkspaceSingleValue_sptr createWorkspaceSingleValueWithError(double value,
return boost::make_shared<WorkspaceSingleValue>(value, error);
Gigg, Martyn Anthony
committed
/** Perform some finalization on event workspace stuff */
void eventWorkspace_Finalize(EventWorkspace_sptr ew) {
// get a proton charge
ew->mutableRun().integrateProtonCharge();
}
Gigg, Martyn Anthony
committed
/** Create event workspace with:
* 500 pixels
* 1000 histogrammed bins.
*/
EventWorkspace_sptr createEventWorkspace() {
return createEventWorkspace(500, 1001, 100, 1000);
Gigg, Martyn Anthony
committed
/** 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) {
Gigg, Martyn Anthony
committed
// add one to the number of bins as this is histogram
numBins++;
auto retVal = boost::make_shared<EventWorkspace>();
retVal->initialize(numPixels, 1, 1);
Janik Zikovsky
committed
// 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));
}
Gigg, Martyn Anthony
committed
}
Gigg, Martyn Anthony
committed
}
Gigg, Martyn Anthony
committed
retVal->setAllX(BinEdges(numBins, LinearGenerator(x0, binDelta)));
Gigg, Martyn Anthony
committed
Gigg, Martyn Anthony
committed
// =====================================================================================
/** 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 (int i = 0; i < numBins; i++)
retVal->getSpectrum(g) += TofEvent((i + 0.5) * binDelta, 1);
retVal->getSpectrum(g).addDetectorID(det);
}
}
Gigg, Martyn Anthony
committed
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)));
// =====================================================================================
/** 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,
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");
Federico Montesino Pouzols
committed
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
Federico Montesino Pouzols
committed
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);
// =====================================================================================
/** 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,
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 * 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.
spec.addDetectorID(g * 9 + i);
spec.setSpectrumNo(g + 1); // Match detector ID and spec NO
Gigg, Martyn Anthony
committed
}
return boost::dynamic_pointer_cast<MatrixWorkspace>(retVal);
}
Gigg, Martyn Anthony
committed
// not strictly creating a workspace, but really helpful to see what one
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] << " ";
Gigg, Martyn Anthony
committed
}
Gigg, Martyn Anthony
committed
}
void displayData(MatrixWorkspace_const_sptr ws) { displayDataX(ws); }
// not strictly creating a workspace, but really helpful to see what one
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] << " ";
Gigg, Martyn Anthony
committed
}
Gigg, Martyn Anthony
committed
}
Gigg, Martyn Anthony
committed
// not strictly creating a workspace, but really helpful to see what one
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] << " ";
Janik Zikovsky
committed
}
Janik Zikovsky
committed
// =====================================================================================
/** 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,
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,
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());
Janik Zikovsky
committed
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);
Janik Zikovsky
committed
}
/// 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++) {
// we just made (in createCylInstrumentWithDetInGivenPosisions) det ID-s to
// start from 1
// and this is absolutely different nummer, corresponding to det ID just by
// chance ? -- some uncertainties remain
// spec->setSpectrumNo(g+1);
// spec->addDetectorID(g*9);
// spec->setSpectrumNo(g+1); // Match detector ID and spec NO
const double dE = (Emax - Emin) / static_cast<double>(numBins);
for (size_t j = 0; j < numPixels; j++) {
std::vector<double> E_transfer;
E_transfer.reserve(numBins);
for (size_t i = 0; i <= numBins; i++) {
E_transfer.push_back(Emin + static_cast<double>(i) * dE);
// set axis, correspondent to the X-values
auto pAxis0 = new NumericAxis(numBins);
const auto &E_transfer = ws->x(0);
for (size_t i = 0; i < numBins; i++) {
double E = 0.5 * (E_transfer[i] + E_transfer[i + 1]);
pAxis0->setValue(i, E);
}
// 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);
Janik Zikovsky
committed
/*
* 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",
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;
const auto &detectorInfo = outputWS->detectorInfo();
for (it = detector_map.begin(); it != detector_map.end(); ++it) {
if (!detectorInfo.isMonitor(detectorInfo.indexOf(it->first))) {
auto &spec = outputWS->getSpectrum(workspaceIndex);
spec.addDetectorID(it->first);
// Start the spectrum number at 1
spec.setSpectrumNo(specnum_t(workspaceIndex + 1));
Janik Zikovsky
committed
RebinnedOutput_sptr createRebinnedOutputWorkspace() {
RebinnedOutput_sptr outputWS =
Mantid::DataObjects::RebinnedOutput_sptr(new RebinnedOutput());
// outputWS->setName("rebinTest");
Mantid::API::AnalysisDataService::Instance().add("rebinTest", outputWS);
// Set Q ('y') axis binning
std::vector<double> qbins{0.0, 1.0, 4.0};
std::vector<double> qaxis;
const int numY =
static_cast<int>(VectorHelper::createAxisFromRebinParams(qbins, qaxis));
// Initialize the workspace
const int numHist = numY - 1;
const int numX = 7;
outputWS->initialize(numHist, numX, numX - 1);
// Set the normal units
outputWS->getAxis(0)->unit() = UnitFactory::Instance().create("DeltaE");
outputWS->setYUnit("Counts");
outputWS->setTitle("Empty_Title");
// Create the i-axis for histogramming.
HistogramData::BinEdges x1{-3.0, -2.0, -1.0, 0.0, 1.0, 2.0, 3.0};
// Create a numeric axis to replace the default vertical one
Axis *const verticalAxis = new NumericAxis(numY);
outputWS->replaceAxis(1, verticalAxis);
// Now set the axis values
for (int i = 0; i < numHist; ++i) {
outputWS->setBinEdges(i, x1);
verticalAxis->setValue(i, qaxis[i]);
}
// One more to set on the 'y' axis
verticalAxis->setValue(numHist, qaxis[numHist]);
// Set the 'y' axis units
verticalAxis->unit() = UnitFactory::Instance().create("MomentumTransfer");
verticalAxis->title() = "|Q|";
// Set the X axis title (for conversion to MD)
outputWS->getAxis(0)->title() = "Energy transfer";
// Now, setup the data
// Q bin #1
// Populates destination starting at index 1 with the following data
// e.g. y(0)[1] = 2.0, y(0)[2] = 3.0 ..etc. as the starting index is 1
// if you change the values in the line below please update this comment!
populateWsWithInitList(outputWS->mutableY(0), 1, {2.0, 3.0, 3.0, 2.0});
populateWsWithInitList(outputWS->mutableE(0), 1, {2.0, 3.0, 3.0, 2.0});
populateWsWithInitList(outputWS->dataF(0), 1, {2.0, 3.0, 3.0, 1.0});
populateWsWithInitList(outputWS->mutableY(1), 1, {1.0, 3.0, 3.0, 2.0, 2.0});
populateWsWithInitList(outputWS->mutableE(1), 1, {1.0, 3.0, 3.0, 2.0, 2.0});
populateWsWithInitList(outputWS->dataF(1), 1, {1.0, 3.0, 3.0, 1.0, 2.0});
populateWsWithInitList(outputWS->mutableY(2), 1, {1.0, 2.0, 3.0, 1.0});
populateWsWithInitList(outputWS->mutableE(2), 1, {1.0, 2.0, 3.0, 1.0});
populateWsWithInitList(outputWS->dataF(2), 1, {1.0, 2.0, 2.0, 1.0});
populateWsWithInitList(outputWS->mutableY(3), 0, {1.0, 2.0, 3.0, 2.0, 1.0});
populateWsWithInitList(outputWS->mutableE(3), 0, {1.0, 2.0, 3.0, 2.0, 1.0});
populateWsWithInitList(outputWS->dataF(3), 0, {1.0, 2.0, 3.0, 2.0, 1.0, 1.0});
// Set representation
outputWS->finalize();
// Make errors squared rooted
for (int i = 0; i < numHist; ++i) {
auto &mutableE = outputWS->mutableE(i);
for (int j = 0; j < numX - 1; ++j) {
mutableE[j] = std::sqrt(mutableE[j]);
/**
* Populates the destination array (usually a mutable histogram)
* starting at the index specified with the doubles provided in an
* initializer list. Note the caller is responsible for ensuring
* the destination has capacity for startingIndex + size(initializer list)
* number of values
*
* @param destination :: The array to populate with data
* @param startingIndex :: The index to start populating data at
* @param values :: The initializer list to populate the array with
* starting at the index specified
*/
template <typename T>
void populateWsWithInitList(T &destination, size_t startingIndex,
const std::initializer_list<double> &values) {
size_t index = 0;
for (const double val : values) {
destination[startingIndex + index] = val;
index++;
}
}
Mantid::DataObjects::PeaksWorkspace_sptr
createPeaksWorkspace(const int numPeaks, const bool createOrientedLattice) {
auto peaksWS = boost::make_shared<PeaksWorkspace>();
Instrument_sptr inst =
ComponentCreationHelper::createTestInstrumentRectangular2(1, 10);
peaksWS->setInstrument(inst);
for (int i = 0; i < numPeaks; ++i) {
Peak peak(inst, i, i + 0.5);
peaksWS->addPeak(peak);
if (createOrientedLattice) {
Mantid::Geometry::OrientedLattice lattice;
peaksWS->mutableSample().setOrientedLattice(&lattice);
}
/** helper method to create preprocessed detector's table workspace */
boost::shared_ptr<DataObjects::TableWorkspace>
createTableWorkspace(const API::MatrixWorkspace_const_sptr &inputWS) {
const size_t nHist = inputWS->getNumberHistograms();
// set the target workspace
auto targWS = boost::make_shared<TableWorkspace>(nHist);
1340
1341
1342
1343
1344
1345
1346
1347
1348
1349
1350
1351
1352
1353
1354
1355
1356
1357
1358
1359
1360
1361
1362
1363
1364
1365
1366
1367
1368
1369
1370
1371
1372
1373
1374
1375
1376
1377
1378
1379
1380
1381
// detectors positions
if (!targWS->addColumn("V3D", "DetDirections"))
throw(std::runtime_error("Can not add column DetDirectrions"));
// sample-detector distance;
if (!targWS->addColumn("double", "L2"))
throw(std::runtime_error("Can not add column L2"));
// Diffraction angle
if (!targWS->addColumn("double", "TwoTheta"))
throw(std::runtime_error("Can not add column TwoTheta"));
if (!targWS->addColumn("double", "Azimuthal"))
throw(std::runtime_error("Can not add column Azimuthal"));
// the detector ID;
if (!targWS->addColumn("int", "DetectorID"))
throw(std::runtime_error("Can not add column DetectorID"));
// stores spectra index which corresponds to a valid detector index;
if (!targWS->addColumn("size_t", "detIDMap"))
throw(std::runtime_error("Can not add column detIDMap"));
// stores detector index which corresponds to the workspace index;
if (!targWS->addColumn("size_t", "spec2detMap"))
throw(std::runtime_error("Can not add column spec2detMap"));
// will see about that
// sin^2(Theta)
// std::vector<double> SinThetaSq;
//,"If the detectors were actually processed from real instrument or generated
// for some fake one ");
return targWS;
}
/** method does preliminary calculations of the detectors positions to convert
results into k-dE space ;
and places the resutls into static cash to be used in subsequent calls to this
algorithm */
void processDetectorsPositions(const API::MatrixWorkspace_const_sptr &inputWS,
DataObjects::TableWorkspace_sptr &targWS,
double Ei) {
Geometry::Instrument_const_sptr instrument = inputWS->getInstrument();
//
Geometry::IComponent_const_sptr source = instrument->getSource();
Geometry::IComponent_const_sptr sample = instrument->getSample();
if ((!source) || (!sample)) {
throw Kernel::Exception::InstrumentDefinitionError(
"Instrubment not sufficiently defined: failed to get source and/or "
"sample");
}
1388
1389
1390
1391
1392
1393
1394
1395
1396
1397
1398
1399
1400
1401
1402
1403
1404
1405
1406
1407
1408
1409
1410
1411
1412
1413
1414
1415
1416
1417
1418
// L1
try {
double L1 = source->getDistance(*sample);
targWS->logs()->addProperty<double>("L1", L1, true);
} catch (Kernel::Exception::NotFoundError &) {
throw Kernel::Exception::InstrumentDefinitionError(
"Unable to calculate source-sample distance for workspace",
inputWS->getTitle());
}
// Instrument name
std::string InstrName = instrument->getName();
targWS->logs()->addProperty<std::string>("InstrumentName", InstrName, true);
targWS->logs()->addProperty<bool>("FakeDetectors", false, true);
targWS->logs()->addProperty<double>("Ei", Ei, true); //"Incident energy for
// Direct or Analysis
// energy for indirect
// instrument");
// get access to the workspace memory
auto &sp2detMap = targWS->getColVector<size_t>("spec2detMap");
auto &detId = targWS->getColVector<int32_t>("DetectorID");
auto &detIDMap = targWS->getColVector<size_t>("detIDMap");
auto &L2 = targWS->getColVector<double>("L2");
auto &TwoTheta = targWS->getColVector<double>("TwoTheta");
auto &Azimuthal = targWS->getColVector<double>("Azimuthal");
auto &detDir = targWS->getColVector<Kernel::V3D>("DetDirections");
//// progress messave appearence
size_t nHist = targWS->rowCount();
//// Loop over the spectra
uint32_t liveDetectorsCount(0);
const auto &spectrumInfo = inputWS->spectrumInfo();
for (size_t i = 0; i < nHist; i++) {
sp2detMap[i] = std::numeric_limits<size_t>::quiet_NaN();
detId[i] = std::numeric_limits<int32_t>::quiet_NaN();
detIDMap[i] = std::numeric_limits<size_t>::quiet_NaN();
L2[i] = std::numeric_limits<double>::quiet_NaN();
TwoTheta[i] = std::numeric_limits<double>::quiet_NaN();
Azimuthal[i] = std::numeric_limits<double>::quiet_NaN();
if (!spectrumInfo.hasDetectors(i) || spectrumInfo.isMonitor(i))
// calculate the requested values;
sp2detMap[i] = liveDetectorsCount;
detId[liveDetectorsCount] = int32_t(spectrumInfo.detector(i).getID());
detIDMap[liveDetectorsCount] = i;
L2[liveDetectorsCount] = spectrumInfo.l2(i);
double polar = spectrumInfo.twoTheta(i);
double azim = spectrumInfo.detector(i).getPhi();
TwoTheta[liveDetectorsCount] = polar;
Azimuthal[liveDetectorsCount] = azim;
double sPhi = sin(polar);
double ez = cos(polar);
double ex = sPhi * cos(azim);
double ey = sPhi * sin(azim);
detDir[liveDetectorsCount].setX(ex);
detDir[liveDetectorsCount].setY(ey);
detDir[liveDetectorsCount].setZ(ez);
// double sinTheta=sin(0.5*polar);
// this->SinThetaSq[liveDetectorsCount] = sinTheta*sinTheta;
targWS->logs()->addProperty<uint32_t>(
"ActualDetectorsNum", liveDetectorsCount,
true); //,"The actual number of detectors receivinv signal");
}
boost::shared_ptr<Mantid::DataObjects::TableWorkspace>
buildPreprocessedDetectorsWorkspace(Mantid::API::MatrixWorkspace_sptr ws) {
Mantid::DataObjects::TableWorkspace_sptr DetPos = createTableWorkspace(ws);
double Ei = ws->run().getPropertyValueAsType<double>("Ei");
processDetectorsPositions(ws, DetPos, Ei);
return DetPos;
}
void create2DAngles(std::vector<double> &L2, std::vector<double> &polar,
std::vector<double> &azim, size_t nPolar, size_t nAzim,
double polStart, double polEnd, double azimStart,
double azimEnd) {
size_t nDet = nPolar * nAzim;
L2.resize(nDet, 10);
polar.resize(nDet);
azim.resize(nDet);
double dPolar = (polEnd - polStart) / static_cast<double>(nDet - 1);
double dAzim = (azimEnd - azimEnd) / static_cast<double>(nDet - 1);
for (size_t i = 0; i < nPolar; i++) {
for (size_t j = 0; j < nAzim; j++) {
polar[i * nPolar + j] = polStart + dPolar * static_cast<double>(i);
azim[i * nPolar + j] = azimStart + dAzim * static_cast<double>(j);
ITableWorkspace_sptr
createEPPTableWorkspace(const std::vector<EPPTableRow> &rows) {
ITableWorkspace_sptr ws = boost::make_shared<TableWorkspace>(rows.size());
auto wsIndexColumn = ws->addColumn("int", "WorkspaceIndex");
auto centreColumn = ws->addColumn("double", "PeakCentre");
auto centreErrorColumn = ws->addColumn("double", "PeakCentreError");
auto sigmaColumn = ws->addColumn("double", "Sigma");
auto sigmaErrorColumn = ws->addColumn("double", "SigmaError");
auto heightColumn = ws->addColumn("double", "Height");
auto heightErrorColumn = ws->addColumn("double", "HeightError");
auto chiSqColumn = ws->addColumn("double", "chiSq");
auto statusColumn = ws->addColumn("str", "FitStatus");
for (size_t i = 0; i != rows.size(); ++i) {
if (row.workspaceIndex < 0) {
wsIndexColumn->cell<int>(i) = static_cast<int>(i);
} else {
wsIndexColumn->cell<int>(i) = row.workspaceIndex;
}
centreColumn->cell<double>(i) = row.peakCentre;
centreErrorColumn->cell<double>(i) = row.peakCentreError;
sigmaColumn->cell<double>(i) = row.sigma;
sigmaErrorColumn->cell<double>(i) = row.sigmaError;
heightColumn->cell<double>(i) = row.height;
heightErrorColumn->cell<double>(i) = row.heightError;
chiSqColumn->cell<double>(i) = row.chiSq;
statusColumn->cell<std::string>(i) =
row.fitStatus == EPPTableRow::FitStatus::SUCCESS ? "success" : "failed";
}
return ws;
Gigg, Martyn Anthony
committed
}
} // namespace WorkspaceCreationHelper