Mantid
Loading...
Searching...
No Matches
ProcessBackground.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/Axis.h"
10#include "MantidAPI/TableRow.h"
19#include "MantidKernel/System.h"
21
22#include <boost/algorithm/string.hpp>
23#include <boost/algorithm/string/predicate.hpp>
24#include <boost/algorithm/string/split.hpp>
25#include <utility>
26
27using namespace Mantid;
28
29using namespace Mantid::API;
30
31using namespace Mantid::Kernel;
32
33using namespace Mantid::DataObjects;
34
35using namespace Mantid::CurveFitting;
36
37using namespace HistogramData;
38
39using namespace std;
40
42
43using namespace CurveFitting;
44
46
47//----------------------------------------------------------------------------------------------
51 : m_dataWS(), m_outputWS(), m_wsIndex(-1), m_lowerBound(DBL_MAX), m_upperBound(DBL_MIN), m_bkgdType(),
52 m_numFWHM(-1.) {}
53
54//----------------------------------------------------------------------------------------------
58 // Input and output Workspace
59 declareProperty(std::make_unique<WorkspaceProperty<Workspace2D>>("InputWorkspace", "Anonymous", Direction::Input),
60 "Name of the output workspace containing the processed background.");
61
62 // Workspace index
63 declareProperty("WorkspaceIndex", 0, "Workspace index for the input workspaces.");
64
65 // Output workspace
66 declareProperty(std::make_unique<WorkspaceProperty<Workspace2D>>("OutputWorkspace", "", Direction::Output),
67 "Output workspace containing processed background");
68
69 // Function Options
70 std::vector<std::string> options{"SelectBackgroundPoints", "RemovePeaks", "DeleteRegion", "AddRegion"};
71
72 auto validator = std::make_shared<Kernel::StringListValidator>(options);
73 declareProperty("Options", "RemovePeaks", validator, "Name of the functionality realized by this algorithm.");
74
75 // Boundary
76 declareProperty("LowerBound", Mantid::EMPTY_DBL(), "Lower boundary of the data to have background processed.");
77 declareProperty("UpperBound", Mantid::EMPTY_DBL(), "Upper boundary of the data to have background processed.");
78
79 auto refwsprop = std::make_unique<WorkspaceProperty<Workspace2D>>("ReferenceWorkspace", "", Direction::Input,
81 declareProperty(std::move(refwsprop), "Name of the workspace containing the data "
82 "required by function AddRegion.");
83 setPropertySettings("ReferenceWorkspace", std::make_unique<VisibleWhenProperty>("Options", IS_EQUAL_TO, "AddRegion"));
84
85 // Optional Function Type
86 std::vector<std::string> bkgdtype{"Polynomial", "Chebyshev"};
87 auto bkgdvalidator = std::make_shared<Kernel::StringListValidator>(bkgdtype);
88 declareProperty("BackgroundType", "Polynomial", bkgdvalidator,
89 "Type of the background. Options include Polynomial and Chebyshev.");
90 setPropertySettings("BackgroundType",
91 std::make_unique<VisibleWhenProperty>("Options", IS_EQUAL_TO, "SelectBackgroundPoints"));
92
93 vector<string> funcoptions{"N/A", "FitGivenDataPoints", "UserFunction"};
94 auto fovalidator = std::make_shared<StringListValidator>(funcoptions);
95 declareProperty("SelectionMode", "N/A", fovalidator,
96 "If choise is UserFunction, background will be selected by "
97 "an input background "
98 "function. Otherwise, background function will be fitted "
99 "from user's input data points.");
100 setPropertySettings("SelectionMode",
101 std::make_unique<VisibleWhenProperty>("Options", IS_EQUAL_TO, "SelectBackgroundPoints"));
102
103 declareProperty("BackgroundOrder", 0, "Order of polynomial or chebyshev background. ");
104 setPropertySettings("BackgroundOrder",
105 std::make_unique<VisibleWhenProperty>("Options", IS_EQUAL_TO, "SelectBackgroundPoints"));
106 setPropertySettings("BackgroundOrder",
107 std::make_unique<VisibleWhenProperty>("SelectionMode", IS_EQUAL_TO, "FitGivenDataPoints"));
108
109 // User input background points for "SelectBackground"
110 auto arrayproperty = std::make_unique<Kernel::ArrayProperty<double>>("BackgroundPoints");
111 declareProperty(std::move(arrayproperty), "Vector of doubles, each of which is the "
112 "X-axis value of the background point "
113 "selected by user.");
114 setPropertySettings("BackgroundPoints",
115 std::make_unique<VisibleWhenProperty>("Options", IS_EQUAL_TO, "SelectBackgroundPoints"));
116 setPropertySettings("BackgroundPoints",
117 std::make_unique<VisibleWhenProperty>("SelectionMode", IS_EQUAL_TO, "FitGivenDataPoints"));
118
119 declareProperty(std::make_unique<WorkspaceProperty<TableWorkspace>>("BackgroundTableWorkspace", "", Direction::Input,
121 "Name of the table workspace containing background "
122 "parameters for mode SelectBackgroundPoints.");
123 setPropertySettings("BackgroundTableWorkspace",
124 std::make_unique<VisibleWhenProperty>("Options", IS_EQUAL_TO, "SelectBackgroundPoints"));
125 setPropertySettings("BackgroundTableWorkspace",
126 std::make_unique<VisibleWhenProperty>("SelectionMode", IS_EQUAL_TO, "UserFunction"));
127
128 // Mode to select background
129 vector<string> pointsselectmode{"All Background Points", "Input Background Points Only"};
130 auto modevalidator = std::make_shared<StringListValidator>(pointsselectmode);
131 declareProperty("BackgroundPointSelectMode", "All Background Points", modevalidator,
132 "Mode to select background points. ");
133 setPropertySettings("BackgroundPointSelectMode",
134 std::make_unique<VisibleWhenProperty>("Options", IS_EQUAL_TO, "SelectBackgroundPoints"));
135 setPropertySettings("BackgroundPointSelectMode",
136 std::make_unique<VisibleWhenProperty>("SelectionMode", IS_EQUAL_TO, "FitGivenDataPoints"));
137
138 // Background tolerance
139 declareProperty("NoiseTolerance", 1.0, "Tolerance of noise range. ");
140 setPropertySettings("NoiseTolerance",
141 std::make_unique<VisibleWhenProperty>("Options", IS_EQUAL_TO, "SelectBackgroundPoints"));
142
143 // Background tolerance
144 declareProperty("NegativeNoiseTolerance", EMPTY_DBL(), "Tolerance of noise range for negative number. ");
145 setPropertySettings("NegativeNoiseTolerance",
146 std::make_unique<VisibleWhenProperty>("Options", IS_EQUAL_TO, "SelectBackgroundPoints"));
147
148 // Optional output workspace
150 std::make_unique<WorkspaceProperty<Workspace2D>>("UserBackgroundWorkspace", "_dummy01", Direction::Output),
151 "Output workspace containing fitted background from points "
152 "specified by users.");
153 setPropertySettings("UserBackgroundWorkspace",
154 std::make_unique<VisibleWhenProperty>("Options", IS_EQUAL_TO, "SelectBackgroundPoints"));
155
156 // Optional output workspace
157 declareProperty(std::make_unique<WorkspaceProperty<TableWorkspace>>("OutputBackgroundParameterWorkspace", "_dummy02",
159 "Output parameter table workspace containing the background fitting "
160 "result. ");
161 setPropertySettings("OutputBackgroundParameterWorkspace",
162 std::make_unique<VisibleWhenProperty>("Options", IS_EQUAL_TO, "SelectBackgroundPoints"));
163
164 // Output background type.
165 std::vector<std::string> outbkgdtype{"Polynomial", "Chebyshev"};
166 auto outbkgdvalidator = std::make_shared<Kernel::StringListValidator>(outbkgdtype);
167 declareProperty("OutputBackgroundType", "Polynomial", outbkgdvalidator,
168 "Type of background to fit with selected background points.");
169 setPropertySettings("OutputBackgroundType",
170 std::make_unique<VisibleWhenProperty>("Options", IS_EQUAL_TO, "SelectBackgroundPoints"));
171
172 // Output background type.
173 declareProperty("OutputBackgroundOrder", 6, "Order of background to fit with selected background points.");
174 setPropertySettings("OutputBackgroundOrder",
175 std::make_unique<VisibleWhenProperty>("Options", IS_EQUAL_TO, "SelectBackgroundPoints"));
176
177 // Peak table workspac for "RemovePeaks"
178 declareProperty(std::make_unique<WorkspaceProperty<TableWorkspace>>("BraggPeakTableWorkspace", "", Direction::Input,
180 "Name of table workspace containing peaks' parameters. ");
181 setPropertySettings("BraggPeakTableWorkspace",
182 std::make_unique<VisibleWhenProperty>("Options", IS_EQUAL_TO, "RemovePeaks"));
183
184 // Number of FWHM to have peak removed
185 declareProperty("NumberOfFWHM", 1.0, "Number of FWHM to as the peak region to have peak removed. ");
186 setPropertySettings("NumberOfFWHM", std::make_unique<VisibleWhenProperty>("Options", IS_EQUAL_TO, "RemovePeaks"));
187}
188
189//----------------------------------------------------------------------------------------------
193 // Process general properties
194 m_dataWS = this->getProperty("InputWorkspace");
195 if (!m_dataWS) {
196 g_log.error() << "Input Workspace cannot be obtained.\n";
197 throw std::invalid_argument("Input Workspace cannot be obtained.");
198 }
199
200 m_bkgdType = getPropertyValue("BackgroundType");
201
202 int intemp = getProperty("WorkspaceIndex");
203 if (intemp < 0)
204 throw std::invalid_argument("WorkspaceIndex is not allowed to be less than 0. ");
205 m_wsIndex = intemp;
206 if (m_wsIndex >= static_cast<int>(m_dataWS->getNumberHistograms()))
207 throw runtime_error("Workspace index is out of boundary.");
208
209 m_lowerBound = getProperty("LowerBound");
210 m_upperBound = getProperty("UpperBound");
212 m_lowerBound = m_dataWS->x(m_wsIndex).front();
214 m_upperBound = m_dataWS->x(m_wsIndex).back();
215
216 // 2. Do different work
217 std::string option = getProperty("Options");
218 if (option == "RemovePeaks") {
219 removePeaks();
220 } else if (option == "DeleteRegion") {
221 deleteRegion();
222 } else if (option == "AddRegion") {
223 addRegion();
224 } else if (option == "SelectBackgroundPoints") {
225
227 } else {
228 g_log.error() << "Option " << option << " is not supported. \n";
229 throw std::invalid_argument("Unsupported option. ");
230 }
231
232 // 3. Set output
233 setProperty("OutputWorkspace", m_outputWS);
234}
235
236//----------------------------------------------------------------------------------------------
240 // Dummy outputs to make it work with python script
241 setPropertyValue("UserBackgroundWorkspace", "dummy0");
242 Workspace2D_sptr dummyws =
243 std::dynamic_pointer_cast<Workspace2D>(WorkspaceFactory::Instance().create("Workspace2D", 1, 1, 1));
244 setProperty("UserBackgroundWorkspace", dummyws);
245
246 setPropertyValue("OutputBackgroundParameterWorkspace", "dummy1");
247 TableWorkspace_sptr dummytbws = std::make_shared<TableWorkspace>();
248 setProperty("OutputBackgroundParameterWorkspace", dummytbws);
249}
250
251//----------------------------------------------------------------------------------------------
255 // Check boundary
257 throw std::invalid_argument("Using DeleteRegion. Both LowerBound and "
258 "UpperBound must be specified.");
259 }
260 if (m_lowerBound >= m_upperBound) {
261 throw std::invalid_argument("Lower boundary cannot be equal or larger than upper boundary.");
262 }
263
264 const auto &dataX = m_dataWS->x(0);
265 const auto &dataY = m_dataWS->y(0);
266 const auto &dataE = m_dataWS->e(0);
267
268 // Find the dimensions of the region excluded by m_lowerBound and m_upperBound
269 std::vector<size_t> incIndexes;
270 for (size_t i = 0; i < dataY.size(); i++) {
271 if (dataX[i] < m_lowerBound || dataX[i] > m_upperBound) {
272 incIndexes.emplace_back(i);
273 }
274 }
275 size_t sizex = incIndexes.size();
276 size_t sizey = sizex;
277 if (dataX.size() > dataY.size()) {
278 sizex++;
279 }
280
281 // Create a new workspace with these dimensions and copy data from the defined
282 // region
283 API::MatrixWorkspace_sptr mws = API::WorkspaceFactory::Instance().create("Workspace2D", 1, sizex, sizey);
284 m_outputWS = std::dynamic_pointer_cast<DataObjects::Workspace2D>(mws);
285 m_outputWS->getAxis(0)->setUnit(m_dataWS->getAxis(0)->unit()->unitID());
286
287 for (size_t i = 0; i < sizey; i++) {
288 size_t index = incIndexes[i];
289 m_outputWS->mutableX(0)[i] = dataX[index];
290 m_outputWS->mutableY(0)[i] = dataY[index];
291 m_outputWS->mutableE(0)[i] = dataE[index];
292 }
293 if (sizex > sizey) {
294 m_outputWS->mutableX(0)[sizex - 1] = dataX.back();
295 }
296
297 // Set up dummies
299}
300
301//----------------------------------------------------------------------------------------------
305 // Check boundary, which are required
307 throw std::invalid_argument("Using AddRegion. Both LowerBound and UpperBound must be specified.");
308 }
309 if (m_lowerBound >= m_upperBound) {
310 throw std::invalid_argument("Lower boundary cannot be equal or larger than upper boundary.");
311 }
312
313 // Copy data to a set of vectors
314 const auto &vecX = m_dataWS->x(0);
315 const auto &vecY = m_dataWS->y(0);
316 const auto &vecE = m_dataWS->e(0);
317
318 std::vector<double> vx, vy, ve;
319 for (size_t i = 0; i < vecY.size(); ++i) {
320 double xtmp = vecX[i];
321 if (xtmp < m_lowerBound || xtmp > m_upperBound) {
322 vx.emplace_back(vecX[i]);
323 vy.emplace_back(vecY[i]);
324 ve.emplace_back(vecE[i]);
325 }
326 }
327
328 // Histogram
329 if (vecX.size() > vecY.size())
330 vx.emplace_back(vecX.back());
331
332 // Get access to reference workspace
333 DataObjects::Workspace2D_const_sptr refWS = getProperty("ReferenceWorkspace");
334 if (!refWS)
335 throw std::invalid_argument("ReferenceWorkspace is not given. ");
336
337 const auto &refX = refWS->x(0);
338 const auto &refY = refWS->y(0);
339 const auto &refE = refWS->e(0);
340
341 // 4. Insert the value of the reference workspace from lowerBoundary to
342 // upperBoundary
343 std::vector<double>::const_iterator refiter;
344 refiter = std::lower_bound(refX.begin(), refX.end(), m_lowerBound);
345 size_t sindex = size_t(refiter - refX.begin());
346 refiter = std::lower_bound(refX.begin(), refX.end(), m_upperBound);
347 size_t eindex = size_t(refiter - refX.begin());
348
349 for (size_t i = sindex; i < eindex; ++i) {
350 const double tmpx = refX[i];
351 const double tmpy = refY[i];
352 const double tmpe = refE[i];
353
354 // Locate the position of tmpx in the array to be inserted
355 auto newit = std::lower_bound(vx.begin(), vx.end(), tmpx);
356 size_t newindex = size_t(newit - vx.begin());
357
358 // insert tmpx, tmpy, tmpe by iterator
359 vx.insert(newit, tmpx);
360
361 newit = vy.begin() + newindex;
362 vy.insert(newit, tmpy);
363
364 newit = ve.begin() + newindex;
365 ve.insert(newit, tmpe);
366 }
367
368 // Check
369 for (auto it = vx.begin() + 1; it != vx.end(); ++it) {
370 if (*it <= *it - 1) {
371 g_log.error() << "The vector X with value inserted is not ordered incrementally\n";
372 throw std::runtime_error("Build new vector error!");
373 }
374 }
375
376 // Construct the new Workspace
377 m_outputWS = std::dynamic_pointer_cast<DataObjects::Workspace2D>(
378 API::WorkspaceFactory::Instance().create("Workspace2D", 1, vx.size(), vy.size()));
379 m_outputWS->getAxis(0)->setUnit(m_dataWS->getAxis(0)->unit()->unitID());
380 m_outputWS->setHistogram(0, Histogram(Points(vx), Counts(vy), CountStandardDeviations(ve)));
381
382 // Write out dummy output workspaces
384}
385
386//----------------------------------------------------------------------------------------------
387// Methods for selecting background points
388//----------------------------------------------------------------------------------------------
389
390//----------------------------------------------------------------------------------------------
394 // Select background points
395 string smode = getProperty("SelectionMode");
396 if (smode == "FitGivenDataPoints") {
398 } else if (smode == "UserFunction") {
400 } else {
401 throw runtime_error("N/A is not supported.");
402 }
403
404 // Fit the background points if output backgrond parameter workspace is given
405 // explicitly
406 string outbkgdparwsname = getPropertyValue("OutputBackgroundParameterWorkspace");
407 if (!outbkgdparwsname.empty() && outbkgdparwsname != "_dummy02") {
408 // Will fit the selected background
409 string bkgdfunctype = getPropertyValue("OutputBackgroundType");
410 fitBackgroundFunction(bkgdfunctype);
411 } else {
413 }
414
415 m_outputWS->getAxis(0)->setUnit(m_dataWS->getAxis(0)->unit()->unitID());
416}
417
418//----------------------------------------------------------------------------------------------
422 // Get special input properties
423 std::vector<double> bkgdpoints = getProperty("BackgroundPoints");
424 string mode = getProperty("BackgroundPointSelectMode");
425
426 // Construct background workspace for fit
427 std::vector<size_t> realIndexes;
428 const auto &vecX = m_dataWS->x(m_wsIndex);
429 const auto &vecY = m_dataWS->y(m_wsIndex);
430 const auto &vecE = m_dataWS->e(m_wsIndex);
431 for (size_t i = 0; i < bkgdpoints.size(); ++i) {
432 // Data range validation
433 double bkgdpoint = bkgdpoints[i];
434 if (bkgdpoint < vecX.front()) {
435 g_log.warning() << "Input background point " << bkgdpoint << " is out of lower boundary. "
436 << "Use X[0] = " << vecX.front() << " instead."
437 << "\n";
438 bkgdpoint = vecX.front();
439 } else if (bkgdpoint > vecX.back()) {
440 g_log.warning() << "Input background point " << bkgdpoint
441 << " is out of upper boundary. Use X[-1] = " << vecX.back() << " instead."
442 << "\n";
443 bkgdpoint = vecX.back();
444 }
445
446 // Find the index in
447 std::vector<double>::const_iterator it;
448 it = std::lower_bound(vecX.begin(), vecX.end(), bkgdpoint);
449 size_t index = size_t(it - vecX.begin());
450
451 g_log.debug() << "DBx502 Background Points " << i << " Index = " << index << " For TOF = " << bkgdpoints[i]
452 << " in [" << vecX[0] << ", " << vecX.back() << "] "
453 << "\n";
454
455 // Add index to list
456 realIndexes.emplace_back(index);
457
458 } // ENDFOR (i)
459
460 size_t wsSize = realIndexes.size();
461 DataObjects::Workspace2D_sptr bkgdWS = std::dynamic_pointer_cast<DataObjects::Workspace2D>(
462 API::WorkspaceFactory::Instance().create("Workspace2D", 1, wsSize, wsSize));
463 for (size_t i = 0; i < wsSize; ++i) {
464 size_t index = realIndexes[i];
465 bkgdWS->mutableX(0)[i] = vecX[index];
466 bkgdWS->mutableY(0)[i] = vecY[index];
467 bkgdWS->mutableE(0)[i] = vecE[index];
468 }
469
470 // Select background points according to mode
471 if (mode == "All Background Points") {
472 // Select (possibly) all background points
474 } else if (mode == "Input Background Points Only") {
475 // Use the input background points only
476 m_outputWS = bkgdWS;
477 } else {
478 stringstream errss;
479 errss << "Background select mode " << mode << " is not supported by ProcessBackground.";
480 g_log.error(errss.str());
481 throw runtime_error(errss.str());
482 }
483}
484
485//----------------------------------------------------------------------------------------------
489 // Process properties
491 TableWorkspace_sptr bkgdtablews = getProperty("BackgroundTableWorkspace");
492
493 // Set up background function from table
494 size_t numrows = bkgdtablews->rowCount();
495 map<string, double> parmap;
496 for (size_t i = 0; i < numrows; ++i) {
497 TableRow row = bkgdtablews->getRow(i);
498 string parname;
499 double parvalue;
500 row >> parname >> parvalue;
501 if (parname[0] == 'A')
502 parmap.emplace(parname, parvalue);
503 }
504
505 auto bkgdorder = static_cast<int>(parmap.size() - 1); // A0 - A(n) total n+1 parameters
506 bkgdfunc->setAttributeValue("n", bkgdorder);
507 for (auto &mit : parmap) {
508 string parname = mit.first;
509 double parvalue = mit.second;
510 bkgdfunc->setParameter(parname, parvalue);
511 }
512
513 // Filter out
515}
516
517//----------------------------------------------------------------------------------------------
521 // Get background type and create bakground function
523
524 int bkgdorder = getProperty("BackgroundOrder");
525 if (bkgdorder == 0)
526 g_log.warning("(Input) background function order is 0. It might not be "
527 "able to give a good estimation.");
528
529 bkgdfunction->setAttributeValue("n", bkgdorder);
530 bkgdfunction->initialize();
531
532 g_log.information() << "Input background points has " << bkgdWS->x(0).size() << " data points for fit " << bkgdorder
533 << "-th order " << bkgdfunction->name() << " (background) function" << bkgdfunction->asString()
534 << "\n";
535
536 // Fit input (a few) background pionts to get initial guess
538 try {
539 fit = this->createChildAlgorithm("Fit", 0.0, 0.2, true);
540 } catch (Exception::NotFoundError &) {
541 g_log.error() << "Requires CurveFitting library.\n";
542 throw;
543 }
544
545 double startx = m_lowerBound;
546 double endx = m_upperBound;
547 fit->setProperty("Function", std::dynamic_pointer_cast<API::IFunction>(bkgdfunction));
548 fit->setProperty("InputWorkspace", bkgdWS);
549 fit->setProperty("WorkspaceIndex", 0);
550 fit->setProperty("MaxIterations", 500);
551 fit->setProperty("StartX", startx);
552 fit->setProperty("EndX", endx);
553 fit->setProperty("Minimizer", "Levenberg-Marquardt");
554 fit->setProperty("CostFunction", "Least squares");
555
556 fit->executeAsChildAlg();
557
558 // Get fit result
559 // a) Status
560 std::string fitStatus = fit->getProperty("OutputStatus");
561 bool allowedfailure = (fitStatus.find("Changes") < fitStatus.size()) && (fitStatus.find("small") < fitStatus.size());
562 if (fitStatus != "success" && !allowedfailure) {
563 g_log.error() << "ProcessBackground: Fit Status = " << fitStatus << ". Not to update fit result\n";
564 throw std::runtime_error("Bad Fit " + fitStatus);
565 }
566
567 // b) check that chi2 got better
568 const double chi2 = fit->getProperty("OutputChi2overDoF");
569 g_log.information() << "Fit background: Fit Status = " << fitStatus << ", chi2 = " << chi2 << "\n";
570
571 // Filter and construct for the output workspace
572 Workspace2D_sptr outws = filterForBackground(bkgdfunction);
573
574 return outws;
575} // END OF FUNCTION
576
577//----------------------------------------------------------------------------------------------
582
583 if (backgroundtype == "Polynomial") {
584 bkgdfunction = std::dynamic_pointer_cast<Functions::BackgroundFunction>(std::make_shared<Functions::Polynomial>());
585 bkgdfunction->initialize();
586 } else if (backgroundtype == "Chebyshev") {
587 Chebyshev_sptr cheby = std::make_shared<Functions::Chebyshev>();
588 bkgdfunction = std::dynamic_pointer_cast<Functions::BackgroundFunction>(cheby);
589 bkgdfunction->initialize();
590
591 g_log.debug() << "[D] Chebyshev is set to range " << m_lowerBound << ", " << m_upperBound << "\n";
592 bkgdfunction->setAttributeValue("StartX", m_lowerBound);
593 bkgdfunction->setAttributeValue("EndX", m_upperBound);
594 } else {
595 stringstream errss;
596 errss << "Background of type " << backgroundtype << " is not supported. ";
597 g_log.error(errss.str());
598 throw std::invalid_argument(errss.str());
599 }
600
601 return bkgdfunction;
602}
603
604//----------------------------------------------------------------------------------------------
608 double posnoisetolerance = getProperty("NoiseTolerance");
609 double negnoisetolerance = getProperty("NegativeNoiseTolerance");
610 if (isEmpty(negnoisetolerance))
611 negnoisetolerance = posnoisetolerance;
612
613 // Calcualte theoretical values
614 const auto &x = m_dataWS->x(m_wsIndex);
615 API::FunctionDomain1DVector domain(x.rawData());
616 API::FunctionValues values(domain);
617 bkgdfunction->function(domain, values);
618
619 g_log.information() << "Function used to select background points : " << bkgdfunction->asString() << "\n";
620
621 // Optional output
622 string userbkgdwsname = getPropertyValue("UserBackgroundWorkspace");
623 if (userbkgdwsname.empty())
624 throw runtime_error("In mode SelectBackgroundPoints, "
625 "UserBackgroundWorkspace must be given!");
626
627 size_t sizex = domain.size();
628 size_t sizey = values.size();
629 MatrixWorkspace_sptr visualws =
630 std::dynamic_pointer_cast<MatrixWorkspace>(WorkspaceFactory::Instance().create("Workspace2D", 4, sizex, sizey));
631 for (size_t i = 0; i < sizex; ++i) {
632 for (size_t j = 0; j < 4; ++j) {
633 visualws->mutableX(j)[i] = domain[i];
634 }
635 }
636 for (size_t i = 0; i < sizey; ++i) {
637 visualws->mutableY(0)[i] = values[i];
638 visualws->mutableY(1)[i] = m_dataWS->y(m_wsIndex)[i] - values[i];
639 visualws->mutableY(2)[i] = posnoisetolerance;
640 visualws->mutableY(3)[i] = -negnoisetolerance;
641 }
642 setProperty("UserBackgroundWorkspace", visualws);
643
644 // Filter for background
645 std::vector<size_t> selectedIndexes;
646 for (size_t i = 0; i < domain.size(); ++i) {
647 double purey = visualws->y(1)[i];
648 if (purey < posnoisetolerance && purey > -negnoisetolerance) {
649 selectedIndexes.emplace_back(i);
650 }
651 }
652 size_t wsSize = selectedIndexes.size();
653 g_log.information() << "Found " << wsSize << " background points out of " << m_dataWS->x(m_wsIndex).size()
654 << " total data points. "
655 << "\n";
656
657 // Build new workspace for OutputWorkspace
658 size_t nspec = 3;
659 Workspace2D_sptr outws = std::dynamic_pointer_cast<DataObjects::Workspace2D>(
660 API::WorkspaceFactory::Instance().create("Workspace2D", nspec, wsSize, wsSize));
661 for (size_t i = 0; i < wsSize; ++i) {
662 size_t index = selectedIndexes[i];
663 for (size_t j = 0; j < nspec; ++j)
664 outws->mutableX(j)[i] = domain[index];
665 outws->mutableY(0)[i] = m_dataWS->y(m_wsIndex)[index];
666 outws->mutableE(0)[i] = m_dataWS->e(m_wsIndex)[index];
667 }
668
669 return outws;
670}
671
672//----------------------------------------------------------------------------------------------
675void ProcessBackground::fitBackgroundFunction(const std::string &bkgdfunctiontype) {
676 // Get background type and create bakground function
677 BackgroundFunction_sptr bkgdfunction = createBackgroundFunction(bkgdfunctiontype);
678
679 int bkgdorder = getProperty("OutputBackgroundOrder");
680 bkgdfunction->setAttributeValue("n", bkgdorder);
681
682 if (bkgdfunctiontype == "Chebyshev") {
683 double xmin = m_outputWS->x(0).front();
684 double xmax = m_outputWS->x(0).back();
685 g_log.information() << "Chebyshev Fit range: " << xmin << ", " << xmax << "\n";
686 bkgdfunction->setAttributeValue("StartX", xmin);
687 bkgdfunction->setAttributeValue("EndX", xmax);
688 }
689
690 g_log.information() << "Fit selected background " << bkgdfunctiontype << " to data workspace with "
691 << m_outputWS->getNumberHistograms() << " spectra."
692 << "\n";
693
694 // Fit input (a few) background pionts to get initial guess
696 try {
697 fit = this->createChildAlgorithm("Fit", 0.9, 1.0, true);
698 } catch (Exception::NotFoundError &) {
699 g_log.error() << "Requires CurveFitting library.\n";
700 throw;
701 }
702
703 g_log.information() << "Fitting background function: " << bkgdfunction->asString() << "\n";
704
705 double startx = m_lowerBound;
706 double endx = m_upperBound;
707 fit->setProperty("Function", std::dynamic_pointer_cast<API::IFunction>(bkgdfunction));
708 fit->setProperty("InputWorkspace", m_outputWS);
709 fit->setProperty("WorkspaceIndex", 0);
710 fit->setProperty("MaxIterations", 500);
711 fit->setProperty("StartX", startx);
712 fit->setProperty("EndX", endx);
713 fit->setProperty("Minimizer", "Levenberg-MarquardtMD");
714 fit->setProperty("CostFunction", "Least squares");
715
716 fit->executeAsChildAlg();
717
718 // Get fit status and chi^2
719 std::string fitStatus = fit->getProperty("OutputStatus");
720 bool allowedfailure = (fitStatus.find("Changes") < fitStatus.size()) && (fitStatus.find("small") < fitStatus.size());
721 if (fitStatus != "success" && !allowedfailure) {
722 g_log.error() << "ProcessBackground: Fit Status = " << fitStatus << ". Not to update fit result\n";
723 throw std::runtime_error("Bad Fit " + fitStatus);
724 }
725
726 const double chi2 = fit->getProperty("OutputChi2overDoF");
727 g_log.information() << "Fit background: Fit Status = " << fitStatus << ", chi2 = " << chi2 << "\n";
728
729 // Get out the parameter names
730 API::IFunction_sptr funcout = fit->getProperty("Function");
731 TableWorkspace_sptr outbkgdparws = std::make_shared<TableWorkspace>();
732 outbkgdparws->addColumn("str", "Name");
733 outbkgdparws->addColumn("double", "Value");
734
735 TableRow typerow = outbkgdparws->appendRow();
736 typerow << bkgdfunctiontype << 0.;
737
738 vector<string> parnames = funcout->getParameterNames();
739 size_t nparam = funcout->nParams();
740 for (size_t i = 0; i < nparam; ++i) {
741 TableRow newrow = outbkgdparws->appendRow();
742 newrow << parnames[i] << funcout->getParameter(i);
743 }
744
745 TableRow chi2row = outbkgdparws->appendRow();
746 chi2row << "Chi-square" << chi2;
747
748 g_log.information() << "Set table workspace (#row = " << outbkgdparws->rowCount()
749 << ") to OutputBackgroundParameterTable. "
750 << "\n";
751 setProperty("OutputBackgroundParameterWorkspace", outbkgdparws);
752
753 // Set output workspace
754 const auto &vecX = m_outputWS->x(0);
755 const auto &vecY = m_outputWS->y(0);
756 FunctionDomain1DVector domain(vecX.rawData());
757 FunctionValues values(domain);
758
759 funcout->function(domain, values);
760
761 auto &dataModel = m_outputWS->mutableY(1);
762 auto &dataDiff = m_outputWS->mutableY(2);
763 for (size_t i = 0; i < dataModel.size(); ++i) {
764 dataModel[i] = values[i];
765 dataDiff[i] = vecY[i] - dataModel[i];
766 }
767}
768
769//----------------------------------------------------------------------------------------------
770// Remove peaks
771//----------------------------------------------------------------------------------------------
772
773//----------------------------------------------------------------------------------------------
777 // Get input
778 TableWorkspace_sptr peaktablews = getProperty("BraggPeakTableWorkspace");
779 if (!peaktablews)
780 throw runtime_error("Option RemovePeaks requires input to BgraggPeaTablekWorkspace.");
781
782 m_numFWHM = getProperty("NumberOfFWHM");
783 if (m_numFWHM <= 0.)
784 throw runtime_error("NumberOfFWHM must be larger than 0. ");
785
786 // Remove peaks
787 RemovePeaks remove;
788 remove.setup(peaktablews);
790
791 // Dummy outputs
793}
794
795//----------------------------------------------------------------------------------------------
798void RemovePeaks::setup(const TableWorkspace_sptr &peaktablews) {
799 // Parse table workspace
801
802 // Check
803 if (m_vecPeakCentre.size() != m_vecPeakFWHM.size())
804 throw runtime_error("Number of peak centres and FWHMs are different!");
805 else if (m_vecPeakCentre.empty())
806 throw runtime_error("There is not any peak entry in input table workspace.");
807}
808
809//----------------------------------------------------------------------------------------------
813 // Check
814 if (m_vecPeakCentre.empty())
815 throw runtime_error("RemovePeaks has not been setup yet. ");
816
817 // Initialize vectors
818 const auto &vecX = dataws->x(wsindex);
819 const auto &vecY = dataws->y(wsindex);
820 const auto &vecE = dataws->e(wsindex);
821
822 size_t sizex = vecX.size();
823 vector<bool> vec_useX(sizex, true);
824
825 // Exclude regions
826 size_t numbkgdpoints = excludePeaks(vecX.rawData(), vec_useX, m_vecPeakCentre, m_vecPeakFWHM, numfwhm);
827 size_t numbkgdpointsy = numbkgdpoints;
828 size_t sizey = vecY.size();
829 if (sizex > sizey)
830 --numbkgdpointsy;
831
832 // Construct output workspace
833 Workspace2D_sptr outws = std::dynamic_pointer_cast<Workspace2D>(
834 WorkspaceFactory::Instance().create("Workspace2D", 1, numbkgdpoints, numbkgdpointsy));
835 outws->getAxis(0)->setUnit(dataws->getAxis(0)->unit()->unitID());
836 auto &outX = outws->mutableX(0);
837 auto &outY = outws->mutableY(0);
838 auto &outE = outws->mutableE(0);
839 size_t index = 0;
840 for (size_t i = 0; i < sizex; ++i) {
841 if (vec_useX[i]) {
842 if (index >= numbkgdpoints)
843 throw runtime_error("Programming logic error (1)");
844 outX[index] = vecX[i];
845 ++index;
846 }
847 }
848 index = 0;
849 for (size_t i = 0; i < sizey; ++i) {
850 if (vec_useX[i]) {
851 if (index >= numbkgdpointsy)
852 throw runtime_error("Programming logic error (2)");
853 outY[index] = vecY[i];
854 outE[index] = vecE[i];
855 ++index;
856 }
857 }
858
859 return outws;
860}
861
862//----------------------------------------------------------------------------------------------
865void RemovePeaks::parsePeakTableWorkspace(const TableWorkspace_sptr &peaktablews, vector<double> &vec_peakcentre,
866 vector<double> &vec_peakfwhm) {
867 // Get peak table workspace information
868 vector<string> colnames = peaktablews->getColumnNames();
869 int index_centre = -1;
870 int index_fwhm = -1;
871 for (int i = 0; i < static_cast<int>(colnames.size()); ++i) {
872 string colname = colnames[i];
873 if (colname == "TOF_h")
874 index_centre = i;
875 else if (colname == "FWHM")
876 index_fwhm = i;
877 }
878
879 if (index_centre < 0 || index_fwhm < 0) {
880 throw runtime_error("Input Bragg peak table workspace does not have TOF_h and/or FWHM");
881 }
882
883 // Get values
884 size_t numrows = peaktablews->rowCount();
885 vec_peakcentre.resize(numrows, 0.);
886 vec_peakfwhm.resize(numrows, 0.);
887
888 for (size_t i = 0; i < numrows; ++i) {
889 double centre = peaktablews->cell<double>(i, index_centre);
890 double fwhm = peaktablews->cell<double>(i, index_fwhm);
891 vec_peakcentre[i] = centre;
892 vec_peakfwhm[i] = fwhm;
893 }
894}
895
896//----------------------------------------------------------------------------------------------
899size_t RemovePeaks::excludePeaks(vector<double> v_inX, vector<bool> &v_useX, vector<double> v_centre,
900 vector<double> v_fwhm, double num_fwhm) {
901 // Validate
902 if (v_centre.size() != v_fwhm.size())
903 throw runtime_error("Input different number of peak centres and fwhm.");
904 if (v_inX.size() != v_useX.size())
905 throw runtime_error("Input differetn number of vec X and flag X.");
906
907 // Flag peak regions
908 size_t numpeaks = v_centre.size();
909 for (size_t i = 0; i < numpeaks; ++i) {
910 // Define boundary
911 double centre = v_centre[i];
912 double fwhm = v_fwhm[i];
913 double xmin = centre - num_fwhm * fwhm;
914 double xmax = centre + num_fwhm * fwhm;
915
916 vector<double>::iterator viter;
917 int i_min, i_max;
918
919 // Locate index in v_inX
920 if (xmin <= v_inX.front())
921 i_min = 0;
922 else if (xmin >= v_inX.back())
923 i_min = static_cast<int>(v_inX.size()) - 1;
924 else {
925 viter = lower_bound(v_inX.begin(), v_inX.end(), xmin);
926 i_min = static_cast<int>(viter - v_inX.begin());
927 }
928
929 if (xmax <= v_inX.front())
930 i_max = 0;
931 else if (xmax >= v_inX.back())
932 i_max = static_cast<int>(v_inX.size()) - 1;
933 else {
934 viter = lower_bound(v_inX.begin(), v_inX.end(), xmax);
935 i_max = static_cast<int>(viter - v_inX.begin());
936 }
937
938 // Flag the excluded region
939 for (int excluded = i_min; excluded <= i_max; ++excluded)
940 v_useX[excluded] = false;
941 }
942
943 return std::count(v_useX.cbegin(), v_useX.cend(), true);
944}
945
946} // namespace Mantid::CurveFitting::Functions
#define DECLARE_ALGORITHM(classname)
Definition: Algorithm.h:576
std::map< DeltaEMode::Type, std::string > index
Definition: DeltaEMode.cpp:19
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
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
void setPropertyValue(const std::string &name, const std::string &value) override
Set the value of a property by string N.B.
Definition: Algorithm.cpp:1975
static bool isEmpty(const NumT toCheck)
checks that the value was not set by users, uses the value in empty double/int.
Implements FunctionDomain1D with its own storage in form of a std::vector.
size_t size() const override
Return the number of arguments in the domain.
A class to store values calculated by a function.
size_t size() const
Return the number of values.
TableRow represents a row in a TableWorkspace.
Definition: TableRow.h:39
A property class for workspaces.
ProcessBackground : Process background obtained from LeBailFit.
DataObjects::Workspace2D_sptr autoBackgroundSelection(const DataObjects::Workspace2D_sptr &bkgdWS)
Select background points automatically.
void fitBackgroundFunction(const std::string &bkgdfunctiontype)
Fit background function.
DataObjects::Workspace2D_const_sptr m_dataWS
void setupDummyOutputWSes()
Set up dummy output optional workspaces.
void removePeaks()
Remove peaks in a certain region.
void selectFromGivenFunction()
Select background points (main)
BackgroundFunction_sptr createBackgroundFunction(const std::string &backgroundtype)
Create a background function from input properties.
double m_numFWHM
Number of FWHM of range of peak to be removed.
DataObjects::Workspace2D_sptr filterForBackground(const BackgroundFunction_sptr &bkgdfunction)
Filter non-background data points out and create a background workspace.
void addRegion()
Add a certain region from a reference workspace.
void deleteRegion()
Remove a certain region from input workspace.
void selectFromGivenXValues()
Select background points (main)
void parsePeakTableWorkspace(const DataObjects::TableWorkspace_sptr &peaktablews, std::vector< double > &vec_peakcentre, std::vector< double > &vec_peakfwhm)
Parse peak centre and FWHM from a table workspace.
void setup(const DataObjects::TableWorkspace_sptr &peaktablews)
Set up: parse peak workspace to vectors.
size_t excludePeaks(std::vector< double > v_inX, std::vector< bool > &v_useX, std::vector< double > v_centre, std::vector< double > v_fwhm, double num_fwhm)
Exclude peak regions.
DataObjects::Workspace2D_sptr removePeaks(const API::MatrixWorkspace_const_sptr &dataws, int wsindex, double numfwhm)
Remove peaks from a input workspace.
Exception for when an item is not found in a collection.
Definition: Exception.h:145
IPropertyManager * setProperty(const std::string &name, const T &value)
Templated method to set the value of a PropertyWithValue.
void setPropertySettings(const std::string &name, std::unique_ptr< IPropertySettings > settings)
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< IAlgorithm > IAlgorithm_sptr
shared pointer to Mantid::API::IAlgorithm
std::shared_ptr< const MatrixWorkspace > MatrixWorkspace_const_sptr
shared pointer to the matrix workspace base class (const version)
std::shared_ptr< IFunction > IFunction_sptr
shared pointer to the function base class
Definition: IFunction.h:732
std::shared_ptr< MatrixWorkspace > MatrixWorkspace_sptr
shared pointer to the matrix workspace base class
std::shared_ptr< Chebyshev > Chebyshev_sptr
Definition: Chebyshev.h:55
std::shared_ptr< BackgroundFunction > BackgroundFunction_sptr
std::shared_ptr< Workspace2D > Workspace2D_sptr
shared pointer to Mantid::DataObjects::Workspace2D
std::shared_ptr< TableWorkspace > TableWorkspace_sptr
shared pointer to Mantid::DataObjects::TableWorkspace
std::shared_ptr< const Workspace2D > Workspace2D_const_sptr
shared pointer to Mantid::DataObjects::Workspace2D (const version)
std::unique_ptr< T > create(const P &parent, const IndexArg &indexArg, const HistArg &histArg)
This is the create() method that all the other create() methods call.
Helper class which provides the Collimation Length for SANS instruments.
constexpr double EMPTY_DBL() noexcept
Returns what we consider an "empty" double within a property.
Definition: EmptyValues.h:43
STL namespace.
@ Input
An input workspace.
Definition: Property.h:53
@ Output
An output workspace.
Definition: Property.h:54