Mantid
Loading...
Searching...
No Matches
FindPeaks.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 +
9
13#include "MantidAPI/TableRow.h"
17#include "MantidIndexing/GlobalSpectrumIndex.h"
18#include "MantidIndexing/IndexInfo.h"
22
25#include <boost/algorithm/string.hpp>
26#include <boost/math/special_functions/round.hpp>
27#include <numeric>
28
29#include <fstream>
30
31using namespace Mantid;
32using namespace Mantid::Kernel;
33using namespace Mantid::API;
34using namespace Mantid::DataObjects;
35using namespace Mantid::HistogramData;
36using namespace Mantid::Indexing;
37
38// const double MINHEIGHT = 2.00000001;
39
40namespace Mantid::Algorithms {
41
42// Register the algorithm into the AlgorithmFactory
43DECLARE_ALGORITHM(FindPeaks)
44
45//----------------------------------------------------------------------------------------------
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) /* for Gaussian */, 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) {}
56
57//----------------------------------------------------------------------------------------------
61 declareProperty(std::make_unique<WorkspaceProperty<>>("InputWorkspace", "", Direction::Input),
62 "Name of the workspace to search");
63
64 auto mustBeNonNegative = std::make_shared<BoundedValidator<int>>();
65 mustBeNonNegative->setLower(0);
66 declareProperty("WorkspaceIndex", EMPTY_INT(), mustBeNonNegative,
67 "If set, only this spectrum will be searched for peaks "
68 "(otherwise all are)");
69
70 auto min = std::make_shared<BoundedValidator<int>>();
71 min->setLower(1);
72 // The estimated width of a peak in terms of number of channels
73 declareProperty("FWHM", 7, min, "Estimated number of points covered by the fwhm of a peak (default 7)");
74
75 // The tolerance allowed in meeting the conditions
76 declareProperty("Tolerance", 4, min,
77 "A measure of the strictness desired in "
78 "meeting the condition on peak "
79 "candidates,\n"
80 "Mariscotti recommends 2 (default 4)");
81
82 declareProperty(std::make_unique<ArrayProperty<double>>("PeakPositions"),
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.");
86
87 declareProperty(std::make_unique<ArrayProperty<double>>("FitWindows"),
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.");
91
92 std::vector<std::string> peakNames = FunctionFactory::Instance().getFunctionNames<API::IPeakFunction>();
93 declareProperty("PeakFunction", "Gaussian", std::make_shared<StringListValidator>(peakNames));
94
95 std::vector<std::string> bkgdtypes{"Flat", "Linear", "Quadratic"};
96 declareProperty("BackgroundType", "Linear", std::make_shared<StringListValidator>(bkgdtypes), "Type of Background.");
97
98 declareProperty("HighBackground", true, "Relatively weak peak in high background");
99
100 auto mustBePositive = std::make_shared<BoundedValidator<int>>();
101 mustBePositive->setLower(1);
102 declareProperty("MinGuessedPeakWidth", 2, mustBePositive,
103 "Minimum guessed peak width for fit. It is in unit of number of pixels.");
104
105 declareProperty("MaxGuessedPeakWidth", 10, mustBePositive,
106 "Maximum guessed peak width for fit. It is in unit of number of pixels.");
107
108 declareProperty("GuessedPeakWidthStep", 2, mustBePositive,
109 "Step of guessed peak width. It is in unit of number of pixels.");
110
111 auto mustBePositiveDBL = std::make_shared<BoundedValidator<double>>();
112 declareProperty("PeakPositionTolerance", EMPTY_DBL(), mustBePositiveDBL,
113 "Tolerance on the found peaks' positions against the input "
114 "peak positions. Non-positive value indicates that this "
115 "option is turned off.");
116
117 // The found peaks in a table
119 "The name of the TableWorkspace in which to store the list "
120 "of peaks found");
121
122 declareProperty("RawPeakParameters", false,
123 "false generates table with effective centre/width/height "
124 "parameters. true generates a table with peak function "
125 "parameters");
126
127 declareProperty("MinimumPeakHeight", DBL_MIN, "Minimum allowed peak height. ");
128
129 declareProperty("MinimumPeakHeightObs", 0.0,
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, "
132 "then "
133 "this peak will not be fit. It is designed for EventWorkspace with "
134 "integer counts.");
135
136 std::array<std::string, 2> costFuncOptions = {{"Chi-Square", "Rwp"}};
137 declareProperty("CostFunction", "Chi-Square",
138 Kernel::IValidator_sptr(new Kernel::ListValidator<std::string>(costFuncOptions)), "Cost functions");
139
140 std::vector<std::string> minimizerOptions = API::FuncMinimizerFactory::Instance().getKeys();
141
142 declareProperty("Minimizer", "Levenberg-MarquardtMD",
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\"");
149
150 declareProperty("StartFromObservedPeakCentre", true, "Use observed value as the starting value of peak centre. ");
151}
152
153//----------------------------------------------------------------------------------------------
157 // Process input
159
160 // Create those functions to fit
162
163 // Set up output table workspace
165
166 // Fit
167 if (!m_vecPeakCentre.empty()) {
168 if (!m_vecFitWindows.empty()) {
169 if (m_vecFitWindows.size() != (m_vecPeakCentre.size() * 2)) {
170 throw std::invalid_argument("Number of FitWindows must be exactly "
171 "twice the number of PeakPositions");
172 }
173 }
174
175 // Perform fit with fixed start positions.
177 } else {
178 // Use Mariscotti's method to find the peak centers
181 }
182
183 // Set output properties
184 g_log.information() << "Total " << m_outPeakTableWS->rowCount() << " peaks found and successfully fitted.\n";
185 setProperty("PeaksList", m_outPeakTableWS);
186} // END: exec()
187
188//----------------------------------------------------------------------------------------------
192 // Input workspace
193 m_dataWS = getProperty("InputWorkspace");
194
195 // WorkspaceIndex
196 int wsIndex = getProperty("WorkspaceIndex");
197 if (!isEmpty(wsIndex)) {
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";
201 throw Kernel::Exception::IndexError(wsIndex, m_dataWS->getNumberHistograms() - 1,
202 "FindPeaks WorkspaceIndex property");
203 }
204 m_indexSet = m_dataWS->indexInfo().makeIndexSet({static_cast<Indexing::GlobalSpectrumIndex>(wsIndex)});
205 } else {
206 m_indexSet = m_dataWS->indexInfo().makeIndexSet();
207 }
208
209 // Peak width
211 int t1 = getProperty("MinGuessedPeakWidth");
212 int t2 = getProperty("MaxGuessedPeakWidth");
213 int t3 = getProperty("GuessedPeakWidthStep");
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;
219 g_log.warning(errss.str());
220 throw std::runtime_error(errss.str());
221 }
222
226
227 m_peakPositionTolerance = getProperty("PeakPositionTolerance");
231
232 // Specified peak positions, which is optional
233 m_vecPeakCentre = getProperty("PeakPositions");
234 if (!m_vecPeakCentre.empty())
235 std::sort(m_vecPeakCentre.begin(), m_vecPeakCentre.end());
236 m_vecFitWindows = getProperty("FitWindows");
237
238 // Peak and ground type
239 m_peakFuncType = getPropertyValue("PeakFunction");
240 m_backgroundType = getPropertyValue("BackgroundType");
241
242 // Fit algorithm
243 m_highBackground = getProperty("HighBackground");
244
245 // Peak parameters are give via a table workspace
246 m_rawPeaksTable = getProperty("RawPeakParameters");
247
248 // Minimum peak height
249 m_minHeight = getProperty("MinimumPeakHeight");
250
251 // About Fit
252 m_minimizer = getPropertyValue("Minimizer");
253 m_costFunction = getPropertyValue("CostFunction");
254
255 m_useObsCentre = getProperty("StartFromObservedPeakCentre");
256
257 m_leastMaxObsY = getProperty("MinimumPeakHeightObs");
258}
259
260//----------------------------------------------------------------------------------------------
264 m_outPeakTableWS = std::make_shared<TableWorkspace>();
265 m_outPeakTableWS->addColumn("int", "spectrum");
266
267 if (m_rawPeaksTable) {
268 // Output raw peak parameters
269 size_t numpeakpars = m_peakFunction->nParams();
270 size_t numbkgdpars = m_backgroundFunction->nParams();
271 m_numTableParams = numpeakpars + numbkgdpars;
272 if (m_peakFuncType == "Gaussian")
273 m_centreIndex = 1;
274 else if (m_peakFuncType == "LogNormal")
275 m_centreIndex = 1;
276 else if (m_peakFuncType == "Lorentzian")
277 m_centreIndex = 1;
278 else if (m_peakFuncType == "PseudoVoigt")
279 m_centreIndex = 2;
280 else
281 m_centreIndex = m_numTableParams; // bad value
282
283 for (size_t i = 0; i < numpeakpars; ++i)
284 m_outPeakTableWS->addColumn("double", m_peakParameterNames[i]);
285 for (size_t i = 0; i < numbkgdpars; ++i)
286 m_outPeakTableWS->addColumn("double", m_bkgdParameterNames[i]);
287 // m_outPeakTableWS->addColumn("double", "f1.A2");
288 } else {
289 // Output centre, weight, height, A0, A1 and A2
291 m_centreIndex = 0;
292 m_outPeakTableWS->addColumn("double", "centre");
293 m_outPeakTableWS->addColumn("double", "width");
294 m_outPeakTableWS->addColumn("double", "height");
295 m_outPeakTableWS->addColumn("double", "backgroundintercept");
296 m_outPeakTableWS->addColumn("double", "backgroundslope");
297 m_outPeakTableWS->addColumn("double", "A2");
298 }
299
300 m_outPeakTableWS->addColumn("double", "chi2");
301}
302
303//----------------------------------------------------------------------------------------------
310void FindPeaks::findPeaksGivenStartingPoints(const std::vector<double> &peakcentres,
311 const std::vector<double> &fitwindows) {
312 bool useWindows = (!fitwindows.empty());
313 std::size_t numPeaks = peakcentres.size();
314
315 // Loop over the spectra searching for peaks
316 m_progress = std::make_unique<Progress>(this, 0.0, 1.0, m_indexSet.size());
317
318 for (const auto spec : m_indexSet) {
319 const auto &vecX = m_dataWS->x(spec);
320
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";
324 {
325 const auto &vecY = m_dataWS->y(spec);
326 const auto &vecE = m_dataWS->e(spec);
327 const size_t numY = vecY.size();
328 size_t i_min = 1;
329 for (; i_min < numY; ++i_min) {
330 if ((vecY[i_min] != 0.) || (vecE[i_min] != 0)) {
331 --i_min; // bring it back one
332 break;
333 }
334 }
335 practical_x_min = vecX[i_min];
336
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)) {
340 ++i_max; // bring it back one
341 break;
342 }
343 }
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];
347 }
348 g_log.information() << "practical x-range = [" << practical_x_min << " -> " << practical_x_max << "]\n";
349
350 for (std::size_t ipeak = 0; ipeak < numPeaks; ipeak++) {
351 // Try to fit at this center
352 double x_center = peakcentres[ipeak];
353
354 std::stringstream infoss;
355 infoss << "Spectrum " << spec << ": Fit peak @ d = " << x_center;
356 if (useWindows) {
357 infoss << " inside fit window [" << fitwindows[2 * ipeak] << ", " << fitwindows[2 * ipeak + 1] << "]";
358 }
359 g_log.information(infoss.str());
360
361 // Check whether it is the in data range
362 if (x_center > practical_x_min && x_center < practical_x_max) {
363 if (useWindows)
364 fitPeakInWindow(m_dataWS, static_cast<int>(spec), x_center, fitwindows[2 * ipeak], fitwindows[2 * ipeak + 1]);
365 else {
366 bool hasLeftPeak = (ipeak > 0);
367 double leftpeakcentre = 0.;
368 if (hasLeftPeak)
369 leftpeakcentre = peakcentres[ipeak - 1];
370
371 bool hasRightPeak = (ipeak < numPeaks - 1);
372 double rightpeakcentre = 0.;
373 if (hasRightPeak)
374 rightpeakcentre = peakcentres[ipeak + 1];
375
376 fitPeakGivenFWHM(m_dataWS, static_cast<int>(spec), x_center, m_inputPeakFWHM, hasLeftPeak, leftpeakcentre,
377 hasRightPeak, rightpeakcentre);
378 }
379 } else {
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";
382 addNonFitRecord(spec, x_center);
383 }
384
385 } // loop through the peaks specified
386
387 m_progress->report();
388
389 } // loop over spectra
390}
391
392//----------------------------------------------------------------------------------------------
396 // At this point the data has not been smoothed yet.
397 auto smoothedData = this->calculateSecondDifference(m_dataWS);
398
399 // The optimum number of points in the smoothing, according to Mariscotti, is
400 // 0.6*fwhm
401 auto w = static_cast<int>(0.6 * m_inputPeakFWHM);
402 // w must be odd
403 if (!(w % 2))
404 ++w;
405
406 smoothData(smoothedData, w, g_z);
407
408 // Now calculate the errors on the smoothed data
409 this->calculateStandardDeviation(m_dataWS, smoothedData, w);
410
411 // Calculate n1 (Mariscotti eqn. 18)
412 const double kz = 1.22; // This kz corresponds to z=5 & w=0.6*fwhm - see Mariscotti Fig. 8
413 const int n1 = boost::math::iround(kz * m_inputPeakFWHM);
414 // Can't calculate n2 or n3 yet because they need i0
415 const int tolerance = getProperty("Tolerance");
416
417 // Loop over the spectra searching for peaks
418 m_progress = std::make_unique<Progress>(this, 0.0, 1.0, m_indexSet.size());
419
420 for (size_t k_out = 0; k_out < m_indexSet.size(); ++k_out) {
421 const size_t k = m_indexSet[k_out];
422 const auto &S = smoothedData[k_out].y();
423 const auto &F = smoothedData[k_out].e();
424
425 // This implements the flow chart given on page 320 of Mariscotti
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) {
428
429 int M = 0;
430 if (S[i] > F[i])
431 M = 1;
432 else {
433 S[i] > 0 ? M = 2 : M = 3;
434 }
435
436 if (S[i - 1] > F[i - 1]) {
437 switch (M) {
438 case 3:
439 i3 = i;
440 /* no break */
441 // intentional fall-through
442 case 2:
443 i2 = i - 1;
444 break;
445 case 1:
446 // do nothing
447 break;
448 default:
449 assert(false);
450 // should never happen
451 break;
452 }
453 } else if (S[i - 1] > 0) {
454 switch (M) {
455 case 3:
456 i3 = i;
457 break;
458 case 2:
459 // do nothing
460 break;
461 case 1:
462 i1 = i;
463 break;
464 default:
465 assert(false);
466 // should never happen
467 break;
468 }
469 } else {
470 switch (M) {
471 case 3:
472 // do nothing
473 break;
474 case 2: // fall through (i.e. same action if M = 1 or 2)
475 case 1:
476 i5 = i - 1;
477 break;
478 default:
479 assert(false);
480 // should never happen
481 break;
482 }
483 }
484
485 if (i5 && i1 && i2 && i3) // If i5 has been set then we should have the
486 // full set and can check conditions
487 {
488 i4 = i3; // Starting point for finding i4 - calculated below
489 double num = 0.0, denom = 0.0;
490 for (int j = i3; j <= i5; ++j) {
491 // Calculate i4 - it's at the minimum value of Si between i3 & i5
492 if (S[j] <= S[i4])
493 i4 = j;
494 // Calculate sums for i0 (Mariscotti eqn. 27)
495 num += j * S[j];
496 denom += S[j];
497 }
498 i0 = static_cast<int>(num / denom);
499
500 // Check we have a correctly ordered set of points. If not, reset and
501 // continue
502 if (i1 > i2 || i2 > i3 || i3 > i4 || i5 <= i4) {
503 i5 = 0;
504 continue;
505 }
506
507 // Check if conditions are fulfilled - if any are not, loop onto the
508 // next i in the spectrum
509 // Mariscotti eqn. (14)
510 if (std::abs(S[i4]) < 2 * F[i4]) {
511 i5 = 0;
512 continue;
513 }
514 // Mariscotti eqn. (19)
515 if (abs(i5 - i3 + 1 - n1) > tolerance) {
516 i5 = 0;
517 continue;
518 }
519 // Calculate n2 (Mariscotti eqn. 20)
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)));
522 if (n2b > n2)
523 n2 = n2b;
524 // Mariscotti eqn. (21)
525 const int testVal = n2 ? n2 : 1;
526 if (i3 - i2 - 1 > testVal) {
527 i5 = 0;
528 continue;
529 }
530 // Calculate n3 (Mariscotti eqn. 22)
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]))));
533 if (n3b < n3)
534 n3 = n3b;
535 // Mariscotti eqn. (23)
536 if (i2 - i1 + 1 < n3) {
537 i5 = 0;
538 continue;
539 }
540
541 // If we get to here then we've identified a peak
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';
544
545 // Use i0, i2 and i4 to find out i_min and i_max, i0: right, i2: left,
546 // i4: centre
547 auto wssize = static_cast<int>(m_dataWS->x(k).size());
548
549 int iwidth = i0 - i2;
550 if (iwidth <= 0)
551 iwidth = 1;
552
553 int i_min = 1;
554 if (i4 > 5 * iwidth)
555 i_min = i4 - 5 * iwidth;
556
557 int i_max = i4 + 5 * iwidth;
558 if (i_max >= wssize)
559 i_max = wssize - 1;
560
561 this->fitSinglePeak(m_dataWS, static_cast<int>(k), i_min, i_max, i4);
562
563 // reset and go searching for the next peak
564 i1 = 0, i2 = 0, i3 = 0, i4 = 0, i5 = 0;
565 }
566
567 } // loop through a single spectrum
568
569 m_progress->report();
570
571 } // loop over spectra
572}
573
574//----------------------------------------------------------------------------------------------
584 std::vector<Histogram> diffed;
585
586 // Loop over spectra
587 for (const auto i : m_indexSet) {
588 diffed.emplace_back(input->histogram(i));
589 diffed.back().mutableY() = 0.0;
590 diffed.back().mutableE() = 0.0;
591
592 const auto &Y = input->y(i);
593 auto &S = diffed.back().mutableY();
594 // Go through each spectrum calculating the second difference at each point
595 // First and last points in each spectrum left as zero (you'd never be able
596 // to find peaks that close to the edge anyway)
597 for (size_t j = 1; j < S.size() - 1; ++j) {
598 S[j] = Y[j - 1] - 2 * Y[j] + Y[j + 1];
599 }
600 }
601
602 return diffed;
603}
604
605//----------------------------------------------------------------------------------------------
612void FindPeaks::smoothData(std::vector<Histogram> &histograms, const int w, const int g_z) {
613 for (auto &histogram : histograms)
614 for (int i = 0; i < g_z; ++i)
615 histogram = smooth(histogram, w);
616}
617
618//----------------------------------------------------------------------------------------------
629 std::vector<HistogramData::Histogram> &smoothed, const int &w) {
630 // Guard against anyone changing the value of z, which would mean different
631 // phi values were needed (see Marriscotti p.312)
632 static_assert(g_z == 5, "Value of z has changed!");
633 // Have to adjust for fact that I normalise Si (unlike the paper)
634 const auto factor = static_cast<int>(std::pow(static_cast<double>(w), g_z));
635
636 const double constant = sqrt(static_cast<double>(this->computePhi(w))) / factor;
637
638 for (size_t i = 0; i < m_indexSet.size(); ++i) {
639 size_t i_in = m_indexSet[i];
640 smoothed[i].mutableE() = input->e(i_in) * constant;
641 }
642}
643
644//----------------------------------------------------------------------------------------------
653long long FindPeaks::computePhi(const int &w) const {
654 const int m = (w - 1) / 2;
655 int zz = 0;
656 int max_index_prev = 1;
657 int n_el_prev = 3;
658 std::vector<long long> previous(n_el_prev);
659 previous[0] = 1;
660 previous[1] = -2;
661 previous[2] = 1;
662
663 // Can't happen at present
664 if (g_z == 0)
665 return std::accumulate(previous.begin(), previous.end(), static_cast<long long>(0),
667
668 std::vector<long long> next;
669 // Calculate the Cij iteratively.
670 do {
671 zz++;
672 int max_index = zz * m + 1;
673 int n_el = 2 * max_index + 1;
674 next.resize(n_el);
675 std::fill(next.begin(), next.end(), 0);
676 for (int i = 0; i < n_el; ++i) {
677 int delta = -max_index + i;
678 for (int l = delta - m; l <= delta + m; l++) {
679 int index = l + max_index_prev;
680 if (index >= 0 && index < n_el_prev)
681 next[i] += previous[index];
682 }
683 }
684 previous.resize(n_el);
685 std::copy(next.begin(), next.end(), previous.begin());
686 max_index_prev = max_index;
687 n_el_prev = n_el;
688 } while (zz != g_z);
689
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";
693 return retval;
694}
695
696//----------------------------------------------------------------------------------------------
702int FindPeaks::getIndex(const HistogramX &vecX, double x) {
703 int index;
704
705 if (x <= vecX.front()) {
706 // Left or equal to lower boundary
707 index = 0;
708 } else if (x >= vecX.back()) {
709 // Right or equal to upper boundary
710 index = static_cast<int>(vecX.size()) - 1;
711 } else {
712 // within the range
713 index = static_cast<int>(std::lower_bound(vecX.begin(), vecX.end(), x) - vecX.begin());
714
715 // check lower boundary
716 if (index == 0) {
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.";
720 g_log.warning(errss.str());
721 throw std::runtime_error(errss.str());
722 } else if (x < vecX[index - 1] || x > vecX[index]) {
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.";
726 g_log.warning(errss.str());
727 throw std::runtime_error(errss.str());
728 }
729
730 // Find the index of the nearest value to return
731 if (x - vecX[index - 1] < vecX[index] - x)
732 --index;
733 }
734
735 return index;
736}
737
738//----------------------------------------------------------------------------------------------
755void FindPeaks::fitPeakGivenFWHM(const API::MatrixWorkspace_sptr &input, const int wsIndex, const double center_guess,
756 const int fitWidth, const bool hasleftpeak, const double leftpeakcentre,
757 const bool hasrightpeak, const double rightpeakcentre) {
758 // The X axis you are looking at
759 const auto &vecX = input->x(wsIndex);
760 const auto &vecY = input->y(wsIndex);
761
762 // Find i_center - the index of the center - The guess is within the X axis?
763 int i_centre = this->getIndex(vecX, center_guess);
764
765 // Set up lower fit boundary
766 int i_min = i_centre - 5 * fitWidth;
767 if (i_min < 1)
768 i_min = 1;
769
770 if (hasleftpeak) {
771 // Use 2/3 distance as the seperation for right peak
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);
776 }
777
778 // Set up upper boundary
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;
782
783 if (hasrightpeak) {
784 // Use 2/3 distance as the separation for right peak
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);
789 }
790
791 // Check
792 if (i_max - i_min <= 0)
793 throw std::runtime_error("Impossible to i_min >= i_max.");
794
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];
798 g_log.information(outss.str());
799
800 fitSinglePeak(input, wsIndex, i_min, i_max, i_centre);
801}
802
803//----------------------------------------------------------------------------------------------
813void FindPeaks::fitPeakInWindow(const API::MatrixWorkspace_sptr &input, const int wsIndex, const double centre_guess,
814 const double xmin, const double xmax) {
815 // Check
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. ");
820 addNonFitRecord(wsIndex, centre_guess);
821 return;
822 }
823
824 // The X axis you are looking at
825 const auto &vecX = input->x(wsIndex);
826
827 // The centre index
828 int i_centre = this->getIndex(vecX, centre_guess);
829
830 // The left index
831 int i_min = getIndex(vecX, xmin);
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";
835 addNonFitRecord(wsIndex, centre_guess);
836 return;
837 }
838
839 // The right index
840 int i_max = getIndex(vecX, xmax);
841 if (i_max < i_centre) {
842 g_log.warning() << "Input peak centre @ " << centre_guess << " is out side of maximum x = " << xmax << "\n";
843 addNonFitRecord(wsIndex, centre_guess);
844 return;
845 }
846
847 // finally do the actual fit
848 fitSinglePeak(input, wsIndex, i_min, i_max, i_centre);
849}
850
851//----------------------------------------------------------------------------------------------
855void FindPeaks::fitSinglePeak(const API::MatrixWorkspace_sptr &input, const int spectrum, const int i_min,
856 const int i_max, const int i_centre) {
857 const auto &vecX = input->x(spectrum);
858 const auto &vecY = input->y(spectrum);
859
860 // Exclude peak with peak counts
861 bool hasHighCounts = false;
862 for (int i = i_min; i <= i_max; ++i)
863 if (vecY[i] > m_leastMaxObsY) {
864 hasHighCounts = true;
865 break;
866 }
867 if (!hasHighCounts) {
868 std::stringstream ess;
869 ess << "Peak supposed at " << vecY[i_centre] << " does not have enough counts as " << m_leastMaxObsY;
870 g_log.debug(ess.str());
871 addNonFitRecord(spectrum, vecY[i_centre]);
872 return;
873 }
874
875 {
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 << "). ";
879 g_log.information(outss.str());
880 }
881
882 // Estimate background
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) {
887 // Estimate background roughly for a failed case
888 estimateBackground(vecX, vecY, i_min, i_max, vecbkgdparvalue);
889 }
890
891 for (size_t i = 0; i < vecbkgdparvalue.size(); ++i)
892 if (i < m_bkgdOrder)
893 m_backgroundFunction->setParameter(i, vecbkgdparvalue[i]);
894
895 // Estimate peak parameters
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()) {
902 // Unable to estimate peak
903 i_obscentre = i_centre;
904 est_fwhm = 1.;
905 est_height = 1.;
906 g_log.warning(errmsg);
907 }
908
909 // Set peak parameters to
910 if (m_useObsCentre)
911 m_peakFunction->setCentre(vecX[i_obscentre]);
912 else
913 m_peakFunction->setCentre(vecX[i_centre]);
914 m_peakFunction->setHeight(est_height);
915 m_peakFunction->setFwhm(est_fwhm);
916
917 if (usefpdresult < 0) {
918 // Estimate peak range based on estimated linear background and peak
919 // parameter estimated from observation
920 if (!m_useObsCentre)
921 i_obscentre = i_centre;
922 estimatePeakRange(vecX, i_obscentre, i_min, i_max, est_leftfwhm, est_rightfwhm, vecpeakrange);
923 }
924
925 //-------------------------------------------------------------------------
926 // Fit Peak
927 //-------------------------------------------------------------------------
928 std::vector<double> fitwindow(2);
929 fitwindow[0] = vecX[i_min];
930 fitwindow[1] = vecX[i_max];
931
932 double costfuncvalue = callFitPeak(input, spectrum, m_peakFunction, m_backgroundFunction, fitwindow, vecpeakrange,
934
935 bool fitsuccess = false;
936 if (costfuncvalue < DBL_MAX && costfuncvalue >= 0. && m_peakFunction->height() > m_minHeight) {
937 fitsuccess = true;
938 }
939 if (fitsuccess && m_usePeakPositionTolerance) {
940 fitsuccess = (fabs(m_peakFunction->centre() - vecX[i_centre]) < m_peakPositionTolerance);
941 }
942
943 //-------------------------------------------------------------------------
944 // Process Fit result
945 //-------------------------------------------------------------------------
946 // Update output
947 if (fitsuccess)
949 else
950 addNonFitRecord(spectrum, m_peakFunction->centre());
951}
952
953//----------------------------------------------------------------------------------------------
957int FindPeaks::findPeakBackground(const MatrixWorkspace_sptr &input, int spectrum, size_t i_min, size_t i_max,
958 std::vector<double> &vecBkgdParamValues, std::vector<double> &vecpeakrange) {
959 const auto &vecX = input->x(spectrum);
960
961 // Call FindPeakBackground
962 auto estimate = createChildAlgorithm("FindPeakBackground");
963 estimate->setLoggingOffset(1);
964 estimate->setProperty("InputWorkspace", input);
965 estimate->setProperty("WorkspaceIndex", spectrum);
966 // estimate->setProperty("SigmaConstant", 1.0);
967 std::vector<double> fwvec;
968 fwvec.emplace_back(vecX[i_min]);
969 fwvec.emplace_back(vecX[i_max]);
970 estimate->setProperty("BackgroundType", m_backgroundType);
971 estimate->setProperty("FitWindow", fwvec);
972 estimate->executeAsChildAlg();
973 // Get back the result
974 Mantid::API::ITableWorkspace_sptr peaklisttablews = estimate->getProperty("OutputWorkspace");
975
976 // Determine whether to use FindPeakBackground's result.
977 int fitresult = -1;
978 if (peaklisttablews->columnCount() < 7)
979 throw std::runtime_error("No 7th column for use FindPeakBackground result or not. ");
980
981 if (peaklisttablews->rowCount() > 0) {
987 const int hiddenFitresult = peaklisttablews->Int(0, 6);
988 g_log.information() << "fitresult=" << hiddenFitresult << "\n";
989 }
990
991 // Local check whether FindPeakBackground gives a reasonable value
992 vecpeakrange.resize(2);
993 if (fitresult > 0) {
994 // Use FitPeakBackgroud's result
995 size_t i_peakmin, i_peakmax;
996 i_peakmin = peaklisttablews->Int(0, 1);
997 i_peakmax = peaklisttablews->Int(0, 2);
998
999 g_log.information() << "FindPeakBackground successful. "
1000 << "iMin = " << i_min << ", iPeakMin = " << i_peakmin << ", iPeakMax = " << i_peakmax
1001 << ", iMax = " << i_max << "\n";
1002
1003 if (i_peakmin < i_peakmax && i_peakmin > i_min + 2 && i_peakmax < i_max - 2) {
1004 // FIXME - It is assumed that there are 3 background parameters set to
1005 // FindPeaksBackground
1006 double bg0, bg1, bg2;
1007 bg0 = peaklisttablews->Double(0, 3);
1008 bg1 = peaklisttablews->Double(0, 4);
1009 bg2 = peaklisttablews->Double(0, 5);
1010
1011 // Set output
1012 vecBkgdParamValues.resize(3, 0.);
1013 vecBkgdParamValues[0] = bg0;
1014 vecBkgdParamValues[1] = bg1;
1015 vecBkgdParamValues[2] = bg2;
1016
1017 g_log.information() << "Background parameters (from FindPeakBackground) A0=" << bg0 << ", A1=" << bg1
1018 << ", A2=" << bg2 << "\n";
1019
1020 vecpeakrange[0] = vecX[i_peakmin];
1021 vecpeakrange[1] = vecX[i_peakmax];
1022 } else {
1023 // Do manual estimation again
1024 g_log.debug("FindPeakBackground result is ignored due to wrong in peak range. ");
1025 }
1026 } else {
1027 g_log.information("Failed to get background estimation\n");
1028 }
1029
1030 std::stringstream outx;
1031 outx << "FindPeakBackground Result: Given window (" << vecX[i_min] << ", " << vecX[i_max]
1032 << "); Determine peak range: (" << vecpeakrange[0] << ", " << vecpeakrange[1] << "). ";
1033 g_log.information(outx.str());
1034
1035 return fitresult;
1036}
1037
1038//----------------------------------------------------------------------------------------------
1055std::string FindPeaks::estimatePeakParameters(const HistogramX &vecX, const HistogramY &vecY, size_t i_min,
1056 size_t i_max, const std::vector<double> &vecbkgdparvalues,
1057 size_t &iobscentre, double &height, double &fwhm, double &leftfwhm,
1058 double &rightfwhm) {
1059 // Search for maximum considering background
1060 const double bg0 = vecbkgdparvalues[0];
1061 double bg1 = 0;
1062 double bg2 = 0;
1063 if (vecbkgdparvalues.size() >= 2) {
1064 bg1 = vecbkgdparvalues[1];
1065 if (vecbkgdparvalues.size() >= 3)
1066 bg2 = vecbkgdparvalues[2];
1067 }
1068
1069 // Starting value
1070 iobscentre = i_min;
1071 const double tmpx = vecX[i_min];
1072 height = vecY[i_min] - (bg0 + bg1 * tmpx + bg2 * tmpx * tmpx);
1073 double lowest = height;
1074
1075 // Extreme case
1076 if (i_max == vecY.size())
1077 i_max = i_max - 1;
1078
1079 // Searching
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);
1083
1084 if (tmpheight > height) {
1085 iobscentre = i;
1086 height = tmpheight;
1087 } else if (tmpheight < lowest) {
1088 lowest = tmpheight;
1089 }
1090 }
1091
1092 // Summarize on peak centre
1093 double obscentre = vecX[iobscentre];
1094 double drop = height - lowest;
1095 if (drop == 0) {
1096 // Flat spectrum. No peak parameter can be estimated.
1097 // FIXME - should have a second notice such that there will be no more
1098 // fitting.
1099 return "Flat spectrum";
1100 } else if (height <= m_minHeight) {
1101 // The peak is not high enough!
1102 // FIXME - should have a second notice such that there will be no more
1103 // fitting.
1104 return "Fluctuation is less than minimum allowed value.";
1105 }
1106
1107 // If maximum point is on the edge 2 points, return false. One side of peak
1108 // must have at least 3 points
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 << ").";
1113 return dbss.str();
1114 }
1115
1116 // Search for half-maximum: no need to very precise
1117
1118 // Slope at the left side of peak.
1119 leftfwhm = -1.;
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) {
1124 // Ideal case
1125 // FIXME - Need a linear interpolation on it!
1126 leftfwhm = obscentre - 0.5 * (vecX[i] + vecX[i + 1]);
1127 break;
1128 }
1129 }
1130
1131 // Slope at the right side of peak
1132 rightfwhm = -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;
1138 break;
1139 }
1140 }
1141
1142 if (leftfwhm <= 0 || rightfwhm <= 0) {
1143 std::stringstream errmsg;
1144 errmsg << "Estimate peak parameters error (FWHM cannot be zero): Input "
1145 "data size = "
1146 << vecX.size() << ", Xmin = " << vecX[i_min] << "(" << i_min << "), Xmax = " << vecX[i_max] << "(" << i_max
1147 << "); "
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();
1152 }
1153
1154 fwhm = leftfwhm + rightfwhm;
1155 if (fwhm < 1.e-200) // very narrow peak
1156 {
1157 std::stringstream errmsg;
1158 errmsg << "Estimate peak parameters error (FWHM cannot be zero): Input "
1159 "data size = "
1160 << vecX.size() << ", Xmin = " << vecX[i_min] << "(" << i_min << "), Xmax = " << vecX[i_max] << "(" << i_max
1161 << "); "
1162 << "Estimated peak centre @ " << vecX[iobscentre] << "(" << iobscentre << ") with height = " << height
1163 << "; Lowest Y value = " << lowest << "; Output error: . fwhm = " << fwhm << ".";
1164 return errmsg.str();
1165 }
1166
1167 g_log.information() << "Estimated peak parameters: Centre = " << obscentre << ", Height = " << height
1168 << ", FWHM = " << fwhm << " = " << leftfwhm << " + " << rightfwhm << ".\n";
1169
1170 return std::string();
1171}
1172
1173//----------------------------------------------------------------------------------------------
1184void FindPeaks::estimateBackground(const HistogramX &X, const HistogramY &Y, const size_t i_min, const size_t i_max,
1185 std::vector<double> &vecbkgdparvalues) {
1186 // Validate input
1187 if (i_min >= i_max)
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.);
1192
1193 // FIXME - THIS IS A MAGIC NUMBER
1194 const size_t MAGICNUMBER = 12;
1195 size_t numavg;
1196 if (i_max - i_min > MAGICNUMBER)
1197 numavg = 3;
1198 else
1199 numavg = 1;
1200
1201 // Get (x0, y0) and (xf, yf)
1202 double x0 = 0.0;
1203 double y0 = 0.0;
1204 double xf = 0.0;
1205 double yf = 0.0;
1206
1207 for (size_t i = 0; i < numavg; ++i) {
1208 x0 += X[i_min + i];
1209 y0 += Y[i_min + i];
1210
1211 xf += X[i_max - i];
1212 // X has one more value than Y
1213 yf += Y[i_max - i - 1];
1214 }
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);
1219
1220 // Esitmate
1221 vecbkgdparvalues[2] = 0.;
1222 if (m_bkgdOrder >= 1) {
1223 // linear background
1224 vecbkgdparvalues[1] = (y0 - yf) / (x0 - xf);
1225 vecbkgdparvalues[0] = (xf * y0 - x0 * yf) / (xf - x0);
1226 } else {
1227 // flat background
1228 vecbkgdparvalues[1] = 0.;
1229 vecbkgdparvalues[0] = 0.5 * (y0 + yf);
1230 }
1231}
1232
1233//----------------------------------------------------------------------------------------------
1237void FindPeaks::estimatePeakRange(const HistogramX &vecX, size_t i_centre, size_t i_min, size_t i_max,
1238 const double &leftfwhm, const double &rightfwhm, std::vector<double> &vecpeakrange) {
1239 // Check
1240 if (vecpeakrange.size() < 2)
1241 vecpeakrange.resize(2, 0.);
1242
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. ");
1245
1246 // Search peak left by using 6 * half of FWHM
1247 double peakleftbound = vecX[i_centre] - 6. * leftfwhm;
1248 double peakrightbound = vecX[i_centre] + 6. * rightfwhm;
1249
1250 // Deal with case the peak boundary is too close to fit window
1251 auto ipeakleft = static_cast<size_t>(getIndex(vecX, peakleftbound));
1252 if (ipeakleft <= i_min) {
1253 size_t numbkgdpts = (i_centre - i_min) / 6;
1254 // FIXME - 3 is a magic number
1255 if (numbkgdpts < 3)
1256 numbkgdpts = 3;
1257 ipeakleft = i_min + numbkgdpts;
1258 if (ipeakleft >= i_centre)
1259 ipeakleft = i_min + 1;
1260
1261 peakleftbound = vecX[ipeakleft];
1262 }
1263
1264 auto ipeakright = static_cast<size_t>(getIndex(vecX, peakrightbound));
1265 if (ipeakright >= i_max) {
1266 size_t numbkgdpts = (i_max - i_centre) / 6;
1267 // FIXME - 3 is a magic number
1268 if (numbkgdpts < 3)
1269 numbkgdpts = 3;
1270 ipeakright = i_max - numbkgdpts;
1271 if (ipeakright <= i_centre)
1272 ipeakright = i_max - 1;
1273
1274 peakrightbound = vecX[ipeakright];
1275 }
1276
1277 // Set result to output vector
1278 vecpeakrange[0] = peakleftbound;
1279 vecpeakrange[1] = peakrightbound;
1280}
1281
1282//----------------------------------------------------------------------------------------------
1291void FindPeaks::addInfoRow(const size_t spectrum, const API::IPeakFunction_const_sptr &peakfunction,
1292 const API::IBackgroundFunction_sptr &bkgdfunction, const bool isoutputraw,
1293 const double mincost) {
1294 // Check input validity
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. ");
1299
1300 // Add fitted parameters to output table workspace
1301 API::TableRow t = m_outPeakTableWS->appendRow();
1302
1303 // spectrum
1304 t << static_cast<int>(spectrum);
1305
1306 // peak and background function parameters
1307 if (isoutputraw) {
1308 // Output of raw peak parameters
1309 size_t nparams = peakfunction->nParams();
1310 size_t nparamsb = bkgdfunction->nParams();
1311
1312 size_t numcols = m_outPeakTableWS->columnCount();
1313 if (nparams + nparamsb + 2 != numcols) {
1314 throw std::runtime_error("Error 1307 number of columns do not matches");
1315 }
1316
1317 for (size_t i = 0; i < nparams; ++i) {
1318 t << peakfunction->getParameter(i);
1319 }
1320 for (size_t i = 0; i < nparamsb; ++i) {
1321 t << bkgdfunction->getParameter(i);
1322 }
1323 } else {
1324 // Output of effective peak parameters
1325 // Set up parameters to peak function and get effective value
1326 double peakcentre = peakfunction->centre();
1327 double fwhm = peakfunction->fwhm();
1328 double height = peakfunction->height();
1329
1330 t << peakcentre << fwhm << height;
1331
1332 // Set up parameters to background function
1333
1334 // FIXME - Use Polynomial for all 3 background types.
1335 double a0 = bkgdfunction->getParameter("A0");
1336 double a1 = 0.;
1337 if (bkgdfunction->name() != "FlatBackground")
1338 a1 = bkgdfunction->getParameter("A1");
1339 double a2 = 0;
1340 if (bkgdfunction->name() != "LinearBackground" && bkgdfunction->name() != "FlatBackground")
1341 a2 = bkgdfunction->getParameter("A2");
1342
1343 t << a0 << a1 << a2;
1344
1345 g_log.debug() << "Peak parameters found: cen=" << peakcentre << " fwhm=" << fwhm << " height=" << height
1346 << " a0=" << a0 << " a1=" << a1 << " a2=" << a2;
1347 }
1348 g_log.debug() << " chsq=" << mincost << "\n";
1349 // Minimum cost function value
1350 t << mincost;
1351}
1352
1353//----------------------------------------------------------------------------------------------
1358void FindPeaks::addNonFitRecord(const size_t spectrum, const double centre) {
1359 // Add a new row
1360 API::TableRow t = m_outPeakTableWS->appendRow();
1361
1362 g_log.information() << "Failed to fit peak at " << centre << "\n";
1363 // 1st column
1364 t << static_cast<int>(spectrum);
1365
1366 // Parameters
1367 for (std::size_t i = 0; i < m_numTableParams; i++) {
1368 if (i == m_centreIndex)
1369 t << centre;
1370 else
1371 t << 0.;
1372 }
1373
1374 // HUGE chi-square
1375 t << DBL_MAX;
1376}
1377
1378//----------------------------------------------------------------------------------------------
1382 // Setup the background
1383 // FIXME (No In This Ticket) Need to have a uniformed routine to name
1384 // background function
1385 std::string backgroundposix;
1386 if (m_backgroundType != "Quadratic") {
1387 // FlatBackground, LinearBackground, Quadratic
1388 backgroundposix = "Background";
1389 }
1390 m_backgroundFunction = std::dynamic_pointer_cast<IBackgroundFunction>(
1391 API::FunctionFactory::Instance().createFunction(m_backgroundType + backgroundposix));
1392 g_log.information() << "Background function (" << m_backgroundFunction->name() << ") has been created. "
1393 << "\n";
1394
1395 m_bkgdParameterNames = m_backgroundFunction->getParameterNames();
1396 // FIXME - Need to add method nOrder to background function;
1397 m_bkgdOrder = m_backgroundFunction->nParams() - 1;
1398
1399 // Set up peak function
1401 std::dynamic_pointer_cast<IPeakFunction>(API::FunctionFactory::Instance().createFunction(m_peakFuncType));
1402 m_peakParameterNames = m_peakFunction->getParameterNames();
1403}
1404
1405//----------------------------------------------------------------------------------------------
1408double FindPeaks::callFitPeak(const MatrixWorkspace_sptr &dataws, int wsindex,
1409 const API::IPeakFunction_sptr &peakfunction,
1410 const API::IBackgroundFunction_sptr &backgroundfunction,
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;
1415 g_log.information(dbss.str());
1416
1417 double userFWHM = m_peakFunction->fwhm();
1418 bool fitwithsteppedfwhm = (guessedFWHMStep > 0);
1419
1420 FitOneSinglePeak fitpeak;
1421 fitpeak.setChild(true);
1422 fitpeak.setWorskpace(dataws, wsindex);
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]);
1428
1429 if (estBackResult == 1) {
1430 g_log.information("simpleFit");
1431 fitpeak.simpleFit();
1432 } else if (m_highBackground) {
1433 g_log.information("highBkgdFit");
1434 fitpeak.highBkgdFit();
1435 } else {
1436 g_log.information("simpleFit");
1437 fitpeak.simpleFit();
1438 }
1439
1440 double costfuncvalue = fitpeak.getFitCostFunctionValue();
1441 std::string dbinfo = fitpeak.getDebugMessage();
1442 g_log.information(dbinfo);
1443
1444 return costfuncvalue;
1445}
1446
1447//----------------------------------------------------------------------------------------------
1453 size_t numpeakpars = m_peakFunction->nParams();
1454 std::vector<double> vec_value(numpeakpars);
1455 for (size_t i = 0; i < numpeakpars; ++i) {
1456 double parvalue = m_peakFunction->getParameter(i);
1457 vec_value[i] = parvalue;
1458 }
1459
1460 return vec_value;
1461}
1462
1463} // namespace Mantid::Algorithms
1464
1465// 0.5044, 0.5191, 0.535, 0.5526, 0.5936, 0.6178, 0.6453, 0.6768, 0.7134,
1466// 0.7566, 0.8089, 0.8737, 0.9571, 1.0701, 1.2356, 1.5133, 2.1401
#define DECLARE_ALGORITHM(classname)
Definition: Algorithm.h:576
CostFunctions::CostFuncFitting & m_costFunction
The cost function.
const double MAGICNUMBER
Definition: FitPeak.cpp:47
double height
Definition: GetAllEi.cpp:155
std::map< DeltaEMode::Type, std::string > index
Definition: DeltaEMode.cpp:19
#define fabs(x)
Definition: Matrix.cpp:22
double tolerance
void declareProperty(std::unique_ptr< Kernel::Property > p, const std::string &doc="") override
Add a property to the list of managed properties.
Definition: Algorithm.cpp:1913
std::string getPropertyValue(const std::string &name) const override
Get the value of a property as a string.
Definition: Algorithm.cpp:2026
TypedValue getProperty(const std::string &name) const override
Get the value of a property.
Definition: Algorithm.cpp:2076
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.
Definition: Algorithm.cpp:842
Kernel::Logger & g_log
Definition: Algorithm.h:451
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.
Definition: Algorithm.cpp:166
An interface to a peak function, which extend the interface of IFunctionWithLocation by adding method...
Definition: IPeakFunction.h:51
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.
Definition: TableRow.h:39
A property class for workspaces.
This algorithm searches for peaks in a dataset.
Definition: FindPeaks.h:50
static const int g_z
The number of smoothing iterations.
Definition: FindPeaks.h:156
Indexing::SpectrumIndexSet m_indexSet
list of workspace indicies to check
Definition: FindPeaks.h:166
API::ITableWorkspace_sptr m_outPeakTableWS
Storage of the peak data.
Definition: FindPeaks.h:159
void init() override
Initialize and declare properties.
Definition: FindPeaks.cpp:60
API::IBackgroundFunction_sptr m_backgroundFunction
Definition: FindPeaks.h:182
bool m_useObsCentre
Start values.
Definition: FindPeaks.h:205
std::vector< double > m_vecPeakCentre
Definition: FindPeaks.h:178
void generateOutputPeakParameterTable()
Generate a table workspace for output peak parameters.
Definition: FindPeaks.cpp:263
std::unique_ptr< API::Progress > m_progress
Progress reporting.
Definition: FindPeaks.h:161
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.
Definition: FindPeaks.cpp:755
int getIndex(const HistogramData::HistogramX &vecX, double x)
needed by FindPeaksBackground
Definition: FindPeaks.cpp:702
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.
Definition: FindPeaks.cpp:957
void processAlgorithmProperties()
Process algorithm's properties.
Definition: FindPeaks.cpp:191
void createFunctions()
Create peak and background functions.
Definition: FindPeaks.cpp:1381
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)
Definition: FindPeaks.cpp:813
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.
Definition: FindPeaks.cpp:1237
API::MatrixWorkspace_sptr m_dataWS
workspace to check for peaks
Definition: FindPeaks.h:164
double m_leastMaxObsY
Minimum value of peak's observed maximum Y value.
Definition: FindPeaks.h:202
double m_minHeight
Minimum peak height.
Definition: FindPeaks.h:200
API::IPeakFunction_sptr m_peakFunction
Definition: FindPeaks.h:183
std::size_t m_numTableParams
parameters or effective (centre, width, height)
Definition: FindPeaks.h:172
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'.
Definition: FindPeaks.cpp:1408
std::vector< double > m_vecFitWindows
Definition: FindPeaks.h:179
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.
Definition: FindPeaks.cpp:1055
void calculateStandardDeviation(const API::MatrixWorkspace_const_sptr &input, std::vector< HistogramData::Histogram > &smoothed, const int &w)
Calculates the statistical error on the smoothed data.
Definition: FindPeaks.cpp:628
long long computePhi(const int &w) const
Calculates the coefficient phi which goes into the calculation of the error on the smoothed data Uses...
Definition: FindPeaks.cpp:653
void exec() override
Execute the findPeaks algorithm.
Definition: FindPeaks.cpp:156
bool m_highBackground
flag for find relatively weak peak in high
Definition: FindPeaks.h:168
void findPeaksGivenStartingPoints(const std::vector< double > &peakcentres, const std::vector< double > &fitwindows)
Find peaks according to given peak positions.
Definition: FindPeaks.cpp:310
void findPeaksUsingMariscotti()
Find peaks by searching peak position using Mariscotti algorithm.
Definition: FindPeaks.cpp:395
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.
Definition: FindPeaks.cpp:855
std::vector< std::string > m_bkgdParameterNames
Definition: FindPeaks.h:151
std::vector< std::string > m_peakParameterNames
Definition: FindPeaks.h:150
void addNonFitRecord(const size_t spectrum, const double centre)
Add the fit record (failure) to output workspace.
Definition: FindPeaks.cpp:1358
std::vector< double > getStartingPeakValues()
Get the peak parameter values from m_peakFunction and output to a list in the same order of m_peakPar...
Definition: FindPeaks.cpp:1452
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.
Definition: FindPeaks.cpp:1184
std::vector< HistogramData::Histogram > calculateSecondDifference(const API::MatrixWorkspace_const_sptr &input)
Methods searving for findPeaksUsingMariscotti()
Definition: FindPeaks.cpp:583
void smoothData(std::vector< HistogramData::Histogram > &histograms, const int w, const int g_z)
Smooth data for Mariscotti.
Definition: FindPeaks.cpp:612
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.
Definition: FindPeaks.cpp:1291
bool m_rawPeaksTable
background
Definition: FindPeaks.h:170
int m_inputPeakFWHM
holder for the requested peak FWHM
Definition: FindPeaks.h:165
FitOneSinglePeak: a class to perform peak fitting on a single peak.
Definition: FitPeak.h:30
void setPeakRange(double xpeakleft, double xpeakright)
Set peak range.
Definition: FitPeak.cpp:112
void setupGuessedFWHM(double usrwidth, int minfwhm, int maxfwhm, int stepsize, bool fitwithsteppedfwhm)
Set peak width to guess.
Definition: FitPeak.cpp:156
std::string getDebugMessage()
Get debug message.
Definition: FitPeak.cpp:276
void setFunctions(const API::IPeakFunction_sptr &peakfunc, const API::IBackgroundFunction_sptr &bkgdfunc)
Set functions.
Definition: FitPeak.cpp:82
void highBkgdFit()
Fit peak first considering high background.
Definition: FitPeak.cpp:459
bool simpleFit()
Fit peak and background together.
Definition: FitPeak.cpp:281
void setFittingMethod(std::string minimizer, const std::string &costfunction)
Set fitting method.
Definition: FitPeak.cpp:129
void setFitWindow(double leftwindow, double rightwindow)
Set fit range.
Definition: FitPeak.cpp:93
void setWorskpace(const API::MatrixWorkspace_sptr &dataws, size_t wsindex)
Set workspaces.
Definition: FitPeak.cpp:65
double getFitCostFunctionValue()
Get cost function value from fitting.
Definition: FitPeak.cpp:974
Support for a property that holds an array of values.
Definition: ArrayProperty.h:28
Exception for index errors.
Definition: Exception.h:284
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...
Definition: ListValidator.h:29
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
void information(const std::string &msg)
Logs at information level.
Definition: Logger.cpp:105
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.
Definition: IValidator.h:26
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.
Definition: EmptyValues.h:25
constexpr double EMPTY_DBL() noexcept
Returns what we consider an "empty" double within a property.
Definition: EmptyValues.h:43
@ Input
An input workspace.
Definition: Property.h:53
@ Output
An output workspace.
Definition: Property.h:54
Functor to accumulate a sum of squares.
Definition: VectorHelper.h:170