12#include "MantidHistogramData/HistogramX.h"
13#include "MantidHistogramData/HistogramY.h"
19#include <gsl/gsl_sf_erf.h>
23using Mantid::HistogramData::HistogramY;
31const double PEAKRANGECONSTANT = 5.0;
33const string CHEBYSHEV_BACKGROUND(
"Chebyshev");
34const string POLYNOMIAL_BACKGROUND(
"Polynomial");
35const string FULLPROF_POLYNOMIAL_BACKGROUND(
"FullprofPolynomial");
44LeBailFunction::LeBailFunction(
const std::string &peaktype) {
59 errss <<
"Input peak type " << peaktype <<
" is not a recoganizable Mantid function.";
60 throw runtime_error(errss.str());
65 errss <<
"Input peak type " << peaktype <<
" is not a IPowderDiffPeakFunction.";
66 throw runtime_error(errss.str());
86LeBailFunction::~LeBailFunction() =
default;
104HistogramY LeBailFunction::function(
const Mantid::HistogramData::HistogramX &xvalues,
bool calpeaks,
105 bool calbkgd)
const {
108 std::vector<double> out(xvalues.size(), 0);
109 const auto &xvals = xvalues.rawData();
113 for (
size_t ipk = 0; ipk <
m_numPeaks; ++ipk) {
115 vector<double> temp(xvalues.size(), 0);
117 peak->function(temp, xvals);
118 transform(out.begin(), out.end(), temp.begin(), out.begin(), ::plus<double>());
125 throw runtime_error(
"Must define background first!");
132 size_t numpts = out.size();
133 for (
size_t i = 0; i < numpts; ++i)
137 return HistogramY(out);
142HistogramY LeBailFunction::calPeak(
size_t ipk,
const std::vector<double> &xvalues,
size_t ySize)
const {
146 errss <<
"Try to calculate peak indexed " << ipk <<
". But number of peaks = " <<
m_numPeaks;
148 throw runtime_error(errss.str());
151 std::vector<double> out(ySize, 0);
153 peak->function(out, xvalues);
154 return HistogramY(out);
161bool LeBailFunction::hasProfileParameter(
const std::string ¶mname) {
170 string matchparname = *fiter;
171 if (matchparname != paramname)
183bool LeBailFunction::isParameterValid(
double maxfwhm)
const {
190 bool arevalid =
true;
193 bool isvalid = peak->isPhysical();
194 if (isvalid && maxfwhm >= 0)
195 isvalid = peak->fwhm() < maxfwhm;
200 peak->getMillerIndex(h, k, l);
201 g_log.
information() <<
"Peak [" << h <<
", " << k <<
", " << l <<
"] @ TOF = " << peak->centre()
202 <<
" has unphysical parameters or unreasonable large FWHM"
214void LeBailFunction::calculatePeakParameterValues()
const {
217 peak->calculateParameters(
false);
229void LeBailFunction::setPeakCentreTolerance(
double peakpostol,
double tofmin,
double tofmax) {
239void LeBailFunction::addPeaks(std::vector<std::vector<int>> peakhkls) {
242 throw runtime_error(
"Client must set up profile parameter vlaues by calling "
243 "setProfileParameterValues() first! ");
246 for (
size_t ipk = 0; ipk < peakhkls.size(); ++ipk) {
247 vector<int> hkl = peakhkls[ipk];
250 if (hkl.size() != 3) {
252 errss <<
"Error of " << ipk <<
"-th input Miller Index. It has " << peakhkls[ipk].size()
253 <<
" items, but not required 3 items.";
255 throw runtime_error(errss.str());
265 throw runtime_error(
"Unable to generate peak.");
268 double tofh = newpeak->centre();
270 g_log.
information() <<
"Peak " << h <<
", " << k <<
", " << l <<
" 's centre is at TOF = " << tofh
275 double dsp = newpeak->getPeakParameter(
"d_h");
300 peak->setMillerIndex(h, k, l);
303 peak->setParameter(parname, parvalue);
324bool LeBailFunction::calculatePeaksIntensities(
const vector<double> &vecX,
const vector<double> &vecY,
325 vector<double> &vec_summedpeaks) {
327 std::fill(vec_summedpeaks.begin(), vec_summedpeaks.end(), 0.0);
330 vector<vector<pair<double, IPowderDiffPeakFunction_sptr>>> peakgroupvec;
331 vector<IPowderDiffPeakFunction_sptr> outboundpeakvec;
332 double xmin = vecX.front();
333 double xmax = vecX.back();
334 groupPeaks(peakgroupvec, outboundpeakvec, xmin, xmax);
337 bool allpeakheightsphysical =
true;
338 for (
size_t ig = 0; ig < peakgroupvec.size(); ++ig) {
339 g_log.
debug() <<
"[Fx351] Calculate peaks heights for (peak) group " << ig
340 <<
" : number of peaks = " << peakgroupvec[ig].size() <<
"\n";
344 if (!peakheightsphysical)
345 allpeakheightsphysical =
false;
349 for (
const auto &peak : outboundpeakvec) {
353 return allpeakheightsphysical;
367bool LeBailFunction::calculateGroupPeakIntensities(vector<pair<double, IPowderDiffPeakFunction_sptr>> peakgroup,
368 const vector<double> &vecX,
const vector<double> &vecY,
369 vector<double> &vec_summedpeaks) {
371 if (peakgroup.empty()) {
372 throw runtime_error(
"Programming error such that input peak group cannot be empty!");
374 g_log.
debug() <<
"[Fx155] Peaks group size = " << peakgroup.size() <<
"\n";
376 if (peakgroup.size() > 1)
377 sort(peakgroup.begin(), peakgroup.end());
380 if (vec_summedpeaks.size() != vecY.size()) {
382 errss <<
"Input vector 'allpeaksvalues' has wrong size = " << vec_summedpeaks.size()
383 <<
" != data workspace Y's size = " << vecY.size();
385 throw runtime_error(errss.str());
390 double leftbound = leftpeak->centre() - PEAKRANGECONSTANT * leftpeak->fwhm();
391 if (leftbound < vecX.front()) {
394 leftpeak->getMillerIndex(h, k, l);
395 msg <<
"Peak group (containing " << peakgroup.size() <<
" peaks) has its left boundary (TOF = " << leftbound
396 <<
") out side of input data workspace's left boundary (" << vecX.front()
397 <<
"). Accuracy of its peak intensity might be affected. "
398 <<
"Group's left boundary is determined by its leftmost peak (" << h <<
", " << k <<
", " << l
399 <<
") at TOF = " << leftpeak->centre() <<
" with FWHM = " << leftpeak->fwhm() <<
". ";
403 leftbound = vecX[0] + 0.1;
406 double rightbound = rightpeak->centre() + PEAKRANGECONSTANT * rightpeak->fwhm();
407 if (rightbound > vecX.back()) {
409 msg <<
"Peak group's right boundary " << rightbound <<
" is out side of "
410 <<
"input data workspace's right bound (" << vecX.back()
411 <<
")! Accuracy of its peak intensity might be affected. ";
415 rightbound = vecX.back() - 0.1;
419 vector<double>::const_iterator cviter;
421 cviter = lower_bound(vecX.begin(), vecX.end(), leftbound);
422 size_t ileft =
static_cast<size_t>(cviter - vecX.begin());
426 cviter = lower_bound(vecX.begin(), vecX.end(), rightbound);
427 size_t iright =
static_cast<size_t>(cviter - vecX.begin());
428 if (iright <= vecX.size() - 1)
431 size_t ndata = iright - ileft;
432 if (ileft >= iright) {
434 errss <<
"[Calcualte Peak Intensity] Group range is unphysical. iLeft = " << ileft <<
", iRight = " << iright
435 <<
"; Number of peaks = " << peakgroup.size() <<
"; Left boundary = " << leftbound
436 <<
", Right boundary = " << rightbound <<
"; Left peak FWHM = " << leftpeak->fwhm()
437 <<
", Right peak FWHM = " << rightpeak->fwhm();
438 for (
size_t ipk = 0; ipk < peakgroup.size(); ++ipk) {
440 errss <<
"Peak " << ipk <<
": d_h = " << peakgroup[ipk].first <<
", TOF_h = " << thispeak->centre()
441 <<
", FWHM = " << thispeak->fwhm() <<
"\n";
442 vector<string> peakparamnames = thispeak->getParameterNames();
443 for (
auto &peakparamname : peakparamnames) {
444 errss <<
"\t" << peakparamname <<
" = " << thispeak->getParameter(peakparamname) <<
"\n";
449 throw runtime_error(errss.str());
453 vector<double> datax(vecX.begin() + ileft, vecX.begin() + iright);
454 vector<double> datay(vecY.begin() + ileft, vecY.begin() + iright);
455 if (datax.size() != ndata) {
457 errmsg <<
"Impossible: Partial peak data size = " << datax.size() <<
" != ndata = " << ndata;
459 throw runtime_error(errmsg.str());
461 g_log.
debug() <<
"[DBx356] Number of data points = " << ndata <<
" index from " << ileft <<
" to " << iright
462 <<
"; Size(datax, datay) = " << datax.size() <<
"\n";
465 vector<double> sumYs(ndata, 0.0);
466 size_t numPeaks(peakgroup.size());
467 vector<vector<double>> peakvalues(numPeaks);
470 bool datavalueinvalid =
false;
471 for (
size_t ipk = 0; ipk < numPeaks; ++ipk) {
475 peak->setHeight(1.0);
476 vector<double> localpeakvalue(ndata, 0.0);
477 peak->function(localpeakvalue, datax);
481 vector<double>::const_iterator localpeakvalue_end = localpeakvalue.end();
482 for (
auto it = localpeakvalue.begin(); it != localpeakvalue_end; ++it) {
483 if ((*it != 0.) && (*it < NEG_DBL_MAX || *it > DBL_MAX)) {
489 if (numbadpts == 0) {
491 for (
size_t i = 0; i < ndata; ++i) {
493 sumYs[i] += localpeakvalue[i];
498 peak->getMillerIndex(h, k, l);
500 warnss <<
"Peak (" << h <<
", " << k <<
", " << l <<
") @ TOF = " << peak->centre() <<
" has " << numbadpts
502 <<
"whose values exceed limit (i.e., not physical). ";
504 datavalueinvalid =
true;
506 peakvalues[ipk].assign(localpeakvalue.begin(), localpeakvalue.end());
510 bool peakheightsphysical = !datavalueinvalid;
511 if (peakheightsphysical) {
512 for (
size_t ipk = 0; ipk < peakgroup.size(); ++ipk) {
514 double intensity = 0.0;
516 for (
size_t i = 0; i < ndata; ++i) {
518 if (sumYs[i] > 1.0E-5) {
520 double peaktogroupratio = peakvalues[ipk][i] / sumYs[i];
521 temp = datay[i] * peaktogroupratio;
528 deltax = datax[1] - datax[0];
530 deltax = datax[i] - datax[i - 1];
532 intensity += temp * deltax;
535 if (intensity != intensity) {
538 peakheightsphysical =
false;
541 peak->getMillerIndex(h, k, l);
542 g_log.
warning() <<
"Peak (" << h <<
", " << k <<
", " << l <<
") has unphysical intensity = NaN!\n";
544 }
else if (intensity <= -DBL_MAX || intensity >= DBL_MAX) {
547 peakheightsphysical =
false;
550 peak->getMillerIndex(h, k, l);
551 g_log.
warning() <<
"Peak (" << h <<
", " << k <<
", " << l <<
") has unphysical intensity = Infty!\n";
552 }
else if (intensity < 0.0) {
554 g_log.
debug() <<
"[Fx134] Set peak @ " << peak->centre() <<
"'s intensity to 0.0 instead of " << intensity
558 g_log.
debug() <<
"[Fx407] Peak @ " << peak->centre() <<
": Set Intensity = " << intensity <<
"\n";
559 peak->setHeight(intensity);
562 for (
size_t i = ileft; i < iright; ++i) {
563 vec_summedpeaks[i] += (intensity * peakvalues[ipk][i - ileft]);
569 return peakheightsphysical;
582 double peakheight,
bool setpeakheight) {
587 throw runtime_error(
"Requiring update flag: peak value changed and etc.");
641void LeBailFunction::setProfileParameterValues(map<std::string, double> parammap) {
642 const double MINDIFF = 1.0E-10;
644 map<std::string, double>::iterator inpiter, curiter;
647 for (
size_t i = 0; i < numpars; ++i) {
651 inpiter = parammap.find(parname);
652 if (inpiter != parammap.end()) {
658 errmsg <<
"Parameter " << parname <<
" is in parameter name list, but not in profile "
659 <<
"parameter map. It violates the programming logic.";
661 throw runtime_error(errmsg.str());
665 double curvalue = curiter->second;
666 double newvalue = inpiter->second;
667 bool localnewvalue =
false;
668 if (
fabs(curvalue - newvalue) > MINDIFF) {
669 curiter->second = newvalue;
671 localnewvalue =
true;
679 for (
size_t ipk = 0; ipk <
m_numPeaks; ++ipk) {
681 peak->setParameter(i, newvalue);
685 g_log.
debug() <<
"Parameter " << parname <<
" is not a profile parameter. Length of string = " << parname.size()
703void LeBailFunction::groupPeaks(vector<vector<pair<double, IPowderDiffPeakFunction_sptr>>> &peakgroupvec,
704 vector<IPowderDiffPeakFunction_sptr> &outboundpeakvec,
double xmin,
double xmax) {
709 std::stringstream errmsg;
710 errmsg <<
"Group peaks: No peak is found in the peak vector. ";
712 throw std::runtime_error(errmsg.str());
716 peakgroupvec.clear();
717 outboundpeakvec.clear();
718 vector<pair<double, IPowderDiffPeakFunction_sptr>> peakgroup;
722 bool outbound =
true;
726 if (peak->centre() <= xmin) {
728 outboundpeakvec.emplace_back(peak);
741 if (thispeak->centre() < xmax) {
753 double thispeak_rightbound = thispeak->centre() + PEAKRANGECONSTANT * thispeak->fwhm();
754 double rightpeak_leftbound = rightpeak->centre() - PEAKRANGECONSTANT * rightpeak->fwhm();
756 if (thispeak_rightbound < rightpeak_leftbound) {
759 peakgroupvec.emplace_back(std::move(peakgroup));
768 peakgroupvec.emplace_back(peakgroup);
776 g_log.
information() <<
"[Fx301] Group peak: peak @ " << thispeak->centre() <<
" causes grouping "
777 <<
"peak over at maximum TOF = " << xmax <<
".\n";
779 if (!peakgroup.empty()) {
780 peakgroupvec.emplace_back(peakgroup);
791 g_log.
debug() <<
"[Calculate Peak Intensity]: Number of Peak Groups = " << peakgroupvec.size() <<
"\n";
805void LeBailFunction::addBackgroundFunction(
const string &backgroundtype,
const unsigned int &order,
806 const std::vector<std::string> &vecparnames,
807 const std::vector<double> &vecparvalues,
double startx,
double endx) {
809 if (backgroundtype != POLYNOMIAL_BACKGROUND && backgroundtype != CHEBYSHEV_BACKGROUND &&
810 backgroundtype != FULLPROF_POLYNOMIAL_BACKGROUND) {
812 warnss <<
"Cliet specified background type " << backgroundtype <<
" may not be supported properly.";
815 if (vecparnames.size() != vecparvalues.size())
816 throw runtime_error(
"Input parameter names and parameter values are not matched. ");
818 g_log.
information() <<
"Add background: type = " << backgroundtype <<
", order = " << order
819 <<
", number of parameters/attributes = " << vecparnames.size() <<
"\n";
823 m_background = std::dynamic_pointer_cast<Functions::BackgroundFunction>(background);
826 m_background->setAttributeValue(
"n",
static_cast<int>(order));
830 size_t numpars = vecparnames.size();
831 for (
size_t i = 0; i < numpars; ++i) {
832 const string &parname = vecparnames[i];
833 if (parname !=
"Bkpos")
835 else if (backgroundtype == FULLPROF_POLYNOMIAL_BACKGROUND)
836 m_background->setAttributeValue(
"Bkpos", vecparvalues[i]);
838 throw runtime_error(
"Bkpos should not be in the parameter list. ");
841 if (backgroundtype == CHEBYSHEV_BACKGROUND) {
855void LeBailFunction::setFitProfileParameter(
const string ¶mname,
double minvalue,
double maxvalue) {
857 for (
size_t ipk = 1; ipk <
m_numPeaks; ++ipk) {
858 stringstream ss1, ss2;
859 ss1 <<
"f" << (ipk - 1) <<
"." << paramname;
860 ss2 <<
"f" << ipk <<
"." << paramname;
861 string tiepart1 = ss1.str();
862 string tiepart2 = ss2.str();
864 g_log.
debug() <<
"LeBailFunction::Fit(Tie) / " << tiepart1 <<
" / " << tiepart2 <<
" /\n";
868 std::stringstream parss;
869 parss <<
"f0." << paramname;
870 string parnamef0 = parss.str();
871 auto bc = std::make_unique<Constraints::BoundaryConstraint>(
m_compsiteFunction.get(), parnamef0, minvalue, maxvalue);
880void LeBailFunction::fixPeakParameter(
const string ¶mname,
double paramvalue) {
881 for (
size_t ipk = 0; ipk <
m_numPeaks; ++ipk) {
882 stringstream ss1, ss2;
883 ss1 <<
"f" << ipk <<
"." << paramname;
885 string tiepart1 = ss1.str();
886 string tievalue = ss2.str();
889 g_log.
debug() <<
"Set up tie | " << tiepart1 <<
" <---> " << tievalue <<
" | \n";
905void LeBailFunction::fixBackgroundParameters() {
908 for (
size_t iparam = 0; iparam < numbkgdparams; ++iparam)
915void LeBailFunction::setFixPeakHeights() {
916 for (
size_t ipk = 0; ipk <
m_numPeaks; ++ipk) {
927void LeBailFunction::setPeakHeights(
const std::vector<double> &inheights) {
929 throw runtime_error(
"It is not implemented properly.");
951 errmsg <<
"Try to access peak " << peakindex <<
" out of range [0, " <<
m_numPeaks <<
").";
953 throw runtime_error(errmsg.str());
964double LeBailFunction::getPeakParameter(std::vector<int> hkl,
const std::string &parname)
const {
970 errss <<
"Peak with Miller index (" << hkl[0] <<
", " << hkl[1] <<
"," << hkl[2]
971 <<
") does not exist in Le Bail function.";
973 throw runtime_error(errss.str());
986double LeBailFunction::getPeakParameter(
size_t index,
const std::string &parname)
const {
989 errss <<
"getPeakParameter() tries to reach a peak with index " <<
index <<
", which is out of range " <<
m_numPeaks
992 throw std::runtime_error(errss.str());
1007 const std::string &parname)
const {
1017 string matchparname = *vsiter;
1018 if (parname != matchparname)
1026 parvalue = peak->getParameter(parname);
1029 parvalue = peak->getPeakParameter(parname);
1038double LeBailFunction::getPeakMaximumValue(std::vector<int> hkl,
const std::vector<double> &xvalues,
size_t &ix) {
1044 errss <<
"Peak with Miller index (" << hkl[0] <<
", " << hkl[1] <<
"," << hkl[2]
1045 <<
") does not exist in Le Bail function.";
1047 throw runtime_error(errss.str());
1052 double maxvalue = peak->getMaximumValue(xvalues, ix);
double value
The value of the point.
std::map< DeltaEMode::Type, std::string > index
#define UNUSED_ARG(x)
Function arguments are sometimes unused in certain implmentations but are required for documentation ...
A composite function is a function containing other functions.
Implements FunctionDomain1D with its own storage in form of a std::vector.
A class to store values calculated by a function.
bool m_isInputValue
Has first value set up.
std::vector< std::string > m_peakParameterNameVec
Name of peak parameter names (be same as the order in IPowderDiffPeakFunction)
void calculatePeakParameterValues() const
Calculate all peaks' parameter value.
void groupPeaks(std::vector< std::vector< std::pair< double, API::IPowderDiffPeakFunction_sptr > > > &peakgroupvec, std::vector< API::IPowderDiffPeakFunction_sptr > &outboundpeakvec, double xmin, double xmax)
Group close peaks together.
size_t m_numPeaks
Number of peaks.
std::vector< std::pair< double, API::IPowderDiffPeakFunction_sptr > > m_dspPeakVec
Vector of pair <peak position in d-space, Peak> sortable.
double m_maxTOFPeakCentre
double m_minTOFPeakCentre
bool calculateGroupPeakIntensities(std::vector< std::pair< double, API::IPowderDiffPeakFunction_sptr > > peakgroup, const std::vector< double > &vecX, const std::vector< double > &vecY, std::vector< double > &vec_summedpeaks)
Calculate the peaks intensities in same group.
Functions::BackgroundFunction_sptr m_background
Background function.
bool m_hasNewPeakValue
Has new peak values.
API::IPowderDiffPeakFunction_sptr generatePeak(int h, int k, int l)
Generate a peak with parameter set by.
std::vector< API::IPowderDiffPeakFunction_sptr > m_vecPeaks
Vector of all peaks.
double getPeakParameterValue(const API::IPowderDiffPeakFunction_sptr &peak, const std::string &parname) const
Retrieve peak's parameter. may be native or calculated.
API::CompositeFunction_sptr m_compsiteFunction
Composite functions for all peaks and background.
std::vector< std::string > m_orderedProfileParameterNames
Ordered profile parameter names for search.
std::map< std::vector< int >, API::IPowderDiffPeakFunction_sptr > m_mapHKLPeak
Vector of all peak's Miller indexes.
std::map< std::string, double > m_functionParameters
Parameters.
std::string m_peakType
Peak type.
The Logger class is in charge of the publishing messages from the framework through various channels.
void debug(const std::string &msg)
Logs at debug level.
void error(const std::string &msg)
Logs at error level.
void warning(const std::string &msg)
Logs at warning level.
void information(const std::string &msg)
Logs at information level.
static T & Instance()
Return a reference to the Singleton instance, creating it if it does not already exist Creation is do...
std::shared_ptr< IPowderDiffPeakFunction > IPowderDiffPeakFunction_sptr
Kernel::Logger g_log("ExperimentInfo")
static logger object
std::shared_ptr< IFunction > IFunction_sptr
shared pointer to the function base class
std::shared_ptr< CompositeFunction > CompositeFunction_sptr
shared pointer to the composite function base class
static constexpr double h
Planck constant in J*s.