Mantid
Loading...
Searching...
No Matches
PeakParameterHelper.cpp
Go to the documentation of this file.
1// Mantid Repository : https://github.com/mantidproject/mantid
2//
3// Copyright © 2021 ISIS Rutherford Appleton Laboratory UKRI,
4// NScD Oak Ridge National Laboratory, European Spallation Source,
5// Institut Laue - Langevin & CSNS, Institute of High Energy Physics, CAS
6// SPDX - License - Identifier: GPL - 3.0 +
7
9
10using namespace Mantid;
11using namespace Mantid::API;
12using Mantid::HistogramData::Histogram;
13
15
19template <typename vector_like> size_t findXIndex(const vector_like &vecx, const double x, const size_t startindex) {
20 size_t index;
21 if (x <= vecx.front()) {
22 index = 0;
23 } else if (x >= vecx.back()) {
24 index = vecx.size() - 1;
25 } else {
26 const auto fiter = std::lower_bound(vecx.cbegin() + startindex, vecx.cend(), x);
27 if (fiter == vecx.cend())
28 throw std::runtime_error("It seems impossible to have this value. ");
29
30 index = static_cast<size_t>(fiter - vecx.cbegin());
31 if (index > 0 && (x - vecx[index - 1] < vecx[index] - x))
32 --index;
33 }
34
35 return index;
36}
37
38//----------------------------------------------------------------------------------------------
51int observePeakCenter(const Histogram &histogram, const FunctionValues &bkgd_values, size_t start_index,
52 size_t stop_index, double &peak_center, size_t &peak_center_index, double &peak_height) {
53 const auto &vector_x = histogram.points();
54
55 // find the original starting point
56 auto peak_center_iter = std::lower_bound(vector_x.begin() + start_index, vector_x.begin() + stop_index, peak_center);
57 if (peak_center_iter == vector_x.begin() + stop_index)
58 return OUTOFBOUND; // suggested center is not in the window
59 peak_center_index = static_cast<size_t>(peak_center_iter - vector_x.begin());
60
61 // initialize the search to that in case something goes wrong below
62 const auto &vector_y = histogram.y();
63 peak_center = *peak_center_iter; // set the value in case something goes wrong below
64 peak_height = vector_y[peak_center_index] - bkgd_values.getCalculated(peak_center_index - start_index);
65
66 // assume that the actual peak is within 40% (in index number)
67 // of the window size of the suggested peak. This assumes that
68 // the minimum search size is 5 bins (arbitrary).
69 const size_t windowSize = stop_index - start_index;
70 const size_t searchBox = std::max(static_cast<size_t>(.3 * static_cast<double>(windowSize)), static_cast<size_t>(5));
71 size_t left = std::max(peak_center_index - searchBox, start_index);
72 if (searchBox > peak_center_index) {
73 // prevent overflow since these are unsigned
74 left = start_index;
75 }
76 const size_t rght = std::min(peak_center_index + searchBox, stop_index);
77
78 for (size_t i = left; i < rght; ++i) {
79 const double y = vector_y[i] - bkgd_values.getCalculated(i - start_index);
80 if (y > peak_height) {
81 peak_height = y;
82 peak_center = vector_x[i];
83 peak_center_index = i;
84 }
85 }
86
87 return GOOD;
88}
89
90//----------------------------------------------------------------------------------------------
101double observePeakFwhm(const Histogram &histogram, const FunctionValues &bkgd_values, size_t ipeak, size_t istart,
102 size_t istop, const EstimatePeakWidth peakWidthEstimateApproach,
103 const double peakWidthPercentage) {
104 double peak_fwhm(-0.);
105
106 if (peakWidthEstimateApproach == EstimatePeakWidth::InstrumentResolution) {
107 // width from guessing from delta(D)/D
108 const double peak_center = histogram.points()[ipeak];
109 peak_fwhm = peak_center * peakWidthPercentage;
110 } else if (peakWidthEstimateApproach == EstimatePeakWidth::Observation) {
111 // estimate fwhm from area assuming Gaussian (more robust than using
112 // moments which overestimates variance by factor ~5 depending on background
113 // estimation over window much wider than peak)
114 const auto &x_vec = histogram.points();
115 const auto &y_vec = histogram.y();
116 const size_t numPoints = std::min(istart - istop, bkgd_values.size()) - 1;
117
118 double area = 0.0;
119 // integrate using Newton's method
120 for (size_t i = 0; i < numPoints; ++i) {
121 // skip counting area if negative to give a better fwhm estimate
122 if (y_vec[istart + i] < 0.0) {
123 continue;
124 }
125 const auto yavg = 0.5 * (y_vec[istart + i] - bkgd_values.getCalculated(i) + y_vec[istart + i + 1] -
126 bkgd_values.getCalculated(i + 1));
127 const auto dx = x_vec[istart + i + 1] - x_vec[istart + i];
128 area += yavg * dx;
129 }
130 peak_fwhm = 2 * std::sqrt(M_LN2 / M_PI) * area / y_vec[ipeak];
131 } else {
132 // get from last peak or from input!
133 throw std::runtime_error("This case for observing peak width is not supported.");
134 }
135
136 return peak_fwhm;
137}
138
139//----------------------------------------------------------------------------------------------
156int estimatePeakParameters(const Histogram &histogram, const std::pair<size_t, size_t> &peak_window,
157 const API::IPeakFunction_sptr &peakfunction,
158 const API::IBackgroundFunction_sptr &bkgdfunction, bool observe_peak_width,
159 const EstimatePeakWidth peakWidthEstimateApproach, const double peakWidthPercentage,
160 const double minPeakHeight) {
161
162 // calculate background
163 const auto &vector_x = histogram.points();
164 FunctionDomain1DVector domain(vector_x.cbegin() + peak_window.first, vector_x.cbegin() + peak_window.second);
165 FunctionValues bkgd_values(domain);
166 bkgdfunction->function(domain, bkgd_values);
167
168 // Estimate peak center
169 double peak_height;
170 double peak_center = peakfunction->centre();
171 size_t peak_center_index;
172
173 int result = observePeakCenter(histogram, bkgd_values, peak_window.first, peak_window.second, peak_center,
174 peak_center_index, peak_height);
175
176 // return if failing to 'observe' peak center
177 if (result != GOOD)
178 return result;
179
180 if (peak_height < minPeakHeight || std::isnan(peak_height))
181 return LOWPEAK;
182
183 // use values from background to locate FWHM
184 peakfunction->setHeight(peak_height);
185 peakfunction->setCentre(peak_center);
186
187 // Estimate FHWM (peak width)
188 if (observe_peak_width && peakWidthEstimateApproach != EstimatePeakWidth::NoEstimation) {
189 const double peak_fwhm = observePeakFwhm(histogram, bkgd_values, peak_center_index, peak_window.first,
190 peak_window.second, peakWidthEstimateApproach, peakWidthPercentage);
191 if (peak_fwhm > 0.0) {
192 peakfunction->setFwhm(peak_fwhm);
193 }
194 }
195
196 return GOOD;
197}
198
199template MANTID_ALGORITHMS_DLL size_t findXIndex(const HistogramData::Points &, const double, const size_t);
200template MANTID_ALGORITHMS_DLL size_t findXIndex(const HistogramData::HistogramX &, const double, const size_t);
201template MANTID_ALGORITHMS_DLL size_t findXIndex(const std::vector<double> &, const double, const size_t);
202
203} // namespace Mantid::Algorithms::PeakParameterHelper
std::map< DeltaEMode::Type, std::string > index
Definition: DeltaEMode.cpp:19
double left
Definition: LineProfile.cpp:80
Implements FunctionDomain1D with its own storage in form of a std::vector.
A class to store values calculated by a function.
size_t size() const
Return the number of values.
double getCalculated(size_t i) const
Get i-th calculated value.
std::shared_ptr< IBackgroundFunction > IBackgroundFunction_sptr
std::shared_ptr< IPeakFunction > IPeakFunction_sptr
MANTID_ALGORITHMS_DLL double observePeakFwhm(const HistogramData::Histogram &histogram, const API::FunctionValues &bkgd_values, size_t ipeak, size_t istart, size_t istop, const EstimatePeakWidth peakWidthEstimateApproach, const double peakWidthPercentage)
Observe peak width.
MANTID_ALGORITHMS_DLL int observePeakCenter(const HistogramData::Histogram &histogram, const API::FunctionValues &bkgd_values, size_t start_index, size_t stop_index, double &peak_center, size_t &peak_center_index, double &peak_height)
observe peak center
MANTID_ALGORITHMS_DLL size_t findXIndex(const vector_like &vecx, const double x, const size_t startindex=0)
Get an index of a value in a sorted vector.
MANTID_ALGORITHMS_DLL int estimatePeakParameters(const HistogramData::Histogram &histogram, const std::pair< size_t, size_t > &peak_window, const API::IPeakFunction_sptr &peakfunction, const API::IBackgroundFunction_sptr &bkgdfunction, bool observe_peak_width, const EstimatePeakWidth peakWidthEstimateApproach, const double peakWidthPercentage, const double minPeakHeight)
Estimate peak parameters by 'observation'.
Helper class which provides the Collimation Length for SANS instruments.