LeanElasticPeak to calculated wavelength from qLab

Framework/DataObjects/inc/MantidDataObjects/BasePeak.h

@@ -118,6 +118,12 @@ public:
   void setAbsorptionWeightedPathLength(double pathLength) override;
   double getAbsorptionWeightedPathLength() const override;
 
+protected:
+  double calculateWavelengthFromQLab(const Mantid::Kernel::V3D qLab);
+
+  // ki-kf for Inelastic convention; kf-ki for Crystallography convention
+  std::string convention;
+
 private:
   /// Name of the parent bank
   std::string m_bankName;

Framework/DataObjects/inc/MantidDataObjects/LeanElasticPeak.h

@@ -8,6 +8,7 @@
 
 #include "MantidDataObjects/BasePeak.h"
 #include "MantidGeometry/Crystal/IPeak.h"
+#include "MantidGeometry/Instrument/ReferenceFrame.h"
 #include "MantidKernel/Logger.h"
 #include "MantidKernel/Matrix.h"
 #include "MantidKernel/PhysicalConstants.h"
@@ -34,12 +35,11 @@ public:
 
   LeanElasticPeak();
   LeanElasticPeak(const Mantid::Kernel::V3D &QSampleFrame);
-  LeanElasticPeak(const Mantid::Kernel::V3D &QSampleFrame,
-                  const Mantid::Kernel::Matrix<double> &goniometer);
-  LeanElasticPeak(const Mantid::Kernel::V3D &QSampleFrame, double wavelength);
   LeanElasticPeak(const Mantid::Kernel::V3D &QSampleFrame,
                   const Mantid::Kernel::Matrix<double> &goniometer,
-                  double wavelength);
+                  boost::optional<std::shared_ptr<Geometry::ReferenceFrame>>
+                      refFrame = boost::none);
+  LeanElasticPeak(const Mantid::Kernel::V3D &QSampleFrame, double wavelength);
 
   /// Copy constructor
   LeanElasticPeak(const LeanElasticPeak &other);
@@ -66,6 +66,9 @@ public:
   void setInstrument(const Geometry::Instrument_const_sptr &) override;
   Geometry::IDetector_const_sptr getDetector() const override;
   Geometry::Instrument_const_sptr getInstrument() const override;
+  void setReferenceFrame(std::shared_ptr<Geometry::ReferenceFrame> frame);
+  std::shared_ptr<const Geometry::ReferenceFrame>
+  getReferenceFrame() const override;
 
   bool findDetector() override;
   bool findDetector(const Geometry::InstrumentRayTracer &) override;
@@ -80,6 +83,8 @@ public:
 
   void setQSampleFrame(const Mantid::Kernel::V3D &QSampleFrame,
                        boost::optional<double> = boost::none) override;
+  void setQSampleFrame(const Mantid::Kernel::V3D &QSampleFrame,
+                       const Mantid::Kernel::Matrix<double> &goniometer);
   void setQLabFrame(const Mantid::Kernel::V3D &qLab,
                     boost::optional<double> = boost::none) override;
 
@@ -108,9 +113,11 @@ private:
   /// Q_sample vector
   Mantid::Kernel::V3D m_Qsample;
 
-  /// Initial energy of neutrons at the peak
+  /// Wavelength of neutrons at the peak
   double m_wavelength;
 
+  std::shared_ptr<Geometry::ReferenceFrame> m_refFrame;
+
   /// Static logger
   static Mantid::Kernel::Logger g_log;
 };

Framework/DataObjects/inc/MantidDataObjects/Peak.h

@@ -85,6 +85,8 @@ public:
   void setInstrument(const Geometry::Instrument_const_sptr &inst) override;
   Geometry::IDetector_const_sptr getDetector() const override;
   Geometry::Instrument_const_sptr getInstrument() const override;
+  std::shared_ptr<const Geometry::ReferenceFrame>
+  getReferenceFrame() const override;
 
   bool findDetector() override;
   bool findDetector(const Geometry::InstrumentRayTracer &tracer) override;
