diff --git a/Code/Mantid/Framework/Algorithms/inc/MantidAlgorithms/Stitch1D.h b/Code/Mantid/Framework/Algorithms/inc/MantidAlgorithms/Stitch1D.h
index c33366a2b7a394916bc67ef705e39d0bf9e91ae4..1c5bf9878a4cd065c8354771e428067d2cb06758 100644
--- a/Code/Mantid/Framework/Algorithms/inc/MantidAlgorithms/Stitch1D.h
+++ b/Code/Mantid/Framework/Algorithms/inc/MantidAlgorithms/Stitch1D.h
@@ -3,6 +3,7 @@
#include "MantidKernel/System.h"
#include "MantidAPI/Algorithm.h"
+#include <boost/tuple/tuple.hpp>
namespace Mantid
{
@@ -47,9 +48,15 @@ namespace Algorithms
void exec();
double getStartOverlap(const double& min, const double& max) const;
double getEndOverlap(const double& min, const double& max) const;
+ bool hasNonzeroErrors(Mantid::API::MatrixWorkspace_sptr& ws) const;
Mantid::MantidVec getRebinParams(Mantid::API::MatrixWorkspace_sptr& lhsWS, Mantid::API::MatrixWorkspace_sptr& rhsWS) const;
Mantid::API::MatrixWorkspace_sptr rebin(Mantid::API::MatrixWorkspace_sptr& input, const Mantid::MantidVec& params);
Mantid::API::MatrixWorkspace_sptr integration(Mantid::API::MatrixWorkspace_sptr& input, const double& start, const double& stop);
+ 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 double& factor);
+ Mantid::API::MatrixWorkspace_sptr singleValueWS(double val);
+ Mantid::API::MatrixWorkspace_sptr weightedMean(Mantid::API::MatrixWorkspace_sptr& inOne, Mantid::API::MatrixWorkspace_sptr& inTwo);
+ boost::tuple<int,int> findStartEndIndexes(double startOverlap, double endOverlap, Mantid::API::MatrixWorkspace_sptr& workspace);
};
diff --git a/Code/Mantid/Framework/Algorithms/src/Stitch1D.cpp b/Code/Mantid/Framework/Algorithms/src/Stitch1D.cpp
index 20f62fff20e7fad7a090333c41568fa3aa0f2077..95b9bc6777e202ce042e95244f0fa4d6a4f6e8a6 100644
--- a/Code/Mantid/Framework/Algorithms/src/Stitch1D.cpp
+++ b/Code/Mantid/Framework/Algorithms/src/Stitch1D.cpp
@@ -31,6 +31,11 @@ namespace
MantidVec rhs_x = rhsWS->readX(0);
return MinMaxTuple(rhs_x.front(), lhs_x.back());
}
+
+ bool isNonzero (double i)
+ {
+ return (0 != i);
+ }
}
namespace Mantid
@@ -159,7 +164,7 @@ namespace Mantid
}
endOverlapVal = intesectionMax;
std::stringstream buffer;
- buffer << "StartOverlap calculated to be: " << endOverlapVal;
+ buffer << "EndOverlap calculated to be: " << endOverlapVal;
g_log.information(buffer.str());
}
return endOverlapVal;
@@ -229,6 +234,75 @@ namespace Mantid
return outWS;
}
+ 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;
+ }
+
+ 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;
+ }
+
+ 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;
+ }
+
+ 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;
+ }
+
+ 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);
+
+ }
+
+ bool Stitch1D::hasNonzeroErrors(MatrixWorkspace_sptr& ws) const
+ {
+ size_t ws_size = ws->getNumberHistograms();
+ for (size_t i = 0; i < ws_size; ++i)
+ {
+ auto e = ws->readE(i);
+ std::vector<double, std::allocator<double> >::iterator error = std::find_if(e.begin(), e.end(), isNonzero);
+ if (error != e.end())
+ {
+ return true;
+ }
+ }
+ return false;
+ }
+
//----------------------------------------------------------------------------------------------
/** Execute the algorithm.
*/
@@ -236,7 +310,7 @@ namespace Mantid
{
MatrixWorkspace_sptr rhsWS = this->getProperty("RHSWorkspace");
MatrixWorkspace_sptr lhsWS = this->getProperty("LHSWorkspace");
- const MinMaxTuple intesectionXRegion = calculateXIntersection(rhsWS, lhsWS);
+ const MinMaxTuple intesectionXRegion = calculateXIntersection(lhsWS, rhsWS);
const double intersectionMin = intesectionXRegion.get<0>();
const double intersectionMax = intesectionXRegion.get<1>();
@@ -283,6 +357,10 @@ namespace Mantid
auto rebinnedLHS = rebin(lhsWS, params);
auto rebinnedRHS = rebin(rhsWS, params);
+ boost::tuple<int,int> startEnd = findStartEndIndexes(startOverlap, endOverlap, rebinnedLHS);
+
+ //manualscalefactor if
+
auto rhsOverlapIntegrated = integration(rebinnedRHS, startOverlap, endOverlap);
auto lhsOverlapIntegrated = integration(rebinnedLHS, startOverlap, endOverlap);
@@ -302,7 +380,43 @@ namespace Mantid
scaleFactor = y2[0]/y1[0];
}
- setProperty("OutputWorkspace", rhsWS);
+ //manualscalefactor end if
+
+ int a1 = boost::tuples::get<0>(startEnd);
+ int a2 = boost::tuples::get<1>(startEnd);
+
+ // Mask out everything BUT the overlap region as a new workspace.
+ MatrixWorkspace_sptr overlap1 = multiplyRange(rebinnedLHS,0,a1,0);
+ overlap1 = multiplyRange(overlap1,a2,0);
+
+ // Mask out everything BUT the overlap region as a new workspace.
+ MatrixWorkspace_sptr overlap2 = multiplyRange(rebinnedRHS,0,a1,0);
+ overlap2 = multiplyRange(overlap2,a2,0);
+
+ // 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);
+
+ // Calculate a weighted mean for the overlap region
+
+ 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;
+
+ setProperty("OutputWorkspace", result);
setProperty("OutScaleFactor", scaleFactor);
}
diff --git a/Code/Mantid/Framework/Algorithms/test/Stitch1DTest.h b/Code/Mantid/Framework/Algorithms/test/Stitch1DTest.h
index 418735d753bb9fea5206e51a39963b28007842b8..1f8aba562cfb47943e6713a1b9ae2754c4f8f936 100644
--- a/Code/Mantid/Framework/Algorithms/test/Stitch1DTest.h
+++ b/Code/Mantid/Framework/Algorithms/test/Stitch1DTest.h
@@ -57,6 +57,21 @@ private:
return outWS;
}
+ MatrixWorkspace_sptr create1DWorkspace(MantidVec& xData, MantidVec& yData)
+ {
+ auto createWorkspace = AlgorithmManager::Instance().create("CreateWorkspace");
+ createWorkspace->setChild(true);
+ createWorkspace->initialize();
+ createWorkspace->setProperty("UnitX", "1/q");
+ createWorkspace->setProperty("DataX", xData);
+ createWorkspace->setProperty("DataY", yData);
+ createWorkspace->setProperty("NSpec", 1);
+ createWorkspace->setPropertyValue("OutputWorkspace", "dummy");
+ createWorkspace->execute();
+ MatrixWorkspace_sptr outWS = createWorkspace->getProperty("OutputWorkspace");
+ return outWS;
+ }
+
MatrixWorkspace_sptr a;
MatrixWorkspace_sptr b;
MantidVec x;
@@ -154,6 +169,22 @@ public:
return ResultType(stitched, scaleFactor);
}
+ ResultType do_stitch1D(MatrixWorkspace_sptr& lhs, MatrixWorkspace_sptr& rhs, const MantidVec& params)
+ {
+ Stitch1D alg;
+ alg.setChild(true);
+ alg.setRethrows(true);
+ alg.initialize();
+ alg.setProperty("LHSWorkspace", lhs);
+ alg.setProperty("RHSWorkspace", rhs);
+ alg.setProperty("Params", params);
+ alg.setPropertyValue("OutputWorkspace", "dummy_value");
+ alg.execute();
+ MatrixWorkspace_sptr stitched = alg.getProperty("OutputWorkspace");
+ double scaleFactor = alg.getProperty("OutScaleFactor");
+ return ResultType(stitched, scaleFactor);
+ }
+
void test_Init()
{
Stitch1D alg;
@@ -172,7 +203,7 @@ public:
{
auto lhs_ws = make_arbitrary_point_ws();
auto rhs_ws = make_arbitrary_histogram_ws();
- TSM_ASSERT_THROWS("LHS WS must be histogram", do_stitch1D(lhs_ws, rhs_ws, -1, 1, MantidVec(0.2)),
+ TSM_ASSERT_THROWS("LHS WS must be histogram", do_stitch1D(lhs_ws, rhs_ws, -1, 1, MantidVec(1, 0.2)),
std::invalid_argument&);
}
@@ -180,7 +211,7 @@ public:
{
auto lhs_ws = make_arbitrary_histogram_ws();
auto rhs_ws = make_arbitrary_point_ws();
- TSM_ASSERT_THROWS("RHS WS must be histogram", do_stitch1D(lhs_ws, rhs_ws, -1, 1, MantidVec(0.2)),
+ TSM_ASSERT_THROWS("RHS WS must be histogram", do_stitch1D(lhs_ws, rhs_ws, -1, 1, MantidVec(1, 0.2)),
std::invalid_argument&);
}
@@ -198,6 +229,52 @@ public:
TS_ASSERT_EQUALS(xMax, 1.0);
}
+ void test_stitching_determines_params()
+ {
+ MantidVec x1 = boost::assign::list_of(-1.0)(-0.8)(-0.6)(-0.4)(-0.2)(0.0)(0.2)(0.4)(0.6)(0.8).convert_to_container<MantidVec>();
+ MantidVec x2 = boost::assign::list_of(0.4)(0.6)(0.8)(1.0)(1.2)(1.4)(1.6).convert_to_container<MantidVec>();
+
+ MatrixWorkspace_sptr ws1 = create1DWorkspace(x1, boost::assign::list_of(1)(1)(1)(1)(1)(1)(1)(1)(1).convert_to_container<MantidVec>());
+ MatrixWorkspace_sptr ws2 = create1DWorkspace(x2, boost::assign::list_of(1)(1)(1)(1)(1)(1).convert_to_container<MantidVec>());
+ double demanded_step_size = 0.2;
+ auto ret = do_stitch1D(ws1,ws2,0.4,1.0,boost::assign::list_of(demanded_step_size).convert_to_container<MantidVec>());
+
+ //Check the ranges on the output workspace against the param inputs.
+ MantidVec out_x_values = ret.get<0>()->readX(0);
+ double x_min = *std::min_element(out_x_values.begin(), out_x_values.end());
+ double x_max = *std::max_element(out_x_values.begin(), out_x_values.end());
+ double step_size = out_x_values[1] - out_x_values[0];
+
+ TS_ASSERT_EQUALS(x_min, -1);
+ TS_ASSERT_DELTA(x_max - demanded_step_size, 1.4, 0.000001);
+ TS_ASSERT_DELTA(step_size, demanded_step_size, 0.000001);
+ }
+
+ void test_stitching_determines_start_and_end_overlap()
+ {
+ MantidVec x1 = boost::assign::list_of(-1.0)(-0.8)(-0.6)(-0.4)(-0.2)(0.0)(0.2)(0.4).convert_to_container<MantidVec>();
+ MantidVec x2 = boost::assign::list_of(-0.4)(-0.2)(0.0)(0.2)(0.4)(0.6)(0.8)(1.0).convert_to_container<MantidVec>();
+ MatrixWorkspace_sptr ws1 = create1DWorkspace(x1, boost::assign::list_of(1)(1)(1)(3)(3)(3)(3).convert_to_container<MantidVec>());
+ MatrixWorkspace_sptr ws2 = create1DWorkspace(x2, boost::assign::list_of(1)(1)(1)(1)(3)(3)(3).convert_to_container<MantidVec>());
+
+ auto ret = do_stitch1D(ws1,ws2,boost::assign::list_of(-1.0)(0.2)(1.0).convert_to_container<MantidVec>());
+
+ MantidVec stitched_y = ret.get<0>()->readY(0);
+ MantidVec stitched_x = ret.get<0>()->readX(0);
+ std::vector<size_t> overlap_indexes = std::vector<size_t>();
+ for (size_t itr = 0; itr < stitched_y.size(); ++itr)
+ {
+ if (stitched_y[itr] >= 1.0009 && stitched_y[itr] <= 3.0001)
+ {
+ overlap_indexes.push_back(itr);
+ }
+ }
+
+ double start_overlap_determined = stitched_x[overlap_indexes[0]];
+ double end_overlap_determined = stitched_x[overlap_indexes[overlap_indexes.size()-1]];
+ TS_ASSERT_DELTA(start_overlap_determined, -0.4, 0.000000001);
+ TS_ASSERT_DELTA(end_overlap_determined, 0.2, 0.000000001);
+ }
};
#endif /* MANTID_ALGORITHMS_STITCH1DTEST_H_ */