Mantid
Loading...
Searching...
No Matches
ConvertSpectrumAxis.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"
15#include "MantidHistogramData/Histogram.h"
16#include "MantidHistogramData/HistogramBuilder.h"
20#include "MantidKernel/Unit.h"
22
23#include <cfloat>
24
25namespace Mantid::Algorithms {
26// Register the algorithm into the AlgorithmFactory
27DECLARE_ALGORITHM(ConvertSpectrumAxis)
28using namespace Kernel;
29using namespace API;
30using namespace Geometry;
31using namespace DataObjects;
32using namespace HistogramData;
33namespace {
34constexpr double rad2deg = 180. / M_PI;
35}
36
38 // Validator for Input Workspace
39 auto wsVal = std::make_shared<CompositeValidator>();
40 wsVal->add<SpectraAxisValidator>();
41 wsVal->add<InstrumentValidator>();
42
43 declareProperty(std::make_unique<WorkspaceProperty<>>("InputWorkspace", "", Direction::Input, wsVal),
44 "The name of the input workspace.");
45 declareProperty(std::make_unique<WorkspaceProperty<>>("OutputWorkspace", "", Direction::Output),
46 "The name to use for the output workspace.");
47 std::vector<std::string> targetOptions = Mantid::Kernel::UnitFactory::Instance().getKeys();
48 targetOptions.emplace_back("theta");
49 targetOptions.emplace_back("signed_theta");
50 declareProperty("Target", "", std::make_shared<StringListValidator>(targetOptions),
51 "The unit to which the spectrum axis should be converted. "
52 "This can be either \"theta\" (for <math>\\theta</math> "
53 "degrees), or any of the IDs known to the [[Unit Factory]].");
54 std::vector<std::string> eModeOptions;
55 eModeOptions.emplace_back("Direct");
56 eModeOptions.emplace_back("Indirect");
57 declareProperty("EMode", "Direct", std::make_shared<StringListValidator>(eModeOptions),
58 "Some unit conversions require this value to be set "
59 "(\"Direct\" or \"Indirect\")");
60 auto mustBePositive = std::make_shared<BoundedValidator<double>>();
61 mustBePositive->setLower(0.0);
62 declareProperty("EFixed", EMPTY_DBL(), mustBePositive,
63 "Value of fixed energy in meV : EI (EMode=Direct) or EF "
64 "(EMode=Indirect))");
65}
66
68 // Get the input workspace
69 MatrixWorkspace_const_sptr inputWS = getProperty("InputWorkspace");
70 std::string unitTarget = getProperty("Target");
71 // Loop over the original spectrum axis, finding the theta (n.b. not 2theta!)
72 // for each spectrum and storing it's corresponding workspace index
73 // Map will be sorted on theta, so resulting axis will be ordered as well
74 std::multimap<double, size_t> indexMap;
75 const size_t nHist = inputWS->getNumberHistograms();
76 const size_t nBins = inputWS->blocksize();
77 const bool isHist = inputWS->isHistogramData();
78 auto &spectrumInfo = inputWS->spectrumInfo();
79 size_t nxBins;
80 if (isHist) {
81 nxBins = nBins + 1;
82 } else {
83 nxBins = nBins;
84 }
85 if (unitTarget != "theta" && unitTarget != "signed_theta") {
86 Kernel::Unit_sptr fromUnit = inputWS->getAxis(0)->unit();
87 Kernel::Unit_sptr toUnit = UnitFactory::Instance().create(unitTarget);
88 std::vector<double> emptyVector;
89 const double l1 = spectrumInfo.l1();
90 const std::string emodeStr = getProperty("EMode");
91 auto emode = Kernel::DeltaEMode::fromString(emodeStr);
92 size_t nfailures = 0;
93 for (size_t i = 0; i < nHist; i++) {
94 std::vector<double> xval{inputWS->x(i).front(), inputWS->x(i).back()};
96
97 double efixedProp = getProperty("Efixed");
98 if (efixedProp != EMPTY_DBL()) {
99 pmap[UnitParams::efixed] = efixedProp;
100 g_log.debug() << "Detector: " << spectrumInfo.detector(i).getID() << " Efixed: " << efixedProp << "\n";
101 }
102
103 spectrumInfo.getDetectorValues(*fromUnit, *toUnit, emode, false, i, pmap);
104 double value = 0.;
105 try {
106 fromUnit->toTOF(xval, emptyVector, l1, emode, pmap);
107 toUnit->fromTOF(xval, emptyVector, l1, emode, pmap);
108 value = (xval.front() + xval.back()) / 2;
109 } catch (std::runtime_error &) {
110 nfailures++;
111 g_log.warning() << "Unable to calculate new spectrum axis value for "
112 "workspace index "
113 << i;
114 value = inputWS->getAxis(1)->getValue(i);
115 }
116 indexMap.emplace(value, i);
117 }
118 if (nfailures == nHist) {
119 throw std::runtime_error("Unable to convert spectrum axis values on all spectra");
120 }
121 } else {
122 // Set up binding to memeber funtion. Avoids condition as part of loop over
123 // nHistograms.
124 using namespace std::placeholders;
125 std::function<double(const IDetector &)> thetaFunction;
126 if (unitTarget == "signed_theta") {
127 thetaFunction = std::bind(&MatrixWorkspace::detectorSignedTwoTheta, inputWS, _1);
128 } else {
129 thetaFunction = std::bind(&MatrixWorkspace::detectorTwoTheta, inputWS, _1);
130 }
131
132 bool warningGiven = false;
133 for (size_t i = 0; i < nHist; ++i) {
134 try {
135 IDetector_const_sptr det = inputWS->getDetector(i);
136 // Invoke relevant member function.
137 indexMap.emplace(thetaFunction(*det) * rad2deg, i);
138 } catch (Exception::NotFoundError &) {
139 if (!warningGiven)
140 g_log.warning("The instrument definition is incomplete - spectra "
141 "dropped from output");
142 warningGiven = true;
143 }
144 }
145 }
146 // Create the output workspace. Can't re-use the input one because we'll be
147 // re-ordering the spectra.
148 HistogramBuilder builder;
149 builder.setX(nxBins);
150 builder.setY(nBins);
151 builder.setDistribution(inputWS->isDistribution());
152 MatrixWorkspace_sptr outputWS = create<MatrixWorkspace>(*inputWS, indexMap.size(), builder.build());
153 // Now set up a new, numeric axis holding the theta values corresponding to
154 // each spectrum
155 auto newAxis = std::make_unique<NumericAxis>(indexMap.size());
156 auto newAxisRaw = newAxis.get();
157 outputWS->replaceAxis(1, std::move(newAxis));
158 // The unit of this axis is radians. Use the 'radians' unit defined above.
159 if (unitTarget == "theta" || unitTarget == "signed_theta") {
160 newAxisRaw->unit() = std::make_shared<Units::Degrees>();
161 } else {
162 newAxisRaw->unit() = UnitFactory::Instance().create(unitTarget);
163 }
164 std::multimap<double, size_t>::const_iterator it;
165 size_t currentIndex = 0;
166 for (it = indexMap.begin(); it != indexMap.end(); ++it) {
167 // Set the axis value
168 newAxisRaw->setValue(currentIndex, it->first);
169 // Now copy over the data
170 outputWS->setHistogram(currentIndex, inputWS->histogram(it->second));
171 // We can keep the spectrum numbers etc.
172 outputWS->getSpectrum(currentIndex).copyInfoFrom(inputWS->getSpectrum(it->second));
173 ++currentIndex;
174 }
175 setProperty("OutputWorkspace", outputWS);
176}
177
179 const MatrixWorkspace_const_sptr &inputWS, int emode) const {
180 double efixed(0);
181 double efixedProp = getProperty("Efixed");
182 if (efixedProp != EMPTY_DBL()) {
183 efixed = efixedProp;
184 g_log.debug() << "Detector: " << detector.getID() << " Efixed: " << efixed << "\n";
185 } else {
186 if (emode == 1) {
187 if (inputWS->run().hasProperty("Ei")) {
188 Kernel::Property *p = inputWS->run().getProperty("Ei");
189 auto *doublep = dynamic_cast<Kernel::PropertyWithValue<double> *>(p);
190 if (doublep) {
191 efixed = (*doublep)();
192 } else {
193 efixed = 0.0;
194 g_log.warning() << "Efixed could not be found for detector " << detector.getID() << ", set to 0.0\n";
195 }
196 } else {
197 efixed = 0.0;
198 g_log.warning() << "Efixed could not be found for detector " << detector.getID() << ", set to 0.0\n";
199 }
200 } else if (emode == 2) {
201 std::vector<double> efixedVec = detector.getNumberParameter("Efixed");
202 if (efixedVec.empty()) {
203 int detid = detector.getID();
204 IDetector_const_sptr detectorSingle = inputWS->getInstrument()->getDetector(detid);
205 efixedVec = detectorSingle->getNumberParameter("Efixed");
206 }
207 if (!efixedVec.empty()) {
208 efixed = efixedVec.at(0);
209 g_log.debug() << "Detector: " << detector.getID() << " EFixed: " << efixed << "\n";
210 } else {
211 efixed = 0.0;
212 g_log.warning() << "Efixed could not be found for detector " << detector.getID() << ", set to 0.0\n";
213 }
214 }
215 }
216 return efixed;
217}
218} // namespace Mantid::Algorithms
#define DECLARE_ALGORITHM(classname)
Definition: Algorithm.h:576
double value
The value of the point.
Definition: FitMW.cpp:51
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
A validator which checks that a workspace has a valid instrument.
double detectorSignedTwoTheta(const Geometry::IDetector &det) const
Returns the signed 2Theta scattering angle for a detector.
double detectorTwoTheta(const Geometry::IDetector &det) const
Returns the 2Theta scattering angle for a detector.
A validator which checks whether the input workspace has the Spectra number in the axis.
A property class for workspaces.
double getEfixed(const Mantid::Geometry::IDetector &detector, const API::MatrixWorkspace_const_sptr &inputWS, int emode) const
Getting Efixed.
void init() override
Initialisation code.
virtual std::vector< double > getNumberParameter(const std::string &pname, bool recursive=true) const =0
Get a parameter defined as a double.
Interface class for detector objects.
Definition: IDetector.h:43
virtual detid_t getID() const =0
Get the detector ID.
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 debug(const std::string &msg)
Logs at debug level.
Definition: Logger.cpp:114
void warning(const std::string &msg)
Logs at warning level.
Definition: Logger.cpp:86
The concrete, templated class for properties.
Base class for properties.
Definition: Property.h:94
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::shared_ptr< MatrixWorkspace > MatrixWorkspace_sptr
shared pointer to the matrix workspace base class
constexpr double rad2deg
Radians to degrees conversion factor.
Definition: AngleUnits.h:23
std::shared_ptr< const Mantid::Geometry::IDetector > IDetector_const_sptr
Shared pointer to IDetector (const version)
Definition: IDetector.h:102
std::unordered_map< UnitParams, double > UnitParametersMap
Definition: Unit.h:30
std::shared_ptr< Unit > Unit_sptr
Shared pointer to the Unit base class.
Definition: Unit.h:229
constexpr double EMPTY_DBL() noexcept
Returns what we consider an "empty" double within a property.
Definition: EmptyValues.h:43
Generate a tableworkspace to store the calibration results.
static Type fromString(const std::string &modeStr)
Returns the emode from the given string.
Definition: DeltaEMode.cpp:69
@ Input
An input workspace.
Definition: Property.h:53
@ Output
An output workspace.
Definition: Property.h:54