19#include <boost/lexical_cast.hpp>
28using namespace Kernel;
30using namespace Geometry;
51 auto val = std::make_shared<CompositeValidator>();
56 "The X units of this workspace must be time of flight with times in\n"
58 auto mustBePositive = std::make_shared<BoundedValidator<int>>();
59 mustBePositive->setLower(0);
61 "The spectrum number of the output of the first monitor, e.g. MAPS\n"
62 "41474, MARI 2, MERLIN 69634");
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);
69 "An approximate value for the typical incident energy, energy of\n"
70 "neutrons leaving the source (meV)");
90 double dist2moni0 = -1, dist2moni1 = -1;
91 getGeometry(inWS, mon1Spec, mon2Spec, dist2moni0, dist2moni1);
99 const double peakLoc0 = 1e6 *
timeToFly(dist2moni0, E_est);
103 "will be searched for at TOF "
104 << peakLoc0 <<
" micro seconds +/-" << boost::lexical_cast<std::string>(100.0 *
HALF_WINDOW)
106 const double peakLoc1 = 1e6 *
timeToFly(dist2moni1, E_est);
108 "will be searched for at TOF "
109 << peakLoc1 <<
" micro seconds +/-" << boost::lexical_cast<std::string>(100.0 *
HALF_WINDOW)
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";
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";
126 double meanSpeed = (dist2moni1 - dist2moni0) / (1e-6 * (t_monitor1 - t_monitor0));
131 g_log.
notice() <<
"The incident energy has been calculated to be " << E_i <<
" meV"
132 <<
" (your estimate was " << E_est <<
" meV)\n";
137 inWS->mutableRun().addProperty(incident_energy,
true);
150 double &monitor0Dist,
double &monitor1Dist)
const {
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();
163 const auto &spectrumInfo = WS->spectrumInfo();
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();
172 monitor0Dist = spectrumInfo.l1() + spectrumInfo.l2(monWI);
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();
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();
189 monitor1Dist = spectrumInfo.l1() + spectrumInfo.l2(monWI);
209 std::vector<specnum_t> specNumTemp(&specNum1, &specNum1 + 1);
210 auto wsInds = WS->getIndicesFromSpectra(specNumTemp);
212 if (wsInds.size() != 1) {
214 g_log.
error() <<
"Couldn't find the first monitor "
221 specNumTemp[0] = specNum2;
222 auto wsIndexTemp = WS->getIndicesFromSpectra(specNumTemp);
223 if (wsIndexTemp.size() != 1) {
226 "monitor spectrum, number "
231 wsInds.emplace_back(wsIndexTemp[0]);
267 const auto ×Array = WS->x(monitIn);
270 double halfWin = (timesArray.back() - timesArray.front()) *
HALF_WINDOW;
271 if (monitIn < std::numeric_limits<int>::max()) {
272 auto ivsInd =
static_cast<int>(monitIn);
276 extractSpec(ivsInd, peakTime - halfWin, peakTime + halfWin);
279 " integer number defined for this OS."
280 " This behaviour is not yet supported");
290 int64_t centreGausInd;
291 double height, backGroundlev;
303 return (lHalf + rHalf) / 2;
323 childAlg->setProperty(
"XMin", start);
324 childAlg->setProperty(
"XMax", end);
325 childAlg->setProperty(
"StartWorkspaceIndex", wsInd);
326 childAlg->setProperty(
"EndWorkspaceIndex", wsInd);
327 childAlg->executeAsChildAlg();
329 m_tempWS = childAlg->getProperty(
"OutputWorkspace");
359 background = std::accumulate(
Y.begin(),
Y.end(), 0.0);
362 auto maxHeight = std::max_element(
Y.begin(),
Y.end());
364 centreInd = std::distance(
Y.begin(), maxHeight);
366 background = background /
static_cast<double>(
Y.size());
368 throw std::invalid_argument(
"No peak was found or its height is less than the threshold of " +
370 " times the mean background, was the energy estimate (" +
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";
393 auto endInd = startInd;
397 while (
Y[endInd] > (
height + noise) / 2.0) {
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 "
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 "
413 if (std::abs(
static_cast<int64_t
>(endInd - startInd)) <
PEAK_THRESH_W) {
414 g_log.
error() <<
"Likely precision problem or error, one half height "
415 "distance is less than the threshold number of bins from "
417 << std::abs(
static_cast<int>(endInd - startInd)) <<
"<" <<
PEAK_THRESH_W
418 <<
". Check the monitor peak\n";
422 double hOverN = (
height - noise) / noise;
424 std::abs(hOverN *
static_cast<double>(endInd - startInd)) <
PEAK_THRESH_A) {
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 (>" +
429 " bins). Was the energy estimate close enough?");
433 double halfTime = (
X[endInd] +
X[endInd + 1]) / 2;
436 if (endInd != startInd) {
441 double gradient = (
Y[endInd] -
Y[endInd - go]) / (
X[endInd] -
X[endInd - go]);
444 double deltaY =
Y[endInd] - (
height + noise) / 2.0;
447 halfTime -= deltaY / gradient;
450 g_log.
debug() <<
"One half height point found at TOF = " << halfTime <<
" microseconds\n";
#define DECLARE_ALGORITHM(classname)
Base class from which all concrete algorithm classes should be derived.
void declareProperty(std::unique_ptr< Kernel::Property > p, const std::string &doc="") override
Add a property to the list of managed properties.
std::string getPropertyValue(const std::string &name) const override
Get the value of a property as a string.
TypedValue getProperty(const std::string &name) const override
Get the value of a property.
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.
void progress(double p, const std::string &msg="", double estimatedTime=0.0, int progressPrecision=0)
Sends ProgressNotification.
void interruption_point()
This is called during long-running operations, and check if the algorithm has requested that it be ca...
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
static const double HALF_WINDOW
the range of TOF X-values over which the peak will be searched is double this value,...
direction
used by the function findHalfLoc to indicate whether to search left or right
@ GO_RIGHT
flag value to search right
@ GO_LEFT
flag value to serch left
double m_fracCompl
An estimate of the percentage of the algorithm runtimes that has been completed.
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
void exec() override
Executes the algorithm.
static const double CROP
fraction of algorithm time taken up with running CropWorkspace
API::MatrixWorkspace_sptr m_tempWS
name of the tempory workspace that we create and use
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...
void init() override
Virtual method - must be overridden by concrete algorithm.
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...
double neutron_E_At(double speed) const
Get the kinetic energy of a neuton in joules given it speed using E=mv^2/2.
void getPeakEstimates(double &height, int64_t ¢reInd, double &background) const
Finds the largest peak by looping through the histogram and finding the maximum value.
static const double FIT_PEAK
single call to ConvertToDistribution
std::vector< size_t > getMonitorWsIndexs(const API::MatrixWorkspace_const_sptr &WS, specnum_t specNum1, specnum_t specNum2) const
Converts detector IDs to spectra indexes.
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.
GetEi()
Empty default constructor algorithm() calls the constructor in the base class.
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 ...
void extractSpec(int wsInd, double start, double end)
Calls CropWorkspace as a Child Algorithm and passes to it the InputWorkspace property.
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,...
void advanceProgress(double toAdd)
Update the percentage complete estimate assuming that the algorithm has completed a task with estimat...
static const double GET_COUNT_RATE
fraction of algorithm taken by a
Exception for when an item is not found in a collection.
Marks code as not implemented yet.
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.
void notice(const std::string &msg)
Logs at notice level.
void error(const std::string &msg)
Logs at error level.
void information(const std::string &msg)
Logs at information level.
The concrete, templated class for properties.
Base class for properties.
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.
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.
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.
@ Output
An output workspace.