Mantid
Loading...
Searching...
No Matches
SofQCommon.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/Run.h"
14#include "MantidKernel/Unit.h"
16
17namespace Mantid::Algorithms {
26 // Retrieve the emode & efixed properties
27 const std::string emode = hostAlgorithm->getProperty("EMode");
28 // Convert back to an integer representation
29 m_emode = 0;
30 if (emode == "Direct")
31 m_emode = 1;
32 else if (emode == "Indirect")
33 m_emode = 2;
34 m_efixed = hostAlgorithm->getProperty("EFixed");
35
36 // Check whether they should have supplied an EFixed value
37 if (m_emode == 1) // Direct
38 {
39 // If GetEi was run then it will have been stored in the workspace, if not
40 // the user will need to enter one
41 if (m_efixed == 0.0) {
42 if (workspace.run().hasProperty("Ei")) {
43 Kernel::Property *p = workspace.run().getProperty("Ei");
44 auto *eiProp = dynamic_cast<Kernel::PropertyWithValue<double> *>(p);
45 if (!eiProp)
46 throw std::runtime_error("Input workspace contains Ei but its "
47 "property type is not a double.");
48 m_efixed = (*eiProp)();
49 } else {
50 throw std::invalid_argument("Input workspace does not contain an "
51 "EFixed value. Please provide one or run "
52 "GetEi.");
53 }
54 } else {
55 m_efixedGiven = true;
56 }
57 } else {
58 if (m_efixed != 0.0) {
59 m_efixedGiven = true;
60 }
61 }
62}
63
73 double efixed(0.0);
74 if (m_emode == 1) // Direct
75 {
77 } else // Indirect
78 {
79 if (m_efixedGiven)
80 efixed = m_efixed; // user provided a value
81 else {
82 const std::vector<double> param = det.getNumberParameter("EFixed");
83 if (param.empty())
84 throw std::runtime_error("Cannot find EFixed parameter for component \"" + det.getName() +
85 "\". This is required in indirect mode. Please check the IDF "
86 "contains these values.");
87 efixed = param[0];
88 }
89 }
90 return efixed;
91}
92
101double SofQCommon::q(const double deltaE, const double twoTheta, const Geometry::IDetector *det) const {
102 if (m_emode == 1) {
103 return directQ(deltaE, twoTheta);
104 }
105 return indirectQ(deltaE, twoTheta, det);
106}
107
115std::pair<double, double> SofQCommon::qBinHints(const API::MatrixWorkspace &ws, const double minE,
116 const double maxE) const {
117 if (m_emode == 1) {
118 return qBinHintsDirect(ws, minE, maxE);
119 }
120 return qBinHintsIndirect(ws, minE, maxE);
121}
122
129double SofQCommon::directQ(const double deltaE, const double twoTheta) const {
131 const double ki = std::sqrt(m_efixed / E_mev_toNeutronWavenumberSq);
132 const double kf = std::sqrt((m_efixed - deltaE) / E_mev_toNeutronWavenumberSq);
133 return std::sqrt(ki * ki + kf * kf - 2. * ki * kf * std::cos(twoTheta));
134}
135
143double SofQCommon::indirectQ(const double deltaE, const double twoTheta, const Geometry::IDetector *det) const {
145 if (!det) {
146 throw std::runtime_error("indirectQ: det is nullptr.");
147 }
148 const auto efixed = getEFixed(*det);
149 const double ki = std::sqrt((efixed + deltaE) / E_mev_toNeutronWavenumberSq);
150 const double kf = std::sqrt(efixed / E_mev_toNeutronWavenumberSq);
151 return std::sqrt(ki * ki + kf * kf - 2. * ki * kf * std::cos(twoTheta));
152}
153
162std::pair<double, double> SofQCommon::qBinHintsDirect(const API::MatrixWorkspace &ws, const double minE,
163 const double maxE) const {
164 using namespace Mantid::PhysicalConstants;
165 if (maxE > m_efixed) {
166 throw std::invalid_argument("Cannot compute Q binning range: maximum "
167 "energy transfer is greater than the incident "
168 "energy.");
169 }
170 auto minTwoTheta = std::numeric_limits<double>::max();
171 auto maxTwoTheta = std::numeric_limits<double>::lowest();
172 const auto &spectrumInfo = ws.spectrumInfo();
173 for (size_t i = 0; i < spectrumInfo.size(); ++i) {
174 if (spectrumInfo.isMasked(i) || spectrumInfo.isMonitor(i)) {
175 continue;
176 }
177 const auto twoTheta = spectrumInfo.twoTheta(i);
178 minTwoTheta = std::min(minTwoTheta, twoTheta);
179 maxTwoTheta = std::max(maxTwoTheta, twoTheta);
180 }
181 if (minTwoTheta == std::numeric_limits<double>::max()) {
182 throw std::runtime_error("Could not determine Q binning: workspace does "
183 "not contain usable spectra.");
184 }
185 std::array<double, 4> q;
186 q[0] = directQ(minE, minTwoTheta);
187 q[1] = directQ(minE, maxTwoTheta);
188 q[2] = directQ(maxE, minTwoTheta);
189 q[3] = directQ(maxE, maxTwoTheta);
190 const auto minmaxQ = std::minmax_element(q.cbegin(), q.cend());
191 return std::make_pair(*minmaxQ.first, *minmaxQ.second);
192}
193
204std::pair<double, double> SofQCommon::qBinHintsIndirect(const API::MatrixWorkspace &ws, const double minE,
205 const double maxE) const {
206 using namespace Mantid::PhysicalConstants;
207 auto minQ = std::numeric_limits<double>::max();
208 auto maxQ = std::numeric_limits<double>::lowest();
209 const auto &detectorInfo = ws.detectorInfo();
210 for (size_t i = 0; i < detectorInfo.size(); ++i) {
211 if (detectorInfo.isMasked(i) || detectorInfo.isMonitor(i)) {
212 continue;
213 }
214 const auto twoTheta = detectorInfo.twoTheta(i);
215 const auto &det = detectorInfo.detector(i);
216 const auto Q1 = indirectQ(minE, twoTheta, &det);
217 const auto Q2 = indirectQ(maxE, twoTheta, &det);
218 if (!std::isfinite(Q1) || !std::isfinite(Q2)) {
219 throw std::invalid_argument("Cannot compute Q binning range: non-finite "
220 "Q found for detector ID " +
221 std::to_string(detectorInfo.detectorIDs()[i]));
222 }
223 const auto minmaxQ = std::minmax(Q1, Q2);
224 minQ = std::min(minQ, minmaxQ.first);
225 maxQ = std::max(maxQ, minmaxQ.second);
226 }
227 if (minQ == std::numeric_limits<double>::max()) {
228 throw std::runtime_error("Could not determine Q binning: workspace does "
229 "not contain usable spectra.");
230 }
231 return std::make_pair(minQ, maxQ);
232}
233} // namespace Mantid::Algorithms
IPeaksWorkspace_sptr workspace
Definition: IndexPeaks.cpp:114
Base class from which all concrete algorithm classes should be derived.
Definition: Algorithm.h:85
TypedValue getProperty(const std::string &name) const override
Get the value of a property.
Definition: Algorithm.cpp:2076
const SpectrumInfo & spectrumInfo() const
Return a reference to the SpectrumInfo object.
const Geometry::DetectorInfo & detectorInfo() const
Return a const reference to the DetectorInfo object.
Base MatrixWorkspace Abstract Class.
virtual std::vector< double > getNumberParameter(const std::string &pname, bool recursive=true) const =0
Get a parameter defined as a double.
virtual std::string getName() const =0
Get the IComponent name.
Interface class for detector objects.
Definition: IDetector.h:43
The concrete, templated class for properties.
Base class for properties.
Definition: Property.h:94
A namespace containing physical constants that are required by algorithms and unit routines.
Definition: Atom.h:14
static constexpr double E_mev_toNeutronWavenumberSq
Transformation coefficient to transform neutron energy into neutron wavevector: K-neutron[m^-10] = sq...
std::string to_string(const wide_integer< Bits, Signed > &n)
std::pair< double, double > qBinHintsDirect(const API::MatrixWorkspace &ws, const double minE, const double maxE) const
Return a pair of (minimum Q, maximum Q) for given direct geometry workspace.
Definition: SofQCommon.cpp:162
double indirectQ(const double deltaE, const double twoTheta, const Geometry::IDetector *det) const
Calculate the Q value for an indirect instrument.
Definition: SofQCommon.cpp:143
std::pair< double, double > qBinHints(const API::MatrixWorkspace &ws, const double minE, const double maxE) const
Estimate minimum and maximum momentum transfer.
Definition: SofQCommon.cpp:115
double directQ(const double deltaE, const double twoTheta) const
Calculate the Q value for a direct instrument.
Definition: SofQCommon.cpp:129
bool m_efixedGiven
EFixed has been provided.
Definition: SofQCommon.h:23
double getEFixed(const Geometry::IDetector &det) const
Get the efixed value for the given detector.
Definition: SofQCommon.cpp:72
std::pair< double, double > qBinHintsIndirect(const API::MatrixWorkspace &ws, const double minE, const double maxE) const
Return a pair of (minimum Q, maximum Q) for given indirect geometry workspace.
Definition: SofQCommon.cpp:204
double q(const double deltaE, const double twoTheta, const Geometry::IDetector *det) const
Calculate the Q value.
Definition: SofQCommon.cpp:101
void initCachedValues(const API::MatrixWorkspace &workspace, API::Algorithm *const hostAlgorithm)
The procedure analyses emode and efixed properties provided to the algorithm and identify the energy ...
Definition: SofQCommon.cpp:25