28 double operator()(
double x)
const {
return x *
x; }
32const double guessSignalToNoiseRatio = 1e15;
52 "Workspace indices for spectra to process. "
53 "If empty smooth all spectra.");
55 "An output workspace.");
64 std::vector<int> wsIndexList =
getProperty(
"WorkspaceIndexList");
67 const size_t nInputSpectra = inputWS->getNumberHistograms();
70 if (wsIndexList.size() > nInputSpectra) {
71 throw std::invalid_argument(
"Workspace index list has more indices than "
72 "there are spectra in the input workspace.");
76 if (wsIndexList.empty()) {
78 wsIndexList.resize(nInputSpectra);
79 wsIndexList.front() = 0;
80 for (
auto index = wsIndexList.begin() + 1;
index != wsIndexList.end(); ++
index) {
86 const size_t nOutputSpectra = wsIndexList.size();
89 auto wsIndex =
static_cast<size_t>(wsIndexList.front());
100 auto inAxis = inputWS->getAxis(1);
101 auto outAxis = std::unique_ptr<API::Axis>(inAxis->clone(nOutputSpectra, outputWS.get()));
103 bool isSpectra = outAxis->isSpectra();
104 bool isNumeric = outAxis->isNumeric();
106 auto outTextAxis =
dynamic_cast<API::TextAxis *
>(outAxis.get());
112 for (
size_t outIndex = 0; outIndex < nOutputSpectra; ++outIndex) {
113 auto inIndex = wsIndexList[outIndex];
117 outputWS->setHistogram(outIndex, next->histogram(0));
121 auto &inSpectrum = inputWS->getSpectrum(inIndex);
122 auto &outSpectrum = outputWS->getSpectrum(outIndex);
123 outSpectrum.setSpectrumNo(inSpectrum.getSpectrumNo());
124 outSpectrum.setDetectorIDs(inSpectrum.getDetectorIDs());
125 }
else if (isNumeric) {
126 outAxis->setValue(outIndex, inAxis->getValue(inIndex));
127 }
else if (inTextAxis && outTextAxis) {
128 outTextAxis->setLabel(outIndex, inTextAxis->label(inIndex));
132 outputWS->replaceAxis(1, std::move(outAxis));
145 size_t dataSize = inputWS->blocksize();
149 g_log.
debug() <<
"No smoothing, spectrum copied.\n";
154 const bool isOddSize = dataSize % 2 != 0;
160 auto &
X = inputWS->x(wsIndex);
161 auto &
Y = inputWS->y(wsIndex);
162 auto &E = inputWS->e(wsIndex);
163 double dx =
X[dataSize - 1] -
X[dataSize - 2];
164 auto histogram = inputWS->histogram(wsIndex);
165 histogram.resize(histogram.size() + 1);
166 histogram.mutableX().back() =
X.back() + dx;
167 histogram.mutableY().back() =
Y.back();
168 histogram.mutableE().back() = E.back();
169 inputWS->setHistogram(wsIndex, histogram);
173 auto &
X = inputWS->x(wsIndex);
174 auto &
Y = inputWS->y(wsIndex);
187 if (nbreak * 3 > dataSize)
188 nbreak = dataSize / 3;
196 g_log.
debug() <<
"Spline break points " << nbreak <<
'\n';
200 auto xInterval =
getStartEnd(
X, inputWS->isHistogramData());
202 spline->setAttributeValue(
"NBreak",
static_cast<int>(nbreak));
204 spline->setParameter(0,
Y.front());
206 size_t lastParamIndex = spline->nParams() - 1;
207 spline->setParameter(lastParamIndex,
Y.back());
208 spline->fix(lastParamIndex);
213 fit->setProperty(
"Function", spline);
214 fit->setProperty(
"InputWorkspace", inputWS);
215 fit->setProperty(
"WorkspaceIndex",
static_cast<int>(wsIndex));
216 fit->setProperty(
"CreateOutput",
true);
225 fourier->initialize();
226 fourier->setProperty(
"InputWorkspace", fitOut);
227 fourier->setProperty(
"WorkspaceIndex", 2);
229 fourier->setProperty(
"IgnoreXBins",
true);
237 auto &powerSpec = fourierOut->mutableY(2);
239 std::transform(powerSpec.begin(), powerSpec.end(), powerSpec.begin(), PowerSpectrum());
242 size_t n2 = powerSpec.size();
243 double noise = std::accumulate(powerSpec.begin() + n2 / 2, powerSpec.end(), 0.0);
244 noise /=
static_cast<double>(n2);
248 static_cast<size_t>(std::distance(powerSpec.begin(), std::max_element(powerSpec.begin(), powerSpec.end())));
251 noise = powerSpec[imax] / guessSignalToNoiseRatio;
254 g_log.
debug() <<
"Maximum signal " << powerSpec[imax] <<
'\n';
258 std::vector<double> wf(n2);
276 for (
size_t i = 0; i < n2; ++i) {
277 double cd1 = powerSpec[i] / noise;
278 if (cd1 < 1.0 && i > imax) {
282 double cd2 = log(cd1);
283 wf[i] = cd1 / (1.0 + cd1);
284 auto j =
static_cast<double>(i + 1);
292 g_log.
debug() <<
"Noise start index " << i0 <<
'\n';
295 auto ri0f =
static_cast<double>(i0 + 1);
296 double xm = (1.0 + ri0f) / 2;
298 double a1 = (xy - ri0f * xm * ym) / (xx - ri0f * xm * xm);
299 double b1 = ym - a1 * xm;
301 g_log.
debug() <<
"(a1,b1) = (" << a1 <<
',' << b1 <<
')' <<
'\n';
303 const double dblev = -20.0;
305 double ri1 = floor((dblev / 4 - b1) / a1);
306 if (ri1 <
static_cast<double>(i0)) {
307 g_log.
warning() <<
"Failed to build Wiener filter: no smoothing.\n";
308 ri1 =
static_cast<double>(i0);
310 auto i1 =
static_cast<size_t>(ri1);
313 for (
size_t i = i0; i < i1; ++i) {
314 double s = exp(a1 *
static_cast<double>(i + 1) + b1);
315 wf[i] = s / (1.0 + s);
319 g_log.
debug() <<
"Cut-off index " << i1 <<
'\n';
321 g_log.
warning() <<
"Power spectrum has an unexpected shape: no smoothing\n";
326 auto &re = fourierOut->mutableY(0);
327 auto &im = fourierOut->mutableY(1);
329 std::transform(re.begin(), re.end(), wf.begin(), re.begin(), std::multiplies<double>());
330 std::transform(im.begin(), im.end(), wf.begin(), im.begin(), std::multiplies<double>());
334 fourier->initialize();
335 fourier->setProperty(
"InputWorkspace", fourierOut);
336 fourier->setProperty(
"IgnoreXBins",
true);
337 fourier->setPropertyValue(
"Transform",
"Backward");
341 auto &background = fitOut->y(1);
342 auto &
y = out->mutableY(0);
344 if (
y.size() != background.size()) {
345 throw std::logic_error(
"Logic error: inconsistent arrays");
349 std::transform(
y.begin(),
y.end(), background.begin(),
y.begin(), std::plus<double>());
355 auto histogram = out->histogram(0);
356 histogram.setSharedX(inputWS->sharedX(wsIndex));
357 histogram.setSharedE(inputWS->sharedE(wsIndex));
359 auto newSize = histogram.y().size() - 1;
360 histogram.resize(newSize);
362 out->setHistogram(0, histogram);
364 out->setSharedX(0, inputWS->sharedX(wsIndex));
365 out->setSharedE(0, inputWS->sharedE(wsIndex));
379 bool isHistogram)
const {
380 const size_t n =
X.size();
383 throw std::runtime_error(
"Number of bins/data points cannot be smaller than 3.");
386 return std::make_pair((
X[0] +
X[1]) / 2, (
X[
n - 1] +
X[
n - 2]) / 2);
389 return std::make_pair(
X.front(),
X.back());
401 alg->setProperty(
"InputWorkspace", inputWS);
402 alg->setProperty(
"WorkspaceIndex",
static_cast<int>(wsIndex));
409 const double currentEndX = spline->getAttribute(
"EndX").asDouble();
410 if (startX >= currentEndX) {
411 spline->setAttributeValue(
"EndX", endX);
412 spline->setAttributeValue(
"StartX", startX);
414 spline->setAttributeValue(
"StartX", startX);
415 spline->setAttributeValue(
"EndX", endX);
#define DECLARE_ALGORITHM(classname)
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.
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.
Helper class for reporting progress from algorithms.
Class to represent a text axis of a workspace.
A property class for workspaces.
void init() override
Initialize the algorithm's properties.
void setSplineDataBounds(API::IFunction_sptr &spline, const double startX, const double endX)
std::pair< double, double > getStartEnd(const Mantid::HistogramData::HistogramX &X, bool isHistogram) const
Get the start and end of the x-interval.
void exec() override
Execute the algorithm.
const std::string category() const override
Algorithm's category for identification.
const std::string summary() const override
Algorithm's summary for use in the GUI and help.
API::MatrixWorkspace_sptr copyInput(const API::MatrixWorkspace_sptr &inputWS, size_t wsIndex)
Exctract the input spectrum into a separate workspace.
API::MatrixWorkspace_sptr smoothSingleSpectrum(API::MatrixWorkspace_sptr inputWS, size_t wsIndex)
Execute smoothing of a single spectrum.
int version() const override
Algorithm's version for identification.
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 warning(const std::string &msg)
Logs at warning 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< IFunction > IFunction_sptr
shared pointer to the function base class
std::shared_ptr< MatrixWorkspace > MatrixWorkspace_sptr
shared pointer to the matrix workspace base class
Describes the direction (within an algorithm) of a Property.
@ Input
An input workspace.
@ Output
An output workspace.