diff --git a/Framework/DataObjects/src/PeaksWorkspace.cpp b/Framework/DataObjects/src/PeaksWorkspace.cpp
index 7cbe6aea457b9b5fb8b05259675843ad0cf027dd..a4708efa97b393b349f6880997cf9c5600f87aa5 100644
--- a/Framework/DataObjects/src/PeaksWorkspace.cpp
+++ b/Framework/DataObjects/src/PeaksWorkspace.cpp
@@ -16,6 +16,7 @@
#include "MantidKernel/PhysicalConstants.h"
#include "MantidKernel/Quat.h"
#include "MantidKernel/Unit.h"
+#include "MantidKernel/UnitConversion.h"
#include "MantidKernel/V3D.h"
#include <algorithm>
@@ -39,6 +40,8 @@ namespace DataObjects {
/// Register the workspace as a type
DECLARE_WORKSPACE(PeaksWorkspace)
+Mantid::Kernel::Logger g_log("PeaksWorkspace");
+
//---------------------------------------------------------------------------------------------
/** Constructor. Create a table with all the required columns.
*
@@ -278,7 +281,29 @@ PeaksWorkspace::createPeak(const Kernel::V3D &position,
*/
Peak *PeaksWorkspace::createPeakQSample(const V3D &position) const {
// Create a peak from QSampleFrame
- const auto goniometer = run().getGoniometer();
+
+ Geometry::Goniometer goniometer;
+
+ LogManager_const_sptr props = getLogs();
+ if (props->hasProperty("wavelength") || props->hasProperty("energy")) {
+ // Assume constant wavelenth
+ // Calculate Q lab from Q sample and wavelength
+ double wavelength;
+ if (props->hasProperty("energy")) {
+ wavelength = Kernel::UnitConversion::run(
+ "Energy", "Wavelength",
+ props->getPropertyValueAsType<double>("energy"), 0, 0, 0,
+ Kernel::DeltaEMode::Elastic, 0);
+ } else {
+ wavelength = props->getPropertyValueAsType<double>("wavelength");
+ }
+ goniometer.calcFromQSampleAndWavelength(position, wavelength);
+ g_log.information() << "Found goniometer rotation to be "
+ << goniometer.getEulerAngles()[0]
+ << " degrees for Q sample = " << position << "\n";
+ } else {
+ goniometer = run().getGoniometer();
+ }
// create a peak using the qLab frame
auto peak = new Peak(getInstrument(), position, goniometer.getR());
// Take the run number from this
diff --git a/Framework/Geometry/inc/MantidGeometry/Instrument/Goniometer.h b/Framework/Geometry/inc/MantidGeometry/Instrument/Goniometer.h
index 3a1b47ce8111db2c41327ff230ef805e114ea4f4..7ffbb2626a3298f3df831123776fdfa312a4e250 100644
--- a/Framework/Geometry/inc/MantidGeometry/Instrument/Goniometer.h
+++ b/Framework/Geometry/inc/MantidGeometry/Instrument/Goniometer.h
@@ -91,6 +91,10 @@ public:
// Set rotation angle for an axis in the units the angle is set (default --
// degrees)
void setRotationAngle(size_t axisnumber, double value);
+ // Calculate goniometer for rotation around y-axis for constant wavelength
+ // from Q Sample
+ void calcFromQSampleAndWavelength(const Mantid::Kernel::V3D &Q,
+ double wavelength);
// Get axis object
const GoniometerAxis &getAxis(size_t axisnumber) const;
// Get axis object
diff --git a/Framework/Geometry/src/Instrument/Goniometer.cpp b/Framework/Geometry/src/Instrument/Goniometer.cpp
index 3bae3f88603880e33666d42fce6e9d27970da9c9..1002fd2ed07877a088b42ea910ef460aa4520cc9 100644
--- a/Framework/Geometry/src/Instrument/Goniometer.cpp
+++ b/Framework/Geometry/src/Instrument/Goniometer.cpp
@@ -6,6 +6,7 @@
#include <vector>
#include <boost/algorithm/string.hpp>
#include <cstdlib>
+#include "MantidKernel/ConfigService.h"
#include "MantidKernel/Strings.h"
#include "MantidKernel/Logger.h"
@@ -178,6 +179,40 @@ void Goniometer::setRotationAngle(size_t axisnumber, double value) {
recalculateR();
}
+/**Calculate goniometer for rotation around y-asix for constant wavelength from
+ * Q Sample
+ * @param position :: Q Sample position in reciprocal space
+ * @param wavelength :: wavelength
+*/
+void Goniometer::calcFromQSampleAndWavelength(
+ const Mantid::Kernel::V3D &position, double wavelength) {
+ V3D Q(position);
+ if (Kernel::ConfigService::Instance().getString("Q.convention") ==
+ "Crystallography")
+ Q *= -1;
+ double wv = 2.0 * M_PI / wavelength;
+ double norm_q2 = Q.norm2();
+ double theta = acos(1 - norm_q2 / (2 * wv * wv)); // [0, pi]
+ double phi = asin(-Q[1] / wv * sin(theta)); // [-pi/2, pi/2]
+ V3D Q_lab(-wv * sin(theta) * cos(phi), -wv * sin(theta) * sin(phi),
+ wv * (1 - cos(theta)));
+
+ // Solve to find rotation matrix, assuming only rotation around y-axis
+ // A * X = B
+ Matrix<double> A({Q[0], Q[2], Q[2], -Q[0]}, 2, 2);
+ A.Invert();
+ std::vector<double> B{Q_lab[0], Q_lab[2]};
+ std::vector<double> X = A * B;
+ double rot = atan2(X[1], X[0]);
+
+ Matrix<double> goniometer(3, 3, true);
+ goniometer[0][0] = cos(rot);
+ goniometer[0][2] = sin(rot);
+ goniometer[2][0] = -sin(rot);
+ goniometer[2][2] = cos(rot);
+ setR(goniometer);
+}
+
/// Get GoniometerAxis obfject using motor number
/// @param axisnumber :: axis number (from 0)
const GoniometerAxis &Goniometer::getAxis(size_t axisnumber) const {
diff --git a/Framework/Geometry/test/GoniometerTest.h b/Framework/Geometry/test/GoniometerTest.h
index 7cfad5d1b3fb4dfaac13bceeaa4b28c1b3f75f3c..afcdf3a9a3a604cb29b4c96ba3e131a8967a1e44 100644
--- a/Framework/Geometry/test/GoniometerTest.h
+++ b/Framework/Geometry/test/GoniometerTest.h
@@ -179,6 +179,56 @@ public:
}
}
+ void test_calcFromQSampleAndWavelength() {
+ Goniometer G;
+ double wavelength = 2 * M_PI; // k=1
+ DblMatrix R;
+ V3D Q;
+
+ // 0 degree rotation
+ Q = V3D(-1, 0, 1);
+ G.calcFromQSampleAndWavelength(Q, wavelength);
+ R = G.getR();
+ TS_ASSERT_DELTA(R[0][0], 1.0, 0.001);
+ TS_ASSERT_DELTA(R[0][1], 0.0, 0.001);
+ TS_ASSERT_DELTA(R[0][2], 0.0, 0.001);
+ TS_ASSERT_DELTA(R[1][0], 0.0, 0.001);
+ TS_ASSERT_DELTA(R[1][1], 1.0, 0.001);
+ TS_ASSERT_DELTA(R[1][2], 0.0, 0.001);
+ TS_ASSERT_DELTA(R[2][0], 0.0, 0.001);
+ TS_ASSERT_DELTA(R[2][1], 0.0, 0.001);
+ TS_ASSERT_DELTA(R[2][2], 1.0, 0.001);
+
+ // -90 degree rotation
+ Q = V3D(1, 0, 1);
+ G.calcFromQSampleAndWavelength(Q, wavelength);
+ R = G.getR();
+ TS_ASSERT_DELTA(R[0][0], 0.0, 0.001);
+ TS_ASSERT_DELTA(R[0][1], 0.0, 0.001);
+ TS_ASSERT_DELTA(R[0][2], -1.0, 0.001);
+ TS_ASSERT_DELTA(R[1][0], 0.0, 0.001);
+ TS_ASSERT_DELTA(R[1][1], 1.0, 0.001);
+ TS_ASSERT_DELTA(R[1][2], 0.0, 0.001);
+ TS_ASSERT_DELTA(R[2][0], 1.0, 0.001);
+ TS_ASSERT_DELTA(R[2][1], 0.0, 0.001);
+ TS_ASSERT_DELTA(R[2][2], 0.0, 0.001);
+
+ // 30 degree rotation
+ wavelength = 1.54;
+ Q = V3D(-0.63523489, -0.12302677, -0.29517982);
+ G.calcFromQSampleAndWavelength(Q, wavelength);
+ R = G.getR();
+ TS_ASSERT_DELTA(R[0][0], 0.866, 0.001);
+ TS_ASSERT_DELTA(R[0][1], 0.0, 0.001);
+ TS_ASSERT_DELTA(R[0][2], 0.5, 0.01);
+ TS_ASSERT_DELTA(R[1][0], 0.0, 0.001);
+ TS_ASSERT_DELTA(R[1][1], 1.0, 0.001);
+ TS_ASSERT_DELTA(R[1][2], 0.0, 0.001);
+ TS_ASSERT_DELTA(R[2][0], -0.5, 0.01);
+ TS_ASSERT_DELTA(R[2][1], 0.0, 0.001);
+ TS_ASSERT_DELTA(R[2][2], 0.866, 0.001);
+ }
+
/** Save and load to NXS file */
void test_nexus() {
NexusTestHelper th(true);
diff --git a/Framework/MDAlgorithms/src/FindPeaksMD.cpp b/Framework/MDAlgorithms/src/FindPeaksMD.cpp
index 9856bc510dc0f0a7615eec1aed5d915113910279..b71032b5f837448c74b838ba9d35c720bdd9bcd2 100644
--- a/Framework/MDAlgorithms/src/FindPeaksMD.cpp
+++ b/Framework/MDAlgorithms/src/FindPeaksMD.cpp
@@ -4,6 +4,7 @@
#include "MantidDataObjects/MDHistoWorkspace.h"
#include "MantidDataObjects/PeaksWorkspace.h"
#include "MantidGeometry/Crystal/EdgePixel.h"
+#include "MantidGeometry/Instrument/Goniometer.h"
#include "MantidGeometry/Objects/InstrumentRayTracer.h"
#include "MantidKernel/BoundedValidator.h"
#include "MantidKernel/EnabledWhenProperty.h"
@@ -295,30 +296,12 @@ FindPeaksMD::createPeak(const Mantid::Kernel::V3D &Q, const double binCount,
if (calcGoniometer) {
// Calculate Q lab from Q sample and wavelength
double wavelength = getProperty("Wavelength");
- double wv = 2.0 * M_PI / wavelength;
- double norm_q2 = Q.norm2();
- double theta = acos(1 - norm_q2 / (2 * wv * wv));
- double phi = asin(-Q[1] / wv * sin(theta));
- V3D Q_lab(-wv * sin(theta) * cos(phi), -wv * sin(theta) * sin(phi),
- wv * (1 - cos(theta)));
-
- // Solve to find rotation matrix, assuming only rotation around y-axis
- // A * X = B
- Matrix<double> A({Q[0], Q[2], Q[2], -Q[0]}, 2, 2);
- A.Invert();
- std::vector<double> B{Q_lab[0], Q_lab[2]};
- std::vector<double> X = A * B;
- double rot = atan2(X[1], X[0]);
+ Geometry::Goniometer goniometer;
+ goniometer.calcFromQSampleAndWavelength(Q, wavelength);
g_log.information() << "Found goniometer rotation to be "
- << rot * 180 / M_PI
- << " degrees for peak at Q sample = " << Q << "\n";
-
- Matrix<double> goniometer(3, 3, true);
- goniometer[0][0] = cos(rot);
- goniometer[0][2] = sin(rot);
- goniometer[2][0] = -sin(rot);
- goniometer[2][2] = cos(rot);
- p = boost::make_shared<Peak>(inst, Q, goniometer);
+ << goniometer.getEulerAngles()[0]
+ << " degrees for Q sample = " << Q << "\n";
+ p = boost::make_shared<Peak>(inst, Q, goniometer.getR());
} else {
p = boost::make_shared<Peak>(inst, Q, m_goniometer);
diff --git a/docs/source/algorithms/FindPeaksMD-v1.rst b/docs/source/algorithms/FindPeaksMD-v1.rst
index 924f356c771cc84af50d15d79f181f67d191252c..48f18985df86a044e7c1f615d2ce8d490ad3e582 100644
--- a/docs/source/algorithms/FindPeaksMD-v1.rst
+++ b/docs/source/algorithms/FindPeaksMD-v1.rst
@@ -85,6 +85,8 @@ First calculate the :math:`\textbf{Q}_{lab}` using
.. math:: \phi = \sin^{-1}(-\textbf{Q}_{sample}^y \sin(\theta)/k)
+where :math:`\theta` is from 0 to :math:`\pi` and :math:`\phi` is from :math:`-\pi/2` to :math:`\pi/2`. This means that it will assume your detector position is on the left of the beam even it it's not.
+
Now you have :math:`\theta`, :math:`\phi` and k you can get :math:`\textbf{Q}_{lab}` using (1).
We need to now solve :math:`G \textbf{Q}_{sample} =
diff --git a/docs/source/concepts/PeaksWorkspace.rst b/docs/source/concepts/PeaksWorkspace.rst
index d6c43dbd29df77dd17d972089f366459c41f7c73..1195279873f749a0475002717924a79335cf1d04 100644
--- a/docs/source/concepts/PeaksWorkspace.rst
+++ b/docs/source/concepts/PeaksWorkspace.rst
@@ -49,6 +49,27 @@ Each Peak object contains a PeakShape. Only the integration algorithms which act
Subtypes of PeakShape will then provide additional information. For example PeakShapeSpherical provides the radius as well as background inner, and background outer radius.
+Calculate Goniometer For Constant Wavelength
+--------------------------------------------
+
+If you set the `wavelength` (in Ã…) or `energy` (in meV) property on a
+PeaksWorkspace when the createPeak method is used the goniometer
+rotation with be calculated. This allows you to use one instrument
+definition for multiple goniometer rotations, for example adding peaks
+in Slice Viewer from multiple combined MD workspaces. It only works
+for a constant wavelength source and only for Q sample workspaces. It
+also assumes the goniometer rotation is around the y-axis only. For
+details on the calculation see "Calculate Goniometer For Constant
+Wavelength" at :ref:`FindPeaksMD <algm-FindPeaksMD>`.
+
+.. code-block:: python
+
+ pws = mtd['name_of_peaks_workspace']
+ pws.run().addProperty('wavelength', 1.54, True)
+ # or
+ pws.run().addProperty('energy', 34.48, True)
+
+
Using PeaksWorkspaces in Python
---------------------------------