@@ -159,9 +161,6 @@ private:
 
   /// Static logger
   static Mantid::Kernel::Logger g_log;
-
-  // ki-kf for Inelastic convention; kf-ki for Crystallography convention
-  std::string convention;
 };
 
 } // namespace DataObjects

Framework/DataObjects/src/BasePeak.cpp

@@ -36,14 +36,15 @@ BasePeak::BasePeak()
       m_GoniometerMatrix(3, 3, true), m_InverseGoniometerMatrix(3, 3, true),
       m_runNumber(0), m_monitorCount(0), m_row(-1), m_col(-1), m_peakNumber(0),
       m_intHKL(V3D(0, 0, 0)), m_intMNP(V3D(0, 0, 0)),
-      m_peakShape(std::make_shared<NoShape>()) {}
+      m_peakShape(std::make_shared<NoShape>()) {
+  convention = Kernel::ConfigService::Instance().getString("Q.convention");
+}
 
 //----------------------------------------------------------------------------------------------
 /** Constructor including goniometer
  *
  * @param goniometer :: a 3x3 rotation matrix
  */
-
 BasePeak::BasePeak(const Mantid::Kernel::Matrix<double> &goniometer)
     : m_H(0), m_K(0), m_L(0), m_intensity(0), m_sigmaIntensity(0),
       m_binCount(0), m_absorptionWeightedPathLength(0),
@@ -51,15 +52,18 @@ BasePeak::BasePeak(const Mantid::Kernel::Matrix<double> &goniometer)
       m_runNumber(0), m_monitorCount(0), m_row(-1), m_col(-1), m_peakNumber(0),
       m_intHKL(V3D(0, 0, 0)), m_intMNP(V3D(0, 0, 0)),
       m_peakShape(std::make_shared<NoShape>()) {
+  convention = Kernel::ConfigService::Instance().getString("Q.convention");
+
   if (fabs(m_InverseGoniometerMatrix.Invert()) < 1e-8)
     throw std::invalid_argument(
         "BasePeak::ctor(): Goniometer matrix must non-singular.");
 }
 
 BasePeak::BasePeak(const BasePeak &other)
-    : m_bankName(other.m_bankName), m_H(other.m_H), m_K(other.m_K),
-      m_L(other.m_L), m_intensity(other.m_intensity),
-      m_sigmaIntensity(other.m_sigmaIntensity), m_binCount(other.m_binCount),
+    : convention(other.convention), m_bankName(other.m_bankName),
+      m_H(other.m_H), m_K(other.m_K), m_L(other.m_L),
+      m_intensity(other.m_intensity), m_sigmaIntensity(other.m_sigmaIntensity),
+      m_binCount(other.m_binCount),
       m_absorptionWeightedPathLength(other.m_absorptionWeightedPathLength),
       m_GoniometerMatrix(other.m_GoniometerMatrix),
       m_InverseGoniometerMatrix(other.m_InverseGoniometerMatrix),
