GetEi.cpp 21.5 KB
Newer Older
1
2
3
4
5
6
7
/*WIKI* 

Uses E= (1/2)mv^2 to calculate the energy of neutrons leaving the source. The velocity is calculated from the time it takes for the neutron pulse to travel between the two monitors whose spectra were specified.

An initial energy guess is required for the algorithm to find the correct peak. The analysis will be done on the highest peak that is within 8% of the estimated TOF given by the estimate.

Not all neutrons arrive at the monitors at the same time because their kinetic energies, and therefore velocities, are all different. The time of arrival of the neutron pulse is taken to be the mean of the two half peak height locations. The half height points are found as follows:
8
9
10
11
# the peak height is the largest number of counts above the background in any bin in the window
# the half height is half the above number
# examine bins to the left of the bin with the highest number of counts looking for a bin with less than half that number above background
# interpolate between this point bin and the one immediately previous to find the first half height location
12
# repeat the steps 3 and 4 looking to the right of the highest point to get the second half height point
13
14
15
# the mean of the X-values of the two half height points is the TOF arrival time of the neutrons

The above process is illustrated on a peak is shown below in the image below
16
[[File:Monitorspect_getei.jpg|Monitor Peak|centre|618px]]
17
18
19
20

The distances between the monitors are read from the instrument definition file. It is assumed that the source and the monitors all lie on one line and that the monitors have the same delay time.

*WIKI*/
21

22
23
#include "MantidAlgorithms/GetEi.h"
#include "MantidKernel/ArrayProperty.h"
24
#include "MantidAPI/FileProperty.h"
25
26
27
28
29
#include "MantidKernel/PhysicalConstants.h"
#include "MantidAPI/WorkspaceValidators.h"
#include <boost/lexical_cast.hpp>
#include "MantidKernel/Exception.h" 
#include <cmath>
30
#include "MantidKernel/BoundedValidator.h"
31
32
33
34
35
36
37
38
39

