15#include "MantidHistogramData/HistogramIterator.h"
16#include "MantidHistogramData/LinearGenerator.h"
23const std::string COMPONENTS{
"Components"};
24const std::string DIAGNOSTICS_WS{
"DiagnosticsWorkspace"};
25const std::string EPP_WS{
"EPPWorkspace"};
26const std::string INPUT_WS{
"InputWorkspace"};
27const std::string OUTPUT_WS{
"OutputWorkspace"};
28const std::string POLYNOMIAL_DEGREE{
"Degree"};
29const std::string SIGMA_MULTIPLIER{
"NonBkgRegionInSigmas"};
39 size_t const firstColumnIndex) {
40 std::vector<double> ranges;
42 bool needsRangeStart{
true};
43 bool needsRangeEnd{
false};
45 if (spectrumInfo.isMasked(i) || statuses.
cell<std::string>(i + firstColumnIndex) !=
"success") {
47 needsRangeStart =
true;
48 needsRangeEnd =
false;
50 ranges.emplace_back(
static_cast<double>(i) - 0.5);
54 if (needsRangeStart) {
55 needsRangeStart =
false;
57 ranges.emplace_back(
static_cast<double>(i) - 0.5);
69 std::vector<double> peakStarts;
71 std::vector<double> peakEnds;
84PeakBounds peakBounds(
size_t const firstIndex,
size_t const lastIndex,
double const sigmaMultiplier,
88 size_t const boundsSize = lastIndex - firstIndex + 1;
89 bounds.peakStarts.resize(boundsSize);
90 bounds.peakEnds.resize(boundsSize);
91 for (
size_t i = firstIndex; i <= lastIndex; ++i) {
92 auto const boundsIndex = i - firstIndex;
93 if (fitStatusColumn.
cell<std::string>(i) ==
"success") {
94 auto const peakCentre = peakCentreColumn.
cell<
double>(i);
95 auto const peakWidth = sigmaMultiplier * sigmaColumn.
cell<
double>(i);
96 bounds.peakStarts[boundsIndex] = peakCentre - peakWidth;
97 bounds.peakEnds[boundsIndex] = peakCentre + peakWidth;
99 bounds.peakStarts[boundsIndex] = std::numeric_limits<double>::lowest();
100 bounds.peakEnds[boundsIndex] = std::numeric_limits<double>::max();
117 throw NotFoundError(
"EPPWorkspace does not contain 'PeakCentre' column.",
"PeakCentre");
120 throw NotFoundError(
"EPPWorkspace does not contain 'Sigma' column",
"Sigma");
123 throw NotFoundError(
"EPPWorkspace does not contain 'FitStatus' column.",
"FitStatus");
137 throw NotFoundError(
"Component not found in InputWorkspace's instrument.", componentName);
174 auto const &ys = componentBkgWS.
y(0);
175 auto const &es = componentBkgWS.
e(0);
176 for (
size_t i = 0; i < ys.size(); ++i) {
177 targetWS.
mutableY(firstTargetWSIndex + i) = ys[i];
178 targetWS.
mutableE(firstTargetWSIndex + i) = es[i];
197 return "CorrectionFunctions\\BackgroundCorrections;ILL\\Direct";
202 return {
"CalculateFlatBackground",
"CalculatePolynomialBackground",
"CreateUserDefinedBackground",
"RemoveBackground",
208 return "Fits polynomial backgrounds over the pixels of position sensitive "
217 "A workspace to fit the backgrounds to.");
220 "The fitted backgrounds.");
222 "A list of component names for which to calculate the backgrounds.");
225 "A table workspace containing results from the FindEPP algorithm.");
226 auto positiveDouble = std::make_shared<Kernel::BoundedValidator<double>>();
227 positiveDouble->setLowerExclusive(0.);
229 "Half width of the range excluded from background around "
230 "the elastic peaks in multiplies of 'Sigma' in the EPP table.'.");
231 auto nonnegativeInt = std::make_shared<Kernel::BoundedValidator<int>>();
232 nonnegativeInt->setLower(0);
233 declareProperty(Prop::POLYNOMIAL_DEGREE, 0, nonnegativeInt,
"The degree of the background polynomial.");
236 "Detector diagnostics workspace for masking.");
241 std::map<std::string, std::string> issues;
244 if (inWS->getNumberHistograms() != eppWS->rowCount()) {
245 issues[Prop::EPP_WS] =
"Wrong EPP workspace? The number of the table rows "
246 "does not match the number of histograms in InputWorkspace.";
250 if (inWS->getNumberHistograms() != maskWS->getNumberHistograms()) {
251 issues[Prop::EPP_WS] =
"Wrong diagnostics workspace? The number of histograms "
252 "does not match with InputWorkspace.";
264 for (
size_t i = 0; i < bkgWS->getNumberHistograms(); ++i) {
265 bkgWS->convertToFrequencies(i);
268 double const sigmaMultiplier =
getProperty(Prop::SIGMA_MULTIPLIER);
269 auto peakCentreColumn = eppWS->getColumn(
"PeakCentre");
270 auto sigmaColumn = eppWS->getColumn(
"Sigma");
271 auto fitStatusColumn = eppWS->getColumn(
"FitStatus");
272 checkEPPColumnsExist(peakCentreColumn, sigmaColumn, fitStatusColumn);
274 std::vector<std::string>
const componentNames =
components(*instrument);
277 for (int64_t i = 0; i < static_cast<int64_t>(componentNames.size()); ++i) {
279 auto const &componentName = componentNames[
static_cast<size_t>(i)];
280 progress.report(
"Processing " + componentName);
281 checkComponentExists(componentName, *instrument);
283 auto const wsIndexRange = componentWSIndexRange(*componentWS, *ws);
284 auto const bkgRanges = bkgFittingRanges(*componentWS, *fitStatusColumn, wsIndexRange.first);
285 if (bkgRanges.empty()) {
288 auto const bounds = peakBounds(wsIndexRange.first, wsIndexRange.last, sigmaMultiplier, *peakCentreColumn,
289 *sigmaColumn, *fitStatusColumn);
292 writeComponentBackgroundToOutput(*fittedComponentBkg, *bkgWS, wsIndexRange.first);
297 for (
size_t i = 0; i < bkgWS->getNumberHistograms(); ++i) {
298 bkgWS->convertToCounts(i);
315 mask->setProperty(
"Workspace", ws);
316 mask->setProperty(
"MaskedWorkspace", diagnosticsWS);
329 std::string
const COMPONENTS_PARAMETER{
"components-for-backgrounds"};
331 throw std::runtime_error(
"Could not find '" + COMPONENTS_PARAMETER +
332 "' in instrument parameters file. Component "
333 "names must be given using the '" +
334 Prop::COMPONENTS +
"' property.");
336 auto const componentList = instrument.
getStringParameter(COMPONENTS_PARAMETER).front();
349 std::string
const &componentName) {
351 crop->setProperty(
"InputWorkspace", ws);
352 crop->setProperty(
"OutputWorkspace",
"_unused");
353 crop->setProperty(
"ComponentNames", std::vector<std::string>{componentName});
355 return crop->getProperty(
"OutputWorkspace");
365 std::vector<double>
const &xRanges) {
366 int const degree =
getProperty(Prop::POLYNOMIAL_DEGREE);
368 calculateBkg->setProperty(
"InputWorkspace", ws);
369 calculateBkg->setProperty(
"OutputWorkspace",
"_unused");
370 calculateBkg->setProperty(
"Degree", degree);
371 calculateBkg->setProperty(
"XRanges", xRanges);
372 calculateBkg->setProperty(
"CostFunction",
"Unweighted least squares");
373 calculateBkg->execute();
374 return calculateBkg->getProperty(
"OutputWorkspace");
385 std::vector<double>
const &peakStarts,
386 std::vector<double>
const &peakEnds) {
387 HistogramData::Points wsIndices{ws.
getNumberHistograms(), HistogramData::LinearGenerator(0., 1.)};
391 HistogramData::Histogram modelHistogram{wsIndices, zeroCounts};
394 auto const peakStart = peakStarts[i];
395 auto const peakEnd = peakEnds[i];
398 double &sum = averageWS->mutableY(0)[i];
399 double &
error = averageWS->mutableE(0)[i];
400 for (
auto &histogramItem : ws.
histogram(i)) {
401 auto const center = histogramItem.center();
402 if (peakStart <= center && center <= peakEnd) {
405 sum += histogramItem.frequency();
406 auto const stdDev = histogramItem.frequencyStandardDeviation();
410 if (itemCount != 0) {
411 sum /=
static_cast<double>(itemCount);
412 error /=
static_cast<double>(itemCount);
#define DECLARE_ALGORITHM(classname)
#define PARALLEL_START_INTERRUPT_REGION
Begins a block to skip processing is the algorithm has been interupted Note the end of the block if n...
#define PARALLEL_END_INTERRUPT_REGION
Ends a block to skip processing is the algorithm has been interupted Note the start of the block if n...
#define PARALLEL_FOR_IF(condition)
Empty definitions - to enable set your complier to enable openMP.
#define PARALLEL_CHECK_INTERRUPT_REGION
Adds a check after a Parallel region to see if it was interupted.
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.
bool isDefault(const std::string &name) const
void setPropertyValue(const std::string &name, const std::string &value) override
Set the value of a property by string N.B.
Column is the base class for columns of TableWorkspace.
T & cell(size_t index)
Templated method for returning a value. No type checks are done.
const SpectrumInfo & spectrumInfo() const
Return a reference to the SpectrumInfo object.
Geometry::Instrument_const_sptr getInstrument() const
Returns the parameterized instrument.
specnum_t getSpectrumNo() const
Base MatrixWorkspace Abstract Class.
virtual ISpectrum & getSpectrum(const size_t index)=0
Return the underlying ISpectrum ptr at the given workspace index.
const HistogramData::HistogramE & e(const size_t index) const
virtual std::size_t getNumberHistograms() const =0
Returns the number of histograms in the workspace.
bool isDistribution() const
Are the Y-values dimensioned?
HistogramData::Histogram histogram(const size_t index) const
Returns the Histogram at the given workspace index.
HistogramData::HistogramE & mutableE(const size_t index) &
size_t getIndexFromSpectrumNumber(const specnum_t specNo) const
Given a spectrum number, find the corresponding workspace index.
HistogramData::HistogramY & mutableY(const size_t index) &
const HistogramData::HistogramY & y(const size_t index) const
Helper class for reporting progress from algorithms.
A property class for workspaces.
DirectILLTubeBackground : Fits polynomial backgrounds over the pixels of position sensitive tubes.
API::MatrixWorkspace_sptr applyDiagnostics(API::MatrixWorkspace_sptr ws)
Apply a mask workspace (if given) to ws.
const std::string summary() const override
Algorithm's summary for use in the GUI and help.
API::MatrixWorkspace_sptr fitComponentBackground(API::MatrixWorkspace_sptr &ws, std::vector< double > const &xRanges)
Fits a polynomial background.
int version() const override
Algorithm's version for identification.
std::map< std::string, std::string > validateInputs() override
Validate input properties.
void exec() override
Execute the algorithm.
void init() override
Initialize the algorithm's properties.
std::vector< std::string > components(Geometry::Instrument const &instrument)
Return a list of component names for the algorithm to process.
API::MatrixWorkspace_sptr peakExcludingAverage(const API::MatrixWorkspace &ws, std::vector< double > const &peakStarts, std::vector< double > const &peakEnds)
Average the histograms of a workspace.
const std::string category() const override
Algorithm's category for identification.
const std::vector< std::string > seeAlso() const override
Return a vector of related algorithms.
API::MatrixWorkspace_sptr cropToComponent(API::MatrixWorkspace_sptr &ws, std::string const &componentName)
Crop a component workspace out of ws.
std::shared_ptr< const IComponent > getComponentByName(const std::string &cname, int nlevels=0) const override
Returns a pointer to the first component of assembly encountered with the given name.
std::vector< std::string > getStringParameter(const std::string &pname, bool recursive=true) const override
Get a parameter defined as a string.
bool hasParameter(const std::string &name, bool recursive=true) const override
Returns a boolean indicating if the component has the named parameter.
Support for a property that holds an array of values.
Exception for when an item is not found in a collection.
IPropertyManager * setProperty(const std::string &name, const T &value)
Templated method to set the value of a PropertyWithValue.
std::shared_ptr< Column > Column_sptr
std::shared_ptr< ITableWorkspace > ITableWorkspace_sptr
shared pointer to Mantid::API::ITableWorkspace
std::shared_ptr< MatrixWorkspace > MatrixWorkspace_sptr
shared pointer to the matrix workspace base class
The exception classes used by Mantid.
std::enable_if< std::is_pointer< Arg >::value, bool >::type threadSafe(Arg workspace)
Thread-safety check Checks the workspace to ensure it is suitable for multithreaded access.
String constants for algorithm's properties.
@ Input
An input workspace.
@ Output
An output workspace.