diff --git a/Code/Mantid/Framework/Algorithms/src/FindPeakBackground.cpp b/Code/Mantid/Framework/Algorithms/src/FindPeakBackground.cpp
index 5068c9a99ef1d36257800c8e6cba6b516e8cf296..fdda95a8c531eca3277e0a6c46aee48a5f577fe8 100644
--- a/Code/Mantid/Framework/Algorithms/src/FindPeakBackground.cpp
+++ b/Code/Mantid/Framework/Algorithms/src/FindPeakBackground.cpp
@@ -76,9 +76,8 @@ namespace Algorithms
auto inwsprop = new WorkspaceProperty<MatrixWorkspace>("InputWorkspace", "Anonymous", Direction::Input);
declareProperty(inwsprop, "Name of input MatrixWorkspace that contains peaks.");
- declareProperty(new ArrayProperty<int>("WorkspaceIndices"), "Optional: enter a comma-separated list of the "
- "workspace indices to have peak and background separated. "
- "Default is to calculate for all spectra.");
+ declareProperty("WorkspaceIndex", EMPTY_INT(), "workspace indices to have peak and background separated. "
+ "No default is taken. ");
declareProperty("SigmaConstant", 1.0, "Multiplier of standard deviations of the variance for convergence of "
"peak elimination. Default is 1.0. ");
@@ -109,32 +108,27 @@ namespace Algorithms
{
// 1. Get input and validate
MatrixWorkspace_const_sptr inpWS = getProperty("InputWorkspace");
- std::vector<int> inpwsindex = getProperty("WorkspaceIndices");
+ int inpwsindex = getProperty("WorkspaceIndex");
std::vector<double> m_vecFitWindows = getProperty("FitWindow");
m_backgroundType = getPropertyValue("BackgroundType");
double k = getProperty("SigmaConstant");
- bool separateall = false;
- if (inpwsindex.size() == 0)
+ if (inpwsindex < 0 || inpwsindex >= static_cast<int>(inpWS->getNumberHistograms()))
{
- separateall = true;
+ stringstream errss;
+ errss << "Input workspace " << inpWS->name() << " has " << inpWS->getNumberHistograms()
+ << " spectra. Input workspace index " << inpwsindex << " is out of boundary. ";
+ throw runtime_error(errss.str());
}
+ if (isEmpty(inpwsindex))
+ throw runtime_error("WorkspaceIndex must be given. ");
// 2. Generate output
- size_t numspec;
- if (separateall)
- {
- numspec = inpWS->getNumberHistograms();
- }
- else
- {
- numspec = inpwsindex.size();
- }
- size_t sizex = inpWS->readX(0).size();
- size_t sizey = inpWS->readY(0).size();
+ const MantidVec& inpX = inpWS->readX(inpwsindex);
+ size_t sizex = inpWS->readX(inpwsindex).size();
+ size_t sizey = inpWS->readY(inpwsindex).size();
size_t n = sizey;
size_t l0 = 0;
- const MantidVec& inpX = inpWS->readX(0);
if (m_vecFitWindows.size() > 1)
{
@@ -152,145 +146,121 @@ namespace Algorithms
m_outPeakTableWS->addColumn("double", "bkg0");
m_outPeakTableWS->addColumn("double", "bkg1");
m_outPeakTableWS->addColumn("double", "bkg2");
- for( size_t i = 0; i < numspec; ++i )
- m_outPeakTableWS->appendRow();
+
+ m_outPeakTableWS->appendRow();
// 3. Get Y values
- Progress prog(this, 0, 1.0, numspec);
- PARALLEL_FOR2(inpWS, m_outPeakTableWS)
- for (int i = 0; i < static_cast<int>(numspec); ++i)
+ Progress prog(this, 0, 1.0, 1);
+
+ // Find background
+
+ const MantidVec& inpY = inpWS->readY(inpwsindex);
+
+ double Ymean, Yvariance, Ysigma;
+ MantidVec maskedY;
+ MantidVec::const_iterator in = std::min_element(inpY.begin(), inpY.end());
+ double bkg0 = inpY[in - inpY.begin()];
+ for (size_t l = l0; l < n; ++l)
+ {
+ maskedY.push_back(inpY[l]-bkg0);
+ }
+ MantidVec mask(n-l0,0.0);
+ double xn = static_cast<double>(n-l0);
+ do
+ {
+ Statistics stats = getStatistics(maskedY);
+ Ymean = stats.mean;
+ Yvariance = stats.standard_deviation * stats.standard_deviation;
+ Ysigma = std::sqrt((moment4(maskedY,n-l0,Ymean)-(xn-3.0)/(xn-1.0) * Yvariance)/xn);
+ MantidVec::const_iterator it = std::max_element(maskedY.begin(), maskedY.end());
+ const size_t pos = it - maskedY.begin();
+ maskedY[pos] = 0;
+ mask[pos] = 1.0;
+ }
+ while (std::abs(Ymean-Yvariance) > k * Ysigma);
+
+ if(n-l0 > 5)
{
- PARALLEL_START_INTERUPT_REGION
- // a) figure out wsindex
- size_t wsindex;
- if (separateall)
+ // remove single outliers
+ if (mask[1] == mask[2] && mask[2] == mask[3])
+ mask[0] = mask[1];
+ if (mask[0] == mask[2] && mask[2] == mask[3])
+ mask[1] = mask[2];
+ for (size_t l = 2; l < n-l0-3; ++l)
{
- // Update wsindex to index in input workspace
- wsindex = static_cast<size_t>(i);
+ if (mask[l-1] == mask[l+1] && (mask[l-1] == mask[l-2] || mask[l+1] == mask[l+2]))
+ {
+ mask[l] = mask[l+1];
+ }
}
- else
+ if (mask[n-l0-2] == mask[n-l0-3] && mask[n-l0-3] == mask[n-l0-4])
+ mask[n-l0-1] = mask[n-l0-2];
+ if (mask[n-l0-1] == mask[n-l0-3] && mask[n-l0-3] == mask[n-l0-4])
+ mask[n-l0-2] = mask[n-l0-1];
+
+ // mask regions not connected to largest region
+ // for loop can start > 1 for multiple peaks
+ vector<cont_peak> peaks;
+ if (mask[0] == 1)
{
- // Use the wsindex as the input
- wsindex = static_cast<size_t>(inpwsindex[i]);
- if (wsindex >= inpWS->getNumberHistograms())
+ peaks.push_back(cont_peak());
+ peaks[peaks.size()-1].start = l0;
+ }
+ for (size_t l = 1; l < n-l0; ++l)
+ {
+ if (mask[l] != mask[l-1] && mask[l] == 1)
+ {
+ peaks.push_back(cont_peak());
+ peaks[peaks.size()-1].start = l+l0;
+ }
+ else if (peaks.size() > 0)
{
- stringstream errmsg;
- errmsg << "Input workspace index " << inpwsindex[i] << " is out of input workspace range = "
- << inpWS->getNumberHistograms() << endl;
+ size_t ipeak = peaks.size()-1;
+ if (mask[l] != mask[l-1] && mask[l] == 0)
+ {
+ peaks[ipeak].stop = l+l0;
+ }
+ if (inpY[l+l0] > peaks[ipeak].maxY) peaks[ipeak].maxY = inpY[l+l0];
}
}
-
- // Find background
- const MantidVec& inpX = inpWS->readX(wsindex);
- const MantidVec& inpY = inpWS->readY(wsindex);
-
- double Ymean, Yvariance, Ysigma;
- MantidVec maskedY;
- MantidVec::const_iterator in = std::min_element(inpY.begin(), inpY.end());
- double bkg0 = inpY[in - inpY.begin()];
- for (size_t l = l0; l < n; ++l)
+ size_t min_peak, max_peak;
+ double a0,a1,a2;
+ if(peaks.size()> 0)
{
- maskedY.push_back(inpY[l]-bkg0);
+ if(peaks[peaks.size()-1].stop == 0) peaks[peaks.size()-1].stop = n-1;
+ std::sort(peaks.begin(), peaks.end(), by_len());
+
+ // save endpoints
+ min_peak = peaks[0].start;
+ // extra point for histogram input
+ max_peak = peaks[0].stop + sizex - sizey;
+ estimateBackground(inpX, inpY, l0, n,
+ peaks[0].start, peaks[0].stop, a0, a1, a2);
}
- MantidVec mask(n-l0,0.0);
- double xn = static_cast<double>(n-l0);
- do
+ else
{
- Statistics stats = getStatistics(maskedY);
- Ymean = stats.mean;
- Yvariance = stats.standard_deviation * stats.standard_deviation;
- Ysigma = std::sqrt((moment4(maskedY,n-l0,Ymean)-(xn-3.0)/(xn-1.0) * Yvariance)/xn);
- MantidVec::const_iterator it = std::max_element(maskedY.begin(), maskedY.end());
- const size_t pos = it - maskedY.begin();
- maskedY[pos] = 0;
- mask[pos] = 1.0;
+ // assume background is 12 first and last points
+ min_peak = l0+12;
+ max_peak = n-13;
+ if (min_peak > sizey)min_peak = sizey-1;
+ a0 = 0.0;
+ a1 = 0.0;
+ a2 = 0.0;
}
- while (std::abs(Ymean-Yvariance) > k * Ysigma);
- if(n-l0 > 5)
- {
- // remove single outliers
- if (mask[1] == mask[2] && mask[2] == mask[3])
- mask[0] = mask[1];
- if (mask[0] == mask[2] && mask[2] == mask[3])
- mask[1] = mask[2];
- for (size_t l = 2; l < n-l0-3; ++l)
- {
- if (mask[l-1] == mask[l+1] && (mask[l-1] == mask[l-2] || mask[l+1] == mask[l+2]))
- {
- mask[l] = mask[l+1];
- }
- }
- if (mask[n-l0-2] == mask[n-l0-3] && mask[n-l0-3] == mask[n-l0-4])
- mask[n-l0-1] = mask[n-l0-2];
- if (mask[n-l0-1] == mask[n-l0-3] && mask[n-l0-3] == mask[n-l0-4])
- mask[n-l0-2] = mask[n-l0-1];
-
- // mask regions not connected to largest region
- // for loop can start > 1 for multiple peaks
- vector<cont_peak> peaks;
- if (mask[0] == 1)
- {
- peaks.push_back(cont_peak());
- peaks[peaks.size()-1].start = l0;
- }
- for (size_t l = 1; l < n-l0; ++l)
- {
- if (mask[l] != mask[l-1] && mask[l] == 1)
- {
- peaks.push_back(cont_peak());
- peaks[peaks.size()-1].start = l+l0;
- }
- else if (peaks.size() > 0)
- {
- size_t ipeak = peaks.size()-1;
- if (mask[l] != mask[l-1] && mask[l] == 0)
- {
- peaks[ipeak].stop = l+l0;
- }
- if (inpY[l+l0] > peaks[ipeak].maxY) peaks[ipeak].maxY = inpY[l+l0];
- }
- }
- size_t min_peak, max_peak;
- double a0,a1,a2;
- if(peaks.size()> 0)
- {
- if(peaks[peaks.size()-1].stop == 0) peaks[peaks.size()-1].stop = n-1;
- std::sort(peaks.begin(), peaks.end(), by_len());
-
- // save endpoints
- min_peak = peaks[0].start;
- // extra point for histogram input
- max_peak = peaks[0].stop + sizex - sizey;
- estimateBackground(inpX, inpY, l0, n,
- peaks[0].start, peaks[0].stop, a0, a1, a2);
- }
- else
- {
- // assume background is 12 first and last points
- min_peak = l0+12;
- max_peak = n-13;
- if (min_peak > sizey)min_peak = sizey-1;
- a0 = 0.0;
- a1 = 0.0;
- a2 = 0.0;
- }
-
- // Add a new row
- API::TableRow t = m_outPeakTableWS->getRow(i);
- t << static_cast<int>(wsindex) << static_cast<int>(min_peak) << static_cast<int>(max_peak) << a0 << a1 <<a2;
- }
+ // Add a new row
+ API::TableRow t = m_outPeakTableWS->getRow(0);
+ t << static_cast<int>(inpwsindex) << static_cast<int>(min_peak) << static_cast<int>(max_peak) << a0 << a1 <<a2;
+ }
- prog.report();
- PARALLEL_END_INTERUPT_REGION
- } // ENDFOR
- PARALLEL_CHECK_INTERUPT_REGION
+ prog.report();
// 4. Set the output
setProperty("OutputWorkspace", m_outPeakTableWS);
return;
}
+
//----------------------------------------------------------------------------------------------
/** Estimate background
* @param X :: vec for X
diff --git a/Code/Mantid/Framework/Algorithms/test/FindPeakBackgroundTest.h b/Code/Mantid/Framework/Algorithms/test/FindPeakBackgroundTest.h
index 670838a0dcfe09d18502551b296636ed39d1b19b..264e07aeb066343a741786e633e1b109773e7328 100644
--- a/Code/Mantid/Framework/Algorithms/test/FindPeakBackgroundTest.h
+++ b/Code/Mantid/Framework/Algorithms/test/FindPeakBackgroundTest.h
@@ -6,6 +6,7 @@
#include "MantidAlgorithms/FindPeakBackground.h"
#include "MantidAPI/WorkspaceFactory.h"
#include "MantidAPI/MatrixWorkspace.h"
+#include "MantidAPI/ITableWorkspace.h"
#include "MantidDataObjects/Workspace2D.h"
#include <cmath>
@@ -30,78 +31,222 @@ public:
delete suite;
}
- void test_Calculation() {
+ void test_Calculation()
+ {
+ // 1. Generate input workspace
+ MatrixWorkspace_sptr inWS = generateTestWorkspace();
- // 1. Generate input workspace
- MatrixWorkspace_sptr inWS = generateTestWorkspace();
+ // 2. Create
+ Algorithms::FindPeakBackground alg;
- // 2. Create
- Algorithms::FindPeakBackground alg;
+ alg.initialize();
+ TS_ASSERT(alg.isInitialized());
- alg.initialize();
- TS_ASSERT(alg.isInitialized());
+ alg.setProperty("InputWorkspace", inWS);
+ alg.setProperty("OutputWorkspace", "Signal");
+ alg.setProperty("WorkspaceIndex", 0);
- alg.setProperty("InputWorkspace", inWS);
- alg.setProperty("OutputWorkspace", "Signal");
- alg.setProperty("WorkspaceIndices", "0");
+ alg.execute();
+ TS_ASSERT(alg.isExecuted());
- alg.execute();
- TS_ASSERT(alg.isExecuted());
+ Mantid::API::ITableWorkspace_sptr peaklist = boost::dynamic_pointer_cast<Mantid::API::ITableWorkspace>
+ (Mantid::API::AnalysisDataService::Instance().retrieve("Signal"));
- Mantid::API::ITableWorkspace_sptr peaklist = boost::dynamic_pointer_cast<Mantid::API::ITableWorkspace>
- (Mantid::API::AnalysisDataService::Instance().retrieve("Signal"));
+ TS_ASSERT( peaklist );
+ TS_ASSERT_EQUALS( peaklist->rowCount() , 1 );
+ TS_ASSERT_DELTA( peaklist->Int(0,1), 4, 0.01 );
+ TS_ASSERT_DELTA( peaklist->Int(0,2), 19, 0.01 );
+ TS_ASSERT_DELTA( peaklist->Double(0,3), 1.2, 0.01 );
+ TS_ASSERT_DELTA( peaklist->Double(0,4), 0.04, 0.01 );
+ TS_ASSERT_DELTA( peaklist->Double(0,5), 0.0, 0.01 );
- TS_ASSERT( peaklist );
- TS_ASSERT_EQUALS( peaklist->rowCount() , 1 );
- TS_ASSERT_DELTA( peaklist->Int(0,1), 4, 0.01 );
- TS_ASSERT_DELTA( peaklist->Int(0,2), 19, 0.01 );
- TS_ASSERT_DELTA( peaklist->Double(0,3), 1.2, 0.01 );
- TS_ASSERT_DELTA( peaklist->Double(0,4), 0.04, 0.01 );
- TS_ASSERT_DELTA( peaklist->Double(0,5), 0.0, 0.01 );
+ AnalysisDataService::Instance().remove("Signal");
- return;
+ return;
}
/** Generate a workspace for test
*/
- MatrixWorkspace_sptr generateTestWorkspace() {
- vector<double> data;
- data.push_back(1);
- data.push_back(2);
- data.push_back(1);
- data.push_back(1);
- data.push_back(9);
- data.push_back(11);
- data.push_back(13);
- data.push_back(20);
- data.push_back(24);
- data.push_back(32);
- data.push_back(28);
- data.push_back(48);
- data.push_back(42);
- data.push_back(77);
- data.push_back(67);
- data.push_back(33);
- data.push_back(27);
- data.push_back(20);
- data.push_back(9);
- data.push_back(2);
-
- MatrixWorkspace_sptr ws = boost::dynamic_pointer_cast<MatrixWorkspace>(
- WorkspaceFactory::Instance().create("Workspace2D", 1,
- data.size(), data.size()));
-
- MantidVec& vecX = ws->dataX(0);
- MantidVec& vecY = ws->dataY(0);
- MantidVec& vecE = ws->dataE(0);
-
- for (size_t i = 0; i < data.size(); ++i) {
- vecX[i] = static_cast<double> (i);
- vecY[i] = data[i];
- vecE[i] = sqrt(data[i]);
- }
-
- return ws;
+ MatrixWorkspace_sptr generateTestWorkspace()
+ {
+ vector<double> data;
+ data.push_back(1);
+ data.push_back(2);
+ data.push_back(1);
+ data.push_back(1);
+ data.push_back(9);
+ data.push_back(11);
+ data.push_back(13);
+ data.push_back(20);
+ data.push_back(24);
+ data.push_back(32);
+ data.push_back(28);
+ data.push_back(48);
+ data.push_back(42);
+ data.push_back(77);
+ data.push_back(67);
+ data.push_back(33);
+ data.push_back(27);
+ data.push_back(20);
+ data.push_back(9);
+ data.push_back(2);
+
+ MatrixWorkspace_sptr ws = boost::dynamic_pointer_cast<MatrixWorkspace>(
+ WorkspaceFactory::Instance().create("Workspace2D", 1,
+ data.size(), data.size()));
+
+ MantidVec& vecX = ws->dataX(0);
+ MantidVec& vecY = ws->dataY(0);
+ MantidVec& vecE = ws->dataE(0);
+
+ for (size_t i = 0; i < data.size(); ++i) {
+ vecX[i] = static_cast<double> (i);
+ vecY[i] = data[i];
+ vecE[i] = sqrt(data[i]);
+ }
+
+ return ws;
+ }
+
+ //--------------------------------------------------------------------------------------------
+ /** Test on a spectrum without peak
+ */
+ void test_FindBackgroundOnFlat()
+ {
+ // Add workspace
+ MatrixWorkspace_sptr testws = generate2SpectraTestWorkspace();
+ AnalysisDataService::Instance().addOrReplace("Test2Workspace", testws);
+
+ // Set up algorithm
+ Algorithms::FindPeakBackground alg;
+
+ alg.initialize();
+ TS_ASSERT(alg.isInitialized());
+
+ alg.setProperty("InputWorkspace", "Test2Workspace");
+ alg.setProperty("OutputWorkspace", "Signal3");
+ alg.setProperty("WorkspaceIndex", 0);
+
+ // Execute
+ alg.execute();
+ TS_ASSERT(alg.isExecuted());
+
+ // Check result
+ ITableWorkspace_sptr outws = boost::dynamic_pointer_cast<ITableWorkspace>(
+ AnalysisDataService::Instance().retrieve("Signal3"));
+ TS_ASSERT(outws);
+ if (!outws)
+ return;
+
+ TS_ASSERT_EQUALS(outws->rowCount(), 1);
+ if (outws->rowCount() < 1)
+ return;
+
+ int ipeakmin = outws->Int(0, 1);
+ int ipeakmax = outws->Int(0, 2);
+ TS_ASSERT(ipeakmin >= ipeakmax);
+
+ return;
+ }
+
+ //--------------------------------------------------------------------------------------------
+ /** Test on a spectrum without peak
+ */
+ void test_FindBackgroundOnSpec1()
+ {
+ // Add workspace
+ MatrixWorkspace_sptr testws = generate2SpectraTestWorkspace();
+ AnalysisDataService::Instance().addOrReplace("Test2Workspace", testws);
+
+ // Set up algorithm
+ Algorithms::FindPeakBackground alg;
+
+ alg.initialize();
+ TS_ASSERT(alg.isInitialized());
+
+ alg.setProperty("InputWorkspace", "Test2Workspace");
+ alg.setProperty("OutputWorkspace", "Signal2");
+ alg.setProperty("WorkspaceIndex", 1);
+
+ // Execute
+ alg.execute();
+ TS_ASSERT(alg.isExecuted());
+
+ // Check result
+ ITableWorkspace_sptr outws = boost::dynamic_pointer_cast<ITableWorkspace>(
+ AnalysisDataService::Instance().retrieve("Signal2"));
+ TS_ASSERT(outws);
+ if (!outws)
+ return;
+
+ TS_ASSERT_EQUALS(outws->rowCount(), 1);
+ if (outws->rowCount() < 1)
+ return;
+
+ TS_ASSERT_EQUALS( outws->rowCount() , 1 );
+ TS_ASSERT_DELTA( outws->Int(0,1), 4, 0.01 );
+ TS_ASSERT_DELTA( outws->Int(0,2), 19, 0.01 );
+ TS_ASSERT_DELTA( outws->Double(0,3), 1.2, 0.01 );
+ TS_ASSERT_DELTA( outws->Double(0,4), 0.04, 0.01 );
+ TS_ASSERT_DELTA( outws->Double(0,5), 0.0, 0.01 );
+
+ return;
+ }
+
+ //--------------------------------------------------------------------------------------------
+ /** Generate a workspace with 2 spectra for test
+ */
+ MatrixWorkspace_sptr generate2SpectraTestWorkspace()
+ {
+ vector<double> data;
+ data.push_back(1);
+ data.push_back(2);
+ data.push_back(1);
+ data.push_back(1);
+ data.push_back(9);
+ data.push_back(11);
+ data.push_back(13);
+ data.push_back(20);
+ data.push_back(24);
+ data.push_back(32);
+ data.push_back(28);
+ data.push_back(48);
+ data.push_back(42);
+ data.push_back(77);
+ data.push_back(67);
+ data.push_back(33);
+ data.push_back(27);
+ data.push_back(20);
+ data.push_back(9);
+ data.push_back(2);
+
+ MatrixWorkspace_sptr ws = boost::dynamic_pointer_cast<MatrixWorkspace>(
+ WorkspaceFactory::Instance().create("Workspace2D", 2,
+ data.size(), data.size()));
+
+ // Workspace index = 0
+ MantidVec& vecX = ws->dataX(0);
+ MantidVec& vecY = ws->dataY(0);
+ MantidVec& vecE = ws->dataE(0);
+ for (size_t i = 0; i < data.size(); ++i)
+ {
+ vecX[i] = static_cast<double> (i);
+ vecY[i] = 0.0;
+ vecE[i] = 1.0;
+ }
+
+ // Workspace index = 1
+ MantidVec& vecX1 = ws->dataX(1);
+ MantidVec& vecY1 = ws->dataY(1);
+ MantidVec& vecE1 = ws->dataE(1);
+ for (size_t i = 0; i < data.size(); ++i)
+ {
+ vecX1[i] = static_cast<double> (i);
+ vecY1[i] = data[i];
+ vecE1[i] = sqrt(data[i]);
+ }
+
+ return ws;
}
};