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");
50const static std::string DETECTOR_ID(
"Detector ID");
52const static std::string SPECTRUM_NUMBER(
"Spectrum Number");
54const static std::string WORKSPACE_INDEX(
"Workspace Index");
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");
65const static std::string L2(
"L2");
67const static std::string REFERENCE_SPECTRA(
"ReferenceSpectra");
68const static std::string REFERENCE_WORKSPACE(
"ReferenceWorkspace");
74const static std::string INCIDENT_ENERGY(
"Ei");
75const static std::string WAVELENGTH(
"wavelength");
87template <
typename Map>
size_t mapIndex(
const int index,
const Map &indexMap) {
89 return indexMap.at(
index);
90 }
catch (std::out_of_range &) {
91 throw std::runtime_error(PropertyNames::REFERENCE_SPECTRA +
" out of range.");
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);
114 throw std::runtime_error(PropertyNames::REFERENCE_SPECTRA +
" out of range.");
116 return static_cast<size_t>(
index);
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) {
130 return indexMap.at(i);
131 } catch (std::out_of_range &) {
132 throw std::runtime_error(PropertyNames::REFERENCE_SPECTRA +
" out of range.");
141const std::string CorrectTOFAxis::name()
const {
return "CorrectTOFAxis"; }
144int CorrectTOFAxis::version()
const {
return 1; }
147const std::string CorrectTOFAxis::category()
const {
return "Inelastic\\Corrections"; }
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.";
158void CorrectTOFAxis::init() {
159 auto tofWorkspace = std::make_shared<Kernel::CompositeValidator>();
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);
168 Direction::Input, tofWorkspace),
169 "An input workspace.");
172 "An output workspace.");
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.");
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 +
"').");
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.",
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);
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 +
".";
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 +
".";
217 if (!m_referenceWs->run().hasProperty(SampleLog::INCIDENT_ENERGY)) {
218 issues[PropertyNames::REFERENCE_WORKSPACE] =
"'Ei' is missing from the sample logs.";
220 if (!m_referenceWs->run().hasProperty(SampleLog::WAVELENGTH)) {
221 issues[PropertyNames::REFERENCE_WORKSPACE] =
"'wavelength' is missing from the sample logs.";
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);
237 issues[PropertyNames::REFERENCE_SPECTRA] =
238 "Either " + PropertyNames::REFERENCE_SPECTRA +
" or " + PropertyNames::L2 +
" has to be specified.";
242 m_eppTable = getProperty(PropertyNames::EPP_TABLE);
244 issues[PropertyNames::EPP_TABLE] =
245 "No EPP table specified nor " + PropertyNames::ELASTIC_BIN_INDEX +
" specified.";
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.";
254 if (spectra.empty()) {
255 issues[PropertyNames::REFERENCE_SPECTRA] =
"No reference spectra selected.";
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.";
268 if (!spectrumInfo.hasDetectors(i)) {
269 issues[PropertyNames::REFERENCE_SPECTRA] =
"No detectors attached to workspace index " +
std::to_string(i) +
".";
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.";
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.";
292void CorrectTOFAxis::exec() {
293 m_inputWs = getProperty(PropertyNames::INPUT_WORKSPACE);
295 if (outputWs != m_inputWs) {
296 outputWs = m_inputWs->clone();
299 useReferenceWorkspace(outputWs);
301 correctManually(outputWs);
303 setProperty(PropertyNames::OUTPUT_WORKSPACE, outputWs);
312 const auto histogramCount =
static_cast<int64_t
>(m_referenceWs->getNumberHistograms());
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());
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));
325 if (outputWs->run().hasProperty(SampleLog::WAVELENGTH)) {
326 outputWs->mutableRun()
327 .getProperty(SampleLog::WAVELENGTH)
328 ->setValueFromProperty(*m_referenceWs->run().getProperty(SampleLog::WAVELENGTH));
340 const auto &spectrumInfo = m_inputWs->spectrumInfo();
341 const double l1 = spectrumInfo.l1();
344 auto const fractionalBinIndex =
static_cast<double>(getProperty(PropertyNames::ELASTIC_BIN_INDEX));
345 auto const elasticBinIndexOffset = fractionalBinIndex - floor(fractionalBinIndex);
346 auto eppOffset = 0.0;
348 averageL2AndEPP(spectrumInfo,
l2, epp);
350 epp = m_inputWs->points(0)[m_elasticBinIndex];
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;
357 double Ei = getProperty(PropertyNames::FIXED_ENERGY);
359 Ei = m_inputWs->run().getPropertyAsSingleValue(SampleLog::INCIDENT_ENERGY);
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);
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());
374 for (int64_t i = 0; i < histogramCount; ++i) {
376 outputWs->mutableX(i) += shift;
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);
395 const auto indexCount =
static_cast<int64_t
>(m_workspaceIndices.size());
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) {
405 const double d = spectrumInfo.
l2(
index);
407 const double epp = (*peakPositionColumn)[
index];
410 g_log.
debug() <<
"Including workspace index " <<
index <<
" - distance: " <<
d <<
" EPP: " << epp <<
".\n";
412 g_log.
debug() <<
"Excluding masked workspace index " <<
index <<
".\n";
415 g_log.
debug() <<
"Excluding detector with unsuccessful fit at workspace index " <<
index <<
".\n";
421 throw std::runtime_error(
"No successful detector fits found in " + PropertyNames::EPP_TABLE);
423 l2Average = l2Sum /
static_cast<double>(
n);
425 eppAverage = eppSum /
static_cast<double>(
n);
432 const auto indexCount =
static_cast<int64_t
>(m_workspaceIndices.size());
434 for (int64_t i = 0; i < indexCount; ++i) {
436 const size_t index = m_workspaceIndices[i];
437 interruption_point();
439 const double d = spectrumInfo.
l2(
index);
443 g_log.
debug() <<
"Excluding masked workspace index " <<
index <<
".\n";
449 throw std::runtime_error(
"No unmasked detectors found in " + PropertyNames::REFERENCE_SPECTRA);
451 const double l2 = l2Sum /
static_cast<double>(indexCount);
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;
#define DECLARE_ALGORITHM(classname)
std::map< DeltaEMode::Type, std::string > index
#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...
#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....
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.
void debug(const std::string &msg)
Logs at debug level.
void information(const std::string &msg)
Logs at information level.
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.
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.
std::string to_string(const wide_integer< Bits, Signed > &n)
Describes the direction (within an algorithm) of a Property.