Mantid
Loading...
Searching...
No Matches
NeutronBk2BkExpConvPVoigt.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 +
11#include "MantidKernel/Logger.h"
13
15
16#include <gsl/gsl_sf_erf.h>
17
18const double PEAKRANGE = 5.0;
19const double NEG_DBL_MAX = -1. * DBL_MAX;
20
21using namespace std;
22
24
25using namespace CurveFitting;
26namespace {
28Kernel::Logger g_log("NeutronBk2BkExpConvPV");
29} // namespace
30
31DECLARE_FUNCTION(NeutronBk2BkExpConvPVoigt)
32
33//----------------------------------------------------------------------------------------------
37 : API::IPowderDiffPeakFunction(), m_Alpha(), m_Beta(), m_Sigma2(), m_Gamma(), m_eta(), m_N() {
38 mHKLSet = false;
39}
40
43
44//----------------------------------------------------------------------------------------------
50 // Peak height (0)
51 declareParameter("Height", 1.0, "Intensity of peak");
52
53 // Instrument geometry related (1 ~ 3)
54 declareParameter("Dtt1", 1.0, "coefficient 1 for d-spacing calculation for epithermal neutron part");
55 declareParameter("Dtt2", 1.0, "coefficient 2 for d-spacing calculation for epithermal neutron part");
56 declareParameter("Zero", 0.0, "Zero shift for epithermal neutron");
57
58 // Peak profile related (4 ~ 7) Back to back Expoential
59 declareParameter("Alph0", 1.6, "exponential constant for rising part of epithermal neutron pulse");
60 declareParameter("Alph1", 1.5, "exponential constant for rising part of expithermal neutron pulse");
61 declareParameter("Beta0", 1.6, "exponential constant of decaying part of epithermal neutron pulse");
62 declareParameter("Beta1", 1.5, "exponential constant of decaying part of epithermal neutron pulse");
63
64 // Pseudo-Voigt (8 ~ 13)
65 declareParameter("Sig0", 1.0,
66 "variance parameter 1 of the Gaussian "
67 "component of the psuedovoigt function");
68 declareParameter("Sig1", 1.0,
69 "variance parameter 2 of the Gaussian "
70 "component of the psuedovoigt function");
71 declareParameter("Sig2", 1.0,
72 "variance parameter 3 of the Gaussian "
73 "component of the psuedovoigt function");
74
75 declareParameter("Gam0", 0.0,
76 "FWHM parameter 1 of the Lorentzian component "
77 "of the psuedovoigt function");
78 declareParameter("Gam1", 0.0,
79 "FWHM parameter 2 of the Lorentzian component "
80 "of the psuedovoigt function");
81 declareParameter("Gam2", 0.0,
82 "FWHM parameter 3 of the Lorentzian component "
83 "of the psuedovoigt function");
84
85 // Lattice parameter (14)
86 declareParameter("LatticeConstant", 10.0, "lattice constant for the sample");
87
88 declareParameter("H", 0.0, "Miller index H. ");
89 declareParameter("K", 0.0, "Miller index K. ");
90 declareParameter("L", 0.0, "Miller index L. ");
91
92 // Some build-in constant
93 LATTICEINDEX = 14;
94 HEIGHTINDEX = 0;
95
96 // Unit cell
97 m_unitCellSize = -DBL_MAX;
98
99 // Set flag
101}
102
103//----------------------------------------------------------------------------------------------
110double NeutronBk2BkExpConvPVoigt::getPeakParameter(const std::string &paramname) {
111 // Calculate peak parameters if required
113 calculateParameters(false);
114 }
115
116 // Get value
117 double paramvalue(EMPTY_DBL());
118
119 if (paramname == "Alpha")
120 paramvalue = m_Alpha;
121 else if (paramname == "Beta")
122 paramvalue = m_Beta;
123 else if (paramname == "Sigma2")
124 paramvalue = m_Sigma2;
125 else if (paramname == "Gamma")
126 paramvalue = m_Gamma;
127 else if (paramname == "d_h")
128 paramvalue = m_dcentre;
129 else if (paramname == "Eta")
130 paramvalue = m_eta;
131 else if (paramname == "TOF_h")
132 paramvalue = m_centre;
133 else if (paramname == "FWHM")
134 paramvalue = m_fwhm;
135 else {
136 stringstream errss;
137 errss << "Parameter " << paramname << " does not exist in peak function " << this->name()
138 << "'s calculated parameters. "
139 << "Candidates are Alpha, Beta, Sigma2, Gamma, d_h and FWHM. ";
140 g_log.error(errss.str());
141 throw runtime_error(errss.str());
142 }
143
144 return paramvalue;
145}
146
147//----------------------------------------------------------------------------------------------
151void NeutronBk2BkExpConvPVoigt::calculateParameters(bool explicitoutput) const {
152 // Obtain parameters (class) with pre-set order
153 double dtt1 = getParameter(1);
154 double dtt2 = getParameter(2);
155 double zero = getParameter(3);
156
157 double alph0 = getParameter(4);
158 double alph1 = getParameter(5);
159 double beta0 = getParameter(6);
160 double beta1 = getParameter(7);
161
162 double sig0 = getParameter(8);
163 double sig1 = getParameter(9);
164 double sig2 = getParameter(10);
165 double gam0 = getParameter(11);
166 double gam1 = getParameter(12);
167 double gam2 = getParameter(13);
168
169 double latticeconstant = getParameter(LATTICEINDEX);
170
171 if (!mHKLSet) {
172 // Miller index can only be set once. Either via parameters or
173 // client-called setMillerIndex
174 double h = getParameter(15);
175 double k = getParameter(16);
176 double l = getParameter(17);
177 mH = static_cast<int>(h);
178 mK = static_cast<int>(k);
179 mL = static_cast<int>(l);
180
181 // Check value valid or not
182 if (mH * mH + mK * mK + mL * mL < 1.0E-8) {
183 stringstream errmsg;
184 errmsg << "H = K = L = 0 is not allowed";
185 g_log.error(errmsg.str());
186 throw std::invalid_argument(errmsg.str());
187 }
188
189 g_log.debug() << "Set (HKL) from input parameter = " << mH << ", " << mK << ", " << mL << "\n";
190
191 mHKLSet = true;
192 }
193
194 double dh, tof_h, eta, alpha, beta, H, sigma2, gamma, N;
195
196 // Calcualte Peak Position d-spacing if necessary
198 // TODO : make it more general to all type of lattice structure
199 m_unitCell.set(latticeconstant, latticeconstant, latticeconstant, 90.0, 90.0, 90.0);
200 dh = m_unitCell.d(mH, mK, mL);
201 m_dcentre = dh;
203 } else {
204 dh = m_dcentre;
205 }
206
207 // Calculate all the parameters
208 /*
209 alpha(d) = alpha0 + alpha1/d_h
210 beta(d) = beta0 + beta1/d_h^4
211 tof(d) = zero + Dtt1*d_h + Dtt2*d_h^2
212 */
213
214 alpha = alph0 + alph1 / dh;
215 beta = beta0 + beta1 / pow(dh, 4.);
216 tof_h = zero + dtt1 * dh + dtt2 * dh * dh;
217
218 sigma2 = sig0 * sig0 + sig1 * sig1 * std::pow(dh, 2) + sig2 * sig2 * std::pow(dh, 4);
219 gamma = gam0 + gam1 * dh + gam2 * std::pow(dh, 2);
220
221 g_log.debug() << "[F001] TOF_h = " << tof_h << ", Alpha = " << alpha << ", Beta = " << beta << ", Gamma = " << gamma
222 << "(Gam-0 = " << gam0 << ", Gam-1 = " << gam1 << ", Gam-2 = " << gam2 << ")."
223 << "\n";
224
225 // Calcualte H for the peak
226 calHandEta(sigma2, gamma, H, eta);
227
228 N = alpha * beta * 0.5 / (alpha + beta);
229
230 // Record recent value
231 m_Alpha = alpha;
232 m_Beta = beta;
233 m_Sigma2 = sigma2;
234 m_Gamma = gamma;
235 m_fwhm = H;
236 m_centre = tof_h;
237 m_N = N;
238 m_eta = eta;
239
240 // Check whether all the parameters are physical
241 if (alpha != alpha || beta != beta || sigma2 != sigma2 || gamma != gamma || H != H || H <= 0.) {
242 m_parameterValid = false;
243 } else {
244 m_parameterValid = true;
245 }
246
247 // Debug output
248 if (explicitoutput) {
249 stringstream msgss;
250 msgss << " dh = " << dh << "; TOF = " << tof_h << ", "
251 << "alpha = " << alpha << ", beta = " << beta;
252 g_log.information(msgss.str());
253 }
254
255 // Reset the flag
257}
258
259//----------------------------------------------------------------------------------------------
262void NeutronBk2BkExpConvPVoigt::setParameter(size_t i, const double &value, bool explicitlySet) {
263 if (i == LATTICEINDEX) {
264 // Lattice parameter
265 if (fabs(m_unitCellSize - value) > 1.0E-8) {
266 // If change in value is non-trivial
268 ParamFunction::setParameter(i, value, explicitlySet);
271 }
272 } else {
273 // Non lattice parameter
274 ParamFunction::setParameter(i, value, explicitlySet);
276 }
277}
278
279//----------------------------------------------------------------------------------------------
282void NeutronBk2BkExpConvPVoigt::setParameter(const std::string &name, const double &value, bool explicitlySet) {
283 if (name == "LatticeConstant") {
284 // Lattice parameter
285 if (fabs(m_unitCellSize - value) > 1.0E-8) {
286 // If change in value is non-trivial
291 }
292 } else {
293 ParamFunction::setParameter(name, value, explicitlySet);
295 }
296}
297
298//----------------------------------------------------------------------------------------------
305void NeutronBk2BkExpConvPVoigt::function(vector<double> &out, const vector<double> &xValues) const {
306 // Calculate peak parameters
308 calculateParameters(false);
309 else
310 g_log.debug("Function() has no new parameters to calculate. ");
311
312 // Get some other parmeters
313 const double HEIGHT = getParameter(HEIGHTINDEX);
314 const double INVERT_SQRT2SIGMA = 1.0 / sqrt(2.0 * m_Sigma2);
315
316 // Calculate where to start and to end calculating
317 const double RANGE = m_fwhm * PEAKRANGE;
318
319 const double LEFT_VALUE = m_centre - RANGE;
320 auto iter = std::lower_bound(xValues.cbegin(), xValues.cend(), LEFT_VALUE);
321
322 const double RIGHT_VALUE = m_centre + RANGE;
323 auto iter_end = std::lower_bound(iter, xValues.cend(), RIGHT_VALUE);
324
325 // Calcualte
326 std::size_t pos(std::distance(xValues.begin(), iter)); // second loop variable
327 for (; iter != iter_end; ++iter) {
328 out[pos] = HEIGHT * calOmega(*iter - m_centre, m_eta, m_N, m_Alpha, m_Beta, m_fwhm, m_Sigma2, INVERT_SQRT2SIGMA);
329 pos++;
330 } // ENDFOR data points
331}
332
333//----------------------------------------------------------------------------------------------
336void NeutronBk2BkExpConvPVoigt::function1D(double *out, const double *xValues, const size_t nData) const {
337 // Calculate peak parameters
339 calculateParameters(false);
340 else
341 g_log.debug("Function() has no new parameters to calculate. ");
342
343 // Get some other parameters
344 const double HEIGHT = getParameter(HEIGHTINDEX);
345 const double INVERT_SQRT2SIGMA = 1.0 / sqrt(2.0 * m_Sigma2);
346
347 // On each data point
348 const double RANGE = m_fwhm * PEAKRANGE;
349 g_log.debug() << "[F002] Peak centre = " << m_centre << "; Calcualtion Range = " << RANGE << ".\n";
350
351 for (size_t i = 0; i < nData; ++i) {
352 if (fabs(xValues[i] - m_centre) < RANGE) {
353 // In peak range
354 out[i] =
355 HEIGHT * calOmega(xValues[i] - m_centre, m_eta, m_N, m_Alpha, m_Beta, m_fwhm, m_Sigma2, INVERT_SQRT2SIGMA);
356 g_log.debug() << "TOF = " << xValues[i] << " = " << out[i] << "\n";
357 } else {
358 // Out of peak range
359 out[i] = 0.0;
360 g_log.debug() << "TOF = " << xValues[i] << " out of calculation range. "
361 << ".\n";
362 }
363 }
364}
365
366//----------------------------------------------------------------------------------------------
369void NeutronBk2BkExpConvPVoigt::calHandEta(double sigma2, double gamma, double &H, double &eta) const {
370 // 1. Calculate H
371 double H_G = sqrt(8.0 * sigma2 * M_LN2);
372 double H_L = gamma;
373
374 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) +
375 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);
376
377 H = std::pow(temp1, 0.2);
378
379 // 2. Calculate eta
380 double gam_pv = H_L / H;
381 eta = 1.36603 * gam_pv - 0.47719 * std::pow(gam_pv, 2) + 0.11116 * std::pow(gam_pv, 3);
382
383 if (eta > 1 || eta < 0) {
384 g_log.warning() << "Calculated eta = " << eta << " is out of range [0, 1].\n";
385 } else {
386 g_log.debug() << "[DBx121] Eta = " << eta << "; Gamma = " << gamma << ".\n";
387 }
388}
389
390//----------------------------------------------------------------------------------------------
394double NeutronBk2BkExpConvPVoigt::calOmega(const double x, const double eta, const double N, const double alpha,
395 const double beta, const double H, const double sigma2,
396 const double invert_sqrt2sigma, const bool explicitoutput) const {
397 // Transform to variable u, v, y, z
398 const double u = 0.5 * alpha * (alpha * sigma2 + 2. * x);
399 const double y = (alpha * sigma2 + x) * invert_sqrt2sigma;
400
401 const double v = 0.5 * beta * (beta * sigma2 - 2. * x);
402 const double z = (beta * sigma2 - x) * invert_sqrt2sigma;
403
404 // Calculate Gaussian part
405 const double erfcy = gsl_sf_erfc(y);
406 double part1(0.);
407 if (fabs(erfcy) > DBL_MIN)
408 part1 = exp(u) * erfcy;
409
410 const double erfcz = gsl_sf_erfc(z);
411 double part2(0.);
412 if (fabs(erfcz) > DBL_MIN)
413 part2 = exp(v) * erfcz;
414
415 const double omega1 = (1. - eta) * N * (part1 + part2);
416
417 // Calculate Lorenzian part
418 double omega2(0.);
419
420 g_log.debug() << "Eta = " << eta << "; X = " << x << " N = " << N << ".\n";
421
422 if (eta >= 1.0E-8) {
423 const double SQRT_H_5 = sqrt(H) * .5;
424 std::complex<double> p(alpha * x, alpha * SQRT_H_5);
425 std::complex<double> q(-beta * x, beta * SQRT_H_5);
426 double omega2a = imag(exp(p) * Mantid::API::E1(p));
427 double omega2b = imag(exp(q) * Mantid::API::E1(q));
428 omega2 = -1.0 * N * eta * (omega2a + omega2b) * M_2_PI;
429
430 g_log.debug() << "Exp(p) = " << exp(p) << ", Exp(q) = " << exp(q) << ".\n";
431
432 if (omega2 != omega2) {
433 g_log.debug() << "Omega2 is not physical. Omega2a = " << omega2a << ", Omega2b = " << omega2b
434 << ", p = " << p.real() << ", " << p.imag() << ".\n";
435 } else {
436 g_log.debug() << "X = " << x << " is OK. Omega 2 = " << omega2 << ", Omega2A = " << omega2a
437 << ", Omega2B = " << omega2b << "\n";
438 }
439 }
440 const double omega = omega1 + omega2;
441
442 if (explicitoutput || omega != omega) {
443 if (omega <= NEG_DBL_MAX || omega >= DBL_MAX || omega != omega) {
444 stringstream errss;
445 errss << "Peak (" << mH << mK << mL << "): TOF = " << m_centre << ", dX = " << x << ", (" << x / m_fwhm
446 << " FWHM) ";
447 errss << "Omega = " << omega << " is infinity! omega1 = " << omega1 << ", omega2 = " << omega2 << "\n";
448 errss << " u = " << u << ", v = " << v << ", erfc(y) = " << gsl_sf_erfc(y) << ", erfc(z) = " << gsl_sf_erfc(z)
449 << "\n";
450 errss << " alpha = " << alpha << ", beta = " << beta << " sigma2 = " << sigma2 << ", N = " << N << "\n";
451 g_log.warning(errss.str());
452 }
453 }
454 g_log.debug() << "[DB] Final Value of Omega = " << omega << ".\n";
455
456 return omega;
457}
458
459} // namespace Mantid::CurveFitting::Functions
double value
The value of the point.
Definition: FitMW.cpp:51
#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
const double PEAKRANGE
const double NEG_DBL_MAX
static Kernel::Logger g_log
Logger instance.
Definition: IFunction1D.h:74
An interface to a peak function, which extend the interface of IFunctionWithLocation by adding method...
bool m_cellParamValueChanged
An indicator to re-calculate peak d-space position.
bool m_parameterValid
Flag to indicate whether peaks' parameters value can generate a valid peak.
bool m_hasNewParameterValue
Flag if any parameter value changed.
double m_dcentre
Centre of the peak in d-space.
void setParameter(size_t, const double &value, bool explicitlySet=true) override
Set i-th parameter.
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.
NeutronBk2BkExpConvPVoigt : Back-to-back exponential function convoluted with pseudo-voigt for epithe...
void calculateParameters(bool explicitoutput) const override
Calculate peak parameters (alpha, beta, sigma2..)
void setParameter(size_t i, const double &value, bool explicitlySet=true) override
Override setting a new value to the i-th parameter.
void function(std::vector< double > &out, const std::vector< double > &xValues) const override
Function (local) of the vector version.
static int s_peakRadius
Default value for the peak radius.
double getPeakParameter(const std::string &) override
Get peak parameters.
std::string name() const override
Overwrite IFunction base class method: name.
void calHandEta(double sigma2, double gamma, double &H, double &eta) const
Calcualte H and Eta.
double m_Alpha
Set 2 functions to be hidden from client.
void function1D(double *out, const double *xValues, const size_t nData) const override
Function you want to fit to.
double calOmega(const double x, const double eta, const double N, const double alpha, const double beta, const double H, const double sigma2, const double invert_sqrt2sigma, const bool explicitoutput=false) const
Calculate peak profile I(TOF) = Omega(TOF)
void set(double _a, double _b, double _c, double _alpha, double _beta, double _gamma, const int angleunit=angDegrees)
Set lattice parameters.
Definition: UnitCell.cpp:303
double d(double h, double k, double l) const
Return d-spacing ( ) for a given h,k,l coordinate.
Definition: UnitCell.cpp:700
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
void warning(const std::string &msg)
Logs at warning level.
Definition: Logger.cpp:86
void information(const std::string &msg)
Logs at information level.
Definition: Logger.cpp:105
std::complex< double > MANTID_API_DLL E1(std::complex< double > z)
Integral for Gamma.
Kernel::Logger g_log("ExperimentInfo")
static logger object
constexpr double EMPTY_DBL() noexcept
Returns what we consider an "empty" double within a property.
Definition: EmptyValues.h:43
STL namespace.