Mantid
Loading...
Searching...
No Matches
BackToBackExponential.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
13#include <cmath>
14#include <gsl/gsl_multifit_nlin.h>
15#include <gsl/gsl_sf_erf.h>
16#include <gsl/gsl_sf_lambert.h>
17
19
20using namespace CurveFitting;
21
22using namespace Kernel;
23
24using namespace API;
25
26DECLARE_FUNCTION(BackToBackExponential)
27
29 // Do not change the order of these parameters!
30 declareParameter("I", 0.0, "integrated intensity of the peak"); // 0
31 declareParameter("A", 1.0,
32 "exponential constant of rising part of neutron pulse"); // 1
33 declareParameter("B", 0.05, "exponential constant of decaying part of neutron pulse"); // 2
34 declareParameter("X0", 0.0, "peak position"); // 3
35 declareParameter("S", 1.0,
36 "standard deviation of gaussian part of peakshape function"); // 4
37}
38
43 // height = I * ln2 /(FWHM - ln2 * S)
44 const double I = getParameter(0);
45 const double s = getParameter(4);
46 const auto width = fwhm();
47 const auto h = I * M_LN2 / (width - M_LN2 * s);
48 return h;
49}
50
56 double h0 = height();
57 if (h0 == 0.0) {
58 setParameter(0, 1e-6);
59 h0 = height();
60 }
61 double area = getParameter(0); // ==I
62 area *= h / h0;
63 if (area <= 0.0) {
64 area = 1e-6;
65 }
66 if (!std::isfinite(area)) {
67 area = std::numeric_limits<double>::max() / 2;
68 }
69 setParameter(0, area);
70}
71
76 // intrinsic width of B2B exp
77 auto s = getParameter("S");
78 auto w0 = expWidth();
79 // Gaussian and B2B exp widths don't add in quadrature. The following
80 // tends to gaussian at large S and at S=0 is equal to the intrinsic width of
81 // the B2B exp (good to <3% for typical params)
82 return w0 * exp(-0.5 * M_LN2 * s / w0) + 2 * sqrt(2 * M_LN2) * s;
83}
84
89void BackToBackExponential::setFwhm(const double w) {
90 // get height so can reset after changing S
91 const auto h0 = height();
92 const auto w0 = expWidth();
93 if (w > w0) {
94 const auto a = 0.5 * M_LN2;
95 const auto b = 2 * sqrt(2 * M_LN2);
96 // calculate new value of S (from solving eq in fwhm func)
97 setParameter("S", w0 * (gsl_sf_lambert_W0(-(a / b) * exp(-(a / b) * (w / w0))) / a + (w / w0) / b));
98 } else {
99 // set to some small number relative to w0
100 setParameter("S", 1e-6);
101 }
102 // reset height
103 setHeight(h0);
104}
105
106void BackToBackExponential::function1D(double *out, const double *xValues, const size_t nData) const {
107 /*
108 const double& I = getParameter("I");
109 const double& a = getParameter("A");
110 const double& b = getParameter("B");
111 const double& x0 = getParameter("X0");
112 const double& s = getParameter("S");
113 */
114
115 const double I = getParameter(0);
116 const double a = getParameter(1);
117 const double b = getParameter(2);
118 const double x0 = getParameter(3);
119 const double s = getParameter(4);
120
121 // find the reasonable extent of the peak ~100 fwhm
122 double extent = expWidth();
123 if (s > extent)
124 extent = s;
125 extent *= 100;
126
127 double s2 = s * s;
128 double normFactor = a * b / (a + b) / 2;
129 // Needed for IntegratePeaksMD for cylinder profile fitted with b=0
130 if (normFactor == 0.0)
131 normFactor = 1.0;
132 for (size_t i = 0; i < nData; i++) {
133 double diff = xValues[i] - x0;
134 if (fabs(diff) < extent) {
135 double val = 0.0;
136 double arg1 = a / 2 * (a * s2 + 2 * diff);
137 val += exp(arg1 + gsl_sf_log_erfc((a * s2 + diff) / sqrt(2 * s2))); // prevent overflow
138 double arg2 = b / 2 * (b * s2 - 2 * diff);
139 val += exp(arg2 + gsl_sf_log_erfc((b * s2 - diff) / sqrt(2 * s2))); // prevent overflow
140 out[i] = I * val * normFactor;
141 } else
142 out[i] = 0.0;
143 }
144}
145
149void BackToBackExponential::functionDeriv1D(Jacobian *jacobian, const double *xValues, const size_t nData) {
150 FunctionDomain1DView domain(xValues, nData);
151 this->calNumericalDeriv(domain, *jacobian);
152}
153
158 const double a = getParameter(1);
159 const double b = getParameter(2);
160 // Needed for IntegratePeaksMD for cylinder profile fitted with b=0
161 if (a * b == 0.0)
162 return M_LN2;
163 return M_LN2 * (a + b) / (a * b);
164}
165
166} // namespace Mantid::CurveFitting::Functions
#define DECLARE_FUNCTION(classname)
Macro for declaring a new type of function to be used with the FunctionFactory.
#define fabs(x)
Definition: Matrix.cpp:22
1D domain - a wrapper around an array of doubles.
void calNumericalDeriv(const FunctionDomain &domain, Jacobian &jacobian)
Calculate numerical derivatives.
Definition: IFunction.cpp:1031
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
double getParameter(size_t i) const override
Get i-th parameter.
Provide BackToBackExponential peak shape function interface to IPeakFunction.
void setFwhm(const double w) override
Set new peak width approximately.
void function1D(double *out, const double *xValues, const size_t nData) const override
General implementation of the method for all peaks.
double height() const override
Get approximate height (maximum value) of the peak (not at X0)
void setHeight(const double h) override
Set new height of the peak.
double fwhm() const override
Get approximate peak width.
void functionDeriv1D(API::Jacobian *jacobian, const double *xValues, const size_t nData) override
Evaluate function derivatives numerically.
double expWidth() const
Calculate contribution to the width by the exponentials.