Mantid
Loading...
Searching...
No Matches
LoadILLPolarizationFactors.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/TextAxis.h"
14#include "MantidHistogramData/Histogram.h"
15#include "MantidHistogramData/Interpolate.h"
16
17#include <fstream>
18
19namespace {
21namespace Prop {
22const static std::string FILENAME{"Filename"};
23const static std::string OUT_WS{"OutputWorkspace"};
24const static std::string REF_WS{"WavelengthReference"};
25} // namespace Prop
26
28struct FactorDefinition {
30 std::vector<double> limits;
32 std::vector<double> fitFactors;
33};
34
36enum class Factor { F1, F2, Phi, P1, P2 };
37
39Factor factor(const std::string &l) {
40 const auto id = l.substr(0, 2);
41 if (id == "F1")
42 return Factor::F1;
43 if (id == "F2")
44 return Factor::F2;
45 if (id == "Ph")
46 return Factor::Phi;
47 if (id == "P1")
48 return Factor::P1;
49 if (id == "P2")
50 return Factor::P2;
51 throw std::runtime_error("Syntax error.");
52}
53
55const std::array<Factor, 5> factor_list() { return {{Factor::F1, Factor::F2, Factor::P1, Factor::P2, Factor::Phi}}; }
56
58std::string to_string(const Factor f) {
59 switch (f) {
60 case Factor::F1:
61 return "F1";
62 case Factor::F2:
63 return "F2";
64 case Factor::P1:
65 return "P1";
66 case Factor::P2:
67 return "P2";
68 case Factor::Phi:
69 return "Phi";
70 }
71 throw std::runtime_error("Unknown polarization correction factor tag.");
72}
73
75std::string cleanse_comments(const std::string &l) {
76 const auto commentBegin = l.find(';');
77 // commentBegin == npos is OK.
78 return l.substr(0, commentBegin);
79}
80
82std::string cleanse_whitespace(const std::string &l) {
83 std::string s;
84 s.reserve(l.size());
85 for (const auto c : l) {
86 if (!isblank(c)) {
87 s.push_back(c);
88 }
89 }
90 return s;
91}
92
94bool contains_limits(const std::string &l) { return l.find("_limits") != std::string::npos; }
95
97std::vector<double> extract_values(const std::string &l) {
98 const auto valBegin = l.find('[');
99 const auto valEnd = l.find(']');
100 if (valBegin == std::string::npos || valEnd == std::string::npos)
101 throw std::runtime_error("Syntax error.");
102 std::vector<double> values;
103 std::string valStr;
104 for (size_t i = valBegin + 1; i != valEnd; i++) {
105 const auto c = l[i];
106 if (c != ',') {
107 valStr.push_back(c);
108 if (i != valEnd - 1) {
109 continue;
110 }
111 }
112 double v;
113 try {
114 v = std::stod(valStr);
115 } catch (std::exception &) {
116 throw std::runtime_error("Syntax error");
117 }
118 values.emplace_back(v);
119 valStr.clear();
120 }
121 return values;
122}
123
125Mantid::HistogramData::Histogram make_histogram(const std::vector<double> &points, const double maxWavelength) {
126 Mantid::HistogramData::Points p(points.size() + 2);
127 p.mutableRawData().front() = 0;
128 p.mutableRawData().back() = maxWavelength > points.back() ? maxWavelength : 2 * points.back();
129 for (size_t i = 1; i != p.size() - 1; ++i) {
130 p.mutableData()[i] = points[i - 1];
131 }
132 const Mantid::HistogramData::Counts c(p.size());
133 return Mantid::HistogramData::Histogram(p, c);
134}
135
137void calculate_factors_in_place(Mantid::HistogramData::Histogram &h, const std::vector<double> &piecewiseFactors) {
138 const auto &xs = h.x();
139 auto &ys = h.mutableY();
140 ys[0] = piecewiseFactors.front();
141 for (size_t i = 1; i != h.size(); ++i) {
142 ys[i] = ys[i - 1] + piecewiseFactors[i] * (xs[i] - xs[i - 1]);
143 }
144}
145
147std::map<Factor, FactorDefinition> parse(std::istream &in) {
148 std::map<Factor, FactorDefinition> factors;
149 std::string l;
150 while (!in.eof()) {
151 try {
152 std::getline(in, l);
153 } catch (std::exception &e) {
154 throw std::runtime_error(std::string("Unknown exception: ") + e.what());
155 }
156 l = cleanse_whitespace(l);
157 l = cleanse_comments(l);
158 if (l.empty())
159 continue;
160 auto &fDef = factors[factor(l)];
161 const auto values = extract_values(l);
162 if (contains_limits(l)) {
163 fDef.limits = values;
164 } else {
165 fDef.fitFactors = values;
166 }
167 }
168 return factors;
169}
170
172void definition_map_sanity_check(const std::map<Factor, FactorDefinition> &m) {
173 const auto factors = factor_list();
174 for (const auto f : factors) {
175 const auto i = m.find(f);
176 if (i == m.end()) {
177 throw std::runtime_error("One of the factors is missing.");
178 }
179 const auto &fDef = (*i).second;
180 if (fDef.limits.empty()) {
181 throw std::runtime_error("No limits defined for a factor.");
182 }
183 if (fDef.fitFactors.empty()) {
184 throw std::runtime_error("No fitting information defined for a factor.");
185 }
186 if (fDef.limits.size() + 2 != fDef.fitFactors.size()) {
187 throw std::runtime_error("Size mismatch between limits and fitting information.");
188 }
189 }
190}
191
193void addErrors(Mantid::HistogramData::Histogram &h, const Factor tag) {
194 // The error estimates are taken from the LAMP/COSMOS software.
195 const auto f = [tag]() {
196 switch (tag) {
197 case Factor::F1:
198 case Factor::F2:
199 return 1. / 3000.;
200 case Factor::P1:
201 case Factor::P2:
202 case Factor::Phi:
203 return 1. / 500.;
204 default:
205 throw std::logic_error("Logic error: unknown efficiency factor tag.");
206 }
207 }();
208 h.mutableE() = (h.y() * f).rawData();
209}
210
212void setUnits(Mantid::API::MatrixWorkspace &ws) {
213 auto xAxis = ws.getAxis(0);
214 xAxis->setUnit("Wavelength");
215 ws.setYUnit("Polarization efficiency");
216}
217} // namespace
218
219namespace Mantid::DataHandling {
220
223
224// Register the algorithm into the AlgorithmFactory
225DECLARE_ALGORITHM(LoadILLPolarizationFactors)
226
227//----------------------------------------------------------------------------------------------
228
229
230const std::string LoadILLPolarizationFactors::name() const { return "LoadILLPolarizationFactors"; }
231
233int LoadILLPolarizationFactors::version() const { return 1; }
234
236const std::string LoadILLPolarizationFactors::category() const { return "DataHandling\\Text;ILL\\Reflectometry"; }
237
239const std::string LoadILLPolarizationFactors::summary() const {
240 return "Loads ILL formatted reflectometry polarization efficiency factors.";
241}
242
243//----------------------------------------------------------------------------------------------
247 declareProperty(std::make_unique<API::FileProperty>(Prop::FILENAME, "", API::FileProperty::Load),
248 "Path to the polarization efficiency file.");
249 const auto refWSValidator = std::make_shared<API::IncreasingAxisValidator>();
251 std::make_unique<WorkspaceProperty<API::MatrixWorkspace>>(Prop::OUT_WS, "", Direction::Output, refWSValidator),
252 "An output workspace containing the efficiencies at the "
253 "reference workspace's wavelength points.");
255 "A reference workspace to get the wavelength axis from.");
256}
257
258//----------------------------------------------------------------------------------------------
262 API::MatrixWorkspace_const_sptr refWS = getProperty(Prop::REF_WS);
263 HistogramData::Histogram tmplHist{refWS->histogram(0).points()};
264 API::MatrixWorkspace_sptr outWS = DataObjects::create<DataObjects::Workspace2D>(5, tmplHist);
265 auto outVertAxis = std::make_unique<API::TextAxis>(5);
266 const auto maxWavelength = tmplHist.x().back();
267
268 const std::string filename = getProperty(Prop::FILENAME);
269 std::ifstream in(filename);
270 if (in.bad()) {
271 throw std::runtime_error("Couldn't open file " + filename);
272 }
273 const auto fittingData = [&in, &filename]() {
274 try {
275 const auto data = parse(in);
276 definition_map_sanity_check(data);
277 return data;
278 } catch (std::exception &e) {
279 throw std::runtime_error("Error while reading " + filename + ": " + e.what());
280 }
281 }();
282 const auto factorTags = factor_list();
284 for (int i = 0; i < static_cast<int>(factorTags.size()); ++i) {
286 const auto tag = factorTags[i];
287 const auto &fDef = fittingData.at(tag);
288 auto source = make_histogram(fDef.limits, maxWavelength);
289 calculate_factors_in_place(source, fDef.fitFactors);
290 auto target = outWS->histogram(i);
291 HistogramData::interpolateLinearInplace(source, target);
292 addErrors(target, tag);
293 HistogramData::Histogram outH{refWS->histogram(0)};
294 outH.setSharedY(target.sharedY());
295 outH.setSharedE(target.sharedE());
296 outWS->setHistogram(i, outH);
297 outVertAxis->setLabel(i, to_string(tag));
299 }
301 outWS->replaceAxis(1, std::move(outVertAxis));
302 setUnits(*outWS);
303 outWS->setTitle("Polarization efficiency factors");
304 setProperty(Prop::OUT_WS, outWS);
305}
306
310std::map<std::string, std::string> LoadILLPolarizationFactors::validateInputs() {
311 std::map<std::string, std::string> issues;
312 API::MatrixWorkspace_const_sptr refWS = getProperty(Prop::REF_WS);
313 if (refWS->getNumberHistograms() == 0) {
314 issues[Prop::REF_WS] = "The reference workspace does not contain any histograms.";
315 return issues;
316 }
317 const auto &xs = refWS->x(0);
318 // A validator should have checked that xs is ordered.
319 if (xs.front() < 0) {
320 issues[Prop::REF_WS] = "The reference workspace contains negative X values.";
321 }
322 return issues;
323}
324
325} // namespace Mantid::DataHandling
#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
TypedValue getProperty(const std::string &name) const override
Get the value of a property.
Definition: Algorithm.cpp:2076
virtual const std::shared_ptr< Kernel::Unit > & setUnit(const std::string &unitName)
Set the unit on the Axis.
Definition: Axis.cpp:39
@ Load
allowed here which will be passed to the algorithm
Definition: FileProperty.h:52
Base MatrixWorkspace Abstract Class.
virtual Axis * getAxis(const std::size_t &axisIndex) const
Get a non owning pointer to a workspace axis.
void setYUnit(const std::string &newUnit)
Sets a new unit for the data (Y axis) in the workspace.
A property class for workspaces.
LoadILLPolarizationFactors : Load reflectometry polarization efficiency correction factors from disk.
int version() const override
Algorithm's version for identification.
const std::string category() const override
Algorithm's category for identification.
void init() override
Initialize the algorithm's properties.
const std::string summary() const override
Algorithm's summary for use in the GUI and help.
std::map< std::string, std::string > validateInputs() override
Validates the algorithm's inputs.
IPropertyManager * setProperty(const std::string &name, const T &value)
Templated method to set the value of a PropertyWithValue.
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
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.
std::string to_string(const wide_integer< Bits, Signed > &n)
Describes the direction (within an algorithm) of a Property.
Definition: Property.h:50
@ Input
An input workspace.
Definition: Property.h:53
@ Output
An output workspace.
Definition: Property.h:54