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 "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"
Gigg, Martyn Anthony
committed
#include "MantidGeometry/Instrument/Detector.h"
#include "MantidGeometry/Instrument/ParameterMap.h"
#include "MantidGeometry/Instrument/ReferenceFrame.h"
#include "MantidGeometry/Instrument/Component.h"
Gigg, Martyn Anthony
committed
#include "MantidGeometry/Objects/ShapeFactory.h"
#include "MantidGeometry/Crystal/OrientedLattice.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;
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);
}
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
Workspace2D_sptr Create1DWorkspaceRand(int size) {
MantidVec y1(size);
MantidVec e1(size);
Federico Montesino Pouzols
committed
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;
Gigg, Martyn Anthony
committed
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;
Janik Zikovsky
committed
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;
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,
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);
Gigg, Martyn Anthony
committed
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);
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;
}
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);
}
Gigg, Martyn Anthony
committed
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);
}
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 += "<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);
Gigg, Martyn Anthony
committed
}
workspace->setInstrument(instrument);
Gigg, Martyn Anthony
committed
}
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);
Gigg, Martyn Anthony
committed
}
}
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,
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());
Gigg, Martyn Anthony
committed
}
Gigg, Martyn Anthony
committed
/** Create a 2D workspace with this many histograms and bins.
*/
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);
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,
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);
/**
* 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;
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)
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);
//================================================================================================================
/** 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
MantidVec x = ws->readX(0);
auto ax0 = new NumericAxis(x.size());
ax0->setUnit("dSpacing");
for (size_t i = 0; i < x.size(); 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);
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();
ws->getSpectrum(wi).clear(true);
ws->getSpectrum(wi).setDetectorID(detectorIds[wi]);
/**
* 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");
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());
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);
workspace->getSpectrum(i).clearDetectorIDs();
workspace->getSpectrum(i).addDetectorID(det->getID());
Gigg, Martyn Anthony
committed
}
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,
double error) {
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,
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");
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,
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);
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 * 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 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] << " ";
Gigg, Martyn Anthony
committed
}
Gigg, Martyn Anthony
committed
}
}
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] << " ";
Gigg, Martyn Anthony
committed
}
Gigg, Martyn Anthony
committed
}
Gigg, Martyn Anthony
committed
// 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] << " ";
}
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,
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());
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
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);
}
// 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;
for (it = detector_map.begin(); it != detector_map.end(); ++it) {
if (!it->second->isMonitor()) {
auto &spec = outputWS->getSpectrum(workspaceIndex);
spec.addDetectorID(it->first);
// Start the spectrum number at 1
spec.setSpectrumNo(specnum_t(workspaceIndex + 1));
Janik Zikovsky
committed
1012
1013
1014
1015
1016
1017
1018
1019
1020
1021
1022
1023
1024
1025
1026
1027
1028
1029
1030
1031
1032
1033
1034
1035
1036
1037
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
MantidVec qbins;
qbins.push_back(0.0);
qbins.push_back(1.0);
qbins.push_back(4.0);
MantidVec 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 x-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);
1047
1048
1049
1050
1051
1052
1053
1054
1055
1056
1057
1058
1059
1060
1061
1062
1063
1064
1065
1066
1067
1068
1069
1070
1071
1072
1073
1074
1075
1076
1077
1078
1079
1080
1081
1082
1083
1084
1085
1086
1087
1088
1089
1090
1091
1092
1093
1094
1095
1096
1097
1098
1099
1100
1101
1102
1103
1104
1105
1106
1107
1108
1109
1110
1111
1112
1113
1114
1115
1116
1117
1118
1119
1120
1121
1122
1123
1124
1125
1126
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
outputWS->dataY(0)[1] = 2.0;
outputWS->dataY(0)[2] = 3.0;
outputWS->dataY(0)[3] = 3.0;
outputWS->dataY(0)[4] = 2.0;
outputWS->dataE(0)[1] = 2.0;
outputWS->dataE(0)[2] = 3.0;
outputWS->dataE(0)[3] = 3.0;
outputWS->dataE(0)[4] = 2.0;
outputWS->dataF(0)[1] = 2.0;
outputWS->dataF(0)[2] = 3.0;
outputWS->dataF(0)[3] = 3.0;
outputWS->dataF(0)[4] = 1.0;
// Q bin #2
outputWS->dataY(1)[1] = 1.0;
outputWS->dataY(1)[2] = 3.0;
outputWS->dataY(1)[3] = 3.0;
outputWS->dataY(1)[4] = 2.0;
outputWS->dataY(1)[5] = 2.0;
outputWS->dataE(1)[1] = 1.0;
outputWS->dataE(1)[2] = 3.0;
outputWS->dataE(1)[3] = 3.0;
outputWS->dataE(1)[4] = 2.0;
outputWS->dataE(1)[5] = 2.0;
outputWS->dataF(1)[1] = 1.0;
outputWS->dataF(1)[2] = 3.0;
outputWS->dataF(1)[3] = 3.0;
outputWS->dataF(1)[4] = 1.0;
outputWS->dataF(1)[5] = 2.0;
// Q bin #3
outputWS->dataY(2)[1] = 1.0;
outputWS->dataY(2)[2] = 2.0;
outputWS->dataY(2)[3] = 3.0;
outputWS->dataY(2)[4] = 1.0;
outputWS->dataE(2)[1] = 1.0;
outputWS->dataE(2)[2] = 2.0;
outputWS->dataE(2)[3] = 3.0;
outputWS->dataE(2)[4] = 1.0;
outputWS->dataF(2)[1] = 1.0;
outputWS->dataF(2)[2] = 2.0;
outputWS->dataF(2)[3] = 2.0;
outputWS->dataF(2)[4] = 1.0;
// Q bin #4
outputWS->dataY(3)[0] = 1.0;
outputWS->dataY(3)[1] = 2.0;
outputWS->dataY(3)[2] = 3.0;
outputWS->dataY(3)[3] = 2.0;
outputWS->dataY(3)[4] = 1.0;
outputWS->dataE(3)[0] = 1.0;
outputWS->dataE(3)[1] = 2.0;
outputWS->dataE(3)[2] = 3.0;
outputWS->dataE(3)[3] = 2.0;
outputWS->dataE(3)[4] = 1.0;
outputWS->dataF(3)[0] = 1.0;
outputWS->dataF(3)[1] = 2.0;
outputWS->dataF(3)[2] = 3.0;
outputWS->dataF(3)[3] = 2.0;
outputWS->dataF(3)[4] = 1.0;
outputWS->dataF(3)[5] = 1.0;
// Set representation
outputWS->finalize();
// Make errors squared rooted
for (int i = 0; i < numHist; ++i) {
for (int j = 0; j < numX - 1; ++j) {
outputWS->dataE(i)[j] = std::sqrt(outputWS->dataE(i)[j]);
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);
1159
1160
1161
1162
1163
1164
1165
1166
1167
1168
1169
1170
1171
1172
1173
1174
1175
1176
1177
1178
1179
1180
1181
1182
1183
1184
1185
1186
1187
1188
1189
1190
1191
1192
1193
1194
1195
1196
1197
1198
1199
1200
// 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");
}
1207
1208
1209
1210
1211
1212
1213
1214
1215
1216
1217
1218
1219
1220
1221
1222
1223
1224
1225
1226
1227
1228
1229
1230
1231
1232
1233
1234
1235
1236
1237
1238
1239
1240
1241
1242
1243
1244
1245
1246
1247
1248
1249
1250
1251
1252
// 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);
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();
// get detector or detector group which corresponds to the spectra i
Geometry::IDetector_const_sptr spDet;
try {
spDet = inputWS->getDetector(i);
} catch (Kernel::Exception::NotFoundError &) {
continue;
}
// Check that we aren't dealing with monitor...
if (spDet->isMonitor())
continue;
// calculate the requested values;
sp2detMap[i] = liveDetectorsCount;
detId[liveDetectorsCount] = int32_t(spDet->getID());
detIDMap[liveDetectorsCount] = i;
L2[liveDetectorsCount] = spDet->getDistance(*sample);
double polar = inputWS->detectorTwoTheta(*spDet);
double azim = spDet->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);