GetEi.cpp 21.7 KB
Newer Older
1
2
3
// Mantid Repository : https://github.com/mantidproject/mantid
//
// Copyright © 2018 ISIS Rutherford Appleton Laboratory UKRI,
4
5
//   NScD Oak Ridge National Laboratory, European Spallation Source,
//   Institut Laue - Langevin & CSNS, Institute of High Energy Physics, CAS
6
// SPDX - License - Identifier: GPL - 3.0 +
7
#include "MantidAlgorithms/GetEi.h"
8
9
#include "MantidAPI/HistogramValidator.h"
#include "MantidAPI/InstrumentValidator.h"
10
#include "MantidAPI/Run.h"
11
#include "MantidAPI/SpectrumInfo.h"
12
#include "MantidAPI/WorkspaceUnitValidator.h"
13
#include "MantidGeometry/Instrument.h"
14
#include "MantidKernel/BoundedValidator.h"
LamarMoore's avatar
LamarMoore committed
15
16
#include "MantidKernel/CompositeValidator.h"
#include "MantidKernel/Exception.h"
17
#include "MantidKernel/PhysicalConstants.h"
18

19
20
#include <boost/lexical_cast.hpp>
#include <cmath>
21
#include <numeric>
22

23
24
namespace Mantid {
namespace Algorithms {
25
26
27
28
29
30
31
32

// Register the algorithm into the algorithm factory
DECLARE_ALGORITHM(GetEi)

using namespace Kernel;
using namespace API;
using namespace Geometry;

33
34
35
36
37
// adjustable fit criteria, increase the first number or reduce any of the last
// three for more promiscuous peak fitting
// from the estimated location of the peak search forward by the following
// fraction and backward by the same fraction
const double GetEi::HALF_WINDOW = 8.0 / 100;
38
39
const double GetEi::PEAK_THRESH_H = 3.0;
const double GetEi::PEAK_THRESH_A = 5.0;
Peterson, Peter's avatar
Peterson, Peter committed
40
const int64_t GetEi::PEAK_THRESH_W = 3;
41
42
43
44
45

// progress estimates
const double GetEi::CROP = 0.15;
const double GetEi::GET_COUNT_RATE = 0.15;
const double GetEi::FIT_PEAK = 0.2;
46
47

/// Empty default constructor algorith() calls the constructor in the base class
48
GetEi::GetEi() : Algorithm(), m_tempWS(), m_fracCompl(0.0) {}
49

50
void GetEi::init() {
51
  // Declare required input parameters for algorithm and do some validation here
52
  auto val = std::make_shared<CompositeValidator>();
53
54
55
  val->add<WorkspaceUnitValidator>("TOF");
  val->add<HistogramValidator>();
  val->add<InstrumentValidator>();
Samuel Jones's avatar
Samuel Jones committed
56
57
58
  declareProperty(std::make_unique<WorkspaceProperty<>>("InputWorkspace", "", Direction::Input, val),
                  "The X units of this workspace must be time of flight with times in\n"
                  "micro-seconds");
59
  auto mustBePositive = std::make_shared<BoundedValidator<int>>();
60
  mustBePositive->setLower(0);
Samuel Jones's avatar
Samuel Jones committed
61
62
63
64
65
66
  declareProperty("Monitor1Spec", -1, mustBePositive,
                  "The spectrum number of the output of the first monitor, e.g. MAPS\n"
                  "41474, MARI 2, MERLIN 69634");
  declareProperty("Monitor2Spec", -1, mustBePositive,
                  "The spectrum number of the output of the second monitor e.g. MAPS\n"
                  "41475, MARI 3, MERLIN 69638");
67
  auto positiveDouble = std::make_shared<BoundedValidator<double>>();
68
  positiveDouble->setLower(0);
Samuel Jones's avatar
Samuel Jones committed
69
70
71
  declareProperty("EnergyEstimate", -1.0, positiveDouble,
                  "An approximate value for the typical incident energy, energy of\n"
                  "neutrons leaving the source (meV)");
72
  declareProperty("IncidentEnergy", -1.0, Direction::Output);
73
  declareProperty("FirstMonitorPeak", -1.0, Direction::Output);
74
75
76
77
78

  m_fracCompl = 0.0;
}

/** Executes the algorithm
LamarMoore's avatar
LamarMoore committed
79
80
81
82
83
84
85
86
 *  @throw out_of_range if the peak runs off the edge of the histogram
 *  @throw NotFoundError if one of the requested spectrum numbers was not found
 * in the workspace
 *  @throw IndexError if there is a problem converting spectra indexes to
 * spectra numbers, which would imply there is a problem with the workspace
 *  @throw invalid_argument if a good peak fit wasn't made or the input
 * workspace does not have common binning
 */
87
void GetEi::exec() {
88
  MatrixWorkspace_sptr inWS = getProperty("InputWorkspace");
89
90
  const specnum_t mon1Spec = getProperty("Monitor1Spec");
  const specnum_t mon2Spec = getProperty("Monitor2Spec");
91
  double dist2moni0 = -1, dist2moni1 = -1;
92
  getGeometry(inWS, mon1Spec, mon2Spec, dist2moni0, dist2moni1);
93

94
95
96
  // the E_i estimate is used to find (identify) the monitor peaks, checking
  // prior to fitting will throw an exception if this estimate is too big or
  // small
97
  const double E_est = getProperty("EnergyEstimate");
98
99
100
101
102
103
  // we're assuming that the time units for the X-values in the workspace are
  // micro-seconds
  const double peakLoc0 = 1e6 * timeToFly(dist2moni0, E_est);
  // write a lot of stuff to the log at user level as it is very possible for
  // fit routines not to the expected thing
  g_log.information() << "Based on the user selected energy the first peak "
LamarMoore's avatar
LamarMoore committed
104
                         "will be searched for at TOF "
Samuel Jones's avatar
Samuel Jones committed
105
                      << peakLoc0 << " micro seconds +/-" << boost::lexical_cast<std::string>(100.0 * HALF_WINDOW)
106
107
108
                      << "%\n";
  const double peakLoc1 = 1e6 * timeToFly(dist2moni1, E_est);
  g_log.information() << "Based on the user selected energy the second peak "
LamarMoore's avatar
LamarMoore committed
109
                         "will be searched for at TOF "
Samuel Jones's avatar
Samuel Jones committed
110
                      << peakLoc1 << " micro seconds +/-" << boost::lexical_cast<std::string>(100.0 * HALF_WINDOW)
111
112
113
                      << "%\n";

  // get the histograms created by the monitors
114
  std::vector<size_t> indexes = getMonitorWsIndexs(inWS, mon1Spec, mon2Spec);
115

Samuel Jones's avatar
Samuel Jones committed
116
  g_log.information() << "Looking for a peak in the first monitor spectrum, spectra index " << indexes[0] << '\n';
117
  double t_monitor0 = getPeakCentre(inWS, indexes[0], peakLoc0);
Samuel Jones's avatar
Samuel Jones committed
118
  g_log.notice() << "The first peak has been found at TOF = " << t_monitor0 << " microseconds\n";
119
120
  setProperty("FirstMonitorPeak", t_monitor0);

Samuel Jones's avatar
Samuel Jones committed
121
  g_log.information() << "Looking for a peak in the second monitor spectrum, spectra index " << indexes[1] << '\n';
122
  double t_monitor1 = getPeakCentre(inWS, indexes[1], peakLoc1);
Samuel Jones's avatar
Samuel Jones committed
123
  g_log.information() << "The second peak has been found at TOF = " << t_monitor1 << " microseconds\n";
124

125
126
  // assumes that the source and the both mintors lie on one straight line, the
  // 1e-6 converts microseconds to seconds as the mean speed needs to be in m/s
Samuel Jones's avatar
Samuel Jones committed
127
  double meanSpeed = (dist2moni1 - dist2moni0) / (1e-6 * (t_monitor1 - t_monitor0));
128

129
130
131
  // uses 0.5mv^2 to get the kinetic energy in joules which we then convert to
  // meV
  double E_i = neutron_E_At(meanSpeed) / PhysicalConstants::meV;
Samuel Jones's avatar
Samuel Jones committed
132
  g_log.notice() << "The incident energy has been calculated to be " << E_i << " meV"
133
                 << " (your estimate was " << E_est << " meV)\n";
134
135

  setProperty("IncidentEnergy", E_i);
136
  // store property in input workspace
Samuel Jones's avatar
Samuel Jones committed
137
  Property *incident_energy = new PropertyWithValue<double>("Ei", E_i, Direction::Input);
138
  inWS->mutableRun().addProperty(incident_energy, true);
139
140
}
/** Gets the distances between the source and detectors whose IDs you pass to it
LamarMoore's avatar
LamarMoore committed
141
142
143
144
145
146
147
148
149
 *  @param WS :: the input workspace
 *  @param mon0Spec :: Spectrum number of the output from the first monitor
 *  @param mon1Spec :: Spectrum number of the output from the second monitor
 *  @param monitor0Dist :: the calculated distance to the detector whose ID was
 * passed to this function first
 *  @param monitor1Dist :: calculated distance to the detector whose ID was
 * passed to this function second
 *  @throw NotFoundError if no detector is found for the detector ID given
 */
Samuel Jones's avatar
Samuel Jones committed
150
void GetEi::getGeometry(const API::MatrixWorkspace_const_sptr &WS, specnum_t mon0Spec, specnum_t mon1Spec,
David Fairbrother's avatar
David Fairbrother committed
151
                        double &monitor0Dist, double &monitor1Dist) const {
152
  const IComponent_const_sptr source = WS->getInstrument()->getSource();
153
154

  // retrieve a pointer to the first detector and get its distance
155
  size_t monWI = 0;
156
  try {
157
    monWI = WS->getIndexFromSpectrumNumber(mon0Spec);
158
  } catch (std::runtime_error &) {
Samuel Jones's avatar
Samuel Jones committed
159
    g_log.error() << "Could not find the workspace index for the monitor at spectrum " << mon0Spec << "\n";
160
    g_log.error() << "Error retrieving data for the first monitor\n";
161
162
163
    throw std::bad_cast();
  }

164
165
166
  const auto &spectrumInfo = WS->spectrumInfo();

  if (!spectrumInfo.hasUniqueDetector(monWI)) {
167
168
169
    g_log.error() << "The detector for spectrum number " << mon0Spec
                  << " was either not found or is a group, grouped monitors "
                     "are not supported by this algorithm\n";
170
    g_log.error() << "Error retrieving data for the first monitor\n";
171
    throw std::bad_cast();
172
  }
173
  monitor0Dist = spectrumInfo.l1() + spectrumInfo.l2(monWI);
174

175
  // repeat for the second detector
176
  try {
Alex Buts's avatar
Alex Buts committed
177
    monWI = WS->getIndexFromSpectrumNumber(mon1Spec);
178
  } catch (std::runtime_error &) {
Samuel Jones's avatar
Samuel Jones committed
179
    g_log.error() << "Could not find the workspace index for the monitor at spectrum " << mon0Spec << "\n";
180
181
182
    g_log.error() << "Error retrieving data for the second monitor\n";
    throw std::bad_cast();
  }
183
  if (!spectrumInfo.hasUniqueDetector(monWI)) {
184
185
186
    g_log.error() << "The detector for spectrum number " << mon1Spec
                  << " was either not found or is a group, grouped monitors "
                     "are not supported by this algorithm\n";
187
188
    g_log.error() << "Error retrieving data for the second monitor\n";
    throw std::bad_cast();
189
  }
190
  monitor1Dist = spectrumInfo.l1() + spectrumInfo.l2(monWI);
191
}
192
/** Converts detector IDs to spectra indexes
LamarMoore's avatar
LamarMoore committed
193
194
195
196
197
198
199
200
 *  @param WS :: the workspace on which the calculations are being performed
 *  @param specNum1 :: spectrum number of the output of the first monitor
 *  @param specNum2 :: spectrum number of the output of the second monitor
 *  @return the indexes of the histograms created by the detector whose ID were
 * passed
 *  @throw NotFoundError if one of the requested spectrum numbers was not found
 * in the workspace
 */
Samuel Jones's avatar
Samuel Jones committed
201
202
203
204
std::vector<size_t>
GetEi::getMonitorWsIndexs(const API::MatrixWorkspace_const_sptr &WS, specnum_t specNum1,
                          specnum_t specNum2) const { // getting spectra numbers from detector IDs is
                                                      // hard because the map works the other way,
Nick Draper's avatar
Nick Draper committed
205
206
  // getting index numbers from spectra numbers has
  // the same problem and we are about to do both
207

208
  // get the index number of the histogram for the first monitor
209

210
  std::vector<specnum_t> specNumTemp(&specNum1, &specNum1 + 1);
211
212
213
  auto wsInds = WS->getIndicesFromSpectra(specNumTemp);

  if (wsInds.size() != 1) { // the monitor spectrum isn't present in the
Hahn, Steven's avatar
Hahn, Steven committed
214
215
    // workspace, we can't continue from here
    g_log.error() << "Couldn't find the first monitor "
LamarMoore's avatar
LamarMoore committed
216
217
                     "spectrum, number "
                  << specNum1 << '\n';
218
    throw Exception::NotFoundError("GetEi::getMonitorWsIndexs()", specNum1);
219
220
221
222
  }

  // nowe the second monitor
  specNumTemp[0] = specNum2;
223
224
  auto wsIndexTemp = WS->getIndicesFromSpectra(specNumTemp);
  if (wsIndexTemp.size() != 1) { // the monitor spectrum isn't present in the
Hahn, Steven's avatar
Hahn, Steven committed
225
226
    // workspace, we can't continue from here
    g_log.error() << "Couldn't find the second "
LamarMoore's avatar
LamarMoore committed
227
228
                     "monitor spectrum, number "
                  << specNum2 << '\n';
229
    throw Exception::NotFoundError("GetEi::getMonitorWsIndexs()", specNum2);
230
  }
231

232
  wsInds.emplace_back(wsIndexTemp[0]);
233
  return wsInds;
234
}
235
/** Uses E_KE = mv^2/2 and s = vt to calculate the time required for a neutron
LamarMoore's avatar
LamarMoore committed
236
237
238
239
240
 *  to travel a distance, s
 * @param s :: ditance travelled in meters
 * @param E_KE :: kinetic energy in meV
 * @return the time to taken to travel that uninterrupted distance in seconds
 */
241
double GetEi::timeToFly(double s, double E_KE) const {
242
243
244
245
246
247
248
  // E_KE = mv^2/2, s = vt
  // t = s/v, v = sqrt(2*E_KE/m)
  // t = s/sqrt(2*E_KE/m)

  // convert E_KE to joules kg m^2 s^-2
  E_KE *= PhysicalConstants::meV;

249
  return s / sqrt(2 * E_KE / PhysicalConstants::NeutronMass);
250
}
251

252
/** Looks for and examines a peak close to that specified by the input
LamarMoore's avatar
LamarMoore committed
253
254
255
256
257
258
259
260
261
262
263
264
265
266
 * parameters and
 *  examines it to find a representative time for when the neutrons hit the
 * detector
 *  @param WS :: the workspace containing the monitor spectrum
 *  @param monitIn :: the index of the histogram that contains the monitor
 * spectrum
 *  @param peakTime :: the estimated TOF of the monitor peak in the time units
 * of the workspace
 *  @return a time of flight value in the peak in microseconds
 *  @throw invalid_argument if a good peak fit wasn't made or the input
 * workspace does not have common binning
 *  @throw out_of_range if the peak runs off the edge of the histogram
 *  @throw runtime_error a Child Algorithm just falls over
 */
Samuel Jones's avatar
Samuel Jones committed
267
double GetEi::getPeakCentre(const API::MatrixWorkspace_const_sptr &WS, const size_t monitIn, const double peakTime) {
268
  const auto &timesArray = WS->x(monitIn);
269
270
271
  // we search for the peak only inside some window because there are often more
  // peaks in the monitor histogram
  double halfWin = (timesArray.back() - timesArray.front()) * HALF_WINDOW;
272
  if (monitIn < std::numeric_limits<int>::max()) {
273
    auto ivsInd = static_cast<int>(monitIn);
274

Alex Buts's avatar
Alex Buts committed
275
    // runs CropWorkspace as a Child Algorithm to and puts the result in a new
276
    // temporary workspace that will be deleted when this algorithm has finished
Alex Buts's avatar
Alex Buts committed
277
278
    extractSpec(ivsInd, peakTime - halfWin, peakTime + halfWin);
  } else {
Samuel Jones's avatar
Samuel Jones committed
279
280
281
    throw Kernel::Exception::NotImplementedError("Spectra number exceeds maximal"
                                                 " integer number defined for this OS."
                                                 " This behaviour is not yet supported");
282
  }
283
284
  // converting the workspace to count rate is required by the fitting algorithm
  // if the bin widths are not all the same
285
  WorkspaceHelpers::makeDistribution(m_tempWS);
286
287
  // look out for user cancel messgages as the above command can take a bit of
  // time
288
289
  advanceProgress(GET_COUNT_RATE);

290
  // to store fit results
Peterson, Peter's avatar
Peterson, Peter committed
291
  int64_t centreGausInd;
292
293
  double height, backGroundlev;
  getPeakEstimates(height, centreGausInd, backGroundlev);
294
295
  // look out for user cancel messgages
  advanceProgress(FIT_PEAK);
296

297
298
  // the peak centre is defined as the centre of the two half maximum points as
  // this is better for asymmetric peaks
299
  // first loop backwards along the histogram to get the first half height point
Samuel Jones's avatar
Samuel Jones committed
300
  const double lHalf = findHalfLoc(centreGausInd, height, backGroundlev, GO_LEFT);
301
  // go forewards to get the half height on the otherside of the peak
Samuel Jones's avatar
Samuel Jones committed
302
  const double rHalf = findHalfLoc(centreGausInd, height, backGroundlev, GO_RIGHT);
303
304
  // the peak centre is defined as the mean of the two half height times
  return (lHalf + rHalf) / 2;
305
}
306
/** Calls CropWorkspace as a Child Algorithm and passes to it the InputWorkspace
LamarMoore's avatar
LamarMoore committed
307
308
309
310
311
312
313
314
315
316
317
 * property
 *  @param wsInd :: the index number of the histogram to extract
 *  @param start :: the number of the first bin to include (starts counting bins
 * at 0)
 *  @param end :: the number of the last bin to include (starts counting bins at
 * 0)
 *  @throw out_of_range if start, end or wsInd are set outside of the valid
 * range for the workspace
 *  @throw runtime_error if the algorithm just falls over
 *  @throw invalid_argument if the input workspace does not have common binning
 */
318
void GetEi::extractSpec(int wsInd, double start, double end) {
319
  auto childAlg = createChildAlgorithm("CropWorkspace", 100 * m_fracCompl, 100 * (m_fracCompl + CROP));
320
  m_fracCompl += CROP;
321

Samuel Jones's avatar
Samuel Jones committed
322
  childAlg->setProperty<MatrixWorkspace_sptr>("InputWorkspace", getProperty("InputWorkspace"));
Alex Buts's avatar
Alex Buts committed
323

324
325
  childAlg->setProperty("XMin", start);
  childAlg->setProperty("XMax", end);
326
327
  childAlg->setProperty("StartWorkspaceIndex", wsInd);
  childAlg->setProperty("EndWorkspaceIndex", wsInd);
328
  childAlg->executeAsChildAlg();
329
330
331

  m_tempWS = childAlg->getProperty("OutputWorkspace");

332
333
334
  // DEBUGGING CODE uncomment out the line below if you want to see the TOF
  // window that was analysed
  // AnalysisDataService::Instance().addOrReplace("croped_dist_del", m_tempWS);
335
336
337
338
  progress(m_fracCompl);
  interruption_point();
}

339
/** Finds the largest peak by looping through the histogram and finding the
LamarMoore's avatar
LamarMoore committed
340
341
342
343
344
345
346
347
348
 * maximum
 *  value
 * @param height :: its passed value ignored it is set to the peak height
 * @param centreInd :: passed value is ignored it will be set to the bin index
 * of the peak center
 * @param background :: passed value ignored set mean number of counts per bin
 * in the spectrum
 * @throw invalid_argument if the peak is not clearly above the background
 */
Samuel Jones's avatar
Samuel Jones committed
349
void GetEi::getPeakEstimates(double &height, int64_t &centreInd, double &background) const {
350
351
352

  const auto &X = m_tempWS->x(0);
  const auto &Y = m_tempWS->y(0);
353
354
  // take note of the number of background counts as error checking, do we have
  // a peak or just a bump in the background
355
356
  background = 0;
  // start at the first Y value
357
  height = Y[0];
358
  centreInd = 0;
359

Dimitar Tasev's avatar
Dimitar Tasev committed
360
  background = std::accumulate(Y.begin(), Y.end(), 0.0);
361
  // then loop through all the Y values and find the tallest peak
Dimitar Tasev's avatar
Dimitar Tasev committed
362
363
364
365
  // Todo use std::max to find max element and record index?
  auto maxHeight = std::max_element(Y.begin(), Y.end());
  height = *maxHeight;
  centreInd = std::distance(Y.begin(), maxHeight);
Dimitar Tasev's avatar
Dimitar Tasev committed
366

367
  background = background / static_cast<double>(Y.size());
368
  if (height < PEAK_THRESH_H * background) {
Samuel Jones's avatar
Samuel Jones committed
369
370
371
372
    throw std::invalid_argument("No peak was found or its height is less than the threshold of " +
                                boost::lexical_cast<std::string>(PEAK_THRESH_H) +
                                " times the mean background, was the energy estimate (" +
                                getPropertyValue("EnergyEstimate") + " meV) close enough?");
373
  }
374

375
376
  g_log.debug() << "Peak position is the bin that has the maximum Y value in "
                   "the monitor spectrum, which is at TOF "
Samuel Jones's avatar
Samuel Jones committed
377
                << (X[centreInd] + X[centreInd + 1]) / 2 << " (peak height " << height << " counts/microsecond)\n";
378
}
379
/** Estimates the closest time, looking either or back, when the number of
LamarMoore's avatar
LamarMoore committed
380
381
382
383
384
385
386
387
388
389
390
391
392
 * counts is
 *  half that in the bin whose index that passed
 *  @param startInd :: index of the bin to search around, e.g. the index of the
 * peak centre
 *  @param height :: the number of counts (or count rate) to compare against
 * e.g. a peak height
 *  @param noise :: mean number of counts in each bin in the workspace
 *  @param go :: either GetEi::GO_LEFT or GetEi::GO_RIGHT
 *  @return estimated TOF of the half maximum point
 *  @throw out_of_range if the end of the histogram is reached before the point
 * is found
 *  @throw invalid_argument if the peak is too thin
 */
Samuel Jones's avatar
Samuel Jones committed
393
double GetEi::findHalfLoc(size_t startInd, const double height, const double noise, const direction go) const {
394
  auto endInd = startInd;
395

396
397
398
  const auto &X = m_tempWS->x(0);
  const auto &Y = m_tempWS->y(0);
  while (Y[endInd] > (height + noise) / 2.0) {
399
    endInd += go;
400
    if (endInd < 1) {
Samuel Jones's avatar
Samuel Jones committed
401
402
403
404
      throw std::out_of_range("Can't analyse peak, some of the peak is outside the " +
                              boost::lexical_cast<std::string>(HALF_WINDOW * 100) +
                              "% window, at TOF values that are too low. Was the energy estimate "
                              "close enough?");
405
    }
406
    if (endInd > Y.size() - 2) {
Samuel Jones's avatar
Samuel Jones committed
407
408
409
410
      throw std::out_of_range("Can't analyse peak, some of the peak is outside the " +
                              boost::lexical_cast<std::string>(HALF_WINDOW * 100) +
                              "% window, at TOF values that are too high. Was the energy estimate "
                              "close enough?");
411
412
    }
  }
413

Samuel Jones's avatar
Samuel Jones committed
414
  if (std::abs(static_cast<int64_t>(endInd - startInd)) < PEAK_THRESH_W) { // we didn't find a significant peak
415
416
417
    g_log.error() << "Likely precision problem or error, one half height "
                     "distance is less than the threshold number of bins from "
                     "the central peak: "
Samuel Jones's avatar
Samuel Jones committed
418
419
                  << std::abs(static_cast<int>(endInd - startInd)) << "<" << PEAK_THRESH_W
                  << ". Check the monitor peak\n";
420
  }
421
422
423
424
  // we have a peak in range, do an area check to see if the peak has any
  // significance
  double hOverN = (height - noise) / noise;
  if (hOverN < PEAK_THRESH_A &&
Samuel Jones's avatar
Samuel Jones committed
425
426
427
428
429
430
      std::abs(hOverN * static_cast<double>(endInd - startInd)) < PEAK_THRESH_A) { // the peak could just be noise on
                                                                                   // the background, ignore it
    throw std::invalid_argument("No good peak was found. The ratio of the height to the background "
                                "multiplied either half widths must be above the threshold (>" +
                                boost::lexical_cast<std::string>(PEAK_THRESH_A) +
                                " bins). Was the energy estimate close enough?");
431
  }
432
433
  // get the TOF value in the middle of the bin with the first value below the
  // half height
434
  double halfTime = (X[endInd] + X[endInd + 1]) / 2;
435
436
437
438
439
440
441
  // interpolate back between the first bin with less than half the counts to
  // the bin before it
  if (endInd != startInd) { // let the bin that we found have coordinates (x_1,
                            // y_1) the distance of the half point (x_2, y_2)
                            // from this is (y_1-y_2)/gradient. Gradient =
                            // (y_3-y_1)/(x_3-x_1) where (x_3, y_3) are the
                            // coordinates of the other bin we are using
Samuel Jones's avatar
Samuel Jones committed
442
    double gradient = (Y[endInd] - Y[endInd - go]) / (X[endInd] - X[endInd - go]);
443
444
    // we don't need to check for a zero or negative gradient if we assume the
    // endInd bin was found when the Y-value dropped below the threshold
445
    double deltaY = Y[endInd] - (height + noise) / 2.0;
446
447
448
    // correct for the interpolation back in the direction towards the peak
    // centre
    halfTime -= deltaY / gradient;
449
450
  }

Samuel Jones's avatar
Samuel Jones committed
451
  g_log.debug() << "One half height point found at TOF = " << halfTime << " microseconds\n";
452
  return halfTime;
453
}
454
/** Get the kinetic energy of a neuton in joules given it speed using E=mv^2/2
LamarMoore's avatar
LamarMoore committed
455
456
457
 *  @param speed :: the instantanious speed of a neutron in metres per second
 *  @return the energy in joules
 */
458
double GetEi::neutron_E_At(double speed) const {
459
  // E_KE = mv^2/2
460
  return PhysicalConstants::NeutronMass * speed * speed / (2);
461
462
}

463
464
465
/// Update the percentage complete estimate assuming that the algorithm has
/// completed a task with estimated RunTime toAdd
void GetEi::advanceProgress(double toAdd) {
466
467
468
469
470
  m_fracCompl += toAdd;
  progress(m_fracCompl);
  // look out for user cancel messgages
  interruption_point();
}
471
472
473

} // namespace Algorithms
} // namespace Mantid