Mantid
Loading...
Searching...
No Matches
GramCharlierComptonProfile.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 <gsl/gsl_errno.h>
14#include <gsl/gsl_sf_gamma.h> // for factorial
15#include <gsl/gsl_spline.h>
16
17#include <boost/math/special_functions/hermite.hpp>
18
19#include <cmath>
20#include <sstream>
21
23
24using namespace CurveFitting;
25// Register into factory
26DECLARE_FUNCTION(GramCharlierComptonProfile)
27
28namespace {
30const char *WIDTH_PARAM = "Width";
31const char *HERMITE_PREFIX = "C_";
32const char *KFSE_NAME = "FSECoeff";
33const char *HERMITE_C_NAME = "HermiteCoeffs";
34const int NFINE_Y = 1000;
35
43double trapzf(const std::vector<double> &xv, const std::vector<double> &yv) {
44 const double stepsize = xv[1] - xv[0];
45 const size_t endpoint = xv.size() - 1;
46 double area(0.0);
47 for (size_t i = 1; i < endpoint; i++) {
48 area += yv[i];
49 }
50 area = stepsize / 2 * (yv[0] + 2 * area + yv[endpoint]); // final step
51 return area;
52}
53
54// Cannot put these inside massProfile definition as you can't use local types
55// as template
56// arguments on some compilers
57
58// massProfile needs to sort 1 vector but keep another in order. Use a struct of
59// combined points
60// and a custom comparator
61struct Point {
62 Point(double yvalue = 0.0, double qvalue = 0.0) : y(yvalue), q(qvalue) {}
63 double y;
64 double q;
65};
66
67struct InY {
69 bool operator()(Point const &a, Point const &b) { return a.y < b.y; }
70};
72} // namespace
73
75 : ComptonProfile(), m_hermite(), m_yfine(), m_qfine(), m_voigt(), m_voigtProfile(), m_userFixedFSE(false) {}
76
80std::string GramCharlierComptonProfile::name() const { return "GramCharlierComptonProfile"; }
81
83 // Base class ones
85 declareParameter(WIDTH_PARAM, 1.0, "Gaussian width parameter");
86 declareParameter(KFSE_NAME, 1.0, "FSE coefficient k");
87 // Other parameters depend on the Hermite attribute...
88}
89
91 // Base class ones
93 declareAttribute(HERMITE_C_NAME,
94 IFunction::Attribute("")); // space-separated string 1/0
95 // indicating which coefficients
96 // are active
97}
98
103void GramCharlierComptonProfile::setAttribute(const std::string &name, const Attribute &value) {
104 if (name == HERMITE_C_NAME)
105 setHermiteCoefficients(value.asString());
107}
108
115 if (coeffs.empty()) {
116 throw std::invalid_argument("GramCharlierComptonProfile - Hermite polynomial string is empty!");
117 }
118 m_hermite.clear();
119 m_hermite.reserve(3); // Maximum guess
120 std::istringstream is(coeffs);
121 while (!is.eof()) {
122 short value;
123 is >> value;
124 if (!is) {
125 throw std::invalid_argument("NCSCountRate - Error reading int from hermite coefficient string: " + coeffs);
126 }
127 m_hermite.emplace_back(value);
128 }
130}
131
138 // Gram-Chelier parameters are the even coefficents of the Hermite
139 // polynomials, i.e
140 // setting hermite coefficients to "1 0 1" uses coefficents C_0,C_4 and C_2 is
141 // skipped.
142 for (size_t i = 0; i < m_hermite.size(); ++i) {
143 if (m_hermite[i] > 0) {
144 std::ostringstream os;
145 os << HERMITE_PREFIX << 2 * i;
146 this->declareParameter(os.str(), 1.0, "Hermite polynomial coefficent");
147 }
148 }
149}
150
152 assert(!m_hermite.empty());
153
154 std::vector<size_t> indices;
155 indices.reserve(m_hermite.size() + 1);
156 for (size_t i = 0; i < m_hermite.size(); ++i) {
157 if (m_hermite[i] > 0) {
158 std::ostringstream os;
159 os << HERMITE_PREFIX << 2 * i; // refactor to have method that produces the name
160 indices.emplace_back(this->parameterIndex(os.str()));
161 }
162 }
163 // Include Kfse if it is not fixed
164 const size_t kIndex = this->parameterIndex(KFSE_NAME);
165 if (isActive(kIndex)) {
166 indices.emplace_back(kIndex);
167 }
168
169 return indices;
170}
171
182 const HistogramData::HistogramE &errors) const {
183 std::vector<double> profile(NFINE_Y, 0.0);
184 const size_t nData(ySpace().size());
185 std::vector<double> result(nData, 0.0);
186
187 // If the FSE term is fixed by user then it's contribution is convoluted with
188 // Voigt and summed with the first column
189 // otherwise it gets a column of it's own at the end.
190 // Either way it needs to be computed, so do this first
191
192 std::vector<double> fse(NFINE_Y, 0.0);
193 std::vector<double> convolvedFSE(nData, 0.0);
194 addFSETerm(fse);
195 convoluteVoigt(convolvedFSE.data(), nData, fse);
196
197 size_t col(0);
198 for (unsigned int i = 0; i < m_hermite.size(); ++i) {
199 if (m_hermite[i] == 0)
200 continue;
201
202 const unsigned int npoly = 2 * i;
203 addMassProfile(profile.data(), npoly);
204 convoluteVoigt(result.data(), nData, profile);
205 if (i == 0 && m_userFixedFSE) {
206 std::transform(result.begin(), result.end(), convolvedFSE.begin(), result.begin(), std::plus<double>());
207 }
208 std::transform(result.begin(), result.end(), errors.begin(), result.begin(), std::divides<double>());
209 cmatrix.setColumn(start + col, result);
210
211 std::fill_n(profile.begin(), NFINE_Y, 0.0);
212 std::fill_n(result.begin(), nData, 0.0);
213 ++col;
214 }
215
216 if (!m_userFixedFSE) // Extra column for He3
217 {
218 std::transform(convolvedFSE.begin(), convolvedFSE.end(), errors.begin(), convolvedFSE.begin(),
219 std::divides<double>());
220 cmatrix.setColumn(start + col, convolvedFSE);
221 ++col;
222 }
223 return col;
224}
225
234void GramCharlierComptonProfile::massProfile(double *result, const size_t nData) const {
235 UNUSED_ARG(nData);
236
237 using namespace Mantid::Kernel;
238
239 // Hermite expansion (only even terms) + FSE term
240 const size_t nhermite(m_hermite.size());
241
242 // Sum over polynomials for each y
243 std::vector<double> sumH(NFINE_Y);
244 for (unsigned int i = 0; i < nhermite; ++i) {
245 if (m_hermite[i] == 0)
246 continue;
247 const unsigned int npoly = 2 * i; // Only even ones
248 addMassProfile(sumH.data(), npoly);
249 }
250 addFSETerm(sumH);
251 convoluteVoigt(result, nData, sumH);
252}
253
263void GramCharlierComptonProfile::addMassProfile(double *result, const unsigned int npoly) const {
264
265 using namespace Mantid::Kernel;
266
267 const double amp(1.0), wg(getParameter(WIDTH_PARAM));
268 const double ampNorm = amp / (std::sqrt(2.0 * M_PI) * wg);
269
270 std::ostringstream os;
271 os << HERMITE_PREFIX << npoly;
272 const double hermiteCoeff = getParameter(os.str());
273 const double factorial = gsl_sf_fact(npoly / 2);
274 // Intel compiler doesn't overload pow for unsigned types
275 const double denom = ((std::pow(2.0, static_cast<int>(npoly))) * factorial);
276
277 for (int j = 0; j < NFINE_Y; ++j) {
278 const double y = m_yfine[j] / M_SQRT2 / wg;
279 const double hermiteI = boost::math::hermite(npoly, y);
280 result[j] += ampNorm * std::exp(-y * y) * hermiteI * hermiteCoeff / denom;
281 }
282}
283
288void GramCharlierComptonProfile::addFSETerm(std::vector<double> &lhs) const {
289 assert(static_cast<size_t>(NFINE_Y) == lhs.size());
290
291 using namespace Mantid::Kernel;
292
293 const double amp(1.0), wg(getParameter(WIDTH_PARAM));
294 const double ampNorm = amp / (std::sqrt(2.0 * M_PI) * wg);
295
296 double kfse(getParameter(KFSE_NAME));
297 if (m_userFixedFSE)
298 kfse *= getParameter("C_0");
299
300 if (m_yfine.empty() || m_qfine.empty()) {
301 throw std::runtime_error("The Y values or Q values have not been set");
302 }
303
304 for (int j = 0; j < NFINE_Y; ++j) {
305 const double y = m_yfine[j] / M_SQRT2 / wg;
306 const double he3 = boost::math::hermite(3, y);
307 lhs[j] += ampNorm * std::exp(-y * y) * he3 * (kfse / m_qfine[j]);
308 }
309}
310
317void GramCharlierComptonProfile::convoluteVoigt(double *result, const size_t nData,
318 const std::vector<double> &profile) const {
319 const auto &modq = modQ();
320 const auto &ei = e0();
321
322 // Now convolute with the Voigt function (pre-calculated in setWorkspace as
323 // its expensive)
324 for (size_t i = 0; i < nData; ++i) {
325 const std::vector<double> &voigt = m_voigt[i];
326 // Multiply voigt with polynomial sum and put result in voigt to save using
327 // another vector
328 std::transform(voigt.begin(), voigt.end(), profile.begin(), m_voigtProfile.begin(), std::multiplies<double>());
329 const double prefactor = std::pow(ei[i], 0.1) * mass() / modq[i];
330 result[i] = prefactor * trapzf(m_yfine, m_voigtProfile);
331 }
332}
333
341void GramCharlierComptonProfile::setMatrixWorkspace(std::shared_ptr<const API::MatrixWorkspace> workspace, size_t wi,
342 double startX, double endX) {
344 endX); // Do base-class calculation first
345}
346
351void GramCharlierComptonProfile::cacheYSpaceValues(const HistogramData::Points &tseconds,
352 const Algorithms::DetectorParams &detpar) {
354 detpar); // base-class calculations
355
356 // Is FSE fixed at the moment?
357 // The ComptonScatteringCountRate fixes it but we still need to know if the
358 // user wanted it fixed
359 m_userFixedFSE = !this->isActive(this->parameterIndex(KFSE_NAME));
360
361 const auto &yspace = ySpace();
362 const auto &modq = modQ();
363
364 // massProfile is calculated over a large range of Y, constructed by
365 // interpolation
366 // This is done over an interpolated range between ymin & ymax and y and hence
367 // q must be sorted
368 const size_t ncoarseY(yspace.size());
369 std::vector<Point> points(ncoarseY);
370 for (size_t i = 0; i < ncoarseY; ++i) {
371 points[i] = Point(yspace[i], modq[i]);
372 }
373 std::sort(points.begin(), points.end(), InY());
374 // Separate back into vectors as GSL requires them separate
375 std::vector<double> sortedy(ncoarseY), sortedq(ncoarseY);
376 for (size_t i = 0; i < ncoarseY; ++i) {
377 const auto &p = points[i];
378 sortedy[i] = p.y;
379 sortedq[i] = p.q;
380 }
381
382 // Generate a more-finely grained y axis and interpolate Q values
383 m_yfine.resize(NFINE_Y);
384 m_qfine.resize(NFINE_Y);
385 const double miny(sortedy.front()), maxy(sortedy.back());
386 const double step = (maxy - miny) / static_cast<double>((NFINE_Y - 1));
387
388 // Set up GSL interpolater
389 gsl_interp_accel *acc = gsl_interp_accel_alloc();
390 gsl_spline *spline = gsl_spline_alloc(gsl_interp_linear, ncoarseY); // Actually a linear interpolater
391 gsl_spline_init(spline, sortedy.data(), sortedq.data(), ncoarseY);
392 for (int i = 0; i < NFINE_Y - 1; ++i) {
393 const double xi = miny + step * i;
394 m_yfine[i] = xi;
395 m_qfine[i] = gsl_spline_eval(spline, xi, acc);
396 }
397 // Final value to ensure it ends at maxy
398 m_yfine.back() = maxy;
399 m_qfine.back() = gsl_spline_eval(spline, maxy, acc);
400 gsl_spline_free(spline);
401 gsl_interp_accel_free(acc);
402
403 // Cache voigt function over yfine
404 std::vector<double> minusYFine(NFINE_Y);
405 using std::placeholders::_1;
406 std::transform(m_yfine.begin(), m_yfine.end(), minusYFine.begin(), std::bind(std::multiplies<double>(), _1, -1.0));
407 std::vector<double> ym(NFINE_Y); // Holds result of (y[i] - yfine) for each original y
408 m_voigt.resize(ncoarseY);
409
410 for (size_t i = 0; i < ncoarseY; ++i) {
411 std::vector<double> &voigt = m_voigt[i];
412 voigt.resize(NFINE_Y);
413
414 const double yi = yspace[i];
415 std::transform(minusYFine.begin(), minusYFine.end(), ym.begin(),
416 std::bind(std::plus<double>(), _1, yi)); // yfine is actually -yfine
417 m_resolutionFunction->voigtApprox(voigt, ym, 0, 1.0);
418 }
419
420 m_voigtProfile.resize(NFINE_Y); // Value holder for later to avoid repeated
421 // memory allocations when creating a new
422 // vector
423}
424
425} // 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.
IPeaksWorkspace_sptr workspace
Definition: IndexPeaks.cpp:114
#define UNUSED_ARG(x)
Function arguments are sometimes unused in certain implmentations but are required for documentation ...
Definition: System.h:64
Attribute is a non-fitting parameter.
Definition: IFunction.h:282
bool isActive(size_t i) const
Check if an active parameter i is actually active.
Definition: IFunction.cpp:160
void declareAttribute(const std::string &name, const API::IFunction::Attribute &defaultValue)
Declare a single attribute.
Definition: IFunction.cpp:1418
virtual void declareAttributes()
Override to declare function attributes.
Definition: IFunction.h:676
virtual void setAttribute(const std::string &name, const Attribute &)
Set a value to attribute attName.
Definition: IFunction.cpp:1409
size_t parameterIndex(const std::string &name) const override
Returns the index of parameter name.
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.
This class serves as a base-class for ComptonProfile type functions.
const std::vector< double > & modQ() const
Access Q values cache.
void declareParameters() override
Declare parameters that will never participate in the fit.
std::shared_ptr< VesuvioResolution > m_resolutionFunction
Vesuvio resolution function.
void setMatrixWorkspace(std::shared_ptr< const API::MatrixWorkspace > workspace, size_t wsIndex, double startX, double endX) override
Cache a copy of the workspace pointer and pull out the parameters.
const std::vector< double > & ySpace() const
Access y-values cache.
const std::vector< double > & e0() const
Access e0 values.
void cacheYSpaceValues(const HistogramData::Points &tseconds, const Algorithms::DetectorParams &detpar, const ResolutionParams &respar)
Pre-calculate the Y-space values with specified resolution parameters.
void declareParameters() override
Declare the function parameters.
std::string name() const override
A string identifier for this function.
void addMassProfile(double *result, const unsigned int npoly) const
Compute the contribution to mass profile nth Hermite polynomial coefficient.
std::vector< double > m_qfine
Interpolated Q values over a finer Y range.
void declareAttributes() override
Declare non-parameter attributes.
std::vector< double > m_yfine
Y values over a finer range.
void massProfile(double *result, const size_t nData) const override
Compute the sum for all Hermite polynomial coefficents.
void declareGramCharlierParameters()
Declare the Gram-Charlier (Hermite) coefficients.
std::vector< double > m_voigtProfile
Holds the result Voigt multiplied by the profile function for the extended Y space range.
void setHermiteCoefficients(const std::string &coeffs)
Parse the active hermite polynomial coefficents.
bool m_userFixedFSE
Flag to hold whether the FSE parameter is fixed by the user.
void setAttribute(const std::string &name, const Attribute &value) override
Set an attribute value (and possibly cache its value)
GramCharlierComptonProfile()
Default constructor required by factory.
void addFSETerm(std::vector< double > &lhs) const
Add FSE term based on current parameter setting.
std::vector< short > m_hermite
The active hermite coefficents.
std::vector< std::vector< double > > m_voigt
Holds the value of the Voigt function for each coarse y-space point as this is an expensive calculati...
void convoluteVoigt(double *result, const size_t nData, const std::vector< double > &profile) const
Convolute with resolution.
void setMatrixWorkspace(std::shared_ptr< const API::MatrixWorkspace > workspace, size_t wi, double startX, double endX) override
Called by the framework when a workspace is set.
size_t fillConstraintMatrix(Kernel::DblMatrix &cmatrix, const size_t start, const HistogramData::HistogramE &errors) const override
Fill in the columns of the matrix for this mass.
std::vector< size_t > intensityParameterIndices() const override
Returns the indices of the intensity parameters.
void cacheYSpaceValues(const HistogramData::Points &tseconds, const Algorithms::DetectorParams &detpar) override
Pre-calculate the Y-space values.
void setColumn(const size_t nCol, const std::vector< T > &newCol)
Definition: Matrix.cpp:675
Simple data structure to store nominal detector values It avoids some functions taking a huge number ...