13#include <gsl/gsl_errno.h>
14#include <gsl/gsl_sf_gamma.h>
15#include <gsl/gsl_spline.h>
17#include <boost/math/special_functions/hermite.hpp>
24using namespace CurveFitting;
31const char *HERMITE_PREFIX =
"C_";
32const char *KFSE_NAME =
"FSECoeff";
33const char *HERMITE_C_NAME =
"HermiteCoeffs";
34const int NFINE_Y = 1000;
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;
47 for (
size_t i = 1; i < endpoint; i++) {
50 area = stepsize / 2 * (yv[0] + 2 * area + yv[endpoint]);
62 Point(
double yvalue = 0.0,
double qvalue = 0.0) :
y(yvalue), q(qvalue) {}
69 bool operator()(Point
const &a, Point
const &b) {
return a.y < b.y; }
75 :
ComptonProfile(), m_hermite(), m_yfine(), m_qfine(), m_voigt(), m_voigtProfile(), m_userFixedFSE(false) {}
104 if (
name == HERMITE_C_NAME)
115 if (coeffs.empty()) {
116 throw std::invalid_argument(
"GramCharlierComptonProfile - Hermite polynomial string is empty!");
120 std::istringstream is(coeffs);
125 throw std::invalid_argument(
"NCSCountRate - Error reading int from hermite coefficient string: " + coeffs);
142 for (
size_t i = 0; i <
m_hermite.size(); ++i) {
144 std::ostringstream os;
145 os << HERMITE_PREFIX << 2 * i;
154 std::vector<size_t> indices;
156 for (
size_t i = 0; i <
m_hermite.size(); ++i) {
158 std::ostringstream os;
159 os << HERMITE_PREFIX << 2 * i;
166 indices.emplace_back(kIndex);
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);
192 std::vector<double> fse(NFINE_Y, 0.0);
193 std::vector<double> convolvedFSE(nData, 0.0);
198 for (
unsigned int i = 0; i <
m_hermite.size(); ++i) {
202 const unsigned int npoly = 2 * i;
206 std::transform(result.begin(), result.end(), convolvedFSE.begin(), result.begin(), std::plus<double>());
208 std::transform(result.begin(), result.end(), errors.begin(), result.begin(), std::divides<double>());
211 std::fill_n(profile.begin(), NFINE_Y, 0.0);
212 std::fill_n(result.begin(), nData, 0.0);
218 std::transform(convolvedFSE.begin(), convolvedFSE.end(), errors.begin(), convolvedFSE.begin(),
219 std::divides<double>());
220 cmatrix.
setColumn(start + col, convolvedFSE);
243 std::vector<double> sumH(NFINE_Y);
244 for (
unsigned int i = 0; i < nhermite; ++i) {
247 const unsigned int npoly = 2 * i;
268 const double ampNorm = amp / (std::sqrt(2.0 * M_PI) * wg);
270 std::ostringstream os;
271 os << HERMITE_PREFIX << npoly;
273 const double factorial = gsl_sf_fact(npoly / 2);
275 const double denom = ((std::pow(2.0,
static_cast<int>(npoly))) * factorial);
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;
289 assert(
static_cast<size_t>(NFINE_Y) == lhs.size());
294 const double ampNorm = amp / (std::sqrt(2.0 * M_PI) * wg);
301 throw std::runtime_error(
"The Y values or Q values have not been set");
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]);
318 const std::vector<double> &profile)
const {
319 const auto &modq =
modQ();
320 const auto &ei =
e0();
324 for (
size_t i = 0; i < nData; ++i) {
325 const std::vector<double> &voigt =
m_voigt[i];
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];
342 double startX,
double endX) {
361 const auto &yspace =
ySpace();
362 const auto &modq =
modQ();
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]);
373 std::sort(points.begin(), points.end(), InY());
375 std::vector<double> sortedy(ncoarseY), sortedq(ncoarseY);
376 for (
size_t i = 0; i < ncoarseY; ++i) {
377 const auto &p = points[i];
385 const double miny(sortedy.front()), maxy(sortedy.back());
386 const double step = (maxy - miny) /
static_cast<double>((NFINE_Y - 1));
389 gsl_interp_accel *acc = gsl_interp_accel_alloc();
390 gsl_spline *spline = gsl_spline_alloc(gsl_interp_linear, ncoarseY);
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;
395 m_qfine[i] = gsl_spline_eval(spline, xi, acc);
399 m_qfine.back() = gsl_spline_eval(spline, maxy, acc);
400 gsl_spline_free(spline);
401 gsl_interp_accel_free(acc);
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);
410 for (
size_t i = 0; i < ncoarseY; ++i) {
411 std::vector<double> &voigt =
m_voigt[i];
412 voigt.resize(NFINE_Y);
414 const double yi = yspace[i];
415 std::transform(minusYFine.begin(), minusYFine.end(), ym.begin(),
416 std::bind(std::plus<double>(), _1, yi));
double value
The value of the point.
#define DECLARE_FUNCTION(classname)
Macro for declaring a new type of function to be used with the FunctionFactory.
IPeaksWorkspace_sptr workspace
#define UNUSED_ARG(x)
Function arguments are sometimes unused in certain implmentations but are required for documentation ...
Attribute is a non-fitting parameter.
bool isActive(size_t i) const
Check if an active parameter i is actually active.
void declareAttribute(const std::string &name, const API::IFunction::Attribute &defaultValue)
Declare a single attribute.
virtual void declareAttributes()
Override to declare function attributes.
virtual void setAttribute(const std::string &name, const Attribute &)
Set a value to attribute attName.
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.
double mass() const
Access the mass.
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)
Simple data structure to store nominal detector values It avoids some functions taking a huge number ...