Commit 9478f730 authored by Anthony Lim's avatar Anthony Lim
Browse files

updated fitting range for PSI background subtraction

parent 5d284699
......@@ -35,8 +35,9 @@ private:
/// mock.
API::IAlgorithm_sptr setupFitAlgorithm(const std::string &wsName);
virtual std::tuple<double, double>
calculateBackgroundFromFit(API::IAlgorithm_sptr &fit, double maxTime,
int workspaceIndex);
calculateBackgroundFromFit(API::IAlgorithm_sptr &fit,
const std::pair<double, double> &range,
const int &workspaceIndex);
/// Perform validation of inputs to the algorithm
std::map<std::string, std::string> validateInputs() override;
......
......@@ -11,6 +11,9 @@
#include "MantidAPI/WorkspaceFactory.h"
#include "MantidDataObjects/WorkspaceSingleValue.h"
#include "MantidKernel/BoundedValidator.h"
#include "MantidAPI/Run.h"
#include <math.h>
using namespace Mantid::API;
using namespace Mantid::Kernel;
......@@ -22,6 +25,28 @@ namespace Muon {
namespace {
constexpr char *MINIMISER = "Levenberg-Marquardt";
constexpr double FIT_TOLERANCE = 10;
const std::string FIRST_GOOD = "First good spectra ";
const std::string LAST_GOOD = "Last good spectra ";
std::pair<double, double> getRange(const MatrixWorkspace_sptr &inputWorkspace,
const size_t &index) {
auto firstGoodIndex = std::stoi(
inputWorkspace->getLog( FIRST_GOOD+ std::to_string(index))
->value());
auto lastGoodIndex = std::stoi(
inputWorkspace->getLog(LAST_GOOD + std::to_string(index))
->value());
auto midGoodIndex =
int(floor(((double(lastGoodIndex) - double(firstGoodIndex)) / 2.)));
double midGood = inputWorkspace->readX(index)[midGoodIndex];
double lastGood = inputWorkspace->readX(index)[lastGoodIndex];
return std::make_pair(midGood, lastGood);
}
} // namespace
// Register the algorithm into the AlgorithmFactory
......@@ -40,6 +65,11 @@ void PSIBackgroundSubtraction::init() {
declareProperty(
"MaxIterations", 500, mustBePositive,
"Stop after this number of iterations if a good fit is not found");
auto mustBeGreater0 = std::make_shared<Kernel::BoundedValidator<int>>();
mustBeGreater0->setLower(1);
declareProperty("Binning", 1, mustBeGreater0,
"Constant sized rebinning of the data");
}
std::map<std::string, std::string> PSIBackgroundSubtraction::validateInputs() {
......@@ -50,95 +80,130 @@ std::map<std::string, std::string> PSIBackgroundSubtraction::validateInputs() {
if (inputWS->YUnit() != "Counts") {
errors["InputWorkspace"] = "Input Workspace should be a counts workspace.";
}
return errors;
}
const Mantid::API::Run &run = inputWS->run();
for (size_t index = 0; index < inputWS->getNumberHistograms(); index++) {
void PSIBackgroundSubtraction::exec() {
MatrixWorkspace_sptr inputWS = getProperty("InputWorkspace");
// Caclulate and subtract background from inputWS
calculateBackgroundUsingFit(inputWS);
}
/**
* Calculate the background of a PSI workspace by performing a fit, comprising
of a FlatBackground and ExpDecayMuon, on the second half of the PSI data.
* @param inputWorkspace :: Input PSI workspace which is modified inplace
*/
void PSIBackgroundSubtraction::calculateBackgroundUsingFit(
MatrixWorkspace_sptr &inputWorkspace) {
IAlgorithm_sptr fit = setupFitAlgorithm(inputWorkspace->getName());
auto numberOfHistograms = inputWorkspace->getNumberHistograms();
std::vector<double> backgroundValues(numberOfHistograms);
for (size_t index = 0; index < numberOfHistograms; ++index) {
auto maxtime = inputWorkspace->x(index).back();
auto [background, fitQuality] =
calculateBackgroundFromFit(fit, maxtime, static_cast<int>(index));
// If fit quality is poor, do not subtract the background and instead log a
// warning
if (fitQuality > FIT_TOLERANCE) {
g_log.warning()
<< "Fit quality obtained in PSIBackgroundSubtraction is poor. "
<< "Skipping background calculation for WorkspaceIndex: "
<< static_cast<int>(index) << "\n";
} else {
backgroundValues[index] = background;
int firstGood = 0;
int lastGood = 0;
try {firstGood = std::stod(
run.getProperty(FIRST_GOOD + std::to_string(index))->value());
} catch (Kernel::Exception::NotFoundError) {
errors["InputWorkspace"] =
"Input Workspace should should contain first food data. ";
}
try {
lastGood = std::stod(
run.getProperty(LAST_GOOD + std::to_string(index))->value());
} catch (Kernel::Exception::NotFoundError) {
errors["InputWorkspace"] +=
"\n Input Workspace should should contain last food data. ";
}
if (lastGood <= firstGood) {
errors["InputWorkspace"] +=
"\n Input Workspace should have last good data > first good data. ";
}
if (firstGood <= 0) {
errors["InputWorkspace"] +=
"\n Input Workspace should have first good data > 0. ";
}
if (lastGood > inputWS->readX(index).size()) {
errors["InputWorkspace"] +=
"\n Input Workspace should have last good data < number of bins. ";
}
}
// Create background workspace
IAlgorithm_sptr wsAlg = createChildAlgorithm("CreateWorkspace", 0.7, 1.0);
wsAlg->setProperty<std::vector<double>>("DataX", std::vector<double>(2, 0.0));
wsAlg->setProperty<std::vector<double>>("DataY", std::move(backgroundValues));
wsAlg->setProperty<int>("NSpec", static_cast<int>(numberOfHistograms));
wsAlg->execute();
MatrixWorkspace_sptr backgroundWS = wsAlg->getProperty("OutputWorkspace");
backgroundWS->setYUnit("Counts");
// Subtract background from inputworkspace
auto minusAlg = createChildAlgorithm("Minus");
minusAlg->setProperty("LHSWorkspace", inputWorkspace);
minusAlg->setProperty("RHSWorkspace", backgroundWS);
minusAlg->setProperty("OutputWorkspace", inputWorkspace);
minusAlg->execute();
}
/**
* Setup the fit algorithm used to obtain the background from PSI data
of a FlatBackground and ExpDecayMuon, on the second half of the PSI data.
* @param wsName :: Input workspace name
* @return Initalised fitting algorithm
*/
IAlgorithm_sptr
PSIBackgroundSubtraction::setupFitAlgorithm(const std::string &wsName) {
std::string functionstring = "name=FlatBackground,A0=0;name=ExpDecayMuon";
IFunction_sptr func =
FunctionFactory::Instance().createInitialized(functionstring);
IAlgorithm_sptr fit = createChildAlgorithm("Fit");
int maxIterations = getProperty("MaxIterations");
fit->initialize();
fit->setProperty("Function", func);
fit->setProperty("MaxIterations", maxIterations);
fit->setPropertyValue("Minimizer", MINIMISER);
fit->setProperty("CreateOutput", false);
fit->setProperty("InputWorkspace", wsName);
return fit;
}
/**
* Runs the fit algorithm used to obtain the background from PSI data.
* @param fit :: Reference to the child fit algorithm used to obtain the
* background
* @param maxTime :: The maximum value of time (x-axis), used for fit range
* @param workspaceIndex :: Index which the fit will be performed on
* @return A tuple of background and fit quality
*/
std::tuple<double, double> PSIBackgroundSubtraction::calculateBackgroundFromFit(
IAlgorithm_sptr &fit, double maxTime, int workspaceIndex) {
fit->setProperty("StartX", maxTime / 2);
fit->setProperty("EndX", maxTime);
fit->setProperty("WorkspaceIndex", workspaceIndex);
fit->execute();
IFunction_sptr func = fit->getProperty("Function");
double flatbackground = func->getParameter("f0.A0");
double chi2 = std::stod(fit->getPropertyValue("OutputChi2overDof"));
return std::make_tuple(flatbackground, chi2);
}
return errors;
}
void PSIBackgroundSubtraction::exec() {
MatrixWorkspace_sptr inputWS = getProperty("InputWorkspace");
// Caclulate and subtract background from inputWS
calculateBackgroundUsingFit(inputWS);
}
/**
* Calculate the background of a PSI workspace by performing a fit, comprising
of a FlatBackground and ExpDecayMuon, on the second half of the PSI data.
* @param inputWorkspace :: Input PSI workspace which is modified inplace
*/
void PSIBackgroundSubtraction::calculateBackgroundUsingFit(
MatrixWorkspace_sptr & inputWorkspace) {
IAlgorithm_sptr fit = setupFitAlgorithm(inputWorkspace->getName());
auto numberOfHistograms = inputWorkspace->getNumberHistograms();
std::vector<double> backgroundValues(numberOfHistograms);
for (size_t index = 0; index < numberOfHistograms; ++index) {
std::pair<double, double> range = getRange(inputWorkspace, index);
auto [background, fitQuality] =
calculateBackgroundFromFit(fit, range, static_cast<int>(index));
// If fit quality is poor, do not subtract the background and instead log
// a warning
if (fitQuality > FIT_TOLERANCE) {
g_log.warning()
<< "Fit quality obtained in PSIBackgroundSubtraction is poor. "
<< "Skipping background calculation for WorkspaceIndex: "
<< static_cast<int>(index) << "\n";
} else {
backgroundValues[index] = background;
}
}
// Create background workspace
IAlgorithm_sptr wsAlg = createChildAlgorithm("CreateWorkspace", 0.7, 1.0);
wsAlg->setProperty<std::vector<double>>("DataX",
std::vector<double>(2, 0.0));
wsAlg->setProperty<std::vector<double>>("DataY",
std::move(backgroundValues));
wsAlg->setProperty<int>("NSpec", static_cast<int>(numberOfHistograms));
wsAlg->execute();
MatrixWorkspace_sptr backgroundWS = wsAlg->getProperty("OutputWorkspace");
backgroundWS->setYUnit("Counts");
// Subtract background from inputworkspace
auto minusAlg = createChildAlgorithm("Minus");
minusAlg->setProperty("LHSWorkspace", inputWorkspace);
minusAlg->setProperty("RHSWorkspace", backgroundWS);
minusAlg->setProperty("OutputWorkspace", inputWorkspace);
minusAlg->execute();
}
/**
* Setup the fit algorithm used to obtain the background from PSI data
of a FlatBackground and ExpDecayMuon, on the second half of the PSI data.
* @param wsName :: Input workspace name
* @return Initalised fitting algorithm
*/
IAlgorithm_sptr PSIBackgroundSubtraction::setupFitAlgorithm(
const std::string &wsName) {
std::string functionstring = "name=FlatBackground,A0=0;name=ExpDecayMuon";
IFunction_sptr func =
FunctionFactory::Instance().createInitialized(functionstring);
IAlgorithm_sptr fit = createChildAlgorithm("Fit");
int maxIterations = getProperty("MaxIterations");
fit->initialize();
fit->setProperty("Function", func);
fit->setProperty("MaxIterations", maxIterations);
fit->setPropertyValue("Minimizer", MINIMISER);
fit->setProperty("CreateOutput", false);
fit->setProperty("InputWorkspace", wsName);
return fit;
}
/**
* Runs the fit algorithm used to obtain the background from PSI data.
* @param fit :: Reference to the child fit algorithm used to obtain the
* background
* @param maxTime :: The maximum value of time (x-axis), used for fit range
* @param workspaceIndex :: Index which the fit will be performed on
* @return A tuple of background and fit quality
*/
std::tuple<double, double>
PSIBackgroundSubtraction::calculateBackgroundFromFit(
IAlgorithm_sptr & fit, const std::pair<double, double> &range,
const int &workspaceIndex) {
fit->setProperty("StartX", range.first);
fit->setProperty("EndX", range.second);
fit->setProperty("WorkspaceIndex", workspaceIndex);
fit->execute();
IFunction_sptr func = fit->getProperty("Function");
double flatbackground = func->getParameter("f0.A0");
double chi2 = std::stod(fit->getPropertyValue("OutputChi2overDof"));
return std::make_tuple(flatbackground, chi2);
}
} // namespace Muon
} // namespace Muon
} // namespace Mantid
......@@ -25,6 +25,13 @@ MatrixWorkspace_sptr createCountsTestWorkspace(const size_t numberOfHistograms,
MatrixWorkspace_sptr ws = WorkspaceCreationHelper::create2DWorkspace(
numberOfHistograms, numberOfBins);
ws->setYUnit("Counts");
for (int index = 0; index < numberOfHistograms; index++) {
ws->mutableRun().addProperty("First good spectra "+std::to_string(index), int(numberOfBins/2.));
ws->mutableRun().addProperty("Last good spectra " +
std::to_string(index), numberOfBins);
}
AnalysisDataService::Instance().addOrReplace(WORKSPACE_NAME, ws);
return ws;
......@@ -58,7 +65,7 @@ public:
private:
std::tuple<double, double> calculateBackgroundFromFit(IAlgorithm_sptr &,
double, int) override {
const std::pair<double,double> &range, const int &index) override {
return std::make_tuple(m_background, m_fitQuality);
}
double m_background{0.00};
......
......@@ -18,11 +18,12 @@ The background is removed from the original data through the following formula,
where :math:`y_{cor}(t)` is the corrected bin data, :math:`y_{org}(t)` the original counts data and :math:`B` is the
flat background.
To obtain the flat background, :math:`B`, the second-half of the raw-data is fitted with the following function,
To obtain the flat background, :math:`B`, the second-half of the good raw-data is fitted with the following function,
.. math:: f(t) = \mbox{A}e^{-\lambda t} + B
where the first term represents a ExpDecay function, see :ref:`func-ExpDecayMuon`.
where the first term represents a ExpDecay function, see :ref:`func-ExpDecayMuon`.
The good raw-data is defined as being from the bin containing the `first good data` to the bin containing the `last good data`.
The algorithm takes in an input workspace and performs the correction inplace on the workspace. The number of iterations
can be specified through an optional parameter. If a poor quality fit is returned, a warning will be displayed in the
......@@ -47,6 +48,9 @@ Usage
# Create workspaces
input_workspace = CreateWorkspace(time, counts)
input_workspace.setYUnit("Counts")
run = input_workspace.getRun()
run.addProperty("First good spectra 0","0","None",True)
run.addProperty("Last good spectra 0","10","None",True)
workspace_copy = input_workspace.clone()
# Run PSIBackgroundSubtraction Algorithm
......
Supports Markdown
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