GetEi.cpp 17.6 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
#include "MantidAlgorithms/GetEi.h"
#include "MantidKernel/ArrayProperty.h"
#include "MantidKernel/FileProperty.h"
#include "MantidKernel/PhysicalConstants.h"
#include "MantidAPI/WorkspaceValidators.h"
#include "MantidAPI/SpectraDetectorMap.h"
#include <boost/lexical_cast.hpp>
#include "MantidKernel/Exception.h" 
#include <cmath>

namespace Mantid
{
namespace Algorithms
{

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

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

// from the estimated location of the peak search forward by the following fraction and backward by the same fraction
const double GetEi::HALF_WINDOW = 4.0/100;
const double GetEi::PEAK_THRESH_H = 3.0;
const int GetEi::PEAK_THRESH_W = 3;

// progress estimates
const double GetEi::CROP = 0.15;
const double GetEi::GET_COUNT_RATE = 0.15;
const double GetEi::FIT_PEAK = 2.0;

/// Empty default constructor algorith() calls the constructor in the base class
GetEi::GetEi() : Algorithm(),
36
  m_tempWS(), m_fracCompl(0.0)
37
38
39
40
41
42
43
44
45
46
47
48
49
50
{
}

void GetEi::init()
{
  // Declare required input parameters for algorithm
  WorkspaceUnitValidator<Workspace2D> *val =
    new WorkspaceUnitValidator<Workspace2D>("TOF");
  declareProperty(new WorkspaceProperty<Workspace2D>(
    "InputWorkspace","",Direction::Input,val),
    "The X units of this workspace must be time of flight with times in\n"
    "micro-seconds");
  BoundedValidator<int> *mustBePositive = new BoundedValidator<int>();
  mustBePositive->setLower(0);
51
52
53
54
  declareProperty("Monitor1Spec", -1, mustBePositive,
    "The spectrum number of the output of the first monitor");
  declareProperty("Monitor2Spec", -1, mustBePositive->clone(),
    "The spectrum number of the output of the second monitor");
55
56
57
58
59
60
61
62
63
64
65
66
  BoundedValidator<double> *positiveDouble = new BoundedValidator<double>();
  positiveDouble->setLower(0);
  declareProperty("EnergyEstimate", -1.0, positiveDouble,
    "An approximate value for the typical incident energy, energy of\n"
    "neutrons leaving the source (meV)");
  declareProperty("IncidentEnergy", -1.0, Direction::Output);

  m_fracCompl = 0.0;
}

/** Executes the algorithm
*  @throw out_of_range if the peak runs off the edge of the histogram
67
*  @throw NotFoundError if one of the requested spectrum numbers was not found in the workspace
68
69
*  @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
70
*  @throw runtime_error if there is a problem with the SpectraDetectorMap or a sub-algorithm falls over
71
72
73
*/
void GetEi::exec()
{
74
75
76
  Workspace2D_const_sptr inWS = getProperty("InputWorkspace");
  const int mon1Spec = getProperty("Monitor1Spec");
  const int mon2Spec = getProperty("Monitor2Spec");
77
  double dist2moni0 = -1, dist2moni1 = -1;
78
  getGeometry(inWS, mon1Spec, mon2Spec, dist2moni0, dist2moni1);
79
80
81
82
83
84
85
86
87
88

  // 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");
  // 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 will be searched for at TOF " << peakLoc0 << " micro seconds +/-" << boost::lexical_cast<std::string>(100.0*HALF_WINDOW) << "%" << std::endl;
  const double peakLoc1 = 1e6*timeToFly(dist2moni1, E_est);
  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) << "%" << std::endl;

89
90
91
    // get the histograms created by the monitors
  std::vector<int> indexes = getMonitorSpecIndexs(inWS, mon1Spec, mon2Spec);

92
  g_log.information() << "Looking for a peak in the first monitor spectrum, spectra index " << indexes[0] << std::endl;
93
  double t_monitor0 = getPeakCentre(inWS, indexes[0], peakLoc0);
94
95
  g_log.information() << "The first peak has been found at TOF = " << t_monitor0 << std::endl;
  g_log.information() << "Looking for a peak in the second monitor spectrum, spectra index " << indexes[1] << std::endl;
96
  double t_monitor1 = getPeakCentre(inWS, indexes[1], peakLoc1);
97
98
99
100
101
102
103
  g_log.information() << "The second peak has been found at TOF = " << t_monitor1 << std::endl;

