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 INPUT_WS_GROUP{"InputWorkspaceGroup"};
33static const std::string OUTPUT_WS{"OutputWorkspace"};
34static const std::string ADD_SPIN_STATE_LOG{"AddSpinStateToLog"};
35} // namespace Prop
36
42void checkInputExists(const Mantid::API::MatrixWorkspace_sptr &ws, const std::string &tag) {
43 if (!ws) {
44 throw std::runtime_error("A workspace designated as " + tag + " is missing in inputs.");
45 }
46}
47
69void fourInputsCorrectedAndErrors(Eigen::Vector4d &corrected, Eigen::Vector4d &errors, const double ppy,
70 const double ppyE, const double pmy, const double pmyE, const double mpy,
71 const double mpyE, const double mmy, const double mmyE, const double f1,
72 const double f1E, const double f2, const double f2E, const double p1,
73 const double p1E, const double p2, const double p2E) {
74 using namespace boost::math;
75 // Note that f1 and f2 correspond to 1-F1 and 1-F2 in [Wildes, 1999].
76 // These are inverted forms of the efficiency matrices.
77 const auto diag1 = 1. / f1;
78 const auto off1 = (f1 - 1.) / f1;
79 Eigen::Matrix4d F1m;
80 F1m << 1., 0., 0., 0., 0., 1., 0., 0., off1, 0., diag1, 0., 0., off1, 0., diag1;
81
82 const auto diag2 = 1. / f2;
83 const auto off2 = (f2 - 1.) / f2;
84 Eigen::Matrix4d F2m;
85 F2m << 1., 0., 0., 0., off2, diag2, 0., 0., 0., 0., 1., 0., 0., 0., off2, diag2;
86 const auto diag3 = p1 / (2. * p1 - 1);
87 const auto off3 = (p1 - 1.) / (2. * p1 - 1.);
88 Eigen::Matrix4d P1m;
89 P1m << diag3, 0, off3, 0, 0, diag3, 0, off3, off3, 0, diag3, 0, 0, off3, 0, diag3;
90 const auto diag4 = p2 / (2. * p2 - 1.);
91 const auto off4 = (p2 - 1.) / (2. * p2 - 1.);
92 Eigen::Matrix4d P2m;
93 P2m << diag4, off4, 0., 0., off4, diag4, 0., 0., 0., 0., diag4, off4, 0., 0., off4, diag4;
94 const Eigen::Vector4d intensities(ppy, pmy, mpy, mmy);
95 const auto FProduct = F2m * F1m;
96 const auto PProduct = P2m * P1m;
97 const auto PFProduct = PProduct * FProduct;
98 corrected = PFProduct * intensities;
99 // The error matrices here are element-wise algebraic derivatives of
100 // the matrices above, multiplied by the error.
101 const auto elemE1 = -1. / pow<2>(f1) * f1E;
102 Eigen::Matrix4d F1Em;
103 F1Em << 0., 0., 0., 0., 0., 0., 0., 0., -elemE1, 0., elemE1, 0., 0., -elemE1, 0., elemE1;
104 const auto elemE2 = -1. / pow<2>(f2) * f2E;
105 Eigen::Matrix4d F2Em;
106 F2Em << 0., 0., 0., 0., -elemE2, elemE2, 0., 0., 0., 0., 0., 0., 0., 0., -elemE2, elemE2;
107 const auto elemE3 = 1. / pow<2>(2. * p1 - 1.) * p1E;
108 Eigen::Matrix4d P1Em;
109 P1Em << elemE3, 0., -elemE3, 0., 0., elemE3, 0., -elemE3, -elemE3, 0., elemE3, 0., 0., -elemE3, 0., elemE3;
110 const auto elemE4 = 1. / pow<2>(2. * p2 - 1.) * p2E;
111 Eigen::Matrix4d P2Em;
112 P2Em << elemE4, -elemE4, 0., 0., -elemE4, elemE4, 0., 0., 0., 0., elemE4, -elemE4, 0., 0., -elemE4, elemE4;
113 const Eigen::Vector4d yErrors(ppyE, pmyE, mpyE, mmyE);
114 const auto e1 = (P2Em * P1m * FProduct * intensities).array();
115 const auto e2 = (P2m * P1Em * FProduct * intensities).array();
116 const auto e3 = (PProduct * F2Em * F1m * intensities).array();
117 const auto e4 = (PProduct * F2m * F1Em * intensities).array();
118 const auto sqPFProduct = (PFProduct.array() * PFProduct.array()).matrix();
119 const auto sqErrors = (yErrors.array() * yErrors.array()).matrix();
120 const auto e5 = (sqPFProduct * sqErrors).array();
121 errors = (e1 * e1 + e2 * e2 + e3 * e3 + e4 * e4 + e5).sqrt();
122}
123
140double twoInputsErrorEstimate01(const double i00, const double e00, const double i11, const double e11, const double p1,
141 const double p1E, const double p2, const double p2E, const double f1, const double f1E,
142 const double f2, const double f2E) {
143 using namespace boost::math;
144 // Derivatives of the equation which solves the I01 intensities
145 // with respect to i00, i11, f1, etc.
146 const auto pmdi00 =
147 -((f1 * (-1. + 2. * p1) * (-f2 * pow<2>(1. - 2. * p2) + pow<2>(f2) * pow<2>(1. - 2. * p2) + (-1. + p2) * p2)) /
148 (f2 * p1 * (-1. + p1 + 2. * p2 - 2. * p1 * p2) +
149 f1 * (-1. + 2. * p1) * ((1. - p2) * p2 + f2 * (-1. + p1 + p2) * (-1. + 2. * p2))));
150 const auto pmdi11 = (f2 * p1 * (-1. + p1 + 2. * p2 - 2. * p1 * p2)) /
151 (f2 * p1 * (-1. + p1 + 2. * p2 - 2. * p1 * p2) +
152 f1 * (-1. + 2. * p1) * ((1. - p2) * p2 + f2 * (-1. + p1 + p2) * (-1. + 2. * p2)));
153 const auto pmdf1 =
154 -(((-1. + 2. * p1) * ((1. - p2) * p2 + f2 * (-1. + p1 + p2) * (-1. + 2. * p2)) *
155 (f2 * i11 * p1 * (-1. + p1 + 2. * p2 - 2. * p1 * p2) -
156 f1 * i00 * (-1. + 2. * p1) *
157 (-f2 * pow<2>(1. - 2. * p2) + pow<2>(f2) * pow<2>(1. - 2. * p2) + (-1. + p2) * p2))) /
158 pow<2>(f2 * p1 * (-1. + p1 + 2. * p2 - 2. * p1 * p2) +
159 f1 * (-1. + 2. * p1) * ((1. - p2) * p2 + f2 * (-1. + p1 + p2) * (-1. + 2. * p2)))) -
160 (i00 * (-1. + 2. * p1) * (-f2 * pow<2>(1. - 2. * p2) + pow<2>(f2) * pow<2>(1. - 2. * p2) + (-1. + p2) * p2)) /
161 (f2 * p1 * (-1. + p1 + 2. * p2 - 2. * p1 * p2) +
162 f1 * (-1. + 2. * p1) * ((1. - p2) * p2 + f2 * (-1. + p1 + p2) * (-1. + 2. * p2)));
163 const auto pmdf2 =
164 -(((f1 * (-1. + 2. * p1) * (-1. + p1 + p2) * (-1 + 2 * p2) + p1 * (-1. + p1 + 2. * p2 - 2. * p1 * p2)) *
165 (f2 * i11 * p1 * (-1. + p1 + 2. * p2 - 2. * p1 * p2) -
166 f1 * i00 * (-1. + 2. * p1) *
167 (-f2 * pow<2>(1. - 2. * p2) + pow<2>(f2) * pow<2>(1. - 2. * p2) + (-1. + p2) * p2))) /
168 pow<2>(f2 * p1 * (-1. + p1 + 2. * p2 - 2. * p1 * p2) +
169 f1 * (-1. + 2. * p1) * ((1. - p2) * p2 + f2 * (-1. + p1 + p2) * (-1. + 2. * p2)))) +
170 (-f1 * i00 * (-1. + 2. * p1) * (-pow<2>(1. - 2. * p2) + 2 * f2 * pow<2>(1. - 2. * p2)) +
171 i11 * p1 * (-1. + p1 + 2. * p2 - 2. * p1 * p2)) /
172 (f2 * p1 * (-1. + p1 + 2. * p2 - 2. * p1 * p2) +
173 f1 * (-1. + 2. * p1) * ((1. - p2) * p2 + f2 * (-1. + p1 + p2) * (-1. + 2. * p2)));
174 const auto pmdp1 =
175 -(((f2 * i11 * p1 * (-1. + p1 + 2. * p2 - 2. * p1 * p2) -
176 f1 * i00 * (-1. + 2. * p1) *
177 (-f2 * pow<2>(1. - 2. * p2) + pow<2>(f2) * pow<2>(1. - 2. * p2) + (-1. + p2) * p2)) *
178 (f2 * p1 * (1. - 2. * p2) + f1 * f2 * (-1. + 2. * p1) * (-1. + 2. * p2) +
179 f2 * (-1. + p1 + 2. * p2 - 2. * p1 * p2) +
180 2. * f1 * ((1. - p2) * p2 + f2 * (-1. + p1 + p2) * (-1. + 2. * p2)))) /
181 pow<2>(f2 * p1 * (-1. + p1 + 2. * p2 - 2. * p1 * p2) +
182 f1 * (-1. + 2. * p1) * ((1. - p2) * p2 + f2 * (-1. + p1 + p2) * (-1. + 2. * p2)))) +
183 (f2 * i11 * p1 * (1. - 2. * p2) + f2 * i11 * (-1. + p1 + 2. * p2 - 2. * p1 * p2) -
184 2. * f1 * i00 * (-f2 * pow<2>(1. - 2. * p2) + pow<2>(f2) * pow<2>(1. - 2. * p2) + (-1. + p2) * p2)) /
185 (f2 * p1 * (-1. + p1 + 2. * p2 - 2. * p1 * p2) +
186 f1 * (-1. + 2. * p1) * ((1. - p2) * p2 + f2 * (-1. + p1 + p2) * (-1. + 2. * p2)));
187 const auto pmdp2 =
188 -(((f2 * i11 * p1 * (-1. + p1 + 2. * p2 - 2. * p1 * p2) -
189 f1 * i00 * (-1. + 2. * p1) *
190 (-f2 * pow<2>(1. - 2. * p2) + pow<2>(f2) * pow<2>(1. - 2. * p2) + (-1. + p2) * p2)) *
191 (f2 * (2. - 2. * p1) * p1 +
192 f1 * (-1. + 2. * p1) * (1. - 2. * p2 + 2. * f2 * (-1. + p1 + p2) + f2 * (-1. + 2. * p2)))) /
193 pow<2>(f2 * p1 * (-1. + p1 + 2. * p2 - 2. * p1 * p2) +
194 f1 * (-1. + 2. * p1) * ((1. - p2) * p2 + f2 * (-1. + p1 + p2) * (-1. + 2. * p2)))) +
195 (f2 * i11 * (2. - 2. * p1) * p1 -
196 f1 * i00 * (-1. + 2. * p1) * (-1. + 4. * f2 * (1. - 2. * p2) - 4. * pow<2>(f2) * (1. - 2. * p2) + 2. * p2)) /
197 (f2 * p1 * (-1. + p1 + 2. * p2 - 2. * p1 * p2) +
198 f1 * (-1. + 2. * p1) * ((1. - p2) * p2 + f2 * (-1. + p1 + p2) * (-1. + 2. * p2)));
199 // Estimate the error components using linearized extrapolation,
200 // sum in squares.
201 const auto e01_I00 = pow<2>(pmdi00 * e00);
202 const auto e01_I11 = pow<2>(pmdi11 * e11);
203 const auto e01_F1 = pow<2>(pmdf1 * f1E);
204 const auto e01_F2 = pow<2>(pmdf2 * f2E);
205 const auto e01_P1 = pow<2>(pmdp1 * p1E);
206 const auto e01_P2 = pow<2>(pmdp2 * p2E);
207 return std::sqrt(e01_I00 + e01_I11 + e01_F1 + e01_F2 + e01_P1 + e01_P2);
208}
209
226double twoInputsErrorEstimate10(const double i00, const double e00, const double i11, const double e11, const double p1,
227 const double p1E, const double p2, const double p2E, const double f1, const double f1E,
228 const double f2, const double f2E) {
229 using namespace boost::math;
230 // Derivatives of the equation which solves the I10 intensities
231 // with respect to i00, i11, f1, etc.
232 const auto a = -1. + p1 + 2. * p2 - 2. * p1 * p2;
233 const auto b = -1. + 2. * p1;
234 const auto c = -1. + 2. * p2;
235 const auto d = -1. + p2;
236 const auto mpdi00 = (-pow<2>(f1) * f2 * pow<2>(b) * c + f1 * f2 * pow<2>(b) * c + f2 * p1 * a) /
237 (f2 * p1 * a + f1 * b * (-d * p2 + f2 * (p1 + d) * c));
238 const auto mpdi11 = -((f1 * b * d * p2) / (f2 * p1 * a + f1 * b * (-d * p2 + f2 * (p1 + d) * c)));
239 const auto mpdf1 =
240 -(((-1. + 2. * p1) * ((1. - p2) * p2 + f2 * (-1. + p1 + p2) * (-1. + 2. * p2)) *
241 (-pow<2>(f1) * f2 * i00 * pow<2>(1. - 2. * p1) * (-1. + 2. * p2) +
242 f2 * i00 * p1 * (-1. + p1 + 2. * p2 - 2. * p1 * p2) +
243 f1 * (-1. + 2. * p1) * (-i11 * (-1. + p2) * p2 + f2 * i00 * (-1. + 2. * p1) * (-1. + 2. * p2)))) /
244 pow<2>(f2 * p1 * (-1. + p1 + 2. * p2 - 2. * p1 * p2) +
245 f1 * (-1. + 2. * p1) * ((1. - p2) * p2 + f2 * (-1. + p1 + p2) * (-1. + 2. * p2)))) +
246 (-2. * f1 * f2 * i00 * pow<2>(1. - 2. * p1) * (-1. + 2. * p2) +
247 (-1. + 2. * p1) * (-i11 * (-1. + p2) * p2 + f2 * i00 * (-1. + 2. * p1) * (-1. + 2. * p2))) /
248 (f2 * p1 * (-1. + p1 + 2. * p2 - 2. * p1 * p2) +
249 f1 * (-1. + 2. * p1) * ((1. - p2) * p2 + f2 * (-1. + p1 + p2) * (-1. + 2. * p2)));
250 const auto mpdf2 = -(((f1 * b * (p1 + d) * c + p1 * a) * (-pow<2>(f1) * f2 * i00 * pow<2>(b) * c + f2 * i00 * p1 * a +
251 f1 * b * (-i11 * d * p2 + f2 * i00 * b * c))) /
252 pow<2>(f2 * p1 * a + f1 * b * (-d * p2 + f2 * (p1 + d) * c))) +
253 (-pow<2>(f1) * i00 * pow<2>(b) * c + f1 * i00 * pow<2>(b) * c + i00 * p1 * a) /
254 (f2 * p1 * a + f1 * b * (-d * p2 + f2 * (p1 + d) * c));
255 const auto mpdp1 =
256 -(((-pow<2>(f1) * f2 * i00 * pow<2>(b) * c + f2 * i00 * p1 * a + f1 * b * (-i11 * d * p2 + f2 * i00 * b * c)) *
257 (f2 * p1 * -c + f1 * f2 * b * c + f2 * a + 2. * f1 * (-d * p2 + f2 * (p1 + d) * c))) /
258 pow<2>(f2 * p1 * a + f1 * b * (-d * p2 + f2 * (p1 + d) * c))) +
259 (f2 * i00 * p1 * -c + 4. * pow<2>(f1) * f2 * i00 * -b * c + 2. * f1 * f2 * i00 * b * c + f2 * i00 * a +
260 2. * f1 * (-i11 * d * p2 + f2 * i00 * b * c)) /
261 (f2 * p1 * a + f1 * b * (-d * p2 + f2 * (p1 + d) * c));
262 const auto mpdp2 =
263 -(((f2 * (2. - 2. * p1) * p1 + f1 * b * (1. - 2. * p2 + 2. * f2 * (p1 + d) + f2 * c)) *
264 (-pow<2>(f1) * f2 * i00 * pow<2>(b) * c + f2 * i00 * p1 * a + f1 * b * (-i11 * d * p2 + f2 * i00 * b * c))) /
265 pow<2>(f2 * p1 * a + f1 * b * (-d * p2 + f2 * (p1 + d) * c))) +
266 (-2. * pow<2>(f1) * f2 * i00 * pow<2>(b) + f2 * i00 * (2. - 2. * p1) * p1 +
267 f1 * b * (2. * f2 * i00 * b - i11 * d - i11 * p2)) /
268 (f2 * p1 * a + f1 * b * (-d * p2 + f2 * (p1 + d) * c));
269 // Estimate the error components using linearized extrapolation,
270 // sum in squares.
271 const auto e10_I00 = pow<2>(mpdi00 * e00);
272 const auto e10_I11 = pow<2>(mpdi11 * e11);
273 const auto e10_F1 = pow<2>(mpdf1 * f1E);
274 const auto e10_F2 = pow<2>(mpdf2 * f2E);
275 const auto e10_P1 = pow<2>(mpdp1 * p1E);
276 const auto e10_P2 = pow<2>(mpdp2 * p2E);
277 return std::sqrt(e10_I00 + e10_I11 + e10_F1 + e10_F2 + e10_P1 + e10_P2);
278}
279
280Mantid::API::MatrixWorkspace_sptr createWorkspaceWithHistory(const Mantid::API::MatrixWorkspace_const_sptr &inputWS) {
281 Mantid::API::MatrixWorkspace_sptr outputWS = Mantid::DataObjects::create<Mantid::DataObjects::Workspace2D>(*inputWS);
282 outputWS->history().addHistory(inputWS->getHistory());
283 return outputWS;
284}
285} // namespace
286
287namespace Mantid::Algorithms {
288
289// Register the algorithm into the AlgorithmFactory
290DECLARE_ALGORITHM(PolarizationCorrectionWildes)
291
292//----------------------------------------------------------------------------------------------
293
294
295const std::string PolarizationCorrectionWildes::name() const { return "PolarizationCorrectionWildes"; }
296
298int PolarizationCorrectionWildes::version() const { return 1; }
299
301const std::string PolarizationCorrectionWildes::category() const { return "Reflectometry"; }
302
304const std::string PolarizationCorrectionWildes::summary() const {
305 return "Corrects a group of polarization analysis workspaces for polarizer "
306 "and analyzer efficiencies.";
307}
308
310const std::vector<std::string> PolarizationCorrectionWildes::seeAlso() const {
311 return {"PolarizationEfficiencyCor", "PolarizationEfficienciesWildes"};
312}
313
318size_t PolarizationCorrectionWildes::WorkspaceMap::size() const noexcept {
319 return (mmWS ? 1 : 0) + (mpWS ? 1 : 0) + (pmWS ? 1 : 0) + (ppWS ? 1 : 0);
320}
321
322//----------------------------------------------------------------------------------------------
325void PolarizationCorrectionWildes::init() {
326 declareProperty(std::make_unique<Kernel::ArrayProperty<std::string>>(
327 Prop::INPUT_WS, "", std::make_shared<API::ADSValidator>(true, true), Kernel::Direction::Input),
328 "A list of workspaces to be corrected corresponding to the "
329 "flipper configurations.");
330 declareProperty(std::make_unique<API::WorkspaceProperty<API::WorkspaceGroup>>(
331 Prop::INPUT_WS_GROUP, "", Kernel::Direction::Input, API::PropertyMode::Optional),
332 "A group of workspaces to be corrected corresponding to the "
333 "flipper configurations.");
334 declareProperty(
335 std::make_unique<API::WorkspaceProperty<API::WorkspaceGroup>>(Prop::OUTPUT_WS, "", Kernel::Direction::Output),
336 "A group of polarization efficiency corrected workspaces.");
337
338 const auto flipperConfigValidator =
339 std::make_shared<Kernel::SpinStateValidator>(std::unordered_set<int>{1, 2, 3, 4}, true);
340 declareProperty(Prop::FLIPPERS,
341 std::string(FlipperConfigurations::OFF_OFF) + ", " + FlipperConfigurations::OFF_ON + ", " +
342 FlipperConfigurations::ON_OFF + ", " + FlipperConfigurations::ON_ON,
343 flipperConfigValidator, "Flipper configurations of the input workspaces.");
344 const auto spinStateValidator =
345 std::make_shared<Kernel::SpinStateValidator>(std::unordered_set<int>{0, 2, 4}, false, "+", "-", true);
346 declareProperty(Prop::SPIN_STATES, "", spinStateValidator, "The order of the spin states in the output workspace.");
347 declareProperty(
348 std::make_unique<API::WorkspaceProperty<API::MatrixWorkspace>>(Prop::EFFICIENCIES, "", Kernel::Direction::Input),
349 "A workspace containing the efficiency factors P1, P2, F1 and F2 as "
350 "histograms");
351 declareProperty(
352 Prop::ADD_SPIN_STATE_LOG, false,
353 "Whether to add the final spin state into the sample log of each child workspace in the output group.");
354}
355
356//----------------------------------------------------------------------------------------------
359void PolarizationCorrectionWildes::exec() {
360 const std::string flipperProperty = getProperty(Prop::FLIPPERS);
361 const auto flippers = Kernel::SpinStateHelpers::splitSpinStateString(flipperProperty);
362 // Check if the input flipper configuration includes an analyser
363 const bool hasAnalyser = flippers.front().size() > 1;
364 const auto inputs = mapInputsToDirections(flippers);
365 checkConsistentNumberHistograms(inputs);
366 const EfficiencyMap efficiencies = efficiencyFactors();
367 checkConsistentX(inputs, efficiencies);
368 WorkspaceMap outputs;
369 switch (inputs.size()) {
370 case 1:
371 outputs = directBeamCorrections(inputs, efficiencies);
372 break;
373 case 2:
374 if (hasAnalyser) {
375 outputs = twoInputCorrections(inputs, efficiencies);
376 } else {
377 outputs = analyzerlessCorrections(inputs, efficiencies);
378 }
379 break;
380 case 3:
381 outputs = threeInputCorrections(inputs, efficiencies);
382 break;
383 case 4:
384 outputs = fullCorrections(inputs, efficiencies);
385 }
386 setProperty(Prop::OUTPUT_WS, groupOutput(outputs, hasAnalyser));
387}
388
393std::map<std::string, std::string> PolarizationCorrectionWildes::validateInputs() {
394 std::map<std::string, std::string> issues;
395 API::MatrixWorkspace_const_sptr factorWS = getProperty(Prop::EFFICIENCIES);
396 if (factorWS) {
397 const auto &factorAxis = factorWS->getAxis(1);
398 if (!factorAxis) {
399 issues[Prop::EFFICIENCIES] = "The workspace is missing a vertical axis.";
400 } else if (!factorAxis->isText()) {
401 issues[Prop::EFFICIENCIES] = "The vertical axis in the workspace is not text axis.";
402 } else if (factorWS->getNumberHistograms() < 4) {
403 issues[Prop::EFFICIENCIES] = "The workspace should contain at least 4 histograms.";
404 } else {
405 std::vector<std::string> tags{{"P1", "P2", "F1", "F2"}};
406 for (size_t i = 0; i != factorAxis->length(); ++i) {
407 const auto label = factorAxis->label(i);
408 auto found = std::find(tags.begin(), tags.end(), label);
409 if (found != tags.cend()) {
410 std::swap(tags.back(), *found);
411 tags.pop_back();
412 }
413 }
414 if (!tags.empty()) {
415 issues[Prop::EFFICIENCIES] = "A histogram labeled " + tags.front() + " is missing from the workspace.";
416 }
417 }
418 }
419 const std::vector<std::string> inputs = getProperty(Prop::INPUT_WS);
420 const API::WorkspaceGroup_sptr inputsGroup = getProperty(Prop::INPUT_WS_GROUP);
421 if (inputs.empty() && !inputsGroup) {
422 issues[Prop::INPUT_WS] = "At least one of the input workspace properties must be provided.";
423 } else if (!inputs.empty() && inputsGroup) {
424 issues[Prop::INPUT_WS] = "Only one of the input workspace properties can be provided.";
425 }
426 const auto flipperConfig = Kernel::SpinStateHelpers::splitSpinStateString(getPropertyValue(Prop::FLIPPERS));
427 const auto flipperCount = flipperConfig.size();
428 const size_t inputsCount = (inputsGroup) ? inputsGroup->getNumberOfEntries() : inputs.size();
429 if (inputsCount != flipperCount) {
430 issues[Prop::FLIPPERS] = "The number of flipper configurations (" + std::to_string(flipperCount) +
431 ") does not match the number of input workspaces (" + std::to_string(inputsCount) + ")";
432 }
433 // SpinStates checks.
434 const auto spinStates = Kernel::SpinStateHelpers::splitSpinStateString(getPropertyValue(Prop::SPIN_STATES));
435 if (inputsCount == 1 && !spinStates.empty()) {
436 issues[Prop::SPIN_STATES] = "Output workspace order cannot be set for direct beam calculations.";
437 } else if (!spinStates.empty()) {
438 if (flipperConfig.front().size() == 1 && spinStates.size() != 2) {
439 issues[Prop::SPIN_STATES] =
440 "Incorrect number of workspaces in output configuration: " + std::to_string(spinStates.size()) +
441 ". Only two output workspaces are produced when an analyzer is not used.";
442 }
443 if (flipperConfig.front().size() == 2 && spinStates.size() != 4) {
444 issues[Prop::SPIN_STATES] =
445 "Incorrect number of workspaces in output configuration: " + std::to_string(spinStates.size()) +
446 ". Four output workspaces are produced by the corrections.";
447 }
448 }
449 return issues;
450}
451
456void PolarizationCorrectionWildes::checkConsistentNumberHistograms(const WorkspaceMap &inputs) {
457 size_t nHist{0};
458 bool nHistValid{false};
459 // A local helper function to check the number of histograms.
460 auto checkNHist = [&nHist, &nHistValid](const API::MatrixWorkspace_sptr &ws, const std::string &tag) {
461 if (nHistValid) {
462 if (nHist != ws->getNumberHistograms()) {
463 throw std::runtime_error("Number of histograms mismatch in " + tag);
464 }
465 } else {
466 nHist = ws->getNumberHistograms();
467 nHistValid = true;
468 }
469 };
470 if (inputs.mmWS) {
471 checkNHist(inputs.mmWS, FlipperConfigurations::ON_ON);
472 }
473 if (inputs.mpWS) {
474 checkNHist(inputs.mpWS, FlipperConfigurations::ON_OFF);
475 }
476 if (inputs.pmWS) {
477 checkNHist(inputs.pmWS, FlipperConfigurations::OFF_ON);
478 }
479 if (inputs.ppWS) {
480 checkNHist(inputs.ppWS, FlipperConfigurations::OFF_OFF);
481 }
482}
483
489void PolarizationCorrectionWildes::checkConsistentX(const WorkspaceMap &inputs, const EfficiencyMap &efficiencies) {
490 // Compare everything to F1 efficiency.
491 const auto &F1x = efficiencies.F1->x();
492 // A local helper function to check a HistogramX against F1.
493 auto checkX = [&F1x](const HistogramData::HistogramX &x, const std::string &tag) {
494 if (x.size() != F1x.size()) {
495 throw std::runtime_error("Mismatch of histogram lengths between F1 and " + tag + '.');
496 }
497 for (size_t i = 0; i != x.size(); ++i) {
498 if (x[i] != F1x[i]) {
499 throw std::runtime_error("Mismatch of X data between F1 and " + tag + '.');
500 }
501 }
502 };
503 const auto &F2x = efficiencies.F2->x();
504 checkX(F2x, "F2");
505 const auto &P1x = efficiencies.P1->x();
506 checkX(P1x, "P1");
507 const auto &P2x = efficiencies.P2->x();
508 checkX(P2x, "P2");
509 // A local helper function to check an input workspace against F1.
510 auto checkWS = [&checkX](const API::MatrixWorkspace_sptr &ws, const std::string &tag) {
511 const auto nHist = ws->getNumberHistograms();
512 for (size_t i = 0; i != nHist; ++i) {
513 checkX(ws->x(i), tag);
514 }
515 };
516 if (inputs.mmWS) {
517 checkWS(inputs.mmWS, FlipperConfigurations::ON_ON);
518 }
519 if (inputs.mpWS) {
520 checkWS(inputs.mpWS, FlipperConfigurations::ON_OFF);
521 }
522 if (inputs.pmWS) {
523 checkWS(inputs.pmWS, FlipperConfigurations::OFF_ON);
524 }
525 if (inputs.ppWS) {
526 checkWS(inputs.ppWS, FlipperConfigurations::OFF_OFF);
527 }
528}
529
538API::WorkspaceGroup_sptr PolarizationCorrectionWildes::groupOutput(const WorkspaceMap &outputs,
539 const bool hasAnalyser) {
540 const auto &outWSName = getPropertyValue(Prop::OUTPUT_WS);
541 auto spinStateOrder = getPropertyValue(Prop::SPIN_STATES);
542 const bool addSpinStateLog = getProperty(Prop::ADD_SPIN_STATE_LOG);
543
544 std::vector<std::string> names;
545 if (!spinStateOrder.empty()) {
546 names.resize(Kernel::SpinStateHelpers::splitSpinStateString(spinStateOrder).size());
547 }
548
549 if (outputs.ppWS) {
550 addSpinStateOutput(names, spinStateOrder, outWSName, outputs.ppWS, SpinStateConfigurationsWildes::PLUS_PLUS,
551 addSpinStateLog, hasAnalyser);
552 }
553 if (outputs.pmWS) {
554 addSpinStateOutput(names, spinStateOrder, outWSName, outputs.pmWS, SpinStateConfigurationsWildes::PLUS_MINUS,
555 addSpinStateLog, hasAnalyser);
556 }
557 if (outputs.mpWS) {
558 addSpinStateOutput(names, spinStateOrder, outWSName, outputs.mpWS, SpinStateConfigurationsWildes::MINUS_PLUS,
559 addSpinStateLog, hasAnalyser);
560 }
561 if (outputs.mmWS) {
562 addSpinStateOutput(names, spinStateOrder, outWSName, outputs.mmWS, SpinStateConfigurationsWildes::MINUS_MINUS,
563 addSpinStateLog, hasAnalyser);
564 }
565
566 auto group = createChildAlgorithm("GroupWorkspaces");
567 group->initialize();
568 group->setProperty("InputWorkspaces", names);
569 group->setProperty("OutputWorkspace", outWSName);
570 group->execute();
571 API::WorkspaceGroup_sptr outWS = group->getProperty("OutputWorkspace");
572 return outWS;
573}
574
587void PolarizationCorrectionWildes::addSpinStateOutput(std::vector<std::string> &names,
588 const std::string &spinStateOrder, const std::string &baseName,
589 const API::MatrixWorkspace_sptr &ws, const std::string &spinState,
590 const bool addSpinStateLog, const bool hasAnalyser) {
591 if (addSpinStateLog) {
592 addSpinStateLogToWs(ws, spinState, hasAnalyser);
593 }
594
595 if (spinStateOrder.empty()) {
596 names.emplace_back(baseName + "_" + spinState);
597 API::AnalysisDataService::Instance().addOrReplace(names.back(), ws);
598 } else {
599 const auto &maybeIndex = Kernel::SpinStateHelpers::indexOfWorkspaceForSpinState(
600 Kernel::SpinStateHelpers::splitSpinStateString(spinStateOrder), spinState);
601 if (!maybeIndex.has_value()) {
602 throw std::invalid_argument("Required spin state (" + spinState + ") not found in spin state order (" +
603 spinStateOrder + ").");
604 }
605 const auto index = maybeIndex.value();
606 names[index] = baseName + "_" + spinState;
607 API::AnalysisDataService::Instance().addOrReplace(names[index], ws);
608 }
609}
610
617void PolarizationCorrectionWildes::addSpinStateLogToWs(const API::MatrixWorkspace_sptr &ws,
618 const std::string &spinState, const bool hasAnalyser) {
619 if (!hasAnalyser) {
620 if (spinState == SpinStateConfigurationsWildes::PLUS_PLUS) {
621 SpinStatesORSO::addORSOLogForSpinState(ws, SpinStateConfigurationsWildes::PLUS);
622 return;
623 }
624
625 if (spinState == SpinStateConfigurationsWildes::MINUS_MINUS) {
626 SpinStatesORSO::addORSOLogForSpinState(ws, SpinStateConfigurationsWildes::MINUS);
627 return;
628 }
629 }
630
631 SpinStatesORSO::addORSOLogForSpinState(ws, spinState);
632}
633
638PolarizationCorrectionWildes::EfficiencyMap PolarizationCorrectionWildes::efficiencyFactors() {
639 EfficiencyMap e;
640 API::MatrixWorkspace_const_sptr factorWS = getProperty(Prop::EFFICIENCIES);
641 const auto &vertAxis = factorWS->getAxis(1);
642 for (size_t i = 0; i != vertAxis->length(); ++i) {
643 const auto label = vertAxis->label(i);
644 if (label == "P1") {
645 e.P1 = &factorWS->getSpectrum(i);
646 } else if (label == "P2") {
647 e.P2 = &factorWS->getSpectrum(i);
648 } else if (label == "F1") {
649 e.F1 = &factorWS->getSpectrum(i);
650 } else if (label == "F2") {
651 e.F2 = &factorWS->getSpectrum(i);
652 }
653 // Ignore other histograms such as 'Phi' in ILL's efficiency ws.
654 }
655 return e;
656}
657
665PolarizationCorrectionWildes::WorkspaceMap
666PolarizationCorrectionWildes::directBeamCorrections(const WorkspaceMap &inputs, const EfficiencyMap &efficiencies) {
667 using namespace boost::math;
668 checkInputExists(inputs.ppWS, FlipperConfigurations::OFF);
669 WorkspaceMap outputs;
670 outputs.ppWS = createWorkspaceWithHistory(inputs.ppWS);
671 const size_t nHisto = inputs.ppWS->getNumberHistograms();
672 // cppcheck-suppress unreadVariable
673 const bool threadSafe = Kernel::threadSafe(*inputs.ppWS, *outputs.ppWS);
674 PARALLEL_FOR_IF(threadSafe)
675 for (int wsIndex = 0; wsIndex < static_cast<int>(nHisto); ++wsIndex) {
676 PARALLEL_START_INTERRUPT_REGION
677 const auto &ppY = inputs.ppWS->y(wsIndex);
678 const auto &ppE = inputs.ppWS->e(wsIndex);
679 auto &ppYOut = outputs.ppWS->mutableY(wsIndex);
680 auto &ppEOut = outputs.ppWS->mutableE(wsIndex);
681 PARALLEL_FOR_IF(threadSafe)
682 for (int binIndex = 0; binIndex < static_cast<int>(ppY.size()); ++binIndex) {
683 const auto P1 = efficiencies.P1->y()[binIndex];
684 const auto P2 = efficiencies.P2->y()[binIndex];
685 const double f = 1. - P1 - P2 + 2. * P1 * P2;
686 ppYOut[binIndex] = ppY[binIndex] / f;
687 const auto P1E = efficiencies.P1->e()[binIndex];
688 const auto P2E = efficiencies.P2->e()[binIndex];
689 const auto e1 = pow<2>(P1E * (2. * P1 - 1.) / pow<2>(f) * ppY[binIndex]);
690 const auto e2 = pow<2>(P2E * (2. * P2 - 1.) / pow<2>(f) * ppY[binIndex]);
691 const auto e3 = pow<2>(ppE[binIndex] / f);
692 const auto errorSum = std::sqrt(e1 + e2 + e3);
693 ppEOut[binIndex] = errorSum;
694 }
695 PARALLEL_END_INTERRUPT_REGION
696 }
697 PARALLEL_CHECK_INTERRUPT_REGION
698 return outputs;
699}
700
709PolarizationCorrectionWildes::WorkspaceMap
710PolarizationCorrectionWildes::analyzerlessCorrections(const WorkspaceMap &inputs, const EfficiencyMap &efficiencies) {
711 using namespace boost::math;
712 checkInputExists(inputs.mmWS, FlipperConfigurations::ON);
713 checkInputExists(inputs.ppWS, FlipperConfigurations::OFF);
714 WorkspaceMap outputs;
715 outputs.mmWS = createWorkspaceWithHistory(inputs.mmWS);
716 outputs.ppWS = createWorkspaceWithHistory(inputs.ppWS);
717 const size_t nHisto = inputs.mmWS->getNumberHistograms();
718 // cppcheck-suppress unreadVariable
719 const bool threadSafe = Kernel::threadSafe(*inputs.mmWS, *inputs.ppWS, *outputs.mmWS, *outputs.ppWS);
720 PARALLEL_FOR_IF(threadSafe)
721 for (int wsIndex = 0; wsIndex < static_cast<int>(nHisto); ++wsIndex) {
722 PARALLEL_START_INTERRUPT_REGION
723 const auto &mmY = inputs.mmWS->y(wsIndex);
724 const auto &mmE = inputs.mmWS->e(wsIndex);
725 const auto &ppY = inputs.ppWS->y(wsIndex);
726 const auto &ppE = inputs.ppWS->e(wsIndex);
727 auto &mmYOut = outputs.mmWS->mutableY(wsIndex);
728 auto &mmEOut = outputs.mmWS->mutableE(wsIndex);
729 auto &ppYOut = outputs.ppWS->mutableY(wsIndex);
730 auto &ppEOut = outputs.ppWS->mutableE(wsIndex);
731 PARALLEL_FOR_IF(threadSafe)
732 for (int binIndex = 0; binIndex < static_cast<int>(mmY.size()); ++binIndex) {
733 const auto F1 = efficiencies.F1->y()[binIndex];
734 const auto P1 = efficiencies.P1->y()[binIndex];
735 Eigen::Matrix2d F1m;
736 F1m << 1., 0., (F1 - 1.) / F1, 1. / F1;
737 const double divisor = (2. * P1 - 1.);
738 const double off = (P1 - 1.) / divisor;
739 const double diag = P1 / divisor;
740 Eigen::Matrix2d P1m;
741 P1m << diag, off, off, diag;
742 const Eigen::Vector2d intensities(ppY[binIndex], mmY[binIndex]);
743 const auto PFProduct = P1m * F1m;
744 const auto corrected = PFProduct * intensities;
745 ppYOut[binIndex] = corrected[0];
746 mmYOut[binIndex] = corrected[1];
747 const auto F1E = efficiencies.F1->e()[binIndex];
748 const auto P1E = efficiencies.P1->e()[binIndex];
749 const auto elemE1 = -1. / pow<2>(F1) * F1E;
750 Eigen::Matrix2d F1Em;
751 F1Em << 0., 0., -elemE1, elemE1;
752 const auto elemE2 = 1. / pow<2>(divisor) * P1E;
753 Eigen::Matrix2d P1Em;
754 P1Em << elemE2, -elemE2, -elemE2, elemE2;
755 const Eigen::Vector2d errors(ppE[binIndex], mmE[binIndex]);
756 const auto e1 = (P1Em * F1m * intensities).array();
757 const auto e2 = (P1m * F1Em * intensities).array();
758 const auto sqPFProduct = (PFProduct.array() * PFProduct.array()).matrix();
759 const auto sqErrors = (errors.array() * errors.array()).matrix();
760 const auto e3 = (sqPFProduct * sqErrors).array();
761 const auto errorSum = (e1 * e1 + e2 * e2 + e3).sqrt();
762 ppEOut[binIndex] = errorSum[0];
763 mmEOut[binIndex] = errorSum[1];
764 }
765 PARALLEL_END_INTERRUPT_REGION
766 }
767 PARALLEL_CHECK_INTERRUPT_REGION
768 return outputs;
769}
770
780PolarizationCorrectionWildes::WorkspaceMap
781PolarizationCorrectionWildes::twoInputCorrections(const WorkspaceMap &inputs, const EfficiencyMap &efficiencies) {
782 using namespace boost::math;
783 checkInputExists(inputs.mmWS, FlipperConfigurations::ON_ON);
784 checkInputExists(inputs.ppWS, FlipperConfigurations::OFF_OFF);
785 WorkspaceMap fullInputs = inputs;
786 fullInputs.mpWS = createWorkspaceWithHistory(inputs.mmWS);
787 fullInputs.pmWS = createWorkspaceWithHistory(inputs.ppWS);
788 twoInputsSolve01And10(fullInputs, inputs, efficiencies);
789 return fullCorrections(fullInputs, efficiencies);
790}
791
801PolarizationCorrectionWildes::WorkspaceMap
802PolarizationCorrectionWildes::threeInputCorrections(const WorkspaceMap &inputs, const EfficiencyMap &efficiencies) {
803 WorkspaceMap fullInputs = inputs;
804 checkInputExists(inputs.mmWS, FlipperConfigurations::ON_ON);
805 checkInputExists(inputs.ppWS, FlipperConfigurations::OFF_OFF);
806 if (!inputs.mpWS) {
807 checkInputExists(inputs.pmWS, FlipperConfigurations::OFF_ON);
808 threeInputsSolve10(fullInputs, efficiencies);
809 } else {
810 checkInputExists(inputs.mpWS, FlipperConfigurations::ON_OFF);
811 threeInputsSolve01(fullInputs, efficiencies);
812 }
813 return fullCorrections(fullInputs, efficiencies);
814}
815
824PolarizationCorrectionWildes::WorkspaceMap
825PolarizationCorrectionWildes::fullCorrections(const WorkspaceMap &inputs, const EfficiencyMap &efficiencies) {
826 using namespace boost::math;
827 checkInputExists(inputs.mmWS, FlipperConfigurations::ON_ON);
828 checkInputExists(inputs.mpWS, FlipperConfigurations::ON_OFF);
829 checkInputExists(inputs.pmWS, FlipperConfigurations::OFF_ON);
830 checkInputExists(inputs.ppWS, FlipperConfigurations::OFF_OFF);
831 WorkspaceMap outputs;
832 outputs.mmWS = createWorkspaceWithHistory(inputs.mmWS);
833 outputs.mpWS = createWorkspaceWithHistory(inputs.mpWS);
834 outputs.pmWS = createWorkspaceWithHistory(inputs.pmWS);
835 outputs.ppWS = createWorkspaceWithHistory(inputs.ppWS);
836 const auto F1 = efficiencies.F1->y();
837 const auto F1E = efficiencies.F1->e();
838 const auto F2 = efficiencies.F2->y();
839 const auto F2E = efficiencies.F2->e();
840 const auto P1 = efficiencies.P1->y();
841 const auto P1E = efficiencies.P1->e();
842 const auto P2 = efficiencies.P2->y();
843 const auto P2E = efficiencies.P2->e();
844 const size_t nHisto = inputs.mmWS->getNumberHistograms();
845 // cppcheck-suppress unreadVariable
846 const bool threadSafe = Kernel::threadSafe(*inputs.mmWS, *inputs.mpWS, *inputs.pmWS, *inputs.ppWS, *outputs.mmWS,
847 *outputs.mpWS, *outputs.pmWS, *outputs.ppWS);
848 PARALLEL_FOR_IF(threadSafe)
849 for (int wsIndex = 0; wsIndex < static_cast<int>(nHisto); ++wsIndex) {
850 PARALLEL_START_INTERRUPT_REGION
851 const auto &mmY = inputs.mmWS->y(wsIndex);
852 const auto &mmE = inputs.mmWS->e(wsIndex);
853 const auto &mpY = inputs.mpWS->y(wsIndex);
854 const auto &mpE = inputs.mpWS->e(wsIndex);
855 const auto &pmY = inputs.pmWS->y(wsIndex);
856 const auto &pmE = inputs.pmWS->e(wsIndex);
857 const auto &ppY = inputs.ppWS->y(wsIndex);
858 const auto &ppE = inputs.ppWS->e(wsIndex);
859 auto &mmYOut = outputs.mmWS->mutableY(wsIndex);
860 auto &mmEOut = outputs.mmWS->mutableE(wsIndex);
861 auto &mpYOut = outputs.mpWS->mutableY(wsIndex);
862 auto &mpEOut = outputs.mpWS->mutableE(wsIndex);
863 auto &pmYOut = outputs.pmWS->mutableY(wsIndex);
864 auto &pmEOut = outputs.pmWS->mutableE(wsIndex);
865 auto &ppYOut = outputs.ppWS->mutableY(wsIndex);
866 auto &ppEOut = outputs.ppWS->mutableE(wsIndex);
867 PARALLEL_FOR_IF(threadSafe)
868 for (int binIndex = 0; binIndex < static_cast<int>(mmY.size()); ++binIndex) {
869 Eigen::Vector4d corrected;
870 Eigen::Vector4d errors;
871 fourInputsCorrectedAndErrors(corrected, errors, ppY[binIndex], ppE[binIndex], pmY[binIndex], pmE[binIndex],
872 mpY[binIndex], mpE[binIndex], mmY[binIndex], mmE[binIndex], F1[binIndex],
873 F1E[binIndex], F2[binIndex], F2E[binIndex], P1[binIndex], P1E[binIndex],
874 P2[binIndex], P2E[binIndex]);
875 ppYOut[binIndex] = corrected[0];
876 pmYOut[binIndex] = corrected[1];
877 mpYOut[binIndex] = corrected[2];
878 mmYOut[binIndex] = corrected[3];
879 ppEOut[binIndex] = errors[0];
880 pmEOut[binIndex] = errors[1];
881 mpEOut[binIndex] = errors[2];
882 mmEOut[binIndex] = errors[3];
883 }
884 PARALLEL_END_INTERRUPT_REGION
885 }
886 PARALLEL_CHECK_INTERRUPT_REGION
887 return outputs;
888}
889
895PolarizationCorrectionWildes::WorkspaceMap
896PolarizationCorrectionWildes::mapInputsToDirections(const std::vector<std::string> &flippers) {
897 const std::vector<std::string> inputNames = getProperty(Prop::INPUT_WS);
898 const API::WorkspaceGroup_sptr inputGroup = getProperty(Prop::INPUT_WS_GROUP);
899 WorkspaceMap inputs;
900 for (size_t i = 0; i < flippers.size(); ++i) {
901 auto ws = (inputGroup) ? dynamic_pointer_cast<API::MatrixWorkspace>(inputGroup->getItem(i))
902 : API::AnalysisDataService::Instance().retrieveWS<API::MatrixWorkspace>(inputNames[i]);
903 if (!ws) {
904 throw std::runtime_error("One of the input workspaces doesn't seem to be a MatrixWorkspace.");
905 }
906 const auto &f = flippers[i];
907 if (f == FlipperConfigurations::ON_ON || f == FlipperConfigurations::ON) {
908 inputs.mmWS = ws;
909 } else if (f == FlipperConfigurations::ON_OFF) {
910 inputs.mpWS = ws;
911 } else if (f == FlipperConfigurations::OFF_ON) {
912 inputs.pmWS = ws;
913 } else if (f == FlipperConfigurations::OFF_OFF || f == FlipperConfigurations::OFF) {
914 inputs.ppWS = ws;
915 } else {
916 throw std::runtime_error(std::string{"Unknown entry in "} + Prop::FLIPPERS);
917 }
918 }
919 return inputs;
920}
921
928void PolarizationCorrectionWildes::threeInputsSolve01(WorkspaceMap &inputs, const EfficiencyMap &efficiencies) {
929 inputs.pmWS = createWorkspaceWithHistory(inputs.mpWS);
930 const auto &F1 = efficiencies.F1->y();
931 const auto &F2 = efficiencies.F2->y();
932 const auto &P1 = efficiencies.P1->y();
933 const auto &P2 = efficiencies.P2->y();
934 const auto nHisto = inputs.pmWS->getNumberHistograms();
935 // cppcheck-suppress unreadVariable
936 const bool threadSafe = Kernel::threadSafe(*inputs.mmWS, *inputs.mpWS, *inputs.pmWS, *inputs.ppWS);
937 PARALLEL_FOR_IF(threadSafe)
938 for (int wsIndex = 0; wsIndex < static_cast<int>(nHisto); ++wsIndex) {
939 PARALLEL_START_INTERRUPT_REGION
940 const auto &I00 = inputs.ppWS->y(wsIndex);
941 auto &I01 = inputs.pmWS->mutableY(wsIndex);
942 const auto &I10 = inputs.mpWS->y(wsIndex);
943 const auto &I11 = inputs.mmWS->y(wsIndex);
944 PARALLEL_FOR_IF(threadSafe)
945 for (int binIndex = 0; binIndex < static_cast<int>(I00.size()); ++binIndex) {
946 const auto f1 = F1[binIndex];
947 const auto f2 = F2[binIndex];
948 const auto p1 = P1[binIndex];
949 const auto p2 = P2[binIndex];
950 const auto i00 = I00[binIndex];
951 const auto i10 = I10[binIndex];
952 const auto i11 = I11[binIndex];
953 I01[binIndex] =
954 (f1 * i00 * (-1. + 2. * p1) - (i00 - i10 + i11) * (p1 - p2) - f2 * (i00 - i10) * (-1. + 2. * p2)) /
955 (-p1 + f1 * (-1. + 2. * p1) + p2);
956 // The errors are left to zero.
957 }
958 PARALLEL_END_INTERRUPT_REGION
959 }
960 PARALLEL_CHECK_INTERRUPT_REGION
961}
962
969void PolarizationCorrectionWildes::threeInputsSolve10(WorkspaceMap &inputs, const EfficiencyMap &efficiencies) {
970 inputs.mpWS = createWorkspaceWithHistory(inputs.pmWS);
971 const auto &F1 = efficiencies.F1->y();
972 const auto &F2 = efficiencies.F2->y();
973 const auto &P1 = efficiencies.P1->y();
974 const auto &P2 = efficiencies.P2->y();
975 const auto nHisto = inputs.mpWS->getNumberHistograms();
976 // cppcheck-suppress unreadVariable
977 const bool threadSafe = Kernel::threadSafe(*inputs.mmWS, *inputs.mpWS, *inputs.pmWS, *inputs.ppWS);
978 PARALLEL_FOR_IF(threadSafe)
979 for (int wsIndex = 0; wsIndex < static_cast<int>(nHisto); ++wsIndex) {
980 PARALLEL_START_INTERRUPT_REGION
981 const auto &I00 = inputs.ppWS->y(wsIndex);
982 const auto &I01 = inputs.pmWS->y(wsIndex);
983 auto &I10 = inputs.mpWS->mutableY(wsIndex);
984 const auto &I11 = inputs.mmWS->y(wsIndex);
985 PARALLEL_FOR_IF(threadSafe)
986 for (int binIndex = 0; binIndex < static_cast<int>(I00.size()); ++binIndex) {
987 const auto f1 = F1[binIndex];
988 const auto f2 = F2[binIndex];
989 const auto p1 = P1[binIndex];
990 const auto p2 = P2[binIndex];
991 const auto i00 = I00[binIndex];
992 const auto i01 = I01[binIndex];
993 const auto i11 = I11[binIndex];
994 I10[binIndex] =
995 (-f1 * (i00 - i01) * (-1. + 2. * p1) + (i00 - i01 + i11) * (p1 - p2) + f2 * i00 * (-1. + 2. * p2)) /
996 (p1 - p2 + f2 * (-1. + 2. * p2));
997 // The errors are left to zero.
998 }
999 PARALLEL_END_INTERRUPT_REGION
1000 }
1001 PARALLEL_CHECK_INTERRUPT_REGION
1002}
1003
1011void PolarizationCorrectionWildes::twoInputsSolve01And10(WorkspaceMap &fullInputs, const WorkspaceMap &inputs,
1012 const EfficiencyMap &efficiencies) {
1013 using namespace boost::math;
1014 const auto &F1 = efficiencies.F1->y();
1015 const auto &F1E = efficiencies.F1->e();
1016 const auto &F2 = efficiencies.F2->y();
1017 const auto &F2E = efficiencies.F2->e();
1018 const auto &P1 = efficiencies.P1->y();
1019 const auto &P1E = efficiencies.P1->e();
1020 const auto &P2 = efficiencies.P2->y();
1021 const auto &P2E = efficiencies.P2->e();
1022 const auto nHisto = inputs.mmWS->getNumberHistograms();
1023 // cppcheck-suppress unreadVariable
1024 const bool threadSafe = Kernel::threadSafe(*inputs.mmWS, *inputs.ppWS, *fullInputs.pmWS, *fullInputs.mpWS);
1025 PARALLEL_FOR_IF(threadSafe)
1026 for (int wsIndex = 0; wsIndex < static_cast<int>(nHisto); ++wsIndex) {
1027 PARALLEL_START_INTERRUPT_REGION
1028 const auto &I00 = inputs.ppWS->y(wsIndex);
1029 const auto &E00 = inputs.ppWS->e(wsIndex);
1030 const auto &I11 = inputs.mmWS->y(wsIndex);
1031 const auto &E11 = inputs.mmWS->e(wsIndex);
1032 auto &I01 = fullInputs.pmWS->mutableY(wsIndex);
1033 auto &E01 = fullInputs.pmWS->mutableE(wsIndex);
1034 auto &I10 = fullInputs.mpWS->mutableY(wsIndex);
1035 auto &E10 = fullInputs.mpWS->mutableE(wsIndex);
1036 PARALLEL_FOR_IF(threadSafe)
1037 for (int binIndex = 0; binIndex < static_cast<int>(I00.size()); ++binIndex) {
1038 const auto i00 = I00[binIndex];
1039 const auto i11 = I11[binIndex];
1040 const auto f1 = F1[binIndex];
1041 const auto f2 = F2[binIndex];
1042 const auto p1 = P1[binIndex];
1043 const auto p2 = P2[binIndex];
1044 const auto a = -1. + p1 + 2. * p2 - 2. * p1 * p2;
1045 const auto b = -1. + 2. * p1;
1046 const auto c = -1. + 2. * p2;
1047 const auto d = -1. + p2;
1048 // Case: 01
1049 const auto divisor = f2 * p1 * a + f1 * b * (-d * p2 + f2 * (p1 + d) * c);
1050 I01[binIndex] = (f2 * i11 * p1 * a - f1 * i00 * b * (-f2 * pow<2>(c) + pow<2>(f2 * c) + d * p2)) / divisor;
1051 E01[binIndex] = twoInputsErrorEstimate01(i00, E00[binIndex], i11, E11[binIndex], p1, P1E[binIndex], p2,
1052 P2E[binIndex], f1, F1E[binIndex], f2, F2E[binIndex]);
1053 // Case: 10
1054 I10[binIndex] =
1055 (-pow<2>(f1) * f2 * i00 * pow<2>(b) * c + f2 * i00 * p1 * a + f1 * b * (-i11 * d * p2 + f2 * i00 * b * c)) /
1056 divisor;
1057 E10[binIndex] = twoInputsErrorEstimate10(i00, E00[binIndex], i11, E11[binIndex], p1, P1E[binIndex], p2,
1058 P2E[binIndex], f1, F1E[binIndex], f2, F2E[binIndex]);
1059 }
1060 PARALLEL_END_INTERRUPT_REGION
1061 }
1062 PARALLEL_CHECK_INTERRUPT_REGION
1063}
1064} // 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:93
Mantid::API::WorkspaceProperty
A property class for workspaces.
Definition WorkspaceProperty.h:58
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:489
Mantid::Algorithms::PolarizationCorrectionWildes::validateInputs
std::map< std::string, std::string > validateInputs() override
Validate the algorithm's input properties.
Definition PolarizationCorrectionWildes.cpp:393
Mantid::Algorithms::PolarizationCorrectionWildes::efficiencyFactors
EfficiencyMap efficiencyFactors()
Make a convenience access object to the efficiency factors.
Definition PolarizationCorrectionWildes.cpp:638
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:969
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:587
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:666
Mantid::Algorithms::PolarizationCorrectionWildes::init
void init() override
Initialize the algorithm's properties.
Definition PolarizationCorrectionWildes.cpp:325
Mantid::Algorithms::PolarizationCorrectionWildes::version
int version() const override
Algorithm's version for identification.
Definition PolarizationCorrectionWildes.cpp:298
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:617
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:928
Mantid::Algorithms::PolarizationCorrectionWildes::fullCorrections
WorkspaceMap fullCorrections(const WorkspaceMap &inputs, const EfficiencyMap &efficiencies)
Correct for non-ideal instrument effects.
Definition PolarizationCorrectionWildes.cpp:825
Mantid::Algorithms::PolarizationCorrectionWildes::seeAlso
const std::vector< std::string > seeAlso() const override
Algorithm's related algorithms.
Definition PolarizationCorrectionWildes.cpp:310
Mantid::Algorithms::PolarizationCorrectionWildes::category
const std::string category() const override
Algorithm's category for identification.
Definition PolarizationCorrectionWildes.cpp:301
Mantid::Algorithms::PolarizationCorrectionWildes::summary
const std::string summary() const override
Algorithm's summary for use in the GUI and help.
Definition PolarizationCorrectionWildes.cpp:304
Mantid::Algorithms::PolarizationCorrectionWildes::checkConsistentNumberHistograms
void checkConsistentNumberHistograms(const WorkspaceMap &inputs)
Check that all workspaces in inputs have the same number of histograms.
Definition PolarizationCorrectionWildes.cpp:456
Mantid::Algorithms::PolarizationCorrectionWildes::threeInputCorrections
WorkspaceMap threeInputCorrections(const WorkspaceMap &inputs, const EfficiencyMap &efficiencies)
Correct for non-ideal instrument effects.
Definition PolarizationCorrectionWildes.cpp:802
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:1011
Mantid::Algorithms::PolarizationCorrectionWildes::twoInputCorrections
WorkspaceMap twoInputCorrections(const WorkspaceMap &inputs, const EfficiencyMap &efficiencies)
Correct for non-ideal instrument effects.
Definition PolarizationCorrectionWildes.cpp:781
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:896
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:538
Mantid::Algorithms::PolarizationCorrectionWildes::analyzerlessCorrections
WorkspaceMap analyzerlessCorrections(const WorkspaceMap &inputs, const EfficiencyMap &efficiencies)
Correct for non-ideal instrument effects.
Definition PolarizationCorrectionWildes.cpp:710
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:588
Mantid::API::MatrixWorkspace_sptr
std::shared_ptr< MatrixWorkspace > MatrixWorkspace_sptr
shared pointer to the matrix workspace base class
Definition MatrixWorkspace.h:586
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:318
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