@@ -86,6 +90,8 @@ BasePeak::BasePeak(const Geometry::IPeak &ipeak)
       m_col(ipeak.getCol()), m_peakNumber(ipeak.getPeakNumber()),
       m_intHKL(ipeak.getIntHKL()), m_intMNP(ipeak.getIntMNP()),
       m_peakShape(std::make_shared<NoShape>()) {
+  convention = Kernel::ConfigService::Instance().getString("Q.convention");
+
   if (fabs(m_InverseGoniometerMatrix.Invert()) < 1e-8)
     throw std::invalid_argument(
         "Peak::ctor(): Goniometer matrix must non-singular.");
@@ -406,6 +412,45 @@ double BasePeak::getAbsorptionWeightedPathLength() const {
   return m_absorptionWeightedPathLength;
 }
 
+double BasePeak::calculateWavelengthFromQLab(const V3D qLab) {
+  /* The q-vector direction of the peak is = goniometer * ub * hkl_vector
+   * The incident neutron wavevector is along the beam direction, ki = 1/wl
+   * (usually z, but referenceframe is definitive).
+   * In the inelastic convention, q = ki - kf.
+   * The final neutron wavector kf = -qx in x; -qy in y; and (-q.beam_dir+1/wl)
+   * in beam direction.
+   * AND: norm(kf) = norm(ki) = 2*pi/wavelength
+   * THEREFORE: 1/wl = norm(q)^2 / (2*q.beam_dir)
+   */
+  const double norm_q = qLab.norm();
+  if (norm_q == 0.0)
+    throw std::invalid_argument("BasePeak::setQLabFrame(): Q cannot be 0,0,0.");
+
+  std::shared_ptr<const ReferenceFrame> refFrame = getReferenceFrame();
+  V3D refBeamDir(0, 0, 1); // default beam direction +Z
+  if (refFrame)
+    refBeamDir = refFrame->vecPointingAlongBeam();
+
+  // Default for ki-kf has -q
+  const double qSign = (convention != "Crystallography") ? 1.0 : -1.0;
+  const double qBeam = qLab.scalar_prod(refBeamDir) * qSign;
+
+  if (qBeam == 0.0)
+    throw std::invalid_argument(
+        "BasePeak::setQLabFrame(): Q cannot be 0 in the beam direction.");
+
+  const double one_over_wl = (norm_q * norm_q) / (2.0 * qBeam);
+  const double wl = (2.0 * M_PI) / one_over_wl;
+  if (wl < 0.0) {
+    std::ostringstream mess;
+    mess << "BasePeak::setQLabFrame(): Wavelength found was negative (" << wl
+         << " Ang)! This Q is not physical.";
+    throw std::invalid_argument(mess.str());
+  }
+
+  return wl;
+}
+
 Mantid::Kernel::Logger BasePeak::g_log("PeakLogger");
 
 } // namespace DataObjects

Framework/DataObjects/src/LeanElasticPeak.cpp

@@ -6,8 +6,6 @@
 // SPDX - License - Identifier: GPL - 3.0 +
 #include "MantidDataObjects/LeanElasticPeak.h"
 #include "MantidDataObjects/NoShape.h"
-#include "MantidGeometry/Instrument/RectangularDetector.h"
-#include "MantidGeometry/Instrument/ReferenceFrame.h"
 #include "MantidGeometry/Objects/InstrumentRayTracer.h"
 #include "MantidGeometry/Surfaces/LineIntersectVisit.h"
 #include "MantidKernel/ConfigService.h"
@@ -42,9 +40,7 @@ LeanElasticPeak::LeanElasticPeak()
  *the sample frame (goniometer rotation accounted for).
  */
 LeanElasticPeak::LeanElasticPeak(const Mantid::Kernel::V3D &QSampleFrame)
