Skip to content

Commit

Permalink
Made the new feature work. Refs #8472.
Browse files Browse the repository at this point in the history
  • Loading branch information
@wdzhou
wdzhou committed Nov 21, 2013
1 parent 7a38ac1 commit 172f805
Show file tree
Hide file tree
Showing 2 changed files with 158 additions and 14 deletions.
8 changes: 7 additions & 1 deletion Code/Mantid/Framework/CurveFitting/inc/MantidCurveFitting/ProcessBackground.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -101,7 +101,7 @@ class DLLExport ProcessBackground : public API::Algorithm
DataObjects::Workspace2D_sptr autoBackgroundSelection(DataObjects::Workspace2D_sptr bkgdWS);

/// Create a background function from input properties
BackgroundFunction_sptr createBackgroundFunction();
BackgroundFunction_sptr createBackgroundFunction(const std::string backgroundtype);

/// Filter non-background data points out and create a background workspace
DataObjects::Workspace2D_sptr filterForBackground(BackgroundFunction_sptr bkgdfunction);
Expand All @@ -114,6 +114,10 @@ class DLLExport ProcessBackground : public API::Algorithm
double m_lowerBound;
double m_upperBound;

std::string m_bkgdType;

bool m_doFitBackground;

// double mTolerance;

/// Number of FWHM of range of peak to be removed.
Expand All @@ -127,6 +131,8 @@ class DLLExport ProcessBackground : public API::Algorithm

/// Add a certain region from a reference workspace
void addRegion();

void fitBackgroundFunction(std::string bkgdfunctiontype);

};

Expand Down
164 changes: 151 additions & 13 deletions Code/Mantid/Framework/CurveFitting/src/ProcessBackground.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -179,20 +179,23 @@ DECLARE_ALGORITHM(ProcessBackground)
new VisibleWhenProperty("Options", IS_EQUAL_TO, "SelectBackgroundPoints"));

// Optional output workspace
// TODO - Implement more upon this option!
declareProperty(new WorkspaceProperty<Workspace2D>("OutputBackgroundParameterTable", "", Direction::Output,
PropertyMode::Optional),
"Output .... ... .");
declareProperty(new WorkspaceProperty<TableWorkspace>("OutputBackgroundParameterTable", "", Direction::Output,
PropertyMode::Optional),
"Output parameter table workspace containing the background fitting result. ");
setPropertySettings("OutputBackgroundParameterTable",
new VisibleWhenProperty("Options", IS_EQUAL_TO, "SelectBackgroundPoints"));

// TODO - Develop from this new option
// Output background type.
declareProperty("OutputBackgroundType", "Polynomial", bkgdvalidator,
"Type of background to fit with selected ... .");
"Type of background to fit with selected background points.");
setPropertySettings("OutputBackgroundType",
new VisibleWhenProperty("Options", IS_EQUAL_TO, "SelectBackgroundPoints"));

// Output background type.
declareProperty("OutputBackgroundOrder", 6,
"Order of background to fit with selected background points.");
setPropertySettings("OutputBackgroundOrder",
new VisibleWhenProperty("Options", IS_EQUAL_TO, "SelectBackgroundPoints"));

