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

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

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 &copy; 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 <http://www.gnu.org/licenses/>.
+
+    File change history is stored at: <https://github.com/mantidproject/mantid>
+    Code Documentation is available at: <http://doxygen.mantidproject.org>
+  */
+  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<const API::Workspace> 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<API::IPeakFunction> m_gauss;
+    /// Voigt function for higher-energy peaks
+    boost::shared_ptr<API::IPeakFunction> 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_ */

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,

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<double>(), 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<IPeakFunction>(FunctionFactory::Instance().createFunction("Gaussian"));
+    m_voigt = boost::dynamic_pointer_cast<IPeakFunction>(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<const API::Workspace> ws)
+  {
+    auto workspace = boost::dynamic_pointer_cast<const API::MatrixWorkspace>(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<int>(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<size_t>(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

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<std::string>(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");

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<std::string>(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