Mantid
Loading...
Searching...
No Matches
CorrectTOFAxis.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
12#include "MantidAPI/Run.h"
23
24namespace Mantid::Algorithms {
25
28
29// Register the algorithm into the AlgorithmFactory
30DECLARE_ALGORITHM(CorrectTOFAxis)
31
32namespace {
36namespace EPPTableLiterals {
38const static std::string FIT_STATUS_COLUMN("FitStatus");
40const static std::string PEAK_CENTRE_COLUMN("PeakCentre");
42const static std::string FIT_STATUS_SUCCESS("success");
43} // namespace EPPTableLiterals
44
48namespace IndexTypes {
50const static std::string DETECTOR_ID("Detector ID");
52const static std::string SPECTRUM_NUMBER("Spectrum Number");
54const static std::string WORKSPACE_INDEX("Workspace Index");
55} // namespace IndexTypes
56
59namespace PropertyNames {
60const static std::string ELASTIC_BIN_INDEX("ElasticBinIndex");
61const static std::string EPP_TABLE("EPPTable");
62const static std::string FIXED_ENERGY("EFixed");
63const static std::string INDEX_TYPE("IndexType");
64const static std::string INPUT_WORKSPACE("InputWorkspace");
65const static std::string L2("L2");
66const static std::string OUTPUT_WORKSPACE("OutputWorkspace");
67const static std::string REFERENCE_SPECTRA("ReferenceSpectra");
68const static std::string REFERENCE_WORKSPACE("ReferenceWorkspace");
69} // namespace PropertyNames
70
73namespace SampleLog {
74const static std::string INCIDENT_ENERGY("Ei");
75const static std::string WAVELENGTH("wavelength");
76} // namespace SampleLog
77
87template <typename Map> size_t mapIndex(const int index, const Map &indexMap) {
88 try {
89 return indexMap.at(index);
90 } catch (std::out_of_range &) {
91 throw std::runtime_error(PropertyNames::REFERENCE_SPECTRA + " out of range.");
92 }
93}
94
105size_t toWorkspaceIndex(const int index, const std::string &indexType, const API::MatrixWorkspace_const_sptr &ws) {
106 if (indexType == IndexTypes::DETECTOR_ID) {
107 const auto indexMap = ws->getDetectorIDToWorkspaceIndexMap();
108 return mapIndex(index, indexMap);
109 } else if (indexType == IndexTypes::SPECTRUM_NUMBER) {
110 const auto indexMap = ws->getSpectrumToWorkspaceIndexMap();
111 return mapIndex(index, indexMap);
112 } else {
113 if (index < 0) {
114 throw std::runtime_error(PropertyNames::REFERENCE_SPECTRA + " out of range.");
115 }
116 return static_cast<size_t>(index);
117 }
118}
119
125template <typename Map>
126void mapIndices(const std::vector<int> &spectra, const Map &indexMap, std::vector<size_t> &workspaceIndices) {
127 auto back = std::back_inserter(workspaceIndices);
128 std::transform(spectra.cbegin(), spectra.cend(), back, [&indexMap](int i) {
129 try {
130 return indexMap.at(i);
131 } catch (std::out_of_range &) {
132 throw std::runtime_error(PropertyNames::REFERENCE_SPECTRA + " out of range.");
133 }
134 });
135}
136} // anonymous namespace
137
138//----------------------------------------------------------------------------------------------
139
141const std::string CorrectTOFAxis::name() const { return "CorrectTOFAxis"; }
142
144int CorrectTOFAxis::version() const { return 1; }
145
147const std::string CorrectTOFAxis::category() const { return "Inelastic\\Corrections"; }
148
150const std::string CorrectTOFAxis::summary() const {
151 return "Corrects the time-of-flight axis with regards to the incident energy "
152 "and the L1+L2 distance or a reference workspace.";
153}
154
155//----------------------------------------------------------------------------------------------
158void CorrectTOFAxis::init() {
159 auto tofWorkspace = std::make_shared<Kernel::CompositeValidator>();
160 tofWorkspace->add<API::WorkspaceUnitValidator>("TOF");
161 tofWorkspace->add<API::InstrumentValidator>();
162 auto mustBePositiveDouble = std::make_shared<Kernel::BoundedValidator<double>>();
163 mustBePositiveDouble->setLower(0);
164 auto mustBePositiveInt = std::make_shared<Kernel::BoundedValidator<int>>();
165 mustBePositiveInt->setLower(0);
166
167 declareProperty(std::make_unique<WorkspaceProperty<API::MatrixWorkspace>>(PropertyNames::INPUT_WORKSPACE, "",
168 Direction::Input, tofWorkspace),
169 "An input workspace.");
170 declareProperty(
171 std::make_unique<WorkspaceProperty<API::MatrixWorkspace>>(PropertyNames::OUTPUT_WORKSPACE, "", Direction::Output),
172 "An output workspace.");
173 declareProperty(std::make_unique<WorkspaceProperty<API::MatrixWorkspace>>(PropertyNames::REFERENCE_WORKSPACE, "",
174 Direction::Input,
175 API::PropertyMode::Optional, tofWorkspace),
176 "A reference workspace from which to copy the TOF axis as "
177 "well as the 'Ei' and 'wavelength' sample logs.");
178 declareProperty(std::make_unique<WorkspaceProperty<API::ITableWorkspace>>(
179 PropertyNames::EPP_TABLE.c_str(), "", Direction::Input, API::PropertyMode::Optional),
180 "An input EPP table.");
181 const std::vector<std::string> indexTypes{IndexTypes::DETECTOR_ID, IndexTypes::SPECTRUM_NUMBER,
182 IndexTypes::WORKSPACE_INDEX};
183 declareProperty(PropertyNames::INDEX_TYPE, IndexTypes::DETECTOR_ID,
184 std::make_shared<Kernel::StringListValidator>(indexTypes),
185 "The type of indices used in " + PropertyNames::REFERENCE_SPECTRA + " or " +
186 PropertyNames::ELASTIC_BIN_INDEX + " (default: '" + IndexTypes::DETECTOR_ID + "').");
187 declareProperty(std::make_unique<Kernel::ArrayProperty<int>>(PropertyNames::REFERENCE_SPECTRA.c_str()),
188 "A list of reference spectra.");
189 declareProperty(PropertyNames::ELASTIC_BIN_INDEX, EMPTY_DBL(), mustBePositiveDouble,
190 "Bin index of the nominal elastic TOF channel.", Direction::Input);
191 declareProperty(PropertyNames::FIXED_ENERGY, EMPTY_DBL(), mustBePositiveDouble,
192 "Incident energy if the 'EI' sample log is not present/incorrect.", Direction::Input);
193 declareProperty(PropertyNames::L2, EMPTY_DBL(), mustBePositiveDouble, "Sample to detector distance, in meters.",
194 Direction::Input);
195}
196
200std::map<std::string, std::string> CorrectTOFAxis::validateInputs() {
201 std::map<std::string, std::string> issues;
202 m_inputWs = getProperty(PropertyNames::INPUT_WORKSPACE);
203 m_referenceWs = getProperty(PropertyNames::REFERENCE_WORKSPACE);
204 if (m_referenceWs) {
205 m_referenceWs = getProperty(PropertyNames::REFERENCE_WORKSPACE);
206 if (m_inputWs->getNumberHistograms() != m_referenceWs->getNumberHistograms()) {
207 issues[PropertyNames::REFERENCE_WORKSPACE] =
208 "Number of histograms don't match with" + PropertyNames::INPUT_WORKSPACE + ".";
209 }
210 for (size_t i = 0; i < m_inputWs->getNumberHistograms(); ++i) {
211 if (m_inputWs->x(i).size() != m_referenceWs->x(i).size()) {
212 issues[PropertyNames::REFERENCE_WORKSPACE] =
213 "X axis sizes don't match with " + PropertyNames::INPUT_WORKSPACE + ".";
214 break;
215 }
216 }
217 if (!m_referenceWs->run().hasProperty(SampleLog::INCIDENT_ENERGY)) {
218 issues[PropertyNames::REFERENCE_WORKSPACE] = "'Ei' is missing from the sample logs.";
219 }
220 if (!m_referenceWs->run().hasProperty(SampleLog::WAVELENGTH)) {
221 issues[PropertyNames::REFERENCE_WORKSPACE] = "'wavelength' is missing from the sample logs.";
222 }
223 // If reference workspace is given, the rest of the properties are
224 // skipped.
225 return issues;
226 }
227 // If no reference workspace, we either use a predefined elastic channel
228 // or EPP tables to declare the elastic TOF.
229 const double elasticBinFullIndex = getProperty(PropertyNames::ELASTIC_BIN_INDEX);
230 const int elasticBinIndex = static_cast<int>(elasticBinFullIndex);
231 const std::vector<int> spectra = getProperty(PropertyNames::REFERENCE_SPECTRA);
232 const double l2 = getProperty(PropertyNames::L2);
233 if (elasticBinFullIndex != EMPTY_DBL()) {
234 const std::string indexType = getProperty(PropertyNames::INDEX_TYPE);
235 m_elasticBinIndex = toWorkspaceIndex(elasticBinIndex, indexType, m_inputWs);
236 if (spectra.empty() && l2 == EMPTY_DBL()) {
237 issues[PropertyNames::REFERENCE_SPECTRA] =
238 "Either " + PropertyNames::REFERENCE_SPECTRA + " or " + PropertyNames::L2 + " has to be specified.";
239 return issues;
240 }
241 } else {
242 m_eppTable = getProperty(PropertyNames::EPP_TABLE);
243 if (!m_eppTable) {
244 issues[PropertyNames::EPP_TABLE] =
245 "No EPP table specified nor " + PropertyNames::ELASTIC_BIN_INDEX + " specified.";
246 return issues;
247 }
248 const auto peakPositionColumn = m_eppTable->getColumn(EPPTableLiterals::PEAK_CENTRE_COLUMN);
249 const auto fitStatusColumn = m_eppTable->getColumn(EPPTableLiterals::FIT_STATUS_COLUMN);
250 if (!peakPositionColumn || !fitStatusColumn) {
251 issues[PropertyNames::EPP_TABLE] = "EPP table doesn't contain the expected columns.";
252 return issues;
253 }
254 if (spectra.empty()) {
255 issues[PropertyNames::REFERENCE_SPECTRA] = "No reference spectra selected.";
256 return issues;
257 }
258 }
259 m_workspaceIndices = referenceWorkspaceIndices();
260 std::sort(m_workspaceIndices.begin(), m_workspaceIndices.end());
261 m_workspaceIndices.erase(std::unique(m_workspaceIndices.begin(), m_workspaceIndices.end()), m_workspaceIndices.end());
262 const auto &spectrumInfo = m_inputWs->spectrumInfo();
263 for (const auto i : m_workspaceIndices) {
264 if (spectrumInfo.isMonitor(i)) {
265 issues[PropertyNames::REFERENCE_SPECTRA] = "Monitor found among the given spectra.";
266 break;
267 }
268 if (!spectrumInfo.hasDetectors(i)) {
269 issues[PropertyNames::REFERENCE_SPECTRA] = "No detectors attached to workspace index " + std::to_string(i) + ".";
270 break;
271 }
272 if (m_eppTable) {
273 const auto peakPositionColumn = m_eppTable->getColumn(EPPTableLiterals::PEAK_CENTRE_COLUMN);
274 if (i >= peakPositionColumn->size()) {
275 issues[PropertyNames::REFERENCE_SPECTRA] =
276 "Workspace index " + std::to_string(i) + " not found in the EPP table.";
277 }
278 }
279 }
280
281 if (getPointerToProperty(PropertyNames::FIXED_ENERGY)->isDefault()) {
282 if (!m_inputWs->run().hasProperty(SampleLog::INCIDENT_ENERGY)) {
283 issues[PropertyNames::INPUT_WORKSPACE] = "'Ei' is missing from the sample logs.";
284 }
285 }
286 return issues;
287}
288
289//----------------------------------------------------------------------------------------------
292void CorrectTOFAxis::exec() {
293 m_inputWs = getProperty(PropertyNames::INPUT_WORKSPACE);
294 API::MatrixWorkspace_sptr outputWs = getProperty(PropertyNames::OUTPUT_WORKSPACE);
295 if (outputWs != m_inputWs) {
296 outputWs = m_inputWs->clone();
297 }
298 if (m_referenceWs) {
299 useReferenceWorkspace(outputWs);
300 } else {
301 correctManually(outputWs);
302 }
303 setProperty(PropertyNames::OUTPUT_WORKSPACE, outputWs);
304}
305
311void CorrectTOFAxis::useReferenceWorkspace(const API::MatrixWorkspace_sptr &outputWs) {
312 const auto histogramCount = static_cast<int64_t>(m_referenceWs->getNumberHistograms());
313 PARALLEL_FOR_IF(threadSafe(*m_referenceWs, *outputWs))
314 for (int64_t i = 0; i < histogramCount; ++i) {
316 std::copy(m_referenceWs->x(i).cbegin(), m_referenceWs->x(i).cend(), outputWs->mutableX(i).begin());
318 }
320 if (outputWs->run().hasProperty(SampleLog::INCIDENT_ENERGY)) {
321 outputWs->mutableRun()
322 .getProperty(SampleLog::INCIDENT_ENERGY)
323 ->setValueFromProperty(*m_referenceWs->run().getProperty(SampleLog::INCIDENT_ENERGY));
324 }
325 if (outputWs->run().hasProperty(SampleLog::WAVELENGTH)) {
326 outputWs->mutableRun()
327 .getProperty(SampleLog::WAVELENGTH)
328 ->setValueFromProperty(*m_referenceWs->run().getProperty(SampleLog::WAVELENGTH));
329 }
330}
331
339void CorrectTOFAxis::correctManually(const API::MatrixWorkspace_sptr &outputWs) {
340 const auto &spectrumInfo = m_inputWs->spectrumInfo();
341 const double l1 = spectrumInfo.l1();
342 auto l2 = 0.0;
343 auto epp = 0.0;
344 auto const fractionalBinIndex = static_cast<double>(getProperty(PropertyNames::ELASTIC_BIN_INDEX));
345 auto const elasticBinIndexOffset = fractionalBinIndex - floor(fractionalBinIndex);
346 auto eppOffset = 0.0;
347 if (m_eppTable) {
348 averageL2AndEPP(spectrumInfo, l2, epp);
349 } else {
350 epp = m_inputWs->points(0)[m_elasticBinIndex];
351 eppOffset =
352 elasticBinIndexOffset * (m_inputWs->points(0)[m_elasticBinIndex + 1] - m_inputWs->points(0)[m_elasticBinIndex]);
353 const double l2Property = getProperty(PropertyNames::L2);
354 l2 = l2Property == EMPTY_DBL() ? averageL2(spectrumInfo) : l2Property;
355 }
356 g_log.information() << "EPP: " << epp << ".\n";
357 double Ei = getProperty(PropertyNames::FIXED_ENERGY);
358 if (Ei == EMPTY_DBL()) {
359 Ei = m_inputWs->run().getPropertyAsSingleValue(SampleLog::INCIDENT_ENERGY);
360 } else {
361 // Save user-given Ei and wavelength to the output workspace.
362 outputWs->mutableRun().addProperty(SampleLog::INCIDENT_ENERGY, Ei, true);
363 const double wavelength =
364 Kernel::UnitConversion::run("Energy", "Wavelength", Ei, l1, l2, 0, Kernel::DeltaEMode::Direct, 0);
365 outputWs->mutableRun().addProperty(SampleLog::WAVELENGTH, wavelength, true);
366 }
367 // In microseconds.
368 const double TOF = (l1 + l2) / std::sqrt(2 * Ei * PhysicalConstants::meV / PhysicalConstants::NeutronMass) * 1e6;
369 g_log.information() << "Calculated TOF for L1+L2 distance of " << l1 + l2 << "m: " << TOF << '\n';
370 const double shift = TOF - epp - eppOffset;
371 g_log.debug() << "TOF shift: " << shift << '\n';
372 const auto histogramCount = static_cast<int64_t>(m_inputWs->getNumberHistograms());
373 PARALLEL_FOR_IF(threadSafe(*m_inputWs, *outputWs))
374 for (int64_t i = 0; i < histogramCount; ++i) {
376 outputWs->mutableX(i) += shift;
378 }
380}
381
389void CorrectTOFAxis::averageL2AndEPP(const API::SpectrumInfo &spectrumInfo, double &l2Average, double &eppAverage) {
390 auto peakPositionColumn = m_eppTable->getColumn(EPPTableLiterals::PEAK_CENTRE_COLUMN);
391 auto fitStatusColumn = m_eppTable->getColumn(EPPTableLiterals::FIT_STATUS_COLUMN);
392 double l2Sum = 0;
393 double eppSum = 0;
394 size_t n = 0;
395 const auto indexCount = static_cast<int64_t>(m_workspaceIndices.size());
396
397 PRAGMA_OMP(parallel for if (m_eppTable->threadSafe())
398 reduction(+: n, l2Sum, eppSum))
399 for (int64_t i = 0; i < indexCount; ++i) {
401 const size_t index = m_workspaceIndices[i];
402 interruption_point();
403 if (fitStatusColumn->cell<std::string>(index) == EPPTableLiterals::FIT_STATUS_SUCCESS) {
404 if (!spectrumInfo.isMasked(index)) {
405 const double d = spectrumInfo.l2(index);
406 l2Sum += d;
407 const double epp = (*peakPositionColumn)[index];
408 eppSum += epp;
409 ++n;
410 g_log.debug() << "Including workspace index " << index << " - distance: " << d << " EPP: " << epp << ".\n";
411 } else {
412 g_log.debug() << "Excluding masked workspace index " << index << ".\n";
413 }
414 } else {
415 g_log.debug() << "Excluding detector with unsuccessful fit at workspace index " << index << ".\n";
416 }
418 }
420 if (n == 0) {
421 throw std::runtime_error("No successful detector fits found in " + PropertyNames::EPP_TABLE);
422 }
423 l2Average = l2Sum / static_cast<double>(n);
424 g_log.information() << "Average L2 distance: " << l2Average << ".\n";
425 eppAverage = eppSum / static_cast<double>(n);
426 g_log.information() << "Average EPP: " << eppAverage << ".\n";
427}
428
429double CorrectTOFAxis::averageL2(const API::SpectrumInfo &spectrumInfo) {
430 double l2Sum = 0;
431 size_t n = 0;
432 const auto indexCount = static_cast<int64_t>(m_workspaceIndices.size());
433 PRAGMA_OMP(parallel for reduction(+: n, l2Sum))
434 for (int64_t i = 0; i < indexCount; ++i) {
436 const size_t index = m_workspaceIndices[i];
437 interruption_point();
438 if (!spectrumInfo.isMasked(index)) {
439 const double d = spectrumInfo.l2(index);
440 ++n;
441 l2Sum += d;
442 } else {
443 g_log.debug() << "Excluding masked workspace index " << index << ".\n";
444 }
446 }
448 if (n == 0) {
449 throw std::runtime_error("No unmasked detectors found in " + PropertyNames::REFERENCE_SPECTRA);
450 }
451 const double l2 = l2Sum / static_cast<double>(indexCount);
452 g_log.information() << "Average L2 distance: " << l2 << ".\n";
453 return l2;
454}
455
459std::vector<size_t> CorrectTOFAxis::referenceWorkspaceIndices() const {
460 const std::vector<int> indices = getProperty(PropertyNames::REFERENCE_SPECTRA);
461 const std::string indexType = getProperty(PropertyNames::INDEX_TYPE);
462 std::vector<size_t> workspaceIndices(indices.size());
463 std::transform(indices.cbegin(), indices.cend(), workspaceIndices.begin(),
464 [&indexType, this](int index) { return toWorkspaceIndex(index, indexType, m_inputWs); });
465 return workspaceIndices;
466}
467
468} // namespace Mantid::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_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 PRAGMA_OMP(expression)
#define PARALLEL_CHECK_INTERRUPT_REGION
Adds a check after a Parallel region to see if it was interupted.
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
double l2(const size_t index) const
Returns L2 (distance from sample to spectrum).
bool isMasked(const size_t index) const
Returns true if the detector(s) associated with the spectrum are masked.
A property class for workspaces.
A validator which checks that the unit of the workspace referred to by a WorkspaceProperty is the exp...
Support for a property that holds an array of values.
Definition: ArrayProperty.h:28
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
Kernel::Logger g_log("ExperimentInfo")
static logger object
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
static const std::string OUTPUT_WORKSPACE
static const std::string INPUT_WORKSPACE
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 NeutronMass
Mass of the neutron in kg.
static constexpr double meV
1 meV in Joules.
constexpr double EMPTY_DBL() noexcept
Returns what we consider an "empty" double within a property.
Definition: EmptyValues.h:43
std::string to_string(const wide_integer< Bits, Signed > &n)
Describes the direction (within an algorithm) of a Property.
Definition: Property.h:50