Mantid
Loading...
Searching...
No Matches
RemoveBackground.cpp
Go to the documentation of this file.
1// Mantid Repository : https://github.com/mantidproject/mantid
2//
3// Copyright © 2018 ISIS Rutherford Appleton Laboratory UKRI,
4// NScD Oak Ridge National Laboratory, European Spallation Source,
5// Institut Laue - Langevin & CSNS, Institute of High Energy Physics, CAS
6// SPDX - License - Identifier: GPL - 3.0 +
8
9#include "MantidAPI/Axis.h"
12#include "MantidAPI/Run.h"
22#include "MantidKernel/Unit.h"
25
26using Mantid::HistogramData::HistogramE;
27using Mantid::HistogramData::HistogramX;
28using Mantid::HistogramData::HistogramY;
29
30namespace Mantid::Algorithms {
31
32// Register the class into the algorithm factory
33DECLARE_ALGORITHM(RemoveBackground)
34
35using namespace Kernel;
36using namespace API;
37
42 auto sourceValidator = std::make_shared<CompositeValidator>();
43 sourceValidator->add<InstrumentValidator>();
44 sourceValidator->add<HistogramValidator>();
45 declareProperty(std::make_unique<WorkspaceProperty<>>("InputWorkspace", "", Direction::Input, sourceValidator),
46 "Workspace containing the input data");
47 declareProperty(std::make_unique<WorkspaceProperty<>>("OutputWorkspace", "", Direction::Output),
48 "The name to give the output workspace");
49
50 auto vsValidator = std::make_shared<CompositeValidator>();
51 vsValidator->add<WorkspaceUnitValidator>("TOF");
52 vsValidator->add<HistogramValidator>();
53 declareProperty(std::make_unique<WorkspaceProperty<>>("BkgWorkspace", "", Direction::Input, vsValidator),
54 "An optional histogram workspace in the units of TOF defining background "
55 "for removal during rebinning."
56 "The workspace has to have single value or contain the same number of "
57 "spectra as the \"InputWorkspace\" and single Y value per each spectra,"
58 "representing flat background in the background time region. "
59 "If such workspace is present, the value of the flat background provided "
60 "by this workspace is removed "
61 "from each spectra of the rebinned workspace. This works for histogram "
62 "and event workspace when events are not retained "
63 "but actually useful mainly for removing background while rebinning an "
64 "event workspace in the units different from TOF.");
65
66 std::vector<std::string> dE_modes = Kernel::DeltaEMode::availableTypes();
68 std::make_shared<Kernel::StringListValidator>(dE_modes),
69 "The energy conversion mode used to define the conversion "
70 "from the units of the InputWorkspace to TOF",
72 declareProperty("NullifyNegativeValues", false,
73 "When background is subtracted, signals in some time "
74 "channels may become negative.\n"
75 "If this option is true, signal in such bins is nullified "
76 "and the module of the removed signal"
77 "is added to the error. If false, the negative signal and "
78 "correspondent errors remain unchanged",
80}
81
88
89 // Get the input workspace
90 MatrixWorkspace_sptr inputWS = getProperty("InputWorkspace");
91 MatrixWorkspace_sptr outputWS = getProperty("OutputWorkspace");
92
93 API::MatrixWorkspace_const_sptr bkgWksp = getProperty("BkgWorkspace");
94 bool nullifyNegative = getProperty("NullifyNegativeValues");
95
96 if (!(bkgWksp->getNumberHistograms() == 1 || inputWS->getNumberHistograms() == bkgWksp->getNumberHistograms())) {
97 throw std::invalid_argument(" Background Workspace: " + bkgWksp->getName() +
98 " should have the same number of spectra as "
99 "source workspace or be a single histogram "
100 "workspace");
101 }
102
103 const std::string emodeStr = getProperty("EMode");
104 auto eMode = Kernel::DeltaEMode::fromString(emodeStr);
105
106 // Removing background in-place
107 bool inPlace = (inputWS == outputWS);
108 // workspace independent determination of length
109 const auto histnumber = static_cast<int>(inputWS->getNumberHistograms());
110
111 if (!inPlace) {
112 // make the copy of the output Workspace from the input. Also copies
113 // X-vectors and axis
114 outputWS = API::WorkspaceFactory::Instance().create(inputWS);
115 }
116
117 int nThreads = PARALLEL_GET_MAX_THREADS;
118 m_BackgroundHelper.initialize(bkgWksp, inputWS, eMode, &g_log, nThreads, inPlace, nullifyNegative);
119
120 Progress prog(this, 0.0, 1.0, histnumber);
121 PARALLEL_FOR_IF(Kernel::threadSafe(*inputWS, *outputWS))
122 for (int hist = 0; hist < histnumber; ++hist) {
124 // get references to output Workspace X-arrays.
125 auto &XValues = outputWS->mutableX(hist);
126 // get references to output workspace data and error. If it is new
127 // workspace, data will be copied there, if old, modified in-place
128 auto &YValues = outputWS->mutableY(hist);
129 auto &YErrors = outputWS->mutableE(hist);
130
131 // output data arrays are implicitly filled by function
132 int id = PARALLEL_THREAD_NUMBER;
133 m_BackgroundHelper.removeBackground(hist, XValues, YValues, YErrors, id);
134
135 prog.report(name());
137 }
139 // Assign it to the output workspace property
140 setProperty("OutputWorkspace", outputWS);
141}
142//-------------------------------------------------------------------------------------------------------------------------------
143//---------------- BACKGROUND HELPER
144//---------------------------------------------------------
145//-------------------------------------------------------------------------------------------------------------------------------
148 : m_WSUnit(), m_bgWs(), m_wkWS(), m_spectrumInfo(nullptr), m_pgLog(nullptr), m_inPlace(true),
149 m_singleValueBackground(false), m_NBg(0), m_dtBg(1), m_ErrSq(0), m_Emode(DeltaEMode::Elastic),
150 m_nullifyNegative(false), m_previouslyRemovedBkgMode(false) {}
151
166 Kernel::Logger *pLog, int nThreads, bool inPlace, bool nullifyNegative) {
167 m_bgWs = bkgWS;
168 m_wkWS = sourceWS;
169 m_Emode = emode;
170 m_pgLog = pLog;
171 m_inPlace = inPlace;
172 m_nullifyNegative = nullifyNegative;
173
174 std::string bgUnits = bkgWS->getAxis(0)->unit()->unitID();
175 if (bgUnits != "TOF")
176 throw std::invalid_argument(" Background Workspace: " + bkgWS->getName() + " should be in the units of TOF");
177
178 if (!(bkgWS->getNumberHistograms() == 1 || sourceWS->getNumberHistograms() == bkgWS->getNumberHistograms()))
179 throw std::invalid_argument(" Background Workspace: " + bkgWS->getName() +
180 " should have the same number of spectra as "
181 "source workspace or be a single histogram "
182 "workspace");
183
184 auto WSUnit = sourceWS->getAxis(0)->unit();
185 if (!WSUnit)
186 throw std::invalid_argument(" Source Workspace: " + sourceWS->getName() + " should have units");
187
188 m_spectrumInfo = &sourceWS->spectrumInfo();
189
190 // allocate the array of units converters to avoid units reallocation within a
191 // loop
192 m_WSUnit.clear();
193 m_WSUnit.reserve(nThreads);
194 for (int i = 0; i < nThreads; i++) {
195 m_WSUnit.emplace_back(std::unique_ptr<Kernel::Unit>(WSUnit->clone()));
196 }
197
199 if (bkgWS->getNumberHistograms() <= 1)
201
202 const auto &dataX = bkgWS->x(0);
203 const auto &dataY = bkgWS->y(0);
204 const auto &dataE = bkgWS->e(0);
205 m_NBg = dataY[0];
206 m_dtBg = dataX[1] - dataX[0];
207 m_ErrSq = dataE[0] * dataE[0]; // needs further clarification
208
211 if (m_NBg < 1.e-7 && m_ErrSq < 1.e-7)
213 }
214}
226void BackgroundHelper::removeBackground(int nHist, HistogramX &x_data, HistogramY &y_data, HistogramE &e_data,
227 int threadNum) const {
228
229 double dtBg, IBg, ErrBgSq;
231 dtBg = m_dtBg;
232 ErrBgSq = m_ErrSq;
233 IBg = m_NBg;
234 } else {
235 auto &dataX = m_bgWs->x(nHist);
236 auto &dataY = m_bgWs->y(nHist);
237 auto &dataE = m_bgWs->e(nHist);
238 dtBg = (dataX[1] - dataX[0]);
239 IBg = dataY[0];
240 ErrBgSq = dataE[0] * dataE[0]; // Needs further clarification
241 }
242
243 try {
244 double L1 = m_spectrumInfo->l1();
245
246 // get access to source workspace in case if target is different from source
247 auto &XValues = m_wkWS->x(nHist);
248 auto &YValues = m_wkWS->y(nHist);
249 auto &YErrors = m_wkWS->e(nHist);
250
251 // use thread-specific unit conversion class to avoid multithreading issues
252 Kernel::Unit *unitConv = m_WSUnit[threadNum].get();
254 m_spectrumInfo->getDetectorValues(*unitConv, Kernel::Units::TOF{}, m_Emode, false, nHist, pmap);
255 unitConv->initialize(L1, m_Emode, pmap);
256
257 x_data[0] = XValues[0];
258 double tof1 = unitConv->singleToTOF(x_data[0]);
259 for (size_t i = 0; i < y_data.size(); i++) {
260 double X = XValues[i + 1];
261 double tof2 = unitConv->singleToTOF(X);
262 double Jack = std::fabs((tof2 - tof1) / dtBg);
263 double normBkgrnd = IBg * Jack;
264 double errBkgSq = ErrBgSq * Jack * Jack;
265 tof1 = tof2;
266 if (m_inPlace) {
267 y_data[i] -= normBkgrnd;
268 // e_data[i] =std::sqrt(ErrSq*Jack*Jack+e_data[i]*e_data[i]); //
269 // needs further clarification -- Gaussian error summation?
270 //--> assume error for background is sqrt(signal):
271 e_data[i] = std::sqrt((errBkgSq + e_data[i] * e_data[i]));
272 } else {
273 x_data[i + 1] = X;
274 y_data[i] = YValues[i] - normBkgrnd;
275 e_data[i] = std::sqrt((errBkgSq + YErrors[i] * YErrors[i]));
276 }
277 if (m_nullifyNegative && y_data[i] < 0) {
279 // background have been removed
280 // and not we remove negative signal and estimate errors
281 double prev_rem_bkg = -y_data[i];
282 e_data[i] = e_data[i] > prev_rem_bkg ? e_data[i] : prev_rem_bkg;
283 } else {
287 e_data[i] = e_data[i] > normBkgrnd ? e_data[i] : normBkgrnd;
288 }
289 y_data[i] = 0;
290 }
291 }
292
293 } catch (...) {
294 // no background removal for this spectra as it does not have a detector or
295 // other reason
296 if (m_pgLog)
297 m_pgLog->debug() << " Can not remove background for the spectra with number (id)" << nHist;
298 }
299}
300
301} // namespace Mantid::Algorithms
#define DECLARE_ALGORITHM(classname)
Definition: Algorithm.h:576
#define PARALLEL_THREAD_NUMBER
#define PARALLEL_START_INTERRUPT_REGION
Begins a block to skip processing is the algorithm has been interupted Note the end of the block if n...
Definition: MultiThreaded.h:94
#define PARALLEL_END_INTERRUPT_REGION
Ends a block to skip processing is the algorithm has been interupted Note the start of the block if n...
#define PARALLEL_FOR_IF(condition)
Empty definitions - to enable set your complier to enable openMP.
#define PARALLEL_GET_MAX_THREADS
#define PARALLEL_CHECK_INTERRUPT_REGION
Adds a check after a Parallel region to see if it was interupted.
void declareProperty(std::unique_ptr< Kernel::Property > p, const std::string &doc="") override
Add a property to the list of managed properties.
Definition: Algorithm.cpp:1913
TypedValue getProperty(const std::string &name) const override
Get the value of a property.
Definition: Algorithm.cpp:2076
Kernel::Logger & g_log
Definition: Algorithm.h:451
A validator which checks that a workspace contains histogram data (the default) or point data as requ...
A validator which checks that a workspace has a valid instrument.
Helper class for reporting progress from algorithms.
Definition: Progress.h:25
double l1() const
Returns L1 (distance from source to sample).
void getDetectorValues(const Kernel::Unit &inputUnit, const Kernel::Unit &outputUnit, const Kernel::DeltaEMode::Type emode, const bool signedTheta, int64_t wsIndex, Kernel::UnitParametersMap &pmap) const
Get the detector values relevant to unit conversion for a workspace index.
A property class for workspaces.
A validator which checks that the unit of the workspace referred to by a WorkspaceProperty is the exp...
Kernel::DeltaEMode::Type m_Emode
const API::SpectrumInfo * m_spectrumInfo
std::vector< std::unique_ptr< Kernel::Unit > > m_WSUnit
void removeBackground(int nHist, HistogramData::HistogramX &x_data, HistogramData::HistogramY &y_data, HistogramData::HistogramE &e_data, int threadNum=0) const
Method removes background from vectors which represent a histogram data for a single spectra.
API::MatrixWorkspace_const_sptr m_wkWS
void initialize(const API::MatrixWorkspace_const_sptr &bkgWS, const API::MatrixWorkspace_sptr &sourceWS, Kernel::DeltaEMode::Type emode, Kernel::Logger *pLog=nullptr, int nThreads=1, bool inPlace=true, bool nullifyNegative=false)
Initialization method:
API::MatrixWorkspace_const_sptr m_bgWs
void init() override
Initialization method.
void exec() override
Executes the rebin algorithm.
const std::string name() const override
Algorithm's name for identification overriding a virtual method.
IPropertyManager * setProperty(const std::string &name, const T &value)
Templated method to set the value of a PropertyWithValue.
The Logger class is in charge of the publishing messages from the framework through various channels.
Definition: Logger.h:52
void debug(const std::string &msg)
Logs at debug level.
Definition: Logger.cpp:114
void report()
Increments the loop counter by 1, then sends the progress notification on behalf of its algorithm.
Definition: ProgressBase.h:51
static T & Instance()
Return a reference to the Singleton instance, creating it if it does not already exist Creation is do...
The base units (abstract) class.
Definition: Unit.h:41
virtual double singleToTOF(const double x) const =0
Convert a single X value to TOF.
void initialize(const double &_l1, const int &_emode, const UnitParametersMap &params)
Initialize the unit to perform conversion using singleToTof() and singleFromTof()
Definition: Unit.cpp:132
Time of flight in microseconds.
Definition: Unit.h:283
std::shared_ptr< const MatrixWorkspace > MatrixWorkspace_const_sptr
shared pointer to the matrix workspace base class (const version)
std::shared_ptr< MatrixWorkspace > MatrixWorkspace_sptr
shared pointer to the matrix workspace base class
std::unordered_map< UnitParams, double > UnitParametersMap
Definition: Unit.h:30
std::enable_if< std::is_pointer< Arg >::value, bool >::type threadSafe(Arg workspace)
Thread-safety check Checks the workspace to ensure it is suitable for multithreaded access.
Definition: MultiThreaded.h:22
Generate a tableworkspace to store the calibration results.
Defines the possible energy transfer modes:
Definition: DeltaEMode.h:23
static const std::vector< std::string > availableTypes()
Returns the string list of available modes.
Definition: DeltaEMode.cpp:35
static Type fromString(const std::string &modeStr)
Returns the emode from the given string.
Definition: DeltaEMode.cpp:69
Type
Define the available energy transfer modes It is important to assign enums proper numbers,...
Definition: DeltaEMode.h:29
@ Input
An input workspace.
Definition: Property.h:53
@ Output
An output workspace.
Definition: Property.h:54