Mantid
Loading...
Searching...
No Matches
GeneratePeaks.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#include "MantidAPI/Column.h"
16#include "MantidHistogramData/HistogramBuilder.h"
17#include "MantidIndexing/IndexInfo.h"
21#include "MantidTypes/SpectrumDefinition.h"
22
23#include <boost/algorithm/string/classification.hpp>
24#include <boost/algorithm/string/split.hpp>
25
26using namespace Mantid::Kernel;
27using namespace Mantid::API;
28using namespace Mantid::DataObjects;
29using namespace Mantid::HistogramData;
30
31namespace Mantid::Algorithms {
32DECLARE_ALGORITHM(GeneratePeaks)
33
34//----------------------------------------------------------------------------------------------
38 : m_peakFunction(), m_bkgdFunction(), m_vecPeakParamValues(), m_vecBkgdParamValues(), m_SpectrumMap(),
39 m_spectraSet(), m_useAutoBkgd(false), m_funcParamWS(), inputWS(), m_newWSFromParent(false), binParameters(),
40 m_genBackground(false), m_useRawParameter(false), m_maxChi2(0.), m_numPeakWidth(0.), m_funcParameterNames(),
41 i_height(-1), i_centre(-1), i_width(-1), i_a0(-1), i_a1(-1), i_a2(-1), m_useFuncParamWS(false), m_wsIndex(-1) {}
42
43//----------------------------------------------------------------------------------------------
48 "PeakParametersWorkspace", "", Direction::Input, PropertyMode::Optional),
49 "Input TableWorkspace for peak's parameters.");
50
51 // Peak function properties
52 std::vector<std::string> peakNames = FunctionFactory::Instance().getFunctionNames<IPeakFunction>();
53 std::vector<std::string> peakFullNames = addFunctionParameterNames(peakNames);
54 declareProperty("PeakType", "", std::make_shared<StringListValidator>(peakFullNames), "Peak function type. ");
55
56 for (size_t i = 0; i < peakFullNames.size(); ++i)
57 g_log.debug() << "Peak function " << i << ": " << peakFullNames[i] << "\n";
58
59 declareProperty(std::make_unique<ArrayProperty<double>>("PeakParameterValues"),
60 "List of peak parameter values. They must have a 1-to-1 "
61 "mapping to PeakParameterNames list. ");
62
63 // Background properties
64 std::vector<std::string> bkgdtypes{"Auto",
65 "Flat (A0)",
66 "Linear (A0, A1)",
67 "Quadratic (A0, A1, A2)",
68 "Flat"
69 "Linear",
70 "Quadratic"};
71 declareProperty("BackgroundType", "Linear", std::make_shared<StringListValidator>(bkgdtypes), "Type of Background.");
72
73 declareProperty(std::make_unique<ArrayProperty<double>>("BackgroundParameterValues"),
74 "List of background parameter values. They must have a "
75 "1-to-1 mapping to PeakParameterNames list. ");
76
79 "InputWorkspace (optional) to take information for the instrument, and "
80 "where to evaluate the x-axis.");
81
82 declareProperty("WorkspaceIndex", 0,
83 "Spectrum of the peak to be generated. "
84 "It is only applied to the case by input parameter values in "
85 "vector format. ");
86
87 declareProperty(std::make_unique<Kernel::ArrayProperty<double>>("BinningParameters",
88 std::make_shared<Kernel::RebinParamsValidator>(true)),
89 "A comma separated list of first bin boundary, width, last "
90 "bin boundary. Optionally\n"
91 "this can be followed by a comma and more widths and last "
92 "boundary pairs.\n"
93 "Negative width values indicate logarithmic binning.");
94
95 declareProperty("NumberWidths", 2., "Number of peak width to evaluate each peak for. Default=2.");
97 std::make_unique<API::WorkspaceProperty<API::MatrixWorkspace>>("OutputWorkspace", "", Direction::Output),
98 "Output Workspace to put the calculated data.");
99
100 declareProperty("GenerateBackground", true, "Whether or not to generate the background");
101
102 declareProperty("MaxAllowedChi2", 100.0, "Maximum chi^2 of the peak allowed to calculate. Default 100.");
103
104 declareProperty("IgnoreWidePeaks", false,
105 "If selected, the peaks that are wider than fit window "
106 "(denoted by negative chi^2) are ignored.");
107
108 declareProperty("IsRawParameter", true,
109 "Flag to show whether the parameter table contains raw parameters. "
110 "In the case that parameter values are input via vector, and this flag "
111 "is set to false, "
112 "the default order of effective peak parameters is centre, height and "
113 "width; "
114 "the default order of effective background parameters is A0, A1 and "
115 "A2. ");
116}
117
118//----------------------------------------------------------------------------------------------
122 // Process input parameters
123 std::string peaktype, bkgdtype;
124 processAlgProperties(peaktype, bkgdtype);
125
126 // Create functions
127 createFunction(peaktype, bkgdtype);
128
129 // Process parameter table
130 if (m_useFuncParamWS) {
133 }
134
135 // Create output workspace
137 setProperty("OutputWorkspace", outputWS);
138
139 // Generate peaks
140 std::map<specnum_t, std::vector<std::pair<double, API::IFunction_sptr>>> functionmap;
142 importPeaksFromTable(functionmap);
143 else {
144 std::vector<std::pair<double, API::IFunction_sptr>> vecpeakfunc;
145 importPeakFromVector(vecpeakfunc);
146 functionmap.emplace(m_wsIndex, vecpeakfunc);
147 }
148
149 generatePeaks(functionmap, outputWS);
150}
151
152//----------------------------------------------------------------------------------------------
155void GeneratePeaks::processAlgProperties(std::string &peakfunctype, std::string &bkgdfunctype) {
156 // Function parameters
157 std::string paramwsname = getPropertyValue("PeakParametersWorkspace");
158 if (!paramwsname.empty()) {
159 // Using parameter table workspace has a higher priority
160 m_useFuncParamWS = true;
161 m_funcParamWS = getProperty("PeakParametersWorkspace");
162 } else {
163 m_useFuncParamWS = false;
164 m_vecPeakParamValues = getProperty("PeakParameterValues");
165 m_vecBkgdParamValues = getProperty("BackgroundParameterValues");
166 }
167
168 peakfunctype = getPropertyValue("PeakType");
169 // Remove extra helping message
170 if (peakfunctype.find('(') != std::string::npos) {
171 std::vector<std::string> strs;
172 boost::split(strs, peakfunctype, boost::is_any_of(" ("));
173 peakfunctype = strs[0];
174 }
175
176 bkgdfunctype = getPropertyValue("BackgroundType");
177 // Remove extra helping message
178 if (bkgdfunctype.find('(') != std::string::npos) {
179 std::vector<std::string> strs;
180 boost::split(strs, bkgdfunctype, boost::is_any_of(" ("));
181 bkgdfunctype = strs[0];
182 }
183
184 if (bkgdfunctype == "Auto") {
185 m_useAutoBkgd = true;
186 bkgdfunctype = "Quadratic";
187 } else if (bkgdfunctype == "None") {
188 m_useAutoBkgd = false;
189 m_genBackground = false;
190 } else if (bkgdfunctype == "Linear" || bkgdfunctype == "Flat") {
191 m_useAutoBkgd = false;
192 bkgdfunctype = bkgdfunctype + "Background";
193 }
194
195 binParameters = this->getProperty("BinningParameters");
196 inputWS = this->getProperty("InputWorkspace");
197
198 m_genBackground = getProperty("GenerateBackground");
199
200 m_useRawParameter = getProperty("IsRawParameter");
201
202 // Special properties related
203 m_maxChi2 = this->getProperty("MaxAllowedChi2");
204 m_numPeakWidth = this->getProperty("NumberWidths");
205
206 // Spectrum set if not using parameter table workspace
207 // One and only one peak
208 if (!m_useFuncParamWS) {
209 m_wsIndex = getProperty("WorkspaceIndex");
210 m_spectraSet.insert(static_cast<specnum_t>(m_wsIndex));
211 m_SpectrumMap.emplace(static_cast<specnum_t>(m_wsIndex), 0);
212 }
213}
214
215//----------------------------------------------------------------------------------------------
221 std::map<specnum_t, std::vector<std::pair<double, API::IFunction_sptr>>> &functionmap) {
222 size_t numpeaks = m_funcParamWS->rowCount();
223 size_t icolchi2 = m_funcParamWS->columnCount() - 1;
224 size_t numpeakparams = m_peakFunction->nParams();
225 size_t numbkgdparams = 0;
226 if (m_bkgdFunction)
227 numbkgdparams = m_bkgdFunction->nParams();
228 else
229 g_log.warning("There is no background function specified. ");
230
231 // Create data structure for all peaks functions
232 std::map<specnum_t, std::vector<std::pair<double, API::IFunction_sptr>>>::iterator mapiter;
233
234 // Go through the table workspace to create peak/background functions
235 for (size_t ipeak = 0; ipeak < numpeaks; ++ipeak) {
236 // Spectrum
237 int wsindex = m_funcParamWS->cell<int>(ipeak, 0);
238
239 // Ignore peak with large chi^2/Rwp
240 double chi2 = m_funcParamWS->cell<double>(ipeak, icolchi2);
241 if (chi2 > m_maxChi2) {
242 g_log.notice() << "Skip Peak " << ipeak << " (chi^2 " << chi2 << " > " << m_maxChi2 << ".) "
243 << "\n";
244 continue;
245 } else if (chi2 < 0.) {
246 g_log.notice() << "Skip Peak " << ipeak << " (chi^2 " << chi2 << " < 0 )"
247 << "\n";
248 continue;
249 } else {
250 g_log.debug() << "[DB] Chi-square = " << chi2 << "\n";
251 }
252
253 // Set up function
254 if (m_useRawParameter) {
255 for (size_t p = 0; p < numpeakparams; ++p) {
256 std::string parname = m_funcParameterNames[p];
257 double parvalue = m_funcParamWS->cell<double>(ipeak, p + 1);
258 m_peakFunction->setParameter(parname, parvalue);
259 }
260 if (m_genBackground) {
261 for (size_t p = 0; p < numbkgdparams; ++p) {
262 std::string parname = m_funcParameterNames[p + numpeakparams];
263 double parvalue = m_funcParamWS->cell<double>(ipeak, p + 1 + numpeakparams);
264 m_bkgdFunction->setParameter(parname, parvalue);
265 }
266 }
267 } else {
268 double tmpheight = m_funcParamWS->cell<double>(ipeak, i_height);
269 double tmpwidth = m_funcParamWS->cell<double>(ipeak, i_width);
270 double tmpcentre = m_funcParamWS->cell<double>(ipeak, i_centre);
271 m_peakFunction->setHeight(tmpheight);
272 m_peakFunction->setCentre(tmpcentre);
273 m_peakFunction->setFwhm(tmpwidth);
274
275 if (m_genBackground) {
276 double tmpa0 = m_funcParamWS->cell<double>(ipeak, i_a0);
277 double tmpa1 = m_funcParamWS->cell<double>(ipeak, i_a1);
278 double tmpa2 = m_funcParamWS->cell<double>(ipeak, i_a2);
279 m_bkgdFunction->setParameter("A0", tmpa0);
280 m_bkgdFunction->setParameter("A1", tmpa1);
281 m_bkgdFunction->setParameter("A2", tmpa2);
282 }
283 }
284
285 double centre = m_peakFunction->centre();
286
287 // Generate function to plot
288 API::CompositeFunction_sptr plotfunc = std::make_shared<CompositeFunction>();
289 plotfunc->addFunction(m_peakFunction);
290 if (m_genBackground)
291 plotfunc->addFunction(m_bkgdFunction);
292
293 // Clone to another function
294 API::IFunction_sptr clonefunction = plotfunc->clone();
295
296 // Existed?
297 mapiter = functionmap.find(wsindex);
298 if (mapiter == functionmap.end()) {
299 std::vector<std::pair<double, API::IFunction_sptr>> tempvector;
300 std::pair<std::map<specnum_t, std::vector<std::pair<double, API::IFunction_sptr>>>::iterator, bool> ret;
301 ret = functionmap.emplace(wsindex, tempvector);
302 mapiter = ret.first;
303 }
304
305 // Generate peak function
306 mapiter->second.emplace_back(centre, clonefunction);
307
308 g_log.information() << "Peak " << ipeak << ": Spec = " << wsindex << " func: " << clonefunction->asString() << "\n";
309
310 } // ENDFOR (ipeak)
311
312 // Sort by peak position
313 for (mapiter = functionmap.begin(); mapiter != functionmap.end(); ++mapiter) {
314 std::vector<std::pair<double, API::IFunction_sptr>> &vec_centrefunc = mapiter->second;
315 std::sort(vec_centrefunc.begin(), vec_centrefunc.end());
316 }
317}
318
319//----------------------------------------------------------------------------------------------
322void GeneratePeaks::importPeakFromVector(std::vector<std::pair<double, API::IFunction_sptr>> &functionmap) {
323 API::CompositeFunction_sptr compfunc = std::make_shared<API::CompositeFunction>();
324
325 // Set up and clone peak function
326 if (m_useRawParameter) {
327 // Input vector of values are for raw parameter name
328 size_t numpeakparams = m_peakFunction->nParams();
329 if (m_vecPeakParamValues.size() == numpeakparams) {
330 for (size_t i = 0; i < numpeakparams; ++i)
331 m_peakFunction->setParameter(i, m_vecPeakParamValues[i]);
332 } else {
333 // Number of input parameter values is not correct. Throw!
334 std::stringstream errss;
335 errss << "Number of input peak parameters' value (" << m_vecPeakParamValues.size()
336 << ") is not correct (should be " << numpeakparams << " for peak of type " << m_peakFunction->name()
337 << "). ";
338 throw std::runtime_error(errss.str());
339 }
340 } else {
341 // Input vector of values are for effective parameter names
342 if (m_vecPeakParamValues.size() != 3)
343 throw std::runtime_error("Input peak parameters must have 3 numbers for "
344 "effective parameter names.");
345
349 }
350 compfunc->addFunction(m_peakFunction->clone());
351
352 // Set up and clone background function
353 if (m_genBackground) {
354 size_t numbkgdparams = m_bkgdFunction->nParams();
355 if (m_useRawParameter) {
356 // Raw background parameters
357 if (m_vecBkgdParamValues.size() != numbkgdparams)
358 throw std::runtime_error("Number of background parameters' value is not correct. ");
359 else {
360 for (size_t i = 0; i < numbkgdparams; ++i)
361 m_bkgdFunction->setParameter(i, m_vecBkgdParamValues[i]);
362 }
363 } else {
364 // Effective background parameters
365 if (m_vecBkgdParamValues.size() < 3 && m_vecBkgdParamValues.size() < numbkgdparams) {
366 throw std::runtime_error("There is no enough effective background parameter values.");
367 }
368
369 // FIXME - Assume that all background functions define parameter i for A_i
370 for (size_t i = 0; i < numbkgdparams; ++i)
371 m_bkgdFunction->setParameter(i, m_vecBkgdParamValues[i]);
372 }
373
374 compfunc->addFunction(m_bkgdFunction->clone());
375 }
376
377 // Set up function map
378 functionmap.emplace_back(m_peakFunction->centre(), compfunc);
379}
380
381//----------------------------------------------------------------------------------------------
388 const std::map<specnum_t, std::vector<std::pair<double, API::IFunction_sptr>>> &functionmap,
389 const API::MatrixWorkspace_sptr &dataWS) {
390 // Calcualte function
391 std::map<specnum_t, std::vector<std::pair<double, API::IFunction_sptr>>>::const_iterator mapiter;
392 for (mapiter = functionmap.begin(); mapiter != functionmap.end(); ++mapiter) {
393 // Get spec id and translated to wsindex in the output workspace
394 specnum_t specid = mapiter->first;
395 specnum_t wsindex;
397 wsindex = specid;
398 else
399 wsindex = m_SpectrumMap[specid];
400
401 const std::vector<std::pair<double, API::IFunction_sptr>> &vec_centrefunc = mapiter->second;
402 size_t numpeaksinspec = mapiter->second.size();
403
404 for (size_t ipeak = 0; ipeak < numpeaksinspec; ++ipeak) {
405 const std::pair<double, API::IFunction_sptr> &centrefunc = vec_centrefunc[ipeak];
406
407 // Determine boundary
408 API::IPeakFunction_sptr thispeak = getPeakFunction(centrefunc.second);
409 double centre = centrefunc.first;
410 double fwhm = thispeak->fwhm();
411
412 //
413 const auto &X = dataWS->x(wsindex);
414 double leftbound = centre - m_numPeakWidth * fwhm;
415 if (ipeak > 0) {
416 // Not left most peak.
417 API::IPeakFunction_sptr leftPeak = getPeakFunction(vec_centrefunc[ipeak - 1].second);
418
419 double middle = 0.5 * (centre + leftPeak->centre());
420 if (leftbound < middle)
421 leftbound = middle;
422 }
423 auto left = std::lower_bound(X.cbegin(), X.cend(), leftbound);
424 if (left == X.end())
425 left = X.begin();
426
427 double rightbound = centre + m_numPeakWidth * fwhm;
428 if (ipeak != numpeaksinspec - 1) {
429 // Not the rightmost peak
430 IPeakFunction_sptr rightPeak = getPeakFunction(vec_centrefunc[ipeak + 1].second);
431 double middle = 0.5 * (centre + rightPeak->centre());
432 if (rightbound > middle)
433 rightbound = middle;
434 }
435 auto right = std::lower_bound(left + 1, X.cend(), rightbound);
436
437 // Build domain & function
439 right); // dataWS->dataX(wsindex));
440
441 // Evaluate the function
442 API::FunctionValues values(domain);
443 centrefunc.second->function(domain, values);
444
445 // Put to output
446 std::size_t offset = (left - X.begin());
447 std::size_t numY = values.size();
448
449 auto &dataY = dataWS->mutableY(wsindex);
450 for (std::size_t i = 0; i < numY; i++) {
451 dataY[i + offset] += values[i];
452 }
453
454 } // ENDFOR(ipeak)
455 }
456}
457
458//----------------------------------------------------------------------------------------------
462void GeneratePeaks::createFunction(std::string &peaktype, std::string &bkgdtype) {
463 // Create peak function
464 m_peakFunction = std::dynamic_pointer_cast<IPeakFunction>(API::FunctionFactory::Instance().createFunction(peaktype));
465
466 // create the background
467 if (m_genBackground)
469 std::dynamic_pointer_cast<IBackgroundFunction>(API::FunctionFactory::Instance().createFunction(bkgdtype));
470}
471
472//----------------------------------------------------------------------------------------------
476 using namespace boost::algorithm;
477
478 // Initial check
479 std::vector<std::string> colnames = m_funcParamWS->getColumnNames();
480
481 if (colnames[0] != "spectrum")
482 throw std::runtime_error("First column must be 'spectrum' in integer. ");
483 if (colnames.back() != "chi2")
484 throw std::runtime_error("Last column must be 'chi2'.");
485
486 // Process column names in case that there are not same as parameter names
487 // fx.name might be available
488 m_funcParameterNames.resize(colnames.size() - 2);
489 for (size_t i = 0; i < colnames.size() - 2; ++i) {
490 std::string str = colnames[i + 1];
491 std::vector<std::string> tokens;
492 split(tokens, str, is_any_of(".")); // here it is
493 m_funcParameterNames[i] = tokens.back();
494 }
495
496 // Check column number
497 size_t numparamcols = colnames.size() - 2;
498 if (m_useRawParameter) {
499 // Number of parameters must be peak parameters + background parameters
500 size_t numpeakparams = m_peakFunction->nParams();
501 size_t numbkgdparams = 0;
502 if (m_genBackground)
503 numbkgdparams = m_bkgdFunction->nParams();
504 size_t numparams = numpeakparams + numbkgdparams;
505 if (numparamcols < numparams)
506 throw std::runtime_error("Parameters number is not correct. ");
507 else if (numparamcols > numparams)
508 g_log.warning("Number of parameters given in table workspace is more than "
509 "number of parameters of function(s) to generate peaks. ");
510
511 // Check column names are same as function parameter naems
512 for (size_t i = 0; i < numpeakparams; ++i) {
514 std::stringstream errss;
515 errss << "Peak function " << m_peakFunction->name() << " does not have paramter " << m_funcParameterNames[i]
516 << "\n"
517 << "Allowed function parameters are ";
518 std::vector<std::string> parnames = m_peakFunction->getParameterNames();
519 for (auto &parname : parnames)
520 errss << parname << ", ";
521 throw std::runtime_error(errss.str());
522 }
523 }
524
525 // Background function
526 for (size_t i = 0; i < numbkgdparams; ++i) {
527 if (!hasParameter(m_bkgdFunction, m_funcParameterNames[i + numpeakparams])) {
528 std::stringstream errss;
529 errss << "Background function does not have paramter " << m_funcParameterNames[i + numpeakparams];
530 throw std::runtime_error(errss.str());
531 }
532 }
533 } else {
534 // Effective parameter names
535 if (numparamcols != 6)
536 throw std::runtime_error("Number of columns must be 6 if not using raw.");
537
538 // Find the column index of height, width, centre, a0, a1 and a2
539 i_height = static_cast<int>(std::find(m_funcParameterNames.begin(), m_funcParameterNames.end(), "height") -
540 m_funcParameterNames.begin()) +
541 1;
542 i_centre = static_cast<int>(std::find(m_funcParameterNames.begin(), m_funcParameterNames.end(), "centre") -
543 m_funcParameterNames.begin()) +
544 1;
545 i_width = static_cast<int>(std::find(m_funcParameterNames.begin(), m_funcParameterNames.end(), "width") -
546 m_funcParameterNames.begin()) +
547 1;
548 i_a0 = static_cast<int>(std::find(m_funcParameterNames.begin(), m_funcParameterNames.end(), "backgroundintercept") -
549 m_funcParameterNames.begin()) +
550 1;
551 i_a1 = static_cast<int>(std::find(m_funcParameterNames.begin(), m_funcParameterNames.end(), "backgroundslope") -
552 m_funcParameterNames.begin()) +
553 1;
554 i_a2 = static_cast<int>(std::find(m_funcParameterNames.begin(), m_funcParameterNames.end(), "A2") -
555 m_funcParameterNames.begin()) +
556 1;
557 }
558}
559
560//----------------------------------------------------------------------------------------------
566 size_t numpeaks = peakParmsWS->rowCount();
567 API::Column_const_sptr col = peakParmsWS->getColumn("spectrum");
568
569 for (size_t ipk = 0; ipk < numpeaks; ipk++) {
570 // Spectrum
571 auto specid = static_cast<specnum_t>((*col)[ipk]);
572 m_spectraSet.insert(specid);
573
574 std::stringstream outss;
575 outss << "Peak " << ipk << ": specid = " << specid;
576 g_log.debug(outss.str());
577 }
578
579 specnum_t icount = 0;
580 for (const auto specnum : m_spectraSet) {
581 m_SpectrumMap.emplace(specnum, icount);
582 ++icount;
583 }
584}
585
586//----------------------------------------------------------------------------------------------
590 // Not a composite function
591 API::CompositeFunction_sptr compfunc = std::dynamic_pointer_cast<API::CompositeFunction>(infunction);
592
593 // If it is a composite function (complete part for return)
594 if (compfunc) {
595 for (size_t i = 0; i < compfunc->nFunctions(); ++i) {
596 API::IFunction_sptr subfunc = compfunc->getFunction(i);
597 API::IPeakFunction_sptr peakfunc = std::dynamic_pointer_cast<API::IPeakFunction>(subfunc);
598 if (peakfunc)
599 return peakfunc;
600 }
601 }
602
603 // Return if not a composite function
604 API::IPeakFunction_sptr peakfunc = std::dynamic_pointer_cast<API::IPeakFunction>(infunction);
605
606 return peakfunc;
607}
608
609//----------------------------------------------------------------------------------------------
612bool GeneratePeaks::hasParameter(const API::IFunction_sptr &function, const std::string &paramname) {
613 std::vector<std::string> parnames = function->getParameterNames();
614 std::vector<std::string>::iterator piter;
615 piter = std::find(parnames.begin(), parnames.end(), paramname);
616 return piter != parnames.end();
617}
618
619//----------------------------------------------------------------------------------------------
623 // Reference workspace and output workspace
625
626 m_newWSFromParent = true;
627 if (!inputWS && binParameters.empty()) {
628 // Error! Neither bin parameters or reference workspace is given.
629 std::string errmsg("Must define either InputWorkspace or BinningParameters.");
630 g_log.error(errmsg);
631 throw std::invalid_argument(errmsg);
632 } else if (inputWS) {
633 // Generate Workspace2D from input workspace
634 if (!binParameters.empty())
635 g_log.notice() << "Both binning parameters and input workspace are given. "
636 << "Using input worksapce to generate output workspace!\n";
637
638 HistogramBuilder builder;
639 builder.setX(inputWS->x(0).size());
640 builder.setY(inputWS->y(0).size());
641
642 builder.setDistribution(inputWS->isDistribution());
643 outputWS = create<MatrixWorkspace>(*inputWS, builder.build());
644 // Only copy the X-values from spectra with peaks specified in the table
645 // workspace.
646 for (const auto &iws : m_spectraSet) {
647 outputWS->setSharedX(iws, inputWS->sharedX(iws));
648 }
649
650 m_newWSFromParent = true;
651 } else {
652 // Generate a one-spectrum Workspace2D from binning
654 m_newWSFromParent = false;
655 }
656
657 return outputWS;
658}
659
660//----------------------------------------------------------------------------------------------
664 // Check validity
665 if (m_spectraSet.empty())
666 throw std::invalid_argument("Input spectra list is empty. Unable to generate a new workspace.");
667
668 if (binparameters.size() < 3) {
669 std::stringstream errss;
670 errss << "Number of input binning parameters are not enough (" << binparameters.size() << "). "
671 << "Binning parameters should be 3 (x0, step, xf).";
672 g_log.error(errss.str());
673 throw std::invalid_argument(errss.str());
674 }
675
676 double x0 = binparameters[0];
677 double dx = binparameters[1]; // binDelta
678 double xf = binparameters[2]; // max value
679 if (x0 >= xf || (xf - x0) < dx || dx == 0.) {
680 std::stringstream errss;
681 errss << "Order of input binning parameters is not correct. It is not "
682 "logical to have "
683 << "x0 = " << x0 << ", xf = " << xf << ", dx = " << dx;
684 g_log.error(errss.str());
685 throw std::invalid_argument(errss.str());
686 }
687
688 std::vector<double> xarray;
689 double xvalue = x0;
690
691 while (xvalue <= xf) {
692 // Push current value to vector
693 xarray.emplace_back(xvalue);
694
695 // Calculate next value, linear or logarithmic
696 if (dx > 0)
697 xvalue += dx;
698 else
699 xvalue += fabs(dx) * xvalue;
700 }
701
702 std::vector<Indexing::SpectrumNumber> specNums;
703 for (const auto &item : m_SpectrumMap) {
704 specnum_t specid = item.first;
705 g_log.debug() << "Build WorkspaceIndex-Spectrum " << specNums.size() << " , " << specid << "\n";
706 specNums.emplace_back(specid);
707 }
708
709 Indexing::IndexInfo indices(specNums.size());
710 // There is no instrument, so the automatic build of a 1:1 mapping would fail.
711 // Need to set empty grouping manually.
712 indices.setSpectrumDefinitions(std::vector<SpectrumDefinition>(specNums.size()));
713 indices.setSpectrumNumbers(std::move(specNums));
714 return create<Workspace2D>(indices, BinEdges(std::move(xarray)));
715}
716
719std::vector<std::string> GeneratePeaks::addFunctionParameterNames(const std::vector<std::string> &funcnames) {
720 std::vector<std::string> vec_funcparnames;
721
722 for (auto &funcname : funcnames) {
723 // Add original name in
724 vec_funcparnames.emplace_back(funcname);
725
726 // Add a full function name and parameter names in
727 IFunction_sptr tempfunc = FunctionFactory::Instance().createFunction(funcname);
728
729 std::stringstream parnamess;
730 parnamess << funcname << " (";
731 std::vector<std::string> funcpars = tempfunc->getParameterNames();
732 for (size_t j = 0; j < funcpars.size(); ++j) {
733 parnamess << funcpars[j];
734 if (j != funcpars.size() - 1)
735 parnamess << ", ";
736 }
737 parnamess << ")";
738
739 vec_funcparnames.emplace_back(parnamess.str());
740 }
741
742 return vec_funcparnames;
743}
744
745} // namespace Mantid::Algorithms
#define DECLARE_ALGORITHM(classname)
Definition: Algorithm.h:576
double left
Definition: LineProfile.cpp:80
double right
Definition: LineProfile.cpp:81
#define fabs(x)
Definition: Matrix.cpp:22
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
std::string getPropertyValue(const std::string &name) const override
Get the value of a property as a string.
Definition: Algorithm.cpp:2026
TypedValue getProperty(const std::string &name) const override
Get the value of a property.
Definition: Algorithm.cpp:2076
Kernel::Logger & g_log
Definition: Algorithm.h:451
Implements FunctionDomain1D with its own storage in form of a std::vector.
A class to store values calculated by a function.
size_t size() const
Return the number of values.
An interface to a peak function, which extend the interface of IFunctionWithLocation by adding method...
Definition: IPeakFunction.h:51
A property class for workspaces.
GeneratePeaks : Generate peaks from a table workspace containing peak parameters.
Definition: GeneratePeaks.h:26
API::IBackgroundFunction_sptr m_bkgdFunction
Background function.
Definition: GeneratePeaks.h:87
std::set< specnum_t > m_spectraSet
Set of spectra (workspace indexes) in the original workspace that contain peaks to generate.
Definition: GeneratePeaks.h:99
int m_wsIndex
Spectrum if only 1 peak is to plot.
int i_height
Indexes of height, centre, width, a0, a1, and a2 in input parameter table.
void createFunction(std::string &peaktype, std::string &bkgdtype)
Create a function for fitting.
void init() override
Define algorithm's properties.
void exec() override
Implement abstract Algorithm methods.
API::MatrixWorkspace_sptr createOutputWorkspace()
Create output workspace.
std::vector< double > m_vecPeakParamValues
Definition: GeneratePeaks.h:90
bool m_useFuncParamWS
Flag to use parameter table workspace.
bool hasParameter(const API::IFunction_sptr &function, const std::string &paramname)
Check whether function has a certain parameter.
void generatePeaks(const std::map< specnum_t, std::vector< std::pair< double, API::IFunction_sptr > > > &functionmap, const API::MatrixWorkspace_sptr &dataWS)
Generate peaks in output data workspaces.
std::map< specnum_t, specnum_t > m_SpectrumMap
Spectrum map from full spectra workspace to partial spectra workspace.
Definition: GeneratePeaks.h:95
API::MatrixWorkspace_const_sptr inputWS
Input workspace (optional)
double m_maxChi2
Maximum chi-square to have peak generated.
API::MatrixWorkspace_sptr createDataWorkspace(std::vector< double > binparameters)
Create a Workspace2D (MatrixWorkspace) with given spectra and bin parameters.
std::vector< std::string > m_funcParameterNames
List of functions' parameters naems.
bool m_useRawParameter
Flag to indicate parameter table workspace containing raw parameters names.
std::vector< std::string > addFunctionParameterNames(const std::vector< std::string > &funcnames)
Add function parameter names to.
void getSpectraSet(const DataObjects::TableWorkspace_const_sptr &peakParmsWS)
Get set of spectra of the input table workspace Spectra is set to the column named 'spectrum'.
std::vector< double > binParameters
Binning parameters.
void processAlgProperties(std::string &peakfunctype, std::string &bkgdfunctype)
Process algorithm properties.
API::IPeakFunction_sptr getPeakFunction(const API::IFunction_sptr &infunction)
Get the IPeakFunction part in the input function.
API::IPeakFunction_sptr m_peakFunction
Peak function.
Definition: GeneratePeaks.h:84
void importPeaksFromTable(std::map< specnum_t, std::vector< std::pair< double, API::IFunction_sptr > > > &functionmap)
Import peak and background functions from table workspace.
void processTableColumnNames()
Process column names with peak parameter names.
bool m_newWSFromParent
Flag whether the new workspace is exactly as input.
bool m_useAutoBkgd
Flag to use automatic background (???)
bool m_genBackground
Flag to generate background.
std::vector< double > m_vecBkgdParamValues
Definition: GeneratePeaks.h:92
double m_numPeakWidth
Number of FWHM for peak to extend.
void importPeakFromVector(std::vector< std::pair< double, API::IFunction_sptr > > &functionmap)
Import peak and background function parameters from vector.
DataObjects::TableWorkspace_sptr m_funcParamWS
Parameter table workspace.
Support for a property that holds an array of values.
Definition: ArrayProperty.h:28
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
void information(const std::string &msg)
Logs at information level.
Definition: Logger.cpp:105
static T & Instance()
Return a reference to the Singleton instance, creating it if it does not already exist Creation is do...
std::shared_ptr< IPeakFunction > IPeakFunction_sptr
std::shared_ptr< IFunction > IFunction_sptr
shared pointer to the function base class
Definition: IFunction.h:732
std::shared_ptr< const Column > Column_const_sptr
Definition: Column.h:229
std::shared_ptr< MatrixWorkspace > MatrixWorkspace_sptr
shared pointer to the matrix workspace base class
std::shared_ptr< CompositeFunction > CompositeFunction_sptr
shared pointer to the composite function base class
std::shared_ptr< const TableWorkspace > TableWorkspace_const_sptr
shared pointer to Mantid::DataObjects::TableWorkspace (const version)
void split(const int A, int &S, int &V)
Split a number into the sign and positive value.
Definition: Acomp.cpp:42
int32_t specnum_t
Typedef for a spectrum Number.
Definition: IDTypes.h:16
@ Input
An input workspace.
Definition: Property.h:53
@ Output
An output workspace.
Definition: Property.h:54