14#include "MantidHistogramData/Histogram.h"
15#include "MantidHistogramData/Interpolate.h"
22const static std::string FILENAME{
"Filename"};
23const static std::string OUT_WS{
"OutputWorkspace"};
24const static std::string REF_WS{
"WavelengthReference"};
28struct FactorDefinition {
30 std::vector<double> limits;
32 std::vector<double> fitFactors;
36enum class Factor { F1, F2, Phi, P1, P2 };
39Factor factor(
const std::string &l) {
40 const auto id = l.substr(0, 2);
51 throw std::runtime_error(
"Syntax error.");
55const std::array<Factor, 5> factor_list() {
return {{Factor::F1, Factor::F2, Factor::P1, Factor::P2, Factor::Phi}}; }
71 throw std::runtime_error(
"Unknown polarization correction factor tag.");
75std::string cleanse_comments(
const std::string &l) {
76 const auto commentBegin = l.find(
';');
78 return l.substr(0, commentBegin);
82std::string cleanse_whitespace(
const std::string &l) {
85 for (
const auto c : l) {
94bool contains_limits(
const std::string &l) {
return l.find(
"_limits") != std::string::npos; }
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;
104 for (
size_t i = valBegin + 1; i != valEnd; i++) {
108 if (i != valEnd - 1) {
114 v = std::stod(valStr);
115 }
catch (std::exception &) {
116 throw std::runtime_error(
"Syntax error");
118 values.emplace_back(v);
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];
132 const Mantid::HistogramData::Counts c(p.size());
133 return Mantid::HistogramData::Histogram(p, c);
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]);
147std::map<Factor, FactorDefinition> parse(std::istream &in) {
148 std::map<Factor, FactorDefinition> factors;
153 }
catch (std::exception &e) {
154 throw std::runtime_error(std::string(
"Unknown exception: ") + e.what());
156 l = cleanse_whitespace(l);
157 l = cleanse_comments(l);
160 auto &fDef = factors[factor(l)];
161 const auto values = extract_values(l);
162 if (contains_limits(l)) {
163 fDef.limits = values;
165 fDef.fitFactors = values;
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);
177 throw std::runtime_error(
"One of the factors is missing.");
179 const auto &fDef = (*i).second;
180 if (fDef.limits.empty()) {
181 throw std::runtime_error(
"No limits defined for a factor.");
183 if (fDef.fitFactors.empty()) {
184 throw std::runtime_error(
"No fitting information defined for a factor.");
186 if (fDef.limits.size() + 2 != fDef.fitFactors.size()) {
187 throw std::runtime_error(
"Size mismatch between limits and fitting information.");
193void addErrors(Mantid::HistogramData::Histogram &h,
const Factor tag) {
195 const auto f = [tag]() {
205 throw std::logic_error(
"Logic error: unknown efficiency factor tag.");
208 h.mutableE() = (h.y() * f).rawData();
215 ws.
setYUnit(
"Polarization efficiency");
240 return "Loads ILL formatted reflectometry polarization efficiency factors.";
248 "Path to the polarization efficiency file.");
249 const auto refWSValidator = std::make_shared<API::IncreasingAxisValidator>();
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.");
263 HistogramData::Histogram tmplHist{refWS->histogram(0).points()};
265 auto outVertAxis = std::make_unique<API::TextAxis>(5);
266 const auto maxWavelength = tmplHist.x().back();
268 const std::string filename =
getProperty(Prop::FILENAME);
269 std::ifstream in(filename);
271 throw std::runtime_error(
"Couldn't open file " + filename);
273 const auto fittingData = [&in, &filename]() {
275 const auto data = parse(in);
276 definition_map_sanity_check(data);
278 }
catch (std::exception &e) {
279 throw std::runtime_error(
"Error while reading " + filename +
": " + e.what());
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));
301 outWS->replaceAxis(1, std::move(outVertAxis));
303 outWS->setTitle(
"Polarization efficiency factors");
311 std::map<std::string, std::string> issues;
313 if (refWS->getNumberHistograms() == 0) {
314 issues[Prop::REF_WS] =
"The reference workspace does not contain any histograms.";
317 const auto &xs = refWS->x(0);
319 if (xs.front() < 0) {
320 issues[Prop::REF_WS] =
"The reference workspace contains negative X values.";
#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 const std::shared_ptr< Kernel::Unit > & setUnit(const std::string &unitName)
Set the unit on the Axis.
@ Load
allowed here which will be passed to the algorithm
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.
void exec() override
Execute the algorithm.
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.
String constants for algorithm's properties.
std::string to_string(const wide_integer< Bits, Signed > &n)
Describes the direction (within an algorithm) of a Property.
@ Input
An input workspace.
@ Output
An output workspace.