PSIBackgroundSubtraction.cpp 7.82 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
// Mantid Repository : https://github.com/mantidproject/mantid
//
// Copyright © 2020 ISIS Rutherford Appleton Laboratory UKRI,
//   NScD Oak Ridge National Laboratory, European Spallation Source,
//   Institut Laue - Langevin & CSNS, Institute of High Energy Physics, CAS
// SPDX - License - Identifier: GPL - 3.0 +
#include "MantidMuon/PSIBackgroundSubtraction.h"
#include "MantidAPI/FunctionFactory.h"
#include "MantidAPI/IFunction.h"
#include "MantidAPI/MatrixWorkspace.h"
11
#include "MantidAPI/Run.h"
12
13
#include "MantidAPI/WorkspaceFactory.h"
#include "MantidDataObjects/WorkspaceSingleValue.h"
14
#include "MantidKernel/BoundedValidator.h"
15
16

#include <math.h>
17
18
19
20
21
22
23
24
25
26
27

using namespace Mantid::API;
using namespace Mantid::Kernel;
using namespace Mantid::DataObjects;

namespace Mantid {
namespace Muon {

namespace {
constexpr char *MINIMISER = "Levenberg-Marquardt";
constexpr double FIT_TOLERANCE = 10;
28
29
30
31
32
33
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(
34
      inputWorkspace->getLog(FIRST_GOOD + std::to_string(index))->value());
35
  auto lastGoodIndex = std::stoi(
36
      inputWorkspace->getLog(LAST_GOOD + std::to_string(index))->value());
37
38
39
40
41
42
43
44
45
46
47

  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);
}

48
49
50
51
52
53
54
55
56
} // namespace

// Register the algorithm into the AlgorithmFactory
DECLARE_ALGORITHM(PSIBackgroundSubtraction)

void PSIBackgroundSubtraction::init() {

  declareProperty(
      std::make_unique<WorkspaceProperty<>>(
57
          "InputWorkspace", "", Direction::InOut, PropertyMode::Mandatory),
58
59
60
61
62
63
64
65
      "Input workspace containing the PSI bin data "
      "which the background correction will be applied to.");

  auto mustBePositive = std::make_shared<Kernel::BoundedValidator<int>>();
  mustBePositive->setLower(0);
  declareProperty(
      "MaxIterations", 500, mustBePositive,
      "Stop after this number of iterations if a good fit is not found");
66
67
68
69
70

  auto mustBeGreater0 = std::make_shared<Kernel::BoundedValidator<int>>();
  mustBeGreater0->setLower(1);
  declareProperty("Binning", 1, mustBeGreater0,
                  "Constant sized rebinning of the data");
71
72
73
74
75
76
77
}

std::map<std::string, std::string> PSIBackgroundSubtraction::validateInputs() {
  std::map<std::string, std::string> errors;

  MatrixWorkspace_sptr inputWS = getProperty("InputWorkspace");

78
  if (inputWS->YUnit() != "Counts") {
79
80
    errors["InputWorkspace"] = "Input Workspace should be a counts workspace.";
  }
81
82
  const Mantid::API::Run &run = inputWS->run();
  for (size_t index = 0; index < inputWS->getNumberHistograms(); index++) {
83

84
85
86
87
88
    int firstGood = 0;
    int lastGood = 0;
    try {
      firstGood = std::stod(
          run.getProperty(FIRST_GOOD + std::to_string(index))->value());
89
90
91
92
93
94
95
96
97
98
99
100
101
102
    } 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. ";
103
    }
104
    if (firstGood <= 0) {
105
106
107
      errors["InputWorkspace"] +=
          "\n Input Workspace should have first good data > 0. ";
    }
108
109
110
    if (lastGood > inputWS->readX(index).size()) {
      errors["InputWorkspace"] +=
          "\n Input Workspace should have last good data < number of bins. ";
111
112
    }
  }
113
114
  return errors;
}
115

116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
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;
144
145
    }
  }
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
  // 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);
}
204
} // namespace Muon
205
} // namespace Mantid