diff --git a/Code/Mantid/Framework/Algorithms/inc/MantidAlgorithms/Qhelper.h b/Code/Mantid/Framework/Algorithms/inc/MantidAlgorithms/Qhelper.h
index 18b68d48f4c210e5e10a506c11394fd19b62822e..1bfb19eec544fb057a2f9f28cd8ee77d8201bd85 100644
--- a/Code/Mantid/Framework/Algorithms/inc/MantidAlgorithms/Qhelper.h
+++ b/Code/Mantid/Framework/Algorithms/inc/MantidAlgorithms/Qhelper.h
@@ -37,6 +37,11 @@ namespace Algorithms {
*/
class Qhelper {
public:
+ void examineInput(API::MatrixWorkspace_const_sptr dataWS,
+ API::MatrixWorkspace_const_sptr binAdj,
+ API::MatrixWorkspace_const_sptr detectAdj,
+ API::MatrixWorkspace_const_sptr qResolution);
+
void examineInput(API::MatrixWorkspace_const_sptr dataWS,
API::MatrixWorkspace_const_sptr binAdj,
API::MatrixWorkspace_const_sptr detectAdj);
diff --git a/Code/Mantid/Framework/Algorithms/src/Q1D2.cpp b/Code/Mantid/Framework/Algorithms/src/Q1D2.cpp
index 4f8fa7643512d41ca56033adf00363d0c5c9d9af..5f1d01f6e86f0068135f906e1c7d33f12964eedf 100644
--- a/Code/Mantid/Framework/Algorithms/src/Q1D2.cpp
+++ b/Code/Mantid/Framework/Algorithms/src/Q1D2.cpp
@@ -86,6 +86,11 @@ void Q1D2::init() {
declareProperty(
"ExtraLength", 0.0, mustBePositive,
"Additional length for gravity correction.");
+
+ declareProperty(
+ new WorkspaceProperty<>("QResolution", "", Direction::Input,
+ PropertyMode::Optional, dataVal),
+ "Workspace to calculate the Q resolution.\n");
}
/**
@ throw invalid_argument if the workspaces are not mututially compatible
@@ -95,12 +100,14 @@ void Q1D2::exec() {
MatrixWorkspace_const_sptr waveAdj = getProperty("WavelengthAdj");
MatrixWorkspace_const_sptr pixelAdj = getProperty("PixelAdj");
MatrixWorkspace_const_sptr wavePixelAdj = getProperty("WavePixelAdj");
+ MatrixWorkspace_const_sptr qResolution = getProperty("QResolution");
+
const bool doGravity = getProperty("AccountForGravity");
m_doSolidAngle = getProperty("SolidAngleWeighting");
// throws if we don't have common binning or another incompatibility
Qhelper helper;
- helper.examineInput(m_dataWS, waveAdj, pixelAdj);
+ helper.examineInput(m_dataWS, waveAdj, pixelAdj, qResolution);
// FIXME: how to examine the wavePixelAdj?
g_log.debug() << "All input workspaces were found to be valid\n";
// normalization as a function of wavelength (i.e. centers of x-value bins)
@@ -122,10 +129,15 @@ void Q1D2::exec() {
MantidVec normSum(YOut.size(), 0.0);
// the error on the normalisation
MantidVec normError2(YOut.size(), 0.0);
+ // the averaged Q resolution
+ MantidVec &qResolutionOut = outputWS->dataDx(0);
const int numSpec = static_cast<int>(m_dataWS->getNumberHistograms());
Progress progress(this, 0.05, 1.0, numSpec + 1);
+ // Flag to decide if Q Resolution is to be used
+ auto useQResolution = qResolution ? true : false;
+
PARALLEL_FOR3(m_dataWS, outputWS, pixelAdj)
for (int i = 0; i < numSpec; ++i) {
PARALLEL_START_INTERUPT_REGION
@@ -160,6 +172,7 @@ void Q1D2::exec() {
continue;
}
+
const size_t numWavbins = m_dataWS->readY(i).size() - wavStart;
// make just one call to new to reduce CPU overhead on each thread, access
// to these
@@ -184,6 +197,12 @@ void Q1D2::exec() {
MantidVec::const_iterator YIn = m_dataWS->readY(i).begin() + wavStart;
MantidVec::const_iterator EIn = m_dataWS->readE(i).begin() + wavStart;
+ // Pointers to the QResolution data. Note that the xdata was initially the same, hence
+ // the same indexing applies to the y values of m_dataWS and qResolution
+ // If we want to use it set it to the correct value, else to YIN, although that does not matter, as
+ // we won't use it
+ MantidVec::const_iterator QResIn = useQResolution ? (qResolution->readY(i).begin() + wavStart) : YIn;
+
// when finding the output Q bin remember that the input Q bins (from the
// convert to wavelength) start high and reduce
MantidVec::const_iterator loc = QOut.end();
@@ -205,13 +224,20 @@ void Q1D2::exec() {
// at the end
EOutTo2[bin] += (*EIn) * (*EIn);
normError2[bin] += *normETo2s;
+ if (useQResolution) {
+ qResolutionOut[bin] += (*YIn) * (*QResIn);
+ }
}
}
+
+ // Increment the QResolution iterator
+ if (useQResolution) {
+ ++QResIn;
+ }
}
PARALLEL_CRITICAL(q1d_spectra_map) {
progress.report("Computing I(Q)");
-
// Add up the detector IDs in the output spectrum at workspace index 0
const ISpectrum *inSpec = m_dataWS->getSpectrum(i);
ISpectrum *outSpec = outputWS->getSpectrum(0);
@@ -222,12 +248,30 @@ void Q1D2::exec() {
}
PARALLEL_CHECK_INTERUPT_REGION
+ if (useQResolution) {
+ // The number of Q (x)_ values is N, while the number of DeltaQ values is N-1,
+ // Richard Heenan suggested to duplicate the last entry of DeltaQ
+ Mantid::MantidVec::const_iterator countsIterator = YOut.begin();
+ Mantid::MantidVec::iterator qResolutionIterator = qResolutionOut.begin();
+ for (;countsIterator!= YOut.end(); ++countsIterator, ++qResolutionIterator) {
+ // Divide by the counts of the Qbin, if the counts are 0, the the qresolution will also be 0
+ if ((*countsIterator) > 0.0) {
+ *qResolutionIterator = (*qResolutionIterator)/(*countsIterator);
+ }
+ }
+ // Now duplicate write the second to last element into the last element of the deltaQ vector
+ if (qResolutionOut.size() > 1) {
+ qResolutionOut.rbegin()[0] = qResolutionOut.rbegin()[1];
+ }
+ }
+
bool doOutputParts = getProperty("OutputParts");
if (doOutputParts) {
MatrixWorkspace_sptr ws_sumOfCounts =
WorkspaceFactory::Instance().create(outputWS);
ws_sumOfCounts->dataX(0) = outputWS->dataX(0);
ws_sumOfCounts->dataY(0) = outputWS->dataY(0);
+ ws_sumOfCounts->dataDx(0) = outputWS->dataDx(0);
for (size_t i = 0; i < outputWS->dataE(0).size(); i++) {
ws_sumOfCounts->dataE(0)[i] = sqrt(outputWS->dataE(0)[i]);
}
@@ -235,6 +279,7 @@ void Q1D2::exec() {
MatrixWorkspace_sptr ws_sumOfNormFactors =
WorkspaceFactory::Instance().create(outputWS);
ws_sumOfNormFactors->dataX(0) = outputWS->dataX(0);
+ ws_sumOfNormFactors->dataDx(0) = outputWS->dataDx(0);
for (size_t i = 0; i < ws_sumOfNormFactors->dataY(0).size(); i++) {
ws_sumOfNormFactors->dataY(0)[i] = normSum[i];
ws_sumOfNormFactors->dataE(0)[i] = sqrt(normError2[i]);
diff --git a/Code/Mantid/Framework/Algorithms/src/Qhelper.cpp b/Code/Mantid/Framework/Algorithms/src/Qhelper.cpp
index 9dff883a551de317696752947e936c7c1b2bd9ca..07df46969aec725828b2fa05f682d488067d9728 100644
--- a/Code/Mantid/Framework/Algorithms/src/Qhelper.cpp
+++ b/Code/Mantid/Framework/Algorithms/src/Qhelper.cpp
@@ -17,6 +17,46 @@ using namespace Kernel;
using namespace API;
using namespace Geometry;
+/** Checks if workspaces input to Q1D or Qxy are reasonable
+ @param dataWS data workspace
+ @param binAdj (WavelengthAdj) workpace that will be checked to see if it has
+ one spectrum and the same number of bins as dataWS
+ @param detectAdj (PixelAdj) passing NULL for this wont raise an error, if set
+ it will be checked this workspace has as many histograms as dataWS each with
+ one bin
+ @param qResolution: the QResolution workspace
+ @throw invalid_argument if the workspaces are not mututially compatible
+*/
+void Qhelper::examineInput(API::MatrixWorkspace_const_sptr dataWS,
+ API::MatrixWorkspace_const_sptr binAdj,
+ API::MatrixWorkspace_const_sptr detectAdj,
+ API::MatrixWorkspace_const_sptr qResolution) {
+
+ // Check the compatibility of dataWS, binAdj and detectAdj
+ examineInput(dataWS, binAdj, detectAdj);
+
+ // Check the compatibility of the QResolution workspace
+ if (qResolution) {
+ // We require the same number of histograms
+ if (qResolution->getNumberHistograms() != dataWS->getNumberHistograms()) {
+ throw std::invalid_argument("The QResolution should have one spectrum"
+ "per spectrum of the input workspace");
+ }
+
+ // We require the same binning for the input workspace and the q resolution
+ // workspace
+ MantidVec::const_iterator reqX = dataWS->readX(0).begin();
+ MantidVec::const_iterator qResX = qResolution->readX(0).begin();
+ for (; reqX != dataWS->readX(0).end(); ++reqX, ++qResX) {
+ if (*reqX != *qResX) {
+ throw std::invalid_argument("The QResolution needs to have the same binning as"
+ "as the input workspace.");
+ }
+ }
+ }
+}
+
+
/** Checks if workspaces input to Q1D or Qxy are reasonable
@param dataWS data workspace
@param binAdj (WavelengthAdj) workpace that will be checked to see if it has
@@ -103,6 +143,7 @@ void Qhelper::examineInput(API::MatrixWorkspace_const_sptr dataWS,
}
}
+
/** Finds the first index number of the first wavelength bin that should
* included based on the
* the calculation: W = Wcut (Rcut-R)/Rcut
diff --git a/Code/Mantid/Framework/Algorithms/test/Q1D2Test.h b/Code/Mantid/Framework/Algorithms/test/Q1D2Test.h
index 1bcbb109d36f30405d56a7c5ee8ec4e9350e231c..c474d081aa436c9e7e58c8822a4241de3c8b0f13 100644
--- a/Code/Mantid/Framework/Algorithms/test/Q1D2Test.h
+++ b/Code/Mantid/Framework/Algorithms/test/Q1D2Test.h
@@ -13,6 +13,8 @@
#include "MantidTestHelpers/WorkspaceCreationHelper.h"
+#include <algorithm>
+
using namespace Mantid::API;
using namespace Mantid::Kernel;
using namespace Mantid::DataHandling;
@@ -24,6 +26,8 @@ static double flat_cell061Es [] = {8.140295E-03, 8.260089E-03, 8.198814E-03, 8.3
/// defined below this creates some input data
void createInputWorkspaces(int start, int end, Mantid::API::MatrixWorkspace_sptr & input, Mantid::API::MatrixWorkspace_sptr & wave, Mantid::API::MatrixWorkspace_sptr & pixels);
void createInputWorkSpacesForMasking ( Mantid::API::MatrixWorkspace_sptr & input, Mantid::API::MatrixWorkspace_sptr & wave, Mantid::API::MatrixWorkspace_sptr & pixels);
+void createQResolutionWorkspace(Mantid::API::MatrixWorkspace_sptr & qResolution, Mantid::API::MatrixWorkspace_sptr & alteredInput, Mantid::API::MatrixWorkspace_sptr & input, double value1, double value2);
+
class Q1D2Test : public CxxTest::TestSuite
{
@@ -341,7 +345,7 @@ public:
Mantid::API::AnalysisDataService::Instance().remove(m_noGrav);
}
-
+
void testInvalidInput()
{
Mantid::Algorithms::Q1D2 Q1D;
@@ -362,8 +366,94 @@ public:
Q1D.execute();
TS_ASSERT( ! Q1D.isExecuted() )
+
+ // Restore the old binning
+ xData[15] -= 0.001;
}
-
+
+ void testRunsWithQResolution() {
+ // Arrange
+ Mantid::API::MatrixWorkspace_sptr qResolution, alteredInput;
+ const double value1 = 1;
+ const double value2 = 2;
+ createQResolutionWorkspace(qResolution, alteredInput, m_inputWS, value1, value2);
+
+ // Act
+ Mantid::Algorithms::Q1D2 Q1D;
+ TS_ASSERT_THROWS_NOTHING( Q1D.initialize() );
+ TS_ASSERT( Q1D.isInitialized() )
+
+ const std::string outputWS("Q1D2Test_result");
+ Q1D.setProperty("DetBankWorkspace", alteredInput);
+ Q1D.setProperty("WavelengthAdj", m_wavNorm);
+ Q1D.setProperty("PixelAdj", m_pixel);
+ Q1D.setPropertyValue("OutputWorkspace", outputWS);
+ Q1D.setPropertyValue("OutputBinning", "0.1,-0.02,0.5");
+ Q1D.setProperty("QResolution", qResolution);
+ Q1D.setProperty("OutputParts", true);
+
+ // Assert
+ TS_ASSERT_THROWS_NOTHING(Q1D.execute());
+
+ TS_ASSERT( Q1D.isExecuted() )
+
+ Mantid::API::MatrixWorkspace_sptr result;
+ TS_ASSERT_THROWS_NOTHING( result = boost::dynamic_pointer_cast<Mantid::API::MatrixWorkspace>
+ (Mantid::API::AnalysisDataService::Instance().retrieve(outputWS)) )
+
+ Mantid::API::MatrixWorkspace_sptr sumOfCounts;
+ TS_ASSERT_THROWS_NOTHING( sumOfCounts = boost::dynamic_pointer_cast<Mantid::API::MatrixWorkspace>
+ (Mantid::API::AnalysisDataService::Instance().retrieve(outputWS+"_sumOfCounts")) )
+
+ // That output will be SUM_i(Yin_i*QRES_in_i)/(SUM_i(Y_in_i)) for each q value
+ // In our test workspace we set QRes_in_1 to 1, this means that all DX values should be SUM_i(Y_in_i)/(SUM_i(Y_in_i))
+ // which is either 1 or 0. It can be 0 if no data falls into this bin. We make sure that there is at least one bin
+ // with a count of 1
+ auto& dataDX = result->dataDx(0);
+ unsigned int counter = 0;
+ for (auto it = dataDX.begin(); it != dataDX.end(); ++it) {
+ if (*it==1.0) {
+ counter++;
+ }
+ // Make sure that the value is either 1 or 0
+ TSM_ASSERT("The DX value should be either 1 or 0", (*it==1.0) || (*it==0.0));
+ }
+ TSM_ASSERT("There should be at least one bin which is not 0", counter > 0);
+
+ // Clean Up
+ Mantid::API::AnalysisDataService::Instance().remove(outputWS);
+ Mantid::API::AnalysisDataService::Instance().remove(outputWS+"_sumOfCounts");
+ Mantid::API::AnalysisDataService::Instance().remove(outputWS+"_sumOfNormFactors");
+ }
+
+ void testThatDxIsNotPopulatedWhenQResolutionIsNotChoosen() {
+ // Arrange
+ Mantid::Algorithms::Q1D2 Q1D;
+ TS_ASSERT_THROWS_NOTHING( Q1D.initialize() );
+ TS_ASSERT( Q1D.isInitialized() )
+
+ const std::string outputWS("Q1D2Test_result");
+ Q1D.setProperty("DetBankWorkspace", m_inputWS);
+ Q1D.setProperty("WavelengthAdj", m_wavNorm);
+ Q1D.setProperty("PixelAdj", m_pixel);
+ Q1D.setPropertyValue("OutputWorkspace", outputWS);
+ Q1D.setPropertyValue("OutputBinning", "0.1,-0.02,0.5");
+ TS_ASSERT_THROWS_NOTHING(Q1D.execute());
+
+ TS_ASSERT( Q1D.isExecuted() )
+
+ Mantid::API::MatrixWorkspace_sptr result;
+ TS_ASSERT_THROWS_NOTHING( result = boost::dynamic_pointer_cast<Mantid::API::MatrixWorkspace>
+ (Mantid::API::AnalysisDataService::Instance().retrieve(outputWS)) )
+
+ // Make sure that the Q resolution is not calculated
+ auto& dataDX = result->dataDx(0);
+ for (auto it = dataDX.begin(); it != dataDX.end(); ++it) {
+ TSM_ASSERT("All Dx values should be 0, as we didn't use the QResolution ooption", (*it==0.0));
+ }
+ Mantid::API::AnalysisDataService::Instance().remove(outputWS);
+ }
+
///stop the constructor from being run every time algorithms test suite is initialised
static Q1D2Test *createSuite() { return new Q1D2Test(); }
static void destroySuite(Q1D2Test *suite) { delete suite; }
@@ -483,4 +573,25 @@ void createInputWorkSpacesForMasking ( Mantid::API::MatrixWorkspace_sptr & input
{
createInputWorkspaces ( 9001, 9030, input, wave, pixels );
}
+
+void createQResolutionWorkspace(Mantid::API::MatrixWorkspace_sptr & qResolution, Mantid::API::MatrixWorkspace_sptr & alteredInput, Mantid::API::MatrixWorkspace_sptr & input, double value1, double value2) {
+ //The q resolution workspace is almost the same to the input workspace, except for the y value, we set all Y values to 1
+ qResolution = Mantid::API::MatrixWorkspace_sptr(input->clone().release());
+ alteredInput = Mantid::API::MatrixWorkspace_sptr(input->clone().release());
+
+ // Populate Y with Value1
+ for (size_t i = 0; i < qResolution->getNumberHistograms(); ++i) {
+ auto& data = qResolution->dataY(i);
+ std::fill(data.begin(), data.end(), value1);
+ }
+
+ // Populate Y with Value2
+ for (size_t i = 0; i < alteredInput->getNumberHistograms(); ++i) {
+ auto& data = alteredInput->dataY(i);
+ std::fill(data.begin(), data.end(), value2);
+ }
+
+
+}
+
#endif /*Q1D2Test_H_*/