Mantid
Loading...
Searching...
No Matches
PolarizationCorrectionWildes.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 +
7#include "MantidAlgorithms/PolarizationCorrectionWildes.h"
8
9#include "MantidAPI/ADSValidator.h"
10#include "MantidAPI/Axis.h"
11#include "MantidAPI/WorkspaceGroup.h"
12#include "MantidAPI/WorkspaceHistory.h"
13#include "MantidAlgorithms/PolarizationCorrections/PolarizationCorrectionsHelpers.h"
14#include "MantidDataObjects/Workspace2D.h"
15#include "MantidDataObjects/WorkspaceCreation.h"
16#include "MantidKernel/ArrayProperty.h"
17#include "MantidKernel/ListValidator.h"
18#include "MantidKernel/SpinStateHelpers.h"
19#include "MantidKernel/SpinStateValidator.h"
20#include "MantidKernel/StringTokenizer.h"
21
22#include <Eigen/Dense>
23#include <boost/math/special_functions/pow.hpp>
24
25namespace {
27namespace Prop {
28static const std::string FLIPPERS{"Flippers"};
29static const std::string SPIN_STATES{"SpinStates"};
30static const std::string EFFICIENCIES{"Efficiencies"};
31static const std::string INPUT_WS{"InputWorkspaces"};
32static const std::string OUTPUT_WS{"OutputWorkspace"};
33static const std::string ADD_SPIN_STATE_LOG{"AddSpinStateToLog"};
34} // namespace Prop
35
41void checkInputExists(const Mantid::API::MatrixWorkspace_sptr &ws, const std::string &tag) {
42 if (!ws) {
43 throw std::runtime_error("A workspace designated as " + tag + " is missing in inputs.");
44 }
45}
46
68void fourInputsCorrectedAndErrors(Eigen::Vector4d &corrected, Eigen::Vector4d &errors, const double ppy,
69 const double ppyE, const double pmy, const double pmyE, const double mpy,
70 const double mpyE, const double mmy, const double mmyE, const double f1,
71 const double f1E, const double f2, const double f2E, const double p1,
72 const double p1E, const double p2, const double p2E) {
73 using namespace boost::math;
74 // Note that f1 and f2 correspond to 1-F1 and 1-F2 in [Wildes, 1999].
75 // These are inverted forms of the efficiency matrices.
76 const auto diag1 = 1. / f1;
77 const auto off1 = (f1 - 1.) / f1;
78 Eigen::Matrix4d F1m;
79 F1m << 1., 0., 0., 0., 0., 1., 0., 0., off1, 0., diag1, 0., 0., off1, 0., diag1;
80
81 const auto diag2 = 1. / f2;
82 const auto off2 = (f2 - 1.) / f2;
83 Eigen::Matrix4d F2m;
84 F2m << 1., 0., 0., 0., off2, diag2, 0., 0., 0., 0., 1., 0., 0., 0., off2, diag2;
85 const auto diag3 = p1 / (2. * p1 - 1);
86 const auto off3 = (p1 - 1.) / (2. * p1 - 1.);
87 Eigen::Matrix4d P1m;
88 P1m << diag3, 0, off3, 0, 0, diag3, 0, off3, off3, 0, diag3, 0, 0, off3, 0, diag3;
89 const auto diag4 = p2 / (2. * p2 - 1.);
90 const auto off4 = (p2 - 1.) / (2. * p2 - 1.);
91 Eigen::Matrix4d P2m;
92 P2m << diag4, off4, 0., 0., off4, diag4, 0., 0., 0., 0., diag4, off4, 0., 0., off4, diag4;
93 const Eigen::Vector4d intensities(ppy, pmy, mpy, mmy);
94 const auto FProduct = F2m * F1m;
95 const auto PProduct = P2m * P1m;
96 const auto PFProduct = PProduct * FProduct;
97 corrected = PFProduct * intensities;
98 // The error matrices here are element-wise algebraic derivatives of
99 // the matrices above, multiplied by the error.
100 const auto elemE1 = -1. / pow<2>(f1) * f1E;
101 Eigen::Matrix4d F1Em;
102 F1Em << 0., 0., 0., 0., 0., 0., 0., 0., -elemE1, 0., elemE1, 0., 0., -elemE1, 0., elemE1;
103 const auto elemE2 = -1. / pow<2>(f2) * f2E;
104 Eigen::Matrix4d F2Em;
105 F2Em << 0., 0., 0., 0., -elemE2, elemE2, 0., 0., 0., 0., 0., 0., 0., 0., -elemE2, elemE2;
106 const auto elemE3 = 1. / pow<2>(2. * p1 - 1.) * p1E;
107 Eigen::Matrix4d P1Em;
108 P1Em << elemE3, 0., -elemE3, 0., 0., elemE3, 0., -elemE3, -elemE3, 0., elemE3, 0., 0., -elemE3, 0., elemE3;
109 const auto elemE4 = 1. / pow<2>(2. * p2 - 1.) * p2E;
110 Eigen::Matrix4d P2Em;
111 P2Em << elemE4, -elemE4, 0., 0., -elemE4, elemE4, 0., 0., 0., 0., elemE4, -elemE4, 0., 0., -elemE4, elemE4;
112 const Eigen::Vector4d yErrors(ppyE, pmyE, mpyE, mmyE);
113 const auto e1 = (P2Em * P1m * FProduct * intensities).array();
114 const auto e2 = (P2m * P1Em * FProduct * intensities).array();
115 const auto e3 = (PProduct * F2Em * F1m * intensities).array();
116 const auto e4 = (PProduct * F2m * F1Em * intensities).array();
117 const auto sqPFProduct = (PFProduct.array() * PFProduct.array()).matrix();
118 const auto sqErrors = (yErrors.array() * yErrors.array()).matrix();
119 const auto e5 = (sqPFProduct * sqErrors).array();
120 errors = (e1 * e1 + e2 * e2 + e3 * e3 + e4 * e4 + e5).sqrt();
121}
122
139double twoInputsErrorEstimate01(const double i00, const double e00, const double i11, const double e11, const double p1,
140 const double p1E, const double p2, const double p2E, const double f1, const double f1E,
141 const double f2, const double f2E) {
142 using namespace boost::math;
143 // Derivatives of the equation which solves the I01 intensities
144 // with respect to i00, i11, f1, etc.
145 const auto pmdi00 =
146 -((f1 * (-1. + 2. * p1) * (-f2 * pow<2>(1. - 2. * p2) + pow<2>(f2) * pow<2>(1. - 2. * p2) + (-1. + p2) * p2)) /
147 (f2 * p1 * (-1. + p1 + 2. * p2 - 2. * p1 * p2) +
148 f1 * (-1. + 2. * p1) * ((1. - p2) * p2 + f2 * (-1. + p1 + p2) * (-1. + 2. * p2))));
149 const auto pmdi11 = (f2 * p1 * (-1. + p1 + 2. * p2 - 2. * p1 * p2)) /
150 (f2 * p1 * (-1. + p1 + 2. * p2 - 2. * p1 * p2) +
151 f1 * (-1. + 2. * p1) * ((1. - p2) * p2 + f2 * (-1. + p1 + p2) * (-1. + 2. * p2)));
152 const auto pmdf1 =
153 -(((-1. + 2. * p1) * ((1. - p2) * p2 + f2 * (-1. + p1 + p2) * (-1. + 2. * p2)) *
154 (f2 * i11 * p1 * (-1. + p1 + 2. * p2 - 2. * p1 * p2) -
155 f1 * i00 * (-1. + 2. * p1) *
156 (-f2 * pow<2>(1. - 2. * p2) + pow<2>(f2) * pow<2>(1. - 2. * p2) + (-1. + p2) * p2))) /
157 pow<2>(f2 * p1 * (-1. + p1 + 2. * p2 - 2. * p1 * p2) +
158 f1 * (-1. + 2. * p1) * ((1. - p2) * p2 + f2 * (-1. + p1 + p2) * (-1. + 2. * p2)))) -
159 (i00 * (-1. + 2. * p1) * (-f2 * pow<2>(1. - 2. * p2) + pow<2>(f2) * pow<2>(1. - 2. * p2) + (-1. + p2) * p2)) /
160 (f2 * p1 * (-1. + p1 + 2. * p2 - 2. * p1 * p2) +
161 f1 * (-1. + 2. * p1) * ((1. - p2) * p2 + f2 * (-1. + p1 + p2) * (-1. + 2. * p2)));
162 const auto pmdf2 =
163 -(((f1 * (-1. + 2. * p1) * (-1. + p1 + p2) * (-1 + 2 * p2) + p1 * (-1. + p1 + 2. * p2 - 2. * p1 * p2)) *
164 (f2 * i11 * p1 * (-1. + p1 + 2. * p2 - 2. * p1 * p2) -
165 f1 * i00 * (-1. + 2. * p1) *
166 (-f2 * pow<2>(1. - 2. * p2) + pow<2>(f2) * pow<2>(1. - 2. * p2) + (-1. + p2) * p2))) /
167 pow<2>(f2 * p1 * (-1. + p1 + 2. * p2 - 2. * p1 * p2) +
168 f1 * (-1. + 2. * p1) * ((1. - p2) * p2 + f2 * (-1. + p1 + p2) * (-1. + 2. * p2)))) +
169 (-f1 * i00 * (-1. + 2. * p1) * (-pow<2>(1. - 2. * p2) + 2 * f2 * pow<2>(1. - 2. * p2)) +
170 i11 * p1 * (-1. + p1 + 2. * p2 - 2. * p1 * p2)) /
171 (f2 * p1 * (-1. + p1 + 2. * p2 - 2. * p1 * p2) +
172 f1 * (-1. + 2. * p1) * ((1. - p2) * p2 + f2 * (-1. + p1 + p2) * (-1. + 2. * p2)));
173 const auto pmdp1 =
174 -(((f2 * i11 * p1 * (-1. + p1 + 2. * p2 - 2. * p1 * p2) -
175 f1 * i00 * (-1. + 2. * p1) *
176 (-f2 * pow<2>(1. - 2. * p2) + pow<2>(f2) * pow<2>(1. - 2. * p2) + (-1. + p2) * p2)) *
177 (f2 * p1 * (1. - 2. * p2) + f1 * f2 * (-1. + 2. * p1) * (-1. + 2. * p2) +
178 f2 * (-1. + p1 + 2. * p2 - 2. * p1 * p2) +
179 2. * f1 * ((1. - p2) * p2 + f2 * (-1. + p1 + p2) * (-1. + 2. * p2)))) /
180 pow<2>(f2 * p1 * (-1. + p1 + 2. * p2 - 2. * p1 * p2) +
181 f1 * (-1. + 2. * p1) * ((1. - p2) * p2 + f2 * (-1. + p1 + p2) * (-1. + 2. * p2)))) +
182 (f2 * i11 * p1 * (1. - 2. * p2) + f2 * i11 * (-1. + p1 + 2. * p2 - 2. * p1 * p2) -
183 2. * f1 * i00 * (-f2 * pow<2>(1. - 2. * p2) + pow<2>(f2) * pow<2>(1. - 2. * p2) + (-1. + p2) * p2)) /
184 (f2 * p1 * (-1. + p1 + 2. * p2 - 2. * p1 * p2) +
185 f1 * (-1. + 2. * p1) * ((1. - p2) * p2 + f2 * (-1. + p1 + p2) * (-1. + 2. * p2)));
186 const auto pmdp2 =
187 -(((f2 * i11 * p1 * (-1. + p1 + 2. * p2 - 2. * p1 * p2) -
188 f1 * i00 * (-1. + 2. * p1) *
189 (-f2 * pow<2>(1. - 2. * p2) + pow<2>(f2) * pow<2>(1. - 2. * p2) + (-1. + p2) * p2)) *
190 (f2 * (2. - 2. * p1) * p1 +
191 f1 * (-1. + 2. * p1) * (1. - 2. * p2 + 2. * f2 * (-1. + p1 + p2) + f2 * (-1. + 2. * p2)))) /
192 pow<2>(f2 * p1 * (-1. + p1 + 2. * p2 - 2. * p1 * p2) +
193 f1 * (-1. + 2. * p1) * ((1. - p2) * p2 + f2 * (-1. + p1 + p2) * (-1. + 2. * p2)))) +
194 (f2 * i11 * (2. - 2. * p1) * p1 -
195 f1 * i00 * (-1. + 2. * p1) * (-1. + 4. * f2 * (1. - 2. * p2) - 4. * pow<2>(f2) * (1. - 2. * p2) + 2. * p2)) /
196 (f2 * p1 * (-1. + p1 + 2. * p2 - 2. * p1 * p2) +
197 f1 * (-1. + 2. * p1) * ((1. - p2) * p2 + f2 * (-1. + p1 + p2) * (-1. + 2. * p2)));
198 // Estimate the error components using linearized extrapolation,
199 // sum in squares.
200 const auto e01_I00 = pow<2>(pmdi00 * e00);
201 const auto e01_I11 = pow<2>(pmdi11 * e11);
202 const auto e01_F1 = pow<2>(pmdf1 * f1E);
203 const auto e01_F2 = pow<2>(pmdf2 * f2E);
204 const auto e01_P1 = pow<2>(pmdp1 * p1E);
205 const auto e01_P2 = pow<2>(pmdp2 * p2E);
206 return std::sqrt(e01_I00 + e01_I11 + e01_F1 + e01_F2 + e01_P1 + e01_P2);
207}
208
225double twoInputsErrorEstimate10(const double i00, const double e00, const double i11, const double e11, const double p1,
226 const double p1E, const double p2, const double p2E, const double f1, const double f1E,
227 const double f2, const double f2E) {
228 using namespace boost::math;
229 // Derivatives of the equation which solves the I10 intensities
230 // with respect to i00, i11, f1, etc.
231 const auto a = -1. + p1 + 2. * p2 - 2. * p1 * p2;
232 const auto b = -1. + 2. * p1;
233 const auto c = -1. + 2. * p2;
234 const auto d = -1. + p2;
235 const auto mpdi00 = (-pow<2>(f1) * f2 * pow<2>(b) * c + f1 * f2 * pow<2>(b) * c + f2 * p1 * a) /
236 (f2 * p1 * a + f1 * b * (-d * p2 + f2 * (p1 + d) * c));
237 const auto mpdi11 = -((f1 * b * d * p2) / (f2 * p1 * a + f1 * b * (-d * p2 + f2 * (p1 + d) * c)));
238 const auto mpdf1 =
239 -(((-1. + 2. * p1) * ((1. - p2) * p2 + f2 * (-1. + p1 + p2) * (-1. + 2. * p2)) *
240 (-pow<2>(f1) * f2 * i00 * pow<2>(1. - 2. * p1) * (-1. + 2. * p2) +
241 f2 * i00 * p1 * (-1. + p1 + 2. * p2 - 2. * p1 * p2) +
242 f1 * (-1. + 2. * p1) * (-i11 * (-1. + p2) * p2 + f2 * i00 * (-1. + 2. * p1) * (-1. + 2. * p2)))) /
243 pow<2>(f2 * p1 * (-1. + p1 + 2. * p2 - 2. * p1 * p2) +
244 f1 * (-1. + 2. * p1) * ((1. - p2) * p2 + f2 * (-1. + p1 + p2) * (-1. + 2. * p2)))) +
245 (-2. * f1 * f2 * i00 * pow<2>(1. - 2. * p1) * (-1. + 2. * p2) +
246 (-1. + 2. * p1) * (-i11 * (-1. + p2) * p2 + f2 * i00 * (-1. + 2. * p1) * (-1. + 2. * p2))) /
247 (f2 * p1 * (-1. + p1 + 2. * p2 - 2. * p1 * p2) +
248 f1 * (-1. + 2. * p1) * ((1. - p2) * p2 + f2 * (-1. + p1 + p2) * (-1. + 2. * p2)));
249 const auto mpdf2 = -(((f1 * b * (p1 + d) * c + p1 * a) * (-pow<2>(f1) * f2 * i00 * pow<2>(b) * c + f2 * i00 * p1 * a +
250 f1 * b * (-i11 * d * p2 + f2 * i00 * b * c))) /
251 pow<2>(f2 * p1 * a + f1 * b * (-d * p2 + f2 * (p1 + d) * c))) +
252 (-pow<2>(f1) * i00 * pow<2>(b) * c + f1 * i00 * pow<2>(b) * c + i00 * p1 * a) /
253 (f2 * p1 * a + f1 * b * (-d * p2 + f2 * (p1 + d) * c));
254 const auto mpdp1 =
255 -(((-pow<2>(f1) * f2 * i00 * pow<2>(b) * c + f2 * i00 * p1 * a + f1 * b * (-i11 * d * p2 + f2 * i00 * b * c)) *
256 (f2 * p1 * -c + f1 * f2 * b * c + f2 * a + 2. * f1 * (-d * p2 + f2 * (p1 + d) * c))) /
257 pow<2>(f2 * p1 * a + f1 * b * (-d * p2 + f2 * (p1 + d) * c))) +
258 (f2 * i00 * p1 * -c + 4. * pow<2>(f1) * f2 * i00 * -b * c + 2. * f1 * f2 * i00 * b * c + f2 * i00 * a +
259 2. * f1 * (-i11 * d * p2 + f2 * i00 * b * c)) /
260 (f2 * p1 * a + f1 * b * (-d * p2 + f2 * (p1 + d) * c));
261 const auto mpdp2 =
262 -(((f2 * (2. - 2. * p1) * p1 + f1 * b * (1. - 2. * p2 + 2. * f2 * (p1 + d) + f2 * c)) *
263 (-pow<2>(f1) * f2 * i00 * pow<2>(b) * c + f2 * i00 * p1 * a + f1 * b * (-i11 * d * p2 + f2 * i00 * b * c))) /
264 pow<2>(f2 * p1 * a + f1 * b * (-d * p2 + f2 * (p1 + d) * c))) +
265 (-2. * pow<2>(f1) * f2 * i00 * pow<2>(b) + f2 * i00 * (2. - 2. * p1) * p1 +
266 f1 * b * (2. * f2 * i00 * b - i11 * d - i11 * p2)) /
267 (f2 * p1 * a + f1 * b * (-d * p2 + f2 * (p1 + d) * c));
268 // Estimate the error components using linearized extrapolation,
269 // sum in squares.
270 const auto e10_I00 = pow<2>(mpdi00 * e00);
271 const auto e10_I11 = pow<2>(mpdi11 * e11);
272 const auto e10_F1 = pow<2>(mpdf1 * f1E);
273 const auto e10_F2 = pow<2>(mpdf2 * f2E);
274 const auto e10_P1 = pow<2>(mpdp1 * p1E);
275 const auto e10_P2 = pow<2>(mpdp2 * p2E);
276 return std::sqrt(e10_I00 + e10_I11 + e10_F1 + e10_F2 + e10_P1 + e10_P2);
277}
278
279Mantid::API::MatrixWorkspace_sptr createWorkspaceWithHistory(const Mantid::API::MatrixWorkspace_const_sptr &inputWS) {
280 Mantid::API::MatrixWorkspace_sptr outputWS = Mantid::DataObjects::create<Mantid::DataObjects::Workspace2D>(*inputWS);
281 outputWS->history().addHistory(inputWS->getHistory());
282 return outputWS;
283}
284} // namespace
285
286namespace Mantid::Algorithms {
287
288// Register the algorithm into the AlgorithmFactory
289DECLARE_ALGORITHM(PolarizationCorrectionWildes)
290
291//----------------------------------------------------------------------------------------------
292
293
294const std::string PolarizationCorrectionWildes::name() const { return "PolarizationCorrectionWildes"; }
295
297int PolarizationCorrectionWildes::version() const { return 1; }
298
300const std::string PolarizationCorrectionWildes::category() const { return "Reflectometry"; }
301
303const std::string PolarizationCorrectionWildes::summary() const {
304 return "Corrects a group of polarization analysis workspaces for polarizer "
305 "and analyzer efficiencies.";
306}
307
309const std::vector<std::string> PolarizationCorrectionWildes::seeAlso() const {
310 return {"PolarizationEfficiencyCor", "PolarizationEfficienciesWildes"};
311}
312
317size_t PolarizationCorrectionWildes::WorkspaceMap::size() const noexcept {
318 return (mmWS ? 1 : 0) + (mpWS ? 1 : 0) + (pmWS ? 1 : 0) + (ppWS ? 1 : 0);
319}
320
321//----------------------------------------------------------------------------------------------
324void PolarizationCorrectionWildes::init() {
325 declareProperty(std::make_unique<Kernel::ArrayProperty<std::string>>(
326 Prop::INPUT_WS, "", std::make_shared<API::ADSValidator>(), Kernel::Direction::Input),
327 "A list of workspaces to be corrected corresponding to the "
328 "flipper configurations.");
329 declareProperty(
330 std::make_unique<API::WorkspaceProperty<API::WorkspaceGroup>>(Prop::OUTPUT_WS, "", Kernel::Direction::Output),
331 "A group of polarization efficiency corrected workspaces.");
332
333 const auto flipperConfigValidator =
334 std::make_shared<Kernel::SpinStateValidator>(std::unordered_set<int>{1, 2, 3, 4}, true);
335 declareProperty(Prop::FLIPPERS,
336 std::string(FlipperConfigurations::OFF_OFF) + ", " + FlipperConfigurations::OFF_ON + ", " +
337 FlipperConfigurations::ON_OFF + ", " + FlipperConfigurations::ON_ON,
338 flipperConfigValidator, "Flipper configurations of the input workspaces.");
339 const auto spinStateValidator =
340 std::make_shared<Kernel::SpinStateValidator>(std::unordered_set<int>{0, 2, 4}, false, "+", "-", true);
341 declareProperty(Prop::SPIN_STATES, "", spinStateValidator, "The order of the spin states in the output workspace.");
342 declareProperty(
343 std::make_unique<API::WorkspaceProperty<API::MatrixWorkspace>>(Prop::EFFICIENCIES, "", Kernel::Direction::Input),
344 "A workspace containing the efficiency factors P1, P2, F1 and F2 as "
345 "histograms");
346 declareProperty(
347 Prop::ADD_SPIN_STATE_LOG, false,
348 "Whether to add the final spin state into the sample log of each child workspace in the output group.");
349}
350
351//----------------------------------------------------------------------------------------------
354void PolarizationCorrectionWildes::exec() {
355 const std::string flipperProperty = getProperty(Prop::FLIPPERS);
356 const auto flippers = Kernel::SpinStateHelpers::splitSpinStateString(flipperProperty);
357 // Check if the input flipper configuration includes an analyser
358 const bool hasAnalyser = flippers.front().size() > 1;
359 const auto inputs = mapInputsToDirections(flippers);
360 checkConsistentNumberHistograms(inputs);
361 const EfficiencyMap efficiencies = efficiencyFactors();
362 checkConsistentX(inputs, efficiencies);
363 WorkspaceMap outputs;
364 switch (inputs.size()) {
365 case 1:
366 outputs = directBeamCorrections(inputs, efficiencies);
367 break;
368 case 2:
369 if (hasAnalyser) {
370 outputs = twoInputCorrections(inputs, efficiencies);
371 } else {
372 outputs = analyzerlessCorrections(inputs, efficiencies);
373 }
374 break;
375 case 3:
376 outputs = threeInputCorrections(inputs, efficiencies);
377 break;
378 case 4:
379 outputs = fullCorrections(inputs, efficiencies);
380 }
381 setProperty(Prop::OUTPUT_WS, groupOutput(outputs, hasAnalyser));
382}
383
388std::map<std::string, std::string> PolarizationCorrectionWildes::validateInputs() {
389 std::map<std::string, std::string> issues;
390 API::MatrixWorkspace_const_sptr factorWS = getProperty(Prop::EFFICIENCIES);
391 if (factorWS) {
392 const auto &factorAxis = factorWS->getAxis(1);
393 if (!factorAxis) {
394 issues[Prop::EFFICIENCIES] = "The workspace is missing a vertical axis.";
395 } else if (!factorAxis->isText()) {
396 issues[Prop::EFFICIENCIES] = "The vertical axis in the workspace is not text axis.";
397 } else if (factorWS->getNumberHistograms() < 4) {
398 issues[Prop::EFFICIENCIES] = "The workspace should contain at least 4 histograms.";
399 } else {
400 std::vector<std::string> tags{{"P1", "P2", "F1", "F2"}};
401 for (size_t i = 0; i != factorAxis->length(); ++i) {
402 const auto label = factorAxis->label(i);
403 auto found = std::find(tags.begin(), tags.end(), label);
404 if (found != tags.cend()) {
405 std::swap(tags.back(), *found);
406 tags.pop_back();
407 }
408 }
409 if (!tags.empty()) {
410 issues[Prop::EFFICIENCIES] = "A histogram labeled " + tags.front() + " is missing from the workspace.";
411 }
412 }
413 }
414 const std::vector<std::string> inputs = getProperty(Prop::INPUT_WS);
415 const auto flipperConfig = Kernel::SpinStateHelpers::splitSpinStateString(getPropertyValue(Prop::FLIPPERS));
416 const auto flipperCount = flipperConfig.size();
417 if (inputs.size() != flipperCount) {
418 issues[Prop::FLIPPERS] = "The number of flipper configurations (" + std::to_string(flipperCount) +
419 ") does not match the number of input workspaces (" + std::to_string(inputs.size()) + ")";
420 }
421 // SpinStates checks.
422 const auto spinStates = Kernel::SpinStateHelpers::splitSpinStateString(getPropertyValue(Prop::SPIN_STATES));
423 if (inputs.size() == 1 && !spinStates.empty()) {
424 issues[Prop::SPIN_STATES] = "Output workspace order cannot be set for direct beam calculations.";
425 } else if (!spinStates.empty()) {
426 if (flipperConfig.front().size() == 1 && spinStates.size() != 2) {
427 issues[Prop::SPIN_STATES] =
428 "Incorrect number of workspaces in output configuration: " + std::to_string(spinStates.size()) +
429 ". Only two output workspaces are produced when an analyzer is not used.";
430 }
431 if (flipperConfig.front().size() == 2 && spinStates.size() != 4) {
432 issues[Prop::SPIN_STATES] =
433 "Incorrect number of workspaces in output configuration: " + std::to_string(spinStates.size()) +
434 ". Four output workspaces are produced by the corrections.";
435 }
436 }
437 return issues;
438}
439
444void PolarizationCorrectionWildes::checkConsistentNumberHistograms(const WorkspaceMap &inputs) {
445 size_t nHist{0};
446 bool nHistValid{false};
447 // A local helper function to check the number of histograms.
448 auto checkNHist = [&nHist, &nHistValid](const API::MatrixWorkspace_sptr &ws, const std::string &tag) {
449 if (nHistValid) {
450 if (nHist != ws->getNumberHistograms()) {
451 throw std::runtime_error("Number of histograms mismatch in " + tag);
452 }
453 } else {
454 nHist = ws->getNumberHistograms();
455 nHistValid = true;
456 }
457 };
458 if (inputs.mmWS) {
459 checkNHist(inputs.mmWS, FlipperConfigurations::ON_ON);
460 }
461 if (inputs.mpWS) {
462 checkNHist(inputs.mpWS, FlipperConfigurations::ON_OFF);
463 }
464 if (inputs.pmWS) {
465 checkNHist(inputs.pmWS, FlipperConfigurations::OFF_ON);
466 }
467 if (inputs.ppWS) {
468 checkNHist(inputs.ppWS, FlipperConfigurations::OFF_OFF);
469 }
470}
471
477void PolarizationCorrectionWildes::checkConsistentX(const WorkspaceMap &inputs, const EfficiencyMap &efficiencies) {
478 // Compare everything to F1 efficiency.
479 const auto &F1x = efficiencies.F1->x();
480 // A local helper function to check a HistogramX against F1.
481 auto checkX = [&F1x](const HistogramData::HistogramX &x, const std::string &tag) {
482 if (x.size() != F1x.size()) {
483 throw std::runtime_error("Mismatch of histogram lengths between F1 and " + tag + '.');
484 }
485 for (size_t i = 0; i != x.size(); ++i) {
486 if (x[i] != F1x[i]) {
487 throw std::runtime_error("Mismatch of X data between F1 and " + tag + '.');
488 }
489 }
490 };
491 const auto &F2x = efficiencies.F2->x();
492 checkX(F2x, "F2");
493 const auto &P1x = efficiencies.P1->x();
494 checkX(P1x, "P1");
495 const auto &P2x = efficiencies.P2->x();
496 checkX(P2x, "P2");
497 // A local helper function to check an input workspace against F1.
498 auto checkWS = [&checkX](const API::MatrixWorkspace_sptr &ws, const std::string &tag) {
499 const auto nHist = ws->getNumberHistograms();
500 for (size_t i = 0; i != nHist; ++i) {
501 checkX(ws->x(i), tag);
502 }
503 };
504 if (inputs.mmWS) {
505 checkWS(inputs.mmWS, FlipperConfigurations::ON_ON);
506 }
507 if (inputs.mpWS) {
508 checkWS(inputs.mpWS, FlipperConfigurations::ON_OFF);
509 }
510 if (inputs.pmWS) {
511 checkWS(inputs.pmWS, FlipperConfigurations::OFF_ON);
512 }
513 if (inputs.ppWS) {
514 checkWS(inputs.ppWS, FlipperConfigurations::OFF_OFF);
515 }
516}
517
526API::WorkspaceGroup_sptr PolarizationCorrectionWildes::groupOutput(const WorkspaceMap &outputs,
527 const bool hasAnalyser) {
528 const auto &outWSName = getPropertyValue(Prop::OUTPUT_WS);
529 auto spinStateOrder = getPropertyValue(Prop::SPIN_STATES);
530 const bool addSpinStateLog = getProperty(Prop::ADD_SPIN_STATE_LOG);
531
532 std::vector<std::string> names;
533 if (!spinStateOrder.empty()) {
534 names.resize(Kernel::SpinStateHelpers::splitSpinStateString(spinStateOrder).size());
535 }
536
537 if (outputs.ppWS) {
538 addSpinStateOutput(names, spinStateOrder, outWSName, outputs.ppWS, SpinStateConfigurationsWildes::PLUS_PLUS,
539 addSpinStateLog, hasAnalyser);
540 }
541 if (outputs.pmWS) {
542 addSpinStateOutput(names, spinStateOrder, outWSName, outputs.pmWS, SpinStateConfigurationsWildes::PLUS_MINUS,
543 addSpinStateLog, hasAnalyser);
544 }
545 if (outputs.mpWS) {
546 addSpinStateOutput(names, spinStateOrder, outWSName, outputs.mpWS, SpinStateConfigurationsWildes::MINUS_PLUS,
547 addSpinStateLog, hasAnalyser);
548 }
549 if (outputs.mmWS) {
550 addSpinStateOutput(names, spinStateOrder, outWSName, outputs.mmWS, SpinStateConfigurationsWildes::MINUS_MINUS,
551 addSpinStateLog, hasAnalyser);
552 }
553
554 auto group = createChildAlgorithm("GroupWorkspaces");
555 group->initialize();
556 group->setProperty("InputWorkspaces", names);
557 group->setProperty("OutputWorkspace", outWSName);
558 group->execute();
559 API::WorkspaceGroup_sptr outWS = group->getProperty("OutputWorkspace");
560 return outWS;
561}
562
575void PolarizationCorrectionWildes::addSpinStateOutput(std::vector<std::string> &names,
576 const std::string &spinStateOrder, const std::string &baseName,
577 const API::MatrixWorkspace_sptr &ws, const std::string &spinState,
578 const bool addSpinStateLog, const bool hasAnalyser) {
579 if (addSpinStateLog) {
580 addSpinStateLogToWs(ws, spinState, hasAnalyser);
581 }
582
583 if (spinStateOrder.empty()) {
584 names.emplace_back(baseName + "_" + spinState);
585 API::AnalysisDataService::Instance().addOrReplace(names.back(), ws);
586 } else {
587 const auto &maybeIndex = Kernel::SpinStateHelpers::indexOfWorkspaceForSpinState(
588 Kernel::SpinStateHelpers::splitSpinStateString(spinStateOrder), spinState);
589 if (!maybeIndex.has_value()) {
590 throw std::invalid_argument("Required spin state (" + spinState + ") not found in spin state order (" +
591 spinStateOrder + ").");
592 }
593 const auto index = maybeIndex.value();
594 names[index] = baseName + "_" + spinState;
595 API::AnalysisDataService::Instance().addOrReplace(names[index], ws);
596 }
597}
598
605void PolarizationCorrectionWildes::addSpinStateLogToWs(const API::MatrixWorkspace_sptr &ws,
606 const std::string &spinState, const bool hasAnalyser) {
607 if (!hasAnalyser) {
608 if (spinState == SpinStateConfigurationsWildes::PLUS_PLUS) {
609 SpinStatesORSO::addORSOLogForSpinState(ws, SpinStateConfigurationsWildes::PLUS);
610 return;
611 }
612
613 if (spinState == SpinStateConfigurationsWildes::MINUS_MINUS) {
614 SpinStatesORSO::addORSOLogForSpinState(ws, SpinStateConfigurationsWildes::MINUS);
615 return;
616 }
617 }
618
619 SpinStatesORSO::addORSOLogForSpinState(ws, spinState);
620}
621
626PolarizationCorrectionWildes::EfficiencyMap PolarizationCorrectionWildes::efficiencyFactors() {
627 EfficiencyMap e;
628 API::MatrixWorkspace_const_sptr factorWS = getProperty(Prop::EFFICIENCIES);
629 const auto &vertAxis = factorWS->getAxis(1);
630 for (size_t i = 0; i != vertAxis->length(); ++i) {
631 const auto label = vertAxis->label(i);
632 if (label == "P1") {
633 e.P1 = &factorWS->getSpectrum(i);
634 } else if (label == "P2") {
635 e.P2 = &factorWS->getSpectrum(i);
636 } else if (label == "F1") {
637 e.F1 = &factorWS->getSpectrum(i);
638 } else if (label == "F2") {
639 e.F2 = &factorWS->getSpectrum(i);
640 }
641 // Ignore other histograms such as 'Phi' in ILL's efficiency ws.
642 }
643 return e;
644}
645
653PolarizationCorrectionWildes::WorkspaceMap
654PolarizationCorrectionWildes::directBeamCorrections(const WorkspaceMap &inputs, const EfficiencyMap &efficiencies) {
655 using namespace boost::math;
656 checkInputExists(inputs.ppWS, FlipperConfigurations::OFF);
657 WorkspaceMap outputs;
658 outputs.ppWS = createWorkspaceWithHistory(inputs.ppWS);
659 const size_t nHisto = inputs.ppWS->getNumberHistograms();
660 // cppcheck-suppress unreadVariable
661 const bool threadSafe = Kernel::threadSafe(*inputs.ppWS, *outputs.ppWS);
662 PARALLEL_FOR_IF(threadSafe)
663 for (int wsIndex = 0; wsIndex < static_cast<int>(nHisto); ++wsIndex) {
664 PARALLEL_START_INTERRUPT_REGION
665 const auto &ppY = inputs.ppWS->y(wsIndex);
666 const auto &ppE = inputs.ppWS->e(wsIndex);
667 auto &ppYOut = outputs.ppWS->mutableY(wsIndex);
668 auto &ppEOut = outputs.ppWS->mutableE(wsIndex);
669 PARALLEL_FOR_IF(threadSafe)
670 for (int binIndex = 0; binIndex < static_cast<int>(ppY.size()); ++binIndex) {
671 const auto P1 = efficiencies.P1->y()[binIndex];
672 const auto P2 = efficiencies.P2->y()[binIndex];
673 const double f = 1. - P1 - P2 + 2. * P1 * P2;
674 ppYOut[binIndex] = ppY[binIndex] / f;
675 const auto P1E = efficiencies.P1->e()[binIndex];
676 const auto P2E = efficiencies.P2->e()[binIndex];
677 const auto e1 = pow<2>(P1E * (2. * P1 - 1.) / pow<2>(f) * ppY[binIndex]);
678 const auto e2 = pow<2>(P2E * (2. * P2 - 1.) / pow<2>(f) * ppY[binIndex]);
679 const auto e3 = pow<2>(ppE[binIndex] / f);
680 const auto errorSum = std::sqrt(e1 + e2 + e3);
681 ppEOut[binIndex] = errorSum;
682 }
683 PARALLEL_END_INTERRUPT_REGION
684 }
685 PARALLEL_CHECK_INTERRUPT_REGION
686 return outputs;
687}
688
697PolarizationCorrectionWildes::WorkspaceMap
698PolarizationCorrectionWildes::analyzerlessCorrections(const WorkspaceMap &inputs, const EfficiencyMap &efficiencies) {
699 using namespace boost::math;
700 checkInputExists(inputs.mmWS, FlipperConfigurations::ON);
701 checkInputExists(inputs.ppWS, FlipperConfigurations::OFF);
702 WorkspaceMap outputs;
703 outputs.mmWS = createWorkspaceWithHistory(inputs.mmWS);
704 outputs.ppWS = createWorkspaceWithHistory(inputs.ppWS);
705 const size_t nHisto = inputs.mmWS->getNumberHistograms();
706 // cppcheck-suppress unreadVariable
707 const bool threadSafe = Kernel::threadSafe(*inputs.mmWS, *inputs.ppWS, *outputs.mmWS, *outputs.ppWS);
708 PARALLEL_FOR_IF(threadSafe)
709 for (int wsIndex = 0; wsIndex < static_cast<int>(nHisto); ++wsIndex) {
710 PARALLEL_START_INTERRUPT_REGION
711 const auto &mmY = inputs.mmWS->y(wsIndex);
712 const auto &mmE = inputs.mmWS->e(wsIndex);
713 const auto &ppY = inputs.ppWS->y(wsIndex);
714 const auto &ppE = inputs.ppWS->e(wsIndex);
715 auto &mmYOut = outputs.mmWS->mutableY(wsIndex);
716 auto &mmEOut = outputs.mmWS->mutableE(wsIndex);
717 auto &ppYOut = outputs.ppWS->mutableY(wsIndex);
718 auto &ppEOut = outputs.ppWS->mutableE(wsIndex);
719 PARALLEL_FOR_IF(threadSafe)
720 for (int binIndex = 0; binIndex < static_cast<int>(mmY.size()); ++binIndex) {
721 const auto F1 = efficiencies.F1->y()[binIndex];
722 const auto P1 = efficiencies.P1->y()[binIndex];
723 Eigen::Matrix2d F1m;
724 F1m << 1., 0., (F1 - 1.) / F1, 1. / F1;
725 const double divisor = (2. * P1 - 1.);
726 const double off = (P1 - 1.) / divisor;
727 const double diag = P1 / divisor;
728 Eigen::Matrix2d P1m;
729 P1m << diag, off, off, diag;
730 const Eigen::Vector2d intensities(ppY[binIndex], mmY[binIndex]);
731 const auto PFProduct = P1m * F1m;
732 const auto corrected = PFProduct * intensities;
733 ppYOut[binIndex] = corrected[0];
734 mmYOut[binIndex] = corrected[1];
735 const auto F1E = efficiencies.F1->e()[binIndex];
736 const auto P1E = efficiencies.P1->e()[binIndex];
737 const auto elemE1 = -1. / pow<2>(F1) * F1E;
738 Eigen::Matrix2d F1Em;
739 F1Em << 0., 0., -elemE1, elemE1;
740 const auto elemE2 = 1. / pow<2>(divisor) * P1E;
741 Eigen::Matrix2d P1Em;
742 P1Em << elemE2, -elemE2, -elemE2, elemE2;
743 const Eigen::Vector2d errors(ppE[binIndex], mmE[binIndex]);
744 const auto e1 = (P1Em * F1m * intensities).array();
745 const auto e2 = (P1m * F1Em * intensities).array();
746 const auto sqPFProduct = (PFProduct.array() * PFProduct.array()).matrix();
747 const auto sqErrors = (errors.array() * errors.array()).matrix();
748 const auto e3 = (sqPFProduct * sqErrors).array();
749 const auto errorSum = (e1 * e1 + e2 * e2 + e3).sqrt();
750 ppEOut[binIndex] = errorSum[0];
751 mmEOut[binIndex] = errorSum[1];
752 }
753 PARALLEL_END_INTERRUPT_REGION
754 }
755 PARALLEL_CHECK_INTERRUPT_REGION
756 return outputs;
757}
758
768PolarizationCorrectionWildes::WorkspaceMap
769PolarizationCorrectionWildes::twoInputCorrections(const WorkspaceMap &inputs, const EfficiencyMap &efficiencies) {
770 using namespace boost::math;
771 checkInputExists(inputs.mmWS, FlipperConfigurations::ON_ON);
772 checkInputExists(inputs.ppWS, FlipperConfigurations::OFF_OFF);
773 WorkspaceMap fullInputs = inputs;
774 fullInputs.mpWS = createWorkspaceWithHistory(inputs.mmWS);
775 fullInputs.pmWS = createWorkspaceWithHistory(inputs.ppWS);
776 twoInputsSolve01And10(fullInputs, inputs, efficiencies);
777 return fullCorrections(fullInputs, efficiencies);
778}
779
789PolarizationCorrectionWildes::WorkspaceMap
790PolarizationCorrectionWildes::threeInputCorrections(const WorkspaceMap &inputs, const EfficiencyMap &efficiencies) {
791 WorkspaceMap fullInputs = inputs;
792 checkInputExists(inputs.mmWS, FlipperConfigurations::ON_ON);
793 checkInputExists(inputs.ppWS, FlipperConfigurations::OFF_OFF);
794 if (!inputs.mpWS) {
795 checkInputExists(inputs.pmWS, FlipperConfigurations::OFF_ON);
796 threeInputsSolve10(fullInputs, efficiencies);
797 } else {
798 checkInputExists(inputs.mpWS, FlipperConfigurations::ON_OFF);
799 threeInputsSolve01(fullInputs, efficiencies);
800 }
801 return fullCorrections(fullInputs, efficiencies);
802}
803
812PolarizationCorrectionWildes::WorkspaceMap
813PolarizationCorrectionWildes::fullCorrections(const WorkspaceMap &inputs, const EfficiencyMap &efficiencies) {
814 using namespace boost::math;
815 checkInputExists(inputs.mmWS, FlipperConfigurations::ON_ON);
816 checkInputExists(inputs.mpWS, FlipperConfigurations::ON_OFF);
817 checkInputExists(inputs.pmWS, FlipperConfigurations::OFF_ON);
818 checkInputExists(inputs.ppWS, FlipperConfigurations::OFF_OFF);
819 WorkspaceMap outputs;
820 outputs.mmWS = createWorkspaceWithHistory(inputs.mmWS);
821 outputs.mpWS = createWorkspaceWithHistory(inputs.mpWS);
822 outputs.pmWS = createWorkspaceWithHistory(inputs.pmWS);
823 outputs.ppWS = createWorkspaceWithHistory(inputs.ppWS);
824 const auto F1 = efficiencies.F1->y();
825 const auto F1E = efficiencies.F1->e();
826 const auto F2 = efficiencies.F2->y();
827 const auto F2E = efficiencies.F2->e();
828 const auto P1 = efficiencies.P1->y();
829 const auto P1E = efficiencies.P1->e();
830 const auto P2 = efficiencies.P2->y();
831 const auto P2E = efficiencies.P2->e();
832 const size_t nHisto = inputs.mmWS->getNumberHistograms();
833 // cppcheck-suppress unreadVariable
834 const bool threadSafe = Kernel::threadSafe(*inputs.mmWS, *inputs.mpWS, *inputs.pmWS, *inputs.ppWS, *outputs.mmWS,
835 *outputs.mpWS, *outputs.pmWS, *outputs.ppWS);
836 PARALLEL_FOR_IF(threadSafe)
837 for (int wsIndex = 0; wsIndex < static_cast<int>(nHisto); ++wsIndex) {
838 PARALLEL_START_INTERRUPT_REGION
839 const auto &mmY = inputs.mmWS->y(wsIndex);
840 const auto &mmE = inputs.mmWS->e(wsIndex);
841 const auto &mpY = inputs.mpWS->y(wsIndex);
842 const auto &mpE = inputs.mpWS->e(wsIndex);
843 const auto &pmY = inputs.pmWS->y(wsIndex);
844 const auto &pmE = inputs.pmWS->e(wsIndex);
845 const auto &ppY = inputs.ppWS->y(wsIndex);
846 const auto &ppE = inputs.ppWS->e(wsIndex);
847 auto &mmYOut = outputs.mmWS->mutableY(wsIndex);
848 auto &mmEOut = outputs.mmWS->mutableE(wsIndex);
849 auto &mpYOut = outputs.mpWS->mutableY(wsIndex);
850 auto &mpEOut = outputs.mpWS->mutableE(wsIndex);
851 auto &pmYOut = outputs.pmWS->mutableY(wsIndex);
852 auto &pmEOut = outputs.pmWS->mutableE(wsIndex);
853 auto &ppYOut = outputs.ppWS->mutableY(wsIndex);
854 auto &ppEOut = outputs.ppWS->mutableE(wsIndex);
855 PARALLEL_FOR_IF(threadSafe)
856 for (int binIndex = 0; binIndex < static_cast<int>(mmY.size()); ++binIndex) {
857 Eigen::Vector4d corrected;
858 Eigen::Vector4d errors;
859 fourInputsCorrectedAndErrors(corrected, errors, ppY[binIndex], ppE[binIndex], pmY[binIndex], pmE[binIndex],
860 mpY[binIndex], mpE[binIndex], mmY[binIndex], mmE[binIndex], F1[binIndex],
861 F1E[binIndex], F2[binIndex], F2E[binIndex], P1[binIndex], P1E[binIndex],
862 P2[binIndex], P2E[binIndex]);
863 ppYOut[binIndex] = corrected[0];
864 pmYOut[binIndex] = corrected[1];
865 mpYOut[binIndex] = corrected[2];
866 mmYOut[binIndex] = corrected[3];
867 ppEOut[binIndex] = errors[0];
868 pmEOut[binIndex] = errors[1];
869 mpEOut[binIndex] = errors[2];
870 mmEOut[binIndex] = errors[3];
871 }
872 PARALLEL_END_INTERRUPT_REGION
873 }
874 PARALLEL_CHECK_INTERRUPT_REGION
875 return outputs;
876}
877
883PolarizationCorrectionWildes::WorkspaceMap
884PolarizationCorrectionWildes::mapInputsToDirections(const std::vector<std::string> &flippers) {
885 const std::vector<std::string> inputNames = getProperty(Prop::INPUT_WS);
886 WorkspaceMap inputs;
887 for (size_t i = 0; i < flippers.size(); ++i) {
888 auto ws = (API::AnalysisDataService::Instance().retrieveWS<API::MatrixWorkspace>(inputNames[i]));
889 if (!ws) {
890 throw std::runtime_error("One of the input workspaces doesn't seem to be a MatrixWorkspace.");
891 }
892 const auto &f = flippers[i];
893 if (f == FlipperConfigurations::ON_ON || f == FlipperConfigurations::ON) {
894 inputs.mmWS = ws;
895 } else if (f == FlipperConfigurations::ON_OFF) {
896 inputs.mpWS = ws;
897 } else if (f == FlipperConfigurations::OFF_ON) {
898 inputs.pmWS = ws;
899 } else if (f == FlipperConfigurations::OFF_OFF || f == FlipperConfigurations::OFF) {
900 inputs.ppWS = ws;
901 } else {
902 throw std::runtime_error(std::string{"Unknown entry in "} + Prop::FLIPPERS);
903 }
904 }
905 return inputs;
906}
907
914void PolarizationCorrectionWildes::threeInputsSolve01(WorkspaceMap &inputs, const EfficiencyMap &efficiencies) {
915 inputs.pmWS = createWorkspaceWithHistory(inputs.mpWS);
916 const auto &F1 = efficiencies.F1->y();
917 const auto &F2 = efficiencies.F2->y();
918 const auto &P1 = efficiencies.P1->y();
919 const auto &P2 = efficiencies.P2->y();
920 const auto nHisto = inputs.pmWS->getNumberHistograms();
921 // cppcheck-suppress unreadVariable
922 const bool threadSafe = Kernel::threadSafe(*inputs.mmWS, *inputs.mpWS, *inputs.pmWS, *inputs.ppWS);
923 PARALLEL_FOR_IF(threadSafe)
924 for (int wsIndex = 0; wsIndex < static_cast<int>(nHisto); ++wsIndex) {
925 PARALLEL_START_INTERRUPT_REGION
926 const auto &I00 = inputs.ppWS->y(wsIndex);
927 auto &I01 = inputs.pmWS->mutableY(wsIndex);
928 const auto &I10 = inputs.mpWS->y(wsIndex);
929 const auto &I11 = inputs.mmWS->y(wsIndex);
930 PARALLEL_FOR_IF(threadSafe)
931 for (int binIndex = 0; binIndex < static_cast<int>(I00.size()); ++binIndex) {
932 const auto f1 = F1[binIndex];
933 const auto f2 = F2[binIndex];
934 const auto p1 = P1[binIndex];
935 const auto p2 = P2[binIndex];
936 const auto i00 = I00[binIndex];
937 const auto i10 = I10[binIndex];
938 const auto i11 = I11[binIndex];
939 I01[binIndex] =
940 (f1 * i00 * (-1. + 2. * p1) - (i00 - i10 + i11) * (p1 - p2) - f2 * (i00 - i10) * (-1. + 2. * p2)) /
941 (-p1 + f1 * (-1. + 2. * p1) + p2);
942 // The errors are left to zero.
943 }
944 PARALLEL_END_INTERRUPT_REGION
945 }
946 PARALLEL_CHECK_INTERRUPT_REGION
947}
948
955void PolarizationCorrectionWildes::threeInputsSolve10(WorkspaceMap &inputs, const EfficiencyMap &efficiencies) {
956 inputs.mpWS = createWorkspaceWithHistory(inputs.pmWS);
957 const auto &F1 = efficiencies.F1->y();
958 const auto &F2 = efficiencies.F2->y();
959 const auto &P1 = efficiencies.P1->y();
960 const auto &P2 = efficiencies.P2->y();
961 const auto nHisto = inputs.mpWS->getNumberHistograms();
962 // cppcheck-suppress unreadVariable
963 const bool threadSafe = Kernel::threadSafe(*inputs.mmWS, *inputs.mpWS, *inputs.pmWS, *inputs.ppWS);
964 PARALLEL_FOR_IF(threadSafe)
965 for (int wsIndex = 0; wsIndex < static_cast<int>(nHisto); ++wsIndex) {
966 PARALLEL_START_INTERRUPT_REGION
967 const auto &I00 = inputs.ppWS->y(wsIndex);
968 const auto &I01 = inputs.pmWS->y(wsIndex);
969 auto &I10 = inputs.mpWS->mutableY(wsIndex);
970 const auto &I11 = inputs.mmWS->y(wsIndex);
971 PARALLEL_FOR_IF(threadSafe)
972 for (int binIndex = 0; binIndex < static_cast<int>(I00.size()); ++binIndex) {
973 const auto f1 = F1[binIndex];
974 const auto f2 = F2[binIndex];
975 const auto p1 = P1[binIndex];
976 const auto p2 = P2[binIndex];
977 const auto i00 = I00[binIndex];
978 const auto i01 = I01[binIndex];
979 const auto i11 = I11[binIndex];
980 I10[binIndex] =
981 (-f1 * (i00 - i01) * (-1. + 2. * p1) + (i00 - i01 + i11) * (p1 - p2) + f2 * i00 * (-1. + 2. * p2)) /
982 (p1 - p2 + f2 * (-1. + 2. * p2));
983 // The errors are left to zero.
984 }
985 PARALLEL_END_INTERRUPT_REGION
986 }
987 PARALLEL_CHECK_INTERRUPT_REGION
988}
989
997void PolarizationCorrectionWildes::twoInputsSolve01And10(WorkspaceMap &fullInputs, const WorkspaceMap &inputs,
998 const EfficiencyMap &efficiencies) {
999 using namespace boost::math;
1000 const auto &F1 = efficiencies.F1->y();
1001 const auto &F1E = efficiencies.F1->e();
1002 const auto &F2 = efficiencies.F2->y();
1003 const auto &F2E = efficiencies.F2->e();
1004 const auto &P1 = efficiencies.P1->y();
1005 const auto &P1E = efficiencies.P1->e();
1006 const auto &P2 = efficiencies.P2->y();
1007 const auto &P2E = efficiencies.P2->e();
1008 const auto nHisto = inputs.mmWS->getNumberHistograms();
1009 // cppcheck-suppress unreadVariable
1010 const bool threadSafe = Kernel::threadSafe(*inputs.mmWS, *inputs.ppWS, *fullInputs.pmWS, *fullInputs.mpWS);
1011 PARALLEL_FOR_IF(threadSafe)
1012 for (int wsIndex = 0; wsIndex < static_cast<int>(nHisto); ++wsIndex) {
1013 PARALLEL_START_INTERRUPT_REGION
1014 const auto &I00 = inputs.ppWS->y(wsIndex);
1015 const auto &E00 = inputs.ppWS->e(wsIndex);
1016 const auto &I11 = inputs.mmWS->y(wsIndex);
1017 const auto &E11 = inputs.mmWS->e(wsIndex);
1018 auto &I01 = fullInputs.pmWS->mutableY(wsIndex);
1019 auto &E01 = fullInputs.pmWS->mutableE(wsIndex);
1020 auto &I10 = fullInputs.mpWS->mutableY(wsIndex);
1021 auto &E10 = fullInputs.mpWS->mutableE(wsIndex);
1022 PARALLEL_FOR_IF(threadSafe)
1023 for (int binIndex = 0; binIndex < static_cast<int>(I00.size()); ++binIndex) {
1024 const auto i00 = I00[binIndex];
1025 const auto i11 = I11[binIndex];
1026 const auto f1 = F1[binIndex];
1027 const auto f2 = F2[binIndex];
1028 const auto p1 = P1[binIndex];
1029 const auto p2 = P2[binIndex];
1030 const auto a = -1. + p1 + 2. * p2 - 2. * p1 * p2;
1031 const auto b = -1. + 2. * p1;
1032 const auto c = -1. + 2. * p2;
1033 const auto d = -1. + p2;
1034 // Case: 01
1035 const auto divisor = f2 * p1 * a + f1 * b * (-d * p2 + f2 * (p1 + d) * c);
1036 I01[binIndex] = (f2 * i11 * p1 * a - f1 * i00 * b * (-f2 * pow<2>(c) + pow<2>(f2 * c) + d * p2)) / divisor;
1037 E01[binIndex] = twoInputsErrorEstimate01(i00, E00[binIndex], i11, E11[binIndex], p1, P1E[binIndex], p2,
1038 P2E[binIndex], f1, F1E[binIndex], f2, F2E[binIndex]);
1039 // Case: 10
1040 I10[binIndex] =
1041 (-pow<2>(f1) * f2 * i00 * pow<2>(b) * c + f2 * i00 * p1 * a + f1 * b * (-i11 * d * p2 + f2 * i00 * b * c)) /
1042 divisor;
1043 E10[binIndex] = twoInputsErrorEstimate10(i00, E00[binIndex], i11, E11[binIndex], p1, P1E[binIndex], p2,
1044 P2E[binIndex], f1, F1E[binIndex], f2, F2E[binIndex]);
1045 }
1046 PARALLEL_END_INTERRUPT_REGION
1047 }
1048 PARALLEL_CHECK_INTERRUPT_REGION
1049}
1050} // namespace Mantid::Algorithms
name
std::string name
Definition Run.cpp:60
DECLARE_ALGORITHM
#define DECLARE_ALGORITHM(classname)
Definition Algorithm.h:538
index
std::map< DeltaEMode::Type, std::string > index
Definition DeltaEMode.cpp:19
PARALLEL_START_INTERRUPT_REGION
#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
PARALLEL_END_INTERRUPT_REGION
#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...
Definition MultiThreaded.h:102
PARALLEL_FOR_IF
#define PARALLEL_FOR_IF(condition)
Empty definitions - to enable set your complier to enable openMP.
Definition MultiThreaded.h:218
PARALLEL_CHECK_INTERRUPT_REGION
#define PARALLEL_CHECK_INTERRUPT_REGION
Adds a check after a Parallel region to see if it was interupted.
Definition MultiThreaded.h:117
NXentrytype::group
@ group
Mantid::API::Algorithm::declareProperty
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:1803
Mantid::API::Algorithm::getPropertyValue
std::string getPropertyValue(const std::string &name) const override
Get the value of a property as a string.
Definition Algorithm.cpp:1916
Mantid::API::Algorithm::getProperty
TypedValue getProperty(const std::string &name) const override
Get the value of a property.
Definition Algorithm.cpp:1966
Mantid::API::Algorithm::createChildAlgorithm
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:860
Mantid::API::ISpectrum::x
const HistogramData::HistogramX & x() const
Definition ISpectrum.h:171
Mantid::API::ISpectrum::e
virtual const HistogramData::HistogramE & e() const
Definition ISpectrum.h:173
Mantid::API::ISpectrum::y
virtual const HistogramData::HistogramY & y() const
Definition ISpectrum.h:172
Mantid::API::MatrixWorkspace
Base MatrixWorkspace Abstract Class.
Definition MatrixWorkspace.h:92
Mantid::API::WorkspaceProperty
A property class for workspaces.
Definition WorkspaceProperty.h:51
Mantid::Algorithms::PolarizationCorrectionWildes
PolarizationCorrectionWildes : This algorithm corrects for non-ideal component efficiencies in polari...
Definition PolarizationCorrectionWildes.h:26
Mantid::Algorithms::PolarizationCorrectionWildes::checkConsistentX
void checkConsistentX(const WorkspaceMap &inputs, const EfficiencyMap &efficiencies)
Check that all workspaces and efficicencies have the same X data.
Definition PolarizationCorrectionWildes.cpp:477
Mantid::Algorithms::PolarizationCorrectionWildes::validateInputs
std::map< std::string, std::string > validateInputs() override
Validate the algorithm's input properties.
Definition PolarizationCorrectionWildes.cpp:388
Mantid::Algorithms::PolarizationCorrectionWildes::efficiencyFactors
EfficiencyMap efficiencyFactors()
Make a convenience access object to the efficiency factors.
Definition PolarizationCorrectionWildes.cpp:626
Mantid::Algorithms::PolarizationCorrectionWildes::threeInputsSolve10
void threeInputsSolve10(WorkspaceMap &inputs, const EfficiencyMap &efficiencies)
Solve in-place the 10 flipper configuration from the assumption that for the corrected intensities R0...
Definition PolarizationCorrectionWildes.cpp:955
Mantid::Algorithms::PolarizationCorrectionWildes::addSpinStateOutput
void addSpinStateOutput(std::vector< std::string > &names, const std::string &spinStateOrder, const std::string &baseName, const API::MatrixWorkspace_sptr &ws, const std::string &spinState, const bool addSpinStateLog, const bool hasAnalyser)
Add an output name in the correct position in the vector and to the ADS.
Definition PolarizationCorrectionWildes.cpp:575
Mantid::Algorithms::PolarizationCorrectionWildes::directBeamCorrections
WorkspaceMap directBeamCorrections(const WorkspaceMap &inputs, const EfficiencyMap &efficiencies)
Correct a direct beam measurement for non-ideal instrument effects.
Definition PolarizationCorrectionWildes.cpp:654
Mantid::Algorithms::PolarizationCorrectionWildes::init
void init() override
Initialize the algorithm's properties.
Definition PolarizationCorrectionWildes.cpp:324
Mantid::Algorithms::PolarizationCorrectionWildes::version
int version() const override
Algorithm's version for identification.
Definition PolarizationCorrectionWildes.cpp:297
Mantid::Algorithms::PolarizationCorrectionWildes::addSpinStateLogToWs
void addSpinStateLogToWs(const API::MatrixWorkspace_sptr &ws, const std::string &spinState, const bool hasAnalyser)
Adds a sample log to the workspace that gives the spin state of the data using Reflectometry ORSO not...
Definition PolarizationCorrectionWildes.cpp:605
Mantid::Algorithms::PolarizationCorrectionWildes::threeInputsSolve01
void threeInputsSolve01(WorkspaceMap &inputs, const EfficiencyMap &efficiencies)
Solve in-place the 01 flipper configuration from the assumption that for the corrected intensities,...
Definition PolarizationCorrectionWildes.cpp:914
Mantid::Algorithms::PolarizationCorrectionWildes::fullCorrections
WorkspaceMap fullCorrections(const WorkspaceMap &inputs, const EfficiencyMap &efficiencies)
Correct for non-ideal instrument effects.
Definition PolarizationCorrectionWildes.cpp:813
Mantid::Algorithms::PolarizationCorrectionWildes::seeAlso
const std::vector< std::string > seeAlso() const override
Algorithm's related algorithms.
Definition PolarizationCorrectionWildes.cpp:309
Mantid::Algorithms::PolarizationCorrectionWildes::category
const std::string category() const override
Algorithm's category for identification.
Definition PolarizationCorrectionWildes.cpp:300
Mantid::Algorithms::PolarizationCorrectionWildes::summary
const std::string summary() const override
Algorithm's summary for use in the GUI and help.
Definition PolarizationCorrectionWildes.cpp:303
Mantid::Algorithms::PolarizationCorrectionWildes::checkConsistentNumberHistograms
void checkConsistentNumberHistograms(const WorkspaceMap &inputs)
Check that all workspaces in inputs have the same number of histograms.
Definition PolarizationCorrectionWildes.cpp:444
Mantid::Algorithms::PolarizationCorrectionWildes::threeInputCorrections
WorkspaceMap threeInputCorrections(const WorkspaceMap &inputs, const EfficiencyMap &efficiencies)
Correct for non-ideal instrument effects.
Definition PolarizationCorrectionWildes.cpp:790
Mantid::Algorithms::PolarizationCorrectionWildes::twoInputsSolve01And10
void twoInputsSolve01And10(WorkspaceMap &fullInputs, const WorkspaceMap &inputs, const EfficiencyMap &efficiencies)
Solve in-place the 01 and 10 flipper configurations from the assumption that for the corrected intens...
Definition PolarizationCorrectionWildes.cpp:997
Mantid::Algorithms::PolarizationCorrectionWildes::twoInputCorrections
WorkspaceMap twoInputCorrections(const WorkspaceMap &inputs, const EfficiencyMap &efficiencies)
Correct for non-ideal instrument effects.
Definition PolarizationCorrectionWildes.cpp:769
Mantid::Algorithms::PolarizationCorrectionWildes::mapInputsToDirections
WorkspaceMap mapInputsToDirections(const std::vector< std::string > &flippers)
Make a set of workspaces to correct from input properties.
Definition PolarizationCorrectionWildes.cpp:884
Mantid::Algorithms::PolarizationCorrectionWildes::groupOutput
API::WorkspaceGroup_sptr groupOutput(const WorkspaceMap &outputs, const bool hasAnalyser)
Make a workspace group out of the given set of workspaces.
Definition PolarizationCorrectionWildes.cpp:526
Mantid::Algorithms::PolarizationCorrectionWildes::analyzerlessCorrections
WorkspaceMap analyzerlessCorrections(const WorkspaceMap &inputs, const EfficiencyMap &efficiencies)
Correct for non-ideal instrument effects.
Definition PolarizationCorrectionWildes.cpp:698
Mantid::Kernel::ArrayProperty
Support for a property that holds an array of values.
Definition ArrayProperty.h:28
Mantid::Kernel::IPropertyManager::setProperty
IPropertyManager * setProperty(const std::string &name, const T &value)
Templated method to set the value of a PropertyWithValue.
Definition IPropertyManager.h:270
Mantid::API::WorkspaceGroup_sptr
std::shared_ptr< WorkspaceGroup > WorkspaceGroup_sptr
shared pointer to Mantid::API::WorkspaceGroup
Definition WorkspaceGroup_fwd.h:20
Mantid::API::MatrixWorkspace_const_sptr
std::shared_ptr< const MatrixWorkspace > MatrixWorkspace_const_sptr
shared pointer to the matrix workspace base class (const version)
Definition MatrixWorkspace.h:587
Mantid::API::MatrixWorkspace_sptr
std::shared_ptr< MatrixWorkspace > MatrixWorkspace_sptr
shared pointer to the matrix workspace base class
Definition MatrixWorkspace.h:585
Mantid::Algorithms::SpinStatesORSO::addORSOLogForSpinState
MANTID_ALGORITHMS_DLL void addORSOLogForSpinState(const Mantid::API::MatrixWorkspace_sptr &ws, const std::string &spinState)
Definition PolarizationCorrectionsHelpers.cpp:87
Mantid::CurveFitting::Functions::f1
double f1(const double x, const double G, const double w0)
Definition DynamicKuboToyabe.cpp:132
Mantid::Kernel::SpinStateHelpers::splitSpinStateString
MANTID_KERNEL_DLL std::vector< std::string > splitSpinStateString(const std::string &spinStates)
Definition SpinStateHelpers.cpp:40
Mantid::Kernel::SpinStateHelpers::indexOfWorkspaceForSpinState
MANTID_KERNEL_DLL std::optional< size_t > indexOfWorkspaceForSpinState(const std::vector< std::string > &spinStateOrder, std::string targetSpinState)
Definition SpinStateHelpers.cpp:24
Mantid::Kernel::threadSafe
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
Prop
String constants for algorithm's properties.
Definition ApplyFloodWorkspace.cpp:21
Prop::ADD_SPIN_STATE_LOG
static const std::string ADD_SPIN_STATE_LOG
Definition PolarizationCorrectionFredrikze.cpp:34
Prop::EFFICIENCIES
static const std::string EFFICIENCIES
Definition PolarizationCorrectionFredrikze.cpp:29
std
STL namespace.
std::to_string
std::string to_string(const wide_integer< Bits, Signed > &n)
Mantid::Algorithms::PolarizationCorrectionWildes::WorkspaceMap
A convenience set of workspaces corresponding flipper configurations.
Definition PolarizationCorrectionWildes.h:36
Mantid::Algorithms::PolarizationCorrectionWildes::WorkspaceMap::size
size_t size() const noexcept
Count the non-nullptr workspaces.
Definition PolarizationCorrectionWildes.cpp:317
Mantid::Kernel::Direction::Input
@ Input
An input workspace.
Definition Property.h:53
Mantid::Kernel::Direction::Output
@ Output
An output workspace.
Definition Property.h:54
Generated by doxygen 1.9.8