Mantid
Loading...
Searching...
No Matches
ParallaxCorrection.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 +
10#include "MantidAPI/Progress.h"
15#include "MantidHistogramData/HistogramMath.h"
19
20#include <math.h>
21#include <muParser.h>
22
23namespace {
24static const std::string PARALLAX_PARAMETER = "parallax";
25static const std::string DIRECTION_PARAMETER = "direction";
26
33std::string validateFormula(const std::string &parallax, const std::string &direction) {
34 if (direction != "x" && direction != "y") {
35 return "Direction must be x or y";
36 }
37 mu::Parser muParser;
38 double t = 0.;
39 muParser.DefineVar("t", &t);
40 muParser.SetExpr(parallax);
41 try {
42 muParser.Eval();
43 } catch (mu::Parser::exception_type &e) {
44 return e.GetMsg();
45 }
46 return "";
47}
48} // namespace
49
50namespace Mantid::Algorithms {
51
52// Register the algorithm into the AlgorithmFactory
53DECLARE_ALGORITHM(ParallaxCorrection)
54
55//----------------------------------------------------------------------------------------------
56
57
58const std::string ParallaxCorrection::name() const { return "ParallaxCorrection"; }
59
61int ParallaxCorrection::version() const { return 1; }
62
64const std::string ParallaxCorrection::category() const { return "SANS"; }
65
67const std::string ParallaxCorrection::summary() const {
68 return "Performs parallax correction for tube based SANS instruments.";
69}
70
72std::map<std::string, std::string> ParallaxCorrection::validateInputs() {
73 std::map<std::string, std::string> results;
74 const std::vector<double> angleOffsets = getProperty("AngleOffsets");
75 const std::vector<std::string> componentNames = getProperty("ComponentNames");
76 if (angleOffsets.size() != componentNames.size() && angleOffsets.size() != 1) {
77 results["AngleOffsets"] = "Angle offsets should have one value or as many "
78 "as there are components";
79 }
80 return results;
81}
82
83//----------------------------------------------------------------------------------------------
87 auto validator = std::make_shared<Kernel::CompositeValidator>();
88 validator->add(std::make_unique<API::InstrumentValidator>());
89 auto lengthValidator = std::make_shared<Kernel::ArrayLengthValidator<std::string>>();
90 lengthValidator->setLengthMin(1);
91 declareProperty(std::make_unique<API::WorkspaceProperty<API::MatrixWorkspace>>("InputWorkspace", "",
92 Kernel::Direction::Input, validator),
93 "An input workspace.");
94 declareProperty(std::make_unique<Kernel::ArrayProperty<std::string>>("ComponentNames", lengthValidator),
95 "List of instrument components to perform the corrections for.");
97 std::make_unique<API::WorkspaceProperty<API::MatrixWorkspace>>("OutputWorkspace", "", Kernel::Direction::Output),
98 "An output workspace.");
99
100 declareProperty(std::make_unique<Kernel::ArrayProperty<double>>("AngleOffsets", std::vector<double>{0.0}),
101 "The values of offset angles [degrees] to be subtracted from "
102 "the scattering angle (per component).");
103}
104
105//----------------------------------------------------------------------------------------------
114void ParallaxCorrection::performCorrection(const API::MatrixWorkspace_sptr &outWS, const std::vector<size_t> &indices,
115 const std::string &parallax, const std::string &direction,
116 const double angleOffset) {
117 double t;
118 mu::Parser muParser;
119 muParser.DefineVar("t", &t);
120 muParser.SetExpr(parallax);
121 const auto &detectorInfo = outWS->detectorInfo();
122 // note that this is intenionally serial
123 for (const auto wsIndex : indices) {
124 const Kernel::V3D pos = detectorInfo.position(wsIndex);
125 if (direction == "y") {
126 t = std::fabs(std::atan2(pos.X(), pos.Z()));
127 } else {
128 t = std::fabs(std::atan2(pos.Y(), pos.Z()));
129 }
130 t -= angleOffset * M_PI / 180.;
131 const double correction = muParser.Eval();
132 if (correction > 0.) {
133 auto &spectrum = outWS->mutableY(wsIndex);
134 auto &errors = outWS->mutableE(wsIndex);
135 spectrum /= correction;
136 errors /= correction;
137 } else {
138 g_log.warning() << "Correction is <=0 for workspace index " << wsIndex << ". Skipping the correction.\n";
139 }
140 }
141}
142
143//----------------------------------------------------------------------------------------------
147 API::MatrixWorkspace_const_sptr inputWorkspace = getProperty("InputWorkspace");
148 API::MatrixWorkspace_sptr outputWorkspace = getProperty("OutputWorkspace");
149 if (inputWorkspace != outputWorkspace) {
150 outputWorkspace = inputWorkspace->clone();
151 }
152 const std::vector<double> angleOffsets = getProperty("AngleOffsets");
153 const std::vector<std::string> componentNames = getProperty("ComponentNames");
154 const auto &instrument = inputWorkspace->getInstrument();
155 const auto &detectorInfo = outputWorkspace->detectorInfo();
156 const auto &allDetIDs = detectorInfo.detectorIDs();
157 const auto &componentInfo = outputWorkspace->componentInfo();
158 auto progress = std::make_unique<API::Progress>(this, 0., 1., componentNames.size());
159 size_t componentIndex = 0;
160 for (const auto &componentName : componentNames) {
161 const double angleOffset = angleOffsets[(angleOffsets.size() == 1) ? 0 : componentIndex];
162 ++componentIndex;
163 progress->report("Performing parallax correction for component " + componentName);
164 const auto component = instrument->getComponentByName(componentName);
165 if (!component) {
166 g_log.error() << "No component defined with name " << componentName << "\n";
167 continue;
168 }
169 if (!component->hasParameter(PARALLAX_PARAMETER) || !component->hasParameter(DIRECTION_PARAMETER)) {
170 g_log.error() << "No parallax correction defined in IPF for component " << componentName << "\n";
171 continue;
172 }
173 const std::string parallax = component->getStringParameter(PARALLAX_PARAMETER)[0];
174 const std::string direction = component->getStringParameter(DIRECTION_PARAMETER)[0];
175 const auto valid = validateFormula(parallax, direction);
176 if (!valid.empty()) {
177 g_log.error() << "Unable to parse the parallax formula and direction "
178 "for component "
179 << componentName << ". Reason: " << valid << "\n";
180 continue;
181 }
182 const auto &detectorIndices = componentInfo.detectorsInSubtree(componentInfo.indexOfAny(componentName));
183 if (detectorIndices.empty()) {
184 g_log.error() << "No detectors found in component " << componentName << "\n";
185 continue;
186 }
187 std::vector<detid_t> detIDs;
188 detIDs.reserve(detectorIndices.size());
189 std::transform(detectorIndices.cbegin(), detectorIndices.cend(), std::back_inserter(detIDs),
190 [&allDetIDs](size_t i) { return allDetIDs[i]; });
191 const auto indices = outputWorkspace->getIndicesFromDetectorIDs(detIDs);
192 performCorrection(outputWorkspace, indices, parallax, direction, angleOffset);
193 }
194 setProperty("OutputWorkspace", outputWorkspace);
195}
196
197} // namespace Mantid::Algorithms
#define DECLARE_ALGORITHM(classname)
Definition: Algorithm.h:576
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
Kernel::Logger & g_log
Definition: Algorithm.h:451
void progress(double p, const std::string &msg="", double estimatedTime=0.0, int progressPrecision=0)
Sends ProgressNotification.
Definition: Algorithm.cpp:231
A property class for workspaces.
ParallaxCorrection : Performs geometrical correction for parallax effect in tube based SANS instrumen...
const std::string summary() const override
Algorithm's summary for use in the GUI and help.
const std::string category() const override
Algorithm's category for identification.
std::map< std::string, std::string > validateInputs() override
Validate inputs.
int version() const override
Algorithm's version for identification.
void performCorrection(const API::MatrixWorkspace_sptr &, const std::vector< size_t > &, const std::string &, const std::string &, const double)
ParallaxCorrection::performCorrection for the given bank.
void exec() override
Execute the algorithm.
void init() override
Initialize the algorithm's properties.
Support for a property that holds an array of values.
Definition: ArrayProperty.h:28
IPropertyManager * setProperty(const std::string &name, const T &value)
Templated method to set the value of a PropertyWithValue.
void error(const std::string &msg)
Logs at error level.
Definition: Logger.cpp:77
void warning(const std::string &msg)
Logs at warning level.
Definition: Logger.cpp:86
Class for 3D vectors.
Definition: V3D.h:34
constexpr double X() const noexcept
Get x.
Definition: V3D.h:232
constexpr double Y() const noexcept
Get y.
Definition: V3D.h:233
constexpr double Z() const noexcept
Get z.
Definition: V3D.h:234
std::shared_ptr< const MatrixWorkspace > MatrixWorkspace_const_sptr
shared pointer to the matrix workspace base class (const version)
std::shared_ptr< MatrixWorkspace > MatrixWorkspace_sptr
shared pointer to the matrix workspace base class
STL namespace.
@ Input
An input workspace.
Definition: Property.h:53
@ Output
An output workspace.
Definition: Property.h:54