-    : BasePeak() {
-  this->setQSampleFrame(QSampleFrame);
-}
+    : BasePeak(), m_Qsample(QSampleFrame), m_wavelength(0.) {}
 
 //----------------------------------------------------------------------------------------------
 /** Constructor that uses the Q position of the peak (in the lab frame).
@@ -52,11 +48,17 @@ LeanElasticPeak::LeanElasticPeak(const Mantid::Kernel::V3D &QSampleFrame)
  *
  * @param QSampleFrame :: Q of the center of the peak, in reciprocal space
  * @param goniometer :: a 3x3 rotation matrix
+ * @param refFrame :: optional reference frame, will default to beam along +Z
  */
 LeanElasticPeak::LeanElasticPeak(
     const Mantid::Kernel::V3D &QSampleFrame,
-    const Mantid::Kernel::Matrix<double> &goniometer)
-    : BasePeak(goniometer), m_Qsample(QSampleFrame) {}
+    const Mantid::Kernel::Matrix<double> &goniometer,
+    boost::optional<std::shared_ptr<ReferenceFrame>> refFrame)
+    : BasePeak() {
+  if (refFrame.is_initialized())
+    setReferenceFrame(refFrame.get());
+  setQSampleFrame(QSampleFrame, goniometer);
+}
 
 //----------------------------------------------------------------------------------------------
 /** Constructor that uses the Q position of the peak (in the sample frame)
@@ -71,26 +73,13 @@ LeanElasticPeak::LeanElasticPeak(const Mantid::Kernel::V3D &QSampleFrame,
                                  double wavelength)
     : BasePeak(), m_Qsample(QSampleFrame), m_wavelength(wavelength) {}
 
-//----------------------------------------------------------------------------------------------
-/** Constructor that uses the Q position of the peak (in the lab frame).
- * No detector ID is set.
- *
- * @param QSampleFrame :: Q of the center of the peak, in reciprocal space
- * @param goniometer :: a 3x3 rotation matrix
- * @param wavelength :: wavelength in Angstroms.
- */
-LeanElasticPeak::LeanElasticPeak(
-    const Mantid::Kernel::V3D &QSampleFrame,
-    const Mantid::Kernel::Matrix<double> &goniometer, double wavelength)
-    : BasePeak(goniometer), m_Qsample(QSampleFrame), m_wavelength(wavelength) {}
-
 /**
  * @brief Copy constructor
  * @param other : Source
  */
 LeanElasticPeak::LeanElasticPeak(const LeanElasticPeak &other)
     : BasePeak(other), m_Qsample(other.m_Qsample),
-      m_wavelength(other.m_wavelength) {}
+      m_wavelength(other.m_wavelength), m_refFrame(other.m_refFrame) {}
 
 //----------------------------------------------------------------------------------------------
 /** Constructor making a LeanElasticPeak from IPeak interface
@@ -149,6 +138,24 @@ Geometry::Instrument_const_sptr LeanElasticPeak::getInstrument() const {
       "LeanElasticPeak::setInstrument(): Has no instrument");
 }
 
+/**
+Getter for the reference frame.
+@return : reference frame.
+*/
+
+std::shared_ptr<const Geometry::ReferenceFrame>
+LeanElasticPeak::getReferenceFrame() const {
+  return m_refFrame;
+}
+
+/**
+Setter for the reference frame.
+@param frame : reference frame object to use.
+*/
+void LeanElasticPeak::setReferenceFrame(std::shared_ptr<ReferenceFrame> frame) {
+  m_refFrame = std::move(frame);
+}
+
 // -------------------------------------------------------------------------------------
 /** Return the neutron wavelength (in angstroms) */
 double LeanElasticPeak::getWavelength() const { return m_wavelength; }
@@ -212,6 +219,23 @@ void LeanElasticPeak::setQSampleFrame(const Mantid::Kernel::V3D &QSampleFrame,
   m_Qsample = QSampleFrame;
 }
 
+void LeanElasticPeak::setQSampleFrame(
+    const Mantid::Kernel::V3D &QSampleFrame,
+    const Mantid::Kernel::Matrix<double> &goniometer) {
+  m_Qsample = QSampleFrame;
+  setGoniometerMatrix(goniometer);
+
+  const V3D qLab = getQLabFrame();
+
+  try {
+    double wl = calculateWavelengthFromQLab(qLab);
+    setWavelength(wl);
+  } catch (std::exception &e) {
+    g_log.warning() << "Unable to determine wavelength from q-lab\n"
+                    << e.what() << '\n';
+  }
+}
+
 //----------------------------------------------------------------------------------------------
 /** Set the peak using the peak's position in reciprocal space, in the lab
  *frame.

Framework/DataObjects/src/Peak.cpp

@@ -31,9 +31,7 @@ namespace DataObjects {
 //----------------------------------------------------------------------------------------------
 /** Default constructor */
 Peak::Peak()
