Mantid
Loading...
Searching...
No Matches
IkedaCarpenterPV.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//----------------------------------------------------------------------
22
23#include <cmath>
24#include <gsl/gsl_math.h>
25#include <gsl/gsl_multifit_nlin.h>
26#include <gsl/gsl_sf_erf.h>
27#include <limits>
28
30
31using namespace CurveFitting;
32
33namespace {
35Kernel::Logger g_log("IkedaCarpenterPV");
36} // namespace
37
38using namespace Kernel;
39
40using namespace CurveFitting::SpecialFunctionSupport;
41
42using namespace Geometry;
43using namespace Constraints;
44
45DECLARE_FUNCTION(IkedaCarpenterPV)
46
47double IkedaCarpenterPV::centre() const { return getParameter("X0"); }
48
49void IkedaCarpenterPV::setHeight(const double h) {
50 // calculate height of peakshape function corresponding to intensity = 1
51 setParameter("I", 1);
52 double h0 = height();
53
54 // to avoid devide by zero and to insane value for I to be set
55 double minCutOff = 100.0 * std::numeric_limits<double>::min();
56 if (h0 > 0 && h0 < minCutOff)
57 h0 = minCutOff;
58 if (h0 < 0 && h0 > -minCutOff)
59 h0 = -minCutOff;
60
61 // The intensity is then estimated to be h/h0
62 setParameter("I", h / h0);
63}
64
66 // return the function value at centre()
67 // using arrays - otherwise coverity warning
68 double h0[1];
69 double toCentre[1];
70 toCentre[0] = centre();
71 constFunction(h0, toCentre, 1);
72 return h0[0];
73}
74
75double IkedaCarpenterPV::fwhm() const {
76 double sigmaSquared = getParameter("SigmaSquared");
77 double gamma = getParameter("Gamma");
78
79 if (sigmaSquared < 0) {
80 g_log.debug() << "SigmaSquared NEGATIVE!.\n"
81 << "Likely due to a fit not converging properly\n"
82 << "If this is frequent problem please report to Mantid team.\n"
83 << "For now to calculate width force SigmaSquared positive.\n";
84 sigmaSquared = -sigmaSquared;
85 }
86 if (gamma < 0) {
87 g_log.debug() << "Gamma NEGATIVE!.\n"
88 << "Likely due to a fit not converging properly\n"
89 << "If this is frequent problem please report to Mantid team.\n"
90 << "For now to calculate width force Gamma positive.\n";
91 gamma = -gamma;
92 ;
93 }
94 return sqrt(8.0 * M_LN2 * sigmaSquared) + gamma;
95}
96
97void IkedaCarpenterPV::setFwhm(const double w) {
98 setParameter("SigmaSquared", w * w / (32.0 * M_LN2)); // used 4.0 * 8.0 = 32.0
99 setParameter("Gamma", w / 2.0);
100}
101
102void IkedaCarpenterPV::setCentre(const double c) { setParameter("X0", c); }
103
105 declareParameter("I", 0.0,
106 "The integrated intensity of the peak. I.e. "
107 "approximately equal to HWHM times height of "
108 "peak");
109 this->lowerConstraint0("I");
110 declareParameter("Alpha0", 1.6, "Used to model fast decay constant");
111 this->lowerConstraint0("Alpha0");
112 declareParameter("Alpha1", 1.5, "Used to model fast decay constant");
113 this->lowerConstraint0("Alpha1");
114 declareParameter("Beta0", 31.9, "Inverse of slow decay constant");
115 this->lowerConstraint0("Beta0");
116 declareParameter("Kappa", 46.0, "Controls contribution of slow decay term");
117 this->lowerConstraint0("Kappa");
118 declareParameter("SigmaSquared", 1.0, "standard deviation squared (Voigt Guassian broadening)");
119 this->lowerConstraint0("SigmaSquared");
120 declareParameter("Gamma", 1.0, "Voigt Lorentzian broadening");
121 this->lowerConstraint0("Gamma");
122 declareParameter("X0", 0.0, "Peak position");
123 this->lowerConstraint0("X0");
124}
125
126void IkedaCarpenterPV::lowerConstraint0(const std::string &paramName) {
127 auto mixingConstraint = std::make_unique<BoundaryConstraint>(this, paramName, 0.0, true);
128 mixingConstraint->setPenaltyFactor(1e9);
129
130 addConstraint(std::move(mixingConstraint));
131}
132
141void IkedaCarpenterPV::calWavelengthAtEachDataPoint(const double *xValues, const size_t &nData) const {
142 // if wavelength vector already have the right size no need for resizing it
143 // further we make the assumption that no need to recalculate this vector if
144 // it already has the right size
145
146 if (m_waveLength.size() != nData) {
147 m_waveLength.resize(nData);
148
149 Mantid::Kernel::Unit_sptr wavelength = Mantid::Kernel::UnitFactory::Instance().create("Wavelength");
150 for (size_t i = 0; i < nData; i++) {
151 m_waveLength[i] = xValues[i];
152 }
153
154 // note if a version of convertValue was added which allows a double* as
155 // first argument then could avoid copying above plus only have to resize
156 // m_wavelength when its size smaller than nData
158 if (mws) {
159 Instrument_const_sptr instrument = mws->getInstrument();
160 Geometry::IComponent_const_sptr sample = instrument->getSample();
161 if (sample != nullptr) {
162 convertValue(m_waveLength, wavelength, mws, m_workspaceIndex);
163 } else {
164 g_log.warning() << "No sample set for instrument in workspace.\n"
165 << "Can't calculate wavelength in IkedaCarpenter.\n"
166 << "Default all wavelengths to one.\n"
167 << "Solution is to load appropriate instrument into workspace.\n";
168 for (size_t i = 0; i < nData; i++)
169 m_waveLength[i] = 1.0;
170 }
171 } else {
172 g_log.warning() << "Workspace not set.\n"
173 << "Can't calculate wavelength in IkedaCarpenter.\n"
174 << "Default all wavelengths to one.\n"
175 << "Solution call setMatrixWorkspace() for function.\n";
176 for (size_t i = 0; i < nData; i++)
177 m_waveLength[i] = 1.0;
178 }
179 }
180}
181
189void IkedaCarpenterPV::convertVoigtToPseudo(const double &voigtSigmaSq, const double &voigtGamma, double &H,
190 double &eta) const {
191 double fwhmGsq = 8.0 * M_LN2 * voigtSigmaSq;
192 double fwhmG = sqrt(fwhmGsq);
193 double fwhmG4 = fwhmGsq * fwhmGsq;
194 double fwhmL = voigtGamma;
195 double fwhmLsq = voigtGamma * voigtGamma;
196 double fwhmL4 = fwhmLsq * fwhmLsq;
197
198 H = pow(fwhmG4 * fwhmG + 2.69269 * fwhmG4 * fwhmL + 2.42843 * fwhmGsq * fwhmG * fwhmLsq +
199 4.47163 * fwhmGsq * fwhmLsq * fwhmL + 0.07842 * fwhmG * fwhmL4 + fwhmL4 * fwhmL,
200 0.2);
201
202 if (H == 0.0)
203 H = std::numeric_limits<double>::epsilon() * 1000.0;
204
205 double tmp = fwhmL / H;
206
207 eta = 1.36603 * tmp - 0.47719 * tmp * tmp + 0.11116 * tmp * tmp * tmp;
208}
209
210void IkedaCarpenterPV::constFunction(double *out, const double *xValues, const int &nData) const {
211 const double I = getParameter("I");
212 const double alpha0 = getParameter("Alpha0");
213 const double alpha1 = getParameter("Alpha1");
214 const double beta0 = getParameter("Beta0");
215 const double kappa = getParameter("Kappa");
216 const double voigtsigmaSquared = getParameter("SigmaSquared");
217 const double voigtgamma = getParameter("Gamma");
218 const double X0 = getParameter("X0");
219
220 // cal pseudo voigt sigmaSq and gamma and eta
221 double gamma = 1.0; // dummy initialization
222 double eta = 0.5; // dummy initialization
223 convertVoigtToPseudo(voigtsigmaSquared, voigtgamma, gamma, eta);
224 double sigmaSquared = gamma * gamma / (8.0 * M_LN2);
225
226 const double beta = 1 / beta0;
227
228 // equations taken from Fullprof manual
229
230 const double k = 0.05;
231
232 // Not entirely sure what to do if sigmaSquared ever negative
233 // for now just post a warning
234 double someConst = std::numeric_limits<double>::max() / 100.0;
235 if (sigmaSquared > 0)
236 someConst = 1 / sqrt(2.0 * sigmaSquared);
237 else if (sigmaSquared < 0) {
238 g_log.warning() << "sigmaSquared negative in functionLocal.\n";
239 }
240
241 // update wavelength vector
242 calWavelengthAtEachDataPoint(xValues, nData);
243
244 for (int i = 0; i < nData; i++) {
245 double diff = xValues[i] - X0;
246
247 double R = exp(-81.799 / (m_waveLength[i] * m_waveLength[i] * kappa));
248 double alpha = 1.0 / (alpha0 + m_waveLength[i] * alpha1);
249
250 double a_minus = alpha * (1 - k);
251 double a_plus = alpha * (1 + k);
252 double x = a_minus - beta;
253 double y = alpha - beta;
254 double z = a_plus - beta;
255
256 double Nu = 1 - R * a_minus / x;
257 double Nv = 1 - R * a_plus / z;
258 double Ns = -2 * (1 - R * alpha / y);
259 double Nr = 2 * R * alpha * alpha * beta * k * k / (x * y * z);
260
261 double u = a_minus * (a_minus * sigmaSquared - 2 * diff) / 2.0;
262 double v = a_plus * (a_plus * sigmaSquared - 2 * diff) / 2.0;
263 double s = alpha * (alpha * sigmaSquared - 2 * diff) / 2.0;
264 double r = beta * (beta * sigmaSquared - 2 * diff) / 2.0;
265
266 double yu = (a_minus * sigmaSquared - diff) * someConst;
267 double yv = (a_plus * sigmaSquared - diff) * someConst;
268 double ys = (alpha * sigmaSquared - diff) * someConst;
269 double yr = (beta * sigmaSquared - diff) * someConst;
270
271 std::complex<double> zs = std::complex<double>(-alpha * diff, 0.5 * alpha * gamma);
272 std::complex<double> zu = (1 - k) * zs;
273 std::complex<double> zv = (1 + k) * zs;
274 std::complex<double> zr = std::complex<double>(-beta * diff, 0.5 * beta * gamma);
275
276 double N = 0.25 * alpha * (1 - k * k) / (k * k);
277
278 out[i] = I * N *
279 ((1 - eta) * (Nu * exp(u + gsl_sf_log_erfc(yu)) + Nv * exp(v + gsl_sf_log_erfc(yv)) +
280 Ns * exp(s + gsl_sf_log_erfc(ys)) + Nr * exp(r + gsl_sf_log_erfc(yr))) -
281 eta * 2.0 / M_PI *
282 (Nu * exponentialIntegral(zu).imag() + Nv * exponentialIntegral(zv).imag() +
283 Ns * exponentialIntegral(zs).imag() + Nr * exponentialIntegral(zr).imag()));
284 }
285}
286
287void IkedaCarpenterPV::functionLocal(double *out, const double *xValues, const size_t nData) const {
288 const double I = getParameter("I");
289 const double alpha0 = getParameter("Alpha0");
290 const double alpha1 = getParameter("Alpha1");
291 const double beta0 = getParameter("Beta0");
292 const double kappa = getParameter("Kappa");
293 const double voigtsigmaSquared = getParameter("SigmaSquared");
294 const double voigtgamma = getParameter("Gamma");
295 const double X0 = getParameter("X0");
296
297 // cal pseudo voigt sigmaSq and gamma and eta
298 double gamma = 1.0; // dummy initialization
299 double eta = 0.5; // dummy initialization
300 convertVoigtToPseudo(voigtsigmaSquared, voigtgamma, gamma, eta);
301 double sigmaSquared = gamma * gamma / (8.0 * M_LN2); // pseudo voigt sigma^2
302
303 const double beta = 1 / beta0;
304
305 // equations taken from Fullprof manual
306
307 const double k = 0.05;
308
309 // Not entirely sure what to do if sigmaSquared ever negative
310 // for now just post a warning
311 double someConst = std::numeric_limits<double>::max() / 100.0;
312 if (sigmaSquared > 0)
313 someConst = 1 / sqrt(2.0 * sigmaSquared);
314 else if (sigmaSquared < 0) {
315 g_log.warning() << "sigmaSquared negative in functionLocal.\n";
316 }
317
318 // update wavelength vector
319 calWavelengthAtEachDataPoint(xValues, nData);
320
321 for (size_t i = 0; i < nData; i++) {
322 double diff = xValues[i] - X0;
323
324 double R = exp(-81.799 / (m_waveLength[i] * m_waveLength[i] * kappa));
325 double alpha = 1.0 / (alpha0 + m_waveLength[i] * alpha1);
326
327 double a_minus = alpha * (1 - k);
328 double a_plus = alpha * (1 + k);
329 double x = a_minus - beta;
330 double y = alpha - beta;
331 double z = a_plus - beta;
332
333 double Nu = 1 - R * a_minus / x;
334 double Nv = 1 - R * a_plus / z;
335 double Ns = -2 * (1 - R * alpha / y);
336 double Nr = 2 * R * alpha * alpha * beta * k * k / (x * y * z);
337
338 double u = a_minus * (a_minus * sigmaSquared - 2 * diff) / 2.0;
339 double v = a_plus * (a_plus * sigmaSquared - 2 * diff) / 2.0;
340 double s = alpha * (alpha * sigmaSquared - 2 * diff) / 2.0;
341 double r = beta * (beta * sigmaSquared - 2 * diff) / 2.0;
342
343 double yu = (a_minus * sigmaSquared - diff) * someConst;
344 double yv = (a_plus * sigmaSquared - diff) * someConst;
345 double ys = (alpha * sigmaSquared - diff) * someConst;
346 double yr = (beta * sigmaSquared - diff) * someConst;
347
348 std::complex<double> zs = std::complex<double>(-alpha * diff, 0.5 * alpha * gamma);
349 std::complex<double> zu = (1 - k) * zs;
350 std::complex<double> zv = (1 + k) * zs;
351 std::complex<double> zr = std::complex<double>(-beta * diff, 0.5 * beta * gamma);
352
353 double N = 0.25 * alpha * (1 - k * k) / (k * k);
354
355 out[i] = I * N *
356 ((1 - eta) * (Nu * exp(u + gsl_sf_log_erfc(yu)) + Nv * exp(v + gsl_sf_log_erfc(yv)) +
357 Ns * exp(s + gsl_sf_log_erfc(ys)) + Nr * exp(r + gsl_sf_log_erfc(yr))) -
358 eta * 2.0 / M_PI *
359 (Nu * exponentialIntegral(zu).imag() + Nv * exponentialIntegral(zv).imag() +
360 Ns * exponentialIntegral(zs).imag() + Nr * exponentialIntegral(zr).imag()));
361 }
362}
363
364void IkedaCarpenterPV::functionDerivLocal(API::Jacobian * /*jacobian*/, const double * /*xValues*/,
365 const size_t /*nData*/) {
366 throw Mantid::Kernel::Exception::NotImplementedError("functionDerivLocal is not implemented for IkedaCarpenterPV.");
367}
368
370 calNumericalDeriv(domain, jacobian);
371}
372
375 auto interval = getDomainInterval(1e-2);
376
378 API::IntegrationResult result = integrator.integrate(*this, interval.first, interval.second);
379
380 return result.result;
381}
382
383void IkedaCarpenterPV::setMatrixWorkspace(std::shared_ptr<const API::MatrixWorkspace> workspace, size_t wi,
384 double startX, double endX) {
385 if (workspace) {
386 // convert inital parameters that depend on x axis to correct units so
387 // inital guess is reasonable
388 auto tof = Mantid::Kernel::UnitFactory::Instance().create("TOF");
389 const auto centre = getParameter("X0");
390 const auto scaleFactor = centre / convertValue(centre, tof, workspace, wi);
391 if (scaleFactor != 0) {
392 if (isActive(parameterIndex("Alpha0")))
393 setParameter("Alpha0", getParameter("Alpha0") * scaleFactor);
394 if (isActive(parameterIndex("Alpha1")))
395 setParameter("Alpha1", getParameter("Alpha1") * scaleFactor);
396 if (isActive(parameterIndex("Beta0")))
397 setParameter("Beta0", getParameter("Beta0") * scaleFactor);
398 }
399 }
401}
402
403} // namespace Mantid::CurveFitting::Functions
gsl_vector * tmp
#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
Base class that represents the domain of a function.
size_t m_workspaceIndex
An index to a spectrum.
Definition: IFunctionMW.h:44
void setMatrixWorkspace(std::shared_ptr< const API::MatrixWorkspace > workspace, size_t wi, double startX, double endX) override
Set MatrixWorkspace.
Definition: IFunctionMW.cpp:22
std::shared_ptr< const API::MatrixWorkspace > getMatrixWorkspace() const
Get shared pointer to the workspace.
Definition: IFunctionMW.cpp:33
bool isActive(size_t i) const
Check if an active parameter i is actually active.
Definition: IFunction.cpp:160
virtual double getParameter(size_t i) const =0
Get i-th parameter.
double convertValue(double value, Kernel::Unit_sptr &outUnit, const std::shared_ptr< const MatrixWorkspace > &ws, size_t wsIndex) const
Convert a value from one unit (inUnit) to unit defined in workspace (ws)
Definition: IFunction.cpp:1267
void calNumericalDeriv(const FunctionDomain &domain, Jacobian &jacobian)
Calculate numerical derivatives.
Definition: IFunction.cpp:1031
virtual size_t parameterIndex(const std::string &name) const =0
Returns the index of parameter name.
virtual void addConstraint(std::unique_ptr< IConstraint > ic)
Add a constraint to function.
Definition: IFunction.cpp:376
virtual void declareParameter(const std::string &name, double initValue=0, const std::string &description="")=0
Declare a new parameter.
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.
virtual std::pair< double, double > getDomainInterval(double level=DEFAULT_SEARCH_LEVEL) const
Get the interval on which the peak has all its values above a certain level.
Represents the Jacobian in IFitFunction::functionDeriv.
Definition: Jacobian.h:22
IntegrationResult integrate(const IPeakFunction &peakFunction, double lowerLimit, double upperLimit) const
Integration of peak function on the interval [a, b].
Provide Ikeda-Carpenter-pseudo-Voigt peak shape function interface to IPeakFunction.
void functionLocal(double *out, const double *xValues, const size_t nData) const override
Function evaluation method to be implemented in the inherited classes.
void constFunction(double *out, const double *xValues, const int &nData) const
calculate the const function
void setHeight(const double h) override
Sets the parameters such that height == h.
void setCentre(const double c) override
Sets the parameters such that centre == c.
double intensity() const override
Returns the integral intensity of the peak.
double centre() const override
overwrite IPeakFunction base class methods
void init() override
overwrite IFunction base class method, which declare function parameters
void setMatrixWorkspace(std::shared_ptr< const API::MatrixWorkspace > workspace, size_t wi, double startX, double endX) override
Set MatrixWorkspace.
void setFwhm(const double w) override
Sets the parameters such that FWHM = w.
void lowerConstraint0(const std::string &paramName)
constrain all parameters to be non-negative
std::vector< double > m_waveLength
container for storing wavelength values for each data point
double height() const override
Returns the height of the function.
void convertVoigtToPseudo(const double &voigtSigmaSq, const double &voigtGamma, double &H, double &eta) const
convert voigt params to pseudo voigt params
void functionDerivLocal(API::Jacobian *out, const double *xValues, const size_t nData) override
Derivative evaluation method. Default is to calculate numerically.
void functionDeriv(const API::FunctionDomain &domain, API::Jacobian &jacobian) override
Derivatives of function with respect to active parameters.
double fwhm() const override
Returns the peak FWHM.
void calWavelengthAtEachDataPoint(const double *xValues, const size_t &nData) const
method for updating m_waveLength
Marks code as not implemented yet.
Definition: Exception.h:138
void debug(const std::string &msg)
Logs at debug level.
Definition: Logger.cpp:114
void warning(const std::string &msg)
Logs at warning level.
Definition: Logger.cpp:86
static T & Instance()
Return a reference to the Singleton instance, creating it if it does not already exist Creation is do...
Kernel::Logger g_log("ExperimentInfo")
static logger object
std::shared_ptr< const MatrixWorkspace > MatrixWorkspace_const_sptr
shared pointer to the matrix workspace base class (const version)
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.
std::shared_ptr< const IComponent > IComponent_const_sptr
Typdef of a shared pointer to a const IComponent.
Definition: IComponent.h:161
std::shared_ptr< const Instrument > Instrument_const_sptr
Shared pointer to an const instrument object.
std::shared_ptr< Unit > Unit_sptr
Shared pointer to the Unit base class.
Definition: Unit.h:229