Commit 6ccdb3ea authored by Gigg, Martyn Anthony's avatar Gigg, Martyn Anthony
Browse files

Added separate Compton profile functions. Refs #6316

These will replace NCSCountRate as they are more fleixble and they
isolate the code more effectively. It also allows a new profile to be
used without too much effort.
parent 5d4d2599
......@@ -58,12 +58,15 @@ public:
///Destructor
virtual ~CompositeFunction();
/* Overriden methods */
/* Overriden methods */
/// Returns the function's name
virtual std::string name()const {return "CompositeFunction";}
/// Writes itself into a string
std::string asString()const;
/// Sets the workspace for each member function
void setWorkspace(boost::shared_ptr<const Workspace> ws);
/// Function you want to fit to.
/// @param domain :: The buffer for writing the calculated values. Must be big enough to accept dataSize() values
virtual void function(const FunctionDomain& domain, FunctionValues& values)const;
......
......@@ -54,12 +54,13 @@ public:
/// Set new peak radius
static void setPeakRadius(const int& r = 5);
protected:
/// Function evaluation method to be implemented in the inherited classes
virtual void functionLocal(double* out, const double* xValues, const size_t nData)const = 0;
/// Derivative evaluation method to be implemented in the inherited classes
virtual void functionDerivLocal(Jacobian* out, const double* xValues, const size_t nData) = 0;
/// Defines the area around the centre where the peak values are to be calculated (in FWHM).
protected:
/// Defines the area around the centre where the peak values are to be calculated (in FWHM).
static int s_peakRadius;
};
......
......@@ -119,6 +119,19 @@ std::string CompositeFunction::asString()const
return ostr.str();
}
/**
* @param ws A pointer to the workspace being fitted
*/
void CompositeFunction::setWorkspace(boost::shared_ptr<const Workspace> ws)
{
// Pass it on to each member
auto iend = m_functions.end();
for(auto it = m_functions.begin(); it != iend; ++it)
{
(*it)->setWorkspace(ws);
}
}
/** Function you want to fit to.
* @param domain :: The buffer for writing the calculated values. Must be big enough to accept dataSize() values
*/
......
......@@ -10,6 +10,7 @@ set ( SRC_FILES
src/Bk2BkExpConvPV.cpp
src/BoundaryConstraint.cpp
src/Chebyshev.cpp
src/ComptonProfile.cpp
src/Convolution.cpp
src/CostFuncFitting.cpp
src/CostFuncLeastSquares.cpp
......@@ -33,7 +34,9 @@ set ( SRC_FILES
src/GausOsc.cpp
src/Gaussian.cpp
src/Gaussian1D.cpp
src/GaussianComptonProfile.cpp
src/GaussianLinearBG1D.cpp
src/GramCharlierComptonProfile.cpp
src/IkedaCarpenterPV.cpp
src/LeBailFit.cpp
src/LeBailFit2.cpp
......@@ -93,6 +96,7 @@ set ( INC_FILES
inc/MantidCurveFitting/BoundaryConstraint.h
inc/MantidCurveFitting/Chebyshev.h
inc/MantidCurveFitting/CompositeValues.h
inc/MantidCurveFitting/ComptonProfile.h
inc/MantidCurveFitting/Convolution.h
inc/MantidCurveFitting/CostFuncFitting.h
inc/MantidCurveFitting/CostFuncLeastSquares.h
......@@ -120,7 +124,9 @@ set ( INC_FILES
inc/MantidCurveFitting/GausOsc.h
inc/MantidCurveFitting/Gaussian.h
inc/MantidCurveFitting/Gaussian1D.h
inc/MantidCurveFitting/GaussianComptonProfile.h
inc/MantidCurveFitting/GaussianLinearBG1D.h
inc/MantidCurveFitting/GramCharlierComptonProfile.h
inc/MantidCurveFitting/IkedaCarpenterPV.h
inc/MantidCurveFitting/Jacobian.h
inc/MantidCurveFitting/LeBailFit.h
......@@ -177,6 +183,7 @@ set ( TEST_FILES
BoundaryConstraintTest.h
ChebyshevTest.h
CompositeFunctionTest.h
ComptonProfileTest.h
ConvolutionTest.h
DampingMinimizerTest.h
DeltaFunctionTest.h
......@@ -193,10 +200,12 @@ set ( TEST_FILES
GausDecayTest.h
GausOscTest.h
Gaussian1DTest.h
GaussianComptonProfileTest.h
GaussianTest.h
GramCharlierComptonProfileTest.h
IkedaCarpenterPVTest.h
LeBailFitTest.h
LeBailFit2Test.h
LeBailFitTest.h
LeBailFunctionTest.h
LeastSquaresTest.h
LevenbergMarquardtMDTest.h
......
#ifndef MANTID_CURVEFITTING_COMPTONPROFILE_H_
#define MANTID_CURVEFITTING_COMPTONPROFILE_H_
#include "MantidCurveFitting/DllConfig.h"
#include "MantidAPI/IPeakFunction.h"
#include "MantidAPI/MatrixWorkspace.h"
#include "MantidAPI/ParamFunction.h"
namespace Mantid
{
namespace CurveFitting
{
/**
This class serves as a base-class for ComptonProfile type functions. @see GaussianComptonProfile, GramCharlierComptonProfile
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 ComptonProfile : public virtual API::ParamFunction, public virtual API::IFunction1D
{
public:
/// Default constructor required for factory
ComptonProfile();
/** @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 pul out the parameters
void setWorkspace(boost::shared_ptr<const API::Workspace> ws);
///@}
protected:
/// Declare parameters that will never participate in the fit
void declareAttributes();
/// Set an attribute value (and possibly cache its value)
void setAttribute(const std::string& name,const Attribute& value);
/// Override to calculate the value of the profile for this mass
virtual void massProfile(std::vector<double> & result,const double lorentzFWHM, const double resolutionFWHM) const = 0;
/// Access current y-values
inline const std::vector<double> & ySpace() const { return m_yspace; }
/// Access current Q values
inline const std::vector<double> & modQ() const { return m_modQ; }
/// Compute Voigt function interpolated around the given values
void voigtApproxDiff(std::vector<double> & voigtDiff, const std::vector<double> & yspace, const double lorentzPos, const double lorentzAmp,
const double lorentzWidth, const double gaussWidth) const;
/// Compute Voigt function
void voigtApprox(std::vector<double> & voigt, const std::vector<double> & yspace, const double lorentzPos, const double lorentzAmp,
const double lorentzWidth, const double gaussWidth) const;
private:
/// Retrieve a component parameter
double getComponentParameter(const Geometry::IComponent & comp,const std::string &name) const;
/// The workspace providing the data
API::MatrixWorkspace_const_sptr m_workspace;
/// Current workspace index, required to access instrument parameters
size_t m_wsIndex;
/// Store the mass values
double m_mass;
/// Source to sample distance
double m_l1;
/// Std dev of l1 distance
double m_sigmaL1;
/// Sample to detector distance
double m_l2;
/// Std dev of l1 distance
double m_sigmaL2;
/// Theta value for this spectrum
double m_theta;
/// Std Dev theta value for this spectrum
double m_sigmaTheta;
/// Final energy
double m_e1;
/// T0 value for this spectrum in seconds
double m_t0;
/// Lorentzian HWHM of the foil analyser energy
double m_hwhmGaussE;
/// Gaussian HWHM of the foil analyser energy
double m_hwhmLorentzE;
/// Voigt function
boost::shared_ptr<API::IPeakFunction> m_voigt;
/** @name Caches for commonly used values*/
///@{
/// Current y-values
mutable std::vector<double> m_yspace;
/// Current Q values
mutable std::vector<double> m_modQ;
///@}
};
} // namespace CurveFitting
} // namespace Mantid
#endif /* MANTID_CURVEFITTING_COMPTONPROFILE_H_ */
#ifndef MANTID_CURVEFITTING_GAUSSIANCOMPTONPROFILE_H_
#define MANTID_CURVEFITTING_GAUSSIANCOMPTONPROFILE_H_
#include "MantidCurveFitting/DllConfig.h"
#include "MantidCurveFitting/ComptonProfile.h"
namespace Mantid
{
namespace CurveFitting
{
/**
Implements a function to calculate the Compton profile of a nucleus using a Gaussian approximation
convoluted with an instrument resolution function that is approximated by a Voigt function. The function calculates
\f[\frac{E_0I(E_0)}{q}A_M J_M(y_M)\otimes R_M(t)\f]
for the given mass where \f$J_M\f$ is approximated using a Gaussian and \f$R_M\f$ is the resolution function
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 GaussianComptonProfile : public ComptonProfile
{
public:
/// Default constructor required for factory
GaussianComptonProfile();
private:
/// A string identifier for this function
std::string name() const;
/// Declare the function parameters
void declareParameters();
/// Compute the function
void massProfile(std::vector<double> & result,const double lorentzFWHM, const double resolutionFWHM) const;
};
} // namespace CurveFitting
} // namespace Mantid
#endif /* MANTID_CURVEFITTING_GAUSSIANCOMPTONPROFILE_H_ */
#ifndef MANTID_CURVEFITTING_GRAMCHARLIERCOMPTONPROFILE_H_
#define MANTID_CURVEFITTING_GRAMCHARLIERCOMPTONPROFILE_H_
#include "MantidCurveFitting/DllConfig.h"
#include "MantidCurveFitting/ComptonProfile.h"
namespace Mantid
{
namespace CurveFitting
{
/**
Implements a function to calculate the Compton profile of a nucleus using a Gram-Charlier approximation
convoluted with an instrument resolution function that is approximated by a Voigt function.
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 DLLExport GramCharlierComptonProfile : public ComptonProfile
{
public:
/// Default constructor required by factory
GramCharlierComptonProfile();
private:
/// A string identifier for this function
std::string name() const;
/// Declare the function parameters
void declareParameters();
/// Declare non-parameter attributes
void declareAttributes();
/// Set an attribute value (and possibly cache its value)
void setAttribute(const std::string& name,const Attribute& value);
/// Parse the active hermite polynomial coefficents
void setHermiteCoefficients(const std::string & coeffs);
/// Declare the Gram-Charlier (Hermite) coefficients
void declareGramCharlierParameters();
/// Compute the function
void massProfile(std::vector<double> & result,const double lorentzFWHM, const double resolutionFWHM) const;
/// The active hermite coefficents
std::vector<short> m_hermite;
};
} // namespace CurveFitting
} // namespace Mantid
#endif /* MANTID_CURVEFITTING_GRAMCHARLIERCOMPTONPROFILE_H_ */
//-----------------------------------------------------------------------------
// Includes
//-----------------------------------------------------------------------------
#include "MantidCurveFitting/ComptonProfile.h"
#include "MantidAPI/FunctionFactory.h"
#include <gsl/gsl_poly.h>
namespace Mantid
{
namespace CurveFitting
{
namespace
{
///@cond
const char * WSINDEX_NAME = "WorkspaceIndex";
const char * MASS_NAME = "Mass";
const double STDDEV_TO_HWHM = std::sqrt(std::log(4));
///@endcond
}
/**
*/
ComptonProfile::ComptonProfile() : API::ParamFunction(), API::IFunction1D(),
m_workspace(), m_wsIndex(0), m_mass(0.0), m_l1(0.0), m_sigmaL1(0.0),
m_l2(0.0), m_sigmaL2(0.0), m_theta(0.0), m_sigmaTheta(0.0), m_e1(0.0),
m_t0(0.0), m_hwhmGaussE(0.0), m_hwhmLorentzE(0.0), m_voigt()
{}
/**
*/
void ComptonProfile::declareAttributes()
{
declareAttribute(WSINDEX_NAME, IFunction::Attribute(static_cast<int>(m_wsIndex)));
declareAttribute(MASS_NAME, IFunction::Attribute(m_mass));
}
/**
* @param name The name of the attribute
* @param value The attribute's value
*/
void ComptonProfile::setAttribute(const std::string& name,const Attribute& value)
{
if(name == WSINDEX_NAME) m_wsIndex = static_cast<size_t>(value.asInt());
else if(name == MASS_NAME) m_mass = value.asDouble();
}
/**
* Also caches parameters from the instrument
* Throws if it is not a MatrixWorkspace
* @param ws The workspace set as input
*/
void ComptonProfile::setWorkspace(boost::shared_ptr<const API::Workspace> ws)
{
m_workspace = boost::dynamic_pointer_cast<const API::MatrixWorkspace>(ws);
if(!m_workspace)
{
throw std::invalid_argument("NCSCountRate expected an object of type MatrixWorkspace, type=" + ws->id());
}
auto inst = m_workspace->getInstrument();
auto sample = inst->getSample();
auto source = inst->getSource();
if(!sample || !source)
{
throw std::invalid_argument("NCSCountRate - Workspace has no source/sample.");
}
Geometry::IDetector_const_sptr det;
try
{
det = m_workspace->getDetector(m_wsIndex);
}
catch (Kernel::Exception::NotFoundError &)
{
throw std::invalid_argument("NCSCountRate - Workspace has not detector attached to histogram at index " + boost::lexical_cast<std::string>(m_wsIndex));
}
m_l1 = sample->getDistance(*source);
m_l2 = det->getDistance(*sample);
m_theta = m_workspace->detectorTwoTheta(det);
// parameters
m_sigmaL1 = getComponentParameter(*det, "sigma_l1");
m_sigmaL2 = getComponentParameter(*det, "sigma_l2");
m_sigmaTheta = getComponentParameter(*det, "sigma_theta");
m_e1 = getComponentParameter(*det,"efixed");
m_t0 = getComponentParameter(*det,"t0")*1e-6; // Convert to seconds
m_hwhmLorentzE = getComponentParameter(*det, "hwhm_energy_lorentz");
m_hwhmGaussE = STDDEV_TO_HWHM*getComponentParameter(*det, "sigma_energy_gauss");
}
/*
* Creates the internal caches
*/
void ComptonProfile::setUpForFit()
{
using namespace Mantid::API;
m_voigt = boost::dynamic_pointer_cast<IPeakFunction>(FunctionFactory::Instance().createFunction("Voigt"));
}
//-------------------------------------- Function evaluation -----------------------------------------
/**
* 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. It is assumed to be times in microseconds
* @param nData The length of the out & xValues arrays
*/
void ComptonProfile::function1D(double* out, const double* xValues, const size_t nData) const
{
std::vector<double> tInSecs(xValues, xValues + nData);
std::transform(tInSecs.begin(), tInSecs.end(), tInSecs.begin(), std::bind2nd(std::multiplies<double>(), 1e-6)); // Convert to seconds
const double mn = PhysicalConstants::NeutronMassAMU;
const double mevToK = PhysicalConstants::E_mev_toNeutronWavenumberSq;
const double massToMeV = 0.5*PhysicalConstants::NeutronMass/PhysicalConstants::meV; // Includes factor of 1/2
const double v1 = std::sqrt(m_e1/massToMeV);
const double k1 = std::sqrt(m_e1/mevToK);
const double l2l1 = m_l2/m_l1;
// -------------------------------------- Resolution dependence -------------------------------------
double x0(0.0),x1(0.0);
gsl_poly_solve_quadratic(m_mass-1.0, 2.0*std::cos(m_theta), -(m_mass+1.0), &x0, &x1);
const double k0k1 = std::max(x0,x1); // K0/K1 at y=0
double qy0(0.0), wgauss(0.0), lorentzFWHM(0.0);
;
if(m_mass > 1.0)
{
qy0 = std::sqrt(k1*k1*m_mass*(k0k1*k0k1 - 1));
double k0k1p3 = std::pow(k0k1,3);
double r1 = -(1.0 + l2l1*k0k1p3);
double r2 = 1.0 - l2l1*k0k1p3 + l2l1*std::pow(k0k1,2)*std::cos(m_theta) - k0k1*std::cos(m_theta);
double factor = (0.2413/qy0)*((m_mass/mn)*r1 - r2);
lorentzFWHM = std::abs(factor*m_hwhmLorentzE*2);
wgauss = std::abs(factor*m_hwhmGaussE*2);
}
else
{
qy0 = k1*std::tan(m_theta);
double factor = (0.2413*2.0/k1)*std::abs((std::cos(m_theta) + l2l1)/std::sin(m_theta));
lorentzFWHM = m_hwhmLorentzE*factor;
wgauss = m_hwhmGaussE*factor;
}
double k0y0 = k1*k0k1; // k0_y0 = k0 value at y=0
double wtheta = 2.0*STDDEV_TO_HWHM*std::abs(k0y0*k1*std::sin(m_theta)/qy0)*m_sigmaTheta;
double common = (m_mass/mn) - 1 + k1*std::cos(m_theta)/k0y0;
double wl1 = 2.0*STDDEV_TO_HWHM*std::abs((std::pow(k0y0,2)/(qy0*m_l1))*common)*m_sigmaL1;
double wl2 = 2.0*STDDEV_TO_HWHM*std::abs((std::pow(k0y0,3)/(k1*qy0*m_l1))*common)*m_sigmaL2;
const double resolutionSigma = std::sqrt(std::pow(wgauss,2) + std::pow(wtheta,2) + std::pow(wl1,2) + std::pow(wl2,2));
// -------------------------------------- Profile factor --------------------------------------------------------------
// Calculate energy dependent factors and transform q to Yspace
std::vector<double> e0(nData,0.0);
m_modQ.resize(nData);
m_yspace.resize(nData);
for(size_t i = 0; i < nData; ++i)
{
const double v0 = m_l1/(tInSecs[i] - m_t0 - (m_l2/v1));
const double ei = massToMeV*v0*v0;
e0[i] = ei;
const double w = ei - m_e1;
const double k0 = std::sqrt(ei/mevToK);
const double q = std::sqrt(k0*k0 + k1*k1 - 2.0*k0*k1*std::cos(m_theta));
m_modQ[i] = q;
m_yspace[i] = 0.2393*(m_mass/q)*(w - mevToK*q*q/m_mass);
}
std::vector<double> profile(nData);
this->massProfile(profile,lorentzFWHM, resolutionSigma);
// Multiply by the mass & final prefactor e0^0.1/q
for(size_t j = 0; j < nData; ++j)
{
out[j] = std::pow(e0[j],0.1)*m_mass*profile[j]/m_modQ[j];
}
// for(size_t i = 0; i < nmasses; ++i)
// {
// const auto & yi = yspace[i];
// auto & j1i = j1[i];
// std::ostringstream os;
// os << WIDTH_PREFIX << i;
// const double gaussWidth(getParameter(os.str()));
// const double lorentzWidth(lorentzW[i]);
// const double gaussRes = sigmaRes[i];
// if(i == 0)
// {
// const double amp(1.0);
// firstMassJ(j1i, yi, modQ, amp, kfse, gaussWidth, lorentzWidth, gaussRes);
// }
// else
// {
// os.str("");
// os << INTENSITY_PREFIX << i;
// double lorentzPos(0.0), amplitude(getParameter(os.str())), lorentzFWHM(lorentzWidth);
// double gaussFWHM = std::sqrt(std::pow(gaussRes,2) + std::pow(2.0*STDDEV_TO_HWHM*gaussWidth,2));
// voigtApprox(j1i, yi,lorentzPos, amplitude, lorentzFWHM, gaussFWHM); // Answer goes into j1i
// voigtApproxDiff(voigtDiffResult, yi, lorentzPos, amplitude, lorentzFWHM, gaussFWHM);
// for(size_t j = 0; j < nData; ++j)
// {
// const double factor = std::pow(gaussWidth,4.0)/(3.0*modQ[j]);
// j1i[j] -= factor*voigtDiffResult[j];
// }
// }
// // Multiply by mass
// std::transform(j1i.begin(), j1i.end(), j1i.begin(), std::bind2nd(std::multiplies<double>(), m_masses[i]));
// }
//
// // Sum over each mass and scale by prefactor to get answer
// for(size_t j = 0; j < nData; ++j)
// {
// for(size_t i = 0; i < nmasses; ++i)
// {
// out[j] += j1[i][j];
// }
// out[j] *= std::pow(e0[j],0.1)/modQ[j];