diff --git a/Code/Mantid/Framework/Algorithms/inc/MantidAlgorithms/Stitch1D.h b/Code/Mantid/Framework/Algorithms/inc/MantidAlgorithms/Stitch1D.h
index 26366341f99cf259b0d2119440f3ac7730ba935b..aa6017f87cc9234fb4fe5cf37dbf431a26ae3196 100644
--- a/Code/Mantid/Framework/Algorithms/inc/MantidAlgorithms/Stitch1D.h
+++ b/Code/Mantid/Framework/Algorithms/inc/MantidAlgorithms/Stitch1D.h
@@ -12,41 +12,59 @@ namespace Mantid
/** Stitch1D : Stitches two Matrix Workspaces together into a single output.
- Copyright © 2014 ISIS Rutherford Appleton Laboratory & NScD Oak Ridge National Laboratory
+ Copyright © 2014 ISIS Rutherford Appleton Laboratory & NScD Oak Ridge National Laboratory
- This file is part of Mantid.
+ This file is part of Mantid.
- Mantid is free software; you can redistribute it and/or modify
- it under the terms of the GNU General Public License as published by
- the Free Software Foundation; either version 3 of the License, or
- (at your option) any later version.
+ Mantid is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 3 of the License, or
+ (at your option) any later version.
- Mantid is distributed in the hope that it will be useful,
- but WITHOUT ANY WARRANTY; without even the implied warranty of
- MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- GNU General Public License for more details.
+ Mantid is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
- You should have received a copy of the GNU General Public License
- along with this program. If not, see <http://www.gnu.org/licenses/>.
+ You should have received a copy of the GNU General Public License
+ along with this program. If not, see <http://www.gnu.org/licenses/>.
- File change history is stored at: <https://github.com/mantidproject/mantid>
- Code Documentation is available at: <http://doxygen.mantidproject.org>
- */
- class DLLExport Stitch1D : public API::Algorithm
+ File change history is stored at: <https://github.com/mantidproject/mantid>
+ Code Documentation is available at: <http://doxygen.mantidproject.org>
+ */
+ class DLLExport Stitch1D: public API::Algorithm
{
public:
/// Default constructor
- Stitch1D(){};
+ Stitch1D()
+ {
+ }
+ ;
/// Destructor
- virtual ~Stitch1D(){};
+ virtual ~Stitch1D()
+ {
+ }
+ ;
/// Algorithm's name for identification. @see Algorithm::name
- virtual const std::string name() const {return "Stitch1D";}
+ virtual const std::string name() const
+ {
+ return "Stitch1D";
+ }
/// Algorithm's version for identification. @see Algorithm::version
- virtual int version() const {return 3;}
+ virtual int version() const
+ {
+ return 3;
+ }
/// Algorithm's category for identification. @see Algorithm::category
- virtual const std::string category() const {return "Reflectometry";}
+ virtual const std::string category() const
+ {
+ return "Reflectometry";
+ }
///Summary of algorithm's purpose
- virtual const std::string summary() const {return "Stitches single histogram matrix workspaces together";}
+ virtual const std::string summary() const
+ {
+ return "Stitches single histogram matrix workspaces together";
+ }
/// Does the x-axis have non-zero errors
bool hasNonzeroErrors(Mantid::API::MatrixWorkspace_sptr ws);
private:
@@ -60,27 +78,38 @@ namespace Mantid
double getEndOverlap(const double& intesectionMin, const double& intesectionMax) const;
/// Get the rebin parameters
- Mantid::MantidVec getRebinParams(Mantid::API::MatrixWorkspace_sptr& lhsWS, Mantid::API::MatrixWorkspace_sptr& rhsWS, const bool scaleRHSWS) const;
+ Mantid::MantidVec getRebinParams(Mantid::API::MatrixWorkspace_sptr& lhsWS,
+ Mantid::API::MatrixWorkspace_sptr& rhsWS, const bool scaleRHSWS) const;
/// Perform rebin
- Mantid::API::MatrixWorkspace_sptr rebin(Mantid::API::MatrixWorkspace_sptr& input, const Mantid::MantidVec& params);
+ Mantid::API::MatrixWorkspace_sptr rebin(Mantid::API::MatrixWorkspace_sptr& input,
+ const Mantid::MantidVec& params);
/// Perform integration
- Mantid::API::MatrixWorkspace_sptr specialIntegration(Mantid::API::MatrixWorkspace_sptr& input, const double& start, const double& stop);
+ Mantid::API::MatrixWorkspace_sptr specialIntegration(Mantid::API::MatrixWorkspace_sptr& input,
+ const double& start, const double& stop);
/// Perform multiplication over a range
- Mantid::API::MatrixWorkspace_sptr multiplyRange(Mantid::API::MatrixWorkspace_sptr& input, const int& startBin, const int& endBin, const double& factor);
+ Mantid::API::MatrixWorkspace_sptr multiplyRange(Mantid::API::MatrixWorkspace_sptr& input,
+ const int& startBin, const int& endBin, const double& factor);
/// Perform multiplication over a range
- Mantid::API::MatrixWorkspace_sptr multiplyRange(Mantid::API::MatrixWorkspace_sptr& input, const int& startBin, const double& factor);
+ Mantid::API::MatrixWorkspace_sptr multiplyRange(Mantid::API::MatrixWorkspace_sptr& input,
+ const int& startBin, const double& factor);
/// Create a single valued workspace
Mantid::API::MatrixWorkspace_sptr singleValueWS(double val);
/// Calclate the weighted mean
- Mantid::API::MatrixWorkspace_sptr weightedMean(Mantid::API::MatrixWorkspace_sptr& inOne, Mantid::API::MatrixWorkspace_sptr& inTwo);
+ Mantid::API::MatrixWorkspace_sptr weightedMean(Mantid::API::MatrixWorkspace_sptr& inOne,
+ Mantid::API::MatrixWorkspace_sptr& inTwo);
/// Find the start and end indexes
- boost::tuple<int,int> findStartEndIndexes(double startOverlap, double endOverlap, Mantid::API::MatrixWorkspace_sptr& workspace);
+ boost::tuple<int, int> findStartEndIndexes(double startOverlap, double endOverlap,
+ Mantid::API::MatrixWorkspace_sptr& workspace);
+ /// Mask out everything but the data in the ranges
+ Mantid::API::MatrixWorkspace_sptr maskAllBut(int a1, int a2,
+ Mantid::API::MatrixWorkspace_sptr & source);
+ /// Mask out everything but the data in the ranges, but do it inplace.
+ void maskInPlace(int a1, int a2, Mantid::API::MatrixWorkspace_sptr source);
/// Range tolerance
static const double range_tolerance;
};
-
} // namespace Algorithms
} // namespace Mantid
diff --git a/Code/Mantid/Framework/Algorithms/src/Stitch1D.cpp b/Code/Mantid/Framework/Algorithms/src/Stitch1D.cpp
index 6cb8e14c76ded914c8da9729d518d101d860a7fd..ddd32fdba88afacef752871599d212293318d9e5 100644
--- a/Code/Mantid/Framework/Algorithms/src/Stitch1D.cpp
+++ b/Code/Mantid/Framework/Algorithms/src/Stitch1D.cpp
@@ -16,7 +16,6 @@
#include <vector>
#include <algorithm>
-
using namespace Mantid::Kernel;
using namespace Mantid::API;
using Mantid::MantidVec;
@@ -50,473 +49,526 @@ namespace Mantid
* inclusive of bin boundaries if they are sitting ontop of the bin boundaries.
*/
const double Stitch1D::range_tolerance = 1e-9;
- // Register the algorithm into the AlgorithmFactory
+// Register the algorithm into the AlgorithmFactory
DECLARE_ALGORITHM(Stitch1D)
- //----------------------------------------------------------------------------------------------
- /** Initialize the algorithm's properties.
- */
- void Stitch1D::init()
- {
- Kernel::IValidator_sptr histogramValidator = boost::make_shared<HistogramValidator>();
-
- declareProperty(
- new WorkspaceProperty<MatrixWorkspace>("LHSWorkspace", "", Direction::Input,
- histogramValidator->clone()), "LHS input workspace.");
- declareProperty(
- new WorkspaceProperty<MatrixWorkspace>("RHSWorkspace", "", Direction::Input,
- histogramValidator->clone()), "RHS input workspace.");
- declareProperty(new WorkspaceProperty<MatrixWorkspace>("OutputWorkspace", "", Direction::Output),
- "Output stitched workspace.");
- declareProperty(
- new PropertyWithValue<double>("StartOverlap", Mantid::EMPTY_DBL(), Direction::Input),
- "Start overlap x-value in units of x-axis. Optional.");
- declareProperty(new PropertyWithValue<double>("EndOverlap", Mantid::EMPTY_DBL(), Direction::Input),
- "End overlap x-value in units of x-axis. Optional.");
- declareProperty(
- new ArrayProperty<double>("Params", boost::make_shared<RebinParamsValidator>(true)),
- "Rebinning Parameters. See Rebin for format. If only a single value is provided, start and end are taken from input workspaces.");
- declareProperty(new PropertyWithValue<bool>("ScaleRHSWorkspace", true, Direction::Input),
- "Scaling either with respect to workspace 1 or workspace 2");
- declareProperty(new PropertyWithValue<bool>("UseManualScaleFactor", false, Direction::Input),
- "True to use a provided value for the scale factor.");
- auto manualScaleFactorValidator = boost::make_shared<BoundedValidator<double> >();
- manualScaleFactorValidator->setLower(0);
- manualScaleFactorValidator->setExclusive(true);
- declareProperty(new PropertyWithValue<double>("ManualScaleFactor", 1.0, manualScaleFactorValidator, Direction::Input),
- "Provided value for the scale factor. Optional.");
- declareProperty(
- new PropertyWithValue<double>("OutScaleFactor", Mantid::EMPTY_DBL(), Direction::Output),
- "The actual used value for the scaling factor.");
- }
-
- /**Gets the start of the overlapping region
- @param intesectionMin :: The minimum possible value for the overlapping region to inhabit
- @param intesectionMax :: The maximum possible value for the overlapping region to inhabit
- @return a double contianing the start of the overlapping region
+ /**
+ * Zero out all y and e data that is not in the region a1 to a2.
+ * @param a1 : Zero based bin index (first one)
+ * @param a2 : Zero based bin index (last one inclusive)
+ * @param source : Workspace providing the source data.
+ * @return Masked workspace.
*/
- double Stitch1D::getStartOverlap(const double& intesectionMin, const double& intesectionMax) const
+ MatrixWorkspace_sptr Stitch1D::maskAllBut(int a1, int a2, MatrixWorkspace_sptr & source)
{
- Property* startOverlapProp = this->getProperty("StartOverlap");
- double startOverlapVal = this->getProperty("StartOverlap");
- startOverlapVal -= this->range_tolerance;
- const bool startOverlapBeyondRange = (startOverlapVal < intesectionMin)
- || (startOverlapVal > intesectionMax);
- if (startOverlapProp->isDefault() || startOverlapBeyondRange)
+ MatrixWorkspace_sptr product = WorkspaceFactory::Instance().create(source);
+ const int histogramCount = static_cast<int>(source->getNumberHistograms());
+ PARALLEL_FOR2(source,product)
+ for (int i = 0; i < histogramCount; ++i)
{
- if (!startOverlapProp->isDefault() && startOverlapBeyondRange)
- {
- char message[200];
- std::sprintf(message,
- "StartOverlap is outside range at %0.4f, Min is %0.4f, Max is %0.4f . Forced to be: %0.4f",
- startOverlapVal, intesectionMin, intesectionMax, intesectionMin);
- g_log.warning(std::string(message));
- }
- startOverlapVal = intesectionMin;
- std::stringstream buffer;
- buffer << "StartOverlap calculated to be: " << startOverlapVal;
- g_log.information(buffer.str());
- }
- return startOverlapVal;
+ PARALLEL_START_INTERUPT_REGION
+ // Copy over the bin boundaries
+ product->setX(i, source->refX(i));
+ // Copy over the data
+ const MantidVec& sourceY = source->readY(i);
+ const MantidVec& sourceE = source->readE(i);
+
+ // initially zero - out the data.
+ product->dataY(i) = MantidVec(sourceY.size(), 0);
+ product->dataE(i) = MantidVec(sourceE.size(), 0);
+
+ MantidVec& newY = product->dataY(i);
+ MantidVec& newE = product->dataE(i);
+
+ // Copy over the non-zero stuff
+ std::copy(sourceY.begin() + a1 + 1, sourceY.begin() + a2, newY.begin() + a1 + 1);
+ std::copy(sourceE.begin() + a1 + 1, sourceE.begin() + a2, newE.begin() + a1 + 1);
+
+ PARALLEL_END_INTERUPT_REGION
}
+ PARALLEL_CHECK_INTERUPT_REGION
+ return product;
+ }
- /**Gets the end of the overlapping region
- @param intesectionMin :: The minimum possible value for the overlapping region to inhabit
- @param intesectionMax :: The maximum possible value for the overlapping region to inhabit
- @return a double contianing the end of the overlapping region
- */
- double Stitch1D::getEndOverlap(const double& intesectionMin, const double& intesectionMax) const
+ void Stitch1D::maskInPlace(int a1, int a2, MatrixWorkspace_sptr source)
+ {
+ MatrixWorkspace_sptr product = WorkspaceFactory::Instance().create(source);
+ const int histogramCount = static_cast<int>(source->getNumberHistograms());
+ PARALLEL_FOR2(source,product)
+ for (int i = 0; i < histogramCount; ++i)
{
- Property* endOverlapProp = this->getProperty("EndOverlap");
- double endOverlapVal = this->getProperty("EndOverlap");
- endOverlapVal += this->range_tolerance;
- const bool endOverlapBeyondRange = (endOverlapVal < intesectionMin)
- || (endOverlapVal > intesectionMax);
- if (endOverlapProp->isDefault() || endOverlapBeyondRange)
+ PARALLEL_START_INTERUPT_REGION
+ // Copy over the data
+ MantidVec& sourceY = source->dataY(i);
+ MantidVec& sourceE = source->dataE(i);
+
+ for (int i = a1; i < a2; ++i)
{
- if (!endOverlapProp->isDefault() && endOverlapBeyondRange)
- {
- char message[200];
- std::sprintf(message,
- "EndOverlap is outside range at %0.4f, Min is %0.4f, Max is %0.4f . Forced to be: %0.4f",
- endOverlapVal, intesectionMin, intesectionMax, intesectionMax);
- g_log.warning(std::string(message));
- }
- endOverlapVal = intesectionMax;
- std::stringstream buffer;
- buffer << "EndOverlap calculated to be: " << endOverlapVal;
- g_log.information(buffer.str());
+ sourceY[i] = 0;
+ sourceE[i] = 0;
}
- return endOverlapVal;
- }
- /**Gets the rebinning parameters and calculates any missing values
- @param lhsWS :: The left hand side input workspace
- @param rhsWS :: The right hand side input workspace
- @param scaleRHS :: Scale the right hand side workspace
- @return a vector<double> contianing the rebinning parameters
- */
- MantidVec Stitch1D::getRebinParams(MatrixWorkspace_sptr& lhsWS, MatrixWorkspace_sptr& rhsWS, const bool scaleRHS) const
- {
- MantidVec inputParams = this->getProperty("Params");
- Property* prop = this->getProperty("Params");
- const bool areParamsDefault = prop->isDefault();
+ PARALLEL_END_INTERUPT_REGION
+ }
+PARALLEL_CHECK_INTERUPT_REGION
+}
- const MantidVec& lhsX = lhsWS->readX(0);
- auto it = std::min_element(lhsX.begin(), lhsX.end());
- const double minLHSX = *it;
+//----------------------------------------------------------------------------------------------
+/** Initialize the algorithm's properties.
+ */
+void Stitch1D::init()
+{
+Kernel::IValidator_sptr histogramValidator = boost::make_shared<HistogramValidator>();
+
+declareProperty(
+ new WorkspaceProperty<MatrixWorkspace>("LHSWorkspace", "", Direction::Input,
+ histogramValidator->clone()), "LHS input workspace.");
+declareProperty(
+ new WorkspaceProperty<MatrixWorkspace>("RHSWorkspace", "", Direction::Input,
+ histogramValidator->clone()), "RHS input workspace.");
+declareProperty(new WorkspaceProperty<MatrixWorkspace>("OutputWorkspace", "", Direction::Output),
+ "Output stitched workspace.");
+declareProperty(new PropertyWithValue<double>("StartOverlap", Mantid::EMPTY_DBL(), Direction::Input),
+ "Start overlap x-value in units of x-axis. Optional.");
+declareProperty(new PropertyWithValue<double>("EndOverlap", Mantid::EMPTY_DBL(), Direction::Input),
+ "End overlap x-value in units of x-axis. Optional.");
+declareProperty(new ArrayProperty<double>("Params", boost::make_shared<RebinParamsValidator>(true)),
+ "Rebinning Parameters. See Rebin for format. If only a single value is provided, start and end are taken from input workspaces.");
+declareProperty(new PropertyWithValue<bool>("ScaleRHSWorkspace", true, Direction::Input),
+ "Scaling either with respect to workspace 1 or workspace 2");
+declareProperty(new PropertyWithValue<bool>("UseManualScaleFactor", false, Direction::Input),
+ "True to use a provided value for the scale factor.");
+auto manualScaleFactorValidator = boost::make_shared<BoundedValidator<double> >();
+manualScaleFactorValidator->setLower(0);
+manualScaleFactorValidator->setExclusive(true);
+declareProperty(
+ new PropertyWithValue<double>("ManualScaleFactor", 1.0, manualScaleFactorValidator,
+ Direction::Input), "Provided value for the scale factor. Optional.");
+declareProperty(new PropertyWithValue<double>("OutScaleFactor", Mantid::EMPTY_DBL(), Direction::Output),
+ "The actual used value for the scaling factor.");
+}
- const MantidVec& rhsX = rhsWS->readX(0);
- it = std::max_element(rhsX.begin(), rhsX.end());
- const double maxRHSX = *it;
+/**Gets the start of the overlapping region
+ @param intesectionMin :: The minimum possible value for the overlapping region to inhabit
+ @param intesectionMax :: The maximum possible value for the overlapping region to inhabit
+ @return a double contianing the start of the overlapping region
+ */
+double Stitch1D::getStartOverlap(const double& intesectionMin, const double& intesectionMax) const
+{
+Property* startOverlapProp = this->getProperty("StartOverlap");
+double startOverlapVal = this->getProperty("StartOverlap");
+startOverlapVal -= this->range_tolerance;
+const bool startOverlapBeyondRange = (startOverlapVal < intesectionMin)
+ || (startOverlapVal > intesectionMax);
+if (startOverlapProp->isDefault() || startOverlapBeyondRange)
+{
+ if (!startOverlapProp->isDefault() && startOverlapBeyondRange)
+ {
+ char message[200];
+ std::sprintf(message,
+ "StartOverlap is outside range at %0.4f, Min is %0.4f, Max is %0.4f . Forced to be: %0.4f",
+ startOverlapVal, intesectionMin, intesectionMax, intesectionMin);
+ g_log.warning(std::string(message));
+ }
+ startOverlapVal = intesectionMin;
+ std::stringstream buffer;
+ buffer << "StartOverlap calculated to be: " << startOverlapVal;
+ g_log.information(buffer.str());
+}
+return startOverlapVal;
+}
- MantidVec result;
- if (areParamsDefault)
- {
- MantidVec calculatedParams;
-
- // Calculate the step size based on the existing step size of the LHS workspace. That way scale factors should be reasonably maintained.
- double calculatedStep =0;
- if(scaleRHS)
- {
- // Calculate the step from the workspace that will not be scaled. The LHS workspace.
- calculatedStep = lhsX[1] - lhsX[0];
- }
- else
- {
- // Calculate the step from the workspace that will not be scaled. The RHS workspace.
- calculatedStep = rhsX[1] - rhsX[0];
- }
-
- calculatedParams.push_back(minLHSX);
- calculatedParams.push_back(calculatedStep);
- calculatedParams.push_back(maxRHSX);
- result = calculatedParams;
- }
- else
- {
- if (inputParams.size() == 1)
- {
- MantidVec calculatedParams;
- calculatedParams.push_back(minLHSX);
- calculatedParams.push_back(inputParams.front()); // Use the step supplied.
- calculatedParams.push_back(maxRHSX);
- result = calculatedParams;
- }
- else
- {
- result = inputParams; // user has provided params. Use those.
- }
- }
- return result;
- }
+/**Gets the end of the overlapping region
+ @param intesectionMin :: The minimum possible value for the overlapping region to inhabit
+ @param intesectionMax :: The maximum possible value for the overlapping region to inhabit
+ @return a double contianing the end of the overlapping region
+ */
+double Stitch1D::getEndOverlap(const double& intesectionMin, const double& intesectionMax) const
+{
+Property* endOverlapProp = this->getProperty("EndOverlap");
+double endOverlapVal = this->getProperty("EndOverlap");
+endOverlapVal += this->range_tolerance;
+const bool endOverlapBeyondRange = (endOverlapVal < intesectionMin) || (endOverlapVal > intesectionMax);
+if (endOverlapProp->isDefault() || endOverlapBeyondRange)
+{
+ if (!endOverlapProp->isDefault() && endOverlapBeyondRange)
+ {
+ char message[200];
+ std::sprintf(message,
+ "EndOverlap is outside range at %0.4f, Min is %0.4f, Max is %0.4f . Forced to be: %0.4f",
+ endOverlapVal, intesectionMin, intesectionMax, intesectionMax);
+ g_log.warning(std::string(message));
+ }
+ endOverlapVal = intesectionMax;
+ std::stringstream buffer;
+ buffer << "EndOverlap calculated to be: " << endOverlapVal;
+ g_log.information(buffer.str());
+}
+return endOverlapVal;
+}
- /**Runs the Rebin Algorithm as a child
- @param input :: The input workspace
- @param params :: a vector<double> containing rebinning parameters
- @return A shared pointer to the resulting MatrixWorkspace
- */
- MatrixWorkspace_sptr Stitch1D::rebin(MatrixWorkspace_sptr& input, const MantidVec& params)
- {
- auto rebin = this->createChildAlgorithm("Rebin");
- rebin->setProperty("InputWorkspace", input);
- rebin->setProperty("Params", params);
- rebin->execute();
- MatrixWorkspace_sptr outWS = rebin->getProperty("OutputWorkspace");
- return outWS;
- }
+/**Gets the rebinning parameters and calculates any missing values
+ @param lhsWS :: The left hand side input workspace
+ @param rhsWS :: The right hand side input workspace
+ @param scaleRHS :: Scale the right hand side workspace
+ @return a vector<double> contianing the rebinning parameters
+ */
+MantidVec Stitch1D::getRebinParams(MatrixWorkspace_sptr& lhsWS, MatrixWorkspace_sptr& rhsWS,
+ const bool scaleRHS) const
+{
+MantidVec inputParams = this->getProperty("Params");
+Property* prop = this->getProperty("Params");
+const bool areParamsDefault = prop->isDefault();
- /**Runs the Integration Algorithm as a child after replacing special values.
- @param input :: The input workspace
- @param start :: a double defining the start of the region to integrate
- @param stop :: a double defining the end of the region to integrate
- @return A shared pointer to the resulting MatrixWorkspace
- */
- MatrixWorkspace_sptr Stitch1D::specialIntegration(MatrixWorkspace_sptr& input, const double& start,
- const double& stop)
- {
- // Effectively ignore values that will trip the integration.
- auto replace = this->createChildAlgorithm("ReplaceSpecialValues");
- replace->setProperty("InputWorkspace", input);
- replace->setProperty("NaNValue", 0.0);
- replace->setProperty("NaNError", 0.0);
- replace->setProperty("InfinityValue", 0.0);
- replace->setProperty("InfinityError", 0.0);
- replace->execute();
- MatrixWorkspace_sptr patchedWS = replace->getProperty("OutputWorkspace");
-
- auto integration = this->createChildAlgorithm("Integration");
- integration->setProperty("InputWorkspace", patchedWS);
- integration->setProperty("RangeLower", start);
- integration->setProperty("RangeUpper", stop);
- integration->execute();
- MatrixWorkspace_sptr outWS = integration->getProperty("OutputWorkspace");
- return outWS;
- }
+const MantidVec& lhsX = lhsWS->readX(0);
+auto it = std::min_element(lhsX.begin(), lhsX.end());
+const double minLHSX = *it;
- /**Runs the MultiplyRange Algorithm as a child defining an end bin
- @param input :: The input workspace
- @param startBin :: The first bin int eh range to multiply
- @param endBin :: The last bin in the range to multiply
- @param factor :: The multiplication factor
- @return A shared pointer to the resulting MatrixWorkspace
- */
- MatrixWorkspace_sptr Stitch1D::multiplyRange(MatrixWorkspace_sptr& input, const int& startBin,
- const int& endBin, const double& factor)
- {
- auto multiplyRange = this->createChildAlgorithm("MultiplyRange");
- multiplyRange->setProperty("InputWorkspace", input);
- multiplyRange->setProperty("StartBin", startBin);
- multiplyRange->setProperty("EndBin", endBin);
- multiplyRange->setProperty("Factor", factor);
- multiplyRange->execute();
- MatrixWorkspace_sptr outWS = multiplyRange->getProperty("OutputWorkspace");
- return outWS;
- }
+const MantidVec& rhsX = rhsWS->readX(0);
+it = std::max_element(rhsX.begin(), rhsX.end());
+const double maxRHSX = *it;
- /**Runs the MultiplyRange Algorithm as a child
- @param input :: The input workspace
- @param startBin :: The first bin int eh range to multiply
- @param factor :: The multiplication factor
- @return A shared pointer to the resulting MatrixWorkspace
- */
- MatrixWorkspace_sptr Stitch1D::multiplyRange(MatrixWorkspace_sptr& input, const int& startBin,
- const double& factor)
- {
- auto multiplyRange = this->createChildAlgorithm("MultiplyRange");
- multiplyRange->setProperty("InputWorkspace", input);
- multiplyRange->setProperty("StartBin", startBin);
- multiplyRange->setProperty("Factor", factor);
- multiplyRange->execute();
- MatrixWorkspace_sptr outWS = multiplyRange->getProperty("OutputWorkspace");
- return outWS;
- }
+MantidVec result;
+if (areParamsDefault)
+{
+ MantidVec calculatedParams;
- /**Runs the WeightedMean Algorithm as a child
- @param inOne :: The first input workspace
- @param inTwo :: The second input workspace
- @return A shared pointer to the resulting MatrixWorkspace
- */
- MatrixWorkspace_sptr Stitch1D::weightedMean(MatrixWorkspace_sptr& inOne, MatrixWorkspace_sptr& inTwo)
- {
- auto weightedMean = this->createChildAlgorithm("WeightedMean");
- weightedMean->setProperty("InputWorkspace1", inOne);
- weightedMean->setProperty("InputWorkspace2", inTwo);
- weightedMean->execute();
- MatrixWorkspace_sptr outWS = weightedMean->getProperty("OutputWorkspace");
- return outWS;
- }
+ // Calculate the step size based on the existing step size of the LHS workspace. That way scale factors should be reasonably maintained.
+ double calculatedStep = 0;
+ if (scaleRHS)
+ {
+ // Calculate the step from the workspace that will not be scaled. The LHS workspace.
+ calculatedStep = lhsX[1] - lhsX[0];
+ }
+ else
+ {
+ // Calculate the step from the workspace that will not be scaled. The RHS workspace.
+ calculatedStep = rhsX[1] - rhsX[0];
+ }
- /**Runs the CreateSingleValuedWorkspace Algorithm as a child
- @param val :: The double to convert to a single value workspace
- @return A shared pointer to the resulting MatrixWorkspace
- */
- MatrixWorkspace_sptr Stitch1D::singleValueWS(double val)
- {
- auto singleValueWS = this->createChildAlgorithm("CreateSingleValuedWorkspace");
- singleValueWS->setProperty("DataValue", val);
- singleValueWS->execute();
- MatrixWorkspace_sptr outWS = singleValueWS->getProperty("OutputWorkspace");
- return outWS;
- }
+ calculatedParams.push_back(minLHSX);
+ calculatedParams.push_back(calculatedStep);
+ calculatedParams.push_back(maxRHSX);
+ result = calculatedParams;
+}
+else
+{
+ if (inputParams.size() == 1)
+ {
+ MantidVec calculatedParams;
+ calculatedParams.push_back(minLHSX);
+ calculatedParams.push_back(inputParams.front()); // Use the step supplied.
+ calculatedParams.push_back(maxRHSX);
+ result = calculatedParams;
+ }
+ else
+ {
+ result = inputParams; // user has provided params. Use those.
+ }
+}
+return result;
+}
- /**finds the bins containing the ends of the overlappign region
- @param startOverlap :: The start of the overlapping region
- @param endOverlap :: The end of the overlapping region
- @param workspace :: The workspace to determine the overlaps inside
- @return a boost::tuple<int,int> containing the bin indexes of the overlaps
- */
- boost::tuple<int, int> Stitch1D::findStartEndIndexes(double startOverlap, double endOverlap,
- MatrixWorkspace_sptr& workspace)
- {
- int a1 = static_cast<int>(workspace->binIndexOf(startOverlap));
- int a2 = static_cast<int>(workspace->binIndexOf(endOverlap));
- if (a1 == a2)
- {
- throw std::runtime_error(
- "The Params you have provided for binning yield a workspace in which start and end overlap appear in the same bin. Make binning finer via input Params.");
- }
- return boost::tuple<int, int>(a1, a2);
+/**Runs the Rebin Algorithm as a child
+ @param input :: The input workspace
+ @param params :: a vector<double> containing rebinning parameters
+ @return A shared pointer to the resulting MatrixWorkspace
+ */
+MatrixWorkspace_sptr Stitch1D::rebin(MatrixWorkspace_sptr& input, const MantidVec& params)
+{
+auto rebin = this->createChildAlgorithm("Rebin");
+rebin->setProperty("InputWorkspace", input);
+rebin->setProperty("Params", params);
+rebin->execute();
+MatrixWorkspace_sptr outWS = rebin->getProperty("OutputWorkspace");
+return outWS;
+}
- }
+/**Runs the Integration Algorithm as a child after replacing special values.
+ @param input :: The input workspace
+ @param start :: a double defining the start of the region to integrate
+ @param stop :: a double defining the end of the region to integrate
+ @return A shared pointer to the resulting MatrixWorkspace
+ */
+MatrixWorkspace_sptr Stitch1D::specialIntegration(MatrixWorkspace_sptr& input, const double& start,
+ const double& stop)
+{
+// Effectively ignore values that will trip the integration.
+auto replace = this->createChildAlgorithm("ReplaceSpecialValues");
+replace->setProperty("InputWorkspace", input);
+replace->setProperty("NaNValue", 0.0);
+replace->setProperty("NaNError", 0.0);
+replace->setProperty("InfinityValue", 0.0);
+replace->setProperty("InfinityError", 0.0);
+replace->execute();
+MatrixWorkspace_sptr patchedWS = replace->getProperty("OutputWorkspace");
+
+auto integration = this->createChildAlgorithm("Integration");
+integration->setProperty("InputWorkspace", patchedWS);
+integration->setProperty("RangeLower", start);
+integration->setProperty("RangeUpper", stop);
+integration->execute();
+MatrixWorkspace_sptr outWS = integration->getProperty("OutputWorkspace");
+return outWS;
+}
- /**Determines if a workspace has non zero errors
- @param ws :: The input workspace
- @return True if there are any non-zero errors in the workspace
- */
- bool Stitch1D::hasNonzeroErrors(MatrixWorkspace_sptr ws)
- {
- int64_t ws_size = static_cast<int64_t>(ws->getNumberHistograms());
- bool hasNonZeroErrors = false;
- PARALLEL_FOR1(ws)
- for (int i = 0; i < ws_size; ++i)
- {
- PARALLEL_START_INTERUPT_REGION
- if (!hasNonZeroErrors) // Keep checking
- {
- auto e = ws->readE(i);
- auto it = std::find_if(e.begin(), e.end(),isNonzero);
- if (it != e.end())
- {
- PARALLEL_CRITICAL(has_non_zero)
- {
- hasNonZeroErrors = true; // Set flag. Should not need to check any more spectra.
- }
- }
- }
- PARALLEL_END_INTERUPT_REGION
- }
- PARALLEL_CHECK_INTERUPT_REGION
+/**Runs the MultiplyRange Algorithm as a child defining an end bin
+ @param input :: The input workspace
+ @param startBin :: The first bin int eh range to multiply
+ @param endBin :: The last bin in the range to multiply
+ @param factor :: The multiplication factor
+ @return A shared pointer to the resulting MatrixWorkspace
+ */
+MatrixWorkspace_sptr Stitch1D::multiplyRange(MatrixWorkspace_sptr& input, const int& startBin,
+ const int& endBin, const double& factor)
+{
+auto multiplyRange = this->createChildAlgorithm("MultiplyRange");
+multiplyRange->setProperty("InputWorkspace", input);
+multiplyRange->setProperty("StartBin", startBin);
+multiplyRange->setProperty("EndBin", endBin);
+multiplyRange->setProperty("Factor", factor);
+multiplyRange->execute();
+MatrixWorkspace_sptr outWS = multiplyRange->getProperty("OutputWorkspace");
+return outWS;
+}
- return hasNonZeroErrors;
- }
+/**Runs the MultiplyRange Algorithm as a child
+ @param input :: The input workspace
+ @param startBin :: The first bin int eh range to multiply
+ @param factor :: The multiplication factor
+ @return A shared pointer to the resulting MatrixWorkspace
+ */
+MatrixWorkspace_sptr Stitch1D::multiplyRange(MatrixWorkspace_sptr& input, const int& startBin,
+ const double& factor)
+{
+auto multiplyRange = this->createChildAlgorithm("MultiplyRange");
+multiplyRange->setProperty("InputWorkspace", input);
+multiplyRange->setProperty("StartBin", startBin);
+multiplyRange->setProperty("Factor", factor);
+multiplyRange->execute();
+MatrixWorkspace_sptr outWS = multiplyRange->getProperty("OutputWorkspace");
+return outWS;
+}
- //----------------------------------------------------------------------------------------------
- /** Execute the algorithm.
- */
- void Stitch1D::exec()
- {
- MatrixWorkspace_sptr rhsWS = this->getProperty("RHSWorkspace");
- MatrixWorkspace_sptr lhsWS = this->getProperty("LHSWorkspace");
- const MinMaxTuple intesectionXRegion = calculateXIntersection(lhsWS, rhsWS);
+/**Runs the WeightedMean Algorithm as a child
+ @param inOne :: The first input workspace
+ @param inTwo :: The second input workspace
+ @return A shared pointer to the resulting MatrixWorkspace
+ */
+MatrixWorkspace_sptr Stitch1D::weightedMean(MatrixWorkspace_sptr& inOne, MatrixWorkspace_sptr& inTwo)
+{
+auto weightedMean = this->createChildAlgorithm("WeightedMean");
+weightedMean->setProperty("InputWorkspace1", inOne);
+weightedMean->setProperty("InputWorkspace2", inTwo);
+weightedMean->execute();
+MatrixWorkspace_sptr outWS = weightedMean->getProperty("OutputWorkspace");
+return outWS;
+}
- const double intersectionMin = intesectionXRegion.get<0>();
- const double intersectionMax = intesectionXRegion.get<1>();
+/**Runs the CreateSingleValuedWorkspace Algorithm as a child
+ @param val :: The double to convert to a single value workspace
+ @return A shared pointer to the resulting MatrixWorkspace
+ */
+MatrixWorkspace_sptr Stitch1D::singleValueWS(double val)
+{
+auto singleValueWS = this->createChildAlgorithm("CreateSingleValuedWorkspace");
+singleValueWS->setProperty("DataValue", val);
+singleValueWS->execute();
+MatrixWorkspace_sptr outWS = singleValueWS->getProperty("OutputWorkspace");
+return outWS;
+}
- const double startOverlap = getStartOverlap(intersectionMin, intersectionMax);
- const double endOverlap = getEndOverlap(intersectionMin, intersectionMax);
+/**finds the bins containing the ends of the overlappign region
+ @param startOverlap :: The start of the overlapping region
+ @param endOverlap :: The end of the overlapping region
+ @param workspace :: The workspace to determine the overlaps inside
+ @return a boost::tuple<int,int> containing the bin indexes of the overlaps
+ */
+boost::tuple<int, int> Stitch1D::findStartEndIndexes(double startOverlap, double endOverlap,
+ MatrixWorkspace_sptr& workspace)
+{
+int a1 = static_cast<int>(workspace->binIndexOf(startOverlap));
+int a2 = static_cast<int>(workspace->binIndexOf(endOverlap));
+if (a1 == a2)
+{
+ throw std::runtime_error(
+ "The Params you have provided for binning yield a workspace in which start and end overlap appear in the same bin. Make binning finer via input Params.");
+}
+return boost::tuple<int, int>(a1, a2);
- if (startOverlap > endOverlap)
+}
+
+/**Determines if a workspace has non zero errors
+ @param ws :: The input workspace
+ @return True if there are any non-zero errors in the workspace
+ */
+bool Stitch1D::hasNonzeroErrors(MatrixWorkspace_sptr ws)
+{
+int64_t ws_size = static_cast<int64_t>(ws->getNumberHistograms());
+bool hasNonZeroErrors = false;
+PARALLEL_FOR1(ws)
+for (int i = 0; i < ws_size; ++i)
+{
+ PARALLEL_START_INTERUPT_REGION
+ if (!hasNonZeroErrors) // Keep checking
+ {
+ auto e = ws->readE(i);
+ auto it = std::find_if(e.begin(), e.end(), isNonzero);
+ if (it != e.end())
{
- std::string message = boost::str(
- boost::format(
- "Stitch1D cannot have a StartOverlap > EndOverlap. StartOverlap: %0.9f, EndOverlap: %0.9f")
- % startOverlap % endOverlap);
- throw std::runtime_error(message);
+ PARALLEL_CRITICAL(has_non_zero)
+ {
+ hasNonZeroErrors = true; // Set flag. Should not need to check any more spectra.
+ }
}
+ }
+PARALLEL_END_INTERUPT_REGION
+}
+PARALLEL_CHECK_INTERUPT_REGION
- const bool scaleRHS = this->getProperty("ScaleRHSWorkspace");
- MantidVec params = getRebinParams(lhsWS, rhsWS, scaleRHS);
+return hasNonZeroErrors;
+}
- const double& xMin = params.front();
- const double& xMax = params.back();
+//----------------------------------------------------------------------------------------------
+/** Execute the algorithm.
+ */
+void Stitch1D::exec()
+{
+MatrixWorkspace_sptr rhsWS = this->getProperty("RHSWorkspace");
+MatrixWorkspace_sptr lhsWS = this->getProperty("LHSWorkspace");
+const MinMaxTuple intesectionXRegion = calculateXIntersection(lhsWS, rhsWS);
- if (startOverlap < xMin)
- {
- std::string message =
- boost::str(
- boost::format(
- "Stitch1D StartOverlap is outside the available X range after rebinning. StartOverlap: %10.9f, X min: %10.9f")
- % startOverlap % xMin);
+const double intersectionMin = intesectionXRegion.get<0>();
+const double intersectionMax = intesectionXRegion.get<1>();
- throw std::runtime_error(message);
- }
- if (endOverlap > xMax)
- {
- std::string message =
- boost::str(
- boost::format(
- "Stitch1D EndOverlap is outside the available X range after rebinning. EndOverlap: %10.9f, X max: %10.9f")
- % endOverlap % xMax);
+const double startOverlap = getStartOverlap(intersectionMin, intersectionMax);
+const double endOverlap = getEndOverlap(intersectionMin, intersectionMax);
- throw std::runtime_error(message);
- }
+if (startOverlap > endOverlap)
+{
+std::string message = boost::str(
+ boost::format(
+ "Stitch1D cannot have a StartOverlap > EndOverlap. StartOverlap: %0.9f, EndOverlap: %0.9f")
+ % startOverlap % endOverlap);
+throw std::runtime_error(message);
+}
- auto rebinnedLHS = rebin(lhsWS, params);
- auto rebinnedRHS = rebin(rhsWS, params);
+const bool scaleRHS = this->getProperty("ScaleRHSWorkspace");
+MantidVec params = getRebinParams(lhsWS, rhsWS, scaleRHS);
- boost::tuple<int, int> startEnd = findStartEndIndexes(startOverlap, endOverlap, rebinnedLHS);
+const double& xMin = params.front();
+const double& xMax = params.back();
- const bool useManualScaleFactor = this->getProperty("UseManualScaleFactor");
- double scaleFactor = 0;
- double errorScaleFactor = 0;
+if (startOverlap < xMin)
+{
+std::string message =
+ boost::str(
+ boost::format(
+ "Stitch1D StartOverlap is outside the available X range after rebinning. StartOverlap: %10.9f, X min: %10.9f")
+ % startOverlap % xMin);
- if (useManualScaleFactor)
- {
- double manualScaleFactor = this->getProperty("ManualScaleFactor");
- MatrixWorkspace_sptr manualScaleFactorWS = singleValueWS(manualScaleFactor);
+throw std::runtime_error(message);
+}
+if (endOverlap > xMax)
+{
+std::string message =
+ boost::str(
+ boost::format(
+ "Stitch1D EndOverlap is outside the available X range after rebinning. EndOverlap: %10.9f, X max: %10.9f")
+ % endOverlap % xMax);
- if (scaleRHS)
- {
- rebinnedRHS = rebinnedRHS * manualScaleFactorWS;
- }
- else
- {
- rebinnedLHS = rebinnedLHS * manualScaleFactorWS;
- }
- scaleFactor = manualScaleFactor;
- errorScaleFactor = manualScaleFactor;
- }
- else
- {
- auto rhsOverlapIntegrated = specialIntegration(rebinnedRHS, startOverlap, endOverlap);
- auto lhsOverlapIntegrated = specialIntegration(rebinnedLHS, startOverlap, endOverlap);
+throw std::runtime_error(message);
+}
- MatrixWorkspace_sptr ratio;
- if (scaleRHS)
- {
- ratio = lhsOverlapIntegrated / rhsOverlapIntegrated;
- rebinnedRHS = rebinnedRHS * ratio;
- }
- else
- {
- ratio = rhsOverlapIntegrated / lhsOverlapIntegrated;
- rebinnedLHS = rebinnedLHS * ratio;
- }
- scaleFactor = ratio->readY(0).front();
- errorScaleFactor = ratio->readE(0).front();
- if(scaleFactor < 1e-2 || scaleFactor > 1e2 || boost::math::isnan(scaleFactor))
- {
- std::stringstream messageBuffer;
- messageBuffer << "Stitch1D calculated scale factor is: " << scaleFactor << ". Check that in both input workspaces the integrated overlap region is non-zero.";
- g_log.warning(messageBuffer.str());
- }
+auto rebinnedLHS = rebin(lhsWS, params);
+auto rebinnedRHS = rebin(rhsWS, params);
- }
+boost::tuple<int, int> startEnd = findStartEndIndexes(startOverlap, endOverlap, rebinnedLHS);
- int a1 = boost::tuples::get<0>(startEnd);
- int a2 = boost::tuples::get<1>(startEnd);
+const bool useManualScaleFactor = this->getProperty("UseManualScaleFactor");
+double scaleFactor = 0;
+double errorScaleFactor = 0;
- // Mask out everything BUT the overlap region as a new workspace.
- MatrixWorkspace_sptr overlap1 = multiplyRange(rebinnedLHS, 0, a1, 0);
- overlap1 = multiplyRange(overlap1, a2, 0);
+if (useManualScaleFactor)
+{
+double manualScaleFactor = this->getProperty("ManualScaleFactor");
+MatrixWorkspace_sptr manualScaleFactorWS = singleValueWS(manualScaleFactor);
+
+if (scaleRHS)
+{
+ rebinnedRHS = rebinnedRHS * manualScaleFactorWS;
+}
+else
+{
+ rebinnedLHS = rebinnedLHS * manualScaleFactorWS;
+}
+scaleFactor = manualScaleFactor;
+errorScaleFactor = manualScaleFactor;
+}
+else
+{
+auto rhsOverlapIntegrated = specialIntegration(rebinnedRHS, startOverlap, endOverlap);
+auto lhsOverlapIntegrated = specialIntegration(rebinnedLHS, startOverlap, endOverlap);
- // Mask out everything BUT the overlap region as a new workspace.
- MatrixWorkspace_sptr overlap2 = multiplyRange(rebinnedRHS, 0, a1, 0);
- overlap2 = multiplyRange(overlap2, a2, 0);
+MatrixWorkspace_sptr ratio;
+if (scaleRHS)
+{
+ ratio = lhsOverlapIntegrated / rhsOverlapIntegrated;
+ rebinnedRHS = rebinnedRHS * ratio;
+}
+else
+{
+ ratio = rhsOverlapIntegrated / lhsOverlapIntegrated;
+ rebinnedLHS = rebinnedLHS * ratio;
+}
+scaleFactor = ratio->readY(0).front();
+errorScaleFactor = ratio->readE(0).front();
+if (scaleFactor < 1e-2 || scaleFactor > 1e2 || boost::math::isnan(scaleFactor))
+{
+ std::stringstream messageBuffer;
+ messageBuffer << "Stitch1D calculated scale factor is: " << scaleFactor
+ << ". Check that in both input workspaces the integrated overlap region is non-zero.";
+ g_log.warning(messageBuffer.str());
+}
- // Mask out everything AFTER the start of the overlap region
- rebinnedLHS = multiplyRange(rebinnedLHS, a1 + 1, 0);
+}
- // Mask out everything BEFORE the end of the overlap region
- rebinnedRHS = multiplyRange(rebinnedRHS, 0, a2 - 1, 0);
+int a1 = boost::tuples::get<0>(startEnd);
+int a2 = boost::tuples::get<1>(startEnd);
- // Calculate a weighted mean for the overlap region
+// Mask out everything BUT the overlap region as a new workspace.
+MatrixWorkspace_sptr overlap1 = maskAllBut(a1, a2, rebinnedLHS);
+// Mask out everything BUT the overlap region as a new workspace.
+MatrixWorkspace_sptr overlap2 = maskAllBut(a1, a2, rebinnedRHS);
+// Mask out everything AFTER the overlap region as a new workspace.
+maskInPlace(a1 + 1, rebinnedLHS->blocksize(), rebinnedLHS);
+// Mask out everything BEFORE the overlap region as a new workspace.
+maskInPlace(0, a2, rebinnedRHS);
- MatrixWorkspace_sptr overlapave;
- if (hasNonzeroErrors(overlap1) && hasNonzeroErrors(overlap2))
- {
- overlapave = weightedMean(overlap1, overlap2);
- }
- else
- {
- g_log.information("Using un-weighted mean for Stitch1D overlap mean");
- MatrixWorkspace_sptr sum = overlap1 + overlap2;
- MatrixWorkspace_sptr denominator = singleValueWS(2.0);
- overlapave = sum / denominator;
- }
+MatrixWorkspace_sptr overlapave;
+if (hasNonzeroErrors(overlap1) && hasNonzeroErrors(overlap2))
+{
+overlapave = weightedMean(overlap1, overlap2);
+}
+else
+{
+g_log.information("Using un-weighted mean for Stitch1D overlap mean");
+MatrixWorkspace_sptr sum = overlap1 + overlap2;
+MatrixWorkspace_sptr denominator = singleValueWS(2.0);
+overlapave = sum / denominator;
+}
- MatrixWorkspace_sptr result = rebinnedLHS + overlapave + rebinnedRHS;
+MatrixWorkspace_sptr result = rebinnedLHS + overlapave + rebinnedRHS;
- // Provide log information about the scale factors used in the calculations.
- std::stringstream messageBuffer;
- messageBuffer << "Scale Factor Y is: " << scaleFactor << " Scale Factor E is: " << errorScaleFactor;
- g_log.notice(messageBuffer.str());
+// Provide log information about the scale factors used in the calculations.
+std::stringstream messageBuffer;
+messageBuffer << "Scale Factor Y is: " << scaleFactor << " Scale Factor E is: " << errorScaleFactor;
+g_log.notice(messageBuffer.str());
- setProperty("OutputWorkspace", result);
- setProperty("OutScaleFactor", scaleFactor);
+setProperty("OutputWorkspace", result);
+setProperty("OutScaleFactor", scaleFactor);
- }
+}
} // namespace Algorithms
} // namespace Mantid