29using namespace Kernel;
48 auto wsValidator = std::make_shared<WorkspaceUnitValidator>(
"Empty");
51 "Name of the input workspace");
53 "Name of the output workspace, can be the same as the input");
55 "A list of spectra workspace indices as a string with ranges; e.g. "
57 "Optional: if not specified, then the Start/EndIndex fields "
59 "If specified, the range and the list are combined (without "
60 "duplicating indices). For example, a range of 10 to 20 and "
61 "a list '12,15,26,28' gives '10-20,26,28'.");
64 "A list of spectra indices as a string with ranges; e.g. "
66 "Optional: if not specified, then the Start/EndIndex fields "
68 "If specified, the range and the list are combined (without "
69 "duplicating indices). For example, a range of 10 to 20 and "
70 "a list '12,15,26,28' gives '10-20,26,28'.");
73 auto mustBePositive = std::make_shared<BoundedValidator<int>>();
74 mustBePositive->setLower(0);
76 "Value of elastic peak position if none of the above are filled in.");
78 "Workspace Index used for elastic peak position above.");
88 std::vector<int> spectraIndices =
getProperty(
"ListOfSpectraIndices");
89 std::vector<int> channelIndices =
getProperty(
"ListOfChannelIndices");
91 const int eppSpectrum =
getProperty(
"ElasticPeakPositionSpectrum");
93 std::vector<double> tofAxis;
94 double channelWidth{0.0};
95 if (
m_inputWS->run().hasProperty(
"channel_width")) {
96 channelWidth =
m_inputWS->run().getPropertyValueAsType<
double>(
"channel_width");
98 const std::string mesg =
"No property channel_width found in the input workspace....";
99 throw std::runtime_error(mesg);
106 "ElasticPeakPositionSpectrum");
108 double wavelength{0.0};
109 if (
m_inputWS->run().hasProperty(
"wavelength")) {
110 wavelength =
m_inputWS->run().getPropertyValueAsType<
double>(
"wavelength");
112 std::string mesg =
"No property wavelength found in the input workspace....";
113 throw std::runtime_error(mesg);
133 if (
g_log.
is(Logger::Priority::PRIO_DEBUG)) {
134 for (
auto &e : eppMap) {
135 g_log.
debug() <<
"Spectra idx =" << e.first <<
", epp=" << e.second <<
'\n';
141 tofAxis =
makeTofAxis(eppAndEpTof.first, eppAndEpTof.second,
m_inputWS->blocksize() + 1, channelWidth);
161 auto nHist =
m_inputWS->getNumberHistograms();
163 g_log.
information(
"No spectrum number given. Using all spectra to calculate "
164 "the elastic peak.");
168 for (
unsigned int i = 0; i < nHist; ++i)
171 const auto found = std::find_if(
172 v.cbegin(), v.cend(), [nHist](
const auto index) { return index < 0 || static_cast<size_t>(index) >= nHist; });
173 if (found != v.cend()) {
174 throw std::runtime_error(
"Spectra index out of limits: " +
std::to_string(*found));
185 auto blockSize =
m_inputWS->blocksize() + 1;
188 "spectrum!) to calculate the elastic peak.");
190 v.reserve(blockSize);
191 for (
unsigned int i = 0; i < blockSize; ++i)
194 const auto found = std::find_if(v.cbegin(), v.cend(), [blockSize](
const auto index) {
195 return index < 0 || static_cast<size_t>(index) >= blockSize;
197 if (found != v.cend()) {
198 throw std::runtime_error(
"Channel index out of limits: " +
std::to_string(*found));
209 const std::vector<int> &channelIndices) {
211 std::map<int, int> eppMap;
214 assert(
static_cast<size_t>(*(channelIndices.end() - 1)) <
m_inputWS->blocksize() + 1);
218 for (
auto spectrumIndex : spectraIndices) {
220 auto &thisSpecY =
m_inputWS->y(spectrumIndex);
222 int minChannelIndex = *(channelIndices.begin());
223 int maxChannelIndex = *(channelIndices.end() - 1);
226 minX =
static_cast<double>(minChannelIndex);
227 maxX =
static_cast<double>(maxChannelIndex);
230 g_log.
debug() <<
"Peak estimate :: center=" << center <<
"\t sigma=" <<
sigma <<
"\t height=" <<
height
231 <<
"\t minX=" << minX <<
"\t maxX=" << maxX <<
'\n';
235 g_log.
error() <<
"doFitGaussianPeak failed...\n";
236 throw std::runtime_error(
"Gaussin Peak Fit failed....");
239 g_log.
debug() <<
"Peak Fitting :: center=" << center <<
"\t sigma=" <<
sigma <<
"\t height=" <<
height
240 <<
"\t minX=" << minX <<
"\t maxX=" << maxX <<
'\n';
243 eppMap[spectrumIndex] =
roundUp(center);
253 double &
height,
double &minX,
double &maxX) {
255 auto maxValueIt = std::max_element(spec.begin() +
static_cast<size_t>(minX),
256 spec.begin() +
static_cast<size_t>(maxX));
258 size_t maxIndex = std::distance(spec.begin(), maxValueIt);
261 std::find_if(MantidVec::const_reverse_iterator(maxValueIt), MantidVec::const_reverse_iterator(spec.cbegin()),
262 [
maxValue](
double value) { return value < 0.5 * maxValue; });
263 auto maxFwhmIt = std::find_if(maxValueIt, spec.end(), [
maxValue](
double value) { return value < 0.5 * maxValue; });
266 double fwhm =
static_cast<double>(std::distance(minFwhmIt.base(), maxFwhmIt) + 1);
269 center =
static_cast<double>(maxIndex);
273 g_log.
debug() <<
"Peak estimate : center=" << center <<
"\t sigma=" <<
sigma <<
"\t h=" <<
height <<
'\n';
276 size_t ipeak_min = std::max(
277 std::size_t{0}, maxIndex -
static_cast<size_t>(2.5 *
static_cast<double>(std::distance(maxValueIt, maxFwhmIt))));
278 size_t ipeak_max = std::min(
279 spec.size(), maxIndex +
static_cast<size_t>(2.5 *
static_cast<double>(std::distance(maxFwhmIt, maxValueIt))));
280 size_t i_delta_peak = ipeak_max - ipeak_min;
282 g_log.
debug() <<
"Peak estimate xmin/max: " << ipeak_min - 1 <<
"\t" << ipeak_max + 1 <<
'\n';
284 minX =
static_cast<double>(ipeak_min - 2 * i_delta_peak);
285 maxX =
static_cast<double>(ipeak_max + 2 * i_delta_peak);
302 double startX,
double endX) {
311 auto gaussianpeak = std::dynamic_pointer_cast<API::IPeakFunction>(temppeak);
312 gaussianpeak->setHeight(
height);
313 gaussianpeak->setCentre(center);
314 gaussianpeak->setFwhm(
sigma);
317 double centerleftend = center -
sigma * 0.5;
318 double centerrightend = center +
sigma * 0.5;
319 std::ostringstream os;
320 os << centerleftend <<
" < PeakCentre < " << centerrightend;
321 auto centerbound = std::unique_ptr<API::IConstraint>(
323 gaussianpeak->addConstraint(std::move(centerbound));
325 g_log.
debug(
"Calling createChildAlgorithm : Fit...");
328 fitalg->initialize();
330 fitalg->setProperty(
"Function", std::dynamic_pointer_cast<API::IFunction>(gaussianpeak));
331 fitalg->setProperty(
"InputWorkspace",
m_inputWS);
332 fitalg->setProperty(
"WorkspaceIndex", workspaceindex);
333 fitalg->setProperty(
"Minimizer",
"Levenberg-MarquardtMD");
334 fitalg->setProperty(
"CostFunction",
"Least squares");
335 fitalg->setProperty(
"MaxIterations", 1000);
336 fitalg->setProperty(
"Output",
"FitGaussianPeak");
337 fitalg->setProperty(
"StartX", startX);
338 fitalg->setProperty(
"EndX", endX);
341 bool successfulfit = fitalg->execute();
342 if (!fitalg->isExecuted() || !successfulfit) {
344 g_log.
warning() <<
"Fitting Gaussian peak for peak around " << gaussianpeak->centre() <<
'\n';
349 center = gaussianpeak->centre();
350 height = gaussianpeak->height();
351 double fwhm = gaussianpeak->fwhm();
365 double wavelength{0.0};
366 if (
m_inputWS->run().hasProperty(
"wavelength")) {
367 wavelength =
m_inputWS->run().getPropertyValueAsType<
double>(
"wavelength");
369 std::string mesg =
"No property wavelength found in the input workspace....";
370 throw std::runtime_error(mesg);
373 std::vector<double> epTofList;
374 std::vector<int> eppList;
377 for (
const auto &epp : eppMap) {
381 g_log.
error() <<
"firstL2=" << firstL2 <<
" , "
382 <<
"l2=" <<
l2 <<
'\n';
383 throw std::runtime_error(
"All the pixels for selected spectra must have "
384 "the same distance from the sample!");
390 eppList.emplace_back(epp.first);
392 g_log.
debug() <<
"WS index = " << epp.first <<
", l1 = " << l1 <<
", l2 = " <<
l2
393 <<
", TOF(l1+l2) = " << *(epTofList.end() - 1) <<
'\n';
396 double averageEpTof =
397 std::accumulate(epTofList.begin(), epTofList.end(), 0.0) /
static_cast<double>(epTofList.size());
398 int averageEpp =
roundUp(
static_cast<double>(std::accumulate(eppList.begin(), eppList.end(), 0)) /
399 static_cast<double>(eppList.size()));
401 g_log.
debug() <<
"Average epp=" << averageEpp <<
" , Average epTof=" << averageEpTof <<
'\n';
402 return std::make_pair(averageEpp, averageEpTof);
406 if (wavelength <= 0) {
407 throw std::runtime_error(
"Wavelenght is <= 0");
412 return distance / velocity;
426 std::vector<double> axis(size);
428 g_log.
debug() <<
"Building the TOF X Axis: epp=" << epp <<
", epTof=" << epTof <<
", Channel Width=" << channelWidth
430 for (
size_t i = 0; i < size; ++i) {
431 axis[i] = epTof + channelWidth *
static_cast<double>(
static_cast<int>(i) - epp) -
434 g_log.
debug() <<
"TOF X Axis: [start,end] = [" << *axis.begin() <<
"," << *(axis.end() - 1) <<
"]\n";
440 const size_t numberOfSpectra =
m_inputWS->getNumberHistograms();
442 g_log.
debug() <<
"Setting the TOF X Axis for numberOfSpectra=" << numberOfSpectra <<
'\n';
444 HistogramData::BinEdges edges(tofAxis);
445 Progress prog(
this, 0.0, 0.2, numberOfSpectra);
447 for (
size_t i = 0; i < numberOfSpectra; ++i) {
449 outputWS->setBinEdges(i, edges);
#define DECLARE_ALGORITHM(classname)
double value
The value of the point.
std::map< DeltaEMode::Type, std::string > index
void declareProperty(std::unique_ptr< Kernel::Property > p, const std::string &doc="") override
Add a property to the list of managed properties.
TypedValue getProperty(const std::string &name) const override
Get the value of a property.
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.
Helper class for reporting progress from algorithms.
double l2(const size_t index) const
Returns L2 (distance from sample to spectrum).
double l1() const
Returns L1 (distance from source to sample).
A property class for workspaces.
bool doFitGaussianPeak(int, double &, double &, double &, double, double)
Fit peak without background i.e, with background removed inspired from FitPowderDiffPeaks....
double calculateTOF(double, double)
API::MatrixWorkspace_sptr m_outputWS
void validateChannelIndices(std::vector< int > &v)
Check if the channel indices are in the limits of the number of the block size in the input workspace...
void validateWorkspaceIndices(std::vector< int > &v)
Check if spectra indices are in the limits of the number of histograms in the input workspace.
void exec() override
Execute the algorithm.
const std::string category() const override
Algorithm's category for identification.
void estimateFWHM(const Mantid::HistogramData::HistogramY &, double &, double &, double &, double &, double &)
Estimated the FWHM for Gaussian peak fitting.
std::map< int, int > findElasticPeakPositions(const std::vector< int > &, const std::vector< int > &)
Looks for the EPP positions in the spectraIndices.
bool areEqual(double, double, double)
Compare two double with a precision epsilon.
int version() const override
Algorithm's version for identification.
std::vector< double > makeTofAxis(int, double, size_t, double)
Builds the X time axis.
void init() override
Initialize the algorithm's properties.
std::pair< int, double > findAverageEppAndEpTof(const std::map< int, int > &)
Finds the TOF for a given epp.
DataObjects::Workspace2D_sptr m_inputWS
void setTofInWS(const std::vector< double > &, const API::MatrixWorkspace_sptr &)
const API::SpectrumInfo * m_spectrumInfo
Support for a property that holds an array of values.
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.
void error(const std::string &msg)
Logs at error level.
void warning(const std::string &msg)
Logs at warning level.
bool is(int level) const
Returns true if at least the given log level is set.
void information(const std::string &msg)
Logs at information level.
void report()
Increments the loop counter by 1, then sends the progress notification on behalf of its algorithm.
static T & Instance()
Return a reference to the Singleton instance, creating it if it does not already exist Creation is do...
std::shared_ptr< MatrixWorkspace > MatrixWorkspace_sptr
shared pointer to the matrix workspace base class
static constexpr double NeutronMass
Mass of the neutron in kg.
static constexpr double h
Planck constant in J*s.
constexpr int EMPTY_INT() noexcept
Returns what we consider an "empty" integer within a property.
std::string to_string(const wide_integer< Bits, Signed > &n)
@ Input
An input workspace.
@ Output
An output workspace.