Mantid
Loading...
Searching...
No Matches
CreateFloodWorkspace.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 +
12#include "MantidAPI/Run.h"
18
19#include <algorithm>
20#include <iterator>
21#include <limits>
22
23using namespace Mantid::Kernel;
24using namespace Mantid::API;
25
26namespace {
27namespace Prop {
28std::string const FILENAME("Filename");
29std::string const OUTPUT_WORKSPACE("OutputWorkspace");
30std::string const START_X("StartSpectrum");
31std::string const END_X("EndSpectrum");
32std::string const EXCLUDE("ExcludeSpectra");
33std::string const BACKGROUND("Background");
34std::string const CENTRAL_PIXEL("CentralPixelSpectrum");
35std::string const RANGE_LOWER("RangeLower");
36std::string const RANGE_UPPER("RangeUpper");
37} // namespace Prop
38
39// Too large number makes mantid crash when trying to open a plot
40double const VERY_BIG_VALUE = 1.0e200;
41std::map<std::string, std::string> const funMap{{"Linear", "name=LinearBackground"}, {"Quadratic", "name=Quadratic"}};
42
43} // namespace
44
45namespace Mantid::Algorithms {
46
47// Register the algorithm into the AlgorithmFactory
48DECLARE_ALGORITHM(CreateFloodWorkspace)
49
50const std::string CreateFloodWorkspace::name() const { return "CreateFloodWorkspace"; }
51
52const std::string CreateFloodWorkspace::summary() const {
53 return "Algorithm to create a flood correction workspace for reflectometry "
54 "data reduction.";
55}
56
57int CreateFloodWorkspace::version() const { return 1; }
58
59const std::vector<std::string> CreateFloodWorkspace::seeAlso() const { return {"ReflectometryReductionOneAuto"}; }
60
61const std::string CreateFloodWorkspace::category() const { return "Reflectometry\\ISIS"; }
62
66 const FacilityInfo &defaultFacility = ConfigService::Instance().getFacility();
67 std::vector<std::string> exts = defaultFacility.extensions();
68
69 declareProperty(std::make_unique<MultipleFileProperty>(Prop::FILENAME, exts),
70 "The name of the flood run file(s) to read. Multiple runs "
71 "can be loaded and added together, e.g. INST10+11+12+13.ext");
72
73 declareProperty(Prop::START_X, EMPTY_DBL(), "Start value of the fitting interval");
74 declareProperty(Prop::END_X, EMPTY_DBL(), "End value of the fitting interval");
75
76 declareProperty(std::make_unique<ArrayProperty<double>>(Prop::EXCLUDE), "Spectra to exclude");
77
78 declareProperty(Prop::RANGE_LOWER, EMPTY_DBL(), "The lower integration limit (an X value).");
79
80 declareProperty(Prop::RANGE_UPPER, EMPTY_DBL(), "The upper integration limit (an X value).");
81
82 declareProperty(Prop::CENTRAL_PIXEL, EMPTY_INT(), "A spectrum number of the central pixel.");
83
84 std::vector<std::string> allowedValues;
85 allowedValues.reserve(funMap.size());
86 std::transform(funMap.cbegin(), funMap.cend(), std::back_inserter(allowedValues),
87 [](const auto &i) { return i.first; });
88 auto backgroundValidator = std::make_shared<ListValidator<std::string>>(allowedValues);
89 declareProperty(Prop::BACKGROUND, "Linear", backgroundValidator, "Background function.");
90
91 declareProperty(std::make_unique<WorkspaceProperty<MatrixWorkspace>>(Prop::OUTPUT_WORKSPACE, "", Direction::Output),
92 "The flood correction workspace.");
93}
94
96 std::string fileName = getProperty(Prop::FILENAME);
97 auto alg = createChildAlgorithm("Load", 0, 0.8);
98 alg->setProperty("Filename", fileName);
99 if (alg->existsProperty("LoadMonitors")) {
100 alg->setProperty("LoadMonitors", false);
101 }
102 alg->setProperty("OutputWorkspace", "dummy");
103 alg->execute();
104 Workspace_sptr ws = alg->getProperty("OutputWorkspace");
105 auto input = std::dynamic_pointer_cast<MatrixWorkspace>(ws);
106 if (!input) {
107 throw std::invalid_argument("Loaded files do not produce a single MatrixWorkspace as expected.");
108 }
109 return input;
110}
111
112std::string CreateFloodWorkspace::getBackgroundFunction() { return funMap.at(getPropertyValue(Prop::BACKGROUND)); }
113
115 auto alg = createChildAlgorithm("Integration");
116 alg->setProperty("InputWorkspace", ws);
117 alg->setProperty("OutputWorkspace", "dummy");
118 if (!isDefault(Prop::RANGE_LOWER)) {
119 alg->setProperty("RangeLower", static_cast<double>(getProperty(Prop::RANGE_LOWER)));
120 }
121 if (!isDefault(Prop::RANGE_UPPER)) {
122 alg->setProperty("RangeUpper", static_cast<double>(getProperty(Prop::RANGE_UPPER)));
123 }
124 alg->execute();
125 return alg->getProperty("OutputWorkspace");
126}
127
129 auto alg = createChildAlgorithm("Transpose");
130 alg->setProperty("InputWorkspace", ws);
131 alg->setProperty("OutputWorkspace", "dummy");
132 alg->execute();
133 return alg->getProperty("OutputWorkspace");
134}
135
136bool CreateFloodWorkspace::shouldRemoveBackground() { return isDefault(Prop::CENTRAL_PIXEL); }
137
139
141 return std::find(m_excludedSpectra.begin(), m_excludedSpectra.end(), spec) != m_excludedSpectra.end();
142}
143
145 g_log.information() << "Remove background " << getPropertyValue(Prop::BACKGROUND) << '\n';
146 auto fitWS = transpose(ws);
147 auto const &x = fitWS->x(0);
148
149 // Define the fitting interval
150 double startX = getProperty(Prop::START_X);
151 double endX = getProperty(Prop::END_X);
152 std::vector<double> excludeFromFit;
153 if (isDefault(Prop::START_X)) {
154 startX = x.front();
155 } else {
156 excludeFromFit.emplace_back(x.front());
157 excludeFromFit.emplace_back(startX);
158 }
159 if (isDefault(Prop::END_X)) {
160 endX = x.back();
161 } else {
162 excludeFromFit.emplace_back(endX);
163 excludeFromFit.emplace_back(x.back());
164 }
165
166 // Exclude any bad detectors.
167 for (auto i : m_excludedSpectra) {
168 excludeFromFit.emplace_back(i);
169 excludeFromFit.emplace_back(i);
170 }
171
172 std::string const function = getBackgroundFunction();
173
174 // Fit the data to determine unwanted background
175 auto alg = createChildAlgorithm("Fit", 0.9, 0.99);
176 alg->setProperty("Function", function);
177 alg->setProperty("InputWorkspace", fitWS);
178 alg->setProperty("WorkspaceIndex", 0);
179 if (!excludeFromFit.empty()) {
180 alg->setProperty("Exclude", excludeFromFit);
181 }
182 alg->setProperty("Output", "fit");
183 alg->execute();
184
185 IFunction_sptr func = alg->getProperty("Function");
186 g_log.information() << "Background function parameters:\n";
187 for (size_t i = 0; i < func->nParams(); ++i) {
188 g_log.information() << " " << func->parameterName(i) << ": " << func->getParameter(i) << '\n';
189 }
190
191 // Divide the workspace by the fitted curve to remove the background
192 // and scale to values around 1
193 MatrixWorkspace_sptr bkgWS = alg->getProperty("OutputWorkspace");
194 auto const &bkg = bkgWS->y(1);
195 auto const nHisto = static_cast<int>(ws->getNumberHistograms());
197 for (int i = 0; i < nHisto; ++i) {
199 auto const xVal = x[i];
200 if (isExcludedSpectrum(xVal)) {
201 ws->mutableY(i)[0] = VERY_BIG_VALUE;
202 ws->mutableE(i)[0] = 0.0;
203 } else if (xVal >= startX && xVal <= endX) {
204 auto const background = bkg[i];
205 if (background <= 0.0) {
206 throw std::runtime_error("Background is expected to be positive, found value " + std::to_string(background) +
207 " at spectrum with workspace index " + std::to_string(i));
208 }
209 ws->mutableY(i)[0] /= background;
210 ws->mutableE(i)[0] /= background;
211 } else {
212 ws->mutableY(i)[0] = 1.0;
213 ws->mutableE(i)[0] = 0.0;
214 }
216 }
218
219 // Remove the logs
220 ws->setSharedRun(make_cow<Run>());
221
222 return ws;
223}
224
226 int const centralSpectrum = getProperty(Prop::CENTRAL_PIXEL);
227 auto const nHisto = static_cast<int>(ws->getNumberHistograms());
228 if (centralSpectrum >= nHisto) {
229 throw std::invalid_argument("Spectrum index " + std::to_string(centralSpectrum) + " passed to property " +
230 Prop::CENTRAL_PIXEL + " is outside the range 0-" + std::to_string(nHisto - 1));
231 }
232 auto const spectraMap = ws->getSpectrumToWorkspaceIndexMap();
233 auto const centralIndex = spectraMap.at(centralSpectrum);
234 auto const scaleFactor = ws->y(centralIndex).front();
235 g_log.information() << "Scale to central pixel, factor = " << scaleFactor << '\n';
236 if (scaleFactor <= 0.0) {
237 throw std::runtime_error("Scale factor muhst be > 0, found " + std::to_string(scaleFactor));
238 }
239 auto const axis = ws->getAxis(1);
240 auto const sa = dynamic_cast<const SpectraAxis *>(axis);
241 double const startX = isDefault(Prop::START_X) ? sa->getMin() : getProperty(Prop::START_X);
242 double const endX = isDefault(Prop::END_X) ? sa->getMax() : getProperty(Prop::END_X);
244 for (int i = 0; i < nHisto; ++i) {
246 auto const spec = ws->getSpectrum(i).getSpectrumNo();
247 if (isExcludedSpectrum(spec)) {
248 ws->mutableY(i)[0] = VERY_BIG_VALUE;
249 ws->mutableE(i)[0] = 0.0;
250 } else if (spec >= startX && spec <= endX) {
251 ws->mutableY(i)[0] /= scaleFactor;
252 ws->mutableE(i)[0] /= scaleFactor;
253 } else {
254 ws->mutableY(i)[0] = 1.0;
255 ws->mutableE(i)[0] = 0.0;
256 }
258 }
260 return ws;
261}
262
266 progress(0.0);
267 auto ws = getInputWorkspace();
268 ws = integrate(ws);
269 progress(0.9);
272 ws = removeBackground(ws);
273 } else {
274 ws = scaleToCentralPixel(ws);
275 }
276 progress(1.0);
277 setProperty(Prop::OUTPUT_WORKSPACE, ws);
278}
279
280} // namespace Mantid::Algorithms
#define DECLARE_ALGORITHM(classname)
Definition: Algorithm.h:576
#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
std::string getPropertyValue(const std::string &name) const override
Get the value of a property as a string.
Definition: Algorithm.cpp:2026
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
bool isDefault(const std::string &name) const
Definition: Algorithm.cpp:2084
Class to represent the spectra axis of a workspace.
Definition: SpectraAxis.h:31
A property class for workspaces.
Algorithm to create a flood correction workspace for reflectometry data reduction.
API::MatrixWorkspace_sptr removeBackground(API::MatrixWorkspace_sptr ws)
API::MatrixWorkspace_sptr transpose(const API::MatrixWorkspace_sptr &ws)
void exec() override
Execute the algorithm.
const std::vector< std::string > seeAlso() const override
Function to return all of the seeAlso algorithms related to this algorithm.
const std::string summary() const override
function returns a summary message that will be displayed in the default GUI, and in the help.
int version() const override
function to return a version of the algorithm, must be overridden in all algorithms
API::MatrixWorkspace_sptr integrate(const API::MatrixWorkspace_sptr &ws)
void init() override
Initialize the algorithm's properties.
API::MatrixWorkspace_sptr scaleToCentralPixel(API::MatrixWorkspace_sptr ws)
const std::string category() const override
function to return a category of the algorithm.
Support for a property that holds an array of values.
Definition: ArrayProperty.h:28
A class that holds information about a facility.
Definition: FacilityInfo.h:36
const std::vector< std::string > extensions() const
Returns a list of file extensions.
Definition: FacilityInfo.h:48
IPropertyManager * setProperty(const std::string &name, const T &value)
Templated method to set the value of a PropertyWithValue.
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< Workspace > Workspace_sptr
shared pointer to Mantid::API::Workspace
Definition: Workspace_fwd.h:20
std::shared_ptr< IFunction > IFunction_sptr
shared pointer to the function base class
Definition: IFunction.h:732
std::shared_ptr< MatrixWorkspace > MatrixWorkspace_sptr
shared pointer to the matrix workspace base class
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
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
String constants for algorithm's properties.
STL namespace.
std::string to_string(const wide_integer< Bits, Signed > &n)
@ Output
An output workspace.
Definition: Property.h:54