-    : BasePeak(), m_detectorID(-1), m_initialEnergy(0.), m_finalEnergy(0.) {
-  convention = Kernel::ConfigService::Instance().getString("Q.convention");
-}
+    : BasePeak(), m_detectorID(-1), m_initialEnergy(0.), m_finalEnergy(0.) {}
 
 //----------------------------------------------------------------------------------------------
 /** Constructor that uses the Q position of the peak (in the lab frame).
@@ -49,7 +47,6 @@ Peak::Peak(const Geometry::Instrument_const_sptr &m_inst,
            const Mantid::Kernel::V3D &QLabFrame,
            boost::optional<double> detectorDistance)
     : BasePeak() {
-  convention = Kernel::ConfigService::Instance().getString("Q.convention");
   this->setInstrument(m_inst);
   this->setQLabFrame(QLabFrame, std::move(detectorDistance));
 }
@@ -72,7 +69,6 @@ Peak::Peak(const Geometry::Instrument_const_sptr &m_inst,
            const Mantid::Kernel::Matrix<double> &goniometer,
            boost::optional<double> detectorDistance)
     : BasePeak(goniometer) {
-  convention = Kernel::ConfigService::Instance().getString("Q.convention");
   this->setInstrument(m_inst);
   this->setQSampleFrame(QSampleFrame, std::move(detectorDistance));
 }
@@ -88,7 +84,6 @@ Peak::Peak(const Geometry::Instrument_const_sptr &m_inst,
 Peak::Peak(const Geometry::Instrument_const_sptr &m_inst, int m_detectorID,
            double m_Wavelength)
     : BasePeak() {
-  convention = Kernel::ConfigService::Instance().getString("Q.convention");
   this->setInstrument(m_inst);
   this->setDetectorID(m_detectorID);
   this->setWavelength(m_Wavelength);
@@ -106,7 +101,6 @@ Peak::Peak(const Geometry::Instrument_const_sptr &m_inst, int m_detectorID,
 Peak::Peak(const Geometry::Instrument_const_sptr &m_inst, int m_detectorID,
            double m_Wavelength, const Mantid::Kernel::V3D &HKL)
     : BasePeak() {
-  convention = Kernel::ConfigService::Instance().getString("Q.convention");
   this->setInstrument(m_inst);
   this->setDetectorID(m_detectorID);
   this->setWavelength(m_Wavelength);
@@ -127,7 +121,6 @@ Peak::Peak(const Geometry::Instrument_const_sptr &m_inst, int m_detectorID,
            double m_Wavelength, const Mantid::Kernel::V3D &HKL,
            const Mantid::Kernel::Matrix<double> &goniometer)
     : BasePeak(goniometer) {
-  convention = Kernel::ConfigService::Instance().getString("Q.convention");
   this->setInstrument(m_inst);
   this->setDetectorID(m_detectorID);
   this->setWavelength(m_Wavelength);
@@ -144,7 +137,6 @@ Peak::Peak(const Geometry::Instrument_const_sptr &m_inst, int m_detectorID,
 Peak::Peak(const Geometry::Instrument_const_sptr &m_inst, double scattering,
            double m_Wavelength)
     : BasePeak() {
-  convention = Kernel::ConfigService::Instance().getString("Q.convention");
   this->setInstrument(m_inst);
   this->setWavelength(m_Wavelength);
   m_detectorID = -1;
@@ -163,7 +155,7 @@ Peak::Peak(const Peak &other)
       m_detectorID(other.m_detectorID), m_initialEnergy(other.m_initialEnergy),
       m_finalEnergy(other.m_finalEnergy), sourcePos(other.sourcePos),
       samplePos(other.samplePos), detPos(other.detPos),
-      m_detIDs(other.m_detIDs), convention(other.convention) {}
+      m_detIDs(other.m_detIDs) {}
 
 //----------------------------------------------------------------------------------------------
 /** Constructor making a Peak from IPeak interface
@@ -329,6 +321,12 @@ Geometry::IDetector_const_sptr Peak::getDetector() const { return m_det; }
 /** Return a shared ptr to the instrument for this peak. */
 Geometry::Instrument_const_sptr Peak::getInstrument() const { return m_inst; }
 
