From 6c9b13e71e5c8b07939990cdf4faec8c1a0ed67f Mon Sep 17 00:00:00 2001 From: Martyn Gigg Date: Fri, 13 Dec 2013 14:27:22 +0000 Subject: [PATCH] Implement a profile function to fit a mass peak in Y-space. Used in neutron compton scattering to Fit a spectrum using only the proton peak. Refs #8598 --- .../Framework/CurveFitting/CMakeLists.txt | 3 + .../MantidCurveFitting/ComptonPeakProfile.h | 88 ++++++++++ .../inc/MantidCurveFitting/ConvertToYSpace.h | 2 + .../CurveFitting/src/ComptonPeakProfile.cpp | 157 ++++++++++++++++++ .../CurveFitting/src/ComptonProfile.cpp | 7 +- .../CurveFitting/src/ConvertToYSpace.cpp | 35 ++++ .../test/ComptonPeakProfileTest.h | 85 ++++++++++ 7 files changed, 371 insertions(+), 6 deletions(-) create mode 100644 Code/Mantid/Framework/CurveFitting/inc/MantidCurveFitting/ComptonPeakProfile.h create mode 100644 Code/Mantid/Framework/CurveFitting/src/ComptonPeakProfile.cpp create mode 100644 Code/Mantid/Framework/CurveFitting/test/ComptonPeakProfileTest.h diff --git a/Code/Mantid/Framework/CurveFitting/CMakeLists.txt b/Code/Mantid/Framework/CurveFitting/CMakeLists.txt index da7f7750ac9c..8a9e2a07462e 100644 --- a/Code/Mantid/Framework/CurveFitting/CMakeLists.txt +++ b/Code/Mantid/Framework/CurveFitting/CMakeLists.txt @@ -13,6 +13,7 @@ set ( SRC_FILES src/BoundaryConstraint.cpp src/CalculateGammaBackground.cpp src/Chebyshev.cpp + src/ComptonPeakProfile.cpp src/ComptonProfile.cpp src/ComptonScatteringCountRate.cpp src/ConvertToYSpace.cpp @@ -112,6 +113,7 @@ set ( INC_FILES inc/MantidCurveFitting/CalculateGammaBackground.h inc/MantidCurveFitting/Chebyshev.h inc/MantidCurveFitting/CompositeValues.h + inc/MantidCurveFitting/ComptonPeakProfile.h inc/MantidCurveFitting/ComptonProfile.h inc/MantidCurveFitting/ComptonScatteringCountRate.h inc/MantidCurveFitting/ConvertToYSpace.h @@ -214,6 +216,7 @@ set ( TEST_FILES CalculateGammaBackgroundTest.h ChebyshevTest.h CompositeFunctionTest.h + ComptonPeakProfileTest.h ComptonProfileTest.h ComptonScatteringCountRateTest.h ConvertToYSpaceTest.h diff --git a/Code/Mantid/Framework/CurveFitting/inc/MantidCurveFitting/ComptonPeakProfile.h b/Code/Mantid/Framework/CurveFitting/inc/MantidCurveFitting/ComptonPeakProfile.h new file mode 100644 index 000000000000..abce8939024a --- /dev/null +++ b/Code/Mantid/Framework/CurveFitting/inc/MantidCurveFitting/ComptonPeakProfile.h @@ -0,0 +1,88 @@ +#ifndef MANTID_CURVEFITTING_COMPTONPEAKPROFILE_H_ +#define MANTID_CURVEFITTING_COMPTONPEAKPROFILE_H_ + +#include "MantidCurveFitting/DllConfig.h" +#include "MantidAPI/IPeakFunction.h" +#include "MantidAPI/MatrixWorkspace.h" +#include "MantidAPI/ParamFunction.h" + +namespace Mantid +{ +namespace CurveFitting +{ + /** + This implements a resolution function for fitting a single mass in a compton scattering spectrum. The + function has two domains defined by the VoigtEnergyCutOff attribute: + + - < VoigtEnergyCutOff: Voigt approximation to spectrum is used + - >= VoigtEnergyCutOff: Gaussian approximation to spectrum is used. + + Both are normalized by the peak area. + + Copyright © 2013 ISIS Rutherford Appleton Laboratory & NScD Oak Ridge National Laboratory + + This file is part of Mantid. + + Mantid is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 3 of the License, or + (at your option) any later version. + + Mantid is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with this program. If not, see . + + File change history is stored at: + Code Documentation is available at: + */ + class MANTID_CURVEFITTING_DLL ComptonPeakProfile : public virtual API::ParamFunction, + public virtual API::IFunction1D + { + public: + /// Default constructor required for factory + ComptonPeakProfile(); + + private: + std::string name() const; + + /** @name Function evaluation */ + ///@{ + /// Calculate the function + void function1D(double* out, const double* xValues, const size_t nData) const; + /// Ensure the object is ready to be fitted + void setUpForFit(); + /// Cache a copy of the workspace pointer and pull out the parameters + void setWorkspace(boost::shared_ptr ws); + ///@} + + void declareParameters(); + void declareAttributes(); + void setAttribute(const std::string& name,const Attribute& value); + + /// WorkspaceIndex attribute + size_t m_wsIndex; + /// Mass of peak + double m_mass; + /// Below this value a Voigt is used for profile approximation + double m_voigtCutOff; + + /// Gaussian function for lower-energy peaks + boost::shared_ptr m_gauss; + /// Voigt function for higher-energy peaks + boost::shared_ptr m_voigt; + + /// Final energy of analyser + double m_efixed; + /// Calculated value of lorentz width + double m_hwhmLorentz; + }; + + +} // namespace CurveFitting +} // namespace Mantid + +#endif /* MANTID_CURVEFITTING_COMPTONPEAKPROFILE_H_ */ diff --git a/Code/Mantid/Framework/CurveFitting/inc/MantidCurveFitting/ConvertToYSpace.h b/Code/Mantid/Framework/CurveFitting/inc/MantidCurveFitting/ConvertToYSpace.h index 6d6c7f1f9049..b824b1c9164c 100644 --- a/Code/Mantid/Framework/CurveFitting/inc/MantidCurveFitting/ConvertToYSpace.h +++ b/Code/Mantid/Framework/CurveFitting/inc/MantidCurveFitting/ConvertToYSpace.h @@ -53,6 +53,8 @@ namespace CurveFitting int version() const; const std::string category() const; + /// Creates a POD struct containing the required detector parameters for this spectrum + static DetectorParams getDetectorParameters(const API::MatrixWorkspace_const_sptr & ws, const size_t index); /// Retrieve a component parameter static double getComponentParameter(const Geometry::IComponent_const_sptr & comp, const Geometry::ParameterMap &pmap, diff --git a/Code/Mantid/Framework/CurveFitting/src/ComptonPeakProfile.cpp b/Code/Mantid/Framework/CurveFitting/src/ComptonPeakProfile.cpp new file mode 100644 index 000000000000..60b9cd4ea1aa --- /dev/null +++ b/Code/Mantid/Framework/CurveFitting/src/ComptonPeakProfile.cpp @@ -0,0 +1,157 @@ +//----------------------------------------------------------------------------- +// Includes +//----------------------------------------------------------------------------- +#include "MantidCurveFitting/ComptonPeakProfile.h" +#include "MantidCurveFitting/ConvertToYSpace.h" +#include "MantidAPI/FunctionFactory.h" + +namespace Mantid +{ +namespace CurveFitting +{ + + DECLARE_FUNCTION(ComptonPeakProfile); + + namespace + { + ///@cond + const char * WSINDEX_NAME = "WorkspaceIndex"; + const char * MASS_NAME = "Mass"; + const char * VOIGT_CUT_OFF = "VoigtEnergyCutOff"; + + const char * AMP_PARAM = "Intensity"; + const char * POS_PARAM = "Position"; + const char * WIDTH_PARAM = "SigmaGauss"; + + /// Conversion constant + const double MASS_TO_MEV = 0.5*PhysicalConstants::NeutronMass/PhysicalConstants::meV; + const double STDDEV_TO_FWHM = 0.5*std::sqrt(std::log(4.0)); + ///@endcond + } + + /** + */ + ComptonPeakProfile::ComptonPeakProfile() : API::ParamFunction(), API::IFunction1D(), + m_wsIndex(0), m_mass(0.0), m_voigtCutOff(5000.), + m_gauss(), m_voigt(), m_efixed(0.0), m_hwhmLorentz(0.0) + {} + + //-------------------------------------- Private methods ----------------------------------------- + + /// A string identifier for this function + std::string ComptonPeakProfile::name() const { return "ComptonPeakProfile"; } + + /** + * Calculates the value of the function for each x value and stores in the given output array + * @param out An array of size nData to store the results + * @param xValues The input X data array of size nData. + * @param nData The length of the out & xValues arrays + */ + void ComptonPeakProfile::function1D(double* out, const double* xValues, const size_t nData) const + { + double amp(getParameter(0)), pos(getParameter(1)), gaussSigma(getParameter(2)); + + if(m_efixed < m_voigtCutOff) + { + double gaussFWHM(getParameter(2)*STDDEV_TO_FWHM), lorentzFWHM(2.0*m_hwhmLorentz); + m_voigt->setParameter(0,amp); + m_voigt->setParameter(1,pos); + m_voigt->setParameter(2,lorentzFWHM); + m_voigt->setParameter(3, gaussFWHM); + m_voigt->functionLocal(out, xValues, nData); + const double norm = 1.0/(0.5*M_PI*lorentzFWHM); + std::transform(out, out + nData, out, std::bind2nd(std::multiplies(), norm)); + } + else + { + double sigmaTotalSq = m_hwhmLorentz*m_hwhmLorentz + gaussSigma*gaussSigma; + // Our gaussian isn't normalised by the width. Here we set the height to + // amp/(2*sigma^2) so that it will be normalised correctly + m_gauss->setParameter(0, 0.5*amp/M_PI/sigmaTotalSq); + m_gauss->setParameter(1, pos); + m_gauss->setParameter(2, sqrt(sigmaTotalSq)); + m_gauss->functionLocal(out, xValues, nData); + } + } + + /* + * Creates the internal caches + */ + void ComptonPeakProfile::setUpForFit() + { + // Voigt & Gaussian + using namespace Mantid::API; + m_gauss = boost::dynamic_pointer_cast(FunctionFactory::Instance().createFunction("Gaussian")); + m_voigt = boost::dynamic_pointer_cast(FunctionFactory::Instance().createFunction("Voigt")); + } + + /** + * Also caches parameters from the instrument + * Throws if it is not a MatrixWorkspace + * @param ws The workspace set as input + */ + void ComptonPeakProfile::setWorkspace(boost::shared_ptr ws) + { + auto workspace = boost::dynamic_pointer_cast(ws); + if(!workspace) + { + throw std::invalid_argument("ComptonPeakProfile expected an object of type MatrixWorkspace, type=" + ws->id()); + } + + DetectorParams detpar = ConvertToYSpace::getDetectorParameters(workspace, m_wsIndex); + m_efixed = detpar.efixed; + + // Recoil time + const double sth = sin(detpar.theta); + const double distFact = (cos(detpar.theta) + sqrt(m_mass*m_mass - sth*sth))/(m_mass+1.0); + const double ei = detpar.efixed/pow(distFact,2.0); + const double v1 = std::sqrt(detpar.efixed/MASS_TO_MEV); + const double k1 = std::sqrt(detpar.efixed/PhysicalConstants::E_mev_toNeutronWavenumberSq); + const double v2 = std::sqrt(ei/MASS_TO_MEV); + const double trec = detpar.l1/v1 + detpar.l2/v2; + + // Compute lorentz width due to in Y due to spread in energy hwhm_lorentz + const double dELorentz = + ConvertToYSpace::getComponentParameter(workspace->getDetector(m_wsIndex), workspace->constInstrumentParameters(), + "hwhm_lorentz"); + double yplus(0.0), yminus(0.0), dummy(0.0); + detpar.efixed += dELorentz; + ConvertToYSpace::calculateY(yplus,dummy,dummy,m_mass,trec*1e6,k1,v1,detpar); + detpar.efixed -= 2.0*dELorentz; + ConvertToYSpace::calculateY(yminus,dummy,dummy,m_mass,trec*1e6,k1,v1,detpar); + // lorentzian width + m_hwhmLorentz = 0.5*(yplus-yminus); + } + + /** + */ + void ComptonPeakProfile::declareParameters() + { + declareParameter(AMP_PARAM, 1.0, "Intensity parameter"); + declareParameter(POS_PARAM, 1.0, "Peak position parameter"); + declareParameter(WIDTH_PARAM, 1.0, "Width parameter"); + } + + /** + */ + void ComptonPeakProfile::declareAttributes() + { + declareAttribute(WSINDEX_NAME, IFunction::Attribute(static_cast(m_wsIndex))); + declareAttribute(MASS_NAME, IFunction::Attribute(m_mass)); + declareAttribute(VOIGT_CUT_OFF, IFunction::Attribute(m_voigtCutOff)); + } + + /** + * @param name The name of the attribute + * @param value The attribute's value + */ + void ComptonPeakProfile::setAttribute(const std::string& name,const Attribute& value) + { + IFunction::setAttribute(name,value); // Make sure the base-class stores it + if(name == WSINDEX_NAME) m_wsIndex = static_cast(value.asInt()); + else if(name == MASS_NAME) m_mass = value.asDouble(); + else if(name == VOIGT_CUT_OFF) m_voigtCutOff = value.asDouble(); + } + +} // namespace CurveFitting +} // namespace Mantid diff --git a/Code/Mantid/Framework/CurveFitting/src/ComptonProfile.cpp b/Code/Mantid/Framework/CurveFitting/src/ComptonProfile.cpp index 6b63700ee92e..45df24eb2235 100644 --- a/Code/Mantid/Framework/CurveFitting/src/ComptonProfile.cpp +++ b/Code/Mantid/Framework/CurveFitting/src/ComptonProfile.cpp @@ -86,13 +86,8 @@ namespace CurveFitting throw std::invalid_argument("ComptonProfile - Workspace has no detector attached to histogram at index " + boost::lexical_cast(m_wsIndex)); } - DetectorParams detpar; + DetectorParams detpar = ConvertToYSpace::getDetectorParameters(workspace, m_wsIndex); const auto & pmap = workspace->constInstrumentParameters(); - detpar.l1 = sample->getDistance(*source); - detpar.l2 = det->getDistance(*sample); - detpar.theta = workspace->detectorTwoTheta(det); - detpar.t0 = ConvertToYSpace::getComponentParameter(det, pmap, "t0")*1e-6; // Convert to seconds - detpar.efixed = ConvertToYSpace::getComponentParameter(det, pmap, "efixed"); ResolutionParams respar; respar.dl1 = ConvertToYSpace::getComponentParameter(det, pmap, "sigma_l1"); diff --git a/Code/Mantid/Framework/CurveFitting/src/ConvertToYSpace.cpp b/Code/Mantid/Framework/CurveFitting/src/ConvertToYSpace.cpp index 1a01b2747ab2..217da6cc14f9 100644 --- a/Code/Mantid/Framework/CurveFitting/src/ConvertToYSpace.cpp +++ b/Code/Mantid/Framework/CurveFitting/src/ConvertToYSpace.cpp @@ -55,6 +55,41 @@ namespace CurveFitting } //---------------------------------------------------------------------------------------------- + /** + * @param ws The workspace with attached instrument + * @param index Index of the spectrum + * @return DetectorParams structure containing the relevant parameters + */ + DetectorParams ConvertToYSpace::getDetectorParameters(const API::MatrixWorkspace_const_sptr & ws, + const size_t index) + { + auto inst = ws->getInstrument(); + auto sample = inst->getSample(); + auto source = inst->getSource(); + if(!sample || !source) + { + throw std::invalid_argument("ConvertToYSpace - Workspace has no source/sample."); + } + Geometry::IDetector_const_sptr det; + try + { + det = ws->getDetector(index); + } + catch (Kernel::Exception::NotFoundError &) + { + throw std::invalid_argument("ConvertToYSpace - Workspace has no detector attached to histogram at index " + \ + boost::lexical_cast(index)); + } + + DetectorParams detpar; + const auto & pmap = ws->constInstrumentParameters(); + detpar.l1 = sample->getDistance(*source); + detpar.l2 = det->getDistance(*sample); + detpar.theta = ws->detectorTwoTheta(det); + detpar.t0 = ConvertToYSpace::getComponentParameter(det, pmap, "t0")*1e-6; // Convert to seconds + detpar.efixed = ConvertToYSpace::getComponentParameter(det, pmap, "efixed"); + return detpar; + } /** * If a DetectorGroup is encountered then the parameters are averaged over the group diff --git a/Code/Mantid/Framework/CurveFitting/test/ComptonPeakProfileTest.h b/Code/Mantid/Framework/CurveFitting/test/ComptonPeakProfileTest.h new file mode 100644 index 000000000000..1e52ccd2d9b8 --- /dev/null +++ b/Code/Mantid/Framework/CurveFitting/test/ComptonPeakProfileTest.h @@ -0,0 +1,85 @@ +#ifndef MANTID_CURVEFITTING_ComptonPeakProfileTEST_H_ +#define MANTID_CURVEFITTING_ComptonPeakProfileTEST_H_ + +#include + +#include "MantidCurveFitting/ComptonPeakProfile.h" +#include "ComptonProfileTestHelpers.h" + +using Mantid::CurveFitting::ComptonPeakProfile; + +class ComptonPeakProfileTest : public CxxTest::TestSuite +{ + +public: + // This pair of boilerplate methods prevent the suite being created statically + // This means the constructor isn't called when running other tests + static ComptonPeakProfileTest *createSuite() { return new ComptonPeakProfileTest(); } + static void destroySuite( ComptonPeakProfileTest *suite ) { delete suite; } + + void test_initialized_object_has_three_parameters() + { + auto profile = createFunction(); + TS_ASSERT_EQUALS(3, profile->nParams()); + } + + void test_initialized_object_has_expected_attributes() + { + auto profile = createFunction(); + static const size_t nattrs(3); + const char * expectedAttrs[nattrs] = {"WorkspaceIndex", "Mass", "VoigtEnergyCutOff"}; + + TS_ASSERT_EQUALS(nattrs, profile->nAttributes()); + + // Test names as they are used in scripts + if(profile->nAttributes() > 0) + { + std::set expectedAttrSet(expectedAttrs, expectedAttrs + nattrs); + std::vector actualNames = profile->getAttributeNames(); + + for(size_t i = 0; i < nattrs; ++i) + { + const std::string & name = actualNames[i]; + size_t keyCount = expectedAttrSet.count(name); + TSM_ASSERT_EQUALS("Expected " + name + " to be found as attribute but it was not.", 1, keyCount); + } + } + } + + void test_function_gives_expected_value_for_given_input_data() + { + using namespace Mantid::API; + + auto peakProfile = createFunction(); + + auto domain = boost::shared_ptr(new FunctionDomain1DVector(-1,1, 3)); + Mantid::API::FunctionValues outputs(*domain); + peakProfile->setParameter(0,0.93); + peakProfile->setParameter(1,0.4); + peakProfile->setParameter(2,4.29); + + peakProfile->function(*domain, outputs); + + TS_ASSERT_DELTA(0.14694800,outputs.getCalculated(0),1e-08); + TS_ASSERT_DELTA(0.34795949,outputs.getCalculated(1),1e-08); + TS_ASSERT_DELTA(0.31659214,outputs.getCalculated(2),1e-08); + } + +private: + + Mantid::API::IFunction_sptr createFunction() + { + Mantid::API::IFunction_sptr profile = boost::make_shared(); + profile->initialize(); + auto paramWS = ComptonProfileTestHelpers::createSingleSpectrumWorkspace(300,351,0.5,true,true); // Only using for parameters + profile->setAttributeValue("Mass", 1.0079); + TS_ASSERT_THROWS_NOTHING(profile->setWorkspace(paramWS)); + profile->setUpForFit(); + + return profile; + } + +}; + + +#endif /* MANTID_CURVEFITTING_COMPTONPEAKPROFILETEST_H_ */