Mantid
Loading...
Searching...
No Matches
CorrectKiKf.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#include "MantidAPI/Run.h"
19
20namespace Mantid::Algorithms {
21
22// Register with the algorithm factory
23DECLARE_ALGORITHM(CorrectKiKf)
24
25using namespace Kernel;
26using namespace API;
27using namespace DataObjects;
28using namespace Geometry;
29using std::size_t;
30
33 auto wsValidator = std::make_shared<WorkspaceUnitValidator>("DeltaE");
34
35 this->declareProperty(
36 std::make_unique<WorkspaceProperty<API::MatrixWorkspace>>("InputWorkspace", "", Direction::Input, wsValidator),
37 "Name of the input workspace");
38 this->declareProperty(
39 std::make_unique<WorkspaceProperty<API::MatrixWorkspace>>("OutputWorkspace", "", Direction::Output),
40 "Name of the output workspace, can be the same as the input");
41
42 std::vector<std::string> propOptions{"Direct", "Indirect"};
43 this->declareProperty("EMode", "Direct", std::make_shared<StringListValidator>(propOptions),
44 "The energy mode (default: Direct)");
45 auto mustBePositive = std::make_shared<BoundedValidator<double>>();
46 mustBePositive->setLower(0.0);
47 this->declareProperty("EFixed", EMPTY_DBL(), mustBePositive,
48 "Value of fixed energy in meV : EI (EMode=Direct) or "
49 "EF (EMode=Indirect) .");
50}
51
53 // Get the workspaces
54 MatrixWorkspace_const_sptr inputWS = this->getProperty("InputWorkspace");
55 MatrixWorkspace_sptr outputWS = this->getProperty("OutputWorkspace");
56
57 // Check if it is an event workspace
58 EventWorkspace_const_sptr eventW = std::dynamic_pointer_cast<const EventWorkspace>(inputWS);
59 if (eventW != nullptr) {
60 this->execEvent();
61 return;
62 }
63
64 // If input and output workspaces are not the same, create a new workspace for
65 // the output using the clone method to account for different workspace types
66 if (outputWS != inputWS) {
67 outputWS = inputWS->clone();
68 }
69
70 const size_t size = inputWS->blocksize();
71 // Calculate the number of spectra in this workspace
72 const auto numberOfSpectra = static_cast<int>(inputWS->size() / size);
73 API::Progress prog(this, 0.0, 1.0, numberOfSpectra);
74 bool negativeEnergyWarning = false;
75
76 const std::string emodeStr = getProperty("EMode");
77 double efixedProp = getProperty("EFixed");
78
79 if (efixedProp == EMPTY_DBL()) {
80 if (emodeStr == "Direct") {
81 // Check if it has been store on the run object for this workspace
82 if (inputWS->run().hasProperty("Ei")) {
83 efixedProp = inputWS->run().getPropertyValueAsType<double>("Ei");
84 g_log.debug() << "Using stored Ei value " << efixedProp << "\n";
85 } else {
86 throw std::invalid_argument("No Ei value has been set or stored within the run information.");
87 }
88 } else {
89 // If not specified, will try to get Ef from the parameter file for
90 // indirect geometry,
91 // but it will be done for each spectrum separately, in case of different
92 // analyzer crystals
93 }
94 }
95
96 // Get the parameter map
97 const ParameterMap &pmap = outputWS->constInstrumentParameters();
98 const auto &spectrumInfo = inputWS->spectrumInfo();
99
100 PARALLEL_FOR_IF(Kernel::threadSafe(*inputWS, *outputWS))
101 for (int64_t i = 0; i < int64_t(numberOfSpectra); ++i) {
103 double Efi = 0;
104 // Now get the detector object for this histogram to check if monitor
105 // or to get Ef for indirect geometry
106 if (emodeStr == "Indirect") {
107 if (efixedProp != EMPTY_DBL())
108 Efi = efixedProp;
109 // If a DetectorGroup is present should provide a value as a property
110 // instead
111 else if (spectrumInfo.hasUniqueDetector(i)) {
112 getEfixedFromParameterMap(Efi, i, spectrumInfo, pmap);
113 } else {
114 g_log.information() << "Workspace Index " << i << ": cannot find detector"
115 << "\n";
116 }
117 }
118
119 auto &yOut = outputWS->mutableY(i);
120 auto &eOut = outputWS->mutableE(i);
121 const auto &xIn = inputWS->points(i);
122 auto &yIn = inputWS->y(i);
123 auto &eIn = inputWS->e(i);
124 // Copy the energy transfer axis
125 outputWS->setSharedX(i, inputWS->sharedX(i));
126 for (unsigned int j = 0; j < size; ++j) {
127 const double deltaE = xIn[j];
128 double Ei = 0.;
129 double Ef = 0.;
130 double kioverkf = 1.;
131 if (emodeStr == "Direct") // Ei=Efixed
132 {
133 Ei = efixedProp;
134 Ef = Ei - deltaE;
135 } else // Ef=Efixed
136 {
137 Ef = Efi;
138 Ei = Efi + deltaE;
139 }
140 // if Ei or Ef is negative, it should be a warning
141 // however, if the intensity is 0 (histogram goes to energy transfer
142 // higher than Ei) it is still ok, so no warning.
143 if ((Ei <= 0) || (Ef <= 0)) {
144 kioverkf = 0.;
145 if (yIn[j] != 0)
146 negativeEnergyWarning = true;
147 } else
148 kioverkf = std::sqrt(Ei / Ef);
149
150 yOut[j] = yIn[j] * kioverkf;
151 eOut[j] = eIn[j] * kioverkf;
152 }
153 prog.report();
155 } // end for i
157
158 if (negativeEnergyWarning)
159 g_log.information() << "Ef <= 0 or Ei <= 0 in at least one spectrum!!!!\n";
160 if ((negativeEnergyWarning) && (efixedProp == EMPTY_DBL()))
161 g_log.information() << "Try to set fixed energy\n";
162 this->setProperty("OutputWorkspace", outputWS);
163}
164
171 g_log.information("Processing event workspace");
172
173 const MatrixWorkspace_const_sptr matrixInputWS = getProperty("InputWorkspace");
174 auto inputWS = std::dynamic_pointer_cast<const EventWorkspace>(matrixInputWS);
175
176 // generate the output workspace pointer
177 API::MatrixWorkspace_sptr matrixOutputWS = getProperty("OutputWorkspace");
178 if (matrixOutputWS != matrixInputWS) {
179 matrixOutputWS = matrixInputWS->clone();
180 setProperty("OutputWorkspace", matrixOutputWS);
181 }
182 auto outputWS = std::dynamic_pointer_cast<EventWorkspace>(matrixOutputWS);
183
184 const std::string emodeStr = getProperty("EMode");
185 double efixedProp = getProperty("EFixed"), efixed;
186
187 if (efixedProp == EMPTY_DBL()) {
188 if (emodeStr == "Direct") {
189 // Check if it has been store on the run object for this workspace
190 if (inputWS->run().hasProperty("Ei")) {
191 efixedProp = inputWS->run().getPropertyValueAsType<double>("Ei");
192 g_log.debug() << "Using stored Ei value " << efixedProp << "\n";
193 } else {
194 throw std::invalid_argument("No Ei value has been set or stored within the run information.");
195 }
196 } else {
197 // If not specified, will try to get Ef from the parameter file for
198 // indirect geometry,
199 // but it will be done for each spectrum separately, in case of different
200 // analyzer crystals
201 }
202 }
203
204 // Get the parameter map
205 const ParameterMap &pmap = outputWS->constInstrumentParameters();
206
207 auto numHistograms = static_cast<int64_t>(inputWS->getNumberHistograms());
208 const auto &spectrumInfo = inputWS->spectrumInfo();
209 API::Progress prog = API::Progress(this, 0.0, 1.0, numHistograms);
211 for (int64_t i = 0; i < numHistograms; ++i) {
213 // Now get the detector object for this histogram to check if monitor
214 // or to get Ef for indirect geometry
215 if (emodeStr == "Indirect") {
216 double Efi = 0;
217 if (efixedProp != EMPTY_DBL()) {
218 Efi = efixedProp;
219 // If a DetectorGroup is present should provide a value as a property
220 // instead
221 } else if (spectrumInfo.hasUniqueDetector(i)) {
222 getEfixedFromParameterMap(Efi, i, spectrumInfo, pmap);
223 } else {
224 g_log.information() << "Workspace Index " << i << ": cannot find detector"
225 << "\n";
226 }
227 efixed = Efi;
228 } else
229 efixed = efixedProp;
230
231 // Do the correction
232 auto &evlist = outputWS->getSpectrum(i);
233 switch (evlist.getEventType()) {
234 case TOF:
235 // Switch to weights if needed.
236 evlist.switchTo(WEIGHTED);
237 /* no break */
238 // Fall through
239
240 case WEIGHTED:
241 correctKiKfEventHelper(evlist.getWeightedEvents(), efixed, emodeStr);
242 break;
243
244 case WEIGHTED_NOTIME:
245 correctKiKfEventHelper(evlist.getWeightedEventsNoTime(), efixed, emodeStr);
246 break;
247 }
248
249 prog.report();
251 }
253
254 outputWS->clearMRU();
255 if (inputWS->getNumberEvents() != outputWS->getNumberEvents()) {
256 g_log.information() << "Ef <= 0 or Ei <= 0 for " << inputWS->getNumberEvents() - outputWS->getNumberEvents()
257 << " events, out of " << inputWS->getNumberEvents() << '\n';
258 if (efixedProp == EMPTY_DBL())
259 g_log.information() << "Try to set fixed energy\n";
260 }
261}
262
263template <class T>
264void CorrectKiKf::correctKiKfEventHelper(std::vector<T> &wevector, double efixed, const std::string &emodeStr) {
265 double Ei, Ef;
266 float kioverkf;
267 typename std::vector<T>::iterator it;
268 for (it = wevector.begin(); it < wevector.end();) {
269 if (emodeStr == "Direct") // Ei=Efixed
270 {
271 Ei = efixed;
272 Ef = Ei - it->tof();
273 } else // Ef=Efixed
274 {
275 Ef = efixed;
276 Ei = Ef + it->tof();
277 }
278 // if Ei or Ef is negative, delete the event
279 if ((Ei <= 0) || (Ef <= 0)) {
280 it = wevector.erase(it);
281 } else {
282 kioverkf = static_cast<float>(std::sqrt(Ei / Ef));
283 it->m_weight *= kioverkf;
284 it->m_errorSquared *= kioverkf * kioverkf;
285 ++it;
286 }
287 }
288}
289
290void CorrectKiKf::getEfixedFromParameterMap(double &Efi, int64_t i, const SpectrumInfo &spectrumInfo,
291 const ParameterMap &pmap) {
292 Efi = 0;
293
294 if (spectrumInfo.isMonitor(i))
295 return;
296
297 const auto &det = spectrumInfo.detector(i);
298 Parameter_sptr par = pmap.getRecursive(&det, "Efixed");
299 if (par) {
300 Efi = par->value<double>();
301 g_log.debug() << "Detector: " << det.getID() << " EFixed: " << Efi << "\n";
302 }
303}
304
305} // namespace Mantid::Algorithms
#define DECLARE_ALGORITHM(classname)
Definition: Algorithm.h:576
#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 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.
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
Helper class for reporting progress from algorithms.
Definition: Progress.h:25
API::SpectrumInfo is an intermediate step towards a SpectrumInfo that is part of Instrument-2....
Definition: SpectrumInfo.h:53
bool isMonitor(const size_t index) const
Returns true if the detector(s) associated with the spectrum are monitors.
const Geometry::IDetector & detector(const size_t index) const
Return a const reference to the detector or detector group of the spectrum with given index.
A property class for workspaces.
void getEfixedFromParameterMap(double &Efi, int64_t i, const Mantid::API::SpectrumInfo &spectrumInfo, const Mantid::Geometry::ParameterMap &pmap)
void execEvent()
Execute CorrectKiKf for event workspaces.
void exec() override
Virtual method - must be overridden by concrete algorithm.
Definition: CorrectKiKf.cpp:52
void correctKiKfEventHelper(std::vector< T > &wevector, double efixed, const std::string &emodeStr)
Execute CorrectKiKf for event lists.
void init() override
Initialisation method.
Definition: CorrectKiKf.cpp:32
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 information(const std::string &msg)
Logs at information level.
Definition: Logger.cpp:105
void report()
Increments the loop counter by 1, then sends the progress notification on behalf of its algorithm.
Definition: ProgressBase.h:51
std::shared_ptr< const MatrixWorkspace > MatrixWorkspace_const_sptr
shared pointer to the matrix workspace base class (const version)
@ WEIGHTED_NOTIME
Definition: IEventList.h:18
std::shared_ptr< MatrixWorkspace > MatrixWorkspace_sptr
shared pointer to the matrix workspace base class
std::shared_ptr< const EventWorkspace > EventWorkspace_const_sptr
shared pointer to a const Workspace2D
std::shared_ptr< Parameter > Parameter_sptr
Typedef for the shared pointer.
Definition: Parameter.h:195
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
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.
@ Input
An input workspace.
Definition: Property.h:53
@ Output
An output workspace.
Definition: Property.h:54