Mantid
Loading...
Searching...
No Matches
PolarizationCorrectionFredrikze.cpp
Go to the documentation of this file.
1// Mantid Repository : https://github.com/mantidproject/mantid
2//
3// Copyright © 2018 ISIS Rutherford Appleton Laboratory UKRI,
4// NScD Oak Ridge National Laboratory, European Spallation Source,
5// Institut Laue - Langevin & CSNS, Institute of High Energy Physics, CAS
6// SPDX - License - Identifier: GPL - 3.0 +
8#include "MantidAPI/Axis.h"
10#include "MantidAPI/TextAxis.h"
16
17#include <memory>
18
19#include <algorithm>
20
21using namespace Mantid::API;
22using namespace Mantid::Kernel;
23using namespace Mantid::Geometry;
24
25namespace {
26
27const std::string pNRLabel("PNR");
28
29const std::string pALabel("PA");
30
31const std::string crhoLabel("Rho");
32
33const std::string cppLabel("Pp");
34
35const std::string cAlphaLabel("Alpha");
36
37const std::string cApLabel("Ap");
38
39const std::string efficienciesLabel("Efficiencies");
40
41std::vector<std::string> modes() {
42 std::vector<std::string> modes;
43 modes.emplace_back(pALabel);
44 modes.emplace_back(pNRLabel);
45 return modes;
46}
47
48Instrument_const_sptr fetchInstrument(WorkspaceGroup const *const groupWS) {
49 if (groupWS->size() == 0) {
50 throw std::invalid_argument("Input group workspace has no children.");
51 }
52 Workspace_sptr firstWS = groupWS->getItem(0);
53 MatrixWorkspace_sptr matrixWS = std::dynamic_pointer_cast<MatrixWorkspace>(firstWS);
54 return matrixWS->getInstrument();
55}
56
57void validateInputWorkspace(WorkspaceGroup_sptr &ws) {
59 for (size_t i = 0; i < ws->size(); ++i) {
60
61 Workspace_sptr item = ws->getItem(i);
62
63 if (MatrixWorkspace_sptr ws2d = std::dynamic_pointer_cast<MatrixWorkspace>(item)) {
64
65 // X-units check
66 auto wsUnit = ws2d->getAxis(0)->unit();
67 auto expectedUnit = Units::Wavelength();
68 if (wsUnit->unitID() != expectedUnit.unitID()) {
69 throw std::invalid_argument("Input workspaces must have units of Wavelength");
70 }
71
72 // More detailed checks based on shape.
73 if (lastWS) {
74 if (lastWS->getNumberHistograms() != ws2d->getNumberHistograms()) {
75 throw std::invalid_argument("Not all workspaces in the "
76 "InputWorkspace WorkspaceGroup have the "
77 "same number of spectrum");
78 }
79 if (lastWS->blocksize() != ws2d->blocksize()) {
80 throw std::invalid_argument("Number of bins do not match between all "
81 "workspaces in the InputWorkspace "
82 "WorkspaceGroup");
83 }
84
85 auto &currentX = ws2d->x(0);
86 auto &lastX = lastWS->x(0);
87 auto xMatches = std::equal(lastX.cbegin(), lastX.cend(), currentX.cbegin());
88 if (!xMatches) {
89 throw std::invalid_argument("X-arrays do not match between all "
90 "workspaces in the InputWorkspace "
91 "WorkspaceGroup.");
92 }
93 }
94
95 lastWS = ws2d; // Cache the last workspace so we can use it for comparison
96 // purposes.
97
98 } else {
99 std::stringstream messageBuffer;
100 messageBuffer << "Item with index: " << i << "in the InputWorkspace is not a MatrixWorkspace";
101 throw std::invalid_argument(messageBuffer.str());
102 }
103 }
104}
105
106using VecDouble = std::vector<double>;
107} // namespace
108
109namespace Mantid::Algorithms {
110
111// Register the algorithm into the AlgorithmFactory
112DECLARE_ALGORITHM(PolarizationCorrectionFredrikze)
113
114//----------------------------------------------------------------------------------------------
116const std::string PolarizationCorrectionFredrikze::name() const { return "PolarizationCorrectionFredrikze"; }
117
120
122const std::string PolarizationCorrectionFredrikze::category() const { return "Reflectometry"; }
123
128 return "Makes corrections for polarization efficiencies of the polarizer and "
129 "analyzer in a reflectometry neutron spectrometer.";
130}
131
139 auto multiply = this->createChildAlgorithm("Multiply");
140 auto rhsWS = std::make_shared<DataObjects::WorkspaceSingleValue>(rhs);
141 multiply->initialize();
142 multiply->setProperty("LHSWorkspace", lhsWS);
143 multiply->setProperty("RHSWorkspace", rhsWS);
144 multiply->execute();
145 MatrixWorkspace_sptr outWS = multiply->getProperty("OutputWorkspace");
146 return outWS;
147}
148
156 auto plus = this->createChildAlgorithm("Plus");
157 auto rhsWS = std::make_shared<DataObjects::WorkspaceSingleValue>(rhs);
158 plus->initialize();
159 plus->setProperty("LHSWorkspace", lhsWS);
160 plus->setProperty("RHSWorkspace", rhsWS);
161 plus->execute();
162 MatrixWorkspace_sptr outWS = plus->getProperty("OutputWorkspace");
163 return outWS;
164}
165
166//----------------------------------------------------------------------------------------------
171 std::make_unique<WorkspaceProperty<Mantid::API::WorkspaceGroup>>("InputWorkspace", "", Direction::Input),
172 "An input workspace to process.");
173
174 auto propOptions = modes();
175 declareProperty("PolarizationAnalysis", "PA", std::make_shared<StringListValidator>(propOptions),
176 "What Polarization mode will be used?\n"
177 "PNR: Polarized Neutron Reflectivity mode\n"
178 "PA: Full Polarization Analysis PNR-PA");
179
181 std::make_unique<API::WorkspaceProperty<API::MatrixWorkspace>>(efficienciesLabel, "", Kernel::Direction::Input),
182 "A workspace containing the efficiency factors Pp, Ap, Rho and Alpha "
183 "as histograms");
184
186 std::make_unique<WorkspaceProperty<Mantid::API::WorkspaceGroup>>("OutputWorkspace", "", Direction::Output),
187 "An output workspace.");
188}
189
191
192 size_t itemIndex = 0;
193 MatrixWorkspace_sptr Ipp = std::dynamic_pointer_cast<MatrixWorkspace>(inWS->getItem(itemIndex++));
194 MatrixWorkspace_sptr Ipa = std::dynamic_pointer_cast<MatrixWorkspace>(inWS->getItem(itemIndex++));
195 MatrixWorkspace_sptr Iap = std::dynamic_pointer_cast<MatrixWorkspace>(inWS->getItem(itemIndex++));
196 MatrixWorkspace_sptr Iaa = std::dynamic_pointer_cast<MatrixWorkspace>(inWS->getItem(itemIndex++));
197
198 Ipp->setTitle("Ipp");
199 Iaa->setTitle("Iaa");
200 Ipa->setTitle("Ipa");
201 Iap->setTitle("Iap");
202
203 const auto rho = this->getEfficiencyWorkspace(crhoLabel);
204 const auto pp = this->getEfficiencyWorkspace(cppLabel);
205 const auto alpha = this->getEfficiencyWorkspace(cAlphaLabel);
206 const auto ap = this->getEfficiencyWorkspace(cApLabel);
207
208 const auto A0 = (Iaa * pp * ap) + (Ipa * ap * rho * pp) + (Iap * ap * alpha * pp) + (Ipp * ap * alpha * rho * pp);
209 const auto A1 = Iaa * pp;
210 const auto A2 = Iap * pp;
211 const auto A3 = Iaa * ap;
212 const auto A4 = Ipa * ap;
213 const auto A5 = Ipp * ap * alpha;
214 const auto A6 = Iap * ap * alpha;
215 const auto A7 = Ipp * pp * rho;
216 const auto A8 = Ipa * pp * rho;
217
218 const auto D = pp * ap * (rho + alpha + 1.0 + (rho * alpha));
219
220 const auto nIpp = (A0 - A1 + A2 - A3 + A4 + A5 - A6 + A7 - A8 + Ipp + Iaa - Ipa - Iap) / D;
221 const auto nIaa = (A0 + A1 - A2 + A3 - A4 - A5 + A6 - A7 + A8 + Ipp + Iaa - Ipa - Iap) / D;
222 const auto nIap = (A0 - A1 + A2 + A3 - A4 - A5 + A6 + A7 - A8 - Ipp - Iaa + Ipa + Iap) / D;
223 const auto nIpa = (A0 + A1 - A2 - A3 + A4 + A5 - A6 - A7 + A8 - Ipp - Iaa + Ipa + Iap) / D;
224
225 WorkspaceGroup_sptr dataOut = std::make_shared<WorkspaceGroup>();
226 dataOut->addWorkspace(nIpp);
227 dataOut->addWorkspace(nIpa);
228 dataOut->addWorkspace(nIap);
229 dataOut->addWorkspace(nIaa);
230 size_t totalGroupEntries(dataOut->getNumberOfEntries());
231 for (size_t i = 1; i < totalGroupEntries; i++) {
232 auto alg = this->createChildAlgorithm("ReplaceSpecialValues");
233 alg->setProperty("InputWorkspace", dataOut->getItem(i));
234 alg->setProperty("OutputWorkspace", "dataOut_" + std::to_string(i));
235 alg->setProperty("NaNValue", 0.0);
236 alg->setProperty("NaNError", 0.0);
237 alg->setProperty("InfinityValue", 0.0);
238 alg->setProperty("InfinityError", 0.0);
239 alg->execute();
240 }
241 // Preserve the history of the inside workspaces
242 nIpp->history().addHistory(Ipp->getHistory());
243 nIaa->history().addHistory(Iaa->getHistory());
244 nIpa->history().addHistory(Ipa->getHistory());
245 nIap->history().addHistory(Iap->getHistory());
246
247 return dataOut;
248}
249
251 size_t itemIndex = 0;
252 MatrixWorkspace_sptr Ip = std::dynamic_pointer_cast<MatrixWorkspace>(inWS->getItem(itemIndex++));
253 MatrixWorkspace_sptr Ia = std::dynamic_pointer_cast<MatrixWorkspace>(inWS->getItem(itemIndex++));
254
255 const auto rho = this->getEfficiencyWorkspace(crhoLabel);
256 const auto pp = this->getEfficiencyWorkspace(cppLabel);
257
258 const auto D = pp * (rho + 1);
259
260 const auto nIp = (Ip * (rho * pp + 1.0) + Ia * (pp - 1.0)) / D;
261 const auto nIa = (Ip * (rho * pp - 1.0) + Ia * (pp + 1.0)) / D;
262
263 // Preserve the history of the inside workspaces
264 nIp->history().addHistory(Ip->getHistory());
265 nIa->history().addHistory(Ia->getHistory());
266
267 WorkspaceGroup_sptr dataOut = std::make_shared<WorkspaceGroup>();
268 dataOut->addWorkspace(nIp);
269 dataOut->addWorkspace(nIa);
270
271 return dataOut;
272}
273
278std::shared_ptr<Mantid::API::MatrixWorkspace>
280 MatrixWorkspace_sptr efficiencies = getProperty(efficienciesLabel);
281 auto const &axis = dynamic_cast<TextAxis &>(*efficiencies->getAxis(1));
282 size_t index = axis.length();
283 for (size_t i = 0; i < axis.length(); ++i) {
284 if (axis.label(i) == label) {
285 index = i;
286 break;
287 }
288 }
289
290 if (index == axis.length()) {
291 // Check if we need to fetch polarization parameters from the instrument's
292 // parameters
293 static std::map<std::string, std::string> loadableProperties{
294 {crhoLabel, "crho"}, {cppLabel, "cPp"}, {cApLabel, "cAp"}, {cAlphaLabel, "calpha"}};
295 WorkspaceGroup_sptr inWS = getProperty("InputWorkspace");
296 Instrument_const_sptr instrument = fetchInstrument(inWS.get());
297 auto vals = instrument->getStringParameter(loadableProperties[label]);
298 if (vals.empty()) {
299 throw std::invalid_argument("Efficiencey property not found: " + label);
300 }
301 auto extract = createChildAlgorithm("CreatePolarizationEfficiencies");
302 extract->initialize();
303 extract->setProperty("InputWorkspace", efficiencies);
304 extract->setProperty(label, vals.front());
305 extract->execute();
306 MatrixWorkspace_sptr outWS = extract->getProperty("OutputWorkspace");
307 return outWS;
308 } else {
309 auto extract = createChildAlgorithm("ExtractSingleSpectrum");
310 extract->initialize();
311 extract->setProperty("InputWorkspace", efficiencies);
312 extract->setProperty("WorkspaceIndex", static_cast<int>(index));
313 extract->execute();
314 MatrixWorkspace_sptr outWS = extract->getProperty("OutputWorkspace");
315 return outWS;
316 }
317}
318
319//----------------------------------------------------------------------------------------------
323 WorkspaceGroup_sptr inWS = getProperty("InputWorkspace");
324 const std::string analysisMode = getProperty("PolarizationAnalysis");
325 const size_t nWorkspaces = inWS->size();
326
327 validateInputWorkspace(inWS);
328
330 if (analysisMode == pALabel) {
331 if (nWorkspaces != 4) {
332 throw std::invalid_argument("For PA analysis, input group must have 4 periods.");
333 }
334 g_log.notice("PA polarization correction");
335 outWS = execPA(inWS);
336 } else if (analysisMode == pNRLabel) {
337 if (nWorkspaces != 2) {
338 throw std::invalid_argument("For PNR analysis, input group must have 2 periods.");
339 }
340 outWS = execPNR(inWS);
341 g_log.notice("PNR polarization correction");
342 }
343 this->setProperty("OutputWorkspace", outWS);
344}
345
346} // namespace Mantid::Algorithms
#define DECLARE_ALGORITHM(classname)
Definition: Algorithm.h:576
const std::vector< double > & rhs
std::map< DeltaEMode::Type, std::string > index
Definition: DeltaEMode.cpp:19
std::vector< double > VecDouble
void declareProperty(std::unique_ptr< Kernel::Property > p, const std::string &doc="") override
Add a property to the list of managed properties.
Definition: Algorithm.cpp:1913
TypedValue getProperty(const std::string &name) const override
Get the value of a property.
Definition: Algorithm.cpp:2076
virtual std::shared_ptr< Algorithm > createChildAlgorithm(const std::string &name, const double startProgress=-1., const double endProgress=-1., const bool enableLogging=true, const int &version=-1)
Create a Child Algorithm.
Definition: Algorithm.cpp:842
Kernel::Logger & g_log
Definition: Algorithm.h:451
Class to represent a text axis of a workspace.
Definition: TextAxis.h:36
Class to hold a set of workspaces.
Workspace_sptr getItem(const size_t index) const
Return the ith workspace.
size_t size() const
Return the size of the group, so it is more like a container.
A property class for workspaces.
PolarizationCorrectionFredrikze : Algorithm to perform polarisation corrections on multi-period group...
std::shared_ptr< Mantid::API::WorkspaceGroup > execPA(const std::shared_ptr< Mantid::API::WorkspaceGroup > &inWS)
std::shared_ptr< Mantid::API::MatrixWorkspace > add(std::shared_ptr< Mantid::API::MatrixWorkspace > &lhsWS, const double &rhs)
Add a constant value to a workspace.
void init() override
Initialize the algorithm's properties.
std::shared_ptr< Mantid::API::WorkspaceGroup > execPNR(const std::shared_ptr< Mantid::API::WorkspaceGroup > &inWS)
std::shared_ptr< Mantid::API::MatrixWorkspace > getEfficiencyWorkspace(const std::string &label)
Extract a spectrum from the Efficiencies workspace as a 1D workspace.
int version() const override
Algorithm's version for identification.
std::shared_ptr< Mantid::API::MatrixWorkspace > multiply(std::shared_ptr< Mantid::API::MatrixWorkspace > &lhsWS, const double &rhs)
Multiply a workspace by a constant value.
const std::string category() const override
Algorithm's category for identification.
IPropertyManager * setProperty(const std::string &name, const T &value)
Templated method to set the value of a PropertyWithValue.
void notice(const std::string &msg)
Logs at notice level.
Definition: Logger.cpp:95
Wavelength in Angstrom.
Definition: Unit.h:302
std::shared_ptr< WorkspaceGroup > WorkspaceGroup_sptr
shared pointer to Mantid::API::WorkspaceGroup
std::shared_ptr< Workspace > Workspace_sptr
shared pointer to Mantid::API::Workspace
Definition: Workspace_fwd.h:20
std::shared_ptr< MatrixWorkspace > MatrixWorkspace_sptr
shared pointer to the matrix workspace base class
std::shared_ptr< const Instrument > Instrument_const_sptr
Shared pointer to an const instrument object.
STL namespace.
std::string to_string(const wide_integer< Bits, Signed > &n)
@ Input
An input workspace.
Definition: Property.h:53
@ Output
An output workspace.
Definition: Property.h:54