diff --git a/Code/Mantid/Framework/DataObjects/inc/MantidDataObjects/Peak.h b/Code/Mantid/Framework/DataObjects/inc/MantidDataObjects/Peak.h index ab55a5f3ae97fbe8d0f26869f834c8f6b291cc5d..38a1b5407ce6d5f7a4b1630585183866a36b687d 100644 --- a/Code/Mantid/Framework/DataObjects/inc/MantidDataObjects/Peak.h +++ b/Code/Mantid/Framework/DataObjects/inc/MantidDataObjects/Peak.h @@ -4,6 +4,7 @@ #include "MantidKernel/System.h" #include "MantidGeometry/Math/Matrix.h" #include "MantidGeometry/IInstrument.h" +#include "MantidKernel/PhysicalConstants.h" namespace Mantid @@ -20,42 +21,52 @@ namespace DataObjects { public: Peak(Mantid::Geometry::IInstrument_sptr m_inst, int m_DetectorID, double m_InitialEnergy); - ~Peak(); + + // Copy constructor is compiler-provided. + // Peak(const Peak & other); + virtual ~Peak(); int getDetectorID() const; + void setDetectorID(int m_DetectorID); + Mantid::Geometry::IDetector_sptr getDetector() const; + Mantid::Geometry::IInstrument_sptr getInstrument() const; + + int getRunNumber() const; + void setRunNumber(int m_RunNumber); double getH() const; double getK() const; double getL() const; - - double getInitialEnergy() const; - double getFinalEnergy() const; - - double getIntensity() const; - double getSigIntensity() const; - - void setDetectorID(int m_DetectorID); - void setH(double m_H); void setK(double m_K); void setL(double m_L); void setHKL(double H, double K, double L); + Mantid::Geometry::V3D getQLabFrame() const; + Mantid::Geometry::V3D getQSampleFrame() const; + + void setWavelength(double wavelength); + double getWavelength() const; + double getDSpacing() const; + double getTOF() const; + + double getInitialEnergy() const; + double getFinalEnergy() const; void setInitialEnergy(double m_InitialEnergy); void setFinalEnergy(double m_FinalEnergy); + double getIntensity() const; + double getSigmaIntensity() const; + void setIntensity(double m_Intensity); - void setSigIntensity(double m_SigIntensity); + void setSigmaIntensity(double m_SigmaIntensity); - Mantid::Geometry::Matrix<double> getOrientationMatrix() const; - void setOrientationMatrix(Mantid::Geometry::Matrix<double> m_OrientationMatrix); + Mantid::Geometry::Matrix<double> getGoniometerMatrix() const; + void setGoniometerMatrix(Mantid::Geometry::Matrix<double> m_GoniometerMatrix); - // ------ Still to implement ------- - double getWavelength(); - double getDSpacing(); - std::string getBankName(); - int getRow(); - int getCol(); + std::string getBankName() const; + int getRow() const; + int getCol() const; protected: /// Shared pointer to the instrument (for calculating some values ) @@ -80,7 +91,7 @@ namespace DataObjects double m_Intensity; /// Error (sigma) on peak intensity - double m_SigIntensity; + double m_SigmaIntensity; /// Initial energy of neutrons at the peak double m_InitialEnergy; @@ -88,8 +99,25 @@ namespace DataObjects /// Final energy of the neutrons at peak (normally same as m_InitialEnergy) double m_FinalEnergy; - /// Sample orientation matrix - Mantid::Geometry::Matrix<double> m_OrientationMatrix; + /// Orientation matrix of the goniometer angles. + Mantid::Geometry::Matrix<double> m_GoniometerMatrix; + + /// Originating run number for this peak + int m_RunNumber; + + /// Cached row in the detector + int m_Row; + + /// Cached column in the detector + int m_Col; + + /// Cached source position + Mantid::Geometry::V3D sourcePos; + /// Cached sample position + Mantid::Geometry::V3D samplePos; + /// Cached detector position + Mantid::Geometry::V3D detPos; + }; diff --git a/Code/Mantid/Framework/DataObjects/src/Peak.cpp b/Code/Mantid/Framework/DataObjects/src/Peak.cpp index c02eb2f9370a09c796617d17f0e075e132b49d23..5e09db1f19fe9aec7b2b432ad661cfd8c103f52a 100644 --- a/Code/Mantid/Framework/DataObjects/src/Peak.cpp +++ b/Code/Mantid/Framework/DataObjects/src/Peak.cpp @@ -17,20 +17,28 @@ namespace DataObjects * * @param m_inst :: Shared pointer to the instrument for this peak detection * @param m_DetectorID :: ID to the detector of the center of the peak - * @param m_InitialEnergy :: incident neutron energy, in mEv (UNITS??) + * @param m_Wavelength :: incident neutron wavelength, in Angstroms * @return */ - Peak::Peak(Mantid::Geometry::IInstrument_sptr m_inst, int m_DetectorID, double m_InitialEnergy) + 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_Intensity(0), m_SigIntensity(0), - m_InitialEnergy(m_InitialEnergy), - m_FinalEnergy(m_InitialEnergy), - m_OrientationMatrix(3,3) + m_Intensity(0), m_SigmaIntensity(0), + m_GoniometerMatrix(3,3), + m_RunNumber(0) { this->setDetectorID(m_DetectorID); + this->setWavelength(m_Wavelength); } +// /** Copy constructor +// * @param other :: Peak to copy. +// */ +// Peak::Peak(const Peak & other) +// : m_inst(other. +// { +// } + //---------------------------------------------------------------------------------------------- /** Destructor */ @@ -38,98 +46,203 @@ namespace DataObjects { } - int Peak::getDetectorID() const + //---------------------------------------------------------------------------------------------- + /** Set the incident wavelength of the neutron. Calculates the energy from this. + * + * @param wavelength :: wavelength in Angstroms. + */ + void Peak::setWavelength(double wavelength) { - return m_DetectorID; + // Velocity of the neutron (non-relativistic) + double velocity = PhysicalConstants::h / (wavelength * 1e-10 * PhysicalConstants::NeutronMass ); + // Energy in J of the neutron + double energy = PhysicalConstants::NeutronMass * velocity * velocity / 2; + // Convert to meV + m_InitialEnergy = energy / PhysicalConstants::meV; + m_FinalEnergy = m_InitialEnergy; } - double Peak::getFinalEnergy() const + //---------------------------------------------------------------------------------------------- + /** 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) { - return m_FinalEnergy; + 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!"); + + + // Cache some positions + const Geometry::IObjComponent_sptr sourceObj = m_inst->getSource(); + if (sourceObj == NULL) + throw Exception::InstrumentDefinitionError("Failed to get source component from instrument"); + const Geometry::IObjComponent_sptr sampleObj = m_inst->getSample(); + if (sampleObj == NULL) + throw Exception::InstrumentDefinitionError("Failed to get sample component from instrument"); + + sourcePos = sourceObj->getPos(); + samplePos = sampleObj->getPos(); + detPos = m_det->getPos(); + + // We now look for the row/column. -1 if not found. + m_Row = -1; + m_Col = -1; + if (!m_det->getParent()) return; + RectangularDetector_const_sptr retDet = boost::dynamic_pointer_cast<const RectangularDetector>(m_det->getParent()); + if (!retDet) return; + std::pair<int,int> xy = retDet->getXYForDetectorID(m_DetectorID); + m_Row = xy.second; + m_Col = xy.first; } - double Peak::getH() const - { - return m_H; - } - double Peak::getInitialEnergy() const + // ------------------------------------------------------------------------------------- + /** Calculate the neutron wavelength (in angstroms) at the peak */ + double Peak::getWavelength() const { - return m_InitialEnergy; + // Energy in J of the neutron + double energy = PhysicalConstants::meV * m_InitialEnergy; + // v = sqrt(E / 2*m) + double velocity = sqrt(energy/(2*PhysicalConstants::NeutronMass)); + // wavelength = h / mv + double wavelength = PhysicalConstants::h / (PhysicalConstants::NeutronMass * velocity); + // Return it in angstroms + return wavelength * 1e10; } - double Peak::getIntensity() const + // ------------------------------------------------------------------------------------- + /** Calculate the time of flight (in microseconds) of the neutrons for this peak, + * using the geometry of the detector */ + double Peak::getTOF() const { - return m_Intensity; + // First, calculate the total neutron traveling distance + double distance = (detPos - samplePos).norm() + (samplePos - sourcePos).norm(); + // Energy in J of the neutron + double energy = PhysicalConstants::meV * m_InitialEnergy; + // v = sqrt(E / 2*m) + double velocity = sqrt(energy/(2*PhysicalConstants::NeutronMass)); + // Time of flight in seconds = distance / speed + double tof = distance / velocity; + // Return in microsecond units + return tof * 1e6; } - double Peak::getK() const + // ------------------------------------------------------------------------------------- + /** Calculate the d-spacing of the peak, in 1/Angstroms */ + double Peak::getDSpacing() const { - return m_K; + // Get the neutron wavelength in angstroms + double wavelength = this->getWavelength(); + // The detector is at 2 theta scattering angle + V3D beamDir = samplePos - sourcePos; + double two_theta = detPos.angle(beamDir); + double sin_theta = sin(two_theta/2.0); + // Bragg condition is n*wavelength = 2 * d * sin(theta) + return wavelength / (2 * sin_theta); } - double Peak::getL() const + //---------------------------------------------------------------------------------------------- + /** Return the Q change (of the lattice, k_i - k_f) for this peak. + * The Q is in the Lab frame: the goniometer rotation was NOT taken out. + * + * Note: There is no 2*pi factor used, so |Q| = 1/wavelength. + * */ + Mantid::Geometry::V3D Peak::getQLabFrame() const { - return m_L; + // Normalized beam direction + V3D beamDir = samplePos - sourcePos; + beamDir /= beamDir.norm(); + // Normalized detector direction + V3D detDir = (detPos - samplePos); + detDir /= detDir.norm(); + // Normalized Q direction + V3D Q_dir = beamDir - detDir; + // Now calculate the wavevector of the incident neutron + double wavevector = 1.0 / this->getWavelength(); + // And Q in the lab frame = this direction * the wave vector. + return Q_dir * wavevector; } - double Peak::getSigIntensity() const + //---------------------------------------------------------------------------------------------- + /** Return the Q change (of the lattice, k_i - k_f) for this peak. + * The Q is in the Sample frame: the goniometer rotation WAS taken out. */ + Mantid::Geometry::V3D Peak::getQSampleFrame() const { - return m_SigIntensity; + V3D Qlab = this->getQLabFrame(); + // Multiply by the goniometer matrix to get the sample frame + // TODO: should I multiply by the inverse??? + V3D Qsample = m_GoniometerMatrix * Qlab; + return Qsample; } - void Peak::setDetectorID(int m_DetectorID) - { - this->m_DetectorID = m_DetectorID; - if (m_inst) - this->m_det = m_inst->getDetector(this->m_DetectorID); - } - void Peak::setFinalEnergy(double m_FinalEnergy) - { - this->m_FinalEnergy = m_FinalEnergy; - } - void Peak::setH(double m_H) - { - this->m_H = m_H; - } + //---------------------------------------------------------------------------------------------- + /** Return a shared ptr to the detector at center of peak. */ + Mantid::Geometry::IDetector_sptr Peak::getDetector() const + { return m_det; } - void Peak::setInitialEnergy(double m_InitialEnergy) - { - this->m_InitialEnergy = m_InitialEnergy; - } + /** Return a shared ptr to the instrument for this peak. */ + Mantid::Geometry::IInstrument_sptr Peak::getInstrument() const + { return m_inst; } - void Peak::setIntensity(double m_Intensity) - { - this->m_Intensity = m_Intensity; - } + //---------------------------------------------------------------------------------------------- + /** Return the run number this peak was measured at. */ + int Peak::getRunNumber() const + { return m_RunNumber; } - void Peak::setK(double m_K) - { - this->m_K = m_K; - } + /** Set the run number that measured this peak + * @param m_RunNumber :: the run number */ + void Peak::setRunNumber(int m_RunNumber) + { this->m_RunNumber = m_RunNumber; } - void Peak::setL(double m_L) - { - this->m_L = m_L; - } + //---------------------------------------------------------------------------------------------- + /** Get the ID of the detector at the center of the peak */ + int Peak::getDetectorID() const + { return m_DetectorID; } - void Peak::setSigIntensity(double m_SigIntensity) - { - this->m_SigIntensity = m_SigIntensity; - } + //---------------------------------------------------------------------------------------------- + /** Get the final neutron energy */ + double Peak::getFinalEnergy() const + { return m_FinalEnergy; } - Mantid::Geometry::Matrix<double> Peak::getOrientationMatrix() const - { - return this->m_OrientationMatrix; - } + /** Get the initial (incident) neutron energy */ + double Peak::getInitialEnergy() const + { return m_InitialEnergy; } - void Peak::setOrientationMatrix(Mantid::Geometry::Matrix<double> m_OrientationMatrix) - { - this->m_OrientationMatrix = m_OrientationMatrix; - } + //---------------------------------------------------------------------------------------------- + /** Get the H index of the peak */ + double Peak::getH() const + { return m_H; } + + /** Get the K index of the peak */ + double Peak::getK() const + { return m_K; } + + /** Get the L index of the peak */ + double Peak::getL() const + { return m_L; } + + //---------------------------------------------------------------------------------------------- + /** Set the H index of this peak + * @param m_H :: index to set */ + void Peak::setH(double m_H) + { this->m_H = m_H; } + + /** Set the K index of this peak + * @param m_K :: index to set */ + void Peak::setK(double m_K) + { this->m_K = m_K; } + + /** Set the L index of this peak + * @param m_L :: index to set */ + void Peak::setL(double m_L) + { this->m_L = m_L; } + /** Set all three H,K,L indices of the peak */ void Peak::setHKL(double H, double K, double L) { m_H = H; @@ -137,26 +250,56 @@ namespace DataObjects m_L = L; } + //---------------------------------------------------------------------------------------------- + /** Return the integrated peak intensity */ + double Peak::getIntensity() const + { return m_Intensity; } + + /** Return the error on the integrated peak intensity */ + double Peak::getSigmaIntensity() const + { return m_SigmaIntensity; } + + /** Set the integrated peak intensity + * @param m_Intensity :: intensity value */ + void Peak::setIntensity(double m_Intensity) + { this->m_Intensity = m_Intensity; } + + /** Set the error on the integrated peak intensity + * @param m_Intensity :: intensity error value */ + void Peak::setSigmaIntensity(double m_SigmaIntensity) + { this->m_SigmaIntensity = m_SigmaIntensity; } + + + void Peak::setFinalEnergy(double m_FinalEnergy) + { this->m_FinalEnergy = m_FinalEnergy; } + + void Peak::setInitialEnergy(double m_InitialEnergy) + { this->m_InitialEnergy = m_InitialEnergy; } + + // ------------------------------------------------------------------------------------- - /** Calculate the neutron wavelength at the peak */ - double Peak::getWavelength() + /** Get the goniometer rotation matrix at which this peak was measured. */ + Mantid::Geometry::Matrix<double> Peak::getGoniometerMatrix() const { - return 0.0; //TODO: + return this->m_GoniometerMatrix; } - // ------------------------------------------------------------------------------------- - /** Calculate the d-spacing of the peak */ - double Peak::getDSpacing() + /** Set the goniometer rotation matrix at which this peak was measured. + * @throw std::invalid_argument if matrix is not 3x3*/ + void Peak::setGoniometerMatrix(Mantid::Geometry::Matrix<double> goniometerMatrix) { - return 0.0; //TODO: + if ((goniometerMatrix.numCols() != 3) || (goniometerMatrix.numRows() != 3)) + throw std::invalid_argument("Goniometer matrix must be 3x3."); + this->m_GoniometerMatrix = goniometerMatrix; } + // ------------------------------------------------------------------------------------- /** Find the name of the bank that is the parent of the detector. This works * best for RectangularDetector instruments (goes up one level) * @return name of the bank. */ - std::string Peak::getBankName() + std::string Peak::getBankName() const { if (!m_det) return ""; if (!m_det->getParent()) return ""; @@ -164,29 +307,16 @@ namespace DataObjects } // ------------------------------------------------------------------------------------- - /** For RectangularDetectors only, returns the row (y) of the pixel of the detector. */ - int Peak::getRow() - { - if (!m_det) throw std::runtime_error("No detector in Peak!"); - if (!m_det->getParent()) throw std::runtime_error("No parent to the detector in Peak!"); - RectangularDetector_const_sptr retDet = boost::dynamic_pointer_cast<const RectangularDetector>(m_det->getParent()); - if (!retDet) throw std::runtime_error("Parent of the detector in Peak is not RectangularDetector: can't find the row!"); - std::pair<int,int> xy = retDet->getXYForDetectorID(m_DetectorID); - return xy.second; - } + /** For RectangularDetectors only, returns the row (y) of the pixel of the detector. + * Returns -1 if it could not find it. */ + int Peak::getRow() const + { return m_Row; } // ------------------------------------------------------------------------------------- - /** For RectangularDetectors only, returns the row of the pixel of the detector. */ - int Peak::getCol() - { - if (!m_det) throw std::runtime_error("No detector in Peak!"); - if (!m_det->getParent()) throw std::runtime_error("No parent to the detector in Peak!"); - RectangularDetector_const_sptr retDet = boost::dynamic_pointer_cast<const RectangularDetector>(m_det->getParent()); - if (!retDet) throw std::runtime_error("Parent of the detector in Peak is not RectangularDetector: can't find the col!"); - std::pair<int,int> xy = retDet->getXYForDetectorID(m_DetectorID); - return xy.first; - } - + /** For RectangularDetectors only, returns the column (x) of the pixel of the detector. + * Returns -1 if it could not find it. */ + int Peak::getCol() const + { return m_Col; } diff --git a/Code/Mantid/Framework/DataObjects/test/PeakTest.h b/Code/Mantid/Framework/DataObjects/test/PeakTest.h index d7c4bbb5dddc2ef4539ae70b0fcc15563be6c4c6..c0071060f0169e61ceb764b31eb04f31e0e9c376 100644 --- a/Code/Mantid/Framework/DataObjects/test/PeakTest.h +++ b/Code/Mantid/Framework/DataObjects/test/PeakTest.h @@ -28,12 +28,43 @@ public: { // detector IDs start at 10000 Peak p(inst, 10000, 2.0); - TS_ASSERT_DELTA(p.getInitialEnergy(), 2.0, 1e-5) - TS_ASSERT_DELTA(p.getFinalEnergy(), 2.0, 1e-5) TS_ASSERT_DELTA(p.getH(), 0.0, 1e-5) TS_ASSERT_DELTA(p.getK(), 0.0, 1e-5) TS_ASSERT_DELTA(p.getL(), 0.0, 1e-5) TS_ASSERT_EQUALS(p.getDetectorID(), 10000) + TS_ASSERT_EQUALS(p.getDetector()->getID(), 10000) + TS_ASSERT_EQUALS(p.getInstrument(), inst) + } + + void test_copyConstructor() + { + Peak p(inst, 10102, 2.0); + p.setHKL(1,2,3); + p.setRunNumber(1234); + // Default (not-explicit) copy constructor + Peak p2(p); + TS_ASSERT_EQUALS(p.getRow(), p2.getRow()); + TS_ASSERT_EQUALS(p.getCol(), p2.getCol()); + TS_ASSERT_EQUALS(p.getH(), p2.getH()); + TS_ASSERT_EQUALS(p.getK(), p2.getK()); + TS_ASSERT_EQUALS(p.getL(), p2.getL()); + TS_ASSERT_EQUALS(p.getGoniometerMatrix(), p2.getGoniometerMatrix()); + TS_ASSERT_EQUALS(p.getRunNumber(), p2.getRunNumber()); + TS_ASSERT_EQUALS(p.getDetector(), p2.getDetector()) + TS_ASSERT_EQUALS(p.getInstrument(), p2.getInstrument()) + } + + /** Set the wavelength and see the other "versions" of it get calculated. */ + void test_wavelength_conversion() + { + //1 angstroms wavelength, and at the opposite corner of the detector + Peak p(inst, 19999, 1.0); + // Energy in meV + TS_ASSERT_DELTA(p.getInitialEnergy(), 81.805, 1e-3) // Conversion table at : www.ncnr.nist.gov/instruments/dcs/dcs_usersguide/Conversion_Factors.pdf + TS_ASSERT_DELTA(p.getFinalEnergy(), p.getInitialEnergy(), 1e-5) + // TODO: Check that the conversion is correct (I just took the results and put them back into the test, they are within a reasonable range though) + TS_ASSERT_DELTA(p.getDSpacing(), 9.0938, 1e-3); + TS_ASSERT_DELTA(p.getTOF(), 7645.996, 1e-2); } void test_badDetectorID_throws() @@ -42,15 +73,57 @@ public: TS_ASSERT_THROWS_ANYTHING( p.setDetectorID(7) ); } + void test_runNumber() + { + Peak p(inst, 10000, 2.0); + p.setRunNumber(12345); + TS_ASSERT_EQUALS( p.getRunNumber(), 12345); + } + + void test_GoniometerMatrix() + { + Peak p(inst, 10000, 2.0); + Matrix<double> mat(3,3); + for (int x=0; x<3; x++) + for (int y=0; y<3; y++) + mat[x][y]=x+y; + p.setGoniometerMatrix(mat); + TS_ASSERT_EQUALS( p.getGoniometerMatrix(), mat); + + // Matrix must be 3x3 + Matrix<double> mat2(4,3); + TS_ASSERT_THROWS_ANYTHING( p.setGoniometerMatrix(mat2) ); + } + + void test_HKL() + { + Peak p(inst, 10000, 2.0); + p.setHKL(1.0, 2.0, 3.0); + TS_ASSERT_EQUALS( p.getH(), 1.0); + TS_ASSERT_EQUALS( p.getK(), 2.0); + TS_ASSERT_EQUALS( p.getL(), 3.0); + p.setH(5); + p.setK(6); + p.setL(7); + TS_ASSERT_EQUALS( p.getH(), 5.0); + TS_ASSERT_EQUALS( p.getK(), 6.0); + TS_ASSERT_EQUALS( p.getL(), 7.0); + } + void test_getBank_and_row() { Peak p(inst, 10000, 2.0); TS_ASSERT_EQUALS(p.getBankName(), "bank1") TS_ASSERT_EQUALS(p.getRow(), 0) TS_ASSERT_EQUALS(p.getCol(), 0) + p.setDetectorID(10050); + TS_ASSERT_EQUALS(p.getRow(), 50) + TS_ASSERT_EQUALS(p.getCol(), 0) + p.setDetectorID(10100); + TS_ASSERT_EQUALS(p.getRow(), 0) + TS_ASSERT_EQUALS(p.getCol(), 1) } - };