Commit 16f43a39 authored by Whitfield, Ross's avatar Whitfield, Ross
Browse files

Allow PeaksWorkspace createPeakQSample to calculate goniometer for CW

Moved code from FindPeaksMD into new method on Goniometer. Then added option onto peaks workspace.
parent 5176c76c
......@@ -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
......
......@@ -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
......
......@@ -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 {
......
......@@ -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);
......
......@@ -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);
......
......@@ -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} =
......
......@@ -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
---------------------------------
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment