Skip to content
Snippets Groups Projects
Commit 441a4bb8 authored by Zhou, Wenduo's avatar Zhou, Wenduo
Browse files

Refactored FindPeakBackground. Refs #9020. 

1. Refactored and thus fixed some issues;
2. Added  more unit tests;
3. Change input property. Algorithm can work with 1 spectrum only.
parent 12827dca
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Code/Mantid/Framework/Algorithms/src/FindPeakBackground.cpp
+ 107
− 137
View file @ 441a4bb8
...@@ -76,9 +76,8 @@ namespace Algorithms ...@@ -76,9 +76,8 @@ namespace Algorithms
auto inwsprop = new WorkspaceProperty<MatrixWorkspace>("InputWorkspace", "Anonymous", Direction::Input); auto inwsprop = new WorkspaceProperty<MatrixWorkspace>("InputWorkspace", "Anonymous", Direction::Input);
declareProperty(inwsprop, "Name of input MatrixWorkspace that contains peaks."); declareProperty(inwsprop, "Name of input MatrixWorkspace that contains peaks.");
declareProperty(new ArrayProperty<int>("WorkspaceIndices"), "Optional: enter a comma-separated list of the " declareProperty("WorkspaceIndex", EMPTY_INT(), "workspace indices to have peak and background separated. "
"workspace indices to have peak and background separated. " "No default is taken. ");
"Default is to calculate for all spectra.");
declareProperty("SigmaConstant", 1.0, "Multiplier of standard deviations of the variance for convergence of " declareProperty("SigmaConstant", 1.0, "Multiplier of standard deviations of the variance for convergence of "
"peak elimination. Default is 1.0. "); "peak elimination. Default is 1.0. ");
...@@ -109,32 +108,27 @@ namespace Algorithms ...@@ -109,32 +108,27 @@ namespace Algorithms
{ {
// 1. Get input and validate // 1. Get input and validate
MatrixWorkspace_const_sptr inpWS = getProperty("InputWorkspace"); MatrixWorkspace_const_sptr inpWS = getProperty("InputWorkspace");
std::vector<int> inpwsindex = getProperty("WorkspaceIndices"); int inpwsindex = getProperty("WorkspaceIndex");
std::vector<double> m_vecFitWindows = getProperty("FitWindow"); std::vector<double> m_vecFitWindows = getProperty("FitWindow");
m_backgroundType = getPropertyValue("BackgroundType"); m_backgroundType = getPropertyValue("BackgroundType");
double k = getProperty("SigmaConstant"); double k = getProperty("SigmaConstant");
bool separateall = false; if (inpwsindex < 0 || inpwsindex >= static_cast<int>(inpWS->getNumberHistograms()))
if (inpwsindex.size() == 0)
{ {
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 // 2. Generate output
size_t numspec; const MantidVec& inpX = inpWS->readX(inpwsindex);
if (separateall) size_t sizex = inpWS->readX(inpwsindex).size();
{ size_t sizey = inpWS->readY(inpwsindex).size();
numspec = inpWS->getNumberHistograms();
}
else
{
numspec = inpwsindex.size();
}
size_t sizex = inpWS->readX(0).size();
size_t sizey = inpWS->readY(0).size();
size_t n = sizey; size_t n = sizey;
size_t l0 = 0; size_t l0 = 0;
const MantidVec& inpX = inpWS->readX(0);
if (m_vecFitWindows.size() > 1) if (m_vecFitWindows.size() > 1)
{ {
...@@ -152,145 +146,121 @@ namespace Algorithms ...@@ -152,145 +146,121 @@ namespace Algorithms
m_outPeakTableWS->addColumn("double", "bkg0"); m_outPeakTableWS->addColumn("double", "bkg0");
m_outPeakTableWS->addColumn("double", "bkg1"); m_outPeakTableWS->addColumn("double", "bkg1");
m_outPeakTableWS->addColumn("double", "bkg2"); m_outPeakTableWS->addColumn("double", "bkg2");
for( size_t i = 0; i < numspec; ++i )
m_outPeakTableWS->appendRow(); m_outPeakTableWS->appendRow();
// 3. Get Y values // 3. Get Y values
Progress prog(this, 0, 1.0, numspec); Progress prog(this, 0, 1.0, 1);
PARALLEL_FOR2(inpWS, m_outPeakTableWS)
for (int i = 0; i < static_cast<int>(numspec); ++i) // 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 // remove single outliers
// a) figure out wsindex if (mask[1] == mask[2] && mask[2] == mask[3])
size_t wsindex; mask[0] = mask[1];
if (separateall) 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 if (mask[l-1] == mask[l+1] && (mask[l-1] == mask[l-2] || mask[l+1] == mask[l+2]))
wsindex = static_cast<size_t>(i); {
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 peaks.push_back(cont_peak());
wsindex = static_cast<size_t>(inpwsindex[i]); peaks[peaks.size()-1].start = l0;
if (wsindex >= inpWS->getNumberHistograms()) }
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; size_t ipeak = peaks.size()-1;
errmsg << "Input workspace index " << inpwsindex[i] << " is out of input workspace range = " if (mask[l] != mask[l-1] && mask[l] == 0)
<< inpWS->getNumberHistograms() << endl; {
peaks[ipeak].stop = l+l0;
}
if (inpY[l+l0] > peaks[ipeak].maxY) peaks[ipeak].maxY = inpY[l+l0];
} }
} }
size_t min_peak, max_peak;
// Find background double a0,a1,a2;
const MantidVec& inpX = inpWS->readX(wsindex); if(peaks.size()> 0)
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)
{ {
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); else
double xn = static_cast<double>(n-l0);
do
{ {
Statistics stats = getStatistics(maskedY); // assume background is 12 first and last points
Ymean = stats.mean; min_peak = l0+12;
Yvariance = stats.standard_deviation * stats.standard_deviation; max_peak = n-13;
Ysigma = std::sqrt((moment4(maskedY,n-l0,Ymean)-(xn-3.0)/(xn-1.0) * Yvariance)/xn); if (min_peak > sizey)min_peak = sizey-1;
MantidVec::const_iterator it = std::max_element(maskedY.begin(), maskedY.end()); a0 = 0.0;
const size_t pos = it - maskedY.begin(); a1 = 0.0;
maskedY[pos] = 0; a2 = 0.0;
mask[pos] = 1.0;
} }
while (std::abs(Ymean-Yvariance) > k * Ysigma);
if(n-l0 > 5) // Add a new row
{ API::TableRow t = m_outPeakTableWS->getRow(0);
// remove single outliers t << static_cast<int>(inpwsindex) << static_cast<int>(min_peak) << static_cast<int>(max_peak) << a0 << a1 <<a2;
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;
}
prog.report(); prog.report();
PARALLEL_END_INTERUPT_REGION
} // ENDFOR
PARALLEL_CHECK_INTERUPT_REGION
// 4. Set the output // 4. Set the output
setProperty("OutputWorkspace", m_outPeakTableWS); setProperty("OutputWorkspace", m_outPeakTableWS);
return; return;
} }
//---------------------------------------------------------------------------------------------- //----------------------------------------------------------------------------------------------
/** Estimate background /** Estimate background
* @param X :: vec for X * @param X :: vec for X
......
Code/Mantid/Framework/Algorithms/test/FindPeakBackgroundTest.h
+ 205
− 60
View file @ 441a4bb8
...@@ -6,6 +6,7 @@ ...@@ -6,6 +6,7 @@
#include "MantidAlgorithms/FindPeakBackground.h" #include "MantidAlgorithms/FindPeakBackground.h"
#include "MantidAPI/WorkspaceFactory.h" #include "MantidAPI/WorkspaceFactory.h"
#include "MantidAPI/MatrixWorkspace.h" #include "MantidAPI/MatrixWorkspace.h"
#include "MantidAPI/ITableWorkspace.h"
#include "MantidDataObjects/Workspace2D.h" #include "MantidDataObjects/Workspace2D.h"
#include <cmath> #include <cmath>
...@@ -30,78 +31,222 @@ public: ...@@ -30,78 +31,222 @@ public:
delete suite; delete suite;
} }
void test_Calculation() { void test_Calculation()
{
// 1. Generate input workspace
MatrixWorkspace_sptr inWS = generateTestWorkspace();
// 1. Generate input workspace // 2. Create
MatrixWorkspace_sptr inWS = generateTestWorkspace(); Algorithms::FindPeakBackground alg;
// 2. Create alg.initialize();
Algorithms::FindPeakBackground alg; TS_ASSERT(alg.isInitialized());
alg.initialize(); alg.setProperty("InputWorkspace", inWS);
TS_ASSERT(alg.isInitialized()); alg.setProperty("OutputWorkspace", "Signal");
alg.setProperty("WorkspaceIndex", 0);
alg.setProperty("InputWorkspace", inWS); alg.execute();
alg.setProperty("OutputWorkspace", "Signal"); TS_ASSERT(alg.isExecuted());
alg.setProperty("WorkspaceIndices", "0");
alg.execute(); Mantid::API::ITableWorkspace_sptr peaklist = boost::dynamic_pointer_cast<Mantid::API::ITableWorkspace>
TS_ASSERT(alg.isExecuted()); (Mantid::API::AnalysisDataService::Instance().retrieve("Signal"));
Mantid::API::ITableWorkspace_sptr peaklist = boost::dynamic_pointer_cast<Mantid::API::ITableWorkspace> TS_ASSERT( peaklist );
(Mantid::API::AnalysisDataService::Instance().retrieve("Signal")); 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 ); AnalysisDataService::Instance().remove("Signal");
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 );
return; return;
} }
/** Generate a workspace for test /** Generate a workspace for test
*/ */
MatrixWorkspace_sptr generateTestWorkspace() { MatrixWorkspace_sptr generateTestWorkspace()
vector<double> data; {
data.push_back(1); vector<double> data;
data.push_back(2); data.push_back(1);
data.push_back(1); data.push_back(2);
data.push_back(1); data.push_back(1);
data.push_back(9); data.push_back(1);
data.push_back(11); data.push_back(9);
data.push_back(13); data.push_back(11);
data.push_back(20); data.push_back(13);
data.push_back(24); data.push_back(20);
data.push_back(32); data.push_back(24);
data.push_back(28); data.push_back(32);
data.push_back(48); data.push_back(28);
data.push_back(42); data.push_back(48);
data.push_back(77); data.push_back(42);
data.push_back(67); data.push_back(77);
data.push_back(33); data.push_back(67);
data.push_back(27); data.push_back(33);
data.push_back(20); data.push_back(27);
data.push_back(9); data.push_back(20);
data.push_back(2); data.push_back(9);
data.push_back(2);
MatrixWorkspace_sptr ws = boost::dynamic_pointer_cast<MatrixWorkspace>(
WorkspaceFactory::Instance().create("Workspace2D", 1, MatrixWorkspace_sptr ws = boost::dynamic_pointer_cast<MatrixWorkspace>(
data.size(), data.size())); WorkspaceFactory::Instance().create("Workspace2D", 1,
data.size(), data.size()));
MantidVec& vecX = ws->dataX(0);
MantidVec& vecY = ws->dataY(0); MantidVec& vecX = ws->dataX(0);
MantidVec& vecE = ws->dataE(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); for (size_t i = 0; i < data.size(); ++i) {
vecY[i] = data[i]; vecX[i] = static_cast<double> (i);
vecE[i] = sqrt(data[i]); vecY[i] = data[i];
} vecE[i] = sqrt(data[i]);
}
return ws;
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;
} }
}; };
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment