Mantid
Loading...
Searching...
No Matches
LeBailFunction.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 +
12#include "MantidHistogramData/HistogramX.h"
13#include "MantidHistogramData/HistogramY.h"
14#include "MantidKernel/System.h"
15
16#include <sstream>
17#include <utility>
18
19#include <gsl/gsl_sf_erf.h>
20
21using namespace Mantid::API;
22using namespace Mantid::Kernel;
23using Mantid::HistogramData::HistogramY;
24
25using namespace std;
26
27const double NEG_DBL_MAX(-1. * DBL_MAX);
28
30namespace {
31const double PEAKRANGECONSTANT = 5.0;
32
33const string CHEBYSHEV_BACKGROUND("Chebyshev");
34const string POLYNOMIAL_BACKGROUND("Polynomial");
35const string FULLPROF_POLYNOMIAL_BACKGROUND("FullprofPolynomial");
36} // namespace
37
38// Get a reference to the logger
39Kernel::Logger g_log("LeBailFunction");
40
41//----------------------------------------------------------------------------------------------
44LeBailFunction::LeBailFunction(const std::string &peaktype) {
45 // Set initial values to some class variables
47 m_compsiteFunction = m_function;
48
49 m_numPeaks = 0;
50
51 m_isInputValue = false;
52 m_hasNewPeakValue = false;
53
54 // Peak type, validate and parameter name vectors
55 m_peakType = peaktype;
56 IFunction_sptr ifunc = FunctionFactory::Instance().createFunction(m_peakType);
57 if (!ifunc) {
58 stringstream errss;
59 errss << "Input peak type " << peaktype << " is not a recoganizable Mantid function.";
60 throw runtime_error(errss.str());
61 }
62 IPowderDiffPeakFunction_sptr peakfunc = std::dynamic_pointer_cast<IPowderDiffPeakFunction>(ifunc);
63 if (!peakfunc) {
64 stringstream errss;
65 errss << "Input peak type " << peaktype << " is not a IPowderDiffPeakFunction.";
66 throw runtime_error(errss.str());
67 }
68
69 m_peakParameterNameVec = peakfunc->getParameterNames();
72
73 // Peak parameter values
74 for (auto parname : m_peakParameterNameVec) {
75 m_functionParameters.emplace(parname, 0.0);
76 }
77
78 // Importing peak position tolerance
80 m_maxTOFPeakCentre = DBL_MAX;
81}
82
83//----------------------------------------------------------------------------------------------
86LeBailFunction::~LeBailFunction() = default;
87
88//----------------------------------------------------------------------------------------------
91API::IFunction_sptr LeBailFunction::getFunction() {
92 return m_compsiteFunction;
93 // return std::dynamic_pointer_cast<IFunction_sptr>(m_compsiteFunction);
94}
95
96//----------------------------------------------------------------------------------------------
104HistogramY LeBailFunction::function(const Mantid::HistogramData::HistogramX &xvalues, bool calpeaks,
105 bool calbkgd) const {
106
107 // Reset output elements to zero
108 std::vector<double> out(xvalues.size(), 0);
109 const auto &xvals = xvalues.rawData();
110
111 // Peaks
112 if (calpeaks) {
113 for (size_t ipk = 0; ipk < m_numPeaks; ++ipk) {
114 // Reset temporary vector for output
115 vector<double> temp(xvalues.size(), 0);
117 peak->function(temp, xvals);
118 transform(out.begin(), out.end(), temp.begin(), out.begin(), ::plus<double>());
119 }
120 }
121
122 // Background if required
123 if (calbkgd) {
124 if (!m_background) {
125 throw runtime_error("Must define background first!");
126 }
127
128 FunctionDomain1DVector domain(xvals);
129 FunctionValues values(domain);
130 g_log.information() << "Background function (in LeBailFunction): " << m_background->asString() << ".\n";
131 m_background->function(domain, values);
132 size_t numpts = out.size();
133 for (size_t i = 0; i < numpts; ++i)
134 out[i] += values[i];
135 }
136
137 return HistogramY(out);
138}
139
142HistogramY LeBailFunction::calPeak(size_t ipk, const std::vector<double> &xvalues, size_t ySize) const {
143
144 if (ipk >= m_numPeaks) {
145 stringstream errss;
146 errss << "Try to calculate peak indexed " << ipk << ". But number of peaks = " << m_numPeaks;
147 g_log.error(errss.str());
148 throw runtime_error(errss.str());
149 }
150
151 std::vector<double> out(ySize, 0);
153 peak->function(out, xvalues);
154 return HistogramY(out);
155}
156
157//----------------------------------------------------------------------------------------------
161bool LeBailFunction::hasProfileParameter(const std::string &paramname) {
162 auto fiter = lower_bound(m_orderedProfileParameterNames.cbegin(), m_orderedProfileParameterNames.cend(), paramname);
163
164 bool found = true;
165 if (fiter == m_orderedProfileParameterNames.end()) {
166 // End of the vector
167 found = false;
168 } else {
169 // Middle of vector
170 string matchparname = *fiter;
171 if (matchparname != paramname)
172 found = false;
173 }
174
175 return found;
176}
177
178//----------------------------------------------------------------------------------------------
183bool LeBailFunction::isParameterValid(double maxfwhm) const {
184 // Re-calculate peak parameter if there is some modification
185 if (m_hasNewPeakValue) {
187 }
188
189 // Check whether each peak has valid value
190 bool arevalid = true;
191 for (size_t i = 0; i < m_numPeaks; ++i) {
193 bool isvalid = peak->isPhysical();
194 if (isvalid && maxfwhm >= 0)
195 isvalid = peak->fwhm() < maxfwhm;
196 if (!isvalid) {
197 arevalid = false;
198
199 int h, k, l;
200 peak->getMillerIndex(h, k, l);
201 g_log.information() << "Peak [" << h << ", " << k << ", " << l << "] @ TOF = " << peak->centre()
202 << " has unphysical parameters or unreasonable large FWHM"
203 << ".\n";
204 break;
205 }
206 }
207
208 return arevalid;
209}
210
211//----------------------------------------------------------------------------------------------
214void LeBailFunction::calculatePeakParameterValues() const {
215 for (size_t i = 0; i < m_numPeaks; ++i) {
217 peak->calculateParameters(false);
218 }
219
220 m_hasNewPeakValue = false;
221}
222
223//----------------------------------------------------------------------------------------------
229void LeBailFunction::setPeakCentreTolerance(double peakpostol, double tofmin, double tofmax) {
230 // m_usePeakPosTol = true;
231 m_minTOFPeakCentre = tofmin - peakpostol;
232 m_maxTOFPeakCentre = tofmax + peakpostol;
233}
234
235//----------------------------------------------------------------------------------------------
239void LeBailFunction::addPeaks(std::vector<std::vector<int>> peakhkls) {
240 // Prerequisit
241 if (!m_isInputValue)
242 throw runtime_error("Client must set up profile parameter vlaues by calling "
243 "setProfileParameterValues() first! ");
244
245 // Add peaks
246 for (size_t ipk = 0; ipk < peakhkls.size(); ++ipk) {
247 vector<int> hkl = peakhkls[ipk];
248
249 // Check input Miller Index
250 if (hkl.size() != 3) {
251 stringstream errss;
252 errss << "Error of " << ipk << "-th input Miller Index. It has " << peakhkls[ipk].size()
253 << " items, but not required 3 items.";
254 g_log.error(errss.str());
255 throw runtime_error(errss.str());
256 }
257
258 // Generate new peak
259 int h = hkl[0];
260 int k = hkl[1];
261 int l = hkl[2];
263 if (!newpeak) {
264 g_log.error("Unable to generate peak. ");
265 throw runtime_error("Unable to generate peak.");
266 }
267
268 double tofh = newpeak->centre();
269 if (tofh < m_minTOFPeakCentre || tofh > m_maxTOFPeakCentre) {
270 g_log.information() << "Peak " << h << ", " << k << ", " << l << " 's centre is at TOF = " << tofh
271 << ", which is out of user specified boundary (" << m_minTOFPeakCentre << ", "
272 << m_maxTOFPeakCentre << "). "
273 << ".\n";
274 } else {
275 double dsp = newpeak->getPeakParameter("d_h");
276
277 // Add new peak to all related data storage
278 m_vecPeaks.emplace_back(newpeak);
279 // FIXME - Refining lattice size is not considered here!
280 m_dspPeakVec.emplace_back(dsp, newpeak);
281 m_mapHKLPeak.emplace(hkl, newpeak);
282 }
283 }
284
285 m_numPeaks = m_vecPeaks.size();
286
287 g_log.information() << "Total " << m_numPeaks << " after trying to add " << peakhkls.size() << " peaks. \n";
288} // END of addPeaks()
289
290//----------------------------------------------------------------------------------------------
296IPowderDiffPeakFunction_sptr LeBailFunction::generatePeak(int h, int k, int l) {
298 IPowderDiffPeakFunction_sptr peak = std::dynamic_pointer_cast<IPowderDiffPeakFunction>(f);
299
300 peak->setMillerIndex(h, k, l);
301 for (const auto &parname : m_peakParameterNameVec) {
302 double parvalue = m_functionParameters[parname];
303 peak->setParameter(parname, parvalue);
304 }
305
306 return peak;
307}
308
309//----------------------------------------------------------------------------------------------
324bool LeBailFunction::calculatePeaksIntensities(const vector<double> &vecX, const vector<double> &vecY,
325 vector<double> &vec_summedpeaks) {
326 // Clear inputs
327 std::fill(vec_summedpeaks.begin(), vec_summedpeaks.end(), 0.0);
328
329 // Divide peaks into groups from peak's parameters
330 vector<vector<pair<double, IPowderDiffPeakFunction_sptr>>> peakgroupvec;
331 vector<IPowderDiffPeakFunction_sptr> outboundpeakvec;
332 double xmin = vecX.front();
333 double xmax = vecX.back();
334 groupPeaks(peakgroupvec, outboundpeakvec, xmin, xmax);
335
336 // Calculate each peak's intensity and set
337 bool allpeakheightsphysical = true;
338 for (size_t ig = 0; ig < peakgroupvec.size(); ++ig) {
339 g_log.debug() << "[Fx351] Calculate peaks heights for (peak) group " << ig
340 << " : number of peaks = " << peakgroupvec[ig].size() << "\n";
341
342 bool peakheightsphysical = calculateGroupPeakIntensities(peakgroupvec[ig], vecX, vecY, vec_summedpeaks);
343
344 if (!peakheightsphysical)
345 allpeakheightsphysical = false;
346 }
347
348 // Set zero to all peaks out of boundary
349 for (const auto &peak : outboundpeakvec) {
350 peak->setHeight(0.);
351 }
352
353 return allpeakheightsphysical;
354}
355
356//----------------------------------------------------------------------------------------------
367bool LeBailFunction::calculateGroupPeakIntensities(vector<pair<double, IPowderDiffPeakFunction_sptr>> peakgroup,
368 const vector<double> &vecX, const vector<double> &vecY,
369 vector<double> &vec_summedpeaks) {
370 // Check input peaks group and sort peak by d-spacing
371 if (peakgroup.empty()) {
372 throw runtime_error("Programming error such that input peak group cannot be empty!");
373 } else {
374 g_log.debug() << "[Fx155] Peaks group size = " << peakgroup.size() << "\n";
375 }
376 if (peakgroup.size() > 1)
377 sort(peakgroup.begin(), peakgroup.end());
378
379 // Check input vector validity
380 if (vec_summedpeaks.size() != vecY.size()) {
381 stringstream errss;
382 errss << "Input vector 'allpeaksvalues' has wrong size = " << vec_summedpeaks.size()
383 << " != data workspace Y's size = " << vecY.size();
384 g_log.error(errss.str());
385 throw runtime_error(errss.str());
386 }
387
388 // Check boundary
389 IPowderDiffPeakFunction_sptr leftpeak = peakgroup[0].second;
390 double leftbound = leftpeak->centre() - PEAKRANGECONSTANT * leftpeak->fwhm();
391 if (leftbound < vecX.front()) {
392 stringstream msg;
393 int h, k, l;
394 leftpeak->getMillerIndex(h, k, l);
395 msg << "Peak group (containing " << peakgroup.size() << " peaks) has its left boundary (TOF = " << leftbound
396 << ") out side of input data workspace's left boundary (" << vecX.front()
397 << "). Accuracy of its peak intensity might be affected. "
398 << "Group's left boundary is determined by its leftmost peak (" << h << ", " << k << ", " << l
399 << ") at TOF = " << leftpeak->centre() << " with FWHM = " << leftpeak->fwhm() << ". ";
400
401 g_log.information(msg.str());
402
403 leftbound = vecX[0] + 0.1;
404 }
405 IPowderDiffPeakFunction_sptr rightpeak = peakgroup.back().second;
406 double rightbound = rightpeak->centre() + PEAKRANGECONSTANT * rightpeak->fwhm();
407 if (rightbound > vecX.back()) {
408 stringstream msg;
409 msg << "Peak group's right boundary " << rightbound << " is out side of "
410 << "input data workspace's right bound (" << vecX.back()
411 << ")! Accuracy of its peak intensity might be affected. ";
412
413 g_log.information(msg.str());
414
415 rightbound = vecX.back() - 0.1;
416 }
417
418 // Determine calculation range to input workspace: [ileft, iright)
419 vector<double>::const_iterator cviter;
420
421 cviter = lower_bound(vecX.begin(), vecX.end(), leftbound);
422 size_t ileft = static_cast<size_t>(cviter - vecX.begin());
423 if (ileft > 0)
424 --ileft;
425
426 cviter = lower_bound(vecX.begin(), vecX.end(), rightbound);
427 size_t iright = static_cast<size_t>(cviter - vecX.begin());
428 if (iright <= vecX.size() - 1)
429 ++iright;
430
431 size_t ndata = iright - ileft;
432 if (ileft >= iright) {
433 stringstream errss;
434 errss << "[Calcualte Peak Intensity] Group range is unphysical. iLeft = " << ileft << ", iRight = " << iright
435 << "; Number of peaks = " << peakgroup.size() << "; Left boundary = " << leftbound
436 << ", Right boundary = " << rightbound << "; Left peak FWHM = " << leftpeak->fwhm()
437 << ", Right peak FWHM = " << rightpeak->fwhm();
438 for (size_t ipk = 0; ipk < peakgroup.size(); ++ipk) {
439 IPowderDiffPeakFunction_sptr thispeak = peakgroup[ipk].second;
440 errss << "Peak " << ipk << ": d_h = " << peakgroup[ipk].first << ", TOF_h = " << thispeak->centre()
441 << ", FWHM = " << thispeak->fwhm() << "\n";
442 vector<string> peakparamnames = thispeak->getParameterNames();
443 for (auto &peakparamname : peakparamnames) {
444 errss << "\t" << peakparamname << " = " << thispeak->getParameter(peakparamname) << "\n";
445 }
446 }
447
448 g_log.error(errss.str());
449 throw runtime_error(errss.str());
450 }
451
452 // Generate a subset of vecX and vecY to calculate peak intensities
453 vector<double> datax(vecX.begin() + ileft, vecX.begin() + iright);
454 vector<double> datay(vecY.begin() + ileft, vecY.begin() + iright);
455 if (datax.size() != ndata) {
456 stringstream errmsg;
457 errmsg << "Impossible: Partial peak data size = " << datax.size() << " != ndata = " << ndata;
458 g_log.error(errmsg.str());
459 throw runtime_error(errmsg.str());
460 }
461 g_log.debug() << "[DBx356] Number of data points = " << ndata << " index from " << ileft << " to " << iright
462 << "; Size(datax, datay) = " << datax.size() << "\n";
463
464 // Prepare to integrate dataY to calculate peak intensity
465 vector<double> sumYs(ndata, 0.0);
466 size_t numPeaks(peakgroup.size());
467 vector<vector<double>> peakvalues(numPeaks);
468
469 // Integrage peak by peak
470 bool datavalueinvalid = false;
471 for (size_t ipk = 0; ipk < numPeaks; ++ipk) {
472 // calculate peak function value. Peak height should be set to a non-zero
473 // value
474 IPowderDiffPeakFunction_sptr peak = peakgroup[ipk].second;
475 peak->setHeight(1.0);
476 vector<double> localpeakvalue(ndata, 0.0);
477 peak->function(localpeakvalue, datax);
478
479 // check data
480 size_t numbadpts(0);
481 vector<double>::const_iterator localpeakvalue_end = localpeakvalue.end();
482 for (auto it = localpeakvalue.begin(); it != localpeakvalue_end; ++it) {
483 if ((*it != 0.) && (*it < NEG_DBL_MAX || *it > DBL_MAX)) {
484 numbadpts++;
485 }
486 }
487
488 // report the problem and/or integrate data
489 if (numbadpts == 0) {
490 // Data is fine. Integrate them all
491 for (size_t i = 0; i < ndata; ++i) {
492 // If value is physical
493 sumYs[i] += localpeakvalue[i];
494 }
495 } else {
496 // Report the problem
497 int h, k, l;
498 peak->getMillerIndex(h, k, l);
499 stringstream warnss;
500 warnss << "Peak (" << h << ", " << k << ", " << l << ") @ TOF = " << peak->centre() << " has " << numbadpts
501 << " data points, "
502 << "whose values exceed limit (i.e., not physical). ";
503 g_log.debug(warnss.str());
504 datavalueinvalid = true;
505 }
506 peakvalues[ipk].assign(localpeakvalue.begin(), localpeakvalue.end());
507 } // For All peaks
508
509 // Calculate intensity of all peaks
510 bool peakheightsphysical = !datavalueinvalid;
511 if (peakheightsphysical) {
512 for (size_t ipk = 0; ipk < peakgroup.size(); ++ipk) {
513 IPowderDiffPeakFunction_sptr peak = peakgroup[ipk].second;
514 double intensity = 0.0;
515
516 for (size_t i = 0; i < ndata; ++i) {
517 double temp;
518 if (sumYs[i] > 1.0E-5) {
519 // Reasonable non-zero value
520 double peaktogroupratio = peakvalues[ipk][i] / sumYs[i];
521 temp = datay[i] * peaktogroupratio;
522 } else {
523 // SumY too smaller
524 temp = 0.0;
525 }
526 double deltax;
527 if (i == 0)
528 deltax = datax[1] - datax[0];
529 else
530 deltax = datax[i] - datax[i - 1];
531
532 intensity += temp * deltax;
533 } // for data points
534
535 if (intensity != intensity) {
536 // Unphysical intensity: NaN
537 intensity = 0.0;
538 peakheightsphysical = false;
539
540 int h, k, l;
541 peak->getMillerIndex(h, k, l);
542 g_log.warning() << "Peak (" << h << ", " << k << ", " << l << ") has unphysical intensity = NaN!\n";
543
544 } else if (intensity <= -DBL_MAX || intensity >= DBL_MAX) {
545 // Unphysical intensity: NaN
546 intensity = 0.0;
547 peakheightsphysical = false;
548
549 int h, k, l;
550 peak->getMillerIndex(h, k, l);
551 g_log.warning() << "Peak (" << h << ", " << k << ", " << l << ") has unphysical intensity = Infty!\n";
552 } else if (intensity < 0.0) {
553 // No negative intensity
554 g_log.debug() << "[Fx134] Set peak @ " << peak->centre() << "'s intensity to 0.0 instead of " << intensity
555 << ".\n";
556 intensity = 0.0;
557 }
558 g_log.debug() << "[Fx407] Peak @ " << peak->centre() << ": Set Intensity = " << intensity << "\n";
559 peak->setHeight(intensity);
560
561 // Add peak's value to peaksvalues
562 for (size_t i = ileft; i < iright; ++i) {
563 vec_summedpeaks[i] += (intensity * peakvalues[ipk][i - ileft]);
564 }
565
566 } // ENDFOR each peak
567 }
568
569 return peakheightsphysical;
570}
571
572//----------------------------------------------------------------------------------------------
581void LeBailFunction::setPeakParameters(const IPowderDiffPeakFunction_sptr &peak, const map<string, double> &parammap,
582 double peakheight, bool setpeakheight) {
583 UNUSED_ARG(peak);
584 UNUSED_ARG(parammap);
585 UNUSED_ARG(peakheight);
586 UNUSED_ARG(setpeakheight);
587 throw runtime_error("Requiring update flag: peak value changed and etc.");
588 /*
589 // FIXME - The best solution for speeding is to have a set of peak
590 parameter listed in the order
591 // of peak function's parameters' indexed. Then no need to do
592 search anymore.
593
594 // 1. Prepare, sort parameters by name
595 std::map<std::string, double>::iterator pit;
596 vector<string> peakparamnames = peak->getParameterNames();
597
598 // 2. Apply parameters values to peak function
599 for (pit = parammap.begin(); pit != parammap.end(); ++pit)
600 {
601 // a) Check whether the parameter is a peak parameter
602 std::string parname = pit->first;
603 std::vector<std::string>::iterator ifind =
604 std::find(peakparamnames.begin(), peakparamnames.end(), parname);
605
606 // b) Set parameter value
607 if (ifind == peakparamnames.end())
608 {
609 // If not a peak profile parameter, skip
610 g_log.debug() << "Parameter '" << parname << "' in input parameter
611 table workspace "
612 << "is not for peak function " << peak->name() << ".\n";
613 }
614 else
615 {
616 // Set value
617 double value = pit->second;
618 peak->setParameter(parname, value);
619 g_log.debug() << "LeBailFit Set " << parname << "= " << value << "\n";
620 }
621 } // ENDFOR: parameter iterator
622
623 // 3. Peak height
624 if (setpeakheight)
625 peak->setHeight(peakheight);
626
627 return;*/
628}
629
630//----------------------------------------------------------------------------------------------
641void LeBailFunction::setProfileParameterValues(map<std::string, double> parammap) {
642 const double MINDIFF = 1.0E-10;
643
644 map<std::string, double>::iterator inpiter, curiter;
645
646 size_t numpars = m_peakParameterNameVec.size();
647 for (size_t i = 0; i < numpars; ++i) {
648 string &parname = m_peakParameterNameVec[i];
649
650 // Find iterator of this parameter in input parammap
651 inpiter = parammap.find(parname);
652 if (inpiter != parammap.end()) {
653 // Find iterator to parameter value in class' parameter map (parameter is
654 // found in input map)
655 curiter = m_functionParameters.find(parname);
656 if (curiter == m_functionParameters.end()) {
657 stringstream errmsg;
658 errmsg << "Parameter " << parname << " is in parameter name list, but not in profile "
659 << "parameter map. It violates the programming logic.";
660 g_log.error(errmsg.str());
661 throw runtime_error(errmsg.str());
662 }
663
664 // Set value if difference is large
665 double curvalue = curiter->second;
666 double newvalue = inpiter->second;
667 bool localnewvalue = false;
668 if (fabs(curvalue - newvalue) > MINDIFF) {
669 curiter->second = newvalue;
670 m_hasNewPeakValue = true;
671 localnewvalue = true;
672 }
673
674 // Set new value to each peak
675 if (!localnewvalue)
676 continue;
677
678 // Set new parameter to each peak
679 for (size_t ipk = 0; ipk < m_numPeaks; ++ipk) {
681 peak->setParameter(i, newvalue);
682 }
683 } // If parameter name is a profile parameter
684 else {
685 g_log.debug() << "Parameter " << parname << " is not a profile parameter. Length of string = " << parname.size()
686 << "\n";
687 }
688 } // ENDFOR [All profile parameter]
689
690 // Set the flag to indicate that client has input parameters
692 m_isInputValue = true;
693}
694
695//----------------------------------------------------------------------------------------------
703void LeBailFunction::groupPeaks(vector<vector<pair<double, IPowderDiffPeakFunction_sptr>>> &peakgroupvec,
704 vector<IPowderDiffPeakFunction_sptr> &outboundpeakvec, double xmin, double xmax) {
705 // Sort peaks
706 if (m_numPeaks > 1) {
707 sort(m_dspPeakVec.begin(), m_dspPeakVec.end());
708 } else if (m_numPeaks == 0) {
709 std::stringstream errmsg;
710 errmsg << "Group peaks: No peak is found in the peak vector. ";
711 g_log.error() << errmsg.str() << "\n";
712 throw std::runtime_error(errmsg.str());
713 }
714
715 // Set up starting value
716 peakgroupvec.clear();
717 outboundpeakvec.clear();
718 vector<pair<double, IPowderDiffPeakFunction_sptr>> peakgroup; // one group of peaks
719 size_t ipk = 0;
720
721 // Group peaks from low-d to high-d
722 bool outbound = true;
723 while (outbound && ipk < m_numPeaks) {
724 // Group peaks out of lower boundary to a separate vector of peaks
726 if (peak->centre() <= xmin) {
727 // Add peak
728 outboundpeakvec.emplace_back(peak);
729 ipk += 1;
730 } else {
731 // Get out of while loop if peak is in bound
732 outbound = false;
733 }
734 }
735
736 bool inbound = true;
737 while (inbound && ipk < m_numPeaks) {
738 // Group peaks in the boundary
739 IPowderDiffPeakFunction_sptr thispeak = m_dspPeakVec[ipk].second;
740
741 if (thispeak->centre() < xmax) {
742 // Peak is in the boundary still
743
744 // add peak to CURRENT peak group
745 peakgroup.emplace_back(m_dspPeakVec[ipk]);
746
747 if (ipk < m_numPeaks - 1) {
748 // Any peak but not the last (rightmost) peak
749
750 // test whether next peak will be in a different group
751 IPowderDiffPeakFunction_sptr rightpeak = m_dspPeakVec[ipk + 1].second;
752
753 double thispeak_rightbound = thispeak->centre() + PEAKRANGECONSTANT * thispeak->fwhm();
754 double rightpeak_leftbound = rightpeak->centre() - PEAKRANGECONSTANT * rightpeak->fwhm();
755
756 if (thispeak_rightbound < rightpeak_leftbound) {
757 // this peak and its right peak are well separated.
758 // finish this group by swapping values
759 peakgroupvec.emplace_back(std::move(peakgroup));
760 peakgroup = {};
761 } else {
762 // this peak and its right peak are close enough to be in same group.
763 // do nothing
764 ;
765 }
766 } else {
767 // Rightmost peak. Finish the current peak
768 peakgroupvec.emplace_back(peakgroup);
769 }
770
771 ++ipk;
772 } // still in bound
773 else {
774 // Peak is get out of boundary
775 inbound = false;
776 g_log.information() << "[Fx301] Group peak: peak @ " << thispeak->centre() << " causes grouping "
777 << "peak over at maximum TOF = " << xmax << ".\n";
778
779 if (!peakgroup.empty()) {
780 peakgroupvec.emplace_back(peakgroup);
781 }
782 } // FIRST out of boundary
783 } // ENDWHILE
784
785 while (ipk < m_numPeaks) {
786 // Group peaks out of uppper boundary to a separate vector of peaks
787 outboundpeakvec.emplace_back(m_dspPeakVec[ipk].second);
788 ipk += 1;
789 }
790
791 g_log.debug() << "[Calculate Peak Intensity]: Number of Peak Groups = " << peakgroupvec.size() << "\n";
792}
793
794//----------------------------------------------------------------------------------------------
805void LeBailFunction::addBackgroundFunction(const string &backgroundtype, const unsigned int &order,
806 const std::vector<std::string> &vecparnames,
807 const std::vector<double> &vecparvalues, double startx, double endx) {
808 // Check
809 if (backgroundtype != POLYNOMIAL_BACKGROUND && backgroundtype != CHEBYSHEV_BACKGROUND &&
810 backgroundtype != FULLPROF_POLYNOMIAL_BACKGROUND) {
811 stringstream warnss;
812 warnss << "Cliet specified background type " << backgroundtype << " may not be supported properly.";
813 g_log.warning(warnss.str());
814 }
815 if (vecparnames.size() != vecparvalues.size())
816 throw runtime_error("Input parameter names and parameter values are not matched. ");
817
818 g_log.information() << "Add background: type = " << backgroundtype << ", order = " << order
819 << ", number of parameters/attributes = " << vecparnames.size() << "\n";
820
821 // Create background function from factory
822 auto background = FunctionFactory::Instance().createFunction(backgroundtype);
823 m_background = std::dynamic_pointer_cast<Functions::BackgroundFunction>(background);
824
825 // Set order and initialize
826 m_background->setAttributeValue("n", static_cast<int>(order));
827 m_background->initialize();
828
829 // Set parameters & attribute
830 size_t numpars = vecparnames.size();
831 for (size_t i = 0; i < numpars; ++i) {
832 const string &parname = vecparnames[i];
833 if (parname != "Bkpos")
834 m_background->setParameter(parname, vecparvalues[i]);
835 else if (backgroundtype == FULLPROF_POLYNOMIAL_BACKGROUND)
836 m_background->setAttributeValue("Bkpos", vecparvalues[i]);
837 else
838 throw runtime_error("Bkpos should not be in the parameter list. ");
839 }
840
841 if (backgroundtype == CHEBYSHEV_BACKGROUND) {
842 if (startx > 0.)
843 m_background->setAttributeValue("StartX", startx);
844 if (endx > 0.)
845 m_background->setAttributeValue("EndX", endx);
846 }
847}
848
849//----------------------------------------------------------------------------------------------
855void LeBailFunction::setFitProfileParameter(const string &paramname, double minvalue, double maxvalue) {
856 // Make ties in composition function
857 for (size_t ipk = 1; ipk < m_numPeaks; ++ipk) {
858 stringstream ss1, ss2;
859 ss1 << "f" << (ipk - 1) << "." << paramname;
860 ss2 << "f" << ipk << "." << paramname;
861 string tiepart1 = ss1.str();
862 string tiepart2 = ss2.str();
863 m_compsiteFunction->tie(tiepart1, tiepart2);
864 g_log.debug() << "LeBailFunction::Fit(Tie) / " << tiepart1 << " / " << tiepart2 << " /\n";
865 }
866
867 // Set contrains of the parameter on any of the tied parameter.
868 std::stringstream parss;
869 parss << "f0." << paramname;
870 string parnamef0 = parss.str();
871 auto bc = std::make_unique<Constraints::BoundaryConstraint>(m_compsiteFunction.get(), parnamef0, minvalue, maxvalue);
872 m_compsiteFunction->addConstraint(std::move(bc));
873}
874
875//----------------------------------------------------------------------------------------------
880void LeBailFunction::fixPeakParameter(const string &paramname, double paramvalue) {
881 for (size_t ipk = 0; ipk < m_numPeaks; ++ipk) {
882 stringstream ss1, ss2;
883 ss1 << "f" << ipk << "." << paramname;
884 ss2 << paramvalue;
885 string tiepart1 = ss1.str();
886 string tievalue = ss2.str();
887 m_compsiteFunction->tie(tiepart1, tievalue);
888
889 g_log.debug() << "Set up tie | " << tiepart1 << " <---> " << tievalue << " | \n";
890
891 // FIXME & TODO: Make a map between peak parameter name and index. And use
892 // fix() to replace tie
893 /*-- Code prepared to replace the existing block
894 ThermalNeutronBk2BkExpConvPVoigt_sptr thispeak = m_dspPeaks[ipk].second;
895 size_t iparam = findIndex(thispeak, funcparam.name);
896 thispeak->fix(iparam);
897 --*/
898
899 } // For each peak
900}
901
902//----------------------------------------------------------------------------------------------
905void LeBailFunction::fixBackgroundParameters() {
906 size_t numbkgdparams = m_background->nParams();
907
908 for (size_t iparam = 0; iparam < numbkgdparams; ++iparam)
909 m_background->fix(iparam);
910}
911
912//----------------------------------------------------------------------------------------------
915void LeBailFunction::setFixPeakHeights() {
916 for (size_t ipk = 0; ipk < m_numPeaks; ++ipk) {
917 // a. Get peak height
918 IPowderDiffPeakFunction_sptr thispeak = m_dspPeakVec[ipk].second;
919 thispeak->fix(0);
920 } // For each peak
921}
922
923//----------------------------------------------------------------------------------------------
927void LeBailFunction::setPeakHeights(const std::vector<double> &inheights) {
928 UNUSED_ARG(inheights);
929 throw runtime_error("It is not implemented properly.");
930 /*
931 if (inheights.size() != heights.size())
932 {
933 g_log.error() << "Input number of peaks (height) is not same as peaks. "
934 << '\n';
935 throw std::logic_error("Input number of peaks (height) is not same as
936 peaks. ");
937 }
938
939 for (size_t ih = 0; ih < inheights.size(); ++ih)
940 heights[ih] = inheights[ih];
941
942 return;*/
943}
944
945//----------------------------------------------------------------------------------------------
948IPowderDiffPeakFunction_sptr LeBailFunction::getPeak(size_t peakindex) {
949 if (peakindex >= m_numPeaks) {
950 stringstream errmsg;
951 errmsg << "Try to access peak " << peakindex << " out of range [0, " << m_numPeaks << ").";
952 g_log.error(errmsg.str());
953 throw runtime_error(errmsg.str());
954 }
955
956 IPowderDiffPeakFunction_sptr rpeak = m_vecPeaks[peakindex];
957
958 return rpeak;
959}
960
961//----------------------------------------------------------------------------------------------
964double LeBailFunction::getPeakParameter(std::vector<int> hkl, const std::string &parname) const {
965 // Search peak in map
966 map<vector<int>, IPowderDiffPeakFunction_sptr>::const_iterator fiter;
967 fiter = m_mapHKLPeak.find(hkl);
968 if (fiter == m_mapHKLPeak.end()) {
969 stringstream errss;
970 errss << "Peak with Miller index (" << hkl[0] << ", " << hkl[1] << "," << hkl[2]
971 << ") does not exist in Le Bail function.";
972 g_log.error(errss.str());
973 throw runtime_error(errss.str());
974 }
975
976 IPowderDiffPeakFunction_sptr peak = fiter->second;
977
978 double parvalue = getPeakParameterValue(peak, parname);
979
980 return parvalue;
981}
982
983//----------------------------------------------------------------------------------------------
986double LeBailFunction::getPeakParameter(size_t index, const std::string &parname) const {
987 if (index >= m_numPeaks) {
988 stringstream errss;
989 errss << "getPeakParameter() tries to reach a peak with index " << index << ", which is out of range " << m_numPeaks
990 << "/" << m_vecPeaks.size() << ".";
991 g_log.error(errss.str());
992 throw std::runtime_error(errss.str());
993 }
994
996 double value = getPeakParameterValue(peak, parname);
997
998 return value;
999}
1000
1001//----------------------------------------------------------------------------------------------
1006double LeBailFunction::getPeakParameterValue(const API::IPowderDiffPeakFunction_sptr &peak,
1007 const std::string &parname) const {
1008 // Locate the category of the parameter name
1009 auto vsiter = lower_bound(m_orderedProfileParameterNames.cbegin(), m_orderedProfileParameterNames.cend(), parname);
1010
1011 bool found = true;
1012 if (vsiter == m_orderedProfileParameterNames.end()) {
1013 // End of vector
1014 found = false;
1015 } else {
1016 // Middle of vector. But no match
1017 string matchparname = *vsiter;
1018 if (parname != matchparname)
1019 found = false;
1020 }
1021
1022 // Get parameter
1023 double parvalue;
1024 if (found) {
1025 // It is a native peak parameter
1026 parvalue = peak->getParameter(parname);
1027 } else {
1028 // It is a calculated peak parameter
1029 parvalue = peak->getPeakParameter(parname);
1030 }
1031
1032 return parvalue;
1033}
1034
1035//----------------------------------------------------------------------------------------------
1038double LeBailFunction::getPeakMaximumValue(std::vector<int> hkl, const std::vector<double> &xvalues, size_t &ix) {
1039 // Search peak in map
1040 map<vector<int>, IPowderDiffPeakFunction_sptr>::const_iterator fiter;
1041 fiter = m_mapHKLPeak.find(hkl);
1042 if (fiter == m_mapHKLPeak.end()) {
1043 stringstream errss;
1044 errss << "Peak with Miller index (" << hkl[0] << ", " << hkl[1] << "," << hkl[2]
1045 << ") does not exist in Le Bail function.";
1046 g_log.error(errss.str());
1047 throw runtime_error(errss.str());
1048 }
1049
1050 IPowderDiffPeakFunction_sptr peak = fiter->second;
1051
1052 double maxvalue = peak->getMaximumValue(xvalues, ix);
1053
1054 return maxvalue;
1055}
1056
1057} // namespace Mantid::CurveFitting::Algorithms
double value
The value of the point.
Definition: FitMW.cpp:51
std::map< DeltaEMode::Type, std::string > index
Definition: DeltaEMode.cpp:19
#define fabs(x)
Definition: Matrix.cpp:22
const double NEG_DBL_MAX
#define UNUSED_ARG(x)
Function arguments are sometimes unused in certain implmentations but are required for documentation ...
Definition: System.h:64
A composite function is a function containing other functions.
Implements FunctionDomain1D with its own storage in form of a std::vector.
A class to store values calculated by a function.
std::vector< std::string > m_peakParameterNameVec
Name of peak parameter names (be same as the order in IPowderDiffPeakFunction)
void calculatePeakParameterValues() const
Calculate all peaks' parameter value.
void groupPeaks(std::vector< std::vector< std::pair< double, API::IPowderDiffPeakFunction_sptr > > > &peakgroupvec, std::vector< API::IPowderDiffPeakFunction_sptr > &outboundpeakvec, double xmin, double xmax)
Group close peaks together.
std::vector< std::pair< double, API::IPowderDiffPeakFunction_sptr > > m_dspPeakVec
Vector of pair <peak position in d-space, Peak> sortable.
bool calculateGroupPeakIntensities(std::vector< std::pair< double, API::IPowderDiffPeakFunction_sptr > > peakgroup, const std::vector< double > &vecX, const std::vector< double > &vecY, std::vector< double > &vec_summedpeaks)
Calculate the peaks intensities in same group.
Functions::BackgroundFunction_sptr m_background
Background function.
API::IPowderDiffPeakFunction_sptr generatePeak(int h, int k, int l)
Generate a peak with parameter set by.
std::vector< API::IPowderDiffPeakFunction_sptr > m_vecPeaks
Vector of all peaks.
double getPeakParameterValue(const API::IPowderDiffPeakFunction_sptr &peak, const std::string &parname) const
Retrieve peak's parameter. may be native or calculated.
API::CompositeFunction_sptr m_compsiteFunction
Composite functions for all peaks and background.
std::vector< std::string > m_orderedProfileParameterNames
Ordered profile parameter names for search.
std::map< std::vector< int >, API::IPowderDiffPeakFunction_sptr > m_mapHKLPeak
Vector of all peak's Miller indexes.
std::map< std::string, double > m_functionParameters
Parameters.
The Logger class is in charge of the publishing messages from the framework through various channels.
Definition: Logger.h:52
void debug(const std::string &msg)
Logs at debug level.
Definition: Logger.cpp:114
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< IPowderDiffPeakFunction > IPowderDiffPeakFunction_sptr
Kernel::Logger g_log("ExperimentInfo")
static logger object
std::shared_ptr< IFunction > IFunction_sptr
shared pointer to the function base class
Definition: IFunction.h:732
std::shared_ptr< CompositeFunction > CompositeFunction_sptr
shared pointer to the composite function base class
static constexpr double h
Planck constant in J*s.
STL namespace.