Mantid
Loading...
Searching...
No Matches
Lorentzian.cpp
Go to the documentation of this file.
1// Mantid Repository : https://github.com/mantidproject/mantid
2//
3// Copyright © 2018 ISIS Rutherford Appleton Laboratory UKRI,
4// NScD Oak Ridge National Laboratory, European Spallation Source,
5// Institut Laue - Langevin & CSNS, Institute of High Energy Physics, CAS
6// SPDX - License - Identifier: GPL - 3.0 +
7//----------------------------------------------------------------------
8// Includes
9//----------------------------------------------------------------------
12#include <cmath>
13
15
16using namespace CurveFitting;
17
18using namespace Kernel;
19
20using namespace API;
21
22DECLARE_FUNCTION(Lorentzian)
23
24Lorentzian::Lorentzian() : m_amplitudeEqualHeight(false) {}
25
27 declareParameter("Amplitude", 1.0, "Intensity scaling");
28 declareParameter("PeakCentre", 0.0, "Centre of peak");
29 declareParameter("FWHM", 0.0, "Full-width at half-maximum");
30}
31
32double Lorentzian::height() const {
33 double gamma = getParameter("FWHM");
34 if (gamma == 0.0) {
35 return getParameter("Amplitude");
36 } else {
37 return getParameter("Amplitude") * 2.0 / (gamma * M_PI);
38 }
39}
40
41void Lorentzian::setHeight(const double h) {
42 double gamma = getParameter("FWHM");
43 if (gamma == 0.0) {
45 setParameter("Amplitude", h);
46 } else {
47 setParameter("Amplitude", h * gamma * M_PI / 2.0);
48 }
49}
50
51void Lorentzian::setFwhm(const double w) {
52 auto gamma = getParameter("FWHM");
53 if (gamma == 0.0 && w != 0.0 && m_amplitudeEqualHeight) {
54 auto h = getParameter("Amplitude");
55 setParameter("Amplitude", h * w * M_PI / 2.0);
56 }
57 if (w != 0.0) {
59 }
60 setParameter("FWHM", w);
61}
62
63void Lorentzian::fixCentre(bool isDefault) { fixParameter("PeakCentre", isDefault); }
64
65void Lorentzian::unfixCentre() { unfixParameter("PeakCentre"); }
66
67void Lorentzian::fixIntensity(bool isDefault) { fixParameter("Amplitude", isDefault); }
68
70
71void Lorentzian::functionLocal(double *out, const double *xValues, const size_t nData) const {
72 const double amplitude = getParameter("Amplitude");
73 const double peakCentre = getParameter("PeakCentre");
74 const double halfGamma = 0.5 * getParameter("FWHM");
75
76 const double invPI = 1.0 / M_PI;
77 for (size_t i = 0; i < nData; i++) {
78 double diff = (xValues[i] - peakCentre);
79 out[i] = amplitude * invPI * halfGamma / (diff * diff + (halfGamma * halfGamma));
80 }
81}
82
83void Lorentzian::functionDerivLocal(Jacobian *out, const double *xValues, const size_t nData) {
84 const double amplitude = getParameter("Amplitude");
85 const double peakCentre = getParameter("PeakCentre");
86 const double gamma = getParameter("FWHM");
87 const double halfGamma = 0.5 * gamma;
88
89 const double invPI = 1.0 / M_PI;
90 for (size_t i = 0; i < nData; i++) {
91 double diff = xValues[i] - peakCentre;
92 const double invDen1 = 1.0 / (gamma * gamma + 4.0 * diff * diff);
93 const double dfda = 2.0 * invPI * gamma * invDen1;
94 out->set(i, 0, dfda);
95
96 double invDen2 = 1 / (diff * diff + halfGamma * halfGamma);
97 const double dfdxo = amplitude * invPI * gamma * diff * invDen2 * invDen2;
98 out->set(i, 1, dfdxo);
99
100 const double dfdg = -2.0 * amplitude * invPI * (gamma * gamma - 4.0 * diff * diff) * invDen1 * invDen1;
101 out->set(i, 2, dfdg);
102 }
103}
104
112void Lorentzian::histogram1D(double *out, double left, const double *right, const size_t nBins) const {
113
114 const double amplitude = getParameter("Amplitude");
115 const double peakCentre = getParameter("PeakCentre");
116 const double gamma = getParameter("FWHM");
117 const double halfGamma = 0.5 * gamma;
118
119 auto cumulFun = [halfGamma, peakCentre](double x) { return atan((x - peakCentre) / halfGamma) / M_PI; };
120 double cLeft = cumulFun(left);
121 for (size_t i = 0; i < nBins; ++i) {
122 double cRight = cumulFun(right[i]);
123 out[i] = amplitude * (cRight - cLeft);
124 cLeft = cRight;
125 }
126}
127
134void Lorentzian::histogramDerivative1D(Jacobian *jacobian, double left, const double *right, const size_t nBins) const {
135 const double amplitude = getParameter("Amplitude");
136 const double c = getParameter("PeakCentre");
137 const double g = getParameter("FWHM");
138 const double g2 = g * g;
139
140 auto cumulFun = [g, c](double x) { return atan((x - c) / g * 2) / M_PI; };
141 auto denom = [g2, c](double x) { return (g2 + 4 * pow(c - x, 2)) * M_PI; };
142
143 double xl = left;
144 double denomLeft = denom(left);
145 double cLeft = cumulFun(left);
146 for (size_t i = 0; i < nBins; ++i) {
147 double xr = right[i];
148 double denomRight = denom(xr);
149 double cRight = cumulFun(xr);
150 jacobian->set(i, 0, cRight - cLeft);
151 jacobian->set(i, 1, -2.0 * (g / denomRight - g / denomLeft) * amplitude);
152 jacobian->set(i, 2, -2.0 * ((xr - c) / denomRight - (xl - c) / denomLeft) * amplitude);
153 denomLeft = denomRight;
154 cLeft = cRight;
155 xl = xr;
156 }
157}
158
159} // namespace Mantid::CurveFitting::Functions
#define DECLARE_FUNCTION(classname)
Macro for declaring a new type of function to be used with the FunctionFactory.
double left
Definition: LineProfile.cpp:80
double right
Definition: LineProfile.cpp:81
void unfixParameter(const std::string &name)
Free a parameter.
Definition: IFunction.cpp:1522
void fixParameter(const std::string &name, bool isDefault=false)
Fix a parameter.
Definition: IFunction.cpp:1515
void setParameter(const std::string &name, const double &value, bool explicitlySet=true) override
Override parent so that we may bust the cache when a parameter is set.
Represents the Jacobian in IFitFunction::functionDeriv.
Definition: Jacobian.h:22
virtual void set(size_t iY, size_t iP, double value)=0
Set a value to a Jacobian matrix element.
void declareParameter(const std::string &name, double initValue=0, const std::string &description="") override
Declare a new parameter.
double getParameter(size_t i) const override
Get i-th parameter.
Provide lorentzian peak shape function interface to IPeakFunction.
Definition: Lorentzian.h:32
void functionDerivLocal(API::Jacobian *out, const double *xValues, const size_t nData) override
Derivative evaluation method. Default is to calculate numerically.
Definition: Lorentzian.cpp:83
void setHeight(const double h) override
Sets the parameters such that height == h.
Definition: Lorentzian.cpp:41
void functionLocal(double *out, const double *xValues, const size_t nData) const override
Function evaluation method to be implemented in the inherited classes.
Definition: Lorentzian.cpp:71
void fixCentre(bool isDefault=false) override
Fix a parameter or set up a tie such that value returned by centre() is constant during fitting.
Definition: Lorentzian.cpp:63
void histogram1D(double *out, double left, const double *right, const size_t nBins) const override
Calculate histogram data.
Definition: Lorentzian.cpp:112
void unfixIntensity() override
Free the intensity parameter.
Definition: Lorentzian.cpp:69
void histogramDerivative1D(API::Jacobian *jacobian, double left, const double *right, const size_t nBins) const override
Devivatives of the histogram.
Definition: Lorentzian.cpp:134
double height() const override
Returns the height of the function.
Definition: Lorentzian.cpp:32
bool m_amplitudeEqualHeight
When Amplitude is set via setHeight() and fwhm() == 0 height is made equal to Amplitude.
Definition: Lorentzian.h:70
void fixIntensity(bool isDefault=false) override
Fix a parameter or set up a tie such that value returned by intensity() is constant during fitting.
Definition: Lorentzian.cpp:67
void init() override
overwrite IFunction base class method, which declare function parameters
Definition: Lorentzian.cpp:26
void unfixCentre() override
Free the centre parameter.
Definition: Lorentzian.cpp:65
void setFwhm(const double w) override
Sets the parameters such that FWHM = w.
Definition: Lorentzian.cpp:51