Get FindPeaksMD working for DEMAND

Framework/Geometry/inc/MantidGeometry/Instrument/Goniometer.h

@@ -79,7 +79,8 @@ public:
   // Calculate goniometer for rotation around y-axis for constant wavelength
   // from Q Sample
   void calcFromQSampleAndWavelength(const Mantid::Kernel::V3D &Q,
-                                    double wavelength);
+                                    double wavelength, bool flip_x = false,
+                                    bool inner = false);
   // Get axis object
   const GoniometerAxis &getAxis(size_t axisnumber) const;
   // Get axis object

Framework/Geometry/src/Instrument/Goniometer.cpp

@@ -192,20 +192,39 @@ void Goniometer::setRotationAngle(size_t axisnumber, double value) {
  * Q Sample
  * @param position :: Q Sample position in reciprocal space
  * @param wavelength :: wavelength
+ * @param flip_x :: option if the q_lab x value should be negative, hence the
+ * detector of the other side of the beam
+ * @param inner :: whether the goniometer is the most inner (phi) or most outer
+ * (omega)
  */
 void Goniometer::calcFromQSampleAndWavelength(
-    const Mantid::Kernel::V3D &position, double wavelength) {
+    const Mantid::Kernel::V3D &position, double wavelength, bool flip_x,
+    bool inner) {
   V3D Q(position);
   if (Kernel::ConfigService::Instance().getString("Q.convention") ==
       "Crystallography")
     Q *= -1;
+
+  Matrix<double> starting_goniometer = getR();
+
+  if (!inner)
+    Q = starting_goniometer * Q;
+
   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]
+  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)));
 
+  if (flip_x)
+    Q_lab[0] *= -1;
+
+  if (inner) {
+    starting_goniometer.Invert();
+    Q_lab = starting_goniometer * Q_lab;
+  }
+
   // 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);
@@ -219,7 +238,12 @@ void Goniometer::calcFromQSampleAndWavelength(
   goniometer[0][2] = sin(rot);
   goniometer[2][0] = -sin(rot);
   goniometer[2][2] = cos(rot);
-  setR(goniometer);
+  if (inner) {
+    starting_goniometer.Invert();
+    setR(starting_goniometer * goniometer);
+  } else {
+    setR(goniometer * starting_goniometer);
+  }
 }
 
 /// Get GoniometerAxis obfject using motor number

Framework/MDAlgorithms/src/FindPeaksMD.cpp

@@ -206,12 +206,28 @@ void FindPeaksMD::init() {
   declareProperty("Wavelength", DBL_MAX, nonNegativeDbl,
                   "Wavelength to use when calculating goniometer angle. If not"
                   "set will use the wavelength parameter on the instrument.");
-
   setPropertySettings("Wavelength",
                       std::make_unique<EnabledWhenProperty>(
                           "CalculateGoniometerForCW",
                           Mantid::Kernel::ePropertyCriterion::IS_NOT_DEFAULT));
 
+  declareProperty("InnerGoniometer", false,
+                  "Whether the goniometer to be calculated is the most inner "
+                  "(phi) or most outer (omega)");
+  setPropertySettings("InnerGoniometer",
+                      std::make_unique<EnabledWhenProperty>(
+                          "CalculateGoniometerForCW",
+                          Mantid::Kernel::ePropertyCriterion::IS_NOT_DEFAULT));
+
+  declareProperty("FlipX", false,
+                  "Used when calculating goniometer angle if the q_lab x value "
+                  "should be negative, hence the detector of the other side "
+                  "(right) of the beam");
+  setPropertySettings("FlipX",
+                      std::make_unique<EnabledWhenProperty>(
+                          "CalculateGoniometerForCW",
+                          Mantid::Kernel::ePropertyCriterion::IS_NOT_DEFAULT));
+
   declareProperty(std::make_unique<WorkspaceProperty<PeaksWorkspace>>(
                       "OutputWorkspace", "", Direction::Output),
                   "An output PeaksWorkspace with the peaks' found positions.");
@@ -314,11 +330,14 @@ FindPeaksMD::createPeak(const Mantid::Kernel::V3D &Q, const double binCount,
         }
       }
 
-      Geometry::Goniometer goniometer;
-      goniometer.calcFromQSampleAndWavelength(Q, wavelength);
-      g_log.information() << "Found goniometer rotation to be "
-                          << goniometer.getEulerAngles()[0]
-                          << " degrees for Q sample = " << Q << "\n";
+      Geometry::Goniometer goniometer(m_goniometer);
+      goniometer.calcFromQSampleAndWavelength(
+          Q, wavelength, getProperty("FlipX"), getProperty("InnerGoniometer"));
+      std::vector<double> angles = goniometer.getEulerAngles("YZY");
+      g_log.information()
+          << "Found goniometer rotation to be in YZY convention [" << angles[0]
+          << ", " << angles[1] << ". " << angles[2]
+          << "] degrees for Q sample = " << Q << "\n";
       p = boost::make_shared<Peak>(inst, Q, goniometer.getR());
 
     } else {

docs/source/algorithms/FindPeaksMD-v1.rst

@@ -62,6 +62,11 @@ works for a constant wavelength source and only for Q sample
 workspaces. It also assumes the goniometer rotation is around the
 y-axis only.
 
+You can specify if the goniometer to calculated is omega or phi, outer
+or inner respectively, by setting the option InnerGoniometer. The
+existing goniometer on the workspace will then be taken into account
+correctly.
+
 The goniometer (:math:`G`) is calculated from
 :math:`\textbf{Q}_{sample}` for a given wavelength (:math:`\lambda`)
 by:
@@ -83,9 +88,12 @@ First calculate the :math:`\textbf{Q}_{lab}` using
 
 .. math:: \theta = \cos^{-1}(1-\frac{|\textbf{Q}_{sample}|^2}{2k^2})
 
-.. math:: \phi = \sin^{-1}(-\textbf{Q}_{sample}^y \sin(\theta)/k)
+.. math:: \phi = \sin^{-1}(-\frac{\textbf{Q}_{sample}^y}{k \sin(\theta)})
 
-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.
+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, unless
+the option FlipX is set to True then the :math:`\textbf{Q}_{lab}^x = -\textbf{Q}_{lab}^x`.
 
 Now you have :math:`\theta`, :math:`\phi` and k you can get :math:`\textbf{Q}_{lab}` using (1).
 
@@ -108,7 +116,7 @@ make
 
 .. math:: A = \begin{bmatrix}
           \textbf{Q}_{sample}^x & \textbf{Q}_{sample}^z \\
-          \textbf{Q}_{sample}^z & \textbf{Q}_{sample}^x
+          \textbf{Q}_{sample}^z & -\textbf{Q}_{sample}^x
           \end{bmatrix}
 
 .. math:: B = \begin{bmatrix}