Mantid
Loading...
Searching...
No Matches
GetEi.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"
18
19#include <boost/lexical_cast.hpp>
20#include <cmath>
21#include <numeric>
22
23namespace Mantid::Algorithms {
24
25// Register the algorithm into the algorithm factory
27
28using namespace Kernel;
29using namespace API;
30using namespace Geometry;
31
32// adjustable fit criteria, increase the first number or reduce any of the last
33// three for more promiscuous peak fitting
34// from the estimated location of the peak search forward by the following
35// fraction and backward by the same fraction
36const double GetEi::HALF_WINDOW = 8.0 / 100;
37const double GetEi::PEAK_THRESH_H = 3.0;
38const double GetEi::PEAK_THRESH_A = 5.0;
39const int64_t GetEi::PEAK_THRESH_W = 3;
40
41// progress estimates
42const double GetEi::CROP = 0.15;
43const double GetEi::GET_COUNT_RATE = 0.15;
44const double GetEi::FIT_PEAK = 0.2;
45
47GetEi::GetEi() : Algorithm(), m_tempWS(), m_fracCompl(0.0) {}
48
50 // Declare required input parameters for algorithm and do some validation here
51 auto val = std::make_shared<CompositeValidator>();
52 val->add<WorkspaceUnitValidator>("TOF");
53 val->add<HistogramValidator>();
54 val->add<InstrumentValidator>();
55 declareProperty(std::make_unique<WorkspaceProperty<>>("InputWorkspace", "", Direction::Input, val),
56 "The X units of this workspace must be time of flight with times in\n"
57 "micro-seconds");
58 auto mustBePositive = std::make_shared<BoundedValidator<int>>();
59 mustBePositive->setLower(0);
60 declareProperty("Monitor1Spec", -1, mustBePositive,
61 "The spectrum number of the output of the first monitor, e.g. MAPS\n"
62 "41474, MARI 2, MERLIN 69634");
63 declareProperty("Monitor2Spec", -1, mustBePositive,
64 "The spectrum number of the output of the second monitor e.g. MAPS\n"
65 "41475, MARI 3, MERLIN 69638");
66 auto positiveDouble = std::make_shared<BoundedValidator<double>>();
67 positiveDouble->setLower(0);
68 declareProperty("EnergyEstimate", -1.0, positiveDouble,
69 "An approximate value for the typical incident energy, energy of\n"
70 "neutrons leaving the source (meV)");
71 declareProperty("IncidentEnergy", -1.0, Direction::Output);
72 declareProperty("FirstMonitorPeak", -1.0, Direction::Output);
73
74 m_fracCompl = 0.0;
75}
76
87 MatrixWorkspace_sptr inWS = getProperty("InputWorkspace");
88 const specnum_t mon1Spec = getProperty("Monitor1Spec");
89 const specnum_t mon2Spec = getProperty("Monitor2Spec");
90 double dist2moni0 = -1, dist2moni1 = -1;
91 getGeometry(inWS, mon1Spec, mon2Spec, dist2moni0, dist2moni1);
92
93 // the E_i estimate is used to find (identify) the monitor peaks, checking
94 // prior to fitting will throw an exception if this estimate is too big or
95 // small
96 const double E_est = getProperty("EnergyEstimate");
97 // we're assuming that the time units for the X-values in the workspace are
98 // micro-seconds
99 const double peakLoc0 = 1e6 * timeToFly(dist2moni0, E_est);
100 // write a lot of stuff to the log at user level as it is very possible for
101 // fit routines not to the expected thing
102 g_log.information() << "Based on the user selected energy the first peak "
103 "will be searched for at TOF "
104 << peakLoc0 << " micro seconds +/-" << boost::lexical_cast<std::string>(100.0 * HALF_WINDOW)
105 << "%\n";
106 const double peakLoc1 = 1e6 * timeToFly(dist2moni1, E_est);
107 g_log.information() << "Based on the user selected energy the second peak "
108 "will be searched for at TOF "
109 << peakLoc1 << " micro seconds +/-" << boost::lexical_cast<std::string>(100.0 * HALF_WINDOW)
110 << "%\n";
111
112 // get the histograms created by the monitors
113 std::vector<size_t> indexes = getMonitorWsIndexs(inWS, mon1Spec, mon2Spec);
114
115 g_log.information() << "Looking for a peak in the first monitor spectrum, spectra index " << indexes[0] << '\n';
116 double t_monitor0 = getPeakCentre(inWS, indexes[0], peakLoc0);
117 g_log.notice() << "The first peak has been found at TOF = " << t_monitor0 << " microseconds\n";
118 setProperty("FirstMonitorPeak", t_monitor0);
119
120 g_log.information() << "Looking for a peak in the second monitor spectrum, spectra index " << indexes[1] << '\n';
121 double t_monitor1 = getPeakCentre(inWS, indexes[1], peakLoc1);
122 g_log.information() << "The second peak has been found at TOF = " << t_monitor1 << " microseconds\n";
123
124 // assumes that the source and the both mintors lie on one straight line, the
125 // 1e-6 converts microseconds to seconds as the mean speed needs to be in m/s
126 double meanSpeed = (dist2moni1 - dist2moni0) / (1e-6 * (t_monitor1 - t_monitor0));
127
128 // uses 0.5mv^2 to get the kinetic energy in joules which we then convert to
129 // meV
130 double E_i = neutron_E_At(meanSpeed) / PhysicalConstants::meV;
131 g_log.notice() << "The incident energy has been calculated to be " << E_i << " meV"
132 << " (your estimate was " << E_est << " meV)\n";
133
134 setProperty("IncidentEnergy", E_i);
135 // store property in input workspace
136 Property *incident_energy = new PropertyWithValue<double>("Ei", E_i, Direction::Input);
137 inWS->mutableRun().addProperty(incident_energy, true);
138}
150 double &monitor0Dist, double &monitor1Dist) const {
151 const IComponent_const_sptr source = WS->getInstrument()->getSource();
152
153 // retrieve a pointer to the first detector and get its distance
154 size_t monWI = 0;
155 try {
156 monWI = WS->getIndexFromSpectrumNumber(mon0Spec);
157 } catch (std::runtime_error &) {
158 g_log.error() << "Could not find the workspace index for the monitor at spectrum " << mon0Spec << "\n";
159 g_log.error() << "Error retrieving data for the first monitor\n";
160 throw std::bad_cast();
161 }
162
163 const auto &spectrumInfo = WS->spectrumInfo();
164
165 if (!spectrumInfo.hasUniqueDetector(monWI)) {
166 g_log.error() << "The detector for spectrum number " << mon0Spec
167 << " was either not found or is a group, grouped monitors "
168 "are not supported by this algorithm\n";
169 g_log.error() << "Error retrieving data for the first monitor\n";
170 throw std::bad_cast();
171 }
172 monitor0Dist = spectrumInfo.l1() + spectrumInfo.l2(monWI);
173
174 // repeat for the second detector
175 try {
176 monWI = WS->getIndexFromSpectrumNumber(mon1Spec);
177 } catch (std::runtime_error &) {
178 g_log.error() << "Could not find the workspace index for the monitor at spectrum " << mon0Spec << "\n";
179 g_log.error() << "Error retrieving data for the second monitor\n";
180 throw std::bad_cast();
181 }
182 if (!spectrumInfo.hasUniqueDetector(monWI)) {
183 g_log.error() << "The detector for spectrum number " << mon1Spec
184 << " was either not found or is a group, grouped monitors "
185 "are not supported by this algorithm\n";
186 g_log.error() << "Error retrieving data for the second monitor\n";
187 throw std::bad_cast();
188 }
189 monitor1Dist = spectrumInfo.l1() + spectrumInfo.l2(monWI);
190}
200std::vector<size_t>
202 specnum_t specNum2) const { // getting spectra numbers from detector IDs is
203 // hard because the map works the other way,
204 // getting index numbers from spectra numbers has
205 // the same problem and we are about to do both
206
207 // get the index number of the histogram for the first monitor
208
209 std::vector<specnum_t> specNumTemp(&specNum1, &specNum1 + 1);
210 auto wsInds = WS->getIndicesFromSpectra(specNumTemp);
211
212 if (wsInds.size() != 1) { // the monitor spectrum isn't present in the
213 // workspace, we can't continue from here
214 g_log.error() << "Couldn't find the first monitor "
215 "spectrum, number "
216 << specNum1 << '\n';
217 throw Exception::NotFoundError("GetEi::getMonitorWsIndexs()", specNum1);
218 }
219
220 // nowe the second monitor
221 specNumTemp[0] = specNum2;
222 auto wsIndexTemp = WS->getIndicesFromSpectra(specNumTemp);
223 if (wsIndexTemp.size() != 1) { // the monitor spectrum isn't present in the
224 // workspace, we can't continue from here
225 g_log.error() << "Couldn't find the second "
226 "monitor spectrum, number "
227 << specNum2 << '\n';
228 throw Exception::NotFoundError("GetEi::getMonitorWsIndexs()", specNum2);
229 }
230
231 wsInds.emplace_back(wsIndexTemp[0]);
232 return wsInds;
233}
240double GetEi::timeToFly(double s, double E_KE) const {
241 // E_KE = mv^2/2, s = vt
242 // t = s/v, v = sqrt(2*E_KE/m)
243 // t = s/sqrt(2*E_KE/m)
244
245 // convert E_KE to joules kg m^2 s^-2
247
248 return s / sqrt(2 * E_KE / PhysicalConstants::NeutronMass);
249}
250
266double GetEi::getPeakCentre(const API::MatrixWorkspace_const_sptr &WS, const size_t monitIn, const double peakTime) {
267 const auto &timesArray = WS->x(monitIn);
268 // we search for the peak only inside some window because there are often more
269 // peaks in the monitor histogram
270 double halfWin = (timesArray.back() - timesArray.front()) * HALF_WINDOW;
271 if (monitIn < std::numeric_limits<int>::max()) {
272 auto ivsInd = static_cast<int>(monitIn);
273
274 // runs CropWorkspace as a Child Algorithm to and puts the result in a new
275 // temporary workspace that will be deleted when this algorithm has finished
276 extractSpec(ivsInd, peakTime - halfWin, peakTime + halfWin);
277 } else {
278 throw Kernel::Exception::NotImplementedError("Spectra number exceeds maximal"
279 " integer number defined for this OS."
280 " This behaviour is not yet supported");
281 }
282 // converting the workspace to count rate is required by the fitting algorithm
283 // if the bin widths are not all the same
285 // look out for user cancel messgages as the above command can take a bit of
286 // time
288
289 // to store fit results
290 int64_t centreGausInd;
291 double height, backGroundlev;
292 getPeakEstimates(height, centreGausInd, backGroundlev);
293 // look out for user cancel messgages
295
296 // the peak centre is defined as the centre of the two half maximum points as
297 // this is better for asymmetric peaks
298 // first loop backwards along the histogram to get the first half height point
299 const double lHalf = findHalfLoc(centreGausInd, height, backGroundlev, GO_LEFT);
300 // go forewards to get the half height on the otherside of the peak
301 const double rHalf = findHalfLoc(centreGausInd, height, backGroundlev, GO_RIGHT);
302 // the peak centre is defined as the mean of the two half height times
303 return (lHalf + rHalf) / 2;
304}
317void GetEi::extractSpec(int wsInd, double start, double end) {
318 auto childAlg = createChildAlgorithm("CropWorkspace", 100 * m_fracCompl, 100 * (m_fracCompl + CROP));
319 m_fracCompl += CROP;
320
321 childAlg->setProperty<MatrixWorkspace_sptr>("InputWorkspace", getProperty("InputWorkspace"));
322
323 childAlg->setProperty("XMin", start);
324 childAlg->setProperty("XMax", end);
325 childAlg->setProperty("StartWorkspaceIndex", wsInd);
326 childAlg->setProperty("EndWorkspaceIndex", wsInd);
327 childAlg->executeAsChildAlg();
328
329 m_tempWS = childAlg->getProperty("OutputWorkspace");
330
331 // DEBUGGING CODE uncomment out the line below if you want to see the TOF
332 // window that was analysed
333 // AnalysisDataService::Instance().addOrReplace("croped_dist_del", m_tempWS);
336}
337
348void GetEi::getPeakEstimates(double &height, int64_t &centreInd, double &background) const {
349
350 const auto &X = m_tempWS->x(0);
351 const auto &Y = m_tempWS->y(0);
352 // take note of the number of background counts as error checking, do we have
353 // a peak or just a bump in the background
354 background = 0;
355 // start at the first Y value
356 height = Y[0];
357 centreInd = 0;
358
359 background = std::accumulate(Y.begin(), Y.end(), 0.0);
360 // then loop through all the Y values and find the tallest peak
361 // Todo use std::max to find max element and record index?
362 auto maxHeight = std::max_element(Y.begin(), Y.end());
363 height = *maxHeight;
364 centreInd = std::distance(Y.begin(), maxHeight);
365
366 background = background / static_cast<double>(Y.size());
367 if (height < PEAK_THRESH_H * background) {
368 throw std::invalid_argument("No peak was found or its height is less than the threshold of " +
369 boost::lexical_cast<std::string>(PEAK_THRESH_H) +
370 " times the mean background, was the energy estimate (" +
371 getPropertyValue("EnergyEstimate") + " meV) close enough?");
372 }
373
374 g_log.debug() << "Peak position is the bin that has the maximum Y value in "
375 "the monitor spectrum, which is at TOF "
376 << (X[centreInd] + X[centreInd + 1]) / 2 << " (peak height " << height << " counts/microsecond)\n";
377}
392double GetEi::findHalfLoc(size_t startInd, const double height, const double noise, const direction go) const {
393 auto endInd = startInd;
394
395 const auto &X = m_tempWS->x(0);
396 const auto &Y = m_tempWS->y(0);
397 while (Y[endInd] > (height + noise) / 2.0) {
398 endInd += go;
399 if (endInd < 1) {
400 throw std::out_of_range("Can't analyse peak, some of the peak is outside the " +
401 boost::lexical_cast<std::string>(HALF_WINDOW * 100) +
402 "% window, at TOF values that are too low. Was the energy estimate "
403 "close enough?");
404 }
405 if (endInd > Y.size() - 2) {
406 throw std::out_of_range("Can't analyse peak, some of the peak is outside the " +
407 boost::lexical_cast<std::string>(HALF_WINDOW * 100) +
408 "% window, at TOF values that are too high. Was the energy estimate "
409 "close enough?");
410 }
411 }
412
413 if (std::abs(static_cast<int64_t>(endInd - startInd)) < PEAK_THRESH_W) { // we didn't find a significant peak
414 g_log.error() << "Likely precision problem or error, one half height "
415 "distance is less than the threshold number of bins from "
416 "the central peak: "
417 << std::abs(static_cast<int>(endInd - startInd)) << "<" << PEAK_THRESH_W
418 << ". Check the monitor peak\n";
419 }
420 // we have a peak in range, do an area check to see if the peak has any
421 // significance
422 double hOverN = (height - noise) / noise;
423 if (hOverN < PEAK_THRESH_A &&
424 std::abs(hOverN * static_cast<double>(endInd - startInd)) < PEAK_THRESH_A) { // the peak could just be noise on
425 // the background, ignore it
426 throw std::invalid_argument("No good peak was found. The ratio of the height to the background "
427 "multiplied either half widths must be above the threshold (>" +
428 boost::lexical_cast<std::string>(PEAK_THRESH_A) +
429 " bins). Was the energy estimate close enough?");
430 }
431 // get the TOF value in the middle of the bin with the first value below the
432 // half height
433 double halfTime = (X[endInd] + X[endInd + 1]) / 2;
434 // interpolate back between the first bin with less than half the counts to
435 // the bin before it
436 if (endInd != startInd) { // let the bin that we found have coordinates (x_1,
437 // y_1) the distance of the half point (x_2, y_2)
438 // from this is (y_1-y_2)/gradient. Gradient =
439 // (y_3-y_1)/(x_3-x_1) where (x_3, y_3) are the
440 // coordinates of the other bin we are using
441 double gradient = (Y[endInd] - Y[endInd - go]) / (X[endInd] - X[endInd - go]);
442 // we don't need to check for a zero or negative gradient if we assume the
443 // endInd bin was found when the Y-value dropped below the threshold
444 double deltaY = Y[endInd] - (height + noise) / 2.0;
445 // correct for the interpolation back in the direction towards the peak
446 // centre
447 halfTime -= deltaY / gradient;
448 }
449
450 g_log.debug() << "One half height point found at TOF = " << halfTime << " microseconds\n";
451 return halfTime;
452}
457double GetEi::neutron_E_At(double speed) const {
458 // E_KE = mv^2/2
459 return PhysicalConstants::NeutronMass * speed * speed / (2);
460}
461
464void GetEi::advanceProgress(double toAdd) {
465 m_fracCompl += toAdd;
467 // look out for user cancel messgages
469}
470
471} // namespace Mantid::Algorithms
#define DECLARE_ALGORITHM(classname)
Definition: Algorithm.h:576
double height
Definition: GetAllEi.cpp:155
Base class from which all concrete algorithm classes should be derived.
Definition: Algorithm.h:85
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
std::string getPropertyValue(const std::string &name) const override
Get the value of a property as a string.
Definition: Algorithm.cpp:2026
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
void interruption_point()
This is called during long-running operations, and check if the algorithm has requested that it be ca...
Definition: Algorithm.cpp:1687
A validator which checks that a workspace contains histogram data (the default) or point data as requ...
A validator which checks that a workspace has a valid instrument.
A property class for workspaces.
A validator which checks that the unit of the workspace referred to by a WorkspaceProperty is the exp...
static const double PEAK_THRESH_H
ignore an peaks that are less than this factor of the background
Definition: GetEi.h:88
static const double HALF_WINDOW
the range of TOF X-values over which the peak will be searched is double this value,...
Definition: GetEi.h:86
direction
used by the function findHalfLoc to indicate whether to search left or right
Definition: GetEi.h:62
@ GO_RIGHT
flag value to search right
Definition: GetEi.h:64
@ GO_LEFT
flag value to serch left
Definition: GetEi.h:63
double m_fracCompl
An estimate of the percentage of the algorithm runtimes that has been completed.
Definition: GetEi.h:59
static const double PEAK_THRESH_A
ignore peaks where the half width times the ratio of the peak height to the background is less this
Definition: GetEi.h:91
void exec() override
Executes the algorithm.
Definition: GetEi.cpp:86
static const double CROP
fraction of algorithm time taken up with running CropWorkspace
Definition: GetEi.h:97
API::MatrixWorkspace_sptr m_tempWS
name of the tempory workspace that we create and use
Definition: GetEi.h:56
double getPeakCentre(const API::MatrixWorkspace_const_sptr &WS, const size_t monitIn, const double peakTime)
Looks for and examines a peak close to that specified by the input parameters and examines it to find...
Definition: GetEi.cpp:266
void init() override
Virtual method - must be overridden by concrete algorithm.
Definition: GetEi.cpp:49
double findHalfLoc(MantidVec::size_type startInd, const double height, const double noise, const direction go) const
Estimates the closest time, looking either or back, when the number of counts is half that in the bin...
Definition: GetEi.cpp:392
double neutron_E_At(double speed) const
Get the kinetic energy of a neuton in joules given it speed using E=mv^2/2.
Definition: GetEi.cpp:457
void getPeakEstimates(double &height, int64_t &centreInd, double &background) const
Finds the largest peak by looping through the histogram and finding the maximum value.
Definition: GetEi.cpp:348
static const double FIT_PEAK
single call to ConvertToDistribution
Definition: GetEi.h:100
std::vector< size_t > getMonitorWsIndexs(const API::MatrixWorkspace_const_sptr &WS, specnum_t specNum1, specnum_t specNum2) const
Converts detector IDs to spectra indexes.
Definition: GetEi.cpp:201
void getGeometry(const API::MatrixWorkspace_const_sptr &WS, specnum_t mon0Spec, specnum_t mon1Spec, double &monitor0Dist, double &monitor1Dist) const
Gets the distances between the source and detectors whose IDs you pass to it.
Definition: GetEi.cpp:149
GetEi()
Empty default constructor algorithm() calls the constructor in the base class.
Definition: GetEi.cpp:47
static const int64_t PEAK_THRESH_W
for peaks where the distance to the half heigth is less than this number of bins in either direction ...
Definition: GetEi.h:94
void extractSpec(int wsInd, double start, double end)
Calls CropWorkspace as a Child Algorithm and passes to it the InputWorkspace property.
Definition: GetEi.cpp:317
double timeToFly(double s, double E_KE) const
Uses E_KE = mv^2/2 and s = vt to calculate the time required for a neutron to travel a distance,...
Definition: GetEi.cpp:240
void advanceProgress(double toAdd)
Update the percentage complete estimate assuming that the algorithm has completed a task with estimat...
Definition: GetEi.cpp:464
static const double GET_COUNT_RATE
fraction of algorithm taken by a
Definition: GetEi.h:98
Exception for when an item is not found in a collection.
Definition: Exception.h:145
Marks code as not implemented yet.
Definition: Exception.h:138
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 notice(const std::string &msg)
Logs at notice level.
Definition: Logger.cpp:95
void error(const std::string &msg)
Logs at error level.
Definition: Logger.cpp:77
void information(const std::string &msg)
Logs at information level.
Definition: Logger.cpp:105
The concrete, templated class for properties.
Base class for properties.
Definition: Property.h:94
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
std::shared_ptr< const IComponent > IComponent_const_sptr
Typdef of a shared pointer to a const IComponent.
Definition: IComponent.h:161
static constexpr double NeutronMass
Mass of the neutron in kg.
static constexpr double meV
1 meV in Joules.
int32_t specnum_t
Typedef for a spectrum Number.
Definition: IDTypes.h:16
static void makeDistribution(const MatrixWorkspace_sptr &workspace, const bool forwards=true)
Divides the data in a workspace by the bin width to make it a distribution.
@ Input
An input workspace.
Definition: Property.h:53
@ Output
An output workspace.
Definition: Property.h:54