31#include <boost/algorithm/string/split.hpp>
42using namespace Kernel;
44using namespace Geometry;
45using namespace DataObjects;
48 auto wsValidator = std::make_shared<CompositeValidator>(CompositeRelation::OR);
49 auto monoValidator = std::make_shared<CompositeValidator>(CompositeRelation::AND);
50 auto monoUnitValidator = std::make_shared<CompositeValidator>(CompositeRelation::OR);
51 auto tofValidator = std::make_shared<CompositeValidator>(CompositeRelation::AND);
56 monoValidator->add(monoUnitValidator);
64 wsValidator->add(monoValidator);
65 wsValidator->add(tofValidator);
68 "Input workspace containing the SANS 2D data");
70 "Workspace that will contain the I(Q) data");
72 "The new bin boundaries in the form: :math:`x_1,\\Delta x_1,x_2,\\Delta "
75 auto positiveInt = std::make_shared<BoundedValidator<int>>();
76 positiveInt->setLower(0);
77 auto positiveDouble = std::make_shared<BoundedValidator<double>>();
78 positiveDouble->setLower(0);
81 "Number of sub-pixels used for each detector pixel in each "
82 "direction.The total number of sub-pixels will be "
83 "NPixelDivision*NPixelDivision.");
86 declareProperty(
"NumberOfWedges", 2, positiveInt,
"Number of wedges to calculate.");
87 declareProperty(
"WedgeAngle", 30.0, positiveDouble,
"Opening angle of the wedge, in degrees.");
88 declareProperty(
"WedgeOffset", 0.0, positiveDouble,
"Wedge offset relative to the horizontal axis, in degrees.");
91 "Name for the WorkspaceGroup containing the wedge I(q) distributions.");
93 declareProperty(
"PixelSizeX", 5.15, positiveDouble,
"Pixel size in the X direction (mm).");
94 declareProperty(
"PixelSizeY", 5.15, positiveDouble,
"Pixel size in the Y direction (mm).");
95 declareProperty(
"ErrorWeighting",
false,
"Choose whether each pixel contribution will be weighted by 1/error^2.");
97 declareProperty(
"AsymmetricWedges",
false,
"Choose to produce the results for asymmetric wedges.");
99 declareProperty(
"AccountForGravity",
false,
"Take the nominal gravity drop into account.");
103 "Table workspace containing the shapes (sectors only) drawn in the "
104 "instrument viewer; if specified, the wedges properties defined above "
105 "are not taken into account.");
135 const std::vector<double> binParams =
getProperty(
"OutputBinning");
141 m_nSpec = inputWS->getNumberHistograms();
144 m_nBins = inputWS->readY(0).size();
155 const double wedgeOffset =
getProperty(
"WedgeOffset");
156 const double wedgeAngle =
getProperty(
"WedgeAngle");
159 for (
size_t iw = 0; iw <
m_nWedges; ++iw) {
163 double outerRadius = -1.;
166 double angleRange = wedgeAngle *
deg2rad;
167 double midAngle = M_PI *
static_cast<double>(iw) /
static_cast<double>(
m_nWedges);
170 midAngle += wedgeOffset *
deg2rad;
176 g_log.
warning(
"This option is still in active development and might be "
177 "subject to changes in the next version.");
188 m_intensities = std::vector<std::vector<std::vector<double>>>(
201 size_t rowCount = shapeWs->rowCount();
203 std::map<std::string, std::vector<double>> viewportParams;
208 for (
size_t i = 0; i < rowCount - 1; ++i) {
209 std::map<std::string, std::vector<std::string>> paramMap;
210 std::vector<std::string> splitParams;
211 boost::algorithm::split(splitParams, shapeWs->String(i, 1), boost::algorithm::is_any_of(
"\n"));
212 std::vector<std::string> params;
213 for (std::string val : splitParams) {
216 boost::algorithm::split(params, val, boost::algorithm::is_any_of(
"\t"));
219 paramMap[params[0]] = params;
221 if (paramMap[
"Type"][1] ==
"sector")
224 g_log.
information() <<
"Shape " << i + 1 <<
" is of type " << paramMap[
"Type"][1]
225 <<
" which is not supported. This shape is ignored." << std::endl;
230 return wedgeA.angleMiddle < wedgeB.angleMiddle;
244 std::map<std::string, std::vector<double>> &viewportParams) {
246 std::vector<std::string> params;
247 boost::algorithm::split(params, viewport, boost::algorithm::is_any_of(
"\t, \n"));
249 if (params[0] !=
"Translation") {
251 "No viewport found in the shape table. Please provide a table using shapes drawn in the Full3D projection.");
255 viewportParams[params[0]] = std::vector<double>(2);
256 viewportParams[params[0]][0] = std::stod(params[1]);
257 viewportParams[params[0]][1] = std::stod(params[2]);
260 viewportParams[params[3]] = std::vector<double>(1);
261 viewportParams[params[3]][0] = std::stod(params[4]);
264 viewportParams[params[5]] = std::vector<double>(4);
265 viewportParams[params[5]][0] = std::stod(params[6]);
266 viewportParams[params[5]][1] = std::stod(params[7]);
267 viewportParams[params[5]][2] = std::stod(params[8]);
268 viewportParams[params[5]][3] = std::stod(params[9]);
270 double epsilon = 1e-10;
272 if (std::fabs(viewportParams[
"Rotation"][0]) > epsilon || std::fabs(viewportParams[
"Rotation"][1]) > epsilon ||
273 std::fabs(viewportParams[
"Rotation"][3]) > epsilon || std::fabs(viewportParams[
"Rotation"][2] - 1) > epsilon) {
274 g_log.
warning(
"The shapes were created using a rotated viewport not using Z- projection, which is not supported. "
275 "Results are likely to be erroneous. Consider freezing the rotation in the instrument viewer.");
285 const std::map<std::string, std::vector<double>> &viewport) {
286 double zoom = viewport.at(
"Zoom")[0];
289 double outerRadius = std::stod(params[2]) / zoom;
291 double startAngle = std::stod(params[3]);
292 double endAngle = std::stod(params[4]);
294 double centerAngle = (startAngle + endAngle) / 2;
295 if (endAngle < startAngle)
296 centerAngle = std::fmod(centerAngle + M_PI, 2 * M_PI);
298 double angleRange = std::fmod(endAngle - startAngle, 2 * M_PI);
299 angleRange = angleRange >= 0 ? angleRange : angleRange + 2 * M_PI;
302 centerAngle = std::fmod(3 * M_PI - centerAngle, 2 * M_PI);
304 double xOffset = viewport.at(
"Translation")[0];
305 double yOffset = viewport.at(
"Translation")[1];
307 double centerX = -(std::stod(params[5]) - xOffset) / zoom;
308 double centerY = (std::stod(params[6]) - yOffset) / zoom;
326 [&wedge](
const auto ¶ms) { return wedge.isSymmetric(params); });
339 <<
" degrees." << std::endl;
354 const auto &spectrumInfo = inputWS->spectrumInfo();
355 const V3D sourcePos = spectrumInfo.sourcePosition();
356 const V3D samplePos = spectrumInfo.samplePosition();
358 const V3D beamLine = samplePos - sourcePos;
360 const auto up = inputWS->getInstrument()->getReferenceFrame()->vecPointingUp();
361 double monoWavelength = 0;
363 if (inputWS->run().hasProperty(
"wavelength")) {
364 monoWavelength = inputWS->run().getPropertyAsSingleValue(
"wavelength");
366 throw std::runtime_error(
"Could not find wavelength in the sample logs.");
373 const auto i =
static_cast<size_t>(
index);
375 if (!spectrumInfo.hasDetectors(i) || spectrumInfo.isMonitor(i) || spectrumInfo.isMasked(i)) {
380 std::vector<size_t> maskedBins;
382 if (inputWS->hasMaskedBins(i)) {
383 maskedBins = inputWS->maskedBinsIndices(i);
387 const auto &XIn = inputWS->x(i);
388 const auto &YIn = inputWS->y(i);
389 const auto &EIn = inputWS->e(i);
392 const V3D pos = spectrumInfo.position(i) - samplePos;
399 for (
size_t j = 0; j <
m_nBins; ++j) {
402 if (std::binary_search(maskedBins.cbegin(), maskedBins.cend(), j)) {
406 const double wavelength =
m_isMonochromatic ? monoWavelength : (XIn[j] + XIn[j + 1]) / 2.;
425 const double sinTheta = sin(0.5 *
position.angle(beamLine));
426 const double q = 4.0 * M_PI * sinTheta / wavelength;
455 m_errors[0][j][k] += w * w * EIn[j] * EIn[j];
459 for (
size_t iw = 0; iw <
m_nWedges; ++iw) {
465 double angle = std::fabs((subPix - center).angle(
V3D(cos(centerAngle), sin(centerAngle), 0.0)));
469 bool isWithinAngularRange =
475 if (isWithinAngularRange && isWithinRadii) {
479 m_errors[iw + 1][j][k] += w * w * EIn[j] * EIn[j];
503 auto wsgroup = std::make_shared<WorkspaceGroup>();
506 for (
size_t iw = 0; iw <
m_nWedges; ++iw) {
509 wsgroup->addWorkspace(wedgeWs);
513 if (outputWSGroupName.empty()) {
515 outputWSGroupName = outputWSName +
"_wedges";
521 for (
size_t iout = 0; iout <
m_nWedges + 1; ++iout) {
523 auto ws = (iout == 0) ? outputWS : std::dynamic_pointer_cast<MatrixWorkspace>(wsgroup->getItem(iout - 1));
542 for (
size_t iSample = 0; iSample <
m_nBins; ++iSample) {
543 auto &YOut = outputWS->mutableY(iSample);
544 auto &EOut = outputWS->mutableE(iSample);
547 for (
int iq = 0; iq < static_cast<int>(
m_nQ); ++iq) {
552 EOut[iq] =
m_errors[iout][iSample][iq] / (norm * norm);
567 auto &YOut = outputWS->mutableY(0);
568 auto &EOut = outputWS->mutableE(0);
570 std::vector<double> normLambda(
m_nQ, 0.0);
572 for (
size_t il = 0; il <
m_nBins; ++il) {
574 for (
int iq = 0; iq < static_cast<int>(
m_nQ); ++iq) {
579 EOut[iq] +=
m_errors[iout][il][iq] / (norm * norm);
580 normLambda[iq] += 1.;
587 for (
size_t i = 0; i <
m_nQ; ++i) {
588 YOut[i] /= normLambda[i];
589 EOut[i] = sqrt(EOut[i]) / normLambda[i];
602 const std::vector<double> &binEdges,
const size_t nSpectra) {
606 for (
size_t iSpectra = 0; iSpectra <
nSpectra; ++iSpectra) {
607 outputWS->setBinEdges(iSpectra, binEdges);
609 outputWS->setYUnitLabel(
"1/cm");
610 outputWS->setDistribution(
true);
#define DECLARE_ALGORITHM(classname)
std::map< DeltaEMode::Type, std::string > index
#define PARALLEL_START_INTERRUPT_REGION
Begins a block to skip processing is the algorithm has been interupted Note the end of the block if n...
#define PARALLEL_CRITICAL(name)
#define PARALLEL_END_INTERRUPT_REGION
Ends a block to skip processing is the algorithm has been interupted Note the start of the block if n...
#define PARALLEL_FOR_IF(condition)
Empty definitions - to enable set your complier to enable openMP.
#define PARALLEL_CHECK_INTERRUPT_REGION
Adds a check after a Parallel region to see if it was interupted.
void declareProperty(std::unique_ptr< Kernel::Property > p, const std::string &doc="") override
Add a property to the list of managed properties.
std::string getPropertyValue(const std::string &name) const override
Get the value of a property as a string.
TypedValue getProperty(const std::string &name) const override
Get the value of a property.
void progress(double p, const std::string &msg="", double estimatedTime=0.0, int progressPrecision=0)
Sends ProgressNotification.
bool isDefault(const std::string &name) const
void setPropertyValue(const std::string &name, const std::string &value) override
Set the value of a property by string N.B.
A validator which checks that a workspace contains histogram data (the default) or point data as requ...
A validator which checks that a workspace has a valid instrument.
Helper class for reporting progress from algorithms.
A property class for workspaces.
A validator which checks that the unit of the workspace referred to by a WorkspaceProperty is the exp...
A helper class for calculating neutron's gravitional drop.
double gravitationalDrop(const double wav) const
void getTableShapes()
Q1DWeighted::getTableShapes if the user provided a shape table, parse the stored values and get the v...
void bootstrap(const API::MatrixWorkspace_const_sptr &)
Q1DWeighted::bootstrap initializes the user inputs.
void fillTOFOutput(API::MatrixWorkspace_sptr &, const size_t)
Q1DWeighted::fillTOFOutput Fill in the output workspace for TOF input, averaging over lambda bins.
std::vector< Wedge > m_wedgesParameters
std::vector< std::vector< std::vector< double > > > m_intensities
bool checkIfSymetricalWedge(Wedge &Wedge)
Q1DWeighted::checkIfSymetricalWedge Check if the symetrical wedge to the one defined by the parameter...
void fillMonochromaticOutput(API::MatrixWorkspace_sptr &, const size_t)
Q1DWeighted::fillMonochromaticOutput Fill the output workspace for monochromatic, kinetic input.
std::vector< std::vector< std::vector< double > > > m_errors
void finalize(const API::MatrixWorkspace_const_sptr &)
Q1DWeighted::finalize performs final averaging and sets the output workspaces.
void init() override
Initialisation code.
void getViewportParams(const std::string &, std::map< std::string, std::vector< double > > &)
Q1DWeighted::getViewportParams get the parameters defining the viewport of the instrument view when t...
std::vector< double > m_qBinEdges
void getWedgeParams(const std::vector< std::string > &, const std::map< std::string, std::vector< double > > &)
Q1DWeighted::getWedgeParams.
std::vector< std::vector< std::vector< double > > > m_normalisation
API::MatrixWorkspace_sptr createOutputWorkspace(const API::MatrixWorkspace_const_sptr &, const size_t, const std::vector< double > &, const size_t)
Create an output workspace.
void checkIfSuperposedWedges()
Q1DWeighted::checkIfSuperposedWedges Check if some wedges ahev the same angleMiddle,...
void calculate(const API::MatrixWorkspace_const_sptr &)
Q1DWeighted::calculate Performs the azimuthal averaging for each wavelength bin.
void exec() override
Execution code.
Support for a property that holds an array of values.
IPropertyManager * setProperty(const std::string &name, const T &value)
Templated method to set the value of a PropertyWithValue.
void error(const std::string &msg)
Logs at error level.
void warning(const std::string &msg)
Logs at warning level.
void information(const std::string &msg)
Logs at information level.
static T & Instance()
Return a reference to the Singleton instance, creating it if it does not already exist Creation is do...
double distance(const V3D &v) const noexcept
Calculates the distance between two vectors.
std::shared_ptr< ITableWorkspace > ITableWorkspace_sptr
shared pointer to Mantid::API::ITableWorkspace
std::shared_ptr< const MatrixWorkspace > MatrixWorkspace_const_sptr
shared pointer to the matrix workspace base class (const version)
std::shared_ptr< MatrixWorkspace > MatrixWorkspace_sptr
shared pointer to the matrix workspace base class
constexpr double deg2rad
Defines units/enum for Crystal work.
size_t MANTID_KERNEL_DLL indexOfValueFromEdges(const std::vector< double > &bin_edges, const double value)
Gets the bin of a value from a vector of bin edges.
int MANTID_KERNEL_DLL createAxisFromRebinParams(const std::vector< double > ¶ms, std::vector< double > &xnew, const bool resize_xnew=true, const bool full_bins_only=false, const double xMinHint=std::nan(""), const double xMaxHint=std::nan(""), const bool useReverseLogarithmic=false, const double power=-1)
Creates a new output X array given a 'standard' set of rebinning parameters.
std::enable_if< std::is_pointer< Arg >::value, bool >::type threadSafe(Arg workspace)
Thread-safety check Checks the workspace to ensure it is suitable for multithreaded access.
@ Input
An input workspace.
@ Output
An output workspace.