Mantid
Loading...
Searching...
No Matches
ConvertToYSpace.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
9#include "MantidAPI/Axis.h"
21#include "MantidKernel/Unit.h"
22
24
25// Register the algorithm into the AlgorithmFactory
26DECLARE_ALGORITHM(ConvertToYSpace)
27
28using namespace API;
29using namespace Kernel;
30
31namespace {
33const double MASS_TO_MEV = 0.5 * PhysicalConstants::NeutronMass / PhysicalConstants::meV;
34} // namespace
35
38ConvertToYSpace::ConvertToYSpace()
39 : Algorithm(), m_inputWS(), m_mass(0.0), m_l1(0.0), m_samplePos(), m_outputWS(), m_qOutputWS() {}
40
42const std::string ConvertToYSpace::name() const { return "ConvertToYSpace"; }
43
45int ConvertToYSpace::version() const { return 1; }
46
48const std::string ConvertToYSpace::category() const { return "Transforms\\Units"; }
49
55DetectorParams ConvertToYSpace::getDetectorParameters(const API::MatrixWorkspace_const_sptr &ws, const size_t index) {
56 auto inst = ws->getInstrument();
57 auto sample = inst->getSample();
58 auto source = inst->getSource();
59 if (!sample || !source) {
60 throw std::invalid_argument("ConvertToYSpace - Workspace has no source/sample.");
61 }
62
63 const auto &spectrumInfo = ws->spectrumInfo();
64 if (!spectrumInfo.hasDetectors(index))
65 throw std::invalid_argument("ConvertToYSpace - Workspace has no detector "
66 "attached to histogram at index " +
68
69 DetectorParams detpar;
70 const auto &pmap = ws->constInstrumentParameters();
71 const auto &det = spectrumInfo.detector(index);
72 detpar.l1 = spectrumInfo.l1();
73 detpar.l2 = spectrumInfo.l2(index);
74 detpar.pos = spectrumInfo.position(index);
75 detpar.theta = spectrumInfo.twoTheta(index);
76 detpar.t0 = getComponentParameter(det, pmap, "t0") * 1e-6; // Convert to seconds
77 detpar.efixed = getComponentParameter(det, pmap, "efixed");
78 return detpar;
79}
80
90double ConvertToYSpace::getComponentParameter(const Geometry::IComponent &comp, const Geometry::ParameterMap &pmap,
91 const std::string &name) {
92 double result(0.0);
93 if (const auto &group = dynamic_cast<const Geometry::DetectorGroup *>(&comp)) {
94 double avg(0.0);
95 for (const auto &det : group->getDetectors()) {
96 auto param = pmap.getRecursive(det->getComponentID(), name);
97 if (param)
98 avg += param->value<double>();
99 else
100 throw std::invalid_argument("ComptonProfile - Unable to find "
101 "DetectorGroup component parameter \"" +
102 name + "\".");
103 }
104 result = avg / static_cast<double>(group->nDets());
105 } else {
106 auto param = pmap.getRecursive(comp.getComponentID(), name);
107 if (param) {
108 result = param->value<double>();
109 } else {
110 throw std::invalid_argument("ComptonProfile - Unable to find component parameter \"" + name + "\".");
111 }
112 }
113 return result;
114}
115
116//----------------------------------------------------------------------------------------------
117
130void ConvertToYSpace::calculateY(double &yspace, double &qspace, double &ei, const double mass, const double tsec,
131 const double k1, const double v1, const DetectorParams &detpar) {
132 const double v0 = detpar.l1 / (tsec - detpar.t0 - (detpar.l2 / v1));
133 ei = MASS_TO_MEV * v0 * v0;
134 const double w = ei - detpar.efixed;
135 const double k0 = std::sqrt(ei / PhysicalConstants::E_mev_toNeutronWavenumberSq);
136 qspace = std::sqrt(k0 * k0 + k1 * k1 - 2.0 * k0 * k1 * std::cos(detpar.theta));
137 const double wreduced = PhysicalConstants::E_mev_toNeutronWavenumberSq * qspace * qspace / mass;
138 yspace = 0.2393 * (mass / qspace) * (w - wreduced);
139}
140
141//----------------------------------------------------------------------------------------------
142
145void ConvertToYSpace::init() {
146 auto wsValidator = std::make_shared<CompositeValidator>();
147 wsValidator->add<HistogramValidator>(false); // point data
148 wsValidator->add<InstrumentValidator>();
149 wsValidator->add<WorkspaceUnitValidator>("TOF");
150 declareProperty(std::make_unique<WorkspaceProperty<>>("InputWorkspace", "", Direction::Input, wsValidator),
151 "The input workspace in Time of Flight");
152
153 auto mustBePositive = std::make_shared<BoundedValidator<double>>();
154 mustBePositive->setLower(0.0);
155 mustBePositive->setLowerExclusive(true); // strictly greater than 0.0
156 declareProperty("Mass", -1.0, mustBePositive, "The mass defining the recoil peak in AMU");
157
158 declareProperty(std::make_unique<WorkspaceProperty<>>("OutputWorkspace", "", Direction::Output),
159 "The output workspace in y-Space");
160
162 "The output workspace in q-Space");
163}
164
165//----------------------------------------------------------------------------------------------
166
169void ConvertToYSpace::exec() {
172
173 const auto nhist = static_cast<int64_t>(m_inputWS->getNumberHistograms());
174 const int64_t nreports = nhist;
175 auto progress = std::make_shared<Progress>(this, 0.0, 1.0, nreports);
176
177 auto &spectrumInfo = m_outputWS->mutableSpectrumInfo();
178 SpectrumInfo *qSpectrumInfo{nullptr};
179 if (m_qOutputWS)
180 qSpectrumInfo = &m_qOutputWS->mutableSpectrumInfo();
181
183 for (int64_t i = 0; i < nhist; ++i) {
185
186 if (!convert(i)) {
187 g_log.warning("No detector defined for index=" + std::to_string(i) + ". Zeroing spectrum.");
188 m_outputWS->getSpectrum(i).clearData();
189 PARALLEL_CRITICAL(setMasked) {
190 spectrumInfo.setMasked(i, true);
191 if (m_qOutputWS) {
192 m_qOutputWS->getSpectrum(i).clearData();
193 qSpectrumInfo->setMasked(i, true);
194 }
195 }
196 }
197 progress->report();
199 }
201
202 setProperty("OutputWorkspace", m_outputWS);
203
204 if (m_qOutputWS)
205 setProperty("QWorkspace", m_qOutputWS);
206}
207
214bool ConvertToYSpace::convert(const size_t index) {
215 try {
217 const double v1 = std::sqrt(detPar.efixed / MASS_TO_MEV);
218 const double k1 = std::sqrt(detPar.efixed / PhysicalConstants::E_mev_toNeutronWavenumberSq);
219
220 auto &outX = m_outputWS->mutableX(index);
221 auto &outY = m_outputWS->mutableY(index);
222 auto &outE = m_outputWS->mutableE(index);
223 const auto &inX = m_inputWS->x(index);
224 const auto &inY = m_inputWS->y(index);
225 const auto &inE = m_inputWS->e(index);
226
227 // The t->y mapping flips the order of the axis so we need to reverse it to
228 // have a monotonically increasing axis
229 const size_t npts = inY.size();
230 for (size_t j = 0; j < npts; ++j) {
231 double ys(0.0), qs(0.0), ei(0.0);
232 calculateY(ys, qs, ei, m_mass, inX[j] * 1e-06, k1, v1, detPar);
233 const size_t outIndex = (npts - j - 1);
234 outX[outIndex] = ys;
235 const double prefactor = qs / pow(ei, 0.1);
236 outY[outIndex] = prefactor * inY[j];
237 outE[outIndex] = prefactor * inE[j];
238
239 if (m_qOutputWS) {
240 m_qOutputWS->mutableX(index)[outIndex] = ys;
241 m_qOutputWS->mutableY(index)[outIndex] = qs;
242 }
243 }
244 return true;
245 } catch (Exception::NotFoundError &) {
246 return false;
247 }
248}
249
253void ConvertToYSpace::retrieveInputs() {
254 m_inputWS = getProperty("InputWorkspace");
255 m_mass = getProperty("Mass");
257}
258
262void ConvertToYSpace::createOutputWorkspace() {
263 // y-Space output workspace
265
266 auto xLabel = std::make_shared<Units::Label>("Momentum", "A^-1");
267 m_outputWS->getAxis(0)->unit() = xLabel;
268 m_outputWS->setYUnit("");
269 m_outputWS->setYUnitLabel("");
270
271 // q-Space output workspace
272 if (!getPropertyValue("QWorkspace").empty()) {
274
275 m_qOutputWS->getAxis(0)->unit() = xLabel;
276 m_qOutputWS->setYUnit("");
277 m_qOutputWS->setYUnitLabel("");
278 }
279}
280
281void ConvertToYSpace::cacheInstrumentGeometry() {
282 auto inst = m_inputWS->getInstrument();
283 auto source = inst->getSource();
284 auto sample = inst->getSample();
285 m_l1 = sample->getDistance(*source);
286 m_samplePos = sample->getPos();
287}
288
289} // namespace Mantid::CurveFitting::Algorithms
#define DECLARE_ALGORITHM(classname)
Definition: Algorithm.h:576
std::map< DeltaEMode::Type, std::string > index
Definition: DeltaEMode.cpp:19
#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_CRITICAL(name)
#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 PARALLEL_FOR_IF(condition)
Empty definitions - to enable set your complier to enable openMP.
#define PARALLEL_CHECK_INTERRUPT_REGION
Adds a check after a Parallel region to see if it was interupted.
Base class from which all concrete algorithm classes should be derived.
Definition: Algorithm.h:85
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
void progress(double p, const std::string &msg="", double estimatedTime=0.0, int progressPrecision=0)
Sends ProgressNotification.
Definition: Algorithm.cpp:231
A validator which checks that a workspace contains histogram data (the default) or point data as requ...
A validator which checks that a workspace has a valid instrument.
API::SpectrumInfo is an intermediate step towards a SpectrumInfo that is part of Instrument-2....
Definition: SpectrumInfo.h:53
A property class for workspaces.
A validator which checks that the unit of the workspace referred to by a WorkspaceProperty is the exp...
API::MatrixWorkspace_sptr m_outputWS
Output workspace.
void retrieveInputs()
Check and store appropriate input data.
const std::string name() const override
Algorithm's name for identification.
void cacheInstrumentGeometry()
Compute & store the parameters that are fixed during the correction.
API::MatrixWorkspace_sptr m_inputWS
Input workspace.
static DetectorParams getDetectorParameters(const API::MatrixWorkspace_const_sptr &ws, const size_t index)
Creates a POD struct containing the required detector parameters for this spectrum.
bool convert(const size_t index)
Perform the conversion to Y-space.
static double getComponentParameter(const Geometry::IComponent &det, const Geometry::ParameterMap &pmap, const std::string &name)
Retrieve a component parameter.
void createOutputWorkspace()
Create the output workspace.
static void calculateY(double &yspace, double &qspace, double &ei, const double mass, const double tsec, const double k1, const double v1, const DetectorParams &detpar)
Convert single time value to Y,Q & Ei values.
Holds a collection of detectors.
Definition: DetectorGroup.h:28
base class for Geometric IComponent
Definition: IComponent.h:51
virtual ComponentID getComponentID() const =0
Returns the ComponentID - a unique identifier of the component.
Exception for when an item is not found in a collection.
Definition: Exception.h:145
IPropertyManager * setProperty(const std::string &name, const T &value)
Templated method to set the value of a PropertyWithValue.
void warning(const std::string &msg)
Logs at warning level.
Definition: Logger.cpp:86
static T & Instance()
Return a reference to the Singleton instance, creating it if it does not already exist Creation is do...
std::shared_ptr< const MatrixWorkspace > MatrixWorkspace_const_sptr
shared pointer to the matrix workspace base class (const version)
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
static constexpr double E_mev_toNeutronWavenumberSq
Transformation coefficient to transform neutron energy into neutron wavevector: K-neutron[m^-10] = sq...
static constexpr double NeutronMass
Mass of the neutron in kg.
static constexpr double meV
1 meV in Joules.
Generate a tableworkspace to store the calibration results.
std::string to_string(const wide_integer< Bits, Signed > &n)
Simple data structure to store nominal detector values It avoids some functions taking a huge number ...
double theta
scattering angle in radians
double l1
source-sample distance in metres
double l2
sample-detector distance in metres
@ Input
An input workspace.
Definition: Property.h:53
@ Output
An output workspace.
Definition: Property.h:54