17#include "MantidIndexing/GlobalSpectrumIndex.h"
18#include "MantidIndexing/IndexInfo.h"
25#include <boost/algorithm/string.hpp>
26#include <boost/math/special_functions/round.hpp>
49 :
API::
ParallelAlgorithm(), m_peakParameterNames(), m_bkgdParameterNames(), m_bkgdOrder(0), m_outPeakTableWS(),
50 m_dataWS(), m_inputPeakFWHM(0), m_highBackground(false), m_rawPeaksTable(false), m_numTableParams(0),
51 m_centreIndex(1) , m_peakFuncType(""), m_backgroundType(""), m_vecPeakCentre(),
52 m_vecFitWindows(), m_backgroundFunction(), m_peakFunction(), m_minGuessedPeakWidth(0), m_maxGuessedPeakWidth(0),
53 m_stepGuessedPeakWidth(0), m_usePeakPositionTolerance(false), m_peakPositionTolerance(0.0), m_fitFunctions(),
54 m_peakLeftIndexes(), m_peakRightIndexes(), m_minimizer("Levenberg-MarquardtMD"),
m_costFunction(),
55 m_minHeight(0.0), m_leastMaxObsY(0.), m_useObsCentre(false) {}
62 "Name of the workspace to search");
64 auto mustBeNonNegative = std::make_shared<BoundedValidator<int>>();
65 mustBeNonNegative->setLower(0);
67 "If set, only this spectrum will be searched for peaks "
68 "(otherwise all are)");
70 auto min = std::make_shared<BoundedValidator<int>>();
73 declareProperty(
"FWHM", 7, min,
"Estimated number of points covered by the fwhm of a peak (default 7)");
77 "A measure of the strictness desired in "
78 "meeting the condition on peak "
80 "Mariscotti recommends 2 (default 4)");
83 "Optional: enter a comma-separated list of the expected "
84 "X-position of the centre of the peaks. Only peaks near "
85 "these positions will be fitted.");
88 "Optional: enter a comma-separated list of the expected "
89 "X-position of windows to fit. The number of values must be "
90 "exactly double the number of specified peaks.");
93 declareProperty(
"PeakFunction",
"Gaussian", std::make_shared<StringListValidator>(peakNames));
95 std::vector<std::string> bkgdtypes{
"Flat",
"Linear",
"Quadratic"};
96 declareProperty(
"BackgroundType",
"Linear", std::make_shared<StringListValidator>(bkgdtypes),
"Type of Background.");
98 declareProperty(
"HighBackground",
true,
"Relatively weak peak in high background");
100 auto mustBePositive = std::make_shared<BoundedValidator<int>>();
101 mustBePositive->setLower(1);
103 "Minimum guessed peak width for fit. It is in unit of number of pixels.");
106 "Maximum guessed peak width for fit. It is in unit of number of pixels.");
109 "Step of guessed peak width. It is in unit of number of pixels.");
111 auto mustBePositiveDBL = std::make_shared<BoundedValidator<double>>();
113 "Tolerance on the found peaks' positions against the input "
114 "peak positions. Non-positive value indicates that this "
115 "option is turned off.");
119 "The name of the TableWorkspace in which to store the list "
123 "false generates table with effective centre/width/height "
124 "parameters. true generates a table with peak function "
127 declareProperty(
"MinimumPeakHeight", DBL_MIN,
"Minimum allowed peak height. ");
130 "Least value of the maximum observed Y value of a peak within "
131 "specified region. If any peak's maximum observed Y value is smaller, "
133 "this peak will not be fit. It is designed for EventWorkspace with "
136 std::array<std::string, 2> costFuncOptions = {{
"Chi-Square",
"Rwp"}};
144 "Minimizer to use for fitting. Minimizers available are "
145 "\"Levenberg-Marquardt\", \"Simplex\","
146 "\"Conjugate gradient (Fletcher-Reeves imp.)\", \"Conjugate "
147 "gradient (Polak-Ribiere imp.)\", \"BFGS\", and "
148 "\"Levenberg-MarquardtMD\"");
150 declareProperty(
"StartFromObservedPeakCentre",
true,
"Use observed value as the starting value of peak centre. ");
170 throw std::invalid_argument(
"Number of FitWindows must be exactly "
171 "twice the number of PeakPositions");
198 if (wsIndex >=
static_cast<int>(
m_dataWS->getNumberHistograms())) {
199 g_log.
warning() <<
"The value of WorkspaceIndex provided (" << wsIndex
200 <<
") is larger than the size of this workspace (" <<
m_dataWS->getNumberHistograms() <<
")\n";
202 "FindPeaks WorkspaceIndex property");
204 m_indexSet =
m_dataWS->indexInfo().makeIndexSet({
static_cast<Indexing::GlobalSpectrumIndex
>(wsIndex)});
214 if (t1 > t2 || t1 <= 0 || t3 <= 0) {
215 std::stringstream errss;
216 errss <<
"User specified wrong guessed peak width parameters (must be "
217 "postive and make sense). "
218 <<
"User inputs are min = " << t1 <<
", max = " << t2 <<
", step = " << t3;
220 throw std::runtime_error(errss.str());
283 for (
size_t i = 0; i < numpeakpars; ++i)
285 for (
size_t i = 0; i < numbkgdpars; ++i)
311 const std::vector<double> &fitwindows) {
312 bool useWindows = (!fitwindows.empty());
313 std::size_t numPeaks = peakcentres.size();
319 const auto &vecX =
m_dataWS->x(spec);
321 double practical_x_min = vecX.front();
322 double practical_x_max = vecX.back();
323 g_log.
information() <<
"actual x-range = [" << practical_x_min <<
" -> " << practical_x_max <<
"]\n";
325 const auto &vecY =
m_dataWS->y(spec);
326 const auto &vecE =
m_dataWS->e(spec);
327 const size_t numY = vecY.size();
329 for (; i_min < numY; ++i_min) {
330 if ((vecY[i_min] != 0.) || (vecE[i_min] != 0)) {
335 practical_x_min = vecX[i_min];
337 size_t i_max = numY - 2;
338 for (; i_max > i_min; --i_max) {
339 if ((vecY[i_max] != 0.) || (vecE[i_max] != 0)) {
344 g_log.
debug() <<
"Finding peaks from giving starting point, with interval i_min = " << i_min
345 <<
" i_max = " << i_max <<
'\n';
346 practical_x_max = vecX[i_max];
348 g_log.
information() <<
"practical x-range = [" << practical_x_min <<
" -> " << practical_x_max <<
"]\n";
350 for (std::size_t ipeak = 0; ipeak < numPeaks; ipeak++) {
352 double x_center = peakcentres[ipeak];
354 std::stringstream infoss;
355 infoss <<
"Spectrum " << spec <<
": Fit peak @ d = " << x_center;
357 infoss <<
" inside fit window [" << fitwindows[2 * ipeak] <<
", " << fitwindows[2 * ipeak + 1] <<
"]";
362 if (x_center > practical_x_min && x_center < practical_x_max) {
366 bool hasLeftPeak = (ipeak > 0);
367 double leftpeakcentre = 0.;
369 leftpeakcentre = peakcentres[ipeak - 1];
371 bool hasRightPeak = (ipeak < numPeaks - 1);
372 double rightpeakcentre = 0.;
374 rightpeakcentre = peakcentres[ipeak + 1];
377 hasRightPeak, rightpeakcentre);
380 g_log.
warning() <<
"Given peak centre " << x_center <<
" is out side of given data's range (" << practical_x_min
381 <<
", " << practical_x_max <<
").\n";
412 const double kz = 1.22;
420 for (
size_t k_out = 0; k_out <
m_indexSet.size(); ++k_out) {
422 const auto &S = smoothedData[k_out].y();
423 const auto &F = smoothedData[k_out].e();
426 int i0 = 0, i1 = 0, i2 = 0, i3 = 0, i4 = 0, i5 = 0;
427 for (
int i = 1; i < static_cast<int>(S.size()); ++i) {
433 S[i] > 0 ? M = 2 : M = 3;
436 if (S[i - 1] > F[i - 1]) {
453 }
else if (S[i - 1] > 0) {
485 if (i5 && i1 && i2 && i3)
489 double num = 0.0, denom = 0.0;
490 for (
int j = i3; j <= i5; ++j) {
498 i0 =
static_cast<int>(num / denom);
502 if (i1 > i2 || i2 > i3 || i3 > i4 || i5 <= i4) {
510 if (std::abs(S[i4]) < 2 * F[i4]) {
520 int n2 = std::abs(boost::math::iround(0.5 * (F[i0] / S[i0]) * (n1 +
tolerance)));
521 const int n2b = std::abs(boost::math::iround(0.5 * (F[i0] / S[i0]) * (n1 -
tolerance)));
525 const int testVal = n2 ? n2 : 1;
526 if (i3 - i2 - 1 > testVal) {
531 int n3 = std::abs(boost::math::iround((n1 +
tolerance) * (1 - 2 * (F[i0] / S[i0]))));
532 const int n3b = std::abs(boost::math::iround((n1 -
tolerance) * (1 - 2 * (F[i0] / S[i0]))));
536 if (i2 - i1 + 1 < n3) {
542 g_log.
debug() <<
"Spectrum=" << k <<
" i0=" << i0 <<
" X=" <<
m_dataWS->x(k)[i0] <<
" i1=" << i1 <<
" i2=" << i2
543 <<
" i3=" << i3 <<
" i4=" << i4 <<
" i5=" << i5 <<
'\n';
547 auto wssize =
static_cast<int>(
m_dataWS->x(k).size());
549 int iwidth = i0 - i2;
555 i_min = i4 - 5 * iwidth;
557 int i_max = i4 + 5 * iwidth;
564 i1 = 0, i2 = 0, i3 = 0, i4 = 0, i5 = 0;
584 std::vector<Histogram> diffed;
588 diffed.emplace_back(input->histogram(i));
589 diffed.back().mutableY() = 0.0;
590 diffed.back().mutableE() = 0.0;
592 const auto &
Y = input->y(i);
593 auto &S = diffed.back().mutableY();
597 for (
size_t j = 1; j < S.size() - 1; ++j) {
598 S[j] =
Y[j - 1] - 2 *
Y[j] +
Y[j + 1];
613 for (
auto &histogram : histograms)
614 for (
int i = 0; i <
g_z; ++i)
615 histogram =
smooth(histogram, w);
629 std::vector<HistogramData::Histogram> &smoothed,
const int &w) {
632 static_assert(
g_z == 5,
"Value of z has changed!");
634 const auto factor =
static_cast<int>(std::pow(
static_cast<double>(w),
g_z));
636 const double constant = sqrt(
static_cast<double>(this->
computePhi(w))) / factor;
638 for (
size_t i = 0; i <
m_indexSet.size(); ++i) {
640 smoothed[i].mutableE() = input->e(i_in) * constant;
654 const int m = (w - 1) / 2;
656 int max_index_prev = 1;
658 std::vector<long long> previous(n_el_prev);
665 return std::accumulate(previous.begin(), previous.end(),
static_cast<long long>(0),
668 std::vector<long long> next;
672 int max_index = zz *
m + 1;
673 int n_el = 2 * max_index + 1;
675 std::fill(next.begin(), next.end(), 0);
676 for (
int i = 0; i < n_el; ++i) {
677 int delta = -max_index + i;
679 int index = l + max_index_prev;
681 next[i] += previous[
index];
684 previous.resize(n_el);
685 std::copy(next.begin(), next.end(), previous.begin());
686 max_index_prev = max_index;
690 const long long retval = std::accumulate(previous.begin(), previous.end(),
static_cast<long long>(0),
692 g_log.
debug() <<
"FindPeaks::computePhi - calculated value = " << retval <<
"\n";
705 if (
x <= vecX.front()) {
708 }
else if (
x >= vecX.back()) {
710 index =
static_cast<int>(vecX.size()) - 1;
713 index =
static_cast<int>(std::lower_bound(vecX.begin(), vecX.end(),
x) - vecX.begin());
717 std::stringstream errss;
718 errss <<
"Returned index = 0 for x = " <<
x <<
" with X[0] = " << vecX[0]
719 <<
". This situation is ruled out in this algorithm.";
721 throw std::runtime_error(errss.str());
723 std::stringstream errss;
724 errss <<
"Returned x = " <<
x <<
" is not between " << vecX[
index - 1] <<
" and " << vecX[
index]
725 <<
", which are returned by lower_bound.";
727 throw std::runtime_error(errss.str());
756 const int fitWidth,
const bool hasleftpeak,
const double leftpeakcentre,
757 const bool hasrightpeak,
const double rightpeakcentre) {
759 const auto &vecX = input->x(wsIndex);
760 const auto &vecY = input->y(wsIndex);
763 int i_centre = this->
getIndex(vecX, center_guess);
766 int i_min = i_centre - 5 * fitWidth;
772 double xmin = vecX[i_min];
773 double peaksepline = center_guess - (center_guess - leftpeakcentre) * 0.66;
774 if (xmin < peaksepline)
775 i_min =
getIndex(vecX, peaksepline);
779 int i_max = i_centre + 5 * fitWidth;
780 if (i_max >=
static_cast<int>(vecX.size()) - 1)
781 i_max =
static_cast<int>(vecY.size()) - 2;
785 double xmax = vecX[i_max];
786 double peaksepline = center_guess + (rightpeakcentre - center_guess) * 0.66;
787 if (xmax > peaksepline)
788 i_max =
getIndex(vecX, peaksepline);
792 if (i_max - i_min <= 0)
793 throw std::runtime_error(
"Impossible to i_min >= i_max.");
795 std::stringstream outss;
796 outss <<
"Fit peak with guessed FWHM: starting center = " << center_guess <<
", FWHM = " << fitWidth
797 <<
". Estimated peak fit window from given FWHM: " << vecX[i_min] <<
", " << vecX[i_max];
814 const double xmin,
const double xmax) {
816 g_log.
information() <<
"Fit Peak with given window: Guessed center = " << centre_guess <<
" x-min = " << xmin
817 <<
", x-max = " << xmax <<
"\n";
818 if (xmin >= centre_guess || xmax <= centre_guess) {
819 g_log.
warning(
"Peak centre is on the edge of Fit window. ");
825 const auto &vecX = input->x(wsIndex);
828 int i_centre = this->
getIndex(vecX, centre_guess);
832 if (i_min >= i_centre) {
833 g_log.
warning() <<
"Input peak centre @ " << centre_guess <<
" is out side of minimum x = " << xmin
834 <<
". Input X ragne = " << vecX.front() <<
", " << vecX.back() <<
"\n";
841 if (i_max < i_centre) {
842 g_log.
warning() <<
"Input peak centre @ " << centre_guess <<
" is out side of maximum x = " << xmax <<
"\n";
856 const int i_max,
const int i_centre) {
857 const auto &vecX = input->x(spectrum);
858 const auto &vecY = input->y(spectrum);
861 bool hasHighCounts =
false;
862 for (
int i = i_min; i <= i_max; ++i)
864 hasHighCounts =
true;
867 if (!hasHighCounts) {
868 std::stringstream ess;
869 ess <<
"Peak supposed at " << vecY[i_centre] <<
" does not have enough counts as " <<
m_leastMaxObsY;
876 std::stringstream outss;
877 outss <<
"Fit single peak in X-range " << vecX[i_min] <<
", " << vecX[i_max] <<
", centre at " << vecX[i_centre]
878 <<
" (index = " << i_centre <<
"). ";
883 std::vector<double> vecbkgdparvalue(3, 0.);
884 std::vector<double> vecpeakrange(3, 0.);
885 int usefpdresult =
findPeakBackground(input, spectrum, i_min, i_max, vecbkgdparvalue, vecpeakrange);
886 if (usefpdresult < 0) {
891 for (
size_t i = 0; i < vecbkgdparvalue.size(); ++i)
896 double est_height(0.0), est_fwhm(0.0);
897 size_t i_obscentre(0);
898 double est_leftfwhm(0.0), est_rightfwhm(0.0);
899 std::string errmsg =
estimatePeakParameters(vecX, vecY, i_min, i_max, vecbkgdparvalue, i_obscentre, est_height,
900 est_fwhm, est_leftfwhm, est_rightfwhm);
901 if (!errmsg.empty()) {
903 i_obscentre = i_centre;
917 if (usefpdresult < 0) {
921 i_obscentre = i_centre;
922 estimatePeakRange(vecX, i_obscentre, i_min, i_max, est_leftfwhm, est_rightfwhm, vecpeakrange);
928 std::vector<double> fitwindow(2);
929 fitwindow[0] = vecX[i_min];
930 fitwindow[1] = vecX[i_max];
935 bool fitsuccess =
false;
958 std::vector<double> &vecBkgdParamValues, std::vector<double> &vecpeakrange) {
959 const auto &vecX = input->x(spectrum);
963 estimate->setLoggingOffset(1);
964 estimate->setProperty(
"InputWorkspace", input);
965 estimate->setProperty(
"WorkspaceIndex", spectrum);
967 std::vector<double> fwvec;
968 fwvec.emplace_back(vecX[i_min]);
969 fwvec.emplace_back(vecX[i_max]);
971 estimate->setProperty(
"FitWindow", fwvec);
972 estimate->executeAsChildAlg();
978 if (peaklisttablews->columnCount() < 7)
979 throw std::runtime_error(
"No 7th column for use FindPeakBackground result or not. ");
981 if (peaklisttablews->rowCount() > 0) {
987 const int hiddenFitresult = peaklisttablews->Int(0, 6);
992 vecpeakrange.resize(2);
995 size_t i_peakmin, i_peakmax;
996 i_peakmin = peaklisttablews->Int(0, 1);
997 i_peakmax = peaklisttablews->Int(0, 2);
1000 <<
"iMin = " << i_min <<
", iPeakMin = " << i_peakmin <<
", iPeakMax = " << i_peakmax
1001 <<
", iMax = " << i_max <<
"\n";
1003 if (i_peakmin < i_peakmax && i_peakmin > i_min + 2 && i_peakmax < i_max - 2) {
1006 double bg0, bg1, bg2;
1007 bg0 = peaklisttablews->Double(0, 3);
1008 bg1 = peaklisttablews->Double(0, 4);
1009 bg2 = peaklisttablews->Double(0, 5);
1012 vecBkgdParamValues.resize(3, 0.);
1013 vecBkgdParamValues[0] = bg0;
1014 vecBkgdParamValues[1] = bg1;
1015 vecBkgdParamValues[2] = bg2;
1017 g_log.
information() <<
"Background parameters (from FindPeakBackground) A0=" << bg0 <<
", A1=" << bg1
1018 <<
", A2=" << bg2 <<
"\n";
1020 vecpeakrange[0] = vecX[i_peakmin];
1021 vecpeakrange[1] = vecX[i_peakmax];
1024 g_log.
debug(
"FindPeakBackground result is ignored due to wrong in peak range. ");
1030 std::stringstream outx;
1031 outx <<
"FindPeakBackground Result: Given window (" << vecX[i_min] <<
", " << vecX[i_max]
1032 <<
"); Determine peak range: (" << vecpeakrange[0] <<
", " << vecpeakrange[1] <<
"). ";
1056 size_t i_max,
const std::vector<double> &vecbkgdparvalues,
1057 size_t &iobscentre,
double &
height,
double &fwhm,
double &leftfwhm,
1058 double &rightfwhm) {
1060 const double bg0 = vecbkgdparvalues[0];
1063 if (vecbkgdparvalues.size() >= 2) {
1064 bg1 = vecbkgdparvalues[1];
1065 if (vecbkgdparvalues.size() >= 3)
1066 bg2 = vecbkgdparvalues[2];
1071 const double tmpx = vecX[i_min];
1072 height = vecY[i_min] - (bg0 + bg1 * tmpx + bg2 * tmpx * tmpx);
1076 if (i_max == vecY.size())
1080 for (
size_t i = i_min + 1; i <= i_max; ++i) {
1081 const double x = vecX[i];
1082 const double tmpheight = vecY[i] - (bg0 + bg1 *
x + bg2 *
x *
x);
1084 if (tmpheight >
height) {
1087 }
else if (tmpheight < lowest) {
1093 double obscentre = vecX[iobscentre];
1094 double drop =
height - lowest;
1099 return "Flat spectrum";
1104 return "Fluctuation is less than minimum allowed value.";
1109 if (iobscentre <= i_min + 1 || iobscentre >= i_max - 1) {
1110 std::stringstream dbss;
1111 dbss <<
"Maximum value on edge. Fit window is between " << vecX[i_min] <<
" and " << vecX[i_max]
1112 <<
". Maximum value " << vecX[iobscentre] <<
" is located on (" << iobscentre <<
").";
1120 for (
int i =
static_cast<int>(iobscentre) - 1; i >= 0; --i) {
1121 double xleft = vecX[i];
1122 double yleft = vecY[i] - (bg0 + bg1 * xleft + bg2 * xleft * xleft);
1123 if (yleft <= 0.5 *
height) {
1126 leftfwhm = obscentre - 0.5 * (vecX[i] + vecX[i + 1]);
1133 for (
size_t i = iobscentre + 1; i <= i_max; ++i) {
1134 double xright = vecX[i];
1135 double yright = vecY[i] - (bg0 + bg1 * xright + bg2 * xright * xright);
1136 if (yright <= 0.5 *
height) {
1137 rightfwhm = 0.5 * (vecX[i] + vecX[i - 1]) - obscentre;
1142 if (leftfwhm <= 0 || rightfwhm <= 0) {
1143 std::stringstream errmsg;
1144 errmsg <<
"Estimate peak parameters error (FWHM cannot be zero): Input "
1146 << vecX.size() <<
", Xmin = " << vecX[i_min] <<
"(" << i_min <<
"), Xmax = " << vecX[i_max] <<
"(" << i_max
1148 <<
"Estimated peak centre @ " << vecX[iobscentre] <<
"(" << iobscentre <<
") with height = " <<
height
1149 <<
"; Lowest Y value = " << lowest <<
"; Output error: . leftfwhm = " << leftfwhm
1150 <<
", right fwhm = " << rightfwhm <<
".";
1151 return errmsg.str();
1154 fwhm = leftfwhm + rightfwhm;
1157 std::stringstream errmsg;
1158 errmsg <<
"Estimate peak parameters error (FWHM cannot be zero): Input "
1160 << vecX.size() <<
", Xmin = " << vecX[i_min] <<
"(" << i_min <<
"), Xmax = " << vecX[i_max] <<
"(" << i_max
1162 <<
"Estimated peak centre @ " << vecX[iobscentre] <<
"(" << iobscentre <<
") with height = " <<
height
1163 <<
"; Lowest Y value = " << lowest <<
"; Output error: . fwhm = " << fwhm <<
".";
1164 return errmsg.str();
1168 <<
", FWHM = " << fwhm <<
" = " << leftfwhm <<
" + " << rightfwhm <<
".\n";
1170 return std::string();
1185 std::vector<double> &vecbkgdparvalues) {
1188 throw std::runtime_error(
"when trying to estimate the background parameter "
1189 "values: i_min cannot larger or equal to i_max");
1190 if (vecbkgdparvalues.size() < 3)
1191 vecbkgdparvalues.resize(3, 0.);
1207 for (
size_t i = 0; i < numavg; ++i) {
1213 yf +=
Y[i_max - i - 1];
1215 x0 = x0 /
static_cast<double>(numavg);
1216 y0 = y0 /
static_cast<double>(numavg);
1217 xf = xf /
static_cast<double>(numavg);
1218 yf = yf /
static_cast<double>(numavg);
1221 vecbkgdparvalues[2] = 0.;
1224 vecbkgdparvalues[1] = (y0 - yf) / (x0 - xf);
1225 vecbkgdparvalues[0] = (xf * y0 - x0 * yf) / (xf - x0);
1228 vecbkgdparvalues[1] = 0.;
1229 vecbkgdparvalues[0] = 0.5 * (y0 + yf);
1238 const double &leftfwhm,
const double &rightfwhm, std::vector<double> &vecpeakrange) {
1240 if (vecpeakrange.size() < 2)
1241 vecpeakrange.resize(2, 0.);
1243 if (i_centre < i_min || i_centre > i_max)
1244 throw std::runtime_error(
"Estimate peak range input centre is out of fit window. ");
1247 double peakleftbound = vecX[i_centre] - 6. * leftfwhm;
1248 double peakrightbound = vecX[i_centre] + 6. * rightfwhm;
1251 auto ipeakleft =
static_cast<size_t>(
getIndex(vecX, peakleftbound));
1252 if (ipeakleft <= i_min) {
1253 size_t numbkgdpts = (i_centre - i_min) / 6;
1257 ipeakleft = i_min + numbkgdpts;
1258 if (ipeakleft >= i_centre)
1259 ipeakleft = i_min + 1;
1261 peakleftbound = vecX[ipeakleft];
1264 auto ipeakright =
static_cast<size_t>(
getIndex(vecX, peakrightbound));
1265 if (ipeakright >= i_max) {
1266 size_t numbkgdpts = (i_max - i_centre) / 6;
1270 ipeakright = i_max - numbkgdpts;
1271 if (ipeakright <= i_centre)
1272 ipeakright = i_max - 1;
1274 peakrightbound = vecX[ipeakright];
1278 vecpeakrange[0] = peakleftbound;
1279 vecpeakrange[1] = peakrightbound;
1293 const double mincost) {
1295 if (mincost < 0. || mincost >= DBL_MAX - 1.0E-10)
1296 throw std::runtime_error(
"Minimum cost indicates that fit fails. This "
1297 "method should not be called "
1298 "under this circumstance. ");
1304 t << static_cast<int>(spectrum);
1309 size_t nparams = peakfunction->nParams();
1310 size_t nparamsb = bkgdfunction->nParams();
1313 if (nparams + nparamsb + 2 != numcols) {
1314 throw std::runtime_error(
"Error 1307 number of columns do not matches");
1317 for (
size_t i = 0; i < nparams; ++i) {
1318 t << peakfunction->getParameter(i);
1320 for (
size_t i = 0; i < nparamsb; ++i) {
1321 t << bkgdfunction->getParameter(i);
1326 double peakcentre = peakfunction->centre();
1327 double fwhm = peakfunction->fwhm();
1328 double height = peakfunction->height();
1330 t << peakcentre << fwhm <<
height;
1335 double a0 = bkgdfunction->getParameter(
"A0");
1337 if (bkgdfunction->name() !=
"FlatBackground")
1338 a1 = bkgdfunction->getParameter(
"A1");
1340 if (bkgdfunction->name() !=
"LinearBackground" && bkgdfunction->name() !=
"FlatBackground")
1341 a2 = bkgdfunction->getParameter(
"A2");
1343 t << a0 << a1 << a2;
1345 g_log.
debug() <<
"Peak parameters found: cen=" << peakcentre <<
" fwhm=" << fwhm <<
" height=" <<
height
1346 <<
" a0=" << a0 <<
" a1=" << a1 <<
" a2=" << a2;
1348 g_log.
debug() <<
" chsq=" << mincost <<
"\n";
1364 t << static_cast<int>(spectrum);
1385 std::string backgroundposix;
1388 backgroundposix =
"Background";
1411 const std::vector<double> &vec_fitwindow,
const std::vector<double> &vec_peakrange,
1412 int minGuessFWHM,
int maxGuessFWHM,
int guessedFWHMStep,
int estBackResult) {
1413 std::stringstream dbss;
1414 dbss <<
"[Call FitPeak] Fit 1 peak at X = " << peakfunction->centre() <<
" of spectrum " << wsindex;
1418 bool fitwithsteppedfwhm = (guessedFWHMStep > 0);
1423 fitpeak.
setFitWindow(vec_fitwindow[0], vec_fitwindow[1]);
1425 fitpeak.
setFunctions(peakfunction, backgroundfunction);
1426 fitpeak.
setupGuessedFWHM(userFWHM, minGuessFWHM, maxGuessFWHM, guessedFWHMStep, fitwithsteppedfwhm);
1427 fitpeak.
setPeakRange(vec_peakrange[0], vec_peakrange[1]);
1429 if (estBackResult == 1) {
1444 return costfuncvalue;
1454 std::vector<double> vec_value(numpeakpars);
1455 for (
size_t i = 0; i < numpeakpars; ++i) {
1457 vec_value[i] = parvalue;
#define DECLARE_ALGORITHM(classname)
CostFunctions::CostFuncFitting & m_costFunction
The cost function.
std::map< DeltaEMode::Type, std::string > index
void declareProperty(std::unique_ptr< Kernel::Property > p, const std::string &doc="") override
Add a property to the list of managed properties.
std::string getPropertyValue(const std::string &name) const override
Get the value of a property as a string.
TypedValue getProperty(const std::string &name) const override
Get the value of a property.
virtual std::shared_ptr< Algorithm > createChildAlgorithm(const std::string &name, const double startProgress=-1., const double endProgress=-1., const bool enableLogging=true, const int &version=-1)
Create a Child Algorithm.
static bool isEmpty(const NumT toCheck)
checks that the value was not set by users, uses the value in empty double/int.
void setChild(const bool isChild) override
To set whether algorithm is a child.
An interface to a peak function, which extend the interface of IFunctionWithLocation by adding method...
Base class for algorithms that can run in parallel on all MPI ranks but not in a distributed fashion.
TableRow represents a row in a TableWorkspace.
A property class for workspaces.
This algorithm searches for peaks in a dataset.
static const int g_z
The number of smoothing iterations.
Indexing::SpectrumIndexSet m_indexSet
list of workspace indicies to check
API::ITableWorkspace_sptr m_outPeakTableWS
Storage of the peak data.
double m_peakPositionTolerance
void init() override
Initialize and declare properties.
API::IBackgroundFunction_sptr m_backgroundFunction
bool m_useObsCentre
Start values.
std::vector< double > m_vecPeakCentre
void generateOutputPeakParameterTable()
Generate a table workspace for output peak parameters.
std::unique_ptr< API::Progress > m_progress
Progress reporting.
void fitPeakGivenFWHM(const API::MatrixWorkspace_sptr &input, const int spectrum, const double center_guess, const int fitWidth, const bool hasleftpeak, const double leftpeakcentre, const bool hasrightpeak, const double rightpeakcentre)
Fit peak by given/guessed FWHM.
int getIndex(const HistogramData::HistogramX &vecX, double x)
needed by FindPeaksBackground
int findPeakBackground(const API::MatrixWorkspace_sptr &input, int spectrum, size_t i_min, size_t i_max, std::vector< double > &vecBkgdParamValues, std::vector< double > &vecpeakrange)
Find peak background.
void processAlgorithmProperties()
Process algorithm's properties.
void createFunctions()
Create peak and background functions.
void fitPeakInWindow(const API::MatrixWorkspace_sptr &input, const int spectrum, const double centre_guess, const double xmin, const double xmax)
Fit peak confined in a given window (x-min, x-max)
int m_maxGuessedPeakWidth
bool m_usePeakPositionTolerance
void estimatePeakRange(const HistogramData::HistogramX &vecX, size_t i_centre, size_t i_min, size_t i_max, const double &leftfwhm, const double &rightfwhm, std::vector< double > &vecpeakrange)
Estimate peak range based on background peak parameter.
API::MatrixWorkspace_sptr m_dataWS
workspace to check for peaks
double m_leastMaxObsY
Minimum value of peak's observed maximum Y value.
double m_minHeight
Minimum peak height.
std::size_t m_centreIndex
API::IPeakFunction_sptr m_peakFunction
std::size_t m_numTableParams
parameters or effective (centre, width, height)
double callFitPeak(const API::MatrixWorkspace_sptr &dataws, int wsindex, const API::IPeakFunction_sptr &peakfunction, const API::IBackgroundFunction_sptr &backgroundfunction, const std::vector< double > &vec_fitwindow, const std::vector< double > &vec_peakrange, int minGuessFWHM, int maxGuessFWHM, int guessedFWHMStep, int estBackResult=0)
Fit peak by calling 'FitPeak'.
std::vector< double > m_vecFitWindows
std::string estimatePeakParameters(const HistogramData::HistogramX &vecX, const HistogramData::HistogramY &vecY, size_t i_min, size_t i_max, const std::vector< double > &vecbkgdparvalues, size_t &iobscentre, double &height, double &fwhm, double &leftfwhm, double &rightfwhm)
Estimate peak parameters.
void calculateStandardDeviation(const API::MatrixWorkspace_const_sptr &input, std::vector< HistogramData::Histogram > &smoothed, const int &w)
Calculates the statistical error on the smoothed data.
long long computePhi(const int &w) const
Calculates the coefficient phi which goes into the calculation of the error on the smoothed data Uses...
void exec() override
Execute the findPeaks algorithm.
bool m_highBackground
flag for find relatively weak peak in high
void findPeaksGivenStartingPoints(const std::vector< double > &peakcentres, const std::vector< double > &fitwindows)
Find peaks according to given peak positions.
void findPeaksUsingMariscotti()
Find peaks by searching peak position using Mariscotti algorithm.
void fitSinglePeak(const API::MatrixWorkspace_sptr &input, const int spectrum, const int i_min, const int i_max, const int i_centre)
Fit peak: this is a basic peak fit function as a root function for all different type of user input.
std::vector< std::string > m_bkgdParameterNames
std::vector< std::string > m_peakParameterNames
void addNonFitRecord(const size_t spectrum, const double centre)
Add the fit record (failure) to output workspace.
std::vector< double > getStartingPeakValues()
Get the peak parameter values from m_peakFunction and output to a list in the same order of m_peakPar...
void estimateBackground(const HistogramData::HistogramX &X, const HistogramData::HistogramY &Y, const size_t i_min, const size_t i_max, std::vector< double > &vecbkgdparvalues)
Estimate background of a given range.
std::vector< HistogramData::Histogram > calculateSecondDifference(const API::MatrixWorkspace_const_sptr &input)
Methods searving for findPeaksUsingMariscotti()
int m_minGuessedPeakWidth
std::string m_peakFuncType
void smoothData(std::vector< HistogramData::Histogram > &histograms, const int w, const int g_z)
Smooth data for Mariscotti.
void addInfoRow(const size_t spectrum, const API::IPeakFunction_const_sptr &peakfunction, const API::IBackgroundFunction_sptr &bkgdfunction, const bool isoutputraw, const double mincost)
Add a new row in output TableWorkspace containing information of the fitted peak+background.
int m_stepGuessedPeakWidth
bool m_rawPeaksTable
background
std::string m_costFunction
std::string m_backgroundType
int m_inputPeakFWHM
holder for the requested peak FWHM
FitOneSinglePeak: a class to perform peak fitting on a single peak.
void setPeakRange(double xpeakleft, double xpeakright)
Set peak range.
void setupGuessedFWHM(double usrwidth, int minfwhm, int maxfwhm, int stepsize, bool fitwithsteppedfwhm)
Set peak width to guess.
std::string getDebugMessage()
Get debug message.
void setFunctions(const API::IPeakFunction_sptr &peakfunc, const API::IBackgroundFunction_sptr &bkgdfunc)
Set functions.
void highBkgdFit()
Fit peak first considering high background.
bool simpleFit()
Fit peak and background together.
void setFittingMethod(std::string minimizer, const std::string &costfunction)
Set fitting method.
void setFitWindow(double leftwindow, double rightwindow)
Set fit range.
void setWorskpace(const API::MatrixWorkspace_sptr &dataws, size_t wsindex)
Set workspaces.
double getFitCostFunctionValue()
Get cost function value from fitting.
Support for a property that holds an array of values.
Exception for index errors.
IPropertyManager * setProperty(const std::string &name, const T &value)
Templated method to set the value of a PropertyWithValue.
ListValidator is a validator that requires the value of a property to be one of a defined list of pos...
void debug(const std::string &msg)
Logs at debug 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...
StartsWithValidator is a validator that requires the value of a property to start with one of the str...
std::shared_ptr< IBackgroundFunction > IBackgroundFunction_sptr
std::shared_ptr< IPeakFunction > IPeakFunction_sptr
std::shared_ptr< ITableWorkspace > ITableWorkspace_sptr
shared pointer to Mantid::API::ITableWorkspace
std::shared_ptr< const IPeakFunction > IPeakFunction_const_sptr
std::shared_ptr< const MatrixWorkspace > MatrixWorkspace_const_sptr
shared pointer to the matrix workspace base class (const version)
std::shared_ptr< MatrixWorkspace > MatrixWorkspace_sptr
shared pointer to the matrix workspace base class
HistogramData::Histogram smooth(const HistogramData::Histogram &histogram, int npts)
std::shared_ptr< IValidator > IValidator_sptr
A shared_ptr to an IValidator.
Helper class which provides the Collimation Length for SANS instruments.
constexpr int EMPTY_INT() noexcept
Returns what we consider an "empty" integer within a property.
constexpr double EMPTY_DBL() noexcept
Returns what we consider an "empty" double within a property.
@ Input
An input workspace.
@ Output
An output workspace.
Functor to accumulate a sum of squares.