Mantid
Loading...
Searching...
No Matches
GetEiMonDet3.cpp
Go to the documentation of this file.
1// Mantid Repository : https://github.com/mantidproject/mantid
2//
3// Copyright © 2019 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/Algorithm.tcc"
11#include "MantidAPI/Run.h"
15#include "MantidIndexing/SpectrumIndexSet.h"
25
26namespace {
30namespace Prop {
32const static std::string DETECTOR_WORKSPACE("DetectorWorkspace");
34const static std::string INCIDENT_ENERGY("IncidentEnergy");
36const static std::string MONITOR("MonitorIndex");
38const static std::string MONITOR_EPP_TABLE("MonitorEPPTable");
40const static std::string MONITOR_WORKSPACE("MonitorWorkspace");
42const static std::string PULSE_INTERVAL("PulseInterval");
44const static std::string MAX_ENERGY("MaximumEnergy");
45} // namespace Prop
46
47std::vector<size_t> toWorkspaceIndices(const Mantid::Indexing::SpectrumIndexSet &indices) {
48 std::vector<size_t> wsIndices;
49 wsIndices.reserve(indices.size());
50 for (size_t i = 0; i < indices.size(); ++i) {
51 wsIndices.emplace_back(indices[i]);
52 }
53 return wsIndices;
54}
55} // namespace
56
57namespace Mantid::Algorithms {
58
59using namespace API;
60
61// Register the algorithm into the algorithm factory.
62DECLARE_ALGORITHM(GetEiMonDet3)
63
64
65const std::string GetEiMonDet3::name() const { return "GetEiMonDet"; }
66
68const std::string GetEiMonDet3::summary() const {
69 return "Calculates the kinetic energy of neutrons leaving the source based "
70 "on the time it takes for them to travel between a monitor and a "
71 "set of detectors.";
72}
73
75int GetEiMonDet3::version() const { return 3; }
76
77const std::vector<std::string> GetEiMonDet3::seeAlso() const { return {"GetEi"}; }
78
80const std::string GetEiMonDet3::category() const { return "Inelastic\\Ei"; }
81
84namespace SampleLogs {
86const static std::string PULSE_INTERVAL("pulse_interval");
87} // namespace SampleLogs
88
93 auto tofWorkspace = std::make_shared<Kernel::CompositeValidator>();
94 tofWorkspace->add<API::WorkspaceUnitValidator>("TOF");
95 tofWorkspace->add<API::InstrumentValidator>();
96 auto mandatoryIntProperty = std::make_shared<Kernel::MandatoryValidator<int>>();
97 auto mustBePositive = std::make_shared<Kernel::BoundedValidator<double>>();
98 mustBePositive->setLower(0);
99
100 declareWorkspaceInputProperties<API::MatrixWorkspace, static_cast<int>(IndexType::SpectrumNum) |
101 static_cast<int>(IndexType::WorkspaceIndex)>(
102 Prop::DETECTOR_WORKSPACE, "A workspace containing the detector spectra.", tofWorkspace);
104 std::make_unique<API::WorkspaceProperty<>>(Prop::MONITOR_WORKSPACE.c_str(), "", Kernel::Direction::Input,
105 API::PropertyMode::Optional, tofWorkspace),
106 "A Workspace containing the monitor spectrum; if empty, " + Prop::DETECTOR_WORKSPACE + " will be used.");
108 Prop::MONITOR_EPP_TABLE, "", Kernel::Direction::Input, API::PropertyMode::Optional),
109 "An EPP table corresponding to " + Prop::MONITOR_WORKSPACE);
110 setPropertySettings(Prop::MONITOR_EPP_TABLE,
111 std::make_unique<Kernel::EnabledWhenProperty>(Prop::MONITOR_WORKSPACE, Kernel::IS_NOT_DEFAULT));
112 declareProperty(Prop::MONITOR, EMPTY_INT(), mandatoryIntProperty, "Usable monitor's workspace index.");
113 declareProperty(Prop::PULSE_INTERVAL, EMPTY_DBL(),
114 "Interval between neutron pulses, in microseconds; taken "
115 "from the sample logs, if not specified.");
116 declareProperty(Prop::MAX_ENERGY, EMPTY_DBL(), mustBePositive,
117 "Multiple pulse intervals will be added to the flight time "
118 "the until final energy is less than this value.");
119 declareProperty(Prop::INCIDENT_ENERGY, EMPTY_DBL(), mustBePositive, "Calculated incident energy, in meV.",
121}
122
127 progress(0.);
128 API::MatrixWorkspace_sptr detectorWs;
129 Indexing::SpectrumIndexSet detectorIndices;
130 std::tie(detectorWs, detectorIndices) = getWorkspaceAndIndices<API::MatrixWorkspace>(Prop::DETECTOR_WORKSPACE);
131
132 API::MatrixWorkspace_sptr monitorWs = getProperty(Prop::MONITOR_WORKSPACE);
133 const int monitorIndex = getProperty(Prop::MONITOR);
134 if (!monitorWs) {
135 monitorWs = detectorWs;
136 if (std::find(detectorIndices.begin(), detectorIndices.end(), monitorIndex) != detectorIndices.end()) {
137 throw std::runtime_error("MonitorIndex is also listed in DetectorWorkspaceIndexSet.");
138 }
139 }
140 if (!monitorWs->spectrumInfo().isMonitor(monitorIndex)) {
141 m_log.warning("The monitor spectrum is not marked as a monitor.");
142 }
143 const auto detWsIndices = toWorkspaceIndices(detectorIndices);
144 auto detectorSumWs = groupSpectra(detectorWs, detWsIndices);
145 progress(0.3);
146 const auto sampleToDetectorDistance = detectorSumWs->spectrumInfo().l2(0);
147 double detectorEPP;
148 try {
149 detectorEPP = peakPosition(detectorSumWs);
150 } catch (std::runtime_error &e) {
151 throw std::runtime_error(std::string("Failed to find detector peak for incident energy: ") + e.what());
152 }
153 progress(0.5);
154 const std::vector<size_t> monWsIndices = {static_cast<size_t>(monitorIndex)};
155 auto monitorSumWs = groupSpectra(monitorWs, monWsIndices);
156 double monitorEPP;
157 try {
158 if (isDefault(Prop::MONITOR_EPP_TABLE)) {
159 monitorEPP = peakPosition(monitorSumWs);
160 } else {
161 monitorEPP = monitorPeakPosition(monWsIndices.front());
162 }
163 } catch (std::runtime_error &e) {
164 throw std::runtime_error(std::string("Failed to find monitor peak for incident energy: ") + e.what());
165 }
166 progress(0.7);
167 // SpectrumInfo returns a negative l2 for monitor.
168 const auto monitorToSampleDistance = std::abs(monitorSumWs->spectrumInfo().l2(0));
169 const double minTOF = minimumTOF(*detectorWs, sampleToDetectorDistance);
170
171 double timeOfFlight = computeTOF(*detectorWs, detectorEPP, monitorEPP, minTOF);
172 const double flightLength = sampleToDetectorDistance + monitorToSampleDistance;
173 const double velocity = flightLength / timeOfFlight * 1e6;
174 using namespace PhysicalConstants;
175 const double energy = 0.5 * NeutronMass * velocity * velocity / meV;
176 progress(1.0);
177 g_log.notice() << "Final time-of-flight:" << timeOfFlight << " which gives " << energy << " as "
178 << Prop::INCIDENT_ENERGY << ".\n";
179 setProperty(Prop::INCIDENT_ENERGY, energy);
180}
181
191double GetEiMonDet3::computeTOF(const API::MatrixWorkspace &detectorWs, const double detectorEPP,
192 const double monitorEPP, const double minTOF) {
193 double timeOfFlight = detectorEPP - monitorEPP;
194 // Check if the obtained time-of-flight makes any sense.
195 while (timeOfFlight <= minTOF) {
196 double pulseInterval = getProperty(Prop::PULSE_INTERVAL);
197 if (pulseInterval == EMPTY_DBL()) {
198 if (detectorWs.run().hasProperty(SampleLogs::PULSE_INTERVAL)) {
199 pulseInterval = detectorWs.run().getPropertyAsSingleValue(SampleLogs::PULSE_INTERVAL);
200 pulseInterval *= 1e6; // To microseconds.
201 } else {
202 throw std::invalid_argument("PulseInterval not explicitly given nor found in the sample logs.");
203 }
204 }
205 g_log.notice() << "Frame delay of " << pulseInterval << " microseconds will be added to the time-of-flight.\n";
206 timeOfFlight += pulseInterval;
207 }
208 g_log.notice() << "Calculated time-of-flight: " << timeOfFlight << ".\n";
209 return timeOfFlight;
210}
211
219 const std::vector<size_t> &wsIndices) {
220 auto group = createChildAlgorithm("GroupDetectors");
221 group->setProperty("InputWorkspace", ws);
222 group->setProperty("OutputWorkspace", "unused");
223 group->setProperty("WorkspaceIndexList", wsIndices);
224 group->execute();
225 return group->getProperty("OutputWorkspace");
226}
227
235double GetEiMonDet3::minimumTOF(const API::MatrixWorkspace &ws, const double sampleToDetectorDistance) {
236 const double maxEnergy = getProperty(Prop::MAX_ENERGY);
237 const auto &spectrumInfo = ws.spectrumInfo();
238 return Kernel::UnitConversion::run("Energy", "TOF", maxEnergy, spectrumInfo.l1(), sampleToDetectorDistance, 0.,
240}
241
247double GetEiMonDet3::monitorPeakPosition(const size_t monitorIndex) {
248 API::ITableWorkspace_sptr monitorEPPWs = getProperty(Prop::MONITOR_EPP_TABLE);
249 const auto &status = monitorEPPWs->getRef<std::string>("FitStatus", monitorIndex);
250 if (status != "success" && status != "narrowPeak") {
251 throw std::runtime_error("Monitor EPP fit status shows a failure.");
252 }
253 return monitorEPPWs->getRef<double>("PeakCentre", monitorIndex);
254}
255
262 auto findEPP = createChildAlgorithm("FindEPP");
263 findEPP->setProperty("InputWorkspace", ws);
264 findEPP->setProperty("OutputWorkspace", "unused");
265 findEPP->execute();
266 API::ITableWorkspace_sptr eppTable = findEPP->getProperty("OutputWorkspace");
267 const auto &status = eppTable->getRef<std::string>("FitStatus", 0);
268 if (status != "success" && status != "narrowPeak") {
269 throw std::runtime_error("Could not fit a Gaussian to the data.");
270 }
271 return eppTable->getRef<double>("PeakCentre", 0);
272}
273
274} // namespace Mantid::Algorithms
#define DECLARE_ALGORITHM(classname)
Definition: Algorithm.h:576
double energy
Definition: GetAllEi.cpp:157
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
virtual std::shared_ptr< Algorithm > createChildAlgorithm(const std::string &name, const double startProgress=-1., const double endProgress=-1., const bool enableLogging=true, const int &version=-1)
Create a Child Algorithm.
Definition: Algorithm.cpp:842
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
bool isDefault(const std::string &name) const
Definition: Algorithm.cpp:2084
Kernel::Logger m_log
Logger for this algorithm.
Definition: Algorithm.h:450
const SpectrumInfo & spectrumInfo() const
Return a reference to the SpectrumInfo object.
const Run & run() const
Run details object access.
A validator which checks that a workspace has a valid instrument.
bool hasProperty(const std::string &name) const
Does the property exist on the object.
Definition: LogManager.cpp:265
double getPropertyAsSingleValue(const std::string &name, Kernel::Math::StatisticType statistic=Kernel::Math::Mean) const
Returns a property as a single double value from its name.
Definition: LogManager.cpp:349
Base MatrixWorkspace Abstract Class.
A property class for workspaces.
A validator which checks that the unit of the workspace referred to by a WorkspaceProperty is the exp...
Estimates the incident neutron energy from the time of flight between a monitor and a set of detector...
Definition: GetEiMonDet3.h:23
double computeTOF(const API::MatrixWorkspace &detectorWs, const double detectorEPP, const double monitorEPP, const double minTOF)
Calculates the time of flight from the monitor to the detectors.
const std::string summary() const override
Returns a summary of algorithm's purpose.
const std::string category() const override
Algorithm's category for identification overriding a virtual method.
void init() override
Initialized the algorithm.
int version() const override
Returns algorithm's version for identification.
double peakPosition(API::MatrixWorkspace_sptr &ws)
Returns the TOF of the grouped detectors' elastic peak.
double monitorPeakPosition(const size_t monitorIndex)
Returns the TOF of the monitor's peak.
const std::vector< std::string > seeAlso() const override
Function to return all of the seeAlso algorithms related to this algorithm.
API::MatrixWorkspace_sptr groupSpectra(API::MatrixWorkspace_sptr &ws, const std::vector< size_t > &wsIndices)
Runs GroupDetectors on given workspace indices.
double minimumTOF(const API::MatrixWorkspace &ws, const double sampleToDetectorDistance)
Computes the minimum TOF between monitor and detectors from maximum energy.
void exec() override
Executes the algorithm.
IPropertyManager * setProperty(const std::string &name, const T &value)
Templated method to set the value of a PropertyWithValue.
void setPropertySettings(const std::string &name, std::unique_ptr< IPropertySettings > settings)
void notice(const std::string &msg)
Logs at notice level.
Definition: Logger.cpp:95
void warning(const std::string &msg)
Logs at warning level.
Definition: Logger.cpp:86
static double run(const std::string &src, const std::string &dest, const double srcValue, const double l1, const double l2, const double theta, const DeltaEMode::Type emode, const double efixed)
Convert a single value between the given units (as strings)
std::shared_ptr< ITableWorkspace > ITableWorkspace_sptr
shared pointer to Mantid::API::ITableWorkspace
std::shared_ptr< MatrixWorkspace > MatrixWorkspace_sptr
shared pointer to the matrix workspace base class
static const std::string PULSE_INTERVAL("pulse_interval")
Name of the pulse interval sample log.
A namespace containing physical constants that are required by algorithms and unit routines.
Definition: Atom.h:14
static constexpr double NeutronMass
Mass of the neutron in kg.
static constexpr double meV
1 meV in Joules.
constexpr int EMPTY_INT() noexcept
Returns what we consider an "empty" integer within a property.
Definition: EmptyValues.h:25
constexpr double EMPTY_DBL() noexcept
Returns what we consider an "empty" double within a property.
Definition: EmptyValues.h:43
String constants for algorithm's properties.
STL namespace.
@ Input
An input workspace.
Definition: Property.h:53
@ Output
An output workspace.
Definition: Property.h:54