Skip to content
Snippets Groups Projects
Commit 70bc199c authored by Verena Reimund's avatar Verena Reimund
Browse files

Simpler interface + docs, separation histograms & point data 

Refs #22197
parent 0540ee1f
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Framework/Algorithms/src/Stitch1D.cpp
+ 53
− 47
View file @ 70bc199c
#include "MantidAlgorithms/Stitch1D.h"
#include "MantidAPI/WorkspaceProperty.h"
#include "MantidAPI/MatrixWorkspace.h"
#include "MantidAPI/HistogramValidator.h"
#include "MantidAPI/WorkspaceFactory.h"
#include "MantidHistogramData/HistogramX.h"
#include "MantidHistogramData/HistogramY.h"
#include "MantidHistogramData/HistogramE.h"
#include "MantidKernel/ArrayProperty.h"
#include "MantidKernel/BoundedValidator.h"
#include "MantidKernel/ListValidator.h"
#include "MantidKernel/MultiThreaded.h"
#include "MantidKernel/PropertyWithValue.h"
#include "MantidKernel/RebinParamsValidator.h"
......@@ -118,17 +116,13 @@ void Stitch1D::maskInPlace(int a1, int a2, MatrixWorkspace_sptr source) {
/** Initialize the algorithm's properties.
*/
void Stitch1D::init() {
Kernel::IValidator_sptr histogramValidator =
boost::make_shared<HistogramValidator>();
declareProperty(
make_unique<WorkspaceProperty<MatrixWorkspace>>(
"LHSWorkspace", "", Direction::Input, histogramValidator->clone()),
"LHS input workspace.");
declareProperty(
make_unique<WorkspaceProperty<MatrixWorkspace>>(
"RHSWorkspace", "", Direction::Input, histogramValidator->clone()),
"RHS input workspace.");
declareProperty(make_unique<WorkspaceProperty<MatrixWorkspace>>(
"LHSWorkspace", "", Direction::Input),
"LHS input workspace.");
declareProperty(make_unique<WorkspaceProperty<MatrixWorkspace>>(
"RHSWorkspace", "", Direction::Input),
"RHS input workspace.");
declareProperty(make_unique<WorkspaceProperty<MatrixWorkspace>>(
"OutputWorkspace", "", Direction::Output),
"Output stitched workspace.");
......@@ -161,12 +155,6 @@ void Stitch1D::init() {
declareProperty(make_unique<PropertyWithValue<double>>(
"OutScaleFactor", Mantid::EMPTY_DBL(), Direction::Output),
"The actual used value for the scaling factor.");
const std::vector<std::string> outputXOption{"WeightedMeanOverlap",
"OriginalOverlap"};
declareProperty("OutputXOption", "WeightedMeanOverlap",
boost::make_shared<StringListValidator>(outputXOption),
"Output option to choose how to overlap x values.");
}
/**Gets the start of the overlapping region
......@@ -395,7 +383,7 @@ MatrixWorkspace_sptr Stitch1D::singleValueWS(double val) {
return outWS;
}
/**finds the bins containing the ends of the overlappign region
/**finds the bins containing the ends of the overlapping 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
......@@ -449,6 +437,7 @@ bool Stitch1D::hasNonzeroErrors(MatrixWorkspace_sptr ws) {
void Stitch1D::exec() {
MatrixWorkspace_sptr rhsWS = this->getProperty("RHSWorkspace");
MatrixWorkspace_sptr lhsWS = this->getProperty("LHSWorkspace");
const MinMaxTuple intesectionXRegion = calculateXIntersection(lhsWS, rhsWS);
const double intersectionMin = intesectionXRegion.get<0>();
......@@ -482,9 +471,9 @@ void Stitch1D::exec() {
boost::format("Stitch1D StartOverlap is outside the available X range "
"after rebinning. StartOverlap: %10.9f, X min: %10.9f") %
startOverlap % xMin);
throw std::runtime_error(message);
}
if (endOverlap > xMax) {
std::string message = boost::str(
boost::format("Stitch1D EndOverlap is outside the available X range "
......@@ -499,40 +488,40 @@ void Stitch1D::exec() {
m_infYIndexes.resize(histogramCount);
m_nanEIndexes.resize(histogramCount);
m_infEIndexes.resize(histogramCount);
auto rebinnedLHS = rebin(lhsWS, params);
auto rebinnedRHS = rebin(rhsWS, params);
boost::tuple<int, int> startEnd =
findStartEndIndexes(startOverlap, endOverlap, rebinnedLHS);
const bool useManualScaleFactor = this->getProperty("UseManualScaleFactor");
double scaleFactor = 0;
double errorScaleFactor = 0;
MatrixWorkspace_sptr lhs, rhs;
// If the input workspaces are histograms ...
if (rhsWS->isHistogramData() && lhsWS->isHistogramData()) {
lhs = rebin(lhsWS, params);
rhs = rebin(rhsWS, params);
}
if (useManualScaleFactor) {
double manualScaleFactor = this->getProperty("ManualScaleFactor");
MatrixWorkspace_sptr manualScaleFactorWS = singleValueWS(manualScaleFactor);
if (scaleRHS) {
rebinnedRHS = rebinnedRHS * manualScaleFactorWS;
} else {
rebinnedLHS = rebinnedLHS * manualScaleFactorWS;
}
if (scaleRHS)
rhs *= manualScaleFactorWS;
else
lhs *= manualScaleFactorWS;
scaleFactor = manualScaleFactor;
errorScaleFactor = manualScaleFactor;
} else {
auto rhsOverlapIntegrated =
integration(rebinnedRHS, startOverlap, endOverlap);
auto lhsOverlapIntegrated =
integration(rebinnedLHS, startOverlap, endOverlap);
auto rhsOverlapIntegrated = integration(rhs, startOverlap, endOverlap);
auto lhsOverlapIntegrated = integration(lhs, startOverlap, endOverlap);
MatrixWorkspace_sptr ratio;
if (scaleRHS) {
ratio = lhsOverlapIntegrated / rhsOverlapIntegrated;
rebinnedRHS = rebinnedRHS * ratio;
rhs = rhs * ratio;
} else {
ratio = rhsOverlapIntegrated / lhsOverlapIntegrated;
rebinnedLHS = rebinnedLHS * ratio;
lhs *= ratio;
}
scaleFactor = ratio->y(0).front();
errorScaleFactor = ratio->e(0).front();
......@@ -545,29 +534,46 @@ void Stitch1D::exec() {
}
}
boost::tuple<int, int> startEnd =
findStartEndIndexes(startOverlap, endOverlap, lhs);
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 = maskAllBut(a1, a2, rebinnedLHS);
MatrixWorkspace_sptr overlap1 = maskAllBut(a1, a2, lhs);
// Mask out everything BUT the overlap region as a new workspace.
MatrixWorkspace_sptr overlap2 = maskAllBut(a1, a2, rebinnedRHS);
MatrixWorkspace_sptr overlap2 = maskAllBut(a1, a2, rhs);
// Mask out everything AFTER the overlap region as a new workspace.
maskInPlace(a1 + 1, static_cast<int>(rebinnedLHS->blocksize()), rebinnedLHS);
maskInPlace(a1 + 1, static_cast<int>(lhs->blocksize()), lhs);
// Mask out everything BEFORE the overlap region as a new workspace.
maskInPlace(0, a2, rebinnedRHS);
maskInPlace(0, a2, rhs);
MatrixWorkspace_sptr overlapave;
if (hasNonzeroErrors(overlap1) && hasNonzeroErrors(overlap2)) {
overlapave = weightedMean(overlap1, overlap2);
// If the input workspaces are histograms ...
if (rhsWS->isHistogramData() && lhsWS->isHistogramData()) {
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;
}
} 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;
// If the input workspaces are point data ...
if (!rhsWS->isHistogramData() && !lhsWS->isHistogramData()) {
g_log.error("Missing implementation");
throw std::invalid_argument("Missing implementation");
} else {
g_log.error("Input workspaces must be both histograms or point data");
throw std::invalid_argument(
"Input workspaces must be both histograms or point data");
}
}
MatrixWorkspace_sptr result = rebinnedLHS + overlapave + rebinnedRHS;
MatrixWorkspace_sptr result = lhs + overlapave + rhs;
reinsertSpecialValues(result);
// Provide log information about the scale factors used in the calculations.
......
Framework/Algorithms/test/Stitch1DTest.h
+ 5
− 4
View file @ 70bc199c
......@@ -236,7 +236,8 @@ public:
void test_startoverlap_greater_than_end_overlap_throws() {
TSM_ASSERT_THROWS("Should have thrown with StartOverlap < x max",
do_stitch1D(this->a, this->b, this->x.back(),
this->x.front(), std::vector<double>(1., 0.2)),
this->x.front(),
std::vector<double>(1., 0.2)),
std::runtime_error &);
}
......@@ -244,7 +245,7 @@ public:
auto lhs_ws = make_arbitrary_point_ws();
auto rhs_ws = make_arbitrary_histogram_ws();
TSM_ASSERT_THROWS(
"LHS WS must be histogram",
"Both ws must be either histogram or point data",
do_stitch1D(lhs_ws, rhs_ws, -1., 1., std::vector<double>(1., 0.2)),
std::invalid_argument &);
}
......@@ -253,12 +254,12 @@ public:
auto lhs_ws = make_arbitrary_histogram_ws();
auto rhs_ws = make_arbitrary_point_ws();
TSM_ASSERT_THROWS(
"RHS WS must be histogram",
"Both ws must be either histogram or point data",
do_stitch1D(lhs_ws, rhs_ws, -1., 1., std::vector<double>(1., 0.2)),
std::invalid_argument &);
}
void test_stitching_uses_suppiled_params() {
void test_stitching_uses_supplied_params() {
std::vector<double> params = {-0.8, 0.2, 1.0};
auto ret = do_stitch1D(this->b, this->a, -0.4, 0.4, params);
......
docs/source/algorithms/Stitch1D-v3.rst
+ 5
− 6
View file @ 70bc199c
......@@ -18,9 +18,7 @@ Users can optionally provide :ref:`algm-Rebin` :literal:`Params`, otherwise they
Likewise, :literal:`StartOverlap` and :literal:`EndOverlap` are optional. If not provided, then these
are taken to be the region of X-axis intersection.
The option :literal:`OutputXOption` defines the calculation of x values in the overlap range of the output workspace.
The default option is :literal:`WeightedMeanOverlap`. If the input workspace contains Dx values, the output workspace will not have dx values defined.
Contrarily, the option :literal:`OriginalOverlap` will propagate the corresponding Dx values to the output workspace.
The type of the input workspaces (histogram or point data) define the calculation of x values in the overlap range of the output workspace.
The algorithm workflow is as follows:
......@@ -46,13 +44,14 @@ The algorithm workflow is as follows:
#. The weighted mean of the two workspaces in range [:literal:`StartOverlap`, :literal:`EndOverlap`]
is calculated. Note that if both workspaces have zero errors, an un-weighted mean will be
performed instead.
#. If :literal:`OutputXOption` is :literal:`WeightedMeanOverlap`, the output workspace will be created by summing the left-hand-side workspace (values in range
#. For histograms, the output workspace will be created by summing the left-hand-side workspace (values in range
[:literal:`StartX`, :literal:`StartOverlap`], where :literal:`StartX` is the minimum X value
specified via :literal:`Params` or calculated from the left-hand-side workspace) + weighted
mean workspace + right-hand-side workspace (values in range [:literal:`EndOverlap`, :literal:`EndX`],
where :literal:`EndX` is the maximum X value specified via :literal:`Params` or calculated
from the right-hand-side workspace) multiplied by the scale factor.
If :literal:`OutputXOption` is :literal:`OriginalOverlap`, the output workspace will propagate x values and Dx values present in the overlap range to the :literal:`OutputWorkspace`.
Dx values will not be present in the output workspace.
For point data, the output workspace will propagate x values and Dx values present in the overlap range to the output workspace.
#. The special values are put back in the output workspace. Note that if both the left-hand-side
workspace and the right-hand-side workspace happen to have a different special value in the same bin, this
bin will be set to infinite in the output workspace.
......@@ -83,7 +82,7 @@ Usage
import numpy as np
def gaussian(x, mu, sigma):
"""Creates a gaussian peak centered on mu and with width sigma."""
"""Creates a Gaussian peak centered on mu and with width sigma."""
return (1/ sigma * np.sqrt(2 * np.pi)) * np.exp( - (x-mu)**2 / (2*sigma**2))
#create two histograms with a single peak in each one
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment