diff --git a/Framework/Algorithms/inc/MantidAlgorithms/SofQCommon.h b/Framework/Algorithms/inc/MantidAlgorithms/SofQCommon.h
index befd505daf40d3e5fc51358c67382bd90a8ab853..60248ddf8ea22992f11d68f63c222e7f7f13e2b5 100644
--- a/Framework/Algorithms/inc/MantidAlgorithms/SofQCommon.h
+++ b/Framework/Algorithms/inc/MantidAlgorithms/SofQCommon.h
@@ -32,6 +32,9 @@ struct DLLExport SofQCommon {
/// Estimate minimum and maximum momentum transfer.
std::pair<double, double> qBinHints(const API::MatrixWorkspace &ws, const double minE, const double maxE) const;
+private:
+ std::pair<double, double> qBinHintsDirect(const API::MatrixWorkspace &ws, const double minE, const double maxE) const;
+ std::pair<double, double> qBinHintsIndirect(const API::MatrixWorkspace &ws, const double minE, const double maxE) const;
};
}
}
diff --git a/Framework/Algorithms/inc/MantidAlgorithms/SofQW.h b/Framework/Algorithms/inc/MantidAlgorithms/SofQW.h
index 2326fc7fb2b37ce8293f76ca4e4bf3fde87f5580..d795b9223dea9e3eb16245bf94d6d863fb41f554 100644
--- a/Framework/Algorithms/inc/MantidAlgorithms/SofQW.h
+++ b/Framework/Algorithms/inc/MantidAlgorithms/SofQW.h
@@ -5,7 +5,6 @@
// Includes
//----------------------------------------------------------------------
#include "MantidAPI/DataProcessorAlgorithm.h"
-#include "MantidAlgorithms/SofQCommon.h"
namespace Mantid {
namespace Algorithms {
@@ -48,6 +47,9 @@ 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 SofQCommon;
+
class DLLExport SofQW : public API::DataProcessorAlgorithm {
public:
/// Algorithm's name
@@ -64,19 +66,16 @@ public:
setUpOutputWorkspace(const API::MatrixWorkspace_const_sptr &inputWorkspace,
const std::vector<double> &qbinParams,
std::vector<double> &qAxis,
- const std::vector<double> &ebinParams);
+ const std::vector<double> &ebinParams,
+ const SofQCommon &emodeProperties);
/// Create the input properties on the given algorithm object
static void createCommonInputProperties(API::Algorithm &alg);
- /// Energy to K constant
- static double energyToK();
private:
/// Initialization code
void init() override;
/// Execution code
void exec() override;
-
- SofQCommon m_EmodeProperties;
};
} // namespace Algorithms
diff --git a/Framework/Algorithms/src/SofQCommon.cpp b/Framework/Algorithms/src/SofQCommon.cpp
index e1fbb0753fc5cfcd9f13ee09185bbcf3af73fd2a..9dd16ab244cebdeb64843905b16f7edf16ef0b54 100644
--- a/Framework/Algorithms/src/SofQCommon.cpp
+++ b/Framework/Algorithms/src/SofQCommon.cpp
@@ -114,42 +114,76 @@ std::pair<double, double> SofQCommon::eBinHints(const API::MatrixWorkspace &ws)
* @return a pair containing global minimun and maximum Q
*/
std::pair<double, double> SofQCommon::qBinHints(const API::MatrixWorkspace &ws, const double minE, const double maxE) const {
- using namespace Mantid::PhysicalConstants;
if (m_emode == 1) {
- // Direct geometry
- auto minTheta = std::numeric_limits<double>::max();
- auto maxTheta = std::numeric_limits<double>::lowest();
- const auto &spectrumInfo = ws.spectrumInfo();
- for (size_t i = 0; i < spectrumInfo.size(); ++i) {
- const auto theta = spectrumInfo.twoTheta(i);
- if (theta < minTheta) {
- minTheta = theta;
- }
- if (theta > maxTheta) {
- maxTheta = theta;
- }
+ return qBinHintsDirect(ws, minE, maxE);
+ }
+ return qBinHintsIndirect(ws, minE, maxE);
+}
+
+/**
+ * Return a pair of (minimum Q, maximum Q) for given
+ * direct geometry workspace.
+ * @param ws a workspace
+ * @param minE minimum energy transfer in ws
+ * @param maxE maximum energy transfer in ws
+ * @return a pair containing global minimun and maximum Q
+ */
+std::pair<double, double> SofQCommon::qBinHintsDirect(const API::MatrixWorkspace &ws, const double minE, const double maxE) const {
+ using namespace Mantid::PhysicalConstants;
+ auto minTheta = std::numeric_limits<double>::max();
+ auto maxTheta = std::numeric_limits<double>::lowest();
+ const auto &spectrumInfo = ws.spectrumInfo();
+ for (size_t i = 0; i < spectrumInfo.size(); ++i) {
+ if (spectrumInfo.isMasked(i) || spectrumInfo.isMonitor(i)) {
+ continue;
}
- const auto incidentKSq = m_efixed / E_mev_toNeutronWavenumberSq;
- const auto incidentK = std::sqrt(incidentKSq);
- const auto minEnergy = m_efixed - maxE;
- const auto maxEnergy = m_efixed - minE;
- const auto minKSq = minEnergy / E_mev_toNeutronWavenumberSq;
- const auto minK = std::sqrt(minKSq);
- const auto maxKSq = maxEnergy / E_mev_toNeutronWavenumberSq;
- const auto maxK = std::sqrt(maxKSq);
- std::array<double, 4> q;
- q[0] = std::sqrt(incidentKSq + minKSq - 2. * incidentK * minK * std::cos(minTheta));
- q[1] = std::sqrt(incidentKSq + minKSq - 2. * incidentK * minK * std::cos(maxTheta));
- q[2] = std::sqrt(incidentKSq + maxKSq - 2. * incidentK * maxK * std::cos(minTheta));
- q[4] = std::sqrt(incidentKSq + maxKSq - 2. * incidentK * maxK * std::cos(maxTheta));
- const auto minmaxQ = std::minmax_element(q.cbegin(), q.cend());
- return std::make_pair(*minmaxQ.first, *minmaxQ.second);
+ const auto theta = spectrumInfo.twoTheta(i);
+ if (theta < minTheta) {
+ minTheta = theta;
+ }
+ if (theta > maxTheta) {
+ maxTheta = theta;
+ }
+ }
+ if (minTheta == std::numeric_limits<double>::max()) {
+ throw std::runtime_error("Could not determine Q binning: workspace does not contain usable spectra.");
}
- // Indirect geometry
+ const auto incidentKSq = m_efixed / E_mev_toNeutronWavenumberSq;
+ const auto incidentK = std::sqrt(incidentKSq);
+ const auto minEnergy = m_efixed - maxE;
+ const auto maxEnergy = m_efixed - minE;
+ const auto minKSq = minEnergy / E_mev_toNeutronWavenumberSq;
+ const auto minK = std::sqrt(minKSq);
+ const auto maxKSq = maxEnergy / E_mev_toNeutronWavenumberSq;
+ const auto maxK = std::sqrt(maxKSq);
+ std::array<double, 4> q;
+ q[0] = std::sqrt(incidentKSq + minKSq - 2. * incidentK * minK * std::cos(minTheta));
+ q[1] = std::sqrt(incidentKSq + minKSq - 2. * incidentK * minK * std::cos(maxTheta));
+ q[2] = std::sqrt(incidentKSq + maxKSq - 2. * incidentK * maxK * std::cos(minTheta));
+ q[4] = std::sqrt(incidentKSq + maxKSq - 2. * incidentK * maxK * std::cos(maxTheta));
+ const auto minmaxQ = std::minmax_element(q.cbegin(), q.cend());
+ return std::make_pair(*minmaxQ.first, *minmaxQ.second);
+}
+
+/**
+ * Return a pair of (minimum Q, maximum Q) for given
+ * indirect geometry workspace. Estimates the Q range from all detectors.
+ * If workspace contains grouped detectors/not all detectors are linked
+ * to a spectrum, the returned interval may be larger than actually needed.
+ * @param ws a workspace
+ * @param minE minimum energy transfer in ws
+ * @param maxE maximum energy transfer in ws
+ * @return a pair containing global minimun and maximum Q
+ */
+std::pair<double, double> SofQCommon::qBinHintsIndirect(const API::MatrixWorkspace &ws, const double minE, const double maxE) const {
+ using namespace Mantid::PhysicalConstants;
auto minQSq = std::numeric_limits<double>::max();
auto maxQSq = std::numeric_limits<double>::lowest();
const auto &detectorInfo = ws.detectorInfo();
for (size_t i = 0; i < detectorInfo.size(); ++i) {
+ if (detectorInfo.isMasked(i) || detectorInfo.isMonitor(i)) {
+ continue;
+ }
const auto costheta = std::cos(detectorInfo.twoTheta(i));
const auto eFixed = getEFixed(detectorInfo.detector(i));
const auto minEnergy = eFixed + minE;
@@ -170,6 +204,9 @@ std::pair<double, double> SofQCommon::qBinHints(const API::MatrixWorkspace &ws,
maxQSq = minmaxQSq.second;
}
}
+ if (minQSq == std::numeric_limits<double>::max()) {
+ throw std::runtime_error("Could not determine Q binning: workspace does not contain usable spectra.");
+ }
return std::make_pair(std::sqrt(minQSq), std::sqrt(maxQSq));
}
diff --git a/Framework/Algorithms/src/SofQW.cpp b/Framework/Algorithms/src/SofQW.cpp
index 4fecf133aee7ee8f0fe6b55405cc973bfd428f6f..af5a52fc79fb612dde76730826c2653223d67663 100644
--- a/Framework/Algorithms/src/SofQW.cpp
+++ b/Framework/Algorithms/src/SofQW.cpp
@@ -8,6 +8,7 @@
#include "MantidAPI/SpectrumDetectorMapping.h"
#include "MantidAPI/WorkspaceFactory.h"
#include "MantidAPI/WorkspaceUnitValidator.h"
+#include "MantidAlgorithms/SofQCommon.h"
#include "MantidAlgorithms/SofQW.h"
#include "MantidDataObjects/Histogram1D.h"
#include "MantidGeometry/Instrument/DetectorGroup.h"
@@ -15,7 +16,6 @@
#include "MantidKernel/BoundedValidator.h"
#include "MantidKernel/CompositeValidator.h"
#include "MantidKernel/ListValidator.h"
-#include "MantidKernel/PhysicalConstants.h"
#include "MantidKernel/RebinParamsValidator.h"
#include "MantidKernel/UnitFactory.h"
#include "MantidKernel/VectorHelper.h"
@@ -26,15 +26,6 @@ namespace Algorithms {
// Register the algorithm into the AlgorithmFactory
DECLARE_ALGORITHM(SofQW)
-/// Energy to K constant
-double SofQW::energyToK() {
- static const double energyToK = 8.0 * M_PI * M_PI *
- PhysicalConstants::NeutronMass *
- PhysicalConstants::meV * 1e-20 /
- (PhysicalConstants::h * PhysicalConstants::h);
- return energyToK;
-}
-
using namespace Kernel;
using namespace API;
@@ -93,10 +84,9 @@ void SofQW::createCommonInputProperties(API::Algorithm &alg) {
make_unique<ArrayProperty<double>>(
"QAxisBinning", boost::make_shared<RebinParamsValidator>()),
"The bin parameters to use for the q axis (in the format used by the "
- ":ref:`algm-Rebin` algorithm, except that bin width-only is "
- "unsupported).");
+ ":ref:`algm-Rebin` algorithm).");
- std::vector<std::string> propOptions{"Direct", "Indirect"};
+ const std::vector<std::string> propOptions{"Direct", "Indirect"};
alg.declareProperty("EMode", "",
boost::make_shared<StringListValidator>(propOptions),
"The energy transfer analysis mode (Direct/Indirect)");
@@ -114,8 +104,7 @@ void SofQW::createCommonInputProperties(API::Algorithm &alg) {
make_unique<ArrayProperty<double>>(
"EAxisBinning", boost::make_shared<RebinParamsValidator>(true)),
"The bin parameters to use for the E axis (optional, in the format "
- "used by the :ref:`algm-Rebin` algorithm, except that bin width-only is "
- "unsupported).");
+ "used by the :ref:`algm-Rebin` algorithm).");
}
void SofQW::exec() {
@@ -153,21 +142,35 @@ void SofQW::exec() {
API::MatrixWorkspace_sptr SofQW::setUpOutputWorkspace(
const API::MatrixWorkspace_const_sptr &inputWorkspace,
const std::vector<double> &qbinParams, std::vector<double> &qAxis,
- const std::vector<double> &ebinParams) {
+ const std::vector<double> &ebinParams,
+ const SofQCommon &emodeProperties) {
+ using Kernel::VectorHelper::createAxisFromRebinParams;
// Create vector to hold the new X axis values
HistogramData::BinEdges xAxis(0);
+ auto eHints = std::make_pair<double, double>(std::nan(""), std::nan(""));
int xLength;
if (ebinParams.empty()) {
xAxis = inputWorkspace->refX(0);
xLength = static_cast<int>(xAxis.size());
+ } else if (ebinParams.size() == 1) {
+ eHints = emodeProperties.eBinHints(*inputWorkspace);
+ xLength = createAxisFromRebinParams(ebinParams, xAxis.mutableRawData(), true, true, eHints.first, eHints.second);
} else {
- xLength = static_cast<int>(VectorHelper::createAxisFromRebinParams(
- ebinParams, xAxis.mutableRawData()));
+ xLength = createAxisFromRebinParams(
+ ebinParams, xAxis.mutableRawData());
}
// Create a vector to temporarily hold the vertical ('y') axis and populate
// that
- const int yLength = static_cast<int>(
- VectorHelper::createAxisFromRebinParams(qbinParams, qAxis));
+ int yLength;
+ if (qbinParams.size() == 1) {
+ if (std::isnan(eHints.first)) {
+ eHints = emodeProperties.eBinHints(*inputWorkspace);
+ }
+ const auto qHints = emodeProperties.qBinHints(*inputWorkspace, eHints.first, eHints.second);
+ yLength = createAxisFromRebinParams(qbinParams, qAxis, true, true, qHints.first, qHints.second);
+ } else {
+ yLength = createAxisFromRebinParams(qbinParams, qAxis);
+ }
// Create the output workspace
MatrixWorkspace_sptr outputWorkspace = WorkspaceFactory::Instance().create(
@@ -194,5 +197,6 @@ API::MatrixWorkspace_sptr SofQW::setUpOutputWorkspace(
return outputWorkspace;
}
+
} // namespace Algorithms
} // namespace Mantid
diff --git a/Framework/Algorithms/src/SofQWCentre.cpp b/Framework/Algorithms/src/SofQWCentre.cpp
index 089e8911eb2d0d1e6b1eb82ee671cc40ddf4b0a6..9e1ccbf9fbae90ee5e880ea35ababd05fee15744 100644
--- a/Framework/Algorithms/src/SofQWCentre.cpp
+++ b/Framework/Algorithms/src/SofQWCentre.cpp
@@ -61,10 +61,13 @@ void SofQWCentre::exec() {
"The input workspace must have common binning across all spectra");
}
+ m_EmodeProperties.initCachedValues(*inputWorkspace, this);
+ const int emode = m_EmodeProperties.m_emode;
+
std::vector<double> verticalAxis;
MatrixWorkspace_sptr outputWorkspace =
SofQW::setUpOutputWorkspace(inputWorkspace, getProperty("QAxisBinning"),
- verticalAxis, getProperty("EAxisBinning"));
+ verticalAxis, getProperty("EAxisBinning"), m_EmodeProperties);
setProperty("OutputWorkspace", outputWorkspace);
const auto &xAxis = outputWorkspace->binEdges(0).rawData();
@@ -72,8 +75,6 @@ void SofQWCentre::exec() {
std::vector<specnum_t> specNumberMapping;
std::vector<detid_t> detIDMapping;
- m_EmodeProperties.initCachedValues(*inputWorkspace, this);
- int emode = m_EmodeProperties.m_emode;
const auto &detectorInfo = inputWorkspace->detectorInfo();
const auto &spectrumInfo = inputWorkspace->spectrumInfo();
diff --git a/Framework/Algorithms/src/SofQWNormalisedPolygon.cpp b/Framework/Algorithms/src/SofQWNormalisedPolygon.cpp
index e5833ee1064a41ecf52d2a2c339723849501f863..162ab0188da371daf8d28963cda68ea1f1c504be 100644
--- a/Framework/Algorithms/src/SofQWNormalisedPolygon.cpp
+++ b/Framework/Algorithms/src/SofQWNormalisedPolygon.cpp
@@ -13,6 +13,7 @@
#include "MantidGeometry/Objects/BoundingBox.h"
#include "MantidGeometry/Objects/IObject.h"
#include "MantidIndexing/IndexInfo.h"
+#include "MantidKernel/PhysicalConstants.h"
#include "MantidKernel/UnitFactory.h"
#include "MantidKernel/VectorHelper.h"
#include "MantidTypes/SpectrumDefinition.h"
@@ -72,6 +73,9 @@ void SofQWNormalisedPolygon::exec() {
"The input workspace must have common binning across all spectra");
}
+ // Compute input caches
+ m_EmodeProperties.initCachedValues(*inputWS, this);
+
RebinnedOutput_sptr outputWS =
this->setUpOutputWorkspace(*inputWS, getProperty("QAxisBinning"), m_Qout,
getProperty("EAxisBinning"));
@@ -88,9 +92,6 @@ void SofQWNormalisedPolygon::exec() {
m_progress = boost::shared_ptr<API::Progress>(
new API::Progress(this, 0.0, 1.0, nreports));
- // Compute input caches
- m_EmodeProperties.initCachedValues(*inputWS, this);
-
std::vector<double> par =
inputWS->getInstrument()->getNumberParameter("detector-neighbour-offset");
if (par.empty()) {
@@ -223,14 +224,15 @@ double SofQWNormalisedPolygon::calculateQ(const double efixed, int emode,
const double deltaE,
const double twoTheta,
const double azimuthal) const {
+ using Mantid::PhysicalConstants::E_mev_toNeutronWavenumberSq;
double ki = 0.0;
double kf = 0.0;
if (emode == 1) {
- ki = std::sqrt(efixed * SofQW::energyToK());
- kf = std::sqrt((efixed - deltaE) * SofQW::energyToK());
+ ki = std::sqrt(efixed / E_mev_toNeutronWavenumberSq);
+ kf = std::sqrt((efixed - deltaE) / E_mev_toNeutronWavenumberSq);
} else if (emode == 2) {
- ki = std::sqrt((deltaE + efixed) * SofQW::energyToK());
- kf = std::sqrt(efixed * SofQW::energyToK());
+ ki = std::sqrt((deltaE + efixed) / E_mev_toNeutronWavenumberSq);
+ kf = std::sqrt(efixed / E_mev_toNeutronWavenumberSq);
}
const double Qx = ki - kf * std::cos(twoTheta);
const double Qy = -kf * std::sin(twoTheta) * std::cos(azimuthal);
@@ -411,18 +413,29 @@ RebinnedOutput_sptr SofQWNormalisedPolygon::setUpOutputWorkspace(
using Kernel::VectorHelper::createAxisFromRebinParams;
HistogramData::BinEdges xAxis(0);
+ auto eHints = std::make_pair<double, double>(std::nan(""), std::nan(""));
// Create vector to hold the new X axis values
if (ebinParams.empty()) {
xAxis = inputWorkspace.binEdges(0);
+ } else if (ebinParams.size() == 1) {
+ eHints = m_EmodeProperties.eBinHints(inputWorkspace);
+ createAxisFromRebinParams(ebinParams, xAxis.mutableRawData(), true, true, eHints.first, eHints.second);
} else {
- static_cast<void>(
- createAxisFromRebinParams(ebinParams, xAxis.mutableRawData()));
+ createAxisFromRebinParams(ebinParams, xAxis.mutableRawData());
}
// Create a vector to temporarily hold the vertical ('y') axis and populate
// that
- const int yLength = static_cast<int>(
- VectorHelper::createAxisFromRebinParams(qbinParams, qAxis));
+ int yLength;
+ if (qbinParams.size() == 1) {
+ if (std::isnan(eHints.first)) {
+ eHints = m_EmodeProperties.eBinHints(inputWorkspace);
+ }
+ const auto qHints = m_EmodeProperties.qBinHints(inputWorkspace, eHints.first, eHints.second);
+ yLength = createAxisFromRebinParams(qbinParams, qAxis, true, true, qHints.first, qHints.second);
+ } else {
+ yLength = createAxisFromRebinParams(qbinParams, qAxis);
+ }
// Create output workspace, bin edges are same as in inputWorkspace index 0
auto outputWorkspace =
diff --git a/Framework/Algorithms/src/SofQWPolygon.cpp b/Framework/Algorithms/src/SofQWPolygon.cpp
index 8fe9b4c015a2c7f98485cd0adb88583ae22dbc39..10c0fcdeedc0e0c4bb11bf36ef851763615ee2c5 100644
--- a/Framework/Algorithms/src/SofQWPolygon.cpp
+++ b/Framework/Algorithms/src/SofQWPolygon.cpp
@@ -4,10 +4,11 @@
#include "MantidAPI/SpectraAxis.h"
#include "MantidAPI/SpectrumDetectorMapping.h"
#include "MantidAPI/SpectrumInfo.h"
+#include "MantidDataObjects/FractionalRebinning.h"
#include "MantidGeometry/Math/PolygonIntersection.h"
#include "MantidGeometry/Math/Quadrilateral.h"
#include "MantidGeometry/Instrument/DetectorGroup.h"
-#include "MantidDataObjects/FractionalRebinning.h"
+#include "MantidKernel/PhysicalConstants.h"
namespace Mantid {
namespace Algorithms {
@@ -39,9 +40,11 @@ void SofQWPolygon::exec() {
"The input workspace must have common binning across all spectra");
}
+ m_EmodeProperties.initCachedValues(*inputWS, this);
+
MatrixWorkspace_sptr outputWS =
SofQW::setUpOutputWorkspace(inputWS, getProperty("QAxisBinning"), m_Qout,
- getProperty("EAxisBinning"));
+ getProperty("EAxisBinning"), m_EmodeProperties);
setProperty("OutputWorkspace", outputWS);
const size_t nenergyBins = inputWS->blocksize();
@@ -159,8 +162,9 @@ void SofQWPolygon::exec() {
double SofQWPolygon::calculateDirectQ(const double efixed, const double deltaE,
const double twoTheta,
const double psi) const {
- const double ki = std::sqrt(efixed * SofQW::energyToK());
- const double kf = std::sqrt((efixed - deltaE) * SofQW::energyToK());
+ using Mantid::PhysicalConstants::E_mev_toNeutronWavenumberSq;
+ const double ki = std::sqrt(efixed / E_mev_toNeutronWavenumberSq);
+ const double kf = std::sqrt((efixed - deltaE) / E_mev_toNeutronWavenumberSq);
const double Qx = ki - kf * std::cos(twoTheta);
const double Qy = -kf * std::sin(twoTheta) * std::cos(psi);
const double Qz = -kf * std::sin(twoTheta) * std::sin(psi);
@@ -168,7 +172,7 @@ double SofQWPolygon::calculateDirectQ(const double efixed, const double deltaE,
}
/**
- * Calculate the Q value for a direct instrument
+ * Calculate the Q value for an indirect instrument
* @param efixed An efixed value
* @param deltaE The energy change
* @param twoTheta The value of the scattering angle
@@ -180,8 +184,9 @@ double SofQWPolygon::calculateIndirectQ(const double efixed,
const double twoTheta,
const double psi) const {
UNUSED_ARG(psi);
- const double ki = std::sqrt((efixed + deltaE) * SofQW::energyToK());
- const double kf = std::sqrt(efixed * SofQW::energyToK());
+ using Mantid::PhysicalConstants::E_mev_toNeutronWavenumberSq;
+ const double ki = std::sqrt((efixed + deltaE) / E_mev_toNeutronWavenumberSq);
+ const double kf = std::sqrt(efixed / E_mev_toNeutronWavenumberSq);
const double Qx = ki - kf * std::cos(twoTheta);
const double Qy = -kf * std::sin(twoTheta);
return std::sqrt(Qx * Qx + Qy * Qy);
diff --git a/Framework/Algorithms/test/SofQWNormalisedPolygonTest.h b/Framework/Algorithms/test/SofQWNormalisedPolygonTest.h
index 26588e4816b39cd7e799f6d47617fdc5da4630ff..fa43c0157944a8c9cc1b0c8e5ca4c2a2bc56bd68 100644
--- a/Framework/Algorithms/test/SofQWNormalisedPolygonTest.h
+++ b/Framework/Algorithms/test/SofQWNormalisedPolygonTest.h
@@ -85,6 +85,111 @@ public:
TS_ASSERT_EQUALS(expectedIDs[i], spectrum.getDetectorIDs());
}
}
+
+ void testEAndQBinningParams() {
+ // SofQWNormalisedPolygon uses it's own setUpOutputWorkspace while
+ // the other SofQW* algorithms use the one in SofQW.
+ auto inWS = SofQWTest::loadTestFile();
+ Mantid::Algorithms::SofQWNormalisedPolygon alg;
+ alg.initialize();
+ alg.setChild(true);
+ alg.setRethrows(true);
+ TS_ASSERT_THROWS_NOTHING(alg.setProperty("InputWorkspace", inWS))
+ TS_ASSERT_THROWS_NOTHING(alg.setProperty("OutputWorkspace", "_unused"))
+ TS_ASSERT_THROWS_NOTHING(alg.setProperty("EMode", "Indirect"))
+ TS_ASSERT_THROWS_NOTHING(alg.setProperty("EFixed", 1.84))
+ const std::vector<double> eBinParams{-0.5, 0.1, -0.1, 0.2, 0.4};
+ const std::vector<double> expectedEBinEdges{-0.5, -0.4, -0.3, -0.2, -0.1, 0.1, 0.3, 0.4};
+ TS_ASSERT_THROWS_NOTHING(alg.setProperty("EAxisBinning", eBinParams))
+ const std::vector<double> qBinParams{0.5, 0.1, 1.0, 0.2, 2.};
+ const std::vector<double> expectedQBinEdges{0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.2, 1.4, 1.6, 1.8, 2.};
+ TS_ASSERT_THROWS_NOTHING(alg.setProperty("QAxisBinning", qBinParams))
+ TS_ASSERT_THROWS_NOTHING(alg.execute())
+ TS_ASSERT(alg.isExecuted())
+ Mantid::API::MatrixWorkspace_sptr outWS = alg.getProperty("OutputWorkspace");
+ TS_ASSERT_EQUALS(outWS->getNumberHistograms(), expectedQBinEdges.size() - 1)
+ for (size_t i = 0; i < outWS->getNumberHistograms(); ++i) {
+ const auto &x = outWS->x(i);
+ for (size_t j = 0; j < x.size(); ++j) {
+ TS_ASSERT_DELTA(x[j], expectedEBinEdges[j], 1e-12)
+ }
+ }
+ const auto axis = outWS->getAxis(1);
+ TS_ASSERT_EQUALS(axis->length(), expectedQBinEdges.size())
+ for (size_t i = 0; i < axis->length(); ++i) {
+ TS_ASSERT_DELTA(axis->getValue(i), expectedQBinEdges[i], 1e-12)
+ }
+ }
+
+ void testEBinWidthAsEAxisBinning() {
+ // SofQWNormalisedPolygon uses it's own setUpOutputWorkspace while
+ // the other SofQW* algorithms use the one in SofQW.
+ auto inWS = SofQWTest::loadTestFile();
+ Mantid::Algorithms::SofQWNormalisedPolygon alg;
+ alg.initialize();
+ alg.setChild(true);
+ alg.setRethrows(true);
+ TS_ASSERT_THROWS_NOTHING(alg.setProperty("InputWorkspace", inWS))
+ TS_ASSERT_THROWS_NOTHING(alg.setProperty("OutputWorkspace", "_unused"))
+ TS_ASSERT_THROWS_NOTHING(alg.setProperty("EMode", "Indirect"))
+ TS_ASSERT_THROWS_NOTHING(alg.setProperty("EFixed", 1.84))
+ const double dE{0.3};
+ const std::vector<double> eBinParams{dE};
+ std::vector<double> expectedEBinEdges;
+ const auto firstEdge = inWS->x(0).front();
+ const auto lastEdge = inWS->x(0).back();
+ auto currentEdge = firstEdge;
+ while (currentEdge < lastEdge) {
+ expectedEBinEdges.emplace_back(currentEdge);
+ currentEdge += dE;
+ }
+ expectedEBinEdges.emplace_back(lastEdge);
+ TS_ASSERT_THROWS_NOTHING(alg.setProperty("EAxisBinning", eBinParams))
+ const std::vector<double> qBinParams{0.5, 0.1, 1.0, 0.2, 2.};
+ const std::vector<double> expectedQBinEdges{0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.2, 1.4, 1.6, 1.8, 2.};
+ TS_ASSERT_THROWS_NOTHING(alg.setProperty("QAxisBinning", qBinParams))
+ TS_ASSERT_THROWS_NOTHING(alg.execute())
+ TS_ASSERT(alg.isExecuted())
+ Mantid::API::MatrixWorkspace_sptr outWS = alg.getProperty("OutputWorkspace");
+ TS_ASSERT_EQUALS(outWS->getNumberHistograms(), expectedQBinEdges.size() - 1)
+ for (size_t i = 0; i < outWS->getNumberHistograms(); ++i) {
+ const auto &x = outWS->x(i);
+ for (size_t j = 0; j < x.size(); ++j) {
+ TS_ASSERT_DELTA(x[j], expectedEBinEdges[j], 1e-12)
+ }
+ }
+ const auto axis = outWS->getAxis(1);
+ TS_ASSERT_EQUALS(axis->length(), expectedQBinEdges.size())
+ for (size_t i = 0; i < axis->length(); ++i) {
+ TS_ASSERT_DELTA(axis->getValue(i), expectedQBinEdges[i], 1e-12)
+ }
+ }
+
+ void testQBinWidthAsQAxisBinning() {
+ // SofQWNormalisedPolygon uses it's own setUpOutputWorkspace while
+ // the other SofQW* algorithms use the one in SofQW.
+ auto inWS = SofQWTest::loadTestFile();
+ Mantid::Algorithms::SofQWNormalisedPolygon alg;
+ alg.initialize();
+ alg.setChild(true);
+ alg.setRethrows(true);
+ TS_ASSERT_THROWS_NOTHING(alg.setProperty("InputWorkspace", inWS))
+ TS_ASSERT_THROWS_NOTHING(alg.setProperty("OutputWorkspace", "_unused"))
+ TS_ASSERT_THROWS_NOTHING(alg.setProperty("EMode", "Indirect"))
+ TS_ASSERT_THROWS_NOTHING(alg.setProperty("EFixed", 1.84))
+ const double dQ{0.023};
+ const std::vector<double> qBinParams{dQ};
+ TS_ASSERT_THROWS_NOTHING(alg.setProperty("QAxisBinning", qBinParams))
+ TS_ASSERT_THROWS_NOTHING(alg.execute())
+ TS_ASSERT(alg.isExecuted())
+ Mantid::API::MatrixWorkspace_sptr outWS = alg.getProperty("OutputWorkspace");
+ const auto axis = outWS->getAxis(1);
+ // Test only the Q bin width, not the actual edges.
+ for (size_t i = 0; i < axis->length() - 1; ++i) {
+ const auto delta = axis->getValue(i + 1) - axis->getValue(i);
+ TS_ASSERT_DELTA(delta, dQ, 1e-12);
+ }
+ }
};
class SofQWNormalisedPolygonTestPerformance : public CxxTest::TestSuite {
diff --git a/Framework/Algorithms/test/SofQWTest.h b/Framework/Algorithms/test/SofQWTest.h
index 695e9b5dbb6361df26bfe01a4be8a52f0a32528a..da4b304a6bff0b0c36cee6755270cad5a0195ece 100644
--- a/Framework/Algorithms/test/SofQWTest.h
+++ b/Framework/Algorithms/test/SofQWTest.h
@@ -3,6 +3,7 @@
#include <cxxtest/TestSuite.h>
#include "MantidAlgorithms/SofQW.h"
+#include "MantidAlgorithms/SofQCommon.h"
#include "MantidAPI/AnalysisDataService.h"
#include "MantidAPI/Axis.h"
#include "MantidAPI/MatrixWorkspace.h"
@@ -16,9 +17,8 @@ using namespace Mantid::API;
class SofQWTest : public CxxTest::TestSuite {
public:
- template <typename SQWType>
- static Mantid::API::MatrixWorkspace_sptr
- runSQW(const std::string &method = "") {
+
+ static Mantid::API::MatrixWorkspace_sptr loadTestFile() {
Mantid::DataHandling::LoadNexusProcessed loader;
loader.initialize();
loader.setChild(true);
@@ -31,15 +31,22 @@ public:
auto inWS =
boost::dynamic_pointer_cast<Mantid::API::MatrixWorkspace>(loadedWS);
WorkspaceHelpers::makeDistribution(inWS);
+ return inWS;
+ }
+
+ template <typename SQWType>
+ static Mantid::API::MatrixWorkspace_sptr
+ runSQW(const std::string &method = "") {
+ auto inWS = loadTestFile();
SQWType sqw;
sqw.initialize();
+ sqw.setRethrows(true);
// Cannot be marked as child or history is not recorded
TS_ASSERT_THROWS_NOTHING(sqw.setProperty("InputWorkspace", inWS));
- std::ostringstream wsname;
- wsname << "_tmp_" << loadedWS;
+ const std::string wsname{"_tmp_"};
TS_ASSERT_THROWS_NOTHING(
- sqw.setPropertyValue("OutputWorkspace", wsname.str()));
+ sqw.setPropertyValue("OutputWorkspace", wsname));
TS_ASSERT_THROWS_NOTHING(
sqw.setPropertyValue("QAxisBinning", "0.5,0.25,2"));
TS_ASSERT_THROWS_NOTHING(sqw.setPropertyValue("EMode", "Indirect"));
@@ -52,8 +59,8 @@ public:
auto &dataStore = Mantid::API::AnalysisDataService::Instance();
auto result =
- dataStore.retrieveWS<Mantid::API::MatrixWorkspace>(wsname.str());
- dataStore.remove(wsname.str());
+ dataStore.retrieveWS<Mantid::API::MatrixWorkspace>(wsname);
+ dataStore.remove(wsname);
return result;
}
@@ -130,8 +137,111 @@ public:
TS_ASSERT(!nanFound);
}
+ void testSetUpOutputWorkspace() {
+ auto inWS = loadTestFile();
+ Mantid::Algorithms::SofQW alg;
+ alg.initialize();
+ alg.setProperty("EMode", "Indirect");
+ alg.setProperty("EFixed", 1.84);
+ Mantid::Algorithms::SofQCommon emodeProperties;
+ emodeProperties.initCachedValues(*inWS, &alg);
+ const std::vector<double> eBinParams{-0.5, 0.1, -0.1, 0.2, 0.4};
+ const std::vector<double> expectedEBinEdges{-0.5, -0.4, -0.3, -0.2, -0.1, 0.1, 0.3, 0.4};
+ const std::vector<double> qBinParams{0.5, 0.1, 1.0, 0.2, 2.};
+ const std::vector<double> expectedQBinEdges{0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.2, 1.4, 1.6, 1.8, 2.};
+ std::vector<double> qAxis;
+ auto outWS = Mantid::Algorithms::SofQW::setUpOutputWorkspace(inWS, qBinParams, qAxis, eBinParams, emodeProperties);
+ TS_ASSERT_EQUALS(outWS->getNumberHistograms(), expectedQBinEdges.size() - 1)
+ for (size_t i = 0; i < outWS->getNumberHistograms(); ++i) {
+ const auto &x = outWS->x(i);
+ for (size_t j = 0; j < x.size(); ++j) {
+ TS_ASSERT_DELTA(x[j], expectedEBinEdges[j], 1e-12)
+ }
+ }
+ const auto axis = outWS->getAxis(1);
+ TS_ASSERT_EQUALS(axis->length(), qAxis.size())
+ TS_ASSERT_EQUALS(qAxis.size(), expectedQBinEdges.size())
+ for (size_t i = 0; i < qAxis.size(); ++i) {
+ TS_ASSERT_EQUALS(axis->getValue(i), qAxis[i])
+ TS_ASSERT_DELTA(qAxis[i], expectedQBinEdges[i], 1e-12)
+ }
+ }
+
+ void testSetUpOutputWorkspaceWithEnergyBinWidthOnly() {
+ auto inWS = loadTestFile();
+
+ Mantid::Algorithms::SofQW alg;
+ alg.initialize();
+ alg.setProperty("EMode", "Indirect");
+ alg.setProperty("EFixed", 1.84);
+ Mantid::Algorithms::SofQCommon emodeProperties;
+ emodeProperties.initCachedValues(*inWS, &alg);
+ const double dE{0.3};
+ const std::vector<double> eBinParams{dE};
+ std::vector<double> expectedEBinEdges;
+ const auto firstEdge = inWS->x(0).front();
+ const auto lastEdge = inWS->x(0).back();
+ auto currentEdge = firstEdge;
+ while (currentEdge < lastEdge) {
+ expectedEBinEdges.emplace_back(currentEdge);
+ currentEdge += dE;
+ }
+ expectedEBinEdges.emplace_back(lastEdge);
+ const std::vector<double> qBinParams{0.5, 0.25, 2.};
+ const std::vector<double> expectedQBinEdges{0.5, 0.75, 1., 1.25, 1.5, 1.75, 2.};
+ std::vector<double> qAxis;
+ auto outWS = Mantid::Algorithms::SofQW::setUpOutputWorkspace(inWS, qBinParams, qAxis, eBinParams, emodeProperties);
+ TS_ASSERT_EQUALS(outWS->getNumberHistograms(), expectedQBinEdges.size() - 1)
+ for (size_t i = 0; i < outWS->getNumberHistograms(); ++i) {
+ const auto &x = outWS->x(i);
+ for (size_t j = 0; j < x.size(); ++j) {
+ TS_ASSERT_DELTA(x[j], expectedEBinEdges[j], 1e-12)
+ }
+ }
+ const auto axis = outWS->getAxis(1);
+ TS_ASSERT_EQUALS(axis->length(), qAxis.size())
+ TS_ASSERT_EQUALS(qAxis.size(), expectedQBinEdges.size())
+ for (size_t i = 0; i < qAxis.size(); ++i) {
+ TS_ASSERT_EQUALS(axis->getValue(i), qAxis[i])
+ TS_ASSERT_DELTA(qAxis[i], expectedQBinEdges[i], 1e-12)
+ }
+ }
+
+ void testSetUpOutputWorkspaceWithQBinWidthOnly() {
+ auto inWS = loadTestFile();
+
+ Mantid::Algorithms::SofQW alg;
+ alg.initialize();
+ alg.setProperty("EMode", "Indirect");
+ alg.setProperty("EFixed", 1.84);
+ Mantid::Algorithms::SofQCommon emodeProperties;
+ emodeProperties.initCachedValues(*inWS, &alg);
+ const std::vector<double> eBinParams{-0.3, 0.2, 0.5};
+ const std::vector<double> expectedEBinEdges{-0.3, -0.1, 0.1, 0.3, 0.5};
+ const double dQ{0.023};
+ const std::vector<double> qBinParams{dQ};
+ std::vector<double> qAxis;
+ auto outWS = Mantid::Algorithms::SofQW::setUpOutputWorkspace(inWS, qBinParams, qAxis, eBinParams, emodeProperties);
+ for (size_t i = 0; i < outWS->getNumberHistograms(); ++i) {
+ const auto &x = outWS->x(i);
+ for (size_t j = 0; j < x.size(); ++j) {
+ TS_ASSERT_DELTA(x[j], expectedEBinEdges[j], 1e-12)
+ }
+ }
+ const auto axis = outWS->getAxis(1);
+ TS_ASSERT_EQUALS(axis->length(), qAxis.size())
+ for (size_t i = 0; i < qAxis.size(); ++i) {
+ TS_ASSERT_EQUALS(axis->getValue(i), qAxis[i])
+ }
+ // Test only the Q bin width, not the actual edges.
+ for (size_t i = 0; i < qAxis.size() - 1; ++i) {
+ const auto delta = qAxis[i + 1] - qAxis[i];
+ TS_ASSERT_DELTA(delta, dQ, 1e-12);
+ }
+ }
+
private:
- bool isAlgorithmInHistory(const Mantid::API::MatrixWorkspace &result,
+ static bool isAlgorithmInHistory(const Mantid::API::MatrixWorkspace &result,
const std::string &name) {
// Loaded nexus file has 13 other entries
const auto &wsHistory = result.getHistory();
diff --git a/Framework/Kernel/inc/MantidKernel/VectorHelper.h b/Framework/Kernel/inc/MantidKernel/VectorHelper.h
index c917ab0f7cc756555f6727cffb3f1cc8a1cae37b..ecfe4b64a09b50e7626cb64ff182bfdff99c864c 100644
--- a/Framework/Kernel/inc/MantidKernel/VectorHelper.h
+++ b/Framework/Kernel/inc/MantidKernel/VectorHelper.h
@@ -44,7 +44,9 @@ 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 bool full_bins_only = false,
+ const double xMinHint = std::nan(""),
+ const double xMaxHint = std::nan(""));
void MANTID_KERNEL_DLL
rebin(const std::vector<double> &xold, const std::vector<double> &yold,
diff --git a/Framework/Kernel/src/VectorHelper.cpp b/Framework/Kernel/src/VectorHelper.cpp
index d44302a366f0a7ff69d92894a9c8f91c5185c3a2..1e2ce79c1c554c57ccf858970d7fd1e17475a5d4 100644
--- a/Framework/Kernel/src/VectorHelper.cpp
+++ b/Framework/Kernel/src/VectorHelper.cpp
@@ -18,23 +18,35 @@ namespace VectorHelper {
*,x_n-1,delta_n-1,x_n]
* @param[out] xnew The newly created axis resulting from the input params
* @param[in] resize_xnew If false then the xnew vector is NOT resized. Useful
- *if on the number of bins needs determining. (Default=True)
+ *if the number of bins needs determining. (Default=True)
* @param[in] full_bins_only If true, bins of the size less than the current
- *step are not included
+ *step are not included. (Default=True)
+ * @param[in] xMinHint x_1 if params contains only delta_1.
+ * @param[in] xMaxHint x_2 if params contains only delta_1.
* @return The number of bin boundaries in the new axis
**/
int DLLExport createAxisFromRebinParams(const std::vector<double> ¶ms,
std::vector<double> &xnew,
const bool resize_xnew,
- const bool full_bins_only) {
- if (params.size() == 1) {
- throw std::runtime_error("Rebinning parameters containing only the bin "
- "width are not supported.");
- }
+ const bool full_bins_only,
+ const double xMinHint,
+ const double xMaxHint) {
+ std::vector<double> tmp;
+ const std::vector<double> &fullParams = [¶ms, &tmp, xMinHint, xMaxHint]() {
+ if (params.size() == 1) {
+ if (std::isnan(xMinHint) || std::isnan(xMaxHint)) {
+ throw std::runtime_error("createAxisFromRebinParams: xMinHint and xMaxHint must be supplied if params contains only the bin width.");
+ }
+ tmp.resize(3);
+ tmp = {xMinHint, params.front(), xMaxHint};
+ return tmp;
+ }
+ return params;
+ }();
double xs;
int ibound(2), istep(1), inew(1);
int ibounds = static_cast<int>(
- params.size()); // highest index in params array containing a bin boundary
+ fullParams.size()); // highest index in params array containing a bin boundary
int isteps = ibounds - 1; // highest index in params array containing a step
xnew.clear();
@@ -49,16 +61,16 @@ int DLLExport createAxisFromRebinParams(const std::vector<double> ¶ms,
lastBinCoef = 1.0;
}
- double xcurr = params[0];
+ double xcurr = fullParams[0];
if (resize_xnew)
xnew.push_back(xcurr);
while ((ibound <= ibounds) && (istep <= isteps)) {
// if step is negative then it is logarithmic step
- if (params[istep] >= 0.0)
- xs = params[istep];
+ if (fullParams[istep] >= 0.0)
+ xs = fullParams[istep];
else
- xs = xcurr * fabs(params[istep]);
+ xs = xcurr * fabs(fullParams[istep]);
if (fabs(xs) == 0.0) {
// Someone gave a 0-sized step! What a dope.
@@ -66,7 +78,7 @@ int DLLExport createAxisFromRebinParams(const std::vector<double> ¶ms,
"Invalid binning step provided! Can't creating binning axis.");
}
- if ((xcurr + xs * (1.0 + lastBinCoef)) <= params[ibound]) {
+ if ((xcurr + xs * (1.0 + lastBinCoef)) <= fullParams[ibound]) {
// If we can still fit current bin _plus_ specified portion of a last bin,
// continue
xcurr += xs;
@@ -80,7 +92,7 @@ int DLLExport createAxisFromRebinParams(const std::vector<double> ¶ms,
else
// For non full_bins_only, finish by adding as mush as is left from the
// range
- xcurr = params[ibound];
+ xcurr = fullParams[ibound];
ibound += 2;
istep += 2;
@@ -88,14 +100,6 @@ int DLLExport createAxisFromRebinParams(const std::vector<double> ¶ms,
if (resize_xnew)
xnew.push_back(xcurr);
inew++;
-
- // if (xnew.size() > 10000000)
- // {
- // //Max out at 1 million bins
- // throw std::runtime_error("Over ten million binning steps created.
- // Exiting to avoid infinite loops.");
- // return inew;
- // }
}
return inew;
diff --git a/Framework/Kernel/test/VectorHelperTest.h b/Framework/Kernel/test/VectorHelperTest.h
index a09ca17d5d2749ff4d2d68e67ff05f914122fa62..1d0313ab9eef9f52a97b359899aeff2f097d2710 100644
--- a/Framework/Kernel/test/VectorHelperTest.h
+++ b/Framework/Kernel/test/VectorHelperTest.h
@@ -160,13 +160,21 @@ public:
TS_ASSERT_EQUALS(axis, expectedAxis);
}
- void test_CreateAxisFromRebinParams_ThrowsIfSingleParam() {
+ void test_CreateAxisFromRebinParams_ThrowsIfSingleParamNoHintsProvided() {
const std::vector<double> rbParams = {1.0};
std::vector<double> axis;
TS_ASSERT_THROWS(VectorHelper::createAxisFromRebinParams(rbParams, axis),
const std::runtime_error)
}
+ void test_CreateAxisFromRebinParams_xMinXMaxHints() {
+ const std::vector<double> rbParams = {1.0};
+ std::vector<double> axis;
+ TS_ASSERT_THROWS_NOTHING(VectorHelper::createAxisFromRebinParams(rbParams, axis, true, true, -5., 3.))
+ const std::vector<double> expectedAxis = {-5, -4, -3, -2, -1, 0, 1, 2, 3};
+ TS_ASSERT_EQUALS(axis, expectedAxis);
+ }
+
void test_ConvertToBinBoundary_EmptyInputVector() {
std::vector<double> bin_centers;
std::vector<double> bin_edges;