Skip to content
Snippets Groups Projects
Commit a8bb1c20 authored by Michael Wedel's avatar Michael Wedel
Browse files

Refs #11210. Adding documentation, more unit tests.

parent f40a3c76
No related branches found
No related tags found
No related merge requests found
......@@ -68,11 +68,6 @@ protected:
void functionDerivLocal(API::Jacobian *out, const double *xValues,
const size_t nData);
// void functionDeriv(const API::FunctionDomain &domain,
// API::Jacobian &jacobian) {
// calNumericalDeriv(domain, jacobian);
// }
void init();
};
......
......@@ -29,7 +29,7 @@ void PseudoVoigt::functionLocal(double *out, const double *xValues,
double xDiffSquared = (xValues[i] - x0) * (xValues[i] - x0);
out[i] = h * (gFraction * exp(-0.5 * xDiffSquared / sSquared) +
(lFraction * gSquared / (xDiffSquared + gSquared)));
(lFraction * gSquared / (xDiffSquared + gSquared)));
}
}
......@@ -59,25 +59,24 @@ void PseudoVoigt::functionDerivLocal(Jacobian *out, const double *xValues,
out->set(i, 0, h * (expTerm - lorentzTerm));
out->set(i, 1, gFraction * expTerm + lFraction * lorentzTerm);
out->set(i, 2, h * (gFraction * expTerm * xDiff / sSquared +
lFraction * lorentzTerm * xDiff * 2.0 /
(xDiffSquared + gSquared)));
out->set(i, 3,
gFraction * h * expTerm * xDiffSquared / sSquared / f +
lFraction * h * (lorentzTerm / g -
lorentzTerm * g / (xDiffSquared + gSquared)));
out->set(i, 2, h * xDiff * (gFraction * expTerm / sSquared +
lFraction * lorentzTerm * 2.0 /
(xDiffSquared + gSquared)));
out->set(i, 3, h * (gFraction * expTerm * xDiffSquared / sSquared / f +
lFraction * lorentzTerm *
(1.0 / g - g / (xDiffSquared + gSquared))));
}
}
void PseudoVoigt::init() {
declareParameter("Mixing");
declareParameter("Mixing", 1.0);
declareParameter("Height");
declareParameter("PeakCentre");
declareParameter("FWHM");
BoundaryConstraint *mixingConstraint =
new BoundaryConstraint(this, "Mixing", 0.0, 1.0, true);
mixingConstraint->setPenaltyFactor(1e7);
mixingConstraint->setPenaltyFactor(1e9);
addConstraint(mixingConstraint);
}
......
......@@ -8,8 +8,17 @@
#include "MantidCurveFitting/Jacobian.h"
#include <boost/make_shared.hpp>
#include "MantidCurveFitting/Gaussian.h"
#include "MantidCurveFitting/Lorentzian.h"
#include "MantidAPI/AlgorithmManager.h"
#include "MantidTestHelpers/WorkspaceCreationHelper.h"
#include "MantidDataObjects/Workspace2D.h"
#include "MantidKernel/MersenneTwister.h"
using namespace Mantid::CurveFitting;
using namespace Mantid::API;
using namespace Mantid::DataObjects;
class PseudoVoigtTest : public CxxTest::TestSuite {
public:
......@@ -182,6 +191,94 @@ public:
}
}
void testGaussianEdge() {
IFunction_sptr pv = getInitializedPV(1.0, 4.78, 0.05, 1.0);
Gaussian gaussian;
gaussian.initialize();
gaussian.setCentre(1.0);
gaussian.setHeight(4.78);
gaussian.setFwhm(0.05);
FunctionDomain1DVector domain(m_xValues);
FunctionValues valuesPV(domain);
FunctionValues valuesGaussian(domain);
pv->function(domain, valuesPV);
gaussian.function(domain, valuesGaussian);
for (size_t i = 0; i < valuesPV.size(); ++i) {
TS_ASSERT_DELTA(valuesPV[i], valuesGaussian[i], 1e-15);
}
}
void testLorentzianEdge() {
IFunction_sptr pv = getInitializedPV(1.0, 4.78, 0.05, 0.0);
Lorentzian lorentzian;
lorentzian.initialize();
lorentzian.setCentre(1.0);
lorentzian.setFwhm(0.05);
lorentzian.setHeight(4.78);
FunctionDomain1DVector domain(m_xValues);
FunctionValues valuesPV(domain);
FunctionValues valuesLorentzian(domain);
pv->function(domain, valuesPV);
lorentzian.function(domain, valuesLorentzian);
for (size_t i = 0; i < valuesPV.size(); ++i) {
TS_ASSERT_DELTA(valuesPV[i], valuesLorentzian[i], 1e-15);
}
}
void testFit() {
// Generating a workspace with function values
Workspace2D_sptr ws =
WorkspaceCreationHelper::Create1DWorkspaceConstant(100, 0.0, 0.0);
std::vector<double> &x = ws->dataX(0);
for (size_t i = 0; i < 100; ++i) {
x[i] = static_cast<double>(i) * 0.01 - 0.5;
}
FunctionDomain1DVector domain(x);
IFunction_sptr generatingPV = getInitializedPV(0.0, 112.78, 0.15, 0.7);
FunctionValues generatingValues(domain);
generatingPV->function(domain, generatingValues);
Mantid::Kernel::MersenneTwister rng(2, -0.5, 0.5);
std::vector<double> &y = ws->dataY(0);
std::vector<double> &e = ws->dataE(0);
for (size_t i = 0; i < 100; ++i) {
y[i] = rng.nextValue() + generatingValues[i];
e[i] = sqrt(fabs(y[i]));
}
// Some starting values for the fit
IFunction_sptr pv = getInitializedPV(0.03, 120.03, 0.1, 0.5);
IAlgorithm_sptr fit = AlgorithmManager::Instance().create("Fit");
fit->setProperty("Function", pv);
fit->setProperty("InputWorkspace", ws);
fit->execute();
TS_ASSERT(fit->isExecuted());
IFunction_sptr fitted = fit->getProperty("Function");
TS_ASSERT_DELTA(fitted->getError(0), 0.0, 1e-6);
TS_ASSERT_DELTA(fitted->getError(2), 0.0, 1e-6);
TS_ASSERT_DELTA(fitted->getError(1), 0.0, 1e-6);
TS_ASSERT_DELTA(fitted->getError(3), 0.0, 1e-6);
TS_ASSERT_DELTA(fitted->getParameter("Mixing"), 0.7, 1e-2);
TS_ASSERT_DELTA(fitted->getParameter("PeakCentre"), 0.0, 1e-4);
TS_ASSERT_DELTA(fitted->getParameter("Height"), 112.78, 0.5);
TS_ASSERT_DELTA(fitted->getParameter("FWHM"), 0.15, 1e-2);
}
private:
IFunction_sptr getInitializedPV(double center, double height, double fwhm,
double mixing) {
......
.. _func-PseudoVoigt:
===========
PseudoVoigt
===========
.. index:: PseudoVoigt
Description
-----------
The Pseudo-Voigt function is an approximation for the Voigt function, which is a convolution of Gaussian and Lorentzian function. It is often used as a peak profile in powder diffraction for cases where neither a pure Gaussian or Lorentzian function appropriately describe a peak.
Instead of convoluting those two functions, the Pseudo-Voigt function is defined as the sum of a Gaussian peak :math:`G(x)` and a Lorentzian peak :math:`L(x)`, weighted by a fourth parameter :math:`\eta` (values between 0 and 1) which shifts the profile more towards pure Gaussian or pure Lorentzian when approaching 1 or 0 respectively:
.. math:: PV(x) = \eta G(x) + (1 - \eta)L(x)
Both functions share three parameters: Height (height of the peak at the maximum), PeakCentre (position of the maximum) and FWHM (full width at half maximum of the peak).
The figure below shows data together with a fitted Pseudo-Voigt function, as well as Gaussian and Lorentzian with equal parameters. The mixing parameter for that example is 0.7, which means that the function is behaving more like a Gaussian.
.. figure:: /images/PseudoVoigt.png
:alt: Comparison of Pseudo-Voigt function with Gaussian and Lorentzian profiles.
.. attributes::
.. properties::
.. categories::
Code/Mantid/docs/source/images/PseudoVoigt.png

39 KiB

0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment