GetEi2.cpp 21.6 KB
Newer Older
1
#include "MantidAlgorithms/GetEi2.h"
2

3
#include "MantidKernel/PhysicalConstants.h"
4
#include "MantidKernel/VectorHelper.h" 
5
#include "MantidAPI/WorkspaceValidators.h"
6
#include "MantidGeometry/Instrument/DetectorGroup.h"
7
#include "MantidGeometry/IObjComponent.h"
8
9

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

14
15
16
using namespace Mantid::Kernel;
using namespace Mantid::API;
using namespace Mantid::Geometry;
17
18
19
20
21
22
23
24

namespace Mantid
{
namespace Algorithms
{
  
  // Register the algorithm into the algorithm factory
  DECLARE_ALGORITHM(GetEi2)
25
26
27
28
29
30
31
32
  
  /// Sets documentation strings for this algorithm
  void GetEi2::initDocs()
  {
    this->setWikiSummary("Calculates the kinetic energy of neutrons leaving the source based on the time it takes for them to travel between two monitors. ");
    this->setOptionalMessage("Calculates the kinetic energy of neutrons leaving the source based on the time it takes for them to travel between two monitors.");
  }
  
33
34
35
36
37
38
39
40
41
42
43
44
45

/** 
* Default contructor
*/
GetEi2::GetEi2() : Algorithm(), m_input_ws(), m_peak1_pos(0, 0.0), m_fixedei(false), m_tof_window(0.1), m_peak_signif(2.0), m_peak_deriv(1.0),
  m_binwidth_frac(1.0/12.0), m_bkgd_frac(0.5)
{
  // Conversion factor common for converting between micro seconds and energy in meV
  m_t_to_mev = 5e11 * PhysicalConstants::NeutronMass / PhysicalConstants::meV;
}

void GetEi2::init()

46
{// Declare required input parameters for algorithm and do some validation here
47
48
49
50
  CompositeValidator<> *validator = new CompositeValidator<>;
  validator->add(new WorkspaceUnitValidator<>("TOF"));
  validator->add(new HistogramValidator<>);
  validator->add(new InstrumentValidator<>);
51

52
  declareProperty(new WorkspaceProperty<>(
53
    "InputWorkspace","",Direction::InOut, validator),
54
55
56
57
58
59
60
61
62
63
    "The X units of this workspace must be time of flight with times in\n"
    "microseconds");
  BoundedValidator<int> *mustBePositive = new BoundedValidator<int>();
  mustBePositive->setLower(0);
  declareProperty("Monitor1Spec", -1, mustBePositive,
    "The spectrum number to use as the first monitor\n");
  declareProperty("Monitor2Spec", -1, mustBePositive->clone(),
    "The spectrum number to use as the second monitor\n");
  BoundedValidator<double> *positiveDouble = new BoundedValidator<double>();
  positiveDouble->setLower(0.0);
64
65
  declareProperty("EnergyEstimate", -1.0 , positiveDouble,
    "An approximate value for the typical incident energy, energy of\n"
66
67
68
    "neutrons leaving the source (meV)");
  declareProperty("FixEi", false, "If true, the incident energy will be set to the value of the \n"
    "EnergyEstimate property.");
69
70
  declareProperty("IncidentEnergy", -1.0, Direction::Output);
  declareProperty("FirstMonitorPeak", -1.0, Direction::Output);
71
  declareProperty("FirstMonitorIndex", (size_t)0, Direction::Output);
72
  declareProperty("Tzero", EMPTY_DBL(), Direction::Output);
73

74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
}

/** Executes the algorithm
*  @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
*  @throw runtime_error if there is a problem with the SpectraDetectorMap or a sub-algorithm falls over
*/
void GetEi2::exec()
{
  m_input_ws = getProperty("InputWorkspace");
  m_fixedei = getProperty("FixEi");
  double initial_guess = getProperty("EnergyEstimate");
  double incident_energy = calculateEi(initial_guess);
  if( !m_fixedei )
  {
    g_log.notice() << "Incident energy = " << incident_energy << " meV from initial guess = " << initial_guess << " meV\n";
  }
  
  storeEi(incident_energy);
  
  setProperty("InputWorkspace", m_input_ws);
97
  // Output properties
98
  setProperty("IncidentEnergy", incident_energy);
99
  setProperty("FirstMonitorIndex", m_peak1_pos.first);
100
101
102
  setProperty("FirstMonitorPeak", m_peak1_pos.second);
}

103
/** Calculate the incident energy of the neutrons for the input workspace on this algorithm
104
 *  @param initial_guess :: A guess for value of the incident energy
105
 *  @return The calculated incident energy
106
107
108
 */
double GetEi2::calculateEi(const double initial_guess)
{
109
110
  const specid_t monitor1_spec = getProperty("Monitor1Spec");
  const specid_t monitor2_spec = getProperty("Monitor2Spec");
111
112

  //Covert spectrum numbers to workspace indices
Janik Zikovsky's avatar
Janik Zikovsky committed
113
  std::vector<specid_t> spec_nums(2, monitor1_spec);
114
  spec_nums[1] = monitor2_spec;
Janik Zikovsky's avatar
Janik Zikovsky committed
115
  std::vector<size_t> mon_indices;
116
117
  mon_indices.reserve(2);
  // get the index number of the histogram for the first monitor
118
  m_input_ws->getIndicesFromSpectra(spec_nums, mon_indices);
119
120
121
122
123
124
125
126
127
128
129
130

  if( mon_indices.size() != 2 )
  {
    g_log.error() << "Error retrieving monitor spectra from input workspace. Check input properties.\n";
    throw std::runtime_error("Error retrieving monitor spectra spectra from input workspace.");
  }

  // Calculate actual peak postion for each monitor peak
  double peak_times[2] = {0.0, 0.0};
  double det_distances[2] = {0.0, 0.0};
  for( unsigned int i = 0; i < 2; ++i )
  { 
Peterson, Peter's avatar
Peterson, Peter committed
131
    int64_t ws_index = mon_indices[i];
132
133
134
    det_distances[i] = getDistanceFromSource(ws_index);
    const double peak_guess = det_distances[i]*std::sqrt(m_t_to_mev/initial_guess);
    if( m_fixedei && i == 0 )
135
    {
136
137
138
      m_peak1_pos = std::make_pair(ws_index, peak_guess);
      g_log.information() << "First monitor peak = " << peak_guess << " microseconds from fixed Ei = " << initial_guess << " meV\n"; 
      break;
139
    }
140
141
    const double t_min = (1.0 - m_tof_window)*peak_guess;
    const double t_max = (1.0 + m_tof_window)*peak_guess;
142
    g_log.information() << "Time-of-flight window for peak " << (i+1) << ": tmin = " << t_min << " microseconds, tmax = " << t_max << " microseconds\n";
143
    peak_times[i] = calculatePeakPosition(ws_index, t_min, t_max);
144
    g_log.information() << "Peak for monitor " << (i+1) << " (at " << det_distances[i] << " metres) = " << peak_times[i] << " microseconds\n";
145
146
147
148
149
150
151
152
153
154
155
156
157
158
    if( i == 0 )
    {  
      //Store for later adjustment of bins
      m_peak1_pos = std::make_pair(ws_index, peak_times[i]);
    }
  }
  
  if( m_fixedei )
  {
    return initial_guess;
  }
  else
  {
    double mean_speed = (det_distances[1] - det_distances[0])/(peak_times[1] - peak_times[0]);
159
160
161
162
    
    double tzero = peak_times[1] - ((1.0/mean_speed)*det_distances[1]);
    setProperty("Tzero", tzero);
    
163
164
    g_log.debug() << "T0 = " << tzero << std::endl;
    
165
166
167
    g_log.debug() << "Mean Speed = " << mean_speed << std::endl;
    g_log.debug() << "Energy (meV) = " << mean_speed*mean_speed*m_t_to_mev << std::endl;
        
168
169
170
171
172
    return mean_speed*mean_speed*m_t_to_mev;
  }
}

/** Gets the distance between the source and detectors whose workspace index is passed
173
 *  @param ws_index :: The workspace index of the detector
174
175
176
 *  @return The distance between the source and the given detector(or DetectorGroup)
 *  @throw runtime_error if there is a problem
 */
Peterson, Peter's avatar
Peterson, Peter committed
177
double GetEi2::getDistanceFromSource(int64_t ws_index) const
178
{
179
180
  const IObjComponent_sptr source = m_input_ws->getInstrument()->getSource();
  // Retrieve a pointer detector
181
182
183
  IDetector_sptr det = m_input_ws->getDetector(ws_index);
  if( !det )
  {
184
185
186
    std::ostringstream msg;
	msg << "A detector for monitor at workspace index " << ws_index << " cannot be found. ";
    throw std::runtime_error(msg.str());
187
188
  }
  return det->getDistance(*source);
189
190
191
192
}

/**
 * Calculate the peak position of a given TOF range within the chosen spectrum
193
194
195
 * @param ws_index :: the wokspace index 
 * @param t_min :: the min time to consider
 * @param t_max :: the max time to consider
Nick Draper's avatar
re #100    
Nick Draper committed
196
 * @return the peak position
197
 */
Peterson, Peter's avatar
Peterson, Peter committed
198
double GetEi2::calculatePeakPosition(int64_t ws_index, double t_min, double t_max)
199
200
201
{
   //Crop out the current monitor workspace to the min/max times defined
  MatrixWorkspace_sptr monitor_ws = extractSpectrum(ws_index, t_min, t_max);
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
  // Workspace needs to be a count rate for the fitting algorithm
  WorkspaceHelpers::makeDistribution(monitor_ws);

  const double prominence(4.0);
  std::vector<double> dummyX, dummyY, dummyE;
  double peak_width = calculatePeakWidthAtHalfHeight(monitor_ws, prominence, dummyX, dummyY, dummyE);
  monitor_ws = rebin(monitor_ws, t_min, peak_width*m_binwidth_frac, t_max);

  double t_mean(0.0);
  try
  {
    t_mean = calculateFirstMoment(monitor_ws, prominence);
  }
  catch(std::runtime_error&)
  {
    //Retry with smaller prominence factor
    t_mean = calculateFirstMoment(monitor_ws, 2.0);
  }

221
222
223
  return t_mean;
}

224
/** Calls CropWorkspace as a sub-algorithm and passes to it the InputWorkspace property
225
226
227
 *  @param ws_index :: 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)
228
229
230
231
232
 *  @return The cropped workspace
 *  @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
 */
Peterson, Peter's avatar
Peterson, Peter committed
233
MatrixWorkspace_sptr GetEi2::extractSpectrum(int64_t ws_index, const double start, const double end)
234
235
{
  IAlgorithm_sptr childAlg = createSubAlgorithm("CropWorkspace");
236
  childAlg->setProperty("InputWorkspace", m_input_ws);
237
238
  childAlg->setProperty<int>("StartWorkspaceIndex", ws_index);
  childAlg->setProperty<int>("EndWorkspaceIndex", ws_index);
239
240
  childAlg->setProperty("XMin", start);
  childAlg->setProperty("XMax", end);
241
  childAlg->executeAsSubAlg();
242
  return childAlg->getProperty("OutputWorkspace");
243
244
245
246
}

/**
 * Calculate the width of the peak within the given region
247
248
249
250
251
 * @param data_ws :: The workspace containg the window around the peak
 * @param prominence :: The factor that the peak must be above the error to count as a peak
 * @param peak_x :: An output vector containing just the X values of the peak data
 * @param peak_y :: An output vector containing just the Y values of the peak data
 * @param peak_e :: An output vector containing just the E values of the peak data
252
253
 * @returns The width of the peak at half height
*/
254
double GetEi2::calculatePeakWidthAtHalfHeight(API::MatrixWorkspace_sptr data_ws, const double prominence, 
255
256
257
258
259
260
261
262
263
                                              std::vector<double> & peak_x, std::vector<double> & peak_y, std::vector<double> & peak_e) const
{
  // First create a temporary vector of bin_centre values to work with
  std::vector<double> Xs(data_ws->readX(0).size());
  VectorHelper::convertToBinCentre(data_ws->readX(0), Xs);
  const MantidVec & Ys = data_ws->readY(0);
  const MantidVec & Es = data_ws->readE(0);

  MantidVec::const_iterator peakIt = std::max_element(Ys.begin(), Ys.end());
264
  MantidVec::difference_type iPeak = peakIt - Ys.begin();
265
266
267
  double peakY = Ys[iPeak];
  double peakE = Es[iPeak];

268
  const std::vector<double>::size_type nxvals = Xs.size();
269
270

  //! Find data range that satisfies prominence criterion: im < ipk < ip will be nearest points that satisfy this
Peterson, Peter's avatar
Peterson, Peter committed
271
  int64_t im = static_cast<int64_t>(iPeak-1);
272
273
274
275
276
277
278
279
280
281
  for( ; im >= 0; --im )
  {
    const double ratio = Ys[im]/peakY;
    const double ratio_err = std::sqrt( std::pow(Es[im],2) + std::pow(ratio*peakE,2) )/peakY;
    if ( ratio < (1.0/prominence - m_peak_signif*ratio_err) )
    {
      break;
    }
  }
  
282
  std::vector<double>::size_type ip = iPeak+1;
283
284
285
286
287
288
289
290
  for( ; ip < nxvals; ip++ )
  {
    const double ratio = Ys[ip]/peakY;
    const double ratio_err =
      std::sqrt( std::pow(Es[ip], 2) + std::pow(ratio*peakE, 2) )/peakY;
    if ( ratio < (1.0/prominence - m_peak_signif*ratio_err) )
    {
      break;
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
  }
  
  if ( ip == nxvals || im < 0 )  
  {
    throw std::invalid_argument("No peak found in data that satisfies prominence criterion");
  }

  // We now have a peak, so can start filling output arguments
  // Determine extent of peak using derivatives
  // At this point 1 =< im < ipk < ip =< size(x)
  // After this section, new values will be given to im, ip that still satisfy these inequalities.
  // 
  // The algorithm for negative side skipped if im=1; positive side skipped if ip=size(x); 
  // if fails derivative criterion -> ip=size(x) (+ve)  im=1 (-ve)
  // In either case, we deem that the peak has a tail(s) that extend outside the range of x

  if ( ip < nxvals )
  {
    double deriv = -1000.0;
    double error = 0.0;
    while ( ( ip < nxvals - 1 ) && ( deriv < -m_peak_deriv*error ) )
    {
      double dtp = Xs[ip+1] - Xs[ip];
      double dtm = Xs[ip] - Xs[ip-1];
      deriv = 0.5*( ((Ys[ip+1] - Ys[ip]) / dtp) + ((Ys[ip] - Ys[ip-1]) / dtm) );
      error = 0.5*std::sqrt( ( (std::pow(Es[ip+1], 2) + std::pow(Es[ip], 2) ) / std::pow(dtp,2) ) + ((std::pow(Es[ip], 2) + std::pow(Es[ip-1], 2) )/std::pow(dtm,2) )
        - 2.0*(std::pow(Es[ip], 2) / (dtp*dtm)) );
      ip++;
    }
    ip--;
322

323
    if (deriv < -error)
324
    {
325
326
327
      ip = nxvals -1;      //        ! derivative criterion not met
    }
  }
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
  if (im > 0)
  {
    double deriv = 1000.0;
    double error = 0.0;
    while ( (im > 0) && (deriv > m_peak_deriv*error) )
    {
      double dtp = Xs[im+1] - Xs[im];
      double dtm = Xs[im] - Xs[im-1];
      deriv = 0.5*( ((Ys[im+1] - Ys[im]) / dtp) + ( (Ys[im] - Ys[im-1]) / dtm) );
      error = 0.5*std::sqrt( ( (std::pow(Es[im+1], 2) + std::pow(Es[im], 2) ) / std::pow(dtp, 2) ) + (( std::pow(Es[im], 2) + std::pow(Es[im-1], 2) ) / std::pow(dtm, 2) )
        - 2.0*std::pow(Es[im], 2)/(dtp*dtm) );
      im--;
    }
    im++;
    if (deriv > error) im = 0; //  derivative criterion not met
  }
  double pk_min = Xs[im];
  double pk_max = Xs[ip];
  double pk_width = Xs[ip] - Xs[im];

  // Determine background from either side of peak.
  // At this point, im and ip define the extreme points of the peak
  // Assume flat background
  double bkgd = 0.0;
  double bkgd_range = 0.0;
  double bkgd_min = std::max(Xs.front(), pk_min - m_bkgd_frac*pk_width);
  double bkgd_max = std::min(Xs.back(), pk_max + m_bkgd_frac*pk_width);

  if (im > 0)
  {
    double bkgd_m, bkgd_err_m;
    integrate(bkgd_m, bkgd_err_m, Xs, Ys, Es, bkgd_min, pk_min);
    bkgd = bkgd + bkgd_m;
    bkgd_range = bkgd_range + (pk_min - bkgd_min);
  }
364

365
366
367
368
369
370
371
  if (ip < nxvals - 1)
  {
    double bkgd_p, bkgd_err_p;
    integrate(bkgd_p, bkgd_err_p, Xs, Ys, Es, pk_max, bkgd_max);
    bkgd = bkgd + bkgd_p;
    bkgd_range = bkgd_range + (bkgd_max - pk_max);
  }
372

373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
  if ( im > 0  || ip < nxvals - 1 ) //background from at least one side
  {
    bkgd = bkgd / bkgd_range;
  }
  // Perform moment analysis on the peak after subtracting the background
  // Fill arrays with peak only:
  const std::vector<double>::size_type nvalues = ip - im + 1;
  peak_x.resize(nvalues);
  std::copy( Xs.begin() + im, Xs.begin() + ip + 1, peak_x.begin());
  peak_y.resize(nvalues);
  std::transform( Ys.begin() + im, Ys.begin() + ip + 1, peak_y.begin(), std::bind2nd(std::minus<double>(),bkgd) );
  peak_e.resize(nvalues);
  std::copy( Es.begin()+im, Es.begin() + ip + 1, peak_e.begin());

  // FWHH:
Peterson, Peter's avatar
Peterson, Peter committed
388
  int64_t ipk_int = static_cast<int64_t>(iPeak) - im;  //       ! peak position in internal array
389
  double hby2 = 0.5*peak_y[ipk_int];
Peterson, Peter's avatar
Peterson, Peter committed
390
  int64_t ip1(0), ip2(0);
391
392
  double xp_hh(0);

Peterson, Peter's avatar
Peterson, Peter committed
393
  int64_t nyvals = static_cast<int64_t>(peak_y.size());
394
395
  if (peak_y[nyvals-1] < hby2)
  {
Peterson, Peter's avatar
Peterson, Peter committed
396
    for( int64_t i = ipk_int; i < nyvals;  ++i )
397
398
    {
      if (peak_y[i] < hby2)
399
      {
400
401
402
        // after this point the intensity starts to go below half-height
        ip1 = i-1;       
        break;
403
      }
404
    }
Peterson, Peter's avatar
Peterson, Peter committed
405
    for ( int64_t i = nyvals-1; i >= ipk_int; --i )
406
407
    {
      if (peak_y[i]>hby2)
408
      {
409
410
        ip2 = i+1;           //   ! point closest to peak after which the intensity is always below half height
        break;
411
      }
412
413
414
415
416
417
418
419
    }
    xp_hh = peak_x[ip2] + (peak_x[ip1]-peak_x[ip2])*((hby2-peak_y[ip2])/(peak_y[ip1]-peak_y[ip2]));
  }
  else
  {
    xp_hh = peak_x[nyvals-1];
  }

Peterson, Peter's avatar
Peterson, Peter committed
420
  int64_t im1(0), im2(0);
421
422
423
  double xm_hh(0);
  if (peak_y[0]<hby2)
  {
Peterson, Peter's avatar
Peterson, Peter committed
424
    for( int64_t i = ipk_int; i >= 0; --i )
425
426
    {
      if (peak_y[i]<hby2)
427
      {
428
429
        im1 = i+1;   // ! after this point the intensity starts to go below half-height
        break;
430
      }
431
    }
Peterson, Peter's avatar
Peterson, Peter committed
432
    for ( int64_t i=0; i <= ipk_int; ++i )
433
434
    {
      if (peak_y[i]>hby2)
435
      {
436
437
        im2 = i-1;   // ! point closest to peak after which the intensity is always below half height
        break;
438
439
      }
    }
440
441
442
443
444
445
446
447
448
    xm_hh = peak_x[im2] + (peak_x[im1]-peak_x[im2])*((hby2-peak_y[im2])/(peak_y[im1]-peak_y[im2]));
  }
  else
  {
    xm_hh = peak_x.front();
  }
  return (xp_hh - xm_hh);
}

449
/** Calculate the first moment of the given workspace
450
451
 *  @param monitor_ws :: The workspace containing a single spectrum to work on
 *  @param prominence :: The factor over the background by which a peak is to be considered a "real" peak
452
 *  @return The calculated first moment
453
 */
454
double GetEi2::calculateFirstMoment(API::MatrixWorkspace_sptr monitor_ws, const double prominence)
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
{
  std::vector<double> peak_x, peak_y, peak_e;
  calculatePeakWidthAtHalfHeight(monitor_ws, prominence, peak_x, peak_y, peak_e);

  // Area
  double area(0.0), dummy(0.0);
  double pk_xmin = peak_x.front();
  double pk_xmax = peak_x.back();
  integrate(area, dummy, peak_x, peak_y, peak_e, pk_xmin, pk_xmax);
  // First moment
  std::transform(peak_y.begin(), peak_y.end(), peak_x.begin(), peak_y.begin(), std::multiplies<double>());
  double xbar(0.0);
  integrate(xbar, dummy, peak_x, peak_y, peak_e, pk_xmin, pk_xmax);
  return xbar / area;
}
470
471
472
473

/**
 * Rebin the workspace using the given bin widths 
 * @param monitor_ws
474
475
476
::  * @param first :: The minimum value for the new bin range
 * @param width :: The new bin width
 * @param end :: The maximum value for the new bin range
477
478
 * @returns The rebinned workspace
*/
479
API::MatrixWorkspace_sptr GetEi2::rebin(API::MatrixWorkspace_sptr monitor_ws, const double first, const double width, const double end)
480
481
482
483
484
485
{
  IAlgorithm_sptr childAlg = createSubAlgorithm("Rebin");
  childAlg->setProperty("InputWorkspace", monitor_ws);
  std::ostringstream binParams;
  binParams << first << "," << width << "," << end;
  childAlg->setPropertyValue( "Params", binParams.str());
486
  childAlg->executeAsSubAlg();
487
488
  return childAlg->getProperty("OutputWorkspace");
}
489
490
491

/**
 * Integrate a point data set
492
493
494
495
496
497
498
 * @param integral_val :: The result of the integral
 * @param integral_err :: The error value
 * @param x :: The X values
 * @param s :: The Y values
 * @param e :: The E values
 * @param xmin :: The minimum value for the integration
 * @param xmax :: The maximum value for the integration
499
 */
500
501
502
503
504
505
506
507
508
void GetEi2::integrate(double & integral_val, double &integral_err, const Mantid::MantidVec &x, const Mantid::MantidVec &s, const Mantid::MantidVec &e, 
                       const double xmin, const double xmax) const
{
  // MG: Note that this is integration of a point data set from libisis
  // @todo: Move to Kernel::VectorHelper and improve performance

  MantidVec::const_iterator lowit = std::lower_bound(x.begin(), x.end(), xmin);
  MantidVec::difference_type ml = std::distance(x.begin(),lowit);
  MantidVec::const_iterator highit = std::upper_bound(lowit,x.end(), xmax);
509
510
511
512
513
514
515
  MantidVec::difference_type mu = std::distance(x.begin(),highit);
  if( mu > 0 ) --mu;

  MantidVec::size_type nx(x.size());
  if( mu < ml )
  {
    //special case of no data points in the integration range
516
517
    unsigned int ilo = std::max<unsigned int>((unsigned int)ml - 1, 0);
    unsigned int ihi = std::min<unsigned int>((unsigned int)mu + 1, (unsigned int)nx);
518
519
520
521
522
    double fraction = (xmax - xmin)/(x[ihi] - x[ilo]);
    integral_val = 0.5 * fraction * 
      ( s[ihi]*((xmax - x[ilo]) + (xmin - x[ilo])) + s[ilo]*((x[ihi] - xmax)+(x[ihi] - xmin)) );
    double err_hi = e[ihi]*((xmax - x[ilo]) + (xmin - x[ilo]));
    double err_lo = e[ilo]*((x[ihi] - xmax)+(x[ihi] - xmin));
523
524
525
526
527
528
529
    integral_err = 0.5*fraction*std::sqrt( err_hi*err_hi  + err_lo*err_lo );
    return;
  }

  double x1eff(0.0), s1eff(0.0), e1eff(0.0);
  if( ml > 0 )
  {
530
531
    x1eff = (xmin*(xmin - x[ml-1]) + x[ml-1]*(x[ml] - xmin))/(x[ml] - x[ml-1]);
    double fraction = (x[ml]-xmin) / ((x[ml]-x[ml-1]) + (xmin - x[ml-1]));
532
533
534
535
536
537
538
539
    s1eff = s[ml-1]* fraction;
    e1eff = e[ml-1]* fraction;
  }
  else
  {
    x1eff = x[ml];
    s1eff = 0.0;
  }
540

541
  double xneff(0.0), sneff(0.0), eneff;
542
  if( mu < static_cast<int>(nx - 1))
543
  {
544
545
    xneff = (xmax*(x[mu+1]-xmax) + x[mu+1]*(xmax - x[mu]))/(x[mu+1] - x[mu]);
    const double fraction = (xmax - x[mu])/((x[mu+1] - x[mu]) + (x[mu+1]-xmax));
546
547
548
549
550
551
552
    sneff = s[mu+1]*fraction;
    eneff = e[mu+1]*fraction;
  }
  else
  {
    xneff = x.back();
    sneff = 0.0;
553
    eneff = 0.0; //Make sure to initialize to a value.
554
555
556
  }

  //xmin -> x[ml]
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
  integral_val = (x[ml] - x1eff)*(s[ml] + s1eff);
  integral_err = e1eff*(x[ml] - x1eff);
  integral_err *= integral_err;

  if( mu == ml )
  {
    double ierr = e[ml]*(xneff - x1eff);
    integral_err += ierr * ierr;
  }
  else if( mu == ml + 1 )
  {
    integral_val += (s[mu] + s[ml])*(x[mu] - x[ml]);
    double err_lo = e[ml]*(x[ml+1] - x1eff);
    double err_hi = e[mu]*(x[mu-1] - xneff);
    integral_err += err_lo*err_lo + err_hi*err_hi;
  }
  else
  {
575
    for( int i = (int)ml; i < mu; ++i )
576
    {
577
578
      integral_val += (s[i+1] + s[i])*(x[i+1] - x[i]);
      if( i < mu - 1 )
579
      {
580
581
        double ierr = e[i+1]*(x[i+2] - x[i]);
        integral_err += ierr * ierr;
582
583
      }
    }
584
585
586
587
588
589
590
591
592
593
594
  }

  //x[mu] -> xmax
  integral_val += (xneff - x[mu])*(s[mu] + sneff);
  double err_tmp = eneff*(xneff - x[mu]);
  integral_err += err_tmp*err_tmp;

  integral_val *= 0.5;
  integral_err = 0.5*sqrt(integral_err);
}

595
596
/**
 * Store the value of Ei wihin the log data on the Sample object
597
 * @param ei :: The value of the incident energy of the neutron
598
599
600
601
 */
void GetEi2::storeEi(const double ei) const
{
  Property *incident_energy = new PropertyWithValue<double>("Ei",ei,Direction::Input);
602
  m_input_ws->mutableRun().addProperty(incident_energy, true);
603
604
605
}

}
606
}