Mantid
Loading...
Searching...
No Matches
GetEi2.cpp
Go to the documentation of this file.
1// Mantid Repository : https://github.com/mantidproject/mantid
2//
3// Copyright © 2018 ISIS Rutherford Appleton Laboratory UKRI,
4// NScD Oak Ridge National Laboratory, European Spallation Source,
5// Institut Laue - Langevin & CSNS, Institute of High Energy Physics, CAS
6// SPDX - License - Identifier: GPL - 3.0 +
8
12#include "MantidAPI/Run.h"
22
23#include <algorithm>
24#include <boost/lexical_cast.hpp>
25#include <cmath>
26#include <sstream>
27#include <utility>
28
29using namespace Mantid::Kernel;
30using namespace Mantid::API;
31using namespace Mantid::Geometry;
32
33namespace Mantid::Algorithms {
34
35// Register the algorithm into the algorithm factory
37
38
42 : Algorithm(), m_input_ws(), m_peak1_pos(0, 0.0), m_fixedei(false), m_tof_window(0.1), m_peak_signif(2.0),
43 m_peak_deriv(1.0), m_binwidth_frac(1.0 / 12.0), m_bkgd_frac(0.5) {
44 // Conversion factor common for converting between micro seconds and energy in
45 // meV
47}
48
50
51{ // Declare required input parameters for algorithm and do some validation here
52 auto validator = std::make_shared<CompositeValidator>();
53 validator->add<WorkspaceUnitValidator>("TOF");
54 validator->add<HistogramValidator>();
55 validator->add<InstrumentValidator>();
56
57 declareProperty(std::make_unique<WorkspaceProperty<>>("InputWorkspace", "", Direction::InOut, validator),
58 "The X units of this workspace must be time of flight with times in\n"
59 "microseconds");
60 auto mustBePositive = std::make_shared<BoundedValidator<int>>();
61 mustBePositive->setLower(0);
62 declareProperty("Monitor1Spec", EMPTY_INT(), mustBePositive,
63 "The spectrum number of the output of the first monitor, "
64 "e.g. MAPS 41474, MARI 2, MERLIN 69634.\n If empty, it will "
65 "be read from the instrument file.\n");
66 declareProperty("Monitor2Spec", EMPTY_INT(), mustBePositive,
67 "The spectrum number of the output of the second monitor "
68 "e.g. MAPS 41475, MARI 3, MERLIN 69638.\n If empty, it will "
69 "be read from the instrument file.\n");
70 auto positiveDouble = std::make_shared<BoundedValidator<double>>();
71 positiveDouble->setLower(0.0);
72 declareProperty("EnergyEstimate", EMPTY_DBL(), positiveDouble,
73 "An approximate value for the typical incident energy, energy of\n"
74 "neutrons leaving the source (meV)\n Can be empty if there is a Sample "
75 "Log called EnergyRequest. Otherwise it is mandatory.");
76 declareProperty("FixEi", false,
77 "If true, the incident energy will be set to the value of the \n"
78 "EnergyEstimate property.");
79
80 declareProperty("IncidentEnergy", -1.0, "The energy of neutron in meV, it is also printed to the Mantid's log",
82
83 declareProperty("FirstMonitorPeak", -1.0,
84 "The time in :math:`\\rm{\\mu s}` when the count rate of the "
85 "first monitor, which defaults to the last monitor the beam "
86 "hits before the sample, is greatest. It is the mean X value "
87 "for the bin with the highest number of counts per second "
88 "and is also writen to Mantid's log.",
90
91 declareProperty("FirstMonitorIndex", 0, "The workspace index of the first monitor in the input workspace.",
93
94 declareProperty("Tzero", 0.0, "", Direction::Output);
95
96 auto inRange0toOne = std::make_shared<BoundedValidator<double>>();
97 inRange0toOne->setLower(0.0);
98 inRange0toOne->setUpper(1.0);
99 declareProperty("PeakSearchRange", 0.1, inRange0toOne,
100 "Specifies the relative TOF range where the algorithm tries "
101 "to find the monitor peak. Search occurs within "
102 "PEAK_TOF_Guess * (1 +/- PeakSearchRange) ranges.\n"
103 "Defaults are almost always sufficient but decrease this "
104 "value for very narrow peaks and increase for wide.",
106}
107
118 m_input_ws = getProperty("InputWorkspace");
119 m_fixedei = getProperty("FixEi");
120 m_tof_window = getProperty("PeakSearchRange");
121
122 double initial_guess = getProperty("EnergyEstimate");
123 // check if incident energy guess is left empty, and try to find it as
124 // EnergyRequest parameter
125 if (initial_guess == EMPTY_DBL()) {
126 if (m_input_ws->run().hasProperty("EnergyRequest")) {
127 initial_guess = m_input_ws->getLogAsSingleValue("EnergyRequest");
128 } else {
129 throw std::invalid_argument("Could not find an energy guess");
130 }
131 }
132 double incident_energy = calculateEi(initial_guess);
133
134 storeEi(incident_energy);
135
136 setProperty("InputWorkspace", m_input_ws);
137 // Output properties
138 setProperty("IncidentEnergy", incident_energy);
139 setProperty("FirstMonitorIndex", m_peak1_pos.first);
140 setProperty("FirstMonitorPeak", m_peak1_pos.second);
141}
142
148double GetEi2::calculateEi(const double initial_guess) {
149 const specnum_t monitor1_spec = getProperty("Monitor1Spec");
150 const specnum_t monitor2_spec = getProperty("Monitor2Spec");
151 specnum_t mon1 = monitor1_spec, mon2 = monitor2_spec;
152
153 const ParameterMap &pmap = m_input_ws->constInstrumentParameters();
154 Instrument_const_sptr instrument = m_input_ws->getInstrument();
155
156 // If we have a t0 formula, calculate t0 and return the initial guess
157 Parameter_sptr par = pmap.getRecursive(instrument->getChild(0).get(), "t0_formula");
158 if (par) {
159 std::string formula = par->asString();
160 std::stringstream guess;
161 guess << initial_guess;
162
163 while (formula.find("incidentEnergy") != std::string::npos) {
164 // check if more than one 'value' in m_eq
165 size_t found = formula.find("incidentEnergy");
166 formula.replace(found, 14, guess.str());
167 }
168
169 try {
170 mu::Parser p;
171 p.SetExpr(formula);
172 double tzero = p.Eval();
173 setProperty("Tzero", tzero);
174
175 g_log.debug() << "T0 = " << tzero << '\n';
176 return initial_guess;
177 } catch (mu::Parser::exception_type &e) {
179 "Equation attribute for t0 formula. Muparser error message is: " + e.GetMsg());
180 }
181 }
182
183 if (mon1 == EMPTY_INT()) {
184 par = pmap.getRecursive(instrument->getChild(0).get(), "ei-mon1-spec");
185 if (par) {
186 mon1 = boost::lexical_cast<specnum_t>(par->asString());
187 }
188 }
189 if (mon2 == EMPTY_INT()) {
190 par = pmap.getRecursive(instrument->getChild(0).get(), "ei-mon2-spec");
191 if (par) {
192 mon2 = boost::lexical_cast<specnum_t>(par->asString());
193 }
194 }
195 if ((mon1 == EMPTY_INT()) || (mon2 == EMPTY_INT())) {
196 throw std::invalid_argument("Could not determine spectrum number to use. Try to set it explicitly");
197 }
198 // Covert spectrum numbers to workspace indices
199 std::vector<specnum_t> spec_nums(2, mon1);
200 spec_nums[1] = mon2;
201 // get the index number of the histogram for the first monitor
202 auto mon_indices = m_input_ws->getIndicesFromSpectra(spec_nums);
203
204 if (mon_indices.size() != 2) {
205 g_log.error() << "Error retrieving monitor spectra from input workspace. "
206 "Check input properties.\n";
207 throw std::runtime_error("Error retrieving monitor spectra spectra from input workspace.");
208 }
209
210 // Calculate actual peak postion for each monitor peak
211 double peak_times[2] = {0.0, 0.0};
212 double det_distances[2] = {0.0, 0.0};
213 auto &spectrumInfo = m_input_ws->spectrumInfo();
214 for (unsigned int i = 0; i < 2; ++i) {
215 size_t ws_index = mon_indices[i];
216 det_distances[i] = getDistanceFromSource(ws_index, spectrumInfo);
217 const double peak_guess = det_distances[i] * std::sqrt(m_t_to_mev / initial_guess);
218 const double t_min = (1.0 - m_tof_window) * peak_guess;
219 const double t_max = (1.0 + m_tof_window) * peak_guess;
220 g_log.information() << "Time-of-flight window for peak " << (i + 1) << ": tmin = " << t_min
221 << " microseconds, tmax = " << t_max << " microseconds\n";
222 try {
223 peak_times[i] = calculatePeakPosition(ws_index, t_min, t_max);
224 g_log.information() << "Peak for monitor " << (i + 1) << " (at " << det_distances[i]
225 << " metres) = " << peak_times[i] << " microseconds\n";
226 } catch (std::invalid_argument &) {
227
228 if (!m_fixedei) {
229 throw std::invalid_argument("No peak found for the monitor with spectra num: " + std::to_string(spec_nums[i]) +
230 " (at " + boost::lexical_cast<std::string>(det_distances[i]) +
231 " metres from source).\n");
232 } else {
233 peak_times[i] = peak_guess;
234 g_log.warning() << "No peak found for monitor with spectra num " << spec_nums[i] << " (at " << det_distances[i]
235 << " metres).\n";
236 g_log.warning() << "Using guess time found from energy estimate of Peak = " << peak_times[i]
237 << " microseconds\n";
238 }
239 }
240 if (i == 0) {
241 m_peak1_pos = std::make_pair(static_cast<int>(ws_index), peak_times[i]);
242 if (m_fixedei)
243 break;
244 }
245 }
246
247 if (m_fixedei) {
248 g_log.notice() << "Incident energy fixed at Ei=" << initial_guess << " meV\n";
249 return initial_guess;
250 } else {
251 double mean_speed = (det_distances[1] - det_distances[0]) / (peak_times[1] - peak_times[0]);
252 double tzero = peak_times[1] - ((1.0 / mean_speed) * det_distances[1]);
253 g_log.debug() << "T0 = " << tzero << '\n';
254 g_log.debug() << "Mean Speed = " << mean_speed << '\n';
255 setProperty("Tzero", tzero);
256
257 const double energy = mean_speed * mean_speed * m_t_to_mev;
258 g_log.notice() << "Incident energy calculated at Ei= " << energy << " meV from initial guess = " << initial_guess
259 << " meV\n";
260 return energy;
261 }
262}
263
273double GetEi2::getDistanceFromSource(size_t ws_index, const SpectrumInfo &spectrumInfo) const {
274 g_log.debug() << "Computing distance between spectrum at index '" << ws_index << "' and the source\n";
275
276 const auto &detector = spectrumInfo.detector(ws_index);
277 const IComponent_const_sptr source = m_input_ws->getInstrument()->getSource();
278 if (!spectrumInfo.hasDetectors(ws_index)) {
279 std::ostringstream msg;
280 msg << "A detector for monitor at workspace index " << ws_index << " cannot be found. ";
281 throw std::runtime_error(msg.str());
282 }
283 if (g_log.is(Logger::Priority::PRIO_DEBUG)) {
284 g_log.debug() << "Detector position = " << spectrumInfo.position(ws_index)
285 << ", Source position = " << source->getPos() << "\n";
286 }
287 const double dist = detector.getDistance(*source);
288 g_log.debug() << "Distance = " << dist << " metres\n";
289 return dist;
290}
291
299double GetEi2::calculatePeakPosition(size_t ws_index, double t_min, double t_max) {
300 // Crop out the current monitor workspace to the min/max times defined
301 MatrixWorkspace_sptr monitor_ws = extractSpectrum(ws_index, t_min, t_max);
302 // Workspace needs to be a count rate for the fitting algorithm
304
305 const double prominence(4.0);
306 std::vector<double> dummyX, dummyY, dummyE;
307 double peak_width = calculatePeakWidthAtHalfHeight(monitor_ws, prominence, dummyX, dummyY, dummyE);
308 monitor_ws = rebin(monitor_ws, t_min, peak_width * m_binwidth_frac, t_max);
309
310 double t_mean(0.0);
311 try {
312 t_mean = calculateFirstMoment(monitor_ws, prominence);
313 } catch (std::runtime_error &) {
314 // Retry with smaller prominence factor
315 t_mean = calculateFirstMoment(monitor_ws, 2.0);
316 }
317
318 return t_mean;
319}
320
334MatrixWorkspace_sptr GetEi2::extractSpectrum(size_t ws_index, const double start, const double end) {
335 auto childAlg = createChildAlgorithm("CropWorkspace");
336 childAlg->setProperty("InputWorkspace", m_input_ws);
337 childAlg->setProperty<int>("StartWorkspaceIndex", static_cast<int>(ws_index));
338 childAlg->setProperty<int>("EndWorkspaceIndex", static_cast<int>(ws_index));
339 childAlg->setProperty("XMin", start);
340 childAlg->setProperty("XMax", end);
341 childAlg->executeAsChildAlg();
342 MatrixWorkspace_sptr monitors = childAlg->getProperty("OutputWorkspace");
343 if (std::dynamic_pointer_cast<API::IEventWorkspace>(m_input_ws)) {
344 // Convert to a MatrixWorkspace representation
345 childAlg = createChildAlgorithm("ConvertToMatrixWorkspace");
346 childAlg->setProperty("InputWorkspace", monitors);
347 childAlg->setProperty("OutputWorkspace", monitors);
348 childAlg->executeAsChildAlg();
349 monitors = childAlg->getProperty("OutputWorkspace");
350 }
351
352 return monitors;
353}
354
368double GetEi2::calculatePeakWidthAtHalfHeight(const API::MatrixWorkspace_sptr &data_ws, const double prominence,
369 std::vector<double> &peak_x, std::vector<double> &peak_y,
370 std::vector<double> &peak_e) const {
371 // Use WS->points() to create a temporary vector of bin_centre values to work
372 // with
373 auto Xs = data_ws->points(0);
374 const auto &Ys = data_ws->y(0);
375 const auto &Es = data_ws->e(0);
376
377 auto peakIt = std::max_element(Ys.cbegin(), Ys.cend());
378 double bkg_val = *std::min_element(Ys.begin(), Ys.end());
379 if (*peakIt == bkg_val) {
380 throw std::invalid_argument("No peak in the range specified as minimal and "
381 "maximal values of the function are equal ");
382 }
383 auto iPeak = peakIt - Ys.begin();
384 double peakY = Ys[iPeak] - bkg_val;
385 double peakE = Es[iPeak];
386
387 const std::vector<double>::size_type nxvals = Xs.size();
388
390 // nearest points that satisfy this
391 auto im = static_cast<int64_t>(iPeak - 1);
392 for (; im >= 0; --im) {
393 const double ratio = (Ys[im] - bkg_val) / peakY;
394 const double ratio_err = std::sqrt(std::pow(Es[im], 2) + std::pow(ratio * peakE, 2)) / peakY;
395 if (ratio < (1.0 / prominence - m_peak_signif * ratio_err)) {
396 break;
397 }
398 }
399
400 std::vector<double>::size_type ip = iPeak + 1;
401 for (; ip < nxvals; ip++) {
402 const double ratio = (Ys[ip] - bkg_val) / peakY;
403 const double ratio_err = std::sqrt(std::pow(Es[ip], 2) + std::pow(ratio * peakE, 2)) / peakY;
404 if (ratio < (1.0 / prominence - m_peak_signif * ratio_err)) {
405 break;
406 }
407 }
408
409 if (ip == nxvals || im < 0) {
410 throw std::invalid_argument("No peak found in data that satisfies prominence criterion");
411 }
412
413 // We now have a peak, so can start filling output arguments
414 // Determine extent of peak using derivatives
415 // At this point 1 =< im < ipk < ip =< size(x)
416 // After this section, new values will be given to im, ip that still satisfy
417 // these inequalities.
418 //
419 // The algorithm for negative side skipped if im=1; positive side skipped if
420 // ip=size(x);
421 // if fails derivative criterion -> ip=size(x) (+ve) im=1 (-ve)
422 // In either case, we deem that the peak has a tail(s) that extend outside the
423 // range of x
424
425 if (ip < nxvals) {
426 double deriv = -1000.0;
427 double error = 0.0;
428 while ((ip < nxvals - 1) && (deriv < -m_peak_deriv * error)) {
429 double dtp = Xs[ip + 1] - Xs[ip];
430 double dtm = Xs[ip] - Xs[ip - 1];
431 deriv = 0.5 * (((Ys[ip + 1] - Ys[ip]) / dtp) + ((Ys[ip] - Ys[ip - 1]) / dtm));
432 error = 0.5 * std::sqrt(((std::pow(Es[ip + 1], 2) + std::pow(Es[ip], 2)) / std::pow(dtp, 2)) +
433 ((std::pow(Es[ip], 2) + std::pow(Es[ip - 1], 2)) / std::pow(dtm, 2)) -
434 2.0 * (std::pow(Es[ip], 2) / (dtp * dtm)));
435 ip++;
436 }
437 ip--;
438
439 if (deriv < -error) {
440 ip = nxvals - 1; // ! derivative criterion not met
441 }
442 }
443
444 if (im > 0) {
445 double deriv = 1000.0;
446 double error = 0.0;
447 while ((im > 0) && (deriv > m_peak_deriv * error)) {
448 double dtp = Xs[im + 1] - Xs[im];
449 double dtm = Xs[im] - Xs[im - 1];
450 deriv = 0.5 * (((Ys[im + 1] - Ys[im]) / dtp) + ((Ys[im] - Ys[im - 1]) / dtm));
451 error = 0.5 * std::sqrt(((std::pow(Es[im + 1], 2) + std::pow(Es[im], 2)) / std::pow(dtp, 2)) +
452 ((std::pow(Es[im], 2) + std::pow(Es[im - 1], 2)) / std::pow(dtm, 2)) -
453 2.0 * std::pow(Es[im], 2) / (dtp * dtm));
454 im--;
455 }
456 im++;
457 if (deriv > error)
458 im = 0; // derivative criterion not met
459 }
460 double pk_min = Xs[im];
461 double pk_max = Xs[ip];
462 double pk_width = Xs[ip] - Xs[im];
463
464 // Determine background from either side of peak.
465 // At this point, im and ip define the extreme points of the peak
466 // Assume flat background
467 double bkgd = 0.0;
468 double bkgd_range = 0.0;
469 double bkgd_min = std::max(Xs.front(), pk_min - m_bkgd_frac * pk_width);
470 double bkgd_max = std::min(Xs.back(), pk_max + m_bkgd_frac * pk_width);
471
472 if (im > 0) {
473 double bkgd_m, bkgd_err_m;
474 integrate(bkgd_m, bkgd_err_m, Xs.rawData(), Ys, Es, bkgd_min, pk_min);
475 bkgd = bkgd + bkgd_m;
476 bkgd_range = bkgd_range + (pk_min - bkgd_min);
477 }
478
479 if (ip < nxvals - 1) {
480 double bkgd_p, bkgd_err_p;
481 integrate(bkgd_p, bkgd_err_p, Xs.rawData(), Ys, Es, pk_max, bkgd_max);
482 bkgd = bkgd + bkgd_p;
483 bkgd_range = bkgd_range + (bkgd_max - pk_max);
484 }
485
486 if (im > 0 || ip < nxvals - 1) // background from at least one side
487 {
488 bkgd = bkgd / bkgd_range;
489 }
490 // Perform moment analysis on the peak after subtracting the background
491 // Fill arrays with peak only:
492 const std::vector<double>::size_type nvalues = ip - im + 1;
493 peak_x.resize(nvalues);
494 std::copy(Xs.begin() + im, Xs.begin() + ip + 1, peak_x.begin());
495 peak_y.resize(nvalues);
496 using std::placeholders::_1;
497 std::transform(Ys.begin() + im, Ys.begin() + ip + 1, peak_y.begin(), std::bind(std::minus<double>(), _1, bkgd));
498 peak_e.resize(nvalues);
499 std::copy(Es.begin() + im, Es.begin() + ip + 1, peak_e.begin());
500
501 // FWHH:
502 int64_t ipk_int = static_cast<int64_t>(iPeak) - im; // ! peak position in internal array
503 double hby2 = 0.5 * peak_y[ipk_int];
504 double xp_hh(0);
505
506 auto nyvals = static_cast<int64_t>(peak_y.size());
507 if (peak_y[nyvals - 1] < hby2) {
508 int64_t ip1(0), ip2(0);
509 for (int64_t i = ipk_int; i < nyvals; ++i) {
510 if (peak_y[i] < hby2) {
511 // after this point the intensity starts to go below half-height
512 ip1 = i - 1;
513 break;
514 }
515 }
516 for (int64_t i = nyvals - 1; i >= ipk_int; --i) {
517 if (peak_y[i] > hby2) {
518 ip2 = i + 1; // ! point closest to peak after which the intensity is
519 // always below half height
520 break;
521 }
522 }
523 // A broad peak with many local maxima on the side can cause the algorithm to give the same indices for the two
524 // points (ip1 and ip2) either side of the half-width point. Noise may also result in equal y values for different
525 // ip1 and ip2 points which would cause a divide-by-zero error in the xp_hh calculation below. We know the
526 // algorithm isn't perfect so move one index until the y values are differnt.
527 if (peak_y[ip1] == peak_y[ip2]) {
528 g_log.warning() << "A peak with a local maxima on the trailing edge has "
529 "been found. The estimation of the "
530 << "half-height point will not be as accurate.\n";
531 // Shift the index ip1 until the heights are no longer equal. This avoids the divide-by-zero.
532 while (peak_y[ip1] == peak_y[ip2]) {
533 ip1--;
534 if (ip1 < ipk_int)
535 throw std::invalid_argument("Trailing edge values are equal until up to the peak centre.");
536 }
537 }
538
539 xp_hh = peak_x[ip2] + (peak_x[ip1] - peak_x[ip2]) * ((hby2 - peak_y[ip2]) / (peak_y[ip1] - peak_y[ip2]));
540 } else {
541 xp_hh = peak_x[nyvals - 1];
542 }
543
544 double xm_hh(0);
545 if (peak_y[0] < hby2) {
546 int64_t im1(0), im2(0);
547 for (int64_t i = ipk_int; i >= 0; --i) {
548 if (peak_y[i] < hby2) {
549 im1 = i + 1; // ! before this point the intensity starts to go below
550 // half-height
551 break;
552 }
553 }
554 for (int64_t i = 0; i <= ipk_int; ++i) {
555 if (peak_y[i] > hby2) {
556 im2 = i - 1; // ! point closest to peak before which the intensity is
557 // always below half height
558 break;
559 }
560 }
561 // A broad peak with many local maxima on the side can cause the algorithm to give the same indices for the two
562 // points (im1 and im2) either side of the half-width point. Noise may also result in equal y values for different
563 // im1 and im2 points which would cause a divide-by-zero error in the xp_hh calculation below. We know the
564 // algorithm isn't perfect so move one index until the y values are differnt.
565 if (peak_y[im1] == peak_y[im2]) {
566 g_log.warning() << "A peak with a local maxima on the rising edge has "
567 "been found. The estimation of the "
568 << "half-height point will not be as accurate.\n";
569 // Shift the index ip1 until the heights are no longer equal. This avoids the divide-by-zero.
570 while (peak_y[im1] == peak_y[im2]) {
571 im1++;
572 if (im1 >= ipk_int)
573 throw std::invalid_argument("The rising edge values are equal up to the peak centre.");
574 }
575 }
576
577 xm_hh = peak_x[im2] + (peak_x[im1] - peak_x[im2]) * ((hby2 - peak_y[im2]) / (peak_y[im1] - peak_y[im2]));
578 } else {
579 xm_hh = peak_x.front();
580 }
581 return (xp_hh - xm_hh);
582}
583
590double GetEi2::calculateFirstMoment(const API::MatrixWorkspace_sptr &monitor_ws, const double prominence) {
591 std::vector<double> peak_x, peak_y, peak_e;
592 calculatePeakWidthAtHalfHeight(monitor_ws, prominence, peak_x, peak_y, peak_e);
593
594 // Area
595 double area(0.0), dummy(0.0);
596 double pk_xmin = peak_x.front();
597 double pk_xmax = peak_x.back();
598 integrate(area, dummy, peak_x, peak_y, peak_e, pk_xmin, pk_xmax);
599 // First moment
600 std::transform(peak_y.begin(), peak_y.end(), peak_x.begin(), peak_y.begin(), std::multiplies<double>());
601 double xbar(0.0);
602 integrate(xbar, dummy, peak_x, peak_y, peak_e, pk_xmin, pk_xmax);
603 return xbar / area;
604}
605
615 const double width, const double end) {
616 auto childAlg = createChildAlgorithm("Rebin");
617 childAlg->setProperty("InputWorkspace", monitor_ws);
618 std::ostringstream binParams;
619 binParams << first << "," << width << "," << end;
620 childAlg->setPropertyValue("Params", binParams.str());
621 childAlg->setProperty("IgnoreBinErrors", true);
622 childAlg->executeAsChildAlg();
623 return childAlg->getProperty("OutputWorkspace");
624}
625
636void GetEi2::integrate(double &integral_val, double &integral_err, const HistogramData::HistogramX &x,
637 const HistogramData::HistogramY &s, const HistogramData::HistogramE &e, const double xmin,
638 const double xmax) const {
639 // MG: Note that this is integration of a point data set from libisis
640 // @todo: Move to Kernel::VectorHelper and improve performance
641
642 auto lowit = std::lower_bound(x.cbegin(), x.cend(), xmin);
643 auto ml = std::distance(x.cbegin(), lowit);
644 auto highit = std::upper_bound(lowit, x.cend(), xmax);
645 auto mu = std::distance(x.cbegin(), highit);
646 if (mu > 0)
647 --mu;
648
649 auto nx(x.size());
650 if (mu < ml) {
651 // special case of no data points in the integration range
652 unsigned int ilo = std::max<unsigned int>(static_cast<unsigned int>(ml) - 1, 0);
653 unsigned int ihi = std::min<unsigned int>(static_cast<unsigned int>(mu) + 1, static_cast<unsigned int>(nx));
654 double fraction = (xmax - xmin) / (x[ihi] - x[ilo]);
655 integral_val =
656 0.5 * fraction * (s[ihi] * ((xmax - x[ilo]) + (xmin - x[ilo])) + s[ilo] * ((x[ihi] - xmax) + (x[ihi] - xmin)));
657 double err_hi = e[ihi] * ((xmax - x[ilo]) + (xmin - x[ilo]));
658 double err_lo = e[ilo] * ((x[ihi] - xmax) + (x[ihi] - xmin));
659 integral_err = 0.5 * fraction * std::sqrt(err_hi * err_hi + err_lo * err_lo);
660 return;
661 }
662
663 double x1eff(0.0), s1eff(0.0), e1eff(0.0);
664 if (ml > 0) {
665 x1eff = (xmin * (xmin - x[ml - 1]) + x[ml - 1] * (x[ml] - xmin)) / (x[ml] - x[ml - 1]);
666 double fraction = (x[ml] - xmin) / ((x[ml] - x[ml - 1]) + (xmin - x[ml - 1]));
667 s1eff = s[ml - 1] * fraction;
668 e1eff = e[ml - 1] * fraction;
669 } else {
670 x1eff = x[ml];
671 s1eff = 0.0;
672 }
673
674 double xneff(0.0), sneff(0.0), eneff;
675 if (mu < static_cast<int>(nx - 1)) {
676 xneff = (xmax * (x[mu + 1] - xmax) + x[mu + 1] * (xmax - x[mu])) / (x[mu + 1] - x[mu]);
677 const double fraction = (xmax - x[mu]) / ((x[mu + 1] - x[mu]) + (x[mu + 1] - xmax));
678 sneff = s[mu + 1] * fraction;
679 eneff = e[mu + 1] * fraction;
680 } else {
681 xneff = x.back();
682 sneff = 0.0;
683 eneff = 0.0; // Make sure to initialize to a value.
684 }
685
686 // xmin -> x[ml]
687 integral_val = (x[ml] - x1eff) * (s[ml] + s1eff);
688 integral_err = e1eff * (x[ml] - x1eff);
689 integral_err *= integral_err;
690
691 if (mu == ml) {
692 double ierr = e[ml] * (xneff - x1eff);
693 integral_err += ierr * ierr;
694 } else if (mu == ml + 1) {
695 integral_val += (s[mu] + s[ml]) * (x[mu] - x[ml]);
696 double err_lo = e[ml] * (x[ml + 1] - x1eff);
697 double err_hi = e[mu] * (x[mu - 1] - xneff);
698 integral_err += err_lo * err_lo + err_hi * err_hi;
699 } else {
700 for (auto i = static_cast<int>(ml); i < mu; ++i) {
701 integral_val += (s[i + 1] + s[i]) * (x[i + 1] - x[i]);
702 if (i < mu - 1) {
703 double ierr = e[i + 1] * (x[i + 2] - x[i]);
704 integral_err += ierr * ierr;
705 }
706 }
707 }
708
709 // x[mu] -> xmax
710 integral_val += (xneff - x[mu]) * (s[mu] + sneff);
711 double err_tmp = eneff * (xneff - x[mu]);
712 integral_err += err_tmp * err_tmp;
713
714 integral_val *= 0.5;
715 integral_err = 0.5 * sqrt(integral_err);
716}
717
722void GetEi2::storeEi(const double ei) const {
723 Property *incident_energy = new PropertyWithValue<double>("Ei", ei, Direction::Input);
724 m_input_ws->mutableRun().addProperty(incident_energy, true);
725}
726} // namespace Mantid::Algorithms
#define DECLARE_ALGORITHM(classname)
Definition: Algorithm.h:576
double energy
Definition: GetAllEi.cpp:157
double error
Definition: IndexPeaks.cpp:133
Base class from which all concrete algorithm classes should be derived.
Definition: Algorithm.h:85
void declareProperty(std::unique_ptr< Kernel::Property > p, const std::string &doc="") override
Add a property to the list of managed properties.
Definition: Algorithm.cpp:1913
TypedValue getProperty(const std::string &name) const override
Get the value of a property.
Definition: Algorithm.cpp:2076
virtual std::shared_ptr< Algorithm > createChildAlgorithm(const std::string &name, const double startProgress=-1., const double endProgress=-1., const bool enableLogging=true, const int &version=-1)
Create a Child Algorithm.
Definition: Algorithm.cpp:842
Kernel::Logger & g_log
Definition: Algorithm.h:451
A validator which checks that a workspace contains histogram data (the default) or point data as requ...
A validator which checks that a workspace has a valid instrument.
API::SpectrumInfo is an intermediate step towards a SpectrumInfo that is part of Instrument-2....
Definition: SpectrumInfo.h:53
bool hasDetectors(const size_t index) const
Returns true if the spectrum is associated with detectors in the instrument.
Kernel::V3D position(const size_t index) const
Returns the position of the spectrum with given index.
const Geometry::IDetector & detector(const size_t index) const
Return a const reference to the detector or detector group of the spectrum with given index.
A property class for workspaces.
A validator which checks that the unit of the workspace referred to by a WorkspaceProperty is the exp...
Requires an estimate for the initial neutron energy which it uses to search for monitor peaks and fro...
Definition: GetEi2.h:48
double m_tof_window
The percentage deviation from the estimated peak time that defines the peak region.
Definition: GetEi2.h:107
const double m_bkgd_frac
The fraction of the peakwidth to use in calculating background points.
Definition: GetEi2.h:115
void integrate(double &integral_val, double &integral_err, const HistogramData::HistogramX &x, const HistogramData::HistogramY &s, const HistogramData::HistogramE &e, const double xmin, const double xmax) const
Integrate the point data.
Definition: GetEi2.cpp:636
double getDistanceFromSource(const size_t ws_index, const API::SpectrumInfo &spectrumInfo) const
Get the distance from the source of the detector at the workspace index given.
Definition: GetEi2.cpp:273
double m_t_to_mev
Conversion factor between time and energy.
Definition: GetEi2.h:104
API::MatrixWorkspace_sptr rebin(const API::MatrixWorkspace_sptr &monitor_ws, const double first, const double width, const double end)
Rebin the given workspace using the given parameters.
Definition: GetEi2.cpp:614
void storeEi(const double ei) const
Store the incident energy within the sample object.
Definition: GetEi2.cpp:722
double calculateFirstMoment(const API::MatrixWorkspace_sptr &monitor_ws, const double prominence)
Calculate the value of the first moment of the given spectrum.
Definition: GetEi2.cpp:590
void exec() override
Execute the algorithm.
Definition: GetEi2.cpp:117
double calculatePeakWidthAtHalfHeight(const API::MatrixWorkspace_sptr &data_ws, const double prominence, std::vector< double > &peak_x, std::vector< double > &peak_y, std::vector< double > &peak_e) const
Calculate peak width.
Definition: GetEi2.cpp:368
API::MatrixWorkspace_sptr m_input_ws
The input workspace.
Definition: GetEi2.h:98
double calculateEi(const double initial_guess)
Calculate Ei from the initial guess given.
Definition: GetEi2.cpp:148
void init() override
Initialize the algorithm.
Definition: GetEi2.cpp:49
const double m_peak_signif
Number of std deviations to consider a peak.
Definition: GetEi2.h:109
std::pair< int, double > m_peak1_pos
The calculated position of the first peak.
Definition: GetEi2.h:100
bool m_fixedei
True if the Ei should be fixed at the guess energy.
Definition: GetEi2.h:102
const double m_peak_deriv
Number of std deviations to consider a peak for the derivative.
Definition: GetEi2.h:111
double calculatePeakPosition(const size_t ws_index, const double t_min, const double t_max)
Calculate the peak position within the given window.
Definition: GetEi2.cpp:299
const double m_binwidth_frac
The fraction of the peak width for the new bins.
Definition: GetEi2.h:113
API::MatrixWorkspace_sptr extractSpectrum(const size_t ws_index, const double start, const double end)
Extract the required spectrum from the input workspace.
Definition: GetEi2.cpp:334
Exception for errors associated with the instrument definition.
Definition: Exception.h:220
IPropertyManager * setProperty(const std::string &name, const T &value)
Templated method to set the value of a PropertyWithValue.
void debug(const std::string &msg)
Logs at debug level.
Definition: Logger.cpp:114
void notice(const std::string &msg)
Logs at notice level.
Definition: Logger.cpp:95
void error(const std::string &msg)
Logs at error level.
Definition: Logger.cpp:77
void warning(const std::string &msg)
Logs at warning level.
Definition: Logger.cpp:86
bool is(int level) const
Returns true if at least the given log level is set.
Definition: Logger.cpp:146
void information(const std::string &msg)
Logs at information level.
Definition: Logger.cpp:105
The concrete, templated class for properties.
Base class for properties.
Definition: Property.h:94
std::shared_ptr< MatrixWorkspace > MatrixWorkspace_sptr
shared pointer to the matrix workspace base class
std::shared_ptr< Parameter > Parameter_sptr
Typedef for the shared pointer.
Definition: Parameter.h:195
std::shared_ptr< const IComponent > IComponent_const_sptr
Typdef of a shared pointer to a const IComponent.
Definition: IComponent.h:161
std::shared_ptr< const Instrument > Instrument_const_sptr
Shared pointer to an const instrument object.
static constexpr double NeutronMass
Mass of the neutron in kg.
static constexpr double meV
1 meV in Joules.
constexpr int EMPTY_INT() noexcept
Returns what we consider an "empty" integer within a property.
Definition: EmptyValues.h:25
int32_t specnum_t
Typedef for a spectrum Number.
Definition: IDTypes.h:16
constexpr double EMPTY_DBL() noexcept
Returns what we consider an "empty" double within a property.
Definition: EmptyValues.h:43
Generate a tableworkspace to store the calibration results.
std::string to_string(const wide_integer< Bits, Signed > &n)
static void makeDistribution(const MatrixWorkspace_sptr &workspace, const bool forwards=true)
Divides the data in a workspace by the bin width to make it a distribution.
@ InOut
Both an input & output workspace.
Definition: Property.h:55
@ Input
An input workspace.
Definition: Property.h:53
@ Output
An output workspace.
Definition: Property.h:54