12#include "MantidHistogramData/HistogramDx.h"
13#include "MantidHistogramData/HistogramE.h"
14#include "MantidHistogramData/HistogramX.h"
15#include "MantidHistogramData/HistogramY.h"
23#include <boost/format.hpp>
24#include <boost/math/special_functions.hpp>
25#include <boost/tuple/tuple.hpp>
30using Mantid::HistogramData::HistogramE;
31using Mantid::HistogramData::HistogramY;
36using MinMaxTuple = boost::tuple<double, double>;
38 return MinMaxTuple(rhsWS->x(0).front(), lhsWS->x(0).back());
42bool isNonzero(
double i) {
return (0 != i); }
63 const auto histogramCount =
static_cast<int>(source->getNumberHistograms());
65 for (
int i = 0; i < histogramCount; ++i) {
68 product->setSharedX(i, source->sharedX(i));
70 const auto &sourceY = source->y(i);
71 const auto &sourceE = source->e(i);
74 product->mutableY(i) = HistogramY(sourceY.size(), 0);
75 product->mutableE(i) = HistogramE(sourceE.size(), 0);
77 auto &newY = product->mutableY(i);
78 auto &newE = product->mutableE(i);
81 std::copy(sourceY.begin() + a1 + 1, sourceY.begin() + a2, newY.begin() + a1 + 1);
82 std::copy(sourceE.begin() + a1 + 1, sourceE.begin() + a2, newE.begin() + a1 + 1);
98 const auto histogramCount =
static_cast<int>(source->getNumberHistograms());
100 for (
int i = 0; i < histogramCount; ++i) {
103 auto &sourceY = source->mutableY(i);
104 auto &sourceE = source->mutableE(i);
105 for (
int binIndex = a1; binIndex < a2; ++binIndex) {
106 sourceY[binIndex] = 0;
107 sourceE[binIndex] = 0;
121 "LHS input workspace.");
123 "RHS input workspace, must be same type as LHSWorkspace "
124 "(histogram or point data).");
126 "Output stitched workspace.");
128 "Start overlap x-value in units of x-axis.");
130 "End overlap x-value in units of x-axis.");
132 "Rebinning Parameters. See Rebin for format. If only a "
133 "single value is provided, start and end are taken from "
134 "input workspaces.");
136 "Scaling either with respect to LHS workspace or RHS workspace");
138 "True to use a provided value for the scale factor.");
139 auto manualScaleFactorValidator = std::make_shared<BoundedValidator<double>>();
140 manualScaleFactorValidator->setLower(0);
141 manualScaleFactorValidator->setExclusive(
true);
144 "Provided value for the scale factor.");
146 "The actual used value for the scaling factor.");
153 std::map<std::string, std::string> issues;
157 issues[
"LHSWorkspace"] =
"The LHSWorkspace must be a MatrixWorkspace.";
160 issues[
"RHSWorkspace"] =
"The RHSWorkspace must be a MatrixWorkspace.";
162 if (!issues.empty()) {
168 if (!compatible.empty()) {
170 if (!(compatible ==
"spectra must have either Dx values or not; ") ||
171 (
rhs->isHistogramData()))
172 issues[
"RHSWorkspace"] =
"Workspace " +
rhs->getName() +
" is not compatible: " + compatible +
"\n";
186 double startOverlapVal = this->
getProperty(
"StartOverlap");
188 const bool startOverlapBeyondRange = (startOverlapVal < intesectionMin) || (startOverlapVal > intesectionMax);
189 if (startOverlapProp->
isDefault() || startOverlapBeyondRange) {
190 if (!startOverlapProp->
isDefault() && startOverlapBeyondRange) {
192 std::sprintf(message,
193 "StartOverlap is outside range at %0.4f, Min is "
194 "%0.4f, Max is %0.4f . Forced to be: %0.4f",
195 startOverlapVal, intesectionMin, intesectionMax, intesectionMin);
198 startOverlapVal = intesectionMin;
199 std::stringstream buffer;
200 buffer <<
"StartOverlap calculated to be: " << startOverlapVal;
203 return startOverlapVal;
215 double endOverlapVal = this->
getProperty(
"EndOverlap");
217 const bool endOverlapBeyondRange = (endOverlapVal < intesectionMin) || (endOverlapVal > intesectionMax);
218 if (endOverlapProp->
isDefault() || endOverlapBeyondRange) {
219 if (!endOverlapProp->
isDefault() && endOverlapBeyondRange) {
221 std::sprintf(message,
222 "EndOverlap is outside range at %0.4f, Min is "
223 "%0.4f, Max is %0.4f . Forced to be: %0.4f",
224 endOverlapVal, intesectionMin, intesectionMax, intesectionMax);
227 endOverlapVal = intesectionMax;
228 std::stringstream buffer;
229 buffer <<
"EndOverlap calculated to be: " << endOverlapVal;
232 return endOverlapVal;
242 const bool scaleRHS)
const {
243 std::vector<double> inputParams = this->
getProperty(
"Params");
245 const bool areParamsDefault = prop->
isDefault();
247 const auto &lhsX = lhsWS->x(0);
248 auto it = std::min_element(lhsX.begin(), lhsX.end());
249 const double minLHSX = *it;
251 const auto &rhsX = rhsWS->x(0);
252 it = std::max_element(rhsX.begin(), rhsX.end());
253 const double maxRHSX = *it;
255 std::vector<double> result;
256 if (areParamsDefault) {
257 std::vector<double> calculatedParams;
261 double calculatedStep = 0;
265 calculatedStep = lhsX[1] - lhsX[0];
269 calculatedStep = rhsX[1] - rhsX[0];
272 calculatedParams.emplace_back(minLHSX);
273 calculatedParams.emplace_back(calculatedStep);
274 calculatedParams.emplace_back(maxRHSX);
275 result = calculatedParams;
277 if (inputParams.size() == 1) {
278 std::vector<double> calculatedParams;
279 calculatedParams.emplace_back(minLHSX);
280 calculatedParams.emplace_back(inputParams.front());
281 calculatedParams.emplace_back(maxRHSX);
282 result = calculatedParams;
284 result = inputParams;
297 rebin->setProperty(
"InputWorkspace", input);
298 rebin->setProperty(
"Params", params);
299 std::stringstream ssParams;
300 ssParams << params[0] <<
"," << params[1] <<
"," << params[2];
305 const auto histogramCount =
static_cast<int>(outWS->getNumberHistograms());
310 for (
int i = 0; i < histogramCount; ++i) {
317 auto &sourceY = outWS->mutableY(i);
318 auto &sourceE = outWS->mutableE(i);
320 for (
size_t j = 0; j < sourceY.size(); ++j) {
321 const double value = sourceY[j];
322 if (std::isnan(
value)) {
323 nanYIndexes.emplace_back(j);
325 }
else if (std::isinf(
value)) {
326 infYIndexes.emplace_back(j);
330 const double eValue = sourceE[j];
331 if (std::isnan(eValue)) {
332 nanEIndexes.emplace_back(j);
334 }
else if (std::isinf(eValue)) {
335 infEIndexes.emplace_back(j);
359 g_log.
information(
"Integration RangeLower: " + boost::lexical_cast<std::string>(start));
360 g_log.
information(
"Integration RangeUpper: " + boost::lexical_cast<std::string>(stop));
362 return integration->getProperty(
"OutputWorkspace");
385 const std::string in1 =
"__Stitch1D_intermediate_workspace_1__";
386 const std::string in2 =
"__Stitch1D_intermediate_workspace_2__";
390 conjoinX->initialize();
391 conjoinX->setProperty(
"InputWorkspaces", std::vector<std::string>{in1, in2});
396 return std::dynamic_pointer_cast<Mantid::API::MatrixWorkspace>(ws);
406 auto a1 =
static_cast<int>(
workspace->yIndexOfX(startOverlap));
407 auto a2 =
static_cast<int>(
workspace->yIndexOfX(endOverlap));
409 throw std::runtime_error(
"The Params you have provided for binning yield a "
410 "workspace in which start and end overlap appear "
411 "in the same bin. Make binning finer via input "
414 return boost::tuple<int, int>(a1, a2);
422 auto ws_size =
static_cast<int64_t
>(ws->getNumberHistograms());
423 bool hasNonZeroErrors =
false;
425 for (int64_t i = 0; i < ws_size; ++i) {
427 if (!hasNonZeroErrors)
429 const auto &e = ws->e(i);
430 auto it = std::find_if(e.begin(), e.end(), isNonzero);
433 hasNonZeroErrors =
true;
441 return hasNonZeroErrors;
454 if (ws->size() == dxWS->size()) {
455 for (
size_t i = 0; i < ws->getNumberHistograms(); ++i) {
456 if (dxWS->hasDx(i) && !ws->hasDx(i) && !ws->isHistogramData()) {
457 ws->setSharedDx(i, dxWS->sharedDx(i));
463 if (m_scaleFactor < 1e-2 || m_scaleFactor > 1e2 || std::isnan(
m_scaleFactor)) {
464 std::stringstream messageBuffer;
465 messageBuffer <<
"Stitch1D calculated scale factor is: " <<
m_scaleFactor
466 <<
". Check the overlap region is non-zero.";
477 const MinMaxTuple intesectionXRegion = calculateXIntersection(lhsWS, rhsWS);
479 const size_t histogramCount = rhsWS->getNumberHistograms();
485 const double intersectionMin = intesectionXRegion.get<0>();
486 const double intersectionMax = intesectionXRegion.get<1>();
487 double startOverlap =
getStartOverlap(intersectionMin, intersectionMax);
488 double endOverlap =
getEndOverlap(intersectionMin, intersectionMax);
489 if (startOverlap > endOverlap) {
490 std::string message = boost::str(boost::format(
"Stitch1D cannot have a StartOverlap > EndOverlap. "
491 "StartOverlap: %0.9f, EndOverlap: %0.9f") %
492 startOverlap % endOverlap);
493 throw std::runtime_error(message);
496 const bool scaleRHS = this->
getProperty(
"ScaleRHSWorkspace");
500 if (lhsWS->isHistogramData()) {
502 const double xMin = params.front();
503 const double xMax = params.back();
505 if (std::abs(xMin - startOverlap) < 1E-6)
508 if (std::abs(xMax - endOverlap) < 1E-6)
511 if (startOverlap < xMin) {
512 std::string message = boost::str(boost::format(
"Stitch1D StartOverlap is outside the available X range. "
513 "StartOverlap: %10.9f, X min: %10.9f") %
514 startOverlap % xMin);
515 throw std::runtime_error(message);
517 if (endOverlap > xMax) {
518 std::string message = boost::str(boost::format(
"Stitch1D EndOverlap is outside the available X range. "
519 "EndOverlap: %10.9f, X max: %10.9f") %
521 throw std::runtime_error(message);
523 lhs =
rebin(lhs, params);
529 const bool useManualScaleFactor = this->
getProperty(
"UseManualScaleFactor");
530 if (useManualScaleFactor) {
536 auto rhsOverlapIntegrated =
integration(
rhs, startOverlap, endOverlap);
537 auto lhsOverlapIntegrated =
integration(lhs, startOverlap, endOverlap);
539 auto scalingFactors = lhsOverlapIntegrated / rhsOverlapIntegrated;
542 auto scalingFactors = rhsOverlapIntegrated / lhsOverlapIntegrated;
547 std::stringstream messageBuffer;
552 if (lhsWS->isHistogramData()) {
554 int a1 = boost::tuples::get<0>(startEnd);
555 int a2 = boost::tuples::get<1>(startEnd);
561 maskInPlace(a1 + 1,
static_cast<int>(lhs->blocksize()), lhs);
570 overlapave = sum * 0.5;
572 result = lhs + overlapave +
rhs;
577 g_log.
error(
"Could not retrieve joined workspace.");
581 childAlg->setProperty(
"InputWorkspace", result);
583 result = childAlg->getProperty(
"OutputWorkspace");
593 auto histogramCount =
static_cast<int>(ws->getNumberHistograms());
595 for (
int i = 0; i < histogramCount; ++i) {
598 auto &sourceY = ws->mutableY(i);
601 sourceY[j] = std::numeric_limits<double>::quiet_NaN();
605 sourceY[j] = std::numeric_limits<double>::infinity();
609 sourceY[j] = std::numeric_limits<double>::quiet_NaN();
613 sourceY[j] = std::numeric_limits<double>::infinity();
#define DECLARE_ALGORITHM(classname)
const std::vector< double > & rhs
double value
The value of the point.
IPeaksWorkspace_sptr workspace
#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.
TypedValue getProperty(const std::string &name) const override
Get the value of a property.
virtual std::shared_ptr< Algorithm > createChildAlgorithm(const std::string &name, const double startProgress=-1., const double endProgress=-1., const bool enableLogging=true, const int &version=-1)
Create a Child Algorithm.
A property class for workspaces.
void setReferenceProperties(const API::MatrixWorkspace_sptr &)
Sets the properties of the reference (usually first) workspace, to later check the compatibility of t...
std::string checkCompatibility(const API::MatrixWorkspace_sptr &, bool checkNumberHistograms=false)
Compares the properties of the input workspace with the reference.
Stitch1D : Stitches two Matrix Workspaces together into a single output.
SpecialTypeIndexes m_infEIndexes
Index per workspace spectra of Infs.
Mantid::API::MatrixWorkspace_sptr weightedMean(Mantid::API::MatrixWorkspace_sptr &inOne, Mantid::API::MatrixWorkspace_sptr &inTwo)
Calculate the weighted mean.
Mantid::API::MatrixWorkspace_sptr maskAllBut(int a1, int a2, Mantid::API::MatrixWorkspace_sptr &source)
Mask out everything but the data in the ranges.
SpecialTypeIndexes m_nanYIndexes
bool hasNonzeroErrors(Mantid::API::MatrixWorkspace_sptr &ws)
Does the x-axis have non-zero errors.
boost::tuple< int, int > findStartEndIndexes(double startOverlap, double endOverlap, Mantid::API::MatrixWorkspace_sptr &workspace)
Find the start and end indexes.
Mantid::API::MatrixWorkspace_sptr rebin(Mantid::API::MatrixWorkspace_sptr &input, const std::vector< double > ¶ms)
Perform rebin.
double m_errorScaleFactor
std::vector< double > getRebinParams(Mantid::API::MatrixWorkspace_const_sptr &lhsWS, Mantid::API::MatrixWorkspace_const_sptr &rhsWS, const bool scaleRHS) const
Get the rebin parameters.
Mantid::API::MatrixWorkspace_sptr integration(Mantid::API::MatrixWorkspace_sptr &input, const double start, const double stop)
Perform integration.
void maskInPlace(int a1, int a2, Mantid::API::MatrixWorkspace_sptr &source)
Mask out everything but the data in the ranges, but do it inplace.
Mantid::API::MatrixWorkspace_sptr conjoinXAxis(Mantid::API::MatrixWorkspace_sptr &inOne, Mantid::API::MatrixWorkspace_sptr &inTwo)
Conjoin x axis.
void scaleWorkspace(API::MatrixWorkspace_sptr &ws, API::MatrixWorkspace_sptr &scaleFactorWS, API::MatrixWorkspace_const_sptr &dxWS)
Scale workspace (left hand side or right hand side)
double getStartOverlap(const double intesectionMin, const double intesectionMax) const
Get the start overlap.
SpecialTypeIndexes m_nanEIndexes
Index per workspace spectra of Nans.
void exec() override
Overwrites Algorithm method.
SpecialTypeIndexes m_infYIndexes
double m_scaleFactor
Scaling factors.
void init() override
Overwrites Algorithm method.
double getEndOverlap(const double intesectionMin, const double intesectionMax) const
Get the end overlap.
static const double range_tolerance
Range tolerance.
void reinsertSpecialValues(const Mantid::API::MatrixWorkspace_sptr &ws)
Add back in any special values.
std::map< std::string, std::string > validateInputs() override
Cross-check properties with each other.
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 notice(const std::string &msg)
Logs at notice level.
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.
The concrete, templated class for properties.
Base class for properties.
virtual bool isDefault() const =0
Overriden function that returns if property has the same value that it was initialised with,...
static T & Instance()
Return a reference to the Singleton instance, creating it if it does not already exist Creation is do...
std::shared_ptr< Workspace > Workspace_sptr
shared pointer to Mantid::API::Workspace
std::shared_ptr< const MatrixWorkspace > MatrixWorkspace_const_sptr
shared pointer to the matrix workspace base class (const version)
std::shared_ptr< Algorithm > Algorithm_sptr
Typedef for a shared pointer to an Algorithm.
std::shared_ptr< MatrixWorkspace > MatrixWorkspace_sptr
shared pointer to the matrix workspace base class
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.
std::vector< double > MantidVec
typedef for the data storage used in Mantid matrix workspaces
constexpr double EMPTY_DBL() noexcept
Returns what we consider an "empty" double within a property.
@ Input
An input workspace.
@ Output
An output workspace.