Unverified Commit fef6ad28 authored by Nick Draper's avatar Nick Draper Committed by GitHub
Browse files

Merge pull request #28798 from mantidproject/PSI_BG

Psi background
parents 3f4c5821 7e5ffdbc
......@@ -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;
......
......@@ -8,10 +8,13 @@
#include "MantidAPI/FunctionFactory.h"
#include "MantidAPI/IFunction.h"
#include "MantidAPI/MatrixWorkspace.h"
#include "MantidAPI/Run.h"
#include "MantidAPI/WorkspaceFactory.h"
#include "MantidDataObjects/WorkspaceSingleValue.h"
#include "MantidKernel/BoundedValidator.h"
#include <math.h>
using namespace Mantid::API;
using namespace Mantid::Kernel;
using namespace Mantid::DataObjects;
......@@ -22,6 +25,26 @@ 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 +63,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,6 +78,38 @@ std::map<std::string, std::string> PSIBackgroundSubtraction::validateInputs() {
if (inputWS->YUnit() != "Counts") {
errors["InputWorkspace"] = "Input Workspace should be a counts workspace.";
}
const Mantid::API::Run &run = inputWS->run();
for (size_t index = 0; index < inputWS->getNumberHistograms(); index++) {
int firstGood = 0;
int lastGood = 0;
try {
firstGood = std::stoi(
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::stoi(
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 >= int(inputWS->readX(index).size())) {
errors["InputWorkspace"] +=
"\n Input Workspace should have last good data < number of bins. ";
}
}
return errors;
}
......@@ -69,11 +129,11 @@ void PSIBackgroundSubtraction::calculateBackgroundUsingFit(
auto numberOfHistograms = inputWorkspace->getNumberHistograms();
std::vector<double> backgroundValues(numberOfHistograms);
for (size_t index = 0; index < numberOfHistograms; ++index) {
auto maxtime = inputWorkspace->x(index).back();
std::pair<double, double> range = getRange(inputWorkspace, index);
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
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. "
......@@ -130,9 +190,10 @@ PSIBackgroundSubtraction::setupFitAlgorithm(const std::string &wsName) {
* @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);
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");
......
......@@ -20,11 +20,22 @@ namespace {
constexpr char *WORKSPACE_NAME = "DummyWS";
MatrixWorkspace_sptr createCountsTestWorkspace(const size_t numberOfHistograms,
const size_t numberOfBins) {
const size_t numberOfBins,
bool addLogs = true) {
MatrixWorkspace_sptr ws = WorkspaceCreationHelper::create2DWorkspace(
numberOfHistograms, numberOfBins);
ws->setYUnit("Counts");
if (addLogs) {
for (size_t index = 0; index < numberOfHistograms; index++) {
ws->mutableRun().addProperty("First good spectra " +
std::to_string(index),
int(double(numberOfBins) / 2.));
ws->mutableRun().addProperty("Last good spectra " + std::to_string(index),
numberOfBins);
}
}
AnalysisDataService::Instance().addOrReplace(WORKSPACE_NAME, ws);
return ws;
......@@ -57,8 +68,12 @@ public:
}
private:
std::tuple<double, double> calculateBackgroundFromFit(IAlgorithm_sptr &,
double, int) override {
std::tuple<double, double>
calculateBackgroundFromFit(IAlgorithm_sptr &,
const std::pair<double, double> &range,
const int &index) override {
(void)range;
(void)index;
return std::make_tuple(m_background, m_fitQuality);
}
double m_background{0.00};
......@@ -93,6 +108,77 @@ public:
clearADS();
}
void test_that_algorithm_does_not_execute_if_no_good_data() {
PSIBackgroundSubtraction alg;
auto ws = createCountsTestWorkspace(2, 100, false);
alg.initialize();
alg.setProperty("InputWorkspace", ws);
TS_ASSERT_THROWS(alg.execute(), const std::runtime_error &);
clearADS();
}
void test_that_algorithm_does_not_execute_if_bad_first_good_data() {
PSIBackgroundSubtraction alg;
int numberOfHistograms = 2;
int numberOfBins = 100;
auto ws =
createCountsTestWorkspace(numberOfHistograms, numberOfBins, false);
for (int index = 0; index < numberOfHistograms; index++) {
ws->mutableRun().addProperty(
"First good spectra " + std::to_string(index), -1);
ws->mutableRun().addProperty("Last good spectra " + std::to_string(index),
numberOfBins - 10);
}
alg.initialize();
alg.setProperty("InputWorkspace", ws);
TS_ASSERT_THROWS(alg.execute(), const std::runtime_error &);
clearADS();
}
void test_that_algorithm_does_not_execute_if_bad_last_good_data() {
PSIBackgroundSubtraction alg;
int numberOfHistograms = 2;
int numberOfBins = 100;
auto ws =
createCountsTestWorkspace(numberOfHistograms, numberOfBins, false);
for (int index = 0; index < numberOfHistograms; index++) {
ws->mutableRun().addProperty(
"First good spectra " + std::to_string(index), 1);
ws->mutableRun().addProperty("Last good spectra " + std::to_string(index),
numberOfBins * 2);
}
alg.initialize();
alg.setProperty("InputWorkspace", ws);
TS_ASSERT_THROWS(alg.execute(), const std::runtime_error &);
clearADS();
}
void test_that_algorithm_does_not_execute_if_last_before_first_good_data() {
PSIBackgroundSubtraction alg;
int numberOfHistograms = 2;
int numberOfBins = 100;
auto ws =
createCountsTestWorkspace(numberOfHistograms, numberOfBins, false);
for (int index = 0; index < numberOfHistograms; index++) {
ws->mutableRun().addProperty(
"First good spectra " + std::to_string(index), 50);
ws->mutableRun().addProperty("Last good spectra " + std::to_string(index),
40);
}
alg.initialize();
alg.setProperty("InputWorkspace", ws);
TS_ASSERT_THROWS(alg.execute(), const std::runtime_error &);
clearADS();
}
void test_background_correctly_removed_from_input_workspace() {
MockPSIBackgroundSubtraction alg;
double background = 20;
......
......@@ -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",10,"None",True)
run.addProperty("Last good spectra 0",99,"None",True)
workspace_copy = input_workspace.clone()
# Run PSIBackgroundSubtraction Algorithm
......@@ -54,8 +58,9 @@ Usage
# Find the difference between the workspaces
workspace_diff = Minus(workspace_copy, input_workspace)
diffs = np.round(workspace_diff.readY(0),4)
# The counts in workspace diff should be a flat line corresponding to the background
print("Differences in counts are: {}".format(workspace_diff.dataY(0)))
print("Differences in counts are: {}".format(diffs))
Output:
......
Markdown is supported
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