Mantid
Loading...
Searching...
No Matches
NormaliseByPeakArea.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/IFunction.h"
17#include "MantidHistogramData/HistogramE.h"
18#include "MantidHistogramData/HistogramY.h"
19
22
23#include <memory>
24
26
28double PEAK_POS_GUESS = -0.1;
30double PEAK_WIDTH_GUESS = 4.0;
32double SUMMEDY_BIN_WIDTH = 0.5;
33
34// Register the algorithm into the AlgorithmFactory
36
37using namespace API;
38using namespace Kernel;
39
40//----------------------------------------------------------------------------------------------
43NormaliseByPeakArea::NormaliseByPeakArea()
44 : API::Algorithm(), m_inputWS(), m_mass(0.0), m_sumResults(true), m_normalisedWS(), m_yspaceWS(), m_fittedWS(),
45 m_symmetrisedWS(), m_progress() {}
46
47//----------------------------------------------------------------------------------------------
49const std::string NormaliseByPeakArea::name() const { return "NormaliseByPeakArea"; }
50
52int NormaliseByPeakArea::version() const { return 1; }
53
55const std::string NormaliseByPeakArea::category() const { return "CorrectionFunctions\\NormalisationCorrections"; }
56
57//----------------------------------------------------------------------------------------------
58
59//----------------------------------------------------------------------------------------------
62void NormaliseByPeakArea::init() {
63 auto wsValidator = std::make_shared<CompositeValidator>();
64 wsValidator->add<HistogramValidator>(false); // point data
65 wsValidator->add<InstrumentValidator>();
66 wsValidator->add<WorkspaceUnitValidator>("TOF");
67 declareProperty(std::make_unique<WorkspaceProperty<>>("InputWorkspace", "", Direction::Input, wsValidator),
68 "An input workspace.");
69
70 auto mustBePositive = std::make_shared<BoundedValidator<double>>();
71 mustBePositive->setLower(0.0);
72 mustBePositive->setLowerExclusive(true); // strictly greater than 0.0
73 declareProperty("Mass", -1.0, mustBePositive, "The mass, in AMU, defining the recoil peak to fit");
74 declareProperty("Sum", true,
75 "If true all spectra on the Y-space, fitted & symmetrised workspaces "
76 "are summed in quadrature to produce the final result");
77
78 declareProperty(std::make_unique<WorkspaceProperty<>>("OutputWorkspace", "", Direction::Output),
79 "Input workspace normalised by the fitted peak area");
80 declareProperty(std::make_unique<WorkspaceProperty<>>("YSpaceDataWorkspace", "", Direction::Output),
81 "Input workspace converted to units of Y-space");
82 declareProperty(std::make_unique<WorkspaceProperty<>>("FittedWorkspace", "", Direction::Output),
83 "Output from fit of the single mass peakin y-space. The output units are "
84 "in momentum (A^-1)");
85 declareProperty(std::make_unique<WorkspaceProperty<>>("SymmetrisedWorkspace", "", Direction::Output),
86 "The input data symmetrised about Y=0. The output units are in momentum "
87 "(A^-1)");
88}
89
90//----------------------------------------------------------------------------------------------
93void NormaliseByPeakArea::exec() {
95 const auto yspaceIn = convertInputToY();
96 createOutputWorkspaces(yspaceIn);
97
98 const auto nhist = static_cast<int64_t>(yspaceIn->getNumberHistograms());
99 const auto nreports = static_cast<int64_t>(yspaceIn->getNumberHistograms() +
100 2 * m_symmetrisedWS->getNumberHistograms() * m_symmetrisedWS->blocksize());
101 m_progress = std::make_unique<API::Progress>(this, 0.10, 1.0, nreports);
102
103 for (int64_t i = 0; i < nhist; ++i) {
104 m_normalisedWS->setSharedX(i, m_inputWS->sharedX(i)); // TOF
105 if (!m_sumResults) // avoid setting multiple times if we are summing
106 {
107 m_yspaceWS->setSharedX(i, yspaceIn->sharedX(i)); // momentum
108 m_fittedWS->setSharedX(i, yspaceIn->sharedX(i)); // momentum
109 m_symmetrisedWS->setSharedX(i, yspaceIn->sharedX(i)); // momentum
110 }
111
112 double peakArea = fitToMassPeak(yspaceIn, static_cast<size_t>(i));
113 normaliseTOFData(peakArea, i);
114 saveToOutput(m_yspaceWS, yspaceIn->sharedY(i), yspaceIn->sharedE(i), i);
115
116 m_progress->report();
117 }
118 // This has to be done after the summation of the spectra
120
121 setProperty("OutputWorkspace", m_normalisedWS);
122 setProperty("YSpaceDataWorkspace", m_yspaceWS);
123 setProperty("FittedWorkspace", m_fittedWS);
124 setProperty("SymmetrisedWorkspace", m_symmetrisedWS);
125}
126
130void NormaliseByPeakArea::retrieveInputs() {
131 m_inputWS = getProperty("InputWorkspace");
132 m_mass = getProperty("Mass");
133 m_sumResults = getProperty("Sum");
134}
135
140void NormaliseByPeakArea::createOutputWorkspaces(const API::MatrixWorkspace_sptr &yspaceIn) {
141 m_normalisedWS = WorkspaceFactory::Instance().create(m_inputWS); // TOF data is not resized
142
143 const size_t nhist = m_sumResults ? 1 : yspaceIn->getNumberHistograms();
144
145 m_yspaceWS = WorkspaceFactory::Instance().create(yspaceIn, nhist);
146 m_fittedWS = WorkspaceFactory::Instance().create(yspaceIn, nhist);
147 m_symmetrisedWS = WorkspaceFactory::Instance().create(yspaceIn, nhist);
148 if (m_sumResults) {
149 // Copy over xvalues & assign "high" initial error values to simplify
150 // symmetrisation calculation
151 constexpr double high(1e6);
152
153 m_yspaceWS->setSharedX(0, yspaceIn->sharedX(0));
154 m_fittedWS->setSharedX(0, yspaceIn->sharedX(0));
155 m_symmetrisedWS->setSharedX(0, yspaceIn->sharedX(0));
156 m_yspaceWS->mutableE(0) = high;
157 m_fittedWS->setSharedE(0, m_yspaceWS->sharedE(0));
158 m_symmetrisedWS->setSharedE(0, m_yspaceWS->sharedE(0));
159 }
160
164}
165
169void NormaliseByPeakArea::setUnitsToMomentum(const API::MatrixWorkspace_sptr &workspace) {
170 // Units
171 auto xLabel = std::make_shared<Units::Label>("Momentum", "A^-1");
172 workspace->getAxis(0)->unit() = xLabel;
173 workspace->setYUnit("");
174 workspace->setYUnitLabel("");
175}
176
177/*
178 * Returns a workspace converted to Y-space coordinates. @see ConvertToYSpace.
179 * If summing is requested
180 * then the output is rebinned to a common grid to allow summation onto a common
181 * grid. The rebin min/max
182 * is found from the converted workspace
183 */
184MatrixWorkspace_sptr NormaliseByPeakArea::convertInputToY() {
185 auto alg = createChildAlgorithm("ConvertToYSpace", 0.0, 0.05, false);
186 alg->setProperty("InputWorkspace", m_inputWS);
187 alg->setProperty("Mass", m_mass);
188 alg->execute();
189 MatrixWorkspace_sptr tofInY = alg->getProperty("OutputWorkspace");
190 if (!m_sumResults)
191 return tofInY;
192
193 // Rebin to common grid
194 double xmin(0.0), xmax(0.0);
195 tofInY->getXMinMax(xmin, xmax);
196 std::vector<double> params(3);
197 params[0] = xmin;
198 params[1] = SUMMEDY_BIN_WIDTH;
199 params[2] = xmax;
200
201 alg = createChildAlgorithm("Rebin", 0.05, 0.1, false);
202 alg->setProperty("InputWorkspace", tofInY);
203 alg->setProperty("Params", params);
204 alg->execute();
205 return alg->getProperty("OutputWorkspace");
206}
207
216double NormaliseByPeakArea::fitToMassPeak(const MatrixWorkspace_sptr &yspace, const size_t index) {
217 auto alg = createChildAlgorithm("Fit");
218 auto func = FunctionFactory::Instance().createFunction("ComptonPeakProfile");
219 func->setAttributeValue("Mass", m_mass);
220 func->setAttributeValue("WorkspaceIndex", static_cast<int>(index));
221
222 // starting guesses based on Hydrogen spectrum
223 func->setParameter("Position", PEAK_POS_GUESS);
224 func->setParameter("SigmaGauss", PEAK_WIDTH_GUESS);
225
226 // Guess at intensity
227 const size_t npts = yspace->blocksize();
228 const auto &yVals = yspace->y(index);
229 const auto &xVals = yspace->x(index);
230 double areaGuess(0.0);
231 for (size_t j = 1; j < npts; ++j) {
232 areaGuess += yVals[j - 1] * (xVals[j] - xVals[j - 1]);
233 }
234 func->setParameter("Intensity", areaGuess);
235 if (g_log.is(Logger::Priority::PRIO_DEBUG)) {
236 g_log.debug() << "Starting values for peak fit on spectrum " << yspace->getSpectrum(index).getSpectrumNo() << ":\n"
237 << "area=" << areaGuess << "\n"
238 << "width=" << PEAK_WIDTH_GUESS << "\n"
239 << "position=" << PEAK_POS_GUESS << "\n";
240 }
241 alg->setProperty("Function", func);
242 alg->setProperty("InputWorkspace", std::static_pointer_cast<Workspace>(yspace));
243 alg->setProperty("WorkspaceIndex", static_cast<int>(index));
244 alg->setProperty("CreateOutput", true);
245 alg->execute();
246
247 MatrixWorkspace_sptr fitOutputWS = alg->getProperty("OutputWorkspace");
248 saveToOutput(m_fittedWS, fitOutputWS->sharedY(1), yspace->sharedE(index), index);
249
250 double area = func->getParameter("Intensity");
251 if (g_log.is(Logger::Priority::PRIO_INFORMATION)) {
252 g_log.information() << "Calculated peak area for spectrum " << yspace->getSpectrum(index).getSpectrumNo() << ": "
253 << area << "\n";
254 }
255 return area;
256}
257
263void NormaliseByPeakArea::normaliseTOFData(const double area, const size_t index) {
264 m_normalisedWS->mutableY(index) = m_inputWS->y(index) / area;
265 m_normalisedWS->mutableE(index) = m_inputWS->e(index) / area;
266}
267
274void NormaliseByPeakArea::saveToOutput(const API::MatrixWorkspace_sptr &accumWS,
276 const Kernel::cow_ptr<HistogramData::HistogramE> &eValues, const size_t index) {
277 assert(yValues->rawData().size() == eValues->rawData().size());
278
279 if (m_sumResults) {
280 const size_t npts(accumWS->blocksize());
281 auto &accumY = accumWS->mutableY(0);
282 auto &accumE = accumWS->mutableE(0);
283 const auto &yValuesRaw = yValues->rawData();
284 const auto &eValuesRaw = eValues->rawData();
285
286 for (size_t j = 0; j < npts; ++j) {
287 double accumYj = accumWS->y(0)[j];
288 double accumEj = accumWS->e(0)[j];
289 double rhsYj = yValuesRaw[j];
290 double rhsEj = eValuesRaw[j];
291 if (accumEj < 1e-12 || rhsEj < 1e-12)
292 continue;
293 double err = 1.0 / (accumEj * accumEj) + 1.0 / (rhsEj * rhsEj);
294 accumY[j] = accumYj / (accumEj * accumEj) + rhsYj / (rhsEj * rhsEj);
295 accumY[j] /= err;
296 accumE[j] = 1.0 / sqrt(err);
297 }
298 } else {
299 accumWS->setSharedY(index, yValues);
300 accumWS->setSharedE(index, eValues);
301 }
302}
303
307void NormaliseByPeakArea::symmetriseYSpace() {
308 // A window is defined the around the Y value of each data point & every other
309 // point is
310 // then checked to see if it falls with in the absolute value +/- window
311 // width.
312 // If it does then the signal is added using the error as a weight, i.e
313 //
314 // yout(j) = yout(j)/(eout(j)^2) + y(j)/(e(j)^2)
315
316 // Symmetrise input data in Y-space
317 const double dy = 0.1;
318 const size_t npts(m_yspaceWS->blocksize());
319 const auto nhist = static_cast<int64_t>(m_symmetrisedWS->getNumberHistograms());
320
322 for (int64_t i = 0; i < nhist; ++i) {
323 const auto &xsym = m_symmetrisedWS->x(i);
324 auto &ySymOut = m_symmetrisedWS->mutableY(i);
325 auto &eSymOut = m_symmetrisedWS->mutableE(i);
326 const auto &yyspace = m_yspaceWS->y(i);
327 const auto &eyspace = m_yspaceWS->e(i);
328
329 for (size_t j = 0; j < npts; ++j) {
330 const double ein = eyspace[j];
331 const double absXj = fabs(xsym[j]);
332
333 double yout(0.0), eout(1e8);
334 for (size_t k = 0; k < npts; ++k) {
335 const double yk = yyspace[k];
336 const double ek = eyspace[k];
337 const double absXk = fabs(xsym[k]);
338 if (absXj >= (absXk - dy) && absXj <= (absXk + dy) && ein != 0.0) {
339
340 if (ein > 1e-12) {
341 double invE2 = 1 / (ek * ek);
342 yout /= eout * eout;
343 yout += yk * invE2;
344 double wt = (1 / (eout * eout)) + invE2;
345 yout /= wt;
346 eout = sqrt(1 / wt);
347 } else {
348 yout = 1e-12;
349 eout = 1e-12;
350 }
351 }
352 }
353 ySymOut[j] = yout;
354 eSymOut[j] = eout;
355 m_progress->report();
356 }
357 }
358}
359
360} // namespace Mantid::CurveFitting::Algorithms
#define DECLARE_ALGORITHM(classname)
Definition: Algorithm.h:576
IPeaksWorkspace_sptr workspace
Definition: IndexPeaks.cpp:114
std::map< DeltaEMode::Type, std::string > index
Definition: DeltaEMode.cpp:19
#define fabs(x)
Definition: Matrix.cpp:22
#define PARALLEL_FOR_IF(condition)
Empty definitions - to enable set your complier to enable openMP.
Base class from which all concrete algorithm classes should be derived.
Definition: Algorithm.h:85
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
virtual std::shared_ptr< Algorithm > createChildAlgorithm(const std::string &name, const double startProgress=-1., const double endProgress=-1., const bool enableLogging=true, const int &version=-1)
Create a Child Algorithm.
Definition: Algorithm.cpp:842
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.
A property class for workspaces.
A validator which checks that the unit of the workspace referred to by a WorkspaceProperty is the exp...
API::MatrixWorkspace_sptr convertInputToY()
Convert input workspace to Y coordinates for fitting.
void setUnitsToMomentum(const API::MatrixWorkspace_sptr &workspace)
Set the units meta-data.
void saveToOutput(const API::MatrixWorkspace_sptr &accumWS, const Kernel::cow_ptr< HistogramData::HistogramY > &yValues, const Kernel::cow_ptr< HistogramData::HistogramE > &eValues, const size_t index)
Stores/accumulates the results.
void createOutputWorkspaces(const API::MatrixWorkspace_sptr &yspaceIn)
Create the output workspaces.
API::MatrixWorkspace_sptr m_normalisedWS
Normalised output in TOF.
void normaliseTOFData(const double area, const size_t index)
Normalise given TOF spectrum.
API::MatrixWorkspace_sptr m_symmetrisedWS
Fitted output.
double fitToMassPeak(const API::MatrixWorkspace_sptr &yspace, const size_t index)
Fit the mass peak & find the area value.
void symmetriseYSpace()
Symmetrises the data in yspace about the origin.
API::MatrixWorkspace_sptr m_yspaceWS
Input data converted (and possible summed) to Y space.
void retrieveInputs()
Check and store appropriate input data.
bool m_sumResults
Flag to indicate if results are to be summed.
std::unique_ptr< API::Progress > m_progress
Reporting.
API::MatrixWorkspace_sptr m_fittedWS
Fitted output.
IPropertyManager * setProperty(const std::string &name, const T &value)
Templated method to set the value of a PropertyWithValue.
void debug(const std::string &msg)
Logs at debug level.
Definition: Logger.cpp:114
bool is(int level) const
Returns true if at least the given log level is set.
Definition: Logger.cpp:146
void information(const std::string &msg)
Logs at information level.
Definition: Logger.cpp:105
static T & Instance()
Return a reference to the Singleton instance, creating it if it does not already exist Creation is do...
Implements a copy on write data template.
Definition: cow_ptr.h:41
std::shared_ptr< MatrixWorkspace > MatrixWorkspace_sptr
shared pointer to the matrix workspace base class
double SUMMEDY_BIN_WIDTH
Bin width for rebinning workspace converted from TOF.
double PEAK_WIDTH_GUESS
Starting value of width in y-space for fit.
double PEAK_POS_GUESS
Starting value of peak position in y-space for fit.
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
@ Input
An input workspace.
Definition: Property.h:53
@ Output
An output workspace.
Definition: Property.h:54