  // 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;
104
  g_log.notice() << "The incident energy has been calculated to be " << E_i << " meV" << " (your estimate was " << E_est << " meV)" << std::endl;
105
106
107
108

  setProperty("IncidentEnergy", E_i);
}
/** Gets the distances between the source and detectors whose IDs you pass to it
109
110
111
*  @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
112
113
114
*  @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
115
*  @throw runtime_error if there is a problem with the SpectraDetectorMap
116
*/
117
void GetEi::getGeometry(Workspace2D_const_sptr WS, int mon0Spec, int mon1Spec, double &monitor0Dist, double &monitor1Dist) const
118
{
119
  const IObjComponent_sptr source = WS->getInstrument()->getSource();
120
121

  // retrieve a pointer to the first detector and get its distance
122
123
124
125
  std::vector<int> dets = WS->spectraMap().getDetectors(mon0Spec);
  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" << std::endl;
126
127
    g_log.error() << "Error retrieving data for the first monitor" << std::endl;
    throw std::bad_cast();
128
129
130
  }
  IDetector_sptr det = WS->getInstrument()->getDetector(dets[0]);
  monitor0Dist = det->getDistance(*(source.get()));
131

132
  // repeat for the second detector
133
134
135
136
  dets = WS->spectraMap().getDetectors(mon1Spec);
  if ( dets.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" << std::endl;
137
138
    g_log.error() << "Error retrieving data for the second monitor" << std::endl;
    throw std::bad_cast();
139
140
141
  }
  det = WS->getInstrument()->getDetector(dets[0]);
  monitor1Dist = det->getDistance(*(source.get()));
142
}
143
144
/** Converts detector IDs to spectra indexes
*  @param WS the workspace on which the calculations are being performed
145
146
*  @param specNum1 spectrum number of the output of the first monitor
*  @param specNum2 spectrum number of the output of the second monitor
147
*  @return the indexes of the histograms created by the detector whose ID were passed
148
*  @throw NotFoundError if one of the requested spectrum numbers was not found in the workspace
149
*/
150
std::vector<int> GetEi::getMonitorSpecIndexs(Workspace2D_const_sptr WS, int specNum1, int specNum2) const
151
152
{// 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
  std::vector<int> specInds;
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
  
  // get the index number of the histogram for the first monitor
  std::vector<int> specNumTemp(&specNum1, &specNum1+1);
  WorkspaceHelpers::getIndicesFromSpectra(WS, specNumTemp, specInds);
  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
  std::vector<int> specIndexTemp;
  specNumTemp[0] = specNum2;
  WorkspaceHelpers::getIndicesFromSpectra(WS, specNumTemp, specIndexTemp);
  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]);
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
  return specInds;
}
/** Uses E_KE = mv^2/2 and s = vt to calculate the time required for a neutron
*  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
*/
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 t^-2
  E_KE *= PhysicalConstants::meV;

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

