From 16f43a3977f689e8293ee7e667a005153b623149 Mon Sep 17 00:00:00 2001
From: Ross Whitfield <whitfieldre@ornl.gov>
Date: Tue, 24 Apr 2018 15:59:08 -0400
Subject: [PATCH] Allow PeaksWorkspace createPeakQSample to calculate
goniometer for CW
Moved code from FindPeaksMD into new method on Goniometer. Then added option onto peaks workspace.
---
Framework/DataObjects/src/PeaksWorkspace.cpp | 22 +++++++-
.../MantidGeometry/Instrument/Goniometer.h | 4 ++
.../Geometry/src/Instrument/Goniometer.cpp | 32 ++++++++++++
Framework/Geometry/test/GoniometerTest.h | 50 +++++++++++++++++++
Framework/MDAlgorithms/src/FindPeaksMD.cpp | 28 ++---------
docs/source/algorithms/FindPeaksMD-v1.rst | 2 +
docs/source/concepts/PeaksWorkspace.rst | 21 ++++++++
7 files changed, 134 insertions(+), 25 deletions(-)
diff --git a/Framework/DataObjects/src/PeaksWorkspace.cpp b/Framework/DataObjects/src/PeaksWorkspace.cpp
index 7cbe6aea457..8f8ecd91dc2 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>
@@ -278,7 +279,26 @@ 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);
+ } 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 3a1b47ce811..7ffbb2626a3 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 3bae3f88603..c672be30d65 100644
--- a/Framework/Geometry/src/Instrument/Goniometer.cpp
+++ b/Framework/Geometry/src/Instrument/Goniometer.cpp
@@ -178,6 +178,38 @@ void Goniometer::setRotationAngle(size_t axisnumber, double value) {
recalculateR();
}
+/**Calculate goniometer for rotation around y-asix for constant wavelength from
+ * Q Sample
+ * @param Q :: Q Sample position in reciprocal space
+ * @param wavelength :: wavelength
+*/
+void Goniometer::calcFromQSampleAndWavelength(const Mantid::Kernel::V3D &Q,
+ double wavelength) {
+ 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]);
+ g_log.information() << "Found goniometer rotation to be " << rot * 180 / M_PI
+ << " degrees for 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);
+ 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 7cfad5d1b3f..afcdf3a9a3a 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 9856bc510dc..7c19297f430 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,9 @@ 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]);
- 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);
+ Geometry::Goniometer goniometer;
+ goniometer.calcFromQSampleAndWavelength(Q, wavelength);
+ 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 924f356c771..48f18985df8 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 d6c43dbd29d..1195279873f 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
---------------------------------
--
GitLab