diff --git a/Code/Mantid/Framework/MDEvents/inc/MantidMDEvents/MDHistoWorkspace.h b/Code/Mantid/Framework/MDEvents/inc/MantidMDEvents/MDHistoWorkspace.h
index ef9e1fb3ef3d11c1480515a422fab94808025f00..0aff81ab35ef56bc0c0da1dcace3436fc731284a 100644
--- a/Code/Mantid/Framework/MDEvents/inc/MantidMDEvents/MDHistoWorkspace.h
+++ b/Code/Mantid/Framework/MDEvents/inc/MantidMDEvents/MDHistoWorkspace.h
@@ -75,6 +75,12 @@ namespace MDEvents
MDHistoWorkspace & operator/=(const MDHistoWorkspace & b);
MDHistoWorkspace & operator/=(const Mantid::DataObjects::WorkspaceSingleValue & b);
+ MDHistoWorkspace & operator&=(const MDHistoWorkspace & b);
+ MDHistoWorkspace & operator|=(const MDHistoWorkspace & b);
+ MDHistoWorkspace & operator^=(const MDHistoWorkspace & b);
+
+ void operatorNot();
+
// --------------------------------------------------------------------------------------------
diff --git a/Code/Mantid/Framework/MDEvents/src/MDHistoWorkspace.cpp b/Code/Mantid/Framework/MDEvents/src/MDHistoWorkspace.cpp
index 8a4badd49f43d10a15e76b14f35318795007553a..b558573ad0edabb3b8ef9ba9a68aa32a33c497c0 100644
--- a/Code/Mantid/Framework/MDEvents/src/MDHistoWorkspace.cpp
+++ b/Code/Mantid/Framework/MDEvents/src/MDHistoWorkspace.cpp
@@ -324,14 +324,22 @@ namespace MDEvents
return out;
}
+ //----------------------------------------------------------------------------------------------
+ /** Check if the two workspace's sizes match (for comparison or element-by-element operation
+ *
+ * @param other :: the workspace to compare to
+ * @param operation :: descriptive string (for the error message)
+ * @throw an error if they don't match
+ */
void MDHistoWorkspace::checkWorkspaceSize(const MDHistoWorkspace & other, std::string operation)
{
if (other.getNumDims() != this->getNumDims())
- throw std::invalid_argument("Cannot perform the " + operation + " operation on this workspace. The number of dimensions does not match.");
+ throw std::invalid_argument("Cannot perform the " + operation + " operation on this MDHistoWorkspace. The number of dimensions does not match.");
if (other.m_length != this->m_length)
- throw std::invalid_argument("Cannot perform the " + operation + " operation on this workspace. The length of the signals vector does not match.");
+ throw std::invalid_argument("Cannot perform the " + operation + " operation on this MDHistoWorkspace. The length of the signals vector does not match.");
}
+
//----------------------------------------------------------------------------------------------
/** Perform the += operation, element-by-element, for two MDHistoWorkspace's
*
@@ -366,6 +374,220 @@ namespace MDEvents
return *this;
}
+ //----------------------------------------------------------------------------------------------
+ /** Perform the -= operation, element-by-element, for two MDHistoWorkspace's
+ *
+ * @param b :: workspace on the RHS of the operation
+ * @return *this after operation */
+ MDHistoWorkspace & MDHistoWorkspace::operator-=(const MDHistoWorkspace & b)
+ {
+ checkWorkspaceSize(b, "-=");
+ for (size_t i=0; i<m_length; ++i)
+ {
+ m_signals[i] -= b.m_signals[i];
+ m_errorsSquared[i] += b.m_errorsSquared[i];
+ }
+ return *this;
+ }
+
+ //----------------------------------------------------------------------------------------------
+ /** Perform the -= operation with a scalar as the RHS argument
+ *
+ * @param b :: WorkspaceSingleValue (signal and error) as the RHS argument
+ * @return *this after operation */
+ MDHistoWorkspace & MDHistoWorkspace::operator-=(const Mantid::DataObjects::WorkspaceSingleValue & b)
+ {
+ signal_t signal = b.dataY(0)[0];
+ signal_t errorSquared = b.dataE(0)[0];
+ errorSquared *= errorSquared;
+ for (size_t i=0; i<m_length; ++i)
+ {
+ m_signals[i] -= signal;
+ m_errorsSquared[i] += errorSquared;
+ }
+ return *this;
+ }
+
+ //----------------------------------------------------------------------------------------------
+ /** Perform the *= operation, element-by-element, for two MDHistoWorkspace's
+ *
+ * Error propagation of \f$ f = a * b \f$ is given by:
+ * \f$ df^2 = f^2 * (da^2 / a^2 + db^2 / b^2) \f$
+ *
+ * @param b_ws :: workspace on the RHS of the operation
+ * @return *this after operation */
+ MDHistoWorkspace & MDHistoWorkspace::operator*=(const MDHistoWorkspace & b_ws)
+ {
+ checkWorkspaceSize(b_ws, "*=");
+ for (size_t i=0; i<m_length; ++i)
+ {
+ signal_t a = m_signals[i];
+ signal_t da2 = m_errorsSquared[i];
+
+ signal_t b = b_ws.m_signals[i];
+ signal_t db2 = b_ws.m_errorsSquared[i];
+
+ signal_t f = a * b;
+ signal_t df2 = (f * f) * ( da2 / (a*a) + db2 / (b*b) );
+
+ m_signals[i] = f;
+ m_errorsSquared[i] = df2;
+ }
+ return *this;
+ }
+
+ //----------------------------------------------------------------------------------------------
+ /** Perform the *= operation with a scalar as the RHS argument
+ *
+ * Error propagation of \f$ f = a * b \f$ is given by:
+ * \f$ df^2 = f^2 * (da^2 / a^2 + db^2 / b^2) \f$
+ *
+ * @param b_ws :: WorkspaceSingleValue (signal and error) as the RHS argument
+ * @return *this after operation */
+ MDHistoWorkspace & MDHistoWorkspace::operator*=(const Mantid::DataObjects::WorkspaceSingleValue & b_ws)
+ {
+ signal_t b = b_ws.dataY(0)[0];
+ signal_t db2 = b_ws.dataE(0)[0] * b_ws.dataE(0)[0];
+ signal_t db2_relative = db2 / (b*b);
+ for (size_t i=0; i<m_length; ++i)
+ {
+ signal_t a = m_signals[i];
+ signal_t da2 = m_errorsSquared[i];
+
+ signal_t f = a * b;
+ signal_t df2 = (f * f) * ( da2 / (a*a) + db2_relative );
+
+ m_signals[i] = f;
+ m_errorsSquared[i] = df2;
+ }
+ return *this;
+ }
+
+ //----------------------------------------------------------------------------------------------
+ /** Perform the /= operation, element-by-element, for two MDHistoWorkspace's
+ *
+ * Error propagation of \f$ f = a / b \f$ is given by:
+ * \f$ df^2 = f^2 * (da^2 / a^2 + db^2 / b^2) \f$
+ *
+ * @param b_ws :: workspace on the RHS of the operation
+ * @return *this after operation */
+ MDHistoWorkspace & MDHistoWorkspace::operator/=(const MDHistoWorkspace & b_ws)
+ {
+ checkWorkspaceSize(b_ws, "/=");
+ for (size_t i=0; i<m_length; ++i)
+ {
+ signal_t a = m_signals[i];
+ signal_t da2 = m_errorsSquared[i];
+
+ signal_t b = b_ws.m_signals[i];
+ signal_t db2 = b_ws.m_errorsSquared[i];
+
+ signal_t f = a / b;
+ signal_t df2 = (f * f) * ( da2 / (a*a) + db2 / (b*b) );
+
+ m_signals[i] = f;
+ m_errorsSquared[i] = df2;
+ }
+ return *this;
+ }
+
+ //----------------------------------------------------------------------------------------------
+ /** Perform the /= operation with a scalar as the RHS argument
+ *
+ * Error propagation of \f$ f = a / b \f$ is given by:
+ * \f$ df^2 = f^2 * (da^2 / a^2 + db^2 / b^2) \f$
+ *
+ * @param b_ws :: WorkspaceSingleValue (signal and error) as the RHS argument
+ * @return *this after operation */
+ MDHistoWorkspace & MDHistoWorkspace::operator/=(const Mantid::DataObjects::WorkspaceSingleValue & b_ws)
+ {
+ signal_t b = b_ws.dataY(0)[0];
+ signal_t db2 = b_ws.dataE(0)[0] * b_ws.dataE(0)[0];
+ signal_t db2_relative = db2 / (b*b);
+ for (size_t i=0; i<m_length; ++i)
+ {
+ signal_t a = m_signals[i];
+ signal_t da2 = m_errorsSquared[i];
+
+ signal_t f = a / b;
+ signal_t df2 = (f * f) * ( da2 / (a*a) + db2_relative );
+
+ m_signals[i] = f;
+ m_errorsSquared[i] = df2;
+ }
+ return *this;
+ }
+
+
+ //----------------------------------------------------------------------------------------------
+ /** A boolean &= (and) operation, element-by-element, for two MDHistoWorkspace's.
+ *
+ * 0.0 is "false", all other values are "true". All errors are set to 0.
+ *
+ * @param b :: workspace on the RHS of the operation
+ * @return *this after operation */
+ MDHistoWorkspace & MDHistoWorkspace::operator&=(const MDHistoWorkspace & b)
+ {
+ checkWorkspaceSize(b, "&= (and)");
+ for (size_t i=0; i<m_length; ++i)
+ {
+ m_signals[i] = m_signals[i] && b.m_signals[i];
+ m_errorsSquared[i] = 0;
+ }
+ return *this;
+ }
+
+ //----------------------------------------------------------------------------------------------
+ /** A boolean |= (or) operation, element-by-element, for two MDHistoWorkspace's.
+ *
+ * 0.0 is "false", all other values are "true". All errors are set to 0.
+ *
+ * @param b :: workspace on the RHS of the operation
+ * @return *this after operation */
+ MDHistoWorkspace & MDHistoWorkspace::operator|=(const MDHistoWorkspace & b)
+ {
+ checkWorkspaceSize(b, "|= (or)");
+ for (size_t i=0; i<m_length; ++i)
+ {
+ m_signals[i] = m_signals[i] || b.m_signals[i];
+ m_errorsSquared[i] = 0;
+ }
+ return *this;
+ }
+
+ //----------------------------------------------------------------------------------------------
+ /** A boolean ^= (xor) operation, element-by-element, for two MDHistoWorkspace's.
+ *
+ * 0.0 is "false", all other values are "true". All errors are set to 0.
+ *
+ * @param b :: workspace on the RHS of the operation
+ * @return *this after operation */
+ MDHistoWorkspace & MDHistoWorkspace::operator^=(const MDHistoWorkspace & b)
+ {
+ checkWorkspaceSize(b, "^= (xor)");
+ for (size_t i=0; i<m_length; ++i)
+ {
+ m_signals[i] = bool(m_signals[i]) ^ bool(b.m_signals[i]);
+ m_errorsSquared[i] = 0;
+ }
+ return *this;
+ }
+
+ //----------------------------------------------------------------------------------------------
+ /** A boolean not operation, performed in-place.
+ * All errors are set to 0.
+ *
+ * 0.0 is "false", all other values are "true". All errors are set to 0.
+ */
+ void MDHistoWorkspace::operatorNot()
+ {
+ for (size_t i=0; i<m_length; ++i)
+ {
+ m_signals[i] = !bool(m_signals[i]);
+ m_errorsSquared[i] = 0;
+ }
+ }
+
} // namespace Mantid
} // namespace MDEvents
diff --git a/Code/Mantid/Framework/MDEvents/test/MDHistoWorkspaceTest.h b/Code/Mantid/Framework/MDEvents/test/MDHistoWorkspaceTest.h
index ffbba22b13423c16c3114a61bb63bd88cb7413f2..c258d400ad4c79c5109afd12b59713170477d51d 100644
--- a/Code/Mantid/Framework/MDEvents/test/MDHistoWorkspaceTest.h
+++ b/Code/Mantid/Framework/MDEvents/test/MDHistoWorkspaceTest.h
@@ -14,6 +14,7 @@
#include <iomanip>
#include <iostream>
#include "MantidAPI/IMDWorkspace.h"
+#include "MantidDataObjects/WorkspaceSingleValue.h"
using namespace Mantid::MDEvents;
using namespace Mantid::Geometry;
@@ -32,16 +33,25 @@ public:
MDHistoWorkspace_sptr two;
MDHistoWorkspace_sptr three;
- MDHistoWorkspace_sptr ma
-
+ /** Create some fake workspace, 5x5, with 2.0 and 3.0 each */
MDHistoWorkspaceTest()
{
- MDEventsTestHelper::makeFakeMDHistoWorkspace(
- two.
+ two = MDEventsTestHelper::makeFakeMDHistoWorkspace(2.0, 2, 5, 10.0, 2.0);
+ three = MDEventsTestHelper::makeFakeMDHistoWorkspace(3.0, 2, 5, 10.0, 3.0);
}
+ /** Check that a workspace has the right signal/error*/
+ void checkWorkspace(MDHistoWorkspace_sptr ws, double expectedSignal, double expectedErrorSquared)
+ {
+ for (size_t i=0; i < ws->getNPoints(); i++)
+ {
+ TS_ASSERT_DELTA( ws->getSignalAt(i), expectedSignal, 1e-5 );
+ TS_ASSERT_DELTA( ws->getErrorAt(i), expectedErrorSquared, 1e-5 );
+ }
+ }
+ //--------------------------------------------------------------------------------------
void test_constructor()
{
MDHistoDimension_sptr dimX(new MDHistoDimension("X", "x", "m", -10, 10, 5));
@@ -364,7 +374,7 @@ public:
void test_getSignalAtCoord()
{
// 2D workspace with signal[i] = i (linear index)
- MDHistoWorkspace_sptr ws = MDEventsTestHelper::makeFakeMDHistoWorkspace(1.0, 1.0, 2, 10);
+ MDHistoWorkspace_sptr ws = MDEventsTestHelper::makeFakeMDHistoWorkspace(1.0, 2, 10);
for (size_t i=0; i<100; i++)
ws->setSignalAt(i, double(i));
IMDWorkspace_sptr iws(ws);
@@ -379,6 +389,133 @@ public:
TS_ASSERT( boost::math::isnan(iws->getSignalAtCoord(VMD(3.5, 10.02)) ) );
}
+ //--------------------------------------------------------------------------------------
+ void test_plus_ws()
+ {
+ MDHistoWorkspace_sptr a = MDEventsTestHelper::makeFakeMDHistoWorkspace(2.0, 2, 5, 10.0, 2.5 /*errorSquared*/);
+ MDHistoWorkspace_sptr b = MDEventsTestHelper::makeFakeMDHistoWorkspace(3.0, 2, 5, 10.0, 3.5 /*errorSquared*/);
+ *a += *b;
+ checkWorkspace(a, 5.0, 6.0);
+ }
+
+ void test_plus_scalar()
+ {
+ MDHistoWorkspace_sptr a = MDEventsTestHelper::makeFakeMDHistoWorkspace(2.0, 2, 5, 10.0, 2.5 /*errorSquared*/);
+ WorkspaceSingleValue b(3.0, sqrt(3.5) );
+ *a += b;
+ checkWorkspace(a, 5.0, 6.0);
+ }
+
+ //--------------------------------------------------------------------------------------
+ void test_minus_ws()
+ {
+ MDHistoWorkspace_sptr a = MDEventsTestHelper::makeFakeMDHistoWorkspace(3.0, 2, 5, 10.0, 2.5 /*errorSquared*/);
+ MDHistoWorkspace_sptr b = MDEventsTestHelper::makeFakeMDHistoWorkspace(2.0, 2, 5, 10.0, 3.5 /*errorSquared*/);
+ *a -= *b;
+ checkWorkspace(a, 1.0, 6.0);
+ }
+
+ void test_minus_scalar()
+ {
+ MDHistoWorkspace_sptr a = MDEventsTestHelper::makeFakeMDHistoWorkspace(3.0, 2, 5, 10.0, 2.5 /*errorSquared*/);
+ WorkspaceSingleValue b(2.0, sqrt(3.5) );
+ *a -= b;
+ checkWorkspace(a, 1.0, 6.0);
+ }
+
+ //--------------------------------------------------------------------------------------
+ void test_times_ws()
+ {
+ MDHistoWorkspace_sptr a = MDEventsTestHelper::makeFakeMDHistoWorkspace(2.0, 2, 5, 10.0, 2.0 /*errorSquared*/);
+ MDHistoWorkspace_sptr b = MDEventsTestHelper::makeFakeMDHistoWorkspace(3.0, 2, 5, 10.0, 3.0 /*errorSquared*/);
+ *a *= *b;
+ checkWorkspace(a, 6.0, 36. * (.5 + 1./3.));
+ }
+
+ //--------------------------------------------------------------------------------------
+ void test_times_scalar()
+ {
+ MDHistoWorkspace_sptr a = MDEventsTestHelper::makeFakeMDHistoWorkspace(2.0, 2, 5, 10.0, 2.0 /*errorSquared*/);
+ WorkspaceSingleValue b(3.0, sqrt(3.0) );
+ *a *= b;
+ checkWorkspace(a, 6.0, 36. * (.5 + 1./3.));
+ // Scalar without error
+ MDHistoWorkspace_sptr d = MDEventsTestHelper::makeFakeMDHistoWorkspace(2.0, 2, 5, 10.0, 2.0 /*errorSquared*/);
+ WorkspaceSingleValue e(3.0, 0);
+ *d *= e;
+ checkWorkspace(d, 6.0, 9 * 2.0);
+ }
+
+ //--------------------------------------------------------------------------------------
+ void test_divide_ws()
+ {
+ MDHistoWorkspace_sptr a = MDEventsTestHelper::makeFakeMDHistoWorkspace(3.0, 2, 5, 10.0, 3.0 /*errorSquared*/);
+ MDHistoWorkspace_sptr b = MDEventsTestHelper::makeFakeMDHistoWorkspace(2.0, 2, 5, 10.0, 2.0 /*errorSquared*/);
+ *a /= *b;
+ checkWorkspace(a, 1.5, 1.5 * 1.5 * (.5 + 1./3.));
+ }
+
+ //--------------------------------------------------------------------------------------
+ void test_divide_scalar()
+ {
+ MDHistoWorkspace_sptr a = MDEventsTestHelper::makeFakeMDHistoWorkspace(3.0, 2, 5, 10.0, 3.0 /*errorSquared*/);
+ WorkspaceSingleValue b(2.0, sqrt(2.0) );
+ *a /= b;
+ checkWorkspace(a, 1.5, 1.5 * 1.5 * (.5 + 1./3.));
+ }
+
+ //--------------------------------------------------------------------------------------
+ void test_boolean_and()
+ {
+ MDHistoWorkspace_sptr a = MDEventsTestHelper::makeFakeMDHistoWorkspace(1.23, 2, 5, 10.0, 3.0);
+ MDHistoWorkspace_sptr b = MDEventsTestHelper::makeFakeMDHistoWorkspace(2.34, 2, 5, 10.0, 2.0);
+ MDHistoWorkspace_sptr c = MDEventsTestHelper::makeFakeMDHistoWorkspace(0.00, 2, 5, 10.0, 2.0);
+ *a &= *b;
+ checkWorkspace(a, 1.0, 0.0);
+ *b &= *c;
+ checkWorkspace(b, 0.0, 0.0);
+ }
+
+ //--------------------------------------------------------------------------------------
+ void test_boolean_or()
+ {
+ MDHistoWorkspace_sptr a = MDEventsTestHelper::makeFakeMDHistoWorkspace(1.23, 2, 5, 10.0, 3.0);
+ MDHistoWorkspace_sptr b = MDEventsTestHelper::makeFakeMDHistoWorkspace(2.34, 2, 5, 10.0, 2.0);
+ MDHistoWorkspace_sptr c = MDEventsTestHelper::makeFakeMDHistoWorkspace(0.00, 2, 5, 10.0, 2.0);
+ *a |= *b;
+ checkWorkspace(a, 1.0, 0.0);
+ *b |= *c;
+ checkWorkspace(b, 1.0, 0.0);
+ *c |= *c;
+ checkWorkspace(c, 0.0, 0.0);
+ }
+
+ //--------------------------------------------------------------------------------------
+ void test_boolean_xor()
+ {
+ MDHistoWorkspace_sptr a = MDEventsTestHelper::makeFakeMDHistoWorkspace(1.23, 2, 5, 10.0, 3.0);
+ MDHistoWorkspace_sptr b = MDEventsTestHelper::makeFakeMDHistoWorkspace(2.34, 2, 5, 10.0, 2.0);
+ MDHistoWorkspace_sptr c = MDEventsTestHelper::makeFakeMDHistoWorkspace(0.00, 2, 5, 10.0, 2.0);
+ *a ^= *b;
+ checkWorkspace(a, 0.0, 0.0);
+ *b ^= *c;
+ checkWorkspace(b, 1.0, 0.0);
+ *c ^= *c;
+ checkWorkspace(c, 0.0, 0.0);
+ }
+
+ //--------------------------------------------------------------------------------------
+ void test_boolean_operatorNot()
+ {
+ MDHistoWorkspace_sptr a = MDEventsTestHelper::makeFakeMDHistoWorkspace(1.23, 2, 5, 10.0, 3.0);
+ MDHistoWorkspace_sptr b = MDEventsTestHelper::makeFakeMDHistoWorkspace(0.00, 2, 5, 10.0, 2.0);
+ a->operatorNot();
+ checkWorkspace(a, 0.0, 0.0);
+ b->operatorNot();
+ checkWorkspace(b, 1.0, 0.0);
+ }
+
+
};
diff --git a/Code/Mantid/Framework/TestHelpers/inc/MantidTestHelpers/MDEventsTestHelper.h b/Code/Mantid/Framework/TestHelpers/inc/MantidTestHelpers/MDEventsTestHelper.h
index b1645e54061fa36251b78edd398a205f162f0277..710ce1068133b8860f6eaa7a8023205b510d7646 100644
--- a/Code/Mantid/Framework/TestHelpers/inc/MantidTestHelpers/MDEventsTestHelper.h
+++ b/Code/Mantid/Framework/TestHelpers/inc/MantidTestHelpers/MDEventsTestHelper.h
@@ -43,7 +43,8 @@ namespace MDEventsTestHelper
DLL_TESTHELPERS Mantid::MDEvents::MDEventWorkspace3Lean::sptr makeFileBackedMDEW(std::string wsName, bool fileBacked);
/// Make a fake n-dimensional MDHistoWorkspace
- DLL_TESTHELPERS Mantid::MDEvents::MDHistoWorkspace_sptr makeFakeMDHistoWorkspace(double signal, size_t numDims, size_t numBins = 10, double max = 10.0);
+ DLL_TESTHELPERS Mantid::MDEvents::MDHistoWorkspace_sptr makeFakeMDHistoWorkspace(double signal, size_t numDims, size_t numBins = 10, double max = 10.0,
+ double errorSquared=1.0);
//-------------------------------------------------------------------------------------
diff --git a/Code/Mantid/Framework/TestHelpers/src/MDEventsTestHelper.cpp b/Code/Mantid/Framework/TestHelpers/src/MDEventsTestHelper.cpp
index da0dc1d0c46e50cecd8f98a2dde7ce3c4c043add..eaeee5883d1c31c9c17e15cd0b5b0530752fb274 100644
--- a/Code/Mantid/Framework/TestHelpers/src/MDEventsTestHelper.cpp
+++ b/Code/Mantid/Framework/TestHelpers/src/MDEventsTestHelper.cpp
@@ -205,7 +205,7 @@ namespace MDEventsTestHelper
* @return the MDHisto
*/
Mantid::MDEvents::MDHistoWorkspace_sptr makeFakeMDHistoWorkspace(double signal, size_t numDims, size_t numBins,
- double max)
+ double max, double errorSquared)
{
Mantid::MDEvents::MDHistoWorkspace * ws = NULL;
if (numDims ==1)
@@ -236,7 +236,7 @@ namespace MDEventsTestHelper
);
}
Mantid::MDEvents::MDHistoWorkspace_sptr ws_sptr(ws);
- ws_sptr->setTo(signal, signal);
+ ws_sptr->setTo(signal, errorSquared);
return ws_sptr;
}