Mantid
Loading...
Searching...
No Matches
DirectILLTubeBackground.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
15#include "MantidHistogramData/HistogramIterator.h"
16#include "MantidHistogramData/LinearGenerator.h"
19
20namespace {
22namespace Prop {
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"};
30} // namespace Prop
31
38std::vector<double> bkgFittingRanges(Mantid::API::MatrixWorkspace const &ws, Mantid::API::Column const &statuses,
39 size_t const firstColumnIndex) {
40 std::vector<double> ranges;
41 auto const &spectrumInfo = ws.spectrumInfo();
42 bool needsRangeStart{true};
43 bool needsRangeEnd{false};
44 for (size_t i = 0; i < ws.getNumberHistograms(); ++i) {
45 if (spectrumInfo.isMasked(i) || statuses.cell<std::string>(i + firstColumnIndex) != "success") {
46 if (needsRangeEnd) {
47 needsRangeStart = true;
48 needsRangeEnd = false;
49 // Current spectrum number is the first to exclude.
50 ranges.emplace_back(static_cast<double>(i) - 0.5);
51 }
52 continue;
53 }
54 if (needsRangeStart) {
55 needsRangeStart = false;
56 needsRangeEnd = true;
57 ranges.emplace_back(static_cast<double>(i) - 0.5);
58 }
59 }
60 if (needsRangeEnd) {
61 ranges.emplace_back(static_cast<double>(ws.getNumberHistograms()) - 0.5);
62 }
63 return ranges;
64}
65
67struct PeakBounds {
69 std::vector<double> peakStarts;
71 std::vector<double> peakEnds;
72};
73
84PeakBounds peakBounds(size_t const firstIndex, size_t const lastIndex, double const sigmaMultiplier,
85 Mantid::API::Column const &peakCentreColumn, Mantid::API::Column const &sigmaColumn,
86 Mantid::API::Column const &fitStatusColumn) {
87 PeakBounds bounds;
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;
98 } else {
99 bounds.peakStarts[boundsIndex] = std::numeric_limits<double>::lowest();
100 bounds.peakEnds[boundsIndex] = std::numeric_limits<double>::max();
101 }
102 }
103 return bounds;
104}
105
113void checkEPPColumnsExist(Mantid::API::Column_sptr const &centreColumn, Mantid::API::Column_sptr const &sigmaColumn,
114 Mantid::API::Column_sptr const &statusColumn) {
115 using namespace Mantid::Kernel::Exception;
116 if (!centreColumn) {
117 throw NotFoundError("EPPWorkspace does not contain 'PeakCentre' column.", "PeakCentre");
118 }
119 if (!sigmaColumn) {
120 throw NotFoundError("EPPWorkspace does not contain 'Sigma' column", "Sigma");
121 }
122 if (!statusColumn) {
123 throw NotFoundError("EPPWorkspace does not contain 'FitStatus' column.", "FitStatus");
124 }
125}
126
133void checkComponentExists(std::string const &componentName, Mantid::Geometry::Instrument const &instrument) {
134 using namespace Mantid::Kernel::Exception;
135 auto const component = instrument.getComponentByName(componentName);
136 if (!component) {
137 throw NotFoundError("Component not found in InputWorkspace's instrument.", componentName);
138 }
139}
140
142struct Range {
144 size_t first{0};
146 size_t last{0};
147};
148
155Range componentWSIndexRange(Mantid::API::MatrixWorkspace const &componentWS,
156 Mantid::API::MatrixWorkspace const &originalWS) {
157 Range range;
158 auto const firstComponentSpectrumNo = componentWS.getSpectrum(0).getSpectrumNo();
159 range.first = originalWS.getIndexFromSpectrumNumber(firstComponentSpectrumNo);
160 auto const nComponentHistograms = componentWS.getNumberHistograms();
161 auto const lastComponentSpectrumNo = componentWS.getSpectrum(nComponentHistograms - 1).getSpectrumNo();
162 range.last = originalWS.getIndexFromSpectrumNumber(lastComponentSpectrumNo);
163 return range;
164}
165
172void writeComponentBackgroundToOutput(Mantid::API::MatrixWorkspace const &componentBkgWS,
173 Mantid::API::MatrixWorkspace &targetWS, size_t const firstTargetWSIndex) {
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];
179 }
180}
181
182} // namespace
183
184namespace Mantid::Algorithms {
185
186// Register the algorithm into the AlgorithmFactory
187DECLARE_ALGORITHM(DirectILLTubeBackground)
188
189
190const std::string DirectILLTubeBackground::name() const { return "DirectILLTubeBackground"; }
191
193int DirectILLTubeBackground::version() const { return 1; }
194
196const std::string DirectILLTubeBackground::category() const {
197 return "CorrectionFunctions\\BackgroundCorrections;ILL\\Direct";
198}
199
201const std::vector<std::string> DirectILLTubeBackground::seeAlso() const {
202 return {"CalculateFlatBackground", "CalculatePolynomialBackground", "CreateUserDefinedBackground", "RemoveBackground",
203 "SplineBackground"};
204}
205
207const std::string DirectILLTubeBackground::summary() const {
208 return "Fits polynomial backgrounds over the pixels of position sensitive "
209 "tubes.";
210}
211
216 std::make_unique<API::WorkspaceProperty<API::MatrixWorkspace>>(Prop::INPUT_WS, "", Kernel::Direction::Input),
217 "A workspace to fit the backgrounds to.");
219 std::make_unique<API::WorkspaceProperty<API::MatrixWorkspace>>(Prop::OUTPUT_WS, "", Kernel::Direction::Output),
220 "The fitted backgrounds.");
221 declareProperty(std::make_unique<Kernel::ArrayProperty<std::string>>(Prop::COMPONENTS, std::vector<std::string>()),
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.);
228 declareProperty(Prop::SIGMA_MULTIPLIER, 6., positiveDouble,
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.");
235 Prop::DIAGNOSTICS_WS, "", Kernel::Direction::Input, API::PropertyMode::Optional),
236 "Detector diagnostics workspace for masking.");
237}
238
240std::map<std::string, std::string> DirectILLTubeBackground::validateInputs() {
241 std::map<std::string, std::string> issues;
242 API::MatrixWorkspace_sptr inWS = getProperty(Prop::INPUT_WS);
243 API::ITableWorkspace_sptr eppWS = getProperty(Prop::EPP_WS);
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.";
247 }
248 if (!isDefault(Prop::DIAGNOSTICS_WS)) {
249 API::MatrixWorkspace_sptr maskWS = getProperty(Prop::DIAGNOSTICS_WS);
250 if (inWS->getNumberHistograms() != maskWS->getNumberHistograms()) {
251 issues[Prop::EPP_WS] = "Wrong diagnostics workspace? The number of histograms "
252 "does not match with InputWorkspace.";
253 }
254 }
255 return issues;
256}
257
261 API::MatrixWorkspace_sptr inWS = getProperty(Prop::INPUT_WS);
262 auto ws = applyDiagnostics(inWS->clone());
263 API::MatrixWorkspace_sptr bkgWS = DataObjects::create<DataObjects::Workspace2D>(*ws);
264 for (size_t i = 0; i < bkgWS->getNumberHistograms(); ++i) {
265 bkgWS->convertToFrequencies(i);
266 }
267 API::ITableWorkspace_sptr eppWS = getProperty(Prop::EPP_WS);
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);
273 auto instrument = ws->getInstrument();
274 std::vector<std::string> const componentNames = components(*instrument);
275 API::Progress progress(this, 0.0, 1.0, componentNames.size());
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);
282 auto componentWS = cropToComponent(ws, componentName);
283 auto const wsIndexRange = componentWSIndexRange(*componentWS, *ws);
284 auto const bkgRanges = bkgFittingRanges(*componentWS, *fitStatusColumn, wsIndexRange.first);
285 if (bkgRanges.empty()) {
286 continue;
287 }
288 auto const bounds = peakBounds(wsIndexRange.first, wsIndexRange.last, sigmaMultiplier, *peakCentreColumn,
289 *sigmaColumn, *fitStatusColumn);
290 auto averageWS = peakExcludingAverage(*componentWS, bounds.peakStarts, bounds.peakEnds);
291 auto fittedComponentBkg = fitComponentBackground(averageWS, bkgRanges);
292 writeComponentBackgroundToOutput(*fittedComponentBkg, *bkgWS, wsIndexRange.first);
294 }
296 if (!ws->isDistribution()) {
297 for (size_t i = 0; i < bkgWS->getNumberHistograms(); ++i) {
298 bkgWS->convertToCounts(i);
299 }
300 }
301 setProperty(Prop::OUTPUT_WS, bkgWS);
302}
303
310 if (isDefault(Prop::DIAGNOSTICS_WS)) {
311 return ws;
312 }
313 API::MatrixWorkspace_sptr diagnosticsWS = getProperty(Prop::DIAGNOSTICS_WS);
314 auto mask = createChildAlgorithm("MaskDetectors");
315 mask->setProperty("Workspace", ws);
316 mask->setProperty("MaskedWorkspace", diagnosticsWS);
317 mask->execute();
318 return ws;
319}
320
327std::vector<std::string> DirectILLTubeBackground::components(Geometry::Instrument const &instrument) {
328 if (isDefault(Prop::COMPONENTS)) {
329 std::string const COMPONENTS_PARAMETER{"components-for-backgrounds"};
330 if (!instrument.hasParameter(COMPONENTS_PARAMETER)) {
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.");
335 }
336 auto const componentList = instrument.getStringParameter(COMPONENTS_PARAMETER).front();
337 setPropertyValue(Prop::COMPONENTS, componentList);
338 }
339 return getProperty(Prop::COMPONENTS);
340}
341
349 std::string const &componentName) {
350 auto crop = createChildAlgorithm("CropToComponent");
351 crop->setProperty("InputWorkspace", ws);
352 crop->setProperty("OutputWorkspace", "_unused");
353 crop->setProperty("ComponentNames", std::vector<std::string>{componentName});
354 crop->execute();
355 return crop->getProperty("OutputWorkspace");
356}
357
365 std::vector<double> const &xRanges) {
366 int const degree = getProperty(Prop::POLYNOMIAL_DEGREE);
367 auto calculateBkg = createChildAlgorithm("CalculatePolynomialBackground");
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");
375}
376
385 std::vector<double> const &peakStarts,
386 std::vector<double> const &peakEnds) {
387 HistogramData::Points wsIndices{ws.getNumberHistograms(), HistogramData::LinearGenerator(0., 1.)};
388 // zeroCounts actually holds the mean frequencies but since its
389 // point data the type has to be Counts.
390 HistogramData::Counts zeroCounts(ws.getNumberHistograms(), 0.);
391 HistogramData::Histogram modelHistogram{wsIndices, zeroCounts};
392 API::MatrixWorkspace_sptr averageWS = DataObjects::create<DataObjects::Workspace2D>(1, modelHistogram);
393 for (size_t i = 0; i < ws.getNumberHistograms(); ++i) {
394 auto const peakStart = peakStarts[i];
395 auto const peakEnd = peakEnds[i];
396 auto histogram = ws.histogram(i);
397 size_t itemCount{0};
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) {
403 continue;
404 }
405 sum += histogramItem.frequency();
406 auto const stdDev = histogramItem.frequencyStandardDeviation();
407 error = std::sqrt(error * error + stdDev * stdDev);
408 ++itemCount;
409 }
410 if (itemCount != 0) {
411 sum /= static_cast<double>(itemCount);
412 error /= static_cast<double>(itemCount);
413 }
414 }
415 return averageWS;
416}
417
418} // namespace Mantid::Algorithms
#define DECLARE_ALGORITHM(classname)
Definition: Algorithm.h:576
double error
Definition: IndexPeaks.cpp:133
#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...
Definition: MultiThreaded.h:94
#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.
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
void progress(double p, const std::string &msg="", double estimatedTime=0.0, int progressPrecision=0)
Sends ProgressNotification.
Definition: Algorithm.cpp:231
bool isDefault(const std::string &name) const
Definition: Algorithm.cpp:2084
void setPropertyValue(const std::string &name, const std::string &value) override
Set the value of a property by string N.B.
Definition: Algorithm.cpp:1975
Column is the base class for columns of TableWorkspace.
Definition: Column.h:35
T & cell(size_t index)
Templated method for returning a value. No type checks are done.
Definition: Column.h:127
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
Definition: ISpectrum.cpp:123
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.
Definition: Progress.h:25
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.
Definition: Component.h:238
bool hasParameter(const std::string &name, bool recursive=true) const override
Returns a boolean indicating if the component has the named parameter.
Definition: Component.cpp:486
Base Instrument Class.
Definition: Instrument.h:47
Support for a property that holds an array of values.
Definition: ArrayProperty.h:28
Exception for when an item is not found in a collection.
Definition: Exception.h:145
IPropertyManager * setProperty(const std::string &name, const T &value)
Templated method to set the value of a PropertyWithValue.
std::shared_ptr< Column > Column_sptr
Definition: Column.h:228
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.
Definition: Exception.h:95
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.
Definition: MultiThreaded.h:22
String constants for algorithm's properties.
STL namespace.
@ Input
An input workspace.
Definition: Property.h:53
@ Output
An output workspace.
Definition: Property.h:54