23using namespace Kernel;
28 "The name of the input workspace.");
30 "The name to use for the output workspace.");
32 auto min = std::make_shared<BoundedValidator<int>>();
36 "The number of points covered, on average, "
37 "by the fwhm of a peak (default 7).\nPassed "
38 "through to [[FindPeaks]].");
41 "A measure of the strictness desired in "
42 "meeting the condition on peak "
44 "Mariscotti recommends 2 (default 4)");
47 "Optional: enter a comma-separated list of the expected "
48 "X-position of the centre of the peaks. Only peaks near "
49 "these positions will be fitted.");
52 "Tolerance on the found peaks' positions against the input "
53 "peak positions. Non-positive value indicates that this "
54 "option is turned off.");
56 std::vector<std::string> bkgdtypes{
"Linear",
"Quadratic"};
57 declareProperty(
"BackgroundType",
"Linear", std::make_shared<StringListValidator>(bkgdtypes),
58 "Type of Background. Present choices include 'Linear' and 'Quadratic'");
61 "Flag to indicate that the peaks are "
62 "relatively weak comparing to "
65 auto mustBePositive = std::make_shared<BoundedValidator<int>>();
66 mustBePositive->setLower(0);
68 "If set, will remove peaks only in the given spectrum of the "
69 "workspace. Otherwise, all spectra will be searched.");
71 auto mustBePositiveDbl = std::make_shared<BoundedValidator<double>>();
72 mustBePositiveDbl->setLower(0.);
74 "The maximum chisq value for fits to remove the peak. Default 100.");
86 if (peakslist->rowCount() > 0) {
105 g_log.
debug(
"Calling FindPeaks as a Child Algorithm");
111 std::string backgroundtype =
getProperty(
"BackgroundType");
112 bool highbackground =
getProperty(
"HighBackground");
113 std::vector<double> peakpositions =
getProperty(
"PeakPositions");
114 double peakpostol =
getProperty(
"PeakPositionTolerance");
121 findpeaks->setProperty(
"InputWorkspace", WS);
122 findpeaks->setProperty<
int>(
"FWHM", fwhm);
123 findpeaks->setProperty<
int>(
"Tolerance",
tolerance);
124 findpeaks->setProperty<
int>(
"WorkspaceIndex", wsindex);
126 findpeaks->setProperty<std::vector<double>>(
"PeakPositions", peakpositions);
127 findpeaks->setProperty<std::string>(
"BackgroundType", backgroundtype);
128 findpeaks->setProperty<
bool>(
"HighBackground", highbackground);
129 findpeaks->setProperty<
double>(
"PeakPositionTolerance", peakpostol);
130 findpeaks->setProperty<
bool>(
"RawPeakParameters",
true);
132 findpeaks->executeAsChildAlg();
136 throw std::runtime_error(
"Algorithm FindPeaks() returned a NULL pointer "
137 "for table workspace 'PeaksList'.");
139 std::stringstream infoss;
140 infoss <<
"Call FindPeaks() to find " << peakpositions.size() <<
" peaks with parameters: \n";
141 infoss <<
"\t FWHM = " << fwhm <<
";\n";
142 infoss <<
"\t Tolerance = " <<
tolerance <<
";\n";
143 infoss <<
"\t HighBackground = " << highbackground <<
";\n";
144 infoss <<
"\t BackgroundType = " << backgroundtype <<
";\n";
145 infoss <<
"\t Peak positions = ";
146 for (
size_t i = 0; i < peakpositions.size(); ++i) {
147 infoss << peakpositions[i];
148 if (i < peakpositions.size() - 1)
152 <<
"Returned number of fitted peak = " << peaklistws->rowCount() <<
".";
171 const size_t hists = input->getNumberHistograms();
174 for (
size_t k = 0; k < hists; ++k) {
175 outputWS->setHistogram(k, input->histogram(k));
176 prg += (0.1 /
static_cast<double>(hists));
183 for (
size_t i = 0; i < peakslist->rowCount(); ++i) {
186 auto &
X = outputWS->x(peakslist->getRef<
int>(
"spectrum", i));
187 auto &
Y = outputWS->mutableY(peakslist->getRef<
int>(
"spectrum", i));
189 const double height = peakslist->getRef<
double>(
"Height", i);
190 const double centre = peakslist->getRef<
double>(
"PeakCentre", i);
191 const double sigma = peakslist->getRef<
double>(
"Sigma", i);
192 const double chisq = peakslist->getRef<
double>(
"chi2", i);
196 g_log.
error() <<
"Find Peak with Negative Height\n";
200 g_log.
information() <<
"Error for fit peak at " << centre <<
" is " << chisq <<
", which exceeds allowed value "
202 if (chisq != DBL_MAX)
203 g_log.
error() <<
"StripPeaks(): Peak Index = " << i <<
" @ X = " << centre
204 <<
" Error: Peak fit with too high of chisq " << chisq <<
" > " <<
m_maxChiSq <<
"\n";
206 }
else if (chisq < 0.) {
207 g_log.
warning() <<
"StripPeaks(): Peak Index = " << i <<
" @ X = " << centre
208 <<
". Error: Peak fit with too wide peak width" <<
sigma <<
" denoted by chi^2 = " << chisq
214 const auto center_iter = lower_bound(
X.begin(),
X.end(), centre);
215 const double bin_width = .5 * (*(center_iter + 1) - (*(center_iter - 1)));
218 const double fwhm =
sigma * 2. * std::sqrt(2. * std::log(2.));
220 if ((fwhm / bin_width) < 1.5) {
221 g_log.
warning() <<
"StripPeaks(): Peak Index = " << i <<
" @ X = " << centre
222 <<
" Error: Peak fit with too narrow of peak"
223 <<
" fwhm = " << fwhm <<
" bin size = " << bin_width
224 <<
" num bins in peak = " << 2. * (fwhm / bin_width) <<
" <3\n";
229 g_log.
information() <<
"Subtracting peak " << i <<
" from spectrum " << peakslist->getRef<
int>(
"spectrum", i)
230 <<
" at x = " << centre <<
" h = " <<
height <<
" s = " <<
sigma <<
" chi2 = " << chisq <<
"\n";
236 const std::vector<std::string> columnNames = peakslist->getColumnNames();
237 if (std::find(columnNames.begin(), columnNames.end(),
"A0") != columnNames.end())
238 a0 = peakslist->getRef<
double>(
"A0", i);
239 if (std::find(columnNames.begin(), columnNames.end(),
"A1") != columnNames.end())
240 a1 = peakslist->getRef<
double>(
"A1", i);
241 if (std::find(columnNames.begin(), columnNames.end(),
"A2") != columnNames.end())
242 a2 = peakslist->getRef<
double>(
"A2", i);
243 g_log.
information() <<
" background = " << a0 <<
" + " << a1 <<
" x + " << a2 <<
" x^2\n";
247 auto binCenters = outputWS->points(peakslist->getRef<
int>(
"spectrum", i));
248 const auto spectrumLength =
static_cast<int>(
Y.size());
249 for (
int j = 0; j < spectrumLength; ++j) {
250 double x = binCenters[j];
252 if (
x < centre - 5.0 *
sigma)
254 if (
x > centre + 5.0 *
sigma)
257 const double funcVal =
height * exp(-0.5 * pow((
x - centre) /
sigma, 2));
261 prg += (0.7 /
static_cast<double>(peakslist->rowCount()));
#define DECLARE_ALGORITHM(classname)
void declareProperty(std::unique_ptr< Kernel::Property > p, const std::string &doc="") override
Add a property to the list of managed properties.
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.
void progress(double p, const std::string &msg="", double estimatedTime=0.0, int progressPrecision=0)
Sends ProgressNotification.
A property class for workspaces.
void exec() override
Execution code.
API::MatrixWorkspace_sptr removePeaks(const API::MatrixWorkspace_const_sptr &input, const API::ITableWorkspace_sptr &peakslist)
If a peak was successfully fitted, it is subtracted from the data.
void init() override
Initialisation code.
API::ITableWorkspace_sptr findPeaks(const API::MatrixWorkspace_sptr &WS)
Calls FindPeaks as a Child Algorithm.
Support for a property that holds an array of values.
IPropertyManager * setProperty(const std::string &name, const T &value)
Templated method to set the value of a PropertyWithValue.
void debug(const std::string &msg)
Logs at debug level.
void error(const std::string &msg)
Logs at error level.
void warning(const std::string &msg)
Logs at warning level.
void information(const std::string &msg)
Logs at information level.
static T & Instance()
Return a reference to the Singleton instance, creating it if it does not already exist Creation is do...
std::shared_ptr< ITableWorkspace > ITableWorkspace_sptr
shared pointer to Mantid::API::ITableWorkspace
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
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.