Mantid
Loading...
Searching...
No Matches
Bk2BkExpConvPV.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#include <cmath>
8
9#include "MantidAPI/Axis.h"
14#include "MantidKernel/System.h"
16
17#include <gsl/gsl_sf_erf.h>
18#include <gsl/gsl_sf_lambert.h>
19
20using namespace Mantid::Kernel;
21
22using namespace Mantid::API;
23
25
26using namespace CurveFitting;
27using namespace CurveFitting::SpecialFunctionSupport;
28namespace {
30Kernel::Logger g_log("Bk2BkExpConvPV");
31} // namespace
32
33DECLARE_FUNCTION(Bk2BkExpConvPV)
34
35
37void Bk2BkExpConvPV::init() {
38 declareParameter("X0", -0.0, "Location of the peak");
39 declareParameter("Intensity", 0.0, "Integrated intensity");
40 declareParameter("Alpha", 0.04, "Exponential rise");
41 declareParameter("Beta", 0.02, "Exponential decay");
42 declareParameter("Sigma2", 1.0, "Sigma squared");
43 declareParameter("Gamma", 0.0);
44}
45
48void Bk2BkExpConvPV::setHeight(const double h) {
49 setParameter("Intensity", 1);
50 double h0 = height();
51
52 // avoid divide by zero
53 double minCutOff = 100.0 * std::numeric_limits<double>::min();
54 if (h0 >= 0 && h0 < minCutOff)
55 h0 = minCutOff;
56 else if (h0 < 0 && h0 > -minCutOff)
57 h0 = -minCutOff;
58
59 setParameter("Intensity", h / h0);
60}
61
65double Bk2BkExpConvPV::height() const {
66 double height[1];
67 double peakCentre[1];
68 peakCentre[0] = this->getParameter("X0");
69 this->functionLocal(height, peakCentre, 1);
70 return height[0];
71}
72
76double Bk2BkExpConvPV::fwhm() const {
77 // get sigma of Gauss with same FWHM as voigt (H)
78 double H, eta;
79 calHandEta(getParameter("Sigma2"), getParameter("Gamma"), H, eta);
80 const auto s = H / (2 * sqrt(2 * M_LN2)); // FWHM = 2*sqrt(2*ln(2))*sigma
81 const auto w0 = expWidth();
82 // Gaussian and B2B exp widths don't add in quadrature. The following
83 // tends to gaussian at large S and at S=0 is equal to the intrinsic width of
84 // the B2B exp (good to <3% for typical params)
85 // (same Eq. used in BackToBackExponential)
86 return w0 * exp(-0.5 * M_LN2 * s / w0) + 2 * sqrt(2 * M_LN2) * s;
87}
88
94void Bk2BkExpConvPV::setFwhm(const double w) {
95 const auto h0 = height();
96 const auto w0 = expWidth();
97 if (w > w0) {
98 const auto a = 0.5 * M_LN2;
99 const auto b = 2 * sqrt(2 * M_LN2);
100 // Can calculate FWHM of voigt (from FWHM of Gauss componeont from Eq. in BackToBackExponential)
101 // From Eq. in calHandEta assuming FWHM of Gauss and Lorz are equal fwhm_pv = 1.6364*fwhm_gauss
102 const auto fwhm_gauss =
103 b * w0 * (gsl_sf_lambert_W0(-(a / b) * exp(-(a / b) * (w / w0))) / a + (w / w0) / b) / 1.6364;
104 setParameter("Sigma2", std::pow(fwhm_gauss / b, 2));
105 setParameter("Gamma", fwhm_gauss / 2);
106 } else {
107 // set to some small number relative to w0
108 setParameter("Sigma2", 1e-6);
109 setParameter("Gamma", 1e-6);
110 }
111 setHeight(h0);
112}
113
116void Bk2BkExpConvPV::setCentre(const double c) { setParameter("X0", c); }
117
120double Bk2BkExpConvPV::centre() const { return getParameter("X0"); }
121
122void Bk2BkExpConvPV::setMatrixWorkspace(std::shared_ptr<const API::MatrixWorkspace> workspace, size_t wi, double startX,
123 double endX) {
125}
126
129void Bk2BkExpConvPV::functionLocal(double *out, const double *xValues, const size_t nData) const {
130 // 1. Prepare constants
131 const double alpha = this->getParameter("Alpha");
132 const double beta = this->getParameter("Beta");
133 const double sigma2 = this->getParameter("Sigma2");
134 const double gamma = this->getParameter("Gamma");
135 const double intensity = this->getParameter("Intensity");
136 const double x0 = this->getParameter("X0");
137
138 double invert_sqrt2sigma = 1.0 / sqrt(2.0 * sigma2);
139 double N = alpha * beta * 0.5 / (alpha + beta);
140
141 double H, eta;
142 calHandEta(sigma2, gamma, H, eta);
143
144 g_log.debug() << "DB1143: nData = " << nData << " From " << xValues[0] << " To " << xValues[nData - 1]
145 << " X0 = " << x0 << " Intensity = " << intensity << " alpha = " << alpha << " beta = " << beta
146 << " H = " << H << " eta = " << eta << '\n';
147
148 // 2. Do calculation for each data point within extent of peak (avoid NaNs)
149 double extent = 10 * fwhm();
150 for (size_t id = 0; id < nData; ++id) {
151 double dT = xValues[id] - x0;
152 if (fabs(dT) < extent) {
153 double omega = calOmega(dT, eta, N, alpha, beta, H, sigma2, invert_sqrt2sigma);
154 out[id] = intensity * omega;
155 } else {
156 out[id] = 0.0;
157 }
158 }
159}
160
163void Bk2BkExpConvPV::functionDerivLocal(API::Jacobian * /*jacobian*/, const double * /*xValues*/,
164 const size_t /*nData*/) {
165 throw Mantid::Kernel::Exception::NotImplementedError("functionDerivLocal is not implemented for Bk2BkExpConvPV.");
166}
167
171 calNumericalDeriv(domain, jacobian);
172}
173
176double Bk2BkExpConvPV::calOmega(double x, double eta, double N, double alpha, double beta, double H, double sigma2,
177 double invert_sqrt2sigma) const {
178 // 1. Prepare
179 std::complex<double> p(alpha * x, alpha * H * 0.5);
180 std::complex<double> q(-beta * x, beta * H * 0.5);
181
182 double u = 0.5 * alpha * (alpha * sigma2 + 2 * x);
183 double y = (alpha * sigma2 + x) * invert_sqrt2sigma;
184
185 double v = 0.5 * beta * (beta * sigma2 - 2 * x);
186 double z = (beta * sigma2 - x) * invert_sqrt2sigma;
187
188 // 2. Calculate
189 double omega1 = (1 - eta) * N * (exp(u + gsl_sf_log_erfc(y)) + exp(v + gsl_sf_log_erfc(z)));
190 double omega2;
191 if (eta < 1.0E-8) {
192 omega2 = 0.0;
193 } else {
194 omega2 = 2 * N * eta / M_PI * (imag(exponentialIntegral(p)) + imag(exponentialIntegral(q)));
195 }
196 double omega = omega1 - omega2;
197
198 return omega;
199}
200
201void Bk2BkExpConvPV::geneatePeak(double *out, const double *xValues, const size_t nData) {
202 this->functionLocal(out, xValues, nData);
203}
204
205void Bk2BkExpConvPV::calHandEta(double sigma2, double gamma, double &H, double &eta) const {
206 // 1. Calculate H
207 double H_G = sqrt(8.0 * sigma2 * M_LN2); // FWHM Gauss
208 double H_L = 2 * gamma; // FWHM lorz
209
210 double temp1 = std::pow(H_L, 5) + 0.07842 * H_G * std::pow(H_L, 4) + 4.47163 * std::pow(H_G, 2) * std::pow(H_L, 3) +
211 2.42843 * std::pow(H_G, 3) * std::pow(H_L, 2) + 2.69269 * std::pow(H_G, 4) * H_L + std::pow(H_G, 5);
212
213 H = std::pow(temp1, 0.2); // FWHM of PV
214
215 // 2. Calculate eta
216 double gam_pv = H_L / H;
217 eta = 1.36603 * gam_pv - 0.47719 * std::pow(gam_pv, 2) + 0.11116 * std::pow(gam_pv, 3);
218
219 if (eta > 1 || eta < 0) {
220 g_log.error() << "Bk2BkExpConvPV: Calculated eta = " << eta << " is out of range [0, 1].\n";
221 }
222}
223
228 const double a = getParameter("Alpha");
229 const double b = getParameter("Beta");
230 return M_LN2 * (a + b) / (a * b);
231}
232
233} // namespace Mantid::CurveFitting::Functions
#define DECLARE_FUNCTION(classname)
Macro for declaring a new type of function to be used with the FunctionFactory.
IPeaksWorkspace_sptr workspace
Definition: IndexPeaks.cpp:114
#define fabs(x)
Definition: Matrix.cpp:22
Base class that represents the domain of a function.
void setMatrixWorkspace(std::shared_ptr< const API::MatrixWorkspace > workspace, size_t wi, double startX, double endX) override
Set MatrixWorkspace.
Definition: IFunctionMW.cpp:22
virtual double getParameter(size_t i) const =0
Get i-th parameter.
void calNumericalDeriv(const FunctionDomain &domain, Jacobian &jacobian)
Calculate numerical derivatives.
Definition: IFunction.cpp:1031
virtual double intensity() const
Returns the integral intensity of the peak.
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
Bk2BkExpConvPV : Peak profile as tback-to-back exponential convoluted with pseudo-Voigt.
void functionDeriv(const API::FunctionDomain &domain, API::Jacobian &jacobian) override
Numerical derivative.
double height() const override
Get approximate height (maximum value) of the peak (not at X0)
double centre() const override
overwrite IPeakFunction base class methods
void geneatePeak(double *out, const double *xValues, const size_t nData)
Set up the range of peak calculation for higher efficiency.
double fwhm() const override
Get approximate peak width.
void setHeight(const double h) override
Set peak height.
void setMatrixWorkspace(std::shared_ptr< const API::MatrixWorkspace > workspace, size_t wi, double startX, double endX) override
Set MatrixWorkspace.
double calOmega(double x, double eta, double N, double alpha, double beta, double H, double sigma2, double invert_sqrt2sigma) const
Calculate Omega(x) = ... ...
void setFwhm(const double w) override
Set new peak width approximately using same mapping in BackToBackExponential assuming fwhm_gaus = fwh...
void setCentre(const double c) override
Set peak center.
void functionLocal(double *out, const double *xValues, const size_t nData) const override
Implement the peak calculating formula.
double expWidth() const
Calculate contribution to the width by the exponentials.
void calHandEta(double sigma2, double gamma, double &H, double &eta) const
void functionDerivLocal(API::Jacobian *out, const double *xValues, const size_t nData) override
Local derivative.
Marks code as not implemented yet.
Definition: Exception.h:138
The Logger class is in charge of the publishing messages from the framework through various channels.
Definition: Logger.h:52
void debug(const std::string &msg)
Logs at debug level.
Definition: Logger.cpp:114
void error(const std::string &msg)
Logs at error level.
Definition: Logger.cpp:77
Kernel::Logger g_log("ExperimentInfo")
static logger object
std::complex< double > MANTID_CURVEFITTING_DLL exponentialIntegral(const std::complex< double > &z)
Compute exp(z)*E1(z) where z is complex and E1(z) is the Exponential Integral.