/** 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
*  @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 half way between the two half height locations on the peak
*  @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 sub-algorithm just falls over
*/
double GetEi::getPeakCentre(Workspace2D_const_sptr WS, const int monitIn, const double peakTime)
{
  // we search for the peak only inside some window because there are often more peaks in the monitor histogram
  double halfWin = WS->readX(monitIn).back()*HALF_WINDOW;
  // runs CropWorkspace as a sub-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
  WorkspaceHelpers::makeDistribution(m_tempWS);
  // look out for user cancel messgages as the above command can take a bit of time
  advanceProgress(GET_COUNT_RATE);

  // to store fit results
  double height, centreGaussian;
  getPeakEstimates(height, centreGaussian);
  // look out for user cancel messgages
  advanceProgress(FIT_PEAK);

  // find the index of the centre point. The centre can't be at index zero as this is at the edge of the spectrum, so centreIndex = 0 is the error value
  MantidVec::size_type cenGausIn = 0;
223
  for ( MantidVec::size_type i = 0; i < m_tempWS->readY(0).size(); ++i )
224
225
226
227
228
229
230
  {// assumes that the bin boundaries are all in order of increasing time
    if ( m_tempWS->readX(0)[i] > centreGaussian )
    {
      cenGausIn = i;
      break;
    }
  }
231

232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
  // 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
  MantidVec::size_type lHalf = findHalfLoc(cenGausIn, height, GO_LEFT);
  // go forewards to get the half height on the otherside of the peak
  MantidVec::size_type rHalf = findHalfLoc(cenGausIn, height, GO_RIGHT);
  
  // the centre is the mean of the two end values
  double centreInd = static_cast<double>(lHalf + rHalf)/2.0;
  // convert the index into a time of flight value
  double tCalcu = m_tempWS->readX(0)[static_cast<int>(floor(centreInd))];
  // centreInd could be an integer or a half integer so return the value of the bin boundry or the mean of two bin boundaries
  tCalcu += m_tempWS->readX(0)[static_cast<int>(ceil(centreInd))];
  return tCalcu/2;
}
/** Calls CropWorkspace as a sub-algorithm and passes to it the InputWorkspace property
*  @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)
*  @throw out_of_range if start, end or specInd are set outside of the vaild 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
*/
void GetEi::extractSpec(int specInd, double start, double end)
{
  IAlgorithm_sptr childAlg =
257
    createSubAlgorithm("CropWorkspace", 100*m_fracCompl, 100*(m_fracCompl+CROP) );
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
  m_fracCompl += CROP;
  
  childAlg->setPropertyValue( "InputWorkspace",
                              getPropertyValue("InputWorkspace") );
  childAlg->setProperty( "XMin", start);
  childAlg->setProperty( "XMax", end);
  childAlg->setProperty( "StartWorkspaceIndex", specInd);
  childAlg->setProperty( "EndWorkspaceIndex", specInd);

  try
  {
    childAlg->execute();
  }
  catch (std::exception&)
  {
    g_log.error("Exception thrown while running CropWorkspace as a sub-algorithm");
    throw;
  }

  if ( ! childAlg->isExecuted() )
  {
    g_log.error("The CropWorkspace algorithm failed unexpectedly, aborting.");
    throw std::runtime_error(name() + " failed trying to run CropWorkspace");
  }
  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);
}

/** Finds the largest peak by looping through the histogram and finding the maximum
*  value 
* @param height this will became the peak height found by the fit
* @param centre will be set to the location of the peak center
* @throw invalid_argument if the peak is not much above the background
*/
void GetEi::getPeakEstimates(double &height, double &centre) const
{
  // take note of the number of background counts as error checking, do we have a peak or just a bump in the background
  double backgroundCounts = 0;
  // start at the first Y value
  height = m_tempWS->readY(0)[0];
  centre = m_tempWS->readX(0)[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(); ++i )
  {
    backgroundCounts += m_tempWS->readY(0)[i];
    if ( m_tempWS->readY(0)[i] > height )
    {
      height = m_tempWS->readY(0)[i];
      centre = m_tempWS->readX(0)[i];
    }
  }
  if ( height < PEAK_THRESH_H*backgroundCounts/m_tempWS->readY(0).size() )
  {
    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");
  }
  g_log.debug() << "Initial guess of peak position, based on the maximum Y value in the monitor spectrum, is at TOF " << centre << " (peak height " << height << " counts/microsecond)" << std::endl;
}
/** Gets the index of the bin that is closest to the bin given and contains a number of
*  counts less half of the number passed to this function, bin indexes start at zero
*  @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 go either GetEi::GO_LEFT or GetEi::GO_RIGHT
*  @return the index number of the first bin found where the counts are less than half
*  @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
*/
int GetEi::findHalfLoc(MantidVec::size_type startInd, double height, direction go) const
{
  MantidVec::size_type endInd = startInd;
  while ( m_tempWS->readY(0)[endInd] >  height/2.0 )
  {
    endInd += go;
    if ( endInd < 0 )
    {
      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");
    }
    if ( endInd >= m_tempWS->readY(0).size())
    {
      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");
    }
  }
  if ( std::abs(static_cast<int>(endInd - startInd)) < PEAK_THRESH_W )
  {// we didn't find a insignificant peak
    throw std::invalid_argument("No peak was found or its width is less than the threshold 2x" + boost::lexical_cast<std::string>(PEAK_THRESH_W) + " bins");
  }
  g_log.debug() << "One half height point found at TOF = " << m_tempWS->readX(0)[endInd] << " microseconds" << std::endl;
  return static_cast<int>(endInd);
}

/** Get the kinetic energy of a neuton in joules given it speed using E=mv^2/2
*  @param speed the instantanious speed of a neutron in metres per second
*  @return the energy in joules
*/
double GetEi::neutron_E_At(double speed) const
{
  // E_KE = mv^2/2
  return PhysicalConstants::NeutronMass*speed*speed/(2);
}

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

} // namespace Algorithms
} // namespace Mantid