+/** Return a shared ptr to the instrument for this peak. */
+std::shared_ptr<const Geometry::ReferenceFrame>
+Peak::getReferenceFrame() const {
+  return m_inst->getReferenceFrame();
+}
+
 // -------------------------------------------------------------------------------------
 /** Calculate the neutron wavelength (in angstroms) at the peak
  * (Note for inelastic scattering - it is the wavelength corresponding to the
@@ -500,38 +498,14 @@ void Peak::setQLabFrame(const Mantid::Kernel::V3D &qLab,
   setCol(-1);
   setBankName("None");
 
-  /* The q-vector direction of the peak is = goniometer * ub * hkl_vector
-   * The incident neutron wavevector is along the beam direction, ki = 1/wl
-   * (usually z, but referenceframe is definitive).
-   * In the inelastic convention, q = ki - kf.
-   * The final neutron wavector kf = -qx in x; -qy in y; and (-q.beam_dir+1/wl)
-   * in beam direction.
-   * AND: norm(kf) = norm(ki) = 2*pi/wavelength
-   * THEREFORE: 1/wl = norm(q)^2 / (2*q.beam_dir)
-   */
-  const double norm_q = qLab.norm();
-  if (norm_q == 0.0)
-    throw std::invalid_argument("Peak::setQLabFrame(): Q cannot be 0,0,0.");
-
-  std::shared_ptr<const ReferenceFrame> refFrame =
-      this->m_inst->getReferenceFrame();
+  const double wl = calculateWavelengthFromQLab(qLab);
+
+  std::shared_ptr<const ReferenceFrame> refFrame = getReferenceFrame();
   const V3D refBeamDir = refFrame->vecPointingAlongBeam();
   // Default for ki-kf has -q
   const double qSign = (convention != "Crystallography") ? 1.0 : -1.0;
   const double qBeam = qLab.scalar_prod(refBeamDir) * qSign;
-
-  if (qBeam == 0.0)
-    throw std::invalid_argument(
-        "Peak::setQLabFrame(): Q cannot be 0 in the beam direction.");
-
-  const double one_over_wl = (norm_q * norm_q) / (2.0 * qBeam);
-  const double wl = (2.0 * M_PI) / one_over_wl;
-  if (wl < 0.0) {
-    std::ostringstream mess;
-    mess << "Peak::setQLabFrame(): Wavelength found was negative (" << wl
-         << " Ang)! This Q is not physical.";
-    throw std::invalid_argument(mess.str());
-  }
+  const double one_over_wl = (2.0 * M_PI) / wl;
 
   // Save the wavelength
   this->setWavelength(wl);

Framework/DataObjects/test/LeanElasticPeakTest.h

@@ -87,21 +87,12 @@ public:
 
     TS_ASSERT_EQUALS(p.getQSampleFrame(), V3D(1, 2, 3))
     TS_ASSERT_EQUALS(p.getQLabFrame(), V3D(2, 1, 3))
-  }
-
-  void test_Qsample_wavelength_constructor() {
-    LeanElasticPeak p(V3D(1, 2, 3), 1.);
-    TS_ASSERT_EQUALS(p.getQSampleFrame(), V3D(1, 2, 3))
-    TS_ASSERT_EQUALS(p.getQLabFrame(), V3D(1, 2, 3))
-
-    TS_ASSERT_DELTA(p.getInitialEnergy(), 81.8042024359, 1e-5)
-    TS_ASSERT_DELTA(p.getFinalEnergy(), 81.8042024359, 1e-5)
-    TS_ASSERT_DELTA(p.getWavelength(), 1., 1e-9)
+    TS_ASSERT_DELTA(p.getWavelength(), M_PI * 6 / 7, 1e-9)
     TS_ASSERT_DELTA(p.getDSpacing(), 1.679251908362714, 1e-9)
-    TS_ASSERT_DELTA(p.getScattering(), 0.6046731932, 1e-9)