Mantid
Loading...
Searching...
No Matches
PolarizationAngleCorrectionMD.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 +
11#include "MantidAPI/Run.h"
23#include <math.h>
24
25using namespace Mantid::Kernel;
26using namespace Mantid::API;
27using namespace Mantid::Geometry;
28using namespace Mantid::DataObjects;
29
30namespace Mantid::MDAlgorithms {
31
32// Register the algorithm into the AlgorithmFactory
33DECLARE_ALGORITHM(PolarizationAngleCorrectionMD)
34
35
36const std::string PolarizationAngleCorrectionMD::name() const { return "PolarizationAngleCorrectionMD"; }
37
40
42const std::string PolarizationAngleCorrectionMD::summary() const {
43 return "Apply polarization angle correction to MDEventWorkspace";
44}
45
47const std::string PolarizationAngleCorrectionMD::category() const { return "MDAlgorithms"; }
48
49//---------------------------------------------------------------------------------------------------------
54
56 std::make_unique<WorkspaceProperty<API::IMDEventWorkspace>>("InputWorkspace", "", Kernel::Direction::Input),
57 "An input MDEventWorkspace. Must be in Q_sample/Q_lab frame. Must have an axis as DeltaE");
58
59 auto anglerange = std::make_shared<BoundedValidator<double>>();
60 anglerange->setLower(-180.);
61 anglerange->setUpper(180.);
62 declareProperty("PolarizationAngle", 0., anglerange,
63 "An in-plane polarization angle​, between -180 and 180 degrees");
64
65 auto precisionrange = std::make_shared<BoundedValidator<double>>();
66 precisionrange->setLower(0.);
67 precisionrange->setUpper(1.);
69 "Precision", 1., precisionrange,
70 "Precision (between 0 and 1). Any event whose absolute value of cosine of 2 of its schaf angle less than this "
71 "precision will be ignored.");
72
74 std::make_unique<WorkspaceProperty<API::IMDEventWorkspace>>("OutputWorkspace", "", Kernel::Direction::Output),
75 "The output MDEventWorkspace with polarization angle correction applied");
76}
77
78//---------------------------------------------------------------------------------------------------------
83 // Get input workspace and other parameters
84 API::IMDEventWorkspace_sptr input_ws = getProperty("InputWorkspace");
85 mPolarizationAngle = getProperty("PolarizationAngle");
86 mPolarizationAngle *= M_PI / 180.; // convert to arcs
87 mPrecision = getProperty("Precision");
88
89 // Process input workspace and create output workspace
90 std::string output_ws_name = getPropertyValue("OutputWorkspace");
91
92 API::IMDEventWorkspace_sptr output_ws(0);
93 if (input_ws->getName() == output_ws_name) {
94 // Calcualte in-place
95 output_ws = input_ws;
96 } else {
97 // Clone input workace to output workspace
98 output_ws = input_ws->clone();
99 }
100
101 // Apply polarization angle correction to MDEvents
103
104 // refresh cache for MDBoxes: set correct Box signal
105 output_ws->refreshCache();
106
107 // Clear masking (box flags) from the output workspace
108 output_ws->clearMDMasking();
109
110 // Set output
111 setProperty("OutputWorkspace", output_ws);
112}
113
114//---------------------------------------------------------------------------------------------------------
125std::map<std::string, std::string> PolarizationAngleCorrectionMD::validateInputs() {
126 std::map<std::string, std::string> output;
127
128 // check input dimension
129 std::string dim_error = checkInputMDDimension();
130 if (dim_error != "") {
131 output["InputWorkspace"] = dim_error;
132 }
133
134 return output;
135}
136
137//---------------------------------------------------------------------------------------------
141template <typename MDE, size_t nd>
144 // Get Box from MDEventWorkspace
145 MDBoxBase<MDE, nd> *box1 = ws->getBox();
146 std::vector<API::IMDNode *> boxes;
147 box1->getBoxes(boxes, 1000, true);
148 auto numBoxes = int(boxes.size());
149
150 // Add the boxes in parallel. They should be spread out enough on each
151 // core to avoid stepping on each other.
152
153 PRAGMA_OMP( parallel for if (!ws->isFileBacked()))
154 for (int i = 0; i < numBoxes; ++i) {
156 auto *box = dynamic_cast<MDBox<MDE, nd> *>(boxes[i]);
157 assert(box);
158 if (!box->getIsMasked()) {
159 // get the MEEvents from box
160 std::vector<MDE> &events = box->getEvents();
161 // Add events, with bounds checking
162 for (auto it = events.begin(); it != events.end(); ++it) {
163 // Modify the event
164
165 // Calculate Gamma
166 double qx(0.), qz(0.);
167 if (!mIsQSample) {
168 // Q-lab: gamma = arctan2(Qx​,Qz​)
169 qx = static_cast<double>(it->getCenter(mQxIndex));
170 qz = static_cast<double>(it->getCenter(mQzIndex));
171 } else {
172 // Q-sample
173 // Qlab = R * QSample
174 std::vector<double> qsample;
175 for (auto d = 0u; d < nd; ++d) {
176 qsample.emplace_back(it->getCenter(d));
177 }
178 std::vector<double> qlab = mRotationMatrixMap[it->getExpInfoIndex()] * qsample;
179 qx = qlab[0];
180 qz = qlab[2];
181 }
182 double gamma(std::atan2(qx, qz)); // unit = arc
183
184 // The Scharpf angle \alphs = \gamma - P_A
185 double alpha = gamma - mPolarizationAngle;
186 // Calculate cosine 2*alpha
187 double cosine2alpha = std::cos(2 * alpha);
188 // If absolute value of consine 2*alpha is larger than Precision
189 float factor(0.);
190 if (fabs(cosine2alpha) > mPrecision) {
191 factor = static_cast<float>(1. / cosine2alpha);
192 }
193
194 // calcalate and set intesity: I *= F
195 auto intensity = it->getSignal() * factor;
196 it->setSignal(intensity);
197
198 // calculate and set error: Err2∗=F^2
199 auto error2 = it->getErrorSquared() * factor * factor;
200 // error2 *= factor * factor;
201 it->setErrorSquared(error2);
202 }
203 }
204 box->releaseEvents();
206 }
208
209 return;
210}
211
212//---------------------------------------------------------------------------------------------------------
221 std::string errormsg("");
222
223 // Check Q-dimension and determine Qx and Qz index
224 API::IMDEventWorkspace_sptr inputws = getProperty("InputWorkspace");
225 size_t numdims = inputws->getNumDims();
226 std::string qxname("Q_lab_x");
227 std::string qzname("Q_lab_z");
228 if (numdims < 4) {
229 errormsg = "Input workspace must have at least 4 dimensions";
230 } else {
231 // Get and check the dimensions: Q3D or Q1D
232 const Mantid::Kernel::SpecialCoordinateSystem coordsys = inputws->getSpecialCoordinateSystem();
234 mIsQSample = false;
236 // q3d
237 mIsQSample = true;
238 // reset Qx and Qz name
239 qxname = "Q_sample_x";
240 qzname = "Q_sample_z";
241 } else {
242 // not supported
243 errormsg = "InputWorkspace is not in Q-Sample or Q-lab frame";
244 }
245
246 // determine Qx and Qz index
247 for (size_t i = 0; i < numdims; ++i) {
248 if (inputws->getDimension(i)->getName() == qxname)
249 mQxIndex = i;
250 else if (inputws->getDimension(i)->getName() == qzname)
251 mQzIndex = i;
252 }
253 // verify and information
254 if (mQxIndex != 0 || mQzIndex != 2)
255 throw std::runtime_error("Qx, Qy and Qz are not in (Qx, Qy, Qz) order");
256 else
257 g_log.information() << "Found " << qxname << " at " << mQxIndex << ", " << qzname << " at " << mQzIndex << "\n";
258
259 // Check DeltaE
260 if (errormsg.size() > 0 && inputws->getDimension(3)->getName() != "DeltaE") {
261 errormsg = "4-th dimension is " + inputws->getDimension(3)->getName() + ". Must be DeltaE";
262 return errormsg;
263 }
264 }
265
266 // Check rotation matrix
267 auto numexpinfo = inputws->getNumExperimentInfo();
268 for (uint16_t i = 0; i < numexpinfo; ++i) {
269 const Kernel::Matrix<double> &rotmatrix = inputws->getExperimentInfo(i)->run().getGoniometerMatrix();
270 mRotationMatrixMap[i] = rotmatrix;
271 }
272
273 return errormsg;
274}
275
276//---------------------------------------------------------------------------------------------------------
283 // Get temperture sample log name
284 std::string Estring("Ei");
285 std::stringstream eiss;
286
287 // the input property could be a valid float; if not must search the experiment info
288 uint16_t numexpinfo = mdws->getNumExperimentInfo();
289
290 for (uint16_t i = 0; i < numexpinfo; ++i) {
291
292 // if user specified is not a valid float
293 ExperimentInfo_const_sptr expinfo = mdws->getExperimentInfo(i);
294 if (expinfo->run().hasProperty(Estring)) {
295 std::string eistr = expinfo->run().getProperty(Estring)->value();
296 try {
297 double ei = boost::lexical_cast<double>(eistr);
298 if (ei <= 0) {
299 // Ei is not greater than 0 and is not allowed
300 eiss << "Experiment Info Ei " << ei << " cannot be zero or less than zero.";
301 }
302 } catch (...) {
303 // unable cast to double
304 eiss << "Experiment Info Ei " << eistr << " cannot be cast to a double number";
305 }
306 } else {
307 // does not have Ei
308 eiss << "Experiment Info " << i << " does not have " << Estring;
309 }
310 }
311
312 // return error string
313 std::string ei_error = eiss.str();
314 return ei_error;
315}
316
317} // namespace Mantid::MDAlgorithms
#define DECLARE_ALGORITHM(classname)
Definition: Algorithm.h:576
#define CALL_MDEVENT_FUNCTION(funcname, workspace)
Macro that makes it possible to call a templated method for a MDEventWorkspace using a IMDEventWorksp...
#define fabs(x)
Definition: Matrix.cpp:22
#define PARALLEL_START_INTERRUPT_REGION
Begins a block to skip processing is the algorithm has been interupted Note the end of the block if n...
Definition: MultiThreaded.h:94
#define PARALLEL_END_INTERRUPT_REGION
Ends a block to skip processing is the algorithm has been interupted Note the start of the block if n...
#define PRAGMA_OMP(expression)
#define PARALLEL_CHECK_INTERRUPT_REGION
Adds a check after a Parallel region to see if it was interupted.
void declareProperty(std::unique_ptr< Kernel::Property > p, const std::string &doc="") override
Add a property to the list of managed properties.
Definition: Algorithm.cpp:1913
std::string getPropertyValue(const std::string &name) const override
Get the value of a property as a string.
Definition: Algorithm.cpp:2026
TypedValue getProperty(const std::string &name) const override
Get the value of a property.
Definition: Algorithm.cpp:2076
Kernel::Logger & g_log
Definition: Algorithm.h:451
virtual void getBoxes(std::vector< IMDNode * > &boxes, size_t maxDepth, bool leafOnly)=0
Fill a vector with all the boxes who are the childred of this one up to a certain depth.
A property class for workspaces.
Templated super-class of a multi-dimensional event "box".
Definition: MDBoxBase.h:50
Templated class for a multi-dimensional event "box".
Definition: MDBox.h:45
std::shared_ptr< MDEventWorkspace< MDE, nd > > sptr
Typedef for a shared pointer of this kind of event workspace.
IPropertyManager * setProperty(const std::string &name, const T &value)
Templated method to set the value of a PropertyWithValue.
void information(const std::string &msg)
Logs at information level.
Definition: Logger.cpp:105
Numerical Matrix class.
Definition: Matrix.h:42
PolarizationAngleCorrection : Perform polarization angle correction to an MDEventWorkspace.
std::string checkEi(const API::IMDEventWorkspace_sptr &mdws)
Get Ei.
std::string checkInputMDDimension()
Check input workspace dimension.
std::map< std::string, std::string > validateInputs() override
Validate inputs.
int version() const override
Algorithm's version for identification.
std::map< uint16_t, Mantid::Kernel::Matrix< double > > mRotationMatrixMap
Map.
void applyPolarizationAngleCorrection(typename Mantid::DataObjects::MDEventWorkspace< MDE, nd >::sptr ws)
Apply polarization angle correction to each MDEvent.
std::shared_ptr< IMDEventWorkspace > IMDEventWorkspace_sptr
Shared pointer to Mantid::API::IMDEventWorkspace.
std::shared_ptr< const ExperimentInfo > ExperimentInfo_const_sptr
Shared pointer to const ExperimentInfo.
SpecialCoordinateSystem
Special coordinate systems for Q3D.
STL namespace.
@ Input
An input workspace.
Definition: Property.h:53
@ Output
An output workspace.
Definition: Property.h:54