GetEi.cpp 21 KB
Newer Older
1
#include "MantidAlgorithms/GetEi.h"
2
3
4
#include "MantidAPI/HistogramValidator.h"
#include "MantidAPI/InstrumentValidator.h"
#include "MantidAPI/WorkspaceUnitValidator.h"
5
#include "MantidGeometry/Instrument.h"
6
7
8
#include "MantidKernel/Exception.h"
#include "MantidKernel/CompositeValidator.h"
#include "MantidKernel/BoundedValidator.h"
9
#include "MantidKernel/PhysicalConstants.h"
10

11
12
13
#include <boost/lexical_cast.hpp>
#include <cmath>

14
15
namespace Mantid {
namespace Algorithms {
16
17
18
19
20
21
22
23

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

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

24
25
26
27
28
// 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;
29
30
const double GetEi::PEAK_THRESH_H = 3.0;
const double GetEi::PEAK_THRESH_A = 5.0;
Peterson, Peter's avatar
Peterson, Peter committed
31
const int64_t GetEi::PEAK_THRESH_W = 3;
32
33
34
35
36

// progress estimates
const double GetEi::CROP = 0.15;
const double GetEi::GET_COUNT_RATE = 0.15;
const double GetEi::FIT_PEAK = 0.2;
37
38

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

41
void GetEi::init() {
42
  // Declare required input parameters for algorithm and do some validation here
43
44
45
46
  auto val = boost::make_shared<CompositeValidator>();
  val->add<WorkspaceUnitValidator>("TOF");
  val->add<HistogramValidator>();
  val->add<InstrumentValidator>();
47
  declareProperty(
48
49
      make_unique<WorkspaceProperty<>>("InputWorkspace", "", Direction::Input,
                                       val),
50
51
52
      "The X units of this workspace must be time of flight with times in\n"
      "micro-seconds");
  auto mustBePositive = boost::make_shared<BoundedValidator<int>>();
53
  mustBePositive->setLower(0);
54
55
56
57
58
59
60
61
62
  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");
  auto positiveDouble = boost::make_shared<BoundedValidator<double>>();
63
  positiveDouble->setLower(0);
64
65
66
67
  declareProperty(
      "EnergyEstimate", -1.0, positiveDouble,
      "An approximate value for the typical incident energy, energy of\n"
      "neutrons leaving the source (meV)");
68
  declareProperty("IncidentEnergy", -1.0, Direction::Output);
69
  declareProperty("FirstMonitorPeak", -1.0, Direction::Output);
70
71
72
73
74
75

  m_fracCompl = 0.0;
}

/** Executes the algorithm
*  @throw out_of_range if the peak runs off the edge of the histogram
76
77
78
79
80
81
*  @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
82
*/
83
void GetEi::exec() {
84
  MatrixWorkspace_const_sptr inWS = getProperty("InputWorkspace");
85
86
  const specnum_t mon1Spec = getProperty("Monitor1Spec");
  const specnum_t mon2Spec = getProperty("Monitor2Spec");
87
  double dist2moni0 = -1, dist2moni1 = -1;
88
  getGeometry(inWS, mon1Spec, mon2Spec, dist2moni0, dist2moni1);
89

90
91
92
  // 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
93
  const double E_est = getProperty("EnergyEstimate");
94
95
96
97
98
99
  // 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 "
Tom Perkins's avatar
Tom Perkins committed
100
101
                         "will be searched for at TOF " << peakLoc0
                      << " micro seconds +/-"
102
103
104
105
                      << boost::lexical_cast<std::string>(100.0 * HALF_WINDOW)
                      << "%\n";
  const double peakLoc1 = 1e6 * timeToFly(dist2moni1, E_est);
  g_log.information() << "Based on the user selected energy the second peak "
Tom Perkins's avatar
Tom Perkins committed
106
107
                         "will be searched for at TOF " << peakLoc1
                      << " micro seconds +/-"
108
109
110
111
                      << boost::lexical_cast<std::string>(100.0 * HALF_WINDOW)
                      << "%\n";

  // get the histograms created by the monitors
Janik Zikovsky's avatar
Janik Zikovsky committed
112
  std::vector<size_t> indexes = getMonitorSpecIndexs(inWS, mon1Spec, mon2Spec);
113

114
115
116
  g_log.information()
      << "Looking for a peak in the first monitor spectrum, spectra index "
      << indexes[0] << std::endl;
117
  double t_monitor0 = getPeakCentre(inWS, indexes[0], peakLoc0);
118
119
  g_log.notice() << "The first peak has been found at TOF = " << t_monitor0
                 << " microseconds\n";
120
121
  setProperty("FirstMonitorPeak", t_monitor0);

122
123
124
  g_log.information()
      << "Looking for a peak in the second monitor spectrum, spectra index "
      << indexes[1] << std::endl;
125
  double t_monitor1 = getPeakCentre(inWS, indexes[1], peakLoc1);
126
127
  g_log.information() << "The second peak has been found at TOF = "
                      << t_monitor1 << " microseconds\n";
128

129
130
131
132
  // 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
  double meanSpeed =
      (dist2moni1 - dist2moni0) / (1e-6 * (t_monitor1 - t_monitor0));
133

134
135
136
137
138
139
  // 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;
  g_log.notice() << "The incident energy has been calculated to be " << E_i
                 << " meV"
                 << " (your estimate was " << E_est << " meV)\n";
140
141
142
143

  setProperty("IncidentEnergy", E_i);
}
/** Gets the distances between the source and detectors whose IDs you pass to it
144
145
146
*  @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
147
148
149
150
*  @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
151
152
*  @throw NotFoundError if no detector is found for the detector ID given
*/
153
154
void GetEi::getGeometry(API::MatrixWorkspace_const_sptr WS, specnum_t mon0Spec,
                        specnum_t mon1Spec, double &monitor0Dist,
155
                        double &monitor1Dist) const {
156
  const IComponent_const_sptr source = WS->getInstrument()->getSource();
157
158

  // retrieve a pointer to the first detector and get its distance
159
  size_t monWI = 0;
160
  try {
161
    monWI = WS->getIndexFromSpectrumNumber(mon0Spec);
162
163
164
165
  } catch (std::runtime_error &) {
    g_log.error()
        << "Could not find the workspace index for the monitor at spectrum "
        << mon0Spec << "\n";
166
167
168
    g_log.error() << "Error retrieving data for the first monitor" << std::endl;
    throw std::bad_cast();
  }
169
  const auto &dets = WS->getSpectrum(monWI)->getDetectorIDs();
170

171
172
173
174
  if (dets.size() != 1) {
    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";
175
176
    g_log.error() << "Error retrieving data for the first monitor" << std::endl;
    throw std::bad_cast();
177
  }
178
  IDetector_const_sptr det = WS->getInstrument()->getDetector(*dets.begin());
179
180
  monitor0Dist = det->getDistance(*(source.get()));

181
  // repeat for the second detector
182
  try {
183
    monWI = WS->getIndexFromSpectrumNumber(mon0Spec);
184
185
186
187
  } catch (std::runtime_error &) {
    g_log.error()
        << "Could not find the workspace index for the monitor at spectrum "
        << mon0Spec << "\n";
188
189
190
    g_log.error() << "Error retrieving data for the second monitor\n";
    throw std::bad_cast();
  }
191
  const auto &dets2 = WS->getSpectrum(monWI)->getDetectorIDs();
192
193
194
195
  if (dets2.size() != 1) {
    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";
196
197
    g_log.error() << "Error retrieving data for the second monitor\n";
    throw std::bad_cast();
198
  }
199
  det = WS->getInstrument()->getDetector(*dets2.begin());
200
  monitor1Dist = det->getDistance(*(source.get()));
201
}
202
/** Converts detector IDs to spectra indexes
203
204
205
*  @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
206
207
208
209
*  @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
210
*/
211
std::vector<size_t> GetEi::getMonitorSpecIndexs(
212
213
    API::MatrixWorkspace_const_sptr WS, specnum_t specNum1,
    specnum_t specNum2) const { // getting spectra numbers from detector IDs is
Nick Draper's avatar
Nick Draper committed
214
215
216
                                // hard because the map works the other way,
  // getting index numbers from spectra numbers has
  // the same problem and we are about to do both
217

218
  // get the index number of the histogram for the first monitor
219

220
  std::vector<specnum_t> specNumTemp(&specNum1, &specNum1 + 1);
221
222
223
  auto wsInds = WS->getIndicesFromSpectra(specNumTemp);

  if (wsInds.size() != 1) { // the monitor spectrum isn't present in the
Nick Draper's avatar
Nick Draper committed
224
                            // workspace, we can't continue from here
225
226
    g_log.error() << "Couldn't find the first monitor spectrum, number "
                  << specNum1 << std::endl;
227
228
229
230
231
    throw Exception::NotFoundError("GetEi::getMonitorSpecIndexs()", specNum1);
  }

  // nowe the second monitor
  specNumTemp[0] = specNum2;
232
233
  auto wsIndexTemp = WS->getIndicesFromSpectra(specNumTemp);
  if (wsIndexTemp.size() != 1) { // the monitor spectrum isn't present in the
Nick Draper's avatar
Nick Draper committed
234
                                 // workspace, we can't continue from here
235
236
    g_log.error() << "Couldn't find the second monitor spectrum, number "
                  << specNum2 << std::endl;
237
238
    throw Exception::NotFoundError("GetEi::getMonitorSpecIndexs()", specNum2);
  }
239

Nick Draper's avatar
Nick Draper committed
240
  wsInds.push_back(specNumTemp[0]);
241
  return wsInds;
242
}
243
244
/** Uses E_KE = mv^2/2 and s = vt to calculate the time required for a neutron
*  to travel a distance, s
245
246
* @param s :: ditance travelled in meters
* @param E_KE :: kinetic energy in meV
247
248
* @return the time to taken to travel that uninterrupted distance in seconds
*/
249
double GetEi::timeToFly(double s, double E_KE) const {
250
251
252
253
254
255
256
  // 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;

257
  return s / sqrt(2 * E_KE / PhysicalConstants::NeutronMass);
258
}
259

260
261
262
263
/** Looks for and examines a peak close to that specified by the input
* parameters and
*  examines it to find a representative time for when the neutrons hit the
* detector
264
*  @param WS :: the workspace containing the monitor spectrum
265
266
267
268
*  @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
269
*  @return a time of flight value in the peak in microseconds
270
271
*  @throw invalid_argument if a good peak fit wasn't made or the input workspace
* does not have common binning
272
*  @throw out_of_range if the peak runs off the edge of the histogram
273
*  @throw runtime_error a Child Algorithm just falls over
274
*/
275
276
277
278
279
280
281
282
283
284
285
double GetEi::getPeakCentre(API::MatrixWorkspace_const_sptr WS,
                            const int64_t monitIn, const double peakTime) {
  const MantidVec &timesArray = WS->readX(monitIn);
  // 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;
  // runs CropWorkspace as a Child Algorithm to and puts the result in a new
  // temporary workspace that will be deleted when this algorithm has finished
  extractSpec(monitIn, peakTime - halfWin, peakTime + halfWin);
  // converting the workspace to count rate is required by the fitting algorithm
  // if the bin widths are not all the same
286
  WorkspaceHelpers::makeDistribution(m_tempWS);
287
288
  // look out for user cancel messgages as the above command can take a bit of
  // time
289
290
  advanceProgress(GET_COUNT_RATE);

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

298
299
  // the peak centre is defined as the centre of the two half maximum points as
  // this is better for asymmetric peaks
300
  // first loop backwards along the histogram to get the first half height point
301
302
  const double lHalf =
      findHalfLoc(centreGausInd, height, backGroundlev, GO_LEFT);
303
  // go forewards to get the half height on the otherside of the peak
304
305
306
307
  const double rHalf =
      findHalfLoc(centreGausInd, height, backGroundlev, GO_RIGHT);
  // the peak centre is defined as the mean of the two half height times
  return (lHalf + rHalf) / 2;
308
}
309
310
/** Calls CropWorkspace as a Child Algorithm and passes to it the InputWorkspace
* property
311
*  @param specInd :: the index number of the histogram to extract
312
313
314
315
316
317
*  @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 specInd are set outside of the vaild
* range for the workspace
318
*  @throw runtime_error if the algorithm just falls over
319
320
*  @throw invalid_argument if the input workspace does not have common binning
*/
321
322
323
void GetEi::extractSpec(int64_t specInd, double start, double end) {
  IAlgorithm_sptr childAlg = createChildAlgorithm(
      "CropWorkspace", 100 * m_fracCompl, 100 * (m_fracCompl + CROP));
324
  m_fracCompl += CROP;
325
326
327
328
329
330
331

  childAlg->setProperty<MatrixWorkspace_sptr>("InputWorkspace",
                                              getProperty("InputWorkspace"));
  childAlg->setProperty("XMin", start);
  childAlg->setProperty("XMax", end);
  childAlg->setProperty("StartWorkspaceIndex", specInd);
  childAlg->setProperty("EndWorkspaceIndex", specInd);
332
  childAlg->executeAsChildAlg();
333
334
335

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

336
337
338
  // 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);
339
340
341
342
  progress(m_fracCompl);
  interruption_point();
}

343
344
345
/** Finds the largest peak by looping through the histogram and finding the
* maximum
*  value
346
* @param height :: its passed value ignored it is set to the peak height
347
348
349
350
* @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
351
* @throw invalid_argument if the peak is not clearly above the background
352
*/
353
354
355
356
void GetEi::getPeakEstimates(double &height, int64_t &centreInd,
                             double &background) const {
  // take note of the number of background counts as error checking, do we have
  // a peak or just a bump in the background
357
358
359
360
361
  background = 0;
  // start at the first Y value
  height = m_tempWS->readY(0)[0];
  centreInd = 0;
  // then loop through all the Y values and find the tallest peak
362
  for (MantidVec::size_type i = 1; i < m_tempWS->readY(0).size() - 1; ++i) {
363
    background += m_tempWS->readY(0)[i];
364
    if (m_tempWS->readY(0)[i] > height) {
365
366
367
368
369
      centreInd = i;
      height = m_tempWS->readY(0)[centreInd];
    }
  }

370
371
372
373
374
375
376
377
  background = background / static_cast<double>(m_tempWS->readY(0).size());
  if (height < PEAK_THRESH_H * background) {
    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?");
  }
378

379
380
381
382
  g_log.debug() << "Peak position is the bin that has the maximum Y value in "
                   "the monitor spectrum, which is at TOF "
                << (m_tempWS->readX(0)[centreInd] +
                    m_tempWS->readX(0)[centreInd + 1]) /
Tom Perkins's avatar
Tom Perkins committed
383
384
                       2 << " (peak height " << height
                << " counts/microsecond)\n";
385
}
386
387
/** Estimates the closest time, looking either or back, when the number of
* counts is
388
*  half that in the bin whose index that passed
389
390
391
392
*  @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
393
394
*  @param noise :: mean number of counts in each bin in the workspace
*  @param go :: either GetEi::GO_LEFT or GetEi::GO_RIGHT
395
*  @return estimated TOF of the half maximum point
396
397
*  @throw out_of_range if the end of the histogram is reached before the point
* is found
398
399
*  @throw invalid_argument if the peak is too thin
*/
400
401
double GetEi::findHalfLoc(MantidVec::size_type startInd, const double height,
                          const double noise, const direction go) const {
402
  MantidVec::size_type endInd = startInd;
403

404
  while (m_tempWS->readY(0)[endInd] > (height + noise) / 2.0) {
405
    endInd += go;
406
407
408
409
410
411
    if (endInd < 1) {
      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?");
412
    }
413
414
415
416
417
418
    if (endInd > m_tempWS->readY(0).size() - 2) {
      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?");
419
420
    }
  }
421

422
423
424
425
426
427
428
  if (std::abs(static_cast<int64_t>(endInd - startInd)) <
      PEAK_THRESH_W) { // we didn't find a significant peak
    g_log.error() << "Likely precision problem or error, one half height "
                     "distance is less than the threshold number of bins from "
                     "the central peak: "
                  << std::abs(static_cast<int>(endInd - startInd)) << "<"
                  << PEAK_THRESH_W << ". Check the monitor peak\n";
429
  }
430
431
432
433
434
435
436
437
438
439
440
441
  // 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 &&
      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?");
442
  }
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
  // get the TOF value in the middle of the bin with the first value below the
  // half height
  double halfTime =
      (m_tempWS->readX(0)[endInd] + m_tempWS->readX(0)[endInd + 1]) / 2;
  // 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
    double gradient =
        (m_tempWS->readY(0)[endInd] - m_tempWS->readY(0)[endInd - go]) /
        (m_tempWS->readX(0)[endInd] - m_tempWS->readX(0)[endInd - go]);
    // 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
    double deltaY = m_tempWS->readY(0)[endInd] - (height + noise) / 2.0;
    // correct for the interpolation back in the direction towards the peak
    // centre
    halfTime -= deltaY / gradient;
463
464
  }

465
466
  g_log.debug() << "One half height point found at TOF = " << halfTime
                << " microseconds\n";
467
  return halfTime;
468
}
469
/** Get the kinetic energy of a neuton in joules given it speed using E=mv^2/2
470
*  @param speed :: the instantanious speed of a neutron in metres per second
471
472
*  @return the energy in joules
*/
473
double GetEi::neutron_E_At(double speed) const {
474
  // E_KE = mv^2/2
475
  return PhysicalConstants::NeutronMass * speed * speed / (2);
476
477
}

478
479
480
/// Update the percentage complete estimate assuming that the algorithm has
/// completed a task with estimated RunTime toAdd
void GetEi::advanceProgress(double toAdd) {
481
482
483
484
485
  m_fracCompl += toAdd;
  progress(m_fracCompl);
  // look out for user cancel messgages
  interruption_point();
}
486
487
488

} // namespace Algorithms
} // namespace Mantid