Commit e953d1f0 authored by Janik Zikovsky's avatar Janik Zikovsky
Browse files

Fixes #3129: Peak constructors giving the scattered Q vector instead of the detector.

parent d68684e8
#ifndef MANTID_DATAOBJECTS_PEAK_H_ #ifndef MANTID_DATAOBJECTS_PEAK_H_
#define MANTID_DATAOBJECTS_PEAK_H_ #define MANTID_DATAOBJECTS_PEAK_H_
#include "MantidKernel/System.h"
#include "MantidGeometry/V3D.h"
#include "MantidGeometry/Math/Matrix.h"
#include "MantidGeometry/IInstrument.h" #include "MantidGeometry/IInstrument.h"
#include "MantidGeometry/Math/Matrix.h"
#include "MantidGeometry/V3D.h"
#include "MantidKernel/PhysicalConstants.h" #include "MantidKernel/PhysicalConstants.h"
#include "MantidKernel/System.h"
namespace Mantid namespace Mantid
...@@ -24,14 +24,18 @@ namespace DataObjects ...@@ -24,14 +24,18 @@ namespace DataObjects
/// Allow PeakColumn class to directly access members. /// Allow PeakColumn class to directly access members.
friend class PeakColumn; friend class PeakColumn;
Peak(Mantid::Geometry::IInstrument_sptr m_inst, Mantid::Geometry::V3D QSampleFrame, double detectorDistance=1.0);
Peak(Mantid::Geometry::IInstrument_sptr m_inst, Mantid::Geometry::V3D QSampleFrame, Mantid::Geometry::Matrix<double> goniometer, double detectorDistance=1.0);
Peak(Mantid::Geometry::IInstrument_sptr m_inst, int m_DetectorID, double m_Wavelength); Peak(Mantid::Geometry::IInstrument_sptr m_inst, int m_DetectorID, double m_Wavelength);
Peak(Mantid::Geometry::IInstrument_sptr m_inst, int m_DetectorID, double m_Wavelength, Mantid::Geometry::V3D HKL); Peak(Mantid::Geometry::IInstrument_sptr m_inst, int m_DetectorID, double m_Wavelength, Mantid::Geometry::V3D HKL);
Peak(Mantid::Geometry::IInstrument_sptr m_inst, int m_DetectorID, double m_Wavelength, Mantid::Geometry::V3D HKL,Mantid::Geometry::Matrix<double> goniometer); Peak(Mantid::Geometry::IInstrument_sptr m_inst, int m_DetectorID, double m_Wavelength, Mantid::Geometry::V3D HKL, Mantid::Geometry::Matrix<double> goniometer);
// Copy constructor is compiler-provided. // Copy constructor is compiler-provided.
// Peak(const Peak & other); // Peak(const Peak & other);
virtual ~Peak(); virtual ~Peak();
void setInstrument(Mantid::Geometry::IInstrument_sptr inst);
int getDetectorID() const; int getDetectorID() const;
void setDetectorID(int m_DetectorID); void setDetectorID(int m_DetectorID);
Mantid::Geometry::IDetector_sptr getDetector() const; Mantid::Geometry::IDetector_sptr getDetector() const;
...@@ -53,6 +57,9 @@ namespace DataObjects ...@@ -53,6 +57,9 @@ namespace DataObjects
Mantid::Geometry::V3D getQLabFrame() const; Mantid::Geometry::V3D getQLabFrame() const;
Mantid::Geometry::V3D getQSampleFrame() const; Mantid::Geometry::V3D getQSampleFrame() const;
void setQSampleFrame(Mantid::Geometry::V3D QSampleFrame, double detectorDistance=1.0);
void setQLabFrame(Mantid::Geometry::V3D QLabFrame, double detectorDistance=1.0);
void setWavelength(double wavelength); void setWavelength(double wavelength);
double getWavelength() const; double getWavelength() const;
double getDSpacing() const; double getDSpacing() const;
......
...@@ -12,6 +12,53 @@ namespace DataObjects ...@@ -12,6 +12,53 @@ namespace DataObjects
{ {
//----------------------------------------------------------------------------------------------
/** Constructor that uses the Q position of the peak (in the lab frame).
* No detector ID is set.
*
* @param m_inst :: Shared pointer to the instrument for this peak detection
* @param QLabFrame :: Q of the center of the peak, in reciprocal space
* @param detectorDistance :: distance between the sample and the detector.
* Used to give a valid TOF. Default 1.0 meters.
*/
Peak::Peak(Mantid::Geometry::IInstrument_sptr m_inst, Mantid::Geometry::V3D QLabFrame, double detectorDistance)
: m_H(0), m_K(0), m_L(0),
m_Intensity(0), m_SigmaIntensity(0), m_BinCount(0),
m_GoniometerMatrix(3,3,true),
m_InverseGoniometerMatrix(3,3,true),
m_RunNumber(0)
{
this->setInstrument(m_inst);
this->setQLabFrame(QLabFrame, detectorDistance);
}
//----------------------------------------------------------------------------------------------
/** Constructor that uses the Q position of the peak (in the sample frame)
* and a goniometer rotation matrix.
* No detector ID is set.
*
* @param m_inst :: Shared pointer to the instrument for this peak detection
* @param QLabFrame :: Q of the center of the peak, in reciprocal space
* @param goniometer :: a 3x3 rotation matrix
* @param detectorDistance :: distance between the sample and the detector.
* Used to give a valid TOF. Default 1.0 meters.
*/
Peak::Peak(Mantid::Geometry::IInstrument_sptr m_inst, Mantid::Geometry::V3D QSampleFrame,
Mantid::Geometry::Matrix<double> goniometer, double detectorDistance)
: m_H(0), m_K(0), m_L(0),
m_Intensity(0), m_SigmaIntensity(0), m_BinCount(0),
m_GoniometerMatrix(goniometer),
m_InverseGoniometerMatrix(goniometer),
m_RunNumber(0)
{
if(fabs(m_InverseGoniometerMatrix.Invert())<1e-8) throw std::invalid_argument("Goniometer matrix must non-singular.");
this->setInstrument(m_inst);
this->setQSampleFrame(QSampleFrame, detectorDistance);
}
//---------------------------------------------------------------------------------------------- //----------------------------------------------------------------------------------------------
/** Constructor /** Constructor
* *
...@@ -21,13 +68,13 @@ namespace DataObjects ...@@ -21,13 +68,13 @@ namespace DataObjects
* @return * @return
*/ */
Peak::Peak(Mantid::Geometry::IInstrument_sptr m_inst, int m_DetectorID, double m_Wavelength) Peak::Peak(Mantid::Geometry::IInstrument_sptr m_inst, int m_DetectorID, double m_Wavelength)
: m_inst(m_inst), : m_H(0), m_K(0), m_L(0),
m_H(0), m_K(0), m_L(0),
m_Intensity(0), m_SigmaIntensity(0), m_BinCount(0), m_Intensity(0), m_SigmaIntensity(0), m_BinCount(0),
m_GoniometerMatrix(3,3,true), m_GoniometerMatrix(3,3,true),
m_InverseGoniometerMatrix(3,3,true), m_InverseGoniometerMatrix(3,3,true),
m_RunNumber(0) m_RunNumber(0)
{ {
this->setInstrument(m_inst);
this->setDetectorID(m_DetectorID); this->setDetectorID(m_DetectorID);
this->setWavelength(m_Wavelength); this->setWavelength(m_Wavelength);
} }
...@@ -42,14 +89,14 @@ namespace DataObjects ...@@ -42,14 +89,14 @@ namespace DataObjects
* @param HKL :: vector with H,K,L position of the peak * @param HKL :: vector with H,K,L position of the peak
* @return * @return
*/ */
Peak::Peak(Mantid::Geometry::IInstrument_sptr m_inst, int m_DetectorID, double m_Wavelength, Mantid::Geometry::V3D HKL) : Peak::Peak(Mantid::Geometry::IInstrument_sptr m_inst, int m_DetectorID, double m_Wavelength, Mantid::Geometry::V3D HKL)
m_inst(m_inst), : m_H(HKL[0]), m_K(HKL[1]), m_L(HKL[2]),
m_H(HKL[0]), m_K(HKL[1]), m_L(HKL[2]),
m_Intensity(0), m_SigmaIntensity(0), m_BinCount(0), m_Intensity(0), m_SigmaIntensity(0), m_BinCount(0),
m_GoniometerMatrix(3,3,true), m_GoniometerMatrix(3,3,true),
m_InverseGoniometerMatrix(3,3,true), m_InverseGoniometerMatrix(3,3,true),
m_RunNumber(0) m_RunNumber(0)
{ {
this->setInstrument(m_inst);
this->setDetectorID(m_DetectorID); this->setDetectorID(m_DetectorID);
this->setWavelength(m_Wavelength); this->setWavelength(m_Wavelength);
} }
...@@ -64,8 +111,7 @@ namespace DataObjects ...@@ -64,8 +111,7 @@ namespace DataObjects
* @param goniometer :: a 3x3 rotation matrix * @param goniometer :: a 3x3 rotation matrix
* @return * @return
*/ */
Peak::Peak(Mantid::Geometry::IInstrument_sptr m_inst, int m_DetectorID, double m_Wavelength, Mantid::Geometry::V3D HKL,Mantid::Geometry::Matrix<double> goniometer) : Peak::Peak(Mantid::Geometry::IInstrument_sptr m_inst, int m_DetectorID, double m_Wavelength, Mantid::Geometry::V3D HKL, Mantid::Geometry::Matrix<double> goniometer) :
m_inst(m_inst),
m_H(HKL[0]), m_K(HKL[1]), m_L(HKL[2]), m_H(HKL[0]), m_K(HKL[1]), m_L(HKL[2]),
m_Intensity(0), m_SigmaIntensity(0), m_BinCount(0), m_Intensity(0), m_SigmaIntensity(0), m_BinCount(0),
m_GoniometerMatrix(goniometer), m_GoniometerMatrix(goniometer),
...@@ -73,6 +119,7 @@ namespace DataObjects ...@@ -73,6 +119,7 @@ namespace DataObjects
m_RunNumber(0) m_RunNumber(0)
{ {
if(fabs(m_InverseGoniometerMatrix.Invert())<1e-8) throw std::invalid_argument("Goniometer matrix must non-singular."); if(fabs(m_InverseGoniometerMatrix.Invert())<1e-8) throw std::invalid_argument("Goniometer matrix must non-singular.");
this->setInstrument(m_inst);
this->setDetectorID(m_DetectorID); this->setDetectorID(m_DetectorID);
this->setWavelength(m_Wavelength); this->setWavelength(m_Wavelength);
} }
...@@ -86,8 +133,10 @@ namespace DataObjects ...@@ -86,8 +133,10 @@ namespace DataObjects
{ {
} }
//---------------------------------------------------------------------------------------------- //----------------------------------------------------------------------------------------------
/** Set the incident wavelength of the neutron. Calculates the energy from this. Assumes elastic scattering. /** Set the incident wavelength of the neutron. Calculates the energy from this.
* Assumes elastic scattering.
* *
* @param wavelength :: wavelength in Angstroms. * @param wavelength :: wavelength in Angstroms.
*/ */
...@@ -102,17 +151,17 @@ namespace DataObjects ...@@ -102,17 +151,17 @@ namespace DataObjects
m_FinalEnergy = m_InitialEnergy; m_FinalEnergy = m_InitialEnergy;
} }
//---------------------------------------------------------------------------------------------- //----------------------------------------------------------------------------------------------
/** Set the detector ID and look up and cache values related to it. /** Set the instrument (and save the source/sample pos).
* Call setDetectorID AFTER this call.
* *
* @param m_DetectorID :: ID of detector at the center of the peak * @param inst :: Instrument sptr to use
*/ */
void Peak::setDetectorID(int m_DetectorID) void Peak::setInstrument(Mantid::Geometry::IInstrument_sptr inst)
{ {
this->m_DetectorID = m_DetectorID; m_inst = inst;
if (!m_inst) throw std::runtime_error("No instrument is set!"); if (!inst) throw std::runtime_error("No instrument is set!");
this->m_det = m_inst->getDetector(this->m_DetectorID);
if (!m_det) throw std::runtime_error("No detector was found!");
// Cache some positions // Cache some positions
const Geometry::IObjComponent_sptr sourceObj = m_inst->getSource(); const Geometry::IObjComponent_sptr sourceObj = m_inst->getSource();
...@@ -124,6 +173,20 @@ namespace DataObjects ...@@ -124,6 +173,20 @@ namespace DataObjects
sourcePos = sourceObj->getPos(); sourcePos = sourceObj->getPos();
samplePos = sampleObj->getPos(); samplePos = sampleObj->getPos();
}
//----------------------------------------------------------------------------------------------
/** Set the detector ID and look up and cache values related to it.
*
* @param m_DetectorID :: ID of detector at the center of the peak
*/
void Peak::setDetectorID(int m_DetectorID)
{
this->m_DetectorID = m_DetectorID;
if (!m_inst) throw std::runtime_error("No instrument is set!");
this->m_det = m_inst->getDetector(this->m_DetectorID);
if (!m_det) throw std::runtime_error("No detector was found!");
detPos = m_det->getPos(); detPos = m_det->getPos();
// We now look for the row/column. -1 if not found. // We now look for the row/column. -1 if not found.
...@@ -192,7 +255,8 @@ namespace DataObjects ...@@ -192,7 +255,8 @@ namespace DataObjects
// In general case (2*pi/d)^2=k_i^2+k_f^2-2*k_i*k_f*cos(two_theta) // In general case (2*pi/d)^2=k_i^2+k_f^2-2*k_i*k_f*cos(two_theta)
// E_i,f=k_i,f^2*hbar^2/(2 m) // E_i,f=k_i,f^2*hbar^2/(2 m)
return 1e10*PhysicalConstants::h/sqrt(2.0*PhysicalConstants::NeutronMass*PhysicalConstants::meV)/sqrt(m_InitialEnergy+m_FinalEnergy-2.0*sqrt(m_InitialEnergy*m_FinalEnergy)*cos(two_theta)); return 1e10*PhysicalConstants::h/sqrt(2.0*PhysicalConstants::NeutronMass*PhysicalConstants::meV)
/ sqrt(m_InitialEnergy+m_FinalEnergy-2.0*sqrt(m_InitialEnergy*m_FinalEnergy)*cos(two_theta));
} }
//---------------------------------------------------------------------------------------------- //----------------------------------------------------------------------------------------------
...@@ -237,6 +301,76 @@ namespace DataObjects ...@@ -237,6 +301,76 @@ namespace DataObjects
return Qsample; return Qsample;
} }
//----------------------------------------------------------------------------------------------
/** Set the peak using the peak's position in reciprocal space, in the sample frame.
* The GoniometerMatrix will be used to find the Q in the lab frame, so it should
* be set beforehand.
*
* @param QSampleFrame :: Q of the center of the peak, in reciprocal space
* This is in inelastic convention: momentum transfer of the LATTICE!
* Also, q does NOT have a 2pi factor = it is equal to 1/wavelength.
* @param detectorDistance :: distance between the sample and the detector.
* Used to give a valid TOF. Default 1.0 meters.
*/
void Peak::setQSampleFrame(Mantid::Geometry::V3D QSampleFrame, double detectorDistance)
{
V3D Qlab = m_GoniometerMatrix * QSampleFrame;
this->setQLabFrame(Qlab, detectorDistance);
}
//----------------------------------------------------------------------------------------------
/** Set the peak using the peak's position in reciprocal space, in the lab frame.
* The detector position will be determined.
* DetectorID, row and column will be set to -1 since they are not (necessarily) found.
*
* @param QLabFrame :: Q of the center of the peak, in reciprocal space.
* This is in inelastic convention: momentum transfer of the LATTICE!
* Also, q does NOT have a 2pi factor = it is equal to 1/wavelength (in Angstroms).
* @param detectorDistance :: distance between the sample and the detector.
* Used to give a valid TOF. Default 1.0 meters.
*/
void Peak::setQLabFrame(Mantid::Geometry::V3D QLabFrame, double detectorDistance)
{
// Clear out the detector = we can't know them
m_DetectorID = -1;
m_det = IDetector_sptr();
m_Row = -1;
m_Col = -1;
m_BankName = "None";
// The q-vector direction of the peak is = goniometer * ub * hkl_vector
V3D q = QLabFrame;
/* The incident neutron wavevector is in the +Z direction, ki = 1/wl (in z direction).
* In the inelastic convention, q = ki - kf.
* The final neutron wavector kf = -qx in x; -qy in y; and (-qz+1/wl) in z.
* AND: norm(kf) = norm(ki) = 1.0/wavelength
* THEREFORE: 1/wl = norm(q)^2 / (2*qz)
*/
double norm_q = q.norm();
if (norm_q == 0.0)
throw std::invalid_argument("Peak::setQLabFrame(): Q cannot be 0,0,0.");
if (q.Z() == 0.0)
throw std::invalid_argument("Peak::setQLabFrame(): Q cannot be 0 in the Z (beam) direction.");
double one_over_wl = (norm_q*norm_q) / (2.0 * q.Z());
double wl = 1.0/one_over_wl;
// This is the scattered direction, kf = (-qx, -qy, 1/wl-qz)
V3D beam = q * -1.0;
beam.setZ(one_over_wl - q.Z());
beam.normalize();
// Save the wavelength
this->setWavelength(wl);
// Use the given detector distance to find the detector position.
detPos = beam * detectorDistance;
}
//---------------------------------------------------------------------------------------------- //----------------------------------------------------------------------------------------------
......
...@@ -172,7 +172,6 @@ public: ...@@ -172,7 +172,6 @@ public:
void test_getQSampleFrame() void test_getQSampleFrame()
{ {
// Peak 3 is phi,chi,omega of 90,0,0; giving this matrix: // Peak 3 is phi,chi,omega of 90,0,0; giving this matrix:
Matrix<double> r2(3,3,false); Matrix<double> r2(3,3,false);
r2[0][2] = 1; r2[0][2] = 1;
...@@ -193,6 +192,74 @@ public: ...@@ -193,6 +192,74 @@ public:
TS_ASSERT_EQUALS(qLab, qSampleRotated); TS_ASSERT_EQUALS(qLab, qSampleRotated);
} }
//------------------------------------------------------------------------------------
/** Can't have Q = 0,0,0 or 0 in the Z direction when creating */
void test_setQLabFrame_ThrowsIfQIsNull()
{
Peak p1(inst, 10000, 2.0);
TS_ASSERT_THROWS_ANYTHING(Peak p2(inst, V3D(0,0,0), 1.0));
TS_ASSERT_THROWS_ANYTHING(Peak p2(inst, V3D(1,2,0), 1.0));
}
/** Compare two peaks, but not the detector IDs etc. */
void comparePeaks(Peak & p1, Peak & p2)
{
TS_ASSERT_EQUALS( p1.getQLabFrame(), p2.getQLabFrame() );
TS_ASSERT_EQUALS( p1.getQSampleFrame(), p2.getQSampleFrame() );
TS_ASSERT_EQUALS( p1.getDetPos(), p2.getDetPos() );
TS_ASSERT_EQUALS( p1.getHKL(), p2.getHKL() );
TS_ASSERT_DELTA( p1.getWavelength(), p2.getWavelength(), 1e-5 );
TS_ASSERT_DELTA( p1.getL1(), p2.getL1(), 1e-5 );
TS_ASSERT_DELTA( p1.getL2(), p2.getL2(), 1e-5 );
TS_ASSERT_DELTA( p1.getTOF(), p2.getTOF(), 1e-5 );
TS_ASSERT_DELTA( p1.getInitialEnergy(), p2.getInitialEnergy(), 1e-5 );
TS_ASSERT_DELTA( p1.getFinalEnergy(), p2.getFinalEnergy(), 1e-5 );
TS_ASSERT( p1.getGoniometerMatrix().equals(p2.getGoniometerMatrix(), 1e-5) );
}
/** Create peaks using Q in the lab frame */
void test_setQLabFrame()
{
Peak p1(inst, 19999, 2.0);
V3D Qlab1 = p1.getQLabFrame();
V3D detPos1 = p1.getDetPos();
// Construct using just Q
Peak p2(inst, Qlab1, detPos1.norm());
comparePeaks(p1, p2);
TS_ASSERT_EQUALS( p2.getBankName(), "None");
TS_ASSERT_EQUALS( p2.getRow(), -1);
TS_ASSERT_EQUALS( p2.getCol(), -1);
TS_ASSERT_EQUALS( p2.getDetectorID(), -1);
}
/** Create peaks using Q in sample frame + a goniometer rotation matrix*/
void test_setQSampleFrame()
{
// A goniometer rotation matrix
Matrix<double> r2(3,3,false);
r2[0][2] = 1;
r2[1][1] = 1;
r2[2][0] = -1;
Peak p1(inst, 19999, 2.0, V3D(1,2,3), r2);
V3D q = p1.getQSampleFrame();
V3D detPos1 = p1.getDetPos();
// Construct using Q + rotation matrix
Peak p2(inst, q, r2, detPos1.norm());
p2.setHKL(V3D(1,2,3)); // Make sure HKL matches too.
comparePeaks(p1, p2);
TS_ASSERT_EQUALS( p2.getBankName(), "None");
TS_ASSERT_EQUALS( p2.getRow(), -1);
TS_ASSERT_EQUALS( p2.getCol(), -1);
TS_ASSERT_EQUALS( p2.getDetectorID(), -1);
}
}; };
......
Supports Markdown
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