namespace Mantid
{
namespace Algorithms
{

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

40
41
42
/// Sets documentation strings for this algorithm
void GetEi::initDocs()
{
Dereck Kachere's avatar
Dereck Kachere committed
43
  this->setWikiSummary("Calculates the kinetic energy of neutrons leaving the source based on the time it takes for them to travel between two monitors.");
44
45
46
47
  this->setOptionalMessage("Calculates the kinetic energy of neutrons leaving the source based on the time it takes for them to travel between two monitors.");
}


48
49
50
51
using namespace Kernel;
using namespace API;
using namespace Geometry;

52
// adjustable fit criteria, increase the first number or reduce any of the last three for more promiscuous peak fitting
53
// from the estimated location of the peak search forward by the following fraction and backward by the same fraction
54
const double GetEi::HALF_WINDOW = 8.0/100;
55
56
const double GetEi::PEAK_THRESH_H = 3.0;
const double GetEi::PEAK_THRESH_A = 5.0;
Peterson, Peter's avatar
Peterson, Peter committed
57
const int64_t GetEi::PEAK_THRESH_W = 3;
58
59
60
61
62

// progress estimates
const double GetEi::CROP = 0.15;
const double GetEi::GET_COUNT_RATE = 0.15;
const double GetEi::FIT_PEAK = 0.2;
63
64
65

/// Empty default constructor algorith() calls the constructor in the base class
GetEi::GetEi() : Algorithm(),
66
  m_tempWS(), m_fracCompl(0.0)
67
68
69
70
{
}

void GetEi::init()
71
{
72
  // Declare required input parameters for algorithm and do some validation here
73
74
75
76
  auto val = boost::make_shared<CompositeValidator>();
  val->add<WorkspaceUnitValidator>("TOF");
  val->add<HistogramValidator>();
  val->add<InstrumentValidator>();
77
  declareProperty(new WorkspaceProperty<>(
78
    "InputWorkspace","",Direction::Input,val),
79
80
    "The X units of this workspace must be time of flight with times in\n"
    "micro-seconds");
81
  auto mustBePositive = boost::make_shared<BoundedValidator<int> >();
82
83
84
85
  mustBePositive->setLower(0);
  declareProperty("Monitor1Spec", -1, mustBePositive,
    "The spectrum number of the output of the first monitor, e.g. MAPS\n"
    "41474, MARI 2, MERLIN 69634");
86
  declareProperty("Monitor2Spec", -1, mustBePositive,
87
88
    "The spectrum number of the output of the second monitor e.g. MAPS\n"
    "41475, MARI 3, MERLIN 69638");
89
  auto positiveDouble = boost::make_shared<BoundedValidator<double> >();
90
  positiveDouble->setLower(0);
91
92
  declareProperty("EnergyEstimate", -1.0, positiveDouble,
    "An approximate value for the typical incident energy, energy of\n"
93
    "neutrons leaving the source (meV)");
94
  declareProperty("IncidentEnergy", -1.0, Direction::Output);
95
  declareProperty("FirstMonitorPeak", -1.0, Direction::Output);
96
97
98
99
100
101

  m_fracCompl = 0.0;
}

/** Executes the algorithm
*  @throw out_of_range if the peak runs off the edge of the histogram
102
*  @throw NotFoundError if one of the requested spectrum numbers was not found in the workspace
103
104
105
106
107
*  @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
*/
void GetEi::exec()
{
108
  MatrixWorkspace_const_sptr inWS = getProperty("InputWorkspace");
109
110
  const specid_t mon1Spec = getProperty("Monitor1Spec");
  const specid_t mon2Spec = getProperty("Monitor2Spec");
111
  double dist2moni0 = -1, dist2moni1 = -1;
112
  getGeometry(inWS, mon1Spec, mon2Spec, dist2moni0, dist2moni1);
113
114
115

  // 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
  const double E_est = getProperty("EnergyEstimate");
116
  // we're assuming that the time units for the X-values in the workspace are micro-seconds
117
  const double peakLoc0 = 1e6*timeToFly(dist2moni0, E_est);
118
119
  // 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 will be searched for at TOF " << peakLoc0 << " micro seconds +/-" << boost::lexical_cast<std::string>(100.0*HALF_WINDOW) << "%\n";
120
  const double peakLoc1 = 1e6*timeToFly(dist2moni1, E_est);
121
  g_log.information() << "Based on the user selected energy the second peak will be searched for at TOF " << peakLoc1 << " micro seconds +/-" << boost::lexical_cast<std::string>(100.0*HALF_WINDOW) << "%\n";
122

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

126
127
  g_log.information() << "Looking for a peak in the first monitor spectrum, spectra index " << indexes[0] << std::endl;
  double t_monitor0 = getPeakCentre(inWS, indexes[0], peakLoc0);
128
129
130
  g_log.notice() << "The first peak has been found at TOF = " << t_monitor0 << " microseconds\n";
  setProperty("FirstMonitorPeak", t_monitor0);

131
132
  g_log.information() << "Looking for a peak in the second monitor spectrum, spectra index " << indexes[1] << std::endl;
  double t_monitor1 = getPeakCentre(inWS, indexes[1], peakLoc1);
133
  g_log.information() << "The second peak has been found at TOF = " << t_monitor1 << " microseconds\n";
134
135
136
137
138
139

  // 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));

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

  setProperty("IncidentEnergy", E_i);
}
/** Gets the distances between the source and detectors whose IDs you pass to it
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
150
151
*  @throw NotFoundError if no detector is found for the detector ID given
*/
Doucet, Mathieu's avatar
Doucet, Mathieu committed
152
void GetEi::getGeometry(API::MatrixWorkspace_const_sptr WS, specid_t mon0Spec, specid_t mon1Spec, double &monitor0Dist, double &monitor1Dist) const
153
{
154
  const IObjComponent_const_sptr source = WS->getInstrument()->getSource();
155
156

  // retrieve a pointer to the first detector and get its distance
157
158
159
160
161
  size_t monWI = 0;
  try
  {
    monWI = WS->getIndexFromSpectrumNumber(mon0Spec);
  }
162
  catch (std::runtime_error &)
163
164
165
166
167
168
169
  {
    g_log.error() << "Could not find the workspace index for the monitor at spectrum " << mon0Spec << "\n";
    g_log.error() << "Error retrieving data for the first monitor" << std::endl;
    throw std::bad_cast();
  }
  const std::set<detid_t> & 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";
    g_log.error() << "Error retrieving data for the first monitor" << std::endl;
    throw std::bad_cast();
175
  }
176
  IDetector_const_sptr det = WS->getInstrument()->getDetector(*dets.begin());
177
178
  monitor0Dist = det->getDistance(*(source.get()));

179
  // repeat for the second detector
180
181
182
183
  try
  {
    monWI = WS->getIndexFromSpectrumNumber(mon0Spec);
  }
184
  catch (std::runtime_error &)
185
186
187
188
189
190
191
  {
    g_log.error() << "Could not find the workspace index for the monitor at spectrum " << mon0Spec << "\n";
    g_log.error() << "Error retrieving data for the second monitor\n";
    throw std::bad_cast();
  }
  const std::set<detid_t> & dets2 = WS->getSpectrum(monWI)->getDetectorIDs();
  if ( dets2.size() != 1 )
192
193
194
195
  {
    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";
    g_log.error() << "Error retrieving data for the second monitor\n";
    throw std::bad_cast();
196
  }
197
  det = WS->getInstrument()->getDetector(*dets2.begin());
198
  monitor1Dist = det->getDistance(*(source.get()));
199
}
200
/** Converts detector IDs to spectra indexes
201
202
203
*  @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
204
*  @return the indexes of the histograms created by the detector whose ID were passed
205
*  @throw NotFoundError if one of the requested spectrum numbers was not found in the workspace
206
*/
Doucet, Mathieu's avatar
Doucet, Mathieu committed
207
std::vector<size_t> GetEi::getMonitorSpecIndexs(API::MatrixWorkspace_const_sptr WS, specid_t specNum1, specid_t specNum2) const
208
{// getting spectra numbers from detector IDs is 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
Janik Zikovsky's avatar
Janik Zikovsky committed
209
  std::vector<size_t> specInds;
210
211
  
  // get the index number of the histogram for the first monitor
Janik Zikovsky's avatar
Janik Zikovsky committed
212
  std::vector<specid_t> specNumTemp(&specNum1, &specNum1+1);
213
  WS->getIndicesFromSpectra(specNumTemp, specInds);
214
215
216
217
218
219
220
  if ( specInds.size() != 1 )
  {// the monitor spectrum isn't present in the workspace, we can't continue from here
    g_log.error() << "Couldn't find the first monitor spectrum, number " << specNum1 << std::endl;
    throw Exception::NotFoundError("GetEi::getMonitorSpecIndexs()", specNum1);
  }

  // nowe the second monitor
Janik Zikovsky's avatar
Janik Zikovsky committed
221
  std::vector<size_t> specIndexTemp;
222
  specNumTemp[0] = specNum2;
223
  WS->getIndicesFromSpectra(specNumTemp, specIndexTemp);
224
225
226
227
228
229
230
  if ( specIndexTemp.size() != 1 )
  {// the monitor spectrum isn't present in the workspace, we can't continue from here
    g_log.error() << "Couldn't find the second monitor spectrum, number " << specNum2 << std::endl;
    throw Exception::NotFoundError("GetEi::getMonitorSpecIndexs()", specNum2);
  }
  
  specInds.push_back(specIndexTemp[0]);
231
232
  return specInds;
}
233
234
/** Uses E_KE = mv^2/2 and s = vt to calculate the time required for a neutron
*  to travel a distance, s
235
236
* @param s :: ditance travelled in meters
* @param E_KE :: kinetic energy in meV
237
238
239
240
241
242
243
244
245
246
247
248
249
* @return the time to taken to travel that uninterrupted distance in seconds
*/
double GetEi::timeToFly(double s, double E_KE) const
{
  // 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;

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

/** 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
253
254
255
*  @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
256
*  @return a time of flight value in the peak in microseconds
257
258
*  @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
259
*  @throw runtime_error a Child Algorithm just falls over
260
*/
Peterson, Peter's avatar
Peterson, Peter committed
261
double GetEi::getPeakCentre(API::MatrixWorkspace_const_sptr WS, const int64_t monitIn, const double peakTime)
262
{
263
  const MantidVec& timesArray = WS->readX(monitIn);
264
  // we search for the peak only inside some window because there are often more peaks in the monitor histogram
265
  double halfWin = ( timesArray.back() - timesArray.front() )*HALF_WINDOW;
266
  // 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
267
  extractSpec(monitIn, peakTime-halfWin, peakTime+halfWin);
268
269
270
271
272
  // converting the workspace to count rate is required by the fitting algorithm if the bin widths are not all the same
  WorkspaceHelpers::makeDistribution(m_tempWS);
  // look out for user cancel messgages as the above command can take a bit of time
  advanceProgress(GET_COUNT_RATE);

273
  // to store fit results
Peterson, Peter's avatar
Peterson, Peter committed
274
  int64_t centreGausInd;
275
276
  double height, backGroundlev;
  getPeakEstimates(height, centreGausInd, backGroundlev);
277
278
  // look out for user cancel messgages
  advanceProgress(FIT_PEAK);
279
280
281

  // the peak centre is defined as the centre of the two half maximum points as this is better for asymmetric peaks
  // first loop backwards along the histogram to get the first half height point
282
  const double lHalf = findHalfLoc(centreGausInd, height, backGroundlev, GO_LEFT);
283
  // go forewards to get the half height on the otherside of the peak
284
285
286
  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;
287
}
288
/** Calls CropWorkspace as a Child Algorithm and passes to it the InputWorkspace property
289
290
291
*  @param specInd :: 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)
292
*  @throw out_of_range if start, end or specInd are set outside of the vaild range for the workspace
293
*  @throw runtime_error if the algorithm just falls over
294
295
*  @throw invalid_argument if the input workspace does not have common binning
*/
Peterson, Peter's avatar
Peterson, Peter committed
296
void GetEi::extractSpec(int64_t specInd, double start, double end)
297
298
{
  IAlgorithm_sptr childAlg =
299
    createChildAlgorithm("CropWorkspace", 100*m_fracCompl, 100*(m_fracCompl+CROP) );
300
301
  m_fracCompl += CROP;
  
302
  childAlg->setProperty<MatrixWorkspace_sptr>("InputWorkspace",getProperty("InputWorkspace") );
303
304
305
306
  childAlg->setProperty( "XMin", start);
  childAlg->setProperty( "XMax", end);
  childAlg->setProperty( "StartWorkspaceIndex", specInd);
  childAlg->setProperty( "EndWorkspaceIndex", specInd);
307
  childAlg->executeAsChildAlg();
308
309
310
311
312
313
314
315
316
317
318

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

//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);
  progress(m_fracCompl);
  interruption_point();
}

/** Finds the largest peak by looping through the histogram and finding the maximum
*  value 
319
320
321
* @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
322
* @throw invalid_argument if the peak is not clearly above the background
323
*/
Peterson, Peter's avatar
Peterson, Peter committed
324
void GetEi::getPeakEstimates(double &height, int64_t &centreInd, double &background) const
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
{
  // take note of the number of background counts as error checking, do we have a peak or just a bump in the background
  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
  for ( MantidVec::size_type i = 1; i < m_tempWS->readY(0).size()-1; ++i )
  {
    background += m_tempWS->readY(0)[i];
    if ( m_tempWS->readY(0)[i] > height )
    {
      centreInd = i;
      height = m_tempWS->readY(0)[centreInd];
    }
  }
  
Peterson, Peter's avatar
Peterson, Peter committed
342
  background = background/static_cast<double>(m_tempWS->readY(0).size());
343
344
345
346
347
348
349
350
  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?");
  }

  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])/2 << " (peak height " << height << " counts/microsecond)\n";

}
351
352
/** Estimates the closest time, looking either or back, when the number of counts is
*  half that in the bin whose index that passed
353
354
355
356
*  @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
357
*  @return estimated TOF of the half maximum point
358
359
360
*  @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
*/
361
double GetEi::findHalfLoc(MantidVec::size_type startInd, const double height, const double noise, const direction go) const
362
363
{
  MantidVec::size_type endInd = startInd;
364

365
  while ( m_tempWS->readY(0)[endInd] >  (height+noise)/2.0 )
366
367
  {
    endInd += go;
368
    if ( endInd < 1 )
369
    {
370
      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?");
371
    }
372
    if ( endInd > m_tempWS->readY(0).size()-2)
373
    {
374
      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?");
375
376
    }
  }
377

Peterson, Peter's avatar
Peterson, Peter committed
378
  if ( std::abs(static_cast<int64_t>(endInd - startInd)) < PEAK_THRESH_W )
379
380
381
382
383
  {// 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";
  }
  // we have a peak in range, do an area check to see if the peak has any significance
  double hOverN = (height-noise)/noise;
Peterson, Peter's avatar
Peterson, Peter committed
384
  if ( hOverN < PEAK_THRESH_A && std::abs(hOverN*static_cast<double>(endInd - startInd)) < PEAK_THRESH_A )
385
386
387
  {// 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?");
  }
388
389
390
  // 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
391
392
  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
393
394
395
396
397
398
    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;
399
400
  }

401
  g_log.debug() << "One half height point found at TOF = " << halfTime << " microseconds\n";
402
  return halfTime;
403
}
404
/** Get the kinetic energy of a neuton in joules given it speed using E=mv^2/2
405
*  @param speed :: the instantanious speed of a neutron in metres per second
406
407
*  @return the energy in joules
*/
408
409
410
411
412
413
double GetEi::neutron_E_At(double speed) const
{
  // E_KE = mv^2/2
  return PhysicalConstants::NeutronMass*speed*speed/(2);
}

414
415
416
417
418
419
420
421
/// Update the percentage complete estimate assuming that the algorithm has completed a task with estimated RunTime toAdd
void GetEi::advanceProgress(double toAdd)
{
  m_fracCompl += toAdd;
  progress(m_fracCompl);
  // look out for user cancel messgages
  interruption_point();
}
422
423
424

} // namespace Algorithms
} // namespace Mantid