// Peak table workspac for "RemovePeaks"
declareProperty(new WorkspaceProperty<TableWorkspace>("BraggPeakTableWorkspace", "", Direction::Input,
Expand Down Expand Up @@ -222,6 +225,8 @@ DECLARE_ALGORITHM(ProcessBackground)
throw std::invalid_argument("Input Workspace cannot be obtained.");
}

m_bkgdType = getPropertyValue("BackgroundType");

int intemp = getProperty("WorkspaceIndex");
if (intemp < 0)
throw std::invalid_argument("WorkspaceIndex is not allowed to be less than 0. ");
Expand Down Expand Up @@ -252,6 +257,17 @@ DECLARE_ALGORITHM(ProcessBackground)
}
else if (option.compare("SelectBackgroundPoints") == 0)
{
string outbkgdparwsname = getPropertyValue("OutputBackgroundParameterTable");
if (outbkgdparwsname.size() > 0)
{
// Will fit the selected background
m_doFitBackground = true;
}
else
{
m_doFitBackground = false;
}

string smode = getProperty("SelectionMode");
if (smode == "FitGivenDataPoints")
{
Expand All @@ -265,6 +281,14 @@ DECLARE_ALGORITHM(ProcessBackground)
{
throw runtime_error("N/A is not supported.");
}

if (m_doFitBackground)
{
// Fit the selected background
string bkgdfunctype = getPropertyValue("OutputBackgroundType");
fitBackgroundFunction(bkgdfunctype);
}

}
else
{
Expand Down Expand Up @@ -648,7 +672,7 @@ DECLARE_ALGORITHM(ProcessBackground)
void ProcessBackground::execSelectBkgdPoints2()
{
// Process properties
BackgroundFunction_sptr bkgdfunc = createBackgroundFunction();
BackgroundFunction_sptr bkgdfunc = createBackgroundFunction(m_bkgdType);
TableWorkspace_sptr bkgdtablews = getProperty("BackgroundTableWorkspace");

// Set up background function from table
Expand Down Expand Up @@ -686,7 +710,7 @@ DECLARE_ALGORITHM(ProcessBackground)
DataObjects::Workspace2D_sptr ProcessBackground::autoBackgroundSelection(Workspace2D_sptr bkgdWS)
{
// Get background type and create bakground function
BackgroundFunction_sptr bkgdfunction = createBackgroundFunction();
BackgroundFunction_sptr bkgdfunction = createBackgroundFunction(m_bkgdType);

int bkgdorder = getProperty("BackgroundOrder");
bkgdfunction->setAttributeValue("n", bkgdorder);
Expand Down Expand Up @@ -759,10 +783,8 @@ DECLARE_ALGORITHM(ProcessBackground)
//----------------------------------------------------------------------------------------------
/** Create a background function from input properties
*/
BackgroundFunction_sptr ProcessBackground::createBackgroundFunction()
BackgroundFunction_sptr ProcessBackground::createBackgroundFunction(const string backgroundtype)
{
std::string backgroundtype = getProperty("BackgroundType");

CurveFitting::BackgroundFunction_sptr bkgdfunction;

if (backgroundtype.compare("Polynomial") == 0)
Expand Down Expand Up @@ -858,12 +880,19 @@ DECLARE_ALGORITHM(ProcessBackground)
<< " total data points. " << "\n";

// Build new workspace
size_t nspec;
if (m_doFitBackground)
nspec = 3;
else
nspec = 1;

Workspace2D_sptr outws = boost::dynamic_pointer_cast<DataObjects::Workspace2D>
(API::WorkspaceFactory::Instance().create("Workspace2D", 1, vecx.size(), vecy.size()));
(API::WorkspaceFactory::Instance().create("Workspace2D", nspec, vecx.size(), vecy.size()));

for (size_t i = 0; i < vecx.size(); ++i)
{
outws->dataX(0)[i] = vecx[i];
for(size_t j = 0; j < nspec; ++j)
outws->dataX(j)[i] = vecx[i];
outws->dataY(0)[i] = vecy[i];
outws->dataE(0)[i] = vece[i];
}
Expand Down Expand Up @@ -1060,6 +1089,115 @@ DECLARE_ALGORITHM(ProcessBackground)
return count;
}

//----------------------------------------------------------------------------------------------
/** Fit background function
*/
void ProcessBackground::fitBackgroundFunction(std::string bkgdfunctiontype)
{
// Get background type and create bakground function
BackgroundFunction_sptr bkgdfunction = createBackgroundFunction(bkgdfunctiontype);

int bkgdorder = getProperty("OutputBackgroundOrder");
bkgdfunction->setAttributeValue("n", bkgdorder);

if (bkgdfunctiontype == "Chebyshev")
{
double xmin = m_outputWS->readX(0).front();
double xmax = m_outputWS->readX(0).back();
g_log.information() << "Chebyshev Fit range: " << xmin << ", " << xmax << "\n";
bkgdfunction->setAttributeValue("StartX", xmin);
bkgdfunction->setAttributeValue("EndX", xmax);
}

g_log.information() << "Fit selected background " << bkgdfunctiontype
<< " to data workspace with " << m_outputWS->getNumberHistograms() << " spectra."
<< "\n";

// Fit input (a few) background pionts to get initial guess
API::IAlgorithm_sptr fit;
try
{
fit = this->createChildAlgorithm("Fit", 0.9, 1.0, true);
}
catch (Exception::NotFoundError &)
{
g_log.error() << "Requires CurveFitting library." << std::endl;
throw;
}

g_log.information() << "Fitting background function: " << bkgdfunction->asString() << "\n";

double startx = m_lowerBound;
double endx = m_upperBound;
fit->setProperty("Function", boost::dynamic_pointer_cast<API::IFunction>(bkgdfunction));
fit->setProperty("InputWorkspace", m_outputWS);
fit->setProperty("WorkspaceIndex", 0);
fit->setProperty("MaxIterations", 500);
fit->setProperty("StartX", startx);
fit->setProperty("EndX", endx);
fit->setProperty("Minimizer", "Levenberg-MarquardtMD");
fit->setProperty("CostFunction", "Least squares");

fit->executeAsChildAlg();

// Get fit status and chi^2
std::string fitStatus = fit->getProperty("OutputStatus");
bool allowedfailure = (fitStatus.find("cannot") < fitStatus.size()) &&
(fitStatus.find("tolerance") < fitStatus.size());
if (fitStatus.compare("success") != 0 && !allowedfailure)
{
g_log.error() << "ProcessBackground: Fit Status = " << fitStatus
<< ". Not to update fit result" << std::endl;
throw std::runtime_error("Bad Fit");
}

const double chi2 = fit->getProperty("OutputChi2overDoF");
g_log.information() << "Fit background: Fit Status = " << fitStatus << ", chi2 = "
<< chi2 << "\n";

// Get out the parameter names
API::IFunction_sptr funcout = fit->getProperty("Function");
TableWorkspace_sptr outbkgdparws = boost::make_shared<TableWorkspace>();
outbkgdparws->addColumn("str", "Name");
outbkgdparws->addColumn("double", "Value");

TableRow typerow = outbkgdparws->appendRow();
typerow << bkgdfunctiontype << 0.;

vector<string> parnames = funcout->getParameterNames();
size_t nparam = funcout->nParams();
for (size_t i = 0; i < nparam; ++i)
{
TableRow newrow = outbkgdparws->appendRow();
newrow << parnames[i] << funcout->getParameter(i);
}

TableRow chi2row = outbkgdparws->appendRow();
chi2row << "Chi-square" << chi2;

g_log.information() << "Set table workspace (#row = " << outbkgdparws->rowCount()
<< ") to OutputBackgroundParameterTable. " << "\n";
setProperty("OutputBackgroundParameterTable", outbkgdparws);

// Set output workspace
const MantidVec& vecX = m_outputWS->readX(0);
const MantidVec& vecY = m_outputWS->readY(0);
FunctionDomain1DVector domain(vecX);
FunctionValues values(domain);

funcout->function(domain, values);

MantidVec& dataModel = m_outputWS->dataY(1);
MantidVec& dataDiff = m_outputWS->dataY(2);
for (size_t i = 0; i < dataModel.size(); ++i)
{
dataModel[i] = values[i];
dataDiff[i] = vecY[i] - dataModel[i];
}

return;
}


} // namespace CurveFitting
} // namespace Mantid
Expand Down

0 comments on commit 172f805

Please sign in to comment.