Mantid
Loading...
Searching...
No Matches
StripPeaks.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 +
8
11#include "MantidAPI/TableRow.h"
17
18namespace Mantid::Algorithms {
19
20// Register the algorithm into the AlgorithmFactory
21DECLARE_ALGORITHM(StripPeaks)
22
23using namespace Kernel;
24using namespace API;
25
27 declareProperty(std::make_unique<WorkspaceProperty<>>("InputWorkspace", "", Direction::Input),
28 "The name of the input workspace.");
29 declareProperty(std::make_unique<WorkspaceProperty<>>("OutputWorkspace", "", Direction::Output),
30 "The name to use for the output workspace.");
31
32 auto min = std::make_shared<BoundedValidator<int>>();
33 min->setLower(1);
34 // The estimated width of a peak in terms of number of channels
35 declareProperty("FWHM", 7, min,
36 "The number of points covered, on average, "
37 "by the fwhm of a peak (default 7).\nPassed "
38 "through to [[FindPeaks]].");
39 // The tolerance allowed in meeting the conditions
40 declareProperty("Tolerance", 4, min,
41 "A measure of the strictness desired in "
42 "meeting the condition on peak "
43 "candidates,\n"
44 "Mariscotti recommends 2 (default 4)");
45
46 declareProperty(std::make_unique<ArrayProperty<double>>("PeakPositions"),
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.");
50
51 declareProperty("PeakPositionTolerance", 0.01,
52 "Tolerance on the found peaks' positions against the input "
53 "peak positions. Non-positive value indicates that this "
54 "option is turned off.");
55
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'");
59
60 declareProperty("HighBackground", true,
61 "Flag to indicate that the peaks are "
62 "relatively weak comparing to "
63 "background.");
64
65 auto mustBePositive = std::make_shared<BoundedValidator<int>>();
66 mustBePositive->setLower(0);
67 declareProperty("WorkspaceIndex", EMPTY_INT(), mustBePositive,
68 "If set, will remove peaks only in the given spectrum of the "
69 "workspace. Otherwise, all spectra will be searched.");
70
71 auto mustBePositiveDbl = std::make_shared<BoundedValidator<double>>();
72 mustBePositiveDbl->setLower(0.);
73 declareProperty("MaximumChisq", 100., mustBePositiveDbl,
74 "The maximum chisq value for fits to remove the peak. Default 100.");
75}
76
78 // Retrieve the input workspace
79 MatrixWorkspace_sptr inputWS = getProperty("InputWorkspace");
80 m_maxChiSq = getProperty("MaximumChisq");
81
82 // Call FindPeaks as a Child Algorithm
83 ITableWorkspace_sptr peakslist = this->findPeaks(inputWS);
84
85 MatrixWorkspace_sptr outputWS;
86 if (peakslist->rowCount() > 0) {
87 // Remove the peaks that were successfully fitted
88 outputWS = this->removePeaks(inputWS, peakslist);
89 } else {
90 // No peaks found!!! - just assign the input workspace pointer to the output
91 // one
92 g_log.warning("No peaks to remove!");
93 outputWS = inputWS;
94 }
95
96 // Set the output workspace
97 setProperty("OutputWorkspace", outputWS);
98}
99
105 g_log.debug("Calling FindPeaks as a Child Algorithm");
106
107 // Read from properties
108 int fwhm = getProperty("FWHM");
109 int tolerance = getProperty("Tolerance");
110 int wsindex = getProperty("WorkspaceIndex");
111 std::string backgroundtype = getProperty("BackgroundType");
112 bool highbackground = getProperty("HighBackground");
113 std::vector<double> peakpositions = getProperty("PeakPositions");
114 double peakpostol = getProperty("PeakPositionTolerance");
115 if (peakpostol < 0.)
116 peakpostol = EMPTY_DBL();
117
118 // Set up and execute algorithm
119 bool showlog = true;
120 auto findpeaks = createChildAlgorithm("FindPeaks", 0.0, 0.2, showlog);
121 findpeaks->setProperty("InputWorkspace", WS);
122 findpeaks->setProperty<int>("FWHM", fwhm);
123 findpeaks->setProperty<int>("Tolerance", tolerance);
124 findpeaks->setProperty<int>("WorkspaceIndex", wsindex);
125 // Get the specified peak positions, which is optional
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);
131
132 findpeaks->executeAsChildAlg();
133
134 ITableWorkspace_sptr peaklistws = findpeaks->getProperty("PeaksList");
135 if (!peaklistws)
136 throw std::runtime_error("Algorithm FindPeaks() returned a NULL pointer "
137 "for table workspace 'PeaksList'.");
138
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)
149 infoss << ", ";
150 }
151 infoss << ")\n"
152 << "Returned number of fitted peak = " << peaklistws->rowCount() << ".";
153 g_log.information(infoss.str());
154
155 return peaklistws;
156}
157
166 const API::ITableWorkspace_sptr &peakslist) {
167 g_log.debug("Subtracting peaks");
168 // Create an output workspace - same size as input one
169 MatrixWorkspace_sptr outputWS = WorkspaceFactory::Instance().create(input);
170 // Copy the data over from the input to the output workspace
171 const size_t hists = input->getNumberHistograms();
172 // progress 0.2 to 0.3 in this loop
173 double prg = 0.2;
174 for (size_t k = 0; k < hists; ++k) {
175 outputWS->setHistogram(k, input->histogram(k));
176 prg += (0.1 / static_cast<double>(hists));
177 progress(prg);
178 }
179
180 // progress from 0.3 to 1.0 here
181 prg = 0.3;
182 // Loop over the list of peaks
183 for (size_t i = 0; i < peakslist->rowCount(); ++i) {
184 // Get references to the data - X can be const
185 // but Y needs to be mutable
186 auto &X = outputWS->x(peakslist->getRef<int>("spectrum", i));
187 auto &Y = outputWS->mutableY(peakslist->getRef<int>("spectrum", i));
188 // Get back the gaussian parameters
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);
193 // These are some heuristic rules to discard bad fits.
194 // Hope to be able to remove them when we have better fitting routine
195 if (height < 0) {
196 g_log.error() << "Find Peak with Negative Height\n";
197 continue; // Height must be positive
198 }
199 if (chisq > m_maxChiSq) {
200 g_log.information() << "Error for fit peak at " << centre << " is " << chisq << ", which exceeds allowed value "
201 << m_maxChiSq << "\n";
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";
205 continue;
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
209 << " <= 0. \n";
210 }
211 {
212 // find the bin width at the center of the peak - average of the bins on
213 // either side
214 const auto center_iter = lower_bound(X.begin(), X.end(), centre);
215 const double bin_width = .5 * (*(center_iter + 1) - (*(center_iter - 1)));
216
217 // seek wikipedia if you don't believe the conversion factor
218 const double fwhm = sigma * 2. * std::sqrt(2. * std::log(2.));
219
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";
225 continue;
226 }
227 }
228
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";
231
232 { // log the background function
233 double a0 = 0.;
234 double a1 = 0.;
235 double a2 = 0.;
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";
244 }
245
246 // Get central bin values using points
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];
251 // Skip if not anywhere near this peak
252 if (x < centre - 5.0 * sigma)
253 continue;
254 if (x > centre + 5.0 * sigma)
255 break;
256 // Calculate the value of the Gaussian function at this point
257 const double funcVal = height * exp(-0.5 * pow((x - centre) / sigma, 2));
258 // Subtract the calculated value from the data
259 Y[j] -= funcVal;
260 }
261 prg += (0.7 / static_cast<double>(peakslist->rowCount()));
262 progress(prg);
263 }
264
265 return outputWS;
266}
267
268} // namespace Mantid::Algorithms
#define DECLARE_ALGORITHM(classname)
Definition: Algorithm.h:576
double sigma
Definition: GetAllEi.cpp:156
double height
Definition: GetAllEi.cpp:155
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
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
void progress(double p, const std::string &msg="", double estimatedTime=0.0, int progressPrecision=0)
Sends ProgressNotification.
Definition: Algorithm.cpp:231
A property class for workspaces.
void exec() override
Execution code.
Definition: StripPeaks.cpp:77
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.
Definition: StripPeaks.cpp:165
void init() override
Initialisation code.
Definition: StripPeaks.cpp:26
API::ITableWorkspace_sptr findPeaks(const API::MatrixWorkspace_sptr &WS)
Calls FindPeaks as a Child Algorithm.
Definition: StripPeaks.cpp:104
Support for a property that holds an array of values.
Definition: ArrayProperty.h:28
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.
Definition: Logger.cpp:114
void error(const std::string &msg)
Logs at error level.
Definition: Logger.cpp:77
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...
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.
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